概率積分法開采沉陷預測的數值計算與分析

魏宗海，熊 偉

(河北省地礦局第三地質大隊，河北 張家口 075000)

地下煤炭資源的開采會改變地質體平衡狀態,進而引發地表沉陷問題[1]。由于煤層開采引起的地表沉陷預測在“三下”采煤、土地復墾治理等生產實踐方面具有非常重要的作用，例如：在建筑物下采煤時，根據預測所得的結果可以判別建筑物是否受開采影響和受開采影響的程度，并依此作為對建筑物采取維修、加固、重建等維護措施的依據[2-4]；在鐵路下開采時根據預計結果判斷鐵路下開采的可能性[5]；在水體下開采時，預計結果可以用來判斷礦井受水患威脅的程度及開采對堤壩等設施的破壞和影響程度；在土地復墾工作中，根據預計結果可以確定地表塌陷的深度、范圍和坡度等，并依此作為編制土地復墾報告的重要技術依據[6-7]。

目前，在我國普遍采用的開采沉陷預測方法是概率積分法，該方法已應用于大量的工程實踐中，并積累了豐富的實踐經驗、形成了相對成熟的預測參數體系，而且有依據規范的指導，是當前常用的地表開采沉陷預測方法。由于該方法需要通過計算復雜的二重積分來獲得地表開采沉陷值，人工計算過程繁瑣且容易出錯，對于多工作面、大范圍地表沉陷預測更是難以實現。隨著計算機技術的迅速發展，國內的不少煤炭研究院所和高校都開發了基于概率積分法的開采沉陷預測軟件，其中具有代表性的有：中國礦業大學開發的MSPS軟件、山東省魯南地勘院開發的SODP軟件、河南理工大學開發的MSDFVS軟件以及山東科技大學、安徽理工大學開發的基于GIS的地表開采沉陷預測程序等[8-9]。這些軟件的出現大大降低了沉陷預測工作的難度；另一方面，開采沉陷預測軟件的編制必然要通過數值積分的方法來實現概率積分法中二重積分的計算，而數值積分計算必然會帶來一定的計算誤差，雖然地表沉陷預測的精度主要取決于預測參數的準確度，然而對于需要多工作面影響疊加來預測地表沉陷的情況而言，由于數值積分計算導致的誤差通過累積會大大增加。因此，選擇合適的數值積分方法來編制地表開采沉陷預測軟件對提高地表沉陷預測的精度和預測效率都具有非常重要的實際意義。

目前對于概率積分法開采沉陷預測數值計算方法和計算精度方面的研究較少，針對這種情況，本文先根據概率積分法地表沉陷預測模型，提出了使用變步長辛卜生二重數值積分的方法對地表沉陷進行預測計算，最后結合一個工程實例對該方法的計算精度進行對比分析。

1 概率積分法開采沉陷預測方法

根據地表移動變形預計方法的建立途徑，可以將其分為3類：基于實測資料的經驗法、理論模擬法和影響函數法。概率積分法是影響函數法的一種，該方法以正態分布函數為影響函數，用積分式表示地表下沉盆地形態特征。在波蘭學者李特威尼申創立的巖層移動隨機介質理論的基礎上，我國學者劉寶琛、廖國華提出了實用的概率積分法開采沉陷模型。

如圖1所示，ABCD為某一緩傾斜煤層走向長壁工作面，工作面長為L1，推進長度為L2，工作面走向、下山和上山拐點偏移距分別為S0、S1、S2。在回采工作面下山建立如圖1所示的計算坐標系，則對于地表任意一點P(x,y)，概率積分法開采沉陷預測值算式為：

(1)

Wc m=mq·cosα,

(2)

(3)

式中：W(x，y)為地表任意點P的下沉值，mm；m為工作面采厚，mm；q為下沉系數；為煤層傾角，(°)；r為主要影響半徑，m；H0為工作面平均采深，m；β為主要影響角，(°)；u，v為開采單元在x，y方向上的坐標值，m。

圖1 概率積分法開采沉陷計算示意圖

2 沉陷預測的數值積分計算方法

2.1 數值積分原理

數值積分是用數值逼近的方法對給定的積分函數進行近似計算的一種方法。如果某一函數f(x)在區間[a,b]上連續且n+1階可導，則f(x)可近似為該區間上n+1個結點的n次插值多項式Pn(x)，即：

(4)

那么，f(x)在區間[a,b]上的積分可近似為函數Pn(x)在區間[a,b]上的積分，即：

(5)

式中，ck為求積系數。式(5)表明，任意一個在區間[a,b]上連續且n+1階可導的函數f(x)都可近似為該區間上某些結點處的函數值的線性組合。顯然，取不同的n+1個結點，可以得到不同的插值求積公式，當n=2時，式(5)對應的求積公式即為辛卜生求積公式，即：

(6)

