趙揚(yáng)鋒,趙 偉,王進(jìn)銘,潘一山,3

(1.遼寧工程技術(shù)大學(xué) 力學(xué)與工程學(xué)院,遼寧 阜新 123000;2.遼寧省礦山沉陷災(zāi)害防治重點(diǎn)實驗室,遼寧 阜新 123000;3.遼寧大學(xué) 災(zāi)害巖體力學(xué)研究所,遼寧 沈陽 110036)

0 引言

隨著我國國民經(jīng)濟(jì)的高速發(fā)展,礦山資源開發(fā)由淺層逐漸轉(zhuǎn)向深部開采,礦山災(zāi)害日益頻發(fā),礦山微震震源定位反演作為微震監(jiān)測系統(tǒng)的核心技術(shù),對礦山災(zāi)害的預(yù)警有重要意義。

微震事件定位方法中,有基于射線追蹤原理的定位方法、線性定位法及非線性定位法。射線追蹤技術(shù)用于微震震源定位主要是Julian等[1]提出2點(diǎn)間射線追蹤方法。Nakanishi等[2]將圖形網(wǎng)絡(luò)理論引入到射線路徑計算中,提出最短路徑射線追蹤法。在此基礎(chǔ),上王輝等[3]提出將單純形法與最短射線路徑追蹤法相結(jié)合的微震震源定位法,利用最短射線路徑追蹤法計算檢波器到潛在震源點(diǎn)的反向射線路徑,再在單純形法計算得到的初始解區(qū)間內(nèi)求節(jié)點(diǎn)走時-到時差之和最小值,該最小值對應(yīng)的節(jié)點(diǎn)即為震源點(diǎn)。Peng等[4]提出基于Log-Cosh的震源定位的目標(biāo)函數(shù),為提高定位準(zhǔn)確度和穩(wěn)定性剔除掉較遠(yuǎn)距離的傳感器P波到時數(shù)據(jù),提高定位準(zhǔn)確度及穩(wěn)定性。線性定位法是針對微震參數(shù)數(shù)據(jù)未知量建立線性方程組的數(shù)學(xué)求解方法,在線性定位法中應(yīng)用較廣的是Geiger法,Geiger法是Geiger[5]將定位非線性問題轉(zhuǎn)化為線性問題求解,利用震源位置坐標(biāo)和發(fā)震時刻建立線性方程組,對方程組使用一階泰勒展開后用最小二乘法求解參數(shù)方程得到震源位置。Waldhause等[6]提出雙差定位法(HYPODD),采用絕對走時和互相關(guān)P波和S波差分走時,當(dāng)所有觀測到的事件與監(jiān)測站臺聯(lián)系在一起時,觀測到時和理論傳播時差之間的雙差最小,通過迭代調(diào)整互相關(guān)震源之間的向量差,得到最小二乘解。線性定位法計算效率高,但可能會陷入局部最優(yōu)解,因此有學(xué)者提出非線性定位法。

非線性定位法可以避免陷入局部最優(yōu)解,其主要是建立震源位置的目標(biāo)函數(shù),使用不同的算法求解目標(biāo)函數(shù)最小值,得出震源位置坐標(biāo)。典型的非線性定位法有網(wǎng)格搜索法[7]、牛頓法[8]、梯度下降法[9]等。其中,網(wǎng)格搜索法引入較復(fù)雜的速度模型,定位精度較高,但網(wǎng)格劃分的疏密程度對計算速度和效率有很大的影響。梯度下降法是一階優(yōu)化算法,初始值的選取對定位精度有極大影響。牛頓法是二階最優(yōu)化算法,迭代次數(shù)少,若初始點(diǎn)選取合理,算法很快收斂。

