薄懷志，宋炳忠，李川建，王猛，盧國(guó)宏，陳宗成，李云偉

(1.山東省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局第三水文地質(zhì)工程地質(zhì)大隊(duì)，山東 兗州 272100；2.山東省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局第二地質(zhì)大隊(duì)，山東 兗州 272100；3.山東省地礦局采煤塌陷地綜合治理工程技術(shù)研究中心，山東 兗州 272100；4.武漢大學(xué) 測(cè)繪學(xué)院， 湖北 武漢 430000)

0 引言

目前概率積分法沉陷預(yù)計(jì)在“三下采煤”、礦區(qū)土地復(fù)墾、采煤塌陷地治理規(guī)劃設(shè)計(jì)、礦區(qū)建筑物穩(wěn)定性評(píng)估等工作中應(yīng)用廣泛[1-4]。近年來(lái)，對(duì)任意多邊形的沉陷預(yù)計(jì)研究較多，吳侃[5]首先將非矩形工作面劃分為若干小矩形工作面，分別進(jìn)行預(yù)計(jì)，然后將各個(gè)小矩形預(yù)計(jì)結(jié)果疊加的方法；許冬等[6]將任意多邊形工作面按照角點(diǎn)劃分為若干梯形，分別預(yù)計(jì)后疊加(減)處理實(shí)現(xiàn)沉陷預(yù)計(jì)；徐孟強(qiáng)等[7]提出了利用三角剖分根據(jù)角點(diǎn)將多邊形工作面劃分為若干三角形、分別預(yù)計(jì)后對(duì)各三角形預(yù)計(jì)結(jié)果疊加(減)的方法?？梢钥闯?，對(duì)于任意多邊形開(kāi)采區(qū)域劃分多從角點(diǎn)往內(nèi)進(jìn)行剖分，容易造成小的剖分區(qū)域重復(fù)計(jì)算或遺漏，難以真正實(shí)現(xiàn)任意多變形開(kāi)采區(qū)域沉陷預(yù)計(jì)。概率積分法預(yù)計(jì)參數(shù)反演常用模矢法和遺傳算法。模矢法是傳統(tǒng)的直接尋優(yōu)算法的典型代表，但該算法對(duì)初值可靠性要求較高，初值選取不當(dāng)，求參結(jié)果易陷入局部極值[8-9]。遺傳算法是典型的智能算法，但算法求解有一定的不穩(wěn)定性，即使求參數(shù)據(jù)完全相同，每次求參結(jié)果都會(huì)有一定的差異，且算法執(zhí)行效率低、后期收斂慢[10-11]。

該文提出了使用Delaunay三角剖分方法，將任意多邊形開(kāi)采區(qū)域劃分為若干三角形，然后對(duì)三角形區(qū)域進(jìn)行沉陷預(yù)計(jì)，累加得到總沉陷值的方法，解決了任意多邊形沉陷預(yù)計(jì)的難題。綜合模矢法和遺傳算法優(yōu)點(diǎn)，提出了組合算法，提高了概率積分法參數(shù)反演精度和效率，開(kāi)發(fā)了集實(shí)測(cè)數(shù)據(jù)處理、預(yù)計(jì)模型參數(shù)反演和移動(dòng)變形預(yù)計(jì)三大模塊于一體的地表移動(dòng)變形數(shù)據(jù)處理軟件。

1 基于用Delaunay三角剖分的沉陷預(yù)計(jì)方法

1.1 概率積分法簡(jiǎn)介

概率積分法是以隨機(jī)介質(zhì)理論為基礎(chǔ)，以正態(tài)分布為影響函數(shù)，用積分式描述下沉盆地的方法，因其預(yù)計(jì)公式中包括概率積分計(jì)算而得名[12]。概率積分法預(yù)計(jì)計(jì)算公式如下：

下沉：

水平移動(dòng)：

