基于狼群算法的概率積分法模型參數(shù)反演方法研究

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

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

我國(guó)近70%的天然能源來(lái)自于煤炭，近年來(lái)煤炭資源在能源消費(fèi)中的占比逐漸降低，但短期內(nèi)仍難以被取代。煤炭資源雖使用廣泛、成本較低，但煤炭開(kāi)采會(huì)造成開(kāi)采區(qū)巖層發(fā)生移動(dòng)、變形和破壞，為提前實(shí)施有效的預(yù)防措施，必須準(zhǔn)確地預(yù)計(jì)開(kāi)采沉陷的影響范圍及變形值大小[1]。我國(guó)礦山開(kāi)采沉陷預(yù)計(jì)的主要方法是概率積分法，由于模型簡(jiǎn)單、結(jié)果精準(zhǔn)等優(yōu)點(diǎn)，得到廣泛應(yīng)用[2]。提高開(kāi)采沉陷預(yù)測(cè)精度[3]主要有兩個(gè)途徑：一是提高預(yù)測(cè)模型精度；二是提高預(yù)測(cè)參數(shù)的準(zhǔn)確性與可靠性。在預(yù)測(cè)模型正確的情況下，概率積分法模型參數(shù)的精確性會(huì)直接影響到預(yù)計(jì)的結(jié)果[4]。因此，研究基于地表移動(dòng)實(shí)測(cè)資料的概率積分參數(shù)穩(wěn)健反演方法對(duì)開(kāi)采沉陷精準(zhǔn)預(yù)計(jì)具有重要意義。

概率積分法模型參數(shù)反演歷程[5]中，經(jīng)歷了直接反演、實(shí)驗(yàn)設(shè)計(jì)方法反演、優(yōu)化算法反演、智能算法反演的過(guò)程。智能算法是現(xiàn)階段概率積分法參數(shù)反演的主要方法，查劍鋒等[6]構(gòu)建了基于遺傳算法的概率積分法預(yù)計(jì)參數(shù)反演，抗干擾能力強(qiáng)，在求取參數(shù)方面具有一定的準(zhǔn)確性、可靠性和穩(wěn)定性，可以避免計(jì)算發(fā)散、陷入局部最優(yōu)解；徐夢(mèng)強(qiáng)等[7]提出采用粒子群算法反演概率積分法參數(shù)，搜索速度較快、效率高且算法簡(jiǎn)單；蘇軍明等[8]利用模擬退火算法反演概率積分法參數(shù)，精度較高，并且計(jì)算過(guò)程簡(jiǎn)單、通用；此后又有專家學(xué)者將果蠅算法[9]、量子遺傳算法[10]引入到概率積分法參數(shù)反演中，取得了良好的效果。上述智能算法具有求參精度高、計(jì)算速度快、可用于非線性問(wèn)題等良好的尋優(yōu)特性。為獲得更好的尋優(yōu)特性，吳虎勝等[11]提出了狼群算法。狼群算法作為一種新智能算法，可以運(yùn)用到復(fù)雜的、非線性的、多參數(shù)的尋優(yōu)過(guò)程中，求取參數(shù)時(shí)能夠保證較高準(zhǔn)確性、可靠性和穩(wěn)定性，同時(shí)能夠以較大的速率找到最優(yōu)解，可以從多個(gè)點(diǎn)出發(fā)進(jìn)行搜索，在避免點(diǎn)與點(diǎn)之間的影響的同時(shí)，防止陷入局部最優(yōu)解。目前，狼群算法被實(shí)際應(yīng)用在求解多維背包問(wèn)題、水電站水庫(kù)優(yōu)化調(diào)度以及旅行商等問(wèn)題中[12-15]，在概率積分法參數(shù)反演中還未得到應(yīng)用。

因此，本文擬將狼群算法引入到概率積分法參數(shù)反演中，構(gòu)建基于WPA的概率積分法模型參數(shù)反演方法，通過(guò)模擬實(shí)驗(yàn)系統(tǒng)的討論方法的可靠性，最后開(kāi)展工程應(yīng)用實(shí)驗(yàn)，并且討論該算法在概率積分法參數(shù)反演中的可靠性和適用性。

