黃宇飛 李智慧 寧慧 胡震岳 曹海翊

(1 北京空間飛行器總體設(shè)計(jì)部，北京 100094)(2 哈爾濱工程大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，哈爾濱 150001)

遙感影像在不同波段范圍內(nèi)反映了地物輻射、反射和發(fā)射電磁波的特性等信息。各種礦產(chǎn)資源的形成與產(chǎn)出都和一定的地質(zhì)體與地質(zhì)現(xiàn)象的空間展布等密切聯(lián)系，不同的礦物、巖石或蝕變礦化等也可通過(guò)遙感影像上不同的波譜曲線反映出來(lái)。因此通過(guò)遙感波譜信息及空間信息的提取，可有針對(duì)性地獲得與礦化有關(guān)的地質(zhì)體、地質(zhì)現(xiàn)象的類別和分布狀況[1]。遙感工作者對(duì)不同礦化蝕變類型采取針對(duì)性的異常提取方法在在不同地區(qū)都開(kāi)展了很好的實(shí)踐。

近幾年國(guó)內(nèi)外遙感找礦的研究主要分為主成分分析方法、復(fù)合處理方法以及基于高級(jí)星載熱發(fā)射反射輻射計(jì)(Advanced Spaceborne Thermal Emission and Reflection Radiometer，ASTER)數(shù)據(jù)的方法。文獻(xiàn)[2]提出了“TM掩膜+主成分變換+分類”識(shí)別提取礦產(chǎn)弱信息的技術(shù)方法，在遼寧二道溝金礦提取三價(jià)鐵蝕變火山碎屑巖信息，在河北華北地臺(tái)北緣提取含金鉀化帶信息，均取得了較好的效果。文獻(xiàn)[3]利用“多元數(shù)據(jù)分析+比值+主成份變換+掩膜+分類(分割)”的方法在新疆、內(nèi)蒙古及江西、云南成功的提取了金礦化蝕變信息。文獻(xiàn)[4]利用TM、SPOT數(shù)據(jù)數(shù)字高程模型及地球化學(xué)信息對(duì)愛(ài)琴海盆地Quatemary島弧中的低溫?zé)嵋盒徒鸬V進(jìn)行研究。文獻(xiàn)[5]利用TM圖像743波段合成圖像，根據(jù)巖石變形變質(zhì)特征及遙感影像特征，對(duì)喀喇昆侖區(qū)域進(jìn)行地層劃分研究并成功地劃分出了該區(qū)的地層系統(tǒng)。文獻(xiàn)[6]利用機(jī)載可見(jiàn)紅外成像光譜儀(AVIRIS)和ASTER數(shù)據(jù)對(duì)成層火山的熱液蝕變進(jìn)行了研究。文獻(xiàn)[7]采用ASTER的近紅外，短波紅外波段數(shù)據(jù)以及變差函數(shù)紋理等特征用于巖性分類。文獻(xiàn)[8]利用ASTER波段數(shù)據(jù)、溫度數(shù)據(jù)、高程數(shù)據(jù)等作為特征，用支持向量機(jī)進(jìn)行預(yù)測(cè)。研究表明ASTER數(shù)據(jù)其波段設(shè)置比其他多光譜數(shù)據(jù)具有更好的波譜連續(xù)性，因此ASTER數(shù)據(jù)在巖性信息提取和礦物識(shí)別方面能達(dá)到更高的準(zhǔn)確率。

ASTER傳感器是搭載在“土”衛(wèi)星(Terra)上的星載熱發(fā)射反射輻射計(jì)，于1999年12月18日發(fā)射升空，其數(shù)據(jù)包括近紅外、短波紅外、熱紅外3個(gè)光學(xué)系統(tǒng)，共計(jì)14個(gè)頻段，空間分辨率分別為15 m、30 m、90 m。ASTER遙感數(shù)據(jù)以其更多的頻段，提取蝕變礦物異常類型更為豐富，在礦化信息識(shí)別取得良好的應(yīng)用效果。

本文以ASTER數(shù)據(jù)為數(shù)據(jù)源，利用熱紅外波段數(shù)據(jù)計(jì)算出其對(duì)應(yīng)的礦物學(xué)指數(shù)作為特征，通過(guò)隨機(jī)森林分類算法，對(duì)巖礦信息提取分類并繪制礦藏分布圖。

