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

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

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

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

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

1 模型原理

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

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

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

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

根據(jù)概率積分法的推導(dǎo)經(jīng)驗(yàn)和彈性力學(xué)知識(shí)，可得單元水平移動(dòng)函數(shù)為

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

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

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

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

式中，

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

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

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

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

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

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

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

2.1 FWA

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

( 1 ) 爆炸算子

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

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

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

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

位移操作是對(duì)煙花的每一維進(jìn)行位移，即

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

( 2 ) 變異算子

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

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

( 3 ) 映射規(guī)則

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

( 4 ) 選擇策略

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

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

2.2 方法構(gòu)建

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

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

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

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

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

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

圖1 BFWA-MSPP技術(shù)路線 Fig. 1 Technology roadmap of BFWA-MSPP

3 仿真試驗(yàn)

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

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

3.2 反演試驗(yàn)

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

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

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

表1 仿真試驗(yàn)BFWA-MSPP求參結(jié)果 Table 1 Solved results of BFWA-MSPP for simulation experiment

3.3 抗差試驗(yàn)

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

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

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

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

4 工程應(yīng)用

4.1 礦區(qū)概況

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

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

4.2 試驗(yàn)設(shè)計(jì)及結(jié)果分析

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

表3 工程應(yīng)用BFWA-MSPP和PIM-MSPP求參結(jié)果比較 Table 3 Comparison between the engineering applications of BFWA-MSPP and PIM-MSPP

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

5 結(jié) 論

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

圖5 BFWA-MSPP和PIM-MSPP預(yù)計(jì)擬合下沉和水平移動(dòng)比較 Fig. 5 Comparison of BFWA-MSPP and PIM-MSPP fitting the subsidence and the horizontal movement

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

( 3 ) 將BFWA-MSPP和PIM-MSPP應(yīng)用在顧橋 煤礦1414( 1 )工作面，對(duì)地表開采穩(wěn)定后工作面的 開采沉陷預(yù)計(jì)參數(shù)進(jìn)行了求取，預(yù)計(jì)擬合下沉和 水平移動(dòng)中誤差分別為±83 mm 和±108 mm；BFWA-MSPP求參結(jié)果為：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。