基于光譜解混的城市地物分類研究

黃作維 胡光偉 謝世雄

(湖南工業大學農牧業廢棄物資源化綜合利用湖南省重點實驗室， 株洲 412000)

0 引言

現代城市是人工環境與自然環境的綜合體，地物種類復雜多變，下墊面的組成成分復雜多樣。傳統的遙感數據分析方法對城市地物光譜特征和區分城市地物所需的光譜信息缺乏，難以反映復雜多樣的城市地物覆蓋類型,不利于精細分類[1]。高光譜遙感技術彌補了這一缺陷，極大地提升了對復雜地物的識別能力,從而能夠區分光譜特征非常相似的城市地物。高光譜遙感能獲取地面目標豐富的空間、輻射和光譜3方面的信息，提供了各類地物完整的連續光譜信息，光譜分辨率優于10 nm，使得準確、精細地解譯地物成為可能[2]。由于高光譜圖像探測距離遠，大氣傳輸過程中存在混合效應，成像空間范圍場景復雜，空間分辨率低以及遙感儀器本身的混合效應這些因素導致混合像元廣泛存在[3-4]?；旌舷裨纸鈫栴}一直以來是當代遙感應用中的一個技術難點。高光譜混合像元分解模型一般分為兩類，即線性光譜混合模型和非線性光譜混合模型，線性模型是迄今為止使用最多的一種模型，其突出優點是算法簡單，物理含義明確。高光譜端元提取是分析和處理高光譜數據的前提條件，BOARDMAN[5]提出了利用凸面幾何學分析的純像元指數提取端元算法，文獻[6-7]提出了N-FINDR算法，NASCIMENTO等[8]提出了頂點成分分析算法，文獻[9-10]提出正交子空間投影方法。目前地物分類中利用的大部分是光譜信息，而對其空間領域信息利用較少，充分利用光譜信息和空域信息的結合來提高分類精度是一個亟待解決的問題[11-13]。本文利用相鄰像素在空間上的相關性，提出一種優化的候選端元判斷方法，運用改進的算法，在保證算法效率的同時可有效提高光譜解混的精度。

1 改進的混合像元分解方法

遙感傳感器在瞬時視場IFOV內搜集目標，對地物的探測是利用光子探測器或熱探測器檢測地物所對應的地表物質光譜信號的綜合，線性光譜解混的分解流程如圖1所示,主要有：數據預處理、端元提取、精度評價、豐度估計。線性光譜解混首先進行端元提取，然后進行混合像元分解[14]。端元是組成混合像元的最基本單位,物理意義是代表某種具有相對固定光譜特征的地物類型，端元提取是混合像元分解的難點[15-17]。通過對高光譜遙感數據分析表明：成像外界因素(噪聲)會引起地物光譜特征發生變化，且大多服從多元正態分布，關鍵圖像中的像素是通過光譜和空間特征來共同描述的,高光譜遙感圖像在空間分布上具有一定的連續性，在光譜特征上混合光譜(含有噪聲)與端元光譜具有相關性(在一定的近鄰范圍內端元與混合光譜具有相似性)[18-19]。實驗表明從物理和數學兩個方面的像元空間信息相關參數進行候選端元判斷更具有可靠性。

圖1 混合像元分解Fig.1 Processing of mixed pixel decomposition

1.1 端元優化方法

端元在物理意義上必須是最或然端元，一定位置的端元與其鄰域內的端元屬于相同地物類別的概率比較大，而噪聲則與其相鄰像元光譜不具相似性。在數學上端元矢量矩陣必須是一個非病態矩陣[20-21]?？梢杂霉庾V角來描述像元光譜矢量間的相似性。對混合像元解混中存在的“噪聲”，利用目標點鄰域空間信息的進一步處理可以對一些誤分結果起到修正的作用，有效提高分類精度[22-23]。

(1)空間領域方法：假設用M表示高光譜端元矩陣，圖像在空間坐標(i,j)處地物位置用M(i,j)來表示。假設以像元為中心，鄰域窗口大小為9×9，光譜波段編號為i，光譜波段數為n，對結果矩陣向四周擴張，它的鄰域像元表示為

(1)

(2)

則以(i，j)為中心的像元點鄰域信息表達為

得到鄰域像元的端元集合為

Gi,j={Mi-1,j-1,Mi-1,j,Mi-1,j+1,Mi,j-1,

Mi,j+1,Mi+1,j-1,Mi+1,j,Mi+1,j+1}

(3)

高光譜像元相似性評價：一般可以用光譜角余弦、圖像處理的窗口尺寸以及相似光譜百分數等參數來衡量，在地學統計分析方法中，可以用變差函數來描述某一波段光譜亮度的空間自相關性，變差函數為