1 研究方法

本文采用的遙感圖像為ASTER遙感圖像，表1中列出了各個(gè)子系統(tǒng)對(duì)應(yīng)波段的相關(guān)參數(shù)。

表1 ASTER圖像的波段

巖礦信息提取通過(guò)機(jī)器學(xué)習(xí)方法實(shí)現(xiàn)，首先提取特征，然后通過(guò)隨機(jī)森林方法實(shí)現(xiàn)巖礦信息分類，按以下四步實(shí)現(xiàn)，即輻射定標(biāo)、歸一化處理、特征提取和分類。

1.1 輻射定標(biāo)

將ASTER圖像上的亮度灰度值轉(zhuǎn)換為絕對(duì)的輻射亮度

L=(d-1)×c

(1)

式中：d為從存儲(chǔ)ASTER圖像的文件中讀取的圖像亮度灰度值；c為每個(gè)波段對(duì)應(yīng)的增益值；L為表面輻射值。

1.2 歸一化處理

TIR頻段基于普朗克定律，將所有像素的第13頻段的亮度溫度轉(zhuǎn)換為固定溫度并將其他4個(gè)波段的數(shù)據(jù)轉(zhuǎn)換，得到傳感器數(shù)據(jù)的歸一化輻射值

(2)

式中：Li是由式(1)給出的頻段i中的傳感器數(shù)據(jù)的ASTER輻射值，i為波段號(hào)。λi是每個(gè)頻段i對(duì)應(yīng)的中心波長(zhǎng)(μm)；ε13是頻段13中的假定發(fā)射率為1.0；T是第13頻段被歸一化的固定溫度(K)，通過(guò)經(jīng)驗(yàn)公式的計(jì)算將其值定為300；c1=3.742×108J·m,c2=1.436 9×104K通過(guò)這種歸一化輻射的過(guò)程，能夠?qū)⑻妓猁}指數(shù)對(duì)表面溫度條件的重度依賴性顯著改善。

1.3 特征提取

通過(guò)不同類別的巖石在對(duì)不同波長(zhǎng)的吸收程度分析，分別得出石英指數(shù)(QI)、碳酸鹽指數(shù)(CI)、鎂鐵質(zhì)指數(shù)(MI)的計(jì)算公式為

(3)

(4)

(5)

式中：Ni為式(1)求得的歸一化后的輻射值。

特征提取步驟：

(1)讀取ASTER圖像的TIR的波段數(shù)據(jù)，共五維，利用式(1)將其圖像亮度灰度值轉(zhuǎn)換成輻射值；

(2)將步驟1中的輻射值利用式(2)進(jìn)行歸一化處理；

(3)用歸一化處理之后的數(shù)據(jù)進(jìn)行指數(shù)計(jì)算分別求出每一個(gè)像素對(duì)應(yīng)的QI、CI、MI的值。

每個(gè)已知數(shù)據(jù)點(diǎn)選取其周?chē)?×7鄰域的像素計(jì)算，每個(gè)像素為一個(gè)樣本，將其計(jì)算出的QI、CI、MI的值作為特征輸出。

1.4 隨機(jī)森林分類算法

1.4.1 分類決策樹(shù)

隨機(jī)森林是一種有監(jiān)督學(xué)習(xí)算法，它創(chuàng)建了一個(gè)森林，并使它擁有某種方式隨機(jī)性[9]。所構(gòu)建的“森林”是決策樹(shù)的集成，分類決策樹(shù)算法通過(guò)熵作為評(píng)價(jià)指標(biāo)對(duì)分類效果進(jìn)行評(píng)估。熵值越大，表明數(shù)據(jù)的純度越低。當(dāng)熵等于0，表明樣本數(shù)據(jù)都是同一個(gè)類別。

假設(shè)D表示樣本集個(gè)數(shù)，屬性a有v個(gè)可能的取值(離散或連續(xù))。進(jìn)行最有劃分屬性時(shí)，比如先找到了屬性a，對(duì)a進(jìn)行評(píng)價(jià)，接下來(lái)對(duì)其他屬性重復(fù)a的過(guò)程，分別得到一個(gè)評(píng)分，選擇評(píng)分最高的那個(gè)，即信息增益最大的作為最有劃分屬性。

