張 兵，崔希民，袁德寶，賀軍亮，郭婭玲

(1.石家莊學(xué)院 資源與環(huán)境科學(xué)學(xué)院，河北 石家莊 050035；2.中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院，北京 100083)

0 引 言

開采沉陷預(yù)計(jì)方法中，無論是概率積分法，還是典型曲線法等，計(jì)算的都是地表穩(wěn)定后的沉降和變形值[1-3]，經(jīng)過多年發(fā)展，我國學(xué)者在靜態(tài)預(yù)計(jì)方面的研究成果較為豐富，且非常成熟。眾所周知，地表移動(dòng)和變形是一個(gè)復(fù)雜的過程，與多因素有關(guān)，如開采速度、覆巖性質(zhì)、煤層傾角等[4-8]。研究地表隨時(shí)間變化的動(dòng)態(tài)過程同樣具有重要的現(xiàn)實(shí)意義，可為采前建筑物保護(hù)、采后損害鑒定、采空區(qū)土地復(fù)墾利用等提供指導(dǎo)。國外動(dòng)態(tài)預(yù)計(jì)研究可追溯到20世紀(jì)20—30年代，我國學(xué)者直到20世紀(jì)80年代，才開始重視動(dòng)態(tài)預(yù)計(jì)相關(guān)研究，并在近年來取得了較大的進(jìn)展。崔希民等[6]對(duì)最常用的動(dòng)態(tài)預(yù)計(jì)“Knothe時(shí)間函數(shù)”進(jìn)行了研究，指出了其在理論上的不足，給出了理想時(shí)間函數(shù)的圖像形態(tài)；常占強(qiáng)等[9]為了改善Knothe時(shí)間函數(shù)在地表下沉初始階段預(yù)計(jì)精度較低的問題，將該函數(shù)進(jìn)行了分段表達(dá)，建立了“分段Knothe時(shí)間函數(shù)模型”。胡青峰等[10]對(duì)Knothe時(shí)間參數(shù)的影響因素進(jìn)行了分析，給出了時(shí)間參數(shù)的直接求取方法，提高了編程計(jì)算的效率。陳磊等[11]在動(dòng)態(tài)預(yù)計(jì)中使用冪指數(shù) Knothe 時(shí)間函數(shù)模型，并利用InSAR與水準(zhǔn)測量數(shù)據(jù)對(duì)該時(shí)間函數(shù)進(jìn)行求參。郭旭煒等[12]對(duì)分段Knothe時(shí)間函數(shù)的參數(shù)進(jìn)行了研究，指出其參數(shù)會(huì)隨開采過程而變化，并給出了動(dòng)態(tài)求參方法。另外，在動(dòng)態(tài)預(yù)計(jì)方法方面，楊澤發(fā)等[13]利用InSAR觀測值與Logistic模型之間的函數(shù)關(guān)系求取后者時(shí)間參數(shù)，進(jìn)而再采用Logistic模型進(jìn)行動(dòng)態(tài)預(yù)計(jì)。李懷展等[14]基于加權(quán)法和地質(zhì)力學(xué)參數(shù)的敏感性，提出了一種地表沉降動(dòng)態(tài)計(jì)算模型，并用實(shí)例說明了該方法具有較高的精度。

通過文獻(xiàn)分析可知，在動(dòng)態(tài)預(yù)計(jì)方面，學(xué)者們提出了各種方法，其中對(duì)Knothe時(shí)間函數(shù)的研究和改進(jìn)涉及較多，說明了該時(shí)間函數(shù)仍然具有很強(qiáng)的適用性；另外，現(xiàn)有方法中，很多時(shí)間函數(shù)參數(shù)的求取依賴大量實(shí)測數(shù)據(jù)，且模型較為復(fù)雜，在實(shí)際應(yīng)用中存在較大困難。為了操作簡便，保證較高的預(yù)計(jì)精度，選擇“分段Knothe時(shí)間函數(shù)”作為動(dòng)態(tài)預(yù)計(jì)時(shí)間函數(shù)—筆者曾對(duì)其進(jìn)行過優(yōu)化，并給出可靠時(shí)間參數(shù)計(jì)算方法[15-16]—研究走向主斷面地表沉陷動(dòng)態(tài)預(yù)計(jì)方法，建立預(yù)計(jì)公式并設(shè)計(jì)相應(yīng)的程序算法。限于篇幅，傾向主斷面和全盆地動(dòng)態(tài)預(yù)計(jì)方法將另文研究。

