楊露露，李春芝，陳曉華，王麗

湖州師范學(xué)院 信息工程學(xué)院,湖州 313000

1 引言

高光譜遙感技術(shù)作為一種新型對地觀測技術(shù)，提供了豐富的地物光譜信息，能夠?qū)Φ匚锬繕?biāo)進(jìn)行識別（Zhang 等，2019）和精細(xì)分類（Ghamisi等，2017；Khan等，2018），逐漸成為當(dāng)前遙感技術(shù)研究前沿。由于受空間分辨率的限制，高光譜圖像中單個像元包含多種地物使得混合像元的現(xiàn)象普遍存在，成為像元級應(yīng)用精度難以提高的主要問題（藍(lán)金輝 等，2018）。因此，對混合像元進(jìn)行分解，通過提取構(gòu)成像元的各類基本物質(zhì)組成成分的光譜特征以及這些物質(zhì)在像元內(nèi)所占的比例可以提高高光譜遙感定量應(yīng)用的精度（童慶禧 等，2016；張兵，2016）。

現(xiàn)有的光譜解混方法通?；诰€性混合模型（LMM）或非線性混合模型（NLMM）（Heylen 等，2014）。其中，線性混合模型由于簡單、易求解得到廣泛應(yīng)用。高光譜盲解混能夠同時進(jìn)行端元提取以及豐度估計，本質(zhì)上屬于信號處理領(lǐng)域中的盲源分離問題（Yao 等，2019）。盲源分離的典型方法有獨(dú)立成分分析（ICA）（Wang和Chang，2006）和非負(fù)矩陣分解（NMF）（Lee 和Seung，1999）。由于ICA對源信號做出的獨(dú)立性假設(shè)與高光譜圖像中的豐度和為一約束矛盾，所以在光譜解混中的應(yīng)用受到了限制（Nascimento 和Dias，2005）。與ICA 相比，NMF 能夠更好地適應(yīng)光譜解混的需求（Tsinos等，2017），但僅有的非負(fù)約束易導(dǎo)致NMF陷入局部最優(yōu)解。因此，需要另外引入合適的約束以獲取更好的解混性能。

考慮到LMM凸幾何的特點(diǎn)，Miao和Qi（2007）引入端元單形體體積約束，提出最小體積約束NMF（MVCNMF）。受Miao 和Qi 的啟發(fā)，后繼幾種基于幾何約束的解混算法被提出（Yu 和Sun，2007；Wang 等，2013；Zhi 等，2017）。這些方法首先假設(shè)高光譜圖像位于一個單形體中，并且其頂點(diǎn)被認(rèn)為是所要提取的端元。Marrinan 和Gillis（2021）為解決端元的不平衡性問題，提出一種基于塊的minimaxNMF 模型來處理含有稀少端元的高光譜圖像。此外，由于自然界中地物一般呈現(xiàn)分段的局域同質(zhì)分布，一個像元不會包含所有端元物質(zhì)，因此，除了端元的物理約束外，豐度的稀疏性也可以被用于構(gòu)建NMF解混框架。Qian等（2011）、Yuan 等（2015）、Feng 等（2018）采用L1/2-NMF，與L1范數(shù)（Bruckstein 等，2008）相比，結(jié)果表明L1/2-NMF更能滿足和為一約束，并且在保證豐度稀疏的同時，降低了計算復(fù)雜度。進(jìn)一步地，為保持?jǐn)?shù)據(jù)的結(jié)構(gòu)信息，He 等（2016）施加L1，2范數(shù)進(jìn)行聯(lián)合稀疏約束，解決了L1/2范數(shù)易受噪聲干擾的問題。Li 等（2016）考慮到同一個端元僅分布在個別像元中，應(yīng)用L2，1范數(shù)對豐度的行和列同時進(jìn)行稀疏約束。Han 等（2020）融合低秩性約束，引入使豐度更稀疏的L2，p范數(shù)。上述方法雖然采取了相關(guān)的約束，但忽略了數(shù)據(jù)的空間信息。Wang等（2017）利用線性迭代聚類分割技術(shù)進(jìn)行空間群稀疏約束，結(jié)合豐度的稀疏性提出群稀疏正則NMF（SGSNMF）。Liu 等（2011）融合豐度的分離性和平滑性約束，在保證相同像元中不同端元相關(guān)性的同時考慮數(shù)據(jù)的空間關(guān)系。Qin 等（2021）引入圖全變差約束（gtv）提高豐度的空間平滑度，其采用交替方向乘子法（ADMM）優(yōu)化求解，并且融合Nystr?m 近似策略和Merriman-Bence-Osher（MBO）方法降低求解復(fù)雜度。Wang 等（2016）考慮圖像的空間和光譜信息，利用超圖建模光譜的空間結(jié)構(gòu)以保證相同超邊內(nèi)像元豐度的一致性。此外，由于高光譜數(shù)據(jù)易受復(fù)雜非線性效應(yīng)影響，F(xiàn)évotte和Dobigeon（2015）提出了魯棒NMF（rNMF），該算法引入群稀疏離群項(xiàng)克服非線性問題。為挖掘高光譜數(shù)據(jù)中的非線性結(jié)構(gòu)，Mei 等（2016）、甘玉泉等（2019）、Rathnayake等（2020）將圖模型理論嵌入到NMF模型中，以保持?jǐn)?shù)據(jù)的局部幾何結(jié)構(gòu)。Lu等（2013）在L1/2-NMF的基礎(chǔ)上，亦構(gòu)建了額外的圖正則項(xiàng)，保證各端元豐度之間的緊密聯(lián)系。

