孫 肖， 徐林林， 王曉陽， 田 野， 王 偉， 張中躍

(1.中國地質(zhì)調(diào)查局廊坊自然資源綜合調(diào)查中心，廊坊 065000； 2.中國地質(zhì)大學(xué)(北京)土地科學(xué)技術(shù)學(xué)院，北京 100083)

0 引言

高光譜遙感具有“圖譜合一”的特點(diǎn)，已廣泛應(yīng)用于礦物識(shí)別、固體礦產(chǎn)和油氣勘探、環(huán)境保護(hù)及環(huán)境監(jiān)測、植被分類及檢測、月球和行星探測等方面[1-2]?，F(xiàn)有高光譜遙感礦物識(shí)別技術(shù)主要采用基于礦物波譜曲線的光譜匹配技術(shù)，如混合像元分解技術(shù)、光譜角技術(shù)和混合調(diào)制匹配濾波技術(shù)等[3-5]?；旌舷裨纸饧夹g(shù)是一種常用的高光譜礦物識(shí)別技術(shù)。由于成像光譜儀的硬件問題，高光譜遙感圖像普遍光譜分辨率高而空間分辨率低，容易導(dǎo)致影像像元出現(xiàn)同譜異物現(xiàn)象，因此必須研究解混方法予以克服。

目前常用的解混方法包括像元純凈指數(shù)、N-Findr、頂點(diǎn)成分分析、最小二乘端元提取算法等[6-9]。在前人工作的基礎(chǔ)上，Roberts等[10]提出了經(jīng)典的迭代光譜解混方法(MESMA)； Asner等[11]引入蒙特卡洛理論，將其與迭代解混理論有效結(jié)合起來，提出了AutoMUC方法； Song[12]介紹了一種基于貝葉斯決策的BSMA方法； 吳柯等[13]提出了基于神經(jīng)網(wǎng)絡(luò)的端元可變解混方法； 林紅磊等[14]提出基于單次散射反照率的礦物高光譜稀疏解混方法。上述高光譜解混方法雖然反演精度有了很大的提高，但是經(jīng)典的方法大都依賴于純凈像元，由于圖像的空間分辨率和地表目標(biāo)的復(fù)雜程度，純凈的像元是很難直接從高光譜圖像中獲取的； 此外，這些方法在解釋混合像元的成像機(jī)理方面也還不完善。

Xu等[15]利用線性光譜混合模型來解釋混合像元的成像機(jī)理，利用凈化像元的概念提出K-P-Means算法來進(jìn)行高光譜混合像元分解。該算法分2步迭代(豐度估計(jì)和端元優(yōu)化)，通過基于高斯混合模型的期望最大(expectation maximization，EM)估計(jì)算法直到最終的端元估計(jì)收斂。由于K-P-Means算法估計(jì)端元和豐度的能力較強(qiáng)，因此，本文將該算法應(yīng)用于礦物識(shí)別。該算法端元的優(yōu)化借助于計(jì)算出的地物類別標(biāo)簽，由于地物類別標(biāo)簽可能不正確，導(dǎo)致計(jì)算出的凈化像元存在多個(gè)類別，這種現(xiàn)象會(huì)導(dǎo)致端元優(yōu)化不理想。Fischler等[16]提出的隨機(jī)抽樣一致性算法(random sample consensus，RANSAC)，可以有效剔除因分類標(biāo)簽不正確導(dǎo)致的凈化像元中的異常值的影響。因此本文提出基于RANSAC的穩(wěn)健的K-P-Means算法(RANSAC based robust K-P-Means，RR-K-P-Means)。在K-P-Means算法的基礎(chǔ)上，利用RANSAC算法剔除異常影響，提高端元優(yōu)化的精度。最終將該算法用于美國內(nèi)華達(dá)州銅礦區(qū)的AVIRIS高光譜傳感器Cuprite數(shù)據(jù)的礦物識(shí)別，將估計(jì)出的光譜曲線與美國地質(zhì)調(diào)查局網(wǎng)站提供的礦物波譜曲線進(jìn)行匹配，確定礦物種類。

1 高光譜混合像元分解方法

1.1 基本原理

1.1.1 線性光譜混合模型

線性光譜混合模型是一種常用的高光譜遙感圖像表達(dá)方式[17]。假設(shè)高光譜圖像像元集X由端元矩陣A和豐度矩陣S以及獨(dú)立同分布的高斯噪聲N組成，即

X=SAT+N，

(1)

(2)

(3)

式中：si(i=1,2,…，m)為E×1的豐度向量，用來表達(dá)第j個(gè)端元光譜向量aj(j=1,2,…,E)在P×1維高光譜像元xi中的貢獻(xiàn)；m為圖像像元個(gè)數(shù)；E為端元個(gè)數(shù)；P為波段數(shù)；n為高斯噪聲向量。

