基于影像特征隨機森林的自然保護區礦區動態監測

張凱選，李修會，余卓淵

(1.遼寧工程技術大學，遼寧 阜新 123000；2.中國科學院 地理科學與資源研究所，北京 100101；3.中國科學院大學，北京 100049)

賀蘭山國家自然保護區跨溫帶草原與荒漠兩大植被區域，其特殊的地理位置形成西北干旱地區生物多樣性中心和生物資源寶庫[1]，并蘊含豐富的礦產資源。近年來，自然保護區內非法開采現象嚴重[2]，造成土地侵占損毀和水體污染等現象[3]，導致賀蘭山生態脆弱性明顯降低[4]，生態環境惡化，生物多樣性降低。鑒于以上問題，2017年內蒙古人民自治區政府和寧夏回族自治區人民政府分別出臺“自然保護區內工礦企業退出方案通知”和“賀蘭山國家級自然保護區生態環境綜合整治推進工作方案”，使合法礦山企業有序退出自然保護區，取締非法開采點，并進行礦山生態修復，從而盡可能恢復礦山原有狀態。針對如何有效監測自然保護區內非法開采活動和礦業權是否有序退出的問題，遙感技術的應用，彌補了以往人工調查方式的不足，為礦產資源監測與調查提供了便利。基于傳統的目視解譯和監督分類的多光譜影像數據對礦山識別與監測時，存在以下問題：分類精度不高，效率低下，人力物力成本高，導致應用效果不佳[5-6]。隨機森林不僅可以結合多個有效特征信息解決“同物異譜”、“異物同譜”的問題，提高影像解譯精度，還能經過同質分類器投票決定最優分類結果，該方法具有原理簡單、易于實現等優點。

近年來，大量學者采用隨機森林對不同的地物覆蓋信息進行提取。王紅等[7]以Landsat-8影像結合溫度反演，提取冬閑田為例，對比隨機森林和支持向量機結果，認為前者的精度優于后者；趙亞杰等[8]構建多特征融合影像從Landsat-8、Sentinel-2數據提取地物信息，實現大范圍的精細分類；周正龍等[9]利用多時相的Landsat-8影像，針對福建平潭實驗區，利用隨機森林算法提取土地利用信息，獲得較好效果。現階段，在礦區提取地物信息應用較廣的有面向對象分類法[10]、支持向量機法[11]、最大似然分類[12]等方法。目前，隨機森林被廣泛用于一般地物提取，在自然保護區礦區使用影像特征結合隨機森林用于礦區地物覆蓋信息提取的研究較少。

本研究采用基于影像特征隨機森林算法分別對賀蘭山自然保護區2016—2019年4個時期的Landsat-8多光譜遙感影像進行礦區地物覆蓋信息提取，通過融合多特征信息的隨機森林分類方法與支持向量機方法進行對比，進一步討論提取自然保護區礦區地物覆蓋信息有效性和可行性，獲取精度較高的動態變化監測結果，以期為礦山監測與生態保護提供技術支持和決策依據。

1 研究區及數據來源

賀蘭山，位于寧夏回族自治區與內蒙古自治區交界處，如圖1所示，賀蘭山地理坐標38°19″～39°27″N，105°41″～106°45″，海拔高度2 000～3 000 m，主峰敖包疙瘩高3 556 m。屬溫帶大陸性氣候，季節變化明顯，干旱少雨。區域內動植物種類豐富，有陸棲脊椎動物約177種，野生維管植物約690種，具有很高經濟和開發價值。賀蘭山自然保護區內礦產資源非常豐富，其中煤炭資源主要分布在北部汝萁溝、石炭井等地，以無煙煤和焦煤為主；非煤炭資源主要分布在正義關、紅果子及大小口子一帶，以石英砂巖和賀蘭山石為主。

