周 棒，康建榮，胡晉山，孫富強(qiáng)

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

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

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

1 果蠅算法

1.1 果蠅算法簡(jiǎn)介

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

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

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

2) 給定隨機(jī)方向及步長(zhǎng)k，見(jiàn)式(1)。

(1)

3) 通過(guò)計(jì)算果蠅個(gè)體距離初始位置的距離Disti與味道濃度的判定值Si來(lái)確定食物的具體方位，見(jiàn)式(2)。

(2)

4) 將氣味濃度值代入適應(yīng)度函數(shù)，計(jì)算當(dāng)前每個(gè)果蠅個(gè)體的氣味濃度值Smelli，見(jiàn)式(3)。

Smelli=Function(Si)

(3)

5) 比較當(dāng)前種群中果蠅個(gè)體的氣味濃度值并標(biāo)記氣味濃度最高(或者最低)的果蠅個(gè)體即最佳果蠅個(gè)體，見(jiàn)式(4)。

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

(4)

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

Smellbest=bestSmell

X-axis=X(bestIndex)

Y-axis=Y(bestIndex)

(5)

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

1.2 變步長(zhǎng)果蠅算法

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

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

(6)

(7)

在文獻(xiàn)[21]中對(duì)自適應(yīng)變步長(zhǎng)果蠅算法進(jìn)行了仿真實(shí)驗(yàn)，證明了此算法在收斂速度和精度上的優(yōu)越性。

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

2.1 概率積分法預(yù)計(jì)模型

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

(8)

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

根據(jù)疊加原理，整個(gè)開(kāi)采引起的A點(diǎn)的下沉量計(jì)算式見(jiàn)式(9)。

(9)

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

對(duì)于矩形工作面，各參數(shù)量計(jì)算見(jiàn)式(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為工作面傾向斜長(zhǎng)；D3為工作面走向長(zhǎng)度；s1、s2為下山和上山方向的拐點(diǎn)偏移距；s3、s4為走向左、右邊界的拐點(diǎn)偏移距；l為走向工作面的計(jì)算長(zhǎng)度；L為傾向工作面的計(jì)算長(zhǎng)度，其他參數(shù)同上。

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

(19)

式中：fp為目標(biāo)函數(shù)；Wi為實(shí)測(cè)下沉值；wi為預(yù)計(jì)下沉值；m為觀測(cè)點(diǎn)數(shù)。

2.2 反演流程

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

圖1 改進(jìn)步長(zhǎng)果蠅算法反演概率積分法預(yù)計(jì)參數(shù)流程

3 實(shí)例分析

3.1 試驗(yàn)區(qū)概況

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

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

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

3.2 反演結(jié)果分析

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

圖2 某礦22101工作面地表移動(dòng)觀測(cè)站布設(shè)

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

表2 參數(shù)反演結(jié)果對(duì)比

圖3 擬合下沉值與實(shí)測(cè)下沉值對(duì)比

表3 反演結(jié)果對(duì)比

點(diǎn)號(hào)實(shí)際下沉值/mm擬合下沉值/mm絕對(duì)誤差/mm相對(duì)誤差/%傳統(tǒng)算法改進(jìn)果蠅算法傳統(tǒng)算法改進(jìn)果蠅算法傳統(tǒng)算法改進(jìn)果蠅算法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 結(jié) 語(yǔ)

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