代 陽

(1.安徽理工大學(xué) 測(cè)繪學(xué)院，安徽 淮南 232001；2.礦山采動(dòng)災(zāi)害空天地協(xié)同監(jiān)測(cè)與預(yù)警安徽省教育廳重點(diǎn)實(shí)驗(yàn)室，安徽 淮南 232001)

煤炭在我國能源供給結(jié)構(gòu)中占有主體地位。在地下開采時(shí)，礦區(qū)巖體受力平衡遭到破壞，引發(fā)上覆巖層移動(dòng)變形，導(dǎo)致地表沉降移動(dòng)。且隨著開采范圍的擴(kuò)大，周邊建筑與環(huán)境勢(shì)必受到影響，嚴(yán)重時(shí)可能會(huì)引發(fā)建筑物塌陷、地下水資源被破壞、地表塌陷等地質(zhì)災(zāi)害和次生災(zāi)害，嚴(yán)重影響人們的生活安全[1-2]。因此，研究礦區(qū)開采引發(fā)的地表移動(dòng)變形規(guī)律，在災(zāi)害治理與防范中有著重要的意義。

概率積分法是預(yù)計(jì)礦區(qū)開采沉陷變化規(guī)律的積分模型，其模型的主要思想是基于LITWINISZYN提出的隨機(jī)介質(zhì)理論，利用采空區(qū)實(shí)測(cè)數(shù)據(jù)擬合得到的預(yù)計(jì)參數(shù)，預(yù)測(cè)礦區(qū)地表移動(dòng)變化。構(gòu)建概率積分法模型的關(guān)鍵，在于保證輸入預(yù)計(jì)參數(shù)的準(zhǔn)確性。由于該模型非線性程度高、參數(shù)之間相關(guān)性強(qiáng)、求解易受迭代初值影響，采用傳統(tǒng)的求解參數(shù)方法往往會(huì)帶來局部解甚至無解[3]。為解決上述問題，學(xué)者們引用了優(yōu)化及智能優(yōu)化算法反演模型參數(shù)。葛家新等以皖北礦區(qū)為例，采用模矢法進(jìn)行反演參數(shù)，擬合預(yù)測(cè)的下沉效果較好，但易受到參數(shù)設(shè)計(jì)初值的影響[4]。查劍鋒等采用遺傳算法進(jìn)行模型參數(shù)估計(jì)，并與最小二乘法和模矢法相比，證明遺傳算法的準(zhǔn)確性與穩(wěn)定性更優(yōu)[5]。蘇軍明等采用模擬退火的方法，在搜索最優(yōu)參數(shù)的過程中加入隨機(jī)因素，更好地得到全局最優(yōu)解[6]。徐孟強(qiáng)等采用粒子群算法反演參數(shù)能夠得到準(zhǔn)確的預(yù)計(jì)參數(shù)，且在一定范圍內(nèi)受初值誤差影響較小[7]。智能優(yōu)化算法通過模仿自然現(xiàn)象及生物進(jìn)化規(guī)律，迭代逐步逼近參數(shù)最優(yōu)值，然而在搜索的過程中，容易早熟且每次反演得到的參數(shù)具有一定的隨機(jī)性，為構(gòu)建概率積分模型帶來影響[8]。因此，有必要引入新的算法反演概率積分法中的參數(shù)。

概率積分法參數(shù)反演可以看作是一種非線性多參數(shù)求解最優(yōu)值的過程，提出采用二次序列規(guī)劃(Sequential Quadratic Programming，SQP)算法反演參數(shù)，以達(dá)到擬合參數(shù)精度最優(yōu)的目的。SQP算法屬于最優(yōu)化算法，通過二次規(guī)劃確定每次迭代的下降方向，并通過價(jià)值函數(shù)獲取步長，多次迭代逼近最優(yōu)值[9]。在處理非線性問題中，SQP算法不依賴于初值，并具有良好的收斂與自我校正能力，結(jié)果的穩(wěn)定性與可靠性強(qiáng)。目前，該算法已經(jīng)被成功運(yùn)用于各種非線性最優(yōu)化問題中，例如在飛機(jī)數(shù)字化裝配[10]、陣列信號(hào)處理[11]、高密度電阻率等領(lǐng)域[12]。

