徐盈盈，黎建華

（臺州學(xué)院 電子與信息工程學(xué)院，浙江 臨海 317000）

0 引言

高光譜影像（hyperspectral imagery）是由衛(wèi)星等遙感設(shè)備上搭載的成像光譜儀獲取的三維數(shù)據(jù)，其中二維圖像描述了地表空間特征，而光譜波段揭示了圖像每一像元的光譜曲線，這些光譜涵蓋了紫外、可見光、近紅外和中紅外等波段［1］。高光譜影像具有很高的光譜分辨率（通常波段寬度＜10 nm），其豐富的光譜信息可充分反映不同地物的特征，在作物估產(chǎn)、礦物勘探、環(huán)境監(jiān)測、地形分類和軍事偵察等領(lǐng)域有著廣泛的應(yīng)用［2］。由于成像傳感器通常空間分辨率較低，數(shù)據(jù)采集過程受大氣、光照等干擾，且地物之間存在自然混合等因素，高光譜影像中普遍存在混合像元（mixed pixel），即一個成像單元的光譜是多種物質(zhì)光譜的混疊信號?；旌舷裨纸?，亦稱解混（unmixing），其目的是從混合光譜信號中分離純凈端元（endmember）光譜并估計各端元所占豐度（abundance）比例［3］?；旌舷裨纸饪梢允垢吖庾V影像的應(yīng)用從像元級提升至亞像元級［4］，對提高遙感應(yīng)用的精準(zhǔn)度、推動定量遙感技術(shù)的深入發(fā)展具有重要的科學(xué)意義。

基于先驗光譜庫的半監(jiān)督線性解混方法是近年來的研究熱點(diǎn)，利用事先獲取的過完備端元光譜庫，解混過程等價于從光譜庫中挑選端元并求其豐度［5］。由于光譜庫所含光譜數(shù)目通常遠(yuǎn)大于待解混的端元個數(shù)，該解混方法一般建模為稀疏誘導(dǎo)約束的回歸模型，并采用稀疏回歸算法進(jìn)行求解，因此通常稱此類方法為稀疏解混。目前稀疏解混的研究重點(diǎn)都在于挖掘待求豐度的先驗信息并設(shè)計相應(yīng)的豐度約束，如分別基于對豐度矩陣空間同質(zhì)性和行稀疏性的考慮，Iordache等人先后提出了經(jīng)典的SUnSAL-TV模型［6］與 CLSUnSAL 模型［7］。最近幾年，Wang 等人［8］提出了 CCSU 模型，通過在協(xié)同稀疏解混框架中引入集中約束來更精確地描述豐度圖的非局部冗余性。Huang等人［9］提出JSpBLRU模型，同時將影像塊結(jié)構(gòu)的聯(lián)合稀疏與低秩作為豐度估計約束。Li等人［10］提出了一種局部光譜相似性保留約束，以維持局部區(qū)域光譜相似性的方式來提高稀疏解混的魯棒性。Wang等人［11］提出DRSU-TV模型，通過使用雙權(quán)值來增強(qiáng)光譜和空間域的豐度稀疏性。Zhang等人［12］提出光譜-空間同時加權(quán)的稀疏解混模型，更好地保留了豐度圖的細(xì)節(jié)信息。

本文利用稀疏解混模型中重建數(shù)據(jù)的低秩特性，在解混過程中通過約束重構(gòu)影像的數(shù)據(jù)特征來減少噪聲等誤差對豐度估計的干擾，提出了重建約束與豐度約束協(xié)同的稀疏解混模型。與經(jīng)典解混模型的對比實(shí)驗表明，稀疏解混模型能有效提高解混精度。

1 稀疏解混建模

現(xiàn)有的稀疏解混方法充分利用了待求豐度的結(jié)構(gòu)特性，但沒有從影像數(shù)據(jù)的角度出發(fā)設(shè)計模型約束。本文發(fā)掘線性解混問題中理想的重建數(shù)據(jù)應(yīng)當(dāng)具有的低秩特性，以此為先驗，構(gòu)造了重建影像的低秩約束項，并考慮豐度矩陣的空間平滑特性，提出結(jié)合重建影像低秩約束和豐度TV約束的稀疏解混模型。

