全變差稀疏約束深度非負(fù)矩陣分解高光譜遙感影像解混方法

趙文君，翟 晗，張洪艷

(1.武漢大學(xué) 測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079；2.中國(guó)地質(zhì)大學(xué)(武漢) 地理與信息工程學(xué)院，湖北 武漢 430074)

高光譜遙感結(jié)合了空間成像技術(shù)與光譜技術(shù)，能以納米級(jí)的光譜分辨率獲取地物近乎連續(xù)的光譜曲線，捕獲地物間的細(xì)微差異，為地物精細(xì)化分類提供可能[1]，促進(jìn)了遙感行業(yè)領(lǐng)域的進(jìn)步與發(fā)展[2]，在農(nóng)業(yè)、林業(yè)、食品安全、礦物勘探、環(huán)境監(jiān)測(cè)等諸多領(lǐng)域得到廣泛應(yīng)用[3-6]。

由于成像光譜儀各波段獲得的輻射能量有限，高光譜影像的空間分辨率通常較低[7]，再加上地物的復(fù)雜性，高光譜影像中存在大量混合像元，限制了高光譜遙感的行業(yè)應(yīng)用水平。對(duì)此，高光譜解混從亞像元的角度出發(fā)，通過(guò)光譜混合模型將影像分解成純凈物質(zhì)組分(端元)及其所占比例(豐度)來(lái)提高信息提取精度。

非負(fù)矩陣分解模型理論嚴(yán)謹(jǐn)，且易于與物理先驗(yàn)相結(jié)合，在高光譜解混中得到廣泛應(yīng)用。但由于目標(biāo)函數(shù)的非凸性，導(dǎo)致解空間不唯一，降低了解混精度。故研究者們?cè)谀繕?biāo)函數(shù)中加入約束來(lái)提高算法精度。文獻(xiàn)[8]將分段平滑度和稀疏性引入非負(fù)矩陣分解，同時(shí)保留了光譜數(shù)據(jù)的不連續(xù)性。文獻(xiàn)[9]考慮本征流形結(jié)構(gòu)，保持了原始影像和豐度圖之間的緊密聯(lián)系。文獻(xiàn)[10]受到深度學(xué)習(xí)的啟發(fā)，考慮了隱含信息的層次特征，解決了神經(jīng)網(wǎng)絡(luò)方法梯度擴(kuò)散的問(wèn)題。文獻(xiàn)[11]引入自適應(yīng)圖來(lái)規(guī)范多層非負(fù)矩陣分解模型。

由于地物分布的復(fù)雜性，高光譜傳感器接收到的信號(hào)通常是高度混合的。傳統(tǒng)深度非負(fù)矩陣分解忽略了地物光譜的實(shí)際混合過(guò)程。為此，本文從地物光譜的傳播出發(fā)，反向模擬光譜混合過(guò)程，提出面向端元矩陣的深度非負(fù)矩陣分解模型，并考慮端元的稀疏分布和空間分段平滑性，對(duì)豐度矩陣施加L1/2約束和全變差約束，提高了解混精度。

1 高光譜影像解混

1.1 線性光譜混合模型

線性光譜混合模型假設(shè)觀測(cè)到的像元光譜是各端元光譜按一定比例線性混合而成的，如圖1所示。線性混合模型可以用下式表示

X=AS+N

(1)

式中，X∈RB×N代表高光譜數(shù)據(jù)；B代表波段數(shù)；N代表像元數(shù)；A∈RB×P是端元光譜矩陣；P代表端元數(shù)；S∈RP×N是豐度矩陣；N∈RB×N代表噪聲。

圖1 線性光譜混合模型Figure 1. Explanation of linear spectral mixing model

由于地表地物分布復(fù)雜，到達(dá)高光譜傳感器的信號(hào)通常是電磁波多次反射的結(jié)果，故高光譜傳感器接收到的信號(hào)通常是鄰近地物反射信號(hào)的深度融合。以電磁波信號(hào)在到達(dá)傳感器前進(jìn)行L次線性混合為例，深度線性光譜混合模型可用圖2簡(jiǎn)要表示。

