徐可心,何 榮,白偉森
(河南理工大學(xué) 測(cè)繪與國(guó)土信息工程學(xué)院,河南 焦作 454000)
煤礦地下開(kāi)采造成地表移動(dòng)形變,引起地表裂縫、建筑物損害等,部分地區(qū)出現(xiàn)積水、水土流失等地質(zhì)災(zāi)害,嚴(yán)重破壞礦區(qū)的生態(tài)環(huán)境。因此,有必要通過(guò)對(duì)礦區(qū)地表形變監(jiān)測(cè),獲取變形數(shù)據(jù)、求算相關(guān)參數(shù),總結(jié)地表移動(dòng)盆地的變形規(guī)律,為維護(hù)礦區(qū)經(jīng)濟(jì)可持續(xù)發(fā)展及生態(tài)文明建設(shè)提供基礎(chǔ)數(shù)據(jù)支撐[1]。
概率積分法是我國(guó)目前開(kāi)采沉陷領(lǐng)域科學(xué)研究常用方法之一[2],構(gòu)建了地表變形與地下開(kāi)采空間的關(guān)系模型。不同的開(kāi)采方法和地質(zhì)采礦條件,預(yù)計(jì)參數(shù)具有較大的差異性,因此,科學(xué)確定預(yù)計(jì)參數(shù)大小對(duì)地表變形量計(jì)算至關(guān)重要。其參數(shù)求算主要有特征點(diǎn)求參、曲線擬合和曲面擬合等方法,曲面擬合法將曲線擬合法求參的原理推廣到整個(gè)下沉盆地,充分利用非主斷面地表變形監(jiān)測(cè)數(shù)據(jù),得到更加精確的地表移動(dòng)變形參數(shù)值[3]。李春意[4]等將曲面擬合法用于求取矩形工作面的概率積分法參數(shù),并驗(yàn)證了非線性曲面擬合法用于開(kāi)采沉陷預(yù)計(jì)的實(shí)用性。該方法不要求觀測(cè)站位于下沉盆地的主斷面上,是非主斷面設(shè)站求取概率積分法參數(shù)的主要方法。但其對(duì)初值選取要求較高,工作量大、求參速度慢。隨著計(jì)算機(jī)技術(shù)的發(fā)展,為提高參數(shù)反演的效率,不少學(xué)者將智能算法應(yīng)用到了概率積分法反演參數(shù)的研究中,優(yōu)化算法有模矢法[5]、遺傳算法[6]、果蠅算法[7]、神經(jīng)網(wǎng)絡(luò)算法[8]以及粒子群算法[9]等。與其他的尋優(yōu)算法相比, 粒子群算法擁有較好的收斂能力,計(jì)算速度更快[10]。徐孟強(qiáng)[11]等首次將粒子群算法引入了概率積分法參數(shù)反演,認(rèn)為粒子群算法可以有效解決隨機(jī)誤差和測(cè)點(diǎn)缺失的問(wèn)題,存在全局搜索能力差、容易提前收斂等問(wèn)題;葉偉[12]等將自適應(yīng)魚(yú)群算法用于曲面擬合求參,得到了較為可靠的擬合下沉曲面,所用方法運(yùn)算復(fù)雜、計(jì)算效率較低。
基于此,針對(duì)傳統(tǒng)粒子群算法存在的局限性,使用改進(jìn)標(biāo)準(zhǔn)粒子群算法,進(jìn)行曲面擬合概率積分法參數(shù)反演,計(jì)算參數(shù)擬合中誤差和下沉擬合誤差,驗(yàn)證改進(jìn)標(biāo)準(zhǔn)粒子群算法用于概率積分法參數(shù)反演的有效性。
相比主斷面設(shè)站法,非主斷面設(shè)站具有布設(shè)靈活的特點(diǎn),既方便觀測(cè),又減少了對(duì)農(nóng)田的影響,被廣泛用于礦區(qū)的地表變形監(jiān)測(cè)。觀測(cè)站位置的選取對(duì)參數(shù)反演精度有較大影響,為保證計(jì)算精度,觀測(cè)站應(yīng)盡量布設(shè)在變形值較大位置,或下沉盆地的特征點(diǎn)處[13]。
試驗(yàn)區(qū)位于河北省邯鄲市某礦,區(qū)內(nèi)地形為低丘陵,地勢(shì)南高北低,有較淺沖溝,部分有基巖出露。2516 工作面位于205 采區(qū)運(yùn)輸下山北翼,北面為充填集中軌道巷,西面為已開(kāi)采的其他工作面,東南方向?yàn)榇迩f。工作面開(kāi)采條件見(jiàn)表1。
表1 工作面開(kāi)采條件Table 1 Mining conditions of working face
為監(jiān)測(cè)礦區(qū)的地表移動(dòng)變形,2516 工作面上方布設(shè)了走向與傾向2 條觀測(cè)線,受地表村莊分布和地形限制,以非主斷面設(shè)站法布設(shè)。觀測(cè)線布設(shè)和井上下對(duì)照?qǐng)D如圖1,走向觀測(cè)線長(zhǎng)671 m,共布設(shè)17 個(gè)水準(zhǔn)觀測(cè)點(diǎn),傾向觀測(cè)線長(zhǎng)1 167 m,共布設(shè)28 個(gè)水準(zhǔn)觀測(cè)點(diǎn)。2516 工作面共進(jìn)行5 次觀測(cè),其中走向最大下沉點(diǎn)下沉值為572 mm,傾向最大下沉點(diǎn)為走向與傾向觀測(cè)線交點(diǎn),下沉值為556 mm。
曲面擬合法不受觀測(cè)站布設(shè)方法限制,但僅適用于矩形工作面。由圖1 可知,2516 工作面為非主斷面設(shè)站的矩形工作面,可使用曲面擬合法反演概率積分法參數(shù)。由概率積分法可知,矩形工作面下沉盆地內(nèi)曲面擬合任意點(diǎn)下沉的函數(shù)W(x,y)為:
圖1 地表移動(dòng)觀測(cè)站布設(shè)與井上下對(duì)照?qǐng)DFig.1 Layout of surface movement observation stations and comparison between upper and lower wells
式中:W0為實(shí)測(cè)最大下沉值,mm;H 為走向主斷面平均采深,m;H1、H2分別為下山和上山方向的采深,m;β、β1、β2為主要影響角,(°);D1、D3為矩形工作面傾向和走向主斷面寬度,m;S1、S2分別為上山和下山方向拐點(diǎn)偏移距,m;S3、S4分別為走向左和走向右拐點(diǎn)偏移距,m;θ0為開(kāi)采影響傳播角,(°);α 為煤層傾角,(°)。
粒子群優(yōu)化算法(Particle Swarm Optimization,PSO)是1 種群智能優(yōu)化算法,源于對(duì)鳥(niǎo)群社會(huì)行為的研究。Reynolds[14]于1986 年設(shè)計(jì)了1 個(gè)再現(xiàn)鳥(niǎo)群聚集行為的人工生命系統(tǒng),隨后Eberhart 和Kennedy[15]基于該系統(tǒng)提出了粒子群算法。Y. Shi[16]等學(xué)者引入了慣性權(quán)重ω,并修正了粒子速度公式,引入ω 的粒子群算法被稱為標(biāo)準(zhǔn)粒子群算法(Standard Particle Swarm Optimization,SPSO)。
假設(shè)在D 維的搜索空間中,N 個(gè)粒子組成1 個(gè)群體,第i 個(gè)粒子的位置表示為xi=[xi1,xi2,…,xiD],其飛行速度表示為vi=[vi1,vi2,…,viD],當(dāng)前搜索到的最優(yōu)位置為pi=[pi1,pi2,…,piD],整個(gè)群體當(dāng)前搜索到的最優(yōu)位置為pg=[pi1,pi2,…,piD]。由此得到粒子速度與位置的更新公式:
慣性權(quán)重ω 是SPSO 算法的重要參數(shù),用于控制算法開(kāi)發(fā)和搜索能力。當(dāng)ω 取值較大,算法全局尋優(yōu)能力較強(qiáng),局部尋優(yōu)能力較弱;當(dāng)ω 取值較小,正好相反。在算法搜索過(guò)程中將ω 設(shè)置為動(dòng)態(tài),可以平衡收斂速度和搜索能力,使算法優(yōu)化結(jié)果更精確。
PSO 算法的優(yōu)勢(shì)在于收斂能力較好,但面對(duì)復(fù)雜多維問(wèn)題,PSO 算法常出現(xiàn)早熟,即提前收斂陷入局部最優(yōu)解[17]。針對(duì)PSO 算法容易早熟的問(wèn)題,在SPSO 算法的速度更新公式中引入漸進(jìn)收斂的方法,改善算法全局收斂性,提高算法搜索能力。改進(jìn)SPSO 算法粒子速度迭代公式為:
由式(1),可將下沉盆地內(nèi)地表任意點(diǎn)的下沉值WPos表示為:
式中:Pos 為概率積分法的參數(shù)矩陣,Pos=[q,tanβ,tanβ1,tanβ2,S1,S2,S3,S4];q 為下沉系數(shù);(x,y)為觀測(cè)點(diǎn)坐標(biāo)。
基于改進(jìn)SPSO 算法的概率積分法參數(shù)反演流程如圖2。
圖2 SPSO 算法概率積分法參數(shù)反演流程Fig.2 SPSO algorithm probability integration method parameter inversion process
在改進(jìn)SPSO 算法中,應(yīng)先設(shè)置算法參數(shù),確定適應(yīng)度函數(shù),然后通過(guò)迭代得到最優(yōu)解,反演流程的具體步驟如下。
1)初始化種群。設(shè)置最大迭代次數(shù)、自變量個(gè)數(shù)、粒子的最大速度以及粒子約束條件為搜索空間,在速度空間和搜索空間上初始化粒子飛行速度和位置,設(shè)置種群規(guī)模,生成初始群體。使用MATLAB 語(yǔ)言對(duì)求參模型進(jìn)行編程,假設(shè)空間內(nèi)共有100 個(gè)粒子,學(xué)習(xí)因子c1=c2=2,nmax=500。粒子維數(shù)由所求公式參數(shù)數(shù)量決定,概率積分法參數(shù)矩陣中有8 個(gè)參數(shù),分別為:下沉系數(shù)q,主要影響角正切tanβ、tanβ1、tanβ2以及拐點(diǎn)偏移距Si(i=1,2,3,4),確定粒子維數(shù)為8。為平衡算法搜索能力,設(shè)置慣性權(quán)重隨迭代次數(shù)增加線性遞減,令ωmax為0.9,ωmin為0.4。鑒于工作面實(shí)際情況,結(jié)合礦區(qū)現(xiàn)有巖移參數(shù),確定粒子約束條件,即參數(shù)取值區(qū)間,參數(shù)取值上下限見(jiàn)表2。
表2 參數(shù)取值上下限Table 2 Upper and lower limits of parameter values
2)定義適應(yīng)度函數(shù)。粒子每更新1 次位置就計(jì)算1 次適應(yīng)度值,適應(yīng)度值反應(yīng)粒子的優(yōu)劣程度。根據(jù)誤差平方和最小原則,預(yù)計(jì)值和實(shí)測(cè)值之差的平方和越小,參數(shù)反演精度越高。構(gòu)建適應(yīng)度函數(shù)為:
式中:RSS 為擬合值與實(shí)測(cè)值之差的平方和;m為觀測(cè)點(diǎn)總個(gè)數(shù);Wi為實(shí)測(cè)下沉值,mm;wi為擬合下沉值,mm。
3)根據(jù)式(3)、式(5)更新粒子的速度和位置,并處理超出約束條件的粒子。迭代過(guò)程中,粒子通過(guò)個(gè)體極值和群體極值更新速度和位置。個(gè)體極值為每個(gè)粒子找到最優(yōu)解,從個(gè)體最優(yōu)解中找到全局最優(yōu)解,與歷史全局最優(yōu)解進(jìn)行比較,不斷更新得到最終結(jié)果。
4)當(dāng)相鄰2 代全局最優(yōu)解之間的偏差最小,或循環(huán)達(dá)到最大迭代次數(shù)時(shí),循環(huán)終止,輸出反演參數(shù)。在本研究構(gòu)建的模型中,當(dāng)相鄰2 代全局最優(yōu)解適應(yīng)度值差值為0 時(shí),認(rèn)為找到最優(yōu)解,或迭代次數(shù)達(dá)到500 次時(shí)循環(huán)停止。一般情況下迭代次數(shù)越多精度越高,但多余迭代次數(shù)會(huì)降低運(yùn)算效率,為避免無(wú)效迭代并保證運(yùn)算精度,通過(guò)試驗(yàn)確定最大迭代次數(shù)為500。循環(huán)停止輸出最優(yōu)解,即反演所得概率積分法參數(shù)。
以2516 工作面實(shí)測(cè)下沉數(shù)據(jù)為基礎(chǔ)(圖1),分別使用改進(jìn)SPSO 算法和PSO 算法進(jìn)行參數(shù)反演,對(duì)比其求參結(jié)果和擬合中誤差。為避免試驗(yàn)偶然性,2 種算法以相同設(shè)置各進(jìn)行10 次試驗(yàn),取試驗(yàn)結(jié)果的平均值為參數(shù)值,所得結(jié)果和礦區(qū)經(jīng)驗(yàn)值見(jiàn)表3。礦區(qū)經(jīng)驗(yàn)值參考205 采區(qū)已有地表移動(dòng)變形參數(shù),結(jié)合2516 工作面開(kāi)采技術(shù)確定。
表3 反演參數(shù)對(duì)比Table 3 Comparison of inversion parameters
由表3 可知,除主要影響角正切tan β 外,改進(jìn)SPSO 算法反演參數(shù)擬合中誤差均小于PSO 算法反演參數(shù)擬合中誤差,即改進(jìn)SPSO 算法反演參數(shù)波動(dòng)性小于PSO 算法,證明改進(jìn)SPSO 算法穩(wěn)定性優(yōu)于PSO 算法。
改進(jìn)SPSO 算法和PSO 算法的迭代次數(shù)對(duì)比如圖3,在算法10 次運(yùn)算過(guò)程中,改進(jìn)SPSO 算法出現(xiàn)1 次收斂過(guò)快未求得最優(yōu)解,PSO 算法出現(xiàn)3 次。且改進(jìn)SPSO 算法的平均迭代次數(shù)為47.3 次,PSO 算法平均迭代次數(shù)為300.4 次??梢钥闯龈倪M(jìn)SPSO 算法收斂速度遠(yuǎn)高于PSO 算法。
圖3 迭代次數(shù)對(duì)比Fig.3 Comparison of iteration times
為分析改進(jìn)SPSO 算法反演概率積分法參數(shù)的效果,將改進(jìn)SPSO 算法和PSO 算法反演參數(shù)值和礦區(qū)經(jīng)驗(yàn)值分別代入概率積分法公式得到擬合下沉曲面,與實(shí)測(cè)數(shù)據(jù)對(duì)比,并計(jì)算其擬合誤差。
下沉擬合誤差計(jì)算公式為:
由式(8)得,改進(jìn)SPSO 算法反演參數(shù)擬合誤差σ=4.66%,PSO 算法擬合誤差σ=4.67%,礦區(qū)經(jīng)驗(yàn)值擬合誤差σ=10.2%。改進(jìn)SPSO 算法反演參數(shù)精度和PSO 算法相當(dāng),高于礦區(qū)經(jīng)驗(yàn)值。改進(jìn)SPSO 算法反演參數(shù)和礦區(qū)經(jīng)驗(yàn)值擬合下沉曲面如圖4。算法擬合下沉值及實(shí)測(cè)值對(duì)比如圖5。
由圖4 可知,改進(jìn)SPSO 算法反演參數(shù)擬合下沉曲面更接近實(shí)際下沉盆地的形狀,擬合效果較好。由圖5 可知,工作面走向方向擬合效果較好,傾向方向的部分監(jiān)測(cè)點(diǎn)實(shí)測(cè)值大于擬合值。由圖1 可知,2516 工作面西側(cè)為2 個(gè)已開(kāi)采工作面,且與2516工作面相距較近,2516 工作面的開(kāi)采勢(shì)必會(huì)引起老采空區(qū)的活化。根據(jù)覆巖“活化”機(jī)理[18],老采空區(qū)會(huì)因?yàn)樯细矌r層的“活化”而產(chǎn)生二次破裂,造成地表下沉值大于預(yù)計(jì)值。試驗(yàn)反演概率積分法參數(shù)僅考慮單一工作面開(kāi)采條件,2516 工作面在實(shí)際開(kāi)采過(guò)程中受到了2314 和2509 已采工作面的影響,使得傾向觀測(cè)線上山方向部分監(jiān)測(cè)點(diǎn)擬合效果不理想。
圖4 下沉擬合曲面及實(shí)測(cè)值對(duì)比Fig.4 Subsidence fitting surface and comparison of measured values
圖5 算法擬合下沉值及實(shí)測(cè)值對(duì)比Fig.5 Algorithm fitting sinking values and comparison of measured values
為進(jìn)一步分析改進(jìn)SPSO 算法反演概率積分法參數(shù)的可靠性,基于改進(jìn)SPSO 算法的反演參數(shù)值和概率積分法原理,對(duì)2516 工作面地表移動(dòng)觀測(cè)站的傾斜和曲率變形進(jìn)行了預(yù)測(cè)。算法預(yù)計(jì)值與實(shí)測(cè)值對(duì)比如圖6。
圖6 算法預(yù)計(jì)值與實(shí)測(cè)值對(duì)比Fig.6 Comparison between the predicted values of the algorithm and the measured values
由圖6 可知,算法傾斜及曲率變形預(yù)計(jì)值和實(shí)測(cè)值相比,整體變化趨勢(shì)基本一致,除少部分異常點(diǎn)外,整體誤差較低,預(yù)計(jì)效果較好。證明改進(jìn)SPSO算法反演參數(shù)具有較高可靠性。
1)針對(duì)非主斷面設(shè)站反演概率積分法參數(shù)的問(wèn)題,將漸進(jìn)收斂引入SPSO 算法,建立了基于改進(jìn)SPSO 算法的概率積分法參數(shù)反演模型。研究表明,該方法參數(shù)反演運(yùn)算效率高、擬合效果較好,通過(guò)該方法反演的參數(shù)具有較高可靠性。
2)反演參數(shù)值擬合中誤差和下沉擬合誤差的對(duì)比分析表明,改進(jìn)SPSO 算法反演參數(shù)擬合中誤差整體優(yōu)于PSO 算法,參數(shù)反演精度高于礦區(qū)經(jīng)驗(yàn)值,有效提高了參數(shù)反演準(zhǔn)確性。和PSO 算法相比,改進(jìn)SPSO 算法具有更好的收斂能力和穩(wěn)定性。
3)由下沉、傾斜、曲率預(yù)計(jì)值與實(shí)測(cè)值的對(duì)比分析可知,改進(jìn)SPSO 算法反演參數(shù)預(yù)計(jì)的擬合下沉盆地比較符合實(shí)際情況,對(duì)于提高礦區(qū)開(kāi)采沉陷預(yù)計(jì)精度有一定的參考價(jià)值。