全變差Cauchy非負張量分解高光譜解混算法

吳新浪，葉 軍

(南京郵電大學 理學院，江蘇 南京 210023)

1 概 述

一幅完整的高光譜圖像(Hyperspectral Image，HSI)包含著豐富的空間信息和光譜信息，被廣泛地應用在多個研究領域[1-3]。然而，受遙感器空間分辨率和自然界地物表面特征復雜性的限制，HSI中存在大量的“混合像元”[4]，這大大降低了高光譜圖像地物識別的精度。因此，高光譜圖像解混(Hyperspectral Image Unmixing，HU)是高光譜圖像應用的一個重要步驟。

高光譜解混，是指從混合像元中分解出它們所包含的純物質特征(端元)并確定這些端元之間的比例(豐度)的過程。目前高光譜解混大多都是基于線性混合模型(Linear Mixed Model，LMM)[5]進行研究的。線性混合模型通常假設混合像元光譜是純物質特征的一個線性組合，并以相應的豐度分數加權，有著結構簡單，物理意義明確，適用性廣的優點。近年來，基于線性混合模型的各類方法中，非負矩陣分解(Nonnegative Matrix Factorization，NMF)方法[6]備受廣大研究者們的關注。NMF能夠將高維HSI數據分解成兩個非負矩陣的乘積，同時獲得端元矩陣和豐度矩陣。由于目標函數的非凸性，直接對HSI進行分解常常會陷入局部最小值。因此，常常在NMF框架中加入合適的約束以減小解空間，首先被添加的約束便是豐度非負約束(Abundance Nonnegativity Constraint，ANC)和豐度和為1約束(Abundance Sum-to-one Constraint，ASC)。為了進一步減小解空間，避免局部最小值問題，更多關于端元矩陣和豐度矩陣的約束被加入到NMF模型中。

根據HSI的幾何特征，Miao和Qi[7]將單純性體積最小化和NMF模型結合起來，提出了最小體積約束NMF(Minimum Volume Constrained NMF，MVCNMF)。Lu等[8]在NMF框架中加入了流形正則化項，提出圖正則化NMF(Graph Regularized L1/2-NMF，GLNMF)方法。孔繁鏘等[9]在解混模型中引入豐度矩陣的空間相關性約束，提出了一種基于空間相關性約束聯合子空間追蹤的HU算法。利用HSI中豐富的空間與光譜相關性，袁博等[10]提出了一種結合空間與譜間相關性分析的NMF解混算法。在促進豐度矩陣稀疏性方面，由于L0范數的非凸性不利于模型求解，L1范數最早被用來代替L0范數促進稀疏[11]。但是，L1范數不滿足完全可加性約束[12]，于是Qian等[13]、EJ等[14]、李璠等[15]分別利用L1/2范數、加權L1范數和變形L1范數代替L1范數促進豐度稀疏，均取得了不錯的解混效果。然而，當噪聲大量存在時，上述方法的解混效果顯著降低，這是因為基于NMF方法中的最小二乘目標函數對噪聲敏感。為了減少噪聲的影響，許多文獻引入了魯棒估計量來代替傳統的最小二乘目標函數，提出了許多魯棒NMF方法，如基于熵損失的魯棒NMF[16]、基于一般損失的魯棒NMF[17]等。文獻[18]利用閾值懲罰函數代替非負矩陣欠逼近中的殘差非負約束，提出了L1/2正則化的逐次高光譜圖像光譜解混，抗噪性能良好。但是，魯棒NMF方法，在解混過程中會丟失一部分空間信息，其主要的原因是，NMF是一種矩陣分解方法。張量提供了一個更自然的HSI立方體表示，它表征了兩個空間維和一個光譜維。矩陣向量非負張量分解算法(Matrix-Vector Nonnegative Tensor Factorization，MV-NTF)[19]將HSI分解為幾個分量張量，每個分量張量都是矩陣和向量的外積，分別代表端元和豐度。然而，MV-NTF缺少一些必要的約束條件。

