基于單純形-最短路徑射線追蹤的微震震源混合定位算法

王 輝，梁 苗，朱夢博

(1.鄂爾多斯職業學院，內蒙古 鄂爾多斯 017000；2.鄂爾多斯生態環境職業學院，內蒙古 鄂爾多斯 017000；3.煤炭科學研究總院，北京 100013)

近20年來，隨著我國經濟快速發展，深部礦產資源開采、深埋隧洞開挖、高陡邊坡治理等巖土工程活動進入井噴階段，由此促進了微震監測技術的研究與應用。震源定位是微震監測技術的基礎工作之一，定位的準確性直接影響后續的各種分析工作[1-3]。圍繞震源高精度定位問題，國內外學者提出了眾多定位方法，包括Geiger法、單純形法、網格搜索法、粒子群法和相對位置法等傳統方法。基于上述經典算法，近幾年來國內學者提出了大量的改進方法，包括混合定位方法[4]、基于到時質量控制的定位方法[2,5]、基于觀測系統優化的定位方法[6-8]等。這些改進后的方法較好地解決了傳統定位方法存在的計算效率低、精度不高、容易陷入局部最優解等問題，但是這些方法仍然基于均勻速度假設。當前，復雜巖土工程對震源定位精度提出了更高的要求，同時監測區域越來越大，地質條件越來越復雜，因此基于均勻速度模型的微震震源定位方法應用效果較差，迫切需要開展復雜速度模型下的震源高精度定位方法研究。

針對復雜速度模型下的震源定位問題，目前主要是通過模擬地震波傳播路徑來計算地震波走時場，主要應用于天然地震震源定位和地震勘探領域。近幾年來，隨著地震波射線追蹤相關理論被引入到微震監測領域及計算機性能的提高，部分研究人員開始研究復雜速度模型下的微震震源定位問題。郭超等[9]基于FSM算法及到時差數據庫技術研究了層狀介質下震源高精度定位方法，并將該方法應用于寧武盆地武2-2井水力壓裂工程。郭亮等[10]針對實際工程中帶有空洞和速度分區等復雜速度模型，通過對FMM進行改進，提出了用于計算地震波走時場的MSFM算法，并將之應用與白鶴灘水電站左岸邊坡微震定位。本研究基于地震波傳播原理，采用最短路徑法計算了地震波射線追蹤路徑，并計算到時場；然后基于觀測系統到時場數據庫和單純形(Simplex)定位方法，建立了單純形-最短路徑射線追蹤定位算法。

1 地震波傳播原理概述

1.1 費馬原理

費馬基于光傳播路徑的研究成果總結提出了費馬原理。對于地震波傳播問題，費馬原理的通俗表述為：地震波沿垂直于波前面的路徑(波的射線方向)傳播所需的時間最短。 費馬原理從運動學的角度描述了波的射線傳播特征，由該原理可以推導出波的直線傳播定律、反射定律和折射定律(斯奈爾定律)。

1.2 惠更斯原理

惠更斯原理是指波在介質中傳播到的每一個點都可以視為一個新的波源，即子波源；每個子波源以所在位置處的速度向各個方向傳播微弱的子波?；诨莞乖?，可以根據當前波前計算一定時間間隔后的波前。如圖1所示二維平面，S0為t0時刻的波前，分析S0波前上的4個子波源R1、R2、R3和R4，子波源從t0時刻開始向外傳播子波，Δt時間后，各子波面為一個半徑為V*Δt(V為波速)的圓，4個子波面的包絡面即為t0+Δt時刻的近似波前S1，加密子波源即可得到精度更高的波前S1。

圖1 地震波波前擴展原理Fig.1 Seismic wavefront propagation principle

2 最短路徑射線追蹤方法

