一種自適應(yīng)魯棒最小體積高光譜解混算法

王天成 劉相振 董澤政 王海波

一種自適應(yīng)魯棒最小體積高光譜解混算法

王天成1劉相振1董澤政1王海波1

對(duì)高光譜圖像解混的目的在于從低空間分辨率的高光譜圖像中找到端元與對(duì)應(yīng)的豐度.本文根據(jù)解混算法中的最小體積準(zhǔn)則,提出了一種自適應(yīng)魯棒最小體積高光譜解混算法(Robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing,RMVHU).本算法通過引入負(fù)數(shù)懲罰正則項(xiàng),替換了同類算法中的豐度非負(fù)性約束(Non-negativity constraint,ANC),使算法對(duì)圖像中的噪聲與異常值具有更強(qiáng)的魯棒性;采用循環(huán)最小化方法,將非凸優(yōu)化問題分解為凸優(yōu)化子問題,然后應(yīng)用交替方向乘子法解決隨著像素點(diǎn)個(gè)數(shù)增大帶來的求解困難問題;對(duì)于正則項(xiàng)系數(shù),本算法提出了一種自適應(yīng)調(diào)整策略,提高了算法的收斂性,并且通過定性分析,說明了該調(diào)整方法的合理性.將算法應(yīng)用于合成數(shù)據(jù)與實(shí)際數(shù)據(jù),實(shí)驗(yàn)結(jié)果表明,與同類算法相比,本文提出的算法能夠取得更為優(yōu)秀的效果.

高光譜解混,交替方向乘子法,凸優(yōu)化,最小體積,自適應(yīng)估參

隨著遙感傳感器的飛速發(fā)展,高光譜圖像由于包含豐富的光譜間信息與空間信息,在分類、探測等方面得到了更加廣泛的應(yīng)用.與普通圖像不同,在高光譜圖像中,每一個(gè)像素點(diǎn)(像元)包含了成百上千個(gè)波段的反射率信息[1?2].高光譜的高譜間分辨率,為像素級(jí)乃至亞像素級(jí)的圖像分類與探測提供了條件,例如文獻(xiàn)[3]就利用高光譜的空譜特性進(jìn)行了有效的異常探測.在光譜圖像中,若某一像元中只存在一種物質(zhì),我們將這樣的像元稱為純像元[4].但是高光譜圖像的空間分辨率較低,圖像中存在著大量混合像元[4].混合像元廣泛存在于高光譜圖像中,是影響遙感分類精度和目標(biāo)探測效果的重要因素之一.為了進(jìn)一步利用高光譜數(shù)據(jù),我們需要通過分解混合像元得到圖像中一系列基本物質(zhì)(這樣的基本物質(zhì)被稱作端元)的光譜信息,同時(shí)還要求取端元在混合像元中的占比(豐度),而上述通過原始高光譜圖像得到端元與豐度的過程就是高光譜圖像的解混過程[5].

近年來,學(xué)者們提出了一系列基于線性混合模型的解混算法.在線性混合模型中,像元可以由一系列不相關(guān)的端元線性表示,且豐度滿足“非負(fù)性”(ANC)與“和為1”(ASC)兩個(gè)約束條件[6].在線性模型的基礎(chǔ)上,基于凸面幾何學(xué)的解混算法被廣泛研究.該類算法的主要思想是:根據(jù)線性模型的全約束條件,所有數(shù)據(jù)點(diǎn)均被包含在某一類單形體中,在這些單形體中,由端元作為頂點(diǎn)構(gòu)成的單形體體積是最小的.同時(shí)也有較多學(xué)者研究基于非負(fù)矩陣的解混算法[7?10]與稀疏解混算法[11?15].相較于其他算法,凸面幾何學(xué)類解混算法運(yùn)算速度快,解混精度較高,因此本文將繼續(xù)深入研究凸面幾何學(xué)的算法.

在基于凸面幾何學(xué)的算法中,有一類基于純像元假設(shè)的算法,如 VCA[16]、PPI[17]、NFINDR[18]等.雖然這些算法在圖像滿足純像元假設(shè)時(shí)[5]有不俗的表現(xiàn),但是若圖像中存在大量高度混合的數(shù)據(jù),解混精度將大大下降.為解決上述問題,學(xué)者們提出了一類基于最小體積變換的算法,該類算法能夠在不滿足純像元假設(shè)的情況下取得較好的端元提取效果,如 MVES[19]、MVSA[20]、SISAL[21]等.但是由于MVES、MVSA這類算法須嚴(yán)格滿足豐度非負(fù)約束,觀測數(shù)據(jù)的異常值將會(huì)對(duì)解混精度造成很大的影響;而SISAL雖然通過懲罰項(xiàng)放寬了約束條件,但是在求解優(yōu)化問題時(shí),SISAL存在兩個(gè)問題:1)將目標(biāo)函數(shù)做了二次近似,將導(dǎo)致近似的目標(biāo)函數(shù)由于高階項(xiàng)的缺失,在某些情況下會(huì)大尺度地偏離原目標(biāo)函數(shù);2)二次項(xiàng)的Hesse矩陣采用對(duì)角矩陣近似而非目標(biāo)函數(shù)的Hesse矩陣.上述問題的存在將給解混精度的提高帶來一定的限制.

為了改善上述算法出現(xiàn)的問題,本文提出了一種自適應(yīng)魯棒最小體積解混算法,該算法有如下優(yōu)點(diǎn):

1)將約束條件放寬,通過引入負(fù)數(shù)懲罰正則項(xiàng)替代豐度非負(fù)性約束,獲得了更強(qiáng)的抗干擾性能.

2)求解的目標(biāo)函數(shù)根據(jù)行列式的形式展開,變換之后的目標(biāo)函數(shù)與原函數(shù)完全等價(jià),不再是二次近似,解決了二次逼近帶來的誤差.

3)為了能夠解決懲罰正則項(xiàng)帶來的大規(guī)模凸優(yōu)化求解困難問題,應(yīng)用交替方向乘子法[22](ADMM),獲得了理想的計(jì)算精度與速度.

4)創(chuàng)造性地提出了一種自適應(yīng)的正則系數(shù)調(diào)整方法,提高了算法的收斂性與穩(wěn)定性,并通過定性的分析說明了此自適應(yīng)參數(shù)調(diào)整方法的合理性.

5)對(duì)于ADMM中的懲罰系數(shù),應(yīng)用文獻(xiàn)[22]中自適應(yīng)調(diào)節(jié)的方法,進(jìn)一步加快了算法的收斂速度.

文章結(jié)構(gòu)主體如下所述.在第1節(jié)中詳細(xì)介紹線性混合模型,凸面幾何學(xué)的數(shù)學(xué)基礎(chǔ)以及最小體積解混模型;在第2節(jié)中講述自適應(yīng)魯棒最小體積高光譜解混算法(RMVHU)模型,參數(shù)的自適應(yīng)調(diào)整方法,RMVHU算法步驟以及算法復(fù)雜度與收斂性;在第3節(jié)中講述實(shí)驗(yàn)用的數(shù)據(jù)來源,并將數(shù)據(jù)應(yīng)用于算法的實(shí)驗(yàn)與分析;在第4節(jié)中,總結(jié)實(shí)驗(yàn)的分析結(jié)果,得出結(jié)論.

1 幾何解混相關(guān)工作

1.1 高光譜線性混合模型

通常情況下,在線性混合模型中,高光譜圖像中的每個(gè)像元都可被近似認(rèn)為是圖像中各個(gè)端元的線性組合:

式中,p為端元個(gè)數(shù),yyy∈RL為任意像元的L維光譜向量(L為圖像波段數(shù)),為大小是L×p的端元矩陣為端元的向量表示.為豐度的向量表示,si表示像元中端元所占的比例,為誤差項(xiàng).

線性混合模型一般可分為3種情形:式(1)為無約束的線性混合模型;加上約束條件(2)則為非負(fù)約束混合模型;再加上約束條件(3)則為全約束混合模型.線性解混就是提取端元,同時(shí)求出各個(gè)端元在像元中所占的比例,得到端元豐度的過程.

在研究幾何學(xué)高光譜解混算法之前,本文將依次介紹相關(guān)數(shù)學(xué)符號(hào)的定義、求解凸優(yōu)化問題的著名算法框架——ADMM,以及仿射集與凸包的數(shù)學(xué)概念.

本文應(yīng)用到的符合及其意義如表1所示:

表1 數(shù)學(xué)符號(hào)及其意義Table 1 Mathematical notations and their meaning

1.1.1 交替方向乘子法(ADMM)

考慮如下線性約束優(yōu)化問題：

式中,μ為與收斂速度相關(guān)的常值,ddd為尺度對(duì)偶變量(Scaled dual variable)[22].

采用ADMM 的框架,能夠在一步步的迭代過程中逼近凸問題的最優(yōu)解.