現(xiàn)存NMF 解混算法針對端元和豐度的物理含義，分別施加相應(yīng)的約束以獲得更精確的解混結(jié)果，但大多基于歐氏距離構(gòu)造目標(biāo)函數(shù)，而高光譜圖像內(nèi)部結(jié)構(gòu)的復(fù)雜性導(dǎo)致簡單的線性表示不能很好地度量重構(gòu)誤差。Li 和Chen（2020）提出相關(guān)熵誘導(dǎo)度量代替歐氏距離構(gòu)造的損失函數(shù)，克服了基于歐氏距離的建模難以降低重構(gòu)誤差且極易受異常值影響的弊端。受最優(yōu)傳輸思想啟發(fā)，Rubner等（2000）提出地球移動距離（Earth Mover’s Distance，EMD）計算兩個分布之間的距離以反映最小需求量。鑒于該最優(yōu)傳輸問題，Sandler 和Lindenbaum（2011）提出基于EMD 的非負(fù)矩陣分解算法（EMDNMF），并證實(shí)EMD對不同維度噪聲不敏感，因此，基于EMD 的NMF 算法可以獲得比傳統(tǒng)基于歐氏距離的NMF算法更魯棒的性能。

由于高光譜圖像內(nèi)部包含大量不同維度的噪聲（Li等，2021），而EMD 具有對不同維度之間噪聲不敏感的優(yōu)勢，因此，經(jīng)過初步理論分析，EMD 比歐氏距離更適用于處理含有復(fù)雜噪聲的高光譜數(shù)據(jù)。由Sandler 和Lindenbaum（2011）提出的EMDNMF 在計算機(jī)視覺中的應(yīng)用可知，基于EMD 的NMF 在高光譜數(shù)據(jù)分析中亦可取得不錯的效果。但實(shí)際應(yīng)用中，由于高光譜圖像屬于復(fù)雜大數(shù)據(jù)，EMDNMF 高額的計算代價導(dǎo)致其在大規(guī)模高光譜數(shù)據(jù)分析中存在一定的局限性。Qian 等（2016）提出了一種基于Sinkhorn 距離的非負(fù)矩陣分解方法，可以更簡單快速地求解該距離下的分解模型。Zhang 等（2020）進(jìn)一步引入對偶圖正則項(xiàng)，同時考慮數(shù)據(jù)和特征的流形結(jié)構(gòu)。上述基于Sinkhorn 距離的非負(fù)矩陣分解算法被廣泛應(yīng)用于圖像聚類并取得了不錯的效果。

本文首次將最優(yōu)傳輸技術(shù)和圖模型理論綜合應(yīng)用于高光譜解混中，提出一種基于Sinkhorn距離和圖正則約束的高效解混算法（SDGNMF）。該算法在充分挖掘EMD 優(yōu)點(diǎn)的基礎(chǔ)上，對EMD 施加熵正則約束，將EMD 改進(jìn)為Sinkhorn 距離，并將其作為度量誤差的標(biāo)準(zhǔn)，有效降低了計算復(fù)雜度。此外，被施加了熵正則約束的EMD，即Sinkhorn距離對模型的表示可以更好地建模不同維度特征之間的關(guān)系，特征的相關(guān)性得以充分利用。本文在Sinkhorn距離的基礎(chǔ)上特別引入圖正則約束以進(jìn)一步刻畫數(shù)據(jù)的流形結(jié)構(gòu)。通過在模擬數(shù)據(jù)集以及真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)驗(yàn)證，以期有效改善基于NMF的高光譜解混算法性能。

