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

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

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

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

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

1 果蠅算法

1.1 果蠅算法簡介

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

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

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

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

(1)

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

(2)

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

Smelli=Function(Si)

(3)

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

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

(4)

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

Smellbest=bestSmell

X-axis=X(bestIndex)

Y-axis=Y(bestIndex)

(5)

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

1.2 變步長果蠅算法

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

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

(6)

(7)

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

2 概率積分法優(yōu)化參數(shù)

2.1 概率積分法預計模型

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

(8)

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

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

(9)

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

對于矩形工作面，各參數(shù)量計算見式(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為傾向工作面的計算長度，其他參數(shù)同上。

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

(19)

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

2.2 反演流程

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

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

3 實例分析

3.1 試驗區(qū)概況

西山煤田某礦22101工作面，東西走向長440 m，南北寬120 m，地表呈西高東低，地面標高為+1 290～+1 305 m，無名溝通過工作面東部，溝內(nèi)出露上石盒子組地層，其余全部被黃土覆蓋，黃土層0～80 m，蓋山厚80～230 m，試驗區(qū)內(nèi)地勢較為開闊，地形傾向西南，傾角4°左右，西南坡底有一深切沖溝，邊坡甚為陡峭。開采煤層為二疊系下統(tǒng)山西組2.3#煤層，工作面內(nèi)的煤層厚度由西部開切眼附近4 m，向東逐漸減小至2.6 m，建筑區(qū)下方采煤厚度按綜合采煤機條件為3.5 m，實際平均采厚為3.18 m，工作面煤層基本為一單斜構造，在實驗區(qū)下方，煤層走向北東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。

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

表1 22101工作面開采技術參數(shù)

3.2 反演結果分析

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

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

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

表2 參數(shù)反演結果對比

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

表3 反演結果對比

點號實際下沉值/mm擬合下沉值/mm絕對誤差/mm相對誤差/%傳統(tǒng)算法改進果蠅算法傳統(tǒng)算法改進果蠅算法傳統(tǒng)算法改進果蠅算法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 結 語

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