基于改進(jìn)狼群算法的概率積分法模型參數(shù)反演方法

李靖宇，王 磊，江克貴，滕超群

( 安徽理工大學(xué) 測繪學(xué)院，安徽 淮南 232001 )

我國能源資源結(jié)構(gòu)具有“富煤、貧油、少氣”的特點(diǎn)，其中煤炭資源在一次能源結(jié)構(gòu)的生產(chǎn)和消費(fèi)中所占比例高達(dá)76%和66%，因此在可再生資源形成規(guī)模之前，煤炭將長期作為主體能源擔(dān)負(fù)著國家能源安全和經(jīng)濟(jì)持續(xù)發(fā)展的重任[1]。高強(qiáng)度的煤炭開采勢必會造成礦區(qū)環(huán)境災(zāi)害問題，如地面塌陷、工程地質(zhì)損害等，開采沉陷預(yù)計(jì)對于“三下”采煤、礦區(qū)環(huán)境災(zāi)害防治可起到指導(dǎo)性作用，因而開采沉陷預(yù)計(jì)一直是我國礦山領(lǐng)域的熱點(diǎn)問題。概率積分法是我國普遍采用的開采沉陷預(yù)計(jì)方法，該方法已積累大量的實(shí)踐經(jīng)驗(yàn)，形成了相對成熟的預(yù)計(jì)參數(shù)體系[2]，如何基于實(shí)測數(shù)據(jù)，精確、快速、可靠地獲得概率積分法模型參數(shù)，一直是該方法應(yīng)用的難點(diǎn)。

概率積分法模型參數(shù)反演主要經(jīng)歷了特征點(diǎn)法、線性近似法[3]、正交實(shí)驗(yàn)法[4]，特征點(diǎn)法因曲線特征點(diǎn)難以精確確定進(jìn)而導(dǎo)致求參誤差大；線性近似法對模型參數(shù)初值要求極高，若初值精度不高，則會導(dǎo)致求參模型發(fā)散；正交實(shí)驗(yàn)法求參速度慢且不易借助計(jì)算機(jī)實(shí)施；針對概率積分法模型參數(shù)反演中存在的問題，一些專家學(xué)者提出采用智能優(yōu)化算法反演概率積分法模型參數(shù)。查劍鋒[5]等利用遺傳算法進(jìn)行概率積分法模型參數(shù)反演，此后又有專家學(xué)者將粒子群算法[6]、人工蜂群算法[7]、果蠅算法[8]引入到了概率積分法模型參數(shù)反演中，為概率積分模型參數(shù)反演提供了新途徑。近年來，吳虎勝[9]等提出了狼群算法，與傳統(tǒng)智能優(yōu)化算法相比，狼群算法具有收斂速度較快、可有效避免局部最優(yōu)解、統(tǒng)計(jì)結(jié)果精度較高等良好的尋優(yōu)特性；文獻(xiàn)[10-12]對狼群算法提出了改進(jìn)方案，進(jìn)一步提高了狼群算法的尋優(yōu)特性。目前，狼群算法已在移動(dòng)機(jī)器人路徑規(guī)劃[13]、橋起主梁優(yōu)化[14]、井下排水采氣超聲速噴嘴優(yōu)化[15]等領(lǐng)域得到了廣泛應(yīng)用，但鮮有概率積分法模型參數(shù)反演方面的案例報(bào)道。

鑒于此，本文在分析狼群算法優(yōu)缺點(diǎn)的基礎(chǔ)上提出二次游走、變異行為改進(jìn)策略，形成了改進(jìn)狼群算法( IWPA )，與狼群算法相比，改進(jìn)狼群算法的搜索精細(xì)程度、跳出局部最優(yōu)解能力得到了有效提高，并將改進(jìn)狼群算法引入概率積分法模型參數(shù)反演中，構(gòu)建改進(jìn)狼群算法的概率積分法模型參數(shù)反演方法( MIWPA )。

1 反演方法

1.1 基本狼群算法

人工狼群狩獵區(qū)域[15]為一個(gè)N×D的歐氏空間( N為人工狼總數(shù)，D為待求參數(shù)維度 )，任意一匹人工狼為Wn＝( wn1，wn2，…，wnd，…，wnD)。在人工狼群中，每一匹人工狼都是隨機(jī)生成的，隨機(jī)生成人工狼的公式為