2 理論介紹

2.1 線性光譜混合模型

線性混合模型假設(shè)入射光在到達(dá)傳感器前只與一種地物發(fā)生作用，不考慮物質(zhì)內(nèi)部反射、散射等作用，混合像元光譜表示為端元光譜與其相應(yīng)豐度的線性組合，數(shù)學(xué)表達(dá)式定義如下：

式中，X∈RL×N表示原始圖像，U∈RL×P表示端元矩陣，V∈RP×N表示豐度矩陣，E∈RL×N表示重構(gòu)誤差矩陣，通常被認(rèn)為是加性噪聲（Akhtar和Mian，2017；Zhang 等，2018）。L為波段數(shù)，P為端元數(shù)量，N為圖像中所含像元總數(shù)。矩陣U，V，E均滿足非負(fù)性約束（ANC），并且豐度V還需滿足和為一約束（ASC）。

2.2 非負(fù)矩陣分解

NMF 是由Lee 和Seung（1999）提出的一種無監(jiān)督特征提取方法。對于給定的非負(fù)矩陣X，NMF將其分解成兩個非負(fù)矩陣U和V，使得分解后兩個矩陣的乘積近似等于原始矩陣，即X≈UV。

為了衡量原始矩陣對分解因子間的重構(gòu)效果，需要構(gòu)造目標(biāo)函數(shù)，定義如下：

式中，Dist?表示度量距離，通常定義為歐氏距離，基于該距離的目標(biāo)函數(shù)表示為

式中，J（·）表示度量函數(shù)，‖ · ‖F(xiàn)表示Frobenius范數(shù)。

對于式（3），一般采用乘性迭代算法進(jìn)行更新，得到U和V如下

式中，（·）T表示矩陣轉(zhuǎn)置。

3 基于Sinkhorn距離和圖正則約束的非負(fù)矩陣解混模型

高光譜圖像內(nèi)部的復(fù)雜性導(dǎo)致數(shù)據(jù)中相鄰元素可能由相似的端元構(gòu)成，并且對應(yīng)的豐度也非常接近，可見，高光譜數(shù)據(jù)中存在著豐富的相關(guān)特征。傳統(tǒng)解混算法中，普遍采用歐氏距離構(gòu)造目標(biāo)函數(shù)，無法準(zhǔn)確描述高光譜圖像中各維度特征間關(guān)系。為充分利用圖像空間中特征的相關(guān)性，本文采用基于Sinkhorn距離測量誤差的方法。

3.1 Sinkhorn距離

在本文所提出的SDGNMF 方法中，采用改進(jìn)的EMD 代替?zhèn)鹘y(tǒng)歐氏距離構(gòu)造目標(biāo)函數(shù)。其中，EMD 定義為兩個分布直方圖之間傳輸元素的最小代價。高光譜的EMD 原理詳見圖1，圖1中，紅框表示選取圖像中某一區(qū)域。該區(qū)域的直方圖表示為x；y為經(jīng)EMD 處理后同一區(qū)域的直方圖表示；q、p分別表示直方圖x、y中的區(qū)段號。選取Urban數(shù)據(jù)集中某一波段的灰度圖像，紅框框選其某一區(qū)域的直方圖表示為x；假設(shè)x分布直方圖和y分布直方圖均含有n個區(qū)段數(shù)，從x分布直方圖中選取數(shù)據(jù)組合為r1的元素，共計個元素，產(chǎn)生轉(zhuǎn)移量可表示為將轉(zhuǎn)移量運(yùn)送至y分布直方圖的各區(qū)段中，生成數(shù)據(jù)組合為r2的元素，即式中，q、p分別表示x、y分布圖中的區(qū)段號；理想情況下根據(jù)上述傳輸方案依此反復(fù)，直至x中所有元素被全部分配至y分布直方圖中。EMDdM(x，y)計算公式如下（Rubner等，2000）：

圖1 高光譜圖像的EMD原理Fig.1 EMD principle diagram of hyperspectral image

式中，Φp表示p區(qū)段內(nèi)元素數(shù)，Φq表示q區(qū)段內(nèi)元素數(shù)。

由于EMD 計算過程中涉及到線性規(guī)劃的復(fù)雜求解問題，其計算復(fù)雜度高達(dá)O(n3logn)，限制了其在大規(guī)模數(shù)據(jù)分析中的適用性（Sandler 和Lindenbaum，2011）。本文在EMD 的基礎(chǔ)上引入熵正則項(xiàng)平滑參數(shù)化模型。將改進(jìn)后的EMD 距離稱為Sinkhorn距離，定義為

