基于改進果蠅算法的概率積分法預計參數反演

周 棒，康建榮，胡晉山，孫富強

(江蘇師范大學地理測繪與城鄉規劃學院，江蘇 徐州 221116)

開采沉陷預計是礦山開采沉陷學的核心內容之一，它對于開采沉陷的理論研究以及對礦區的綠色發展、煤礦的安全開采和塌陷區的土地復墾等具有非常重要的意義。基于隨機介質理論的概率積分法，是我國目前較成熟應用最為廣泛的預計方法之一，該方法理論成熟，便于計算機實現，被廣泛應用于礦區地表沉陷預計[1]。

概率積分法本身的模型誤差以及各參數的取值誤差是開采沉陷預測過程中主要的誤差存在[2]，由于概率積分法著重于對開采沉陷現象幾何表達方式缺少力學機理的解釋[3]，使得模型誤差的缺陷難以修正，因此減小預計參數的取值誤差是目前提高開采沉陷預計精度的主要手段。參數反演方法經歷了用線性最小二乘近似法直接反演[4]、用模矢法實驗設計反演[5]、用遺傳算法[6-8]、粒子群算法[9-10]、蟻群算法[11]、神經網絡算法[12-14]、果蠅算法[15-16]等智能算法反演的過程。而智能算法普遍存在容易陷入局部最優值、過早收斂等問題，因此學者們對這些算法加以改進及優化。本文將參數反演問題轉化為多維參數的組合優化問題，用改進的果蠅算法對概率積分法預計參數進行反演。

1 果蠅算法

1.1 果蠅算法簡介

果蠅優化算法FOA(fruit fly optimization algorithm)是由臺灣華夏科技大學潘文超教授于2011年提出，果蠅優化算法是基于果蠅群體在覓食過程中的合作與競爭行為，提出的一種新群智能優化算法[17]。

果蠅群體在尋找食物的過程中，憑借其敏銳的嗅覺和視覺感官能夠快速定位并找到食物，根據果蠅通過感官覓食的過程得到果蠅優化算法的基本原理[18]：首先是嗅覺搜索粗略定位階段，利用嗅覺充分感知空氣中彌漫的各種氣味信息并區分出食物氣味，根據食物氣味判斷食物的大約位置，并向目標飛向；然后是視覺查找精確定位階段，在飛到食物周邊且能夠進行視覺查找的范圍內，準確判斷食物的具體位置，最終飛向食物。算法主要流程如下所述。

1) 初始化果蠅算法相關參數。果蠅種群位置(X-axis,Y-axis)，迭代次數maxgen,種群規模sizepop。

2) 給定隨機方向及步長k，見式(1)。

(1)

3) 通過計算果蠅個體距離初始位置的距離Disti與味道濃度的判定值Si來確定食物的具體方位，見式(2)。

(2)

4) 將氣味濃度值代入適應度函數，計算當前每個果蠅個體的氣味濃度值Smelli，見式(3)。

Smelli=Function(Si)

(3)

5) 比較當前種群中果蠅個體的氣味濃度值并標記氣味濃度最高(或者最低)的果蠅個體即最佳果蠅個體，見式(4)。

[bestSmell,bestIndex]=min/max(Smelli)

(4)

6) 記錄并保留最佳味道濃度值bestSmell以及該最優果蠅個體坐標位置(X,Y)，利用敏銳的視覺飛往最優個體位置。

Smellbest=bestSmell

X-axis=X(bestIndex)

Y-axis=Y(bestIndex)

(5)

7) 迭代尋優。判斷當前迭代次數是否等于最大迭代次數，若當前迭代次數等于最大迭代次數，則結束程序，并輸出果蠅優化算法的結果，否則步驟2開始繼續執行。

1.2 變步長果蠅算法

為了增加果蠅優化算法的普適性，學者們提出了許多改進的方法使得該算法在計算量、復雜度、穩定性以及精度方面都有所改善[19-20]。目前對于果蠅優化算法的改進，主要從兩方面進行研究：一是通過對原始果蠅算法自身的參數或者推演公式進行完善和改進；二是將果蠅算法與其他智能算法進行結合。

