陳善學(xué)，董桓宇，陳雯雯

(1.重慶郵電大學(xué) 通信與信息工程學(xué)院,重慶 400065；2.重慶郵電大學(xué) 移動通信技術(shù)重慶市重點實驗室,重慶 400065)

0 引 言

高光譜圖像具有較高的空間分辨率與豐富的光譜信息的特點使其在多個領(lǐng)域有著卓越的表現(xiàn)，例如在軍事、環(huán)境、海洋、地理、農(nóng)業(yè)等領(lǐng)域，高光譜遙感技術(shù)都發(fā)揮著重要的作用，因此，高光譜技術(shù)的研究在國防、生產(chǎn)、生活中顯得極為重要。尤其高光譜圖像解混技術(shù)，其對地物成分和光譜特性的體現(xiàn)，使其成為高光譜技術(shù)研究的新秀和寵兒，成為高光譜遙感技術(shù)研究最為重要的內(nèi)容之一。

由于高光譜圖像具有豐富的光譜信息和空間特征，使得解混方法蓬勃發(fā)展。端元提取算法主要有純像元指數(shù)法(pixel purity index, PPI)[1]，最大體積法N-Finder[2]，頂點成分分析(vertex component analysis, VCA)[3]，正交基算法(orthogonal bases algorithm, OBA)[4]，曲線成分分析(curvilinear component analysis, CCA)[5]等，豐度估計較常用的算法是全約束最小二乘法(fully constrained least squares, FCLS)[6-7]。

文獻(xiàn)[8]和文獻(xiàn)[9]分別詳細(xì)地總結(jié)了國內(nèi)外基于線性模型和非線性模型的高光譜圖像解混的研究情況，并針對近幾年的研究成果提出了問題與展望。約束非負(fù)矩陣分解是高光譜解混算法中的一個熱點：文獻(xiàn)[10]提出最小體積約束非負(fù)矩陣分解(minimum volume constrained nonnegative matrix factorization,MVC-NMF),用最小體積項(minimum volume term)作為目標(biāo)函數(shù)的約束條件，構(gòu)成新的目標(biāo)函數(shù)，這一算法被后來的一些研究人員進(jìn)一步研究改進(jìn)[11-12]。文獻(xiàn)[13-14]進(jìn)行了關(guān)于L1/2稀疏性約束非負(fù)矩陣分解(L1/2-NMF)和平滑性約束非負(fù)矩陣分解(smoothness constraint nonnegative matrix factorization, CNMF)的研究。近年來國內(nèi)外越來越多的研究人員關(guān)注于約束非負(fù)矩陣分解的算法研究。文獻(xiàn)[15]提出一種了基于神經(jīng)網(wǎng)絡(luò)的稀疏性約束非負(fù)矩陣分解方法：L1-DNMF算法。文獻(xiàn)[16]也研究了結(jié)合深度學(xué)習(xí)的總體差值稀疏性約束非負(fù)矩陣分解(sparsity-constrained deep nonnegative matrix factorization with total variation, SDNMF-TV)的方法。文獻(xiàn)[17]提出基于非線性模型，對L1/2稀疏性約束和平滑性約束的結(jié)合方法，用核方法進(jìn)行優(yōu)化的算法，即基于豐度約束核非負(fù)矩陣分解的非線性光譜解混算法(abundance sparseness and smoothness constrained KNMF, ASSKNMF)。文獻(xiàn)[18]通過結(jié)合Euclidean范數(shù)稀疏約束和Frobenius范數(shù)稀疏約束的方法，提出了一種強(qiáng)魯棒性約束非負(fù)矩陣分解的解混方法。還有其他許多關(guān)于約束非負(fù)矩陣分解的解混方法，這里就不做介紹了。

