江克貴，王 磊，池深深，魏 濤，蔣 創(chuàng)

( 1. 安徽理工大學 測繪學院，安徽 淮南 232001；2. 中國礦業(yè)大學 環(huán)境與測繪學院，江蘇 徐州 221116 )

我國能源資源稟賦具有“富煤、貧油、少氣”的特點。據(jù)預測，即使到2030年煤炭資源仍占中國能源需求的50%以上[1]。大規(guī)模高強度的煤炭資源開發(fā)會引發(fā)一系列礦山地質(zhì)環(huán)境災害問題，如含水層破壞、山體滑坡、地表塌陷和建構筑物損毀等。開采沉陷預計理論在評估潛在地質(zhì)災害，分析礦區(qū)沉降機理方面具有至關重要的作用，因而開采沉陷預計及防治的相關技術與理論一直是我國礦山領域的研究熱點，基于實測資料求取精準、可靠的開采沉陷預計參數(shù)也是預計模型應用過程中的難 點。

常規(guī)概率積分法存在以下幾個問題：沉降變形預計曲線邊緣收斂過快，與實際變形情況有較大差別；不能較好地適用于復雜多樣地質(zhì)采礦條件下的地表沉陷預計，不具有廣適應性；聯(lián)合時間函數(shù)的動態(tài)沉陷預計模型仍然存在一些問題。解決上述問題可以從2個途徑入手：① 構建精準廣適應的開采沉陷預計方法；② 提高預計參數(shù)的準確性和可靠性。對于途徑①，相關學者提出了新的開采沉陷預計方法或進行一定程度的改進。一類是針對不同地質(zhì)采礦條件而提出的沉陷預計方法，戴華陽[2]等提出了地表移動的矢量預計法，為解決大傾向煤層開采沉陷問題提供了預計依據(jù)；郭文兵[3]提出了全采工作面疊加預計方法，提高了深部大采寬條帶開采預計精度；文獻[4-6]提出了形式各異的開采沉陷預計方法，在一定程度上解決了多種不同地質(zhì)采礦條件下的開采沉陷預計問題。另一類是針對改進時間函數(shù)的動態(tài)開采沉陷預計方法，根據(jù)地表移動變形規(guī)律和采動程度，文獻[7]建立了一種新的Knothe時間函數(shù)求參模型，提高了地表動態(tài)下沉預測的準確性和可靠性；文獻[8]采用整體偏移修正和生長函數(shù)模型相結合的方法對時間函數(shù)進行優(yōu)化，進一步提高了地表下沉動態(tài)預計的精度。近年來，一些學者基于Boltzmann函數(shù)提出了新的開采沉陷預計模型[9]，通過分析函數(shù)曲線形態(tài)和參數(shù)特點發(fā)現(xiàn)，與概率積分法相比，Boltzmann函數(shù)在曲線邊緣部分收斂緩慢，且其預計模型參數(shù)與概率積分法模型參數(shù)可以相互轉(zhuǎn)化；另外，研究表明，該模型能夠廣適應于不同開采條件下的地表沉陷預計[10]。鑒于此，在Boltzmann函數(shù)的基礎上，筆者借鑒概率積分法模型的推導經(jīng)驗，通過理論分析，探討拐點偏移修正方法[11]，構建了一種耦合拐點偏移修正的Boltzmann函數(shù)預計模型。

針對途徑②，大量學者把智能優(yōu)化算法引入 到開采沉陷預計參數(shù)求取中來，相對于傳統(tǒng)方 法[12-14]，避免了求參速度慢、對初值依賴性強、易陷入局部最優(yōu)值、抗差能力差等缺陷，極大地提高了求取參數(shù)的準確性和可靠性。其中比較典型的有模矢法[15]、遺傳算法[16]、粒子群算法[17-18]等。筆者擬將FWA( Fireworks Algorithm )應用于反演預計模型參 數(shù)，提出了一種耦合拐點偏移修正的Boltzmann函數(shù)模型參數(shù)反演方法，并對該方法的抗差性能進行討論。

1 模型原理

玻爾茲曼( Boltzmann )函數(shù)模型在應用實踐中得到了廣泛的發(fā)展[19-20]。分析發(fā)現(xiàn)，Boltzmann函數(shù)可以生成一個S形曲線，跟半無限開采地表變形盆地主斷面的下沉曲線形態(tài)相似。對非充分開采、充分開采以及超充分開采的實測結果分析表明，Boltzmann函數(shù)可以較好地預計不同開采程度下的地表開采沉陷?；贐oltzmann函數(shù)的半無限地表變形盆地主斷面的下沉預計公式為

式中，W0=mq cos α，其中，m為煤層采厚，q為下沉系數(shù)，α為煤層傾角；R為新的主要影響半徑，與主要影響半徑r的物理意義相同，其主要影響下沉曲線的分布形態(tài)；S為拐點偏移距。

