基于TFA-DC到時拾取的雙震相微震震源定位方法

賈寶新，李 峰

(遼寧工程技術大學 土木工程學院，遼寧 阜新 123000)

隨著深部礦產資源開采的快速發展，針對深部區域開采環境的安全性研究日益深入。由于蘊含巨大的地應力與構造應力，深部區域的礦產資源在開采的同時會伴隨著一些諸如沖擊地壓、采空區垮落、巖爆、巷道片幫冒頂等環境與安全隱患。

為了監測與預防這些安全隱患，如今礦產資源開采過程中普遍采用微震監測技術對災害源進行監測定位，其中微震源定位方法是微震監測的核心內容，也是提高災害預測精確度的關鍵。提高微震源定位方法的定位精度一直是微震監測技術研究的重點與難點。

微震定位方法主要是利用檢波器收集到的微震信號進行微震波初至到時拾取，并將其與各檢波器空間坐標一同代入優化算法中，從而反演出微震事件發生位置及時刻。因P波波速高于其他震相，所以當P波到達檢波器時，其微震信號波形不易被干擾，進而有助于其初至時刻的精確拾取，故如今國內外學者主要是通過拾取P波初至到時來進行反演定位或者基于P波初至到時對震源定位精度進行相關研究。如王輝等提出的基于單純形-最短路徑射線追蹤的微震震源混合定位算法；李楠等提出的基于L1范數統計的單純形定位方法；王云宏等針對傳統定位方法過度依賴走時精度問題所提出的改進定位目標函數的差分進化微震定位方法；姜天琪等提出的基于網格搜索-牛頓迭代法的微震震源定位算法；賈寶新等提出的小尺度下基于高密度臺陣和粒子群算法的微震源定位方法；LI等通過研究傳感器陣列的空間布局與放置方向提出的可提高定位精度的傳感器網絡優化原則；JENNIFER等通過優化臺站空間布局提出的改進多陣列震源反演定位方法；SOLEDAD等利用方位角優化結合模擬退火與粒子群算法提出的改進全局搜索算法。在此基礎上，也有對雙震相震源反演定位進行研究的，如TIAN等利用P波、S波初至到時時差提出的實時更新速度模型的微震定位事件交叉雙差反演方法；王志剛通過研究不同臺陣分布對應的不同雙曲面場提出的雙波定位原理。

以上研究成果在一定程度上提高了震源定位的精確度并闡明了雙震相定位的作用，但依然存在著一些問題：由于其他震相的干擾導致S波難以識取到時的問題；P波單震相與S波單震相在震源定位方面的區別與聯系沒有闡明的問題；單獨使用S波到時的震源定位效果不明的問題。

針對以上問題，筆者利用基于時頻分析的下山比較法(Time-Frequency Analysis-Downhill Comparison Method，TFA-DC方法)對微震事件信號進行處理，將其所能拾取到的P波精確到時與S波峰值代入整合波速后的雙震相目標函數中，并使用粒子群算法進行目標函數最小值搜索，以此提出了同時使用P波精確到時與S波峰值到時的雙震相微震定位方法——基于TFA-DC到時拾取的雙震相定位方法(TFA-DC-double seismic phases location method，TD-DL方法)。

1 微震定位基本原理

目前，微震源反演定位廣泛使用基于走時殘差理論的震源定位方法。理論計算走時與實際觀測到時的差即走時殘差，利用計算走時殘差的不同方法，可構造出不同的微震定位目標函數。對應震相應選取相應的目標函數，筆者選取分別應用于單震相定位與雙震相定位的2種目標函數。

1.1 單震相目標函數

現有微震定位方法通常采用P波單震相進行到時拾取及震源反演定位。與其他震相相比，P波震相初至響應特征明顯，更利于精確拾取到時，所以其定位精度也較高。筆者采用理論計算走時與實際觀測到時差值的數學期望作為發震時刻的估計值，從而得到P波單震相的走時殘差計算公式：

(1)

(2)

(3)

同理可得S波單震相的走時殘差計算公式：

(4)

(5)

(6)

1.2 雙震相目標函數

