基于改進SPSO 算法的開采沉陷預測參數反演研究

徐可心，何 榮，白偉森

（河南理工大學 測繪與國土信息工程學院，河南 焦作 454000）

煤礦地下開采造成地表移動形變，引起地表裂縫、建筑物損害等，部分地區出現積水、水土流失等地質災害，嚴重破壞礦區的生態環境。因此，有必要通過對礦區地表形變監測，獲取變形數據、求算相關參數，總結地表移動盆地的變形規律，為維護礦區經濟可持續發展及生態文明建設提供基礎數據支撐[1]。

概率積分法是我國目前開采沉陷領域科學研究常用方法之一[2]，構建了地表變形與地下開采空間的關系模型。不同的開采方法和地質采礦條件，預計參數具有較大的差異性，因此，科學確定預計參數大小對地表變形量計算至關重要。其參數求算主要有特征點求參、曲線擬合和曲面擬合等方法，曲面擬合法將曲線擬合法求參的原理推廣到整個下沉盆地，充分利用非主斷面地表變形監測數據，得到更加精確的地表移動變形參數值[3]。李春意[4]等將曲面擬合法用于求取矩形工作面的概率積分法參數，并驗證了非線性曲面擬合法用于開采沉陷預計的實用性。該方法不要求觀測站位于下沉盆地的主斷面上，是非主斷面設站求取概率積分法參數的主要方法。但其對初值選取要求較高，工作量大、求參速度慢。隨著計算機技術的發展，為提高參數反演的效率，不少學者將智能算法應用到了概率積分法反演參數的研究中，優化算法有模矢法[5]、遺傳算法[6]、果蠅算法[7]、神經網絡算法[8]以及粒子群算法[9]等。與其他的尋優算法相比, 粒子群算法擁有較好的收斂能力,計算速度更快[10]。徐孟強[11]等首次將粒子群算法引入了概率積分法參數反演，認為粒子群算法可以有效解決隨機誤差和測點缺失的問題，存在全局搜索能力差、容易提前收斂等問題；葉偉[12]等將自適應魚群算法用于曲面擬合求參，得到了較為可靠的擬合下沉曲面，所用方法運算復雜、計算效率較低。

基于此，針對傳統粒子群算法存在的局限性，使用改進標準粒子群算法，進行曲面擬合概率積分法參數反演，計算參數擬合中誤差和下沉擬合誤差，驗證改進標準粒子群算法用于概率積分法參數反演的有效性。

1 基于非主斷面設站的曲面擬合模型

相比主斷面設站法，非主斷面設站具有布設靈活的特點，既方便觀測，又減少了對農田的影響，被廣泛用于礦區的地表變形監測。觀測站位置的選取對參數反演精度有較大影響，為保證計算精度，觀測站應盡量布設在變形值較大位置，或下沉盆地的特征點處[13]。

1.1 工程實例概況

試驗區位于河北省邯鄲市某礦，區內地形為低丘陵，地勢南高北低，有較淺沖溝，部分有基巖出露。2516 工作面位于205 采區運輸下山北翼，北面為充填集中軌道巷，西面為已開采的其他工作面，東南方向為村莊。工作面開采條件見表1。

表1 工作面開采條件Table 1 Mining conditions of working face

為監測礦區的地表移動變形，2516 工作面上方布設了走向與傾向2 條觀測線，受地表村莊分布和地形限制，以非主斷面設站法布設。觀測線布設和井上下對照圖如圖1，走向觀測線長671 m，共布設17 個水準觀測點，傾向觀測線長1 167 m，共布設28 個水準觀測點。2516 工作面共進行5 次觀測，其中走向最大下沉點下沉值為572 mm，傾向最大下沉點為走向與傾向觀測線交點，下沉值為556 mm。

1.2 曲面擬合法函數模型

曲面擬合法不受觀測站布設方法限制，但僅適用于矩形工作面。由圖1 可知，2516 工作面為非主斷面設站的矩形工作面，可使用曲面擬合法反演概率積分法參數。由概率積分法可知，矩形工作面下沉盆地內曲面擬合任意點下沉的函數W（x，y）為：

圖1 地表移動觀測站布設與井上下對照圖Fig.1 Layout of surface movement observation stations and comparison between upper and lower wells

式中：W0為實測最大下沉值，mm；H 為走向主斷面平均采深，m；H1、H2分別為下山和上山方向的采深，m；β、β1、β2為主要影響角，（°）；D1、D3為矩形工作面傾向和走向主斷面寬度，m；S1、S2分別為上山和下山方向拐點偏移距，m；S3、S4分別為走向左和走向右拐點偏移距，m；θ0為開采影響傳播角，（°）；α 為煤層傾角，（°）。

