基于核方法的高光譜遙感圖像混合像元分解

林娜， 楊武年， 王斌

(1.重慶交通大學土木工程學院測繪系,重慶 400074； 2.成都理工大學地學空間信息技術國土資源部重點實驗室,成都 610059； 3.重慶市地理信息中心,重慶 401121)

基于核方法的高光譜遙感圖像混合像元分解

林娜1,2， 楊武年2， 王斌3

(1.重慶交通大學土木工程學院測繪系,重慶 400074； 2.成都理工大學地學空間信息技術國土資源部重點實驗室,成都 610059； 3.重慶市地理信息中心,重慶 401121)

為了提高高光譜遙感圖像混合像元分解的精度，提出基于核方法的高光譜線性解混算法。采用正交子空間投影(orthogonal subspace projection，OSP)算子、最小二乘正交子空間投影(least squares OSP， LSOSP)算子、非負約束最小二乘(nonnegative constrained least-squares，NCLS)算子和全約束最小二乘(fully constrained least-squares，FCLS)算子等方法分別構建核正交子空間投影(kernel OSP，KOSP)、核最小二乘正交子空間投影(kernel LSOSP，KLSOSP)、核非負約束最小二乘(kernel NCLS，KNCLS)和核全約束最小二乘(kernel FCLS，KFCLS)高光譜圖像混合像元解混模型； 對美國內華達州CUPRITE礦區AVIRIS數據進行KLSOSP，KNCLS和KFCLS與LSOSP，NCLS和FCLS豐度反演對比實驗。結果表明： 對于混合像元廣泛存在的高光譜遙感圖像來說，基于核方法的KLSOSP，KNCLS和KFCLS的解混精度優于LSOSP，NCLS和FCLS，其中又以KFCLS解混的精度最高； 附加約束條件有利于提高豐度反演的精度。

高光譜遙感； 核方法； 混合像元分解； 正交子空間投影(OSP)

0 引言

高光譜遙感是遙感技術的重要發展方向之一[1]，其光譜分辨率高達nm數量級，可獲得地物的連續光滑的光譜曲線[2]，因此得到廣泛應用[3]。但因高光譜成像光譜儀側重于高光譜分辨率，其空間分辨率相對較低； 加之由于地面的復雜多樣性，使高光譜遙感圖像中廣泛存在混合像元[4]，制約了高光譜遙感技術的應用。高光譜解混是解決混合像元問題的一個有效方法，其研究內容主要是提取遙感圖像基本組成成分(端元)的特征光譜和有效估計混合像元中各個端元之間的混合比例(豐度)[5]。

線性光譜混合模型(linear spectral mixture model，LSMM)在混合像元分解中得到了廣泛應用，它假定混合像元中各端元之間的相互影響可以忽略不計，每一個觀測像元矢量都可以表述為端元和它們相應的豐度的線性組合[6]。LSMM的豐度反演主要應用最小二乘方法[7]； 為了提高反演的精度，可加入各種約束條件，根據不同約束條件可分為： 無約束最小二乘法(unconstrained least-squares，ULS)、非負約束最小二乘法(nonnegative constrained least-squares，NCLS)和全約束最小二乘法(fully constrained least-squares，FCLS)。近年來由于核函數方法具有強大的數據映射能力，可有效提高混合像元分解的精度，因而在高光譜遙感圖像解混中得到廣泛應用[8]。國外學者Camps-Valls等[9]、Hosseini等[10]和Broadwater等[11]以及國內學者劉婷婷等[12]、王曉飛等[13]和譚熊等[14]利用支持向量機(support vector machine，SVM)及其改進算法進行混合像元分解。Kwon等[15]在2005年首次提出核正交子空間投影(kernel orthogonal subspace projection，KOSP)算法； 此后，很多學者將該方法應用于圖像分類、目標檢測、混合像元分解等領域。Liu 等[16]在KOSP的基礎上發展了基于核方法的線性光譜分析。王挺等[17]考慮到局域信息，對KOSP方法進行改進，并應用到高光譜圖像目標探測領域。趙春暉等[18]則提出了一種基于自適應核方法的正交子空間投影高光譜圖像異常檢測算法。國內外學者利用核方法進行高光譜混合像元分解的研究集中在對SVM和KOSP算法的改進及應用上，但對核方法與LSMM模型的豐度反演結合的研究還比較少，尚不成體系。