以上介紹的關(guān)于約束非負(fù)矩陣分解的解混算法，都有著明確的數(shù)學(xué)意義和一定程度上的物理意義，并在真實數(shù)據(jù)的實驗中，有著優(yōu)秀性能。但是這些方法或基于圖像的馬卡洛夫模型的隨機(jī)特性，或基于數(shù)據(jù)的稀疏特性，亦或是基于數(shù)據(jù)的空間連續(xù)性、平滑性，都是考慮到數(shù)據(jù)的數(shù)學(xué)特性，而在真實場景的物理背景上考慮的程度較為有限。STC-NMF算法在前人的線性模型上進(jìn)行約束非負(fù)矩陣分解算法研究的基礎(chǔ)上，受Hapke輻射傳輸理論[19]的啟發(fā)。根據(jù)輻射傳輸理論，在大氣的輻射傳輸中，大氣中存在的膠著物和懸浮物對光線的輻射傳輸有著不可忽視的作用。由大氣中的微粒，如煙、塵埃、小水滴及氣溶膠等會引起形成近似夫瑯禾費衍射現(xiàn)象，其方向性比較明確。表現(xiàn)在傳感信號上，由接收儀器接收，就會形成在本像素點的射入射出信號，即被鄰域的散射信號干擾，同時又干擾鄰域信號，因此，STC-NMF算法對物理層面有較多的考慮，這是不同于以往算法的。這一算法在目標(biāo)函數(shù)中，將干擾看作噪聲和具有現(xiàn)實物理意義的散射，并以散射項作為約束項參與約束算法。也正是因為如此，STC-NMF算法雖然還存在著許多不足，但在實驗中仍舊具有非常優(yōu)秀的解混性能。此外，由于約束項在數(shù)學(xué)形式與數(shù)學(xué)特性上與一些較為知名且被較多引用的算法，如CNMF，L-SNMF，SDNMF-TV等都具有較高的相似性，所以，STC-NMF算法在實驗仿真中也有著不錯的表現(xiàn)。

1 相關(guān)工作

1.1 線性混合模型

線性混合模型(linear mixing model, LMM)是一種重要的混合模型，具有簡單明了的數(shù)學(xué)特性和一定的物理意義，在高光譜解混領(lǐng)域擁有舉足輕重的地位，很多研究都建立在線性混合模型之上。在一般情況下，像元的反射率可以近似地看作端元的線性混合

R=AS+n

(1)

(1)式中：R為像元反射率；A為端元反射率矩陣；S為豐度矩陣，地物豐度滿足非負(fù)性與全加性；n為噪聲干擾。

1.2 非負(fù)矩陣分解

在線性解混模型中，將像元的反射率看作端元豐度線性組合，若在理想條件下，符合公式(1)，各種各樣的影響因素都可看作噪聲，那么有

n=R-AS

(2)

若使得‖n‖最小，采用文獻(xiàn)[20]提出的非負(fù)矩陣分解。構(gòu)建目標(biāo)函數(shù)為

(3)

則求min‖n‖的問題將等同于求arg minf(A,S)問題，根據(jù)梯度下降算法，其乘法迭代規(guī)則為

A←A.*(RST)./(ASST)

(4)

S←S.*(ATR)./(ATAS)

(5)

每個符號代表的意義同公式(1)，.*和./分別代表矩陣的點乘和點除。

2 散射項約束非負(fù)矩陣分解

散射相函數(shù)即為介質(zhì)相位函數(shù)，在輻射傳輸理論中常用以描述米氏散射(Mie scattering)，是對于處于均勻介質(zhì)的各向同性的單個介質(zhì)球在單色光照射下，基于麥克斯韋方程邊界條件下的嚴(yán)格數(shù)學(xué)解。其中最簡單最常用的是Henyey-greenstein項函數(shù)

(6)

但由于Henyey-greenstein項函數(shù)的形式較為復(fù)雜，因此本文選用相對簡明的，在文獻(xiàn)[19]中提出的勒讓德多項式形式的散射相函數(shù)。在文獻(xiàn)[19]中，Hapke提出勒讓德多項式形式的散射相函數(shù)為