1.1.2 仿射集與凸包

仿射包是一個(gè)仿射集,所以也能被表示為:

在文獻(xiàn)'[23]中,通過“求解式(10)描述的問題,根據(jù)向量集與k,得到一組參數(shù)(C,d)來表示仿射集.

式中,約束CTC=Ik是為了滿足rank(C)=k的秩約束條件,是投影到仿射集上的誤差,定義為:

問題(10)有如下閉合解的形式[23]:

凸包中的點(diǎn)x若不能被表示為嚴(yán)格凸組合的形式,即: 若且θ/=eeei,?i=1,···,p,稱 x為凸包的頂點(diǎn).若L=p?1且是仿射無關(guān)的,則稱凸包為單形體,向量集為單形體頂點(diǎn)的集合.

仿射集、凸包的概念如下圖所示:

圖1 頂點(diǎn)個(gè)數(shù)為3時(shí),凸包與仿射集的概念說明Fig.1 Illustration of convex hull and aきne hull when the number of vertices is three

從圖1中可以看到,由三個(gè)頂點(diǎn)描述的仿射集是一個(gè)平面,而凸包則是在仿射集平面上,由頂點(diǎn)組成的三角形內(nèi)的區(qū)域,這個(gè)三角形也被稱為二維單形體.是

1.2 最小體積解混模型

一般來說,式(1)中的端元矩陣是列滿秩的,即各個(gè)端元線性無關(guān).又由于式(3)“和為1”的假設(shè),文獻(xiàn)[19]指出,由觀測數(shù)據(jù)集組成的仿射包與由端元組成的仿射包是同一個(gè).因此與都能由同一組參數(shù)(C,ddd)來表示,其中:

考慮如下優(yōu)化問題:

式中,V(β1,···, βp) 代表由β1,···, βp作為頂點(diǎn)構(gòu)成的單形體的體積,代表光譜數(shù)據(jù)的低維表示.文獻(xiàn)[19]指出,在純像元假設(shè)下,問題(18)的最優(yōu)解如下:

同時(shí)文獻(xiàn)[19]指出,在純像元假設(shè)不嚴(yán)格成立的情況下,采用式(19)表示的可行解也接近問題(18)的最優(yōu)解.根據(jù)文獻(xiàn)[24],單形體的體積能夠被表示為:

將式(20)與式(21)代入式(18),式(18)描述的優(yōu)化問題能夠被等價(jià)表示為:

式(23)即為最小體積解混算法模型的數(shù)學(xué)描述.

2 自適應(yīng)魯棒最小體積高光譜解混算法

2.1 RMVHU算法模型

在凸面幾何學(xué)的最小體積類算法中,我們通常要求單形體能夠包圍所有的數(shù)據(jù)點(diǎn)集,即必須滿足豐度系數(shù)非負(fù)性條件.然而,在實(shí)際情況下,在圖像中很有可能存在異常點(diǎn),或者由于噪聲的影響導(dǎo)致像元分布在單行體之外.若在這些情況下仍強(qiáng)行滿足豐度非負(fù)性條件,則可能出現(xiàn)為了將這些單形體之外的點(diǎn)包入估計(jì)的單形體中,所估計(jì)的單形體頂點(diǎn)與真實(shí)的頂點(diǎn)相差甚遠(yuǎn)的情況.因此,為了容忍異常點(diǎn)給端元估計(jì)帶來的影響,我們需要構(gòu)造一個(gè)負(fù)數(shù)懲罰正則項(xiàng)替換非負(fù)性約束,從而在盡可能保證單形體體積最小的情況下,增強(qiáng)算法對(duì)異常值、噪聲值的魯棒性.同時(shí)通過引入此正則項(xiàng),減少在求解非凸優(yōu)化問題時(shí),陷入局部最小解的概率.

考慮引入負(fù)數(shù)懲罰正則項(xiàng):式中,sij為矩陣S的元素.該正則項(xiàng)能夠?qū)ω?fù)系數(shù)進(jìn)行懲罰,而對(duì)于非負(fù)系數(shù)則沒有任何影響.

式(24)中的h(sij),可以寫成如下形式:

將此正則項(xiàng)代替優(yōu)化問題(23)中的非負(fù)約束項(xiàng),建立如下優(yōu)化模型:

將式(25)代入式(26),上述問題等價(jià)為如下無約束非凸最優(yōu)化問題:

式(27)描述的問題即為RMVHU算法的問題模型.但是式(27)中的目標(biāo)函數(shù)非凸,很難求解得到全局最優(yōu)解,因此,可以采用循環(huán)最小化思想,逐行更新矩陣變量H與的值,以將此問題轉(zhuǎn)換為帶有絕對(duì)值的凸優(yōu)化問題,隨后通過代數(shù)余子式展開行列式,且脫去|det(H)|項(xiàng)的絕對(duì)值,將凸優(yōu)化問題轉(zhuǎn)換為兩個(gè)等價(jià)的凸優(yōu)化子問題.

具體的,考慮更新矩陣變量H與的第i行值,得到RMVHU算法的凸優(yōu)化子問題p?與q?:

式(31)中,Hi,j為H中元素hi,j的代數(shù)余子式.

2.2 RMVHU算法子問題的求解

問題(28)與(29)為帶有l(wèi)1范數(shù)的最優(yōu)化問題.下面介紹一種求解此類的方法.

定理 1.問題(33)的最優(yōu)解與問題(34)的最優(yōu)解相同:

則目標(biāo)函數(shù)值有:

只有當(dāng)y0的第i個(gè)元素y0(i)=0或者z0的第i個(gè)元素z0(i)=0時(shí),能夠滿足在當(dāng)前可行解 x0下,目標(biāo)函數(shù)取得最小值.根據(jù)此結(jié)論,結(jié)合約束條件y ? z =A x + b ,有:

式中ai為矩陣A的行向量,bi為向量b的第i個(gè)元素.所以,問題(34)可以等價(jià)為問題(33)的形式.□

利用定理1,可以將凸優(yōu)化子問題(28)與(29)轉(zhuǎn)換為線性規(guī)劃問題求解.

但是此種解法僅僅適用于圖像像素點(diǎn)較少的情況,原因在于A ∈R2N×p,轉(zhuǎn)換后新增的變量y∈R2N,z∈R2N,變量的個(gè)數(shù)為p+4N,約束條件的個(gè)數(shù)為1+6N.當(dāng)像素點(diǎn)個(gè)數(shù)達(dá)到N=10000的時(shí)候,所需解決的就是一個(gè)中規(guī)模的線性規(guī)劃問題.而這僅僅是算法中更新矩陣變量H與ggg第i行元素的子步驟,若要完成算法的一個(gè)完整的迭代過程,將帶來很大的時(shí)間開銷.因此需要能夠快速求解子問題(28)與(29)的新方法.

2010年后,ADMM 在求解大型凸優(yōu)化問題的領(lǐng)域展現(xiàn)了其強(qiáng)大的能力,也給求解此類凸優(yōu)化問題帶來了契機(jī).

以子問題(28)為例,給出采用ADMM算法解決此問題的步驟.

問題(28)有如下等價(jià)形式:

此問題的ADMM形式為:

根據(jù)ADMM算法的框架,采用表2所描述的算法流程,依次更新x,z1, z2,d1,d2,求解子問題(28)中的p?.

表2 RMVHU中求解子問題p?的算法步驟Table 2 Steps for solving subproblem p?in RMVHU

下面將介紹如何求解步驟2至步驟4的優(yōu)化子問題.

步驟2的子問題為:

式(39)的問題為簡單的無約束二次規(guī)劃問題,目標(biāo)函數(shù)梯度為0的極值點(diǎn)即為最優(yōu)解.令問題(39)目標(biāo)函數(shù)的梯度為0,得到步驟2的解析最優(yōu)解:

步驟3的子問題為:

該問題是對(duì)變量z1解耦的不等式約束最優(yōu)化問題,因此我們能夠直接給出問題(41)的最優(yōu)解形式.更方便的是,z1其實(shí)就是一個(gè)實(shí)數(shù).式(41)最優(yōu)解的解析形式表示為:

步驟4描述的問題為:

因此,問題的解可以采用著名的軟閾值[25]給出:

式中,soft??,ν￠即為著名的軟閾值函數(shù).soft??,ν￠表示對(duì)矩陣?中每個(gè)元素?i,j作如下映射:

在矩陣的運(yùn)算中,最為耗時(shí)的是求解矩陣的逆.我們可以發(fā)現(xiàn),在問題(28)的求解中,除了采用式(40)更新 xk時(shí)需要進(jìn)行矩陣的逆運(yùn)算,其余變量的更新只要進(jìn)行簡單的矩陣加減乘除與元素大小的比較操作,因此著重關(guān)注式(40)中矩陣Ψ的大小.根據(jù)式(30)與(31),A∈R2N×p, c∈Rp,因此Ψ∈Rp×p,p為圖像中的端元個(gè)數(shù).在高光譜圖像中,端元個(gè)數(shù)p幾乎都在幾十以內(nèi),因此,求解Ψ?1的運(yùn)算量很小,并且在子問題的求解中,Ψ?1為常值,意味著在迭代過程中只需要計(jì)算一次Ψ?1,這意味著從時(shí)間與空間開銷上考慮,求逆運(yùn)算對(duì)子問題的求解并不會(huì)產(chǎn)生很大的影響.求解問題(29)的方法與表2中的步驟類似,下面直接給出應(yīng)用表2中步驟求解時(shí),問題(29)對(duì)應(yīng)的步驟2至步驟4的最優(yōu)解表達(dá)形式:

關(guān)于表2步驟7中的迭代停止條件,當(dāng)滿足如下兩種情形中的一種時(shí),停止迭代.

情形1.迭代步數(shù)k到達(dá)最大迭代步數(shù).

情形2.原始?xì)埐钆c對(duì)偶?xì)埐?見下文第2.3.2節(jié)中定義)小于預(yù)設(shè)閾值.

特別要說明的是關(guān)于兩個(gè)優(yōu)化問題的優(yōu)化值取舍的問題.如果p?＜ q?,則采用子問題(28)的解,否則采用子問題(29)的解.

最后,假設(shè)RMVHU算法最終求解問題(27)得到的最終解為H?與ggg?,則可由式(48)得到端元矩陣A:

式中,C與ddd由式(14)得到.

根據(jù)式(21),可得到第n個(gè)像元的豐度

式中,p為端元個(gè)數(shù).

2.3 RMVHU算法參數(shù)的自適應(yīng)調(diào)節(jié)策略

2.3.1 正則項(xiàng)系數(shù)的自適應(yīng)調(diào)節(jié)策略

正則項(xiàng)系數(shù)的大小控制著正則項(xiàng)在整個(gè)優(yōu)化函數(shù)中的權(quán)值,也影響著整個(gè)算法的解混效果,當(dāng)正則項(xiàng)所占的比例較小時(shí),其起的作用相應(yīng)較小,反之亦然.在RMVHU算法中,若正則項(xiàng)系數(shù)λ不變,則在每一步的迭代中,由于 c Tx項(xiàng)每次都會(huì)變化,其所起的作用將會(huì)隨每一次的更新而改變,不利于整體優(yōu)化問題的收斂.因此有必要固定其在每一步子優(yōu)化問題中的占比.受文獻(xiàn)[26]啟發(fā),采用如下策略自適應(yīng)調(diào)整參數(shù)λ:

式中, x0為 x在迭代之前的初值.ω為常值比例系數(shù),實(shí)驗(yàn)中發(fā)現(xiàn)其值在30～50之間算法效果較好.

下面通過定性的分析,表明采用上述的自適應(yīng)參數(shù)調(diào)整法是行之有效的.

以上兩點(diǎn)原因,解釋了為何在正則項(xiàng)的系數(shù)與體積約束項(xiàng)大小成正比,與負(fù)數(shù)懲罰正則項(xiàng)大小成反比時(shí),RMVHU算法能夠取得很好的效果.

2.3.2 懲罰系數(shù)μ的自適應(yīng)調(diào)整

根據(jù)文獻(xiàn)[22],子問題ADMM 形式(38)中的懲罰項(xiàng)系數(shù)μ在如下簡單自適應(yīng)調(diào)節(jié)的策略下,使算法的收斂性得到提高:

式中,rk為原始?xì)埐?sk為對(duì)偶?xì)埐?根據(jù)原始?xì)埐钆c對(duì)偶?xì)埐钤谖墨I(xiàn)[22]中的定義,在RMVHU算法中,原始?xì)埐钆c對(duì)偶?xì)埐钅軌蛴梢韵鹿竭M(jìn)行計(jì)算:

算法的主要流程如下表所示:

表3RMVHU算法步驟Table 3 Procedure for solving RMVHU

在表3描述的算法流程中,在更新矩陣變量H與g 時(shí),一次完整的迭代步驟為:采用步驟3至步驟5將矩陣變量的每一行元素都更新一遍.也就是說,在算法的整個(gè)求解過程中,分別最多求解M×(p?1)個(gè)子問題(28)與(29).算法在采用表2描述的ADMM算法求解子問題后是相對(duì)高效的,因此,RMVHU算法必然能夠在有限時(shí)間內(nèi)給出可行解.

2.4 RMVHU算法時(shí)間復(fù)雜度與收斂性分析

為進(jìn)一步研究RMVHU算法的效率,本節(jié)將詳細(xì)分析算法的復(fù)雜度.由于算法的本質(zhì)為求解M×(p?1)個(gè)子問題(28)與(29),因此,將以求解子問題(28)的時(shí)間復(fù)雜度為分析的著眼點(diǎn).

當(dāng)采用表2所描述的算法求解子問題(28)時(shí),更新xxx需做2Np+4N+p次加法或者乘法運(yùn)算,更新z1需做p+2次加法或者乘法運(yùn)算,更新 zz2需做2Np+12N 次加法或者乘法運(yùn)算,更新d1需做p+2次加法或者乘法運(yùn)算,更新dd2需做2Np+6N 次加法或者乘法運(yùn)算,在計(jì)算時(shí),需做次加法或者乘法運(yùn)算,在計(jì)算時(shí),需做2Np+4N+p+8次加法或者乘法運(yùn)算.因此,在表2描述的算法中,循環(huán)一次需要做10Np+32N+5p+8次加法或者乘法運(yùn)算,算法的時(shí)間復(fù)雜度為相應(yīng)的為O(Np+N+p),可以發(fā)現(xiàn),時(shí)間復(fù)雜度與圖像像素點(diǎn)個(gè)數(shù)N 成正比,與端元個(gè)數(shù)p成正比,當(dāng)p?N 時(shí),算法的時(shí)間復(fù)雜度能夠被近似表示為O(N).

當(dāng)采用定理1的方法,將問題(28)轉(zhuǎn)化為線性規(guī)劃問題求解時(shí),變量個(gè)數(shù)為4N+p.采用Karmarkar算法求解時(shí),時(shí)間復(fù)雜度[27]為O(N3.5).

顯然,采用ADMM 算法求解在時(shí)間效率上具有顯然的優(yōu)勢,尤其當(dāng)N 很大的時(shí)候.

為了分析RMVHU算法的收斂性,本文將從如下兩個(gè)方面進(jìn)行闡述.

1)在式(27)描述的優(yōu)化目標(biāo)函數(shù)中,當(dāng)固定正則項(xiàng)系數(shù)λ不變,采用循環(huán)最小化的方法,逐行更新H與g時(shí),所求解的兩個(gè)子問題都是凸問題,并且ADMM算法在求解凸問題時(shí)具有收斂性,因此,式(27)的目標(biāo)函數(shù)值必將逐漸減少,說明在采用循環(huán)最小化的方法,執(zhí)行步驟4至步驟5時(shí),算法是收斂的.

2)當(dāng)執(zhí)行完步驟4與步驟5,算法仍未收斂時(shí),將根據(jù)式(50)的正則系數(shù)自適應(yīng)調(diào)整策略,重新計(jì)算λ,這時(shí),根據(jù)第2.3.1節(jié)中的定性分析,算法將在此種策略下加速收斂.

3 仿真校驗(yàn)

本部分由兩塊內(nèi)容組成,首先介紹實(shí)驗(yàn)數(shù)據(jù)的來源以及算法評(píng)價(jià)準(zhǔn)則,隨后給出實(shí)驗(yàn)結(jié)果并進(jìn)行分析.

3.1 算法實(shí)驗(yàn)數(shù)據(jù)來源與評(píng)價(jià)準(zhǔn)則

算法的實(shí)驗(yàn)數(shù)據(jù)分為合成光譜數(shù)據(jù)與真實(shí)光譜數(shù)據(jù).

3.1.1 合成高光譜數(shù)據(jù)

在USGS光譜庫[28]中選擇p個(gè)光譜向量作為端元.圖2顯示了庫中的3條光譜曲線.

圖2 USGS庫中不同物質(zhì)的光譜曲線Fig.2 USGS library spectra of diあerent materials

根據(jù)文獻(xiàn)[16]的方法,產(chǎn)生像元的豐度信息.像元的豐度信息滿足狄利克雷分布,狄利克雷分布表示為:

式中,第i個(gè)端元的豐度 αi滿足0≤ αi≤ 1,且豐度期望Γ(·)表示Gamma函數(shù).狄利克雷的密度分布不僅能夠保證豐度符合非負(fù)與和為1約束,還能夠根據(jù)不同的參數(shù)μi,產(chǎn)生不同分布特征的豐度系數(shù).為了產(chǎn)生高度混合的數(shù)據(jù),對(duì)于豐度αi＞ρ的像元,將其豐度重新匹配,改為1/p.ρ為表示像元混合度的參數(shù)且滿足ρ∈(0,1),ρ越小,圖像混合度越高.

為了使數(shù)據(jù)更具有真實(shí)性,加入獨(dú)立同分布的零均值高斯加性噪聲.信噪比SNR表示為:

式中，y[n]為真實(shí)的光譜數(shù)據(jù),σ2代表高斯噪聲的方差,L為波段維數(shù),N為像素點(diǎn)個(gè)數(shù).

在仿真數(shù)據(jù)中,還應(yīng)考慮異常點(diǎn)對(duì)算法的影響,因此需要生成異常點(diǎn).假設(shè)異常點(diǎn)的個(gè)數(shù)為k,則重新生成第1至k個(gè)像素點(diǎn)的豐度.第i(0≤i≤k)個(gè)點(diǎn)的豐度由以下公式重新產(chǎn)生:

式中,δ為常數(shù),k=randperm(p),代表在1與p之間滿足均勻分布的隨機(jī)整數(shù),ss(i)為第i個(gè)像元豐度的向量表示,sm(i)代表ss(i)的第m個(gè)元素.

3.1.2 實(shí)際數(shù)據(jù)的來源

本文使用的實(shí)際數(shù)據(jù)是美國內(nèi)華達(dá)州Cuprite地區(qū)1997年成像的224波段0.4μm 至2.5μm,大小為250×191的子圖[29].這組數(shù)據(jù)在高光譜圖像端元提取的研究中廣泛應(yīng)用,具有很好的代表性.由于水汽吸收的干擾和低信噪比的原因,第1～2,104～113,148～167以及221～224波段的數(shù)據(jù)被剔除.偽彩色子圖如圖3所示.

該地區(qū)的礦物分布圖由Tricorder 3.3軟件[30]提供,如圖4所示.需要注意的是此物質(zhì)分布在1995年就已問世,而Cuprite地區(qū)的數(shù)據(jù)是在1997年采集.因此該物質(zhì)分布圖僅作為解混效果好壞的參考.

圖3 Cuprite地區(qū)高光譜圖像偽彩色子圖,R:2.109μm,G:2.209μm,B:2.308μmFig.3 The Pseudo color subimage of the AVIRIS Cuprite dataset,R:2.109μm,G:2.209μm,B:2.308μm

圖4 Cuprite地區(qū)礦物分布圖Fig.4 The mineral distribution map of Cuprite

3.1.3 算法評(píng)價(jià)準(zhǔn)則

為了描述解混算法的效果好壞,通常需要研究解混后端元與真實(shí)端元的相似度以及解混后估計(jì)豐度與真實(shí)豐度的相似程度.可以定義光譜角距離(SAD)和誤差平方根(RMSE)來評(píng)估實(shí)驗(yàn)所得的端元和豐度.

對(duì)于第i個(gè)端元,假設(shè)實(shí)際的端元的向量表示形式為 ai,采用解混算法得到的端元估計(jì)值的向量表示為第i個(gè)端元的光譜角距離表示為φi:

對(duì)于第i個(gè)端元,假設(shè)端元對(duì)應(yīng)的真實(shí)豐度表示為si,采用解混算法得到與端元對(duì)應(yīng)的豐度估計(jì)值表示為N 為圖像像素點(diǎn)個(gè)數(shù),第i個(gè)端元的均方根誤差表示為θi:

為全局地評(píng)價(jià)解混算法的好壞,將采用如下兩個(gè)平均指標(biāo):平均光譜角距離εφ與平均誤差平方根εθ進(jìn)行算法效果的衡量:

式(59)與(60)中,p為圖像中端元個(gè)數(shù).

3.2 實(shí)驗(yàn)及結(jié)果分析

實(shí)驗(yàn)分為合成數(shù)據(jù)解混與真實(shí)數(shù)據(jù)解混兩組實(shí)驗(yàn).在采用合成數(shù)據(jù)的實(shí)驗(yàn)時(shí),首先研究τ與γ的取值對(duì)ADMM 算法收斂的影響,研究第2.3.1節(jié)中正則系數(shù)自適應(yīng)調(diào)整的策略的有效性,隨后在端元數(shù)p=3且存在異常點(diǎn)的數(shù)據(jù)集中,將RMVHU算法與MVES算法進(jìn)行比較,闡述RMVHU算法的優(yōu)點(diǎn);其次將RMVHU算法與被廣泛使用的VCA、MVES、MVSA、SISAL四種凸面幾何學(xué)算法進(jìn)行比較,分別研究真實(shí)端元個(gè)數(shù),噪聲大小,數(shù)據(jù)混合度ρ,異常點(diǎn)個(gè)數(shù),像元個(gè)數(shù)以及錯(cuò)估端元個(gè)數(shù)時(shí)對(duì)不同算法解混效果的影響.最后,對(duì)于真實(shí)數(shù)據(jù),給出RMVHU算法的解混結(jié)果.

3.2.1 合成數(shù)據(jù)的仿真

實(shí)驗(yàn)1.RMVHU懲罰系數(shù)的調(diào)整實(shí)驗(yàn)

本實(shí)驗(yàn)將進(jìn)一步研究第2.3.2節(jié)中懲罰系數(shù)γ與τ的取值對(duì)算法收斂性的影響.實(shí)驗(yàn)中,γ取值為1,10,100,200,τ取值為0.5,1.5,2,4.

圖5與圖6為求解子問題(28)時(shí),目標(biāo)函數(shù)值在不同γ與τ下的變化趨勢.從圖5可以發(fā)現(xiàn),當(dāng)γ較小時(shí),函數(shù)收斂較快,但波動(dòng)較大,當(dāng)γ增大時(shí),目標(biāo)函數(shù)雖然波動(dòng)較小,但收斂速度變慢.產(chǎn)生上述結(jié)果的原因在于,γ較小時(shí),只要原始?xì)埐钆c對(duì)偶?xì)埐畲嬖诓顒e,那么懲罰系數(shù)μ便會(huì)頻繁的變化,雖然這樣的變化能使收斂速度提高,但會(huì)帶來收斂穩(wěn)定性的不足;而當(dāng)γ較大時(shí),原始?xì)埐钆c對(duì)偶?xì)埐畹牟町愋詫⒑茈y影響到懲罰系數(shù)μ的改變,而合適的μ將對(duì)收斂速度產(chǎn)生較大的影響,因此,雖然ADMM算法具有收斂性,但是收斂速度將放緩.在圖5中可以發(fā)現(xiàn),γ=10時(shí),函數(shù)收斂較快,波動(dòng)較小.

圖5 γ變化時(shí)目標(biāo)函數(shù)值Fig.5 The values of object function when γ changes

圖6 τ變化時(shí)目標(biāo)函數(shù)值Fig.6 The values of object function when τ changes

從圖6可以發(fā)現(xiàn),當(dāng)τ＜1時(shí),函數(shù)不收斂,隨著τ的增大時(shí),函數(shù)收斂速度加快,將很快收斂于極值,但是波動(dòng)將會(huì)增加,其原因在于,當(dāng)對(duì)偶?xì)埐钆c原始?xì)埐钶^大時(shí),若τ＜1,根據(jù)調(diào)整策略,懲罰系數(shù)μ將增大,這將導(dǎo)致對(duì)偶?xì)埐罡?勢必造成算法的不收斂;當(dāng)τ＞1,則根據(jù)調(diào)整策略,懲罰系數(shù)μ將減小,這將縮小對(duì)偶?xì)埐钆c原始?xì)埐畹拇笮?使得算法收斂性得到提升,但是過大的τ將使對(duì)偶?xì)埐畈▌?dòng)較大,帶來收斂穩(wěn)定性的不足.在圖6中可以發(fā)現(xiàn)τ=2時(shí),函數(shù)收斂較快,波動(dòng)較小.

因此在接下來的實(shí)驗(yàn)中,取γ=10,τ=2.

實(shí)驗(yàn)2.正則項(xiàng)系數(shù)的自適應(yīng)調(diào)節(jié)對(duì)算法的影響

為了進(jìn)一步研究第2.3.1節(jié)的正則項(xiàng)系數(shù)的自適應(yīng)調(diào)節(jié)策略對(duì)算法的影響,將進(jìn)行相關(guān)實(shí)驗(yàn)進(jìn)行對(duì)比驗(yàn)證.實(shí)驗(yàn)中將比較隨著端元個(gè)數(shù)的變化時(shí),兩種方法的解混精度.固定系數(shù)的算法取參數(shù)λ=0.002,采用第2.3.1節(jié)中描述的策略的實(shí)驗(yàn),取參數(shù)ω=40.