式中：W(x,y)為地表任意點(diǎn)的下沉值(mm)；Ux(x,y)為走向方向的水平移動(dòng)值(mm)；Uy(x,y)為傾斜方向的水平移動(dòng)值(mm)；Wem，Uem為充分采動(dòng)時(shí)地表最大下沉值和最大水平移動(dòng)值(mm)；r為走向主斷面上采空區(qū)邊界(m)；μ,λ分別為x,y方向的變形值(mm)；x,y為計(jì)算點(diǎn)相對(duì)坐標(biāo)；D為開(kāi)采煤層區(qū)域。

1.2 技術(shù)路線

基于用Delaunay三角剖分的概率積分法，沉陷預(yù)計(jì)算法流程如下：

(1)對(duì)于任意多邊形工作面進(jìn)行Delaunay三角形剖分。

(2)計(jì)算每一小三角單元的開(kāi)采對(duì)地表點(diǎn)P的下沉影響理論值。

(3)判斷每一小三角單元距離煤壁的距離，根據(jù)距離和該小三角形開(kāi)采區(qū)域的采深，將乘以一個(gè)影響系數(shù)ρ(ρ≤1)，獲得該小三角形區(qū)域開(kāi)采對(duì)P點(diǎn)的實(shí)際下沉影響值。

(4)對(duì)各小三角單元開(kāi)采影響進(jìn)行疊加，求得開(kāi)采引起的P點(diǎn)的沉降量。

(5)匯總計(jì)算得到預(yù)計(jì)區(qū)域的沉陷值。

1.3 Delaunay三角剖分

Delaunay三角剖分具有以下特性：①最接近：以最近的三點(diǎn)形成三角形，各邊皆不相交；②唯一性：從任何位置開(kāi)始構(gòu)建，剖分結(jié)果均相同；③最優(yōu)性：任意兩個(gè)相鄰三角形組成的凸四邊形對(duì)角線互換，內(nèi)角中最小角不會(huì)變大；④區(qū)域性：新增、刪除、移動(dòng)某一個(gè)頂點(diǎn)時(shí)只會(huì)影響臨近的三角形[13]。Delaunay三角剖分可以最大限度地避免狹長(zhǎng)三角形的產(chǎn)生，采用Delaunay三角剖分，對(duì)任意多邊形開(kāi)采區(qū)域進(jìn)行剖分非常合適。

首先對(duì)開(kāi)采區(qū)域進(jìn)行插點(diǎn)，插點(diǎn)距離為預(yù)計(jì)范圍最短邊長(zhǎng)的1/10，且不小于20m時(shí)能夠獲得比較好的效果。對(duì)開(kāi)采區(qū)域插點(diǎn)完畢后，即可使用三角網(wǎng)生長(zhǎng)算法來(lái)對(duì)開(kāi)采區(qū)域進(jìn)行三角剖分，三角網(wǎng)生長(zhǎng)算法構(gòu)網(wǎng)的實(shí)現(xiàn)過(guò)程如圖1所示。

圖1 Delaunay三角形剖分算法構(gòu)網(wǎng)過(guò)程

該算法對(duì)凸多邊形和凹多邊形剖分結(jié)果如圖2所示，可以看出，該種算法剖分的多邊形區(qū)域單元三角格網(wǎng)比較規(guī)則，沒(méi)有出現(xiàn)狹長(zhǎng)邊。

圖2 任意多邊形的Delaunay三角剖分

1.4 單元三角區(qū)域引起地面沉陷計(jì)算

Delaunay三角剖分后，計(jì)算每一個(gè)小三角區(qū)域的開(kāi)采對(duì)地表點(diǎn)P的影響，疊加即可得到開(kāi)采區(qū)域?qū)Φ乇睃c(diǎn)P的下沉影響。由于單元三角形的范圍較小，可以將單元區(qū)域的煤層視為水平煤層來(lái)計(jì)算單元區(qū)域開(kāi)采對(duì)P點(diǎn)的影響。