P(g)=1+bcosg+c(3cos2g-1)/2

(7)

本文采取將勒讓德多項式形式的散射相函數(shù)作為目標(biāo)函數(shù)中的約束條件對圖像進(jìn)行解混，目標(biāo)函數(shù)變化為

(8)

(8)式中：λ定義下降速度；g是非對稱因子；b與c是可調(diào)節(jié)參數(shù)。由于入射波的相位在粒子上是不均勻的，造成了各子波在空間和時間上的相位差，因此，本文用地物在像素點的豐度與其四鄰域相同地物的豐度的平均值求差，將差值的絕對值作為非對稱因子描述由于大氣散射造成不同像素空間中相同地物間的散射項貢獻(xiàn)，則定義非對稱因子g為

(9)

(9)式中：Spk為k地物在p像素點的豐度；SNpk為這一地物在該像素點四鄰域豐度均值；p是像素點的四鄰域?！础当硎揪仃囋氐暮蜑?/p>

(10)

將散射項約束作為約束條件，其物理意義為，將干擾因素看成2部分組成，分別是噪聲和鄰近區(qū)域米氏散射影響，并將米氏散射對鄰域貢獻(xiàn)看成是線性的。以豐度差值作為非對稱因子，在宏觀上以單個像元所處空間作為單位，以此計算像素總體相位差異，求相同地物在鄰域具有不同空間占有率情況下，因米氏散射對該像元區(qū)域造成的貢獻(xiàn)。在數(shù)學(xué)上，文獻(xiàn)[14]采用MRF模型，將平滑性約束作為約束項，以此達(dá)到對數(shù)據(jù)進(jìn)行分段平滑性約束的目的。平滑性約束為

g(x)=-e-x2/γ+1

(11)

將約束項泰勒展開為麥克勞林級數(shù)有

(12)

而散射項約束項展開為

(13)

(13)式中，o(x),o(g),o(2g)為二階余項。在形式上g(x)與p(g)二階泰勒展開式相似，如果繼續(xù)展開到高階依舊相似，且對豐度約束而言，x和g的定義域相同，說明兩者有著相近的數(shù)學(xué)意義，即勒讓德多項式形式的散射項函數(shù)作為約束項也具有一定程度的分段平滑性約束的能力。由于相近的數(shù)學(xué)意義，則其梯度下降算法的新乘法迭代規(guī)則也可借鑒文獻(xiàn)[14,20-21]確定，即

A←A.*(RST)./(ASST)

(14)