1 動(dòng)態(tài)預(yù)計(jì)基本原理

1.1 動(dòng)態(tài)單元?jiǎng)澐?/h3>

在動(dòng)態(tài)預(yù)計(jì)實(shí)踐中，國內(nèi)外的普遍做法是通過將地下開采單元所引起的沉陷靜態(tài)預(yù)計(jì)結(jié)果乘以對(duì)應(yīng)的時(shí)間函數(shù)值來實(shí)現(xiàn)。因此，在動(dòng)態(tài)預(yù)計(jì)中除了選擇合理的時(shí)間函數(shù)外，還需要根據(jù)開采情況對(duì)工作面進(jìn)行合理的劃分，通常是將工作面劃分為距離相等或不相等的多個(gè)單元，稱其為“動(dòng)態(tài)開采單元”，劃分方法主要有“周期來壓步距法”和“有效尺寸分割法”[17-18]。

在實(shí)踐中，為了方便操作，通常采用有效尺寸分割法，如圖1所示(X為工作面到開切眼的距離，W為地表下沉量)。

H—采深；wi—某點(diǎn)地表下沉量圖1 動(dòng)態(tài)單元開采對(duì)地表點(diǎn)下沉的動(dòng)態(tài)影響Fig.1 Dynamic influence when mining dynamic unit on surface subsidence

圖1中v1t1表示第1個(gè)動(dòng)態(tài)單元的長度，v1表示第1個(gè)動(dòng)態(tài)單元的開采速度，t1表示第1個(gè)動(dòng)態(tài)單元采完所經(jīng)歷的時(shí)間；同理，v2t2表示第2個(gè)動(dòng)態(tài)單元的長度，vntn表示第n個(gè)動(dòng)態(tài)單元的長度。w1表示第n個(gè)單元?jiǎng)偛赏陼r(shí)，第1個(gè)動(dòng)態(tài)單元開采對(duì)地表下沉的影響，w2表示第n個(gè)動(dòng)態(tài)單元?jiǎng)偛赏陼r(shí)，第2個(gè)動(dòng)態(tài)單元對(duì)地表下沉的影響；w3～wn的意義以此類比。對(duì)比可知，第1個(gè)動(dòng)態(tài)單元的影響最大，第n個(gè)單元的影響最小，原因是，當(dāng)?shù)趎個(gè)動(dòng)態(tài)單元?jiǎng)偛赏陼r(shí)，第1個(gè)動(dòng)態(tài)單元對(duì)地表影響的總時(shí)間是t1、t2、…、tn之和，而第2個(gè)動(dòng)態(tài)單元對(duì)地表影響的總時(shí)間則是t2、t3、…、tn之和，比第1個(gè)單元的影響時(shí)間短，因此其影響也小，同理，由于影響時(shí)間不同，其他單元的影響會(huì)依次減小。在t1+t2+…+tn時(shí)刻，將工作面1～n個(gè)單元開采對(duì)地表的影響進(jìn)行疊加求和，即可得到此時(shí)刻的地表下沉，圖1中Wt即t1+t2+…+tn時(shí)刻的動(dòng)態(tài)下沉曲線。因此，要想求某時(shí)刻T的地表下沉量，需首先計(jì)算各單元在T時(shí)刻的動(dòng)態(tài)影響w1、w2、…、wn，這就需要通過采用時(shí)間函數(shù)，將各動(dòng)態(tài)單元的靜態(tài)影響與時(shí)間系數(shù)聯(lián)系起來，用系數(shù)對(duì)其影響進(jìn)行調(diào)整。

1.2 動(dòng)態(tài)預(yù)計(jì)時(shí)間函數(shù)及其計(jì)算