統(tǒng)計分析發(fā)現(xiàn)[14]，R≈r/4.13，其中，r=H/tan β，H為煤層采深，tan β為主要影響角正切。并且當x→∞時，有W( x )→W0；當x=s時，W( x )=W0/2；這些結果都符合開采沉陷一般規(guī)律。

對式( 1 )求微分，可得Boltzmann函數(shù)預計模型的單元影響函數(shù)，即

根據(jù)概率積分法的推導經(jīng)驗和彈性力學知識，可得單元水平移動函數(shù)為

對單元水平移動函數(shù)進行積分，可得地表沿x方向的水平移動為

令b'=B/R，統(tǒng)計分析可得b'=4.13b，則

式中，B為常數(shù)；b為水平移動系數(shù)。

煤壁附近采空區(qū)上方頂板會出現(xiàn)懸頂或過度破壞垮落，其作用效果相當于工作面采空區(qū)兩端煤壁向內(nèi)或向外分別平移了S的距離。文獻[11]表 明，概率積分法將拐點偏移距的影響全部加在了終采線和傾向上山處；由此，根據(jù)實測數(shù)據(jù)擬合求取的開采沉陷預計參數(shù)會因為擬合度低而出現(xiàn)偏差較大的現(xiàn)象。針對此問題，在Boltzmann函數(shù)預計模型的基礎上對拐點偏移進行修正，走向和傾向主斷面有限開采下沉和水平移動計算公式分別為

式中，

其中，D3和D1分別為工作面走向和傾向開采長度； S1，S2，S3，S4分別為工作面傾向下山、傾向上山、開切眼和終采線拐點偏移距；θ為開采影響傳播角。

井下開采尺寸為D3×D1的工作面引起地面任意一點( x，y )沿φ方向( x軸的正向逆時針到指定方向的角值 )的下沉和水平移動分別為

綜上所述，模型預計參數(shù)為

地質(zhì)采礦參數(shù)為

式中，?為工作面傾向方位角，直接參與工作面坐標轉(zhuǎn)換的過程。

構建模型的概化如式( 9 )所示，工作面地質(zhì)采礦參數(shù)和模型預計參數(shù)作為輸入變量控制地表相應預計點的形變；相反，當?shù)乇肀O(jiān)測點的形變獲取后，在工作面地質(zhì)采礦參數(shù)已知的情況下，求取精確可靠的模型預計參數(shù)顯得尤為重要；同時構建模型屬于高度非線性模型。因此，本文引入一種FWA。

2 模型參數(shù)反演方法

2.1 FWA

FWA是受到夜空中煙花爆炸的啟發(fā)而提出的一種具有局部搜索能力和全局搜索能力自調(diào)節(jié)機制的群體智能算法[21]。FWA的性能取決于爆炸算子( 即爆炸產(chǎn)生火花操作 )、變異算子( 高斯變異產(chǎn)生火花操作 )、映射規(guī)則( 映射到指定范圍操作 )和選擇策略( 選擇下一代火花操作 )[22-23]。

( 1 ) 爆炸算子

爆炸算子是FWA的核心組成部分，它由爆炸數(shù)目、爆炸半徑、位移操作3部分組成。假設煙花種群大小為N，對于煙花xi，其爆炸數(shù)目和爆炸半徑的計算公式分別為

式中，Di為第i個煙花爆炸產(chǎn)生的火花數(shù)目；D為常 數(shù)，用來調(diào)整產(chǎn)生的爆炸火花數(shù)目；f( xi)為第i個煙花的適應度值；ymax為當前種群中最大適應度值，ymax=max( f( xi) )；ξ為一個機器最小量，用來避免除零操作；Ai為第i個煙花的爆炸半徑；A為常數(shù)，用來調(diào)整爆炸半徑的大?。粂min為當前種群中最小適應度值，ymin=min( f( xi) )。

為了避免適應度較好的煙花位置產(chǎn)生過多的爆炸火花，同時適應度較差的煙花位置產(chǎn)生過少的爆炸火花，對爆炸的火花數(shù)目進行如下限制：

式中，a，b為常數(shù)；round( )為根據(jù)四舍五入原則的取整函數(shù)。

位移操作是對煙花的每一維進行位移，即

式中，rand( 0，Ai)表示在幅度Ai內(nèi)生成的均勻隨機 數(shù)；xik為第i個煙花在第k個維度上的位置信息。

( 2 ) 變異算子

為了進一步增加種群的多樣性，F(xiàn)WA在隨機選擇的煙花中執(zhí)行高斯變異操作，即

式中，Gaussian( 1，1 )表示高斯分布。

( 3 ) 映射規(guī)則

爆炸算子和變異算子產(chǎn)生的火花可能會超出邊界，利用映射規(guī)則將超出域值的火花映射到一個新的位置。

