徐 君，王彩玲

(1. 西安航空學(xué)院電子工程學(xué)院，陜西 西安 710077; 2. 西安石油大學(xué)計(jì)算機(jī)學(xué)院，陜西 西安 710065)

高光譜圖像雖然具有很高的光譜分辨率，但空間分辨率卻普遍不高，這就導(dǎo)致高光譜圖像中的每個(gè)像元可能不僅包含一種地物，而是多種地物覆蓋類型的混合。對(duì)于這種光譜混合像元，無法根據(jù)提取的光譜特性判斷其物質(zhì)屬性，也不能在分類過程中簡(jiǎn)單地將其歸于某一類地物?；旌舷裨诟吖庾V圖像中的普遍存在使得傳統(tǒng)的像元級(jí)圖像分類和分割精度難以提高，是高光譜遙感技術(shù)向定量化方向發(fā)展的嚴(yán)重障礙。混合像元分解是解決該問題的有效方法，它能夠在亞像元尺度上獲取像元內(nèi)部地物分布的真實(shí)信息，提高圖像分類的精度。

根據(jù)混合像元的產(chǎn)生機(jī)理，高光譜圖像混合模型可以分為線性混合和非線性混合兩種模型。在傳統(tǒng)的線性光譜混合像元分解算法中，一般假設(shè)圖像中的像元都包含全部的端元光譜，對(duì)每個(gè)端元給出一個(gè)對(duì)應(yīng)的豐度估計(jì)，但這往往與實(shí)際情況并不相符[1-4]。

實(shí)際上混合像元并不一定包含所有端元，而是由其中的一小部分端元混合而成，即混合像元只包含一個(gè)端元子集?；旌舷裨袑?shí)際所包含的端元也互不相同，因此混合像元所包含的端元是可變的。常用的線性混合模型在混合像元分解過程中將所有的端元都納入考慮范圍，對(duì)每個(gè)端元向量都進(jìn)行豐度估計(jì)，因此存在過擬合的分解誤差。為了解決混合像元端元可變問題，國(guó)內(nèi)外一些學(xué)者對(duì)如何選擇混合像元的最優(yōu)端元子集作了一些研究。國(guó)外學(xué)者主要提出了多端元光譜混合分析(MESMA)[5-6]、迭代光譜混合分析(ISMA)[7]和端元可變線性混合模型(EVLMM)[8]等最優(yōu)端元子集選擇方法。MESMA算法從光譜庫(kù)中選擇端元，假設(shè)每個(gè)像元最多由3個(gè)端元混合，嘗試對(duì)所有端元的組合進(jìn)行光譜分解，運(yùn)算效率較低，光譜分解的誤差較大。ISMA采用迭代算法去除豐度最小的端元，采用非限制性光譜分解方法進(jìn)行豐度估計(jì)，根據(jù)混合像元光譜與剩余端元重建光譜之間均方根誤差RMSE的變化情況確定最優(yōu)端元子集。EVLMM基于端元組合數(shù)量由少到多進(jìn)行搜索，直到采用更多的端元進(jìn)行光譜重建后所得RMSE比上一步大時(shí)停止。文獻(xiàn)[9]提出的OESSA算法首先利用所有端元進(jìn)行非限制性光譜分解，然后根據(jù)豐度值量化重要程度的大小，對(duì)端元進(jìn)行排序，由豐度最大的端元開始，逐漸增加端元個(gè)數(shù)解混，根據(jù)RMSE的統(tǒng)計(jì)結(jié)果確定最優(yōu)端元子集。文獻(xiàn)[10]提出了一種端元子集的局域確定方法，首先利用所有端元對(duì)混合像元進(jìn)行分解，根據(jù)該像元鄰域內(nèi)各端元的豐度分布狀況直接判定是否為純像元，如果不是純像元，則根據(jù)鄰域內(nèi)端元的豐度累加值大小排序，去掉累加值最小的端元，進(jìn)行精度更高的二次解混。混合像元的端元可變性在實(shí)際中對(duì)應(yīng)于混合像元的豐度矩陣具有稀疏特性。文獻(xiàn)[11—15]嘗試把稀疏原理與高光譜混合像元分解相結(jié)合，提出了基于稀疏約束的混合像元分解算法，基于稀疏約束的混合像元分解算法是在已知一個(gè)過完備的光譜庫(kù)的情況下，尋找方程的最簡(jiǎn)單的解。

