量子粒子群算法反演概率積分法參數(shù)

朱尚軍 王 磊， 魏 濤 蔣 創(chuàng) 江克貴 查劍鋒 孔 川4

（1.安徽理工大學(xué)空間信息與測繪工程學(xué)院，安徽淮南232001；2.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇徐州221116；3.江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室，江蘇徐州221116；4.山東裕隆礦業(yè)集團(tuán)有限公司單家村煤礦，山東曲阜273100）

概率積分法參數(shù)預(yù)計(jì)是我國《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采指南》中常用的開采沉陷預(yù)計(jì)方法，其中概率積分法參數(shù)反演是地表移動觀測數(shù)據(jù)處理的關(guān)鍵環(huán)節(jié)之一［1-2］。概率積分法是我國研究開采沉陷較為成熟且應(yīng)用較為廣泛的預(yù)計(jì)方法之一［3］。該方法是以正態(tài)函數(shù)為影響函數(shù)，采用概率積分方法表示地表移動盆地的下沉、傾斜、曲率、移動和變形［4］。精確反演概率積分法參數(shù)是提高開采沉陷預(yù)計(jì)精度的基礎(chǔ)，但由于概率積分模型是多參數(shù)非線性模型，從而導(dǎo)致參數(shù)反演過程異常復(fù)雜和困難［5］。

反演概率積分法參數(shù)常用的方法有直接求參、擬合求參和智能優(yōu)化算法求參，其中智能優(yōu)化算法可以彌補(bǔ)反演概率積分法參數(shù)過程中計(jì)算量大與過程復(fù)雜等不足。近年來，不少學(xué)者將智能優(yōu)化算法應(yīng)用于概率積分法參數(shù)反演研究，取得了一定的進(jìn)展。吳侃等［6］將模矢法引入概率積分法參數(shù)反演中，在此基礎(chǔ)上開發(fā)了礦區(qū)沉陷預(yù)測預(yù)報(bào)系統(tǒng)，并結(jié)合無人機(jī)、三維激光掃描儀等新技術(shù)建立了一系列變形監(jiān)測模型；查劍鋒等［7］通過研究論證了遺傳算法在概率積分法參數(shù)反演中的適用性與穩(wěn)定性，但該算法易陷入局部最優(yōu)解；徐夢強(qiáng)等［5］將算法簡單、精度高的粒子群優(yōu)化算法引入概率積分法參數(shù)反演中，利用模擬試驗(yàn)和工程應(yīng)用實(shí)例證明了該算法的有效性，但該算法存在易陷入早熟收斂、粒子全局搜索效果較差、收斂速度較慢等不足；陳濤等［8］提出了基于果蠅算法的概率積分法參數(shù)反演方法，通過礦山實(shí)測數(shù)據(jù)驗(yàn)證分析，證明了該算法有助于提高概率積分法參數(shù)反演精度。通過分析研究成果可以發(fā)現(xiàn)：①當(dāng)算法搜索領(lǐng)域不夠充分時(shí)，上述遺傳算法、粒子群優(yōu)化算法、果蠅算法等都會在一定程度上陷入局部最優(yōu)，從而導(dǎo)致算法早熟收斂［2］；②粒子群優(yōu)化算法在參數(shù)反演過程中存在參數(shù)選擇困難、早熟收斂、全局搜索效果較差、算法運(yùn)行耗時(shí)較多等不足［9-13］。通過查閱相關(guān)成果［9-10］發(fā)現(xiàn)，通過加入量子化引力場，有助于提高智能優(yōu)化算法的運(yùn)行效率。魏濤等［9］將量子算法即量子旋轉(zhuǎn)門引入遺傳算法中，提高了算法的準(zhǔn)確性與穩(wěn)定性，并克服了傳統(tǒng)遺傳算法易陷入局部最優(yōu)的缺陷。借鑒上述研究思路，如果將量子化引力場引入粒子群優(yōu)化算法中，利用該算法的運(yùn)行效率優(yōu)勢，將搜索域擴(kuò)大為全局搜索且降低陷入早熟收斂的概率，那么優(yōu)化后的算法（即QPSO算法）在概率積分法參數(shù)反演中相對于傳統(tǒng)粒子群優(yōu)化算法而言，勢必會具有一定的優(yōu)勢。為此，本研究提出了一種基于QPSO算法的概率積分法參數(shù)反演方法，并對其適用性進(jìn)行討論。

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

