淺埋采空區(qū)頂板應力變形與地表沉降的動態(tài)關系分析

任 敏 朱萬成 賈瀚文 程關文 李斕堃 代 風

（東北大學資源與土木工程學院，遼寧沈陽110819）

目前，地下采空區(qū)引起的地表沉降逐漸成為了一種常見地質(zhì)災害。以鞍鋼弓長嶺露天鐵礦為例，早在日偽時期就進行過地下富礦開采，造成地下多處采空區(qū)。隨著露天采場開采范圍拓展與延伸及深部鐵礦資源的開采和利用，地表沉降嚴重威脅著礦山大型設備的安全，給正常采礦生產(chǎn)帶來了極大困難。為了能夠?qū)崟r掌握弓長嶺露天鐵礦淺埋采空區(qū)頂板巖體的穩(wěn)定性和地表沉降的發(fā)展演化，研究人員在采空區(qū)頂板和露天坑底的地表建立了集錨桿應力計、多點位移計、液壓水準儀為一體的聯(lián)合監(jiān)測系統(tǒng)。因此，基于該聯(lián)合監(jiān)測系統(tǒng)的監(jiān)測信息，分析淺埋采空區(qū)頂板巖體內(nèi)應力、位移的變化與地表沉降之間的影響關系與響應機制，對于多傳感器監(jiān)測信息有效融合和精準的預測預警模型構建意義重大。

為避免地表沉降等災害的發(fā)生，越來越多的監(jiān)測技術與監(jiān)測設備應運而生，為地表沉降監(jiān)測提供了多種多樣的監(jiān)測手段和監(jiān)測方法。目前普遍常用的地表沉降監(jiān)測手段主要有全站儀［1］、水準儀［2］、GPS站點［3］等定點式監(jiān)測方法，以及三維激光掃描［4］、雷達掃描［5］等區(qū)域式監(jiān)測方法。在獲取現(xiàn)場實時監(jiān)測數(shù)據(jù)的前提下，為避免地表沉降變形引起災害及人員、財產(chǎn)損失，開展地表沉降預測研究十分必要。近年來，眾多研究者致力于利用不同的方法手段實現(xiàn)這一目標，主要包括人工智能法、灰色預測法、數(shù)學函數(shù)法等。如李勇發(fā)等［6］基于PS-InSAR 雷達干涉測量技術對地表沉降進行了監(jiān)測，對獲取的實時沉降數(shù)據(jù)利用遺傳算法優(yōu)化的反饋神經(jīng)網(wǎng)絡算法進行了地表沉降預測，獲得了良好的預測效果；張濤［7］將DInSAR 技術與小基線集技術進行有機結合對地表沉降進行了監(jiān)測，并使用粒子群算法優(yōu)化支持向量機模型對沉降數(shù)據(jù)進行了擬合與預測；黨星海等［8］采用卡爾曼濾波方法對三維激光掃描沉降監(jiān)測獲取的初始監(jiān)測數(shù)據(jù)進行了去噪處理，而后利用GM（1，1）灰色預測模型進行了沉降變形預測，得到了優(yōu)于傳統(tǒng)灰色預測模型的預測結果；李篷［9］采用GPS方法進行了煤礦地表非采動沉降監(jiān)測，并利用GM（1，1）灰色預測模型對煤礦主井、副井、煤倉附近區(qū)域開展了地表沉降預測；Nie 等［10-11］基于采空區(qū)上覆巖體不同破裂帶的“S”型變形特征，以及塌陷坑GPS 監(jiān)測獲取的“S”型地表沉降曲線，提出了采空區(qū)的反正切函數(shù)預測模型；袁興明等［12］利用數(shù)字水準儀進行了沉降觀測，并分別利用BP 神經(jīng)網(wǎng)絡算法、改進BP 神經(jīng)網(wǎng)絡算法對沉降數(shù)據(jù)進行了預測。就目前研究來看，基于單一監(jiān)測信息的地表沉降預測存在一定的缺陷，是因為數(shù)據(jù)源本身既無法充分反映地表沉降的多種誘導因素，也無法全面體現(xiàn)沉降過程中的多種響應信息，從而導致預測效果不穩(wěn)定且與實際存在一定的偏差。