根據(jù)概率積分模型構(gòu)建目標(biāo)函數(shù)，引用SQP算法進(jìn)行參數(shù)反演，驗(yàn)證SQP算法運(yùn)用在開采沉陷預(yù)測(cè)參數(shù)的可能性。通過模擬監(jiān)測(cè)數(shù)據(jù)受到不同程度的誤差及粗差的影響，進(jìn)行大量仿真實(shí)驗(yàn)，將SQP算法是反演結(jié)果與常用的優(yōu)化及智能優(yōu)化算法作比較，驗(yàn)證SQP算法抗隨機(jī)誤差及粗差的性能。

1 概率積分法預(yù)計(jì)模型

概率積分法是一種預(yù)計(jì)由于礦區(qū)開采引起的地表沉陷與水平移動(dòng)的模型[13](P116-120)。其中，預(yù)測(cè)地表任意點(diǎn)下沉值的函數(shù)模型為：

(1)

式中，(x,y)為地表任意點(diǎn)平面坐標(biāo)；W(x,y)為對(duì)應(yīng)坐標(biāo)下的地表下沉值；W0(x)為傾向方向充分采動(dòng)時(shí)，橫坐標(biāo)值為x時(shí)在走向面的下沉值；W0(y)為走向方向充分采動(dòng)時(shí)，縱坐標(biāo)為y時(shí)在傾向面的下沉值；W0為地表最大下沉值，限于篇幅，W0、W0計(jì)算公式見參考文獻(xiàn)[13]，不再累述。

預(yù)測(cè)地表任意點(diǎn)的水平移動(dòng)值的函數(shù)模型為：

(2)

式中，φ為由x軸正向逆時(shí)針旋轉(zhuǎn)到指定方向的方位角；U(x,y,φ)為對(duì)應(yīng)坐標(biāo)下的地表水平移動(dòng)值；U0(x)為走向主斷面上的水平移動(dòng)值；U0(y)為傾向主斷面上的水平移動(dòng)值。限于篇幅，U0詳細(xì)計(jì)算見參考文獻(xiàn)[13]。

概率積分法模型反演參數(shù)屬于求解多參數(shù)最優(yōu)化極值問題，目的是尋求一組變量，使預(yù)測(cè)的地表下沉值與水平移動(dòng)值與實(shí)測(cè)數(shù)據(jù)的差值平方和最小。其目標(biāo)函數(shù)的具體表現(xiàn)形式為：

(3)

式中，i為測(cè)站數(shù)據(jù)編號(hào)；N為測(cè)站的總數(shù)目；WS為測(cè)站實(shí)測(cè)下沉觀測(cè)值；US為測(cè)站實(shí)測(cè)水平移動(dòng)值。

2 SQP算法的理論基礎(chǔ)

SQP算法作為一種應(yīng)用于非線性函數(shù)參數(shù)反演的最優(yōu)化算法，應(yīng)用于概率積分法反演求參具體表現(xiàn)形式為[14]：

(4)

式中，x=[q,b,θ,tanβ,S1,S2,S3,S4]T為待求優(yōu)化參數(shù)，其中q為下沉系數(shù)，b為水平移動(dòng)系數(shù)，θ為開采影響傳播角，tanβ為主要影響角正切，S1,S2,S3,S4為對(duì)應(yīng)的下山方向、上山方向、左方向和右方向的拐點(diǎn)偏移距。f(x)對(duì)應(yīng)上文的公式(3)，Gj(x)對(duì)應(yīng)的限制條件，根據(jù)礦區(qū)情況和需求進(jìn)行改變。該算法采用拉格朗日求極值法將參數(shù)反演轉(zhuǎn)化為求解二次近似函數(shù)最優(yōu)值的問題：

