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

（1.西北大學(xué) 地質(zhì)學(xué)系，陜西 西安 710069；2.大陸動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710069）

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

1 礦區(qū)塌陷預(yù)測(cè)評(píng)價(jià)模型

BP神經(jīng)網(wǎng)絡(luò)算法是一種信息處理體系，利用神經(jīng)網(wǎng)絡(luò)能夠進(jìn)行復(fù)雜的邏輯運(yùn)算，實(shí)現(xiàn)多因子相互作用的非線性關(guān)系。BP神經(jīng)網(wǎng)絡(luò)由輸入層、隱含層和輸出層3部分組成[10]。網(wǎng)絡(luò)訓(xùn)練過(guò)程主要信號(hào)的正向傳播和誤差反向傳播[11]，通過(guò)BP算法的迭代依次對(duì)隱含層到輸出層、輸入層到隱含層的權(quán)重和偏置進(jìn)行調(diào)節(jié)，使權(quán)重系數(shù)和閾值參數(shù)逐漸收斂于合適的值，當(dāng)相鄰2次誤差值達(dá)到目標(biāo)精度時(shí)，模型趨于合理，訓(xùn)練完成。

1.1 信號(hào)的正向傳播

設(shè)xi為人工神經(jīng)網(wǎng)絡(luò)隱含層的輸入值，則隱含層輸出值Hk為：

式中：n為輸入節(jié)點(diǎn)數(shù)；wik為隱含層間連接權(quán)值；ak為隱含層閾值；g為隱含層激勵(lì)函數(shù)。

激活函數(shù)的選取一般認(rèn)為理想的激活函數(shù)是階躍函數(shù)，但階躍函數(shù)具有不連續(xù)、不可導(dǎo)的缺點(diǎn)。塌陷預(yù)測(cè)采用Sigmoid函數(shù)：

Sigmoid函數(shù)的優(yōu)勢(shì)在于其不同于階躍函數(shù)的二值化輸出結(jié)果，其輸出值為（0，1）區(qū)間內(nèi)的可控?cái)?shù)值，該數(shù)值可以用于表示塌陷易發(fā)性的大小。

1.2 誤差反向傳播

通過(guò)連接權(quán)值和偏置的更新使誤差函數(shù)取值逐漸減小，當(dāng)誤差取最小值時(shí)，認(rèn)為這時(shí)的權(quán)重參數(shù)趨于最優(yōu)解[12]。BP神經(jīng)網(wǎng)絡(luò)算法基于梯度下降法來(lái)求解修正權(quán)值。進(jìn)行參數(shù)的傳遞分析可以發(fā)現(xiàn)存在這樣的影響傳遞鏈條關(guān)系，即wkj影響輸出層輸出值，最終影響誤差的大小。反向傳播中權(quán)值更新公式：

式中：wkj為連接權(quán)值；ej為網(wǎng)絡(luò)預(yù)測(cè)誤差；p為輸入節(jié)點(diǎn)數(shù)；η為學(xué)習(xí)速率。

當(dāng)相鄰2次之間的誤差值小于目標(biāo)值時(shí)認(rèn)為算法收斂，算法迭代結(jié)束[14]。

2 礦區(qū)概況及塌陷影響因子分析

2.1 礦區(qū)地質(zhì)概況

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

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

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

2.2 塌陷發(fā)育情況

地質(zhì)災(zāi)害類型主要為滑坡、崩塌、塌陷和泥石流，其中因煤礦開(kāi)采，該地區(qū)塌陷顯著增多。截至2019年10月，研究區(qū)塌陷已增至10余處，目前塌陷有形成貫通的趨勢(shì)。其中塌陷面積最大達(dá)到3 km2，嚴(yán)重影響到礦區(qū)煤礦的繼續(xù)開(kāi)采以及當(dāng)?shù)鼐用袢粘Ｉ睿V區(qū)塌陷分布如圖1。

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

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

2.3 塌陷影響因子選取

2.3.1 數(shù)據(jù)源

研究區(qū)塌陷預(yù)測(cè)評(píng)價(jià)建模的數(shù)據(jù)源主要有：①塌陷編錄信息和野外調(diào)查獲取的相關(guān)資料，包括研究區(qū)1∶2 000井田勘查報(bào)告、礦井建設(shè)工程地質(zhì)災(zāi)害易發(fā)性評(píng)估報(bào)告、礦井設(shè)計(jì)書(shū)等；②DEM高程數(shù)據(jù)，用于提取地形地貌等基礎(chǔ)環(huán)境因子；③研究區(qū)地質(zhì)調(diào)查資料。