除了提高算法本身的定位準(zhǔn)確度之外,有學(xué)者將算法與其他算法相結(jié)合進(jìn)行聯(lián)合定位,提高算法定位的準(zhǔn)確度及穩(wěn)定性。林峰等[10]提出線性定位和Geiger定位相結(jié)合的聯(lián)合定位法。用線性定位進(jìn)行初步定位,再以線性定位解作為Geiger定位法的迭代初值進(jìn)行求解。姜天琪等[11]提出1種將網(wǎng)格搜索法與牛頓法相結(jié)合的震源定位算法,利用網(wǎng)格搜索法為牛頓法提供迭代初值進(jìn)行震源定位,能夠有效避免由于迭代初值選取不當(dāng)造成定位失敗,提高了牛頓法定位準(zhǔn)確度和穩(wěn)定性。此外,微震監(jiān)測臺網(wǎng)的布置對定位精度有較大影響,因此劉曉明等[12]提出礦山微震監(jiān)測臺網(wǎng)的綜合評價指標(biāo),為監(jiān)測臺網(wǎng)的布置提供一定的理論支持。

多數(shù)震源定位方法需要預(yù)先給定速度模型,速度模型精度直接影響定位結(jié)果準(zhǔn)確性。針對速度模型對定位精度造成的不利影響,本文采用Adam算法[13]對震源進(jìn)行反演定位,通過設(shè)置不同的P波波速來驗證波速對定位結(jié)果的影響,結(jié)果顯示該算法能夠有效避免測速誤差對定位精度造成不利影響。此外,為提升算法定位的穩(wěn)定性,將線性方程法求得的初次定位解作為Adam算法的初始迭代值。

1 線性方程-Adam算法原理以及實現(xiàn)步驟

對礦山微震事件進(jìn)行定位,利用微震波走時進(jìn)行定位反演。在井下合理布置監(jiān)測臺站形成監(jiān)測結(jié)構(gòu),通過各個監(jiān)測臺站檢波器以及被激發(fā)微震信號的P波到時,建立走時方程如式(1)所示:

(1)

式中:ti是礦震微震信號P波到時時刻,s;t0是發(fā)震時刻,s;(xi,yi,zi)是第i個站臺坐標(biāo),m;(x0,y0,z0)是震源坐標(biāo),m;V是P波波速,m/s。

根據(jù)式(1)建立目標(biāo)函數(shù),如式(2)所示:

(2)

式中:n為站臺傳感器數(shù)量。

1.1 線性方程定位法的定位原理

線性方程定位法忽略了復(fù)雜的速度模型,用對確定區(qū)域震源進(jìn)行簡化的方法,從未知的震中(x0,y0,z0)到站臺檢波器(xi,yi,zi)(i=1,…,n)之間距離的斜直線來計算微震走時T,如式(3)所示:

(3)

式中:T為微震信號P波計算走時,s。

將式(3)進(jìn)行變換,如式(4)所示:

(4)

其中震源參數(shù)(x0,y0,z0,t0)是以下n-1個線性方程組的最小二乘解,如式(5)所示:

A·θ=r

(5)

式中:θ為震源參數(shù)向量(x0,y0,z0,t0)T;A是震源參數(shù)的系數(shù)矩陣;r是作差后的向量。

A和r如式(6)～(7)所示:

(6)

(7)

該方法至少需要5個站臺的P波到時,無需反復(fù)迭代,在本文中用于Adam算法的初始迭代值。

1.2 Adam算法原理

Adam算法由梯度下降算法改進(jìn)得到,相比于標(biāo)準(zhǔn)的梯度下降算法,Adam算法引入指數(shù)加權(quán)移動平均值對梯度值進(jìn)行估計同時又增加了偏置修正,修正從初始化原點(diǎn)一階矩和非中心的二階矩估計,使迭代更快速收斂,同時避免陷入局部極小值。

假設(shè)有目標(biāo)函數(shù)f(x1,x2,x3,…,xn),其中(x1,x2,x3,…,xn)是待更新參數(shù),Adam算法的目的是求解該目標(biāo)函數(shù)最小值;f1(x),f2(x)…fT(x)表示在后續(xù)的時間步1,2…T上的迭代更新。用gt=?xft(x)表示目標(biāo)函數(shù)f(x)的梯度。標(biāo)準(zhǔn)梯度下降算法按照式(8)來更新參數(shù):

xt=xt-1-αgt

(8)

式中:xt為待更新參數(shù)向量;α為步長因子。

Adam算法引入指數(shù)加權(quán)移動平均值,如式(9)所示:

vt=βvt-1+(1-β)gt

(9)