1.1 理論基礎(chǔ)

粒子群優(yōu)化（PSO）算法作為一種基于群體智能的優(yōu)化算法，是由1995年美國心理學(xué)家James Kennedy和電氣工程師Russell Eberhart提出［14］。PSO算法屬于群體智能優(yōu)化算法范疇，源于模擬鳥群的覓食過程，通過個(gè)體之間相互合作最終找到食物最多最近的點(diǎn)，即優(yōu)化后的搜索結(jié)果［5］。孫俊等［11］提出了量子粒子群（QPSO）算法，該算法取消了速度變量，并且位置更新是完全隨機(jī)迭代，彌補(bǔ)了PSO算法的不足，是一個(gè)有效且較為完善的群體智能優(yōu)化算法。

有別于傳統(tǒng)PSO算法，QPSO算法并非讓粒子按照某一隨機(jī)確定的軌道移動，而是建立一個(gè)量子化引力勢場來束縛粒子運(yùn)動。在建立量子粒子群模型時(shí)，用向量X表示粒子個(gè)體的當(dāng)前位置，即粒子個(gè)體的思維狀態(tài)；隨后經(jīng)過適應(yīng)度函數(shù)判定找到個(gè)體的歷史最優(yōu)位置，即粒子個(gè)體的最優(yōu)歷史經(jīng)驗(yàn)；gbest表示由N個(gè)粒子組成的粒子群中適應(yīng)度函數(shù)最優(yōu)的粒子當(dāng)前位置。在物理力學(xué)中用束縛態(tài)來描述聚集性，粒子運(yùn)動中心存在某種引力勢場是粒子產(chǎn)生束縛態(tài)的主要原因［11］。因此建立一個(gè)量子化引力勢場是量子粒子群優(yōu)化算法的基礎(chǔ)，處于引力勢場的粒子會依據(jù)一定的概率出現(xiàn)在解空間的任何一個(gè)位置，當(dāng)趨近無窮遠(yuǎn)時(shí)，粒子出現(xiàn)的概率趨于0。可見，QPSO算法的隨機(jī)性遠(yuǎn)大于PSO算法，如何建立一個(gè)引力勢場至關(guān)重要［15］。

1.1.1 勢場模型建立

式中，i、j分別為粒子個(gè)數(shù)及維數(shù)；t為迭代次數(shù)；N為粒子總數(shù)；c1、c2為加速常數(shù)；r1、r2為區(qū)間(0,1)上的隨機(jī)數(shù)；Gij(t)為歷史全局最優(yōu)位置。

由式（1）和式（2）可知，當(dāng)c1=c2時(shí)，φij(t)即為在區(qū)間(0,1)上隨機(jī)且均勻分布的一個(gè)數(shù)。因此可將式（2）中φij(t)表示為直接由隨機(jī)數(shù)產(chǎn)生，可將式（2）表示為

在量子粒子群算法中采用該式表示歷史最優(yōu)位置與群體最優(yōu)位置之間的某個(gè)隨機(jī)多維點(diǎn)的位置。

由上述分析并依據(jù)文獻(xiàn)［14］提出在pij點(diǎn)建立δ勢阱引力場。由于粒子的位置和速度在量子化空間中無法同時(shí)確定，所以需要利用波函數(shù)ψ(x,t)（x為粒子的位置，t為粒子的速度）來描述粒子的運(yùn)動狀態(tài)，粒子的空間三維坐標(biāo)向量為X=(x,y,z)。波函數(shù)表述的物理意義為粒子在解空間某點(diǎn)出現(xiàn)的概率等于波函數(shù)模的平方，可表示為