圖3 單元三角網(wǎng)格開(kāi)采沉陷影響計(jì)算示意圖

如圖3所示，以△OAB的最長(zhǎng)邊為x軸，建立如圖所示的局部坐標(biāo)系，設(shè)O點(diǎn)位原點(diǎn)，A點(diǎn)坐標(biāo)為(x2,0)，B點(diǎn)坐標(biāo)為(x1，y1)。則OB，AB的直線方程可表示為：

采用先積分μ再積分λ進(jìn)行計(jì)算，公式為：

由于單元三角區(qū)域較小，開(kāi)采影響半徑r可以采用單元三角中心點(diǎn)的采深進(jìn)行計(jì)算。令：

使用無(wú)窮級(jí)數(shù)將冪級(jí)數(shù)展開(kāi)至二次項(xiàng)，則有：

將冪指數(shù)積分轉(zhuǎn)換為二元二次函數(shù)積分，通過(guò)變量替換，可以相對(duì)簡(jiǎn)單對(duì)上式進(jìn)行二次積分。

1.5 單元開(kāi)采區(qū)域影響系數(shù)ρ的計(jì)算

改進(jìn)的算法可以實(shí)現(xiàn)任意多邊形開(kāi)采的工作面沉陷預(yù)計(jì)，由于頂板的懸臂作用，距離開(kāi)采邊界越近的單元，區(qū)域頂板下沉量越小，當(dāng)超過(guò)開(kāi)采邊界一定距離后，頂板才能到達(dá)最大下沉量。需要引入一個(gè)修正系數(shù)，對(duì)單元三角區(qū)域下沉值進(jìn)行修正。三角單元開(kāi)采引起地表點(diǎn)實(shí)際沉降量與該單元距離采空區(qū)邊界的距離成正比，不同采深的三角單元開(kāi)采對(duì)地面的影響也不同，這種性質(zhì)與Maxwell-Boltzmann分布的積分函數(shù)變化形態(tài)非常相似(圖4)[14]。

圖4 Maxwell-Boltzmann分布函數(shù)的積分函數(shù)圖像

三角單元影響系數(shù)ρ可使用該函數(shù)表達(dá)：

式中：x為三角單元至采空區(qū)的距離；a為與采礦地質(zhì)條件有關(guān)的系數(shù)。

2 改進(jìn)的概率積分法參數(shù)反演方法

2.1 常用概率積分法參數(shù)反演方法簡(jiǎn)介

目前開(kāi)采沉陷預(yù)計(jì)時(shí)多使用概率積分法，初期概率積分法參數(shù)反演多使用線性近似法，后逐步發(fā)展到實(shí)驗(yàn)設(shè)計(jì)法和智能算法[15]。模矢法是直接尋優(yōu)算法的典型代表，應(yīng)用廣泛，但該算法對(duì)初值可靠性要求較高，初值選取不當(dāng)求參結(jié)果易陷入局部極值[16]。近年來(lái)，概率積分法參數(shù)反演智能算法發(fā)展迅速，典型代表有遺傳算法[10]、神經(jīng)網(wǎng)絡(luò)算法[17]、退火算法[18]、果蠅算法[19]、隨機(jī)粒子群算法[20]等，其中遺傳算法是智能算法的典型代表。

