崔穎， 王恒， 朱海峰

(哈爾濱工程大學(xué) 信息與通信工程學(xué)院，黑龍江 哈爾濱 150001)

光譜解混在高光譜圖像應(yīng)用中扮演著重要角色。解混算法依賴于預(yù)先假設(shè)的線性或非線性的光譜混合模型，由于線性混合模型在實際應(yīng)用中具有靈活、易于實現(xiàn)等優(yōu)點，當(dāng)前大部分光譜解混算法均基于線性混合模型(linear mixing model, LMM)[1]。稀疏解混算法將高光譜圖像中每個混合像元建模為預(yù)先可獲得的光譜庫中地物光譜的線性組合，其目的是在光譜庫中找到最優(yōu)的端元子集并求解該端元子集所對應(yīng)的豐度[2]。

由于參與混合像元的端元數(shù)量遠小于庫中涉及的原子數(shù)量，在求解稀疏解混欠定方程時，需要基于稀疏誘導(dǎo)正則項的有效稀疏回歸技術(shù)。Bioucas-Dias等[3]提出變量分離與增廣拉格朗日(sparse unmixing algorithm via variable splitting and augmented Lagrangian, SUnSAL)算法，該算法通常假設(shè)參與每個像元的端元數(shù)量較少，盡管該方法具有低復(fù)雜度，但它沒有在解混模型中考慮端元的分布情況，解混性能較差。Iordache等[4]提出的協(xié)同變量分離與增廣拉格朗日(the collaborative SUnSAL, CLSUnSAL)算法是在假設(shè)高光譜圖像中的所有像元共享相同的活越端元的情況下開發(fā)的。這意味著如果光譜庫原子的豐度被收集在矩陣中，則豐度矩陣中應(yīng)該只有少數(shù)幾個非零行。但CLSUnSAL算法卻沒有考慮到端元總是出現(xiàn)在空間均勻區(qū)域而并非均勻分布在整個高光譜圖像場景中，根據(jù)這一分析，Zhang等[4]等提出了局部協(xié)同稀疏回歸算法。研究表明，稀疏解混中包含的空間信息對估計的豐度的準(zhǔn)確性有積極的影響[5]。Iordache等[6]提出的全變差正則化變量分離與增廣拉格朗日(sparse unmixing via variable splitting augmented Lagrangian and total variation, SUnSAL-TV)算法將空間信息作為全變差(total variation, TV)正則項應(yīng)用于稀疏解混數(shù)學(xué)模型。 該TV項假設(shè)2個空間相鄰的像元對同一種地物端元具有相似的豐度。然而，地物分布總是充滿了復(fù)雜性，對于那些地物分布邊界上的像元來說，在TV正則項的作用下，SUnSAL-TV算法可能會產(chǎn)生過平滑與邊緣模糊的現(xiàn)象。

Lefkimmiatis等[7]提出的結(jié)構(gòu)張量全變差(structure tensor total variation, STV)正則項族能夠很好地描述圖像內(nèi)局部鄰域周圍的變化度量，將該正則項族應(yīng)用于灰度和矢量值圖像的去噪與去模糊中已取得了很顯著的效果。受此啟發(fā)，本文提出了一種結(jié)構(gòu)張量全變差再優(yōu)化變量分離與增廣拉格朗日(SUnSAL-TV-STV)的稀疏解混算法。該算法通過將STV正則項引入到SUnSAL-TV模型中約束求解的豐度矩陣，自適應(yīng)地校正相鄰像元所對應(yīng)的豐度向量之間的平滑性，該算法具有良好的魯棒性與解混性能。

1 線性稀疏解混模型

LMM中每個像元的光譜由端元光譜的線性混合來近似表達。在這種情況下，令L維向量y表示具有L個光譜帶的高光譜圖像的像元向量。觀測像元yi可以以光譜庫A中的光譜特征的線性組合的形式給出：

yi=Axi+ni

(1)

由于在光譜庫中只有少數(shù)的端元原子參與觀測像元光譜的混合，式(1)中的豐度向量xi是稀疏的。通過上面的這些定義，可以將解混問題寫成P0問題：

(2)

