基于到時時空分布的微震震源定位算法

李明虎，陳 衛(wèi)，蔡 珣

（1.山東大學(xué) 軟件學(xué)院，山東 濟南 250101；2.北京安科興業(yè)科技股份有限公司，北京 102200）

我國沖擊地壓煤礦較多，隨著煤礦開采強度和深度的增加，沖擊地壓礦井數(shù)量和沖擊危險程度顯著增加[1]。礦山微震監(jiān)測預(yù)警系統(tǒng)的使用對于保證安全生產(chǎn)、提高生產(chǎn)效率意義重大，是沖擊地壓災(zāi)害精準防控的關(guān)鍵。微震震源定位，即根據(jù)檢波器陣列采集的微震信號，判斷礦山微震發(fā)生的時空位置，是監(jiān)測預(yù)警工作中的核心問題之一。礦山微震的震相主要是體波，當前我國礦山大多使用單軸檢波器，監(jiān)測豎直方向的振動，拾取的主要是P 波。

微震震源定位算法起源于天然地震的震源定位算法研究，屬于被動地震技術(shù)。當前主流的微震定位算法有3 大類：基于到時的定位法，如經(jīng)典蓋革（Geiger）算法[2]、單純形算法[3]、雙差定位法[4]等；基于波形偏移定位法，如逆時成像法[5]、震源掃描算法[6]（SSA Source Scanning Algorithm）、聯(lián)合微震映射法[7]（CMM Coalescence Microseismic Mapping）等；基于三分量檢波器定位法，如單臺站法[8]、P 波偏振法[9]等?；诘綍r的定位算法提出較早，發(fā)展充分，計算簡單，定位精度較高，是當前生產(chǎn)中應(yīng)用的主流算法。

根據(jù)發(fā)展的驅(qū)動力，震源定位算法經(jīng)歷了4 個階段：①理論分析階段（1912—1970 年）：經(jīng)典蓋革法于1912 年提出，本質(zhì)上是求解超定方程組的最小二乘解，計算量巨大，到20 世紀70 年代計算機技術(shù)普及，求解蓋革法的數(shù)值分析算法才得以在工程上有效實現(xiàn)，這一階段的主要工作是初步的理論研究；②算法驅(qū)動階段（1970—2005 年）：20 世紀70 年代及隨后的30 多年中，一系列的算法進入實際應(yīng)用并獲得反饋和提升，如單純形算法[3]、聯(lián)合反演法[10]、主事件定位法[11]等，蓋革法的思想得到改進和發(fā)展，定位精度大幅提高，震源定位算法的核心概念和基本理念大多產(chǎn)生于這個階段；③計算驅(qū)動階段（2005—2018 年）：進入21 世紀，計算機、物聯(lián)網(wǎng)、數(shù)據(jù)分析等技術(shù)發(fā)展提升，天然地震和礦山微震工作積累了更多的數(shù)據(jù)數(shù)量，數(shù)據(jù)質(zhì)量也大幅提升，其間更多數(shù)據(jù)得到分析處理，更加復(fù)雜的算法得到實施應(yīng)用，與其他學(xué)科的交叉研究，如二次優(yōu)化法[12]、貝葉斯理論差分進化[13]、粒子群差分進化[14]等取得進展，這段時間是震源定位技術(shù)深入發(fā)展的階段；④數(shù)據(jù)驅(qū)動階段（2018 年至今）：2018 年以后，人工智能技術(shù)逐步應(yīng)用于地球物理學(xué)科的各個領(lǐng)域，基于大量數(shù)據(jù)的有監(jiān)督學(xué)習(xí)和無監(jiān)督學(xué)習(xí)已經(jīng)在自動到時拾取[15]，自動高分信號提取[16]等方面進行應(yīng)用，直接或間接提升地震和微震定位的精度，人工智能的發(fā)展與地球物理學(xué)的應(yīng)用交叉融合，提升和突破值得期待。綜上，首先簡要介紹蓋革法，說明傳統(tǒng)目標函數(shù)的全局最優(yōu)值無法有效指向?qū)嶋H震源點（或其附近）；然后通過到時的分布信息和潛在震源點的區(qū)域信息構(gòu)建合理的目標函數(shù)，形成新的D4DA 算法；最后在生產(chǎn)數(shù)據(jù)中對新算法進行驗證。