本文在KOSP算法的基礎上，將核方法與LSMM模型的豐度反演方法相結合，采用最小二乘正交子空間投影(least squares orthogonal subspace projection，LSOSP)算子，構建核最小二乘正交子空間投影(kernel LSOSP，KLSOSP)高光譜混合像元解混方法； 然后用KLSOSP導出核非負約束最小二乘(kernel nonnegative constrained least-squares，KNCLS)，進一步導出核全約束最小二乘(kernel fully constrained least-squares，KFCLS)； 并進行KLSOSP,KNCLS和KFCLS與LSOSP,NCLS和FCLS的高光譜豐度反演對比實驗。

1 線性光譜分解模型

1.1 正交子空間投影(OSP)

線性光譜混合模型假設光譜向量r由一組圖像端元線性混合而成，即

r=Mα+n，

(1)

式中：r為混合像元光譜，r∈Rl，R為端元，l為波段數；M為端元矩陣，M∈Rl×p，p為端元數目；α為豐度向量，α∈Rp；n為模型誤差。

正交子空間投影(orthogonal subspace projection，OSP)方法在進行混合像元分解時，把M分為感興趣的目標端元d與不感興趣的背景端元U這2部分[19-20]。為不失一般性，設d=mp為第p個端元的光譜信號，U=[m1，m2，…，mp-1]為其他端元的光譜信號，則式(1)可改寫為

r=dαp+Uγ+n，

(2)

式中：αp為第p個端元的豐度；U為背景信號矩陣，U=[m1，m2，…，mp-1]；γ為U對應的端元豐度。根據最大化信噪比原則得到未知的豐度αp(即OSP探測算子)為

(3)

1.2 最小二乘正交子空間投影(LSOSP)

利用式(3)中的OSP探測算子求解端元豐度時，需要已知先驗的豐度矩陣，這通常是很難獲得的，但可利用最小二乘原理估計未知豐度αp，即

(4)

(5)

1.3 全約束最小二乘(FCLS)

2 基于核方法的高光譜豐度反演

通常核方法的應用是將所有的訓練樣本通過核函數映射到高維特征空間，再在特征空間進行線性分析，如核主成分分析[21]、核最小噪聲分離變換[22]、核RX探測[23]等。這類方法因所有的數據樣本都需要核化，計算復雜度較高。本文對核方法的應用是通過核化OSP，LSOSP，NCLS和FCLS算子來進行豐度反演，不僅可提高線性解混的精度，而且可大大降低計算復雜度。首先核化OSP探測算子，構建KOSP模型，利用最小二乘法求解導出KLSOSP； 然后加入非負約束條件導出KNCLS，加入全約束條件導出KFCLS。

2.1 KOSP與KLSOSP

(6)

其中

(7)

將式(7)帶入式(6)則得KOSP關系式，即

(8)

因[K(d，Iφ)，K(Iφ，r)]=K(d，r)，利用核方法將式(8)進一步簡化為

(9)

相似地，由式(5)得到

(10)

再將式(7)帶入式(10)，得到

(11)

2.2 核非負約束最小二乘(KNCLS)

因豐度值不能為負值，所有像元必須滿足豐度非負約束條件(ANC)，故可得到方程組

(12)

對于ANC約束下的最小二乘法也就變成了如下最優化問題，即

Min(r-Mα)T(r-Mα)αj≥0。

(13)

構造一個Lagrange乘數，則有

(14)

(15)

(16)

式(15)和(16)的核形式為

(17)

(18)

2.3 核全約束最小二乘(KFCLS)

(19)

(20)

式中λ和δ這2個參數可分別實現“豐度非負”和“豐度總和為1”這2個約束條件。

3 實驗與分析

3.1 實驗數據