1.1.2 K-P-Means模型

由于同一類的混合像元允許多種端元存在，為了更進(jìn)一步移除豐度值較低的端元的影響，使用豐度值最大的端元(而不是用整個(gè)xi)來估計(jì)aj是很合理的。Xu等[15]把移除豐度值較低端元貢獻(xiàn)后的xi稱為凈化像元。

K-P-Means模型是傳統(tǒng)K-Means模型的衍生，考慮到實(shí)際中豐度非負(fù)，因此算法可以表述為：

(4)

式中k為地物類別數(shù)。

相比K-Means考慮物理過程的整體效果，K-P-Means探究影響觀測的物理過程的源頭，其目標(biāo)函數(shù)可以表述為:

(5)

(6)

式中：l為類別；yi為第k類的凈化像元向量。

基于上面描述的模型，由于未知數(shù)遠(yuǎn)遠(yuǎn)多于觀測方程的個(gè)數(shù)，這是一個(gè)嚴(yán)重奇異的問題。EM算法對(duì)于奇異問題的參數(shù)求解具有較好的效果，該算法通過迭代的方法尋找統(tǒng)計(jì)模型的最大似然估計(jì)[18]。K-P-Means模型依據(jù)EM算法分2步迭代進(jìn)行豐度估計(jì)和端元優(yōu)化： 首先,通過非負(fù)最小二乘法進(jìn)行豐度的估計(jì); 然后,通過計(jì)算的豐度對(duì)端元進(jìn)行優(yōu)化。

1.1.3 豐度估計(jì)

給定{aj}，假設(shè)噪聲滿足高斯分布模型，則其噪聲的概率密度函數(shù)p(n)為：

(7)

(8)

估計(jì)豐度的目標(biāo)函數(shù)可以表達(dá)為：

(9)

由于方差可以通過均勻區(qū)域法獲取，并作為權(quán)重計(jì)算豐度，地物類別li通過sik最大確定，即

(10)

因此，式(9)中豐度的估計(jì)本質(zhì)上就是一個(gè)加權(quán)非負(fù)最小二乘(weight nonnegative least square， WNNLS)問題[19]。

1.1.4 端元優(yōu)化

(11)

1.2 RR-K-P-Means算法

從上述K-P-Means算法的原理很容易發(fā)現(xiàn)，端元的優(yōu)化完全依賴于由豐度確定的標(biāo)簽。類別標(biāo)簽錯(cuò)誤或者存在異常值都會(huì)導(dǎo)致端元的優(yōu)化精度降低。實(shí)踐過程中，從凈化像元值的直方圖可以明顯發(fā)現(xiàn)存在多個(gè)類別(圖1)。

圖1 凈化像元直方圖Fig.1 Histogram of purified pixels

凈化像元值存在多個(gè)類別的現(xiàn)象降低了利用K-P-Means算法進(jìn)行混合像元分解的精度。在計(jì)算機(jī)視覺以及其他很多研究領(lǐng)域，RANSAC算法對(duì)錯(cuò)誤率超過50%的數(shù)據(jù)仍然能夠得到理想的處理結(jié)果，是最有效的穩(wěn)健估計(jì)算法之一[20]。因此，本文利用RANSAC算法對(duì)K-P-Means算法混合像元分解過程中產(chǎn)生的凈化像元進(jìn)行優(yōu)化，剔除錯(cuò)誤標(biāo)簽的影響，實(shí)現(xiàn)對(duì)K-P-Means算法的優(yōu)化。

基于RANSAC算法的穩(wěn)健的K-P-Means算法基本過程可以描述為： ①利用式(1)—(10)原理獲取帶有錯(cuò)誤標(biāo)簽的凈化像元； ②尋找一個(gè)模型(一般為線性模型)適應(yīng)于假設(shè)的正確的凈化像元(初始采用隨機(jī)點(diǎn))，利用尋找的凈化像元估算該模型的參數(shù)(圖2)； ③用第二步得到的模型去測試所有的其他端元，若某個(gè)端元適用于估計(jì)的模型，則認(rèn)為它也是該類別的端元，其他端元滿足的條件設(shè)置為大于最小端元差值(設(shè)定閾值)，如果有足夠多的端元被歸類為假設(shè)的正確類別的端元，那么估計(jì)的模型就足夠合理，用所有假設(shè)的正確的端元去重新估計(jì)選取的模型，直到估計(jì)的正確的端元數(shù)量滿足給定條件(設(shè)置最低錯(cuò)誤率)為止； ④重復(fù)第二步和第三步過程來估計(jì)更加穩(wěn)健的參數(shù)，選擇更加合理的端元，將獲取的合理的端元帶入式(11)，獲取優(yōu)化后的端元。