式中，wnd為第n匹人工狼在第d維空間中的位置( 即第n組隨機(jī)參數(shù)中的第d個(gè)參數(shù) )；wmax，wmin分別為wnd中的最大值、最小值；rand為-1至1之間均勻分布的隨機(jī)數(shù)。

( 1 ) 頭狼產(chǎn)生規(guī)則

在人工狼群中，頭狼感知?dú)馕稘舛茸畲?，頭狼不斷變化且不參與游行、奔襲、圍攻行為。

( 2 ) 游走行為

探狼是除頭狼外、感知?dú)馕稘舛容^大的人工狼。探狼每次游走向不同方向前進(jìn)，每次前進(jìn)后返回原來位置，若探狼前進(jìn)后感知?dú)馕稘舛瘸^頭狼感知?dú)馕稘舛龋瑒t更新頭狼并進(jìn)行奔襲行為，若不是，則繼續(xù)游走直至更新頭狼或達(dá)到最大游走次數(shù)；游走結(jié)束后，若探狼感知最大氣味濃度大于當(dāng)前位置感知?dú)馕稘舛?，則更新探狼。第k次前進(jìn)后，探狼的位置為

式中，wnd為開始搜尋獵物前，探狼在第d維空間中的位置；step ad為探狼在第d維空間中的游走步長。

( 3 ) 奔襲行為

猛狼是距頭狼較近的人工狼。猛狼在奔襲過程中，若感知?dú)馕稘舛瘸^頭狼感知?dú)馕稘舛?，則更新頭狼且繼續(xù)奔襲行為，直至猛狼與頭狼之間的距離滿足要求。猛狼奔襲后的位置為

式中，wnd為猛狼在第d維空間中的位置；step bd為猛狼在第d維空間中的奔襲步長為頭狼在第d維空間中的位置。

( 4 ) 圍攻行為

猛狼、探狼接收頭狼的信號后，開始圍攻獵物，若每次攻擊后感知?dú)馕稘舛瘸^原位置感知?dú)馕稘舛?，則更新該人工狼。第k+1次攻擊后的位置為

( 5 ) 狼群更新規(guī)則

圍攻結(jié)束后，淘汰感知?dú)馕稘舛容^小的人工狼，同時(shí)根據(jù)式( 1 )，( 2 )隨機(jī)生成相同數(shù)量的人工狼。

1.2 狼群算法改進(jìn)策略

( 1 ) 二次游走

在狼群算法游走行為中，存在下述問題：首先，對于某一探狼而言，其游走方向數(shù)固定不變；其次，當(dāng)探狼游走至獵物附近時(shí)，游走步長可能會超出探狼與頭狼間的距離，游走后無法感知更大的氣味濃度，此時(shí)探狼與獵物間雖具有一定距離，但無法進(jìn)一步靠近獵物；最后，當(dāng)探狼與獵物間的初始距離較大時(shí)，游走后感知更大的氣味濃度并向獵物前進(jìn)，但游走步長可能略小，直至游走結(jié)束，距離獵物仍有一定距離。假如探狼A，B游走方向數(shù)為4，則狼群算法游走方式如圖1所示。

針對上述問題，為提高狼群算法的搜索精細(xì)程度，提出了二次游走改進(jìn)策略，其思想如下：

① 將狼群算法探狼游走視為一次游走，探狼選擇以一次游走感知?dú)馕稘舛茸畲蟮姆轿粸槠瘘c(diǎn)，進(jìn)行二次游走。

圖1 WPA探狼游走Fig.1 WPA spies the wolf

② 探狼二次游走時(shí)，根據(jù)一次游走結(jié)果，做出不同的反應(yīng)：對于探狼A而言，其二次游走與一次游走相比，方向數(shù)相同但步長較長；對于探狼B而言，其二次游走與一次游走相比，方向數(shù)較多且步長較短。

假如探狼A，B游走方向數(shù)為4，則改進(jìn)狼群算法游走方式如圖2所示。

圖2 IWPA探狼游走Fig.2 IWPA spies away

③ 探狼二次游走中，每一次前進(jìn)后，返回二次游走前位置，第m次前進(jìn)后的位置為

式中，wnd1為探狼一次游走感知?dú)馕稘舛容^大的方位為探狼第d維空間中的二次游走步長。

( 2 ) 變異行為