(4)

式中h——像元之間的歐氏距離

N(h)——距離h之間的像元對個數

位置為xi的像元在j波段的光譜角為dnj(·)，位置為xi+h的像元在k波段的光譜角為dnk(·)。

(2)為充分利用混合光譜的近鄰特征，區分端元光譜和噪聲的敏感度，用K1表示像元相似性閾值，根據M中每個窗口的像元與中心像元的光譜角余弦進行排序，然后對每個窗口中光譜角余弦取中值，最后對窗口所有中值求平均值，作為像元相似性閾值

(5)

式中k——窗口中光譜編號

C——平均值D——中值

S——排序ri——光譜角余弦

優化準則: 根據高光譜像元分解后誤差最大(空間差異性最大)的像元如果是非噪聲像元，則其是候選端元。在端元集合Gi，j中，假設候選端元與M中各個端元的夾角余弦最小值為Kmin，利用光譜角閾值判斷最大誤差的像元矢量是否為候選端元方法是：根據閾值K1,若K1

1.2 改進算法

混合像元分解過程中不可避免地涉及大量矩陣運算和迭代過程，采用單機串行方法處理消耗大量時間，限制了像元分解的效率。如何快速地進行海量數據運算是目前遙感技術發展中亟待解決的問題。隨著網絡傳輸性能的提高和圖形處理器(Graphic processing unit，GPU)的快速發展，針對高光譜數據計算密集，數據處理過程復雜的特點，GPU能夠充分利用NVIDIA的圖形處理器的并行計算引擎，由于GPU有強大的浮點計算能力，采用粗粒度的線程塊間的并行計算和細粒度的線程間的雙層并行計算，隨著CPU技術的快速發展，目前對稱多處理器(Symmetrical multi-processing，SMP)已經成為最流行的并行計算系統之一。并行計算中的GPU 函數稱為內核函數，它可以分為兩個方面的并行，就是每個Grid 里有多個塊(block)并行，每個block又有多個線程(thread)并行。

在SMP集群的并行結構中，具有分布式內存和共享內存兩種并行計算的優勢，基于SMP集群下的并行優勢，提出了一種改進的端元提取算法。在SMP處理器共享內存模型中，為了提高ICA計算過程中每個并行處理單元的執行效率，算法中把每個處理單元稱為線程單元(Thread element，TE)。在實際集群運算中，并行能力能夠延伸到各計算節點之間(計算節點產生進程進行節點間的并行計算)，這種實體被稱為處理單元(Process element，PE)，為保證運算速度，采用多個TE進行并行計算。

改進的端元提取算法如下：

(1)從高光譜圖像中隨機抽取像元組成初始端元矩陣，利用PCA變換對圖像進行去相關處理，得到新圖像Y，并進行數據中心化與白化處理，處理的變換矩陣如下

(6)

其中D=diag(d1,d2,…,dn)E=(e1,e2,…,en)

式中D——協方差矩陣特征值作為對角元素的對角矩陣

E——協方差矩陣單位特征向量為列的矩陣

(2)對矩陣W進行初始化分離，假設一共有a個TE，基于劃分原則把圖像劃分為{Y1,Y2, …,Ya}，第j個TE負責對Yj(j=1,2,…,a)的劃分，第j個TE計算公式為

W(j)(k+1)=[E{Yjg(WTX)}-E{g′(WTYj)}W(k)]

(7)

(3)利用式(7)對每一個TE進行計算，當所有TE完成運算后，對W(k+1)進行規約為

(8)

(4)根據以上優化的候選端元判斷方法，在端元集合Gi，j中選取候選端元，完成后得到第j個TE(j=1,2，…,a)的豐度估計值為

(9)

(5)對SMP集群的并行結構中的同構的計算節點(PE)進行統計，對主PE圖像進行粗劃分，假設一共有b個PE，每個PE含有a個TE，對每個PE進行圖像劃分Y={Y1，Y2，…，Yb}，然后執行步驟(2)、(3)，把第p(p=1, 2，…,b)個從PE細劃分為Yp={Yp,1,Yp, 2, …,Yp, a}，運算后得到W(k+1)，并把結果傳輸給主PE。

(6) 對主PE進行規約得到W(k+1)，進行歸一化和對稱正交化后傳送到各從PE，迭代完成后，從PE獲得分離矩陣，根據各從PE豐度最大絕對值進行最大值規約，同時得到所對應的端元ei, j，其中i=1,2，…,m;j=1,2，…,a。

