基于線性方程-Adam算法的微震震源定位研究*

趙揚鋒,趙 偉,王進銘,潘一山,3

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

0 引言

隨著我國國民經濟的高速發展,礦山資源開發由淺層逐漸轉向深部開采,礦山災害日益頻發,礦山微震震源定位反演作為微震監測系統的核心技術,對礦山災害的預警有重要意義。

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

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

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

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

1 線性方程-Adam算法原理以及實現步驟

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

(1)

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

根據式(1)建立目標函數,如式(2)所示:

(2)

式中:n為站臺傳感器數量。

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

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

(3)

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

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

(4)

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

A·θ=r

(5)

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

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

(6)

(7)

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

1.2 Adam算法原理

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

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

xt=xt-1-αgt

(8)

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

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

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

(9)

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

當β取值接近0時,歷史信息對當前梯度影響較小,這導致優化曲線中包含更多的噪聲,會導致算法震蕩或者不穩定的更新;當β取值接近1時,歷史信息對當前梯度影響較大,算法優化更加穩定,但收斂速度較慢,適用于包含噪聲較大的數據集。

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

(10)

Adam算法的更新規則如下:

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

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

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

(11)

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

(12)

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

考慮到礦山開采環境的復雜性,β1和β2的選取按照Kingma[13]設置的默認值進行計算,即β1=0.9,β2=0.99。

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

(13)

(14)

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

(15)

式中:δ是用于數值穩定的小常數,取10-8。

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

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

(16)

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

1)首先讀取各個檢波器的位置數據及相應的P波到時數據。

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

5)在每次的迭代循環內設置迭代終止條件,即迭代到某次時目標函數的值不再下降時跳出循環,輸出迭代后的震源位置坐標。

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

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

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

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

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

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

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

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

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

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

在實際工程中,地下煤巖的各向異性使P波在傳播過程中的波速是不斷發生變化的,測得的波速會存在誤差,在一些優化算法中,如網格搜索法、模擬退火法等,速度參數對定位結果有顯著影響。因此本文將介質的各向異性體現在速度參數上,具體做法參考李健等[14]關于速度各向異性模型處理,在平均速度的基礎上隨機分配3%和5%的速度波動范圍,表4為在速度各向異性條件下的P波初至時間。

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

2.2 定位結果分析

2.2.1 步長因子對定位結果影響

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

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

(17)

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

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

2.2.2 迭代初值對定位結果的影響

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

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

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

2.2.3 速度參數對定位結果的影響

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

在速度各向異性的干擾下,測速誤差分別為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)速度各向異性對定位結果的影響。圖3(a)中的定位準確度折線展示了在無測速誤差情況下,速度各向異性對定位結果的影響。圖3中的6條折線代表6個微震事件在速度各向異性條件下,使用不同定位方法得到的定位結果。圖3中的橫坐標表示各種方法對應的速度各向異性范圍,以微震事件1為例,利用牛頓迭代法定位,D、E、F分別表示牛頓迭代法在速度各向異性為0、3%和5%時,微震事件1的定位總誤差。

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

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

各向異性范圍為5%時,線性方程-Adam算法對各微震事件仍然能準確定位,而其余3種方法對各微震事件定位結果誤差均較大。綜合上述定位方法與線性方程-Adam算法的對比,在速度各向異性的條件下,只有線性方程-Adam法能夠實現準確定位。

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

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

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

3 工程實例驗證

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

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

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

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

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

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

4 結論

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

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

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

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