基于影像特征隨機森林的自然保護區(qū)礦區(qū)動態(tài)監(jiān)測

張凱選，李修會，余卓淵

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

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

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

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

1 研究區(qū)及數(shù)據(jù)來源

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

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

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

表1 研究區(qū)域Landsat-8影像信息

2 研究方法

2.1 技術(shù)路線

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

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

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

2.2 特征選擇

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

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

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

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

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

2.3 分類方法

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

(4)

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

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

3 結(jié)果與分析

3.1 影像分類結(jié)果

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

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

通過混淆矩陣(表2)可以對不同特征組合策略的分類結(jié)果進行精度評價。采用同質(zhì)分類器，不同組合策略可發(fā)現(xiàn)，分類的整體精度的關(guān)系為RF-4>RF-3>RF-2>RF-1，Kappa系數(shù)關(guān)系為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和紋理特征，采用隨機森林分類器分類效果達到預(yù)期，得到較好的分類效果。

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

3.2 賀蘭山礦業(yè)權(quán)及地物變化分析

1)賀蘭山2016—2019年礦業(yè)權(quán)數(shù)據(jù)變化

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

表3 礦業(yè)權(quán)信息

2)賀蘭山礦山分布規(guī)律

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

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

礦山活動是否劇烈體現(xiàn)在礦業(yè)權(quán)內(nèi)地物覆蓋變化幅度。因此，提取研究區(qū)2016—2019地物覆蓋信息，繪制變化趨勢圖如圖6所示；結(jié)合礦業(yè)權(quán)變更信息數(shù)據(jù)(表3)綜合分析變化趨勢，在2017年開采場的面積為37.63 km2，較2016年33.12 km2增加了4.51 km2(圖6-a)，主要原因是在2017年1月至5月先后新設(shè)立9個礦業(yè)權(quán)，排土場和建筑用地面積進而有所增加(圖6-b，圖6-c)，主要由植被和其他用地轉(zhuǎn)化為開采場、排土場和建筑用地；2018年開采場面積大幅度減少，較2017年少了14.14 km2，減少的面積占2017年開采場總面積的37.58%，同時，排土場和建筑用地面積分別減少了1.05 km2、2.9 km2(圖6-b，圖6-c)。通過對開采場和排土場進行覆土回填、種植植物和建筑物拆除，開采場、排土場和建筑用地減少的面積主要向其他用地轉(zhuǎn)化；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 研究區(qū)地物覆蓋變化(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研究區(qū)分類地物覆蓋影像與礦業(yè)權(quán)范圍矢量圖進行疊加并結(jié)合表3數(shù)據(jù)信息分析可以得出如下問題，礦業(yè)權(quán)在有效期到期未及時變更狀態(tài)下存在非法開采行為；在礦業(yè)權(quán)集中的區(qū)域，部分開采活動存在越界行為；存在侵占損毀土地的堆放區(qū)，為防止衛(wèi)星監(jiān)測將其偽裝成植被等行為；在礦業(yè)權(quán)范圍以外存在小而多的開采圖斑，疑似無證開采；隨著“方案”的實施，礦區(qū)變化明顯，開采場開采和礦區(qū)建筑面積等逐年縮減，植被和其他用地面積逐年增加，說明“方案”的實施得到了很好的效果。

4 結(jié)論

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

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

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

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

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