式中:vt為指數(shù)移動平均值,初始化值為0;β∈[0,1]是控制指數(shù)衰減率的超參數(shù)。

當(dāng)β取值接近0時,歷史信息對當(dāng)前梯度影響較小,這導(dǎo)致優(yōu)化曲線中包含更多的噪聲,會導(dǎo)致算法震蕩或者不穩(wěn)定的更新;當(dāng)β取值接近1時,歷史信息對當(dāng)前梯度影響較大,算法優(yōu)化更加穩(wěn)定,但收斂速度較慢,適用于包含噪聲較大的數(shù)據(jù)集。

在式(9)中計算指數(shù)加權(quán)移動平均值時需初始化v0=0,這導(dǎo)致第一項為零,使v1的估計值偏于零。因此在計算指數(shù)加權(quán)移動平均值時,為使得結(jié)果更加準(zhǔn)確,引入偏差修正,如式(10)所示:

(10)

Adam算法的更新規(guī)則如下:

1)初始化一階矩、二階矩變量s=0,u=0,以及矩估計的指數(shù)衰減率,β1,β2∈[0,1]。

2)根據(jù)式(9)計算梯度的有偏一階矩估計和梯度的二階矩估計,如式(11)～(12)所示:

st=β1st-1+(1-β1)gt

(11)

式中:st為梯度的一階矩估計。

(12)

式中:ut為梯度的二階矩估計。

考慮到礦山開采環(huán)境的復(fù)雜性,β1和β2的選取按照Kingma[13]設(shè)置的默認(rèn)值進(jìn)行計算,即β1=0.9,β2=0.99。

3)根據(jù)式(10)修正一階和二階矩偏差,如式(13)～(14)所示:

(13)

(14)

4)基于梯度下降更新參數(shù),如式(15)所示:

(15)

式中:δ是用于數(shù)值穩(wěn)定的小常數(shù),取10-8。

1.3 線性方程-Adam算法用于震源定位步驟

由式(2)給定初值(x1,y1,z1,t1)進(jìn)行迭代使目標(biāo)函數(shù)達(dá)到最小值,得到震源空間位置。根據(jù)式(2)將目標(biāo)函數(shù)做如下變化,如式(16)所示:

(16)

線性方程-Adam算法用于微震事件定位時的具體實現(xiàn)步驟如下:

1)首先讀取各個檢波器的位置數(shù)據(jù)及相應(yīng)的P波到時數(shù)據(jù)。

3)設(shè)定Adam算法的步長因子α,再將步驟2)計算得到的初步定位解作為Adam算法的初始迭代值進(jìn)行求解計算。

5)在每次的迭代循環(huán)內(nèi)設(shè)置迭代終止條件,即迭代到某次時目標(biāo)函數(shù)的值不再下降時跳出循環(huán),輸出迭代后的震源位置坐標(biāo)。

6)重復(fù)步驟4)～5),直至達(dá)到迭代終止條件。

線性方程-Adam迭代法的流程圖如圖1所示。

2 線性方程-Adam算法模擬實驗及結(jié)果

2.1 線性方程-Adam算法模擬實驗方案

線性方程-Adam算法模擬實驗如圖2所示,設(shè)置的監(jiān)測區(qū)域范圍是2 000 m×2 000 m×1 500 m,圖2中正方體表示檢波器,球體表示微震事件,1～4是位于監(jiān)測區(qū)域內(nèi)的震源,5～6是監(jiān)測區(qū)域外部的震源,各點(diǎn)的坐標(biāo)如表1～2所示。

表1 仿真模型各檢波器坐標(biāo)Table 1 The coordinates of each detector in the simulation model 單位:m

表2 仿真模型各震源坐標(biāo)Table 2 The source coordinates of the simulation model 單位:m

圖2 檢波器與微震事件分布Fig.2 Detector and microseismic event distribution

設(shè)現(xiàn)在微震P波平均波速為4 500 m/s,發(fā)震時刻0 s,按各個震源到檢波器的距離計算得到P波到時,如表3所示。

表3 仿真模型P波初至到時Table 3 Simulation model P wave first arrival time 單位:ms