1 基于到時的震源定位算法

1.1 經(jīng)典蓋革法

基于到時的定位算法的基本思想是利用直達波初至和地層速度模型來反演震源位置和激發(fā)時刻[17]。到時類算法大多以蓋革法為基礎(chǔ)，基本原理如圖1。

圖1 到時不同定位法原理示意圖[18]Fig.1 Locating hypocenter by arrival time

S（x0，y0，z0，t0）和Ti（xi，yi，zi，ti）分別為實際震源和第i 個檢波器的時空位置，v 為地震波的波速。蓋革法的介紹文獻較多，如[2，8]，簡要描述如下：

1）進行到時提取。分析檢波器記錄信號，確定P波初至?xí)r刻。當前，到時拾取工作仍然需要一定程度的人工參與。

2）確定目標函數(shù)。計算潛在震源到各檢波器的距離di，結(jié)合波速v 計算地震波走時tgi，tgi與起震時間ts的和記為計算到時tci，即tci=ts+tgi，理想情況下，計算到時tci應(yīng)當與觀測到時toi相等。設(shè)所有檢波器計算到時向量tc與觀測到時向量to的差向量為a（記為到時偏差），即a=tc-to，使用a 的函數(shù)，如1-范數(shù)和2-范數(shù)。

3）優(yōu)化目標函數(shù)。選擇合適的反演策略找到最優(yōu)的起震時間和震源位置。一般采用數(shù)值分析的方法，針對所使用的觀測數(shù)據(jù)類型、物理方程及約束條件等因素，選擇最優(yōu)反演策略以避免收斂到局部極值，并且快速、準確地找到最優(yōu)解。

實際微震波傳播過程中，會受自然地質(zhì)、工程采動等因素的影響，傳播介質(zhì)并非均勻，更非各向同性，因此，實際微震波的速度場是非均勻且各向異性的，為了聚焦在問題的主要方面，這里采用均勻且各向同性的速度場。當前，完全依賴算法的定位結(jié)果精度一般不會很高，誤差多在20 m 量級或以上，會出現(xiàn)明顯違背常理的結(jié)果。生產(chǎn)中進行震源定位時大多需要人工干預(yù)，如選擇合適的監(jiān)測信號、進行適合的濾波、選擇合理的通道數(shù)或去除不合理結(jié)果等。提高定位工作的精度和智能化程度一直是震源定位工作的目標。

1.2 L1 和L2 目標函數(shù)

設(shè)起震時刻為ts，第i 個檢波器的觀測到時、走時、計算到時為toi、tgi、tci，有tci=ts+tgi。

定義L1目標函數(shù)（1-范數(shù)）為：

L2目標函數(shù)（2－范數(shù)）為：

式中：L1、L2為目標函數(shù)，ms；n 為檢波器個數(shù)；abs（·）為取絕對值；tci為計算到時，ms；toi為觀測到時，ms。

2 個目標函數(shù)性能優(yōu)良，在人工智能、模式識別等學(xué)科中應(yīng)用廣泛，基于到時的定位算法大多使用2 種目標函數(shù)或其變體。由二者的定義可知，L1對每個檢波器的偏重是一樣的，而L2更加偏重于誤差較大的檢波器。L2在定義域上連續(xù)可導(dǎo)，L1連續(xù)，但存在不可導(dǎo)的點。L2在優(yōu)化方法適應(yīng)性上優(yōu)于L1。

2 D4DA 算法的思路與實現(xiàn)

2.1 傳統(tǒng)震源定位算法中的偏離問題

蓋革法的基本假設(shè)是目標函數(shù)最優(yōu)值在實際震源點附近。然而，目標函數(shù)全局最優(yōu)值可能離實際震源點很遠，即使只是在理想情況下加上非常小的擾動，全局最優(yōu)值也可能以較高的概率對實際震源點有非常大的偏離；實驗表明偏離是普遍的，在實際生產(chǎn)量級的擾動下，偽震源點占比可能達到10%甚至40%以上。

在目標函數(shù)普遍偏離的情況下，優(yōu)化算法失去真正的目標，難以取得好的定位結(jié)果。需要技術(shù)人員結(jié)合現(xiàn)場情況進行經(jīng)驗判斷才能進行較好的震源定位。

