基于Matlab的山區急傾斜煤層開采沉陷預計系統

趙 博 梁乃森 吳 初 武 雄

（中國地質大學（北京）水資源與環境學院，北京100083）

煤層開采引起的地表移動變形與地表傾角大小、地形條件密切相關［1-4］。不同的煤層傾角、地形條件適用于不同的開采沉陷預計方法［5-7］，為計算方便，通常借助不同的計算機編程語言開發出多功能的可視化開采沉陷預計系統［8-10］，也有學者將機器學習算法應用于開采沉陷預計方面［11-12］。近年來，對于急傾斜煤層開采引起的地表移動變形規律和預計方面的研究，成果豐碩［13-16］，但對于山區條件下的急傾斜煤層開采引起的地表移動變形預計方面的研究深度稍有欠缺［17-18］。相關研究表明［7］：皮爾森Ⅲ型公式法適用于平地條件下急傾斜煤層開采地表移動變形預計，山區地表移動預計模型主要適用于山坡平均坡度小于30°的水平和緩傾斜煤層開采條件下的地表移動變形預計。本研究將皮爾森Ⅲ型公式法與山區地表移動預計模型［7］進行疊加，實現優勢互補，在Matlab平臺上開發預計系統，對山區急傾斜煤層開采沉陷進行預計。

1 系統設計及編程實現

本研究開采沉陷預計系統分為五大模塊，即坐標轉換、山地地表曲面擬合、地面移動變形計算、地面移動變形繪圖、地表任意點移動變形值計算，通過代碼編寫、整合及界面開發，便于人機交互操作，實現數據的快速處理與分析。

1.1 計算原理

本研究系統開發是建立在對皮爾森Ⅲ型公式法和山區地表移動預計模型疊加的基礎上，因此實現兩者參數意義匹配較為重要。通過對公式分析可知，兩者參數匹配不存在重大矛盾，皮爾森Ⅲ型公式法可直接應用，將山區地表移動預計模型坐標方向和部分參數意義與前者調整對應即可。

皮爾森Ⅲ型公式法的傾向主斷面預計公式為［7］

皮爾森Ⅲ型公式法采用的走向主斷面概率積分法公式［19］如下：

水平地面任意點( x,y )的下沉值W( x,y )和沿φ方向的水平移動值U( x,y)φ的計算公式為［7］

已進行參數調整的山區地表移動預計模型［7］計算公式為

式中，W( x,y)、U( x,y)φ分別為相同地質、開采技術條件下，水平地面任意點( x,y )的下沉量和該點沿φ方向的水平移動量，按式（3）計算，m；Ws( x,y)、Us( x,y)φ分別為疊加山區地表移動預計模型后最終的下沉值和水平移動值，m；Dx,y為點( x,y )的地表特性系數，根據地表類型（表層土與地面植被特征、地貌）取值；α地x,y為點( x,y )的傾角，（°）；φ為點( x,y )的傾斜方向角，（°）；φ為計算方向角，（°）；φ、φ由x軸正向按逆時針方向計算；P[ x]、P[ y ]分別為傾向和走向主斷面的滑移影響函數，按下式計算：

式中，S 為工作面傾向方向的計算長度，即傾向主斷面下沉盆地全長，m；L 為工作面走向方向的計算長度，即走向主斷面計算長度，m；Wm為理論最大下沉值，參考皮爾森Ⅲ型公式法計算，m；A、p、t 為滑移影響參數，參考值分別為2π、2 和π；r 為主要影響半徑，計算公式為

式中，H( x,y )為地表點( x,y )的高程，m；H均為工作面底板平均高程，m；tanβ為主要影響角正切。

1.2 系統開發

結合皮爾森Ⅲ型公式和山區地表移動預計模型的特點，本研究在Matlab 環境下設計了山區急傾斜煤層開采沉陷預計系統，系統主界面如圖1 所示。

1.2.1 統一計算坐標系

如圖2所示，皮爾森Ⅲ型公式法的原始坐標系分別為傾向O1X 與走向O2Y 兩個單軸坐標系，為符合平面計算習慣，系統統一計算坐標系為由原始主斷面坐標系中O1X 與O2Y 分別向下、向左平移至O1、O2重合形成新的坐標系X1OY1，方便理解計算，本質上沒有改變任何變量，其中S3、S4為工作面走向左、右拐點偏移距。

