趙 博 梁乃森 吳 初 武 雄

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

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

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

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

1.1 計算原理

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

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

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

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

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

式中，W( x,y)、U( x,y)φ分別為相同地質(zhì)、開采技術(shù)條件下，水平地面任意點( x,y )的下沉量和該點沿φ方向的水平移動量，按式（3）計算，m；Ws( x,y)、Us( x,y)φ分別為疊加山區(qū)地表移動預計模型后最終的下沉值和水平移動值，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ū)地表移動預計模型的特點，本研究在Matlab 環(huán)境下設計了山區(qū)急傾斜煤層開采沉陷預計系統(tǒng)，系統(tǒng)主界面如圖1 所示。

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

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

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

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

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

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

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

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

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

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

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ù)生成各個移動變形預計函數(shù)，為下一模塊移動變形圖繪制提供計算基礎，圖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ū)地表移動預計

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

曲面法向量采用偏導數(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)，可在任意點獲取對應值。

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

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

1.2.7 地面移動變形圖形繪制

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

根據(jù)繪圖采樣間隔可得計算坐標列表并生成網(wǎng)格采樣點，可將網(wǎng)格點直接代入函數(shù)即可得到對應網(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ù)生成對應的二維或三維圖形，代碼為

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

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

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

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

2 實例分析

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

2.1 預計參數(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ù)值，確定的本研究開采沉陷預計參數(shù)如表1所示。

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

2.2 開采沉陷預計

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

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

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

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

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

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

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

3 結(jié) 論

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

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