1.1 線性解混模型

線性混合假設(shè)指像元光譜是若干端元按照各自的豐度比例線性疊加得到的混合信號，該假設(shè)符合大部分應(yīng)用場景且易于建模，其矩陣形式可以表示為：

其中H∈?b×n表示高光譜影像數(shù)據(jù)的矩陣形式，b是光譜的波段數(shù)，n是影像包含的像元個數(shù)；W∈?b×m表示包含了m條純凈物質(zhì)光譜的端元矩陣或光譜庫；X∈?m×n是豐度矩陣，N∈?b×n為表示模型誤差與噪聲的擾動項。

在基于先驗光譜庫的解混問題中，給定的過完備光譜庫W一般包含成百上千種物質(zhì)的光譜，而實(shí)測影像H內(nèi)包含的e種端元通常只有幾個到幾十個，遠(yuǎn)遠(yuǎn)小于光譜庫的列維數(shù)，即e?m。因此，與冗余光譜庫對應(yīng)的豐度矩陣X理論上僅有e行非零元素，即X具有稀疏性。此外，鑒于豐度的物理屬性，X是滿足每列元素之和為1的非負(fù)矩陣?；诰€性混合模型（1），光譜庫先驗的解混通常建模為如下的稀疏回歸問題：

其中‖X‖0表示用l0范數(shù)統(tǒng)計X中的零元素個數(shù)，‖·‖F(xiàn)表示用于度量回歸誤差的Frobenius范數(shù)，?表示誤差容忍值。式（2）可以用正交匹配追蹤等貪婪算法直接求解［13］，或者將l0范數(shù)近似替換為更易求解的平滑連續(xù)函數(shù)或lp范數(shù)（0＜p≤1）后再優(yōu)化求解［14］。

1.2 基于低秩重建和TV正則的解混模型

不同于僅含豐度約束的稀疏解混模型，文中根據(jù)重建影像的數(shù)據(jù)特點(diǎn)構(gòu)造新的約束項，在SUnSALTV的基礎(chǔ)上提出了更魯棒的擴(kuò)展模型。Zhao等人［15］已經(jīng)論證了在無噪聲的理想情況下，線性混合的高光譜影像是低秩數(shù)據(jù)，而實(shí)際觀測影像由于噪聲等其他干擾往往是滿秩數(shù)據(jù)。模型中重建影像WX，本質(zhì)上正是由光譜庫W中的很少量光譜列以豐度X為權(quán)重線性混疊而成，因此理想的WX是具有低秩特性的矩陣，且WX的秩通常等于激活端元的個數(shù)。據(jù)此，將稀疏解混問題構(gòu)造為如下優(yōu)化方程：

其中rank(·)表示矩陣的秩，TV(·)表示矩陣的全變分，α,β,γ均為非負(fù)參數(shù)。

顯然式（3）是NP-hard的組合優(yōu)化問題，為了便于極小化求解，我們對目標(biāo)約束其進(jìn)行凸松弛近似，并將模型（3）改寫為如下無約束的凸優(yōu)化問題：

其中‖·‖*表示的核范數(shù)是矩陣秩的最小凸包絡(luò)函數(shù)，定義為‖WX‖*=∑iσi(WX)，σi(WX)表示W(wǎng)X矩陣經(jīng)降序排列的第i個奇異值；‖·‖1是矩陣l0范數(shù)的最緊凸松弛l1范數(shù)；τ+(X)表示豐度矩陣的非負(fù)約束，即X的元素值均為非負(fù)時有τ+(X)=0，否則τ+(X)=+∞。PX≡[PhX;PvX]是豐度矩陣的全變分正則的線性算子表示，其中PhX和PvX分別表示豐度圖的橫向和縱向差分。