圖2 深度線性光譜混合模型Figure 2. Explanation of deep linear spectral mixing model

則電磁波經(jīng)過(guò)多次反射深度線性混合后得到的數(shù)據(jù)矩陣X如式(2)所示。

X=A1S1S2…SL+N

(2)

1.2 非負(fù)矩陣分解模型

非負(fù)矩陣分解[12](Nonnegative Matrix Factorization，NMF)將數(shù)據(jù)矩陣近似為兩個(gè)非負(fù)矩陣的乘積，其代價(jià)函數(shù)如下式所示

(3)

式中，‖·‖F(xiàn)代表Frobenius范數(shù)。

1.3 深度非負(fù)矩陣分解模型

為研究高光譜遙感影像中的深度空譜信息，現(xiàn)有的深度非負(fù)矩陣分解從數(shù)學(xué)理論的角度出發(fā)，對(duì)高光譜數(shù)據(jù)進(jìn)行深度分析，通過(guò)對(duì)分解過(guò)程中的系數(shù)矩陣(即豐度矩陣)進(jìn)行深度分解，來(lái)探索高光譜遙感影像中隱含信息的特征。然而，若從數(shù)學(xué)理論層面對(duì)高光譜遙感影像解混進(jìn)行分析，則會(huì)忽略高光譜遙感影像的成像過(guò)程及其所表征的物理意義。

從光譜混合的物理過(guò)程及其所表征的物理意義等方面考慮，由于地物光譜混合后的結(jié)果所表征的物理意義依舊為光譜，故對(duì)豐度矩陣進(jìn)行深度分解是不合理的。這是由于其所表示的物理意義并非光譜，而是端元所占各像元的比例值。因此，本文從地物光譜的物理傳播過(guò)程出發(fā)，對(duì)光譜深度線性混合的物理過(guò)程進(jìn)行反向模擬，提出了面向端元矩陣的深度非負(fù)矩陣分解模型。

2 面向端元矩陣的深度非負(fù)矩陣分解

2.1 面向端元矩陣的深度非負(fù)矩陣分解模型

由于達(dá)到高光譜傳感器的信號(hào)是電磁波在地表多次反射作用后的結(jié)果，故對(duì)原始數(shù)據(jù)進(jìn)行一次分解所得到的端元矩陣不夠純凈。針對(duì)該問(wèn)題，本文從地物光譜的傳播過(guò)程出發(fā)，對(duì)光譜深度線性混合過(guò)程進(jìn)行反向模擬，提出了一種面向端元矩陣的深度非負(fù)矩陣分解模型。通過(guò)對(duì)每一層分解所得的端元矩陣進(jìn)行分解，得到下一層的端元矩陣和豐度矩陣，直至得到最終的端元矩陣和豐度矩陣。

假定到達(dá)高空傳感器的光譜從地物反射出發(fā)共經(jīng)歷了L次線性混合，設(shè)深度非負(fù)矩陣分解的層數(shù)為L(zhǎng)，則高光譜遙感影像數(shù)據(jù)矩陣X的深度非負(fù)矩陣分解模型可表示為如圖3所示的形式。

圖3 深度非負(fù)矩陣分解模型示意圖Figure 3. Explanation of deep NMF model

在第1層中，高光譜遙感影像數(shù)據(jù)矩陣X被分解為較純凈的端元矩陣A1及其豐度矩陣S1。在第l層中，上一層的端元矩陣Al-1被分解為更純凈的端元矩陣Al及其豐度矩陣Sl，則最終對(duì)高光譜遙感影像數(shù)據(jù)進(jìn)行解混后所得到的純凈端元矩陣為A=AL，其豐度矩陣為S=SL…S2S1，則高光譜數(shù)據(jù)矩陣X可用下式表示

X=ALSL…S2S1+N

(4)

基于歐氏距離的代價(jià)函數(shù)可用式(5)進(jìn)行表示。

(5)

