基于豐度分布約束的NMF端元生成方法*

石悅，王宏琦，郭新毅

(1 中國科學(xué)院電子學(xué)研究所， 北京 100190； 2 中國科學(xué)院網(wǎng)絡(luò)信息體系技術(shù)科技創(chuàng)新重點實驗室， 北京 100190；3 中國科學(xué)院大學(xué)， 北京 100049； 4 國電聯(lián)合動力技術(shù)有限公司， 北京 100039)(2018年11月12日收稿； 2019年3月18日收修改稿)

由于高光譜傳感器空間分辨率的限制以及地物的復(fù)雜多樣性，混合像元普遍存在于高光譜遙感圖像中，因此混合像元分析受到越來越廣泛的關(guān)注。線性混合模型(the linear mixture model，LMM)是國內(nèi)外混合像元分解方法中應(yīng)用最廣泛且最簡單的光譜混合模型。該模型把遙感影像的像元光譜看作是各種不同類型的純像元(稱為端元)光譜的線性組合?；贚MM的端元提取方法主要分為端元識別類方法和端元生成類方法。端元識別類方法包括純像元指數(shù)法(pixel purity index, PPI)[1]、正交子空間投影法(orthogonal subspace projection, OSP)[2]、頂點成分分析法(vertex component analysis, VCA)[3]、高斯消元法(Gaussian elimination method,GEM)[4]、N-FINDR[5]、基于Gram行列式的快速端元提取算法(FGDA)[6]、單純形生長算法(simplex generation algorithm, SGA)[7]和基于單形體的端元提取方法(simplex volume algorithm, SVA)[8]等。此種方法需要純像元假設(shè)，即假設(shè)一個場景中每個端元至少有一個純像元。端元生成類方法不需要純像元假設(shè)，主要包括最小體積變換[9]、迭代約束端元法(iterated constrained endmember, ICE)[10]、稀疏優(yōu)化ICE(sparse ICE, SPICE)[11]和幾何優(yōu)化模型(geometry optimization model, GOM)[12]等。

近年來，非負矩陣分解NMF(non-negative matrix factorization)的端元生成方法受到越來越多的關(guān)注。NMF端元生成方法可以同時估計出端元矩陣和豐度矩陣。Lee和Seung[13]證明NMF端元生成方法可利用乘式迭代高效快速地實現(xiàn)NMF目標函數(shù)的優(yōu)化迭代。由于NMF目標函數(shù)是非凸的，NMF方法容易陷入局部極值[14]。為了緩解NMF方法固有的局部極值問題，通常采用增加約束來優(yōu)化NMF端元生成方法。常用的優(yōu)化約束包括最小體積約束(minimum volume constraint, MVC)[15]、豐度分段平滑約束[16]]、端元不相似性約束[17]、豐度稀疏約束[18-19]等?；谝陨霞s束的NMF方法在一定程度上緩解了NMF方法的局部極值問題，但由于約束的引入，如經(jīng)典MVC-NMF方法，乘式迭代規(guī)則被破壞，算法計算復(fù)雜度大幅增加。我們提出一種基于豐度分布約束的NMF方法，通過矩陣跡運算表達豐度分布離散度約束優(yōu)化端元提取精度。同時可利用乘式迭代實現(xiàn)目標函數(shù)的迭代優(yōu)化，算法計算復(fù)雜度小。

1 基于豐度分布約束的NMF端元生成方法

1.1 線性混合模型

假設(shè)線性混合模型成立，圖像中每個像元向量可以表示為幾個端元(純像元)向量的線性組合，而混合系數(shù)為各個端元的豐度。線性混合模型表示如下：

(1)

(2)

0≤cij≤1.

(3)

式中：N為像素數(shù)，p為端元數(shù)，R=[r1,r2,…,rN]是大小為L×N的高光譜圖像(L為波段數(shù))。ci=[ci1,ci2,…,cip]T和cij分別為像元ri的豐度向量以及端元ej在像元ri處的豐度值。E=[e1,e2,…,ep]是L×p的端元矩陣，C=[c1,c2,…,cN]是p×N的豐度矩陣，N=[n1,n2,…,nN]是噪聲矩陣。由于像元的物理意義約束，豐度滿足2個約束條件：豐度和為一約束(abundance sum-to-one constraint，ASC)、豐度非負約束(abundance non-negative constraint，ANC)。

