徐孟強　查劍鋒　王金濤

(1.中國礦業(yè)大學環(huán)境與測繪學院　江蘇徐州 221116；2.國土環(huán)境與災害監(jiān)測國家測繪地理信息局重點實驗室　江蘇徐州 221116)

?

基于鉆孔數(shù)據(jù)煤礦開采地表沉陷預計

徐孟強1，2查劍鋒1，2王金濤1，2

(1.中國礦業(yè)大學環(huán)境與測繪學院江蘇徐州 221116；2.國土環(huán)境與災害監(jiān)測國家測繪地理信息局重點實驗室江蘇徐州 221116)

開采沉陷引起的環(huán)境損害評價是環(huán)評工作中的重要內(nèi)容，現(xiàn)有預計軟件基本不適應在大面積環(huán)評工作中開展沉陷預計工作。針對現(xiàn)有問題，以礦區(qū)鉆孔數(shù)據(jù)為基礎(chǔ)，基于Delaunay三角形剖分算法建立了煤層開采地表沉陷預計模型。通過設計對比試驗，驗證了新建模型預計結(jié)果較傳統(tǒng)的預計方法更加符合實際地表沉陷情況，預計結(jié)果的可靠性較高。將該模型運用到工程實例中，可以節(jié)省人力資源，使煤礦生態(tài)環(huán)評從原來以定性分析、半定量預計為主的狀況，提高到以定量預測為主的高度，更好地為礦區(qū)生態(tài)環(huán)境保服務。

環(huán)境評價Delaunay三角形三角形分割概率積分法開采沉陷預計實驗設計

0　引言

開采沉陷引起的環(huán)境損害評價是環(huán)評工作中的重要內(nèi)容。礦區(qū)開展環(huán)評工作，預計煤礦開采引起的地表沉陷通常采用以下方法：根據(jù)鉆孔數(shù)據(jù)揭示的煤層厚度求取礦區(qū)內(nèi)的平均煤層厚度進行預計或者將礦區(qū)人為地分割成多個等厚薄煤層進行預計[1]。前者處理簡單，預計結(jié)果精度較差；后者預計結(jié)果精度較高，但是費時費力，同時存在多次分割后工作面邊界角點坐標難以確定等問題。為了更好地處理礦區(qū)環(huán)評工作中開采沉陷預計，以概率積分法原理為基礎(chǔ)，結(jié)合Delaunay三角形剖分算法和三角形分割模型設計了基于鉆孔數(shù)據(jù)煤礦開采地表沉陷預計程序。

1　Delaunary 三角形剖分算法原理介紹

Delaunay三角網(wǎng)是在泰森多邊形的基礎(chǔ)上發(fā)展而來的。1908年，G Voronoi首先在數(shù)學上限定了每個離散點數(shù)據(jù)的有效作用范圍，并定義了二維平面上的泰森多邊形[2]。經(jīng)過不斷發(fā)展，1934年，B Delaunay由泰森多邊形演化出了更易于分析、應用更加廣泛的Delaunay三角網(wǎng)[3]。目前Delaunay三角網(wǎng)已經(jīng)發(fā)展成為最為高效的離散數(shù)據(jù)處理方法，廣泛應用于計算機圖形學、自然科學、地學等眾多科學領(lǐng)域?；贒elaunay三角網(wǎng)剖分提出了諸多的剖分算法，主要包括：分割歸并算法，逐點插入算法、三角網(wǎng)增長算法[4]。因逐點插入算法理論嚴密、唯一性好，算法實現(xiàn)比較簡單，基于鉆孔數(shù)據(jù)的煤層開采地表沉陷預計，借助逐點插入算法根據(jù)工作面角點和鉆孔數(shù)據(jù)生成Delaunay三角網(wǎng)。具體構(gòu)網(wǎng)步驟可參考文獻[5-6]。

2　三角形分割算法原理

根據(jù)工作面角點和鉆孔數(shù)據(jù)，基于逐點插入算法生成Delaunay三角網(wǎng)，此時，工作面區(qū)域被劃分為多個獨立的三角形單元，預計工作面開采地表移動變形就是將預計每個三角形單元開采引起的地表移動變形值的疊加。根據(jù)概率積分法基本原理及其計算公式[1]可知，預計工作面開采引起的地表移動變形主要是進行二重積分計算。所以，地表移動變形預計可以采用二重變步長辛卜生算法[7-8]積分直接求取。因此確定積分限就成了關(guān)鍵的問題。如圖1所示，從生成的Delaunay三角網(wǎng)中任意取出一個三角形單元，過角點B作垂線交AC邊于點D，三角形ABC被分為2個小三角形單元，三角形ABD，DBC。此時，積分上下限的確定就比較簡單了。

