王 輝，梁 苗，朱夢博

(1.鄂爾多斯職業(yè)學(xué)院，內(nèi)蒙古 鄂爾多斯 017000；2.鄂爾多斯生態(tài)環(huán)境職業(yè)學(xué)院，內(nèi)蒙古 鄂爾多斯 017000；3.煤炭科學(xué)研究總院，北京 100013)

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

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

1 地震波傳播原理概述

1.1 費(fèi)馬原理

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

1.2 惠更斯原理

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

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

2 最短路徑射線追蹤方法

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

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

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

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

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

3.1 定位原理

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

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 定位數(shù)值模擬實驗

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

3.2.1 層狀速度模型

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

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

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

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

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

3.2.2 帶空區(qū)速度模型

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

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

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

4 工程應(yīng)用實例分析

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

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

圖11 微震觀測系統(tǒng)及標(biāo)定炮定位結(jié)果Fig.11 Microseismic observation system andblast location results

表1 標(biāo)定炮定位誤差分析Table 1 Location error of two methods

5 結(jié) 論

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

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

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