基于SVM 的地下采煤區沉陷災害發育敏感性分區研究

薛永安，鄒友峰，張文志，張明媚，柳廣春，宋明偉

(1.河南理工大學 測繪與國土信息工程學院，河南 焦作 454000；2.太原理工大學 礦業工程學院，山西太原 030024；3.河南省采空區場地生態修復與建設技術工程研究中心，河南 焦作 454000；4.山西能源學院地質與測繪工程系，山西 晉中 030600；5.中國科學院地理科學與資源研究所 中國科學院陸地表層格局與模擬重點實驗室，北京 100101)

地下礦產采出后，開采區周圍巖土體的原始應力平衡狀態被破壞，巖土體上出現位移與變形[1]，這是一個十分復雜的過程，是采礦、地質等許多因素及多場耦合綜合影響的結果，在地表上表現為塌陷坑與地裂縫[2]，屬于礦山地質災害的范疇，具有漸進性特點，本文稱為沉陷災害。沉陷災害是地下采煤擾動導致地表快速或緩慢沉降的一種地質災害，在礦區分布廣泛，既造成大量土地資源損毀，又嚴重破壞當地的地形地貌景觀，還會嚴重影響周邊植物生長[2]，也是誘發斜坡地質災害和影響區域生態地質環境的重要因素。因此，開展地下采煤區沉陷災害發育敏感性分區研究對國土綜合整治與生態修復、下伏老采空區重大工程項目設計與穩定性監測、地質災害防治與預警等至關重要，同時對地下采煤區生態環境修復戰略[3]具有助力作用。

災害敏感性是指在特殊地形或某些因素作用下發生災害的可能性[4]，災害敏感性評價則是針對某一地區現存或潛在災害空間分布概率的定性或定量分析。目前，地質災害敏感性評價工作已經取得眾多研究成果[5]，針對地表沉陷災害敏感性評價，主要集中在開采地下水引起的地表沉降[6-8]和廢棄礦區地表沉降的災害敏感性評價研究[9-11]，提出了眾多評價模型并被應用和開展對比分析[12-14]，為沉陷災害發育敏感性分區研究提供了評價模型選擇參考。敏感性評價模型一般可分為定性模型和定量模型，隨著研究的不斷深入，更多評價模型被引入該領域[15-16]，為災害風險評估與定量評價提供了良好的模型基礎。目前，常用的評價模型主要包括：頻率比(Frequency Ratio,FR)[17]、確定性系數(Certainty Factor,CF)[18]、支持向量機(Support Vector Machine,SVM)[19]、隨機森林(Random Forest,RF)[7-8]和自適應神經模糊推理系統(Adaptive-Network-Based Fuzzy Inference System,ANFIS)[20]等。在模型應用中，確定地表沉陷影響因子是沉陷災害發育敏感性分區評價的前提，如有研究者認為覆巖巖性、地表坡度、采空區范圍、巖體力學性質、地下水等是誘發地表沉陷的主要影響因素[21-22]，沉陷災害敏感性評價則是在上述基礎上所做的進一步研究工作[23]。評價因子的可靠性選擇，評價模型的適用性對比，為開展地下采煤區沉陷災害發育敏感性分區研究提供了豐富的理論與實踐參考[21-23]。

山西省太原市西山地區地下采煤強度高，沉陷災害及次生地質災害分布廣泛，給人民生活安全及經濟發展帶來較大影響。文獻[24]開展了崩塌、滑坡、不穩定斜坡和地面塌陷4 種地質災害與發育特征因子之間關聯性的空間統計分析，為太原西山地區沉陷災害發育敏感性評價提供了特征因子選擇參考。文獻[25]針對斜坡地質災害敏感性評價中地勢起伏度提取最佳尺度開展了研究，為太原西山地區地質災害敏感性評價時地勢起伏度因子可靠提取提供了依據。

針對地下采煤區沉陷災害發育敏感性分區，筆者選擇在小樣本情況下應用廣泛的SVM 模型為敏感性評價模型，分別以不同年度核查編錄的沉陷災害作為模型建立與預測能力驗證數據，開展太原市西山前山區沉陷災害發育敏感性分區研究，以期為地下采煤區沉陷災害及次生地質災害防治提供科學的方法與數據支持。

