0 引言
    目前,由于煤礦盜采具有隱蔽性強(qiáng)、不易被發(fā)現(xiàn)的特點(diǎn)，煤礦監(jiān)管部門迫切需要有一套檢測(cè)精度高、檢測(cè)周期短的監(jiān)測(cè)設(shè)備來(lái)對(duì)地下巷道進(jìn)行實(shí)時(shí)監(jiān)測(cè)和定位。微地震監(jiān)測(cè)是目前礦區(qū)防盜采中最精確、最及時(shí)、信息最豐富的監(jiān)測(cè)手段之一[1]。微地震監(jiān)測(cè)技術(shù)最主要的任務(wù)就是確定震源的位置[2]。而微震監(jiān)測(cè)系統(tǒng)中的拾震器陣列的空間布置是微震監(jiān)測(cè)系統(tǒng)能否有效并可靠運(yùn)行的關(guān)鍵。國(guó)內(nèi)外關(guān)于震源定位技術(shù)中拾震器陣列的空間布置對(duì)系統(tǒng)精度和可靠性的影響及其優(yōu)化的研究的相關(guān)文獻(xiàn)都非常少[3，4]，因此，非常有必要對(duì)微震監(jiān)測(cè)中拾震器陣列的空間布置進(jìn)行研究[5]。本文提出利用神經(jīng)網(wǎng)絡(luò)遺傳算法和微變網(wǎng)格射線追蹤正演模擬方法在微震定位系統(tǒng)中進(jìn)行拾震器布陣研究來(lái)為微震監(jiān)測(cè)系統(tǒng)選出科學(xué)合理的拾震器布陣方案。
1 微震波震源定位正反演方法和布陣方案設(shè)計(jì)
1.1 PSO-BP震源定位反演
    實(shí)際上求解震源定位問(wèn)題就是解非線性方程組的問(wèn)題，本文通過(guò)PSO-BP算法求解，利用PSO算法進(jìn)行BP網(wǎng)絡(luò)參數(shù)的優(yōu)化，不僅避免了BP網(wǎng)絡(luò)原有算法在訓(xùn)練時(shí)易陷入局部極小的問(wèn)題，而且提高了BP網(wǎng)絡(luò)的訓(xùn)練速成度。所謂參數(shù)的優(yōu)化，就是尋找最優(yōu)的BP網(wǎng)絡(luò)的網(wǎng)絡(luò)權(quán)值和閾值，使得BP網(wǎng)絡(luò)的全局誤差最小化，從而尋找震源位置。
    PSO-BP 混合算法具體步驟如下：
    (1)創(chuàng)建三層BP網(wǎng)絡(luò),構(gòu)造BP網(wǎng)絡(luò)的輸入層、隱含層和輸出層，并隨機(jī)產(chǎn)生N組連接權(quán)值和閾值分布作為PSO算法的初始種群，評(píng)價(jià)每個(gè)微粒的適應(yīng)度。
    (2)對(duì)每個(gè)微粒，將其適應(yīng)度值與其經(jīng)過(guò)的最好位置p_best作比較，如果較好，則將其作為當(dāng)前的最好位置p_best。
    (3)對(duì)每個(gè)微粒，將其當(dāng)前的最好位置p_best與當(dāng)前整個(gè)群體中所有微粒發(fā)現(xiàn)的最好位置g_best作比較，如果較好，則將其作為當(dāng)前群體的最好位置g_best。
    (4)用PSO公式(1)、(2)進(jìn)行粒子的速度和位移的更新：
    v_i(n+1)=w·v_i(n)+c₁·r₁(p_i-x_i(n))+c₂·r₂(p_g-x_i(n))                  (1)
    x_i(n+1)=x_i(n)+v_i(n)                                                           (2)
    (5)用種群最優(yōu)位置gbest更新BP網(wǎng)絡(luò)的權(quán)值和閾值，PSO運(yùn)行結(jié)束。迭代終止條件根據(jù)具體問(wèn)題選為最大迭代次數(shù)或(和)粒子群迄今為止搜索到的最優(yōu)位置滿足預(yù)定目標(biāo)最小適應(yīng)閾值。
    (6)BP算法利用權(quán)值和閥值等網(wǎng)絡(luò)參數(shù)震源位置進(jìn)行計(jì)算。