本文選擇的時(shí)間函數(shù)是“優(yōu)化分段Knothe時(shí)間函數(shù)”，簡稱：“優(yōu)化分段時(shí)間函數(shù)”，筆者曾改進(jìn)了原分段函數(shù)的不足[15]，并通過實(shí)測數(shù)據(jù)證明了其在動(dòng)態(tài)預(yù)計(jì)中的可靠性。優(yōu)化分段時(shí)間函數(shù)理論模型如式(1)所示：

(1)

式中：τ為地表點(diǎn)出現(xiàn)最大下沉速度時(shí)刻；T為地表移動(dòng)變形總時(shí)間；c為與地質(zhì)、采礦條件有關(guān)的時(shí)間系數(shù)；t為給定的預(yù)計(jì)時(shí)刻。

在動(dòng)態(tài)預(yù)計(jì)時(shí)，首先要計(jì)算每個(gè)動(dòng)態(tài)單元所對(duì)應(yīng)的時(shí)間函數(shù)值，函數(shù)值可理解為開采影響調(diào)整系數(shù)。如果在t時(shí)刻進(jìn)行預(yù)計(jì)，則第1個(gè)動(dòng)態(tài)單元的調(diào)整系數(shù)可用Φt表示；第2個(gè)動(dòng)態(tài)單元的調(diào)整系數(shù)用Φ(t-t1)表示；第n個(gè)動(dòng)態(tài)單元的調(diào)整系數(shù)則用Φ(t-t1-…-tn-1)表示，根據(jù)時(shí)間函數(shù)的意義，第1個(gè)動(dòng)態(tài)單元時(shí)間函數(shù)值可采用式(1)直接求得，第2到第n個(gè)動(dòng)態(tài)單元的時(shí)間調(diào)整系數(shù)可用式(2)、式(3)等進(jìn)行計(jì)算。

第2個(gè)動(dòng)態(tài)開采單元時(shí)間函數(shù)值計(jì)算公式為

(2)

同理，可求第n個(gè)單元時(shí)間函數(shù)式，即

(3)

2 走向主斷面動(dòng)態(tài)預(yù)計(jì)理論模型

2.1 走向動(dòng)態(tài)單元有限開采原理

對(duì)于每個(gè)動(dòng)態(tài)單元，其開采影響符合有限開采原理，可采用概率積分模型，采用有限開采相應(yīng)計(jì)算公式，具體如圖2所示。

圖2 有限開采地表走向主斷面預(yù)計(jì)疊加原理Fig.2 Superposition principle in strike main section prediction when finite mining

圖2中，s3、s4為工作面左右拐點(diǎn)偏移距；AB為實(shí)采邊界，計(jì)算時(shí)坐標(biāo)原點(diǎn)位于O處，W(x)表示將x>-s3的所有煤層采出引起的地表下沉曲線；而W(x-l)對(duì)應(yīng)的是AB煤層未采，只采x>(l+S4)的煤層時(shí)的地表下沉，那么，如果僅開采AB煤層，其對(duì)地表下沉量W0(x)就可用W0(x)=W(x)-W(x-l)求解，其結(jié)果可用圖中W0(x)表示。如果不考慮拐點(diǎn)偏距的影響，將坐標(biāo)原點(diǎn)設(shè)在P處，則AB煤層開采影響可用W0(x)=W(x)-W(x-D3)來計(jì)算。

2.2 走向主斷面動(dòng)態(tài)預(yù)計(jì)模型構(gòu)建

在建立模型時(shí)，將坐標(biāo)原點(diǎn)設(shè)在圖2中P點(diǎn)處，縱軸表示地表下沉，x軸表示地表點(diǎn)位置。結(jié)合圖1，第1個(gè)動(dòng)態(tài)單元開采后，該單元的起點(diǎn)橫坐標(biāo)為0，終點(diǎn)橫坐標(biāo)為v1t1，設(shè)其下沉影響用W1(x)表示，考慮拐點(diǎn)偏距影響，參照有限開采原理[1]，其終態(tài)影響可用式(4)計(jì)算，即

W1(x)=W(x-s3)-W[x-s3-(v1t1-s3)]=
W(x-s3)-W(x-v1t1)

(4)

