基于Matlab的山區(qū)急傾斜煤層開采沉陷預(yù)計系統(tǒng)

趙 博 梁乃森 吳 初 武 雄

（中國地質(zhì)大學(xué)（北京）水資源與環(huán)境學(xué)院，北京100083）

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

1 系統(tǒng)設(shè)計及編程實現(xiàn)

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

1.1 計算原理

本研究系統(tǒng)開發(fā)是建立在對皮爾森Ⅲ型公式法和山區(qū)地表移動預(yù)計模型疊加的基礎(chǔ)上，因此實現(xiàn)兩者參數(shù)意義匹配較為重要。通過對公式分析可知，兩者參數(shù)匹配不存在重大矛盾，皮爾森Ⅲ型公式法可直接應(yīng)用，將山區(qū)地表移動預(yù)計模型坐標(biāo)方向和部分參數(shù)意義與前者調(diào)整對應(yīng)即可。

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

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

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

已進行參數(shù)調(diào)整的山區(qū)地表移動預(yù)計模型［7］計算公式為

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

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

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

1.2 系統(tǒng)開發(fā)

結(jié)合皮爾森Ⅲ型公式和山區(qū)地表移動預(yù)計模型的特點，本研究在Matlab 環(huán)境下設(shè)計了山區(qū)急傾斜煤層開采沉陷預(yù)計系統(tǒng)，系統(tǒng)主界面如圖1 所示。

1.2.1 統(tǒng)一計算坐標(biāo)系

如圖2所示，皮爾森Ⅲ型公式法的原始坐標(biāo)系分別為傾向O1X 與走向O2Y 兩個單軸坐標(biāo)系，為符合平面計算習(xí)慣，系統(tǒng)統(tǒng)一計算坐標(biāo)系為由原始主斷面坐標(biāo)系中O1X 與O2Y 分別向下、向左平移至O1、O2重合形成新的坐標(biāo)系X1OY1，方便理解計算，本質(zhì)上沒有改變?nèi)魏巫兞?，其中S3、S4為工作面走向左、右拐點偏移距。

1.2.2 坐標(biāo)轉(zhuǎn)換工具箱

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

坐標(biāo)軸旋轉(zhuǎn)公式為

式中，x1，y1為原始坐標(biāo)，x2，y2為坐標(biāo)軸繞原點逆向旋轉(zhuǎn)θ角后的坐標(biāo)。

1.2.3 山區(qū)地表曲面函數(shù)擬合

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

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

C=importdata（［path，file］，'，'）；%將X，Y，Z格式的三維高程數(shù)據(jù)文本導(dǎo)入

X=C（:，1）；

Y=C（:，2）；

Z=C（:，3）；

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

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.；%根據(jù)函數(shù)對象f 的屬性（f.p00 等）獲取系數(shù)值并生成曲面符號函數(shù)。

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

1.2.4 地面移動變形函數(shù)計算

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

1.2.5 皮爾森Ⅲ型公式法

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

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

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

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

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 即得到平地急傾斜煤層開采地表移動值任意點的匿名計算函數(shù)。

1.2.6 山區(qū)地表移動預(yù)計

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

曲面法向量采用偏導(dǎo)數(shù)計算得到，程序代碼為

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軸正向逆時針旋轉(zhuǎn)夾角是否大于180度

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

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

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

1.2.7 地面移動變形圖形繪制

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

根據(jù)繪圖采樣間隔可得計算坐標(biāo)列表并生成網(wǎng)格采樣點，可將網(wǎng)格點直接代入函數(shù)即可得到對應(yīng)網(wǎng)格點上的地表移動變形值，程序代碼為

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

y=Ymin:interval:Ymax；

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

Z=Ws（X，Y）；

Z=us（X，Y）.

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

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

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

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

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

2 實例分析

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

2.1 預(yù)計參數(shù)

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

注：Kα、a1、a2、a3、b0、b1、b2、b3等參數(shù)參考文獻(xiàn)［7］取值。

2.2 開采沉陷預(yù)計

根據(jù)計算原點和采區(qū)范圍將研究區(qū)的地表三維高程數(shù)據(jù)轉(zhuǎn)換為符合系統(tǒng)統(tǒng)一坐標(biāo)系的數(shù)據(jù)格式，而后通過Matlab 的fit 函數(shù)對該地區(qū)復(fù)雜地形進行五次多項式曲面擬合，擬合效果見圖5。

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

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

2.3 數(shù)據(jù)驗證分析

礦區(qū)在地表布設(shè)有若干觀測點，并已獲取到一定時間段內(nèi)的經(jīng)緯度和高程觀測數(shù)據(jù)，通過坐標(biāo)轉(zhuǎn)換工具箱將原始數(shù)據(jù)轉(zhuǎn)為系統(tǒng)內(nèi)對應(yīng)坐標(biāo)，位于默認(rèn)繪圖范圍內(nèi)的觀測點有38個。將觀測點在統(tǒng)一坐標(biāo)系內(nèi)的坐標(biāo)值輸入系統(tǒng)進行批量查詢，并輸出系統(tǒng)預(yù)計值。預(yù)計值與實際觀測值的對比結(jié)果如圖12、圖13、圖14所示。

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

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

3 結(jié) 論

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

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