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

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

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

0 引言

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

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

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

1 研究區概況與數據源

1.1 研究區概況

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

圖1 研究區遙感影像

1.2 數據源及預處理

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

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

2 研究方法

2.1 特征提取及樣本選擇

花崗偉晶巖(圖2(a))在分布位置及產出形態上受構造控制作用明顯,研究區內斷裂沿NW-SE向分布,因此花崗偉晶巖大多沿NW-SE向展布,呈條狀、脈狀產出; 在影像色調上,其多呈灰白色或略帶土黃色; 由于研究區內海拔較高、風沙較多、風化能力較強,且花崗偉晶巖晶體較大,具有較強的抗風化性,因而在高度上要突出于周圍地物,并且其表面紋理粗糙,邊界較明顯。研究區內其他地物類型主要分為道路(圖2(b))、干涸河流(圖2(c))、冰雪(圖2(d))、陽坡裸地(圖2(e))、陰坡裸地(圖2(f))5類,其中陽坡裸地與陰坡裸地以研究區內最高海拔線為界線(圖1),海拔沿界線南北兩側逐漸降低,加上衛星拍攝角度影響,使得二者地物的色調差別較大。陽坡裸地坡向朝南,色調以土黃色、棕色為主; 陰坡裸地坡向朝北,色調以棕黑色為主; 干涸河流的形狀與邊界較為明顯,寬度較窄; 冰雪主要分布在雪線較低的陰坡裸地上,色調以亮白色為主,形狀為長條狀,展布大多為南北向; 道路色調為墨青色、邊界明顯?；诓煌匚锏奶攸c,本文選擇提取地物的光譜特征、紋理特征、指數特征、地形特征、邊緣特征作為隨機森林分類所需的特征變量。光譜特征包括GF-2高空間分辨率影像數據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) 陰坡裸地樣本

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

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

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

,

(1)

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

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

。

(2)

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

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

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

2.2 隨機森林算法及參數組合優化

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

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

表2 最優參數組合

2.3 特征重要性評估與特征優選

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

圖3 特征重要性排序

2.4 形態學運算

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

3 實驗結果與精度分析

3.1 特征優選評價

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

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

3.2 特征變量分析

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

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

3.3 分類結果與精度分析

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

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

表3 分類結果的精度統計

圖6 研究區地物分類結果

3.3.2 CLAHE特征層討論

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

表4 14個特征的分類結果精度統計

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

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

4 結論

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

2)特征變量個數持續增多會產生數據冗余及計算量過大等問題,導致過擬合現象發生。此次研究在特征優化中確定特征變量個數為16個時分類精度最高。

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