基于GF-2影像和隨機森林算法的花崗偉晶巖提取

杜曉川, 婁德波, 徐林剛, 范瑩琳, 張琳, 李婉悅

(1.中國地質大學(北京)地球科學與資源學院,北京 100083; 2.中國地質科學院礦產(chǎn)資源研究所,自然資源部成礦作用與資源評價重點實驗室,北京 100037)

0 引言

隨著近年來新能源、新興材料的發(fā)展,鋰等關鍵金屬越來越受到重視[1],金屬鋰由于其特殊的物理化學性質,被廣泛應用到諸如電池、玻璃、陶瓷等各個工業(yè)領域,具有極高的戰(zhàn)略價值,其不僅是高科技產(chǎn)業(yè)必需材料,也逐漸成為日常生活中廣泛使用的常規(guī)材料。當前我國國內鋰業(yè)發(fā)展高度依賴礦物型鋰礦特別是花崗偉晶巖型鋰礦的供應[2],此類型鋰礦床通常位于大型花崗巖侵入體邊緣,受地質構造控制明顯,通常以偉晶巖脈集群的形式產(chǎn)出,總體呈帶狀分布,單條偉晶巖脈形態(tài)多為脈狀或板狀[3],因此基于高分辨率影像的遙感技術是尋找和發(fā)現(xiàn)偉晶巖脈的有效方法[4]。

遙感領域相關技術不斷發(fā)展成熟,使得遙感技術在地質信息提取的研究中得到了廣泛應用,有許多利用遙感技術進行含鋰偉晶巖探測的實例。如代晶晶等[5-6]在川西甲基卡地區(qū)進行了鋰輝石等偉晶巖稀有金屬礦床相關礦石礦物的波譜測試研究和偉晶巖遙感地質填圖; Cardoso-Fernandes等[7]利用多源遙感影像數(shù)據(jù),采用RGB彩色合成、波段比值分析、主成分分析等方法識別含鋰偉晶巖的相關異常,完成蝕變填圖; 姚佛軍等[4]通過增強ASTER遙感圖像中的差異來增強地質體的巖性差異,并突出含鋰偉晶巖的遙感特征; 王海宇[8]以WorldView-3為數(shù)據(jù)源,并基于全卷積網(wǎng)絡(fully convolutional networks, FCN)網(wǎng)絡和U-net網(wǎng)絡2種深度學習語義分割模型識別偉晶巖脈信息。前人通過遙感影像數(shù)據(jù)在含鋰偉晶巖提取方面開展了大量的理論研究和實踐,但大多數(shù)研究者側重于使用低空間分辨影像數(shù)據(jù),提取方法也較為傳統(tǒng),盡管已有基于高空間分辨率影像數(shù)據(jù)與深度學習網(wǎng)絡模型提取偉晶巖巖脈的先例,但該方法耗時較長,效率較低,且對樣本數(shù)量有一定要求。因此本文在眾多機器學習算法中選用隨機森林算法,其算法具有如下優(yōu)點: ①數(shù)據(jù)兼容性強,數(shù)據(jù)集無須規(guī)范化,且可以直接處理高維數(shù)據(jù); ②抗噪聲能力強,受異常值影響較小; ③魯棒性強,不容易發(fā)生過擬合現(xiàn)象; ④計算效率高、速度快、泛化能力強[9-11],且該方法結合高空間分辨率影像數(shù)據(jù)在提取復雜地物的效果較好。

據(jù)此,以青海省海西蒙古族藏族自治州天峻縣扎卡東南部區(qū)域為研究區(qū),提取預處理后的GF-2高空間分辨率影像數(shù)據(jù)中各類型地物的光譜特征、紋理特征、指數(shù)特征、地形特征及邊緣特征等25個特征變量構建特征子集,對子集中的特征變量進行特征重要性分析,依據(jù)特征重要性得分進行特征優(yōu)選,確定最優(yōu)的提取花崗偉晶巖的特征組合,并進行隨機森林分類,對分類結果進行精度評定,為后續(xù)該地區(qū)的含鋰花崗偉晶巖找礦勘查工作提供依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