在三角形ABD中，已知坐標A(xA，yA)，B(xB，yB)，D(xD，yD)，則x軸積分上下限分別為xA和xB；y軸積分上下限為：

(1)

(2)

同理可以得到三角形DBC在y軸方向上積分上下限。

圖1　三角形分割

3　系統(tǒng)實現(xiàn)流程

以概率積分法原理為基礎(chǔ)，結(jié)合Delaunay三角形剖分算法和三角形分割模型設計了基于鉆孔數(shù)據(jù)的煤層開采地表沉陷預計程序。系統(tǒng)實現(xiàn)流程見圖2。通過參數(shù)輸入模塊將地質(zhì)采礦條件，概率計分方法預計參數(shù)以及相應的鉆孔數(shù)據(jù)輸入到軟件系統(tǒng)中，進行數(shù)據(jù)預處理?？紤]到礦區(qū)使用的坐標系往往是大地坐標系或者地方坐標系，而開采沉陷預計使用的是工作面坐標系，所以需要將礦區(qū)坐標系下的工作點角點坐標和鉆孔坐標轉(zhuǎn)換成工作面坐標系下的坐標。根據(jù)轉(zhuǎn)換后的角點坐標，尋找多邊形重心，利用文獻[9]的計算公式對工作面計算邊界進行調(diào)整。根據(jù)調(diào)整后的工作面邊界將工作面外圍的鉆孔數(shù)據(jù)進行剔除，利用逐點插入算法生成Delaunay三角網(wǎng)。為便于積分運算，將Delaunay三角形劃分為2個獨立的小三角形，采用二重變步長辛卜生算法逐一計算每個小三角形單元開采引起的地表移動變形值，最終將所有三角形單元的預計結(jié)果進行疊加，獲得整個工作面開采地表移動變形值。系統(tǒng)可以將預計結(jié)果以文本形式輸出，利用Visual Basic和surfer混合編程，通過Visual Basic調(diào)用surfer的ActiveX自動化對象，實現(xiàn)預計結(jié)果的可視化輸出，為采煤沉沉陷分析、地表損害評估提供更好地技術(shù)支持。需要指出的是，非特殊地質(zhì)采礦條件下，煤層厚度具有漸變性，相鄰鉆孔揭示的煤層厚度差異較小，因此，煤礦開采沉陷預計時，取3個角點厚度的均值作為單個三角形單元開采的等效開采厚度，在保證計算結(jié)果精度和可靠性的同時，有效降低計算難度。

圖2　系統(tǒng)實現(xiàn)流程

4　預計結(jié)果精度與可靠性檢驗

為了測試系統(tǒng)的精度與可靠性，設計了3組試驗。將3組試驗的預計結(jié)果進行對比分析，驗證系統(tǒng)對于基于鉆孔數(shù)據(jù)的煤層開采地表移動變形預計結(jié)果的準確性及可靠性。

4.1試驗設計及成果處理

試驗1：假設2個不規(guī)則工作面A，B。在工作面A，B上方相同位置每間隔100 m布設1個監(jiān)測點，共布設13個監(jiān)測點(見圖3)。工作面A，B角點坐標相同，地質(zhì)采礦條件和概率積分法預計參數(shù)相同，詳細參數(shù)見表1。工作面B開采范圍內(nèi)額外存在12個鉆孔點。利用編制的預計軟件直接預計工作面A開采引起的地表沉陷以及工作面上方監(jiān)測點的下沉值。首先根據(jù)工作面角點坐標和鉆孔數(shù)據(jù)生成Delaunay三角網(wǎng)，然后預計工作面B開采引起的地表移動變形。