1 研究區與數據源

1.1 研究區概況

太原市西山前山區位于山西省太原市西部山區(圖1)，北西以尖草坪區、萬柏林區和晉源區行政界線為邊界，南東以山前坡腳線為界，總面積441.063 km2。區內采煤活動歷史悠久，已近百年，煤礦開采主要為地下開采方式，既有國有大礦，包括西銘礦、杜兒坪礦、官地礦、白家莊礦和西峪煤礦，還有分布廣泛且數量龐大的鄉辦、村辦、聯營小型煤礦，逐漸在空間上形成了大礦與小礦既有鄰近、又有重疊的地下采煤區，地表沉陷情況極其復雜，沉陷災害及次生地質災害分布廣泛，部分地裂縫發育于山頂基巖或山體出露巖層(圖2)。

圖1 太原市西山前山區地理位置Fig.1 Geographical map of Xishan area of Taiyuan City

圖2 太原市西山前山區沉陷災害Fig.2 Subsidence disaster in Xishan area of Taiyuan City

1.2 數據源

1) 沉陷災害數據

收集到研究區2012 年、2014 年核查編錄的地面塌陷與地裂縫災害共113 處，提取地裂縫的中心位置，以點形式與地面塌陷構建沉陷災害點圖層，其空間分布如圖1 所示。其中，2012 年79 處，2014 年34 處，數據來源于山西省太原市規劃和自然資源局。

2) DEM 數據

在地理空間數據云網站下載了研究區的ASTER GDEM V2 數據，并重采樣為30 m 分辨率的DEM 數據，作為研究區的數字地貌因子提取基礎數據，高程中誤差為±22.723 m[25]。

3) 地質數據

研究區1∶200 000 數字地質圖來源于山西省太原市規劃和自然資源局。

2 研究方法

沉陷災害發育敏感性分區步驟如下：

(1) 收集研究區實際核查編錄的2012 年和2014年的沉陷災害空間分布數據，以2012 年發育的沉陷災害為訓練樣本建立模型，以2014 年發育的沉陷災害為驗證樣本評價模型預測精度。

(2) 從地形地貌和地質因素兩個方面出發建立評價因子組合。

(3) 以評價因子提取結果輸入基于機器學習思想的SVM 模型開展沉陷災害發育敏感性評價，以受試者特征曲線(Receiver Operating Characteristic,ROC)下面積(Area Under the ROC Curve,AUC)作為定量評價指標對分區結果進行精度驗證。

(4) 以自然間斷點法進行等級劃分并制作研究區沉陷災害發育敏感性分區圖。

總體技術路線如圖3 所示。

圖3 總體技術路線Fig.3 Technology roadmap for the study

2.1 評價因子

針對地下采煤活動對礦區微地貌、土地利用等帶來的擾動影響，文獻[26-27]選擇研究區中西部的杜兒坪煤礦南三盤區為對象，其中，文獻[26]采用1∶10 000地形圖所生產的3 期DEM 數據分別代表礦區采前、采中和采后3 個時期的地形表面，開展了數字地貌時空演變特征分析，結果表明坡度與地勢起伏度統計結果隨時間推進而變化明顯，坡向也有一定的影響。因此，DEM 數據的獲取時間對數字地貌因子分析結果是否可靠具有重要約束，災害敏感性預測研究中應盡可能采用與沉陷災害發育時間相匹配的DEM 數據進行數字地貌因子提取。文獻[27]則以采前、采中和采后的5 期遙感影像為數據源，采用最小距離分類法進行了礦區土地利用信息提取與演變分析，結果表明，土地利用類型轉移過程和方向與煤炭開采擾動過程相吻合。因此，遙感影像時相不同，所提取的土地利用信息則不同，以特定年份的遙感影像所提取的土地利用數據和植被覆蓋數據與歷史編錄的沉陷災害數據進行相關分析，其結果不確定性極大。