模型（4）的核心思想就是通過增加對重建影像的結(jié)構(gòu)約束實(shí)現(xiàn)更合理的豐度估計。從理論上分析，目標(biāo)方程F(X)中對待估豐度矩陣X的全變分約束促使各端元對應(yīng)的豐度圖在相鄰像元之間為平滑過渡，符合自然物質(zhì)分布的局部平滑特性；對重建影像WX的低秩約束，使得該解混模型在數(shù)據(jù)重建時對噪聲有更好的魯棒性，進(jìn)一步引導(dǎo)豐度估計朝著有利于無噪聲影像重建的方向進(jìn)行優(yōu)化。

2 基于ADMM的數(shù)值優(yōu)化算法

優(yōu)化方程（4）可以采用交替方向乘子法（Alternating Direction Method of Multipliers,ADMM）［16］進(jìn)行求解。ADMM算法框架的基本思想是通過變量替換的方法，把難以整體優(yōu)化的原問題分裂為多個易于求解的子問題。我們引入V1～V5共五個輔助變量，將式（4）等價改寫為如下帶約束的優(yōu)化問題：

可以將式（5）寫成如下的增廣拉格朗日形式：

其中參數(shù)ρ是大于0的懲罰系數(shù)，D/ρ表示與約束GX+BV=0對應(yīng)的拉格朗日乘子。文獻(xiàn)［17］論證了對于形如式（6）的優(yōu)化問題，在矩陣G滿足列滿秩且方程g是閉凸函數(shù)的條件下，對任意ρ＞0，通過ADMM迭代極小化式（6）可以收斂到存在的解。

為便于迭代公式推導(dǎo)，我們將式（6）所示的增廣拉格朗日方程展開如下：

其中D1～D5分別是對應(yīng)于式（5）中五個等式約束的拉格朗日乘子。

根據(jù)ADMM方法，我們將對式（7）中的每個變量分別優(yōu)化，即在極小化某一個變量相關(guān)的方程時固定其余變量的值，依次交替迭代至達(dá)到算法的終止條件。令t表示迭代次數(shù)，t=0時所有待優(yōu)化的變量均初始化為全零矩陣。以t+1輪迭代為例，本文的模型優(yōu)化被分為以下幾個子問題的求解：

1）式（7）對豐度矩陣X求導(dǎo)并令導(dǎo)數(shù)等于0，可得X的迭代優(yōu)化公式為：

其中I表示單位矩陣。

2）重建影像的替換變量V1的優(yōu)化目標(biāo)是核范數(shù)約束的極小化問題，我們可以套用奇異值閾值（singular value thresholding，SVT）算法［18］，直接給出V1的求解公式：

記矩陣Z的奇異值分解為Z=UΣS*，則根據(jù)文獻(xiàn)［18］可以定義：

其中矩陣U和S的列分別是矩陣Z經(jīng)過奇異值分解后的左奇異向量和右奇異向量，S*表示酉矩陣S的共軛轉(zhuǎn)置；矩陣Σ的對角元素σi是矩陣Z的奇異值。

3）通過求導(dǎo)可得用于空間約束的豐度矩陣替換變量V2的更新公式：

其中卷積算子P可以利用離散傅里葉變換進(jìn)行快速計算。

4）豐度矩陣全變分約束的替換變量V3的目標(biāo)方程可以用軟閾值（soft thresholding）算法［19］進(jìn)行求解，其更新公式如下：

其中soft(·,η)表示以η為參數(shù)的軟閾值函數(shù)。

5）豐度稀疏約束的替換變量V4的迭代公式同樣可用軟閾值方法直接得到：

6）豐度非負(fù)約束的替換變量V4按如下公式將結(jié)果投影至非負(fù)象限：

7）拉格朗日乘子D1～D5的更新公式如下：

算法的終止條件是t達(dá)到最大迭代次數(shù)tmax或者相鄰兩次迭代的相對誤差滿在下文的實(shí)驗中，我們設(shè)定tmax=200,ξ=10-4。

3 實(shí)驗結(jié)果與分析