式中，1 表示全一向量，⊙表示逐元素乘法；ST表示矩陣S的轉(zhuǎn)置，H(S)=-∑p，qSlogS表示矩陣S的熵，λ表示正則化參數(shù)，用于控制熵正則化作用的強(qiáng)度，當(dāng)λ足夠大時，式（7）相當(dāng)于EMD。

對于上述模型的求解，Cuturi（2013）提出了一種Sinkhorn-Knopps 矩陣縮放算法，即通過生成一系列行和列被交替歸一化的矩陣來進(jìn)行求解。Genevay 等（2019）證明了基于該縮放算法的模型求解復(fù)雜度下降為O(n2)，因此，相較于EMD，Sinkhorn 距離可以更快速有效地進(jìn)行大規(guī)模數(shù)據(jù)分析。Sinkhorn-Knopps矩陣K定義為

然而，模型（7）的優(yōu)化需嚴(yán)格滿足歸一化條件，事實(shí)上，在高光譜圖像中，考慮到各像元的不規(guī)則性且相鄰像元間具有相似性，輸出圖像無法被準(zhǔn)確表示為歸一化的直方圖形式，不能嚴(yán)格滿足式（7）中的約束條件。因此，本文對式（7）進(jìn)行松弛，利用KL散度（Kullback-Leibler divergence）對原等式約束進(jìn)行軟懲罰，將模型簡化為一個無約束的問題。上述方法的有效性已在Frogner 等（2015）的研究中得到證明，本文將其推廣應(yīng)用到高光譜解混模型中。Sinkhorn距離最終優(yōu)化為

本文以兩個直方圖之間傳輸元素的最小成本為原則構(gòu)造目標(biāo)函數(shù)，進(jìn)而采用Sinkhorn距離作為誤差的度量標(biāo)準(zhǔn)，其優(yōu)勢具有如下兩方面：一方面，基于Sinkhorn距離的模型采用Sinkhorn-Knopps矩陣縮放算法進(jìn)行求解，相比于EMD，可以更快速有效地分析高光譜大規(guī)模數(shù)據(jù)；另一方面，考慮到高光譜數(shù)據(jù)中存在著豐富的相關(guān)特征，Sinkhorn距離可以克服噪聲影響，建模不同維度特征之間的關(guān)系，能夠充分利用特征間的相關(guān)性。因此，相較于傳統(tǒng)基于歐氏距離的NMF，基于Sinkhorn距離的NMF 模型，能夠更好地提取數(shù)據(jù)內(nèi)部結(jié)構(gòu)信息，提高解混精度。

3.2 光譜特征圖正則化

高光譜圖像中數(shù)據(jù)通常位于一個低維子流形上。為考慮空間信息并有效挖掘數(shù)據(jù)的流形結(jié)構(gòu)，本文在利用Sinkhorn距離構(gòu)造目標(biāo)函數(shù)的基礎(chǔ)上對豐度矩陣V施加圖正則約束，進(jìn)一步捕捉數(shù)據(jù)間的有效特征，從而增強(qiáng)圖正則約束描述幾何流形結(jié)構(gòu)的有效性。一般地，圖正則項(xiàng)表述如下：

式中，Wij表示圖正則項(xiàng)權(quán)重矩陣，可以通過熱核函數(shù)計算得到（Lu等，2013），具體計算公式如下：

式中，xi和xj表示數(shù)據(jù)空間中兩個鄰近點(diǎn)，τ表示核參數(shù)，通常取值為［0，1］。

3.3 SDGNMF模型

綜合式（9）、（10）、（11），由Sinkhorn距離構(gòu)造的解混目標(biāo)函數(shù)定義如下：

式中，Xj表示矩陣X的第j列向量；ξ表示維持特征相關(guān)性和數(shù)據(jù)流形之間平衡的正則化參數(shù)。

由式（12）可知，所提出的模型是非凸問題，同時優(yōu)化U和V是一個NP 難問題。參考Lee 和Seung（2000），Qian 等（2016）將乘性迭代算法應(yīng)用求解Sinkhorn 距離并證明了該算法在Sinkhorn距離上的收斂性。為將變量交替最小化，本文同樣采用乘性迭代算法對式（12）進(jìn)行求解，當(dāng)一個變量迭代時，將其他變量固定為常數(shù)，每個變量在當(dāng)前固定的變量上有條件地更新，從而將所提出的模型計算問題分解成兩個凸問題。

給定非負(fù)矩陣Θ∈RL×P和Ψ∈RP×N作為U和V的拉格朗日乘子，定義拉格朗日函數(shù)Γ為