沉陷災害影響因子眾多，在開采沉陷預計中多采用地形地貌、地質因素、采煤方法、煤層厚度、埋深、工作面尺寸等因子[1,28]，專家學者基于不同目的建立了不同的影響因子序列。受地下采煤相關數據獲取困難、定量表述不足制約，并考慮方法的通用性，本文在現有研究的基礎上[24-25]，綜合考慮地下采煤區數字地貌特征和土地利用時空演變的影響[26-27]，從地形地貌與地質因素出發，選取高程、坡度、坡向、地勢起伏度、地面曲率、地層巖組、地質構造共7 個因子作為沉陷災害發育敏感性評價因子組合(表1)。

表1 研究區沉陷災害敏感性評價因子Table 1 Assessment factors of subsidence disaster sensitivity in the study area

以ASTER GDEM V2 數據在ArcGIS 軟件中提取高程、坡度、坡向、地勢起伏度和地面曲率，基于1∶200 000 地質圖提取地層巖組與地質構造信息，利用ArcGIS 軟件制作研究區2012 年沉陷災害發育敏感性評價因子圖(圖4)。

經圖4 統計，沉陷災害在各行政區之間的空間分布規模差異性較大，2012 年的平均發育密度為0.179 處/km2，萬柏林區的發育密度最高，達到0.336 處/km2，高于研究區平均水平。其次是晉源區，為0.034 處/km2，發育密度相對較低。

圖4 太原市西山前山區沉陷災害評價因子Fig.4 Assessment factors of the subsidence disaster in Xishan area of Taiyuan City

2.2 支持向量機(SVM)

支持向量機是一種基于最大間隔的線性判別分類方法，其功能強大，可以最大化模型的預測準確度，被認為是目前針對小樣本統計估計和預測學習的最佳理論，最早被用于研究小樣本情況下的機器學習[29]，在地質災害敏感性評價研究中應用較多[13,30]。

SVM 通過引入核函數有效地解決了非線性分類問題，使敏感性評價中的非線性分類計算復雜度不再取決于空間維數，而是取決于用于建模的沉陷災害樣本數量，尤其是其中支持向量的災害點數量。顯然，核函數是SVM 克服維數災難實現沉陷災害敏感性評價有效分類的關鍵。SVM 常用核函數(表2)主要有線性核函數(Linear Kernel,LN)、多項式核函數(Polynomial Kernel,PL)、徑向基核函數(Radial Basis Function,RBF)和Sigmoid 核函數(Sigmoid Kernel,SIG)[13]。

表2 SVM常用核函數[31]Table 2 Common kernel functions of SVM[31]

目前，地質災害敏感性評價中徑向基核函數應用較多[30]，但研究區內沉陷災害樣本數量較少，應對比4種核函數SVM 模型的預測精度，選擇最優核函數建立SVM 預測模型。

沉陷災害敏感性分區SVM 模型預測步驟如下：

(1) 建立研究區沉陷災害數據集，包括訓練樣本集和驗證樣本集。

(2) 提取沉陷災害點在評價因子圖中的值，作為屬性值記錄下來，包括：高程、坡度、坡向、地勢起伏度、地面曲率、地層巖組和距地質構造的距離。

(3) 建立4 種SVM 模型，包括LN-SVM 模型、PL -SVM 模型、RBF-SVM 模型和SIG-SVM 模型。

(4) 采用訓練樣本對4 種SVM 模型進行訓練，統計評價因子在每一種模型中的權重排序，分析模型預測正確率與錯誤率，避免出現過擬合和欠擬合。

(5) 選擇預測正確率最高的SVM 模型作為沉陷災害敏感性分區預測模型。

(6) 采用驗證樣本評價最優SVM 模型的預測精度，對比訓練樣本的預測精度，分析預測結果的合理性。

2.3 精度評價

敏感性評價是一個二值分類問題，一般將樣本分為正類和負類，現有研究中通常采用真正類率和真負類率構建受試者特征曲線(ROC)，利用ROC 曲線下的面積(AUC 值)作為敏感性分區預測精度的度量[13,31]，AUC 值的范圍為[0,1]，值越大越好，表明模型預測能力越強，災害敏感性分區預測結果精度越高。

根據AUC 值的大小可以劃分敏感性分區預測模型的精度等級(表3)[13]。