在實際工程中,地下煤巖的各向異性使P波在傳播過程中的波速是不斷發(fā)生變化的,測得的波速會存在誤差,在一些優(yōu)化算法中,如網(wǎng)格搜索法、模擬退火法等,速度參數(shù)對定位結(jié)果有顯著影響。因此本文將介質(zhì)的各向異性體現(xiàn)在速度參數(shù)上,具體做法參考李健等[14]關(guān)于速度各向異性模型處理,在平均速度的基礎(chǔ)上隨機(jī)分配3%和5%的速度波動范圍,表4為在速度各向異性條件下的P波初至?xí)r間。

表4 速度各向異性下的P波初至?xí)r間Table 4 P wave first arrival time under velocity anisotropy 單位:ms

2.2 定位結(jié)果分析

2.2.1 步長因子對定位結(jié)果影響

Adam算法在對目標(biāo)函數(shù)優(yōu)化時,初始迭代值選定后,利用空間距離公式計算得到定位坐標(biāo)與震源坐標(biāo)之間的總誤差,如式(17)所示,不同步長因子定位結(jié)果見表5。步長因子對定位結(jié)果的影響由算法本身的性質(zhì)和目標(biāo)函數(shù)決定,與震源位置無關(guān),因此表5選擇微震事件1計算得到的定位誤差,其他微震事件同樣適用于此規(guī)律。

表5 步長因子對定位結(jié)果的影響Table 5 The influence of step factor on the positioning results 單位:m

(17)

式中:d為定位總誤差,m;(x,y,z)為定位坐標(biāo),m;(x0,y0,z0)為震源坐標(biāo),m。

由表5可知,隨著步長因子從0.1逐漸減小到0.001時,定位誤差在減少。但當(dāng)步長因子從0.001減小到0.000 1時,定位誤差在增加。這是由于單步迭代算法在迭代計算時,計算機(jī)字長有限,每一步運(yùn)算均需四舍五入,因此迭代的每一步會產(chǎn)生舍入誤差。從收斂性結(jié)果來看,步長因子越小即步長越小,結(jié)果越準(zhǔn)確;但從舍入誤差來看,計算步數(shù)的增多導(dǎo)致舍入誤差積累,最終使得計算結(jié)果失真,因此,選擇步長因子10-3進(jìn)行迭代。

2.2.2 迭代初值對定位結(jié)果的影響

對于大多數(shù)優(yōu)化算法而言,迭代初值的選取會直接影響算法能否尋找到全局最優(yōu)解,進(jìn)而影響到定位結(jié)果。同樣,Adam迭代算法的結(jié)果也會受到迭代初值的影響。在速度參數(shù)確定后,不同的迭代初值對定位結(jié)果的影響如表6所示,表6中每個微震事件的第1個迭代初值為線性方程法的定位解,其余2個迭代初值為隨機(jī)初始值。

表6 不同迭代初值對定位結(jié)果的影響Table 6 The influence of different iterative initial values on the positioning results

觀察表6中每個微震事件不同的迭代初值的總誤差,選擇線性定位法的定位結(jié)果作為Adam算法的初始迭代值時,每個事件的總誤差均小于1 m,由隨機(jī)初始化得到的迭代初值的定位總誤差均在100 m以上,甚至達(dá)到300 m。因此每個微震事件的定位結(jié)果,迭代初值的不同對定位結(jié)果影響較大,而選取線性定位法的結(jié)果作為迭代初值,其定位準(zhǔn)確度較高。

2.2.3 速度參數(shù)對定位結(jié)果的影響

為驗證和比較線性方程-Adam定位法的準(zhǔn)確性以及穩(wěn)定性,本文分別選取線性方程-Adam算法、牛頓迭代法、模擬退火法、Geiger定位法4種方法對表3～4數(shù)據(jù)進(jìn)行反演定位。此外,為更貼合實際的工作環(huán)境,加入測速誤差對總定位結(jié)果的影響。

在速度各向異性的干擾下,測速誤差分別為0,-3%,3%,即速度分別是4 500,4 365,4 635 m/s,4種方法的定位誤差如圖3所示。圖3中編號見表7。

表7 圖例編號Table 7 Legend number

圖3 速度各向異性和測速誤差影響下的定位誤差對比Fig.3 Comparison of positioning errors under influence of velocity anisotropy and velocity measurement error