實(shí)驗(yàn)中,模擬數(shù)據(jù)的像素點(diǎn)個(gè)數(shù)N=10000,波段維數(shù)L=220,端元數(shù)p從3至8變化,SNR=30dB,數(shù)據(jù)混合度ρ=0.8,異常點(diǎn)個(gè)數(shù)為25,估計(jì)端元個(gè)數(shù)無錯(cuò)估.實(shí)驗(yàn)結(jié)果如圖7與圖8所示.

圖7 正則項(xiàng)系數(shù)固定或自適應(yīng)變化時(shí),各算法平均光譜角距離Fig.7 The average spectral angle distance of diあerent algorithms when the regularization parameter is fi xed or changes automatically

圖8 正則項(xiàng)系數(shù)固定或自適應(yīng)變化時(shí),各算法平均均方根誤差Fig.8 The average root mean square error of diあerent algorithms when the regularization parameter is fi xed or changes automatically

圖7與圖8表明,當(dāng)正則項(xiàng)系數(shù)固定時(shí),實(shí)驗(yàn)效果將會(huì)隨著端元個(gè)數(shù)的改變發(fā)生很大的變化.而采用第2.3.1節(jié)策略的算法則能夠穩(wěn)定地估計(jì)真實(shí)端元與豐度.由此可以說明,采用第2.3.1節(jié)策略將使算法具有更強(qiáng)的魯棒性.

實(shí)驗(yàn)3.RMVHU的魯棒性驗(yàn)證實(shí)驗(yàn)

本實(shí)驗(yàn)主要通過與一種流行的凸面幾何學(xué)解混算法MVES進(jìn)行比較,直觀地表明RMVHU算法在存在異常點(diǎn)情況下的魯棒性.

模擬數(shù)據(jù)的像素點(diǎn)個(gè)數(shù)N=10000,波段維數(shù)L=220,端元個(gè)數(shù)p=3(選取的端元光譜曲線見圖2),SNR=30dB,數(shù)據(jù)混合度ρ=0.8,異常點(diǎn)個(gè)數(shù)為25(像素點(diǎn)總數(shù)的0.25%),估計(jì)端元個(gè)數(shù)為3(端元個(gè)數(shù)無錯(cuò)估).

RMVHU與MVES端元估計(jì)結(jié)果如圖9所示.在圖9中,實(shí)線代表真實(shí)的端元光譜曲線,虛線代表解混算法提取出的端元曲線.可以發(fā)現(xiàn),在少量異常點(diǎn)的干擾下,RMVHU算法的端元結(jié)果遠(yuǎn)遠(yuǎn)優(yōu)于MVES算法.其主要原因在于,RMVHU算法引入了正則項(xiàng),該正則項(xiàng)能夠容忍豐度值為負(fù)的情形,使得噪聲與異常點(diǎn)對(duì)最小體積逼近的影響不那么明顯,而MVES算法需要嚴(yán)格滿足非負(fù)的條件,因此MVES計(jì)算的單形體頂點(diǎn)需嚴(yán)格包含所有觀測數(shù)據(jù),使異常點(diǎn)對(duì)最小體積逼近的影響尤為突出,導(dǎo)致了端元提取效果的劇烈下降.

圖10顯示了數(shù)據(jù)在二維空間的分布以及RMVHU、MVES算法提取的端元在二維空間的分布,此圖亦印證了上述緣由.

表4與表5分別記錄了RMVHU與MVES算法的光譜角距離,平均光譜角距離與均方根誤差,平均均方根誤差值.結(jié)果最優(yōu)值在表中加粗.從表4與表5中可以發(fā)現(xiàn),若圖像數(shù)據(jù)中存在異常值,則在進(jìn)行圖像的解混時(shí),RMVHU算法無論是在端元提取精度上,還是豐度反演精度上,都遠(yuǎn)遠(yuǎn)優(yōu)于MVES算法.隨后進(jìn)行的合成數(shù)據(jù)仿真實(shí)驗(yàn)中,將分別研究真實(shí)端元個(gè)數(shù),噪聲大小,數(shù)據(jù)混合度,有無異常點(diǎn),像素點(diǎn)個(gè)數(shù)以及錯(cuò)估端元個(gè)數(shù)對(duì)不同解混算法的影響.

表4 RMVHU與MVES算法的光譜角距離與平均光譜角距離Table 4 Spectral angle distance and average spectral angle distance via RMVHU and MVES

表5 RMVHU與MVES算法的均方根誤差與平均均方根誤差Table 5 Root mean square error and average root mean square error via RMVHU and MVES

圖9 RMVHU與MVES端元估計(jì)結(jié)果,(a)、(b)、(c)為RMVHU的端元估計(jì)結(jié)果,(d)、(e)、(f)為MVES的端元估計(jì)結(jié)果Fig.9 Estimated results of endmembers via two algorithms:RMVHU and MVES((a)、(b)、(c)are the results of RMVHU and(d)、(e)、(f)are the results of MVES)

圖10 數(shù)據(jù)集與端元在p?1(二)維子空間的投影Fig.10 Two dimensional subspace projection of datasets,estimated endmembers and real endmembers

實(shí)驗(yàn)4.真實(shí)端元個(gè)數(shù)變化的實(shí)驗(yàn)

在每一幅高光譜圖像中,端元個(gè)數(shù)會(huì)根據(jù)不同的拍攝場景而變化,因此,研究端元個(gè)數(shù)對(duì)解混結(jié)果的影響很有必要.

本實(shí)驗(yàn)中,生成端元個(gè)數(shù)p從3～10變化的8組數(shù)據(jù),其他參數(shù)設(shè)置為:模擬數(shù)據(jù)的像素點(diǎn)個(gè)數(shù)N=10000,波段維數(shù)L=220,SNR=30dB,數(shù)據(jù)混合度ρ=0.8,異常點(diǎn)個(gè)數(shù)為25(像素點(diǎn)總數(shù)的0.25%),端元個(gè)數(shù)無錯(cuò)估.

平均光譜角距離結(jié)果如圖11所示.平均均方根誤差結(jié)果如圖12所示.

圖11與圖12表明,相較于其他算法,端元個(gè)數(shù)的變化對(duì)RMVHU算法的影響并不很明顯.RMVHU算法的解混效果相對(duì)于基于純像元假設(shè)的算法—VCA得到的結(jié)果優(yōu)秀很多.相比于需要嚴(yán)格滿足非負(fù)性條件的最小體積變化算法:MVES與MVSA,RMVHU提取的端元的精度要優(yōu)秀很多.對(duì)于同樣可以不滿足非負(fù)條件,但是采用二次逼近體積項(xiàng)的SISAL算法來說,由于RMVHU算法子問題中的體積項(xiàng)采用等價(jià)的行列式展開公式,因此其求解的精度也在大多數(shù)情況下比SISAL算法高.

圖11 端元個(gè)數(shù)變化時(shí),各算法平均光譜角距離Fig.11 The average spectral angle distance of diあerent algorithms when the number of endmembers changes

圖12 端元個(gè)數(shù)變化時(shí),各算法平均均方根誤差Fig.12 The average root mean square error of diあerent algorithms when the number of endmembers changes

實(shí)驗(yàn)5.抗噪能力的測試實(shí)驗(yàn)

設(shè)計(jì)本實(shí)驗(yàn)的目的在于測試在不同大小的噪聲下,RMVHU算法的穩(wěn)定性.

本實(shí)驗(yàn)中,信噪比SNR 分別設(shè)置為15dB、20dB、25dB、30dB、35dB、40dB,獲得6組數(shù)據(jù),其他參數(shù)設(shè)置為:模擬數(shù)據(jù)的像素點(diǎn)個(gè)數(shù)N=10000,波段維數(shù)L=220,端元個(gè)數(shù)p=6,數(shù)據(jù)混合度ρ=0.8,異常點(diǎn)個(gè)數(shù)為25(像素點(diǎn)總數(shù)的0.25%),端元個(gè)數(shù)無錯(cuò)估.

在不同噪聲大小下,平均光譜角距離與平均均方根誤差結(jié)果如圖13與圖14所示.

從圖13與圖14中可以發(fā)現(xiàn),信噪比越高,RMVHU算法的解混效果越好.由于VCA中采用了根據(jù)不同噪聲大小進(jìn)行不同的數(shù)據(jù)降維方法,因此VCA在低信噪比下有較突出的表現(xiàn).而本文提出的RMVHU算法效果在低信噪比下的端元提取精度僅次于VCA,側(cè)面印證了RMVHU算法正則項(xiàng)及自適應(yīng)的正則算子調(diào)節(jié)方法有較強(qiáng)的抗噪能力.相較于其他幾何學(xué)解混算法,在信噪比大于20dB的情況下,RMVHU算法的表現(xiàn)出色,能夠得到最低的平均光譜角距離與平均均方根誤差值.RMVHU算法較強(qiáng)的抗噪聲能力得到了實(shí)驗(yàn)結(jié)果的支持.