狼群算法為跳出局部最優(yōu)解，對狼群進(jìn)行狼群更新：淘汰部分感知?dú)馕稘舛容^低的人工狼并隨機(jī)生成等數(shù)量的人工狼；但生成的大部分人工狼適應(yīng)度函數(shù)較低或距離頭狼較遠(yuǎn)，無法成為探狼或猛狼，因此沒有參與到狼群游走、奔襲、圍攻行為中，并在下一次迭代前的狼群更新中被淘汰，從而導(dǎo)致狼群算法跳出局部最優(yōu)解能力受到限制。針對上述問題，提出變異行為改進(jìn)策略，其思想如下。

① 召喚行為后，狼群并不立刻進(jìn)行圍攻行為，而是將頭狼與其他人工狼按比例混合，如式( 7 )所示。

② 混合后，若人工狼感知的氣味濃度變大，則更新其位置。

③ 人工狼更新后，重新選擇頭狼、探狼、猛狼并進(jìn)行圍攻行為。

1.3 反演方法構(gòu)建

根據(jù)1.2節(jié)改進(jìn)的狼群算法，結(jié)合開采沉陷概率積分預(yù)計(jì)模型，提出了基于改進(jìn)狼群算法的概率積分法模型參數(shù)反演方法，其主要過程如下：

Step 1：地面實(shí)測下沉值、水平移動(dòng)值分別記為Ss，Ms，利用概率積分法模型參數(shù)計(jì)算的預(yù)計(jì)下沉值、水平移動(dòng)值分別記為Sy，My，以實(shí)測值與預(yù)計(jì)值之差平方和最小為準(zhǔn)則，構(gòu)建適應(yīng)度函數(shù)為

式中，Ssi，Msi，Syi，Myi分別為地面觀測點(diǎn)( Xi，Yi)對應(yīng)的Ss，Ms，Sy，My。

Step 2：設(shè)置狼群數(shù)量N、最大迭代次數(shù)Dmax、最大游走次數(shù)Tmax等參數(shù)；給定概率積分法模型參數(shù)中心值B0＝[q0，tan β0，b0，θ0，S10，S20，S30，S40]，各概率積分法模型參數(shù)的波動(dòng)范圍ΔB＝[Δq，Δtan β，Δb，Δθ，ΔS1，ΔS2，ΔS3，ΔS4]，生成人工狼群Bi＝[qi，tan βi，bi，θi，S1i，S2i，S3i，S4i]( i＝1，…，N )。

Step 3：以適應(yīng)度函數(shù)值為判斷人工狼感知?dú)馕稘舛萗n大小的依據(jù)，適應(yīng)度函數(shù)值越小則人工狼感知?dú)馕稘舛仍酱?，選擇感知?dú)馕稘舛茸畲蟮娜斯だ菫轭^狼，記其位置、感知?dú)馕稘舛确謩e為wlead，Ylead。

Step 4：選擇除頭狼外、感知?dú)馕稘舛容^大的人工狼為探狼，探狼根據(jù)式( 3 )進(jìn)行一次游走后，根據(jù)式( 6 )進(jìn)行二次游走，直至二次游走后探狼感知?dú)馕稘舛萗n2＞Ylead或游走次數(shù)T＝Tmax后，進(jìn)行奔襲行為，若Yn2大于一次游走后該探狼感知?dú)馕稘舛萗n，則更新該探狼。

Step 5：選擇距離頭狼較近的人工狼為猛狼，猛狼根據(jù)式( 4 )向頭狼奔襲，若奔襲過程中猛狼感知?dú)馕稘舛萗n＞Ylead，則更新頭狼，繼續(xù)召喚行為；若Yn＜Ylead，則猛狼繼續(xù)奔襲，直至猛狼與頭狼間距離Dn小于要求距離D。

Step 6：根據(jù)式( 7 )，對人工狼進(jìn)行變異，若變異后感知?dú)馕稘舛萗nB大于變異前感知?dú)馕稘舛萗n，則更新該人工狼，變異結(jié)束后重新選擇頭狼、探狼、猛狼。

Step 7：根據(jù)式( 5 )，探狼聯(lián)合猛狼對獵物進(jìn)行圍攻。

Step 8：對頭狼進(jìn)行更新，淘汰部分感知?dú)馕稘舛容^小的人工狼，并根據(jù)式( 1 )，( 2 )隨機(jī)生成等數(shù)量人工狼。

Step 9：判斷是否達(dá)到精度要求或最大迭代次數(shù)Dmax，若達(dá)到，則輸出頭狼位置，即概率積分法模型參數(shù)，否則轉(zhuǎn)至Step 3。