研究區(qū)位于青海省海西蒙古族藏族自治州天峻縣扎卡東南部,地理坐標為E99°27′16″～99°30′32″, N36°53′36″～36°54′43″(圖1)。天峻縣地處青海省東北部,青海湖西側、祁連山南麓,柴達木盆地東部。境內高山縱橫,山脈呈東南西北走向,地形以山地為主,高山、中低山、山谷和山間盆地相間分布。平均海拔在4 000 m以上,地處高原寒帶氣候,年平均氣溫-1.5 ℃,年降水量360 mm,干旱少雨,多風沙。

圖1 研究區(qū)遙感影像

1.2 數(shù)據(jù)源及預處理

GF-2衛(wèi)星于2014年8月19日發(fā)射,是我國自主研制的首顆空間分辨率優(yōu)于1 m的民用光學遙感衛(wèi)星。搭載有2臺高分辨率的1 m全色與4 m多光譜相機,傳感器的波段參數(shù)如表1所示。對大氣校正后的多光譜影像以及定標后的全色影像做正射校正處理及融合處理,獲得空間分辨率為1 m的彩色合成影像,再通過研究區(qū)域的矢量范圍對影像數(shù)據(jù)進行裁剪,得到研究區(qū)遙感影像數(shù)據(jù)。

表1 GF-2全色多光譜相機波段參數(shù)

2 研究方法

2.1 特征提取及樣本選擇

花崗偉晶巖(圖2(a))在分布位置及產(chǎn)出形態(tài)上受構造控制作用明顯,研究區(qū)內斷裂沿NW-SE向分布,因此花崗偉晶巖大多沿NW-SE向展布,呈條狀、脈狀產(chǎn)出; 在影像色調上,其多呈灰白色或略帶土黃色; 由于研究區(qū)內海拔較高、風沙較多、風化能力較強,且花崗偉晶巖晶體較大,具有較強的抗風化性,因而在高度上要突出于周圍地物,并且其表面紋理粗糙,邊界較明顯。研究區(qū)內其他地物類型主要分為道路(圖2(b))、干涸河流(圖2(c))、冰雪(圖2(d))、陽坡裸地(圖2(e))、陰坡裸地(圖2(f))5類,其中陽坡裸地與陰坡裸地以研究區(qū)內最高海拔線為界線(圖1),海拔沿界線南北兩側逐漸降低,加上衛(wèi)星拍攝角度影響,使得二者地物的色調差別較大。陽坡裸地坡向朝南,色調以土黃色、棕色為主; 陰坡裸地坡向朝北,色調以棕黑色為主; 干涸河流的形狀與邊界較為明顯,寬度較窄; 冰雪主要分布在雪線較低的陰坡裸地上,色調以亮白色為主,形狀為長條狀,展布大多為南北向; 道路色調為墨青色、邊界明顯?；诓煌匚锏奶攸c,本文選擇提取地物的光譜特征、紋理特征、指數(shù)特征、地形特征、邊緣特征作為隨機森林分類所需的特征變量。光譜特征包括GF-2高空間分辨率影像數(shù)據(jù)4個波段的灰度特征(B,G,R,NIR)、影像灰度平均值特征(Img Mean)、影像灰度標準差特征(Img Std),以及對影像原始4個波段進行限制對比度自適應直方圖均衡化的CLAHE特征(CLAHE B,CLAHE G,CLAHE R,CLAHE NIR)。

(a) 花崗偉晶巖樣本 (b) 道路樣本 (c) 干涸河流樣本 (d) 冰雪樣本 (e) 陽坡裸地樣本 (f) 陰坡裸地樣本

紋理特征是一種全局特征,對噪聲具有較強抵御能力,能刻畫出圖像區(qū)域所對應地物表面的性質,具有旋轉不變的特性[12],且已有研究表明紋理特征有助于改善地表信息提取的精度[13-16]。本文在提取紋理特征前先對影像數(shù)據(jù)進行主成分分析變換,該方法通過應用正交變換減少相關數(shù)據(jù)的冗余,減少各組分間的相關性,盡可能將有用信息集中到新的組分中,且各新組分間互不相關,達到信息增強、降維及減少計算量的目的[17]。其中影像數(shù)據(jù)的第一主成分(PC1)包含所有波段中98.5%的特征信息,因此將PC1作為影像的光譜特征的同時,僅對PC1采用灰度共生矩陣方法提取紋理特征,選取均值(Mean)、方差(Variance)、均質性(Homogeneity)、對比度(Contrast)、相異性(Dissimilarity)、熵(Entropy)、相關性(Correlation)、角二階矩(angular second moment,ASM)等8種彼此相關性弱的紋理信息特征統(tǒng)計量作為紋理特征。