單一微震事件通常同時包含P波、S波、面波、反射波等多種震相信息，而除P波外的其他震相因重疊、干擾等問題，通常難以用于定位計算。相關研究表明，對于提高震源定位精度而言，同時使用2種震相的目標函數定位精度要高于只使用單一震相的目標函數。

故此筆者引入同時使用P波單震相與S波單震相的雙震相目標函數，并對其中雙震相的波速進行整合，其目標函數如下：

(7)

式中，(,,,)為雙震相的走時殘差；為震相傳播速度，m/s；為均衡系數，取值為0～1。

因為P波、S波波速未知且時刻發生變化，該變化在傳播過程中的變化幅度也不相同，所以獨立使用2種波速代入計算后很可能導致算法無法收斂。對于定位精度而言，因波速結構的不確定性，同時代入2種單震相的平均波速進行定位反演計算會使定位精度降低。鑒于以上問題，筆者將P波與S波的平均波速設為相同值——震相傳播速度，并與所求震源空間坐標，，一同代入粒子群算法中計算。一則可以避免預先測速，達到實時監測的效果；二則可以降低波速結構對定位精度的影響，使得目標函數盡可能趨于最小值；三則可以提高目標函數的收斂性，因為不同發震條件下，巖土體破裂程度不同，從而使擠壓產生的P波與剪切產生的S波的波速均非恒定值，而在小尺度微震監測下，P波與S波的傳播距離短，其速度變化量較小，故單次微震事件中P波與S波的波速比值可視為常數，即震相傳播速度可表征為二者的統一傳播速度。

改進雙震相目標函數的作用原理是根據已有檢波器坐標及雙震相到時形成含有4個自變量的方程。在定義域內，利用全局搜索算法對目標函數最小值進行搜索，最小值對應的4個自變量數值即為所求震源反演坐標及其對應的震相傳播速度。

1.3 P波和S波到時聯合拾取方法

(1)TFA-DC方法簡介。TFA-DC方法是一種可同時獲取微震信號P波精確到時與S波峰值到時的到時拾取方法。該方法首先通過分析原始信號的語譜圖與功率密度譜圖獲取背景噪聲的位置和規律以及P波、S波初至前后微震信號的頻率、振幅、能量變化。然后根據相關算法進行第1次帶通濾波，過濾掉規律且功率大于P波、S波信號的高、低頻背景噪聲，使功率密度譜圖中S波的主頻清晰顯現。

隨后以S波主頻為圓點，規定濾波范圍為半徑進行第2次帶通濾波，使信號圖像更加平滑且適宜迭代比較。由于此半徑的選擇會涵蓋P波以及S波的主頻范圍，并且P波波速要大于S波，所以濾波后波形中P波精確到時所對應的P波初至波形以及S波峰值到時所對應的全信號振幅波形均可清晰顯現，且在大多數情況下不會產生干擾或遮蓋。

最后通過將全子波振幅的數學期望設置為閾值，并遵循P波和S波功率大小、到達先后、波形重疊三大關系對濾波后信號進行下山比較以此獲得P波精確到時，以及P波震相右側第1個波峰所代表的S波峰值到時。

(2)TFA-DC方法的有效性證明。現利用連續小波變換理論(Continuous Wavelet Transform，CWT)，并以圖5所示模型試驗信號為例，對TFA-DC方法所拾取的P波精確到時與S波峰值到時加以說明，原始微震信號及尺度如圖1所示。

圖1 原始微震信號及尺度

排除500 Hz以下的規律且全程持續的背景噪聲，通過圖1中信號圖與尺度圖對應分析可知，該微震信號能量變化主要可分為4個過程：31.983～31.992 s為第1次能量突變，因其最早到達，故可視為由P波到達所引起的振動響應；31.992～31.999 s為第2次能量突變，因其能量大于P波并且僅次于P波到達，故可視為S波到達所引起的振動響應；31.999～32.005，32.005～32.019 s分別為第3，4次能量突變，因其到達晚于P波與S波，故可判斷為反射波、面波等其余震相。

綜上可證明TFA-DC方法所拾取2種到時分別為的微震信號的P波精確到時與S波峰值到時。