(5)

通過線性化轉(zhuǎn)化后得到QP子問題，其目標(biāo)函數(shù)變換為：

(6)

xk+1=xk+akdk

(7)

式(6)中，d為全變量搜索方向，為梯度，k為迭代次數(shù)，Bk為拉格朗日函數(shù)的海森矩陣，通過半正定牛頓近似法(Dense quasi-newton approximation，BFGS)進(jìn)行計(jì)算得到[15]。式(7)中，ak為第k次迭代時(shí)的步長。

SQP算法反演概率積分模型預(yù)計(jì)參數(shù)的主要步驟為：

1)給定SQP算法反演參數(shù)x的初值、為算法反演提供初始解，同時(shí)令迭代次數(shù)k=0；

2)將反演參數(shù)值帶入，計(jì)算目標(biāo)函數(shù)f(x)及約束函數(shù)Gj(x)，差分計(jì)算梯度向量值，求解二次規(guī)劃子問題，計(jì)算搜索方向dk；

3)如果dk=0或者xk已經(jīng)符合精度檢驗(yàn)要求，則迭代終止，得到的xk即為該非線性函數(shù)的最優(yōu)解；否則，從點(diǎn)xk沿著dk方向搜索下一步迭代值，并確定步長ak，令xk+1=xk+akdk；

4)通過BFGS算法確定海森矩陣Bk，并轉(zhuǎn)化為正定矩陣。令k=k+1，轉(zhuǎn)至步驟(2)。

3 仿真實(shí)驗(yàn)與結(jié)果統(tǒng)計(jì)分析

3.1 設(shè)計(jì)工作面概況

依據(jù)工作面觀測(cè)點(diǎn)的設(shè)計(jì)原則，參考于文獻(xiàn)[8]中的工作面，制定模擬礦區(qū)的觀測(cè)線。具體設(shè)計(jì)為，煤層開采厚度m為4 m、煤層傾角α為3°、工作面的走向長L和傾向長D分別為670 m和300 m、開采深度為260 m，管理頂板為全部垮落法進(jìn)行管理[8]。概率積分參數(shù)的真值分別設(shè)計(jì)為：下沉系數(shù)q為0.64，水平移動(dòng)系數(shù)b為0.3，開采影響角正切值tanβ為2.37，開采影響傳播角θ為88.5°，上山方向、下山方向、左方向與右方向的拐點(diǎn)偏移距S1、S2、S3、S4分別設(shè)計(jì)為52 m。并在礦區(qū)上方沿著走向線布設(shè)25個(gè)監(jiān)測(cè)點(diǎn)(S1～S25)，沿著傾向線布設(shè)44個(gè)監(jiān)測(cè)點(diǎn)(E1～E44)，各個(gè)監(jiān)測(cè)點(diǎn)設(shè)計(jì)的點(diǎn)位布設(shè)，如圖1所示。

圖1 工作面及對(duì)應(yīng)點(diǎn)位模擬圖

3.2 SQP算法參數(shù)反演準(zhǔn)確性分析

以上述S1～S25、E1～E44共69個(gè)監(jiān)測(cè)點(diǎn)的下沉與水平移動(dòng)值作為反演數(shù)據(jù)的原始數(shù)據(jù)，結(jié)合礦區(qū)的基本信息，本文設(shè)計(jì)了五組不同的參數(shù)初始值，將反演結(jié)果與真值對(duì)比，對(duì)算法的穩(wěn)定性進(jìn)行評(píng)估。同時(shí)，通過比較不同初始參數(shù)設(shè)置下反演結(jié)果的差異，驗(yàn)證SQP算法預(yù)計(jì)參數(shù)的結(jié)果是否受到初始設(shè)計(jì)值的影響。

表1 SQP算法反演概率積分法準(zhǔn)確性分析

