非凸稀疏低秩約束的高光譜解混方法

孔繁鏘，卞陳鼎，李云松，郭文駿

(1. 南京航空航天大學 航天學院，江蘇 南京 210016; 2. 西安電子科技大學 綜合業務網理論及關鍵技術國家重點實驗室，陜西 西安 710071)

?

非凸稀疏低秩約束的高光譜解混方法

孔繁鏘1，卞陳鼎1，李云松2，郭文駿1

(1. 南京航空航天大學 航天學院，江蘇 南京 210016; 2. 西安電子科技大學 綜合業務網理論及關鍵技術國家重點實驗室，陜西 西安 710071)

針對高光譜混合像元的豐度矩陣具有行稀疏特性，提出一種非凸稀疏低秩約束的高光譜解混方法．首先，建立高光譜圖像非凸稀疏低秩約束模型，將豐度系數矩陣的非凸p范數作為稀疏約束，并將豐度系數矩陣奇異值的非凸p范數作為低秩約束；其次，構建聯合低秩性先驗與稀疏性先驗的非凸極小化模型，并提出求解的增廣拉格朗日交替極小化算法，將復合正則化問題分解成多個單一正則化問題，交替迭代求解．實驗仿真結果表明，該算法比貪婪算法和凸優化算法能獲得更高的解混精度，并且適用于信噪比較高的高光譜數據．

圖像處理;稀疏解混;稀疏表示;低秩;凸優化

高光譜遙感圖像是由幾十乃至數百個連續波段圖像組成的三維圖像數據，有豐富的地物空間特性和光譜特性，可定性、定量地分析和識別被測對象，在環境監測、土地調查、軍事探測以及海洋遙感等領域有著廣泛的應用．但受光學器件性能的限制和地物復雜多樣性的影響，使得高光譜遙感圖像單個像元中往往包含不止一種地物覆蓋類型，從而形成混合像元．混合像元在很大程度上影響了地物的識別和分類精度，因此，如何有效地解決混合像元分解問題，已成為高光譜遙感圖像定量分析的關鍵．

在基于線性混合模型的高光譜混合像元分解中，一般采用基于幾何的方法和基于統計的方法等[1]．近幾年來，隨著壓縮傳感和稀疏表示理論的快速發展，學者們將稀疏性約束加入到高光譜解混模型中，形成了第3種方法——稀疏解混方法[2]．稀疏解混方法主要有凸優化算法和貪婪算法等，在凸優化算法方面有SUnSAL(Sparse Unmixing by variable Splitting and Augmented Lagrangian)[3]、CLSUnSAL[4]和SUnSAL-TV[5]等，這些方法用L1范數來表示信號的稀疏性，能夠進行高效的求解，然而數據在低噪聲情況下，這些方法對噪聲是不魯棒的．在貪婪算法方面，其代表性算法有匹配追蹤算法(Matching Pursuit，MP)和正交基追蹤算法(Orthogonal Matching Pursuit，OMP)等[2]，然而匹配追蹤算法及正交基追蹤算法容易陷入局部最優問題．低秩表示方法[6-7]能在全局性的低秩約束下求解系數矩陣，來減小局部最優問題的影響．文獻[8]中提出一種低秩約束的高光譜數據解混方法，對豐度系數矩陣進行低秩約束，其低秩約束解混模型為凸優化模型，能夠進行高效的求解．但近幾年來，非凸優化理論[7，9-10]證明在某些條件下相較于凸優化方法，非凸優化方法可進一步表征系數矩陣的稀疏性，在稀疏信號重建中對噪聲具有更強的魯棒性．

基于以上分析，為了提高高光譜解混的精度，筆者提出了一種新的高光譜圖像解混方法，建立聯合低秩性先驗與稀疏性先驗的非凸極小化模型，并使用增廣拉格朗日乘子法對解混模型進行快速求解，提高求解精度．

1 聯合稀疏解混模型

線性高光譜混合模型假設混合像元為光譜特性穩定的純像元(端元)的線性組合．假設觀測值包含L個波段，那么線性光譜混合模型可以表示為

由于高光譜圖像中包含的端元數一般遠小于光譜庫中的端元數，而且高光譜數據中相鄰像素所含的物質相似，因此豐度系數矩陣中非零元素所在位置具有結構化特征，即行稀疏特性．由此文獻[4]中提出了聯合稀疏解混模型:

其中，Y∈RL×K，為高光譜圖像數據，包含L個波段，K個像元；X∈ RM×K，為豐度系數矩陣；N為噪聲矩陣．那么，聯合稀疏解混問題可以表示為如下優化問題:

其中，Xrow-0代表豐度矩陣X中非零元素的行數，·F表示Frobenius范數．

2 非凸稀疏低秩約束的高光譜解混方法

2.1 非凸稀疏低秩約束的高光譜解混模型

