基于下沉等影響原理的傾斜煤層開采沉陷預計方法及應用

田國燦，徐乃忠，高 超

(1.煤炭科學研究總院 開采研究分院，北京 100013；2.天地科技股份有限公司 開采設計事業(yè)部，北京 100013)

基于下沉等影響原理的傾斜煤層開采沉陷預計方法及應用

田國燦1,2，徐乃忠1,2，高 超1,2

(1.煤炭科學研究總院 開采研究分院，北京 100013；2.天地科技股份有限公司 開采設計事業(yè)部，北京 100013)

通過對以往開采沉陷的預計算法進行分析，采用下沉等影響原理建立沉陷預計的模型，設計適用于任意形狀工作面的開采沉陷預計程序。該程序采用C#語言編寫，進一步優(yōu)化了工作面各角點坐標和最終計算結果的處理，并通過對計算區(qū)域采用三角形劃分法的對比得知本模型可以進一步消除傾斜煤層在開采下沉計算中由采深的變化所帶來的誤差，為傾斜煤層的預計提供了更為準確的計算方法。

下沉等影響原理；傾斜煤層；開采沉陷；C#語言；概率積分法

地表沉陷預計是礦山開采沉陷學的重要內容之一，對開采沉陷的理論研究和生產實踐均具有重要意義[1]。在理論上，利用預計的結果可以定量地分析研究受開采影響的地表在時間上和空間上的分布規(guī)律。為了提高開采沉陷預計的準確性，有必要對預計方法所采用的理論模型、預計參數與地質采礦條件之間的定量、定性關系進行深入地研究，從而能更進一步地加深對開采沉陷基本規(guī)律的認識。在生產實踐上，利用開采沉陷預計結果可以對“三下”(建筑物下、鐵路下和水體下)的開采活動進行技術指導[2]。

因此，在借鑒已有理論體系的基礎上，應用下沉等影響原理和概率積分法的基本理論，開發(fā)一個能對任意形狀的傾斜工作面進行處理并具有足夠高精度的開采沉陷預計程序十分必要[3]。

1 程序設計技術路線

在開發(fā)環(huán)境方面，C#作為目前主流的編程語言之一，既繼承了C++語言的強大功能，又有VB的快速編寫方式，對于構建系統界面具有無可比擬的優(yōu)勢。既可以保證快速、穩(wěn)定地構建應用程序，又能在進行開發(fā)的時候調用現有的、無論是什么語言編寫的COM對象，具有極強的交互性。因此決定采用C#語言來完成系統框架的構建[4]。系統整體構架見圖1。

如圖1所示，程序的運行首先需要獲取工作面的相關信息，工作面各角點的坐標不再需要手工錄入，可通過對CAD組件的調用實現直接對任意形狀的多邊形進行角點坐標捕捉。

地質采礦條件參數的輸入包含煤層采高、工作面最大采深、工作面煤層傾斜方向、煤層傾角、下沉系數、主要影響角正切、水平移動系數、開采影響傳播角以及上山、下山、開切眼邊界、停采線邊界的拐點偏移距。

初始參數設置完之后對所有角點進行沿重心的坐標旋轉，使走向方向沿X軸正方向，傾向沿Y軸正方向，以便于隨后求得沿走向或傾向的地表移動變形值。然后通過設計好的預計模型對所要預計的影響范圍、任意剖面或任意點的移動變形進行計算，并對計算結果的坐標點再次進行坐標旋轉，以達到與原坐標系相匹配的目的。最后對數據進行后處理，生成剖面圖或者等值線圖。

2 等影響原理的算法實現

目前我國應用廣泛且較為成熟的預計方法主要是概率積分法，但由于工作面多邊形的不規(guī)則性，使得二重積分的上下限難以界定，所以大部分預計模型都是將工作面的不規(guī)則多邊形切割劃分為規(guī)則的矩形來進行疊加運算。即使有國內學者提出從工作面坐標原點出發(fā)(如圖2中的1角點)分別連接工作面的其他各點，將整個區(qū)域劃分為若干個三角形并對每個三角形進行二重積分的三角形劃分法，以此來確定積分上下限并消除由矩形疊加算法引起的計算誤差[5]，但該解決方法依然沒有對地表沉陷主要影響半徑r值的變化因素進行考慮，只是采用平均采深值將r值設定為一個常量。但在實際開采過程中，傾斜煤層條件下的采深是隨著工作面位置的移動在不斷變化的，工作面的傾角越大，采用平均采深算法的誤差也就越大。因此，有必要建立一個可以對任意形狀工作面以及任意煤層傾角進行沉陷預計的模型。

圖2 三角形劃分法示意

2.1 等影響工作面及相關變量的轉換

由沉陷預計概率積分法計算公式，地下坐標點U(x,y)處的單元開采引起地表坐標P(x0,y0)處的下沉值為：