1 基于WPA的概率積分法模型參數(shù)反演方法

1.1 概率積分法模型及其參數(shù)體系

如圖1所示，考慮單元開(kāi)采引起變形的三維情況，整個(gè)工作面開(kāi)采引起坐標(biāo)為(x,y)的地表任意點(diǎn)A的下沉值的計(jì)算過(guò)程見(jiàn)式(1)～(4)。

(1)

W0=mqcosα

(2)

W0(x)=

(3)

W0(y)=

(4)

地表任意點(diǎn)A沿φ方向的水平移動(dòng)值，其計(jì)算公式見(jiàn)式(5)和式(6)。

U(x,y,φ)=bri(x,y,φ)

(5)

(6)

式中：φ為x軸的正向逆時(shí)針到指定方向的角值；r為主要影響半徑；b為水平移動(dòng)系數(shù)。

圖1 單元開(kāi)采地表任意點(diǎn)下沉預(yù)計(jì)原理圖Fig.1 Schematic diagram of surface subsidenceprediction at any point of unit mining

綜上可知，概率積分法模型的參數(shù)體系為B=[q,tanβ,b,θ,S1,S2,S3,S4]，參數(shù)體系B的解算精度直接決定開(kāi)采沉陷的預(yù)計(jì)精度。

1.2 基于WPA的概率積分法模型參數(shù)反演方法

1.2.1 WPA理論基礎(chǔ)

人工狼群狩獵區(qū)域是一個(gè)N×D的歐氏空間(N為人工狼總數(shù)，D為待求參數(shù)維度)，任意一匹人工狼的位置為Wn=(wn1,wn2,…,wnd,…,wnD)。在人工狼群中，每一匹人工狼都是隨機(jī)生成的，隨機(jī)生成人工狼群的公式見(jiàn)(7)和式(8)。

Wn=(wn1,wn2,…,wnd,…,wnD),

(1≤n≤N,1≤d≤D)

(7)

wnd=wmin+rand×(wmax-wmin)

(8)

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

頭狼產(chǎn)生規(guī)則：在人工狼群中，頭狼感知?dú)馕稘舛茸畲螅^狼不斷變化且不參與游行、奔襲、圍攻行為。

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

(9)

式中：wnd為開(kāi)始搜尋獵物前，探狼在第d維空間中的位置；stepad為探狼在第d維空間中的游走步長(zhǎng)。

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

(10)

圍攻行為：猛狼、探狼接收頭狼的信號(hào)后，開(kāi)始圍攻獵物，若每次攻擊后感知?dú)馕稘舛瘸^(guò)原位置感知?dú)馕稘舛龋瑒t更新該人工狼。第k+1次攻擊后的位置見(jiàn)式(11)。

(11)

狼群更新規(guī)則：圍攻結(jié)束后，淘汰感知?dú)馕稘舛容^小的人工狼，同時(shí)隨機(jī)生成相同數(shù)量的人工狼。

1.2.2 參數(shù)反演方法構(gòu)建

設(shè)開(kāi)采區(qū)任意一點(diǎn)(x,y)的實(shí)測(cè)下沉值為Wxy實(shí)、水平移動(dòng)值為Uxy實(shí)，將狼群算法反演參數(shù)帶入概率積分法模型，得到的該點(diǎn)下沉值為Wxy預(yù)、水平移動(dòng)值為Uxy預(yù)，以實(shí)測(cè)值與預(yù)計(jì)值之差平方和最小為準(zhǔn)則，構(gòu)建適應(yīng)度函數(shù)為式(12)。

F=∑[(Wxy實(shí)-Wxy預(yù))2+(Uxy實(shí)-Uxy預(yù))2]

(12)

令概率積分法參數(shù)B=[q,tanβ,b,θ,S1,S2,S3,S4]為人工狼狩獵區(qū)域，以適應(yīng)度函數(shù)值作為判斷人工狼感知?dú)馕稘舛却笮〉囊罁?jù)。在每一次迭代過(guò)程中，人工狼群在狩獵區(qū)域中感知?dú)馕稘舛炔?duì)獵物(概率積分法參數(shù))進(jìn)行捕獲，達(dá)到終止條件后，輸出最優(yōu)概率積分法參數(shù)。 求參主要流程如下所述。