通過分析表1中的反演結(jié)果可以得到：1)對(duì)比五組不同參數(shù)設(shè)置下的SQP算法反演結(jié)果，隨著設(shè)置的初始值不斷增大，q、b、tanβ、θ與四個(gè)方向的拐點(diǎn)偏移距的預(yù)計(jì)參數(shù)并沒有發(fā)生明顯的變化。2)將反演結(jié)果與真值對(duì)比，參數(shù)q、b、tanβ和θ與四個(gè)方向的拐點(diǎn)偏移距的預(yù)計(jì)結(jié)果均與真值基本保持一致，相對(duì)誤差不到1%。根據(jù)上述實(shí)驗(yàn)結(jié)果，驗(yàn)證了SQP算法受到初始值設(shè)置的影響很小，且可靠性與穩(wěn)定性高。

3.3 SQP算法對(duì)隨機(jī)誤差、粗差的抗干擾能力

在礦區(qū)觀測(cè)中，觀測(cè)數(shù)據(jù)必然會(huì)受到隨機(jī)噪聲和粗差的影響，因此測(cè)試算法抗隨機(jī)誤差與粗差的干擾能力，是評(píng)價(jià)算法性能的重要依據(jù)。本文模擬多組在不同程度的隨機(jī)誤差及粗差影響下的監(jiān)測(cè)數(shù)據(jù)，將SQP算法、模矢法與遺傳算法進(jìn)行對(duì)比實(shí)驗(yàn)，具體步驟如下：

1)設(shè)計(jì)SQP算法、模矢法與遺傳算法的實(shí)驗(yàn)參數(shù)。模矢法的參數(shù)設(shè)計(jì)為：參數(shù)終止步長為0.0001；遺傳算法的參數(shù)設(shè)計(jì)為：種群個(gè)體數(shù)為100，遺傳代數(shù)為200代，交叉概率為0.9，變異概率為0.001；SQP算法的參數(shù)設(shè)計(jì)為：最大迭代次數(shù)200，終止誤差為0.0001。每個(gè)算法的初始參數(shù)的取值范圍均為：q∈[0.61,0.67]；b∈[0.25,0.35]；tanβ∈[1.97,2.77]；θ∈[80°,90°]；S1、S2、S3、S4的取值范圍均為S∈[32,72]。

2)將設(shè)計(jì)的采礦條件與沉陷預(yù)計(jì)參數(shù)真值帶入概率積分模型中，得到設(shè)計(jì)工作面觀測(cè)點(diǎn)的下沉與水平移動(dòng)值的真值，在真值的基礎(chǔ)上分別加入中誤差為5 mm、10 mm與20 mm的隨機(jī)誤差，并在不同的隨機(jī)誤差情況下，各仿真20組實(shí)驗(yàn)數(shù)據(jù)，并統(tǒng)計(jì)由SQP算法、模矢法與遺傳算法反演的各組參數(shù)的均值。

3)在實(shí)驗(yàn)步驟(1)中的誤差為10 mm的仿真數(shù)據(jù)的基礎(chǔ)上，隨機(jī)選取0～5個(gè)觀測(cè)站的下沉與水平觀測(cè)值，在此基礎(chǔ)上添加5σ-20σ的粗差，并統(tǒng)計(jì)由SQP算法、模矢法與遺傳算法反演參數(shù)的均值。

SQP算法、模矢法與遺傳算法反演的參數(shù)均值統(tǒng)計(jì)值結(jié)果，如表2所示。

表2 SQP算法、模矢法與遺傳算法反演參數(shù)計(jì)算結(jié)果