根 據(jù) Karush-Kuhn-Tucker 條 件，即Θ.×U=0，Ψ.×V=0，可以得出U和V的更新規(guī)則：

式中，S*表示最優(yōu)傳輸矩陣，定義為Sinkhorn-Knopps矩陣K=exp(-λM-1)的對角縮放：

式中，A和B表示兩個縮放矩陣，diag(·)表示對角矩陣。Frogner 等（2015）定義了上述兩個縮放矩陣的更新規(guī)則，具體表達(dá)式描述如下：

3.4 算法流程

本文參照虛擬維度算法（Bajorski，2011）計算端元數(shù)量，同時采用頂點(diǎn)成分分析（VCA）初始化端元矩陣U，完全約束最小二乘法（FCLS）初始化豐度矩陣V，為了計算圖正則項(xiàng)權(quán)重矩陣，將核參數(shù)τ設(shè)置為0.2，并將K-近鄰（KNN）的數(shù)目設(shè)置為5。本文設(shè)置最大迭代次數(shù)Imax，如果迭代次數(shù)大于定義的Imax，則迭代結(jié)束。具體算法流程如下：

輸入：原始數(shù)據(jù)矩陣X，參數(shù)λ，γ，ξ，τ，KNN數(shù)目，最大迭代次數(shù)Imax

步驟1：通過VCA-FCLS算法初始化U和V

步驟2：根據(jù)式（11）計算權(quán)重矩陣W

步驟3：根據(jù)式（6）計算距離矩陣M

步驟4：根據(jù)式（8）計算矩陣K

步驟5：循環(huán)

根據(jù)式（17）更新得到A

根據(jù)式（18）更新得到B

根據(jù)式（16）更新得到S*

根據(jù)式（14）更新得到U

根據(jù)式（15）更新得到V

直到滿足停止條件

輸出：U和V

4 結(jié)果與分析

為了證明提出算法的有效性，本文分別在模擬高光譜數(shù)據(jù)集和真實(shí)高光譜數(shù)據(jù)集上進(jìn)行性能驗(yàn)證，并與GNMF（Lu等，2013），SGSNMF（Wang等，2017），rNMF（Févotte 和Dobigeon，2015），minimaxNMF（Marrinan 和Gillis，2021），gtvMBO（Qin 等，2021），VCA-FCLS 等6 種算法進(jìn)行對比分析。所有對比算法均設(shè)置為原文獻(xiàn)中參數(shù)。實(shí)驗(yàn)平臺為CPU 2.9 GHz，內(nèi)存8 G 的Windows10電腦上進(jìn)行，使用Matlab2019a進(jìn)行代碼實(shí)現(xiàn)。

4.1 選定評價指標(biāo)

為定量評估所提出SDGNMF 解混算法的性能，本文采用高光譜解混評估中常用的兩種評價指標(biāo)，即光譜角距離（SAD）和均方根誤差（RMSE）。SAD 用來評估端元估計值的準(zhǔn)確性，RMSE 用來評估豐度估計值的準(zhǔn)確性。SAD和RMSE的表達(dá)式定義如下：

式中，ui和分別表示第i個端元真實(shí)值和估計值；vi和分別表示第i個端元對應(yīng)豐度真實(shí)值和估計值。

4.2 算法收斂性

對SDGNMF 算法的收斂性能的驗(yàn)證結(jié)果如圖2所示。可見在所有的測試數(shù)據(jù)集上，算法的重構(gòu)誤差隨著迭代次數(shù)的增加而減小，當(dāng)達(dá)到一定的迭代次數(shù)時，誤差達(dá)到一個穩(wěn)定的值。為了保證一致性，在所有數(shù)據(jù)集上迭代次數(shù)均設(shè)置為1000次。

圖2 誤差隨迭代次數(shù)變化Fig.2 Figure of error changing with iteration number

4.3 模擬數(shù)據(jù)集實(shí)驗(yàn)

該模擬數(shù)據(jù)集是根據(jù)Hendrix 等（2012）提出的模擬程序生成的，即從美國地質(zhì)調(diào)查局（USGS）光譜庫中選擇9 種物質(zhì)221 個波段的譜線作為端元，生成大小為100×100 的圖像，圖3 展示了9 種端元光譜特征曲線。為滿足和為一約束以及非負(fù)性約束，通過K均值聚類和高斯濾波器生成相應(yīng)的豐度圖，每個端元對應(yīng)的豐度值隨機(jī)分布范圍為0.004—0.900，最大豐度值為0.900。