(1)

r=H/tanβ

(2)

式中，r為主要影響半徑；H為采空區(qū)坐標點(x,y)處的采深；tanβ為地表沉陷預計參數中的主要影響角正切。

在開采區(qū)域較大，并且煤層傾角不為0的情況下，采深H值隨著煤層的傾向和走向在不斷地變化。若圖3中A點為最大采深位置，其坐標為(xa,ya)，工作面ABCD上任意一點U(x,y)的采深值則需要由煤層傾角、煤層傾向以及最大采深來計算求得。

圖3 實際工作面與投影工作面

若煤層傾角為α，煤層傾向為θ，最大采深為H0，則開采區(qū)域任意點(x,y)的采深H為：

(3)

若將所得的各角點坐標沿逆時針方向旋轉(90°-θ)，則新生成的工作面傾向沿y軸正方向，此時θ值為90°，進而消去式(3)中的參數cotθ，同時規(guī)范了坐標系，并為之后的預計計算做準備。旋轉之后的任意點(x,y)的采深H為：

H=H0-(ya-y)tanα

(4)

將地表點P(x0,y0)與工作面上任意點U(x,y)相連接并延長至最大點采深所在的水平面M得到點Q(x1,y1)，并由圖4可知等比關系

圖4 工作面點與投影點的坐標關系

(5)

同理

(6)

如果用r1來表示投影至水平面的等價工作面開采主要影響半徑，則

r1=H0/tanβ

(7)

將式(4)、(5)、(6)、(7)帶入(1)可得

(8)

可以看出投影面上的Q點(x1,y1)與工作面上的U點(x,y)對地表點(x0,y0)的影響是等價的[6]。然后對公式(1)進行面積分得到整個工作面的開采對任意地表點(x0,y0)產生的下沉值：

(9)

式中，wmax為該工作面開采后地表最大下沉值。

對式(9)進行換元并積分可得

w(x0,y0)=

(10)

同時由式(5)、(6)可得

(11)

(12)

最終將式(11)和式(12)帶入式(10)中可得

(13)

該式即為轉換之后的等價投影面的下沉預計面積分公式，D1為地表任意點(x0,y0)與工作面各角點鏈接之后投影至最大采深平面的多邊形(如圖3中多邊形ABC1D1)，在新生成的多邊形中主要影響半徑值為常量，同時每一個要計算的地表點都有與之對應的各不相同的等價投影多邊形。

2.2 線積分轉化

投影面轉換完成D1之后，D1依然是不規(guī)則圖形，對面積分的上下限進行判定并求取二重積分的算法比較復雜，因此將式(13)轉化為線積分對其進行求解。

由數學分析可知，當區(qū)域D是一個單連通域，且函數P(X,Y)，Q(X,Y)在D上具有一階連續(xù)偏導數的條件下，則下面等式成立：

這就是線積分和面積分轉換的格林公式，其中S為區(qū)域D的邊界曲線[7-9]。

由式(13)可知r1為常量的情況下對x進行積分相對比較簡單，又由高斯誤差函數

可知

其中區(qū)域D1是投影在平面M上的由直線(S1，S2…Sn)所圍成的閉合區(qū)域，高斯誤差函數可以單獨進行計算。至此將r值不斷變化的二重積分，轉換為r值為常量的定積分，在計算過程中，只需將地表點所對應的投影面的各角點坐標計算出來，然后對每一條邊進行上式的線積分計算，即可得到與原開采區(qū)域D相同的沉陷預計結果。該方法使工作面邊界得以確定，并考慮了開采影響范圍r值不斷變化對地表下沉產生的影響，同時將二重積分等價轉化為定積分減少了由計算機計算二次積分所帶來的誤差[10]。對該式進行相應的求導變換，即可得到傾斜、曲率、水平變形等地表移動變形函數，由于這些函數均為定積分，不存在上下限難以界定的問題，只需沿每條邊逐一積分即可得出計算結果，水平移動的計算函數也可通過上述的換元方法最終劃為相應的線積分形式，在此不再贅述。

3 程序結構

本程序主要由如下幾個模塊所組成：

準備模塊 錄入地質條件參數以及工作面各角點坐標。

地表影響范圍計算模塊 自動計算開采沉陷影響范圍，以2倍的主要影響半徑為擴展邊界，并確定要計算的地表點坐標。

坐標轉換模塊 將坐標系統轉換為以煤層走向沿x軸正方向，傾向(上山方向)沿y軸正方向。

變形計算模塊 結合編寫完成的預計模型，將經過左邊轉換后的各個參數導入，得到任意地表點的移動變形值。

輸出模塊 生成影響范圍內各地表點的地表移動變形值表格并繪制等值線圖。