1.2 豐度分布約束

由于觀測像元的物理意義約束，任意一個像元ri的豐度值ci=[ci1,ci2,…,cip]T非負，且滿足sum(ci1,ci2,…,cip)=1。實際應(yīng)用中，并非所有的端元都對某個混合像元有貢獻，且每個端元的貢獻系數(shù)也分主次，換句話說，混合像元往往是由某個端元子集混合而成，其他端元的豐度值基本接近于0。Cichocki等[20]提出的HALS盲信號分離方法也進一步印證，在豐度矩陣C分布離散程度高的情況下，可以獲取較好的盲信號分離結(jié)果。

圖1 豐度分布示意圖Fig.1 Sketch map of abundance distribution

假設(shè)所有像元由3個端元構(gòu)成，如圖1所示。端元構(gòu)成的單形體頂點可以是{A,B,C}，也可以是{A′,B′,C′}。在端元單形體ABC中，像元P的豐度向量cABC=[0.8,0.1,0.1]T，豐度向量cABC的距離為0.66；相應(yīng)地，P在端元單形體ABC中，對應(yīng)的豐度向量cA′B′C′=[0.3,0.3,0.4]T，豐度向量的距離cA′B′C′等于0.34。因此，豐度向量分布約束可由豐度向量的分布離散程度，即豐度向量的距離表示。很明顯，豐度向量cABC的分布離散程度優(yōu)于豐度向量cA′B′C′。因此豐度向量的分布離散程度越高，即豐度向量的距離越大，解混效果越好。

豐度分布離散度可以由豐度矩陣C的跡運算表示，即

JC(C)=trace(CTC).

(4)

因此，我們將豐度分布離散度作為一項約束應(yīng)用于NMF端元生成方法。

1.3 基于豐度分布約束的NMF端元生成方法

在NMF標準方法的基礎(chǔ)上，增加豐度分布約束，本方法形成新的目標函數(shù)：

E≥0，C≥0.

(5)

(6)

(7)

通過計算，可知

(8)

因此，

(9)

因此，目標函數(shù)的迭代優(yōu)化可通過乘式迭代實現(xiàn)，迭代規(guī)則可以表示為：

(10)

本方法可利用乘式迭代規(guī)則實現(xiàn)端元目標函數(shù)迭代，處理效率較高，且得到的單形體滿足以下兩個特點：第一是基于線性混合模型下估計和真實數(shù)據(jù)間的誤差最小，類似于外力，使得端元單形體包含盡量多像元；第二是豐度矩陣離散度最高，類似于內(nèi)力，使得端元單形體盡可能緊地包裹像元。

2 實驗與分析

利用仿真數(shù)據(jù)和真實數(shù)據(jù)開展了一系列的實驗。為驗證本方法的有效性，與3種比較有代表性的方法(NFINDR-FCLS[5]、MVC-NMF[15]和NMF-ICE[21])進行比較。本方法采用NFINDR的端元提取結(jié)果作為初始值。

2.1 仿真數(shù)據(jù)

利用美國地質(zhì)調(diào)查(USGS)數(shù)字光譜庫中4種礦物的光譜作為端元。光譜數(shù)據(jù)共包含224個光譜波段，覆蓋波長范圍為0.38～2.5 μm，光譜分辨率為10 nm。仿真圖像數(shù)據(jù)大小為58×58，利用Dirichlet分布生成豐度系數(shù)。為進一步去除純像元，豐度系數(shù)均小于0.9。此外，仿真數(shù)據(jù)還添加了均值為零的白高斯噪聲，以模擬可能存在的誤差和傳感器噪聲。為了公正地評價所提出方法的有效性，在SNR=20 的情況下開展對比實驗。其中NFINDR-FCLS、MVC-NMF、NMF-ICE分別依據(jù)參考文獻[5,15,21]進行參數(shù)選取，本方法選擇迭代次數(shù)為150次，α=0.01，σ=20。

