張 兵，崔希民，趙玉玲,3，李春意

(1.石家莊學院 資源與環(huán)境科學學院,河北 石家莊 050035; 2.中國礦業(yè)大學(北京) 地球科學與測繪工程學院,北京 100083; 3.河北工程大學 礦業(yè)與測繪工程學院,河北 邯鄲 056038; 4.河南理工大學 測繪與國土信息工程學院，河南 焦作 454000)

采動地表建、構物的損害主要歸因于地表點的移動和變形，經(jīng)典的開采沉陷理論與方法主要集中在靜態(tài)預計研究領域[1-3]，但由于地表移動、變形是一個復雜的動態(tài)過程，與開采速度，煤層采深和傾角、頂板管理方法，上覆巖層力學性質等因素密切相關[4-7]，掌握地表點隨時間的動態(tài)變化規(guī)律同樣具有重要意義，能夠更準確的分析采動地表移動變形與采礦地質因素的關系，可為采后損害鑒定、采前建筑物保護、采空區(qū)土地再利用等供技術支撐。在實踐中，國內外最主要的動態(tài)預計方法多是通過用靜態(tài)預計公式乘以相應的時間函數(shù)來實現(xiàn)的[8-10]，我國開采沉陷動態(tài)預計始于20世紀80年代，許多學者做了大量工作。崔希民、彭小沾等[11-12]在動態(tài)預計研究中，對Knothe時間函數(shù)進行了深入的研究，給出了5種時間系數(shù)的確定方法，同時指出了該時間函數(shù)的不足，并給出了理想的時間函數(shù)分布形態(tài)。LUO Yi和CHENG Jianwei[13]通過研究美國長壁開采地表變形規(guī)律，對Knothe時間函數(shù)進行改進，使其適應傾斜煤層開采時的動態(tài)預計。胡青峰等[14]進一步闡述了Knothe時間系數(shù)的影響因素，探討了時間系數(shù)的直接求取方法。張凱等[15]采用生長函數(shù)模型對正態(tài)分布時間函數(shù)進行了優(yōu)化，解決了其密度函數(shù)的有效積分域虧損隨時間參數(shù)c而減小的問題。郭旭煒等[16]給出了分段Konthe時間函數(shù)的動態(tài)求參方法，解決了其參數(shù)不隨時間變化的不足。陳磊等[17]采用水準與 InSAR技術相結合估算冪指數(shù) Knothe 時間函數(shù)模型的參數(shù)，并用實例說明使用該方法可提高動態(tài)預計精度。王軍保等[18]采用巖石力學的非定常流變模型，將Knothe時間系數(shù)看作非常量，構建了一種新的動態(tài)沉陷預計公式。李懷展等[19]提出了一種基于加權法和地質力學參數(shù)敏感性的地表沉降動態(tài)預測方法，并用實例證明了該方法可以提高地表動態(tài)沉降的預測精度。楊澤發(fā)等[20]建立了Logistic模型與InSAR觀測值之間的函數(shù)關系，然后利用InSAR觀測值估計Logistic模型參數(shù)，進而再采用Logistic模型進行動態(tài)預計。很多研究者還給出了不同的動態(tài)預計方法[21-23]，如地表動態(tài)沉降預測的Richards模型、基于雙因素時間函數(shù)的地表點動態(tài)沉降預計模型等。

通過對已有文獻分析可知，現(xiàn)有動態(tài)沉陷預計模型多是純理論的，沒有結合足夠的工程應用，且大多數(shù)預計模型的參數(shù)求取非常不易，推廣應用較為困難。通過分析還可知，Knothe時間函數(shù)在動態(tài)預計研究中占有重要地位，同時，不少學者也指出了該函數(shù)的應用局限性。為此，很多學者對其進行過改進[24-26]，筆者也曾在Knothe時間函數(shù)的基礎上構建了“優(yōu)化分段Knothe時間函數(shù)模型”[27](以下簡稱:“優(yōu)化分段時間函數(shù)”)，提高了動態(tài)預計精度、擴展了該函數(shù)的適應性;另外，在走向主斷面動態(tài)預計研究方面，很多學者做過研究，給出了相應的方法，但在傾向主斷面動態(tài)沉陷預計方面，現(xiàn)有的研究涉及不多，但筆者認為對其進行研究很有必要，相較于對整個下沉盆地的計算而言，研究主斷面上的動態(tài)預計模型及算法，可以大大節(jié)省計算時間，快速得到預計結果。