(3)到時拾取特點。因帶通濾波對突變信號有一定的平滑作用，信號突變處會產生子波振幅遞增或遞減的平滑過渡現象，以此導致TFA-DC方法中P波精確到時拾取實則早于其真實到達時刻；因檢波器所采集的S波信號自到達至達到峰值需要一定的時間過程，并且有時會疊加其他震相于其上，故S波峰值到時實則晚于其真實到達時刻。所以單獨使用上述P波或S波單震相到時數據進行震源反演定位，可能會造成定位結果相較真實震源有一定的偏移。

1.4 粒子群算法

粒子群算法(Particle Swarm Optimization,PSO)是一種常用的廣域隨機搜索算法，其核心思想是共享各個子個體所歷經的目標函數值信息，使得全部子個體利用所有已知的目標函數值信息進行最優化的定向移動，直到某個子個體找到目標函數的全局最優解。其速度和位置進化方程如下：

(8)

(9)

選定目標函數并導入相關到時及檢波器空間坐標信息后，使用粒子群算法進行目標函數最小值搜索。通過迭代搜索所得的目標函數最小值所對應的自變量數值即所求震源的空間坐標。

1.5 TD-DL方法

筆者提出的TD-DL方法是包括一整套微震定位解決方案的新型定位方法，其中包含可同時拾取P波精確到時與S波峰值到時的TFA-DC雙震相到時聯合拾取方法、整合波速后的雙震相目標函數、粒子群全局搜索算法，以上方法相結合便構成了可自動拾取到時并定位計算的TD-DL方法。

2 模型試驗

利用超高頻多通道構造活動監測儀對擺錘撞擊產生的微震信號進行處理分析，以此探究單震相定位方法與TD-DL方法在震源反演定位方面的區別與聯系，并對TD-DL方法的原理及實際震源定位效果進行驗證。

2.1 模型建立

本試驗模型長為2.3 m、寬為1 m、高為1.3 m，采用石英砂、石膏、石灰、水混合后分層堆砌而成，其質量配比約為6∶3∶1∶2。22個檢波器在模型脫模后鉆孔嵌入墻體，檢波器三維空間坐標見表1。

表1 各檢波器三維空間坐標

圖2中黑色球體為單分向微震檢波器，其周圍數字為檢波器序號；黑色立方體為撞擊位置，箭頭指向黑色立方體的字母為震源序號。10～18號檢波器按規定位置設置在模型的=0面，此面記作面；19～22號檢波器按規定位置設置在模型的=230 cm面，此面記作面；1～9號檢波器按規定位置設置在模型的=100 cm面，此面記作面。

圖2 試驗布局

2.2 監測設備

本試驗采用采樣頻率為100 kHz、可監測頻帶寬度為0～50 kHz的超高頻多通道構造活動監測儀(Antenna-Ⅲ)記錄單擺撞擊模型表面產生的微震信號，設備參數見表2，超高頻多通道構造活動監測儀各組成部分如圖3所示。

表2 超高頻多通道構造活動監測儀設備參數

圖3 超高頻多通道構造活動監測儀

2.3 試驗過程

本試驗采用球體定長單擺撞擊模型指定位置，連接鋼球的連接線采用長為30 cm的棉質線，撞擊球體采用半徑為1 cm的鋼球。模型制作好后帶模具養護7 d，以防開裂。卸下模具后，保持開窗通風且每天至少5 h日照。30 d后檢查時，已可確保模型完全干透。因制作材料本身的脆性，鋼球撞擊模型表面時可使表面下陷0.3 cm左右，所以此過程可以保證模型受擠壓破壞產生P波，也可保證模型擠壓破壞區同時也受剪切破壞，從而產生S波。

(1)撞擊位置。試驗設計撞擊位置共3個，分別命名為“震源”、“震源”、“震源”，其三維空間坐標見表3。

表3 各震源三維空間坐標

(2)撞擊方式。在試驗準備階段將3個撞擊位置用記號筆做好標記，試驗時先將鋼球自然懸掛并剛好接觸模型表面的指定撞擊位置。隨后將鋼球拉至垂直于模型表面5 cm后放手，使之自由下落并撞擊指定位置，表3中的3個位置各撞擊10次，試驗現場照片如圖4所示。