如果其他動(dòng)態(tài)單元均未開采，只開采第2個(gè)動(dòng)態(tài)單元，則其起點(diǎn)橫坐標(biāo)x為v1t1，終點(diǎn)橫坐標(biāo)為v2t2，依照式(4)，其對(duì)地表下沉的終態(tài)影響可用式(5)計(jì)算，即

W2(x)=W(x-v1t1)-W(x-v1t1-v2t2)

(5)

同理，可推導(dǎo)第n個(gè)動(dòng)態(tài)單元單獨(dú)開采后Wn(x)的計(jì)算式，即

Wn(x)=W(x-v1t1-…-vn-1tn-1)-W[x

-(v1t1+…+vn-1tn-1+vntn-s4)]

(6)

根據(jù)概率積分法原理，上述公式中W(x)按照式(7)進(jìn)行計(jì)算。即

(7)

式中，W0為最大下沉量。

式(4)—式(6)計(jì)算的是各動(dòng)態(tài)單元獨(dú)立開采對(duì)地表下沉的終態(tài)影響，按照動(dòng)態(tài)預(yù)計(jì)原理，如果將各單元的終態(tài)影響乘以對(duì)應(yīng)的時(shí)間函數(shù)值，即可得指定預(yù)計(jì)時(shí)刻的地表動(dòng)態(tài)下沉預(yù)計(jì)公式，見式(8)。

W(x,t)=Φ(t)[W(x-s3)-W(x-v1t1)]+
Φ(t-t1)[W(x-v1t1)-W(x-v1t1-1v2t2)]+
Φ(t-t1-t2)[W(x-v1t1-v2t2)-W(x-v1t1-
v2t2-v3t3)]+…+Φ(t-t1-t2-…-tn-1)
{W(x-v1t1-…-vn-1tn-1)-W[x-
(v1t1+…+vn-1tn-1+vntn-s4)]}

(8)

如果每個(gè)單元的開采速度和開采時(shí)間均相等，即v1=v2=…=vn=vd，t1=t2=…=tn=td，式(8)可進(jìn)一步化簡。同理，地表水平移動(dòng)、水平變形、傾斜、曲率等，可根據(jù)相應(yīng)概率積分公式參照求得。

3 走向主斷面動(dòng)態(tài)預(yù)計(jì)算法

3.1 算法基本思想

動(dòng)態(tài)預(yù)計(jì)需考慮的問題比靜態(tài)預(yù)計(jì)更為復(fù)雜，首先是動(dòng)態(tài)單元的劃分問題，其次是，在預(yù)計(jì)時(shí)刻各單元開采狀態(tài)的確定問題(存在3種情況：動(dòng)態(tài)單元已開采完畢、部分開采完畢、沒有開采)。時(shí)間函數(shù)選定后，動(dòng)態(tài)單元長度的大小對(duì)預(yù)計(jì)精度也有較大影響，計(jì)算模型中通常以0.1H0(H0為平均采深)作為動(dòng)態(tài)單元的長度[17]。如果以平均速度代入計(jì)算，則每個(gè)動(dòng)態(tài)單元的長度除了最后一個(gè)，其余全部相等。計(jì)算步驟如下：

1)步驟1：在給定的預(yù)計(jì)時(shí)刻T，判斷1～n個(gè)動(dòng)態(tài)單元的開采狀態(tài)，即：已經(jīng)開采、部分開采或均未開采。當(dāng)T比開采所有煤層所需的總時(shí)間TZ大時(shí)，則表明所有單元均已采完。

2)步驟2：對(duì)于已采單元，計(jì)算其在T時(shí)刻對(duì)應(yīng)的時(shí)間函數(shù)值，由于模型采用的是分段時(shí)間函數(shù)，需判斷T與參數(shù)τ的大小，如后者大，采用分段函數(shù)第一段計(jì)算，如前者大，則采用第2段計(jì)算。

3)步驟3：采用式(6)計(jì)算各單元開采對(duì)地表移動(dòng)變形的靜態(tài)影響，再將其乘以對(duì)應(yīng)的時(shí)間函數(shù)值即可得到T時(shí)刻各單元開采的下沉預(yù)計(jì)值。