2.3.2 影響因子

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

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

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

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

結(jié)合研究區(qū)的實(shí)際情況，本研究選取研究區(qū)3個(gè)基礎(chǔ)環(huán)境因子和5個(gè)工程影響因子，其中基礎(chǔ)環(huán)境因子包括水文特征、地層構(gòu)造和覆巖強(qiáng)度[17]，工程影響因子包括終采時(shí)間、頂板跨度、開(kāi)采深度、采煤高度以及空間迭置層數(shù)。礦區(qū)塌陷影響因子量化方法及易發(fā)性分級(jí)見(jiàn)表1。

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

2.4 塌陷易發(fā)性單元?jiǎng)澐?/h3>

利用預(yù)測(cè)模型開(kāi)展塌陷區(qū)易發(fā)性預(yù)測(cè)的1個(gè)關(guān)鍵問(wèn)題是如何進(jìn)行預(yù)測(cè)單元的劃分[18]。常用的預(yù)測(cè)單元多為柵格[19]、同期開(kāi)采區(qū)域、行政區(qū)劃單元等。其中，柵格具有表達(dá)目標(biāo)明確，模型計(jì)算效率高，數(shù)據(jù)后處理易操作等優(yōu)點(diǎn)[20]。GIS系統(tǒng)具有范圍廣闊且功能強(qiáng)大的空間分析能力，使用其核心工具集可以將數(shù)據(jù)快速處理為柵格化，為數(shù)據(jù)空間分析提供了豐富、便利的環(huán)境和算子[21]。

研究利用GIS的數(shù)據(jù)管理將離散型、連續(xù)型數(shù)據(jù)整合為面數(shù)據(jù)，利用GIS的空間分析功能對(duì)各影響因子進(jìn)行柵格化，根據(jù)礦區(qū)現(xiàn)有資料及GIS系統(tǒng)有效精度選擇適宜的網(wǎng)格大小為20 m×20 m，將48.9 km2研究區(qū)按大小劃分為119 762個(gè)單元，通過(guò)最鄰近分配法分別對(duì)各影響因子單元進(jìn)行指標(biāo)采樣，采樣后的柵格數(shù)據(jù)中每1個(gè)像元點(diǎn)為包含8個(gè)塌陷影響因子的樣本組。礦區(qū)的8個(gè)影響因子分區(qū)如圖3。

圖3 影響因子分區(qū)圖Fig.3 Impact factor classification

3 礦區(qū)塌陷預(yù)測(cè)

采用GIS系統(tǒng)對(duì)礦區(qū)影響因子提取466×257個(gè)單元，對(duì)某些數(shù)據(jù)偏離平均值過(guò)大的畸點(diǎn)取相鄰柵格單元平均值，同時(shí)為防止由于數(shù)據(jù)數(shù)量級(jí)的不統(tǒng)一，對(duì)各類評(píng)價(jià)指標(biāo)數(shù)據(jù)進(jìn)行歸一化處理[20]。通過(guò)Matlab軟件利用礦區(qū)中部及南部466×173個(gè)單元作為樣本進(jìn)行訓(xùn)練構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型，并保留礦區(qū)北部466×84個(gè)單元檢驗(yàn)預(yù)測(cè)結(jié)果。

3.1 神經(jīng)網(wǎng)絡(luò)模型塌陷預(yù)測(cè)

神經(jīng)網(wǎng)絡(luò)隱含層的確定是建模的1個(gè)關(guān)鍵問(wèn)題，隱含層節(jié)點(diǎn)過(guò)少，網(wǎng)絡(luò)難以識(shí)別樣本，從而無(wú)法建立復(fù)雜的映射關(guān)系，導(dǎo)致模型訓(xùn)練能力降低，預(yù)測(cè)效果差；隱含層節(jié)點(diǎn)過(guò)多則會(huì)導(dǎo)致網(wǎng)絡(luò)的過(guò)擬合以及訓(xùn)練效率的降低。目前常用的方法是通過(guò)多次實(shí)驗(yàn)以確定隱含層的最優(yōu)節(jié)點(diǎn)數(shù)[22]。參考式（5）～式（7）進(jìn)行隱含層節(jié)點(diǎn)選?。?/p>

式中：l為隱含層節(jié)點(diǎn)數(shù)；n為輸入層節(jié)點(diǎn)數(shù)；m為輸出層節(jié)點(diǎn)數(shù)；a為0～10之間的常數(shù)。