1.2.2 坐標轉換工具箱

如圖3所示，工作面的計算坐標系根據煤層產狀而定（X2OY2），既有的經緯度及高程數據或投影直角平面坐標系以正東、正北為軸（X1OY1），因此需要將后者進行坐標轉換以匹配前者，方便計算。若后者為經緯度數據，工具箱將會自動按照經緯度與實際距離長度（m）的換算關系進行近似換算。

坐標軸旋轉公式為

式中，x1，y1為原始坐標，x2，y2為坐標軸繞原點逆向旋轉θ角后的坐標。

1.2.3 山區地表曲面函數擬合

山區地表移動預計時需要對坡面任意點的傾斜方向角、傾角、地表特征等進行取值，但由于地貌特性不一，人工取值的數據量巨大，難以實現。本研究對采集的地表三維數據利用Matlab 軟件的fit 函數進行多項式擬合，利用高程數據擬合出二元曲面函數方程，模擬山區的地形起伏形態，便于后續計算；利用三維向量夾角計算方法計算坡面任意點傾斜方向角和傾角，通過判斷任意點的凹凸性從而給定地表特性系數值。Matlab軟件中fit函數提供有三次、四次和五次多項式擬合方法，在本研究系統中分別對應的地形復雜程度為簡單、一般、復雜，主要參考地形相對高差、起伏復雜程度進行選擇。

部分函數代碼如下（“%”之后為代碼釋義）：

C=importdata（［path，file］，'，'）；%將X，Y，Z格式的三維高程數據文本導入

X=C（:，1）；

Y=C（:，2）；

Z=C（:，3）；

f=fit（［X，Y］，Z，'poly55'）；%利用fit函數將三維數據進行五次多項式擬合

syms x y；%定義x和y符號變量

z=f.p00+f.p10.*x+f.p01.*y+f.p20.*x.+f.p11.*x.*y+f.p02.*y.+f.p30.*x.+f.p21.*x..*y+f.p12.*x.*y.+f.p03.*y.+f.p40.*x.+f.p31.*x..*y+f.p22.*x..*y.+f.p13.*x.*y.+f.p04.*y.+f.p50.*x.+f.p41.*x..*y+f.p32.*x..*y.+f.p23.*x..*y.+f.p14.*x.*y.+f.p05.*y.+f.p23.*x..*y.+f.p14.*x.*y.+f.p05.*y.；%根據函數對象f 的屬性（f.p00 等）獲取系數值并生成曲面符號函數。

H=matlabFunction（z）；%將符號函數轉換為匿名函數，即地表擬合函數，為山區地表任意點的傾斜方向角、傾角計算提供方便.

1.2.4 地面移動變形函數計算

地表移動變形函數計算是整個系統實現中最重要的部分，該模塊將上述提及的各計算公式進行耦合編程，在上一模塊地表函數擬合完畢后，利用曲面函數和手動輸入的計算參數生成各個移動變形預計函數，為下一模塊移動變形圖繪制提供計算基礎，圖4所示為主要計算流程。

1.2.5 皮爾森Ⅲ型公式法

皮爾森Ⅲ型公式法主要由三部分組成：傾向主斷面計算、走向主斷面計算、地面任意點計算。傾向主斷面計算依據式（1）以匿名函數進行計算，走向主斷面計算依據式（2）概率積分法以匿名函數進行計算，地面任意點計算依據式（3）將前兩者組合而成，現僅給出下沉值和水平移動值函數計算代碼。

傾向主斷面下沉值與水平移動值函數部分代碼如下：

走向主斷面下沉值與水平移動值函數代碼如下：

地面任意點下沉值與水平移動值函數代碼如下：

syms x y；

Cx=Wx（x）./Wmax；

Cy=Wy（y）./Wmax；

w=Cx.*Cy.*Wmax；

u=Cy.*ux（x）.*cosd（Phi）+Cx.*uy（y）.*sind（Phi）；

Wxy=matlabFunction（w，'Vars'，［x y］）；

uxy=matlabFunction（u，'Vars'，［x y］）；%由Wxy 與uxy 即得到平地急傾斜煤層開采地表移動值任意點的匿名計算函數。