需要注意的是，對于取n+1個結點構造的求積公式有n次代數精度，即：對于所有次數不超過n的多項式來說，求積結果均準確可靠；而對于n+1階多項式來說則不一定準確。

在計算二重積分時，可以將二重積分轉化為兩個一重積分的乘積，然后分別使用式(6)對一重積分進行計算，最后相乘即可得到二重積分結果。

2.2 開采沉陷預測的變步長辛卜生二重積分數值計算方法

根據上節介紹的數值積分原理，將含有二重積分的概率積分法開采沉陷預測公式表示為兩個一重積分G(u)、T(v)的乘積，即：

(7)

(8)

由回采工作面開采范圍可知，一重積分G(u)、T(v)的積分區間分別為[S1,L1-S2]、[S0,L2-S0]，由于直接使用式(6)在整個區間內對函數進行積分的誤差比較大，因此，需要將整個積分區間劃分為若干小區間，再在各個區間上分別使用式(6)進行求積，最后將各個區間的求積結果進行疊加，即可得到函數積分在整個積分區間的結果。在數值積分中，積分精度是一個非常重要的問題。一般來說，整個積分區間中劃分的各小區間長度(積分步長)越小，積分精度越高；但是，如果步長取的太小，計算量則會大大增加，而且由于累加積分項數的增加，累積誤差也會增大，進而會導致沉陷預測誤差過大。因此，在實際計算中，通常根據計算精度要求逐步將積分區間進行二等分(即，變步長)。

以計算一重積分G(u)為例，在使用變步長辛卜生求積法時，首先將積分區間[S1,L1-S2]分為n(n=2)個小區間，則求積結果為：

(9)

式中，h為積分步長。然后，再分別將每個小區間劃分為n個子區間，則有：

(10)

根據式(9)和式(10)構造辛卜生求積公式的積分結果為：

(11)

上述過程即為變步長辛卜生求積分方法的一個子計算過程。由于一個子計算過程往往不能對積分結果進行充分近似，通常情況下，需要不斷重復上述過程，直到前后兩次的積分值之差小于某一限制精度為止，即：

|S2n-Sn|<ε.

(12)

另外需要注意的是，為了防止積分步長無休止地劃分下去，從而導致累積誤差過大的現象發生，一般應對積分的最小步長加以限制；也就是說，當積分步長小于限制步長h0時，無論積分結果是否小于限制精度，都應停止運算，返回積分結果。

使用辛卜生變步長求積分方法計算開采沉陷預測值的算法可歸納為：

1)確定預測及計算參數。根據工作面地質采礦條件，確定積分區間范圍、概率積分法預測參數；根據預測需要，確定預測點的平面坐標P(x,y)；根據計算精度要求，確定積分限制精度和限制步長。

2)計算預測點P(x,y)在x方向上的開采沉陷預測函數g(u)的積分值。

①根據x方向的積分區間劃分數n計算步長h和積分值Zn；

②計算步長折半后的積分值Z2n，并構造辛卜生求積公式的計算結果Sn；

③重復步驟①、②，并判斷前后兩次積分結果(Sn，S2n)之差是否符合限制精度要求，判斷當前積分步長d是否達到限制步長要求；

④若前后兩次積分結果之差符合限制精度要求或積分步長達到限制步長要求，則返回積分結果。否則重復①～③的計算步驟。

3)計算預測點P(x,y)在y方向上的開采沉陷預測函數t(v)的積分值，計算過程與g(u)函數積分計算過程相同。

4)將x方向、y方向的積分結果G(u)、T(v)代入式(7)，獲得P(x，y)點的開采沉陷預測值WP(x,y)。

該算法的流程如圖2所示。

2.3 非標準坐標系的概率積分法開采沉陷預測

在實際開采沉陷預測中，往往需要使用礦山采掘工程平面圖(如Autocad格式)作為底圖進行地表沉陷的預計。一般來說，回采工作面在電子版采工圖坐標系中是非標準存在的(如圖3(Ⅰ)所示)，即：坐標原點不處于回采工作面傾斜方向與走向方向的交點；傾斜方向非x軸方向；走向方向非y軸方向。此時，為了使用前文所述的數值方法進行開采沉陷預測，就必須首先進行坐標系的變換，通過坐標平移和旋轉的方法將工作面及預測點坐標轉換到標準的計算坐標系統中，如圖3所示。

圖2 辛卜生二重數值積分開采沉陷預測算法流程

圖3 坐標系變換示意圖

設工作面傾斜方向與走向方向交點在原始坐標系O(xy)中的坐標為A(xA,yA)，則回采工作面角點坐標及預測點坐標(x,y)在平移后的坐標系O′(x′y′)中的坐標(x′,y′)為：

(13)

(14)