本方法利用文獻[14]中的rmsSAD(光譜角度散度)、rmsAAD(豐度角度散度)、rmsSID(頻譜信息散度)、rmsAID(豐度信息散度)評估方法的性能。從圖2和圖3可以看出，在NFINDR端元提取結(jié)果作為初始值的條件下，本文方法具有明顯優(yōu)勢。

圖2 誤差值對比結(jié)果Fig.2 Comparison among different methods

圖3 本方法端元生成結(jié)果與真實端元比較Fig.3 Comparison of the extracted endmembers with real endmembers

此外, 還通過豐度RMSE分析噪聲對本文方法的影響，在模擬數(shù)據(jù)中添加不同的高斯噪聲(SNR=10、15、20、25、30 dB)。從圖4可以看出，本文方法生成的豐度RMSE值較小，且在SNR=20 dB以后，基本不受噪聲影響。

圖4 不同噪聲條件下豐度RMSE結(jié)果Fig.4 RMSE results of abundances at different SNR values

我們提出的端元生成方法支持利用乘式迭代進行端元優(yōu)化迭代。由表1可知，在4核處理器主頻2.3 GHz的條件下，本方法在圖1所示的端元提取精度條件下，與其他方法相比處理時間大幅減少。

表1 模擬數(shù)據(jù)算法處理時間Table 1 Processing time for different methods with the synthetic image

同時我們也分析了本方法在不同初始值條件下的端元提取效果。從表2可以看出，由于本文方法未考慮現(xiàn)有方法提出的約束，需要較好的初始值才能達到較好的端元生成結(jié)果。

表2 本方法在不同初始值條件下的性能評估Table 2 Performance under different initializations

2.2 真實數(shù)據(jù)

這里使用的高光譜圖像數(shù)據(jù)，為ENVI軟件提供的機載可見光紅外成像光譜儀(AVIRIS)赤銅礦數(shù)據(jù)集。該數(shù)據(jù)集是一個400×350圖像，由1.990 8～2.479 0 μm的短波紅外線(SWIR)50波段組成。該場景已被廣泛用于驗證端元提取算法的性能。這個圖像數(shù)據(jù)在第44個波段存在一個異常像素，反射率值為0.8左右，是由于噪聲破壞而造成，如圖5所示。

圖5 真實實驗數(shù)據(jù)Fig.5 Real image data

將提取的端元與美國地質(zhì)勘察(USGS)數(shù)字光譜庫的對應(yīng)真實光譜進行比較，提取8個端元，并利用端元rmsSAD[14]平均值和處理時間評價本文方法的處理精度和處理效率。

設(shè)定提取的端元數(shù)為8，4類方法提取出的8個端元分別與光譜庫中的真實數(shù)據(jù)進行比較。從表3可以看出，與其他方法相比，本文方法rmsSAD值優(yōu)于或接近其他方法。

在4核處理器主頻2.3 GHz的條件下，本次實驗對各方法在以上處理精度下的處理時間進行了統(tǒng)計分析。從表4可以看出，本文方法既能有效提取端元，且處理效率優(yōu)于其他方法。

表3 不同算法處理精度Table 3 Performance evaluation for different methods

表4 真實數(shù)據(jù)算法處理時間Table 4 Processing time for different methods with the real image

3 結(jié)束語

基于約束的NMF方法緩解了NMF目標函數(shù)非凸帶來的局部極值問題，往往也造成乘式迭代規(guī)則被破壞，端元生成效率降低等問題。本文提出一種基于豐度分布約束的NMF方法，首次利用矩陣跡運算表征豐度分布離散度，可利用乘式迭代實現(xiàn)目標函數(shù)的迭代優(yōu)化，處理效率較高。通過模擬數(shù)據(jù)和真實數(shù)據(jù)從處理精度、處理時效性和噪聲魯棒性等方面進行實驗分析，本文方法可提高端元估計結(jié)果的精度和算法處理效率。由于未考慮現(xiàn)有方法在端元方面的約束，本方法對初始值依賴較高，希望在以后的研究中進一步改進。