指數(shù)特征是通過一定公式的波段運算所獲取的灰度影像。指數(shù)特征可以消除地形差異影響,突出地物特征[18]?；陉幤侣愕卮嬖陉幱暗奶攸c,本文提取用于突出陰影區(qū)域與明亮區(qū)域對比的陰影指數(shù)(shadow index,SI)特征,其計算公式為:

SI=(B+G+R+NIR)/4

,

(1)

式中B,G,R,NIR分別為藍光、綠光、紅光及近紅外波段的反射率?；谘芯繀^(qū)內裸地范圍較大,本文提取突出裸地區(qū)域亮度的裸地區(qū)域指數(shù)(bare area index,BAI)特征,其計算公式為:

BAI=(B-NIR)/(B+NIR)

。

(2)

地形特征是描述地勢變化的特征,由于研究區(qū)內地勢起伏較大,且地物與地形因素有一定關聯(lián),因此有必要提取坡度(Slope)、坡向(Aspect)、高程等地形特征。本文從美國地質調查局獲取研究區(qū)30 m空間分辨率的ASTER GDEM數(shù)字高程模型(digital elevation model,DEM)數(shù)據(jù),并分別提取坡度、坡向,DEM地形特征。

邊緣特征是影像中特性(如像素灰度、紋理等)分布不連續(xù)處,影像周圍特性有階躍變化的像素集合,其中研究區(qū)內花崗偉晶巖、干涸河流、冰雪、道路的邊緣特征較明顯,邊緣灰度梯度變化較大。Canny邊緣檢測能較精確估算出圖像邊緣的強度、梯度方向,具有定位準確、單邊響應和信噪比高等優(yōu)勢[19]。將從GF-2高空間分辨率影像中提取出的光譜特征、紋理特征、指數(shù)特征、地形特征與邊緣特征等25個特征變量構成初始特征空間。

各類型地物選取25個影像樣本構成訓練樣本集,圖2為32像素×32像素的影像樣本塊展示。為提高分類檢測精確度,選取原則依據(jù)樣本要盡可能包含各個類別地物的各種不同類型特征,且單個樣本盡量僅包含單個類別地物。根據(jù)隨機森林實驗需求,將樣本集按照2∶8的比例隨機劃分為測試集和訓練集。

2.2 隨機森林算法及參數(shù)組合優(yōu)化

隨機森林算法最早是由Breiman[20]提出的一種機器學習算法,本文依據(jù)Gini系數(shù)最小原則構建分類和決策樹(classification and regression tree,CART),并基于多棵決策樹對樣本進行訓練,根據(jù)訓練得到的模型,對未知待測樣本類別進行預測的一種監(jiān)督學習分類算法[21]。

本文使用Python語言編程調用sklearn庫進行隨機森林模型訓練,采用隨機搜索交叉驗證(RandomizedSearchCV)與網(wǎng)格搜索交叉驗證(GridSearchCV)[22]的雙重方式進行參數(shù)優(yōu)化。單一的GridSearchCV是對每一組排列組合的參數(shù)進行遍歷,最后選取出最好的一組參數(shù),但由于參數(shù)排列組合數(shù)量過多,逐個遍歷面臨相當耗時與維度災難等問題,因此在進行GridSearchCV前先通過RandomizedSearchCV(設置隨機搭配超參數(shù)組合次數(shù)(n_iter)為200次,交叉驗證折數(shù)(cv)為5次)對取值范圍內的參數(shù)進行有限次的隨機排列組合,并對其進行遍歷,從中選取出一組較優(yōu)參數(shù)組合,再通過進行10折的GridSearchCV遍歷選取出最優(yōu)參數(shù)組合。本文最優(yōu)參數(shù)組合設置如表2所示。

表2 最優(yōu)參數(shù)組合

2.3 特征重要性評估與特征優(yōu)選