式中，Q(X,t)為概率密度分布函數(shù)。

易知，式（4）滿足歸一化條件

式（5）即為Schrodinger公式，是粒子在量子空間中的動力學(xué)方程［11］。

為求得多維δ勢阱下粒子的隨機(jī)位置方程，先在一維空間中建立一個(gè)一維勢阱，表達(dá)式為

式中，x,y為一維條件下點(diǎn)的坐標(biāo)；p為δ勢阱引力場一維坐標(biāo)；γ為常數(shù)。

令y=x-p，Schrodinger公式可改寫為

Schrodinger公式在一維條件下的解為

式中，L為勢阱的特征長度；β為勢阱的特征長度的倒數(shù)。

通過式（6）、（7）、（8）求得粒子出現(xiàn)在解空間的概率分布函數(shù)。為求得粒子在每點(diǎn)的適應(yīng)度值，需求出粒子在解空間的確切位置x，進(jìn)而得到粒子的位置更新公式。

1.1.2 QPSO算法粒子位置進(jìn)化方程

在上文已求得一維條件下的波函數(shù)，即得到了粒子相對于pij點(diǎn)的位置。為求得粒子的確切位置，本研究采用蒙特卡羅隨機(jī)模型（逆變換法）方法［11］推導(dǎo)出一維δ勢阱下粒子的隨機(jī)位置方程：

式中，u為區(qū)間(0,1)上的隨機(jī)數(shù)；p為δ勢阱引力場一維坐標(biāo)。

將一維拓展為多維，即在每一維上都建立一個(gè)一維δ勢阱，可得粒子i的隨機(jī)位置方程（即位置進(jìn)化方程）：

式中，θ1為區(qū)間 (0,1)上的隨機(jī)數(shù)；t為迭代次數(shù)；xij(t+1)為第t+1次迭代中粒子i的位置；PGij(t)為第t次迭代中粒子的歷史最優(yōu)位置pij與群體最優(yōu)位置pg之間的某個(gè)隨機(jī)多維點(diǎn)的位置；Lij(t)為多維條件下的勢阱的特征長度。

對函數(shù)Lij(t)進(jìn)行控制時(shí)的運(yùn)算公式為

式中，D為粒子維數(shù)；mbest(t)為所有粒子個(gè)體的歷史最優(yōu)位置的均值；N為粒子總數(shù)；α為收縮因子是量子粒子群中唯一需要確定的參數(shù)，是用來控制量子粒子群的收斂速度，即擴(kuò)張收縮速度［16-17］。

式（11）中參數(shù)α控制公式為

式中，T為算法的總迭代次數(shù)［4］；t為算法的當(dāng)前迭代次數(shù)；αmin與αmax為常數(shù)，通常αmax取1、αmin取0.5。本研究αmax取1，αmin取0.5。

綜上分析，可得粒子的進(jìn)化公式為

通過上述分析可知：QPSO引入平均個(gè)體歷史最優(yōu)位置，取消了速度向量，將粒子搜索區(qū)域擴(kuò)大為全局唯一的解空間；在平均位置的作用下每個(gè)粒子在收斂過程中不得不考慮其他粒子的運(yùn)動狀態(tài)而獨(dú)立地向最優(yōu)點(diǎn)gbest靠近，遠(yuǎn)離gbest的粒子會將平均位置拉向自己，此時(shí)聚集于gbest附近的粒子大范圍搜索最終的群體最優(yōu)位置，等遠(yuǎn)離的點(diǎn)靠近gbest時(shí)再共同縮小搜索區(qū)域，故而加強(qiáng)了粒子的全局搜索能力，降低了粒子陷入局部最優(yōu)的概率［15-21］。

1.2 基于QPSO算法的概率積分法參數(shù)反演模型構(gòu)建

根據(jù)概率積分法原理［1-5］，可將地表任意一點(diǎn)的下沉值與水平移動表示為