圖2 基于RANSAC算法的端元優(yōu)化原理示意圖Fig.2 Schematic diagram of the principle of end-memberoptimization based on RANSAC

2 實(shí)驗(yàn)與結(jié)果

2.1 仿真數(shù)據(jù)實(shí)驗(yàn)

從美國地質(zhì)調(diào)查局公布的地物光譜庫中，隨機(jī)選擇4條地物光譜曲線，按照如下程序混合成64×64大小的圖像。利用這4個(gè)端元，首先把整個(gè)圖像分成8×8的同質(zhì)塊； 然后通過7×7窗口大小的空間低通濾波對(duì)這些同質(zhì)塊降級(jí)； 再通過添加零均值高斯噪聲對(duì)圖像再降級(jí)，從而獲取和真實(shí)情況更為接近的仿真數(shù)據(jù)。

豐度和端元的估計(jì)實(shí)驗(yàn)中，如果2個(gè)端元估計(jì)值的光譜角度距離(spectral angle distance，SAD)小于一個(gè)給定的數(shù)值τ或者預(yù)先規(guī)定迭代的次數(shù)iters，那么2步迭代就停止。本次實(shí)驗(yàn)中τ=0.01，iters=50。每次迭代中采用頂點(diǎn)成分分析的結(jié)果作為端元的初值[8]。利用RR-K-P-Means算法對(duì)端元進(jìn)行優(yōu)化時(shí)，模型采用線性模型，其他端元距離該直線的最小端元差值設(shè)置為5，最低錯(cuò)誤率設(shè)置為15%。

分別利用K-P-Means算法和RR-K-P-Means算法進(jìn)行豐度和端元的估計(jì)。用廣泛使用的SAD和光譜信息散度(spectral information divergence, SID)來評(píng)價(jià)估計(jì)端元和真實(shí)端元的一致性[21]。為了便于表示，將SID的數(shù)值統(tǒng)一擴(kuò)大10 000倍。利用結(jié)構(gòu)相似性(structural similarity，SSIM)和峰值信噪比(peak signal-to-noise ratio，PSNR)來衡量利用估計(jì)的豐度和端元重新獲得的圖像與原圖像的相似性[22-23]。

SAD和SID值越小，估計(jì)的光譜曲線與真實(shí)光譜曲線越相似。SSIM和PSNR越大，表明利用估計(jì)的豐度和端元重新獲得的圖像與原始圖像越相似。從表1中數(shù)據(jù)可以看出，利用RR-K-P-Means算法估計(jì)的端元和豐度SAD達(dá)到0.73，SID達(dá)到3.1，估計(jì)值和真實(shí)值的一致性分別提高8.8%和13.89%。

表1 優(yōu)化前后圖像各指標(biāo)對(duì)比Tab.1 Comparison of image parametersbefore and after optimization

利用估計(jì)的端元和豐度重新獲取的圖像和原始圖像相似性SSIM達(dá)到0.997，PSNR達(dá)到35.67，相似性分別提高10.17%和62.80%。RR-K-P-Means算法明顯優(yōu)于傳統(tǒng)的K-P-Means算法。從PSNR來看，RR-K-P-Means算法能更好地減弱噪聲對(duì)估計(jì)的影響。

RR-K-P-Means算法可以剔除個(gè)數(shù)較少的異常類別的影響，而且選擇的數(shù)據(jù)符合數(shù)據(jù)的分布形式(圖3)。盡管RR-K-P-Means算法不能在大范圍上改善分類結(jié)果，但是在細(xì)節(jié)上相比K-P-Means算法有了很大的改善。RR-K-P-Means算法估計(jì)的端元和豐度如圖4所示。圖4(a)—(d)為4種地物端元的估計(jì)值和真實(shí)值，兩者吻合程度比較高； 圖4(e)—(h)和(i)—(l)分別為4種地物的真實(shí)豐度和估計(jì)豐度，說明本文算法可以很好地計(jì)算豐度。較高精度的端元和豐度估計(jì)結(jié)果，為后續(xù)的礦物識(shí)別做好了準(zhǔn)備。

圖3 真實(shí)標(biāo)簽與優(yōu)化前后2種算法確定的地物標(biāo)簽

圖4 優(yōu)化后估計(jì)的端元、豐度與真實(shí)值對(duì)比

2.2 真實(shí)高光譜數(shù)據(jù)

