黃作維 胡光偉 謝世雄

(湖南工業(yè)大學(xué)農(nóng)牧業(yè)廢棄物資源化綜合利用湖南省重點(diǎn)實(shí)驗(yàn)室， 株洲 412000)

0 引言

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

1 改進(jìn)的混合像元分解方法

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

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

1.1 端元優(yōu)化方法

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

(1)空間領(lǐng)域方法：假設(shè)用M表示高光譜端元矩陣，圖像在空間坐標(biāo)(i,j)處地物位置用M(i,j)來表示。假設(shè)以像元為中心，鄰域窗口大小為9×9，光譜波段編號(hào)為i，光譜波段數(shù)為n，對(duì)結(jié)果矩陣向四周擴(kuò)張，它的鄰域像元表示為

(1)

(2)

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

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

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)

高光譜像元相似性評(píng)價(jià)：一般可以用光譜角余弦、圖像處理的窗口尺寸以及相似光譜百分?jǐn)?shù)等參數(shù)來衡量，在地學(xué)統(tǒng)計(jì)分析方法中，可以用變差函數(shù)來描述某一波段光譜亮度的空間自相關(guān)性，變差函數(shù)為

(4)

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

N(h)——距離h之間的像元對(duì)個(gè)數(shù)

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

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

(5)

式中k——窗口中光譜編號(hào)

C——平均值D——中值

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

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

1.2 改進(jìn)算法

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

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

改進(jìn)的端元提取算法如下：

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

(6)

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

式中D——協(xié)方差矩陣特征值作為對(duì)角元素的對(duì)角矩陣

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

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

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

(7)

(3)利用式(7)對(duì)每一個(gè)TE進(jìn)行計(jì)算，當(dāng)所有TE完成運(yùn)算后，對(duì)W(k+1)進(jìn)行規(guī)約為

(8)

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

(9)

(5)對(duì)SMP集群的并行結(jié)構(gòu)中的同構(gòu)的計(jì)算節(jié)點(diǎn)(PE)進(jìn)行統(tǒng)計(jì)，對(duì)主PE圖像進(jìn)行粗劃分，假設(shè)一共有b個(gè)PE，每個(gè)PE含有a個(gè)TE，對(duì)每個(gè)PE進(jìn)行圖像劃分Y={Y1，Y2，…，Yb}，然后執(zhí)行步驟(2)、(3)，把第p(p=1, 2，…,b)個(gè)從PE細(xì)劃分為Yp={Yp,1,Yp, 2, …,Yp, a}，運(yùn)算后得到W(k+1)，并把結(jié)果傳輸給主PE。

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

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

2 實(shí)驗(yàn)

2.1 實(shí)驗(yàn)數(shù)據(jù)與處理

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

圖2 研究區(qū)影像圖Fig.2 Image of study area

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

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

2.2 端元提取實(shí)驗(yàn)

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

圖3 端元提取結(jié)果Fig.3 Result of endmember extraction

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

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

2.3 地物分類實(shí)驗(yàn)及分析

參考數(shù)據(jù)是指真實(shí)地物圖，為考察準(zhǔn)確性而選擇的像素，本次精度的檢驗(yàn)數(shù)據(jù)是利用了野外調(diào)繪的真實(shí)數(shù)據(jù)，共包含4個(gè)真實(shí)地物類別，標(biāo)定樣本13 219個(gè)，隨機(jī)選擇標(biāo)定樣本的5%～20%作為訓(xùn)練樣本，剩余作為測(cè)試樣本，利用地面真實(shí)數(shù)據(jù)與分類結(jié)果圖進(jìn)行比較處理。當(dāng)遙感圖像有較多的混合像元時(shí)，這種分類方法所得結(jié)果的精度(特別是數(shù)量精度)會(huì)下降?；煜仃囍袛?shù)據(jù)來源有：已知的典型區(qū)域的地物類型圖，類前選擇的訓(xùn)練區(qū)和訓(xùn)練樣本時(shí)確定的各個(gè)類別及其空間分布圖和實(shí)地調(diào)查的分類結(jié)果圖。本研究采用基I-ME算法進(jìn)行混合像元分解后再分類，利用 ENVI 軟件的Confusion Matrix 功能對(duì)上述的分類結(jié)果圖進(jìn)行精度評(píng)價(jià)，精度結(jié)果是可靠的。

圖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

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

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

圖6 不同分類方法的對(duì)比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 基于改進(jìn)算法光譜解混后SVM分類的混淆矩陣Tab.4 Confusion matrix based on improved spectral unmixing SVM classification

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

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

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

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

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

3 結(jié)束語

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