1) 設(shè)置狼群數(shù)量N、最大迭代次數(shù)Dmax、最大游走次數(shù)Tmax等參數(shù)；輸入觀測(cè)點(diǎn)實(shí)測(cè)下沉值Wxy實(shí)、水平移動(dòng)值Uxy實(shí)及概率積分法參數(shù)的波動(dòng)范圍ΔB=[Δq,Δtanβ,Δb,Δθ,ΔS1,ΔS2,ΔS3,ΔS4]，根據(jù)式(7)和式(8)生成人工狼群Bi=[qi,tanβi,bi,θi,S1i,S2i,S3i,S4i](i=1…N)。

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

3) 選擇除頭狼外、感知?dú)馕稘舛容^大的人工狼為探狼，探狼根據(jù)式(9)進(jìn)行游走行為，直至達(dá)到最大游走次數(shù)Tmax或更新頭狼。

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

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

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

7) 判斷是否達(dá)到精度要求或最大迭代次數(shù)Dmax，若達(dá)到，則輸出頭狼位置，即概率積分法預(yù)計(jì)參數(shù)，否則轉(zhuǎn)至第二步驟。

基于WPA的概率積分模型參數(shù)反演方法程序流程圖如圖2所示。

圖2 WPA參數(shù)反演程序流程圖Fig.2 Flow chart of WPA parameter inversion program

2 模擬實(shí)驗(yàn)

2.1 模擬實(shí)驗(yàn)工作面概況

以淮南礦區(qū)煤系地層為背景，設(shè)計(jì)模擬實(shí)驗(yàn)工作面參數(shù)為：煤層采厚m=3.0 m，煤層傾角α=5°，走向長(zhǎng)度D3=800 m，傾向長(zhǎng)度D1=500 m，平均采深H=400 m，工作面采用全部垮落法管理頂板。假設(shè)模擬工作面的概率積分法模型參數(shù)為：下沉系數(shù)q=0.8，主要影響正切角tanβ=2.5，水平移動(dòng)系數(shù)b=0.25，開(kāi)采影響傳播角θ=85°，上拐點(diǎn)、下拐點(diǎn)、左拐點(diǎn)、右拐點(diǎn)偏移距分別為S1=60 m、S2=60 m、S3=60 m、S4=60 m。本次模擬實(shí)驗(yàn)在移動(dòng)盆地內(nèi)設(shè)計(jì)了沿走向、傾向的兩條主斷面的移動(dòng)和變形的觀測(cè)線，相鄰兩個(gè)監(jiān)測(cè)點(diǎn)之間相距20 m，走向觀測(cè)線1 180 m，共60個(gè)觀測(cè)點(diǎn)，傾向觀測(cè)線880 m，共45個(gè)觀測(cè)點(diǎn)。模擬工作面設(shè)計(jì)圖如圖3所示。

2.2 狼群算法參數(shù)反演準(zhǔn)確性

使用WPA進(jìn)行概率積分參數(shù)反演，并進(jìn)行模擬實(shí)驗(yàn)參數(shù)設(shè)計(jì)值與模擬實(shí)驗(yàn)參數(shù)反演值的對(duì)比，通過(guò)參數(shù)平均值、反演參數(shù)相對(duì)誤差以及反演參數(shù)中誤差來(lái)體現(xiàn)狼群算法反演參數(shù)的準(zhǔn)確性。智能算法在一定范圍內(nèi)具有隨機(jī)性，為避免偶然誤差的影響，在相同的實(shí)驗(yàn)條件下連續(xù)進(jìn)行了10次實(shí)驗(yàn)。設(shè)置人工狼數(shù)量為50，最大迭代次數(shù)為50，實(shí)驗(yàn)結(jié)果見(jiàn)表1。