圖4 試驗現場

(3)數據收集。準備好鋼球和連接線并將其安排在指定位置后，開啟超高頻多通道構造活動監測儀并對其進行測試，完畢后進行撞擊試驗，完成后將所得試驗數據保存至硬盤指定位置。

2.4 預先波速分析

通過真實震源位置與P波精確到時、S波峰值到時正演P波、S波的平均波速，結果顯示：P波波速為5 600～7 400 m/s，S波波速為5 100～6 800 m/s，震相傳播速度為5 100～7 400 m/s。震相傳播速度的變化范圍涵蓋了P波、S波2種單震相波速的變化范圍，說明震相傳播速度可以很好的表征2種震相的傳播。P波平均波速大于S波平均波速，并且P波波速的變化范圍略大于S波波速，說明本試驗模型的密實度較高，波速變化規律也與巖石中傳播相近。

3 單雙震相定位分析

通過對模型試驗所得微震信號進行處理并獲取到時數據后，利用單、雙震相目標函數以及粒子群算法獲取3種到時組合的反演定位結果。最后通過3種到時組合所得震源反演坐標與定位誤差的對比分析，對TD-DL方法的有效性進行驗證。

3.1 到時拾取

利用TFA-DC方法處理共得到震源、震源、震源的各10組到時數據，每組到時數據包含P波精確到時22個、S波峰值到時22個，最后一共得到3處震源的1 320個到時數據。TFA-DC方法拾取雙震相到時如圖5所示(對應原始信號如圖1所示)。圖5中S波信號振幅約為P波的6倍，而圖1中該比值為3倍，這是由于TFA-DC方法對原始信號進行了連續2次帶通濾波，故此信號中P波與S波的信號振幅與特征會較原始信號中的略有不同，但此區別并不會影響到時的拾取，該過程產生的光滑信號反而更有利于TFA-DC方法中下山比較法的執行，使其中迭代比較子波振幅的過程容錯率更高。根據1.3.3節所述并對比圖1與圖5可知，TFA-DC方法所得雙震相到時中P波精確到時略早于真實到時，而S波峰值到時略晚于真實到時。故以下著重分析在利用上述到時的前提下，單、雙震相定位結果之間的聯系與區別。

圖5 雙震相到時拾取示意

3.2 單、雙震相定位結果分析

模型試驗定位結果分析的目的是找到P波精確到時與S波峰值到時單獨或聯合使用時，所得定位結果之間的區別與聯系。為方便分析敘述，下述內容只對定位結果的軸與軸坐標進行分析，軸坐標同理，不再贅述。

..震源定位結果

使用P波、S波單震相反演定位計算時，將到時數據與檢波器空間坐標數據分別應用于式(3)與式(6)，可得二者的定位結果。TD-DL方法需要將雙震相到時與檢波器空間坐標數據同時應用于式(7)，可得其定位結果。上述反演定位計算結果如圖6所示。

圖6 震源I，J，K反演定位結果

..定位結果圖像分析

據圖6可知，3種到時組合的定位結果大體可描述為：以TD-DL方法定位結果為中心，P波、S波單震相定位結果以相近的距離分布于相反方向的兩側。

相較真實震源而言，TD-DL方法定位結果表現為圍繞或鄰近于真實震源的四周，而P波、S波單震相的定位結果與真實震源之間呈方向相反但距離同步變化的趨勢。

..TDDL方法原理分析