試驗2：假設兩個不規(guī)則工作面A，B。工作面A，B角點坐標、工作面范圍內(nèi)鉆孔數(shù)量、鉆孔坐標等基本信息同試驗1。鉆孔數(shù)據(jù)揭示工作面A，B煤層厚度在1.6～2.2 m范圍內(nèi)變化，沿著工作面走向方向(圖上從左向右方向)，煤層厚度逐漸減少。為了預計工作面A，B開采地表移動變形值，對于工作面A，以鉆孔數(shù)據(jù)揭示的工作面A開采范圍內(nèi)的平均開采厚度(1.9 m)為基礎(chǔ)進行沉陷預計。對于工作面B，先根據(jù)工作面角點坐標和鉆孔數(shù)據(jù)生成Delaunay三角網(wǎng)，取三角形單元3個角點煤層厚度的均值，作為三角形單元的等效開采厚度來求解每個三角形開采引起的地表移動變形，將所有三角形單元開采預計結(jié)果進行疊加，得到工作面B開采引起的地表移動變形。

(a)工作面A　　　　　　　　　　　　　　　　(b)工作面B

試驗3：在試驗2的基礎(chǔ)上，將工作面A根據(jù)煤層采厚相近原則將開采區(qū)域劃分成多個小的計算域，本次試驗將工作面A劃分成7個計算區(qū)域。利用開采沉陷預計軟件預計工作面A開采地表下沉。

根據(jù)試驗1～3的設計方案，采用5種不同的預計方案，預計煤礦開采工作面上方地表下沉情況，根據(jù)預計結(jié)果，生成相應的下沉等值線圖，見圖4。將3種試驗方案工作面上方布設的監(jiān)測點預計的下沉值進行匯總，以工作面A的預計結(jié)果為標準，計算絕對誤差和相對誤差。

(a)試驗1　無Delaunay三角形分割

(b)試驗1　Delaunay三角形分割

(c)試驗2　無Delaunay三角形分割

(d)試驗2　Delaunay三角形分割

(e)試驗3　人為劃分

(f)試驗3　Delaunay三角形分割

4.2結(jié)果分析

從形態(tài)上分析，試驗1兩種預計方法預計的地表下沉值幾乎沒有差別。將工作面上方地表監(jiān)測點的下沉值進行對比分析，兩種預計方法最小絕對誤差為0 mm，最大絕對誤差僅為-6.6 mm，僅占相應點處下沉值的0.39%。在進行沉陷預計時，工作面A被分割成5個三角形單元，工作面B被分割成49個三角形單元。工作面劃分上的差異性，導致兩種預計結(jié)果存在較小誤差。上述結(jié)果表明等厚開采，是否將工作面進行Delaunay三角形分割對預計結(jié)果都不會產(chǎn)生明顯影響。

分析圖4(c)、圖4(d)可以明顯看出，兩種預計結(jié)果存在差異。圖4(c)表明地表沉陷極值范圍出現(xiàn)在走向方向中心處，圖4(d)表明地表沉陷的極值范圍出現(xiàn)在走向方向上方。試驗2煤層厚度在工作面走向方向上逐漸減少，所以地表沉陷的極值范圍應該出現(xiàn)在走向方向上方。將工作面上方地表監(jiān)測點的下沉值進行對比分析，最大絕對誤差達到了-131.6 mm，占相應點處下沉值的10.8%。剔除工作面分割的差異性對預計結(jié)果的影響，兩組結(jié)果差異仍較大。同時，相鄰鉆孔處揭示的煤層厚度差異較小，工作面B在進行開采沉陷預計時，根據(jù)工作面角點信息和鉆孔數(shù)據(jù)將工作面進行Delaunay三角形分割，取三角形單元3個角點煤層采厚的均值作為三角形單元開采的等效開采厚度進行預計，預計結(jié)果的準確性得到進一步提高。

圖4(e)、圖4(f)在形態(tài)上具有一定的相似性，地表沉陷的極值范圍都出現(xiàn)在走向方向上方，符合實際開采地表沉陷形態(tài)。運用Delaunary三角形分割算法，將預計工作面分割成49個預計單元，較人為劃分成7個預計單元相比，預計結(jié)果可靠性更高。將工作面上方地表監(jiān)測點的下沉值進行對比分析，最大絕對誤差為-65.8 mm，占相應點處下沉值的7.8%。對比可知，試驗3中工作面A的預計結(jié)果明顯優(yōu)于試驗2中工作面A的預計結(jié)果。實際工程應用中，對于大范圍、基于鉆孔數(shù)據(jù)的煤層開采，基于Delaunay三角形分割預計，較人為劃分預計在資源損耗等方面具有明顯的優(yōu)越性。

綜上分析可以認為，對于基于鉆孔數(shù)據(jù)的煤層開采，將工作面進行Delaunay分割后進行預計，可以有效提高預計結(jié)果的精度和可靠性。