模矢法又稱(chēng)為步長(zhǎng)加速法，是求取非線性無(wú)約束最優(yōu)解問(wèn)題的算法。模矢法將待求參數(shù)看作一組矢量，通過(guò)改變尋優(yōu)方向和步長(zhǎng)尋找一組最優(yōu)參數(shù)，具有方便計(jì)算機(jī)編程和加速向最優(yōu)點(diǎn)移動(dòng)的優(yōu)點(diǎn)，且可根據(jù)任意形狀工作面的實(shí)測(cè)地表移動(dòng)觀測(cè)數(shù)據(jù)進(jìn)行概率積分法參數(shù)求取。遺傳算法是一種通過(guò)模擬自然進(jìn)化過(guò)程來(lái)搜索最優(yōu)解的方法，通過(guò)編碼將求解參數(shù)解映射為種群個(gè)體的染色體，使用交叉、變異算子維持種群的多樣性，使用選擇算子使種群向最優(yōu)方向進(jìn)化。遺傳算法具有隱含的并行搜索能力，能夠使求解參數(shù)跨過(guò)局部極值的障礙向全局最優(yōu)值進(jìn)化。

2.2 組合算法

首先利用遺傳算法的全局尋優(yōu)能力求得全局近似最優(yōu)參數(shù)，然后將近似最優(yōu)參數(shù)適當(dāng)取舍后作為模矢法探索的起點(diǎn)，最后利用模矢法的局部探索能力來(lái)獲得預(yù)計(jì)參數(shù)的精確解，具體步驟為：①使用遺傳算法反演得到預(yù)計(jì)參數(shù)的全局近似最優(yōu)值；②取近似最優(yōu)值小數(shù)點(diǎn)后一位作為初始模矢點(diǎn)的參數(shù)值；③使用模矢法反演出精確的最優(yōu)參數(shù)，具體流程如圖5所示。

圖5 組合算法流程圖

3 系統(tǒng)實(shí)現(xiàn)

3.1 系統(tǒng)簡(jiǎn)介

按照改進(jìn)的概率積分法算法，開(kāi)發(fā)了地表移動(dòng)變形數(shù)據(jù)處理軟件。使用C#語(yǔ)言開(kāi)發(fā)的基于Autocad.net的插件，可以加載在Autocad2006及以上的CAD版本中(圖6)。軟件縱向采用了分層設(shè)計(jì)的思想，開(kāi)采沉陷數(shù)據(jù)處理與預(yù)計(jì)軟件縱向結(jié)構(gòu)分為用戶交互層和數(shù)據(jù)處理層。橫向采用了模塊化設(shè)計(jì)思想，主要包括實(shí)測(cè)數(shù)據(jù)處理模塊、移動(dòng)變形預(yù)計(jì)模塊和最優(yōu)參數(shù)反演模塊。

圖6 軟件主界面

3.2 系統(tǒng)主要功能模塊

3.2.1 實(shí)測(cè)數(shù)據(jù)處理模塊

將Excel格式的沉降、水平移動(dòng)觀測(cè)數(shù)據(jù)導(dǎo)入軟件后，使用數(shù)據(jù)處理功能處理即可以直接繪制地表沉降過(guò)程中的各種動(dòng)態(tài)變形曲線圖。

可以根據(jù)實(shí)測(cè)數(shù)據(jù)求取超前影響距(角)與滯后距(角)等動(dòng)態(tài)地表移動(dòng)變形參數(shù)和邊界角、移動(dòng)角、地表移動(dòng)變形持續(xù)時(shí)間、移動(dòng)變形最大值等穩(wěn)態(tài)地表移動(dòng)變形參數(shù)。

3.2.2 最優(yōu)參數(shù)反演模塊

軟件提供模矢法和遺傳算法兩種概率積分法最優(yōu)參數(shù)反演方法，既可使用單一方法求參，也可以使用組合方法。如果沒(méi)有可供參考的相鄰采區(qū)或類(lèi)似地質(zhì)采礦條件下的預(yù)計(jì)參數(shù)作為反演初值，不宜采用模矢法進(jìn)行求參。而遺傳算法、組合算法則基本不受這種因素的影響，但是從算法的準(zhǔn)確性和效率方面來(lái)說(shuō)，組合算法又優(yōu)于遺傳算法，該算法比較適用于沒(méi)有相關(guān)預(yù)計(jì)參數(shù)可供參考情況下的概率積分法預(yù)計(jì)參數(shù)反演。軟件可以以觀測(cè)數(shù)據(jù)的時(shí)間為依據(jù)，進(jìn)行動(dòng)態(tài)或靜態(tài)求參(圖7、圖8)。