圖3 模擬數(shù)據(jù)端元光譜曲線Fig.3 Endmember spectral curve of simulated data

4.3.1 參數(shù)選擇

參數(shù)λ、γ和ξ分別為正則化參數(shù)、松弛參數(shù)和平衡參數(shù)，分別用于調(diào)節(jié)熵正則作用的強(qiáng)度、調(diào)節(jié)模型的松弛程度以及維持特征相關(guān)性和數(shù)據(jù)流形之間的平衡。為研究參數(shù)λ、γ和ξ對算法性能的影響，在模擬數(shù)據(jù)集上，采用固定初始化的方式對3 個參數(shù)進(jìn)行測試。首先，固定γ=ξ=10以測試參數(shù)λ對算法性能的影響，參考Qian 等（2016）的研究，參數(shù)λ取值從［10，50，100，200，300，400，500］中選取，結(jié)果見圖4（a）?？梢?，當(dāng)λ從10 增加到200 時，SAD 和RMSE 值不斷減小并達(dá)到最小值，而隨著λ繼續(xù)增加，SAD和RMSE 值呈上升趨勢，因此，設(shè)置λ=200。其次，固定λ=200，ξ仍設(shè)置為10，設(shè)置γ取值范圍為［10，80］，步長為10 遞增，結(jié)果見圖4（b）。可見，RMSE 值在這區(qū)間保持穩(wěn)定，而當(dāng)γ=50 時，SAD 達(dá)最小值，因此，設(shè)置γ=50。最后，固定λ=200，γ=50，設(shè)置ξ取值為［1E-6，10］，以10 倍步長增加，結(jié)果見圖4（c）?？梢?，參數(shù)ξ在1E-6—1E-2變化時，SAD 與RMSE 值均保持穩(wěn)定。由于參數(shù)ξ用于維持特征相關(guān)性和數(shù)據(jù)流形之間的平衡，當(dāng)ξ＞0.1 時，重構(gòu)項(xiàng)過于松弛，SAD 和RMSE 估計精度呈明顯下降趨勢。因此，為了獲得更精確的解混結(jié)果，將ξ的最優(yōu)值設(shè)置在［1E-6，1E-2］。在以下模擬數(shù)據(jù)集實(shí)驗(yàn)中，參數(shù)ξ均設(shè)置為1E-4，同時分別設(shè)置λ=200和γ=50。

圖4 SDGNMF在不同參數(shù)下的SAD和RMSE值Fig.4 SAD and RMSE values of SDGNMF with different parameters

4.3.2 魯棒性分析

在本實(shí)驗(yàn)中，旨在對比和分析噪聲對算法的影響。在原始無噪聲模擬數(shù)據(jù)集中加入零均值高斯白噪聲，信噪比（SNR）定義為式中，X和e分別表示原始信號和其相應(yīng)的噪聲，E[·]表示期望運(yùn)算符。

為了測試七種算法的魯棒性，在模擬數(shù)據(jù)集上添加加性噪聲，信噪比分別設(shè)置為：10、20、30、40、50、60 dB。每個信噪比下的柱狀圖從左至右依次為算法SDGNMF、GNMF、SGSNMF、rNMF、minimaxNMF、gtvMBO 和VCA，所有方法的SAD值和RMSE 值均取10 次實(shí)驗(yàn)結(jié)果的平均值，結(jié)果見圖5?？梢姡弘S著信噪比的降低，所有算法的解混精度都隨之降低；在不同的信噪比情況下，所提出的SDGNMF 算法獲得的SAD 值和RMSE 值比其他算法更小。因此，實(shí)驗(yàn)結(jié)果表明，噪聲會導(dǎo)致解混結(jié)果變差，但所提出的算法相比于其他算法更具魯棒性。

圖5 不同SNR下幾種算法的SAD和RMSE對比Fig.5 Several algorithms under different SNR with the comparison of SAD and RMSE

4.3.3 不同算法豐度圖對比分析

在本實(shí)驗(yàn)中，設(shè)置參數(shù)λ=200、γ=50、ξ=1E-4，對比分析了模擬數(shù)據(jù)集上其他幾種算法提取的豐度效果圖，結(jié)果如圖6所示。

圖6 各算法在模擬數(shù)據(jù)集上的豐度對比Fig.6 Comparison of the abundance maps of each algorithm on simulated dataset