5　工程實例

某礦目前為基建礦井，井下生產(chǎn)系統(tǒng)尚未形成。根據(jù)前期鉆探結(jié)果，初步掌握了礦區(qū)煤層厚度分布情況，如圖5(a)。現(xiàn)欲對地下煤礦開采引起的地表沉陷情況進行定量預計，根據(jù)預計結(jié)果進一步指導煤礦生態(tài)環(huán)境影響評價和分析工作。根據(jù)要求，利用前期鉆孔揭示的煤層厚度對該礦煤礦開采引起的地表沉陷進行預計，同時不考慮保護煤柱留設。因該礦井為新建礦井，無法根據(jù)實測數(shù)據(jù)進行參數(shù)反演，所以，參考臨近礦區(qū)開采沉陷預計參數(shù)。預計參數(shù)見表2，預計結(jié)果見圖5(c)。將預計結(jié)果同土地利用現(xiàn)狀圖、地形現(xiàn)狀圖和地面設施布置圖、地質(zhì)災害現(xiàn)狀圖等相關(guān)圖件進行疊加便可展開相關(guān)的損害評估工作。

(a)

(b)

(c)

6　結(jié)語

以礦區(qū)鉆孔數(shù)據(jù)為基礎(chǔ)，建立的基于Delaunay三角形剖分的的煤層開采沉陷預計模型運用到大范圍煤礦開采環(huán)評工作中，可以更加準確的掌握大面積煤礦開采可能造成的沉陷面積以及沉陷程度以及煤礦開采所引發(fā)的一系列問題，使煤礦生態(tài)環(huán)評從原來以定性分析、半定量預計為主的狀況，提高到以定量預測為主的高度，為提早介入煤礦沉陷區(qū)域的土地整治和土地復墾提供更加科學的依據(jù)，為礦區(qū)生態(tài)環(huán)境保護做出重要貢獻。

[1]鐔志偉.煤層群開采的地面沉陷評價研究[D].北京：中國地質(zhì)大學(北京),2007.

[2]武曉波,王世新,肖春生. Delaunay三角網(wǎng)的生成算法研究[J].測繪學報,1999,28(1):30-37.

[3]DELAUNAY B，SUR LA SPHERE VIDE. Bulletin of the Academy of Sciences of the USSR [J].Classe des Sciences Mathematiques et Naturelles, 1934(8):793-800.

[4]田錦州.開采沉陷智能化預計系統(tǒng)的研發(fā)[J].礦山測量,2012(1):72-74.

[5]邵春麗,胡鵬,黃承義,等.DELAUNAY三角網(wǎng)的算法詳述及其應用發(fā)展前景[J].測繪科學,2004,29(6):68-71.

[6]羅小華.Delaunay 三角剖分算法研究[D]. 廣東：暨南大學,2011.

[7]毛慧娟.關(guān)于辛卜生公式[J].教材通訊,1990(4):34-35.

[8]蔣和理.二重積分優(yōu)化復化梯形與辛卜生算法[J].合肥工業(yè)大學學報(自然科學版),1987(6):34-43.

[9]朱劉娟,陳俊杰,鄒友峰.任意形狀工作面開采地表移動變形預計的算法實現(xiàn)[J].遼寧工程技術(shù)大學學報,2005,24(3):337-340.

Surface Subsidence Prediction of Mining Based on Borehole Datas

XU Mengqiang1,2ZHA Jianfeng1,2WANG Jintao1,2

(1.SchoolofEnvironmentScienceandSpatialInformatics,ChinaUniversityofMiningandTechnologyXuzhou,Jiangsu221116)

The environmental impact assessment caused by mining subsidence is an important aspect of EIA, and the existing software is not fit for the work of the large area EIA. Aimed at the present problems, the mining subsidence prediction model based on the Delaunay triangulation is established. Through the contrast test, this new model is more consistent with the actual surface subsidence and the reliability is higher compared to the traditional one. By using this model in the practice, the human resources can be saved, making the coal mine ecological environmental impact assessment improved from the original qualitative analysis and semi quantitative analysis to the high degree of quantitative prediction, and better to protect the ecological environment.

environmental assessmentDelaunay triangletriangular partitionprobability integral methodmining subsidence predictionexperimental design

2015-12-07)

徐孟強，男，1989年生，安徽亳州人，碩士研究生，主要從事變形監(jiān)測與開采沉陷的研究。