圖13 在不同噪聲大小下,各算法平均光譜角距離Fig.13 The average spectral angle distance of diあerent algorithms when SNR changes

圖14 在不同噪聲大小下,各算法平均均方根誤差Fig.14 The average root mean square error of diあerent algorithms when SNR changes

實(shí)驗(yàn)6.數(shù)據(jù)混合度變化的實(shí)驗(yàn)

設(shè)計(jì)本實(shí)驗(yàn)的目的在于:驗(yàn)證在不同混合度的場景下,RMVHU算法也有較為突出的表現(xiàn),即能夠勝任不同混合程度的高光譜圖像的解混.

在本實(shí)驗(yàn)中,數(shù)據(jù)混合度ρ分別設(shè)置為0.65、0.7、0.75、0.8、0.85、0.9、0.95、1,獲得 8組數(shù)據(jù).其他參數(shù)設(shè)置如下:模擬數(shù)據(jù)的像素點(diǎn)個(gè)數(shù)N=10000,波段維數(shù)L=220,真實(shí)端元個(gè)數(shù)p為6,SNR為30dB,異常點(diǎn)個(gè)數(shù)為25(像素點(diǎn)總數(shù)的0.25%),端元個(gè)數(shù)無錯(cuò)估.

在不同混合度下,平均光譜角距離與平均均方根誤差結(jié)果如圖15與圖16所示.

由圖15與圖16發(fā)現(xiàn),RMVHU在不同的混合度下,都能取得最為優(yōu)異的結(jié)果.同時(shí),由于異常值的存在,使得本應(yīng)該在混合度較低,即ρ接近1的情況下能夠獲得更好效果的純像元提取類算法VCA變得精度很差.但是本文提出的解混算法RMVHU在異常點(diǎn)存在的情況下,對(duì)數(shù)據(jù)的混合度大小并不敏感.同時(shí)由于實(shí)驗(yàn)4中所述的原因,RMVHU相比于SISAL也更加優(yōu)秀.本實(shí)驗(yàn)的成功意味著該算法在對(duì)不同混合度的高光譜圖像解混時(shí),能夠取得很好的效果.

圖15 在不同混合度下,各算法平均光譜角距離Fig.15 The average spectral angle distance of diあerent algorithms when the purity of pixels changes

圖16 在不同混合度下,各算法平均均方根誤差Fig.16 The average root mean square error of diあerent algorithms when the purity of pixels changes

實(shí)驗(yàn)7.高于異常點(diǎn)個(gè)數(shù)變化的實(shí)驗(yàn)

設(shè)計(jì)本實(shí)驗(yàn)的目的在于:驗(yàn)證在異常點(diǎn)個(gè)數(shù)發(fā)生變化時(shí),RMVHU算法有較強(qiáng)的穩(wěn)定性,并且與其他算法一起,比較異常點(diǎn)的有無及數(shù)量對(duì)算法解混精度的影響.

實(shí)驗(yàn)中, 異常點(diǎn)個(gè)數(shù)分別設(shè)置為像素點(diǎn)總數(shù)的 0% (無異常數(shù)據(jù)),0.25%、0.5%、0.75%、1%、1.25%、1.5%、1.75%、2%,獲得9組數(shù)據(jù).其他參數(shù)設(shè)置如下:模擬數(shù)據(jù)的像素點(diǎn)個(gè)數(shù)N=10000,波段維數(shù)L=220,真實(shí)端元個(gè)數(shù)p為6,SNR為30dB,數(shù)據(jù)混合度ρ設(shè)置為0.8,端元個(gè)數(shù)無錯(cuò)估.

在異常點(diǎn)個(gè)數(shù)改變的情形下,平均光譜角距離與平均均方根誤差結(jié)果如圖17與圖18所示.

圖17 在不同異常點(diǎn)個(gè)數(shù)下,各算法平均光譜角距離Fig.17 The average spectral angle distance of diあerent algorithms when the number of outliers changes

圖18 在不同異常點(diǎn)個(gè)數(shù)下,各算法平均均方根誤差Fig.18 The average root mean square error of diあerent algorithms when the number of outliers changes

由圖17與圖18可以發(fā)現(xiàn),圖像中無異常數(shù)據(jù)時(shí),基于最小體積變化的算法RMVHU、MVES、MVSA及SISAL算法都能夠取得較好端元提取與解混的效果;但是當(dāng)出現(xiàn)異常數(shù)據(jù)時(shí),除了RMVHU算法能夠取得最為優(yōu)異的解混結(jié)果,SISAL算法勉強(qiáng)能夠較為準(zhǔn)確地獲得圖像端元與豐度信息外,其他算法獲得的結(jié)果都與真實(shí)值有較大偏差.同時(shí),相比于SISAL算法的結(jié)果隨著異常點(diǎn)個(gè)數(shù)增加而變差,RMVHU算法的解混誤差基本不變,維持在最低的水平.本實(shí)驗(yàn)進(jìn)一步說明了RMVHU算法的穩(wěn)定性.

實(shí)驗(yàn)8.關(guān)于像素點(diǎn)個(gè)數(shù)的變化的實(shí)驗(yàn)

設(shè)計(jì)本實(shí)驗(yàn)的目的在于研究圖像像素點(diǎn)個(gè)數(shù)對(duì)RMVHU算法解混精度的影響以及算法的時(shí)間復(fù)雜度(即能否在較短的時(shí)間內(nèi)得到較為精確的解混結(jié)果).

在本實(shí)驗(yàn)中, 設(shè)置像素點(diǎn)個(gè)數(shù)為2000、4000、6000、8000、10000、12000、14000、16000、18000、20000,共得到10組數(shù)據(jù).其他參數(shù)設(shè)置如下:真實(shí)端元個(gè)數(shù)p為6,波段維數(shù)L=220,SNR為30dB,數(shù)據(jù)混合度ρ設(shè)置為0.8,異常點(diǎn)個(gè)數(shù)為像素點(diǎn)總數(shù)的0.25%,端元個(gè)數(shù)無錯(cuò)估.計(jì)算機(jī)為聯(lián)想S20工作站.

在不同大小的高光譜圖像中,平均光譜角距離與平均均方根誤差結(jié)果如圖19與圖20所示,運(yùn)行時(shí)間如圖21所示.

圖19 在不同像素點(diǎn)個(gè)數(shù)下,各算法平均光譜角距離Fig.19 The average spectral angle distance of diあerent algorithms when the number of pixels changes

由圖19與圖20可以發(fā)現(xiàn),相較于其他算法,RMVHU算法解混偏差并不會(huì)像SISAL、MVES、MVSA、VCA 等算法隨著像素點(diǎn)的個(gè)數(shù)變化出現(xiàn)巨大的波動(dòng),它的表現(xiàn)一直很穩(wěn)定.當(dāng)像素點(diǎn)個(gè)數(shù)超過4000時(shí),RMVHU算法的解混結(jié)果是所有算法中最優(yōu)秀的.此實(shí)驗(yàn)驗(yàn)證了本文提出的RMVHU算法能夠在各種大小的高光譜圖像解混中,取得穩(wěn)定且優(yōu)秀的解混效果的說法.

在圖21中,可以發(fā)現(xiàn)隨著像素點(diǎn)個(gè)數(shù)的增多,RMVHU算法運(yùn)行的時(shí)間線性增長,這也與第2.4節(jié)中算法的復(fù)雜度分析相吻合.在像素點(diǎn)個(gè)數(shù)較大時(shí),由于不用求解大規(guī)模的線性規(guī)劃問題,其時(shí)間消耗反而要比MVES算法少.所以,在有限時(shí)間內(nèi),應(yīng)用RMVHU算法,定能得到理想的解混效果.相較于SISAL與MVSA等算法,雖然其時(shí)間消耗較多,但若更注重解混的精度與結(jié)果的穩(wěn)定性,那么在時(shí)間上的一些犧牲是值得的.

圖20 在不同像素點(diǎn)個(gè)數(shù)下,各算法平均均方根誤差Fig.20 The average root mean square error of diあerent algorithms when the number of pixels changes

圖21 在不同像素點(diǎn)個(gè)數(shù)下,各算法運(yùn)行時(shí)間Fig.21 Computational time of diあerent algorithms when the number of pixels changes

實(shí)驗(yàn)6.關(guān)于錯(cuò)誤估計(jì)端元個(gè)數(shù)的實(shí)驗(yàn)