將RR-K-P-Means算法用于真實(shí)高光譜遙感圖像數(shù)據(jù)的礦物提取。實(shí)驗(yàn)中采用美國內(nèi)華達(dá)州銅礦區(qū)的AVIRIS高光譜傳感器的Cuprite數(shù)據(jù)集(圖5)，數(shù)據(jù)獲取網(wǎng)站地址： http: //www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes。Cuprite數(shù)據(jù)集除了少數(shù)植被覆蓋外大部分為含指示礦物的蝕變帶，是進(jìn)行高光譜遙感地質(zhì)研究的典型樣區(qū)。該影像數(shù)據(jù)為1997年6月19日機(jī)載AVIRIS獲得，包括224個(gè)波段，地面瞬時(shí)視場約為20 m。剔除水汽影響嚴(yán)重的波段，共選取159個(gè)波段進(jìn)行實(shí)驗(yàn)。

圖5 內(nèi)華達(dá)州銅礦區(qū)的Cuprite數(shù)據(jù)集Fig.5 Cuprite data sets of Nevadacopper mining area

利用RR-K-P-Means算法對(duì)Cuprite數(shù)據(jù)集進(jìn)行端元和豐度的估計(jì)，通過計(jì)算估計(jì)出的端元和美國地質(zhì)調(diào)查局礦物波譜庫的波譜曲線的相關(guān)性確定端元所屬的礦物類別。由于實(shí)際地物比較復(fù)雜，因此實(shí)驗(yàn)中的端元數(shù)設(shè)置為30，多于本地礦物的種類。選擇相關(guān)系數(shù)大于0.75的端元進(jìn)行研究。

通過相關(guān)系數(shù)的匹配得出8個(gè)類別的主要礦物，分別為： 黃鉀鐵礬、水銨長石、玉髓、葉蠟石、綠脫石、綠泥石、蒙脫石和白云母。研究區(qū)主要礦物綠脫石的端元估計(jì)結(jié)果如圖6所示，波長在400～600 nm的端元估計(jì)結(jié)果和真實(shí)值匹配度一致，其他部分基本相同。

圖6 綠脫石的端元估計(jì)值與真實(shí)值Fig.6 Estimated and true endmember values of nontronite

研究區(qū)主要礦物黃鉀鐵礬的端元估計(jì)結(jié)果如圖7所示，估計(jì)端元波譜曲線與真實(shí)值基本相同。在實(shí)際應(yīng)用中，利用RR-K-P-Means算法估計(jì)的端元可以有效提取礦物的波譜曲線。

圖7 黃鉀鐵礬的端元估計(jì)值與真實(shí)值Fig.7 Estimated and true endmember values of jarosite

8種主要礦物的豐度圖見圖8。由圖8可見，礦物聚集性比較好，與已知樣本分布相比較，利用RR-K-P-Means算法獲取的8種主要礦物基本與實(shí)際情況一致。

圖8-2 識(shí)別出的8種礦物的豐度圖

3 結(jié)論

針對(duì)K-P-Means算法中端元優(yōu)化受異常值影響的問題，提出基于RANSAC的穩(wěn)健的K-P-Means算法(RR-K-P-Means算法)，并將該算法用于礦物識(shí)別，通過仿真和真實(shí)數(shù)據(jù)實(shí)驗(yàn)，都取得了理想的效果。得出如下結(jié)論：

1)仿真實(shí)驗(yàn)證明，利用RR-K-P-Means算法估計(jì)的端元和豐度與真實(shí)值的一致性提高明顯。利用估計(jì)的端元和豐度重新獲取的圖像和原始圖像相似性提高較多。從PSNR來看，RR-K-P-Means算法能更好地減弱噪聲對(duì)估計(jì)的影響。本文算法可以剔除個(gè)數(shù)較少的異常類別的影響，而且選擇的數(shù)據(jù)符合數(shù)據(jù)的分布形式。

2)通過真實(shí)數(shù)據(jù)驗(yàn)證，利用RR-K-P-Means算法可以較好地識(shí)別研究區(qū)主要礦物，和美國地質(zhì)調(diào)查局地物波譜庫提供的標(biāo)準(zhǔn)地物波譜相比具有較高的一致性。在高光譜遙感圖像礦物識(shí)別中，RR-K-P-Means算法可以有效提取多種地物光譜，礦物識(shí)別效果較好。

本文端元優(yōu)化采用算數(shù)平均值，豐度較小的標(biāo)簽未參與計(jì)算，未來可以將加權(quán)平均的方法應(yīng)用到本研究中。除此之外，RR-K-P-Means不僅在高光譜遙感圖像礦物識(shí)別中可以得到很好的應(yīng)用，同時(shí)可以推廣到分類、去噪、超分辨率重建等研究中。