1.4.2 模型建立

將樣本集D隨機(jī)分成訓(xùn)練集S與測(cè)試集T，特征維數(shù)為F。決策樹(shù)的數(shù)量為t，每個(gè)節(jié)點(diǎn)使用到的特征數(shù)量f。終止條件：所有的節(jié)點(diǎn)都訓(xùn)練過(guò)了或者被標(biāo)記為葉子節(jié)點(diǎn)。

(1)從S中有放回的抽取大小和S一樣的訓(xùn)練集S(xm)，作為根節(jié)點(diǎn)的樣本，從根節(jié)點(diǎn)開(kāi)始訓(xùn)練。

(2)如果當(dāng)前節(jié)點(diǎn)上達(dá)到終止條件，則設(shè)置當(dāng)前節(jié)點(diǎn)為葉子節(jié)點(diǎn)，該葉子節(jié)點(diǎn)的預(yù)測(cè)輸出為當(dāng)前節(jié)點(diǎn)樣本集合中數(shù)量最多的那一類c(n)，概率p為c(n)占當(dāng)前樣本集的比例，繼續(xù)訓(xùn)練其他節(jié)點(diǎn)；如果當(dāng)前節(jié)點(diǎn)沒(méi)有達(dá)到終止條件，則從F維特征中無(wú)放回的隨機(jī)選取f維特征。利用這f維特征，尋找分類效果最好的一維特征k及其閾值h，當(dāng)前節(jié)點(diǎn)上樣本第k維特征小于h的樣本被劃分到左節(jié)點(diǎn)，其余的被劃分到右節(jié)點(diǎn)，繼續(xù)訓(xùn)練其他節(jié)點(diǎn)。

(3)重復(fù)(1)(2)直到所有節(jié)點(diǎn)都訓(xùn)練過(guò)了或者被標(biāo)記為葉子節(jié)點(diǎn)。

(4)重復(fù)(1)(2)(3)直到t棵樹(shù)都被訓(xùn)練過(guò)。

1.4.3 預(yù)測(cè)分類

在得到隨機(jī)森林訓(xùn)練的分類器模型之后，輸入未經(jīng)訓(xùn)練得巖礦指數(shù)特征，讓森林中的每一棵決策樹(shù)分別進(jìn)行分類：