據1.3.3節可知，TFA-DC方法所拾取的P波精確到時較真實到時略早，而所拾取的S波峰值到時較真實到時略晚。這就導致通過二者所得的定位結果會以相反的方向偏移真實震源。通過模型試驗證明，二者相對于真實震源所偏移的距離大致相同，可見其到時誤差也大致相同。若到時誤差完全相同，因目標函數中2種震相所占權重相同，故TD-DL方法定位結果應與真實震源一起分布于2種單震相定位結果的中間位置。但由于二者拾取過程不同，加之現場其余震相的干擾，導致其到時誤差不完全相同。在此情況下，得益于增加了S波到時的約束信息，TD-DL方法提高了對檢波器坐標、雙震相到時等已知信息的統合能力，減小了到時誤差不同的影響。通常對于不同的傳播路徑，每個傳感器接收的S波信號自到達至峰值的上升時間是不一樣的。各檢波器中S波到時誤差的不同會因誤差疊加效應而增加其定位誤差，從而使得單獨使用S波峰值到時的定位精度要低于單獨使用P波精確到時的定位精度，但目標函數求最小值的過程會使TD-DL方法定位結果趨于2種單震相到時誤差更小的一方，使其定位結果可以在一定程度上偏向真實震源，從而提高了定位精度。

綜上所述，在使用TFA-DC方法拾取P波、S波單震相到時的基礎上，TD-DL方法可在一定程度上抵消P波、S波單震相計算所得的定位誤差，并使定位結果趨于2種單震相定位精度高的一方，從而達到提高震源反演定位精度的效果。

3.3 單、雙震相定位誤差分析

對震源、震源、震源的各10組到時數據進行反演定位計算，共處理得到P波、S波單震相以及TD-DL方法定位結果各30個。現將以上各定位結果距真實震源的距離作為定位誤差，并將定位誤差按照震源、震源、震源分組作圖，所得結果如圖7所示。

自圖7可知，TD-DL方法定位精度與穩定性要高于P波、S波單震相定位。定位精度可用定位誤差平均值表示，定位穩定性可用定位誤差標準差表示，各震相在不同震源下的定位精度與穩定性見表4。

圖7 震源I，J，K反演定位誤差

由表4可知，震源下TD-DL方法定位誤差平均值分別為P波、S波的23.0%和18.0%；震源下TD-DL方法定位誤差平均值分別為P波、S波的24.6% 和18.8%；震源下TD-DL方法定位誤差平均值分別為P波、S波的24.5%和20.1%。震源下TD-DL方法定位誤差標準值分別為P波、S波的50.4% 和32.8%；震源下TD-DL方法定位誤差標準值分別為P波、S波的42.6%和44.5%；震源下TD-DL方法定位誤差標準值分別為P波、S波的67.7% 和35.4%。

表4 各震相在不同震源下的定位效果分析

由此可見，同時使用P波精確到時與S波峰值到時的TD-DL方法的定位精確度與穩定性要優于P波、S波單震相定位方法。

3.4 檢波器布置對定位精度的影響

為了探究檢波器布置對定位精度的影響，首先將定位誤差分別在各坐標軸上做投影，獲得各軸的軸向定位誤差，隨后通過作圖分析二者之間的對應關系。

..軸向定位誤差

軸向定位誤差見表5，其中軸向定位誤差為此震源下對應震相的定位結果與真實震源的之間的距離在對應坐標軸上的投影的平均值。

表5 各震相在不同震源下的軸向定位誤差

..檢波器布置與定位精度的關系

以定位精度最高的TD-DL方法所得的定位結果為例說明檢波器布置與定位精度的關系，其他方法所得結論相同，以下不再贅述。

由圖2可知，10～18號檢波器位于模型面，19～22號檢波器位于模型面，1～9號檢波器位于模型面。已知在模型面上軸坐標固定等于0，所以此面的檢波器對于軸坐標沒有變化范圍，即此面對于軸與軸是自由面。同理可知模型面為軸與軸的自由面，模型面為軸與軸的自由面。

故軸上坐標可變的檢波器為1～18號，共18個；軸上坐標可變的檢波器為19～22號，共4個；軸上坐標可變的檢波器為1～22號，共22個。各坐標軸上的平均定位誤差與對應坐標軸上坐標可變檢波器數量之間的關系如圖8所示。

圖8 軸向定位誤差與對應軸坐標可變檢波器數量的關系

根據相關理論，檢波器空間密度越高，則震源反演定位誤差越小，而對于垂直于各坐標軸的面而言，通常將檢波器布置于此類面上相較于布置于監測范圍內部來的容易，故此若找到此類面上檢波器布置個數與定位誤差之間的關系，便可優化檢波器空間布置并減小人力物力浪費。