本文選擇機載可見/近紅外成像光譜儀(airborne visible/infrared imaging spectrometer，AVIRIS)1997年6月19日獲取的美國內華達州CUPRITE礦區的成像光譜數據為實驗數據。AVIRIS數據是典型的高光譜數據，它覆蓋全反射光譜區域(0.4～2.5 μm)，共有224個波段，平均光譜分辨率為10 nm，地面分辨率約為20 m。從所獲取的高光譜遙感圖像數據中截取548像元×511像元的子區作為研究區數據，該地區礦物類型主要有明礬石、高嶺石、方解石、云母、玉髓和水銨長石等[25-26]。圖1為研究區高光譜圖像立方體及礦物蝕變分布圖[27]。

(a) 高光譜圖像立方體(b) 礦物蝕變分布圖

圖1 研究區高光譜圖像立方體及礦物蝕變分布圖

Fig.1 Maps of hyperspectral image cube and mineral alteration distribution

由于數據的空間分辨率較低，且地面礦物復雜多樣，故圖像中混合像元較多，適合于進行高光譜解混實驗。實驗之前先對其進行了大氣校正，選用短波紅外的50個連續波段(波段172(1.990 8 μm)―波段221(2.479 0 μm))進行實驗。

3.2 實驗方法及結果分析

實驗中的端元采用美國地質調查局(USGS)權威網站公布的CUPRITE礦區的參考端元光譜[28](圖2)，此端元比較準確，有利于后續對豐度反演方法的分析。

(為清晰表示反射率變化，對縱坐標進行了平移)

用MATLAB編程實現KLSOSP,KNCLS和KFCLS豐度反演算法，其中的核函數都選擇使用廣泛的徑向基函數(radial basis function,RBF)，其計算公式為

(21)

式中：xi和xj為不同的樣本；σ為函數的寬度。

核函數中對參數σ的選擇需要經過多次實驗才能確定較好的值，因KFCLS算法中通過λ和δ這2個參數實現ASC和ANC約束條件，δ是一個小值(本文中取1×10-5)，導致KFCLS中的參數σ與KLSOSP和KNCLS中的參數σ相比要小得多。最終，確定KLSOSP,KNCLS和KFCLS中的參數σ分別為2 000,2 000和0.5。

進行KLSOSP,KNCLS和KFCLS與LSOSP,NCLS和FCLS豐度反演對比實驗。對實驗結果從解混礦物豐度圖和解混均方根誤差圖像2個方面進行分析，二者互為補充。

3.2.1 解混礦物豐度圖分析

限于篇幅，僅列出明礬石和高嶺石2種礦物的豐度圖。圖3―5分別是用LSOSP和KLSOSP，NCLS和KNCLS，FCLS和KFCLS解混得到的明礬石和高嶺石豐度圖。豐度圖中色調越亮的部分，表明礦物在該位置的比例越大。

(a) LSOSP (b) KLSOSP

圖3 利用LSOSP和KLSOSP解混得到的明礬石和高嶺石豐度圖

Fig.3 Alunite and kaolinite abundance maps produced by using LSOSP and KLSOSP

(a) NCLS (b) KNCLS

圖4 利用NCLS和KNCLS解混得到的明礬石和高嶺石豐度圖

Fig.4 Alunite and kaolinite abundance maps produced by using NCLS and KNCLS

(a) FCLS (b) KFCLS

圖5 利用FCLS和KFCLS解混得到的明礬石和高嶺石豐度圖

Fig.5 Alunite and kaolinite abundance maps produced by using FCLS and KFCLS

從圖3―5可以清楚地看出，用KLSOSP，KNCLS和KFCLS方法解混得到的明礬石和高嶺石豐度圖較清晰，與礦物分布圖的吻合度也更好。這3種核化方法解混的效果明顯優于其對應的LSOSP，NCLS和FCLS方法。

對豐度圖進行統計分析的結果表明： LSOSP和KLSOSP解混得到的礦物豐度有小于0或大于1的情況，這與實際情況是不符的； 附加豐度非負約束條件之后，NCLS和KNCLS解混得到的礦物豐度都大于等于0； 附加全約束條件之后，FCLS和KFCLS的豐度都在0和1之間。由此可見，附加約束條件有利于提高礦物豐度反演的精度。

3.2.2 解混均方根誤差分析

均方根誤差(root mean square error，RMSE)計算公式為