近幾年來非凸正則化理論[7，9-10]說明: 在一定條件下，非凸優化在理論和實際上優于凸優化，能進一步表征矩陣或向量的稀疏性，在信號重構過程中具有更強的抗噪性能．因此，針對式(3)的線性高光譜混合像元模型，提出一種非凸稀疏低秩約束的高光譜解混模型:

其中，非負約束項lR+(·)定義為

?i，j ．

通過重寫最優化問題式(5)，利用拉格朗日乘子法將約束項作為二次罰函數加入目標函數:

其中，μ為懲罰因子，V1=X，V2=X．式(6)將無約束問題轉化為多項無約束最優化問題，當懲罰因子選擇適當時，能夠較好地逼近最優解．首先求解參數X，然后利用X依次對每個單一正則項進行求解更新．假設當前為第 k+1 次迭代，則對X的更新是求解最優化問題:

然后，利用X(k+1)，將變量V1和V2的最優化問題分別表示為

由此，將求解式(6)的復合正則化問題分解為式(7)～(9)各個獨立的、求解相對簡單的單一正則化問題．

2.2 算法實現與步驟

f(σ)=f(σ(k))+f(σ(k))，σ-σ(k)=p

針對式(9)，可利用硬閾值法進行迭代更新:

對于D1和D2的迭代更新，由下式可得

3 實驗結果與分析

實驗1 模擬數據1是由50×50個混合像元構成的，分別包括5，7，9，11，13，15端元的6組高光譜模擬數據，豐度系數服從Dirichlet隨機分布，每個混合像元滿足非負性與和為一性．在生成的模擬圖像中加入零均值高斯白噪聲，信噪比分別為 20 dB，25 dB，30 dB，35 dB，40 dB，45 dB 和 50 dB．模擬數據實驗從不同端元數解混性能和抗噪聲性能兩個方面分別對以上算法的解混結果進行對比分析．

不同端元數時解混性能的結果如圖1所示，對 30 dB 高斯噪聲下不同端元數時模擬數據解混性能進行比較．在不同端元數的情況下，筆者提出的算法和SOMP算法性能最好，SUnSAL-TV算法的性能要高于CLSUnSAL算法的，但都低于筆者提出的算法的．筆者提出的算法性能要高于SOMP算法的，在 30 dB 高斯噪聲下不同端元數時平均提高了 4.042 dB．抗噪聲性能的結果如圖2所示，是在信噪比為 20～ 50 dB 的高斯噪聲下，9端元模擬數據1解混結果比較．與其他算法相比可以看出，筆者提出的算法的信號重構誤差值在大部分情況下都是最高的．與SOMP算法相比，筆者提出的算法的信號重構誤差值均高于SOMP算法的，在 20～ 50 dB 高斯噪聲下平均提高了 5.35 dB．因此，筆者提出的算法將聯合稀疏特性和低秩特性用于稀疏解混中，相比CLSUnSAL和SOMP算法能有效地提高解混精度．

圖1 30dB高斯噪聲下模擬數據1解混結果比較圖2 9端元時的模擬數據1解混結果比較

圖4 不同信噪比高斯噪聲的模擬數據2解混結果比較

實驗2 模擬數據2采用文獻[6]中的模擬數據實驗2中的數據，由 100× 100個混合像元構成，并由USGS光譜庫中隨機選取的9個光譜曲線線性組合而成，其中3個原始豐度圖像如圖3(a)中原始豐度圖像所示．在生成模擬圖像后，加入零均值高斯白噪聲，信噪比分別為 20 dB，25 dB，30 dB，35 dB，40 dB，45 dB 和 50 dB．

圖4是在信噪比為20～50 dB的高斯噪聲下，模擬數據2解混結果比較．可看出，隨著信噪比的提高，各種算法的信號重構誤差逐步提高．與其他算法相比，CLSUnSAL的信號重構誤差最低，筆者提出的算法和SOMP算法的信號重構誤差值最高．與SOMP算法相比，筆者提出的算法的信號重構誤差值均高于SOMP算法的，在 20～ 50 dB 噪聲下平均提高了 1.45 dB．圖3為 30 dB 高斯噪聲下各種算法的豐度重建圖像比較，從圖中可以看出，筆者提出的算法與SUnSAL-TV算法相比，其豐度重建圖像包含的噪聲點更多，但保留了更多的邊緣區域，更接近真實豐度圖像．與SOMP和CLSUnSAL算法相比，筆者提出的算法的豐度重建圖像噪聲點較少，并更多地保留了邊緣區域，具有更好的視覺效果．

4 總 結