式中，P為概率積分法參數(shù)矩陣，P=[q,b,tanβ,θ,S1,S2,S3,S4]；q為下沉系數(shù)；b為水平移動系數(shù)；tanβ為主要影響角正切值；θ為開采影響傳播角；S1，S2為上下拐點(diǎn)偏移距；S3，S4為左右拐點(diǎn)偏移距；（x，y）為觀測站點(diǎn)的坐標(biāo)

假設(shè)地表任意一點(diǎn)M的實(shí)測下沉值與水平移動值分別為WM實(shí)，UM實(shí)，該點(diǎn)通過QPSO算法得出的下沉與水平移動預(yù)計(jì)值分別為WM實(shí)，UM實(shí).則該點(diǎn)的下沉與水平移動殘差vM為

式中，abs為絕對值運(yùn)算方式；sum為求和運(yùn)算方式。

根據(jù)誤差平方和最小原則，本研究構(gòu)建的概率積分法求參準(zhǔn)則為

式中，m為測點(diǎn)數(shù)，i＜m?？筛鶕?jù)該式采用QPSO算法求解概率積分法參數(shù)。

綜上分析，本研究采用MATLAB語言對QPSO算法進(jìn)行編程實(shí)現(xiàn)。假設(shè)在解空間有100個(gè)粒子，每個(gè)粒子有8維，在第t次迭代中粒子的當(dāng)前位置可表示為xi8=[xi1(t),xi2(t),…,xi8(t)],i=1,2,…,100；粒子的歷史最優(yōu)位置可表示為pi8(t)=[pi1(t),pi2(t),…,pi8(t)]，群體最優(yōu)位置可表示為pg8(t)=[pg1(t),pg2(t),…,pg8(t)]。將下沉和移動變形實(shí)測值與預(yù)計(jì)值之差的絕對值累加值作為適應(yīng)度函數(shù)f：

適應(yīng)度函數(shù)越小說明粒子位置越好，即參數(shù)反演結(jié)果越好。

本研究基于QPSO算法的概率積分法參數(shù)反演流程如圖1所示。

算法實(shí)現(xiàn)步驟為：

（1）初始化種群，將粒子當(dāng)前位置初始化為個(gè)體歷史最優(yōu)位置，計(jì)算適應(yīng)度找到群體最優(yōu)位置。

（2）通過式（10）得到粒子i(1≤i≤100)的介于個(gè)體歷史最優(yōu)與群體最優(yōu)之間的位置PGi。

（3）依據(jù)式（11）計(jì)算mbest，即個(gè)體歷史最優(yōu)位置的均值。

（4）根據(jù)式（13）更新粒子的位置。

（5）計(jì)算當(dāng)前迭代次數(shù)下的粒子適應(yīng)度值，并與前一次迭代比較，如果適應(yīng)度值小，將粒子歷史最優(yōu)位置更換為當(dāng)前粒子的位置；否則，不變。找到群體最優(yōu)位置與前一次迭代結(jié)果比較，適應(yīng)度值小，則進(jìn)行替換；否則，不變。

（6）重復(fù)步驟（1）至步驟（5），若循環(huán)達(dá)到最大迭代次數(shù)或滿足精度，則跳出循環(huán)，輸出最終反演所得的8個(gè)概率積分法參數(shù)。

2 模擬試驗(yàn)

2.1 工作面概況

本研究所模擬的工作面煤層采厚為3.0 m，煤層傾角為5°，走向線長800 m，傾向線長500 m，平均采深為400 m，采用全部垮落法頂板管理，達(dá)到充分開采。地表沉陷預(yù)計(jì)的概率積分法參數(shù)設(shè)計(jì)為：下沉系數(shù)q為0.8；主要影響正切tanβ為2.5；水平移動系數(shù)b為0.25，開采影響傳播角θ為85°；拐點(diǎn)偏距（S1=S2=S3=S4）為0.15倍的平均采深，即60 m。在試驗(yàn)中，在移動盆地內(nèi)設(shè)計(jì)了沿走向和傾向2條主斷面的移動和變形觀測線，每個(gè)監(jiān)測點(diǎn)間距為30 m，其中走向線（a線）長度為800 m，共布設(shè)了45個(gè)監(jiān)測點(diǎn)；傾向線（b線）長度為500 m，共布設(shè)了35個(gè)監(jiān)測點(diǎn)。模擬工作上的測線及測點(diǎn)分布如圖2所示。