RMSE=|∑(n2/N)|1/2,

(22)

式中N為像元總數。

圖6示出用KLSOSP，KNCLS和KFCLS解混之后的均方根誤差圖像。誤差圖像中的色調越亮，誤差就越大； 色調暗并且均勻分布時，便達到了RMSE盡可能小的要求。

(a) KLSOSP(b) KNCLS(c) KFCLS

圖6 解混后的均方根誤差圖像

Fig.6 Root mean square error images after pixel un-mixing

從圖6可以很明顯地看出，KFCLS的誤差圖像中亮像素最少，灰度分布最均勻； KLSOSP的誤差圖像中亮像素最多； KNCLS的誤差圖像中亮像素的數量介于KLSOSP和KFCLS之間。由此可以看出，KFCLS解混的精度最高，其次是KNCLS，而KLSOSP解混的精度最低。

表1列出KLSOSP,KNCLS和KFCLS與LSOSP, NCLS和FCLS解混的均方根誤差對比情況。可以看出，KLSOSP，KNCLS和KFCLS解混的RMSE的最大值與均值分別小于對應的LSOSP，NCLS和FCLS。這進一步說明了KLSOSP，KNCLS和KFCLS方法分別優于其對應的LSOSP，NCLS和FCLS方法，且附加約束條件有利于提高豐度反演的精度。

表1 不同方法解混均方根誤差

4 結論

本文利用核方法，分別核化OSP，LSOSP，NCLS和FCLS算子，構建KOSP，KLSOSP，KNCLS和KFCLS高光譜圖像混合像元分解模型。KLSOSP，KNCLS和KFCLS與LSOSP，NCLS和FCLS高光譜遙感圖像解混對比實驗的結果表明： 對于混合像元廣泛存在的高光譜遙感圖像來說，利用KLSOSP，KNCLS和KFCLS進行豐度反演的效果要優于對應的LSOSP，NCLS和FCLS方法，其中又以KFCLS解混的精度最高； 附加約束條件則有利于提高豐度反演的精度。

[1] 李德仁,童慶禧,李榮興,等.高分辨率對地觀測的若干前沿科學問題[J].中國科學:地球科學,2012,42(6):805-813. Li D R,Tong Q X,Li R X,et al.Current issues in high-resolution Earth observation technology[J].Science China Earth Sciences,2012,55(7):1043-4051.

[2] Chang C I.Hyperspectral Data Processing:Algorithm Design and Analysis[M].Hoboken,NJ,USA:John Wiley and Sons Inc,2013.

[3] Zhang L P,Huang X.Advanced processing techniques for remotely sensed imagery[J].Journal of Remote Sensing,2009,13(4):559-569.