為了獲取更全面、精準的監(jiān)測數(shù)據(jù)，從多種信息融合角度判斷監(jiān)測對象的發(fā)展狀態(tài)，采用多種監(jiān)測手段針對地表沉降變形建立聯(lián)合監(jiān)測系統(tǒng)，開展多手段聯(lián)合監(jiān)測已經(jīng)得到廣泛研究。如何倩等［13］采用三維激光掃描與D-InSAR 聯(lián)合監(jiān)測礦區(qū)地表動態(tài)沉降；臧?。?4］將測量機器人與電子水準儀相結合應用于隧道監(jiān)測工程中；郭利民等［15］利用GPS 和PS-In-SAR 聯(lián)合監(jiān)測西秦嶺斷裂北緣地表運動；李楊楊等［16］將GPS-RTK聯(lián)合水準儀應用于充填工作面地表高精度動態(tài)監(jiān)測中。目前，學術界對于多手段聯(lián)合監(jiān)測數(shù)據(jù)的處理及應用方面存在兩種研究方向：一是將多種監(jiān)測技術進行有機結合，得到精度較高的地表變形監(jiān)測結果，如余禮仁［17］在對淮北袁店二礦煤礦區(qū)地表沉降進行監(jiān)測時，結合了D-InSAR 和GIS技術得到監(jiān)測數(shù)據(jù)，而后采用三次樣條插值與灰色預測組合模型GM（1，1）對該礦7255工作面中心區(qū)域的下沉趨勢進行了擬合并預測；邢學敏［18］將CR-In-SAR 與PS-InSAR 相結合，對礦區(qū)時序地表變形進行了監(jiān)測，提高了礦區(qū)地表變形監(jiān)測的精度和可靠性；二是根據(jù)模糊系統(tǒng)理論，將多個傳感器獲取的信息進行融合，進而對地表沉降狀態(tài)進行綜合估計與評價，如Wu 等［19］采用野外調(diào)查、地球物理勘探、鉆井勘探和D-InSAR 技術等聯(lián)合監(jiān)測手段，應用模糊綜合評價法對地表沉降進行了風險分區(qū)。由此可見，對于監(jiān)測數(shù)據(jù)而言，雖然多傳感器聯(lián)合監(jiān)測總體上提高了數(shù)據(jù)樣本的多樣性和精度，但是在基于多傳感器聯(lián)合監(jiān)測結果的預測方法研究與應用中，對多傳感器監(jiān)測獲得的地表變形時序數(shù)據(jù)未能進行有效融合，對于提升預測預警效果作用不顯著。

程關文等［20-21］認為深部巖體變形是形成地表變形區(qū)域的原因，并采用數(shù)值模擬方法分析了隨著不同開挖步的增加，深部巖體變形產(chǎn)生地表沉降變形的整個過程。宋長根等［22］對深部巖體變形對地表沉降變形的影響進行了研究，認為當深部巖體發(fā)生破壞后，地表會進入加速變形階段，地表與構筑物不久便出現(xiàn)裂縫，嚴重影響了礦區(qū)安全生產(chǎn)與構筑物正常使用。本研究以弓長嶺露天鐵礦某空區(qū)頂板為例，根據(jù)錨桿應力計、多點位移計、液壓水準儀對空區(qū)頂板巖體及地表沉降的聯(lián)合監(jiān)測結果，采用自回歸分布滯后模型（ARDL）分析淺埋采空區(qū)頂板巖體內(nèi)的應力、位移變化對地表沉降的影響關系和響應機制，確定該空區(qū)頂板巖體內(nèi)發(fā)生應力、位移變化影響到地表沉降的滯后時間，這種發(fā)展關系的探究不僅可以為多傳感器的信息融合提供有效途徑，還可以為地表沉降災害的預測預警提供支撐。