1)速度各向異性對定位結(jié)果的影響。圖3(a)中的定位準(zhǔn)確度折線展示了在無測速誤差情況下,速度各向異性對定位結(jié)果的影響。圖3中的6條折線代表6個微震事件在速度各向異性條件下,使用不同定位方法得到的定位結(jié)果。圖3中的橫坐標(biāo)表示各種方法對應(yīng)的速度各向異性范圍,以微震事件1為例,利用牛頓迭代法定位,D、E、F分別表示牛頓迭代法在速度各向異性為0、3%和5%時,微震事件1的定位總誤差。

各向異性范圍為0時,圖中的附圖結(jié)果顯示,線性方程-Adam法、牛頓迭代法、模擬退火法以及Geiger法的定位結(jié)果相對準(zhǔn)確。

當(dāng)各向異性范圍為3%時,圖3(a)中的附圖結(jié)果顯示,線性方程-Adam法的定位結(jié)果沒有受到影響,而其他方法均受到不同程度的影響。觀察橫坐標(biāo)E代表的牛頓迭代法,在各向異性范圍為3%時,定位總誤差明顯增加。模擬退火法在各向異性范圍為3%時,對微震事件的定位誤差雖有增加,但誤差在可接受范圍內(nèi)。而Geiger法在速度各向異性范圍的影響下,對微震事件的定位誤差最大。

各向異性范圍為5%時,線性方程-Adam算法對各微震事件仍然能準(zhǔn)確定位,而其余3種方法對各微震事件定位結(jié)果誤差均較大。綜合上述定位方法與線性方程-Adam算法的對比,在速度各向異性的條件下,只有線性方程-Adam法能夠?qū)崿F(xiàn)準(zhǔn)確定位。

2)測速誤差對定位結(jié)果的影響。速度各向異性范圍為0的條件下,對比分析測速誤差對定位結(jié)果的影響。觀察對比圖3(a)～圖3(c)中橫坐標(biāo)A、D、G、J對應(yīng)的定位總誤差,對比圖3(a)～圖3(c)的附圖橫坐標(biāo)A,線性方程-Adam算法對各微震事件定位的結(jié)果顯示,在各向異性范圍為0時,不同的測速誤差對定位結(jié)果不會產(chǎn)生影響。對比圖3(a)～圖3(c)的橫坐標(biāo)D,結(jié)果顯示隨測速誤差的增加,牛頓迭代法對各微震事件定位的誤差在增大。對比圖3(a)～圖3(c)的橫坐標(biāo)G,結(jié)果顯示隨測速誤差的增加,模擬退火法對各微震事件定位的誤差相差較小,因此測速誤差對模擬退火的定位影響較小,但模擬退火法的定位準(zhǔn)確度和穩(wěn)定性要稍劣于線性方程-Adam算法。對比圖3(a)～圖3(c)的橫坐標(biāo)J,結(jié)果顯示隨著測速誤差的增加,Geiger法對各微震事件的定位誤差在增大,測速誤差對Geiger法的影響比對牛頓迭代法的影響還要大。線性方程-Adam算法不會受測速誤差的影響,測速誤差對模擬退火的影響較小,對牛頓迭代法和Geiger法的影響較大。

3)速度各向異性和測速誤差的綜合影響。圖3(a)～圖3(c)中橫坐標(biāo)B,C,E,F,H,I,K,L表示當(dāng)速度各向異性范圍達(dá)到3%和5%時,各定位方法在測速誤差分別為-3%和3%時的定位誤差,可以觀察到線性方程-Adam算法的定位仍然沒有發(fā)生變化,而牛頓迭代法的定位誤差有顯著的增大,其定位結(jié)果不理想,定位穩(wěn)定性變差,模擬退火法的定位誤差沒有太大的變化,Geiger法的定位誤差有明顯的增大。