本文影像數據使用Landsat陸地資源衛星數據，選取2016—2019年4～5月份Landsat-8 OLI無云影像數據(表1)，該數據可從USGS GloVis(https：//glovis.usgs.gov/app)獲取；使用賀蘭山自然保護區與礦業權范圍矢量數據對影像進行裁剪和顯示；利用天地圖高分影像作為分類參考；自然資源部信息中心提供2019年礦業權信息登記庫的礦業權信息數據，結合分類變化數據進行輔助分析。

利用ENVI軟件對影像進行輻射校正、大氣改正，以此減少和消除畸變，增強影像，提高解譯精度，并對圖像進行裁剪。

表1 研究區域Landsat-8影像信息

2 研究方法

2.1 技術路線

本文對影像數據進行預處理，為了提高分類精度，提取影像紋理特征和計算相關指數，構建多特征信息融合影像，利用隨機森林[13-14]與SVM[15]算法對融合后的影像進行分類；根據《土地利用現狀分類》(GB/T21010—2007)國家標準和研究區的實際情況，本研究制定的分類系統共2級6類型，即：建設用地、水體、植被、采礦場、排土場和其他土地；結合歸一化提取值與天地圖高清影像，選取實驗和驗證樣本進行分類精度評價，并對分類結果進行分類后處理(合并和平滑)，如圖2所示。

圖2 地物覆蓋提取路線Fig.2 Ground feature cover extraction route

RF-1采用預處理后未添加特征信息的影像；RF-2采用紋理構建的影像；RF-3采用光譜特征構建的影像；RF-4采用紋理特征和光譜特征構建的影像；SVM采用光譜特征和紋理特征構建的影像。

2.2 特征選擇

賀蘭山地理位置特殊，在北部山群由于植被覆蓋較少，多以裸露巖石或土為主，植被在光譜顯示與裸巖相似，呈現褐色；礦區內部分水體處在特殊位置，如在山體的陰影處，其光譜顯示與開采場光譜近似呈現黑色；部分排土場在光譜顯示上與開采場近似，通過目視難區分；精確的實驗樣本是準確分類的前提，由于存在以上光譜近似顯示的問題，憑借經驗直接提取樣本容易誤提取，從而降低分類的準確性。為了提高樣本的準確度，分別添加歸一化差分植被指數NDVI；改進歸一化差異水體指數MNDWI；歸一化差值土壤指數NDSI，以上指數不僅可以增加分類的準確性，而且可以為精確提取實驗樣本提供輔助，計算公式如下：

NDVI=(ρNIR-ρR)/(ρNIR+ρR)

(1)

MNDWI=(ρG-ρMIR)/(ρG+ρMIR)

(2)

NDSI=(ρR-ρB)/(ρR+ρB)

(3)

式中：ρR、ρG、ρB、ρNIR和ρMIR是遙感影像中的紅波段、綠波段、藍波段、近紅外波段和中紅外波段。

在NDVI、MNDWI和NDSI指數提取圖中(圖3)可以清晰地觀察到自然保護區植被、水體和土壤大致分布。通過多光譜圖像與NDVI圖像結合可以提取出植被樣本，區分土壤(裸巖)；多光譜圖像與MNDWI圖像可以提取陰影下存在的水體；多光譜圖像與NDSI圖像可以提取其他土地類型，區分植被類型；多光譜圖像結合NDVI或MNDWI中灰度級分別提取開采場和排土場。

圖3 影像歸一化指數提取(a—多光譜；b—NDVI；c—MNDWI；d—NDSI)Fig.3 Image normalization index extraction(a—Multispectral；b—NDVI；c—MNDWI；d—NDSI)

根據遙感影像中的地物獨特的紋理結構，提取灰度共生矩陣(The Gray Level Co-occurrence Matrix，GLCM)，選取對比度(Contrast)、熵(Entropy)、角二階矩和相關性(Correlation)4個紋理特征，如圖4所示。

圖4 影像紋理特征提取(a—對比度；b—熵；c—角二階矩；d—相關性)Fig.4 Image texture feature extraction(a—Contrast；b—Entropy；c—Angularsecond moment；d—Correlation)

2.3 分類方法

隨機森林[16]是由LEO Breiman 等提出，由一系列CART決策樹組合的集成分類器，分類結果是由基分類器采取投票形式，遵從少數服從多數原則得出。在遙感影像分類和變化監測領域均獲得較好成果。隨機森林方法主要過程分為訓練與分類，訓練過程通過 Bootstrap自抽樣方法從原始訓練樣本中隨機抽取K個樣本子集，未選中的訓練樣本視為袋外數據(Out of Bag，OOB)；然后對于K個樣本子集分別構建對應的決策樹，組成隨機森林模型。分類過程是對每一棵決策樹的決策結果使用簡單多數投票法進行綜合，得到最終的分類結果。

(4)

式中：H(x)—隨機森林最終分類結果；hi(x)—單一決策樹模型分類結果；Y—輸出變量(目標變量)；I(?)—示性函數。

在運行隨機森林分類之前需要對決策樹棵數、最小樣本節點和最小雜質等參數設置，通過實驗，各實驗設置400、1和0作為分類參數。

3 結果與分析

3.1 影像分類結果

完成不同組合分類輸出結果(圖5)發現，存在錯分情況如下：其他用地被錯分，原因是在運輸礦石時礦石散落對道路及開采區周圍產生污染，在遙感影像上的光譜特征與開采場或者排土場相似，因此被解譯成開采場或者排土場，通過與天地圖高清圖像對比進行修正。在4種RF不同組合中，整體分類效果較好，其中RF-1是經過預處理，未加入特征信息，分類效果在細節方面較差，植被和其他用地存在碎點多，部分被錯分，且在地物邊緣契合度低；分別在預處理后影像的融合NDVI、MNDWI、NDSI和紋理特征的RF-2、RF-3，植被的破碎點減少，分類地物的整體性提升，分類效果均優于RF-1；融合NDVI、MNDWI、NDSI和紋理特征的RF-4，分類地物破碎點進一步減少，地物邊緣契合度提高，與相同組合策略采用SVM分類器相比，分類效果均有提升。

圖5 不同組合策略分類方法對比Fig.5 Comparison of classification methods under different feature combination strategies

通過混淆矩陣(表2)可以對不同特征組合策略的分類結果進行精度評價。采用同質分類器，不同組合策略可發現，分類的整體精度的關系為RF-4>RF-3>RF-2>RF-1，Kappa系數關系為RF-4>RF-3=RF-2>RF-1；采用相同特征信息的隨機森林的分類整體精度為92.35%、0.91，較SVM分類器提高了1.45%、0.04。SVM分類器分類整體精度較未添加特征信息的RF-1僅提升了0.06%，效果提升有限。在本文中，融合NDVI、MNDWI、NDSI和紋理特征，采用隨機森林分類器分類效果達到預期，得到較好的分類效果。

表2 不同特征組合策略的混淆矩陣Table 2 Confusion matrix of different feature combination strategies

3.2 賀蘭山礦業權及地物變化分析

1)賀蘭山2016—2019年礦業權數據變化