S←S.*(ATR)./(ATAS+λS.*

[-bsinΨ-3ccosΨsinΨ]

(15)

(15)式中，矩陣Ψ為像素點位置對應(yīng)的非對稱因子g所組成的矩陣。

證明借鑒文獻(xiàn)[14,20]，定義有輔助函數(shù)G(S,St)使得

St+1=arg minG(S,St)

(16)

(16)式中：上標(biāo)t表示當(dāng)前狀態(tài)；t+1表示下一狀態(tài)，滿足迭代要求的關(guān)系

F(St+1)≤G(St+1,St)≤G(St,St)=F(St)

(17)

不妨假設(shè)輔助函數(shù)為

G(S,St)=F(St)+(S-St)T(?F(St))+

(18)

(18)式中：

K(St)=diag((ATASt)./St)+λP′(Ψ)

(19)

那么對比目標(biāo)函數(shù)為

F(S)=F(St)+(S-St)T?F(St)+

(20)

根據(jù)(20)式，有G(S,St)≥F(h)，即為

(S-St)T(ATA-K(ht))(S-St)≥0

(21)

可以推出?SG(S,St)趨向于零

(22)

通過化簡與簡化，關(guān)于豐度矩陣的乘法迭代規(guī)則可以確定為

(23)

STC-NMF算法如下。

輸入：高光譜圖像R，像素點數(shù)P，I行J列，波段數(shù)L;

步驟1：使用N-Finder或Hysime初始化端元矩陣A，然后用Fcls初始化豐度矩陣S;

步驟2：使用公式(14)更新端元矩陣A;

步驟3：使用公式(12)更新豐度矩陣S;

步驟4：一直重復(fù)步驟2，步驟3更新，直到達(dá)到最大迭代次數(shù)，或目標(biāo)函數(shù)，即公式(8)滿足fnew(A,S)-fold(A,S)

輸出：端元矩陣A與豐度矩陣S.

3 實驗與討論

3.1 數(shù)據(jù)集

本文使用的數(shù)據(jù)集是用由朱飛云博士整理并在科研在線公開的數(shù)據(jù)，與文獻(xiàn)[14]相同，數(shù)據(jù)是美國地質(zhì)調(diào)查局的urban場景，307×307像素的數(shù)據(jù)區(qū)域，去除1-4, 76, 87, 101-111, 136-153和 198-210波段的162個波段。圖1為該場景第9個波段的成像圖。

圖1 Urban場景第9個波段的成像圖Fig.1 9th band image of urban

3.2 參數(shù)設(shè)置與算法對比

仿真實驗在主頻3.7 GHz，內(nèi)存8 GB的PC機(jī)上進(jìn)行，操作系統(tǒng)為Window(64位)，仿真平臺為Matlab2014a(64位)。散射項約束非負(fù)矩陣分解中有3個參數(shù)λ,b,c可以設(shè)定,其中,λ調(diào)節(jié)乘法迭代規(guī)則的梯度下降速度，b,c作為調(diào)節(jié)參數(shù)，可根據(jù)對數(shù)據(jù)的適應(yīng)性做調(diào)整，根據(jù)(12),(13)式，這2個參數(shù)滿足特定條件時，可以看做進(jìn)行平滑性約束的塑形參數(shù)。這里將這3個參數(shù)依次設(shè)置為0.01,0.2和0.1進(jìn)行實驗。實驗L1/2-NMF算法設(shè)置調(diào)節(jié)參數(shù)λ為0.01，CNM算法設(shè)置調(diào)節(jié)參數(shù)λ為0.01，γ為0.01。實驗結(jié)果，評價指標(biāo)最優(yōu)結(jié)果用下橫線標(biāo)出。關(guān)于實驗結(jié)果的評價標(biāo)準(zhǔn)有：光譜信息散度(special information divergence, SID)，均方根誤差(root mean square error, RMSE)。本次實驗與還原的L1/2-NMF算法、CNMF算法以及VCA+FCLS算法進(jìn)行對比，實驗結(jié)果如表1和表2。

表1 Urban場景的4端元光譜信息散度

表2 Urban場景的4端元均方根誤差

以上端元對比，通過對比可以得到，在上述實驗條件下，與L1/2-NMF和CNMF相比較，STC-NMF算法性能有一定的提升，同時這3種算法相比較VCA+FCLS都有較為明顯的提升，雖然VCA+FCLS在瀝青路和屋頂2個端元解混效果也不錯，但相比較其他算法并沒有明顯優(yōu)勢，所以暫時不做過多討論。在表1、表2的具體數(shù)值比較，STC-NMF算法解混的4種地物端元都擁有著優(yōu)秀的性能，但與L1/2-NMF和CNMF相比較并沒有特別明顯的優(yōu)勢。從端元的波形圖來看，這3種算法解混出來的波形與真實地物的波形都很接近，且大致輪廓基本一致，其中STC-NMF算法計算出來的端元與真實端元最為接近；但是這3種算法在關(guān)于屋頂?shù)亩嗽饣?，效果都不太好，有較為明顯的差異，但其中STC-NMF算法仍優(yōu)于另外2種算法解混出來的，在波形上與真實地物也更為貼合。