( 4 ) 選擇策略

在候選者集合中( 包括煙花、爆炸火花和高斯變異火花 )，適應度最好的個體會首先被確定性地選擇到下一代作為煙花，剩下的后代煙花采用輪盤賭法進行選擇。候選者xi被選擇的概率計算公式 為

式中，R( xi)為當前個體xi與其他個體的距離之和； d( xi，xj)表示任意2個個體xi和xj之間的歐氏距離；xj∈K表示第 j個位置屬于集合K。

2.2 方法構建

綜合地表移動變形、預計模型以及FWA原理，提出了一種耦合拐點偏移修正的Boltzmann函數(shù)模型參數(shù)反演方法( BFWA-MSPP )，具體步驟如下：

( 1 ) 數(shù)據(jù)準備。首先，利用地表移動觀測站觀測目標區(qū)域地表進入穩(wěn)定期后的下沉和水平移動量，獲取工作面地質(zhì)采礦條件參數(shù)；然后，選取適當?shù)腇WA參數(shù)，確定開采沉陷預計參數(shù)取值空間，進而生成煙花初始種群。

( 2 ) 適應度評價。煙花初始種群用作迭代種 群；基于構建的Boltzmann函數(shù)預計模型，結合獲取的地質(zhì)采礦條件，第i組開采沉陷參數(shù)預計目標區(qū)域任意點 j的下沉和水平移動分別為 Wpi,j和 Upi,j；假設地表點 j的實測下沉和水平移動值分別為 Wmj和Umj，以預計值和實測值之差絕對值的和為準則進行適應度評價，適應度的構造如式( 17 )所示；然后判斷是否滿足迭代終止條件，如果滿足，輸出最優(yōu)解，否則執(zhí)行步驟( 3 )。

( 3 ) 生成火花。根據(jù)煙花迭代種群的適應度，結合FWA參數(shù)，計算每個煙花的爆炸半徑與爆炸火花數(shù)目；利用式( 13 )對各個煙花進行位移偏移，生成多個爆炸火花；另外，在煙花種群中，利用變異算子生成多個變異火花。

( 4 ) 映射及后代選擇。利用映射規(guī)則把超出預計參數(shù)范圍的火花映射到可行域范圍內(nèi)。根據(jù)選擇策略，在候選者集合中選擇下一代。

( 5 ) 循環(huán)執(zhí)行步驟( 2 )～( 4 )，直到滿足終止條件，輸出最優(yōu)開采沉陷預計參數(shù)，模型參數(shù)反演方法技術 路線如圖1所示。

圖1 BFWA-MSPP技術路線 Fig. 1 Technology roadmap of BFWA-MSPP

3 仿真試驗

3.1 數(shù)據(jù)模擬

以淮南礦區(qū)煤系地層為背景，模擬工作面上覆巖層巖性為中硬，工作面地質(zhì)采礦條件為：煤層采高m=3 m，煤層采深H=400 m，煤層傾角α=5°，工作面傾向方位角?=0°，工作面開采尺寸D3×D1=800 m×200 m。分別沿工作面走向和傾向主斷面 模擬布設了2條監(jiān)測線E和N，測點間距為30 m，共 布設88個監(jiān)測點，模擬的工作面和布設的監(jiān)測線 如圖2所示。基于構建模型和模擬的開采沉陷預計參數(shù)( 表1 )預計開采沉陷穩(wěn)定后的下沉和水平移 動。

3.2 反演試驗

根據(jù)模擬的地質(zhì)采礦條件和預計的形變數(shù)據(jù)，利用本文構建的BFWA-MSPP進行開采沉陷預計參數(shù)反演試驗，求參結果見表1。

從表1可以看出，反演參數(shù)q，tan β，b，θ的相對誤差最大不超過1.2%，中誤差最大不超過1.2；反演拐點偏移距S的相對誤差和中誤差分別控制在6.2%和6.5 m以內(nèi)；此外，BFWA-MSPP擬合下沉和水平移動誤差為-11.6～5.7 mm，中誤差為±6.1 mm。試驗結果表明，構建的方法可以準確可靠地求取開采沉陷預計參數(shù)。

圖2 模擬工作面和監(jiān)測線示意 Fig. 2 Schematic diagram of simulated the working face and the observation lines

表1 仿真試驗BFWA-MSPP求參結果 Table 1 Solved results of BFWA-MSPP for simulation experiment

3.3 抗差試驗

在實際測量中，觀測值會不可避免地摻入一些隨機誤差和粗差[24-25]，隨機誤差和粗差的存在勢必影響求參結果的質(zhì)量。為了探討B(tài)FWA-MSPP抵抗觀測誤差的能力，基于3.1節(jié)的模擬數(shù)據(jù)，設計試驗如下：在所有監(jiān)測點的下沉值和水平移動值中分別加入均值為0、標準差為30 mm和10 mm的高斯噪 聲；在此基礎上，將所有監(jiān)測值再加入-200 ～200 mm( 服從均勻分布 )的粗差；然后，以包含測量誤差的數(shù)據(jù)為觀測值，利用BFWA-MSPP進行求參，求參結果見表2。