表3 敏感性預測模型精度等級[13]Table 3 Accuracy levels of sensitivity prediction model[13]

對于二值分類問題，假設C1為正類，C2為負類，對應的2×2 混淆矩陣見表4。其中，TP 表示分類方法準確預測為C1的點的數目；FP 表示分類方法預測為C1，但實際屬于C2的點的數目；FN 表示分類方法預測為C2，但實際屬于C1的點的數目；TN 表示分類方法準確預測為C2的點的數目。

表4 二值分類混淆矩陣Table 4 Confusion matrix of binary classification

敏感性Sp，也稱為真正類率，是C1中所有點被正確預測的比例，計算公式[32]如下：

特異性(Xp)，也稱為真負類率，計算公式[33]如下：

AUC 值計算公式[33]如下：

3 沉陷災害發育敏感性分區結果與討論

3.1 結果

本文按照圖3 的技術流程開展研究區沉陷災害發育敏感性分區預測。

訓練樣本和驗證樣本分別隨機生成同等數量的非沉陷災害樣本，建立由79(沉陷災害)+79(非沉陷災害)個樣本點組成的訓練樣本數據集，共有訓練樣本點158 個。建立由34(沉陷災害)+34(非沉陷災害)個樣本點組成的驗證樣本數據集，共有驗證樣本點68 個。其中，沉陷災害樣本點與非沉陷災害樣本點之間的距離為500 m。

3.1.1 評價因子權重

利用SPSS 軟件統計4 種常用核函數SVM 模型評價因子的重要性排序，結果對比如圖5 所示。

圖5 顯示，4 種模型中評價因子的權重順序并不相同，其權重從高到低的排列順序分別為：

圖5 研究區不同核函數SVM 模型沉陷災害敏感性評價因子權重對比Fig.5 Comparison of the weight of subsidence disaster sensitivity evaluation factors in different kernel function SVM models in the study area

(1) LN-SVM 模型：地層巖組，地質構造，高程，坡向，地面曲率，坡度，地勢起伏度。

(2) PL-SVM 模型：地層巖組，坡向，地質構造，坡度，地面曲率，地勢起伏度，高程。

(3) RBF-SVM 模型：地層巖組，地質構造，坡向，坡度，地面曲率，高程，地勢起伏度。

(4) SIG-SVM 模型：地層巖組，地面曲率，坡度，地勢起伏度，坡向，高程，地質構造。

3.1.2 SVM 模型優選

4 種核函數SVM 模型預測正確率在SPSS 軟件中的統計結果見表5。

由表5 可以看出，4 種核函數SVM 模型的預測正確率由高到低排列為：PL-SVM 模型、LN-SVM 模型、RBF-SVM 模型、SIG-SVM 模型，其中，PL-SVM 模型表現最優。

表5 研究區不同核函數SVM 模型預測正確率Table 5 Prediction accuracy of SVM models with different kernel functions in the study area

3.1.3 評價結果

通過因子權重分析及預測正確率對比，PL-SVM模型的因子重要性排序更合理，預測正確率最高，明顯優于其他3 種模型。因此，本文采用PL-SVM 模型進行研究區沉陷災害發育敏感性分區預測，并在ArcGIS 中采用自然間斷點法將研究區劃分為4 個敏感性等級，分別為：極高敏感區、高敏感區、中等敏感區和低敏感區，得到研究區沉陷災害敏感性分區結果(圖6a)，驗證結果如圖6b 所示，對訓練樣本和驗證樣本的分區統計結果見表6。

圖6 研究區沉陷災害敏感性分區與驗證Fig.6 Sensitivity zoning and verification of subsidence disaster in the study area

表6 研究區沉陷災害敏感性分區統計Table 6 Statistics on sensitivity zoning of subsidence disaster in the study area

由表6 可以看出，PL-SVM 模型通過訓練樣本和驗證樣本分別統計得到的頻率比值(頻率比=災害點比例/分區柵格比例)均隨敏感性等級升高而遞增，經對比其趨勢均符合線性函數，擬合方程分別為：

(1) 訓練樣本頻率比趨勢：y=-0.791x+3.24，R2=0.999 7。