以上方法從不同思路對(duì)混合像元的可變性進(jìn)行了研究，但實(shí)際上混合像元一般分布于不同地物的交界地帶，不同地物交界處的混合像元所包含的端元最有可能與周圍的地物分布狀況相關(guān)。針對(duì)這種情況，本文提出一種利用圖像的空間信息選取最優(yōu)端元組合的混合像元分解方法，在混合像元周圍局部鄰域提取純凈像元光譜，與從整幅圖像中提取到的全部端元進(jìn)行對(duì)比，根據(jù)分解后RMSE的變化情況優(yōu)化端元組合，對(duì)混合像元進(jìn)行分解，從而獲得最佳的光譜分解結(jié)果。

1 理論基礎(chǔ)

1.1 線性混合模型

線性混合模型(如圖1所示)是假定混合像元光譜由各端元的光譜通過線性疊加混合而成，由于其原理簡(jiǎn)單，物理含義明確，運(yùn)算效率高，分解精度基本滿足實(shí)際應(yīng)用的需求，因而獲得了廣泛的應(yīng)用。

圖1形象地闡述了線性光譜混合過程，在地面的3種物質(zhì)分別為a1、a2和a3，所占面積的比例分別為s1、s2和s3，陽光照射在這3種地物上，互不干擾，而在成像光譜儀上形成一個(gè)混合像素點(diǎn)，混合像元得到的反射率為x=a1s1+a2s2+a3s3，它是各個(gè)光譜特征的線性疊加。

線性光譜混合模型的數(shù)學(xué)表達(dá)式為

(1)

式中，x為混合像元的反射率；Ai為像元的第i個(gè)端元向量；si為該像元內(nèi)第i個(gè)端元的反射率比例(即豐度)；E為模型誤差；p為端元數(shù)目。

1.2 混合像元分解精度評(píng)價(jià)

對(duì)混合像元的分解結(jié)果所進(jìn)行的評(píng)價(jià)分為定性評(píng)價(jià)和定量評(píng)價(jià)兩種。定性評(píng)價(jià)主要是通過目視的方式對(duì)分解后的端元分布豐度圖進(jìn)行觀察，并與實(shí)際地物的分布狀況進(jìn)行對(duì)比。為了對(duì)混合像元的分解結(jié)果進(jìn)行準(zhǔn)確客觀的評(píng)價(jià)，需要建立相關(guān)的評(píng)價(jià)指標(biāo)進(jìn)行定量評(píng)價(jià)，定量評(píng)價(jià)通過嚴(yán)格的數(shù)學(xué)計(jì)算分解，給出精確的對(duì)比評(píng)價(jià)結(jié)果，因而獲得了廣泛的應(yīng)用。對(duì)真實(shí)高光譜圖像進(jìn)行混合像元分解時(shí)，端元光譜和端元對(duì)應(yīng)的豐度圖像信息一般是未知的，均方根誤差是一種比較簡(jiǎn)單、常用的混合像元分解精度定量評(píng)價(jià)指標(biāo)。

假設(shè)混合像元原始光譜向量為X，光譜分解所得的端元矩陣為A，豐度向量為S，則光譜重構(gòu)后的光譜向量為

(2)

該像元的余差為

(3)

對(duì)于n維的光譜向量，均方根誤差為

(4)

式中,εj為余差向量中的元素。獲得的RMSE值越小說明算法的混合像元分解精度越高，算法越有效。