2.2 QPSO及PSO算法參數(shù)反演準(zhǔn)確性及效率比較

本研究使用QPSO及PSO兩種智能優(yōu)化算法對概率積分法參數(shù)進(jìn)行反演，將下沉與水平移動殘差絕對值之和作為適應(yīng)度函數(shù)；再將參數(shù)設(shè)計(jì)值與參數(shù)反演值進(jìn)行比較，通過比較求取參數(shù)的相對誤差和中誤差的大小反映參數(shù)反演的準(zhǔn)確性，并進(jìn)一步分析比較兩種算法的運(yùn)行時(shí)間。為了避免求參偶然誤差，在相同條件下分別進(jìn)行了10次試驗(yàn)，結(jié)果如表1所示。

注：PSO與QPSO算法參數(shù)取值是通過10次試驗(yàn)取均值得到；10次試驗(yàn)中，PSO與QPSO算法運(yùn)行耗時(shí)分別為1 665.43 s和197.79 s。

由表1可知：①Q(mào)PSO算法反演的概率積分法參數(shù)中除了θ值中誤差略大外，其余參數(shù)的中誤差均小于PSO算法反演的參數(shù)中誤差，表明QPSO算法的穩(wěn)定性優(yōu)于PSO算法；②QPSO算法反演的概率積分法參數(shù)除了部分拐點(diǎn)偏距（S1、S2）的參數(shù)相對誤差略大外，其余參數(shù)的相對誤差都小于PSO算法，反映出QPSO算法反演精度優(yōu)于PSO算法；此外，QPSO算法耗時(shí)也遠(yuǎn)小于PSO算法，表明QPSO算法的運(yùn)行效率優(yōu)于PSO算法。

為了進(jìn)一步分析兩種算法的參數(shù)反演效果，根據(jù)QPSO與PSO算法的概率積分法參數(shù)反演結(jié)果，繪制了兩者的下沉曲線及水平移動曲線，如圖3至圖6所示。

由圖3至圖6分析可知：PSO及QPSO算法的下沉及水平移動擬合效果均較好；QPSO算法的絕對誤差波動小于PSO算法。通過模擬試驗(yàn)數(shù)據(jù)計(jì)算得PSO的參數(shù)取均值后求得的下沉值與水平移動值的擬合中誤差為8.50 mm；QPSO的參數(shù)取均值后求得的下沉值與水平移動值擬合結(jié)果中誤差為3.19 mm。可見，QPSO算法概率積分法參數(shù)反演精度優(yōu)于PSO算法，并且QPSO算法的運(yùn)行效率遠(yuǎn)高于PSO算法。

2.3 QPSO算法可靠性研究

對表1進(jìn)一步分析可知：①利用QPSO算法進(jìn)行10次概率積分法參數(shù)反演后，參數(shù)的最大相對誤差不超過±4%，說明采用QPSO算法反演概率積分法參數(shù)精度較高；②與PSO算法相比，QPSO算法反演的概率積分法參數(shù)除了影響傳播角θ擬合中誤差略大外，其余參數(shù)的擬合中誤差均小于PSO算法，說明QPSO算法穩(wěn)定性優(yōu)于PSO算法；③QPSO算法運(yùn)行效率也明顯優(yōu)于PSO算法。模擬試驗(yàn)結(jié)果反映出QPSO算法可靠性較好。

為了進(jìn)一步分析QPSO算法的穩(wěn)定性，本研究對PSO與QPSO算法反演的參數(shù)分別進(jìn)行了波動性分析，結(jié)果如圖7所示。