為了能夠充分利用空間信息并同時減少噪聲的影響，該文在非負張量分解的框架上利用Cauchy損失函數來代替最小二乘損失函數，考慮HSI的稀疏性和豐度張量的分段光滑性，提出了全變差Cauchy非負張量分解高光譜解混算法(TV-CNTF)，該方法能夠保證HSI空間信息不會丟失，同時又能減小噪聲對HSI解混的影響。具體貢獻如下：

(1)利用張量分解的框架處理解混問題，充分保留了HSI的空間結構信息，由于HSI中相鄰的像素更可能有相近的物質構成比例，因此在張量分解模型下引入TV正則項來保證豐度的分段光滑性。

(2)利用Cauchy損失函數代替最小二乘損失函數，與最小二乘損失函數相比，Cauchy損失函數的影響函數有上界，這意味著HSI中的大噪聲點對解混結果的影響有限。因此，Cauchy損失函數可以一定程度上抑制噪聲對解混結果的影響。

(3)通過交替方向乘子法(ADMM)，推出了TV-CNTF的迭代更新算法，大量實驗結果表明了所提TV-CNTF解混方法的有效性。

2 背 景

2.1 線性混合模型

線性混合模型(Linear Mixture Model，LMM)通常假設觀測到的HSIy∈RI×J×L可以表示成對應的豐度系數和端元光譜的模積形式，即：

y=x×3A+N

(1)

其中，A∈RL×K是端元矩陣，x∈RI×J×K表示對應的豐度圖像，N∈RI×J×L為高斯噪聲和模型重構誤差，I、J和L分別是HSI的寬度、高度和波段數量，K是端元數。通常豐度還應滿足兩個約束，豐度非負約束(ANC)和豐度和為1約束(ASC)。分別表示為：

x≥0

(2)

x×311×K=1I×J×1

(3)

其中，11×K為全1行向量，1I×J×1為I×J×1的全1三階張量。

2.2 稀疏張量分解

由于LMM中假設HSI是由多個端元光譜特征按照它們對應的豐度分數線性組合而成，因此端元矩陣A和豐度張量x表現出非負性。建立非負張量分解模型如下：

(4)

其中，‖?‖F表示張量的Frobenius范數，三階張量x∈RI×J×K的Frobenius范數為：

(5)

由于HSI圖像中的混合像元往往由有限的幾個端元混合而成，因此豐度張量表現出稀疏性，建立非負稀疏張量分解模型如下：

(6)

其中，λ為控制稀疏程度的參數，‖x‖0表示張量x中非零元素的個數。

由于(6)是一個非凸問題，通常把L0范數松弛為L1范數，如下：

(7)

(8)

其中，⊙表示張量元素間的乘法，權值張量M的更新如下:

M(k+1)=1./max(|x(k)|,eps)

(9)

其中，eps為一個很小的正值，該文取eps=1e-12。

2.3 Cauchy損失函數

Cauchy損失函數(Cauchy Loss Function，CLF)定義為：

φ(x)=log(1+(x/c)2)

(10)

其中，c是常數。

對于一個損失函數φ(x)，影響函數ψ(x)是用來衡量φ(x)中元素變化帶來的影響。其中ψ(x)定義如下[20]：

(11)

那么CLF的影響函數為：

(12)

可以看出，隨著元素x的增大，ψ(x)存在上界，并最終趨于0，這意味著隨著噪聲的不斷增大，其對CLF的影響是有限的，最終達到上限。因此Cauchy損失具有魯棒性，它可以有效地抑制大噪聲。

3 全變差Cauchy非負張量分解模型

與最小二乘損失函數相比，Cauchy損失函數能夠減輕單個元素的影響，特別是噪聲較大的元素，因此Cauchy損失函數比最小二乘損失函數更能適應于噪聲環境下的解混。在NTF框架下，該文利用Cauchy損失函數來代替最小二乘損失函數，同時在模型中加入了TV算子來保證豐度張量的分段平滑結構，提出了一種新的解混算法。

3.1 Cauchy損失函數下的重構誤差

基于張量分解框架下的重構誤差的最小二乘損失函數形式如下：

(13)

其中，y::l表示張量y沿第三維(光譜維)的切片，el=‖Y::l-(x×3A)::l‖F表示第l片的片段誤差。

考慮到Cauchy損失函數對噪聲的魯棒性，將(13)的最小二乘損失替換成Cauchy損失，得到：