本文提出的面向端元矩陣的深度非負(fù)矩陣分解模型包括預(yù)訓(xùn)練和微調(diào)兩個(gè)階段，如圖4所示。預(yù)訓(xùn)練階段通過(guò)最小化式(3)中的代價(jià)函數(shù)，用經(jīng)典的非負(fù)矩陣分解對(duì)所有層進(jìn)行逐層優(yōu)化。微調(diào)階段通過(guò)式(5)將預(yù)訓(xùn)練階段所得的所有因素作為整體進(jìn)行微調(diào)，從而減少總重建誤差。

圖4 深度非負(fù)矩陣分解流程圖Figure 4. Illustration of deep NMF

2.2 基于全變差稀疏約束的深度非負(fù)矩陣分解模型

由于非負(fù)矩陣分解的代價(jià)函數(shù)非凸，存在大量極小值，因此有必要根據(jù)高光譜數(shù)據(jù)特點(diǎn)，對(duì)非負(fù)矩陣分解模型施加正則化約束。

研究表明，高光譜遙感影像的大多數(shù)像元僅由場(chǎng)景中某幾個(gè)端元混合而成，其包含的端元數(shù)目通常遠(yuǎn)低于光譜庫(kù)的端元總數(shù)，故豐度矩陣具有稀疏性。文獻(xiàn)[13]證明了L1/2稀疏約束優(yōu)于其它Lq正則化約束?；诖?，本文將L1/2稀疏約束作為豐度矩陣的稀疏先驗(yàn)約束引入模型，得到預(yù)訓(xùn)練階段和微調(diào)階段的代價(jià)函數(shù)分別如式(6)和式(7)所示。

(6)

(7)

式中，λ為豐度稀疏約束項(xiàng)的系數(shù)，用于調(diào)節(jié)豐度稀疏約束的程度；‖·‖1/2代表L1/2稀疏約束，其計(jì)算式為

(8)

式中，si，j是影像中第i個(gè)端元占第j個(gè)像元的豐度；|·|代表絕對(duì)值。

高光譜影像具有空間平滑性。研究表明，全變差正則化能很好地提取影像的空間一階鄰域信息，可增加豐度圖的分段平滑性并減少噪聲的負(fù)面影響[14-16]。基于此，本文采用各向異性全變差正則化約束[17]來(lái)加強(qiáng)豐度圖的分段平滑性，如下式所示。

(9)

綜上，本文提出的全變差稀疏約束深度非負(fù)矩陣分解算法(Total Variation and Sparsity Regularized Deep Nonnegative Matrix Factorization，TVSDNMF)的預(yù)訓(xùn)練階段和微調(diào)階段的代價(jià)函數(shù)為

(10)

α‖Sl‖TV

(11)

式中，α為豐度全變差正則化項(xiàng)的系數(shù)，用于控制豐度圖分段平滑程度。

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

3.1 實(shí)驗(yàn)設(shè)置

本文通過(guò)模擬實(shí)驗(yàn)和真實(shí)實(shí)驗(yàn)分別驗(yàn)證所提出的全變差稀疏約束的深度非負(fù)矩陣分解算法對(duì)高光譜遙感影像解混的有效性。本文選擇頂點(diǎn)成分分析-全約束最小二乘(Vertex Component Analysis-Fully Constrained Least Squares，VCA-FCLS)[18]、L1/2稀疏約束非負(fù)矩陣分解(L1/2Sparsity-constrained Nonnegative Matrix Factorization，L1/2-NMF)[13]、圖正則化L1/2稀疏約束非負(fù)矩陣分解(Graph Regularization of Nonnegative Matrix Factorization，GNMF)[9]、多層非負(fù)矩陣分解(Multilayer Nonnegative Matrix Factorization， MLNMF)[19]和全變差稀疏約束的深度非負(fù)矩陣分解(Sparsity-constrained Deep Nonnegative Matrix Factorization with Total Variation， SDNMF-TV)[10]5種高光譜解混算法作為對(duì)比，并采用光譜角距離和均方根誤差分別對(duì)光譜提取及豐度估計(jì)結(jié)果進(jìn)行定量評(píng)價(jià)。

