基于數值積分法的地下開采地表變形預測算法

孫士國，文小艷，胡業浩

(1. 華北科技學院，河北 三河065201； 2. 臨沂大學機械工程學院，山東 臨沂276005)

1 引言

隨著國民經濟的高速發展，各種能源逐漸被消耗，其中，礦產資源消耗量較高，合理開采與利用礦產資源，對社會與人類未來發展具有重大意義。由于大量的礦山開采活動，開采區域范圍內與礦區內部的穩定性遭到破壞，地表沉陷問題日漸發生，巖層存在移動、變形等情況，隨著開采范圍逐漸擴大，巖層移動變形情況加劇，最終出現地表塌陷、崩塌等嚴重問題。

所謂的地下開采地表變形，即為地下資源被開采結束后，開采作業面與采空區邊緣的原始應力平衡均被打破，因此，巖體或地表均會出現移動變形等問題，最終反映至地表，即構成了地表變形。在19世紀時，就出現過地表變形與沉陷最終使建筑物受到損傷的現象，這種現象由地下開采引起，為此，各開采單位開始重視地下開采造成的負面影響。由于地表變形的影響十分復雜，且不同地質與地下水分布均會導致變形情況的不同，所以，選取適當的預測方法，才能使地表變形預測更加精準。許多學者對變形預測方法進行研究，例如童凱等，研究基于鹽巖腔體收斂過程的地表變形動態預測算法，但該算法預測不同開采深度時的地表水平位移誤差較大。同時還有鐘國強等，研究基于SFLA-GRNN模型的基坑地表最大沉降預測算法，但該算法針對相應開采寬度情況下地表豎向位移預測精度不高，因此，本文研究基于數值積分法的地下開采地表變形預測算法。

數值積分法是現階段礦產開采過程中應用較為廣泛的預測算法。該算法的應用能夠使地表沉陷變形得到有效預測，并改善地下開采對地表變形的影響。本文依據數值積分法的基本原理，構建地表變形預測函數，使地表下沉值與地表相關變形值得到計算，通過該算法旋轉坐標系與處理拐點偏距，實現地表變形預測，并通過預測程序，輸出最終的預測結果。

2 基于數值積分法的地下開采地表變形預測算法

2.1 數值積分法基本原理

在隨機介質原理的基礎上，構建新型的開采地表變形預測函數，即為數值積分法。該算法是較為成熟的變形預測算法，通過式(1)計算地表下沉值

(1)

式(1)中，坐標(，)的下沉量由(，)表示；地表最大下沉值由表示；地下開采范圍由表示；開采影響傳播角由表示，表示下沉量，、分別表示地下開采寬度與深度，表示某段開采巖石。

通過式(2)計算地表某一點沿方向與方向的傾斜

(2)

式(2)中，(，)、(，)分別表示地表點(，)隨著，方向的傾斜，?表示所沿方向。

通過式(3)計算地表變形曲率

(3)

式(3)中，(，)、(，)依次表示地表點(，)隨著，方向的曲率，表示曲率值。

通過式(4)計算地表變形水平位移

(4)

式(4)中，(，)、(，)依次表示地表點(，)隨著，方向的水平位移，表示水平位移值。

通過式(5)計算地表變形的水平變形與扭曲值

(5)

式(5)中，(，)、(，)依次表示地表點(，)隨著，方向的水平變形與扭曲值，表示水平變形與扭曲值。

求解地表變形的剪切值，該值可通過下沉與水平位移的求導獲得，通過式(6)表示

(6)

式(6)中，(，)、(，)依次為地表點(，)的扭曲與剪切，表示扭曲值，表示剪切值。同時，該公式依據直角減小的剪切值為正進行設定。通過以上公式，即能夠實現地下開采地表變形各參數的計算。

2.2 基于數值積分法的算法實現

221 坐標系旋轉調整

通過以上數值積分法在計算地表變形參數時，通常采用軸平行于走向主斷面、軸平行于傾向主斷面的(，)坐標系，但現實情況下的坐標數據會出現不同的坐標形式，因此，需在預測地表變形前轉換坐標位置，使作業面中的角點坐標調整為相對坐標，再進行地表變形預測。通過式(7)表示坐標旋轉變換

=′cos+′sin

=′sin-′cos

(7)

式(7)中，新坐標點由、表示，而原始坐標點則由′、′表示；原坐標系′軸逆時針旋轉至新坐標系軸的角度值則由表示，同時還表示走向方向方位角。通過Rotate函數實現作業面走向方向坐標系的轉換，完成坐標系旋轉的調整，為后續地表變差預測提供全面的坐標系數值。

222 拐點偏距處理

當懸臂作用力下引起的拐點出現偏移，偏移的長度即為拐點偏距，但地下開采地表變形情況下，通常存在任意多邊形變化，因此拐點偏距會出現較大影響。所以，根據數值積分法計算得到作業面數據，能夠得出角點沿作業面走向、傾角方向中的拐點偏距、。

通過加減的方式即能夠計算走向方向內的拐點偏距，但傾斜方向存在資源層傾斜影響，因此需將角點的坐標調整為式(8)。

=-cos-(+sin)·cot

(8)

式(8)中，調整后的坐標由表示；調整前的坐標由表示；該角點區域的采深由表示；煤層傾角由表示。

223 地表變形預測