果蠅優化算法設計的參數主要有果蠅初始位置，果蠅規模以及搜索步長，本文主要對果蠅的搜索步長進行改進，用變步長果蠅算法對概率積分法預計參數進行反演。當果蠅群體初始位置與規模確定時，果蠅個體在搜索食物的過程中，搜索步長越大，相應的搜索范圍越大，提高了全局尋優性能，此時，局部搜索的能力也會隨之下降；反之，若搜索步長越小，則尋優范圍同步減小，此時果蠅個體能夠進行更為精細準確的搜索，改善了局部搜索能力，相應地全局尋優性能變弱甚至由于步長過小導致陷入局部最優的狀況。而固定的步長也會使果蠅算法無法快速收斂到最優解。針對果蠅算法以上的不足，在果蠅優化算法的基礎上引入了自適應變步長，每次迭代搜索步長的取值與當代和上一代的最優味道濃度值的比值有關。對果蠅優化算法改進的搜索步長如下(以最優味道濃度最小為例)。

(6)

(7)

在文獻[21]中對自適應變步長果蠅算法進行了仿真實驗，證明了此算法在收斂速度和精度上的優越性。

2 概率積分法優化參數

2.1 概率積分法預計模型

假設煤層傾斜，根據概率積分法的基本原理，開采單元B(s,t)進行開采引起的地表點A(x,y)的下沉值計算公式見式(8)。

(8)

式中：主要影響半徑r=H/tanβ;H為單元B處的采深,m;tanβ為主要影響角正切；θ0為開采影響傳播角,°。

根據疊加原理，整個開采引起的A點的下沉量計算式見式(9)。

(9)

式中：w0=mcqcosα為充分采動時地表最大下沉值;mc為煤層法向厚度,m;q為下沉系數;α為煤層平均傾角,°；S為考慮拐點偏移距s后的開采區域。

對于矩形工作面，各參數量計算見式(10)～(18)[22]。

(10)

w0(x)=w(x)-w(x-l)

(11)

(12)

(13)

w0(y)=w(y)-w(y-L)

(14)

(15)

(16)

l=D3-s3-s4

(17)

(18)

式中：D1為工作面傾向斜長；D3為工作面走向長度；s1、s2為下山和上山方向的拐點偏移距；s3、s4為走向左、右邊界的拐點偏移距；l為走向工作面的計算長度；L為傾向工作面的計算長度，其他參數同上。

基于下沉值的概率積分法預計參數有：下沉系數q，主要影響角正切tanβ，開采影響傳播角θ0，矩形工作面拐點偏移距si(i=1,2,3,4)。設p=[qtanβθ0s1s2s3s4]，P為p的搜索空間(預計參數的取值范圍)，以改進步長的果蠅算法尋優得到的概率積分法預計參數進行開采沉陷預計的預計值與對應點的實測值的均方差最小為原則，則根據適應度函數fp將概率積分法參數反演的過程表示為向量p的約束優化問題，即在給定空間P中找到一個向量p0，使適應度函數值最小，見式(19)。

(19)

式中：fp為目標函數；Wi為實測下沉值；wi為預計下沉值；m為觀測點數。

2.2 反演流程

根據該礦區的地質采礦資料并參照《建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采指南》內容，設置參數向量p的尋優范圍及步長，由果蠅算法的原理結合概率積分法的基本理論可得基于改進果蠅算法的概率積分法預計參數反演的步驟：①設置初始化種群數、最大循環次數等果蠅算法參數，根據預計參數向量p的范圍生成初始種群，種群初始化的位置即N維向量的初始解，N為待優化參數的個數(此時N=7)；②計算適應度函數值，結合工作面信息進行概率積分法的開采沉陷預計并更新步長；③在每次迭代過程中，判斷并選擇最優適應度值的種群，如此循環直到找到最優食物源，即解向量，其具體實現流程見圖1。

圖1 改進步長果蠅算法反演概率積分法預計參數流程

3 實例分析

3.1 試驗區概況