2 基于改進SPSO 算法的概率積分法參數反演

2.1 改進SPSO 算法

粒子群優化算法（Particle Swarm Optimization，PSO）是1 種群智能優化算法，源于對鳥群社會行為的研究。Reynolds[14]于1986 年設計了1 個再現鳥群聚集行為的人工生命系統，隨后Eberhart 和Kennedy[15]基于該系統提出了粒子群算法。Y. Shi[16]等學者引入了慣性權重ω，并修正了粒子速度公式，引入ω 的粒子群算法被稱為標準粒子群算法（Standard Particle Swarm Optimization，SPSO）。

假設在D 維的搜索空間中，N 個粒子組成1 個群體，第i 個粒子的位置表示為xi＝［xi1，xi2，…，xiD］，其飛行速度表示為vi＝［vi1，vi2，…，viD］，當前搜索到的最優位置為pi＝［pi1，pi2，…，piD］，整個群體當前搜索到的最優位置為pg＝［pi1，pi2，…，piD］。由此得到粒子速度與位置的更新公式：

慣性權重ω 是SPSO 算法的重要參數，用于控制算法開發和搜索能力。當ω 取值較大，算法全局尋優能力較強，局部尋優能力較弱；當ω 取值較小，正好相反。在算法搜索過程中將ω 設置為動態，可以平衡收斂速度和搜索能力，使算法優化結果更精確。

PSO 算法的優勢在于收斂能力較好，但面對復雜多維問題，PSO 算法常出現早熟，即提前收斂陷入局部最優解[17]。針對PSO 算法容易早熟的問題，在SPSO 算法的速度更新公式中引入漸進收斂的方法，改善算法全局收斂性，提高算法搜索能力。改進SPSO 算法粒子速度迭代公式為：

2.2 基于改進SPSO 算法的概率積分法求參模型

由式（1），可將下沉盆地內地表任意點的下沉值WPos表示為：

式中：Pos 為概率積分法的參數矩陣，Pos＝［q，tanβ，tanβ1，tanβ2，S1，S2，S3，S4］；q 為下沉系數；（x，y）為觀測點坐標。

基于改進SPSO 算法的概率積分法參數反演流程如圖2。

圖2 SPSO 算法概率積分法參數反演流程Fig.2 SPSO algorithm probability integration method parameter inversion process

在改進SPSO 算法中，應先設置算法參數，確定適應度函數，然后通過迭代得到最優解，反演流程的具體步驟如下。

1）初始化種群。設置最大迭代次數、自變量個數、粒子的最大速度以及粒子約束條件為搜索空間，在速度空間和搜索空間上初始化粒子飛行速度和位置，設置種群規模，生成初始群體。使用MATLAB 語言對求參模型進行編程，假設空間內共有100 個粒子，學習因子c1=c2=2，nmax=500。粒子維數由所求公式參數數量決定，概率積分法參數矩陣中有8 個參數，分別為：下沉系數q，主要影響角正切tanβ、tanβ1、tanβ2以及拐點偏移距Si（i＝1，2，3，4），確定粒子維數為8。為平衡算法搜索能力，設置慣性權重隨迭代次數增加線性遞減，令ωmax為0.9，ωmin為0.4。鑒于工作面實際情況，結合礦區現有巖移參數，確定粒子約束條件，即參數取值區間，參數取值上下限見表2。

表2 參數取值上下限Table 2 Upper and lower limits of parameter values

2）定義適應度函數。粒子每更新1 次位置就計算1 次適應度值，適應度值反應粒子的優劣程度。根據誤差平方和最小原則，預計值和實測值之差的平方和越小，參數反演精度越高。構建適應度函數為：

式中：RSS 為擬合值與實測值之差的平方和；m為觀測點總個數；Wi為實測下沉值，mm；wi為擬合下沉值，mm。

3）根據式（3）、式（5）更新粒子的速度和位置，并處理超出約束條件的粒子。迭代過程中，粒子通過個體極值和群體極值更新速度和位置。個體極值為每個粒子找到最優解，從個體最優解中找到全局最優解，與歷史全局最優解進行比較，不斷更新得到最終結果。

4）當相鄰2 代全局最優解之間的偏差最小，或循環達到最大迭代次數時，循環終止，輸出反演參數。在本研究構建的模型中，當相鄰2 代全局最優解適應度值差值為0 時，認為找到最優解，或迭代次數達到500 次時循環停止。一般情況下迭代次數越多精度越高，但多余迭代次數會降低運算效率，為避免無效迭代并保證運算精度，通過試驗確定最大迭代次數為500。循環停止輸出最優解，即反演所得概率積分法參數。

