概率積分法沉陷預計與參數反演優(yōu)化算法及實現

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

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

0 引言

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

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

1 基于用Delaunay三角剖分的沉陷預計方法

1.1 概率積分法簡介

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

下沉：

水平移動：

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

1.2 技術路線

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

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

(2)計算每一小三角單元的開采對地表點P的下沉影響理論值。

(3)判斷每一小三角單元距離煤壁的距離，根據距離和該小三角形開采區(qū)域的采深，將乘以一個影響系數ρ(ρ≤1)，獲得該小三角形區(qū)域開采對P點的實際下沉影響值。

(4)對各小三角單元開采影響進行疊加，求得開采引起的P點的沉降量。

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

1.3 Delaunay三角剖分

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

首先對開采區(qū)域進行插點，插點距離為預計范圍最短邊長的1/10，且不小于20m時能夠獲得比較好的效果。對開采區(qū)域插點完畢后，即可使用三角網生長算法來對開采區(qū)域進行三角剖分，三角網生長算法構網的實現過程如圖1所示。

圖1 Delaunay三角形剖分算法構網過程

該算法對凸多邊形和凹多邊形剖分結果如圖2所示，可以看出，該種算法剖分的多邊形區(qū)域單元三角格網比較規(guī)則，沒有出現狹長邊。

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

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

Delaunay三角剖分后，計算每一個小三角區(qū)域的開采對地表點P的影響，疊加即可得到開采區(qū)域對地表點P的下沉影響。由于單元三角形的范圍較小，可以將單元區(qū)域的煤層視為水平煤層來計算單元區(qū)域開采對P點的影響。

圖3 單元三角網格開采沉陷影響計算示意圖

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

采用先積分μ再積分λ進行計算，公式為：

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

使用無窮級數將冪級數展開至二次項，則有：

將冪指數積分轉換為二元二次函數積分，通過變量替換，可以相對簡單對上式進行二次積分。

1.5 單元開采區(qū)域影響系數ρ的計算

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

圖4 Maxwell-Boltzmann分布函數的積分函數圖像

三角單元影響系數ρ可使用該函數表達：

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

2 改進的概率積分法參數反演方法

2.1 常用概率積分法參數反演方法簡介

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

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

2.2 組合算法

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

圖5 組合算法流程圖

3 系統(tǒng)實現

3.1 系統(tǒng)簡介

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

圖6 軟件主界面

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

3.2.1 實測數據處理模塊

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

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

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

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

圖7 沉降動態(tài)求參結果圖

圖8 水平移動動態(tài)求參結果圖

3.2.3 變形預計模塊

基于Delaunay三角剖分算法開發(fā)了地表移動變形預計模塊，實現了任意多邊形工作面沉陷預測，輸出預計結果、各類等值線圖及三維沉陷模型，并可以通過旋轉從多角度查看三維模型(圖9、圖10)。

圖9 等沉陷曲線圖

圖10 沉陷三維模型圖

4 結論

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

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

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

[1] 趙博,李愛軍,劉術明．昭陽煤礦礦山地質環(huán)境影響評估及監(jiān)測建議[J].山東國土資源,2017,33(5):35-41.

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

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

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

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

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

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

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

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

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

[11] 潘國榮，谷川.改進的遺傳算法用于工業(yè)測量數據處理[J].大地測量與地球動力學，2008，28(1)：55-58.

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

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

[14] 陳永良，周斌，李學斌.基于Boltzmann機的礦產靶區(qū)預測[J].地球物理學進展,2012,27(1):179-185．

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

[16] 劉寶琛,戴華陽.概率積分法的由來與研究進展[J].煤礦開采,2016(2):1-3.

[17] 郭文兵,鄧喀中,鄒友峰.概率積分法預計參數選取的神經網絡模型[J].中國礦業(yè)大學學報,2004(3):88-92.

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

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

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