趙文君，翟 晗，張洪艷

(1.武漢大學(xué) 測繪遙感信息工程國家重點實驗室，湖北 武漢 430079；2.中國地質(zhì)大學(xué)(武漢) 地理與信息工程學(xué)院，湖北 武漢 430074)

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

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

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

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

1 高光譜影像解混

1.1 線性光譜混合模型

線性光譜混合模型假設(shè)觀測到的像元光譜是各端元光譜按一定比例線性混合而成的，如圖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ù)雜，到達高光譜傳感器的信號通常是電磁波多次反射的結(jié)果，故高光譜傳感器接收到的信號通常是鄰近地物反射信號的深度融合。以電磁波信號在到達傳感器前進行L次線性混合為例，深度線性光譜混合模型可用圖2簡要表示。

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

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

X=A1S1S2…SL+N

(2)

1.2 非負矩陣分解模型

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

(3)

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

1.3 深度非負矩陣分解模型

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

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

2 面向端元矩陣的深度非負矩陣分解

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

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

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

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

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

X=ALSL…S2S1+N

(4)

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

(5)

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

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

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

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

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

(6)

(7)

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

(8)

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

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

(9)

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

(10)

α‖Sl‖TV

(11)

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

3 實驗結(jié)果與分析

3.1 實驗設(shè)置

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

(12)

(13)

3.2 模擬實驗

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

圖5 模擬實驗端元光譜圖Figure 5. Spectral signatures of synthetic experiments

圖6 模擬實驗豐度圖Figure 6. Abundance maps of synthetic experiments

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

表1 不同算法模擬實驗定量評價結(jié)果

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

(a)

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

(a) (b) (c)

3.3 真實實驗

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

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

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

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

(a)

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

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

(a) (b) (c) (d) (e) (f) (g) 圖11 不同方法提取的Samson數(shù)據(jù)集端元“水”的豐度圖(a)地面真實值 (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ù)集端元“樹”的豐度圖(a)地面真實值 (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)地面真實值 (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é)束語

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

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