筆者以“優(yōu)化分段時間函數(shù)”為基礎，以概率積分模型為影響函數(shù)，重點探討緩傾斜矩形或近矩形工作面開采地表傾向主斷面動態(tài)預計方法問題，并設計相應算法，為動態(tài)預計的多方面工程應用提供解決方案。

1 傾向主斷面動態(tài)預計模型構建

1.1 傾向主斷面沉陷預計基本原理

設煤層在走向上達到了充分采動，傾向上為有限開采，傾向主斷面上的地表沉降變形可采用等影響原理來計算[1]，具體可采用圖1進行輔助說明。圖1中，AB為實采邊界;s1為下山拐點偏距;s2為上山拐點偏距;CD為AB煤層開采時的計算邊界。

在計算公式建立的過程中，選擇圖1中過A點的直線和地表交點(P)作為坐標原點，垂直指向下方表示地表下沉(W),而地表點的坐標可用y軸表示，y軸為沿地表水平指向右方,Wy為y軸方向的下沉。考慮到煤炭開采拐點偏移距的影響，如果開采AB煤層，根據(jù)等影響原理，AB開采所引起的地表傾向主斷面上點的下沉可用W0(y)表示，采用式(1)進行計算[1]:

圖1 有限開采地表傾向主斷面預計等影響計算原理Fig.1 Equal effect principle in inclined principal section prediction when finite mining

W0(y)=W(y-LAC′;t1)-W(y-LAC′-L;t2)

(1)

LAC′=(H1-s1sinα)cotθ0-s1cosα

(2)

(3)

(4)

其中，LAC′為AC′的距離；t1，t2為計算時所對應的上下山邊界參數(shù);α為煤層傾角;θ0為開采影響傳播角;H1和H2為下山和上山采深;D1為AB煤層長度;W0為最大下沉值;y為地表點坐標。式(1)即為AB煤層采動引起的地表下沉靜態(tài)(終態(tài))計算式。

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

時間函數(shù)是動態(tài)預計的重要組成部分，可用于動態(tài)預計的時間函數(shù)有很多，選擇什么樣的時間函數(shù)，在一定程度上決定著動態(tài)預計精度的高低，本文采用 “優(yōu)化分段時間函數(shù)”[27]，該時間函數(shù)是筆者曾改進過的，并用實踐證明了它的可靠性，改進后的時間函數(shù)提高了原函數(shù)的預測精度及對不同地質采礦條件的適應性。另外，筆者還曾對其參數(shù)的確定方法進行了詳細討論[28]，推導了其時間參數(shù)的直接法計算公式，方便了計算編程，“優(yōu)化分段時間函數(shù)”模型如式(5)所示，其函數(shù)圖像如圖2所示:

(5)

式中，t為任意給定的預計時刻；τ為最大下沉速度出現(xiàn)的時刻；c為時間系數(shù)，與覆巖力學性質相關;T為下沉總時間。

圖2 優(yōu)化分段時間函數(shù)圖像形態(tài)Fig.2 Segmented optimization time function

1.3 傾向主斷面動態(tài)預計模型建立