(7)當各TE搜索完候選端元后，根據各端元的豐度最大絕對值進行規約，主PE同步等待各從PE得到的端元提取結果，完畢后獲得最終端元集合，算法結束。

2 實驗

2.1 實驗數據與處理

研究區位于四川省成都市郊區，高光譜數據采集時間是2014年5月，傳感器采用上海技術物理研究所研制的PHI。在數據獲取期間，對實驗區進行了詳細的地面調查，選擇了具有代表性的1塊數據進行實驗，其中主要地物類型為道路、水體、建筑物、植被等，實驗中選擇了0.38～0.96 μm的112個波段，該數據的第3～6、23～34、66、77～89、99～108波段由于為水氣體吸收波段和信噪比低而被剔除，其余的72個波段被用于進一步處理，圖像大小為223 MB。經過圖像配準校正和輻射校正后，選取第7、45、90個波段分別作為R、G、B進行數據合成，合成后的假彩色圖像如圖2所示。

圖2 研究區影像圖Fig.2 Image of study area

為了定量評價算法的性能，在分解之前要獲得參考地物真實分布情況的各端元的參考光譜，首先從高光譜圖像上用 ROI 方法挑選出一片純凈地物(代表一種地物的均勻的區域)，并求出所有像元的均值，然后以此均值為標準建立起來光譜庫，該實驗的參考光譜是根據圖像上的 ROI 區域按照上述方法選取的。在生成過程中一共有92個較純的像元被選取出來，其中對應植被、水體、道路、建筑的像元數分別為20、12、28、32，取其均值作為該端元的參考光譜。

為了進一步定量評價算法，采用經典算法N-FINDR、ICA 和 OSP等方法在同條件下進行分析。實驗運行環境：計算機操作系統為 Windows 7(R)操作系統，CPU型號是intel core(TM)i7，內存4 GB，主頻為2.80 GHz，GPU 程序開發環境為 Microsoft Visual Studio.net 2008，并行處理的編程軟件采MPICH2以及Microsoft Visual C++2008 OpenMP并行程序設計，能有效降低并行計算的難度和復雜度。

2.2 端元提取實驗

根據優化端元提取方法，利用優化后的改進算法(I-ME)，參數的選取與仿真數據相同，不同地物種類植被、水體、道路、建筑物分別用數字編號 1、2、3、4 表示。圖3給出了在不同波段下4種端元的光譜曲線比較結果。從圖3可以看出，端元各不相同, 另外還采用和參考光譜相比較的方法，提取結果和參考光譜之間的光譜角余弦都較大，各提取端元與參考光譜的比較如圖4所示。

圖3 端元提取結果Fig.3 Result of endmember extraction

為了定量評價所提取的端元的準確度，本實驗利用光譜角距離(SAD)和光譜散度(SID)參數來衡量，表1給出了所提取端元的精度，結果表明，相比其他提取方法該算法(I-ME)有更好的精度，算法的精度由高到低依次為：I-ME、ICA、OSP、N-FINDR。

針對高光譜數據量大的特點，系統運用GPU高性能并行改進算法，隨著計算節點的增加，各算法的運行時間逐步減少。圖5為幾種不同方法的加速比曲線。由圖5可知，采用并行算法能獲得更高的效率，改進算法的加速比高于其他的方法(ICA、N-FINDR)，比N-FINDR提高了8.5%，充分說明了改進算法的有效性。

2.3 地物分類實驗及分析

參考數據是指真實地物圖，為考察準確性而選擇的像素，本次精度的檢驗數據是利用了野外調繪的真實數據，共包含4個真實地物類別，標定樣本13 219個，隨機選擇標定樣本的5%～20%作為訓練樣本，剩余作為測試樣本，利用地面真實數據與分類結果圖進行比較處理。當遙感圖像有較多的混合像元時，這種分類方法所得結果的精度(特別是數量精度)會下降?；煜仃囍袛祿碓从校阂阎牡湫蛥^域的地物類型圖，類前選擇的訓練區和訓練樣本時確定的各個類別及其空間分布圖和實地調查的分類結果圖。本研究采用基I-ME算法進行混合像元分解后再分類，利用 ENVI 軟件的Confusion Matrix 功能對上述的分類結果圖進行精度評價，精度結果是可靠的。

圖4 提取端元與參考光譜的比較Fig.4 Comparison result of extracted endmember and reference spectra

E1(SAD/SID)E2(SAD/SID)E3(SAD/SID)E4(SAD/SID)N-FINDER0.0921/0.03310.1364/0.05360.0862/0.00690.2313/0.0479ICA0.0903/0.03260.1371/0.04980.0812/0.00640.2325/0.0475OSP0.0910/0.03340.1358/0.05080.0834/0.00670.2337/0.0438I-ME0.0876/0.03010.1204/0.04310.0778/0.00680.2066/0.0416