1.2 微變網(wǎng)格射線追蹤正演模擬
    射線理論正演能夠?qū)崿F(xiàn)選擇最佳的微震波在介質(zhì)中的運(yùn)動(dòng)軌跡[6]，因此它是進(jìn)行層析反演的理論基礎(chǔ),本文采用射線追蹤理論正演模擬計(jì)算。
    本文結(jié)合射線追蹤方法中的試射法和微變網(wǎng)格法，建立復(fù)雜介質(zhì)模型以進(jìn)行微震正演模擬。真實(shí)的地質(zhì)情況是紛繁復(fù)雜的，對(duì)于這種界面變化復(fù)雜并且速度梯度變化較劇烈的模型，可通過(guò)將復(fù)雜介質(zhì)模型離散化成小網(wǎng)格，應(yīng)用微變網(wǎng)格法進(jìn)行射線追蹤的方法來(lái)求解。
    結(jié)合微震監(jiān)測(cè)作業(yè)特點(diǎn)，建立復(fù)雜介質(zhì)模型：
    (1)首先建立初始模型，輸入檢波器陣列和微震源位置坐標(biāo)，等間距網(wǎng)格化介質(zhì)模型，結(jié)合模型特點(diǎn)設(shè)定出射角范圍。
    (2)從微震源點(diǎn)所在網(wǎng)格出發(fā)，將各網(wǎng)格端點(diǎn)作為出射角,追蹤各網(wǎng)格出射邊，得出交點(diǎn)位置。
    (3)進(jìn)行網(wǎng)格入射邊微擾動(dòng)，獲得下個(gè)網(wǎng)格，如此反復(fù)，直到監(jiān)測(cè)井與射線相交。同網(wǎng)格時(shí)進(jìn)入步驟(4)；否則再次劃分原區(qū)間，轉(zhuǎn)入步驟(3)。
    (4)判斷監(jiān)測(cè)點(diǎn)與監(jiān)測(cè)井交點(diǎn)位置的關(guān)系，符合微震精度則轉(zhuǎn)入步驟(7)；否則重新設(shè)定區(qū)間，并利用二分法劃分區(qū)間，轉(zhuǎn)入步驟(2)。
    (5)記錄射線路徑及出射角，計(jì)算射線旅行時(shí)并記錄輸出。
    假定微震波在介質(zhì)離散化單元內(nèi)進(jìn)行傳播，已知離散化單元的位置、慢度和慢度梯度以及在射線入射口位置和此點(diǎn)方向，此問(wèn)題可歸結(jié)為求取其出口點(diǎn)位置和方向。其數(shù)學(xué)描述為：非均勻介質(zhì)模型，假設(shè)介質(zhì)離散化單元內(nèi)波速梯度為微常量，程函方程為：

1.3 布陣方案
    目前，微震監(jiān)測(cè)技術(shù)在震監(jiān)測(cè)傳感器布陣的可行性方面一直是一個(gè)難題。本文設(shè)計(jì)了幾種布陣方案，如圖1。為提高震源定位精度，合理布陣，并合理放置適量的拾震器，使用MATLAB編程工具實(shí)現(xiàn)對(duì)震源點(diǎn)的反演定位。用軟件計(jì)算出震源點(diǎn)的定位誤差并繪制出其坐標(biāo)誤差圖，比較呈十字形狀、米字形狀、四邊形狀、八邊形狀的傳感器布陣方案的震源定位誤差圖，觀察其震源定位精度，選出更合理的拾震器布陣方案，實(shí)現(xiàn)對(duì)微震監(jiān)測(cè)中拾震器布陣的評(píng)估。