[4] Rosin P L.Robust pixel unmixing[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(9):1978-1983.

[5] 普晗曄,王斌,張立明.基于單形體幾何的高光譜遙感圖像解混算法[J].中國科學:信息科學,2012,42(8):1019-1033. Pu H Y,Wang B,Zhang L M.Simplex geometry-based abundance estimation algorithm for hyperspectral unmixing[J].Scientia Sinica Informationis,2012,42(8):1019-1033.

[6] Wirasakti S,Zein R A,Mafazi F.Comparative study of land cover linear spectral mixture analysis(LSMA)model on multispectral and hyperspectral imagery[C]//34th Asian Conference on Remote Sensing.Bali,Indonesia:[s.n.],2013.

[7] 趙春暉,肖健鈺.一種利用互信息加權的最小二乘法豐度反演算法[J].沈陽大學學報:自然科學版,2014,26(1):45-49. Zhao C H,Xiao J Y.An abundance inversion algorithm based on mutual information weighted least square error[J].Journal of Shenyang University:Natural Science,2014,26(1):45-49.

[8] 唐曉燕,高昆,倪國強.高光譜圖像非線性解混方法的研究進展[J].遙感技術與應用,2013,28(4):731-738. Tang X Y,Gao K,Ni G Q.Nonlinear spectral unmixing of hyperspectral images[J].Remote Sensing Technology and Application,2013,28(4):731-738.

[9] Camps-Valls G,Bruzzone L.Kernel Methods for Remote Sensing Data Analysis[M].Chichester:John Wiley and Sons Ltd,2009.

[10]Hosseini S A,Ghassemian H.A new fast algorithm for multiclass hyperspectral image classification with SVM[J].International Journal of Remote Sensing,2011,32(23):8657-8683.

[11]Broadwater J,Banerjee A.Mapping intimate mixtures using an adaptive kernel-based technique[C]//Proceeding of the 3rd Workshop on Hyperspectral Image and Signal Processing:Evolution in Remote Sensing.Lisbon:IEEE,2011:1-4.

[12]劉婷婷,林琿,張良培,等.利用SVM相關反饋和語義挖掘的遙感影像檢索[J].武漢大學學報:信息科學版,2012,37(4):415-418. Liu T T,Lin H,Zhang L P,et al.SVM-relevance-feedback and semantic-extraction-based RS image retrieval[J].Geomatics and Information Science of Wuhan University,2012,37(4):415-418.

[13]王曉飛,張鈞萍,張曄.高光譜圖像混合像元分解算法[J].紅外與毫米波學報,2010,29(3):210-215,229. Wang X F,Zhang J P,Zhang Y.Unmixing algorithm of hyperspectral images[J].Journal of Infrared and Millimeter Waves,2010,29(3):210-215,229.

[14]譚熊,余旭初,張鵬強,等.基于多核支持向量機的高光譜影像非線性混合像元分解[J].光學精密工程,2014,22(7):1912-1920. Tan X,Yu X C,Zhang P Q,et al.Nonlinear mixed pixel decomposition of hyperspectral imagery based on multiple kernel SVM[J].Optics and Precision Engineering,2014,22(7):1912-1920.

[15]Kwon H,Nasrabadi N M.Kernel orthogonal subspace projection for hyperspectral signal classification[J].IEEE Transactions on Geoscience and Remote Sensing,2005,43(12):2952-2962.

[16]Liu K H,Wong E,Du E Y,et al.Kernel-based linear spectral mixture analysis[J].IEEE Geoscience and Remote Sensing Letters,2012,9(1):129-133.

[17]王挺,杜博,張良培.顧及局域信息的核化正交子空間投影目標探測方法[J].武漢大學學報:信息科學版,2013,38(2):200-203,239. Wang T,Du B,Zhang L P.A local information-based kernelized OSP method for target detection[J].Geomatics and Information Science of Wuhan University,2013,38(2):200-203,239.

[18]趙春暉,尤佳,李曉慧.基于自適應核方法的正交子空間投影異常檢測算法[J].黑龍江大學自然科學學報,2012,29(2):254-258,272. Zhao C H,You J,Li X H.An orthogonal subspace projection anomaly detection algorithm based on adaptive kernel method[J].Journal of Natural Science of Heilongjiang University,2012,29(2):254-258,272.

[19]Bajorski P.Analytical comparison of the matched filter and orthogonal subspace projection detectors for hyperspectral images[J].IEEE Transactions on Geoscience and Remote Sensing,2007,45(7):2394-2402.

[20]Capobianco L,Garzelli A,Camps-Valls G.Target detection with semisupervised kernel orthogonal subspace projection[J].IEEE Transactions on Geoscience and Remote Sensing,2009,47(11):3822-3833.

[21]Fauvel M,Chanussot J,Benediktsson J A.Kernel principal component analysis for the classification of hyperspectral remote sensing data over urban areas[J].EURASIP Journal on Advances in Signal Processing,2009,2009:783194.

[22]林娜,楊武年,王斌.高光譜遙感影像核最小噪聲分離變換特征提取[J].武漢大學學報:信息科學版,2013,38(8):988-992. Lin N,Yang W N,Wang B.Hyperspectral image feature extraction via kernel minimum noise fraction transform[J].Geomatics and Information Science of Wuhan University,2013,38(8):988-992.

[23]Molero J M,Garzón E M,García I,et al.Anomaly detection based on a parallel kernel RX algorithm for multicore platforms[J].Journal of Applied Remote Sensing,2012,6(1):061503.

[24]Broadwater J,Chellappa R,Banerjee A,et al.Kernel fully constrained least squares abundance estimates[C]//Proceedings of IEEE International Geoscience and Remote Sensing Symposium.Barcelona,Spain:IEEE,2007:4041-4044.

[25]林娜,楊武年,王斌.基于KMNF和BP神經網絡的高光譜遙感影像分類[J].計算機工程與設計,2013,34(8):2774-2777,2782. Lin N,Yang W N,Wang B.Hyperspectral image classification on KMNF and BP neural network[J].Computer Engineering and Design,2013,34(8):2774-2777,2782.

[26]林娜,楊武年,王斌.基于核最小噪聲分離變換的高光譜遙感影像多類SVM分類[J].計算機應用與軟件,2014,31(6):116-119. Lin N,Yang W N,Wang B.Multi-class SVM classification for hyperspectral remote sensing image based on kernel minimum noise fraction transform[J].Computer Applications and Software,2014,31(6):116-119.

[27]Swayze G A,Clark R N,Goetz A F H,et al.Mapping advanced argillic alteration at Cuprite,Nevada,using imaging spectroscopy[J].Economic Geology,2014,109(5):1179-1221.

[28]林娜,楊武年,王斌.基于FLAASH的AVIRIS高光譜影像大氣校正[J].地理空間信息,2013,11(4):49-50,54. Lin N,Yang W N,Wang B.Atmospheric correction of AVIRIS hyperspectral image based on FLAASH[J].Geospatial Information,2013,11(4):49-50,54.

(責任編輯： 劉心季)

Pixel un-mixing for hyperspectral remote sensing image based on kernel method

LIN Na1,2， YANG Wunian2， WANG Bin3

(1.SchoolofCivilEngineering，ChongqingJiaotongUniversity，Chongqing400074，China； 2.KeyLaboratoryofGeoscienceSpatialInformationTechnology，MinistryofLandandResources，ChengduUniversityofTechnology，Chengdu610059，China； 3.ChongqingGeomaticsCenter,Chongqing401121，China)

In order to improve the accuracy of hyperspectral pixel un-mixing，the authors proposed a kernel based pixel un-mixing method in this paper. By adopting orthogonal subspace projection(OSP)operator, least squares OSP(LSOSP)operator, nonnegative constrained least squares(NCLS)operator and fully constrained least squares(FCLS)operator respectively, the authors established kernel OSP(KOSP)，kernel LSOSP(KLSOSP)，kernel NCLS(KNCLS)and kernel FCLS(KFCLS)for hyperspectral imagery pixel un-mixing. The comparative experiments of abundance inversion by applying KLSOSP, KNCLS, KFCLS and LSOSP, NCLS, FCLS to CUPRITE AVIRIS data were carried out，and the results show that, for heavily mixed hyperspectral images, the pixel un-mixing accuracy of kernels based KLSOSP，KNCLS and KFCLS is higher than that of LSOSP, NCLS and FCLS. Meanwhile，the constraint conditions can improve the accuracy of abundance estimates.

hyperspectral remote sensing； kernel trick； pixel un-mixing； orthogonal subspace projection (OSP)

10.6046/gtzyyg.2017.01.03

林娜,楊武年,王斌.基于核方法的高光譜遙感圖像混合像元分解[J].國土資源遙感,2017,29(1)：14-20.(Lin N,Yang W N,Wang B.Pixel un-mixing for hyperspectral remote sensing image based on kernel method[J].Remote Sensing for Land and Resources,2017,29(1)：14-20.)

2015-07-29；

2015-09-20

重慶市教委科技項目“基于核方法的高光譜遙感影像非線性特征提取及混合像元非線性分解”(編號： KJ1400325)、測繪遙感信息工程國家重點實驗室開放基金項目“基于核方法的高光譜遙感影像特征提取及混合像元分解研究”(編號： 13R03)和重慶交通大學博士基金項目“基于核方法的高光譜遙感信息提取研究”(編號： 2012KJC2-011)共同資助。

林娜(1981-),女,博士,副教授,主要從事高光譜遙感圖像處理研究。Email： linnawb@126.com。

TP 751.1

A

1001-070X(2017)01-0014-07