2 光譜分解誤差隨端元子集變化情況

現(xiàn)有的研究表明，利用圖像中的所有端元對(duì)混合像元進(jìn)行分解，得到各端元的豐度值，然后從豐度最大的端元開始，按照豐度從大到小的順序，逐漸添加端元，建立并擴(kuò)充端元子集，用這些端元子集分別對(duì)該混合像元進(jìn)行二次分解，則混合像元分解的誤差總體上趨于減小[9]。

如對(duì)于整幅圖像包含10個(gè)端元、實(shí)際上由3個(gè)端元組成的混合像元，圖2給出了RMSE隨端元子集中端元數(shù)目增加的變化趨勢(shì)。

在混合像元分解時(shí)，如果端元子集中漏選實(shí)際存在的端元，那么光譜重構(gòu)圖像與原始圖像的RMSE會(huì)很大；當(dāng)端元子集中恰好包含實(shí)際存在的端元時(shí)，RMSE最??；而當(dāng)端元子集中多選了冗余端元進(jìn)行分解時(shí)為無偏估計(jì)，RMSE變化很小。

圖2中，按照豐度排序，在利用其中3個(gè)豐度最大的端元進(jìn)行混合像元二次分解后，相比上一步2個(gè)端元分解時(shí)重構(gòu)光譜與原始光譜之間的RMSE急劇下降，再向端元子集中增加端元數(shù)量進(jìn)行分解，RMSE的變化也不大，而是趨于平緩，由此可以判斷曲線“拐點(diǎn)”對(duì)應(yīng)的端元集就可以作為組成該混合像元的最優(yōu)端元子集，該混合像元實(shí)際包含的端元數(shù)為3。

3 結(jié)合圖像空間信息的最優(yōu)端元子集選擇

當(dāng)圖像中地物種類較多、存在較多端元的時(shí)候，如果通過全局搜索的方式在全部端元中搜索混合像元的最優(yōu)端元組合效率較低，因自然界地物往往具有一定的空間結(jié)構(gòu)，如農(nóng)田、森林、草地、建筑物，湖泊等，這些不同種類地物區(qū)塊之間光譜差異明顯。光譜混合像元一般位于地物區(qū)塊的邊緣交界地帶，在混合像元的周圍存在純凈像元的可能性較大，而這些純凈像元的光譜往往就是組成混合像元的端元光譜。充分利用圖像的空間信息，從混合像元的就近鄰域開始尋找純凈像元，搜索確定混合像元的最優(yōu)端元組合，可以提高搜尋混合像元最優(yōu)端元子集的效率。如圖3所示，本文提出一種結(jié)合圖像的空間信息尋找混合像元最優(yōu)端元子集的方法。具體步驟如下：

(1) 定義一個(gè)結(jié)構(gòu)元素K，K是一個(gè)中心像元坐標(biāo)為(x,y)、邊長(zhǎng)為k的正方形。讓結(jié)構(gòu)元素K在圖像上滑動(dòng)，在結(jié)構(gòu)元素覆蓋的區(qū)域內(nèi)計(jì)算所有像元光譜向量均值

(5)

計(jì)算結(jié)構(gòu)元素覆蓋區(qū)域內(nèi)k2個(gè)像元的光譜向量到均值光譜向量Ck的歐氏距離，記為dk。

(2) 根據(jù)凸面幾何學(xué)理論，最大dk對(duì)應(yīng)的像元必然是結(jié)構(gòu)元素覆蓋的局部區(qū)域內(nèi)“最純”的像元，將其記為第1個(gè)“端元”e1；然后將結(jié)構(gòu)元素k2個(gè)像元中距離e1最遠(yuǎn)的像元記為第2個(gè)“端元”e2；距離e1和e2所構(gòu)成直線最遠(yuǎn)的像元(這個(gè)點(diǎn)向e1和e2所構(gòu)成的直線作垂線距離最遠(yuǎn))記為第3個(gè)“端元”e3；距離e1、e2和e3所構(gòu)成平面最遠(yuǎn)的像元記為第4個(gè)“端元”e4；以此類推，最多得到k2個(gè)“端元”e1,e2,…,ek。

