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

(1. 南京航空航天大學(xué) 航天學(xué)院，江蘇 南京 210016; 2. 西安電子科技大學(xué) 綜合業(yè)務(wù)網(wǎng)理論及關(guān)鍵技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710071)

?

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

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

(1. 南京航空航天大學(xué) 航天學(xué)院，江蘇 南京 210016; 2. 西安電子科技大學(xué) 綜合業(yè)務(wù)網(wǎng)理論及關(guān)鍵技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710071)

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

圖像處理;稀疏解混;稀疏表示;低秩;凸優(yōu)化

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

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

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

1 聯(lián)合稀疏解混模型

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

由于高光譜圖像中包含的端元數(shù)一般遠(yuǎn)小于光譜庫(kù)中的端元數(shù)，而且高光譜數(shù)據(jù)中相鄰像素所含的物質(zhì)相似，因此豐度系數(shù)矩陣中非零元素所在位置具有結(jié)構(gòu)化特征，即行稀疏特性．由此文獻(xiàn)[4]中提出了聯(lián)合稀疏解混模型:

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

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

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

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

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

其中，非負(fù)約束項(xiàng)lR+(·)定義為

?i，j ．

通過重寫最優(yōu)化問題式(5)，利用拉格朗日乘子法將約束項(xiàng)作為二次罰函數(shù)加入目標(biāo)函數(shù):

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

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

由此，將求解式(6)的復(fù)合正則化問題分解為式(7)～(9)各個(gè)獨(dú)立的、求解相對(duì)簡(jiǎn)單的單一正則化問題．

2.2 算法實(shí)現(xiàn)與步驟

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

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

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

3 實(shí)驗(yàn)結(jié)果與分析

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

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

圖1 30dB高斯噪聲下模擬數(shù)據(jù)1解混結(jié)果比較圖2 9端元時(shí)的模擬數(shù)據(jù)1解混結(jié)果比較

圖4 不同信噪比高斯噪聲的模擬數(shù)據(jù)2解混結(jié)果比較

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

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

4 總 結(jié)

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

[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]. 西安電子科技大學(xué)學(xué)報(bào), 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

時(shí)間：2016-04-01

國(guó)家自然科學(xué)基金資助項(xiàng)目(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