可見，顯示SDGNMF、SGSNMF、minimaxNMF以及VCA-FCLS 可以較好地恢復(fù)豐度。SGSNMF采用超像素分割法實(shí)現(xiàn)空間群稀疏約束，可以保留空間局部信息，但由于沒有準(zhǔn)確考慮相似像元區(qū)域內(nèi)的不同端元結(jié)構(gòu)，SGSNMF 在得到Kaolinite Kgal端元的豐度圖時，誤提取了Dumortierite 端元。VCA-FCLS 算法也同樣存在上述Kaolinite Kgal端元和Dumortierite端元提取不準(zhǔn)確的問題。minimaxNMF在Kaolinite CM9 豐度圖中出現(xiàn)了Dumortierite 端元輪廓，且Kaolinite Kgal端元提取精度較低。gtvMBO無法準(zhǔn)確區(qū)分具有相似結(jié)構(gòu)的端元，導(dǎo)致Kaolinite Kgal 和Kaolinite CM9 端元提取出現(xiàn)混疊。而本文提出的SDGNMF 算法優(yōu)勢在于能夠保留迭代后的相似結(jié)構(gòu)，各端元特征之間的相關(guān)性得到了充分考慮，因而相鄰區(qū)域分布的相似物質(zhì)能夠分離出來。相比于其他幾種算法，SDGNMF 可以更好地展示局部豐度的細(xì)節(jié)，獲得更真實(shí)、完善的豐度圖。GNMF 和rNMF 表現(xiàn)欠佳。GNMF 由于無法給出正確的稀疏表示，幾乎提取不到有價值的豐度信息，rNMF 獲得的豐度圖像素值為0.3—0.5，端元提取不純。綜上，解混結(jié)果表明所提出的SDGNMF 算法能夠更好地恢復(fù)模擬數(shù)據(jù)的豐度圖，相比于其他幾種對比算法，可以獲得更加精確的解混結(jié)果。

4.4 真實(shí)數(shù)據(jù)集實(shí)驗(yàn)

4.4.1 Urban數(shù)據(jù)集

Urban 數(shù)據(jù)集是高光譜解混研究中使用最廣泛的高光譜數(shù)據(jù)之一，可以在網(wǎng)站（http：//lesun.weebly.com/hyperspectral-data-set.html［2021-03-11］）上下載。其子圖像大小為307×307，共包含210個波段，波段范圍為400—2500 nm，空間分辨率和光譜分辨率分別為2 m 和10 nm?？紤]到密集水蒸氣和大氣效應(yīng)的影響，去除波段1—4、76、87、101—111、136—153 和198—210 后，保留了162 個波段進(jìn)行解混。根據(jù)實(shí)際地面參考圖，選取數(shù)據(jù)集中的6 種端元進(jìn)行實(shí)驗(yàn)，分別是瀝青、草地、樹木、屋頂、金屬和土壤。Urban 數(shù)據(jù)立方體和6種端元光譜曲線如圖7所示。

在本實(shí)驗(yàn)中，設(shè)置參數(shù)λ=300、γ=40、ξ=1E-4，得到的豐度圖和地面參考豐度圖如圖8 所示，表1列出了所有方法的SAD 值，其中粗體表示最優(yōu)值，下劃線表示次優(yōu)值。從圖中可以看出，對比參考豐度圖，本文提出的算法恢復(fù)了幾乎所有的端元成分。VCA-FCLS、gtvMBO 均在瀝青豐度圖中出現(xiàn)了小范圍的草地輪廓，并且VCA-FCLS 未提取到較純的草地豐度圖。GNMF僅提取到了泥土，其他5種端元顯示的豐度圖顏色均為綠色，即提取純度較低。對比GNMF，rNMF 豐度顏色值更小，位于0.2—0.3，可見，端元提取很不準(zhǔn)確。相較于上述幾種對比算法，SGSNMF 和minimaxNMF 得到了相對較好的豐度圖，但均由于提取金屬時像素值太小，導(dǎo)致準(zhǔn)確識別該端元無法實(shí)現(xiàn)。此外，表1還表明，SDGNMF 估計出的大部分端元的SAD 值最優(yōu)，雖然個別提取端元的SAD 值次優(yōu)，但平均SAD 比其他算法都低。因此，從純度和準(zhǔn)確性的角度來看，本文提出的算法較其余算法能夠取得更精確的解混效果。

表1 Urban數(shù)據(jù)集各算法的SAD對比值Table 1 SAD ratio of algorithms in Urban dataset

圖8 各算法在Urban數(shù)據(jù)集上的豐度對比Fig.8 Comparison of abundance maps of each algorithm on Urban dataset

4.4.2 Jasper數(shù)據(jù)集