(2) 驗證樣本頻率比趨勢：y=-0.514x+2.47，R2=0.988 7。

其中，y為頻率比值；x為敏感性等級，極高為1，高為2，中為3，低為4。

可以看出，2 個擬合方程的決定系數R2的值均接近于1，表明PL-SVM 模型通過訓練樣本和驗證樣本所計算的頻率比值與敏感性等級之間呈良好的正相關，均以更小的極高與高敏感性分區面積分布了更多的沉陷災害。訓練樣本頻率比的趨勢較驗證樣本頻率比的趨勢相關性更好，表明2012 年與2014 年發育的沉陷災害在研究區內的空間分布略有不同。

3.1.4 精度評價

利用式(3) 分別計算訓練樣本與驗證樣本在PLSVM 模型敏感性分區預測結果中的AUC 值，其中，以2014 年沉陷災害作為驗證樣本的PL-SVM 預測模型的AUC 值為0.755，低于以2012 年沉陷災害作為訓練樣本的PL-SVM 預測模型的AUC 值(0.854)，表明PL-SVM 模型在驗證樣本集的預測精度略低于訓練樣本集，符合預測模型在測試數據集上的精度一般低于訓練數據集精度的規律。但2 種模型預測精度均較高，訓練模型精度等級為非常好，預測模型精度等級為好。預測精度表明，模型在訓練樣本與驗證樣本上既沒有過擬合，也沒有欠擬合，結果較為合理。

3.2 討論

3.2.1 評價因子權重

由圖5 可以看出，地層巖組在4 種模型中均是最重要的評價因子，表明沉陷災害發育與地層巖組分布關系密切，是地面塌陷與地裂縫發育的最重要影響因素。盡管在SIG-SVM 模型中地質構造表現為影響最弱的因子，但從4 種模型中前3 位因子出現的頻率可以看出，地質構造對沉陷災害發育的影響僅次于地層巖組，表明在構造斷裂附近，巖層破碎，地下采煤擾動再次打破斷裂附近的巖體力學平衡，誘發地表塌陷或開裂。

由圖5 可知，地勢起伏度因子在LN-SVM 和RBFSVM 模型中均是沉陷災害發育影響最弱的因子，在PL-SVM 模型中為次最弱因子，表明區域范圍內高程最大值與最小值之間的差值大小對沉陷災害發育影響較小。一般情況下，當沉陷災害處于溝谷邊緣時極易誘發崩塌與滑坡災害，而沉陷災害僅包括地面塌陷與地裂縫，因此，地勢起伏度因子在敏感性評價中貢獻較小。

由圖5 同時可以看出，LN-SVM 和RBF-SVM 模型中地層巖組的權重遠超其他因子，SIG -SVM 模型中各因子的權重相對平衡，地層巖組略高于其他因子，而PL-SVM 模型中地層巖組與坡向因子的權重明顯高于其他因子，在評價因子組合分析中，既沒有過分依賴某一個因子，也沒有分散因子權重，其合理性高于其他3 種核函數模型。

3.2.2 預測模型優選

SVM 預測模型的性能優劣主要取決于核函數的選擇，本文分別對比了LN-SVM 模型、PL-SVM 模型、RBF-SVM 模型和SIG-SVM 模型在研究區沉陷災害敏感性分區中的預測精度，預測正確率分別為75.32%、92.41%、72.15%、62.66%，其中，PL-SVM 模型預測正確率明顯高于其他3 種核函數SVM 模型，分別比LNSVM 模型、RBF-SVM 模型和SIG-SVM 模型的預測正確率高22.69%、28.08%和47.48%。SIG-SVM 模型的預測正確率是4 種模型中最低的，分別比LN-SVM模型、PL-SVM 模型和RBF -SVM 模型的預測正確率低16.81%、32.19%和13.15%。

由圖5 可知，PL-SVM 模型的評價因子權重排序合理性最好，而SIG-SVM 模型則最差，與預測正確率完全一致，選擇PL-SVM 模型作為研究區沉陷災害發育敏感性分區預測模型是可靠的。

3.2.3 沉陷災害敏感性分區預測結果