(3) 對(duì)步驟(2)中獲得的“端元”ei，與從整幅圖像中提取的所有端元ai(i=1,2,…,m)進(jìn)行對(duì)比，計(jì)算它們之間的光譜夾角距離

(6)

如果光譜夾角SADi小于設(shè)定的閾值，光譜高度相似，則可以將對(duì)應(yīng)的端元ai加入結(jié)構(gòu)元素K覆蓋區(qū)域中心位置處混合像元的最優(yōu)端元子集中。此步驟本質(zhì)上是通過結(jié)構(gòu)元素K在混合像元附近鄰域?qū)ふ摇凹兿裨惫庾V作為混合像元的端元光譜，這符合混合像元一般位于地物區(qū)塊邊緣交界地帶的規(guī)律，實(shí)際上是利用了圖像的空間信息就近尋找混合像元的最優(yōu)端元子集。

(4) 最優(yōu)端元子集建立后，應(yīng)用全限制線性混合模型對(duì)混合像元進(jìn)行分解，然后計(jì)算重構(gòu)光譜與原始光譜的RMSE。如果該RMSE近似等于全端元a1,a2,…,am分解時(shí)與原始光譜之間的RMSE，則當(dāng)前的端元子集即為最優(yōu)；如果該RMSE遠(yuǎn)大于全端元a1,a2,…,am分解時(shí)的RMSE，則說明端元子集中缺少實(shí)際存在的端元。

(5) 增大結(jié)構(gòu)元素K(x,y)的尺寸(增大邊長(zhǎng)k)，擴(kuò)大搜索范圍，返回步驟(1)重復(fù)執(zhí)行，直到RMSE減小到近似全端元分解時(shí)的RMSE為止，最終得到混合像元的最優(yōu)端元子集，進(jìn)而得到最優(yōu)的光譜分解結(jié)果。

圖3為擴(kuò)大結(jié)構(gòu)元素尋找最優(yōu)端子集示意圖。圖3(a)中結(jié)構(gòu)元素K1邊長(zhǎng)為3，結(jié)構(gòu)元素內(nèi)只包含兩種物質(zhì)的純像元，因此只能找出結(jié)構(gòu)元素中心位置處混合像元(x,y)包含的兩個(gè)端元e1和e2；擴(kuò)大結(jié)構(gòu)元素，使其邊長(zhǎng)為5(如圖3(b)所示)，結(jié)構(gòu)元素內(nèi)涵蓋3種物質(zhì)的純像元，即可找到混合像元(x,y)的全部端元e1、e2和e3。

4 試驗(yàn)與討論

4.1 模擬圖像試驗(yàn)

在USGS[16]礦物光譜庫(kù)中選擇4種礦物光譜作為端元，其光譜曲線如圖4所示。

添加SNR=30 dB的高斯噪聲，構(gòu)建大小為200×200的模擬圖像。如圖5所示， 模擬圖像4個(gè)角上分布著4個(gè)端元的純像元，陰影部分是由4種端元混合而成的混合像元。

以4個(gè)端元的其中一個(gè)進(jìn)行分析，圖6為端元①的實(shí)際豐度圖及分別采用UCLS、FCLS、迭代光譜混合分析(ISMA)及本文所提出的端元可變混合像元分解方法解混后的豐度圖。圖中顏色由淺到深表示豐度逐漸增加。