最終確定隱含層節(jié)點(diǎn)數(shù)為5。

利用上文生成的輸入、輸出訓(xùn)練樣本在Matlab環(huán)境下進(jìn)行神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練和檢驗(yàn)工作。訓(xùn)練使用newff神經(jīng)網(wǎng)絡(luò)函數(shù)。為保證預(yù)測(cè)誤差較小，隱含層節(jié)點(diǎn)傳遞函數(shù)選用S型函數(shù)tansig函數(shù)，輸出層節(jié)點(diǎn)傳遞選用purelin函數(shù)，控制學(xué)習(xí)速率為0.1，選擇Levenberg_Marquardt的BP算法訓(xùn)練函數(shù)trainlm，性能分析選擇均方差性能分析函數(shù)mse，誤差允許限度設(shè)置為1×10-2。訓(xùn)練過(guò)程發(fā)現(xiàn)迭代次數(shù)為2 253時(shí)誤差達(dá)到預(yù)期目標(biāo)，可以認(rèn)為此時(shí)模型訓(xùn)練完成。

將研究區(qū)柵格化數(shù)據(jù)組代入訓(xùn)練完成的神經(jīng)網(wǎng)絡(luò)模型，將模型輸出結(jié)果使用image函數(shù)進(jìn)行圖像化，礦區(qū)塌陷預(yù)測(cè)分區(qū)圖如圖4。

圖4 礦區(qū)塌陷預(yù)測(cè)分區(qū)圖Fig.4 Mine collapse prediction classification

3.2 塌陷預(yù)測(cè)結(jié)果檢驗(yàn)

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

為檢驗(yàn)BP神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)結(jié)果，根據(jù)礦區(qū)實(shí)際發(fā)生的地面塌陷作為檢驗(yàn)標(biāo)準(zhǔn)，采用ROC曲線完成神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型的結(jié)果驗(yàn)證[23]。大量研究表明，ROC曲線可以作為驗(yàn)證預(yù)測(cè)模型的良好指標(biāo)[24]。選取礦區(qū)北部466×84個(gè)單元采用繪制ROC曲線，計(jì)算得到曲線下面積AUC=0.902，預(yù)測(cè)精度可達(dá)90.2%，BP模型ROC曲線測(cè)試圖如圖5。

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

驗(yàn)證結(jié)果表明：采用BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型對(duì)礦區(qū)地下開(kāi)采引起的地面塌陷可以很好地達(dá)到預(yù)期的預(yù)測(cè)結(jié)果。

4 結(jié)語(yǔ)

1）采礦區(qū)是塌陷地質(zhì)災(zāi)害高發(fā)區(qū)之一，通過(guò)選取3個(gè)基礎(chǔ)環(huán)境因子和5個(gè)工程影響因子構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型，BP神經(jīng)網(wǎng)絡(luò)具有很強(qiáng)的主動(dòng)學(xué)習(xí)、適應(yīng)能力，反映出礦區(qū)地面塌陷發(fā)育過(guò)程中各種復(fù)雜的非線性關(guān)系。

2）結(jié)合GIS技術(shù)將礦區(qū)數(shù)據(jù)柵格化和人工神經(jīng)網(wǎng)絡(luò)模型應(yīng)用于塌陷易發(fā)性分區(qū)，可以很好地滿足預(yù)測(cè)精度，同時(shí)在訓(xùn)練計(jì)算過(guò)程中的迭代次數(shù)和耗時(shí)較少，可以達(dá)到快速、優(yōu)良的收斂效果。

3）將預(yù)測(cè)模型與實(shí)際調(diào)查資料進(jìn)行對(duì)比，BP神經(jīng)網(wǎng)絡(luò)模型用于礦區(qū)易發(fā)性預(yù)測(cè)的效果較好，適用于礦區(qū)開(kāi)采結(jié)束后地面塌陷的預(yù)測(cè)，為環(huán)境保護(hù)、土地復(fù)墾方案以及礦區(qū)工程建設(shè)提供可行性、科學(xué)性及準(zhǔn)確性的防治建議。

4）目前建立塌陷預(yù)測(cè)模型所采集的塌陷實(shí)際發(fā)生樣本均發(fā)生在終采6年之內(nèi)，時(shí)間尺度較短，采空區(qū)頂板的殘余變形、破裂并未完全停止，后期深入研究可根據(jù)時(shí)間尺度更長(zhǎng)的塌陷樣本完成礦區(qū)地面塌陷預(yù)測(cè)模型。