通過將原始作業面調整為計算作業面，實現作業面的分割預計。本文選取剖面線填充算法分割作業面，該算法能夠快速有效實現作業面分割，并將分割后的作業面轉化為一系列的小矩形以供后續計算，之后將小矩形疊加形成地表變形預計值。通過圖1表示預測的步驟。

圖1 地表變形預測流程

2.3 地下開采地表變形預測程序設計

依據數值積分法的計算原理，通過VC6.0設計地下開采地表變形計算分析程序，通過圖2表示該程序設計結果。

圖2 地表變形預測程序圖

在圖2中，前處理單元負責編寫TXT格式的計算文件，輸入地表變形預測所需條件等參數，之后啟動程序，通過數值積分法計算地表變形參數并開始變形預測。預測結束后進入后處理單元，依據所測地表變形數值，保存成為Matlab圖形文件。同時還可依據預測變形結果繪制等值線圖，并通過AutoCAD軟件內的DXF文件形式輸出，還可以采用Matlab三維圖形文件輸出地表變形結果。

3 實例分析

通過仿真形式驗證本文算法的預測性能，并選取文獻[5]基于鹽巖腔體收斂過程的地表變形動態預測算法、文獻[6]基于SFLA-GRNN模型的基坑地表最大沉降預測算法與本文算法進行對比。

選取10個地表測試點，分析應用本文算法后地表變形數據的實測值與預測值誤差，分析結果如表1所示。

表1 預測誤差分析

根據表1可知，經本文算法預測后的下沉值中，僅有位置4、位置9存在較小誤差，且最高誤差僅為2mm；而經預測的曲率值僅在位置4、位置5與位置10存在誤差；在傾斜值的預測中，出現誤差的情況依然相對較小，最高僅為0.002m，由此說明本文算法的預測結果具有較小的誤差。

分析不同開采深度情況下，地表豎向位移曲線，分析應用本文算法后的預測效果，實驗結果如圖3所示。

圖3 不同開采深度豎向位移情況

根據圖3可知，開采深度為20m時，豎向位移的變化要高于40m與60m的變化，而開采深度為60m時，豎向位移的變化最低，說明當開采深度越大，豎向位移的變化越低。同時，通過本文算法預測地表豎向位移變形過程中，本文預測值與豎向位移實測值十分接近，未出現較大變化，因此，本文算法對豎向位移變形預測效果較高。

分析不同開采深度情況下，地表水平位移曲線，分析應用本文算法后的預測效果，實驗結果如圖4所示。

圖4 不同開采深度水平位移情況

根據圖4可知，開采深度越高，地表的水平位移變形波動越小，當開采深度為20m時，水平位移約在0.2m～-0.2m之內波動，當開采深度達到60m時，水平位移約在0.05m～-0.05m之內波動，同時，經本文算法預測后與實際開采情況十分符合，依然未出現較大誤差，因此，本文算法預測地表水平位移時誤差較低。

分析不同開采寬度情況下，本文算法對地表豎向位移的預測情況，分析結果如圖5所示。

圖5 不同開采寬度豎向位移預測情況

根據圖5可知，當開采寬度處于10m時，地表豎向位移變形最大，而開采寬度為30m時，地表豎向位移變形最小。采用本文算法對地表豎向位移進行預測，可知預測結果與實測結果均保持一致，所以，本文算法能夠實現不同開采寬度的豎向位移變形值預測。

分析開采寬度為20m時，三種算法對地表豎向位移的預測情況，分析結果如圖6所示。

圖6 20m開采寬度時不同方法預測結果

根據圖6可知，在處于20m的開采寬度時，實際豎向位移變形約在-0.2m～-0.8m之間波動，當三種算法開始預測后，文獻[6]算法的預測結果與實測值存在較大偏離，該算法預測后，地表豎向位移約在-0.2m～-0.7m之間波動，而文獻[5]算法地表豎向位移變形波動約在-0.2m～-0.65m以內，預測結果明顯與實際結果不符，經本文算法預測后，地表的豎向位移未出現較大差別，不同開采區域的變形波動預測值均與實測值相符，因此，本文算法預測結果較為精準。

分析開采寬度為20m情況下，三種算法對地表水平位移的預測情況，分析結果如圖7所示。

圖7 20m寬度下不同算法水平位移預測結果

根據圖7可知，開采深度為20m時，實際地表水平位移變化約在0.2m～-0.15m之間，而文獻[6]算法的地表水平位移預測值約在-0.02m～0.12m之間波動，而文獻[5]算法的預測結果也與實際結果存在較大誤差，經本文算法預測后，得到的地表水平位移結果與實測結果并未存在較大誤差，因此，本文預測結果較為精準。

4 結論

本文研究基于數值積分法的地下開采地表變形預測算法，通過數值積分法的基本原理計算地表變形的下沉值以及地表下沉后的變化值，并通過該算法調整坐標系位置，實現旋轉坐標系，并完成拐點偏距處理，以調整后的坐標系與拐點偏距以及計算后的地表變形值為基礎，完成地下開采地表變形的預測，并構建地下開采地表變形預測程序，實現地表變形預測結果的輸出。經實驗驗證該算法的預測性能可知，在不同開采深度、開采寬度情況下，本文算法均能夠實現精準的地表豎向、水平位移預測，預測結果要高于其它算法。在未來研究階段，可依據現有研究，對地表變形預測算法繼續優化，實現精度更高的變形預測。