(12)

(13)

3.2 模擬實(shí)驗(yàn)

本文采用的模擬數(shù)據(jù)集由美國(guó)地質(zhì)調(diào)查局光譜庫(kù)中的數(shù)據(jù)生成。該庫(kù)包含近500種典型礦物，有224個(gè)波段，均勻分布在0.4～2.5 μm之間。本文選擇光鹵石、胺基明礬石、黑云母和陽(yáng)起石的188個(gè)低噪聲波段進(jìn)行模擬實(shí)驗(yàn)，其光譜曲線如圖5所示。模擬數(shù)據(jù)集由48×48像元組成[20]。數(shù)據(jù)是使用線性光譜混合模型生成的，同時(shí)在各模擬像元中施加豐度和唯一約束。在模擬數(shù)據(jù)中，存在純凈區(qū)域和混合區(qū)域兩種不同的區(qū)域，第1行為純凈區(qū)域，僅由一種端元生成，第2～4行為混合區(qū)域，分別使用兩個(gè)端元、3個(gè)端元和4個(gè)端元混合生成。各區(qū)域以方形分布，其具體分布如圖6所示。

圖5 模擬實(shí)驗(yàn)端元光譜圖Figure 5. Spectral signatures of synthetic experiments

圖6 模擬實(shí)驗(yàn)豐度圖Figure 6. Abundance maps of synthetic experiments

表1給出了使用不同解混方法得到的模擬數(shù)據(jù)集平均光譜角距離值和均方根誤差的定量評(píng)價(jià)，其中加粗字體為每列最優(yōu)值。從表中可以看到，TVSDNMF的平均光譜角距離值遠(yuǎn)低于其他方法。相比于面向豐度矩陣的深度非負(fù)矩陣分解方法，平均光譜角距離減小了0.007 8，均方根誤差降低了0.006 8，證明了面向端元矩陣的深度非負(fù)矩陣分解方法的優(yōu)越性與合理性。

表1 不同算法模擬實(shí)驗(yàn)定量評(píng)價(jià)結(jié)果

圖7示出了使用不同解混方法提取到的端元光譜，其中實(shí)線為端元真實(shí)光譜，虛線為端元估計(jì)光譜?？梢钥吹?，除了MLNMF，其余各方法端元提取結(jié)果精度較高。相比于面向豐度矩陣的深度非負(fù)矩陣分解方法，面向端元矩陣的深度非負(fù)矩陣分解方法的吻合度更高，再次證明了所提出的面向端元矩陣的深度非負(fù)矩陣分解方法的優(yōu)越性。

(a)

圖8則顯示了使用不同解混方法提取到的端元1所對(duì)應(yīng)的豐度圖。可以看出，兩種深度非負(fù)矩陣分解的豐度圖的噪聲點(diǎn)相較其他算法更少。通過(guò)對(duì)兩種方法進(jìn)行比較，可以發(fā)現(xiàn)，所提出的TVSDNMF算法得到的豐度圖的背景噪聲更小。

(a) (b) (c)

3.3 真實(shí)實(shí)驗(yàn)

本文采用了Samson高光譜數(shù)據(jù)集進(jìn)行真實(shí)實(shí)驗(yàn)，進(jìn)一步驗(yàn)證了所提算法的有效性。Samson數(shù)據(jù)集包含了952×952個(gè)像元，每個(gè)像元有156個(gè)波段，波段范圍為0.401～0.889 μm，其光譜分辨率高達(dá)3.13 nm。由于原始影像范圍較大，這會(huì)使得計(jì)算時(shí)間過(guò)長(zhǎng)，故本文參考之前的研究選擇了(252，332)處95×95個(gè)像元的區(qū)域來(lái)減少計(jì)算量，因?yàn)榇藚^(qū)域的數(shù)據(jù)不會(huì)由于空白光譜波段或嚴(yán)重噪聲的影響而降低精度[21]。該選區(qū)影像中主要含有3個(gè)端元，分別是土壤、樹(shù)木和水，所用選區(qū)影像圖如圖9所示。