在動態(tài)預計中，通常要按照一定的原則將煤層劃分為不同的動態(tài)開采單元[29]。傾向主斷面動態(tài)開采單元的劃分和地表沉陷的動態(tài)發(fā)展規(guī)律可用圖3進行說明，各單元的長度用開采速度v乘以相應的開采時間t表示。由于各單元開采順序不同，在預計時刻其對地表點的影響時間也不同，先開采的單元影響大，后開采的單元影響小。第n個單元開采后，第1個單元采后所經(jīng)歷的時間是t1,t2，…，tn之和，其對地表的影響最大，第2個單元采后所經(jīng)歷的時間是t2,t3,…,tn之和，其影響次之，最后1個單元由于剛采完畢，其影響未波及地表。假設第1個動態(tài)單元對地表的下沉影響用W1(y)表示，第2個單元影響用W2(y)表示，第n個單元影響用Wn(y)表示。

圖3 傾向動態(tài)單元開采對地表點下沉的動態(tài)影響Fig.3 Dynamic influence when mining strike section dynamic unit on surface subsidence

各動態(tài)單元開采對地表沉陷的影響，也符合有限開采原理，可按有限開采相關公式計算，但需要改變相應的概率積分參數(shù)。根據(jù)有限開采傾向主斷面下沉預計公式，下面給出從第1個單元到第n個單元開采的地表下沉計算公式

(6)

(7)

(8)

其中，b1為下山方向水平移動系數(shù);b2為上山方向水平移動系數(shù);tanβ1為下山方向主要影響角正切;tanβ2為上山方向主要影響角正切。以上均指第1個動態(tài)開采單元的相關參數(shù)。

第2個動態(tài)單元開采時，參照式(6)～(8)，可給出W2(y)為

(9)

式(9)中涉及到的各參數(shù)計算方法為

(10)

(11)

同理Wn(y)可由上述推導類比得到，具體如下:

(12)

第n個動態(tài)單元相關參數(shù)可由式(13)和(14)求取:

(13)

W1(y)，W2(y)，…，Wn(y)是計算1～n個動態(tài)單元開采后單獨對地表的終態(tài)(靜態(tài))影響，如果計算某個時刻T的地表動態(tài)下沉，需要將每個單元的靜態(tài)影響值乘以各單元在T時刻對應的時間函數(shù)值，再進行疊加求和，即可得到T時刻的地表下沉值，計算方法如式(15)所示。需要強調的是:在計算W1(y)，W2(y)，…，Wn(y)過程中均考慮了拐點偏距的影響，坐標原點設在了圖1中的P點處。

W(y,t)=Φ(t)W1(y)+Φ(t-t1)W2(y)+…+

Φ(t-t1-t2-…-tn-1)Wn(y)=

(15)

預計T時刻的地表變形如:曲率、傾斜、水平變形和水平移動，可根據(jù)相應的概率積分法計算公式，參照上述下沉公式的推導過程得出，限于篇幅，這里不再重復。

2 傾向主斷面動態(tài)預計算法

地表移動變形計算公式建立之后，為了能夠用計算機語言來實現(xiàn)計算過程，編制預計程序，則需要研究相應的算法，在傾向主斷面動態(tài)預計中要考慮每個動態(tài)單元的相應參數(shù)，如:各單元水平移動系數(shù)，上下山采深及主要影響半徑等。當預計精度不高時，可帶入各參數(shù)的平均值進行計算，如果直接將原有的概率積分參數(shù)代入計算，預計精度則會大大降低。另外，分別計算各動態(tài)單元參數(shù)不但能提高預計精度，還可提高本文模型的應用范圍，擴展其對于較大傾角煤層開采的適用性。

2.1 算法基本思想

計算各個單元的y坐標是算法中較為復雜的部分，開采傳播角的影響使得下沉和變形拐點向下山方向移動。如圖1所示，可以將各動態(tài)單元的計算長度以L為基礎按比例內插求取，動態(tài)移動變形計算步驟如下:

Step1:計算走向方向充分采動程度系數(shù)(Cxm)，在動態(tài)計算結果中乘以該系數(shù)，可以使得計算結果更為精確，也與實際情況相符合。

Step2:計算預計T時刻各單元所對應的時間函數(shù)值(PHT)。這是計算各單元對地表影響值大小的關鍵參數(shù)，表明了各單元的終態(tài)影響在T時刻對應的時間調整系數(shù)。