1 弓長嶺鐵礦某采空區(qū)概況

近年來，弓長嶺露采分公司采用地球物理探測方法研究采空區(qū)的空間分布和大小情況，取得了顯著成效，基本查清了一定深度空區(qū)的分布、大小等信息。但是，僅僅研究采空區(qū)的分布、大小、深度遠遠無法滿足礦山安全生產(chǎn)的需要。因此，結合弓長嶺鐵礦的工程實際情況，建立了一套先進的監(jiān)測系統(tǒng)，包括錨桿測力計、多點位移計、液體靜力水準儀等，實現(xiàn)對地表沉降變形和采空區(qū)頂板巖體的實時在線監(jiān)測，掌握空區(qū)頂板巖體內(nèi)應力、位移變化及地表沉降的發(fā)展，以確定機械設備、車輛能否安全通過，對于露天礦安全開采至關重要。系統(tǒng)測點布置如圖1所示。

本研究選擇2018-06-11—2018-07-01獲取的時間序列監(jiān)測數(shù)據(jù)進行分析，測點為研究區(qū)內(nèi)（圖1）地表沉降液壓水準儀測點S4、A 孔多點位移計位于135.2 m 水平的測點uA4和位于129.2 m 水平的測點uA1，A 孔錨桿測力計位于138.7 m 水平的測點σA3和位于131.7 m 水平的測點σA1。選擇該類測點進行分析，原因在于：①3 種傳感器測點基本處于同一地質(zhì)環(huán)境，上覆巖體中積聚產(chǎn)生的應力和變形會以較單一的形式影響到地表；②各個測點空間距離較近，且在垂直方向上多點位移計測點與錨桿測力計測點在采空區(qū)頂板巖體中交錯布置，基本能夠監(jiān)測到巖體內(nèi)發(fā)生的應力和變形進而影響地表產(chǎn)生沉降的全過程。因此，可根據(jù)該類測點的數(shù)據(jù)，研究上覆巖體內(nèi)的應力、位移變化對地表沉降的影響機制。

2 模型構建及實證分析

2.1 ARDL模型設定

本研究采用自回歸分布滯后模型（ARDL）分析由錨桿應力計、多點位移計監(jiān)測的淺埋采空區(qū)頂板巖體內(nèi)的應力、位移對液壓水準儀監(jiān)測獲得的地表沉降之間的影響關系、響應機制和滯后時間。Pesaran等［23］認為ADRL 模型的優(yōu)勢在于不需要變量同階單整，也可以檢驗變量之間的長期關系，當解釋變量為內(nèi)生變量時，ARDL 模型的估計也不會受到影響，并且在小樣本情況下ARDL的估計結果也更加穩(wěn)定。基于此，本研究構建的地表沉降與頂板巖體應力、位移之間的ARDL模型為

式中，S 為地表沉降，為被解釋變量；u 和σ 分別為采空區(qū)頂板巖體的位移和應力，為解釋變量；q1、q2、q3分別為各變量對地表沉降長期影響的滯后期數(shù)；α1、α2、α3為各變量及其滯后期數(shù)對被解釋變量的影響程度，若解釋變量的滯后期t 值顯著，則該變量的滯后期對地表沉降有顯著影響；β0為常數(shù)項；ut為殘差項。

2.2 基于ARDL模型實證分析

2.2.1 單位根檢驗

為了避免在實證檢驗過程中發(fā)生偽回歸，需要首先對數(shù)據(jù)進行單位根檢驗。盡管ARDL 模型不要求變量具有同階單整，但應用該模型的前提條件需確保變量單整階數(shù)不能超過1。本研究采用最常用的ADF 檢驗統(tǒng)計量對模型變量進行單位根檢驗，結果如表1所示。

注：單位根檢驗形式（C，T，L）中，C為截距項，T為趨勢項，L為滯后階數(shù)。