1.2.6 山區地表移動預計

由上述分析得到了地表曲面擬合函數，利用曲面法向量與X 軸和Z 軸的方向向量夾角計算方法，可以得到公式所需的坡面任意點的傾斜方向角和傾角，以此來代替對地表坡面的數據采集。

曲面法向量采用偏導數計算得到，程序代碼為

syms x y

vx=matlabFunction（-diff（H，x），'Vars'，［x y］）

vy=matlabFunction（-diff（H，y），'Vars'，［x y］）

vz=1.

坡面任意點的傾角計算代碼為

坡面任意點的傾斜方向角計算代碼為

dirAngle=@（x，y）Angle（x，y，vx，vy）；

function dirAngle=Angle（x，y，vx，vy）

if vy（x，y）＞0 %用于判斷傾斜方向角與x軸正向逆時針旋轉夾角是否大于180度

地表特性系數在凹凸處不同，根據曲面函數可計算任意點的凹凸性從而對該點賦值，計算函數為SurfaceValue( x,y)，可在任意點獲取對應值。

最后，山區地表移動預計模型[ 式 (4)]與皮爾森Ⅲ型公式法[ 式 (3)]疊加得到最終計算函數，下沉值與水平移動值計算的函數代碼如下：

曲率值、水平變形值、傾斜值等函數計算可采用對函數Ws 或us 求導的辦法得出。至此，根據所需要輸入的參數類型在系統用戶界面輸入相應的數值，點擊“計算”便可得到山區急傾斜煤層開采地表移動變形函數。

1.2.7 地面移動變形圖形繪制

在完成了地面移動變形函數計算后，本模塊“計算”按鈕自動開啟，可對地表下沉值、地表水平移動值、地表傾斜值、地表曲率值、地表水平變形值等進行二維及三維圖形繪制。為控制圖形繪制精度和時長，系統默認繪圖采樣間隔為1 m，隨著繪圖間隔增大，繪圖時長減小。默認繪圖范圍由上一模塊函數計算時自動生成，可在系統默認的繪圖范圍內進行修改。

根據繪圖采樣間隔可得計算坐標列表并生成網格采樣點，可將網格點直接代入函數即可得到對應網格點上的地表移動變形值，程序代碼為

x=Xmin:interval:Xmax；%interval為采樣間隔

y=Ymin:interval:Ymax；

［X，Y］=meshgrid（x，y）；

Z=Ws（X，Y）；

Z=us（X，Y）.

根據得到的X,Y,Z 值，便可通過圖形繪制函數生成對應的二維或三維圖形，代碼為

surf（X，Y，Z，'EdgeColor'，'interp'）%三維圖形

contour（X，Y，Z，v，'ShowText'，'on'，'LabelSpacing'，6000）%平面等值線圖。

1.2.8 地表任意點移動變形值計算

本研究開采沉陷預計系統可對計算范圍內任意點的地表移動變形值進行計算和查詢，可以單點查詢或批量查詢，單點查詢直接輸入相應點位的( )x,y坐標，批量查詢時可導入點位( )x,y 坐標的txt文檔，選擇計算值類型，點擊“計算”。單點查詢時可以直接顯示，批量查詢時則輸出為新的txt文件。

2 實例分析

根據甘肅省某礦區山區地下急傾斜煤層開采參數，采用本研究預計系統（登記號：2019SR0590923）對其進行地表移動變形預計和分析。

2.1 預計參數

該采區共有上下兩層煤，矩形工作面采用綜采放頂煤開采方法，煤層厚度分別為7.2～8.2 m 和29.87～44.82 m，平均厚度分別為7.7 m 和35.6 m，工作面走向長800～1 200 m，煤層平均傾角為52°～53°。根據實際礦區開采和觀測資料分析，結合“三下”開采規范［7］中提供的地表移動經驗參數值，確定的本研究開采沉陷預計參數如表1所示。

注：Kα、a1、a2、a3、b0、b1、b2、b3等參數參考文獻［7］取值。

2.2 開采沉陷預計

根據計算原點和采區范圍將研究區的地表三維高程數據轉換為符合系統統一坐標系的數據格式，而后通過Matlab 的fit 函數對該地區復雜地形進行五次多項式曲面擬合，擬合效果見圖5。