圖9 Samson高光譜數(shù)據(jù)集Figure 9. Samson dataset used in real data experiment

表2給出了使用不同解混方法得到的Samson數(shù)據(jù)集平均光譜角距離值和均方根誤差的定量評(píng)價(jià)， 其中加粗字體為每列最優(yōu)值。從表中可以看到，TVSDNMF方法的平均光譜角距離值遠(yuǎn)低于其他方法，相比于面向豐度矩陣的深度非負(fù)矩陣分解方法，其平均光譜角距離即均方根誤差均有所降低，這證明了所提方法的優(yōu)越性。

圖10示出了使用不同解混方法提取的端元光譜曲線，其中實(shí)線為端元真實(shí)光譜，虛線為端元估計(jì)光譜。從圖中可以看到，兩種深度非負(fù)矩陣分解方法提取到的端元估計(jì)光譜與端元真實(shí)光譜更為一致。相較于SDNMF-TV，TVSDNMF估計(jì)所得的水的光譜曲線與真實(shí)光譜更為貼合。

(a)

圖11～圖13顯示了使用不同解混方法提取到的不同端元所對(duì)應(yīng)的豐度圖?？梢钥闯?，整體而言，兩種深度非負(fù)矩陣分解方法SDNMF-TV與TVSDNMF估計(jì)所得的豐度圖與真實(shí)值更為接近。在土壤的豐度圖中，除了兩種深度非負(fù)矩陣分解方法，其余方法所得豐度值均較小。相較兩種深度非負(fù)矩陣分解方法，本文所提出的TVSDNMF豐度估計(jì)結(jié)果更為細(xì)致，證明了所提方法的優(yōu)越性。

表2 Samson數(shù)據(jù)集不同算法的定量評(píng)價(jià)結(jié)果

(a) (b) (c) (d) (e) (f) (g) 圖11 不同方法提取的Samson數(shù)據(jù)集端元“水”的豐度圖(a)地面真實(shí)值 (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMFFigure 11. Abundance map of water in Samson dataset extracted by different methods(a)Ground truth (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMF

(a) (b) (c) (d) (e) (f) (g) 圖12 不同方法提取的Samson數(shù)據(jù)集端元“樹(shù)”的豐度圖(a)地面真實(shí)值 (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMFFigure 12. Abundance map of tree of Samson dataset extracted by different methods(a)Ground truth (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMF

圖13 不同方法提取的Samson數(shù)據(jù)集端元“土壤”豐度圖(a)地面真實(shí)值 (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e) MLNMF (f)SDNMF-TV (g)TVSDNMFFigure.13 Abundance map of soil of Samson dataset extracted by different methods(a)Ground Truth(b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMF

4 結(jié)束語(yǔ)

本文針對(duì)現(xiàn)有深度非負(fù)矩陣分解方法未考慮地物光譜實(shí)際混合過(guò)程的問(wèn)題，從地物反射電磁波的傳播過(guò)程出發(fā)，對(duì)光譜深度線性混合過(guò)程進(jìn)行反向建模，提出了面向端元矩陣的深度非負(fù)矩陣分解模型。此外，本文引入L1/2稀疏約束和全變差正則化來(lái)增強(qiáng)豐度稀疏性和空間平滑性，進(jìn)一步提高了解混精度，建立了全變差稀疏約束的深度非負(fù)矩陣分解模型。與5種解混方法(VCA-FCLS、L1/2-NMF、GNMF、MLNMF和SDNMF-TV)的對(duì)比結(jié)果表明，本文所提出的全變差稀疏約束的深度非負(fù)矩陣分解方法能更好地提取端元，具有良好的性能。

今后的工作將集中于將本文模型擴(kuò)展到能夠應(yīng)對(duì)混合噪聲(稀疏噪聲和高斯噪聲)的情形，以提高模型對(duì)噪聲的魯棒性。