特征集中包含著不同屬性、不同尺度的多樣化的眾多特征變量,過多的特征變量會造成數(shù)據(jù)的冗余、分類時間過長、以及分類精度下降等問題[23],因此需要通過特征選擇篩選出對模型結果影響較大的特征變量構成新的特征子集。特征選擇主要是通過一定的準則對特征重要性排序,或者在分類精度保證的情況下選擇一組數(shù)量最小,且最具有代表性的特征子集[24]。本文先通過調用sklearn庫中基于Gini 系數(shù)原理進行特征重要性計算的feature_importances參數(shù)對每個特征進行重要性評估,對得到的特征重要性得分按降序排列(圖3),然后采用順序前向特征選擇法[25],在每一次迭代中選擇當前特征集中特征重要性得分最高的特征變量進行模型訓練,并獲取相應模型的測試集精度,依次增加特征變量數(shù)量直至全部納入,從中選取出測試集精度最高、特征變量數(shù)量相對最少的特征組合。

圖3 特征重要性排序

2.4 形態(tài)學運算

基于最佳參數(shù)組合與最優(yōu)特征變量組合的隨機森林模型對研究區(qū)各地物類型進行分類,對分類結果進行先膨脹后腐蝕的形態(tài)學閉運算操作,閉運算在填充同類型地物間的細微裂隙,使地物間更連接且更平滑的同時,還可以消除孤立噪聲點的影響。

3 實驗結果與精度分析

3.1 特征優(yōu)選評價

分別對特征空間中的25個特征變量通過順序向前特征選擇法逐一訓練后的結果進行評估,生成特征變量個數(shù)與測試集精度關系圖(圖4),由圖4可知: ①由首個特征變量逐個增加至16個特征變量參與模型訓練期間,測試集精度整體呈現(xiàn)迅速上升至平穩(wěn)上升趨勢,由最初的84.01%上升至96.24%; ②特征變量數(shù)量在17～25個期間,由于特征變量數(shù)量過多,產(chǎn)生數(shù)據(jù)冗余、訓練耗時增加等問題,導致分類器出現(xiàn)一定的過擬合現(xiàn)象,造成測試集精度有些許下降,由96.24%下降至95.94%,由此可知并不是特征變量數(shù)量越多測試集精度就越高。根據(jù)測試集精度結果,選取前16個特征變量組合作為最優(yōu)特征集,并基于該特征集進行隨機森林分類。

圖4 特征變量個數(shù)與測試集精度關系

3.2 特征變量分析

遙感影像中不同波段對于分類不同地物具有重要意義[24],波段間反射率值的差異有助于區(qū)分不同類型地物,通過計算不同特征變量在不同類型地物的反射率值,并對反射率值進行歸一化處理,得到各類型地物波譜曲線(圖5)。光譜特征(圖5(a))中的NIR波段的地物波譜曲線反射率值差異較大,可以很好地區(qū)分各類型地物,并且NIR波段的特征重要性得分最高; 紋理特征(圖5(b))中Mean波段的各類型地物反射率值有所差異,具有區(qū)分地物的能力,花崗偉晶巖的紋理特征在大多數(shù)波段反射率值都較高,原因為風化剝蝕等作用使得花崗偉晶巖巖石表面產(chǎn)生較多紋理; 指數(shù)特征(圖5(c))中陰坡裸地的陰影指數(shù)反射率值較低,與其他類型地物有明顯區(qū)分; 地形特征(圖5(d))中,由于研究區(qū)內高程起伏較大,高程范圍由3 671～4 029 m,高差相差358 m,使得DEM波段中不同類型地物反射率值差異較大。

(a) 光譜特征波譜曲線 (b) 紋理特征波譜曲線

3.3 分類結果與精度分析

3.3.1 隨機森林分類結果與精度分析