基于改進(jìn)狼群算法的概率積分法模型參數(shù)反演方法( MIWPA )的運(yùn)行流程如圖3所示。

圖3 MIWPA流程Fig.3 Flow chart of MIWPA

2 模擬試驗(yàn)

2.1 工作面概況

以淮南礦區(qū)地質(zhì)采礦條件為背景，模擬工作面走向開采長度D3＝800 m，傾向開采長度D1＝300 m，煤層平均采厚m＝3 m，煤層傾角α＝5°，平均采深400 m，采用全部垮落法管理頂板開采。在工作面走向每隔20 m布設(shè)1個(gè)監(jiān)測點(diǎn)，共51個(gè)監(jiān)測點(diǎn)；沿工作面傾向每隔20 m布設(shè)1個(gè)監(jiān)測點(diǎn)，共36個(gè)監(jiān)測點(diǎn)，模擬工作面如圖4所示。礦區(qū)概率積分法模型參數(shù)為：下沉系數(shù)q＝0.8，水平移動(dòng)系數(shù)b＝0.25，主要影響角正切tan β＝2.5，最大下沉角θ＝85°，上、下、左、右拐點(diǎn)偏移距為S1＝S2＝S3＝S4＝60 m。

圖4 模擬工作面及地表監(jiān)測點(diǎn)布置Fig.4 Schematic diagram of simulated working face and surface monitoring points

2.2 改進(jìn)狼群算法參數(shù)估計(jì)準(zhǔn)確性

分別利用MWPA，MIWPA進(jìn)行概率積分法模型參數(shù)反演，試驗(yàn)設(shè)計(jì)方案如下：參數(shù)q，tan β，b，θ，S1，S2，S3，S4波動(dòng)范圍分別為0.5～1.1，2.0～3.0，0.05～0.45，82～88，40～80，40～80，40～80，40～80，為了防止試驗(yàn)結(jié)果出現(xiàn)偶然性，進(jìn)行10次概率積分法模型參數(shù)反演試驗(yàn)，計(jì)算參數(shù)中誤差及參數(shù)相對誤差，試驗(yàn)結(jié)果見表1。

表1 MWPA 和MIWPA 參數(shù)估計(jì)值與設(shè)計(jì)值對比Table 1 Estimated values of MWPA and MIWPA parameters were compared with the design values

由表1可知：在參數(shù)相對誤差方面，MWPA反演參數(shù)相對誤差最大不超過6.3%，MIWPA反演參數(shù)相對誤差最大控制在3.4%以內(nèi)。在參數(shù)擬合中誤差方面，MIWPA反演參數(shù)擬合中誤差均小于MWPA，MWPA反演參數(shù)中誤差最大不超過9.39，MIWPA反演參數(shù)中誤差最大不超過4.02。

試驗(yàn)結(jié)果表明，與MWPA的反演參數(shù)結(jié)果相比，MIWPA方法更加穩(wěn)健，且參數(shù)反演參數(shù)精度大幅提高。

2.3 改進(jìn)狼群算法的可靠性

為進(jìn)一步證明MIWPA的可靠性，基于模擬試驗(yàn)，分別利用MWPA和MIWPA反演概率積分模型參數(shù)10次，反演參數(shù)結(jié)果的波動(dòng)情況如圖5所示。由圖5可知，與MWPA相比較，MIWPA反演參數(shù)波動(dòng)范圍較小，反演參數(shù)更接近于真實(shí)值。

綜上所述，MIWPA反演參數(shù)的穩(wěn)定性、精度均優(yōu)于MWPA，MIWPA反演概率積分法模型參數(shù)具有更好的可靠性。

圖5 10次模擬試驗(yàn)反演參數(shù)波動(dòng)情況Fig.5 Fluctuation of parameters was inversed in 10 simulation experiments

3 工程實(shí)例

3.1 礦區(qū)概況