2.2 目標函數(shù)Ls 的構(gòu)建

新目標函數(shù)構(gòu)建的思路：使用到時偏差的分布信息確定單個潛在震源點的基礎(chǔ)目標函數(shù)；使用潛在震源點鄰近的時空區(qū)域信息構(gòu)建最終目標函數(shù)。

具體構(gòu)建包含3 個步驟，①確定到時偏差的尺度；②計算到時偏差的似然函數(shù)作為單點目標函數(shù)；③使用區(qū)域信息計算目標函數(shù)。

2.2.1 確定到時偏差的尺度。

由ai＝tci－toi＝ts＋tgi－toi可知，到時偏差ai的方差可由式（3）進行計算。

式中：σi為到時偏差ai的標準差，ms；σs為起震時刻標準差，ms；σgi為走時標準差，ms；σoi為觀測到時標準差，ms。

檢波器的空間位置固定且到時拾取算法確定時，通過隨機變量S（x0，y0，z0，t0）計算σs、σgi、σoi，進而計算σi，需知道S 在4 維時空的分布，這需要精準的地質(zhì)、采動的模型進行演繹或大量的S 和觀測到時向量to的對應(yīng)數(shù)據(jù)進行歸納，二者都是無法做到的。

分析σi的使用需求和σs、σgi、σoi的特點來尋求σi的合理近似。各個檢波器到時偏差的標準差σi是為了讓不同檢波器通道上的到時偏差可以相互比較，σi的相對大小更有意義，各σi同步的縮放不會對后續(xù)目標函數(shù)的相對大小產(chǎn)生影響，也就不會改變目標函數(shù)最優(yōu)值點的坐標即震源；起震時刻對各檢波器是獨立的，從而各個檢波器的σs是相同的，這即是σs沒有腳標i 的原因；觀測到時的偏差來源于到時拾取的方法和檢波器的位置，到時拾取的方法與檢波器獨立可忽略，檢波器的位置對觀測到時標準差σoi和走時標準差σgi的影響是一樣的。綜合上述，在式（3）中忽略σs，把σoi替換為σgi，得σi≈2 σgi，進一步忽略系數(shù)2，得式（4）。

式中：σi為到時偏差的標準差，ms；i 為檢波器編號；σgi為走時的標準差，ms；std（·）為取標準差；tgi為走時，ms。

式（4）表明，如果求得隨機變量tgi，則到時偏差的尺度問題即告解決。理想情況下求tgi需要震源在三維空間的分布，實際計算時也是需要合理的簡化。檢波器臺網(wǎng)布設(shè)主要關(guān)注微震易發(fā)區(qū)域，危險區(qū)域周邊應(yīng)盡量在空間上被候選點均勻包圍[18]，這就提示，可以構(gòu)建以各檢波器的形心為中心的適當間距擴展的柵格，通過柵格點來構(gòu)建隨機變量tgi。這種方法另外的2 個優(yōu)點：①潛在震源位置受各檢波器的影響相似；②相對于邊緣點，形心附近的潛在震源位置變化時tgi會有合理的波動。

實際計算σi所使用的形心法如圖2。

圖2 形心法計算到時偏差的標準差Fig.2 Calculating standard deviation of arrival time deviations by the centroid-method

如圖2，以所有檢波器的形心為中心，以Δd=Δt×v（Δt，v 分別是采樣時間間隔和波速）為間隔，沿著3 軸6 方向分別對外擴展n 個步長，形成（2n＋1）3的網(wǎng)格，計算每個節(jié)點到各個檢波器的走時tgi，計算tgi的標準差σgi，依據(jù)式（4）將σgi當作σi來使用。可適當調(diào)整n 的值，使得各通道tgi的標準差的平均值與Δt×v 相當，實際應(yīng)用時，n 的值一般取5～25，對應(yīng)的空間立方區(qū)域的邊長大約是40～200 m，范圍過小則σi之間的比例波動劇烈，范圍過大則各σi趨于相等。

2.2.2 計算到時偏差的似然函數(shù)

對傳統(tǒng)定位算法分析發(fā)現(xiàn)，L2的表現(xiàn)好于L1，說明到時偏差大的分量應(yīng)當賦予更大的權(quán)重。對a=tc-to，假設(shè)a 各分量間相互獨立，且到時偏差的每個分量ai=tci-toi符合正態(tài)分布，則可用式（5）的似然函數(shù)描述到時偏差的不確定性，并記為Lp目標函數(shù)。

