崔穎， 王恒， 朱海峰

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

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

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

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

1 線性稀疏解混模型

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

yi=Axi+ni

(1)

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

(2)

式中：‖xi‖0代表豐度向量xi中非零數(shù)值的個(gè)數(shù)；δ是根據(jù)噪聲和建模誤差設(shè)置的容忍誤差。由于P0問題是非凸優(yōu)化問題，正交匹配追蹤算法(orthogonal matching pursuit, OMP)是一種成熟的貪婪算法來解決P0問題[9]。字典A的稀疏度定義為A中最小線性相關(guān)的原子個(gè)數(shù)，它給定一個(gè)簡單而直接的方式計(jì)算P0問題的解的唯一性，然而計(jì)算一個(gè)矩陣的稀疏度與求解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中任意像元點(diǎn)x處沿著方向n的方向?qū)?shù)定義為：

(5)

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

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

(6)

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

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

(7)

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

(8)

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

(9)

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

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

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

2.2.1 SUnSAL-TV-STV解混模型

SUnSAL-TV-STV屬于一種“兩步法”算法，先通過SUnSAL-TV算法求解豐度矩陣，再利用STV正則項(xiàng)來探索豐度的分片平滑結(jié)構(gòu)。這個(gè)再優(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)化問題的一個(gè)簡單而有效的算法。對(duì)于給定的目標(biāo)函數(shù)(11)，有如下約束優(yōu)化問題：

(12)

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

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

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

(13)

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

3 合成數(shù)據(jù)實(shí)驗(yàn)

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

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

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

3.1 合成數(shù)據(jù)實(shí)驗(yàn)1

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

表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ù)實(shí)驗(yàn)2

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

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

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

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

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

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

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

表3 真實(shí)數(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。相對(duì)于SUnSAL-TV模型，該算法引入了結(jié)構(gòu)張量全變差空間正則項(xiàng)，該空間正則項(xiàng)可根據(jù)豐度向量所對(duì)應(yīng)的像元的空間幾何特征，自適應(yīng)地對(duì)SUnSAL-TV算法所求解的豐度矩陣進(jìn)行降噪處理，有效地改善其邊緣模糊與過平滑問題。

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

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