淮南顧橋南礦1414( 1 )為顧橋煤礦南區(qū)的第1個(gè)回采工作面，工作面從礦床邊界開始向巷道往回開采，采用機(jī)械化掘進(jìn)，垮落法管理頂板，一次采全高，工作面地質(zhì)采礦條件為：走向開采長度2 120 m，傾向開采長度251 m，平均采深735 m，煤層平均采厚3 m，煤層傾角5°。主斷面傾向觀測線布置在距開切眼約1 144 m、終采線約976 m處，傾向線長度為1 500 m，相鄰點(diǎn)之間點(diǎn)間距為30 m，共布設(shè)3個(gè)控制點(diǎn)和50個(gè)監(jiān)測點(diǎn)。主斷面走向觀測線布置在下山方向偏離工作面中心線39 m的方向上，走向線長度約3 480 m，相鄰點(diǎn)之間點(diǎn)間距為30 m/60 m，共布設(shè)3個(gè)控制點(diǎn)和95個(gè)監(jiān)測點(diǎn)。觀測站布設(shè)如圖6所示。

圖6 觀測站布設(shè)Fig.6 Observation station layout

3.2 試驗(yàn)結(jié)果及分析

將顧橋南礦1414( 1 )的實(shí)測數(shù)據(jù)經(jīng)過預(yù)處理，在走向線上選取75個(gè)監(jiān)測點(diǎn)，在傾向線上選取47個(gè)監(jiān)測點(diǎn)，利用MWPA，MIWPA進(jìn)行概率積分法模型參數(shù)反演。為避免試驗(yàn)的偶然性，進(jìn)行了10次試驗(yàn)，然后分別計(jì)算參數(shù)平均值與參數(shù)中誤差。其結(jié)果見表2。

表2 MWPA 和MIWPA 在工程應(yīng)用中求參結(jié)果比較Table 2 Comparison of parameter calculation results between MWPA and MIWPA in engineering application

由表2可知：從參數(shù)中誤差方面來看，MIWPA反演參數(shù)q，tan β，b，θ中誤差控制在0.88以內(nèi)，拐點(diǎn)偏移距S1，S2，S3，S4中誤差最大不超過4.65，MWPA反演參數(shù)q，tan β，b，θ的中誤差控制在1.08以內(nèi)，拐點(diǎn)偏移距S1，S2，S3，S4中誤差最大不超過5.34；表明MIWPA穩(wěn)定性優(yōu)于MWPA。利用MIWPA反演顧橋礦1414( 1 )工作面概率積分法模型參數(shù)為：q＝0.93，tan β＝1.98，b＝0.42，θ＝84.53°，S1＝-12.44 m，S2＝-18.80 m，S3＝55.06 m，S4＝33.98 m。

將MIWPA反演的概率積分法模型參數(shù)代入概率積分法模型，計(jì)算出的下沉曲線、水平移動(dòng)曲線分別如圖7，8所示。

圖7 MIWPA工程實(shí)例參數(shù)估計(jì)下沉值擬合情況Fig.7 Parameter estimation subsidence value fitting of MIWPA project example

圖8 MIWPA工程實(shí)例參數(shù)估計(jì)水平移動(dòng)值擬合情況Fig.8 Fitting condition of parameter estimation horizontal moving value of MIWPA project example

由圖7，8可知，下沉值、水平移動(dòng)值的實(shí)測曲線與擬合曲線基本一致，實(shí)測值與擬合值間的絕對值誤差均不超過300 mm，下沉與水平移動(dòng)擬合中誤差為114.88 mm?？梢姡琈IWPA在反演概率積分法模型參數(shù)中具有實(shí)際應(yīng)用價(jià)值。

4 結(jié) 論

( 1 ) 在分析WPA優(yōu)缺點(diǎn)的基礎(chǔ)上，提出二次游走、改進(jìn)變異行為策略，形成了改進(jìn)的WPA( IWPA )，并將IWPA引入概率積分法模型參數(shù)反演中，構(gòu)建了基于改進(jìn)狼群算法的概率積分法模型參數(shù)反演方法( MIWPA )，并編制了工程應(yīng)用程序。

( 2 ) 通過模擬試驗(yàn)論證了MIWPA的精確性與可靠性。模擬試驗(yàn)結(jié)果表明：MIWPA反演參數(shù)相對誤差、參數(shù)中誤差分別不超過3.4%，4.02，滿足工程要求；MIWPA的準(zhǔn)確性、可靠性均優(yōu)于MWPA。

( 3 ) 利用MIWPA對顧橋礦1414( 1 )工作面概率積分法模型參數(shù)反演。概率積分法模型參數(shù)結(jié)果為q ＝0.93，tan β ＝1.98，b ＝0.42，θ ＝84.53°，S1＝-12.44 m，S2＝-18.80 m，S3＝55.06 m，S4＝33.98 m，下沉值與水平移動(dòng)值擬合中誤差為114.88 mm。