從自然資源部信息中心礦業權信息登記庫(2019)中獲取賀蘭山國家自然保護區(內蒙古和寧夏)累計可查詢到礦業權數量共有79個，其中寧夏賀蘭山國家自然保護區內70個，內蒙古賀蘭山國家自然保護區9個，礦業權信息如表3所示，以時間為節點，截至2019礦業權到期64個，占礦業權總數的81.01%；礦業權已注銷35個，占礦業權總數的44.30%(其中煤礦22個，非煤礦13個)；2017年新設采礦權9個，累計采礦權79個，之后未設置新的采礦權，與2017出臺“方案”不允許在自然保護區設立新礦業權要求吻合；以2020年為節點，礦業權未到期的15個，其中有4個礦業權已注銷，實際具備開采條件礦業權降至11個，占礦業權總數的13.92%。

表3 礦業權信息

2)賀蘭山礦山分布規律

賀蘭山自然保護區礦山主要分布在平均海拔1 100～2 600 m，礦山開采區域集中在賀蘭山北部地勢較緩地區，大中型采礦集群如：汝萁溝、石炭井、石嘴山、呼嚕斯太和沙巴臺等采礦權共有59個，占總體的74.68%，其中煤礦48個，非煤礦11個，分別占60.76%、13.92%。其余散落分布在賀蘭山南部區域；汝萁溝、石炭井、呼魯斯太和沙巴臺采礦集群主要以煤為主，石嘴山采礦集群主要以建筑石料為主；已注銷的礦業權48(含礦業權未到期申請注銷)個，占比60.76%，其中石嘴山和沙巴臺采礦集群礦業權已全部注銷。