式中：‖xi‖0代表豐度向量xi中非零數(shù)值的個數(shù)；δ是根據(jù)噪聲和建模誤差設(shè)置的容忍誤差。由于P0問題是非凸優(yōu)化問題，正交匹配追蹤算法(orthogonal matching pursuit, OMP)是一種成熟的貪婪算法來解決P0問題[9]。字典A的稀疏度定義為A中最小線性相關(guān)的原子個數(shù)，它給定一個簡單而直接的方式計算P0問題的解的唯一性，然而計算一個矩陣的稀疏度與求解P0問題同樣困難。已經(jīng)證明，在受限等距特性(restricted isometry property, RIP)的某種條件下，豐度向量的0范數(shù)即‖xi‖0可以松弛到1范數(shù)[10]。松弛到1范數(shù)的P1問題為凸優(yōu)化問題，可表述為：

(3)

Y=AX+N

(4)

2 結(jié)構(gòu)張量全變差正則化解混模型

2.1 結(jié)構(gòu)張量全變差

這種方法將豐度矩陣X轉(zhuǎn)換到Sobolev空間W1,2(R2，Rm)。令n為任意二維空間內(nèi)的單位方向向量，豐度X中任意像元點x處沿著方向n的方向?qū)?shù)定義為：

(5)

式中：JX是X的雅克比(Jacobian)矩陣，定義為：

JX(x)=[X1(x),…,Xm(x)]

(6)

方向?qū)?shù)的大小‖?X/?n‖2描述了豐度X中像元點x沿n的變化量。為了更加有效地捕獲與x相鄰的像元的豐度差異，本文計算3×3窗口中‖?X/?n‖2的加權(quán)均方根(root mean square, RMS)，其中像素點x位于該窗口的中心。RMSK{‖?X/?n‖2}就是所謂的方向變化量：

式中：K(x)代表非負旋轉(zhuǎn)對稱的高斯卷積核；SKX代表像元點x所對應(yīng)的結(jié)構(gòu)張量，可表示為：

(7)

式中：SKX(x)是一個實對稱2×2矩陣，令λ+=λ+(SKX(x))，λ-=λ-(SKX(x))為SKX(x)的2個特征值且有λ+≥λ-，θ+和θ-為2個特征值所對應(yīng)的單位特征向量。令ω∈(-π,π]表示方向向量n與特征向量θ+之間的夾角，利用矩陣的特征分解對式(7)進行展開分析，可以得出2個重要結(jié)論：

(8)

在式(8)的基礎(chǔ)上，引入基于補丁的Jacobian算子[11], 則STV正則項的離散形式可以表示為：

(9)

式中：JK是扮演著線性映射JK:RNx×Ny×m→χ(Nx×Ny=N代表圖像平面)的基于補丁的Jacobian算子，χ?RNx×Ny×(G×m)×2，G是包含在高斯卷積核內(nèi)的元素的總數(shù)，Sp代表p維的Schatten范數(shù)[12]：

式中σn是矩陣M(M∈RN1×N2)的第n個奇異值。

2.2 解混模型與ADMM優(yōu)化算法

2.2.1 SUnSAL-TV-STV解混模型

SUnSAL-TV-STV屬于一種“兩步法”算法，先通過SUnSAL-TV算法求解豐度矩陣，再利用STV正則項來探索豐度的分片平滑結(jié)構(gòu)。這個再優(yōu)化的過程可表示為豐度降噪模型[13]：

(10)

綜上，本文提出SUnSAL-TV-STV求解優(yōu)化問題：

(11)

2.2.2 ADMM優(yōu)化算法

乘子交替法(the alternating direction method of multipliers, ADMM)是求解上述優(yōu)化問題的一個簡單而有效的算法。對于給定的目標(biāo)函數(shù)(11)，有如下約束優(yōu)化問題：

(12)

μ是一個正數(shù)，D/μ是關(guān)聯(lián)于限制GU+BV=0的拉格朗日乘子。V、G、B可表示為：

V=(V1,V2,V3,V4,V5)

至于優(yōu)化問題：

(13)