由表1可知，在5%的顯著性水平下，包括采空區(qū)頂板巖體的位移、應力以及地表沉降在內(nèi)的所有變量均為非平穩(wěn)序列。對這3 個變量進行一階差分后的ADF 檢驗，在5%的顯著性水平上都是平穩(wěn)序列，說明這3 個變量均是一階單整（即I(1)）序列，符合ARDL模型的應用要求。

2.2.2 協(xié)整檢驗

ARDL 模型的第一步是檢驗變量的協(xié)整關系。由于協(xié)整關系對滯后階數(shù)的選擇非常敏感，因此，需要估計向量自回歸模型（Vector Autoregression，VAR），確定滯后階數(shù)。通過信息準則和似然比（LR）檢驗法得到VAR 模型滯后階數(shù)的確定結果如表2 所示。在滯后階數(shù)為1 時，最終預測誤差FPE 和3 個信息準則AIC、SC、HQ 均最小，由4 個參考值確定1 階滯后為最優(yōu)，因此可得滯后階數(shù)為1。

注：LogL 表示對數(shù)似然函數(shù)；表中用“*”表示從每一列標準中選的滯后階數(shù)；

滯后階數(shù)確定后，為方便起見，本研究采用Johansen 協(xié)整檢驗方法進行協(xié)整檢驗，結果如表3 所示。

?

由表3可知：地表沉降S4與采空區(qū)頂板巖體內(nèi)的位移uA1、uA4和應力σA1、σA3之間跡統(tǒng)計量和最大特征值統(tǒng)計量的檢驗結果均存在協(xié)整關系，因而可得出地表沉降與淺埋采空區(qū)頂板巖層中的位移、應力之間存在長期關系。

2.2.3 長期均衡系數(shù)估計

確定各變量之間存在協(xié)整關系后，可以用ARDL模型進行長期關系和短期關系的估計。在估計長期關系的變量系數(shù)之前需要確定模型中的滯后階數(shù)，本研究使用SC準則確定模型中各變量的最優(yōu)滯后階數(shù)。估計結果顯示：地表沉降S4與位移uA1、uA4，應力σA1、σA3以ARDL( 1,1 ,1,1,1 )模型最為合適，長期系數(shù)如表4所示。

注：D 表示差分，是因為單位根檢驗已得到的各變量均為一階差分平穩(wěn)，故所有變量均采用一階差分形式；表示一階差分后的各變量由ARDL 模型選擇的滯后一階形式表示一階差分的原變量；C 為截距項；@TREND 為時間趨勢項；AIC 為赤池AIC信息準則（Akaike Information Criterion），SC為施瓦茨SIC信息準側（Schwarz Criterion），兩者均是衡量回歸模型優(yōu)良性的標準，數(shù)值較小意味著滯后階數(shù)較為合適；R2 統(tǒng)計量用于衡量模型的估計值對實際值擬合的好壞，數(shù)值越大，表明模型對因變量的擬合效果越好。

由表4 可知：淺埋采空區(qū)頂板巖體內(nèi)的應力、位移滯后量對地表沉降影響顯著，說明上覆巖體內(nèi)的應力、變形影響傳遞到地表需要時間，以使得滯后量與地表沉降最相關。變量系數(shù)有正有負，變量系數(shù)為正，說明采空區(qū)頂板巖體內(nèi)產(chǎn)生的應力、位移使得地表沉降加劇；變量系數(shù)為負則較為復雜，存在水平地應力增加的情形，使得上覆巖體內(nèi)的應力、變形對地表呈現(xiàn)出與塌陷相反方向變形的影響。估計得到的地表沉降S4與位移uA1、uA4，應力σA1、σA3的ARDL( 1,1 ,1,1,1 )模型中，滯后一期的位移uA4、應力σA1和σA3對地表沉降的影響均顯著，即前1 d 采空區(qū)頂板內(nèi)產(chǎn)生的應力、位移與當前發(fā)生的地表沉降最相關，其中應力σA1對地表沉降S4的影響最大，其系數(shù)為0.529 29。