綜合上述分析,Adam算法由線性方程確定迭代初值后,對微震事件的定位結(jié)果不受速度各向異性和測速誤差的影響。這是由目標(biāo)函數(shù)和算法參數(shù)更新策略決定,為使式(2)的目標(biāo)函數(shù)達(dá)到最小值,可以認(rèn)為將P波波速和P波走時作為1個整體參數(shù),而震源位置參數(shù)作為另外1個整體參數(shù)。Adam算法在搜索目標(biāo)函數(shù)最小值更新參數(shù)時看作更新2個整體參數(shù)即震源位置參數(shù)和由P波波速及P波走時組成的整體參數(shù),若調(diào)整P波波速,為保證目標(biāo)函數(shù)的最小值,算法更新P波走時抵消P波波速的變化導(dǎo)致目標(biāo)函數(shù)最小值的變化,最終表現(xiàn)為P波波速不會影響定位精確度。但這會導(dǎo)致P波走時誤差增大,而檢波器收集到的P波到時是固定的,因此最后的誤差就會累積到發(fā)震時刻。本質(zhì)上是以發(fā)震時刻的準(zhǔn)確度為代價。本文研究的是震后反演定位,因此對過去的發(fā)震時刻沒有太高的要求。

3 工程實例驗證

本文利用林峰等[10]在柿竹園礦區(qū)爆破實驗以及劉超等[15]人工爆破實驗。在柿竹園礦區(qū)內(nèi)測得的爆破坐標(biāo)為(8 732.70,6 570.60,511.30),此次爆破實驗共計8個傳感器接收到微震信號,實驗測得的P波波速為5 000 m/s,8個臺站位置信息以及P波到時如表8所示。

表8 柿竹園臺站坐標(biāo)及P波初至到時Table 8 The coordinates of Shizhuyuan station and the first arrival time of P wave

劉超等[15]進(jìn)行的人工爆破實驗使用5個檢波器收集爆破信號,在監(jiān)測區(qū)域煤巖體內(nèi)最優(yōu)化的P波傳播速度為2 572.5 m/s,測得的爆破坐標(biāo)(4 290,7 255,-676),臺站坐標(biāo)和P波到時如表9所示。

表9 人工爆破臺站坐標(biāo)及P波初至到時Table 9 The coordinates of artificial blasting station and the first arrival time of P wave

將相應(yīng)的定位參數(shù)帶入線性方程-Adam算法進(jìn)行定位,與模擬退火-單純形法、Geiger法、牛頓迭代法定位結(jié)果進(jìn)行比較見表10和表11。2類爆破實驗的定位結(jié)果顯示,線性方程-Adam定位法相比于其他定位方法,定位準(zhǔn)確度有一定的提高。模擬退火法每次定位都能靠近真實震源空間位置,但因其采用空間隨機(jī)搜索模式,每次定位結(jié)果都不同,導(dǎo)致定位結(jié)果不穩(wěn)定;Geiger法的定位結(jié)果受速度參數(shù)影響較大;牛頓法定位需要求二階導(dǎo)數(shù)的Hessian矩陣增加了計算量;線性方程-Adam算法在迭代過程中,確定迭代初值后定位結(jié)果不受速度參數(shù)的影響,且在相同的初始條件下,每次的定位結(jié)果都是相同的。

表11 人工爆破實驗定位結(jié)果對比分析Table 11 Comparative analysis of positioning results in artificial blasting test 單位:m

4 結(jié)論

1)線性方程-Adam震源定位法,利用線性方程的初次定位解作為Adam算法的迭代初值提升了Adam算法定位的穩(wěn)定性,并有效避免迭代初值選擇不當(dāng)對定位準(zhǔn)確度的影響。

2)利用柿竹園礦爆破數(shù)據(jù)對算法進(jìn)行驗證,結(jié)果顯示線性方程-Adam定位法相比于牛頓迭代法無需計算二階導(dǎo)數(shù),節(jié)省計算資源,且定位精度提高了54%;相比于模擬退火法-單純形法定位精度提高了67%;相比于Geiger法定位精度提高了82%。

3)利用劉超等人工爆破實驗對算法驗證,結(jié)果顯示線性方程-Adam定位法相比于牛頓迭代法定位精度提升了33%,相比于模擬退火-單純形法定位精度提升了66%,相比于Geiger法定位精度提升了74%。

4)線性方程-Adam定位法相比于其他算法,無需預(yù)先獲取復(fù)雜的速度模型,迭代初值確定后,其定位結(jié)果不會受測速誤差影響。