(14)

由于φ(x)是非凸和非線性的，很難直接求解，因此把φ(x)在x=0處進行泰勒展開，得到：

(15)

(16)

(17)

其中：

W=diag(w1,w2,…,wL)

(18)

(19)

3.2 TV-CNTF算法解混模型

為了減少噪聲的影響，并考慮豐度的分段光滑性，該文將Cauchy損失函數下的重構誤差式(17)及TV項重新引入模型(8)中，提出了全變差Cauchy非負張量分解模型，如下：

(20)

3.3 ADMM算法求解

針對模型(20)，為了便于求解，引入輔助變量L，模型(20)變為：

(21)

將約束L=x代入目標函數，模型(21)變為：

(22)

模型(22)轉化為增廣Lagrange函數如下：

(23)

其中，Ψ∈RL×K、Γ∈RK×(I*J)為矩陣形式的Lagrange算子。

最小化增廣Lagrange函數(23)，通過交替方向乘子法(ADMM)，尋找每個變量的解，同時修正其他變量。

(24)

子問題(24)可以通過KKT條件解出：

(25)

(26)

整理(26)，得：

(27)

(28)

其中，δ是ASC約束的參數，該文取δ=15。子問題(27)可以通過KKT條件解出：

(29)

(30)

問題(30)求解等價于問題(31)求解：

(31)

k=1,2,…,K

(32)

問題(32)可以通過FGP算法[21]解出：

(33)

綜合以上推導，得到全變差Cauchy非負張量分解算法(TV-CNTF)如下：

算法1：TV-CNTF算法。

輸入：HSIy，初始端元矩陣A0和豐度張量x0，端元數K，參數λ、μ、τ；

輸出：端元矩陣A和豐度張量x。

Step1初始化A,x,L；

Step2通過(18)計算W；

Step5通過(29)更新x；

Step6通過(33)更新L；

Step7重復Step4-Step6，直到滿足終止條件；

4 數值實驗和結果分析

為了驗證所提TV-CNTF方法的有效性，在模擬HSI數據集和真實數據集進行了實驗。選取沒有添加TV項的Cauchy非負張量分解算法(CNTF)(只在模擬數據實驗下對比)和解混方法GLNMF[8]、MV-NTF[19]、TV-RSNMF[12]進行對比實驗，并對實驗結果進行分析。選取光譜角距離(SAD)和均方根誤差(RMSE)對實驗結果進行評價分析，定義如下：

(34)

(35)

4.1 模擬數據實驗

從美國地質調查局(USGS)數字光譜庫中選取4個光譜特征，包含224個光譜波段，取其中的187個低噪波段使用LMM來生成模擬數據[12]。生成后的數據包含187個波段，像元數為48×48。模擬實驗中，對生成的HSI添加信噪比(Signal Noise Ratio，SNR)為[10,20,30,40](dB)的高斯噪聲，用這個數據來完成對比實驗。

表1和表2分別展示了不同噪聲水平下不同解混方法的SAD值和RMSE值。從表1和表2中可以看出，所提方法對不同SNR下的模擬數據均有很好的解混效果。綜合考慮RMSE值和SAD值，在SNR較低即SNR=10 dB時，TV-CNTF的解混效果明顯高于其他對比方法，解混效果依次比方法GLNMF、MV-NTF、TV-RSNMF和CNTF提升了約17.3%、26.9%、12.8%和12.5%，這說明引入的Cauchy損失函數在噪聲更大時解混性能更突出。與CNTF對比可以看出，TV項的引入可以提高TV-CNTF方法的解混性能，但在表1中SNR=10 dB時，CNTF的SAD值比TV-CNTF要小，這說明TV項的引入并不能很好地降低噪聲對解混的影響。在其他信噪比下，TV-CNTF方法的解混效果與其他方法相比也有一定的優勢。

表1 不同方法進行模擬數據實驗得到的SAD值

表2 不同方法進行模擬數據實驗得到的RMSE值

4.2 真實數據實驗