本文采用文獻[7]中介紹的迭代方案求解優(yōu)化問題(13)。ADMM算法中的停止迭代條件設(shè)置為‖GU(k)+BV(k)‖F(xiàn)≤ε，該優(yōu)化算法的收斂速度依賴于選擇一個合適的參數(shù)μ。試驗仿真將采用文獻[14-15]中描述的一個自適應(yīng)方案來自動調(diào)整參數(shù)μ。

3 合成數(shù)據(jù)實驗

在2個不同的合成數(shù)據(jù)集上，我們通過實驗來主要對比SUnSAL、CLSUnSAL、SUnSAL-TV算法以及本文提出的SUnSAL-TV-STV算法的解混性能。所有對比算法根據(jù)在調(diào)整正則項參數(shù)后各算法的信號重建誤差(signal-to-reconstruction error, SRE)值改變量不超過0.001時來選擇最佳參數(shù)。所有對比試驗均使用MATLAB R2012b在配備G2030 CPU(3.2 GHz)和2 GB RAM臺式機上實現(xiàn)。

SRE的值越大，解混結(jié)果越精確。此外，合成數(shù)據(jù)試驗還采用了成功概率(Ps)的指標(biāo)來評估解混性能。Ps定義為：

這是對相對誤差功率小于等于某個閾值T的概率估計[16]。為了縮短算法的運行時間，本文采用HySime(hyperspectral signal subspace estimation)算法[17]來估計試驗數(shù)據(jù)的信號子空間，在此基礎(chǔ)上采用多信號分類(multi-signal classification, MUSIC) 算法來修剪光譜庫以識別光譜庫的子集，該子集包含端元信號[18]。

3.1 合成數(shù)據(jù)實驗1

表1給出的是不同SNR環(huán)境下各解混算法獲得的SRE(dB)值。圖1對30 dB(DC1)下各解混算法的解混成功率進行了可視化對比。

表1 DC1下各算法所得最優(yōu)SRE值及參數(shù)Table 1 Optimal SRE values obtained by various algorithms under DC1 and the parameters

圖1 30 dB DC1下不同算法的解混成功率曲線Fig.1 Unmixing successful probability for DC1 with SNR=30 dB

3.2 合成數(shù)據(jù)實驗2

第2個實驗使用的光譜庫是由USGS庫中隨機選擇的498種地物光譜生成(A2∈R224×498)。本實驗從光譜庫A2中隨機選擇9個光譜原子作為端元信號并根據(jù)Drichlet分布模擬豐度矩陣來合成具有75×75個像元的合成數(shù)據(jù)集DC2。在生成的DC2中同樣加入了不同SNR(30、40、50 dB)的高斯白噪聲。此外，該實驗修剪光譜庫后保留了20個光譜原子。表2給出了不同SNR環(huán)境下各解混算法獲得的SRE(dB)值。圖2描繪了30 dB(DC2)下各解混算法的解混成功率曲線。

3.3 合成數(shù)據(jù)實驗分析與算法復(fù)雜度分析

本文統(tǒng)計了SNR=20 dB的DC1實驗中所有解混算法的計算處理時間，對應(yīng)于SUnSAL、CLSUnSAL、SUnSAL-TV以及SUnSAL-TV-STV算法的處理時間分別為0.590 2、15.317 9、15.103 8、53.650 0 s。

SUnSAL算法的計算復(fù)雜度是最低的，但由于該算法并沒有考慮到圖像所包含的空間先驗信息，解混精度也是所有算法中最低的。CLSUnSAL算法只是簡單地假設(shè)所有像元共享光譜庫中少數(shù)活躍的端元集，對空間信息的利用仍然不夠充分。SUnSAL-TV-STV算法由于引入了新的空間正則項而提高了計算復(fù)雜度，但該算法采用STV正則項約束SUnSAL-TV算法所求解的豐度矩陣，自適應(yīng)地修正了豐度矩陣的過平滑與邊緣模糊問題，因此該算法能夠以更高的可信度估計端元子集以及豐度矩陣。對比表1和表2中給出的SRE值可有力支持這一理論分析，同時也驗證了本文所提出算法的優(yōu)越性。此外，從圖1與圖2中可以觀察到，在同一閾值T下SUnSAL-TV-STV算法具有最高的解混成功率。

表2 DC2下各算法所得最優(yōu)SRE值及參數(shù)Table 2 Optimal SRE values obtained by various algorithms under DC2 and the parameters