由表1可知,在參數(shù)相對(duì)誤差方面，q、tanβ、b、θ的參數(shù)相對(duì)誤差均不超過(guò)0.3%，S1、S2、S3、S4對(duì)求參敏感性較低,但參數(shù)相對(duì)誤差均小于2.1%；WPA反演概率比積分法模型參數(shù)具有較高的準(zhǔn)確度。在參數(shù)中誤差方面，q、tanβ、b、θ的參數(shù)中誤差均不超過(guò)1.03，S1、S2、S3、S4的參數(shù)中誤差均小于8.9 m；WPA反演概率比積分法模型參數(shù)具有較高的精度。綜上所述，WPA反演概率積分法模型參數(shù)具有較好的穩(wěn)健性。

圖3 模擬工作面設(shè)計(jì)圖Fig.3 Simulated working face design drawing

表1 WPA反演參數(shù)值與設(shè)計(jì)值對(duì)比Table 1 Comparison of WPA inversion parameter valueswith design values

為更加直觀地看出下沉與水平移動(dòng)的擬合情況，得到10次模擬實(shí)驗(yàn)反演參數(shù)均值計(jì)算出的下沉值曲線和水平移動(dòng)值曲線見(jiàn)圖4和圖5。

由圖4和圖5可以看出，下沉值和水平移動(dòng)值的實(shí)測(cè)曲線與擬合曲線基本一致，點(diǎn)位實(shí)測(cè)值與擬合值間的絕對(duì)值誤差均不超過(guò)5 mm，達(dá)到預(yù)期精度。

2.3 基于WPA的概率積分法模型參數(shù)反演方法的可靠性研究

為進(jìn)一步驗(yàn)證基于WPA算法的概率積分法模型參數(shù)反演方法的可靠性，對(duì)模擬實(shí)驗(yàn)10次反演結(jié)果中的概率積分法模型參數(shù)分別進(jìn)行波動(dòng)性分析，對(duì)比WPA反演參數(shù)與設(shè)計(jì)參數(shù)，結(jié)果見(jiàn)圖6。

圖4 模擬實(shí)驗(yàn)WPA擬合下沉值與實(shí)測(cè)下沉值對(duì)比圖Fig.4 The comparison diagram of the simulated WPAfitted subsidence value and the measuredsubsidence value

圖5 模擬實(shí)驗(yàn)WPA擬合水平移動(dòng)值與實(shí)測(cè)水平移動(dòng)值對(duì)比圖Fig.5 The comparison diagram of the simulated WPAfitted horizontal movement value and the measuredhorizontal movement value

圖6 模擬實(shí)驗(yàn)10次反演參數(shù)的波動(dòng)情況Fig.6 The fluctuation of inversion parameters was simulated for 10 times

從圖6可以看出，WPA反演參數(shù)波動(dòng)范圍較小，反演參數(shù)較為接近于真值，基于WPA算法的概率積分法模型參數(shù)反演方法具有良好的可靠性。

3 工程實(shí)例

3.1 礦區(qū)概況

淮南顧橋礦南二采區(qū)1414(1)工作面是顧橋煤礦南區(qū)的第一個(gè)采面，工作面從礦床邊界開(kāi)始向巷道往回開(kāi)采，采用機(jī)械化掘進(jìn)，垮落法管理頂板，一次采全高。工作面長(zhǎng)約2 120 m，寬約251 m，走向?yàn)槌浞植蓜?dòng)，傾向?yàn)榉浅浞植蓜?dòng)，總體為非充分采動(dòng)。平均采高為3.0 m，煤層傾角在平均為5°，工作面平均深度為735 m，在距離切眼1 144 m處以及停采線976 m處布置傾向觀測(cè)線，傾向線長(zhǎng)度為1 500 m，相鄰點(diǎn)之間點(diǎn)間距為30 m，共布設(shè)了3個(gè)控制點(diǎn)和50個(gè)監(jiān)測(cè)點(diǎn)。 在下山方向偏離工作面中心線的距離39 m處設(shè)置走向觀測(cè)線，長(zhǎng)度為3 480 m，相鄰點(diǎn)之間點(diǎn)間距為30 m/60 m，共布設(shè)了3個(gè)控制點(diǎn)和95個(gè)監(jiān)測(cè)點(diǎn)。