式中：Lp為似然函數(shù)即單點目標函數(shù)；p（·）為取概率密度；a 為到時偏差，ms；n 為檢波器個數(shù)；tci為計算到時，ms；toi為觀測到時，ms；K 為歸一化因子；ai為到時偏差a 的第i 個分量，ms；σi為到時偏差的標準差，ms。

相比于L1、L2目標函數(shù)，因為σi的引入，Lp目標函數(shù)考慮了偏差的尺度信息，不同檢波器之間的差異有了可比性，使得到時偏差的意義更加明確。

2.2.3 使用時空區(qū)域信息計算目標函數(shù)

對于實際震源點，如果時空坐標發(fā)生較小的擾動，目標函數(shù)會以相對和緩的方式進行減少。對于偽震源點，其時空坐標的目標函數(shù)只是“恰好”在某一個位置較小，如果時空坐標發(fā)生較小的擾動，那么目標函數(shù)將會發(fā)生比較大的變動。表明將某個時空坐標附近多個點的目標函數(shù)進行聚合，以聚合值作為總體目標函數(shù)將會使算法有更好的魯棒性。

具體地，Ls目標函數(shù)采用如下的方式進行構(gòu)建，設(shè)檢波器采樣間隔為Δt，以潛在震源點S（xq，yq，zq，tq）為中心，空間上以Δd＝Δt×v 為間隔沿3 軸6 向各外擴展n 個步長，時間上以Δt 為間隔向前后各擴展n 個步長，形成（2n＋1）4的4 維時空網(wǎng)格，對每個網(wǎng)格節(jié)點，計算其Lp目標函數(shù)。聚合權(quán)重可為等權(quán)相加，如式（6），也可取4 維高斯卷積核。

式中：Ls為4 維區(qū)域目標函數(shù)；q 為4 維區(qū)域內(nèi)網(wǎng)格節(jié)個數(shù)；Lpi為單點目標函數(shù)。

需要說明的是：①形心法確定偏離方差時使用3 維空間區(qū)域，且只需要計算1 次，Ls使用4 維時空區(qū)域，且需要對每個潛在震源點進行計算；②Lp目標函數(shù)比L2目標函數(shù)略復(fù)雜，但都是O（n）級別，Ls目標函數(shù)的復(fù)雜度為O（n×m4），其中m 為時空區(qū)域核的寬度，由于冪次較高，m 取值不可太大，當n 取2 時m 為5，Ls的計算量是Lp的625 倍；③在形心法和Ls計算時，都使用了Δd=Δt×v 計算空間柵格的間距，空間間距本質(zhì)上是由微震信號的采樣時間間隔Δt 確定的，微震波在此期間行進了Δt×v 的距離，即速度是時間與空間的紐帶。

2.3 D4DA 算法

D4DA 算法流程圖如圖3。

圖3 D4DA 算法流程圖Fig.3 Flowchart of D4DA algorithm

算法輸入包括各檢波器位置Ti及其觀測到時toi、波速v、搜索區(qū)域R。首先，用形心法計算到時偏差標準差σi作為各檢波器到時偏差的尺度；然后，按2.2 節(jié)“使用時空區(qū)域信息計算目標函數(shù)”方法確定潛在時空震源點組，按式（5）確定單點目標函數(shù)Lp并按式（6）匯總得到算法的目標函數(shù)Ls最后，使用擬牛頓法在潛在搜索區(qū)域R 中得到使得目標函數(shù)最優(yōu)值的時空坐標Sc，輸出Sc作為定位結(jié)果。

3 數(shù)值實驗與結(jié)果分析

3.1 數(shù)據(jù)設(shè)定

使用文獻［19］中第五章的空間數(shù)據(jù)，包括16 個檢波器和12 個炮點的位置坐標，劃定搜索空間為可以包含所有檢波器和炮點的、最小的長方體向3 軸6 向外側(cè)各擴展200 m 的長方體的范圍。選取這里的數(shù)據(jù)是因為其搜索空間較大，有較好的深度跨度；有12 組標定炮數(shù)據(jù)，16 個檢波器，數(shù)據(jù)充分，容易體現(xiàn)出統(tǒng)計特征，便于發(fā)現(xiàn)規(guī)律。檢波器和標定炮坐標可見原文獻。潛在震源點的搜索范圍見表1。