圖2 30 dB DC2下不同算法的解混成功率曲線Fig.2 Unmixing successful probability for DC2 with SNR=30 dB

4 真實高光譜數(shù)據(jù)實驗

使用AVIRIS Cuprite開源數(shù)據(jù)集來進行真實高光譜數(shù)據(jù)試驗。該實驗使用的是一個非常具有代表性的250×191的像元子集，224個光譜帶的光譜波長均勻分布在0.4～2.5 μm內(nèi)。本實驗移除了受吸水性與較低信噪比干擾的36個光譜帶。該實驗使用3.1節(jié)的數(shù)字光譜庫A1作為本實驗的光譜庫，并相應(yīng)的移除所對應(yīng)的36個光譜帶。此外，為提高解混精度，避免計算微小的誤差值，將大于0.001的豐度估計值設(shè)為非零豐度并采用HySime-MUSIC算法人工修剪光譜庫A1以保留54個光譜原子。最后要說明的是，在實際高光譜圖像應(yīng)用中，測試圖像與預(yù)先獲得的光譜庫中的端元光譜之間由于成像條件不同總是存在一些差異[20]，但本文將其視作統(tǒng)一控制因素，不予討論上述差異對解混結(jié)果造成的影響。

真實數(shù)據(jù)試驗采用原始高光譜圖像與重建圖像(通過被修剪光譜庫及其相應(yīng)的豐度矩陣重建)之間的絕對重建誤差(ARE)來評估解混性能。ARE定義為：

圖3 由SUnSAL，CLSUnSAL，SUnSAL-TV和SUnSAL-TV-STV算法估計的AVIRIS Cuprite 礦區(qū)250×191像元子集的豐度Fig.3 Fractional abundance maps estimated by SUnSAL, CLSUnSAL, SUnSAL-TV, and SUnSAL-TV-STV(from top to bottom)for the considered 250×191 pixel subset of the AVIRIS Cuprite scene

從表3中各解混算法計算所得ARE值的對比中可以看出，SUnSAL-TV-STV算法具有最低的重建誤差。此外，上述各解混算法估計的每個像元中豐度值大于0.01的平均端元數(shù)分別為5.041 8、5.091 4、5.437 3、5.520 9。這些微小的差異反映出SUnSAL-TV-STV能夠使用光譜庫中少量的端元光譜來解釋每個混合像元，從而強調(diào)了解的稀疏性。定量地評估真實高光譜數(shù)據(jù)的解混結(jié)果是不充分的，還可以定性地觀察到在明礬巖的豐度估計圖中SUnSAL-TV-STV算法自適應(yīng)地修正了由SUnSAL-TV算法產(chǎn)生的過平滑與邊緣模糊問題。綜上所述，可以得出結(jié)論，SUnSAL-TV-STV算法是一種包含空間信息的有效的高光譜圖像稀疏解混算法。

表3 真實數(shù)據(jù)下各算法所得最優(yōu)ARE值及參數(shù)Table 3 Optimal ARE values obtained by various algorithms under the real hyperspectral data and the parameters

5 結(jié)論

1)本文提出了一種包含空間信息的高光譜圖像稀疏解混算法SUnSAL-TV-STV。相對于SUnSAL-TV模型，該算法引入了結(jié)構(gòu)張量全變差空間正則項，該空間正則項可根據(jù)豐度向量所對應(yīng)的像元的空間幾何特征，自適應(yīng)地對SUnSAL-TV算法所求解的豐度矩陣進行降噪處理，有效地改善其邊緣模糊與過平滑問題。

2)與現(xiàn)階段的其他稀疏解混算法相比，該算法在合成數(shù)據(jù)與真實高光譜數(shù)據(jù)的試驗仿真中均獲得了較好的解混性能。

雖然該算法在解混性能上得到了提高，但由于添加了新的空間正則項而犧牲了算法復(fù)雜度，未來工作將進一步改進算法以降低算法復(fù)雜度。此外，受空間非局部相似性的啟發(fā)，未來工作也將考慮在解混模型中引入能夠更充分挖掘圖像空間信息的非局部結(jié)構(gòu)張量全變差[21]正則項來有效提高解混性能。