Step3:確定動態(tài)預計的計算坐標系，即確定坐標原點、下沉軸和y軸(用以表示地表點的位置);計算上山和下山方向各單元對應的y值大小，這是概率積分法的主要參數(shù)之一。

Step4:計算預計T時刻各單元的上下山采深、主要影響半徑、主要影響角正切，以便代入各動態(tài)單元對應的公式進行計算。

Step5:對每個動態(tài)單元按照公式W1(y)，W2(y)，…，Wn(y)進行計算，得到各單元影響靜態(tài)下沉值，計算完成后分別乘以對應的時間函數(shù)值PHIT)，再求和，即可得到預計T時刻的動態(tài)下沉。

2.2 算法的實現(xiàn)

計算模型和算法的難點在于:考慮拐點偏移距時，1～n個動態(tài)單元對應的計算長度Li如何確定？它是概率積分模型主要參數(shù)之一，要計算Li，需先確定各單元所對應的上、下山方向的y值;另外，在工作面概率積分參數(shù)已知時，各動態(tài)單元相應參數(shù)如何計算，也是需要重點注意的問題。傾向主斷面動態(tài)預計算法偽代碼見表1。

3 傾向主斷面動態(tài)預計實例

3.1 動態(tài)預計實例1

以文獻[1]中介紹的某工作面為例，對其進行傾向主斷面動態(tài)預計，該工作面具體情況如下:走向長1 000 m,傾向長250 m，上山采深H2=279 m，下山采深H1=322 m，煤層傾角=12°，下沉系數(shù)q=0.78，采厚m=1.45 m，主要影響角正切tanβ1=2.2，tanβ2=2.0，開采影響傳播角θ0=81.6°，水平移動系數(shù)b1=0.36，b2=0.30，拐點偏距s1=30 m，s2=28 m;覆巖巖性為中硬，用垮落法管理頂板。動態(tài)預計參數(shù)為:動態(tài)單元長度采用有效尺寸分割法[29]確定，確定為0.1H0(H0為工作面平均采深)，開采速度v為5.0 m/d，時間函數(shù)參數(shù)τ和c按照筆者在文獻[28]中所介紹的方法自動計算，動態(tài)預計程序按本文模型及算法編制。

表1 傾向主斷面動態(tài)預計算法偽代碼Table 1 Pseudocode of dynamic subsidence prediction algorithm for inclined main section

如果工作面是“從上山向下山方向開采”，在不同的預計時刻對工作面傾向主斷面進行動態(tài)預計，根據(jù)預計結果，可繪制如圖4所示的動態(tài)下沉和傾斜曲線。

由圖4可知，由于開采速度較小，采深又較大，當T=200 d時，工作面實際開采了100 m，由于頂板的懸臂支撐，此時的下沉和傾斜都很小。隨著T值的增大，已開采的動態(tài)單元數(shù)越來越多，地表的沉陷范圍越來越大，受開采影響的地表點的下沉和傾斜也逐漸增加。如果給定時間足夠大，如在T=800 d時進行動態(tài)預計，地表的下沉和傾斜與T=1 200 d預計的下沉曲線非常接近，因此可認為T=1 200 d時的預計結果為靜態(tài)結果，即地表趨于穩(wěn)定時的終態(tài)下沉曲線，也說明T=800 d時地表移動變形基本就已經(jīng)趨于了穩(wěn)定。

圖4的“動態(tài)預計-下沉”中，黑色曲線是地表下沉的終態(tài)曲線，其拐點不在計算邊界垂直上方，也不在開切眼的垂直上方，而是距離開切眼有一定的距離，向著下山方向偏移，位于圖中O點處，這是因為受到了開采影響傳播角的影響。

為了驗證程序的理論正確性，需要進行對比研究，假設“工作面再從下山向上山方向開采”，在對應的時刻再次進行動態(tài)預計，根據(jù)預計結果，則可得到如圖5所示的地表下沉和傾斜圖。