4)步驟4：將第3步各單元預(yù)計(jì)值進(jìn)行疊加求和，可得到預(yù)計(jì)時(shí)刻已采單元對(duì)地表下沉的總影響。

5)步驟5：繪制走向主斷面移動(dòng)變形圖，為了對(duì)比，可改變T值重新計(jì)算，將下沉，傾斜等分別繪制在同一幅圖中，觀察地表移動(dòng)變形的動(dòng)態(tài)發(fā)展規(guī)律。

3.2 算法的實(shí)現(xiàn)

在將預(yù)計(jì)模型轉(zhuǎn)換為計(jì)算機(jī)程序的過程中，需考慮各動(dòng)態(tài)單元時(shí)間函數(shù)值的計(jì)算方法問題，在計(jì)算時(shí)，如何判斷各單元的時(shí)間函數(shù)值應(yīng)采用第1或第2段函數(shù)進(jìn)行計(jì)算是難點(diǎn)，具體算法見表1。

表1 優(yōu)化分段時(shí)間函數(shù)編程算法Table 1 Programming algorithm of optimized segmented Knothe time function

上述算法解決了各動(dòng)態(tài)單元時(shí)間函數(shù)值的求取問題。除了計(jì)算時(shí)間函數(shù)值，在走向主斷面的動(dòng)態(tài)預(yù)計(jì)模型中，還需計(jì)算各單元對(duì)地表移動(dòng)的靜態(tài)影響，難點(diǎn)在于考慮拐點(diǎn)偏移距后，各單元起始坐標(biāo)如何確定，相應(yīng)算法見表2。

表2 走向主斷面動(dòng)態(tài)預(yù)計(jì)算法Table 2 Dynamic prediction algorithm of strike main section

4 走向主斷面動(dòng)態(tài)預(yù)計(jì)實(shí)例

4.1 預(yù)計(jì)實(shí)例1

據(jù)文獻(xiàn)[19]，某礦工作面1002的靜態(tài)預(yù)計(jì)參數(shù)如下：D3=1 000 m，D1=250 m，m=3 m，q=0.78，H0=500 m，tanβ=2.0，s3=s4=50 m，b=0.28。動(dòng)態(tài)參數(shù)v=5.0 m/d，動(dòng)態(tài)單元長度L=0.1H0，時(shí)間參數(shù)c和τ，按文獻(xiàn)[16]中方法求得。用本文模型與算法編制的程序，在不同的預(yù)計(jì)時(shí)刻T，對(duì)1002工作面開采時(shí)地表走向主斷面下沉和傾斜進(jìn)行動(dòng)態(tài)預(yù)計(jì)，結(jié)果如圖3所示。由于在傾向上為非充分開采，需考慮非充分采動(dòng)因素，故地表實(shí)際下沉比充分采動(dòng)最大下沉量W0(2.38 m)小很多，從圖中還可看出地表下沉的動(dòng)態(tài)變化規(guī)律。當(dāng)T=500 d時(shí)，地表下沉和傾斜很小，這是由于此時(shí)工作面僅開采了210 m，并且采深較大的緣故。

圖3 工作面1002走向主斷面下沉與傾斜動(dòng)態(tài)預(yù)計(jì)Fig.3 Dynamic prediction of surface subsidence and inclination in strick principal section of No.1002 working face

隨著T增加，地表下沉和傾斜不斷增大，影響范圍也逐漸擴(kuò)大。當(dāng)T=900 d時(shí)，距離工作面開采結(jié)束又歷時(shí)700 d，地表移動(dòng)已趨于穩(wěn)定，再增加T，地表下沉和傾斜曲線不會(huì)再發(fā)生任何改變，此時(shí)，所有動(dòng)態(tài)單元對(duì)應(yīng)的時(shí)間函數(shù)值均為1，代表著各單元開采影響達(dá)到了最大值。由于計(jì)算模型中考慮了拐點(diǎn)偏移距的影響，圖3中的下沉曲線拐點(diǎn)出現(xiàn)在實(shí)際拐點(diǎn)的正上方，而預(yù)計(jì)的傾斜值，在拐點(diǎn)處也達(dá)到了最大值，這與理論所揭示的規(guī)律相吻合。