圖5 并行算法的加速比Fig.5 Speedup ratio of parallel algorithm

針對傳統的地物分類方法(MLC、SVM)，一個重要環節就是特征提取和選擇，即為分類尋找最優特征，實現最大可分性。為了表明改進的混合像元解混方法對地物分類精度的影響，將原始圖像進行降維采樣后，進行對比實驗分析。第1組實驗是高光譜圖像進行預處理后直接利用傳統的方法(SVM、MLC)進行分類；第2組實驗是先將原始圖像進行預處理后利用本文的改進方法進行混合像元分解，再用SVM進行分類(I-ME)；第3組實驗是將原始圖像進行預處理后利用N-FINDR方法進行混合像元分解后，再利用分類器SVM進行分類。通過比較評價指標來間接評價解混的效果,采用了混淆矩陣、Kappa 系數、總體分類精度、識別精度、虛警率等指標對結果進行評價。

表2～4分別給出了基于SVM分類結果的混淆矩陣、基于N-FINDER算法的混合像元分解后SVM分類的混淆矩陣、基于改進算法的混合像元分解后SVM分類的混淆矩陣。表5～7給出了經過SVM分類后各類地物分類精度、基于N-FINDR混合像元分解后SVM分類的各類地物分類精度、基于改進算法混合像元分解后SVM的分類精度?？梢钥闯?，經過解混后(圖像空間信息特征的引入)可以改善分類精度。表5給出了用傳統方法(SVM)的分類結果，與解混后的表6、7對比，從Kappa系數、虛警率和識別精度來看，解混后的精度有一定提高。

圖6 不同分類方法的對比Fig.6 Comparison result of different classification methods

建筑物水體植被道路總和建筑物3653405333731水體 1102382831622737植被 58815638355811道路 1168130668937總和 38322571585690113216

表3 基于N-FINDR光譜解混后SVM分類的混淆矩陣Tab.3 Confusion matrix after spectral unmixing based on N-FINDR

表4 基于改進算法光譜解混后SVM分類的混淆矩陣Tab.4 Confusion matrix based on improved spectral unmixing SVM classification

表5 SVM分類后各類地物分類精度Tab.5 Classification accuracy based on SVM classification %

基于改進算法分解后SVM分類結果與基于N-FINDER算法的混合像元分解后相比，除了道路的Kappa系數從79.34%降到了78.53%，基于改進算法比基于N-FINDER算法的分類指標有較大改善，特別是虛警率指標有明顯降低。

表6 解混后各類地物分類精度(N-FINDR)Tab.6 Classification accuracy after spectral unmixing (N-FINDR) %

表7 解混后各類地物分類精度 (改進算法)Tab.7 Classification accuracy after spectral unmixing (improved method) %

實驗定量分析如圖6所示，圖6a為經典SVM分類的結果圖，圖6b為N-FINDR混合像元分解后SVM分類的結果圖，圖6c為改進算法混合像元分解后SVM分類的結果?？梢钥闯龈倪M算法能有效減少端元“噪聲”像元的數目，與基于N-FINDER算法的混合像元分解相比，改善了地物分類效果，可以看出絕大部分的地物(如建筑物、植被、道路、水體)都能被正確地分類，與實際情況大體相符，證明該分類方法的精度較高。但是由于部分地物之間復雜的空間拓撲關系以及其他不可預見的因素也導致了一些誤分的情況，如在植被和裸土相連的部分，裸土上的植被非常少，各種特征不明顯而被錯誤的分類為裸土，同時一部分裸土上有少量的植物被錯分為植被，又如當兩個中間有植被的建筑物相距很近，而且建筑物和植被高度非常接近時，兩者很難分類，但總體的分類精度還是比較理想的。

3 結束語

高光譜圖像分類是高光譜遙感應用領域的基礎，也是遙感圖像處理的核心之一。光譜解混是實現高光譜圖像高精度分類的前提，準確地分解混合像元是高光譜遙感技術得以廣泛應用的關鍵性問題。本文通過優化的候選端元判斷方法提出了一種改進的混合像元分解方法，相關實驗表明解混后得到的分類精度得到了明顯改善，證明了改進算法對混合像元分解的有效性。今后還要進一步改進對各種特征因子的空間分析性能，把CPU+GPU異構系統并行優化方法應用于遙感圖像處理應用中,提取更優的端元來對數據進行描述，從而獲得更高的分類精度。