圖1  4種拾震器布陣圖

2 計(jì)算機(jī)數(shù)值模擬實(shí)驗(yàn)及分析
2.1 建立復(fù)雜介質(zhì)數(shù)學(xué)模型
    設(shè)Rj和Si分別表示拾震器和微震源，S_i(x，y，z)為第i個(gè)微震源的實(shí)際空間坐標(biāo)，R_j(x，y，z)為第j個(gè)檢波器監(jiān)測(cè)點(diǎn)空間坐標(biāo)。微震監(jiān)測(cè)方式選擇井地聯(lián)合監(jiān)測(cè)[7]，采用復(fù)雜介質(zhì)模型，復(fù)雜地質(zhì)模型為1 km×1 km×1 km 的立方體，如圖2所示。

圖2  三維復(fù)雜地質(zhì)模型

    該三維模型由兩層均勻介質(zhì)與一層非均勻介質(zhì)組成，其中第一個(gè)界面為傾斜界面R1，第二個(gè)界面為水平界面R2：z=-800 m，其空間位置滿足關(guān)系14x-14y-5z=1 000；在這兩個(gè)界面之間嵌有一橢球體E，中心位于模型中心(0.5 km，0.5 km，-0.5 km)，其形狀滿足橢球面方程：。第一層界面以上的介質(zhì)I的縱波速度vp、橫波速度vs分別為2 km/s和1.2 km/s，除橢球體之外兩個(gè)界面之間介質(zhì)II的地震波速度分布為：
    v_p(x，y，z)=
    2 500+0.2(x-500)+0.1(y-500)-0.6(z+500)   m/s(7)
    v_s(x，y，z)=
    1 500+0.1(x-500)+0.05(y-500)-0.3(z+500)  m/s(8)
    第二界面以下介質(zhì)III為均勻介質(zhì)，其速度為：v_p=3.5 km/s，v_s=2.1 km/s；橢球體內(nèi)的介質(zhì)速度為：v_p=3 km/s；v_s=1.7 km/s。震源位于(56 m，80 m，-20 m)。對(duì)三維介質(zhì)模型，3個(gè)方向上均以10 m的間隔將模型剖分成100×100×100個(gè)立方體單元。在微機(jī)上利用微變網(wǎng)格射線追蹤算法計(jì)算初至波走時(shí)。
2.2 正演模擬
    在復(fù)雜介質(zhì)模型中，采用微變網(wǎng)格射線路徑追蹤算法進(jìn)行正演模擬，輸出微震波射線追蹤路徑，計(jì)算震源點(diǎn)至各檢波器的準(zhǔn)確初至到時(shí)，即各檢波器到時(shí)差序列△t_i。表1給出了各種布陣方案中震源點(diǎn)至各檢波器到時(shí)差序列結(jié)果(以第一道檢波器到時(shí)為基準(zhǔn)[8])。