圖2 4種算法與真實端元的端元對比Fig.2 Endmembers comparison with real endmembers of the four algorithms

圖3 真實地物四端元豐度圖Fig.3 Abundance for 4 endmembers of ground truth

圖4 STC-NMF算法四端元豐度圖Fig.4 Abundance for 4 endmembers of STC-NMF

圖5 L1/2-NMF算法四端元豐度圖Fig.5 Abundance for 4 endmembers of L1/2-NMF

圖6 CNMF算法四端元豐度圖Fig.6 Abundance for 4 endmembers of CNMF

圖7 VCA算法四端元豐度圖Fig.7 Abundance for 4 endmembers of VCA

從豐度圖直觀的視覺感受來看，3種算法皆明顯優(yōu)于VCA，因此，在此處也暫時對VCA不做過多討論。STC-NMF對于瀝青路的解混效果略優(yōu)于CNMF，但略差于L1/2-NMF；對草地的解混效果3種算法無明顯的感覺差異；對于樹木的解混效果STC-NMF略優(yōu)于另2種算法；對于屋頂?shù)慕饣煨Ч麤]有造成感官上的明顯優(yōu)勢，3種算法都有將瀝青路出現(xiàn)在屋頂?shù)牡匚镓S度之中，但相比之下，在圖像左部區(qū)域，STC-NMF樹木與屋頂?shù)呢S度差異較為明顯一些，也相比之下略微貼近真實地物的豐度圖。

特別需要注意的是，不論在端元對比亦或是豐度對比，關(guān)于屋頂?shù)慕饣煨Ч琕CA略優(yōu)于這3種算法，這是優(yōu)于為了對比更加公正、更加明顯，實驗中STC-NMF，CNMF，L1/2-NMF 3種算法初始化，都是用N-Finder初始化端元，用FCLS初始化豐度；而VCA本身作為聚類方法可以為算法進(jìn)行端元初始化。在實驗中經(jīng)過多次對比，此場景的四端元解混中，進(jìn)行屋頂端元解混，大約在第50個波段開始N-Finder初始化端元開始有較為明顯的向上趨勢，而VCA接近于真實地物，有向下趨勢；而經(jīng)過約束非負(fù)矩陣分解的梯度下降，有效抑制了這種向上的趨勢，特別是STC-NMF最為有效，最為接近真實地物，但這種差異是由于初始化算法造成的，因此，STC-NMF也明顯受到初始化算法聚類效果的影響。而實驗中初始化選擇N-Finder而不是VCA是因為在實驗中，對于端元的初始化效果，不論是四端元或五端元亦或是六端元，N-Finder總體是優(yōu)于VCA的，同時也是為了減小初始化對算法的影響。

4 結(jié)束語

本文提出一種將米氏散射看作干擾的約束非負(fù)矩陣分解算法，散射項約束非負(fù)矩陣分解，并將其與幾種較為知名且被較多引用的算法進(jìn)行對比。通過實驗比較可以看出，這種算法在本組實驗的既定條件下，性能要略微具有優(yōu)勢，即是說，還可以通過對調(diào)節(jié)參數(shù)的深入研究進(jìn)一步使解混能力變優(yōu)。而這種算法從理論角度出發(fā)，STC-NMF將約束項由抽象的數(shù)學(xué)符號到具象的物理意義，在梯度下降過程中，使得所求的端元與豐度更貼近于現(xiàn)實的情況，且其勒讓德函數(shù)形式的數(shù)學(xué)解使得其在數(shù)學(xué)角度，與統(tǒng)計性、稀疏性等基于數(shù)學(xué)特性的約束也高度耦合，因此，STC-NMF應(yīng)具有優(yōu)秀的解混性能。通過實驗仿真，可以證明，STC-NMF也確實有著不錯的解混性能，相較于L1/2-NMF和CNMF，解混性能也具有一定的優(yōu)勢。