西山煤田某礦22101工作面，東西走向長440 m，南北寬120 m，地表呈西高東低，地面標高為+1 290～+1 305 m，無名溝通過工作面東部，溝內出露上石盒子組地層，其余全部被黃土覆蓋，黃土層0～80 m，蓋山厚80～230 m，試驗區內地勢較為開闊，地形傾向西南，傾角4°左右，西南坡底有一深切沖溝，邊坡甚為陡峭。開采煤層為二疊系下統山西組2.3#煤層，工作面內的煤層厚度由西部開切眼附近4 m，向東逐漸減小至2.6 m，建筑區下方采煤厚度按綜合采煤機條件為3.5 m，實際平均采厚為3.18 m，工作面煤層基本為一單斜構造，在實驗區下方，煤層走向北東10°，傾向南西，傾角1～ 3°，煤層底板標高為1 070～1 085 m，平均開采深度約220 m。在工作面上方設置了走向和傾向兩條主斷面地表移動觀測線，觀測線總長為1 060 m，測點45個，共進行了14次觀測，獲得了預期的實測資料。其中走向觀測線上513.5 m，設置測點26個，走向26個點的點號為1～26;傾向觀測線387.8 m，設置測點19個,傾向19個點的點號為27～45。

工作面的地質條件為西部簡單，中部復雜，東部較復雜。本試驗區下方為工作面西部，地質條件屬簡單類型,其地表地形以及點位布設見圖2。根據22101工作面觀測站設計報告可知該工作面開采技術參數，見表1。

表1 22101工作面開采技術參數

3.2 反演結果分析

根據反演流程，22101工作面開采沉陷預計參數的取值下沉系數q為0.3～1，主要影響角正切tanβ為1～3，開采影響傳播角θ0為80～90°，拐點偏移距s為-30～30 m。對變步長果蠅優化算法初始化參數取值分別為種群規模sizepop5、10、20，迭代次數 50、100、500。下沉系數q的步長取值0.1、0.2，主要影響角正切tanβ的步長取值0.2、0.3，開采影響傳播角θ0步長取值5、8，拐點偏移距s為10、20、30。經過調試分析，確定22101工作面開采沉陷預計參數的取值見表2。表2中傳統方法反演的預計參數是根據22101工作面地表移動觀測站成果分析綜合表中分析經驗公式得到的。下沉系數q=0.79±0.1，主要影響角正切tanβ=2.4，開采影響傳播角θ0=90°-kα(k為系數，α為煤層傾角)，拐點偏移距s1、s2為0.05H，s3、s4為0.1H，H為平均采深。

圖2 某礦22101工作面地表移動觀測站布設

利用改進的變步長果蠅算法擬合下沉值與實測下沉值的對比見圖3。考慮到文章篇幅，選取了改進果蠅算法較傳統算法擬合精度有明顯提高的部分點羅列于表3中。從圖3和表3中可以看到大部分觀測點的擬合值與實測值誤差較小，表明變步長果蠅算法在概率積分法預計參數反演中的可行性和有效性。但從圖3中可以看到在測點17～20處擬合下沉值較實測下沉值偏小，而測點23～25處擬合下沉值較實測下沉值偏大，且兩處擬合值與實測值結果有較大差異，其原因是17～20號測點位于最大下沉區域附近，23～25號測點又位于地形條件復雜的工作面中部，這兩處區域可能受對應地表本身復雜的地質條件以及開采影響下的斷層、裂縫、坡體滑移等非連續變形的影響使得擬合下沉值較實測下沉值出現較大偏差，從而影響用概率積分法進行開采沉陷預計的精度。

表2 參數反演結果對比

圖3 擬合下沉值與實測下沉值對比

表3 反演結果對比

點號實際下沉值/mm擬合下沉值/mm絕對誤差/mm相對誤差/%傳統算法改進果蠅算法傳統算法改進果蠅算法傳統算法改進果蠅算法121 5751 751.961 477.74176.96-97.2711.2-6.1131 6681 850.041 571.78182.04-96.2210.9-5.7141 7411 897.561 623.33156.56-117.679.06.7151 7931 917.221 648.64124.22-144.366.9-8.0161 7391 922.221 657.04183.22-81.9610.5-4.7411 4041 976.971 511.88572.97107.8840.87.7421 6961 906.791 658.087210.79-37.9112.4-2.2431 5221 272.411 485.011-249.59-36.99-16.4-2.4

4 結 語

針對傳統方法在求取概率積分法預計參數時存在的易發散且容易陷入局部最優，難以得到最優解等問題，提出了一種新的全局尋優的智能算法，即改進步長的果蠅優化算法，該算法在算法設計、參數設置、尋優精度上具有一定的優勢，以某礦22101工作面為例，實驗結果表明，采用改進步長的果蠅算法反演概率積分法參數效率可觀，預計精度較高，對于大幅度提高概率積分法開采沉陷預計參數的反演精度有一定的借鑒價值。