表2給出了不同組別的仿真實(shí)驗(yàn)，分析反演參數(shù)的結(jié)果可知：1)在不同誤差情況下，SQP算法對(duì)于q、b、tanβ和θ的反演結(jié)果并未有明顯變化，且與真值相比反演結(jié)果基本一致，最大誤差不超過1.2%；而對(duì)于S1、S2、S3、S4反演結(jié)果具有一定的隨機(jī)性，但與真值相比最大誤差仍不超過8%，仍具有一定的可靠性。2)在隨機(jī)誤差10 mm的基礎(chǔ)上加入了一定數(shù)目的粗差時(shí)，SQP算法反演結(jié)果并未發(fā)生改變，與之前相比變化不超過2.8%。3)隨著誤差的不斷增大，SQP算法、遺傳算法與模矢法反演參數(shù)的整體精度隨著誤差的增長有所降低，但結(jié)果與真值相差不大，仍具有一定的可靠性。4)對(duì)于q、b、tanβ和θ參數(shù)的反演結(jié)果，SQP算法受誤差因素影響較小，結(jié)果與模矢法和遺傳算法相比，均更加靠近真值。且在5 mm誤差的情況下，SQP算法優(yōu)勢(shì)更加明顯，反演結(jié)果幾乎可以認(rèn)為是真值。5)對(duì)于S1、S2、S3、S4參數(shù)的反演結(jié)果，SQP算法的精度隨著誤差的增大而有所降低。在5 mm與10 mm隨機(jī)誤差的影響下，SQP算法對(duì)于S1、S2、S3、S4的反演結(jié)果均優(yōu)于遺傳算法與模矢法；在20 mm隨機(jī)誤差的影響下，SQP算法對(duì)于S1、S2的反演結(jié)果均優(yōu)于模矢法與遺傳算法，但對(duì)于S3、S4反演的結(jié)果精度稍差于遺傳算法與模矢法；在10 mm粗差的影響下，SQP算法對(duì)于S1、S2、S3、S4的反演結(jié)果均優(yōu)于模矢法，但在S1、S2、S3方向上預(yù)計(jì)結(jié)果與遺傳算法相比稍差，但考慮到算法的穩(wěn)健性，其仍具有一定的可靠性。6)在監(jiān)測(cè)數(shù)據(jù)受到不同程度的隨機(jī)誤差與粗差的影響下，SQP算法反演參數(shù)的整體精度最高，遺傳算法次之，而模矢法較差。因此，表2證明了SQP算法對(duì)觀測(cè)數(shù)據(jù)中含有隨機(jī)誤差和粗差具有良好的抗干擾能力。

4 結(jié)論

通過分析常用反演參數(shù)算法存在的缺陷，提出使用SQP算法反演概率積分法模型中的參數(shù)，進(jìn)行大量的仿真實(shí)驗(yàn)，得到如下結(jié)論：

1)SQP算法不受初始值設(shè)置影響，且反演參數(shù)的結(jié)果與真值基本一致，準(zhǔn)確度高。

2)SQP算法對(duì)隨機(jī)誤差和粗差具有一定的抗干擾能力，其反演參數(shù)結(jié)果在隨機(jī)誤差與粗差的影響下，并未產(chǎn)生明顯的偏差，依然可靠。

3)在不同的誤差環(huán)境的影響下，SQP算法反演參數(shù)的精度基本高于模矢法與遺傳算法，但在少數(shù)情況下，個(gè)別方向上的拐點(diǎn)偏移距的預(yù)測(cè)精度可能稍低于遺傳算法，考慮到SQP算法的穩(wěn)健性，其結(jié)果依然可靠。

概率積分法預(yù)測(cè)結(jié)果受參數(shù)誤差與模型誤差影響，本文針對(duì)如何降低參數(shù)誤差影響進(jìn)行了改進(jìn)。但在現(xiàn)實(shí)礦區(qū)中，實(shí)測(cè)的下沉值和水平移動(dòng)值還會(huì)受到松散層厚度、斷層等其他因素影響，從而影響到模型自身精度。如何精化概率積分法的數(shù)學(xué)模型，使得函數(shù)預(yù)計(jì)曲線更符合實(shí)際的下沉曲線和水平移動(dòng)曲線，值得進(jìn)一步研究和探索。