由圖8可知，各坐標軸上坐標可變檢波器的數量與定位結果在相應坐標軸上的誤差呈負相關關系。換言之，高精度定位對檢波器空間高密度布置的要求可簡化為檢波器在各軸對應自由面上的高密度布置，或者說在一定監測范圍內，某一坐標軸上坐標可變的檢波器越多，定位結果在此軸上的表現便越為精準。

4 工程驗證

為了進一步證實TD-DL方法在震源反演定位方面的優越性和可靠性，選取某礦微震監測工程中的礦震信號進行到時拾取及反演定位計算，并通過與P波、S波單震相的定位計算結果進行比較，以此對TD-DL方法進行驗證。

4.1 到時拾取

選取不同時期發生的3次礦震所產生的微震信號，并將其分別命名為KZ-1，KZ-2，KZ-3。使用TFA-DC方法對以上3組微震信號進行到時拾取，濾波處理后所得到時拾取如圖9所示，所得P波精確到時與S波峰值到時數據見表6。

圖9 礦震信號到時拾取示意

表6 礦震信號P波與S波到時數據

4.2 檢波器布置

礦山微震監測工程檢波器布置在地表以下2 m深處，僅用于二維監測，即監測震源位置的軸與軸坐標，故檢波器的坐標變化范圍主要分布在軸與軸，其坐標見表7。

表7 礦山各檢波器三維空間坐標

4.3 定位結果分析

為增加雙震相到時的約束信息，提高定位算法的收斂性與定位精度，在計算時除檢波器軸、軸坐標外同時導入軸坐標協同定位，定位結果只采用軸、軸坐標。雖檢波器軸坐標變化范圍遠小于軸與軸，但由于軸坐標的加入，算法可以增加收斂容錯率并減小收斂后的目標函數值。由此說明加入了軸坐標后，可達到提高定位精度的效果。導入檢波器坐標及到時數據后，利用粒子群算法進行反演定位計算，其結果見表8，各震源真實坐標見表9。

表8 礦震信號各震相定位結果

表9 礦震各震源真實坐標

將上述反演定位結果距離真實震源的距離定義為定位誤差，各震相的定位誤差見表10。由表10可知，在震源KZ-1下，TD-DL方法定位誤差分別為P波、S波單震相的14.8%和17.1%；在震源KZ-2下，TD-DL方法定位誤差分別為P波、S波單震相的7.3% 和8.9%；在震源KZ-3下，TD-DL方法定位誤差分別為P波、S波單震相定位誤差的5.2%和15.9%。

表10 礦震信號各震相定位結果分析

由此可見，TD-DL方法對比P波、S波單震相定位方法，在工程現場數據處理方面更為精準且穩定，說明該方法是一種有效的震源反演定位方法。

5 結 論

(1)TD-DL方法由于S波到時約束信息的增加以及目標函數的最小值趨向作用，提高了對檢波器坐標、雙震相到時等已知信息的統合能力，減小了2種單震相到時誤差不同的影響，并使其定位結果趨于2種單震相定位結果較真實震源誤差更小的一方。故在TFA-DC方法拾取P波、S波單震相到時的基礎上，TD-DL方法可在一定程度上抵消P波、S波單震相計算所得的定位誤差，從而達到提高震源反演定位精度的效果。

(2)TD-DL方法的定位精確度與穩定性要優于P波、S波單震相定位方法。模型試驗下，TD-DL方法定位誤差平均值分別為P波、S波單震相定位的23.9%和18.9%；TD-DL方法定位誤差標準差分別為P波、S波的50.9%和36.9%。工程數據下，TD-DL方法定位誤差平均值分別為P波、S波單震相定位的8.1%和12.6%。

(3)各坐標軸上坐標可變檢波器的數量與定位結果在相應坐標軸上的誤差呈負相關關系。高精度定位對檢波器空間高密度布置的要求可簡化為檢波器在各軸對應自由面上高密度的布置，即在一定監測范圍內，某一坐標軸上坐標可變的檢波器越多，定位結果在此軸上的表現便越為精準。