3 基于改進SPSO 算法的模型在工程中的應用

3.1 參數反演結果

以2516 工作面實測下沉數據為基礎（圖1），分別使用改進SPSO 算法和PSO 算法進行參數反演，對比其求參結果和擬合中誤差。為避免試驗偶然性，2 種算法以相同設置各進行10 次試驗，取試驗結果的平均值為參數值，所得結果和礦區經驗值見表3。礦區經驗值參考205 采區已有地表移動變形參數，結合2516 工作面開采技術確定。

表3 反演參數對比Table 3 Comparison of inversion parameters

由表3 可知，除主要影響角正切tan β 外，改進SPSO 算法反演參數擬合中誤差均小于PSO 算法反演參數擬合中誤差，即改進SPSO 算法反演參數波動性小于PSO 算法，證明改進SPSO 算法穩定性優于PSO 算法。

3.2 可靠性分析及收斂速率比較

改進SPSO 算法和PSO 算法的迭代次數對比如圖3，在算法10 次運算過程中，改進SPSO 算法出現1 次收斂過快未求得最優解，PSO 算法出現3 次。且改進SPSO 算法的平均迭代次數為47.3 次，PSO 算法平均迭代次數為300.4 次。可以看出改進SPSO 算法收斂速度遠高于PSO 算法。

圖3 迭代次數對比Fig.3 Comparison of iteration times

為分析改進SPSO 算法反演概率積分法參數的效果，將改進SPSO 算法和PSO 算法反演參數值和礦區經驗值分別代入概率積分法公式得到擬合下沉曲面，與實測數據對比，并計算其擬合誤差。

下沉擬合誤差計算公式為：

由式（8）得，改進SPSO 算法反演參數擬合誤差σ=4.66%，PSO 算法擬合誤差σ=4.67%，礦區經驗值擬合誤差σ=10.2%。改進SPSO 算法反演參數精度和PSO 算法相當，高于礦區經驗值。改進SPSO 算法反演參數和礦區經驗值擬合下沉曲面如圖4。算法擬合下沉值及實測值對比如圖5。

由圖4 可知，改進SPSO 算法反演參數擬合下沉曲面更接近實際下沉盆地的形狀，擬合效果較好。由圖5 可知，工作面走向方向擬合效果較好，傾向方向的部分監測點實測值大于擬合值。由圖1 可知，2516 工作面西側為2 個已開采工作面，且與2516工作面相距較近，2516 工作面的開采勢必會引起老采空區的活化。根據覆巖“活化”機理[18]，老采空區會因為上覆巖層的“活化”而產生二次破裂，造成地表下沉值大于預計值。試驗反演概率積分法參數僅考慮單一工作面開采條件，2516 工作面在實際開采過程中受到了2314 和2509 已采工作面的影響，使得傾向觀測線上山方向部分監測點擬合效果不理想。

圖4 下沉擬合曲面及實測值對比Fig.4 Subsidence fitting surface and comparison of measured values

圖5 算法擬合下沉值及實測值對比Fig.5 Algorithm fitting sinking values and comparison of measured values

為進一步分析改進SPSO 算法反演概率積分法參數的可靠性，基于改進SPSO 算法的反演參數值和概率積分法原理，對2516 工作面地表移動觀測站的傾斜和曲率變形進行了預測。算法預計值與實測值對比如圖6。

圖6 算法預計值與實測值對比Fig.6 Comparison between the predicted values of the algorithm and the measured values

由圖6 可知，算法傾斜及曲率變形預計值和實測值相比，整體變化趨勢基本一致，除少部分異常點外，整體誤差較低，預計效果較好。證明改進SPSO算法反演參數具有較高可靠性。

4 結 語

1）針對非主斷面設站反演概率積分法參數的問題，將漸進收斂引入SPSO 算法，建立了基于改進SPSO 算法的概率積分法參數反演模型。研究表明，該方法參數反演運算效率高、擬合效果較好，通過該方法反演的參數具有較高可靠性。

2）反演參數值擬合中誤差和下沉擬合誤差的對比分析表明，改進SPSO 算法反演參數擬合中誤差整體優于PSO 算法，參數反演精度高于礦區經驗值，有效提高了參數反演準確性。和PSO 算法相比，改進SPSO 算法具有更好的收斂能力和穩定性。

3）由下沉、傾斜、曲率預計值與實測值的對比分析可知，改進SPSO 算法反演參數預計的擬合下沉盆地比較符合實際情況，對于提高礦區開采沉陷預計精度有一定的參考價值。