圖7 沉降動(dòng)態(tài)求參結(jié)果圖

圖8 水平移動(dòng)動(dòng)態(tài)求參結(jié)果圖

3.2.3 變形預(yù)計(jì)模塊

基于Delaunay三角剖分算法開(kāi)發(fā)了地表移動(dòng)變形預(yù)計(jì)模塊，實(shí)現(xiàn)了任意多邊形工作面沉陷預(yù)測(cè)，輸出預(yù)計(jì)結(jié)果、各類(lèi)等值線圖及三維沉陷模型，并可以通過(guò)旋轉(zhuǎn)從多角度查看三維模型(圖9、圖10)。

圖9 等沉陷曲線圖

圖10 沉陷三維模型圖

4 結(jié)論

(1)開(kāi)采沉陷數(shù)據(jù)處理與預(yù)計(jì)軟件開(kāi)采區(qū)域使用Delaunay三角剖分，由于Delaunay三角剖分具有最接近性、唯一性、最優(yōu)性、最規(guī)則、區(qū)域性的特點(diǎn)，適合用于不規(guī)則任意多邊形采煤沉陷預(yù)測(cè)，真正實(shí)現(xiàn)了任意形狀工作面的地表移動(dòng)變形預(yù)計(jì)。傳統(tǒng)的概率積分法根據(jù)拐點(diǎn)偏移距確定計(jì)算范圍，認(rèn)為拐點(diǎn)偏移距至煤壁的距離上覆煤層都是假設(shè)不發(fā)生下沉的。該算法使用Maxwell-Boltzmann函數(shù)，每個(gè)剖分三角單元的沉陷影響系數(shù)ρ(ρ≤1)，來(lái)解決傳統(tǒng)概率積分法拐點(diǎn)偏移距引起的誤差。

(2)將遺傳算法和模矢法用于概率積分法預(yù)計(jì)模型的最優(yōu)參數(shù)反演，優(yōu)劣互補(bǔ)，充分發(fā)揮兩者的優(yōu)點(diǎn)。如果礦區(qū)有比較可靠的參數(shù)參考值，直接以該參考值作為初值使用模矢法進(jìn)行解算，節(jié)省解算時(shí)間。如果礦區(qū)沒(méi)有可靠參數(shù)值，可以給參數(shù)值設(shè)定一個(gè)范圍，使用遺傳算法解算概率積分法最優(yōu)參數(shù)值。也可以將模矢法、遺傳算法進(jìn)行融合，首先利用遺傳算法的全局尋優(yōu)能力求得全局近似最優(yōu)參數(shù)，然后將近似最優(yōu)參數(shù)適當(dāng)取舍后作為模矢法探索的起點(diǎn)，最后利用模矢法的局部探索能力來(lái)獲得預(yù)計(jì)參數(shù)的精確解。

(3)開(kāi)發(fā)了集實(shí)測(cè)數(shù)據(jù)處理、最優(yōu)參數(shù)反演、采煤沉陷預(yù)計(jì)三大模塊于一體的開(kāi)采沉陷數(shù)據(jù)處理與預(yù)計(jì)軟件，即可處理礦山巖移觀測(cè)數(shù)據(jù)，繪制各類(lèi)變形曲線圖，也可實(shí)現(xiàn)概率積分法最優(yōu)參數(shù)反演和采煤沉陷預(yù)測(cè)。

[1] 趙博,李?lèi)?ài)軍,劉術(shù)明．昭陽(yáng)煤礦礦山地質(zhì)環(huán)境影響評(píng)估及監(jiān)測(cè)建議[J].山東國(guó)土資源,2017,33(5):35-41.