3.2 實(shí)驗(yàn)結(jié)果及分析

以淮南顧橋礦南二采區(qū)1414(1)工作面的實(shí)測(cè)數(shù)據(jù)為基礎(chǔ)，將WPA反演參數(shù)均值計(jì)算所得的下沉值、水平移動(dòng)值與監(jiān)測(cè)點(diǎn)的下沉值、水平移動(dòng)值進(jìn)行對(duì)比,為避免實(shí)驗(yàn)的偶然性，在相同條件下進(jìn)行10次實(shí)驗(yàn)，然后分別計(jì)算反演參數(shù)平均值、中誤差，結(jié)果見(jiàn)表2。

由表2可知，在參數(shù)中誤差方面，q、tanβ、b、θ的參數(shù)中誤差均不超過(guò)0.89，S1、S2、S3、S4的參數(shù)中誤差均不超過(guò)4.5 m；WPA反演顧橋南礦1414(1)工作面的概率積分法參數(shù)為：q=0.935 7，tanβ=1.988 1，b=0.422 8，θ=84.138 3°，S1=-9.136 6 m，S2=-21.283 9 m，S3=55.954 4 m，S4=34.111 3 m。

將10次工程實(shí)例反演參數(shù)均值帶入概率積分法模型，計(jì)算出的下沉值曲線和水平移動(dòng)值曲線見(jiàn)圖7和圖8。

表2 WPA反演參數(shù)工程實(shí)例結(jié)果Table 2 Engineering example results of WPAinversion parameters

圖7 工程實(shí)例WPA擬合下沉值與實(shí)測(cè)下沉值對(duì)比圖Fig.7 The comparison diagram of WPA fitted subsidencevalue and the measured subsidence value

圖8 工程實(shí)例WPA擬合水平移動(dòng)值與實(shí)測(cè)水平移動(dòng)值對(duì)比圖Fig.8 The comparison diagram of WPA fitted horizontalmoving value and measured horizontal moving value

從圖7和圖8可以看出，下沉值和水平移動(dòng)值的實(shí)測(cè)曲線與擬合曲線基本一致，點(diǎn)位實(shí)測(cè)值與擬合值間的絕對(duì)值誤差均不超過(guò)300 mm。WPA反演概率積分法參數(shù)取均值后，其下沉與水平移動(dòng)值的擬合中誤差為114.81 mm，其中，下沉值擬合中誤差約為最大下沉值的2.7%，水平移動(dòng)值擬合中誤差約為最大水平移動(dòng)值的8.6%，滿足工程應(yīng)用精度[18]。

4 結(jié) 論

1) 模擬實(shí)驗(yàn)結(jié)果表明，WPA反演參數(shù)q、tanβ、b、θ的參數(shù)相對(duì)誤差均不超過(guò)0.3%，參數(shù)中誤差均不超過(guò)1.03；S1、S2、S3、S4對(duì)求參敏感性較低，但參數(shù)相對(duì)誤差均小于2.1%，S1、S2、S3、S4的參數(shù)中誤差均不超過(guò)8.9 m。WPA反演參數(shù)波動(dòng)范圍較小，反演參數(shù)較為接近于真值。以上證明了WPA在概率積分法模型參數(shù)反演中，具有良好的穩(wěn)健性，且反演參數(shù)具有較高的準(zhǔn)確性、精度。

2) 利用本文提出的基于WPA的概率積分法模型參數(shù)反演方法對(duì)淮南顧橋礦南二采區(qū)1414(1)工作面概率積分法參數(shù)反演，將由反演結(jié)果計(jì)算的下沉值、水平移動(dòng)擬合值與真實(shí)值對(duì)比，獲得的下沉值與水平移動(dòng)擬合中誤差為114.81 mm，下沉值擬合中誤差約為最大下沉值的2.7%，水平移動(dòng)值擬合中誤差約為最大水平移動(dòng)值的8.6%，滿足工程應(yīng)用精度。