試驗(yàn)數(shù)據(jù)見表1，可以看出，UCLS由于沒有添加豐度約束條件，因而產(chǎn)生了最小的RMSE，由于對(duì)豐度和沒有進(jìn)行約束，豐度和在0.95～1.05之間的像元僅占圖像中像元的69.35%；FCLS嚴(yán)格滿足和為1的約束條件，但RMSE較大，分解精度較低。ISMA算法的RMSE也比FCLS大。本文提出的算法在RMSE上優(yōu)于ISMA算法，并且最后的分解結(jié)果豐度和滿足和為1的要求。

表1 模擬圖像試驗(yàn)結(jié)果

4.2 真實(shí)數(shù)據(jù)試驗(yàn)

美國(guó)內(nèi)華達(dá)州Cuprite地區(qū)礦物品種多樣，各種礦物露頭良好，該地區(qū)的AVIRIS數(shù)據(jù)包含224個(gè)波段(0.4～2.5 μm)，空間分辨率和光譜分辨率分別為20 m和10 nm，適合檢驗(yàn)混合像元分解算法的性能。去除其中受水蒸氣吸收影響較大和低信噪比的波段(1～3、104～113、148～167、221～224)，保留187個(gè)波段，經(jīng)過輻射校正得到反射率數(shù)據(jù)，并選取其中200×200大小的區(qū)域用于試驗(yàn)。圖7為該區(qū)域的假彩色合成圖與執(zhí)行UCLS、FCLS ISMA及本文算法所得到的其中一個(gè)端元豐度的分布圖。

表2是對(duì)試驗(yàn)所得結(jié)果的數(shù)據(jù)分析。結(jié)果表明，UCLS算法的RMES最小，但圖像中僅有11.8%的像元豐度和在0.8～1.2之間。FCLS算法雖然滿足豐度和為1的約束條件，但RMES最大。相比本文算法，ISMA算法的RMSE較大，豐度和也不滿足和為1的約束條件。本文方法在選定最優(yōu)端元子集進(jìn)行二次分解時(shí)，采用的是FCLS算法，豐度和為1。因此本文所提算法在性能上是優(yōu)于ISMA的。

表2 AVIRIS數(shù)據(jù)試驗(yàn)結(jié)果

本文方法利用圖像的空間信息，用結(jié)構(gòu)元素搜索周圍鄰域的“純像元”光譜，優(yōu)先組成最優(yōu)端元子集，進(jìn)行光譜分解。表2中真實(shí)高光譜數(shù)據(jù)的試驗(yàn)結(jié)果再次表明，采用本文方法優(yōu)化后的端元子集進(jìn)行光譜分解結(jié)果在RMSE上優(yōu)于FCLS和SALSA；且在二次分解時(shí)采用FCLS算法，因此豐度和為1，符合實(shí)際情況，總體上達(dá)到了較好的光譜分解效果。

5 結(jié) 語

混合像元分解是高光譜遙感信息提取的常用方法。傳統(tǒng)的分解方法一般假設(shè)混合像元包含圖像中的所有端元，但這并不符合實(shí)際情況，因此導(dǎo)致分解結(jié)果對(duì)噪聲敏感，分解精度不高。實(shí)際的混合像元只包含少數(shù)端元，由于地物空間分布上的連續(xù)性，混合像元所包含的端元往往與附近區(qū)域分布的地物有關(guān)。

本文結(jié)合圖像空間信息，用一個(gè)空間結(jié)構(gòu)元素在混合像元附近區(qū)域就近確定選取最優(yōu)端元子集。相比其他方法，該方法既不需要遍歷圖像中所有端元,且通過迭代運(yùn)算來尋找最優(yōu)端元子集，也不需要按照初次解混后的豐度大小進(jìn)行排序后對(duì)端元進(jìn)行取舍，而是將附近鄰域的“純像元”光譜優(yōu)先加入最優(yōu)端元子集，從最有可能屬于該混合像元的端元開始搜索添加，這符合現(xiàn)實(shí)地物分布的實(shí)際情況，有利于提高最優(yōu)端元子集的搜索效率，最終獲得較好的分解效果。