選取的真實數據是AVIRIS Cuprite數據集。該數據集由224個波段組成，波長范圍覆蓋了0.4 μm～2.5 μm，像元數為250×191。本次實驗中，去除了Cuprite數據集中的高噪聲波段和水吸波段(第1-3，104-113，148-167，221-224波段)，保留了剩下的187個波段進行真實數據實驗。同時，從USGS數字光譜庫中選取參考端元光譜特征。對真實數據實驗的所有方法，都使用頂點成分分析法(VCA)來初始化端元矩陣，然后使用最小二乘法(FCLS)來初始化豐度矩陣。為了保證實驗結果的可靠性，所有實驗都將重復10次。

表3中展示了各種常見解混方法的SAD值。從表中可以看出，提出的TV-CNTF方法的SAD值更小，與其他方法相比，解混效果提升了10%～30%，因此證明了該方法對真實數據實驗解混的有效性。圖1是真實數據實驗解混得到的端元光譜與USGS光譜庫中真實端元的對比圖。圖2～圖5展示了真實數據實驗下各方法解混得到的端元光譜與USGS光譜庫中真實端元光譜的對比圖和豐度圖，實驗中選擇了三種典型端元進行展示，分別是：Alunite、Buddingtonite和Chalcedony。

表3 不同方法在AVIRIS Cuprite數據集中提取不同物質的SAD值

從圖2～圖5可以看出，TV-CNTF和GLNMF得到的端元光譜圖與USGS光譜庫中的真實端元光譜擬合程度更好，從豐度圖上看，TV-CNTF方法比其他方法保留了更多的細節，圖像表現更加豐富。因此，可以得出結論：提出的TV-CNTF的解混效果更好。

圖1 真實數據實驗解混得到的端元光譜與USGS光譜庫中真實端元的對比圖

圖2 真實數據實驗下TV-CNTF得到的端元光譜圖和豐度圖

圖3 真實數據實驗下GLNMF得到的端元光譜圖和豐度圖

圖4 真實數據實驗下MV-NTF得到的端元光譜圖和豐度圖

圖5 真實數據實驗下TV-RSNMF得到的端元光譜圖和豐度圖

4.3 參數分析

首先對模擬數據實驗SNR=20 dB時的參數λ、μ和τ進行分析。給定λ的取值范圍為 [5e-4,1e-3,5e-3,0.01,0.05,0.1],μ的取值范圍為[1,10,1e2,1e3,1e4,1e5]，τ的取值范圍為[1e-4,1e-3,1e-2,1e-1]。固定μ=1e2，對參數λ和τ的不同取值進行解混實驗，結果如圖6所示。從圖6中可以看出，當兩個參數的值接近于0時，SAD和RMSE值顯著增加，這表明稀疏正則項和TV正則項的加入對實驗結果有著積極的影響。當λ=0.01，τ=0.01時，SAD和RMSE值同時達到一個很小的值，綜合考慮取λ=0.01，τ=0.01。再固定λ=0.01，對參數μ和τ的不同取值組合進行解混實驗，結果如圖7所示。從圖7中可以看出，當μ=1e2，τ=0.01時，SAD和RMSE值同時達到一個很小的值，綜合考慮取μ=1e2。因此，不同SNR下的模擬數據實驗的參數確定為λ=0.01，μ=1e2，τ=0.01。

圖6 模擬數據實驗在參數λ和τ的不同組合下的解混結果SAD和 RMSE

圖7 模擬數據實驗在參數μ和τ的不同組合下的解混結果SAD和 RMSE

通過多次真實數據實驗的結果，取μ=500，τ=0.01，λ=0.1，其他參數設置與模擬實驗相同。

5 結束語

該文提出了TV-CNTF算法用于HSI解混，充分保留了HSI的空間結構信息，TV項的引入保證了豐度張量的分段光滑性，同時使用Cauchy損失函數來代替最小二乘損失函數，通過減小噪聲點在解混模型中的權重，來降低噪聲對解混結果的影響。最后的模擬數據和真實數據實驗結果表明，無論是在目視效果還是定量評價指標上，該方法對圖像解混效果的提升都很顯著。然而，在HSI解混過程中，張量分解的方法遠不如非負矩陣分解的方法成熟，如何改進張量分解的方法使其更適合HSI圖像解混將會是未來重點研究方向。