根據地表曲面擬合函數和預計參數進行地表移動變形函數計算并繪圖，由本研究預計系統繪制的地表下沉值的三維、二維平面圖如圖6、圖7 所示，繪制的地表傾向、走向移動值的三維、二維平面圖見圖8～圖11。

通過觀察地表下沉值圖形可看出等值線僅有輕微折曲現象，等值線表現較規則，整體來看地表下沉值受山區條件的影響較小；對地表傾向、走向移動值圖形觀察發現，等值線出現明顯的折曲變化，尤其對于走向移動方向的影響較強，而對于傾向方向移動影響較弱，反映出水平移動值受山區條件的影響較大。

2.3 數據驗證分析

礦區在地表布設有若干觀測點，并已獲取到一定時間段內的經緯度和高程觀測數據，通過坐標轉換工具箱將原始數據轉為系統內對應坐標，位于默認繪圖范圍內的觀測點有38個。將觀測點在統一坐標系內的坐標值輸入系統進行批量查詢，并輸出系統預計值。預計值與實際觀測值的對比結果如圖12、圖13、圖14所示。

本研究預計系統可針對地表移動變形值的變化提供較好的趨勢發展預計結果。圖12 中，點1002～1006 間預計下沉值較小，下沉趨勢呈現同步水平下沉，與實際觀測值呈現的凹形下沉趨勢差別較大，其他點位間的預計值與實際觀測數據之間呈現出良好的凹形發展趨勢。圖13 中，預計傾向移動值的變化趨勢與實際觀測值相似，但移動方向相背，趨勢表現不同步，表明在本例中對傾向移動值的預計沒有意義，但是否對所有的礦區都不適用，需要進一步研究。圖14 中，預計走向移動值與實際觀測值間的發展趨勢與圖12 的下沉值相似，在點1002～1006 間預計值與實測值的凹形發展趨勢相差較大，在點1007 之后的預計值與實測值整體表現出良好的凹形發展趨勢，僅在個別點位出現偏差。

通過進一步分析可知：下沉值與走向移動值的預計結果較好地符合了實測值的發展趨勢，符合發展趨勢的點位占全部觀測點的86.8%，而對傾向移動值的預計背離了實際，基本與實測值呈現相反的移動方向，無法提供有效的預計結果。在本例中所體現的數據誤差來自幾個方面，一是預計系統中地表曲面擬合產生的誤差，地形越復雜誤差越大，可能對水平移動值產生影響；二是采用的礦區參數不能完全反映實際開采情況，本例中工作面鄰近區內已經存在有其他開采工作面，疊加效應將影響到實測值的變化，且并未排除，因此出現了實測下沉值大于預計下沉值的現象；三是受地質條件、數據采集精度等因素影響。

3 結 論

（1）將皮爾森Ⅲ型公式法與山區地表移動預計模型進行疊加計算，可實現對山區急傾斜煤層開采地表移動變形的近似預計。實現方法主要是通過Matlab語言對計算公式封裝開發出可視化預計系統，系統將山區地表擬合生成曲面函數，為計算地表坡面任意點的傾斜方向角、傾角等參數提供了方便。另外，預計系統可以對最終生成的移動變形函數進行采樣繪圖，實現在預計范圍內的二維、三維可視化分析，從而可以更加直觀地展示預計結果和地表變形規律。

（2）通過結合實際礦區的開采條件和觀測數據，確定了采區相關參數，利用本研究預計系統進行地表移動變形預計，得到了觀測點的預計值，并將預計值與實測值進行了對比分析。結果表明：實測下沉值、走向水平移動值與預計值的整體發展趨勢基本符合，僅個別點位存在誤差，傾向水平移動值的預計效果較差，在本例中沒有參考意義，可能受到系統地表曲面擬合精度、采區地質條件、數據采集誤差等因素的影響，仍需要進一步驗證分析并繼續改進。通過分析可知：本研究開發的預計系統可以對山區急傾斜煤層開采引起的地表移動變形發展趨勢進行近似預計，但由于受到誤差影響，會出現偏離預計的現象，因此應盡可能保證各項開采沉陷預計參數精確取值。