最短路徑射線追蹤方法的理論基礎是費馬原理和惠更斯原理，具體利用Dijkstra最短路徑算法進行實現。最短路徑射線追蹤算法流程如圖2所示，各步驟作用及實現方法如下所述。步驟1：根據監測區域大小及設置的網格步長，建立三維網格速度模型；構建兩個節點集合A和Q，A為已經計算過走時的節點集合，Q為作過子波源的節點集合；將震源節點初始化到集合A和Q中，將震源節點走時初始化為0，其余節點的走時初始化為無窮大；步驟2：確定與子波源m相連的所有節點集合V，并將相應節點加入到集合A中；步驟3：確定節點集合S(S∈A，且SQ)中節點n到子波源m之間的走時tmn；步驟4：取{tm，tm+tmn}的最小值更新節點n的走時tn，并記錄使節點n形成最小旅行時的前一個節點號；步驟5：取節點集合S中走時值最小的節點更新子波源m，并將該節點加入到集合Q中；步驟6：判斷集合Q是否包括所有模型節點，若否則執行步驟2～步驟5，若是則完成走時計算；步驟7：根據步驟4中記錄的前節點號，推算從源點到接收點的射線路徑。若只是計算地震波走時場，不計算射線路徑，則無需執行步驟7。

圖2 最短路徑射線追蹤算法流程Fig.2 Shortest path ray tracing flow

需要指出的是，對于三維模型，波源的子波面包括4種情況：①波源在3個邊界面交點上，其子波源包含7個節點(圖3(a))；②波源在兩個邊界面交線上(除端部兩個節點)，其子波源包含11個節點(圖3(b))；③波源在單個邊界面的內部，其子波源包含17個節點(圖3(c))；④波源在模型內部，其子波源包含26個節點(圖3(d))。

圖3 三維模型下波源與子波面空間關系Fig.3 Spatial relationship between wave sourceand wavelet surface

3 單純形-最短路徑射線追蹤定位算法

3.1 定位原理

對于微震震源定位問題，可以將檢波器先視為震源點，模型中其他節點為檢波點，采用最短路徑射線追蹤方法，即可計算出觀測系統中每個檢波器到任意節點的地震波走時及其射線路徑。根據檢源互換原理[11-12]，可以得到監測區域內任意潛在震源節點到各檢波器的走時及其射線路徑，因此震源定位問題就可以轉化為以下數學問題：搜索速度模型中滿足式(1)的節點，該節點即為震源的近似位置。

i=1,2,…,l;j=1,2,…,m;k=1,2,…,n

(1)

圖4 單純形-最短路徑射線追蹤定位算法
Fig.4 A new source location algorithm based on simplex and shortest path ray tracing

3.2 定位數值模擬實驗

為了論述單純形-最短路徑射線追蹤定位算法的具體實施步驟及其優勢，本文建立了礦山常見的層狀速度模型和帶空區速度模型，分別進行了震源定位數值模擬實驗，并將計算結果與常規單純形定位方法的結果進行對比。

3.2.1 層狀速度模型

含煤地層是典型的層狀速度模型，一般煤層是低速層，頂底板巖層是高速層。建立如圖5所示的300 m×300 m×300 m的含煤層狀地質模型，設定頂板、煤層和底板的P波波速分別為3 500 m/s、1 500 m/s和3 500 m/s，設定網格步長為5 m，震源點坐標為(0，0，0)。圖6為震源點到的射線路徑對比圖，設定檢波點(300，300，300)，其中，灰色虛線為最短路徑射線追蹤方法正演結果，黑色實線為直射線。由圖6可知，最短路徑射線追蹤方法計算得到的射線路徑在波速界面有折射現象，滿足斯奈爾定律。圖7(a)和圖7(b)分別為最短路徑射線追蹤方法和均勻速度直射線方法計算得到的走時場，其中均勻速度直射線方法中取均勻速度為2 833 m/s。由圖7可知，在速度分界面處地震波波陣面形狀明顯發生改變。

