王為家，耿修瑞

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

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

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

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

1 背景知識

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

(1)

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

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

(2)

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

(3)

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

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

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

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

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

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

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

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

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

(4)

其中：

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

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

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

(5)

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

(6)

于是重建誤差可以表示為

(7)

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

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

2.2 EIC-OSV的建立及求解

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

(8)

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

(9)

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

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

(10)

式(10)等價于

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

(11)

其中：

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

(12)

其中：

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

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

(13)

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

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

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

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

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

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

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

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

3 實驗部分

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

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

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

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

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

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

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

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

圖5 實驗用數(shù)據(jù)

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

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

(14)

其中

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

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

圖6 豐度誤差對比

表1 端元提取光譜角度差

3.3 實驗結(jié)果討論

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

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

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

4 結(jié)論

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