表1 仿真數(shù)據(jù)的搜索區(qū)域Table 1 Search area in simulation data

對于遍歷搜索，各坐標間隔取7 m，立方對角間隔為12.12 m 可以滿足絕大多數(shù)生產(chǎn)環(huán)境中的定位粒度要求。此時每次需要計算30 945 110 個候選震源點。使用Numpy 包進行計算，如果有16 個通道拾取了到時，一次定位需要約20 GB 內(nèi)存，在i7-8700K 處理器上運行約60 s。

到時數(shù)據(jù)的生成：取起震時刻ts為1 000 ms，取各向同性的均勻波速場且波速vp為3.75 m/ms，根據(jù)炮點和各檢波器的坐標計算距離di，得走時tgi=di/vp，由式（7）生成觀測到時。

式中：toi為觀測到時，ms；i 為檢波器編號；ts為起 震 時 刻，ms；tgi為 走 時，ms；r 為 信 號 干 擾 強 度；norm（0，1）為標準正態(tài)分布發(fā)生器。

理論觀測到時為起震時刻加走時，實際觀測到時為理論觀測到時加入強度為r 的標準正態(tài)分布隨機擾動。偏離存在性驗證時r 取1‰，偏離普遍性驗證時r 分別取1‰、1%、10%及50%。一般認為，生產(chǎn)環(huán)境中，干擾強度至少在10%以上，1‰和1%的干擾屬于理想情況，生產(chǎn)中不可能做到。

3.2 定位算法的求解

擬牛頓法是高效的代數(shù)類優(yōu)化算法，適應(yīng)性強且收斂速度快，其基本原理可參考數(shù)值分析的書籍，主流的數(shù)值分析工具如Matlab、Scipy 中都有擬牛頓法的實現(xiàn)。實際生產(chǎn)中，典型的經(jīng)典蓋革法，使用擬牛頓法對L2目標函數(shù)進行優(yōu)化。

遍歷搜索法以一定的粒度在搜索區(qū)域內(nèi)遍歷離散網(wǎng)格，對每個網(wǎng)格點的目標函數(shù)進行計算，搜索最優(yōu)值。遍歷搜索法存在著粒度與計算能力的矛盾，粒度過粗無法滿足生產(chǎn)要求，粒度過細則過于消耗算力，無法在合理的時間給出計算結(jié)果。對于長寬深各3 000 m 的空間，如果做到直線10 m 的精度，在實驗室環(huán)境的主流服務(wù)器上運行1 次L2目標函數(shù)的遍歷搜索是百秒量級，生產(chǎn)環(huán)境下遍歷搜索法顯然不能使用。遍歷搜索法可以提供空間的整體信息，生成用于對比和驗證的基線，防止算法收斂于局部最優(yōu)。相比于其它優(yōu)化算法，遍歷搜索法更多用于研究探索，不太適合于生產(chǎn)實踐。

3.3 偏離存在性的驗證

對第1 次標定炮，干擾強度r 取1‰，以不同的隨機種子按式（7）生成3 組到時數(shù)據(jù)，分別取L1、L2和Ls目標函數(shù)，擬牛頓法定位情況見表2，遍歷搜索法定位情況見表3。

表2 擬牛頓法定位情況（1‰正態(tài)擾動）Table 2 Locating results of quasi-Newton method（1‰ normal disturbance）

表2、表3 的數(shù)據(jù)很好地支持了全局最優(yōu)值可能會離實際震源點非常的遠的結(jié)論。對第2～第12次標定炮，進行了同樣的實驗，可以發(fā)現(xiàn)：

1）L1、L2和Ls目標函數(shù)的偽震源點是存在的。

2）遍歷搜索可以發(fā)揮基線的作用。遍歷搜索法不存在陷入局部最優(yōu)的問題，可以保證給出“不會太差”的優(yōu)化結(jié)果，遍歷搜索法在各種情況下都給出了偽震源點。