圖5 頂板-煤層-底板速度模型Fig.5 Velocity model of roof -coal-floor

圖6 兩種方法下震源點到檢波點的射線路徑Fig.6 Ray paths from source point to detection pointunder two methods

圖7 兩種方法走時場計算結果(單位：s)Fig.7 Travel-time fields of two methods (Unit:s)

圖8 含空區速度模型Fig.8 Velocity model with empty zone

3.2.2 帶空區速度模型

微震監測工程中空區速度模型較為常見，包括金屬礦的采空區、煤礦的采空區和地下廠房及隧道的開挖區等。建立如圖8所示的300 m×300 m×300 m的含空區速度模型，空區范圍為{135

圖9 兩種方法下震源點到檢波點的射線路徑Fig.9 Ray paths from source point to detectionpoint under two methods

圖10 兩種方法走時場計算結果(單位：s)Fig.10 Travel-time fields of two methods (Unit：s)

4 工程應用實例分析

本文以神東某礦31101工作面堅硬頂板壓裂為例，該工作面直接頂巖性為粉-細砂巖，老頂為致密堅硬的細粒砂巖，厚度為1.2～18.8 m，大部分頂板屬于中等冒落～難冒落頂板。此次壓裂工程中，采用微震監測技術評價壓裂效果，微震觀測系統為井地聯合方式。監測區下覆地層依次為：黃土層，平均厚度35 m，速度，P波波速1 553 m/s；風化-粗粒砂巖層，平均厚度25 m，P波波速2 931 m/s；細粒砂巖層，平均厚度20 m，P波波速3 374 m/s；泥巖層，平均厚度8 m，P波波速2 693 m/s；粉砂巖層，平均厚度25 m，P波波速3 655 m/s；細粒砂巖層，平均厚度8 m，P波波速3 312 m/s；開采煤層，平均厚度4.8 m，P波波速2 287 m/s。

本次監測過程中，采用放標定炮的方式試驗微震監測系統的精度，標定炮震源空間坐標為(8 542.5，2 137.8，482.2)。根據地層信息，按照地層厚度與速度加權方式計算得到平均波速為2 764 m/s。分別采用本文方法和平均速度模型單純形方法進行標定炮位置定位，與該標定炮位置進行比較，驗證本文定位方法的有效性。其中，最短路徑射線追蹤計算階段中，設定網格間距5 m。微震觀測系統及兩種方法下標定炮定位結果如圖11所示，表1為標定炮定位誤差分析。由圖11和表1可知，兩種方法均較為準確地定位到標定炮位置，但單純形-最短路徑射線追蹤方法精度高很多。 實驗表明，本文提出的單純形-最短路徑射線追蹤方法能實現震源高精度定位。

圖11 微震觀測系統及標定炮定位結果Fig.11 Microseismic observation system andblast location results

表1 標定炮定位誤差分析Table 1 Location error of two methods

5 結 論

1) 本研究基于檢源互換原理，采用三維最短路徑射線追蹤方法計算了地震波走時場及射線路徑，包括常見的層狀速度模型和含空區速度模型等兩種情況。層狀速度模型下，檢源路徑在波速分界面發生折射；空區速度模型下，檢源路徑能繞開空區。兩個數值模擬實驗表明，最短路徑射線追蹤方法計算得到的走時場及射線路徑滿足地震波傳播的基本規律。

2) 采用最短路徑射線追蹤方法、單純形方法建立了針對復雜速度模型下的微震震源定位方法，分兩步驟定位微震震源位置。相比于基于走時場的直接搜索方法，新方法計算效率高；相比于基于均勻速度模型的傳統單純形方法，新方法定位精度高。

3) 基于以神東某礦31101工作面堅硬頂板壓裂案例，以標定炮為參考，常規單純形方法定位誤差為19.63 m，單純形-最短路徑射線追蹤算法定位誤差4.35 m，震源定位精度大幅度提高，滿足現場需求。