(1)從當(dāng)前樹(shù)的根節(jié)點(diǎn)開(kāi)始，根據(jù)當(dāng)前節(jié)點(diǎn)的閾值h，判斷是進(jìn)入左節(jié)點(diǎn)(

(2)重復(fù)執(zhí)行(1)直到所有t棵樹(shù)都輸出了預(yù)測(cè)值，輸出為所有樹(shù)中預(yù)測(cè)概率總和最大的那一個(gè)類。

2 試驗(yàn)驗(yàn)證

2.1 試驗(yàn)數(shù)據(jù)

首先將文獻(xiàn)[10]公布的巖礦標(biāo)記圖像通過(guò)谷歌地球軟件導(dǎo)出為GEOTIFF格式的圖片，利用ENVI軟件讀取其中已知真實(shí)數(shù)據(jù)點(diǎn)，包括其經(jīng)緯度以及巖礦類別。待研究的樣本是帕米爾高原地區(qū)的，已知的巖礦數(shù)據(jù)點(diǎn)共50個(gè)，分布在編號(hào)為1～10的十張圖像中，每個(gè)數(shù)據(jù)點(diǎn)為某點(diǎn)附近的一小塊區(qū)域。

所判別的巖礦礦物類型共計(jì)五種，分別為硅質(zhì)、硫酸鹽、超鎂鐵質(zhì)、石英以及變質(zhì)硅酸鹽。采用3種測(cè)試方法：①交叉驗(yàn)證測(cè)試，是對(duì)所有已知點(diǎn)數(shù)據(jù)進(jìn)行交叉驗(yàn)證測(cè)試；②新圖像測(cè)試，即測(cè)試集與訓(xùn)練集選取完全不同的圖像進(jìn)行測(cè)試；③整張圖像類別測(cè)試，將所有圖像的已知數(shù)據(jù)點(diǎn)作為訓(xùn)練集，然后對(duì)整張圖像的類別進(jìn)行預(yù)測(cè)。

2.2 實(shí)驗(yàn)方法

2.2.1 交叉驗(yàn)證測(cè)試

鑒于已知數(shù)據(jù)點(diǎn)個(gè)數(shù)有限，所以第一項(xiàng)測(cè)試方法采用K-Fold交叉驗(yàn)證，K值取4，將原始數(shù)據(jù)分為4個(gè)子集，每個(gè)子集分別作為一次測(cè)試集，其他3個(gè)子集作為訓(xùn)練集，進(jìn)行循環(huán)交叉驗(yàn)證，并計(jì)算得到交叉驗(yàn)證正確率。

巖礦分類評(píng)價(jià)指標(biāo)正確率為

(6)

式中：Rj表示用于交叉驗(yàn)證的子集j分類正確的樣本數(shù)，Aj表示子集j的總樣本數(shù)。

測(cè)試方法是將以上是一張圖像已知數(shù)據(jù)點(diǎn)提取整合，打亂順序，然后分組進(jìn)行交叉驗(yàn)證，將所有ACC加和平均作為最終評(píng)價(jià)結(jié)果。圖1顯示K-Fold交叉驗(yàn)證結(jié)果，顯示測(cè)試集編號(hào)分別為1、2、3、4的準(zhǔn)確率。

圖1 帕米爾高原圖像K-Fold交叉驗(yàn)證結(jié)果

結(jié)果顯示，此種方式對(duì)于訓(xùn)練數(shù)據(jù)點(diǎn)與測(cè)試數(shù)據(jù)點(diǎn)在同一張圖像上的樣本分類能夠得到較高的正確率，但其并不能夠反映該分類器能夠以上表得出的正確率去分類一張沒(méi)有進(jìn)行訓(xùn)練的圖像。

2.2.2 新圖像測(cè)試

為了測(cè)試所訓(xùn)練的模型對(duì)不同圖像的分類性能、反映不同圖像間的分類準(zhǔn)確性，現(xiàn)采用第1項(xiàng)測(cè)試，測(cè)試圖像同表1，第2種測(cè)試方法根據(jù)待分類中數(shù)據(jù)點(diǎn)存在的分布選取6張圖像(編號(hào)為3、4、6、7、9、10，這6張圖像中包含所有的數(shù)據(jù)點(diǎn))的數(shù)據(jù)點(diǎn)進(jìn)行訓(xùn)練，剩下4張圖像(1、2、5、8)分別用來(lái)測(cè)試已知數(shù)據(jù)點(diǎn)的正確率，正確率計(jì)算公式同測(cè)試一所述，每張圖像每個(gè)類別對(duì)應(yīng)分類準(zhǔn)確率如表2所示。

表2 帕米爾高原圖像分類結(jié)果

結(jié)果顯示，對(duì)于訓(xùn)練數(shù)據(jù)點(diǎn)與測(cè)試數(shù)據(jù)點(diǎn)取自完全不同的圖像的情況下，平均正確率能夠達(dá)到0.849 7，對(duì)于變質(zhì)硅酸鹽的分類效果與其他類別的分類效果相比較弱，但平均正確率也能達(dá)到0.5以上，對(duì)硫酸鹽和硅質(zhì)的分類效果較好，平均正確率都在0.9以上。

2.2.3 整張圖像分類測(cè)試

將所有圖片中的已知數(shù)據(jù)點(diǎn)作為訓(xùn)練集，訓(xùn)練隨機(jī)森林的分類器，對(duì)整張圖像所有像素進(jìn)行測(cè)試，并將圖像中的不同類別用不同顏色標(biāo)注出來(lái)，與文獻(xiàn)[10]的預(yù)測(cè)方法進(jìn)行對(duì)比。

圖2，圖3分別是對(duì)編號(hào)為8和10的圖像進(jìn)行全圖預(yù)測(cè)的結(jié)果，圖中白色的標(biāo)記點(diǎn)是已知數(shù)據(jù)點(diǎn)，本文預(yù)測(cè)結(jié)果與文獻(xiàn)[10]的顏色不是完全一致，可通過(guò)白色的標(biāo)記點(diǎn)找到對(duì)應(yīng)關(guān)系。通過(guò)同文獻(xiàn)[10]結(jié)果的對(duì)比，可以看出本文預(yù)測(cè)得到的巖礦分類圖與文獻(xiàn)[10]的一致性較高。由于圖像的質(zhì)量，如云層的干擾，積雪的干擾等都對(duì)數(shù)據(jù)產(chǎn)生較大的影響，所以對(duì)全圖預(yù)測(cè)的結(jié)果是有偏差的。

圖3 編號(hào)為10圖像的預(yù)測(cè)類別圖

2.2.4 指標(biāo)分析

在利用ASTER數(shù)據(jù)進(jìn)行巖礦分類的文獻(xiàn)中，文獻(xiàn)[8]利用ASTER波段數(shù)據(jù)、溫度數(shù)據(jù)、高程數(shù)據(jù)等作為特征，用支持向量機(jī)進(jìn)行預(yù)測(cè)，共預(yù)測(cè)9類，其中包括五種巖礦類型，如圖4所示。分類精度能夠達(dá)到0.792 8，對(duì)超鎂鐵質(zhì)的分類精度為0.672 5，但本方法對(duì)超鎂鐵質(zhì)的分類精度能夠達(dá)到0.907 8。

圖4 文獻(xiàn)[8]分類精度

綜合同以上文獻(xiàn)的對(duì)比，本文對(duì)已知數(shù)據(jù)點(diǎn)分類的準(zhǔn)確率相對(duì)較高，能夠識(shí)別硅質(zhì)、硫酸鹽、石英、超鎂鐵質(zhì)以及石英鎂鐵質(zhì)這5個(gè)類別，平均分類正確率能夠達(dá)到0.849 7以上。對(duì)所有研究圖像進(jìn)行了類別預(yù)測(cè)，主觀評(píng)價(jià)結(jié)果也較好。

3 結(jié)論

本文以紅外多光譜數(shù)據(jù)為數(shù)據(jù)源，選取帕米爾高原地區(qū)的巖礦數(shù)據(jù)作為研究對(duì)象，計(jì)算石英指數(shù)，碳酸鹽指數(shù)，鎂鐵質(zhì)指數(shù)為特征值，采用隨機(jī)森林的分類方法訓(xùn)練出分類模型，再利用該分類模型去判斷該地區(qū)的巖礦種類，根據(jù)判斷出的結(jié)果得出已知數(shù)據(jù)點(diǎn)的正確率。

從測(cè)試結(jié)果可以看出，本文所述方法能夠通過(guò)遙感圖像識(shí)別出5種巖礦類型，對(duì)于訓(xùn)練數(shù)據(jù)點(diǎn)與測(cè)試數(shù)據(jù)點(diǎn)取自完全不同的圖像的情況下，平均正確率相對(duì)其他文獻(xiàn)來(lái)說(shuō)較好，對(duì)整張遙感圖像繪制礦藏分布圖，主觀判斷其效果較好。

本文主要用到ASTER遙感圖像數(shù)據(jù)，按照ASTER項(xiàng)目的計(jì)劃，其數(shù)據(jù)應(yīng)用于全球變化研究中，在科研工作中起到了很好的促進(jìn)作用。ASTER圖像質(zhì)量較高，其近紅外和短波紅外波段平均信噪比分別為186.73，196.33。但是在下載ASTER數(shù)據(jù)時(shí)發(fā)現(xiàn)，在2008年后的ASTER-L1T數(shù)據(jù)缺少短波紅外的數(shù)據(jù)，雖然本項(xiàng)目的研究主要用到的是熱紅外的5個(gè)頻段，但對(duì)數(shù)據(jù)點(diǎn)進(jìn)行篩選的過(guò)程還是需要借助于短波紅外的頻段數(shù)據(jù)的，除此之外ASTER的圖像時(shí)間分布不均，很多月份的數(shù)據(jù)不存在或者文件中數(shù)據(jù)丟失，這種狀況對(duì)圖像分析存在影響。對(duì)于巖礦信息提取，采用的主要是熱紅外波段，波長(zhǎng)范圍在8.125～11.650 μm之間(共計(jì)5個(gè)頻段)，由于其熱紅外波段數(shù)較多，對(duì)巖礦信息提取有很大的幫助，但熱紅外波段的空間分辨率較低，應(yīng)該適當(dāng)提高其空間分辨率。