3）L1目標函數(shù)不適合擬牛頓法，多收斂到局部最優(yōu)。L2和Ls目標函數(shù)有更好的一致性和穩(wěn)定性，L2和Ls目標函數(shù)在兩個優(yōu)化算法中都找到了偽震源點。

4）表2 中的L2目標函數(shù)被擬牛頓法優(yōu)化即是典型蓋革法，Ls對應(yīng)的即是D4DA 算法，在仿真數(shù)據(jù)上D4DA 算法的表現(xiàn)略好。

3.4 偏離普遍性的驗證

使用上節(jié)的3 組數(shù)據(jù)相同的隨機種子，干擾信號強度r 分別取1‰、1%、10%及50%時進行遍歷搜索。不同等級擾動下的偽震源點占比（搜索點數(shù)：30 945 110 個）見表4。

表4 不同等級擾動下的偽震源點占比（搜索點數(shù)：30 945 110 個）Table 4 Proportion of pseudo hypocenters under different disturbances（searched in 30 945 110 points）

可以看出：僅1‰的正態(tài)擾動下，L1、L2和Ls目標函數(shù)的偽震源點數(shù)量即非?？捎^，隨著干擾的增大，偽震源點占比急劇上升，偽震源點是普遍現(xiàn)象。Ls的偽震源占比明顯小于L1和L2，與Ls目標函數(shù)構(gòu)建時對時空區(qū)域信息的使用的預(yù)期相符。另外，由表2、表3 可知：Ls6 次實驗都收斂到偽震源，且定位誤差波動較小，說明Ls目標函數(shù)魯棒性較好。

4 現(xiàn)場應(yīng)用

某礦采用北京安科興業(yè)技術(shù)有限公司的微震監(jiān)控系統(tǒng)，系統(tǒng)正常運行后，選擇在1604 工作面兩巷道的生產(chǎn)幫煤體中打孔，2 次炮孔位置分別選在回風(fēng)巷和運輸巷，炮點及檢波器位置見表5。

表5 某礦1604 工作面炮點及檢波器位置Table 5 Blast positions and detector positions of the 1604 working face in a certain mine

運輸巷標定炮的微震信號被1#～6#和8#檢波器有效拾取，回風(fēng)巷標定炮的微震信號被1#～6#和8#～9#檢波器有效拾取。

分別采用典型蓋革法和D4DA 算法對2 次標定炮進行定位，2 種算法分別使用L2目標函數(shù)和Ls目標函數(shù)，都使用擬牛頓法進行優(yōu)化。定位結(jié)果如下：

1）回風(fēng)巷。①典型蓋革法：定位震源點坐標（38 448 402.45，3 911 206.75，-438.36），發(fā)震時刻580.43 ms，偏離18.13 m；②D4DA 法：定位震源點坐標（38 448 402.38，3 911 207.05，-438.64），發(fā)震時刻580.41 ms，偏離15.61 m。

2）運輸巷。①典型蓋革法：定位震源點坐標（38 448 388.53，3 911 035.21，-436.97），發(fā)震時刻556.62 ms，偏離11.38 m；②D4DA 法：定位震源點坐標（38 448 388.33，3 911 036.00，-437.62），發(fā)震時刻556.76 ms，偏離8.41 m。

可以看出：Ls目標函數(shù)較普通的L2目標函數(shù)有更好的定位結(jié)果。

5 結(jié) 語

1）基于到時的震源定位算法理論成熟、定位精度高，是業(yè)內(nèi)使用的主流算法?；诘綍r的震源定位算法中偽震源現(xiàn)象是存在的，而且是普遍的。提出了偽震源點占比的計算方法，并以此對定位算法中的目標函數(shù)進行定量評估。構(gòu)建了基于分布信息和區(qū)域信息的Ls目標函數(shù)，進一步形成D4DA 算法，并在生產(chǎn)數(shù)據(jù)中對D4DA 算法進行了驗證。

2）分布信息和區(qū)域信息的使用是D4DA 算法的核心，在這2 種信息的使用上進行了初步的嘗試，相信更深入的研究會給震源定位工作帶來新的思路。到時類定位算法研究充分，有許多目標函數(shù)的構(gòu)建思路和優(yōu)化算法，如雙差法、單純形法等，這些方法與分布信息和區(qū)域信息相結(jié)合，有望提升震源定位精度，更好地保障煤礦生產(chǎn)安全。