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

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

0 引 言

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

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

1 動態(tài)預(yù)計基本原理

1.1 動態(tài)單元劃分

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

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

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

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

1.2 動態(tài)預(yù)計時間函數(shù)及其計算

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

(1)

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

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

第2個動態(tài)開采單元時間函數(shù)值計算公式為

(2)

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

(3)

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

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

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

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

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

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

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

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

(4)

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

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

(5)

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

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

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

(6)

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

(7)

式中，W0為最大下沉量。

式(4)—式(6)計算的是各動態(tài)單元獨立開采對地表下沉的終態(tài)影響，按照動態(tài)預(yù)計原理，如果將各單元的終態(tài)影響乘以對應(yīng)的時間函數(shù)值，即可得指定預(yù)計時刻的地表動態(tài)下沉預(yù)計公式，見式(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)

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

3 走向主斷面動態(tài)預(yù)計算法

3.1 算法基本思想

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

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

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

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

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

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

3.2 算法的實現(xiàn)

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

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

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

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

4 走向主斷面動態(tài)預(yù)計實例

4.1 預(yù)計實例1

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

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

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

4.2 預(yù)計實例2

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

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

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

(9)

(10)

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

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

5 結(jié) 論

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

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

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