Jasper 數(shù)據(jù)集可以在網(wǎng)站（http：//lesun.weebly.com/hyperspectral-data-set.html［2021-03-11］）上下載。其圖像大小為512×614，共包含224個波段，范圍為380—2500 nm。數(shù)據(jù)集中含4 種端元，即樹木，水，土壤和道路?？紤]到水汽波段的影響，去除波段1—3、108—112、154—166 和220—224 之后，保留了198個通道來進(jìn)行光譜解混。由于此高光譜圖像的復(fù)雜性，無法精確獲取地面實(shí)際情況，因此本文采用100×100 大小的子圖像進(jìn)行實(shí)驗(yàn)。Jasper數(shù)據(jù)立方體和4種端元光譜曲線如圖9所示。

圖9 Jasper高光譜數(shù)據(jù)Fig.9 Jasper hyperspectral data

設(shè)置參數(shù)λ=300、γ=10、ξ=1E -4，實(shí)驗(yàn)得到的豐度圖和地面參考豐度圖如圖10 所示，表2列出了所有方法的SAD值，其中粗體表示最優(yōu)值，下劃線表示次優(yōu)值。從圖中可以看出，本文提出的算法比其他方法更好地恢復(fù)了豐度圖，與參考地物分布十分吻合。VCA-FCLS 和SGSNMF算法均在泥土豐度圖中錯誤地提取了道路，導(dǎo)致道路無法準(zhǔn)確被識別。minimaxNMF 在提取水流端元時不準(zhǔn)確，豐度圖中出現(xiàn)了道路的輪廓。而其他幾種算法均存在無法識別端元，提取精度不純的問題。此外，如表2 所示，SDGNMF 大多可以獲得比其他方法更小的SAD 值，且平均SAD 值最小。上述結(jié)果表明，本文提出的算法可以更準(zhǔn)確地提取端元，并且獲得的豐度圖精度高于其他方法。

表2 Jasper數(shù)據(jù)集各算法的SAD對比值Table 2 SAD ratio of algorithms in Jasper dataset

圖10 各算法在Jasper數(shù)據(jù)集上的豐度對比Fig.10 Comparison of abundance maps of each algorithm on Jasper dataset

最后，表3給出了本文算法以及其他幾種對比算法在兩種真實(shí)數(shù)據(jù)集上的運(yùn)行時間。盡管根據(jù)3.1 節(jié)中理論分析可知，本文提出的SDGNMF 相比于EMD，一定程度上減少了時間復(fù)雜度。但從表3中可以看出SDGNMF 相對于其他解混算法并不是最節(jié)省時間的算法。這是因?yàn)镾DGNMF 相對于SGSNMF、rNMF、minimaxNMF 這3 種算法，額外引入的圖正則項(xiàng)對圖的構(gòu)造付出了很高的計算成本，而對于同樣需要構(gòu)造圖的GNMF 和gtvMBO 算法，額外的縮放矩陣迭代求解使得算法迭代次數(shù)增大，也在一定程度上降低了計算效率。

表3 真實(shí)數(shù)據(jù)集上各算法運(yùn)行時間比較Table 3 Runtime of algorithms on real datasets

5 結(jié)論

本文提出一種新的光譜解混模型，即SDGNMF。在模擬數(shù)據(jù)集、Urban 數(shù)據(jù)集以及Jasper 數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果均表明所提出方法的有效性和優(yōu)越性。不同于傳統(tǒng)解混模型，SDGNMF 模型能夠有效克服噪聲，并同時考慮特征間的相關(guān)性以及數(shù)據(jù)間的流形結(jié)構(gòu)。與傳統(tǒng)的歐氏距離、EMD 相比較，SDGNMF模型采用Sinkhorn距離來度量重構(gòu)誤差具有如下優(yōu)勢：（1）相比于EMD，Sinkhorn距離可以更快速有效地分析高光譜圖像中的大規(guī)模數(shù)據(jù)；（2）相比于歐氏距離，Sinkhorn 距離對高光譜數(shù)據(jù)中的噪聲相對不敏感，同時對不同維度上的特征能夠分別進(jìn)行有效建模，特征間的相關(guān)性被充分挖掘利用。此外，鑒于SDGNMF 模型的低維表示空間，圖正則約束進(jìn)一步增強(qiáng)了對數(shù)據(jù)幾何流形結(jié)構(gòu)的有效描述。

但是，本文所提出的算法仍然存在計算復(fù)雜度較高的問題，因此，后續(xù)工作將嘗試融合豐度以及端元相關(guān)的先驗(yàn)信息于同一框架下，開發(fā)更高效的光譜解混算法。