量子粒子群算法反演概率積分法參數

朱尚軍 王 磊， 魏 濤 蔣 創 江克貴 查劍鋒 孔 川4

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

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

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

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

1.1 理論基礎

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

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

1.1.1 勢場模型建立

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

由式（1）和式（2）可知，當c1=c2時，φij(t)即為在區間(0,1)上隨機且均勻分布的一個數。因此可將式（2）中φij(t)表示為直接由隨機數產生，可將式（2）表示為

在量子粒子群算法中采用該式表示歷史最優位置與群體最優位置之間的某個隨機多維點的位置。

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

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

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

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

為求得多維δ勢阱下粒子的隨機位置方程，先在一維空間中建立一個一維勢阱，表達式為

式中，x,y為一維條件下點的坐標；p為δ勢阱引力場一維坐標；γ為常數。

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

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

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

通過式（6）、（7）、（8）求得粒子出現在解空間的概率分布函數。為求得粒子在每點的適應度值，需求出粒子在解空間的確切位置x，進而得到粒子的位置更新公式。

1.1.2 QPSO算法粒子位置進化方程

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

式中，u為區間(0,1)上的隨機數；p為δ勢阱引力場一維坐標。

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

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

對函數Lij(t)進行控制時的運算公式為

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

式（11）中參數α控制公式為

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

綜上分析，可得粒子的進化公式為

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

1.2 基于QPSO算法的概率積分法參數反演模型構建

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

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

假設地表任意一點M的實測下沉值與水平移動值分別為WM實，UM實，該點通過QPSO算法得出的下沉與水平移動預計值分別為WM實，UM實.則該點的下沉與水平移動殘差vM為

式中，abs為絕對值運算方式；sum為求和運算方式。

根據誤差平方和最小原則，本研究構建的概率積分法求參準則為

式中，m為測點數，i＜m。可根據該式采用QPSO算法求解概率積分法參數。

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

適應度函數越小說明粒子位置越好，即參數反演結果越好。

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

算法實現步驟為：

（1）初始化種群，將粒子當前位置初始化為個體歷史最優位置，計算適應度找到群體最優位置。

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

（3）依據式（11）計算mbest，即個體歷史最優位置的均值。

（4）根據式（13）更新粒子的位置。

（5）計算當前迭代次數下的粒子適應度值，并與前一次迭代比較，如果適應度值小，將粒子歷史最優位置更換為當前粒子的位置；否則，不變。找到群體最優位置與前一次迭代結果比較，適應度值小，則進行替換；否則，不變。

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

2 模擬試驗

2.1 工作面概況

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

2.2 QPSO及PSO算法參數反演準確性及效率比較

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

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

由表1可知：①QPSO算法反演的概率積分法參數中除了θ值中誤差略大外，其余參數的中誤差均小于PSO算法反演的參數中誤差，表明QPSO算法的穩定性優于PSO算法；②QPSO算法反演的概率積分法參數除了部分拐點偏距（S1、S2）的參數相對誤差略大外，其余參數的相對誤差都小于PSO算法，反映出QPSO算法反演精度優于PSO算法；此外，QPSO算法耗時也遠小于PSO算法，表明QPSO算法的運行效率優于PSO算法。

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

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

2.3 QPSO算法可靠性研究

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

為了進一步分析QPSO算法的穩定性，本研究對PSO與QPSO算法反演的參數分別進行了波動性分析，結果如圖7所示。

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

3 工程實例

3.1 礦區概況

淮南顧橋礦南二采區1414（1）是該礦南區的首采工作面，該工作面采用后退式開采，機械化掘進，一次采全高，全部垮落法管理頂板。該工作面沿煤層走向布置，工作面開采尺寸為走向長度×傾向長度為2 120 m×251 m，工作面走向方向為充分采動，傾向方向為非充分采動，總體為非充分采動［22］。工作面平均采高為3.0 m，煤層傾角平均為5°，為近水平煤層。工作面平均深度為735 m，傾向觀測線布置在距離切眼和停采線1 144 m和976 m處，共布設了3個控制點和50個監測點，點間距為30 m，傾向線長度為1 500 m。

3.2 試驗結果及分析

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

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

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

4 結語

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