王為家，耿修瑞

(中國(guó)科學(xué)院空天信息創(chuàng)新研究院 中國(guó)科學(xué)院空間信息處理與應(yīng)用系統(tǒng)技術(shù)重點(diǎn)實(shí)驗(yàn)室， 北京 100094； 中國(guó)科學(xué)院大學(xué)， 北京 100049) (2020年1月13日收稿； 2020年3月26日收修改稿)

由于傳感器空間分辨率和地面物體多樣性的限制，混合像元問題廣泛存在于高光譜遙感數(shù)據(jù)中。因此，混合像元分解已成為高光譜圖像分析的最重要步驟之一。在線性混合模型下，混合像元是相應(yīng)的豐度系數(shù)加權(quán)的幾個(gè)單一光譜特征(稱為端元)的線性組合。

在線性混合模型的假設(shè)下，端元提取算法可以分為端元選擇和端元生成兩類[1]。端元選擇算法假設(shè)數(shù)據(jù)集中存在恰好為端元的像元，通過一定方法將最有可能成為單形體頂點(diǎn)的像元選擇出來。比如純像元指數(shù)算法(pixel purity index，PPI)[2]，NFINDR[3]，以光譜信息熵改進(jìn)的NFINDR[4]，單形體膨脹算法(simplex growing algorithm, SGA)[5]，頂點(diǎn)分析法(vertex component analysis, VCA)[6]，高斯消去法(Gaussian elimination method, GEM)[7]，基于Gram行列式的快速端元確定法(fast gram determinant based algorithm, FGDA)[8]，基于顯著性的端元檢測(cè)(saliency-based endmember detection, SED)[9]，基于K-means聚類的端元提取[10]等。除此之外高光譜圖像的聚類也可用于端元選擇，譬如譜聚類[11]。當(dāng)然對(duì)于實(shí)際的高光譜數(shù)據(jù)集，純像元可能是不存在的，為了解決這種情況新的端元生成算法應(yīng)運(yùn)而生。端元生成算法在本質(zhì)上就是找到盡可能包含所有像元的體積最小的單形體。譬如最小體積轉(zhuǎn)換算法(minimum volume transformation, MVT)[12]，迭代約束端元算法(the iterative constrained endmember method, ICE)[13]，稀疏性ICE(sparsity-promoting ICE, SPICE)[14]，基于非負(fù)矩陣分解的最小體積約束算法(minimum volume constraint NMF, MVC-NMF)[15],最小封閉單形體體積(minimum volume enclosing simplex, MVES)[16-17]等，以及本文提到的幾何優(yōu)化模型(the geometric optimization model, GOM)[18]。