2.3 加入隨機(jī)噪聲擾動(dòng)的PSO-BP反演
    檢波器檢波時(shí)受到噪聲的干擾，這些干擾波無(wú)法過(guò)濾得十分干凈，從而引入了時(shí)間誤差，因此定位誤差是由時(shí)間誤差而來(lái)的。由于這種時(shí)間誤差是隨機(jī)的，用隨機(jī)時(shí)間擾動(dòng)產(chǎn)生模擬定位擾動(dòng)，用其估計(jì)定位誤差，進(jìn)而估計(jì)和比較各種布陣方式的定位精度。本文采用正態(tài)分布作為模擬隨機(jī)噪音分布的函數(shù)形式。其中模擬過(guò)程中擾動(dòng)次數(shù)要盡可能多，而在在這個(gè)模型試驗(yàn)中發(fā)現(xiàn)，擾動(dòng)150次后誤差不再繼續(xù)減小。因此，采用標(biāo)準(zhǔn)偏差為1 ms的隨機(jī)正態(tài)分布模擬時(shí)間擾動(dòng)，噪聲擾動(dòng)150次后加在走時(shí)上進(jìn)行反演。
    反演中目標(biāo)函數(shù)為非線性回歸模型(走時(shí)方程)：
   

    由拾震器坐標(biāo)(Xi，Yi，Zi)得到待估計(jì)的未知參數(shù)(X₀，Y₀，Z₀，V_p，T₀)，使N個(gè)走時(shí)方程殘差平方和達(dá)到最小。
    通過(guò)微變網(wǎng)格射線追蹤正演程序得出了各檢波器監(jiān)測(cè)點(diǎn)理論上的到時(shí)差序列，將它們構(gòu)成學(xué)習(xí)樣本，進(jìn)而進(jìn)行PSO算法，從而優(yōu)化BP神經(jīng)網(wǎng)絡(luò)；在優(yōu)化后得出的最優(yōu)權(quán)值和閾值分布約束下進(jìn)行BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)，輸出結(jié)果即為所求微震源點(diǎn)空間坐標(biāo)。結(jié)合微震源反演模型，采用單隱層前饋網(wǎng)絡(luò)，三層拓?fù)浣Y(jié)構(gòu)分別為輸入層、隱含層和輸出層。輸入向量(到時(shí)差序列)個(gè)數(shù)n=20，輸出向量(微震源點(diǎn)坐標(biāo))個(gè)數(shù)m=3。隨機(jī)產(chǎn)生N組介于[-1，1]之間的權(quán)值和閾值分布。根據(jù)Kolmogorov定理，隱含層節(jié)點(diǎn)數(shù)理論上為2n+1個(gè)，即最大隱節(jié)點(diǎn)數(shù)Hmax=41，最小隱節(jié)點(diǎn)數(shù)Hmin=3，神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)為 20—41—3結(jié)構(gòu)，共有20×41+41×3=943個(gè)權(quán)值，41+3=44個(gè)閾值。BP網(wǎng)絡(luò)預(yù)測(cè)結(jié)果即為所求微震源空間坐標(biāo)，4種傳感器方案分別得出的反演定位結(jié)果與真實(shí)位置的標(biāo)準(zhǔn)差結(jié)果如表2所示。

    通過(guò)表2可以看出PSO-BP混合反演算法在最優(yōu)權(quán)值和閾值參數(shù)約束條件下的預(yù)測(cè)結(jié)果，即反演定位坐標(biāo)(x，y，z)及定位誤差?駐?啄；由誤差值分布情況可以看出，網(wǎng)絡(luò)預(yù)測(cè)誤差取值分布在范圍0～0.6之間，符合微震源反演定位誤差要求，加入隨機(jī)擾動(dòng)噪聲，反演預(yù)測(cè)結(jié)果可靠。
2.4 布陣方案仿真分析
    拾震器呈十字、米字、四邊形和八邊形布陣方案的震源定位誤差圖如圖3所示，圖中分別示出震源點(diǎn)在不同深度的定位誤差，深度按100 m步進(jìn)。

圖3  震源定位誤差圖

    由圖3可以發(fā)現(xiàn)，十字布陣在0～400 m范圍內(nèi)定位誤差小于1，米字布陣在0～700 m范圍內(nèi)定位誤差小于1，在700～1 300 m范圍內(nèi)米字布陣和四邊布陣定位誤差小于2，在1 300 m～1 800 m范圍內(nèi)四邊布陣定位誤差小于1，在1 800 m～2 450 m范圍內(nèi)八邊布陣定位誤差小于1。
2.5 拾震器數(shù)量與定位精度關(guān)系的仿真分析
    定位誤差隨拾震器數(shù)增加的變化曲線如圖4所示。圖4中隨拾震器個(gè)數(shù)的增加，定位誤差不斷減小，但當(dāng)拾震器增加到一定數(shù)量時(shí)，定位誤差趨于平緩，反應(yīng)了增加拾震器對(duì)定位誤差減小的效率問(wèn)題。