基于隨機森林算法對研究區(qū)的GF-2高空間分辨率影像數(shù)據(jù)進行花崗偉晶巖提取,其各類型地物分類結果如圖6所示。此次實驗結果共提取出84處花崗偉晶巖,其形狀大小不一,以脈狀產(chǎn)出為主,成群出現(xiàn),且輪廓較為明顯,與其他類型地物區(qū)分較好。在研究區(qū)內隨機均勻生成1 000個地面驗證點,其中花崗偉晶巖驗證點104個,道路驗證點19個,冰雪驗證點51個,陽坡裸地驗證點414個,干涸河流驗證點87個,陰坡裸地驗證點325個,通過計算地面驗證點與真實地物類型之間的混淆矩陣獲得分類結果的精度統(tǒng)計(表3)評價,其中總體精度達到93.1%; Kappa系數(shù)達到0.902; 花崗偉晶巖、冰雪、陽坡裸地及陰坡裸地共4種地物的用戶精度和生產(chǎn)者精度都高于93%、道路高于81%、干涸河流高于73%?；◢弬ゾr由于其自身的色調、紋理,以及分布地形海拔都較高等優(yōu)勢,在分類結果上效果較好,用戶精度高達94.24%、生產(chǎn)者精度高達98%,以上結果說明基于多特征的隨機森林算法在提取花崗偉晶巖信息中有較好的適用性。

表3 分類結果的精度統(tǒng)計

圖6 研究區(qū)地物分類結果

3.3.2 CLAHE特征層討論

本文新引入的4個CLAHE特征層中,CLAHE B與CLAHE R作為特征變量參與了隨機森林分類,為驗證CLAHE B和CLAHE R對分類結果產(chǎn)生的影響,故對優(yōu)選的16個特征變量中去除CLAHE B和CLAHE R這2個特征變量,只保留余下的14個特征變量參與隨機森林分類,并通過計算驗證相同的1 000個地面驗證點與真實地物類型之間的混淆矩陣,并獲得分類結果的精度(表4)評價,其中總體精度為90.4%,較最佳分類結果的總體精度下降2.7百分點; Kappa系數(shù)為0.864,較最佳分類結果的Kappa系數(shù)下降0.038; 各地物類型的用戶精度與生產(chǎn)者精度都有不同程度的下降。

表4 14個特征的分類結果精度統(tǒng)計

CLAHE方法的優(yōu)越性在于可以通過改變不同地物間的對比度,以突出色調上的差異,相比普通的直方圖均衡化方法只增強全局影像對比度,CLAHE方法在此基礎上引入了自適應性與限制性,自適應性改進了增強全局對比度中只考慮整體影像區(qū)域,而忽略局部影像區(qū)域的情況,進而注重影像的局部像素差異,在突出局部影像對比度的同時獲得更多細節(jié)特征; 限制性是在自適應性基礎上的進一步改進,目的是為了避免均衡化過程中產(chǎn)生的色調不連續(xù)與對比度過度增強等問題。以精度下降最為明顯的干涸河流為例,且為了效果更加直觀,故采用真彩色的方式展示(圖7),CLAHE真彩色樣本影像中(圖7(b)),干涸河流的亮度要明顯亮于陽坡裸地,且對比度更加突出,邊緣輪廓更加清晰,2類地物在色調上有明顯差異,有利于分類精度的提升,證實了CLAHE方法在區(qū)分色調較為相近地物時的有效性,同時CLAHE特征增強了花崗偉晶巖與其他類型地物間亮度與對比度的差異,有助于提升花崗偉晶巖的提取效果,使其用戶精度和生產(chǎn)者精度分別上升0.97百分點與1百分點。

(a) 真彩色影像 (b) CLAHE真彩色影像

4 結論

1)本文基于GF-2高空間分辨率影像數(shù)據(jù)結合特征優(yōu)選的隨機森林算法提出花崗偉晶巖的提取方法。地物分類結果總體精度為93.10%,Kappa系數(shù)為0.902; 花崗偉晶巖的用戶精度及生產(chǎn)者精度分別為94.24%和98.00%,證實了該方法的有效性。

2)特征變量個數(shù)持續(xù)增多會產(chǎn)生數(shù)據(jù)冗余及計算量過大等問題,導致過擬合現(xiàn)象發(fā)生。此次研究在特征優(yōu)化中確定特征變量個數(shù)為16個時分類精度最高。

3)本文新引入的限制對比度自適應直方圖均衡化(CLAHE)特征層中的CLAHE B和CLAHE R的特征重要性得分較高,并參與隨機森林分類。CLAHE特征變量具有提升影像對比度、突出不同類型地物間色調差異等優(yōu)勢,有助于區(qū)分色調相近的地物,有利于地物分類精度的提升。