3)賀蘭山礦山地物變化分析

礦山活動是否劇烈體現在礦業權內地物覆蓋變化幅度。因此，提取研究區2016—2019地物覆蓋信息，繪制變化趨勢圖如圖6所示；結合礦業權變更信息數據(表3)綜合分析變化趨勢，在2017年開采場的面積為37.63 km2，較2016年33.12 km2增加了4.51 km2(圖6-a)，主要原因是在2017年1月至5月先后新設立9個礦業權，排土場和建筑用地面積進而有所增加(圖6-b，圖6-c)，主要由植被和其他用地轉化為開采場、排土場和建筑用地；2018年開采場面積大幅度減少，較2017年少了14.14 km2，減少的面積占2017年開采場總面積的37.58%，同時，排土場和建筑用地面積分別減少了1.05 km2、2.9 km2(圖6-b，圖6-c)。通過對開采場和排土場進行覆土回填、種植植物和建筑物拆除，開采場、排土場和建筑用地減少的面積主要向其他用地轉化；2019年開采場面積進一步減少至16.43 km2，減少了7.06 km2(圖6-a)，截至2019年減少的總面積占最高值的56.34%，植被、水體和其他用地面積增加至790.74 km2、1 884.08 km2(圖6-d，圖6-e，圖6-f)。

圖6 研究區地物覆蓋變化(a—開采場；b—排土場；c—建筑用地；d—植被；e—水體；f—其他用地)Fig.6 Land cover changes in the study area(a—Mining site；b—Dump site；c—Construction land；d—Vegetation；e—Water body；f—Other land)

將2016—2019研究區分類地物覆蓋影像與礦業權范圍矢量圖進行疊加并結合表3數據信息分析可以得出如下問題，礦業權在有效期到期未及時變更狀態下存在非法開采行為；在礦業權集中的區域，部分開采活動存在越界行為；存在侵占損毀土地的堆放區，為防止衛星監測將其偽裝成植被等行為；在礦業權范圍以外存在小而多的開采圖斑，疑似無證開采；隨著“方案”的實施，礦區變化明顯，開采場開采和礦區建筑面積等逐年縮減，植被和其他用地面積逐年增加，說明“方案”的實施得到了很好的效果。

4 結論

本文以賀蘭山國家自然保護區為研究區，選取了Landsat-8遙感影像作為數據源，結合多種影響因子，利用隨機森林有效地實現了對賀蘭山國家自然保護區礦業權內的地物覆蓋信息提取，并對賀蘭山地區2016—2019年礦業權內地物覆蓋進行動態監測，得出以下結論：

1)采用融合多特征的RF-4組合在隨機森林算法提取礦山地物覆蓋信息中分類整體精度和Kappa系數分別為 92.35%、0.91。優于SVM分類算法并取得良好的分類效果。

2)在礦業權界內將提取地物覆蓋信息疊加，存在越界開采、無序開采、亂采濫采等非法開采活動，特點是開采持續周期較長，且未對其進行修復。

3)研究區內的礦山開采場、排土場和建筑面積由于2017年設立新的礦業權，導致面積增加，2018—2019年連續下降，減少幅度達56.34%，在礦區地物覆蓋面積變化上能很好地體現“方案”實施起到的重要作用，并在很大程度上抑制了越界開采、無序開采和亂采濫挖的現象。

本實驗存在一些不足之處，由于選用Landsat-8遙感影像，在空間分辨率上有所限制，僅對開采場的面積進行監測，未能對不同礦區類型進行細分，進而分別進行動態監測，需要進一步加強研究。