[2] 隋建紅.山東省黃河北煤炭礦區(qū)地表沉陷預(yù)測(cè)與綜合治理研究[J].山東國(guó)土資源,2017,33(3):59-63.

[3] 孔凡杜.滕縣煤礦區(qū)礦山地質(zhì)環(huán)境綜合評(píng)價(jià)研究[J].山東國(guó)土資源,2017,33(1):59-64.

[4] 隋建紅,郝啟勇.山東省黃河北煤炭礦區(qū)地質(zhì)環(huán)境承載力評(píng)價(jià)[J].山東國(guó)土資源,2017,33(4):34-40.

[5] 吳侃，葛家新，王玲丁，等．開(kāi)采沉陷預(yù)計(jì)一體化方法[M]．徐州:中國(guó)礦業(yè)大學(xué)出版社，1998：1-20.

[6] 許冬，王臨清，吳侃．任意形狀工作面沉陷預(yù)計(jì)計(jì)算方法[J]．金屬礦山，2014(5):55-59．

[7] 徐孟強(qiáng),查劍鋒.基于三角形剖分算法的開(kāi)采沉陷預(yù)計(jì)系統(tǒng)設(shè)計(jì)[J].金屬礦山,2016(2):128-131.

[8] 朱曉峻,郭廣禮,方齊.概率積分法預(yù)計(jì)參數(shù)反演方法研究進(jìn)展[J].金屬礦山,2015(4):173-177.

[9] 冉典,吳捷.模矢法在厚松散層地表移動(dòng)參數(shù)解算中的應(yīng)用[J].地礦測(cè)繪,2014(2):25-27.

[10] 查劍鋒,馮文凱,朱曉峻.基于遺傳算法的概率積分法預(yù)計(jì)參數(shù)反演[J].采礦與安全工程學(xué)報(bào),2011,28(4):655-659.

[11] 潘國(guó)榮，谷川.改進(jìn)的遺傳算法用于工業(yè)測(cè)量數(shù)據(jù)處理[J].大地測(cè)量與地球動(dòng)力學(xué)，2008，28(1)：55-58.

[12] 國(guó)家煤炭工業(yè)局．建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開(kāi)采規(guī)程[S].北京:煤炭工業(yè)出版社，2004.

[13] 劉建耀,劉保順.基于C#和AutoCAD的Delaunay三角剖分算法的實(shí)現(xiàn)[J].有色礦冶,2014,30(5):9-14.

[14] 陳永良，周斌，李學(xué)斌.基于Boltzmann機(jī)的礦產(chǎn)靶區(qū)預(yù)測(cè)[J].地球物理學(xué)進(jìn)展,2012,27(1):179-185．

[15] 朱曉峻,郭廣禮,方齊.概率積分法預(yù)計(jì)參數(shù)反演方法研究進(jìn)展[J].金屬礦山,2015(4):173-177.

[16] 劉寶琛,戴華陽(yáng).概率積分法的由來(lái)與研究進(jìn)展[J].煤礦開(kāi)采,2016(2):1-3.

[17] 郭文兵,鄧喀中,鄒友峰.概率積分法預(yù)計(jì)參數(shù)選取的神經(jīng)網(wǎng)絡(luò)模型[J].中國(guó)礦業(yè)大學(xué)學(xué)報(bào),2004(3):88-92.

[18] 蘇軍明,朱建軍,伍雅晴,等.模擬退火算法在概率積分法參數(shù)反演中的應(yīng)用[J].工程勘察,2015(12):76-79.

[19] 陳濤,郭廣禮,朱曉峻,等.利用果蠅算法反演概率積分法開(kāi)采沉陷預(yù)計(jì)參數(shù)[J].金屬礦山,2016(6):185-188.

[20] 王正帥,鄧喀中,康建榮.概率積分法參數(shù)反演的文化-隨機(jī)粒子群優(yōu)化算法[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào)(自然科學(xué)版),2013(3):311-315.