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

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

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

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

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

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

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

1.1 經(jīng)典蓋革法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2.1 確定到時偏差的尺度。

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

式中：σi為到時偏差ai的標(biāo)準(zhǔn)差，ms；σs為起震時刻標(biāo)準(zhǔn)差，ms；σgi為走時標(biāo)準(zhǔn)差，ms；σoi為觀測到時標(biāo)準(zhǔn)差，ms。

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

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

式中：σi為到時偏差的標(biāo)準(zhǔn)差，ms；i 為檢波器編號；σgi為走時的標(biāo)準(zhǔn)差，ms；std（·）為取標(biāo)準(zhǔn)差；tgi為走時，ms。

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

實(shí)際計(jì)算σi所使用的形心法如圖2。

圖2 形心法計(jì)算到時偏差的標(biāo)準(zhǔn)差Fig.2 Calculating standard deviation of arrival time deviations by the centroid-method

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

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

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

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

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

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

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

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

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

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

2.3 D4DA 算法

D4DA 算法流程圖如圖3。

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

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

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

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

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

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

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

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

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

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

3.2 定位算法的求解

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

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

3.3 偏離存在性的驗(yàn)證

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

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

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

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

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

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

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

3.4 偏離普遍性的驗(yàn)證

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

表4 不同等級擾動下的偽震源點(diǎn)占比（搜索點(diǎn)數(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目標(biāo)函數(shù)的偽震源點(diǎn)數(shù)量即非常可觀，隨著干擾的增大，偽震源點(diǎn)占比急劇上升，偽震源點(diǎn)是普遍現(xiàn)象。Ls的偽震源占比明顯小于L1和L2，與Ls目標(biāo)函數(shù)構(gòu)建時對時空區(qū)域信息的使用的預(yù)期相符。另外，由表2、表3 可知：Ls6 次實(shí)驗(yàn)都收斂到偽震源，且定位誤差波動較小，說明Ls目標(biāo)函數(shù)魯棒性較好。

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

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

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

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

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

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

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

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

5 結(jié) 語

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

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