通過式(13)及式(14)將非標準坐標系下的工作面角點坐標及預測點坐標系轉換為標準計算坐標系下的相應坐標后，根據工作面角點及預測點在標準坐標系中坐標、概率積分法預測參數及采礦地質參數，利用前文所述的變步長辛卜生數值積分方法進行開采沉陷預測。標準計算坐標下任意預測點P″(x″y″)的沉陷值WP″(x″,y″)即為原始坐標系下對應點P(x,y)的沉陷預計值。

3 實例與分析

某煤礦1300工作面為該礦一采區首采走向長壁工作面，工作面所采煤層為下二疊系山西組底部之3煤層；工作面于2006-11-21開始回采，2007-03-16停采；該工作面長246 m，實際推進長度606 m，煤層傾角4.8°，平均采厚5.5 m；工作面上山采深389 m，下山采深410 m，平均采深399.5 m。為研究該工作面所在采區地表移動變形規律，根據工作面井上下實際情況，在該工作面上方布設了2條地表移動變形觀測線：一條全傾向觀測線和一條半走向觀測線，如圖4所示。2006-09-05對兩條觀測線進行了首次全面觀測，2008-11-27進行了最后一次全面觀測；期間對2條觀測線共進行了19次高程日常觀測和5次平面日常觀測。觀測站全面觀測、日常觀測次數和精度要求均符合《煤礦測量規程》中的有關規定要求。

圖4 工作面位置及測線布設圖

根據該礦相鄰礦井已有的地表移動監測資料，結合回采工作面地質采礦條件，確定該工作面概率積分法預測參數如表1所示。將變步長辛卜生數值積分的限制精度、限制步長分別設置為0.01和0.1，根據該工作面地質采礦參數和概率積分法預測參數對觀測線各點進行沉陷預測。另外，為了驗證文中提出的開采沉陷預測計算方法的準確性，同時使用MSPS軟件對測線各點進行預測。文章提出的方法計算結果、MSPS軟件計算結果和測線最后一次全面觀測的實測數據對比如圖5所示。需要注意的是，考慮到走向觀測線各點和傾向觀測線靠近1301工作面的Q9-Q13點受相鄰工作面重復采動的影響，文中只對傾向測線上Q1-Q8點的實測下沉值與預測值進行比較。

表1 概率積分法開采沉陷預測參數表

圖5 開采沉陷預測結果比對圖

從圖5可以看出，文中提出的方法和MSPS軟件預測結果是非常接近的，走向測線、傾斜測線的提出方法與MSPS軟件計算結果的標準差分別為12.3 mm、21.1 mm。造成這種差異的原因可能是MSPS軟件中概率積分函數的數值計算實現方法與文中提出方法不同；也可能是由于MSPS軟件也使用了文中提出的數值積分方法計算概率積分函數，但是軟件中設置的限制步長和限制精度與文中設置的精度不同造成的。從傾向觀測線上8個測點的實測數據來看，兩種方法與地表實測下沉值之間存在較大的差異，這種情況主要是由于所取的概率積分法預測參數與該區域最優預測參數之間的差異造成的。文中提出方法的預測結果、MSPS軟件預測結果相對于實測值的平均誤差、標準差和均方根誤差如表2所示。可以看出，提出方法的開采沉陷預測結果略優于MSPS軟件的計算結果。另外，從上述結果來看，雖然由于預測參數不足夠精確造成的誤差占沉陷預測誤差的大部分(最大誤差28.5 cm)，但是由于積分方法造成的計算誤差仍是一個不可忽視的因素(約1～3 cm)；特別是在大范圍、多采區、多工作面預測時，往往是通過單獨計算每個工作面的地表沉陷預測值，然后通過疊加的方法計算最終地表沉陷值，在這種情況下，各工作面的計算誤差累積后可能會遠超過由于預測參數不精確導致的誤差。

表2 辛卜生方法、MSPS軟件開采沉陷預測誤差表 mm

4 結 論

文章提出了使用變步長辛卜生二重數值積分對概率積分函數進行計算，從而實現地表開采沉陷預測的方法；針對非標準坐標系下回采工作面的沉陷預測問題，文章提出了通過坐標系平移、旋轉將非標準坐標系轉換為標準計算坐標系下工作面開采沉陷預測的問題。通過實例計算和分析表明：

1)變步長辛卜生數值積分方法可以準確、可靠地對地表開采沉陷預測值進行計算；

2)基于變步長辛卜生數值計算方法的開采沉陷預測結果略優于MSPS軟件計算結果；

3)雖然由于預測參數不足夠準確導致的誤差為沉陷預測誤差的主要誤差，但是由于積分運算導致的誤差仍不可忽視，在開采沉陷預測中應盡量選擇合適的數值積分方法，減小由于計算導致的沉陷預測誤差。

雖然文中只探討了基于概率積分法的地表沉陷數值積分方法，由于其它地表變形值(如傾斜、曲率等)的預測公式與概率積分法地表下沉值預測公式類似，文中提出的數值積分方法同樣適用于地表傾斜、曲率、水平移動和水平變形預測值的數值計算。