在模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)上分別測試文中模型與算法的解混效果。通過對比稀疏解混算法與經(jīng)典解混算法SUnSAL-TV和CLSUnSAL的實(shí)驗結(jié)果，驗證了結(jié)合重建數(shù)據(jù)低秩約束與豐度空間約束模型的有效性。實(shí)驗采用固定其他參數(shù)，單獨(dú)調(diào)節(jié)其中一個參數(shù)的方式交替調(diào)參，所有實(shí)驗方法的模型參數(shù)均選自統(tǒng)一參數(shù)集合Γ={ }0,0.0005,0.001,0.005,0.01,0.05,0.1,0.5,1,1.5,3,5,10,20，最終取使各算法達(dá)到最佳效果的參數(shù)組合，以實(shí)現(xiàn)公平比較。代碼的實(shí)現(xiàn)基于MATLAB R2018b平臺，硬件配置為英特爾i7-8750H處理器（2.20GHz）和16GB內(nèi)存。

3.1 模擬數(shù)據(jù)實(shí)驗

模擬實(shí)驗采用的光譜庫來自USGS公開的splib06a數(shù)字光譜庫［20］，該庫包含了在波長0.4～2.5微米范圍內(nèi)采集的240種純凈物質(zhì)的光譜信號，每條光譜記錄了224個波段的電磁波信號值。利用開源的模擬豐度圖生成代碼［21］，可以得到不同混合程度且滿足豐度值非負(fù)與豐度矩陣列和為1等物理屬性的豐度數(shù)據(jù)集。實(shí)驗以尺寸為64×64個像素且包含5個端元的模擬豐度為例，隨機(jī)生成的豐度數(shù)據(jù)如圖1所示，其中每個像素點(diǎn)的數(shù)值范圍為［0，1］。根據(jù)線性混合模型，先從光譜庫中隨機(jī)挑選若干端元光譜，再以隨機(jī)匹配的豐度圖作為權(quán)重比例線性疊加，即可合成三維的高光譜模擬混合數(shù)據(jù)。

圖1 模擬混合數(shù)據(jù)中5個端元對應(yīng)的真實(shí)豐度圖

由于模擬數(shù)據(jù)實(shí)驗已知真實(shí)的豐度圖，那么可以通過比較解混算法得到的估計豐度與真實(shí)豐度的差距，設(shè)計定量指標(biāo)，給出對解混效果的客觀評判。信號與重建比誤差（signal-to-reconstruction error，SRE）作為最常用的定量指標(biāo)，以分貝（dB）為度量單位的SRE定義如下：

其中Xtrue表示真實(shí)豐度矩陣，X是解混算法估計的豐度矩陣，E[·]表示數(shù)學(xué)期望。SRE指標(biāo)的數(shù)值越大，表示解混的豐度估計越精確。

為了模擬不同信噪比的高光譜影像，我們對混合數(shù)據(jù)添加了不同強(qiáng)度的高斯噪聲，得到信噪比分別為SNR=20，30，40dB的三組模擬數(shù)據(jù)。通過實(shí)驗參數(shù)調(diào)優(yōu)，分別為這三組數(shù)據(jù)設(shè)定模型參數(shù)三元組(α,β,γ)為(5,0.05,0.005),(1.5,0.01,0.005),(0.5,0.001,0.001)。三種實(shí)驗方法在三組模擬數(shù)據(jù)上得到的解混定量指標(biāo)結(jié)果，如圖2所示，可以看到本文方法對不同信噪比數(shù)據(jù)的解混均取得了高于對比方法的SRE值。對于SNR=20dB這樣的大噪聲數(shù)據(jù)，本文方法比SUnSAL-TV算法的結(jié)果高出了約3.5dB的SRE值，充分說明了引入低秩重建約束的有效性。在模擬數(shù)據(jù)的解混實(shí)驗中，定量指標(biāo)SRE值的高低已能夠客觀反映解混結(jié)果的好差。如圖3所示，以SNR=30dB的模擬數(shù)據(jù)的解混結(jié)果為例，直觀展示了文中的豐度估計結(jié)果與圖1所示真實(shí)豐度的相似性。

圖2 本文方法與對比算法對不同信噪比模擬數(shù)據(jù)解混得到的SRE（dB）值