4.2 預(yù)計(jì)實(shí)例2

為了驗(yàn)證模型精度的可靠性，另對(duì)官地礦29401工作面開采走向主斷面進(jìn)行動(dòng)態(tài)預(yù)計(jì)，工作面相關(guān)概率積分參數(shù)見參考文獻(xiàn)[10]，動(dòng)態(tài)預(yù)計(jì)時(shí)間參數(shù)的計(jì)算方法同實(shí)例1。預(yù)計(jì)時(shí)間節(jié)點(diǎn)與實(shí)際觀測時(shí)間相對(duì)應(yīng)，實(shí)測9次即進(jìn)行了9次預(yù)計(jì)。由于走向監(jiān)測點(diǎn)較多，預(yù)計(jì)后，抽樣對(duì)其第5次和第9次預(yù)計(jì)結(jié)果與實(shí)測結(jié)果進(jìn)行對(duì)比分析，并統(tǒng)計(jì)精度，見表3。

提取最大下沉點(diǎn)的9次預(yù)計(jì)結(jié)果，將其與實(shí)測結(jié)果進(jìn)行對(duì)比分析，并統(tǒng)計(jì)其精度，見表4。

在開采沉陷研究中，預(yù)計(jì)結(jié)果精度常采用標(biāo)準(zhǔn)差m和相對(duì)誤差f來衡量[10,20]，可分別用式(9)和式(10)進(jìn)行計(jì)算。

(9)

(10)

通過表4可知，通過對(duì)走向監(jiān)測點(diǎn)的抽樣分析，其第5次預(yù)計(jì)的標(biāo)準(zhǔn)差為271.8 mm，相對(duì)誤差為5.6%，第9次預(yù)計(jì)標(biāo)準(zhǔn)差為270.4 mm，相對(duì)誤差為5.4%。

通過表3可知，在對(duì)最大下沉點(diǎn)的動(dòng)態(tài)預(yù)計(jì)中，前3期所得到的預(yù)測與實(shí)測結(jié)果較為接近，相差最大約186 mm，預(yù)測誤差較小，后6期的預(yù)測結(jié)果小于實(shí)測結(jié)果，預(yù)測誤差相對(duì)較大。經(jīng)統(tǒng)計(jì)，本次預(yù)計(jì)的最大下沉點(diǎn)預(yù)計(jì)中誤差約為±303 mm，預(yù)計(jì)相對(duì)誤差則在5.7%左右，預(yù)計(jì)精度較為穩(wěn)定。

5 結(jié) 論

1)優(yōu)化分段Knothe時(shí)間函數(shù)完善了原函數(shù)的不足，使其具有更好的適應(yīng)性。以該函數(shù)為基礎(chǔ)，推導(dǎo)了各動(dòng)態(tài)單元時(shí)間函數(shù)值的計(jì)算公式，結(jié)合概率積分法，建立了適應(yīng)于水平和緩傾斜煤層開采的地表走向主斷面沉降變形動(dòng)態(tài)預(yù)計(jì)模型。

2)根據(jù)動(dòng)態(tài)預(yù)計(jì)原理和預(yù)計(jì)模型，探討了動(dòng)態(tài)預(yù)計(jì)計(jì)算的詳細(xì)流程，設(shè)計(jì)了時(shí)間函數(shù)值和走向主斷面沉降變形計(jì)算的編程算法，并編制了程序。

3)1002工作面預(yù)計(jì)實(shí)踐表明：當(dāng)給定的預(yù)計(jì)時(shí)間足夠大，動(dòng)態(tài)預(yù)計(jì)結(jié)果與靜態(tài)預(yù)計(jì)結(jié)果相吻合，并在拐點(diǎn)偏移距處的傾斜值達(dá)到了理論最大值。29401工作面預(yù)計(jì)結(jié)果表明：其走向主斷面地表點(diǎn)動(dòng)態(tài)預(yù)計(jì)精度約在6%以內(nèi)，證明了本文預(yù)計(jì)模型和算法具有較高的精度。