此外，還進行了數(shù)值模擬分析，即保證模擬的弓長嶺露天鐵礦淺埋采空區(qū)頂板巖體所受的豎向應力基本控制在監(jiān)測結果范圍內(nèi)，測點所受應力為0.38～0.74 MPa 時，提取地表測點處的沉降值為0.005 97～0.009 39 m（圖2）。由數(shù)值模擬結果分析可知：σA1的改變對地表沉降的影響較大，系數(shù)達0.572 81 左右，與ARDL 模型的估計結果較吻合。故可將淺埋采空區(qū)頂板巖體內(nèi)的應力、位移作為地表沉降的長期影響指標，度量地表沉降。

2.2.4 短期動態(tài)調(diào)整關系分析

在估計得到長期關系系數(shù)后，可進一步求得基于ARDL模型的誤差修正項（ECM）的系數(shù)，分析各個解釋變量的短期變動對被解釋變量的動態(tài)調(diào)整關系。ARDL的誤差修正模型具體估計結果見表5。

注：EC M( -1 )為滯后一期的誤差修正項。

由表5 可知：在10%的顯著性水平下，滯后一期的地表沉降與滯后一期的位移uA1以及當期的應力σA1對地表沉降均存在顯著影響，其中應力的影響最為顯著，該變量的正向變動會導致地表沉降加劇。由誤差修正項ECM 估計結果表明存在反向調(diào)整機制，意味著當受到短期沖擊后，系統(tǒng)會較快恢復至長期均衡狀態(tài)。

因此，基于上述協(xié)整檢驗、ARDL 模型的長期均衡系數(shù)估計和短期動態(tài)調(diào)整關系的分析，證實了地下采空區(qū)頂板巖體內(nèi)產(chǎn)生的應力、位移與地表沉降之間存在長期和短期影響。

3 結論與建議

在地下開采等多因素的影響下，淺埋采空區(qū)頂板巖體會產(chǎn)生應力、變形，作為直接影響因素導致地表沉降，體現(xiàn)在傳感器上能夠捕捉到的信息即為頂板巖體內(nèi)的應力、位移以及地表沉降。為了給地表沉降多傳感器監(jiān)測信息的有效融合及預測預警模型構建提供支撐，本研究利用ARDL模型系統(tǒng)分析了淺埋采空區(qū)頂板巖體內(nèi)應力、位移的變化和地表沉降之間的影響關系。研究表明：地表沉降與淺埋采空區(qū)頂板巖體內(nèi)的應力、位移之間存在長期和短期的影響；利用SC 準則確定的最優(yōu)模型ARDL( )1,1,1,1,1分析表明，作為地表沉降解釋變量的頂板巖體內(nèi)的應力、位移滯后期均為1，即當前產(chǎn)生的地表沉降與前1 d 頂板巖體內(nèi)產(chǎn)生的應力、位移有關，且滯后的解釋變量的影響均很顯著。

針對地表沉降災害錨桿應力計、多點位移計、液壓水準儀的聯(lián)合監(jiān)測，本研究提出如下建議：已證實淺埋采空區(qū)頂板巖體內(nèi)的應力、位移對地表沉降具有長期影響，且滯后的應力、位移影響顯著，故而可以采用各種預測模型（如人工智能算法等）將地表沉降時間序列數(shù)據(jù)作為輸出值，采空區(qū)頂板巖體內(nèi)滯后的應力、位移時間序列數(shù)據(jù)作為輸入值，根據(jù)監(jiān)測數(shù)據(jù)建立變量之間的線性、非線性關系進行預測?；诙鄠鞲衅鞅O(jiān)測數(shù)據(jù)融合構建的預測模型，可以有效體現(xiàn)出地表沉降過程中淺埋采空區(qū)頂板巖體的響應機理，優(yōu)于具有固定理論模型的數(shù)值模擬方法，因而具有進一步深入研究的價值和應用前景。