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

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

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

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

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

本文以ASTER數據為數據源，利用熱紅外波段數據計算出其對應的礦物學指數作為特征，通過隨機森林分類算法，對巖礦信息提取分類并繪制礦藏分布圖。

1 研究方法

本文采用的遙感圖像為ASTER遙感圖像，表1中列出了各個子系統(tǒng)對應波段的相關參數。

表1 ASTER圖像的波段

巖礦信息提取通過機器學習方法實現，首先提取特征，然后通過隨機森林方法實現巖礦信息分類，按以下四步實現，即輻射定標、歸一化處理、特征提取和分類。

1.1 輻射定標

將ASTER圖像上的亮度灰度值轉換為絕對的輻射亮度

L=(d-1)×c

(1)

式中：d為從存儲ASTER圖像的文件中讀取的圖像亮度灰度值；c為每個波段對應的增益值；L為表面輻射值。

1.2 歸一化處理

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

(2)

式中：Li是由式(1)給出的頻段i中的傳感器數據的ASTER輻射值，i為波段號。λi是每個頻段i對應的中心波長(μm)；ε13是頻段13中的假定發(fā)射率為1.0；T是第13頻段被歸一化的固定溫度(K)，通過經驗公式的計算將其值定為300；c1=3.742×108J·m,c2=1.436 9×104K通過這種歸一化輻射的過程，能夠將碳酸鹽指數對表面溫度條件的重度依賴性顯著改善。

1.3 特征提取

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

(3)

(4)

(5)

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

特征提取步驟：

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

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

(3)用歸一化處理之后的數據進行指數計算分別求出每一個像素對應的QI、CI、MI的值。

每個已知數據點選取其周圍7×7鄰域的像素計算，每個像素為一個樣本，將其計算出的QI、CI、MI的值作為特征輸出。

1.4 隨機森林分類算法

1.4.1 分類決策樹

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

假設D表示樣本集個數，屬性a有v個可能的取值(離散或連續(xù))。進行最有劃分屬性時，比如先找到了屬性a，對a進行評價，接下來對其他屬性重復a的過程，分別得到一個評分，選擇評分最高的那個，即信息增益最大的作為最有劃分屬性。

1.4.2 模型建立

將樣本集D隨機分成訓練集S與測試集T，特征維數為F。決策樹的數量為t，每個節(jié)點使用到的特征數量f。終止條件：所有的節(jié)點都訓練過了或者被標記為葉子節(jié)點。

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

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

(3)重復(1)(2)直到所有節(jié)點都訓練過了或者被標記為葉子節(jié)點。

(4)重復(1)(2)(3)直到t棵樹都被訓練過。

1.4.3 預測分類

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

(1)從當前樹的根節(jié)點開始，根據當前節(jié)點的閾值h，判斷是進入左節(jié)點(

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

2 試驗驗證

2.1 試驗數據

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

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

2.2 實驗方法

2.2.1 交叉驗證測試

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

巖礦分類評價指標正確率為

(6)

式中：Rj表示用于交叉驗證的子集j分類正確的樣本數，Aj表示子集j的總樣本數。

測試方法是將以上是一張圖像已知數據點提取整合，打亂順序，然后分組進行交叉驗證，將所有ACC加和平均作為最終評價結果。圖1顯示K-Fold交叉驗證結果，顯示測試集編號分別為1、2、3、4的準確率。

圖1 帕米爾高原圖像K-Fold交叉驗證結果

結果顯示，此種方式對于訓練數據點與測試數據點在同一張圖像上的樣本分類能夠得到較高的正確率，但其并不能夠反映該分類器能夠以上表得出的正確率去分類一張沒有進行訓練的圖像。

2.2.2 新圖像測試

為了測試所訓練的模型對不同圖像的分類性能、反映不同圖像間的分類準確性，現采用第1項測試，測試圖像同表1，第2種測試方法根據待分類中數據點存在的分布選取6張圖像(編號為3、4、6、7、9、10，這6張圖像中包含所有的數據點)的數據點進行訓練，剩下4張圖像(1、2、5、8)分別用來測試已知數據點的正確率，正確率計算公式同測試一所述，每張圖像每個類別對應分類準確率如表2所示。

表2 帕米爾高原圖像分類結果

結果顯示，對于訓練數據點與測試數據點取自完全不同的圖像的情況下，平均正確率能夠達到0.849 7，對于變質硅酸鹽的分類效果與其他類別的分類效果相比較弱，但平均正確率也能達到0.5以上，對硫酸鹽和硅質的分類效果較好，平均正確率都在0.9以上。

2.2.3 整張圖像分類測試

將所有圖片中的已知數據點作為訓練集，訓練隨機森林的分類器，對整張圖像所有像素進行測試，并將圖像中的不同類別用不同顏色標注出來，與文獻[10]的預測方法進行對比。

圖2，圖3分別是對編號為8和10的圖像進行全圖預測的結果，圖中白色的標記點是已知數據點，本文預測結果與文獻[10]的顏色不是完全一致，可通過白色的標記點找到對應關系。通過同文獻[10]結果的對比，可以看出本文預測得到的巖礦分類圖與文獻[10]的一致性較高。由于圖像的質量，如云層的干擾，積雪的干擾等都對數據產生較大的影響，所以對全圖預測的結果是有偏差的。

圖3 編號為10圖像的預測類別圖

2.2.4 指標分析

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

圖4 文獻[8]分類精度

綜合同以上文獻的對比，本文對已知數據點分類的準確率相對較高，能夠識別硅質、硫酸鹽、石英、超鎂鐵質以及石英鎂鐵質這5個類別，平均分類正確率能夠達到0.849 7以上。對所有研究圖像進行了類別預測，主觀評價結果也較好。

3 結論

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

從測試結果可以看出，本文所述方法能夠通過遙感圖像識別出5種巖礦類型，對于訓練數據點與測試數據點取自完全不同的圖像的情況下，平均正確率相對其他文獻來說較好，對整張遙感圖像繪制礦藏分布圖，主觀判斷其效果較好。

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