基于GIS與BP神經網絡的礦區塌陷易發性預測

王紳皓，謝婉麗，奚家米，井 旭

（1.西北大學 地質學系，陜西 西安 710069；2.大陸動力學國家重點實驗室，陜西 西安 710069）

地面塌陷是礦山常見的地質災害之一，礦區的地面塌陷為采礦工作及當地居民帶來了嚴重的安全隱患及經濟損失[1]。因此，礦區塌陷的預測是保障生產工作正常運行，保護居民生命財產安全的重要工作[2]。針對塌陷易發性的預測，已經有眾多學者進行了大量的研究，并提出了各種判別方法。劉寶琛等[3]提出將概率積分法運用到地下開采等實際工程引發地面塌陷的預測工作中；謝和平等[4]使用相似材料作了巷道和采場的巖層塌陷實驗，并基于非線性的連續介質模型對地下工程的塌陷進行了非線性力學有限元分析。孔憲立等[5]論述了采空區地表為斜坡條件下，上覆巖層發生多種變形、破壞的力學特征；陳學軍[6]、包惠明[7]等分別將GIS技術、模糊數學理論應用在巖溶塌陷評價模型中，對其易發性進行定量評價。以上方法在塌陷易發性預測中均有較好的效果，但由于巖土體的空間特性特殊，建立力學模型的方法復雜且難以簡化；數值模擬難以全面考慮礦山塌陷中的復雜因素，選取合理參數，模糊評判存在超模糊現象。為補充以往研究不足，為此選用BP神經網絡進行礦區塌陷的預測研究。BP神經網絡具有主動學習、多層次，可以通過學習大量的樣本來獲得適合的連接權值，反映出各種復雜的非線性關系[8]。綜合GIS系統空間分析功能在礦區塌陷預測突出的優越性，采用柵格為評價單元，提高預測精度，最大限度地體現出礦區內部塌陷易發性的空間差異[9]，為模型預測提供足夠的樣本數。

1 礦區塌陷預測評價模型

BP神經網絡算法是一種信息處理體系，利用神經網絡能夠進行復雜的邏輯運算，實現多因子相互作用的非線性關系。BP神經網絡由輸入層、隱含層和輸出層3部分組成[10]。網絡訓練過程主要信號的正向傳播和誤差反向傳播[11]，通過BP算法的迭代依次對隱含層到輸出層、輸入層到隱含層的權重和偏置進行調節，使權重系數和閾值參數逐漸收斂于合適的值，當相鄰2次誤差值達到目標精度時，模型趨于合理，訓練完成。

1.1 信號的正向傳播

設xi為人工神經網絡隱含層的輸入值，則隱含層輸出值Hk為：

式中：n為輸入節點數；wik為隱含層間連接權值；ak為隱含層閾值；g為隱含層激勵函數。

激活函數的選取一般認為理想的激活函數是階躍函數，但階躍函數具有不連續、不可導的缺點。塌陷預測采用Sigmoid函數：

Sigmoid函數的優勢在于其不同于階躍函數的二值化輸出結果，其輸出值為（0，1）區間內的可控數值，該數值可以用于表示塌陷易發性的大小。

1.2 誤差反向傳播

通過連接權值和偏置的更新使誤差函數取值逐漸減小，當誤差取最小值時，認為這時的權重參數趨于最優解[12]。BP神經網絡算法基于梯度下降法來求解修正權值。進行參數的傳遞分析可以發現存在這樣的影響傳遞鏈條關系，即wkj影響輸出層輸出值，最終影響誤差的大小。反向傳播中權值更新公式：

式中：wkj為連接權值；ej為網絡預測誤差；p為輸入節點數；η為學習速率。

當相鄰2次之間的誤差值小于目標值時認為算法收斂，算法迭代結束[14]。

2 礦區概況及塌陷影響因子分析

2.1 礦區地質概況

研究礦區地處陜西省廟哈孤礦區東南部，該礦區為典型的黃土梁峁地貌，溝壑縱橫，梁面寬緩平坦，場地內地勢較為平坦。區內地層主要為三迭系細粒石英砂巖，厚度約80～200 m；侏羅系粉砂巖，泥巖夾煤層，厚度約89.78 m，新近系紅色黏土，平均厚度3.1 m；第四系黃土覆蓋大部分山梁、緩坡，平均厚度24.15 m。

礦區降水較少，年內降水量變化大，水系較發育，流經礦區的水系主要為2條：孤山川河自井田西側邊界穿過，府谷縣內流長57 km，流域面積1 018 km2，年平均流量3.48 m3/s，溝谷寬緩；流小南川溝流域面積約6 km2，河流受區域獨特的地質地貌條件及氣象因素的控制，雨季洪峰多。

井田內設計開采煤層為22、31、42、51、52上、52煤，煤層平均可采總厚度12.22 m。以52煤為例，煤層厚度1.47～3.32 m，平均2.23 m，煤厚變化較小，埋深淺，煤層埋藏深度在0～225.66 m，煤層頂板巖性為1套厚層粉砂巖、局部夾薄層粉砂質泥巖，厚16 m。