將不同開采方式的2圖進行對比可知，圖4中地表點的下沉和傾斜是從左向右進行傳播的，而圖5則正好相反，這說明，開采順序不同，地表各點出現(xiàn)最大下沉和傾斜值的時刻是不同的,在對地表建筑物進行保護的過程中，判斷地表點移動和變形最大值出現(xiàn)的時刻和位置至關重要。通過對比還可知，無論開采方式是選擇上山開采還是選擇下山開采，傾向主斷面上各點的終態(tài)下沉是完全相同的，即，圖4,5中，T=1 200 d時的曲線是完全重合的，此時無論給定的預計時刻T如何增加，受影響地表點的范圍和移動變形值都不會再發(fā)生任何改變，在理論上說明了模型和算法的科學性。

圖4 上山開采傾向主斷面下沉和傾斜動態(tài)預計Fig.4 Dymamic subsidence and inclination when mining along downhill direction

圖5 下山開采傾向主斷面下沉和傾斜動態(tài)預計Fig.5 Dymamic surface subsidence and inclination when mining along uphill direction

3.2 動態(tài)預計實例2

為驗證本文模型和算法精度的可靠性，以官地礦 29401采煤工作面傾向主斷面為研究對象，對其開采進行動態(tài)預計，對比預計結果和實測數(shù)據(jù)，統(tǒng)計預計精度。

3.2.1動態(tài)預計

官地29401工作面基本情況見文獻[28]，在此不再列舉。為簡化編程算法，在預計前先將實際工作面和地表監(jiān)測點轉換為以工作面走向為橫軸的坐標系，如圖6所示。

圖6 工作面與測點分布Fig.6 Distribution of monitoring points and working face

根據(jù)29401工作面傾向主斷面9次觀測結果繪制的下沉曲線如圖7所示。

圖7 傾向觀測線實測下沉量Fig.7 Measured subsidence of tend towards line

利用本文算法編制的動態(tài)預計程序，在與實際觀測對應的時刻，進行了9次動態(tài)預計，得到的下沉預計結果如圖8所示。

3.2.2預計精度分析

在開始階段，地表的移動變形值都很小，本文抽樣選擇了第4,5,7,9次觀測成果，將其與動態(tài)預計結果進行對比來統(tǒng)計預測精度。實測與預測結果對比如圖9所示。

圖8 傾向觀測線預計下沉量Fig.8 Prediction subsidence of tend towards line

采用文獻[14]和[27]所介紹的方法，對抽樣的4次預測與實測結果進行統(tǒng)計分析，得到的統(tǒng)計數(shù)據(jù)見表2。由表2可知，預測絕對中誤差最小為±54 mm，最大為±452 mm，預測相對中誤差除了第7次較大外，其他相對誤差分布在8.5%左右，在動態(tài)預計中精度相對較高，能夠滿足絕大多數(shù)工程需要。

圖9 第4,5,7,9次預測下沉與實測下沉對比Fig.9 Comparison between predicted subsidence and measured subsidence at the 4th,5th,7th and 9th

表2 動態(tài)預計精度統(tǒng)計Table 2 Statistical table of dynamic prediction accuracy

4 結 論

(1)以“優(yōu)化分段Knothe時間函數(shù)”為基礎，結合概率積分法，建立了適應于緩傾斜、矩形或近似矩形工作面開采時的傾向主斷面地表沉陷動態(tài)預計模型。

(2)根據(jù)動態(tài)預計原理和所建模型，給出了傾向主斷面動態(tài)預計詳細流程，設計了對應的計算機編程算法，計算流程清晰明了，算法簡潔。

(3)預計實踐1表明:采用本文模型及算法，預計的地表動態(tài)移動變形規(guī)律與理論所揭示的移動過程相符合，證明了模型及算法的科學性;預計實踐2表明:29401工作面傾向觀測線動態(tài)預計最大相對誤差在整體上可以控制在8.5%左右，證明了模型在預計精度上的可靠性。