圖3 本文算法對模擬數(shù)據(jù)解混得到的估計豐度圖

3.2 真實(shí)數(shù)據(jù)實(shí)驗

真實(shí)數(shù)據(jù)實(shí)驗采用公開的常用解混數(shù)據(jù)集：城市影像（Urban data）及其光譜庫［22］。實(shí)測城市影像數(shù)據(jù)的第80個光譜波段的場景，如圖4所示，該場景內(nèi)主要包含四種物質(zhì)，分別為屋頂、草地、樹木和瀝青路面。該影像數(shù)據(jù)的空間分辨率約為2 m×2 m，含有307×307個像元，每個像元記錄了162個有效的光譜波段值，實(shí)驗用的先驗光譜庫是包含上述四種物質(zhì)端元光譜的204條光譜集合。

圖4 實(shí)測高光譜數(shù)據(jù)

由于真實(shí)影像數(shù)據(jù)中各物質(zhì)的豐度未知，實(shí)驗結(jié)果可以根據(jù)豐度圖的視覺效果結(jié)合輔助定量指標(biāo)進(jìn)行評判。我們用影像的重建誤差（reconstruction error,RE）和豐度的稀疏度（sparsity of abundance,SA）作為評價指標(biāo)，定義：

其中影像重建誤差RE能反映用估計豐度和對應(yīng)端元光譜重建的混合數(shù)據(jù)與輸入影像的誤差；稀疏度SA的本意是計算豐度矩陣中0元素的占比，但理論0值在實(shí)際優(yōu)化中通常對應(yīng)于極小的非零數(shù)值解，因此我們用NUM(X＜θ)統(tǒng)計矩陣X中數(shù)值小于θ的元素個數(shù)，實(shí)驗中取θ=0.001。好的解混結(jié)果應(yīng)當(dāng)能同時滿足較低的重建誤差與較高的豐度稀疏度。由于真實(shí)數(shù)據(jù)實(shí)驗的兩個輔助指標(biāo)相互影響，那么以RE為首要調(diào)參指標(biāo)，即優(yōu)先選擇使RE值更小的參數(shù)組合，在有多組參數(shù)使得RE值相同的情況下，再選擇使SA值更大的為最優(yōu)參數(shù)。以此為參數(shù)調(diào)優(yōu)準(zhǔn)則，實(shí)驗中設(shè)定的模型參數(shù)(α,β,γ)為（0.01，0.005，0.01）。

三種實(shí)驗方法在城市影像數(shù)據(jù)上的解混指標(biāo)結(jié)果見表1所示?？傮w來看，本文方法取得了最低的重建誤差值，豐度稀疏度則介于兩個對比方法之間。CLSUnSAL方法的估計豐度雖然有最大的稀疏度，但其影像重建誤差遠(yuǎn)遠(yuǎn)大于本文方法。如圖5所示，從估計豐度的視覺效果上可以看到該城市影像中的四種主要物質(zhì)的豐度分布，從左到右分別對應(yīng)屋頂、草地、樹木和路面，與圖4所示的觀測場景基本一致。結(jié)合輔助評價指標(biāo)，可以說明本文方法對提升實(shí)測高光譜影像解混是有效性的。

表1 本文方法與對比算法對城市影像解混得到的重建誤差（RE）與豐度稀疏度（SA）指標(biāo)值

4 結(jié)語

在基于過完備光譜庫的解混問題中，我們利用重建影像具有低秩特性的先驗信息，引入重建數(shù)據(jù)的低秩約束來更好地引導(dǎo)豐度估計，再結(jié)合豐度矩陣的空間全變分約束構(gòu)建了具有良好抗噪性能的稀疏解混模型，在數(shù)值優(yōu)化上采用了基于ADMM的算法實(shí)現(xiàn)快速求解。通過與兩個經(jīng)典稀疏解混方法的對比，我們在高光譜模擬數(shù)據(jù)和實(shí)測影像上的解混結(jié)果均證實(shí)了本文方法能有效提高豐度估計的精度。