2.2 塌陷發育情況

地質災害類型主要為滑坡、崩塌、塌陷和泥石流，其中因煤礦開采，該地區塌陷顯著增多。截至2019年10月，研究區塌陷已增至10余處，目前塌陷有形成貫通的趨勢。其中塌陷面積最大達到3 km2，嚴重影響到礦區煤礦的繼續開采以及當地居民日常生活，礦區塌陷分布如圖1。

圖1 礦區塌陷分布圖Fig.1 Distribution map of mining area collapse

以二采區形成的地面塌陷為例，其形成時間為2013—2018年，位于井田中部，走向長度1.91 km，傾斜長度1.06 km，該采區52號煤層平均厚度2.38 m，煤層結構較簡單。該地面塌陷區面積較大，塌陷區東側較寬，向西逐漸變窄，地面裂縫發育，裂縫長約7～50 m。裂縫成組平行發育，間距5～20 m不等，部分裂縫可見深度達1.5 m，兩側錯落10～50 cm，在斜坡地帶局部發生了滑塌，地面塌陷裂縫方向與采取工作面方向相一致。

2.3 塌陷影響因子選取

2.3.1 數據源

研究區塌陷預測評價建模的數據源主要有：①塌陷編錄信息和野外調查獲取的相關資料，包括研究區1∶2 000井田勘查報告、礦井建設工程地質災害易發性評估報告、礦井設計書等；②DEM高程數據，用于提取地形地貌等基礎環境因子；③研究區地質調查資料。

2.3.2 影響因子

礦區地面塌陷大多數是由于頂板的破壞而引起的，即礦體被開挖、采掘后，失去下部巖層支持的頂板在自重和次生應力作用下發生剪切或拉裂破壞[13]，該礦區采用走向長壁后退式綜合機械化采煤方法，全部垮落法管理頂板。研究區垂直剖面由上至下可簡化為松散覆蓋物、堅硬巖層、含水層等，區內分布構造裂隙，采空區垂直剖面如圖2。

圖2 采空區垂直剖面圖Fig.2 Vertical section view of goaf

在煤層回采后，上覆巖體完全垮落，垮落帶可達松散覆蓋物底部，在松散覆蓋層形成裂隙，使松散層本身發生垮落，導致地表塌陷[14]。

影響礦區塌陷的因素多種多樣，各個采區情況具有獨特性[15]，眾多研究表明塌陷的發生是基礎環境和工程影響共同作用的結果。采煤方法和頂板管理方法在該礦區不同工作面中相同，在預測模型中不予考慮。塌陷預測模型易發性因素包括：①覆巖強度是反映巖體工程地質特征的最根本因素，直接影響著巖體的物理力學性質及其受力變形的全過程；②水文地質特征，地表水、地下水對巖土體的作用是一種地質應力作用，水與巖土體相互作用改變著巖土體的物理、化學及力學性質；③地質構造主要包括不連續結構面、構造應力，地質構造活動使得巖體中聚集大量的彈性變形能，產生構造應力，受到開挖擾動的影響，構造應力釋放使巖體中出現局部應力集中現象，誘發地面塌陷[16]；④巖層的變形和破壞是從直接頂開始，自下而上擴散直至地面的漸變過程，終采時間越長，塌陷的易發性越高；⑤采礦工程將巖體中原始應力平衡狀態破壞，并使其周圍一定范圍的巖層應力發生重新分布，而開采深度和跨度是影響著擾動后原巖應力的大小、方向與分布狀態的重要因素；⑥采煤高度影響著塌陷的易發程度，當采煤高度較高或多煤層迭置開采時，覆巖經多次塌陷分級傳至地面造成地表沉陷，產生大量裂縫。

結合研究區的實際情況，本研究選取研究區3個基礎環境因子和5個工程影響因子，其中基礎環境因子包括水文特征、地層構造和覆巖強度[17]，工程影響因子包括終采時間、頂板跨度、開采深度、采煤高度以及空間迭置層數。礦區塌陷影響因子量化方法及易發性分級見表1。

表1 礦區塌陷影響因子量化方法及易發性分級Table 1 Quantification method and stability classification of mining area collapse influence factor

2.4 塌陷易發性單元劃分

利用預測模型開展塌陷區易發性預測的1個關鍵問題是如何進行預測單元的劃分[18]。常用的預測單元多為柵格[19]、同期開采區域、行政區劃單元等。其中，柵格具有表達目標明確，模型計算效率高，數據后處理易操作等優點[20]。GIS系統具有范圍廣闊且功能強大的空間分析能力，使用其核心工具集可以將數據快速處理為柵格化，為數據空間分析提供了豐富、便利的環境和算子[21]。