程序結構框圖見圖5。

圖5 程序結構

4 應用實例

設定一不規(guī)則形狀工作面，走向長度約1000m，傾向長度約450m，下山方向采深500m，煤層采出厚度3000mm，煤層傾角13°，傾向沿y軸正方向，下沉系數為0.75，主要影響角正切為2.0，最大下沉角90°，水平移動系數0.3，拐點偏移距均為0。將其輸入到編寫好的沉陷預計程序界面中，并在CAD中獲取工作面各角點坐標值。然后用建立好的模型對其進行下沉計算，可得下沉等值線圖形，其中多段線為設定的工作面邊界，最大下沉值為2146mm。采用將工作面劃分為多個三角形的三角形剖分法獲得的下沉等值線圖形，最大下沉值為2095mm。

將2個等值線圖進行疊加對比得到圖6?？梢钥闯霎斂紤]煤層傾角所帶來的采深變化因素時，相對于三角形劃分算法，開采損害影響的邊界范圍變得較小，但是最大下沉值較三角形劃分法的計算結果增大了51mm。

當下沉值越大的時候，下沉等影響算法的影響范圍就越為明顯地大于三角形剖分法的影響范圍。煤層傾角越大，這種變化就會越為明顯。這種整體影響范圍變小但是一定范圍內下沉值變大的變化特點符合傾斜煤層所引起的實際下沉情況，同時也滿足主要影響角不變，采深越小影響范圍越小的理論。

5 結 論

(1)通過下沉等影響原理，將主要影響半徑化為常量，推導出傾斜煤層開采沉陷預計的精確計算公式。

(2)將二重積分形式轉換為線積分形式，確定了積分的上下限，進一步簡化了積分的過程。

(3)通過編程語言建立了開采沉陷預計模型，并與傳統的三角形劃分模型進行對比，分析了2種模型得出的不同預計結果。得知本模型更為符合實際情況，同時也使得概率積分法能夠更為靈活地應用在傾斜煤層開采沉陷的相關預計與計算當中。

[1]何國清，楊 倫，凌賡娣，等.礦山開采沉陷學[M].徐州：中國礦業(yè)大學出版社，1991.

[2]黃樂亭.村莊壓煤開采的技術途徑[J].中國煤炭，1996(8)：33-35.

[3]王 玲，吳 侃，奚新麗.基于ArcGis Engine的開采沉陷預計系統[J].測繪科學，2008(4)：153-154.

[4]李培現，譚志祥，齊公玉.基于MATLAB的開采沉陷預計系統[J].中國礦業(yè)，2008(11)：72-76.

[5]康建榮.任意形多邊形工作面多段線開采沉陷預計系統 (MSPS)[J].礦山測量，2000(1)：24-26.

[6]國家煤炭局.建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采規(guī)程[S].北京：煤炭工業(yè)出版社，2000.

[7]張華興，焦傳武.地表移動精確預計法——線積分運算[J].礦山測量，1998(4)：21-29.

[8]周萬茂，張華興，何瑞華.任意形狀工作面拐點移動距求取方法[J].煤礦開采，2000(4)：13-16.

[9]張華興.地表移動與變形的等價轉換線積分計算[J].煤礦開采，1991(2)：14-19.

[10]田錦州，徐乃忠，李鳳明.誤差函數erf(x)近似計算及其在開采沉陷預計中的應用[J].煤礦開采，2009，14(2)：33-35.

[責任編輯：李 青]

Prediction Method and Application of Incline Coal Seam Mining Subsidence Based on Subsidence Principle

TIAN Guo-can1,2,XU Nai-zhong1,2，GAO Chao1,2

(1.Mining Institute，China Coal Research Institute，Beijing 100013，China;2.Coal Mining & Designing Department，Tiandi Science & Technology Co.，Ltd.，Beijing 100013，China)

The past mining subsidence prediction method was analyzed，subsidence prediction model was built based on subsidence principle，mining subsidence prediction program that could be applied any shape working face was designed.The program was designed by C# language，every corner point coordinate of working face and the final results were optimized，the error that induced by mining depth variation in incline coal seam mining subsidence prediction could be eliminated in the model by triangle division method compared in calculated area，and then more accurate prediction method could be provided for incline coal seam prediction.

subsidence principle; incline coal seam; mining subsidence; C# language; probability integrity method

2016-10-27

10.13532/j.cnki.cn11-3677/td.2017.02.017

國家自然科學基金(51604151)

田國燦(1991-)，男，河北衡水人，在讀碩士研究生，研究方向為開采沉陷預計。

田國燦，徐乃忠，高 超.基于下沉等影響原理的傾斜煤層開采沉陷預計方法及應用[J].煤礦開采，2017，22(2)：65-68.

TD327

A

1006-6225(2017)02-0065-04