在解混算法的實(shí)際應(yīng)用中,大多數(shù)解混算法需要預(yù)先估計(jì)端元的個(gè)數(shù),RMVHU算法也不例外.但是可能由于混合模型的非線性偏差,端元個(gè)數(shù)的估計(jì)有可能發(fā)生錯(cuò)誤,因此,設(shè)計(jì)實(shí)驗(yàn)研究算法在錯(cuò)誤估計(jì)端元個(gè)數(shù)時(shí)的表現(xiàn)也很有必要.

在本實(shí)驗(yàn)中,真實(shí)端元的個(gè)數(shù)p在3～8間變化.估計(jì)的端元個(gè)數(shù)為p?1,模擬端元個(gè)數(shù)錯(cuò)誤估計(jì)的情況.其他實(shí)驗(yàn)參數(shù)設(shè)置如下:模擬數(shù)據(jù)的像素點(diǎn)個(gè)數(shù)N=10000,波段維數(shù)L=220,SNR為30dB,數(shù)據(jù)混合度ρ設(shè)置為0.8,異常點(diǎn)個(gè)數(shù)設(shè)置為0.

在對(duì)被錯(cuò)估端元個(gè)數(shù)的高光譜圖像解混后,平均光譜角距離與平均均方根誤差結(jié)果如圖22與圖23所示.

圖22 在錯(cuò)估端元個(gè)數(shù)時(shí),各算法平均光譜角距離Fig.22 The average spectral angle distance of diあerent algorithms when the number of endmembers is incorrectly estimated

由圖22發(fā)現(xiàn)在絕大多數(shù)端元個(gè)數(shù)錯(cuò)估的情況下,RMVHU算法相比其他的最小體積變換算法:MVSA、MVES及SISAL,能夠取得更為優(yōu)秀的端元估計(jì)結(jié)果.

在對(duì)豐度值的估計(jì)方面,由圖23可以發(fā)現(xiàn),RMVHU算法在不同的端元個(gè)數(shù)下都獲得了不錯(cuò)的效果,且隨著真實(shí)端元個(gè)數(shù)的增加,精度也有提高的趨勢.對(duì)RMVHU算法而言,端元提取誤差與豐度反演的誤差隨著真實(shí)端元個(gè)數(shù)的增加,也有降低的趨勢,這也說明了RMVHU算法在錯(cuò)誤估計(jì)端元個(gè)數(shù)的情況下仍然具有較強(qiáng)的魯棒性.

3.2.2 真實(shí)數(shù)據(jù)的實(shí)驗(yàn)

采用Cuprite子圖,剔除干擾波段后,文獻(xiàn)[8]的研究結(jié)果表明,這片區(qū)域中有14種礦物.由于有些同種類不同化學(xué)成分的礦物的光譜只存在很微小的差別,因此在解混時(shí),將端元的個(gè)數(shù)減少至11.解混后各物質(zhì)豐度圖如圖24所示.

圖23 在錯(cuò)估端元個(gè)數(shù)時(shí),各算法平均均方根誤差Fig.23 The average root mean square error of diあerent algorithms when the number of endmembers is incorrectly estimated

表6 RMVHU、VCA、MVES、MVSA、SISAL算法提取端元與真實(shí)端元的光譜角距離Table 6 The spectral angle distance of real datasets via unmixing algorithms:RMVHU,VCA,MVES,MVSA and SISAL

圖24 Cuprite地區(qū)估計(jì)豐度圖Fig.24 Estimated abundance maps of AVIRIS Cuprite

在圖24的Cuprite地區(qū)估計(jì)豐度圖中,(a)為Desert Varnish的豐度圖,(b)為Sphene的豐度圖,(c)為Nontronite的豐度圖,(d)為Kaolinite的豐度圖,(e)為Dumortierite的豐度圖,(f)為Chalcedony的豐度圖,(g)為Pyrope的豐度圖,(h)為Andradite的豐度圖,(i)為Montmorillonite的豐度圖,(j)為Muscovite的豐度圖,(k)為Alunite的豐度圖.

同時(shí)為了比較端元提取效果,將RMVHU、VCA、MVES、MVSA、SISAL算法的端元提取結(jié)果與光譜庫中的標(biāo)準(zhǔn)數(shù)據(jù)對(duì)比,得到各物質(zhì)的光譜角距離,如表6所示.針對(duì)所有物質(zhì),各算法得到的光譜角距離的最小值已在表中加粗.

將圖24與圖4相比,可以發(fā)現(xiàn),RMVHU算法確實(shí)能夠?qū)⒏吖庾V圖像中的物質(zhì)信息有效地提取出來.

從表6中可以發(fā)現(xiàn),RMVHU算法在對(duì)如下6種物質(zhì):Desert Varnish、Dumortierite、Chalcedony、Pyrope、Andrad ite以及Muscovite的端元提取中,表現(xiàn)最為優(yōu)秀.對(duì)于其他的物質(zhì),除了Sphene,另外四種物質(zhì)的端元提取效果也較為優(yōu)秀,RMVHU算法的端元提取結(jié)果也非常接近由其他算法所獲得的最優(yōu)秀的結(jié)果.最后可以發(fā)現(xiàn),RMVHU算法的平均光譜角誤差是所有算法中最小的,證明了本算法在實(shí)際解混應(yīng)用中性能的優(yōu)越性.

4 結(jié)論

本文通過引入負(fù)數(shù)懲罰正則項(xiàng),提出了自適應(yīng)魯棒最小體積高光譜解混算法(RMVHU),在求解時(shí)采用循環(huán)最小化方法拆分成可以求解的凸優(yōu)化子問題,并成功應(yīng)用ADMM框架對(duì)子問題進(jìn)行求解.此外,文章對(duì)正則項(xiàng)的系數(shù)λ提出了一種自適應(yīng)的調(diào)節(jié)策略,加快了收斂速度,提高了算法的穩(wěn)定性.

仿真實(shí)驗(yàn)表明相較于其他解混算法,RMVHU算法在任意端元個(gè)數(shù)下均擁有優(yōu)秀且穩(wěn)定的解混表現(xiàn);此外,算法對(duì)數(shù)據(jù)噪聲大小,數(shù)據(jù)異常點(diǎn)個(gè)數(shù)以及圖像像素點(diǎn)個(gè)數(shù)亦很不敏感,尤其在端元數(shù)估計(jì)錯(cuò)誤的情況下也有比較穩(wěn)定的表現(xiàn);在不同混合度的場景下RMVHU也有很優(yōu)秀的表現(xiàn).

對(duì)Cuprite真實(shí)數(shù)據(jù)的解混結(jié)果表明,RMVHU算法提取出的端元及估算出的豐度圖與實(shí)際地物組成基本一致.

需要注意的是,RMVHU算法的計(jì)算時(shí)間相對(duì)于諸如VCA等純像元提取算法較長,因此在下一步的研究中將嘗試采用并行化技術(shù),以提高運(yùn)算的效率.

致謝

作者感謝審稿專家在論文修改過程中提出的寶貴建議和意見.作者對(duì)哈爾濱工業(yè)大學(xué)航天學(xué)院檢測實(shí)驗(yàn)室的研究支持表示感謝.

1 Pan Zong-Xu,Yu Jing,Xiao Chuang-Bai,Sun Wei-Dong.Spectral similarity-based super resolution for hyperspectral images.Acta Automatica Sinica,2014,40(12):2797?2807(潘宗序,禹晶,肖創(chuàng)柏,孫衛(wèi)東.基于光譜相似性的高光譜圖像超分辨率算法.自動(dòng)化學(xué)報(bào),2014,40(12):2797?2807)

2 Ni Ding,Ma Hong-Bing.Spectral-spatial classi fi cation of hyperspectral images based on neighborhood collaboration.Acta Automatica Sinica,2015,41(2):273?284(倪鼎,馬洪兵.基于近鄰協(xié)同的高光譜圖像譜-空聯(lián)合分類.自動(dòng)化學(xué)報(bào),2015,41(2):273?284)

3 Du B,Zhang L P.A discriminative metric learning based anomaly detection method.IEEE Transactions on Geoscience and Remote Sensing,2014,52(11):6844?6857

4 Lin C H,Chi C Y,Wang Y H,Chan T H.A fast hyperplanebased minimum-volume enclosing simplex algorithm for blind hyperspectralunmixing.IEEE Transactions on Signal Processing,2016,64(8):1946?1961

5 Bioucas-Dias J M,Plaza A,Dobigeon N,Parente M,Du Q,Gader P,Chanussot J.Hyperspectralunmixing overview:geometrical,statistical,and sparse regression-based approaches.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2012,5(2):354?379

6 Lillesand T,Kiefer R W,Chipman J.Remote Sensing and Image Interpretation(Seventh Edition).New York:John Wiley and Sons,2014.23?61

7 Yuan Y,Fu M,Lu X Q.Substance dependence constrained sparse NMF for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2015,53(6):2975?2986

8 Wang W H,Qian Y T,Tang Y Y.Hypergraph-regularized sparse NMF for hyperspectralunmixing.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2016,9(2):681?694