研究利用GIS的數據管理將離散型、連續型數據整合為面數據，利用GIS的空間分析功能對各影響因子進行柵格化，根據礦區現有資料及GIS系統有效精度選擇適宜的網格大小為20 m×20 m，將48.9 km2研究區按大小劃分為119 762個單元，通過最鄰近分配法分別對各影響因子單元進行指標采樣，采樣后的柵格數據中每1個像元點為包含8個塌陷影響因子的樣本組。礦區的8個影響因子分區如圖3。

圖3 影響因子分區圖Fig.3 Impact factor classification

3 礦區塌陷預測

采用GIS系統對礦區影響因子提取466×257個單元，對某些數據偏離平均值過大的畸點取相鄰柵格單元平均值，同時為防止由于數據數量級的不統一，對各類評價指標數據進行歸一化處理[20]。通過Matlab軟件利用礦區中部及南部466×173個單元作為樣本進行訓練構建BP神經網絡模型，并保留礦區北部466×84個單元檢驗預測結果。

3.1 神經網絡模型塌陷預測

神經網絡隱含層的確定是建模的1個關鍵問題，隱含層節點過少，網絡難以識別樣本，從而無法建立復雜的映射關系，導致模型訓練能力降低，預測效果差；隱含層節點過多則會導致網絡的過擬合以及訓練效率的降低。目前常用的方法是通過多次實驗以確定隱含層的最優節點數[22]。參考式（5）～式（7）進行隱含層節點選取：

式中：l為隱含層節點數；n為輸入層節點數；m為輸出層節點數；a為0～10之間的常數。

最終確定隱含層節點數為5。

利用上文生成的輸入、輸出訓練樣本在Matlab環境下進行神經網絡模型的訓練和檢驗工作。訓練使用newff神經網絡函數。為保證預測誤差較小，隱含層節點傳遞函數選用S型函數tansig函數，輸出層節點傳遞選用purelin函數，控制學習速率為0.1，選擇Levenberg_Marquardt的BP算法訓練函數trainlm，性能分析選擇均方差性能分析函數mse，誤差允許限度設置為1×10-2。訓練過程發現迭代次數為2 253時誤差達到預期目標，可以認為此時模型訓練完成。

將研究區柵格化數據組代入訓練完成的神經網絡模型，將模型輸出結果使用image函數進行圖像化，礦區塌陷預測分區圖如圖4。

圖4 礦區塌陷預測分區圖Fig.4 Mine collapse prediction classification

3.2 塌陷預測結果檢驗

從BP神經網絡模型的預測分區圖來看，預測區塌陷情況整體符合礦區采空塌陷的實際情況：礦區塌陷主要分布于下部有采空的地區，同時塌陷受到地下巷道支護控制，具備煤柱支護的地區塌陷易發性最低，且支護具備一定區域范圍內的塌陷控制能力，將塌陷分割為多個區域且塌陷易發性明顯降低。地層巖性變化較小時，上覆巖層硬質巖層占比越大，采空區塌陷的易發性越低。采空區頂板越厚、跨度越小，越不容易發生離層破壞或頂板折斷，對應區域塌陷易發性越低。塌陷的易發性在一定程度上受到水文特征的影響，其中，中南部礦區內有水系發育的地區塌陷易發性較高。

為檢驗BP神經網絡的預測結果，根據礦區實際發生的地面塌陷作為檢驗標準，采用ROC曲線完成神經網絡預測模型的結果驗證[23]。大量研究表明，ROC曲線可以作為驗證預測模型的良好指標[24]。選取礦區北部466×84個單元采用繪制ROC曲線，計算得到曲線下面積AUC=0.902，預測精度可達90.2%，BP模型ROC曲線測試圖如圖5。

圖5 BP模型ROC曲線測試圖Fig.5 ROC curve test diagram of BP model

驗證結果表明：采用BP神經網絡預測模型對礦區地下開采引起的地面塌陷可以很好地達到預期的預測結果。

4 結語

1）采礦區是塌陷地質災害高發區之一，通過選取3個基礎環境因子和5個工程影響因子構建的BP神經網絡預測模型，BP神經網絡具有很強的主動學習、適應能力，反映出礦區地面塌陷發育過程中各種復雜的非線性關系。

2）結合GIS技術將礦區數據柵格化和人工神經網絡模型應用于塌陷易發性分區，可以很好地滿足預測精度，同時在訓練計算過程中的迭代次數和耗時較少，可以達到快速、優良的收斂效果。

3）將預測模型與實際調查資料進行對比，BP神經網絡模型用于礦區易發性預測的效果較好，適用于礦區開采結束后地面塌陷的預測，為環境保護、土地復墾方案以及礦區工程建設提供可行性、科學性及準確性的防治建議。

4）目前建立塌陷預測模型所采集的塌陷實際發生樣本均發生在終采6年之內，時間尺度較短，采空區頂板的殘余變形、破裂并未完全停止，后期深入研究可根據時間尺度更長的塌陷樣本完成礦區地面塌陷預測模型。