表2 BFWA-MSPP抗差求參結果 Table 2 Anti-error result of BFWA-MSPP

從表2可以看出，在所有觀測點加入觀測誤差后，反演q，tan β，b，θ的相對誤差不超過13%，反演 拐點偏移距S的相對誤差控制在32%以內(nèi)；就擬合效果來看( 圖3 )，在含有觀測誤差干擾的情況下，擬合曲線和模擬觀測曲線較吻合。需要說明的是，本 文為了極限測試算法的穩(wěn)定性加入了較大的誤差，但在工程實踐中很少有如此大的觀測誤差( 通過精密測量儀器和數(shù)據(jù)處理方法，觀測誤差可以得到控制 )。試驗結果表明，BFWA-MSPP可以有效抵抗觀測誤差的影響，在監(jiān)測實施困難區(qū)域，可適當降低地表形變監(jiān)測精度要求。

圖3 BFWA-MSPP抗差擬合下沉和水平移動效果 Fig. 3 Effect diagram of fitting the subsidence and the horizontal movement by anti-error of BFWA-MSPP

4 工程應用

4.1 礦區(qū)概況

以淮南礦區(qū)顧橋煤礦1414( 1 )工作面為例，該工作面采用綜合機械化采煤，一次采全高，全部垮落法管理頂板。工作面地質(zhì)采礦條件為：平均采高3 m，工作面傾角5°，平均采深735 m，工作面開采尺寸為2 115.433 m×241.210 m，走向為超充分采動，傾向為非充分采動，總體為非充分采動。地表移動觀測站共布設6個控制點和145個監(jiān)測點，點間距 為30/60 m，傾向觀測線長度為1 500 m，走向觀測 線長度為3 480 m，工作面和觀測線的布設如圖4所 示。

圖4 工作面和觀測線實地布設 Fig. 4 Field layout of the working face and the observation lines

4.2 試驗設計及結果分析

分別利用BFWA-MSPP和PIM-MSPP對1414( 1 )工作面地表移動觀測站最后一期觀測數(shù)據(jù)( 地表進入穩(wěn)定期 )進行開采沉陷預計參數(shù)反演，反演結果見表3。

表3 工程應用BFWA-MSPP和PIM-MSPP求參結果比較 Table 3 Comparison between the engineering applications of BFWA-MSPP and PIM-MSPP

從表3中可以看出，2種方法反演參數(shù)的差值主要集中在下沉系數(shù)q和拐點偏移距S處，就差值所占總體比例來說，下沉系數(shù)的差值百分比超過了 10%，拐點偏移距的差值百分比平均約為50%；此 外，開切眼和終采線處的拐點偏移距過大，主要由于工作面走向兩側(cè)形成較大砌體梁結構所致。從預計擬合效果來看( 圖5 )，無論是曲線邊緣部分，還是曲線拐點附近，BFWA-MSPP擬合值更加靠近實測值，2種方法預計擬合下沉和水平移動中誤差分別為±83 mm和±108 mm。

5 結 論

( 1 ) 基于Boltzmann函數(shù)，融入拐點偏移修正方法，構建了耦合拐點偏移修正的Boltzmann函數(shù)預計模型，由于該模型屬于高度非線性模型，通過引入FWA，建立了耦合拐點偏移修正的Boltzmann函數(shù)模型參數(shù)反演方法( BFWA-MSPP )。

圖5 BFWA-MSPP和PIM-MSPP預計擬合下沉和水平移動比較 Fig. 5 Comparison of BFWA-MSPP and PIM-MSPP fitting the subsidence and the horizontal movement

( 2 ) 仿真試驗結果表明，BFWA-MSPP求參的平均相對誤差約為2.1%，平均中誤差為2.3，BFWAMSPP預計擬合下沉和水平移動中誤差為±6.1 mm。在觀測值包含一定的隨機誤差和粗差后，反演參數(shù)的平均相對誤差為13.8%，BFWA-MSPP可以有效抵抗觀測誤差的影響。

( 3 ) 將BFWA-MSPP和PIM-MSPP應用在顧橋 煤礦1414( 1 )工作面，對地表開采穩(wěn)定后工作面的 開采沉陷預計參數(shù)進行了求取，預計擬合下沉和 水平移動中誤差分別為±83 mm 和±108 mm；BFWA-MSPP求參結果為：q=1.07，tan β=1.91，b=0.41，θ=88.32°，S1=-18.89 m，S2=14.74 m，S3=48.89 m，S4=64.42 m。