9 Qu Q,Nasrabadi N M,Tran T D.Subspace vertex pursuit:a fast and robust near-separable nonnegative matrix factorization method for hyperspectralunmixing.IEEE Journal of Selected Topics in Signal Processing,2015,9(6):1142?1155

10 Li J,Bioucas-Dias J M,Plaza A,Liu L.Robust collaborative nonnegative matrix factorization for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2016,54(10):6076?6090

11 Feng R Y,Zhong Y F,Zhang L P.An improved nonlocal sparse unmixing algorithm for hyperspectral imagery.IEEE Geoscience and Remote Sensing Letters,2015,12(4):915?919

12 Iordache M D,Bioucas-Dias J M,Plaza A.Collaborative sparse regression for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2014,52(1):341?354

13 Meyer T R,Drumetz L,Chanussot J,Bertozzi A L,Jutten C.Hyperspectralunmixing with material variability using social sparsity.In:Proceedings of 2016 IEEE International Conference on Image Processing.Phoenix,USA:IEEE,2016.2187?2191

14 Giampouras P V,Themelis K E,Rontogiannis A A,Koutroumbas K D.Simultaneously sparse and low-rank abundance matrix estimation for hyperspectral image unmixing.IEEE Transactions on Geoscience and Remote Sensing,2016,54(8):4775?4789

15 Zheng C Y,Li H,Wang Q,Chen C L P.Reweighted sparse regression for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2016,54(1):479?488

16 Nascimento J M P,Dias J M B.Vertex component analysis:a fast algorithm to unmixhyperspectral data.IEEE Transactions on Geoscience and Remote Sensing,2005,43(4):898?910

17 Boardman J W.Automating spectral unmixing of AVIRIS data using convex geometry concepts.In:Proceedings of Summaries of the 4th Annual JPL Airborne Geoscience Workshop.Arlington,Virginia:JPL,1993.11?14

18 Winter M E.N-FINDR:an algorithm for fast autonomous spectral end-member determination in hyperspectral data.In:Proceedings of 1999 SPIE′s International Symposium on Optical Science,Engineering,and Instrumentation.Denver,USA:SPIE,1999.266?275

19 Chan T H,Chi C Y,Huang Y M,Ma W K.A convex analysis-based minimum-volume enclosing simplex algorithm for hyperspectralunmixing.IEEE Transactions on Signal Processing,2009,57(11):4418?4432

20 Li J,Agathos A,Zaharie D,Bioucas-Dias J M,Plaza A,Li X.Minimum volume simplex analysis:a fast algorithm for linear hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2015,53(9):5067?5082

21 Bioucas-Dias J M.A variable splitting augmented Lagrangian approach to linear spectral unmixing.In:Proceedings of the 1st Workshop on Hyperspectral Image and Signal Processing:Evolution in Remote Sensing.Grenoble,France:IEEE,2009.1?4

22 Boyd S,Parikh N,Chu E,Peleato B,Eckstein J.Distributed optimization and statistical learning via the alternating direction method of multipliers.Foundations and Trends in Machine Learning,2010,3(1):1?122

23 Chan T H,Ma W K,Chi C Y,Wang Y.A convex analysis framework for blind separation of non-negative sources.IEEE Transactions on Signal Processing,2008,56(10):5120?5134

24 StrangG.Linear Algebra and Its Applications.San Diego,CA:Thomson,2005.

25 Chen S S,Donoho D L,Saunders M A.Atomic decomposition by basis pursuit.SIAM Review,2001,43(1):129?159

26 Shi Z W,An Z Y,Tan X Y,Zhu Z X,Jiang Z G.Hyperspectralunmixing using non-negative matrix factorization with automatically estimating regularization parameters.In:Proceedings of the 7th International Conference on Natural Computation.Shanghai,China:IEEE,2011.1836?1840

27 Chen Bao-Lin.Optimization Theories and Algorithms(Second Edition).Beijing:Tsinghua University Press,2005.180?193(陳寶林.最優(yōu)化理論與算法.第2版.北京:清華大學(xué)出版社,2005.180?193)

28 Clark R N,Swayze G A,Wise R,Livo E,Hoefen T,Kokaly R,Sutley S J.USGS digital spectral library[Online],available:http://speclab.cr.usgs.gov/spectral-lib.html,September 13,2016.

29 AVIRIS.AVIRIS data-ordering free AVIRIS standard data products[Online],available:http://aviris.jpl.nasa.gov/html/aviris.freedata.html,September 13,2016.

30 Clark R N,Swayze G A,Livo K E,Kokaly R F,Sutley S J,Dalton J B,McDougal R R,Gent C A.Imaging spectroscopy: earth and planetary remote sensing with the usgstetracorder and expert systems[Online],available:http://speclab.cr.usgs.gov/PAPERS/tetracorder,September 13,2016.

A Robust Minimum Volume Based Algorithm with Automatically Estimating Regularization Parameters for Hyperspectral Unmixing

WANG Tian-Cheng1LIU Xiang-Zhen1DONG Ze-Zheng1WANG Hai-Bo1

Hyperspectral unmixing aims at fi nding hidden endmembers and their corresponding abundances from hyperspectral images with low spatial resolution.Based on the well-known minimum volume(MV)rule in geometrical based approaches,a robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing(RMVHU)is proposed.In this algorithm,the ANC constraint is replaced with a negative number punished regularizer which may lead to a more robust result to outliers and noise.A cyclic minimization algorithm is used to split the nonconvex RMVHU problem into convex subproblems,and ADMM is referred to sovle the large scale optimization problem with the increasing number of pixels in the image.To improve the convergence of the algorithm,a strategy to estimate the regularization parameters of the regularizer automatically is proposed.Compared with some existing geometrical based methods,experimental results show the superiority of the RMVHU algorithm on both synthetic datasets and real datasets.

Hyperspectral unmixing,alternating direction method of multipliers(ADMM),convex optimization,minimum volume,automatically estimating regularization parameters

Wang Tian-Cheng,Liu Xiang-Zhen,Dong Ze-Zheng,Wang Hai-Bo.A robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing.Acta Automatica Sinica,2017,43(12):2141?2159

2016-09-13 錄用日期2016-12-10

September 13,2016;accepted December 10,2016

本文責(zé)任編委胡清華

Recommended by Associate Editor HU Qing-Hua

1.上海衛(wèi)星工程研究所上海201109

1.Shanghai Institute of Satellite Engineering,Shanghai 201109

王天成,劉相振,董澤政,王海波.一種自適應(yīng)魯棒最小體積高光譜解混算法.自動(dòng)化學(xué)報(bào),2017,43(12):2141?2159

DOI10.16383/j.aas.2017.c160653

王天成 上海衛(wèi)星工程研究所助理工程師,2016年獲哈爾濱工業(yè)大學(xué)工學(xué)碩士學(xué)位.主要研究方向?yàn)楦吖庾V解混、系統(tǒng)建模與仿真.本文通信作者.

E-mail:wangtcsa@163.com

(WANG Tian-Cheng Assistant engineer at Shanghai Institute of Satellite Engineering.He received his master degree from Harbin Institute of Technology in 2016.His research interest covers hyperspectral unmixing,system modeling and simulating.Corresponding author of this paper.)

劉相振 上海衛(wèi)星工程研究所高級(jí)工程師,2006年獲上海航天技術(shù)研究院工學(xué)碩士學(xué)位.上海衛(wèi)星工程研究所信息與仿真中心主任.主要研究方向?yàn)閺?fù)雜系統(tǒng)的建模與仿真.

E-mail:sirruslxz@163.com

(LIUXiang-Zhen Senior engineer at Shanghai Institute of Satellite Engineering.He received his master degree from Shanghai Academy of Space fl ight Technology in 2006.Director of the Information and Simulation Center for Satellite Engineering,Shanghai Institute of Satellite Engineering.His research interest covers modeling and simulating of complex system.)

董澤政 上海衛(wèi)星工程研究所工程師,2010年獲南京航空航天大學(xué)工學(xué)碩士學(xué)位.主要研究方向?yàn)橄到y(tǒng)建模與仿真.

E-mail:dzzh520@hotmail.com

(DONG Ze-ZhengEngineer at Shanghai Institute of Satellite Engineering.He received his master degree from Nanjing University of Aeronautics and Astronautics in 2010.His research interest covers system modeling and simulating.)

王海波 上海衛(wèi)星工程研究工程師,2012年獲哈爾濱工業(yè)大學(xué)工學(xué)博士學(xué)位.主要研究方向?yàn)橄到y(tǒng)建模與仿真.

E-mail:13766898363@163.com

(WANG Hai-BoEngineer at Shanghai Institute of Satellite Engineering.He received his Ph.D.degree from Harbin Institute of Technology in 2012.His research interest covers system modeling and simulating.)