敏感性評價研究中通常采用總樣本的70%或80%作為訓練點建立預測模型，剩余30% 或20%作為驗證點驗證預測模型精度[30-31]。本文以研究區2012 年和2014 年實際發育的沉陷災害分別作為訓練樣本和驗證樣本，基于PL-SVM 模型進行沉陷災害發育敏感性分區預測。結果表明，2012 年沉陷災害發育在極高、高、中、低4 個等級的數量百分比分別為49.37%、29.11%、17.72%、3.80%，極高和高敏感區分布了共78.48%的沉陷災害。2014 年沉陷災害發育在極高、高、中、低4 個等級的數量百分比分別為38.24%、26.47%、20.59%、14.71%，極高和高敏感區分布了共64.71%的沉陷災害。

由圖6 可知，2012 年沉陷災害在研究區內有明顯的聚集特征，2014 年則演變為隨機分布，部分沉陷災害出現在低敏感區，導致敏感性分區預測結果對驗證樣本(AUC=0.755)的精度略低于訓練樣本(AUC=0.854)，但圖6b 可以看出，模型的預測結果良好，2014 年新發育的沉陷災害大部分位于極高和高敏感區，為沉陷災害核查提供了重點區預測，可極大減少野外調查的盲目性與野外工作量。

3.2.4 綜合分析

目前，地下采煤區地質災害敏感性評價研究中一般將多年災害點混合作為樣本進行模型訓練，但忽視了受采煤擾動影響的礦區地貌特征和土地利用時空演變[26-27]，采用特定年份的災害數據進行建模又需要面對樣本數量過少對敏感性評價結果所帶來的不確定性影響。今后應就此問題繼續展開研究，利用遙感解譯手段增加特定年份災害樣本數量，提升區域災害點小樣本情況下沉陷災害發育敏感性分區預測精度。

研究區含煤地層出露，極高和高敏感區主要分布于王封煤業、西銘煤礦、杜兒坪煤礦、西峪煤礦、白家莊煤礦等多個大型煤礦區域，相比較地形地貌與地質因素，大、小礦地下采煤擾動是研究區沉陷災害發育的主控因素。本文基于方法適用性考慮，未將地下采礦擾動影響列為敏感性評價指標，今后應繼續研究如何定量化引入地下采礦相關指標參與沉陷災害敏感性評價。

太原市西山地區近年所開展的采煤沉陷區治理效果明顯，玉泉山城郊森林公園、長風城郊森林公園、四達溝生態恢復景區等分布于研究區東緣，生態地質環境恢復較好，但沉陷災害極高和高敏感區與植被覆蓋區重疊度較高，對區內地下水資源保護形成潛在威脅，應防止水資源破壞導致植被衰退而引起礦區生態環境再次惡化[28]，也是礦山生態環境治理需要重點關注的區域。同時，分布于林地內的沉陷災害點給遙感解譯和野外調查帶來困難，災害數據集完備性受限，給敏感性分區評價結果的精準性與可靠性帶來一定的影響。

4 結論

a.PL-SVM 模型是研究區沉陷災害敏感性分區預測的最優模型，分別以2012 年和2014 年實際發育沉陷災害點為訓練模型和驗證模型，其訓練精度與驗證精度均較高，且訓練精度較驗證精度高，對地下采煤區沉陷災害發育敏感性分區具有良好的適用性和預測性。

b.PL-SVM 模型劃分了研究區沉陷災害極高、高、中和低敏感區4 個等級的面積比例，極高與高敏感區主要分布于萬柏林區和晉源區邊山一帶，與研究區內大礦與小礦的重合帶高度吻合。災害點比例隨敏感性等級升高而遞增，與敏感性等級之間呈現良好的線性正相關。分區結果合理，對研究區沉陷災害普查與防治重點區確定具有指導價值，對開展同類研究具有參考意義。

c.應進一步改進、完善評價因子組合，并選擇更多評價模型開展對比，組建更可靠的沉陷災害發育敏感性預測體系，同時收集研究區最新的沉陷災害數據進行模型驗證，為地下采煤區沉陷災害發育敏感性分區預測提供最適用的結果。