針對高光譜數據矩陣的秩與所觀測場景的地物類別之間存在近似相等關系，即高光譜豐度系數矩陣存在低秩屬性，筆者提出一種非凸稀疏低秩約束的高光譜解混方法．該方法將非凸稀疏約束和非凸低秩約束融入稀疏解混模型中，并采用增廣拉格朗日交替極小化方法對解混模型進行快速求解．使用兩組高光譜模擬數據進行仿真實驗，用信號重構誤差評價指標做比較，筆者提出的算法在性能和魯棒性方面均優于貪婪算法和凸優化算法的．

[1] BIOUCAS-DIAS J M, PLAZA A, DOBIGEON N, et al. Hyperspectral Unmixing Overview: Geometrical, Statistical, and Sparse Regression-based Approaches[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(2): 354-379.

[2]IORDACHE M D, BIOUCAS-DIASJ M, PLAZA A. Sparse Unmixing of Hyperspectral Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(6): 2014-2039.

[3]AFONSO M V, BIOUCAS-DIAS J M, FIGUEIREDO M A. An Augmented Lagrangian Approach to the Constrained Optimization Formulation of Imaging Inverse Problems[J]. IEEE Transactions on Image Processing, 2011, 20(3): 681-695.

[4]IORDACHE M D, BIOUCAS-DIAS J M, PLAZA A. Collaborative Sparse Regression for Hyperspectral Unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(1): 341-354.

[5]IORDACHE M D, BIOUCAS-DIAS J M, PLAZA A. Total Variation Spatial Regularization for Sparse Hyperspectral Unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(11): 4484-4502.

[6]LIU G C, LIN Z C, YAN S C, et al. Robust Recovery of Subspace Structures by Low-rank Representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2013, 35(1): 171-184.

[7]MAJUMDAR A, WARD R K, ABOULNASR T. Non-convex Algorithm for Sparse and Low-rank Recovery: Application to Dynamic MRI Reconstruction[J]. Magnetic Resonance Imaging, 2013, 31(3): 448-455.

[8]QU Q, NASRABADI N M, TRAN T D. Abundance Estimation for Bilinear Mixture Models via Joint Sparse and Low-rank Representation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(7): 4404-4423.

[9]張文娟, 馮象初. 非凸低秩稀疏約束的圖像超像素分割方法[J]. 西安電子科技大學學報, 2013, 40(5): 86-91.

ZHANG Wenjuan, FENG Xiangchu. Image Super-pixels Segmentation Method Based on the Non-convex Low-rank and Sparse Constraint[J]. Journal of Xidian University, 2013, 40(5): 86-91.

[10]CHARTRAND R, STANEVA V. Restricted Isometry Properties and Nonconvex Compressive Sensing[J]. Inverse Problems, 2008, 24(3): 1-14.

[11]DONG W S, SHI G M, LI X, et al. Compressive Sensing via Nonlocal Low-rank Regularization[J]. IEEE IEEE Transactions on Image Processing, 2014, 23(8): 3618-3632.

[12]SHI Z W, TANG W, DUREN Z, et al. Subspace Matching Pursuit for Sparse Unmixing of Hyperspectral Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(6): 3256-3274.

[13]CLARK R N, SWAYZE G A, WISE R A, et al. USGS Digital Spectral Library Splib06a [EB/OL]. [2014-12-28]. http://speclab.cr.usgs.gov/spectral.lib06.

(編輯：郭 華)

Hyperspectral unmixing method based on the non-convex sparse and low-rank constraints

KONGFanqiang1，BIANChending1，LIYunsong2，GUOWenjun1

(1. College of Astronautics, Nanjing Univ. of Aeronautics and Astronautics, Nanjing 210016, China; 2. State Key Lab. of Integrated Service Networks, Xidian Univ., Xi’an 710071, China)

Aiming at the row sparse feature of the abundance matrix of hyperspectral mixed pixels, a hyperspectral unmixing method based on the non-convex sparse and low-rank constraints is presented. The non-convex sparse representation and non-convex low-rank models are first constructed, which take the non-convex p-norm of the abundance matrix as the sparse constraint and the non-convex p-norm of the singular values of the abundance matrix as the low-rank constraint. Then the low-rank prior and the sparse prior are jointly utilized to construct a non-convex minimization model. An augmented lagrange alternating minimization method is proposed to solve the unmixing model, the compound regularization problem is decomposed into multiple single regularization problems solved by the variable separation method. Experimental results demonstrate that the proposed method outperforms the greedy algorithm and the convex algorithms with a better spectral unmixing accuracy, and is suitable for high signal-to-noise ratio hyperspectral data.

image processing; sparse unmixing; sparse representation; low-rank; convex optimization

2015-11-09

時間：2016-04-01

國家自然科學基金資助項目(61401200，61201365)

孔繁鏘(1980-)，男，講師， E-mail: kongfq@nuaa.edu.cn．

http://www.cnki.net/kcms/detail/61.1076.tn.20160401.1622.040.html

10.3969/j.issn.1001-2400.2016.06.020

TP751.1

A

1001-2400(2016)06-0116-06