由圖7可知：QPSO算法參數(shù)反演結(jié)果中除了影響傳播角θ波動略大外，其余參數(shù)的波動性均小于PSO算法反演的各個(gè)參數(shù)。

3 工程實(shí)例

3.1 礦區(qū)概況

淮南顧橋礦南二采區(qū)1414（1）是該礦南區(qū)的首采工作面，該工作面采用后退式開采，機(jī)械化掘進(jìn)，一次采全高，全部垮落法管理頂板。該工作面沿煤層走向布置，工作面開采尺寸為走向長度×傾向長度為2 120 m×251 m，工作面走向方向?yàn)槌浞植蓜樱瑑A向方向?yàn)榉浅浞植蓜?，總體為非充分采動［22］。工作面平均采高為3.0 m，煤層傾角平均為5°，為近水平煤層。工作面平均深度為735 m，傾向觀測線布置在距離切眼和停采線1 144 m和976 m處，共布設(shè)了3個(gè)控制點(diǎn)和50個(gè)監(jiān)測點(diǎn)，點(diǎn)間距為30 m，傾向線長度為1 500 m。

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

選取顧橋南礦1414（1）工作面122個(gè)監(jiān)測點(diǎn)（走向線上75個(gè)監(jiān)測點(diǎn)，傾向線上47個(gè)監(jiān)測點(diǎn)）的實(shí)測數(shù)據(jù)為基礎(chǔ)，分別采用基于QPSO、PSO算法構(gòu)建參數(shù)反演模型對該工作面進(jìn)行概率積分法參數(shù)反演，并對結(jié)果進(jìn)行對比分析。兩種模型分別進(jìn)行了10次試驗(yàn)，將試驗(yàn)結(jié)果取平均值，并計(jì)算參數(shù)的擬合中誤差，結(jié)果如表2所示，擬合結(jié)果分別如圖8至圖11所示。

由表2以及圖8至圖11分析可知：QPSO算法反演參數(shù)得到的擬合中誤差均優(yōu)于PSO算法，可認(rèn)為QPSO算法穩(wěn)定性較好，并且QPSO算法效率明顯優(yōu)于PSO算法；PSO與QPSO算法的擬合結(jié)果都較為精確，將兩種智能優(yōu)化算法反演的參數(shù)取平均值，計(jì)算出PSO算法下沉值、水平移動值與實(shí)測下沉值、水平移動值擬合中誤差為70.93 mm，QPSO算法下沉值、水平移動值與實(shí)測下沉值、水平移動值擬合中誤差為72.04 mm，兩種算法的下沉值與水平移動值擬合精度相當(dāng)。根據(jù)文獻(xiàn)［9］分析可知，盡管二者的擬合精度都符合工程應(yīng)用標(biāo)準(zhǔn)，但QPSO算法運(yùn)算效率相對于PSO算法有明顯優(yōu)勢，故QPSO算法對于實(shí)現(xiàn)概率積分法參數(shù)高效率、高精度反演有較好的適用性。

注：概率積分法參數(shù)取值區(qū)間是根據(jù)1414（1）工作面已有資料得出；10次試驗(yàn)中，PSO與QPSO算法運(yùn)行耗時(shí)分別為1639.43 s和197.72 s。

4 結(jié)語

為實(shí)現(xiàn)對概率積分法預(yù)計(jì)參數(shù)的精確反演，提出了基于QPSO算法的參數(shù)反演模型。模擬試驗(yàn)以及顧橋南礦1414（1）工作面實(shí)例分析表明：QPSO算法模型反演參數(shù)的波動性略小于PSO算法模型，并且QPSO算法模型的運(yùn)算效率明顯優(yōu)于PSO算法模型，對于實(shí)現(xiàn)概率積分法開采沉陷預(yù)計(jì)參數(shù)高效、精準(zhǔn)反演有一定的參考價(jià)值。