圖4  定位誤差隨拾震器數(shù)增加的變化曲線

    仿真實(shí)驗(yàn)中按4種布陣方式布陣，從5個(gè)拾震器逐個(gè)增加到90個(gè)，得出定位誤差隨拾震器數(shù)增加的變化曲線，縱坐標(biāo)是E(i)誤差,橫坐標(biāo)為拾震器個(gè)數(shù)。從實(shí)驗(yàn)結(jié)果可以看出，當(dāng)按十字方式布陣的拾震器數(shù)量從5增加到70個(gè)時(shí)，定位誤差成倍的減小，通過(guò)增加拾震器數(shù)量來(lái)提高定位精度的效果明顯；但是當(dāng)數(shù)量達(dá)到70個(gè)之后趨于平緩，減少誤差的效率比較低，這時(shí)通過(guò)增加拾震器數(shù)量來(lái)提高精度的效果非常不明顯。出于成本考慮，拾震器數(shù)量采用50～70個(gè)的性價(jià)比較高。按米字方式布陣的拾震器數(shù)量增加到70個(gè)時(shí)，定位誤差趨于平緩；按四邊方式布陣的拾震器數(shù)量增加到66個(gè)時(shí)，定位誤差趨于平緩；按八邊方式布陣的拾震器數(shù)量增加到57個(gè)時(shí)，定位誤差趨于平緩。
3 結(jié)束語(yǔ)
    文中4種拾震器布陣方案的仿真比較結(jié)果和定位精度比較的仿真結(jié)果表明，利用微變網(wǎng)格射線追蹤法建立微震地質(zhì)模型，對(duì)不同拾震器陣列方案的模型正演模擬研究得出的正演數(shù)據(jù)，在施加隨機(jī)擾動(dòng)的情況下，利用PSO-BP算法反演推算震源位置，最后通過(guò)分析對(duì)比震源位置的誤差來(lái)證明各個(gè)布陣方案的特點(diǎn)，這可以十分有效地解決微震監(jiān)測(cè)系統(tǒng)中拾震器如何布陣的問(wèn)題。同時(shí)，根據(jù)得到的震源定位誤差圖分析，可以在拾震器利用率最高的情況下，選出準(zhǔn)確、合理、有效的拾震器布陣方案。本研究對(duì)礦區(qū)微震定位監(jiān)測(cè)系統(tǒng)中的拾震器布陣具有指導(dǎo)意義。
參考文獻(xiàn)
[1] 萬(wàn)聚金，陳仲東.淺析煤礦超層越界開采的預(yù)防與治理[J].國(guó)土資源導(dǎo)刊，2012(5)：41-42.
[2] 平健，李仕雄，陳虹，等.微震定位原理與實(shí)現(xiàn)[J].金屬礦山，2010(1)：1-2.
[3] MENDECKI A J.Seismic monitoring in mines[M].London：Chapman and Hall ss，2010：12-15.
[4] KIJKO A，SCIOCATTI M.Optimal spatial distribution of seismic stations in mines[J].Int J.Rock Mech.and Mining Sci.and Geomechanics Abstracts，2012，32(6)：607-615.
[5] 牟宗龍，竇林名，鞏思園，等.礦井SOS微震監(jiān)測(cè)網(wǎng)絡(luò)優(yōu)化設(shè)計(jì)及震源定位誤差數(shù)值分析[J].煤礦開采，2009，14(3)：2-4.
[6] 張釙，劉洪，李幼銘.射線追蹤方法的發(fā)展現(xiàn)狀[J].地球物理學(xué)進(jìn)展，2010，15(1)：36-45.
[7] Song Fuxian，TOKS?魻Z M N.Full-waveform based complete moment tensorInversion and source parameter estimation from downhole microseismic data for hydrofracture monitoring[J].Geophysics，2011，40(4)：13-16.
[8] 楊曉東.基于射線理論二維正演模型的微地震震源反演[D].北京：中國(guó)石油大學(xué)，2013.