在不滿足純像元條件下，現(xiàn)有的端元生成方法一般通過使重建誤差最小化來得到最終的端元。但是在實(shí)際的高光譜圖像中一般是存在噪聲的，噪聲的存在會(huì)污染原有的像元，使得最小化重建誤差的提取結(jié)果過擬合，失去其應(yīng)有的物理意義。因此很多端元生成方法都引入了一些約束(例如體積約束)來限制過擬合。事實(shí)上，單形體體積的最小化本就是端元特征的體現(xiàn)。我們發(fā)現(xiàn)可以通過約束重建誤差不變以得到更小的單形體體積。重建誤差的大小與噪聲強(qiáng)度是正相關(guān)的，以重建誤差作為約束會(huì)有更好的魯棒性?；诖颂岢龌诠潭Ｐ驼`差的最優(yōu)單形體體積模型，旨在對(duì)端元提取的結(jié)果進(jìn)一步優(yōu)化。

1 背景知識(shí)

線性混合模型認(rèn)為高光譜圖像可以看成豐度矩陣與端元的線性組合加上誤差，線性混合模型滿足以下關(guān)系

(1)

其中：X為l(l為影像波段數(shù))維混合像元光譜，是已知觀測(cè)量；A為l×p(p為端元數(shù)目)端元矩陣或源矩陣，其中每一列為一個(gè)端元的光譜向量；豐度矩陣S為該像元對(duì)應(yīng)各個(gè)端元的豐度；n為l維高斯隨機(jī)噪聲或模型誤差。對(duì)于線性光譜混合模型，豐度矩陣S中的每個(gè)向量滿足2個(gè)限制條件：非負(fù)性限制(abundance non-negativity constraint, ANC)、和為1限制(abundance sum-to-one constraint, ASC),即每個(gè)像元的豐度系數(shù)均大于0且和為1。

給定端元矩陣A與豐度矩陣S，重建誤差可以表示為

(2)

使用V表示單形體體積，統(tǒng)一地看，端元提取的最終目的均是優(yōu)化如下目標(biāo)函數(shù)

(3)

這里介紹一下凸面幾何體理論，凸面幾何體理論是高光譜影像端元提取的理論基礎(chǔ)。如果空間中任意2點(diǎn)的連線仍然在該幾何體內(nèi)，那么該幾何體就是凸的。n維凸面幾何體是由n-1維的凸集構(gòu)成，而單形體是指n維空間中只有n+1個(gè)頂點(diǎn)的凸面幾何體。圖1是在2個(gè)波段的情況下，所有像元組成了一個(gè)三角形。A、B、C 3點(diǎn)是該凸面幾何體的頂點(diǎn)，即該影像的純像元位置；混合像元?jiǎng)t分布在三角形的內(nèi)部，每個(gè)像元均可以由3點(diǎn)線性表示。光譜矢量位于單形體內(nèi)，如果線性混合模型的所有假設(shè)成立，那么其頂點(diǎn)就是要找的端元，端元提取相當(dāng)于識(shí)別該單形體的頂點(diǎn)。

圖1 單形體(頂點(diǎn)代表端元)

最小化重建誤差的意義是為了滿足和為1約束與非負(fù)性約束；體積約束是為了單形體更緊密地包含像元，從而使得端元提取結(jié)果有物理意義。在理想條件下，純像元存在、噪聲忽略不計(jì)，那么重建誤差項(xiàng)為0，這種情況下端元選擇例如NFINDR、PPI等算法就可以準(zhǔn)確定位端元位置。在非理想條件下需要考慮噪聲對(duì)高光譜圖像的污染與純像元不存在的情況。而端元選擇算法無法兼顧所有的約束條件，僅僅將已有的像元作為提取結(jié)果，導(dǎo)致端元提取結(jié)果不準(zhǔn)確。

端元生成算法(例如MVC-NMF)為了彌補(bǔ)上述算法的缺陷，同時(shí)考慮重建誤差與體積的約束，但是對(duì)于式(3)中參數(shù)λ的選擇沒有判據(jù)，生成結(jié)果的好壞極大程度地依賴λ的選擇。

2 固定模型誤差優(yōu)化單形體體積

2.1 幾何優(yōu)化模型(GOM)

幾何優(yōu)化模型利用一種新的衡量重建誤差的方式，它放棄了線性混合模型中同時(shí)使用端元矩陣與豐度系數(shù)表示重建誤差的方式，而是采用只有單一變量端元矩陣A來表示重建誤差。

如圖2所示，一幅高光譜圖像端元矩陣為A，端元數(shù)量為p，圖中任一像元xi的豐度系數(shù)si可以寫成

圖2 幾何優(yōu)化模型

(4)

其中：

Aij=[a1,…,aj-1,xi,aj+1,…,ap]，

其幾何意義為將端元矩陣A中第i個(gè)端元換為xi后重新計(jì)算的端元矩陣的體積。

可以證明，當(dāng)像元xi位于由端元構(gòu)造的單形體內(nèi)時(shí)，從式(4)推導(dǎo)出的點(diǎn)xi的體積比的系數(shù)總和恰好為1。實(shí)際上，可以發(fā)現(xiàn)對(duì)于任意單形體體內(nèi)的像元xi，均滿足

(5)

反之，如果像元xi在單形體外部，那么應(yīng)該滿足

(6)

于是重建誤差可以表示為

(7)

其中：M是像元的數(shù)量，p是端元數(shù)量。

GOM的關(guān)鍵優(yōu)勢(shì)在于它只有1個(gè)變量，即端元矩陣A，因此它避免了由于豐度矩陣S引起的所有問題。

2.2 EIC-OSV的建立及求解

單形體體積在端元提取中具有重要的物理意義。端元提取實(shí)際上是在單形體體積的最小化與重建誤差的最小化之間取得平衡。本文選擇約束重建誤差不變來獲取最優(yōu)的體積，這樣有2個(gè)優(yōu)勢(shì)：一是重建誤差有理想值0，規(guī)避了初值選取問題；二是保證了單形體體積不被約束，從而使求解結(jié)果有物理意義?；谶@個(gè)思路可以建立如下EIC-OSV模型：

(8)

該模型根據(jù)拉格朗日乘子法可以轉(zhuǎn)化為如下無約束優(yōu)化問題：

(9)

其中：ε0為重建誤差的初值，在無噪聲時(shí)為0，考慮噪聲時(shí)為某一固定值。為保證重建誤差約束，需要滿足目標(biāo)函數(shù)H的導(dǎo)數(shù)與重建誤差函數(shù)的導(dǎo)數(shù)互相垂直，即兩者的內(nèi)積為0：

〈H′A(A,λ),ε′(A)〉=0，

(10)

式(10)等價(jià)于

〈V′(A)+λε′(A),ε′(A)〉=0，

(11)

其中：

那么λ就可以根據(jù)重建誤差約束求解：

(12)

其中：

求解方法選用梯度下降法。考慮到求解過程需要一直保證重建誤差空間與目標(biāo)函數(shù)空間正交，因此每一步迭代都需要求解一次λ，每次更新λ之后，利用下式迭代更新端元矩陣A：

A=A-ηH′A(A,λ).

(13)

其中η為步長(zhǎng)。步長(zhǎng)的選擇在該優(yōu)化模型中至關(guān)重要，如果步長(zhǎng)足夠小，則可以保證每一次迭代都減少，但可能導(dǎo)致收斂太慢；如果步長(zhǎng)太大，則不能保證每一次迭代都減少，也不能保證收斂。最終的求解結(jié)果需要多次迭代直至目標(biāo)函數(shù)的變化小于一定閾值。EIC-OSV的本質(zhì)是在重建誤差梯度的正交補(bǔ)空間進(jìn)一步縮小端元體積。只要兩者梯度不是嚴(yán)格平行的，端元體積存在可以下降的梯度方向，就可以利用EIC-OSV進(jìn)行迭代優(yōu)化。

2.3 信噪比與重建誤差的關(guān)系

在不考慮噪聲時(shí)，如果能夠準(zhǔn)確地找到端元，那么也可以準(zhǔn)確地實(shí)現(xiàn)重建誤差為0。然而一般情況下，完全理想的高光譜圖像是不存在的，在有噪聲干擾時(shí)，即使找到的端元是完美的真實(shí)端元，圖像的重建誤差仍然存在。

圖3列舉了添加4種不同幅度噪聲(信噪比依次為40、20、15、10 dB)后的像元分布情況。理想數(shù)據(jù)在添加不同程度的噪聲之后，各個(gè)像元的分布在各個(gè)方向有不同的波動(dòng)，很多像元落在真實(shí)端元圍成的單形體外，因此幾何優(yōu)化模型的重建誤差不為0。噪聲越大，落在單形體外的像元數(shù)量越多，距離越遠(yuǎn)。

接下來討論對(duì)于給定信噪比的高光譜圖像，如何在端元提取中降低噪聲的影響。對(duì)無噪聲的數(shù)據(jù)，可以約束重建誤差為0；對(duì)有噪聲的數(shù)據(jù)，可以采取如下措施：

1) 根據(jù)給定的高光譜圖像的信噪比、波段數(shù)、分辨率、像元個(gè)數(shù)、端元數(shù)量等已有條件，判斷在這一系列條件下真實(shí)端元的重建誤差。

2) 選擇合適的端元矩陣，使其重建誤差等于步驟1)中的重建誤差，將該端元矩陣作為初始端元矩陣，利用EIC-OSV模型進(jìn)行端元優(yōu)化。下面給出EIC-OSV實(shí)現(xiàn)的具體流程：

EIC-OSV:step1. 確定待提取端元數(shù)量p。step2. 將高光譜數(shù)據(jù)降維至p-1維。step3. 選擇重建誤差初值ε0。step4. 根據(jù)ε0確定合適的端元矩陣初值A(chǔ),可以人為選擇初始端元,也可以使用其他已有方法的端元提取結(jié)果作為初值。step5. 計(jì)算λ=-<ε'(A),ε'(A)>.step6. 計(jì)算A=A-ηH'A(A,λ).step7. 判斷V(A)是否收斂,若是,則結(jié)束;否則返回step5。

3 實(shí)驗(yàn)部分

3.1 模擬數(shù)據(jù)驗(yàn)證

在該測(cè)試中，使用來自美國(guó)地質(zhì)調(diào)查局?jǐn)?shù)字光譜庫的3個(gè)參考圖：明礬石、方解石和高嶺石，生成具有1 000個(gè)像元的混合數(shù)據(jù)。通過Dirichlet分布(參數(shù)[1/3,1/3,1/3])產(chǎn)生相應(yīng)的豐度系數(shù)。為了驗(yàn)證不同條件下的算法提取效果，對(duì)這組數(shù)據(jù)分別作2種處理：舍棄豐度系數(shù)大于0.9的像元(為了保證不滿足純像元條件)；在前者基礎(chǔ)上增加SNR=15 dB的高斯白噪聲(為了評(píng)估算法在噪聲方面的性能)。

如圖4(a)，對(duì)于不存在噪聲的理想情況，可以用一個(gè)足夠大的單形體作為初值，重建誤差約束為0，EIC-OSV就可以精準(zhǔn)地生成端元，優(yōu)化結(jié)果如圖4(b)所示。

不同信噪比的不同分布情況，已經(jīng)在圖3中給出。隨著噪聲強(qiáng)度的擴(kuò)大，真正的端元圍成的單形體舍棄的噪聲點(diǎn)逐漸增多?？梢约s束重建誤差等于該信噪比下的某一固定值，優(yōu)化單形體體積。這樣得到的單形體體積具有物理意義。

以加入信噪比為15 dB的噪聲為例，假設(shè)端元已知，多次生成固定信噪比的隨機(jī)數(shù)據(jù)，統(tǒng)計(jì)信噪比與模型誤差的關(guān)系。統(tǒng)計(jì)了1 000組數(shù)據(jù)下的重建誤差，其均值為0.015 1，標(biāo)準(zhǔn)差為0.001 3，重建誤差對(duì)于該組數(shù)據(jù)分布較為穩(wěn)定。事實(shí)上如果像元數(shù)量足夠大，重建誤差會(huì)趨于某個(gè)固定的值。如圖4(c)所示，給定一組已知信噪比為15 dB的待測(cè)數(shù)據(jù)。根據(jù)NFINDR的端元提取結(jié)果作為初始單形體，維持該單形體質(zhì)心不變進(jìn)行膨脹或者縮小，利用二分法多次確定某一初始端元矩陣使得其模型誤差約為0.015 1，利用得到的端元矩陣作為初值進(jìn)行迭代優(yōu)化。優(yōu)化后的端元提取結(jié)果如圖4(d) 所示?？梢钥闯鰞?yōu)化結(jié)果幾乎完全地契合了真實(shí)端元。圖4(e)反映了單形體體積與迭代次數(shù)的關(guān)系。隨著迭代次數(shù)的增加，單形體體積也在不斷減小，符合我們的期望。

(a)(b)為“純像元不存在”情況；(c)(d)(e)為“加入SNR=15 dB噪聲”情況

3.2 真實(shí)數(shù)據(jù)驗(yàn)證

本節(jié)使用真實(shí)的高光譜數(shù)據(jù)集進(jìn)一步驗(yàn)證EIC-OSV的有效性。該數(shù)據(jù)是1997年6月19日在內(nèi)華達(dá)州的銅礦開采場(chǎng)由機(jī)載可見光紅外光譜儀(AVIRIS)拍攝。本次實(shí)驗(yàn)中使用的子數(shù)據(jù)如圖5所示。該數(shù)據(jù)包含370～2 500 nm范圍內(nèi)的總共224個(gè)波段?？紤]到存在水蒸氣吸收或低信噪比的情況，刪除了波段1～3、107～113、153～169、221～224。

圖5 實(shí)驗(yàn)用數(shù)據(jù)

根據(jù)虛擬維度(virtual dimension，VD)[19]以及實(shí)地的地面調(diào)查[20],提取端元數(shù)量設(shè)置為14。真實(shí)端元的光譜通過ENVI軟件中USGS數(shù)字光譜庫中對(duì)相應(yīng)的光譜重采樣獲得。實(shí)驗(yàn)中，首先使用PCA將數(shù)據(jù)降至13維，分別利用NFINDR、ICE、MVC-NMF提取14個(gè)端元，同時(shí)分別將這3種算法提取的端元矩陣作為EIC-OSV的初值，計(jì)算各自的重建誤差并進(jìn)行優(yōu)化。這樣可以大大減少迭代次數(shù)，方便選取迭代步長(zhǎng)。

使用均方根角度誤差(rmsSAE)來衡量端元提取結(jié)果的準(zhǔn)確性。rmsSAE反映了提取端元與真正端元的均方角度差異，定義如下：

(14)

其中

其物理意義為單一提取端元與真實(shí)端元之間的角度。

表1給出了端元提取的結(jié)果與真實(shí)端元之間的角度、均方根角度誤差。利用已經(jīng)提取到的端元，可以分別計(jì)算對(duì)應(yīng)的豐度矩陣以驗(yàn)證精度，其豐度誤差在圖6中給出。

圖6 豐度誤差對(duì)比

表1 端元提取光譜角度差

3.3 實(shí)驗(yàn)結(jié)果討論

對(duì)比以上實(shí)驗(yàn)結(jié)果，模擬數(shù)據(jù)的端元提取結(jié)果符合預(yù)期：當(dāng)不考慮噪聲或者噪聲的影響極小時(shí)，約束重建誤差為0，EIC-OSV可以作為一種理想數(shù)據(jù)的端元生成算法；在有噪聲影響時(shí)，相同端元矩陣對(duì)于固定信噪比的隨機(jī)數(shù)據(jù)，重建誤差在統(tǒng)計(jì)上是大致穩(wěn)定的，將該統(tǒng)計(jì)值作為EIC-OSV的初值進(jìn)行優(yōu)化，提取結(jié)果幾乎完全吻合于真實(shí)端元。這也體現(xiàn)了EIC-OSV的一個(gè)顯著優(yōu)勢(shì)：它的約束值——重建誤差與噪聲強(qiáng)度是正相關(guān)的，當(dāng)確定這個(gè)約束值之后進(jìn)行優(yōu)化，能極大限度地規(guī)避噪聲帶來的影響。

真實(shí)數(shù)據(jù)的實(shí)驗(yàn)中，對(duì)比表1的各組數(shù)據(jù)?？梢园l(fā)現(xiàn)3種初始方案中NFINDR的提取結(jié)果較差，這可能是圖像中不存在純像元導(dǎo)致的，而ICE、MVC-NMF均考慮了純像元不存在的情況。對(duì)于這3種不同的初始方案，無論是rmsSAE還是豐度誤差，EIC-OSV均實(shí)現(xiàn)了優(yōu)于初值的提取結(jié)果。而且對(duì)于更好的初值，最終的優(yōu)化結(jié)果也更接近真實(shí)端元，這也意味著該初值對(duì)應(yīng)的ε0更接近真實(shí)值。值得注意的是，EIC-OSV的優(yōu)化目標(biāo)為整個(gè)端元矩陣的單形體體積，因此并不是所有端元的提取結(jié)果都會(huì)得到改善。

由于儀器參數(shù)的不確定，可用數(shù)據(jù)集受限以及實(shí)地考察結(jié)果的真實(shí)性，并不能完美地?cái)M合初始重建誤差。但考慮到EIC-OSV是在重建誤差梯度的正交補(bǔ)空間進(jìn)一步迭代縮小端元體積，其優(yōu)化目標(biāo)總是迫使迭代結(jié)果更加逼近真實(shí)端元的方向。因此理論上EIC-OSV可以對(duì)現(xiàn)有的所有端元提取算法進(jìn)行優(yōu)化，只要目標(biāo)函數(shù)梯度與模型誤差梯度不平行，就可以通過EIC-OSV優(yōu)化得到更小的單形體體積，使得噪聲對(duì)端元提取結(jié)果的影響降到最低，從而保證端元提取結(jié)果更有物理意義。

4 結(jié)論

隨著遙感技術(shù)的快速發(fā)展，高光譜圖像的應(yīng)用范圍越來越廣，混合像元問題普遍存在于高光譜數(shù)據(jù)中。作為混合像元分析的重要步驟，端元提取算法側(cè)重于最小化重建誤差來提升端元提取的精度，無法解決噪聲帶來的過擬合問題。本文提出的EIC-OSV模型沿用幾何優(yōu)化模型的重建誤差，避開了求解豐度系數(shù)矩陣的繁瑣步驟。針對(duì)重建誤差的約束來最小化單形體體積，將信噪比與重建誤差建立關(guān)聯(lián)，保障在噪聲條件下求解的結(jié)果仍然具有物理意義。在模擬數(shù)據(jù)及真實(shí)數(shù)據(jù)的提取結(jié)果均表現(xiàn)良好。對(duì)于真實(shí)數(shù)據(jù)，本文實(shí)現(xiàn)了提取結(jié)果的優(yōu)化，如果要獲得更理想的結(jié)果需要定量探討信噪比與重建誤差的關(guān)系，以及其他因素比如端元數(shù)量、像元數(shù)量、豐度分布等因素的影響，這也是將來工作致力的方向。