蔣若辰 ，戴峰 ,*，劉燚 ，李昂

a State Key Laboratory of Hydraulics and Mountain River Engineering, College of Water Resource and Hydropower, Sichuan University, Chengdu 610065, China

b School of Highway, Chang’an University, Xi’an 710064, China

1. 引言

地下空間開挖活動(dòng)無可避免地引起工程圍巖的應(yīng)力重分布，進(jìn)而誘發(fā)裂紋的產(chǎn)生與巖石斷裂能的釋放。作為一種有效的實(shí)時(shí)監(jiān)測方法，微震（microseismic, MS）監(jiān)測技術(shù)通過埋設(shè)于工程巖體內(nèi)部的微震傳感器，探測巖石破裂釋放的應(yīng)力波，并基于不同的地球物理反演方法，獲取巖石破裂特征。這一監(jiān)測技術(shù)，可對(duì)即將發(fā)生的地質(zhì)災(zāi)害，諸如巖爆災(zāi)害，提供提前預(yù)警，并采用特定的實(shí)時(shí)支護(hù)措施來保證工程建設(shè)的安全[1-3]。當(dāng)下，微震監(jiān)測這項(xiàng)三維實(shí)時(shí)監(jiān)測技術(shù)，已經(jīng)被廣泛地用于多個(gè)巖石工程領(lǐng)域，如礦業(yè)工程[4-6]、深埋隧洞[1,3]、巖質(zhì)邊坡[7-9]，以及地下廠房[10-12]。

微震震源定位是這項(xiàng)監(jiān)測技術(shù)的重要部分，是其在工程項(xiàng)目中實(shí)現(xiàn)應(yīng)用的基礎(chǔ)[13,14]。對(duì)巖石工程中微震活動(dòng)的解譯，很大程度上依賴于微震震源定位的準(zhǔn)確性。精確且快速的定位方法能夠指引工程圍巖裂紋網(wǎng)絡(luò)分布圈定工作的進(jìn)行。在獲取定位結(jié)果的基礎(chǔ)上，進(jìn)一步的破裂源反演工作才能夠進(jìn)行，用于揭示圍巖變形或者破壞模式，為巖石失穩(wěn)提供提前預(yù)警，減少人員傷亡與設(shè)備損失。因此，微震破裂源定位是微震監(jiān)測技術(shù)領(lǐng)域一個(gè)受到極大關(guān)注的課題。

在微震破裂源定位中，巖體波速模型的確定非常重要。在監(jiān)測范圍較小且?guī)r性均一的情況下，均一波速模型能夠保證微震定位工作快速穩(wěn)定地實(shí)現(xiàn)，因而在巖石工程中得到了廣泛引用。最早的基于數(shù)學(xué)計(jì)算的微震震源定位方法是所謂的Giger法，即將定位問題轉(zhuǎn)化為求解基于微震信號(hào)到時(shí)的線性方程組[15]。在此之后，隨著計(jì)算方法與技術(shù)的發(fā)展，基于這一思路的新方法得以提出與發(fā)展[16-18]。正如文獻(xiàn)[19]所總結(jié)，這些經(jīng)典的定位方法，諸如相對(duì)定位法和聯(lián)合震源定位法，進(jìn)一步提升了微震震源定位的準(zhǔn)確性。然而，這些方法在運(yùn)用到較為復(fù)雜的地質(zhì)環(huán)境中，會(huì)造成相對(duì)較大的誤差。

近年來，一些新方法被提出，用以滿足受復(fù)雜地質(zhì)環(huán)境影響的巖石工程項(xiàng)目中對(duì)微震震源定位的需求。Dong等[20,21]提出了多種應(yīng)用于地下礦井的微震震源定位方法，諸如多步定位法與三維綜合解析法。之后，Dong等[22]與Hu和Dong [23]基于A*搜索算法，提出了無需測速的微震/聲發(fā)射定位方法，用于不規(guī)則復(fù)雜地質(zhì)結(jié)構(gòu)。Feng等[24]開發(fā)了分層波速模型，并系統(tǒng)測試了兩種不同的方法對(duì)于微震震源定位的可行性，結(jié)果展示了所提出的方法能夠減小平均定位誤差[25]。Castellanos與Van der Baan [26]基于礦井中相似波形，提出了一種互相關(guān)的定位方法。Gong等[27]提出了一種應(yīng)用于煤礦中的各向異性波速模型，與均一波速模型相比，該模型能夠提供更準(zhǔn)確的微震震源定位結(jié)果。此外，一些學(xué)者也將微震震源定位問題，轉(zhuǎn)化為高維優(yōu)化問題，并且采用等效波速模型來實(shí)現(xiàn)微震定位。這類方法會(huì)將額外的參數(shù)引入優(yōu)化搜索過程，并且，啟發(fā)式算法，諸如遺傳算法[28]和重力搜索法[29]，也會(huì)被引入以獲取微震破裂源與等效波速模型[30]。然而，在大跨度地下廠房建設(shè)中實(shí)現(xiàn)微震震源定位，仍然面臨諸多挑戰(zhàn)。首先，現(xiàn)有的方法無法實(shí)現(xiàn)在任意復(fù)雜的巖體波速模型中的定位。其次，地下洞室的開挖，特別是大型地下洞室的開挖，會(huì)在巖體中產(chǎn)生形狀不規(guī)則的大空洞區(qū)域。簡單的等效巖體波速模型無法解決這些問題，因此會(huì)干擾最終的定位結(jié)果。

本研究開發(fā)了一種新的微震定位方法，以實(shí)現(xiàn)在含空洞復(fù)雜巖體波速模型中的微震事件定位。建立了基于網(wǎng)格的巖體波速模型，其中不同的網(wǎng)格節(jié)點(diǎn)被賦予不同的縱波速度以反映復(fù)雜的巖體縱波波速分布。之后，二階差分形式的快速行進(jìn)法被引入，用于計(jì)算從破裂源產(chǎn)生的波到達(dá)各個(gè)節(jié)點(diǎn)的理論到時(shí)?；诂F(xiàn)場采集的微震破裂信號(hào)，通過搜索理論到達(dá)時(shí)間與實(shí)際到達(dá)時(shí)間之間殘差最小的最優(yōu)網(wǎng)格節(jié)點(diǎn)，實(shí)現(xiàn)微震破裂源定位。本文所用的方法將在第2節(jié)中介紹，并且其性能表現(xiàn)將在第3節(jié)的一系列數(shù)值仿真實(shí)驗(yàn)中進(jìn)行分析。猴子巖水電站地下洞室開挖過程中所采集的、具有代表性的爆破與微震事件將由本方法來實(shí)現(xiàn)定位，進(jìn)一步展示本方法潛在的應(yīng)用價(jià)值。本文第5節(jié)為主要的結(jié)論。

2. 方法

這一節(jié)將展示利用二階差分形式的快速行進(jìn)法計(jì)算理論到時(shí)的過程。一個(gè)目標(biāo)函數(shù)將被用于計(jì)算理論到時(shí)與實(shí)際到時(shí)之間的殘差值，進(jìn)而實(shí)現(xiàn)微震震源定位。此外，線性插值法與龍格-庫塔法將被用于獲取從破裂源到微震傳感器之間的傳播路徑。

2.1. 二階差分形式的快速行進(jìn)法

對(duì)于地下洞室的微震監(jiān)測，P波到時(shí)通常被用于定位微震事件。在本文中，二階差分快速行進(jìn)法作為一種在離散區(qū)域內(nèi)，采用有限差分方式，以網(wǎng)格為基礎(chǔ)的實(shí)用計(jì)算方法，被用于計(jì)算P波理論初至到時(shí)[31,32]。這一方法能夠有效地避免很多其他射線追蹤方法（如打靶法[33]與彎曲法[34]）所固有的問題。在各向同性的介質(zhì)中，P波的彈性波方程可表示為：

式中， 為標(biāo)量勢函數(shù)；v為P波波速；t為時(shí)間。式（1）的通解可被表示為：

將式[3]與式[4]導(dǎo)入式[1]中，可得式[5]：

式[5]的左側(cè)包含了虛數(shù)項(xiàng)與實(shí)數(shù)項(xiàng)，但式[5]的右側(cè)僅包含實(shí)數(shù)項(xiàng)。故可獲得式[6]與式[7]：

式[6]為傳輸方程，用于計(jì)算A。式[7]在“高頻近似”假設(shè)的基礎(chǔ)上，即時(shí)，可被簡化為程函方程：

式中，T為P波初至到時(shí)；v為微震/地震波的P波波速，它們都是與位置(x,y,z)有關(guān)的函數(shù)。

目前，對(duì)于程函方程，雖然有一些數(shù)值解的方法，但是無法獲得解析解。快速行進(jìn)法將式（8）的微分方程形式轉(zhuǎn)化為差分方程的形式進(jìn)行求解，來獲取P波初至到時(shí)T。原始的計(jì)算區(qū)域?qū)⒈晦D(zhuǎn)化為離散的網(wǎng)格與節(jié)點(diǎn)?；谙嚓P(guān)的地質(zhì)調(diào)查數(shù)據(jù)，巖體與空洞區(qū)域相應(yīng)的網(wǎng)格節(jié)點(diǎn)被賦予相應(yīng)的P波波速值。Sethian與Popovici[35,36]，以及Chopp [37]進(jìn)一步引入熵迎風(fēng)格式（entropy satisfying upwind scheme）來處理梯度不連續(xù)問題，這使得快速行進(jìn)法能夠模擬P波傳播。程函方程的差分形式可被表示為：

同理，相同的差分方式也可用于沿y軸與z軸的差分計(jì)算[圖1（a）]。對(duì)二階快速行進(jìn)法而言，不同階數(shù)的差分算子的使用依賴于二階差分算子是否可用。如果網(wǎng)格節(jié)點(diǎn)的T無法獲得，諸如靠近破裂源點(diǎn)與計(jì)算邊界，那么二階差分的計(jì)算形式將被轉(zhuǎn)化為一階差分的形式。

式（9）描述了在特定網(wǎng)格節(jié)點(diǎn)，采用差分的方式計(jì)算初至到時(shí)。這一計(jì)算過程在不同的待計(jì)算節(jié)點(diǎn)處的連續(xù)實(shí)施，需要正確地確定待計(jì)算節(jié)點(diǎn)的計(jì)算順序，以保證與P波傳播的方向一致。在快速行進(jìn)法中，P波的傳播可以使用窄帶技術(shù)來計(jì)算從迎風(fēng)區(qū)到下風(fēng)區(qū)的初至到時(shí)，其概念如圖1（b）所示。所有的網(wǎng)格節(jié)點(diǎn)被標(biāo)記為激活點(diǎn)、近點(diǎn)與遠(yuǎn)點(diǎn)三類點(diǎn)中的某一類。激活點(diǎn)位于窄帶的迎風(fēng)區(qū)，其初至到時(shí)計(jì)算已經(jīng)完成。近點(diǎn)位于窄帶中，其初至到時(shí)采用式（9）獲得了一個(gè)嘗試值。遠(yuǎn)點(diǎn)位于窄帶的下風(fēng)區(qū)，并未進(jìn)行初至到時(shí)計(jì)算。窄帶的擴(kuò)展過程首先搜索近點(diǎn)中初至到時(shí)嘗試值最小的節(jié)點(diǎn)，將其標(biāo)記為激活點(diǎn)，并將其周圍的所有鄰點(diǎn)中的遠(yuǎn)點(diǎn)標(biāo)記為近點(diǎn)。之后，采用式（9）計(jì)算該新激活點(diǎn)周圍的近點(diǎn)的初至到時(shí)值。窄帶的形狀可被視為P波傳播的初至波陣面，并且上述計(jì)算過程不斷重復(fù)，直到所有網(wǎng)格節(jié)點(diǎn)被標(biāo)記為激活點(diǎn)，完成計(jì)算。

圖1. 二階差分快速行進(jìn)法實(shí)現(xiàn)過程。（a）二階差分形式示意圖；（b）快速行進(jìn)法采用窄帶技術(shù)實(shí)現(xiàn)整體擴(kuò)展。

2.2. 射線路徑追蹤

快速行進(jìn)法有效地將P波傳播簡化為基于網(wǎng)格的射線追蹤[32]。在獲取計(jì)算區(qū)域內(nèi)的所有網(wǎng)格節(jié)點(diǎn)的初至到時(shí)之后，從微震破裂源點(diǎn)到計(jì)算區(qū)域內(nèi)的其他任意節(jié)點(diǎn)的傳播路徑（也就是射線路徑），可以基于初至到時(shí)梯度來獲得。如果將射線傳播路徑視作由多個(gè)直線段組成，并令rn(xn,yn,zn)為第n步計(jì)算之后的射線路徑位置，則在射線路徑上，下一步計(jì)算所達(dá)到的位置如下：

圖2. 基于線性插值的方式計(jì)算到時(shí)梯度。

式中，h為網(wǎng)格間距。同理，B、E、F點(diǎn)的到時(shí)梯度可被表示為：

O點(diǎn)與D點(diǎn)的到時(shí)梯度可以用相同的方式表示為：

因此， 可通過如下計(jì)算獲得：

為了獲得更準(zhǔn)確的rn+1，龍格-庫塔法被用于計(jì)算式（12），其被轉(zhuǎn)化為如下形式：

2.3. 定位方法

微震震源定位通過搜索理論到時(shí)與實(shí)際到時(shí)之間的最小殘差來實(shí)現(xiàn)。假設(shè)誘發(fā)的微震事件發(fā)生在時(shí)刻t0，經(jīng)過Δti的時(shí)長，在時(shí)刻ti達(dá)到第i個(gè)微震傳感器，則對(duì)于第i個(gè)傳感器的理論到時(shí)與實(shí)際到時(shí)之間的殘差ξi滿足如下關(guān)系：

式中，ti可通過提取在巖石工程中采集到的微震信號(hào)的實(shí)際到時(shí)而獲得；Δti可用二階快速行進(jìn)法計(jì)算獲得；而僅使用一個(gè)微震傳感器，無法獲得t0。為了消除t0對(duì)定位結(jié)果的影響，多個(gè)微震傳感器被采用，并且不同傳感器之間的殘差可通過下式獲得：

式中，ξi與ξj分別表示第i個(gè)與第j個(gè)微震傳感器所對(duì)應(yīng)的殘差。雖然在理想條件下，ξi,j等于0，但在實(shí)際中無法實(shí)現(xiàn)（因受多因素的限制，如P波到時(shí)提取，以及Δti和Δtj的計(jì)算）。然而，相比于其他點(diǎn)，微震破裂源點(diǎn)的ξi,j的絕對(duì)值最小。因此，微震震源定位的實(shí)現(xiàn)可通過搜索節(jié)點(diǎn)，使得與其對(duì)應(yīng)的多傳感器的殘差值最小。換言之，在給定波速模型中，首先使用二階快速行進(jìn)法計(jì)算所有網(wǎng)格節(jié)點(diǎn)到每一個(gè)微震傳感器的Δt，之后微震震源點(diǎn)的確定，可通過搜索使得不同傳感器的殘差值達(dá)到最小的網(wǎng)格節(jié)點(diǎn)來實(shí)現(xiàn)。此處，我們基于最小二乘的思路，使用目標(biāo)函數(shù)f來量化殘差，如式（23）所示。

式中，(x,y,z)是波速模型中任意網(wǎng)格節(jié)點(diǎn)的坐標(biāo)。微震破裂源點(diǎn)可使得目標(biāo)函數(shù)f達(dá)到最小值。

至此，整個(gè)定位過程可被總結(jié)如下：

步驟1：建立波速模型，并對(duì)相應(yīng)的網(wǎng)格節(jié)點(diǎn)賦予相應(yīng)的P波波速值。

步驟2：使用二階快速行進(jìn)法計(jì)算所有網(wǎng)格節(jié)點(diǎn)到每一個(gè)微震傳感器的理論到時(shí)。

步驟3：基于實(shí)際到時(shí)，計(jì)算目標(biāo)函數(shù)值。

步驟4：搜索幾個(gè)最小目標(biāo)函數(shù)f值所對(duì)應(yīng)的網(wǎng)格節(jié)點(diǎn)。其坐標(biāo)的平均值被視作微震震源點(diǎn)的位置。在本文中，取10個(gè)最小目標(biāo)函數(shù)值及其所對(duì)應(yīng)的網(wǎng)格節(jié)點(diǎn)來實(shí)現(xiàn)定位。

步驟5：從微震破裂源點(diǎn)到每一個(gè)微震傳感器之間的射線路徑，通過2.2節(jié)中的方法獲得。

3. 數(shù)值實(shí)驗(yàn)

本節(jié)將進(jìn)行多組數(shù)值實(shí)驗(yàn)來展示二階快速行進(jìn)法和射線路徑追蹤的性能。此外，數(shù)值定位實(shí)驗(yàn)被用于驗(yàn)證定位方法的合理性。

3.1. 二階快速行進(jìn)法

在這一部分，將進(jìn)行一階與二階快速行進(jìn)法的對(duì)比實(shí)驗(yàn)，以展示初至到時(shí)的表現(xiàn)。所采用的波速模型為100 m × 100 m × 100 m（長×寬×高），網(wǎng)格間距為1 m。微震破裂源點(diǎn)位于(0, 0, 0)，并且所有網(wǎng)格節(jié)點(diǎn)的P波波速值為4000 m·s-1。兩節(jié)點(diǎn)之間的距離除以P波波速，所得的值可作為兩節(jié)點(diǎn)之間到時(shí)的解析解。數(shù)值解通過一階與二階快速行進(jìn)法計(jì)算獲得。此外，ξti被用于量化計(jì)算誤差，其定義如式（24）所示：

式中，ξti為第i個(gè)節(jié)點(diǎn)處的誤差；tNi與tAi分別為第i個(gè)節(jié)點(diǎn)處的數(shù)值解與解析解。

計(jì)算結(jié)果被展示在圖3中。圖3（a）展示了在波速模型中解析解的初至到時(shí)。計(jì)算誤差被分別展示在圖3（b）、（c）中。顯然，由于采用了高階差分的方式，二階快速行進(jìn)法有更好的表現(xiàn)，并且能夠有效地減小計(jì)算誤差。在圖4中，盒形統(tǒng)計(jì)圖被進(jìn)一步用于量化計(jì)算誤差分布。與一階快速行進(jìn)法相比，二階快速行進(jìn)法計(jì)算誤差更小，且分布在較小的范圍內(nèi)。一階快速行進(jìn)法的誤差中值（4.1 × 10-4s）是二階快速行進(jìn)法相應(yīng)的誤差中值（1.0 × 10-4s）的4倍，并且一階快速行進(jìn)法的下四分位數(shù)誤差（3.15 × 10-4s）也明顯大于二階快速行進(jìn)法。由此可見，二階快速行進(jìn)法能夠獲得更準(zhǔn)確的計(jì)算結(jié)果。

3.2. 射線路徑追蹤

在這一小節(jié)中，建立一個(gè)雙層波速模型與一個(gè)含空洞波速模型來驗(yàn)證基于二階快速行進(jìn)法計(jì)算結(jié)果的射線路徑追蹤方法的合理性。

圖3. 一階與二階快速行進(jìn)法在100 m × 100 m × 100 m的均勻波速模型中的計(jì)算結(jié)果。（a）解析解結(jié)果；（b）一階快速行進(jìn)法結(jié)果；（c）二階快速行進(jìn)法結(jié)果。

圖4. 采用一階與二階快速行進(jìn)法獲取結(jié)果誤差的盒形統(tǒng)計(jì)圖。一階與二階快速行進(jìn)法計(jì)算誤差的離群值分別為9247和2701。Q1：上四分位誤差值；Q3：下四分位誤差值；IQR：Q3～Q1。

雙層的波速模型網(wǎng)格間距為1 m，模型大小為200 m × 200 m × 200 m（ 長×寬×高），0～100.5 m高的P波波速為6000 m·s-1，101～200 m高的P波波速為4000 m·s-1。微震破裂源點(diǎn)位于(100, 100, 0)，并且10個(gè)微震傳感器被安裝，其坐標(biāo)如表1所示。破裂源點(diǎn)與微震傳感器之間的傳播路徑如圖5所示。根據(jù)地震學(xué)的射線理論，在兩個(gè)區(qū)域交界面上，射線傳播應(yīng)該滿足斯奈爾定律（Snell law），如式（25）所示。

圖5. 在兩層波速模型中，從破裂源點(diǎn)到不同傳感器的射線路徑。

式中，θ1和θ2分別為入射角與折射角；v1和v2為兩個(gè)區(qū)域相應(yīng)的P波波速值。此處，v1和v2值分別取為6000 m·s-1和4000 m·s-1，故v1/v1=6000/4000=1.5。這里，Ri被用于量化射線路徑計(jì)算誤差，如式（26）所示。

式中，Ri對(duì)應(yīng)第i個(gè)微震傳感器的誤差；θ1i和θ2i分別為第i個(gè)微震傳感器所對(duì)應(yīng)的入射角與折射角。在表1中，R的最大值不超過1.5%，且大部分R值沒有超過1%，可忽略不計(jì)。可見，基于2.2節(jié)中的思路所獲取的射線路徑，在界面處滿足斯奈爾定律。

表1 Ri的誤差結(jié)果

圖6展示了含空洞的波速模型，其尺寸為200 m ×200 m × 200 m（長×寬×高），網(wǎng)格間距為1 m?？斩磪^(qū)域的P波波速設(shè)為340 m·s-1，模型剩余部分的波速設(shè)為5000 m·s-1?？斩磪^(qū)尺寸為25 m × 65 m（半徑×長度），圓柱中軸線經(jīng)過(50, 35, 50)與(50, 100, 50)兩點(diǎn)。三個(gè)微震破裂源點(diǎn)分別位于(45, 5, 50)、(45, 55, 95)、(70, 70,20)，四個(gè)微震傳感器分別布置在(25, 45, 65)、(25, 45,35)、(75, 45, 65)、(75, 45, 35)。所獲得的射線路徑如圖6所示。在離空洞較遠(yuǎn)的區(qū)域，射線路徑沿直線傳播，如同在均一波速模型中傳播；而在靠近空洞的區(qū)域，射線傳播路徑緊貼空洞區(qū)外側(cè)進(jìn)行傳播。這些結(jié)果展示了采用本文中的方法在含空洞區(qū)域依然能夠獲得合理的射線傳播路徑。

3.3. 數(shù)值定位實(shí)驗(yàn)

圖6. 三個(gè)微震破裂源點(diǎn)到四個(gè)微震傳感器的射線傳播路徑。（a）三維視圖；（b）右視圖；（c）前視圖；（d）俯視圖。

建立一個(gè)大小為383 m ×100 m × 121 m、網(wǎng)格間距為1 m并含空洞的波速模型，用于測試本文的定位方法（圖7）。這一模型基于地球物理學(xué)領(lǐng)域的Marmousi模型所建立。三個(gè)隧洞模型沿Y軸平行布置，尺寸為15 m × 100 m（半徑×長度），其中軸線分別過(75, 50,50)、(176, 50, 50)、(330, 50, 65)三點(diǎn)。由一階快速行進(jìn)法計(jì)算破裂源點(diǎn)到每個(gè)傳感器的到時(shí)，將其視作從P波中提取的實(shí)際到時(shí)，即式（23）中的ti與tj。所有的理論到時(shí)，也就是式（23）中的Δti與Δtj，采用二階快速行進(jìn)法計(jì)算獲得。微震震源定位基于式（23）實(shí)現(xiàn)，并且最終的定位結(jié)果為10個(gè)最小殘差值對(duì)應(yīng)的網(wǎng)格節(jié)點(diǎn)坐標(biāo)的平均值（圖7）。表2展示了所有的定位結(jié)果，這些結(jié)果相應(yīng)的誤差均小于4 m，表明該方法能夠在較為復(fù)雜的波速模型中實(shí)現(xiàn)微震破裂源點(diǎn)定位。

表2 在數(shù)值定位實(shí)驗(yàn)中的定位誤差

4. 在巖石工程中的應(yīng)用

4.1. 工程項(xiàng)目簡介

將本文的定位方法用于定位在猴子巖水電站建設(shè)期間開挖過程中所發(fā)生的爆破與微震事件。猴子巖水電站是典型的大跨度水電工程項(xiàng)目，該水電站建設(shè)在大渡河上，位于中國四川省成都市西南方約450 km。猴子巖水電站的地下廠房洞室群主要由主廠房、主變室、尾水調(diào)壓室以及母線洞組成。三個(gè)主要的地下洞室（包括主廠房、主變室、尾水調(diào)壓室）沿軸向方位北偏西61°（N61°W）平行布置，以高邊墻大跨度為主要特征。主廠房的開挖規(guī)模約為219.5 m × 29.2 m × 68.7 m（長×寬×高），其最小垂直與水平深度約為380 m與250 m。主變室長141.1 m、寬約18.8 m、高約25.2 m，而尾水調(diào)壓室長約60 m、寬約23.5 m、高約73.98 m。多種工程地質(zhì)勘查結(jié)果表明，一些結(jié)構(gòu)面包括小斷層與層間錯(cuò)動(dòng)帶穿過工程開挖區(qū)[38]。最大主應(yīng)力方向大致沿東西向，且地下廠房周圍圍巖主要為完整堅(jiān)硬的變質(zhì)灰?guī)r[39]。更多關(guān)于工程項(xiàng)目與現(xiàn)場工程地質(zhì)情況的詳細(xì)信息，請(qǐng)參考文獻(xiàn)[39]。

圖7. 使用二階快速行進(jìn)法在含空洞復(fù)雜波速模型中的微震破裂源定位。

4.2. 微震監(jiān)測系統(tǒng)構(gòu)建

為了監(jiān)測由于開挖所誘發(fā)的巖體中的微震活動(dòng)，一套由加拿大Engineering Seismology Group（ESG）公司生產(chǎn)的高精度微震監(jiān)測系統(tǒng)被安裝在地下洞室中（圖8）。ESG微震監(jiān)測系統(tǒng)主要由paladin數(shù)字信號(hào)采集系統(tǒng)、Hyperion數(shù)字信號(hào)處理系統(tǒng)，以及多個(gè)微震加速度傳感器所組成。所安裝的微震加速度傳感器的頻響范圍為50～5000 Hz，被布置在邊墻上的鉆孔末端，如圖8所示。獲取的微震信號(hào)通過Hyperion處理系統(tǒng)，以10 000 Hz的采樣頻率處理為數(shù)字信號(hào)。微震傳感器的坐標(biāo)分別為(41.70, 62.00, 1706)、(71.55, 62.20, 1706)、(97.80, 62.30, 1706)、(128.95, 56.80, 1703.5)?；诒茰y試與數(shù)字聲波測試的聯(lián)合反演結(jié)果，圍巖P波波速為5239 m·s-1。監(jiān)測信號(hào)數(shù)據(jù)通過ESG Wavevis軟件以“.txt”的形式導(dǎo)出，使用Python進(jìn)行處理分析。典型的微震波形通過本研究參與者手動(dòng)識(shí)別提取，并采用赤池信息準(zhǔn)則（Akaike information criterion, AIC）[40,41]進(jìn)行微震信號(hào)P波的初至到時(shí)提取。

4.3. 開挖所誘發(fā)的微震活動(dòng)

2013年12月5日至2014年1月16日，在主廠房與母線洞之間安裝的ESG微震監(jiān)測系統(tǒng)捕獲了大量微震事件。在此期間，開挖活動(dòng)采用鉆爆法結(jié)合機(jī)械開挖的形式進(jìn)行。圖9展示了在采集微震監(jiān)測數(shù)據(jù)之前，所完成的開挖的區(qū)域。開挖所誘發(fā)的微震事件主要分布于主廠房的下游邊墻一側(cè)，位于2#與3#母線洞之間?；诰徊ㄋ倌Ｐ偷亩ㄎ唤Y(jié)果被獲取，該結(jié)果忽略了空洞區(qū)域（即開挖區(qū)域）的影響。結(jié)果，一些微震事件被定位于空洞區(qū)域（圖9），這對(duì)開挖損傷區(qū)的圈定造成消極影響。

4.4. 基于二階快速行進(jìn)法的微震事件定位

4.4.1. 含空洞復(fù)雜波速模型構(gòu)建

地下洞室的含空洞波速模型由圖10所示的方式建立。首先，采用三維設(shè)計(jì)軟件（如Blender和AutoCAD）建立三維洞室模型，并且根據(jù)工程項(xiàng)目進(jìn)展，確定計(jì)算區(qū)域，確保模型包含所有的空洞區(qū)域。建立笛卡爾坐標(biāo)系，如圖10所示，使得Y軸方向平行于所有洞室的中軸線，并且Z軸的方向垂直向上。整個(gè)計(jì)算區(qū)域大小為162 m × 220 m × 86 m [圖10（a）]，并且被劃分成網(wǎng)格間距為1 m的正方體網(wǎng)格。所有網(wǎng)格節(jié)點(diǎn)的坐標(biāo)被記錄，并且基于現(xiàn)場爆破與聲波測試聯(lián)合反演的結(jié)果，將整個(gè)計(jì)算區(qū)域的P波波速設(shè)為5239 m·s-1。接下來，進(jìn)行布爾運(yùn)算（Boolean operation）以獲取開挖區(qū)域與整個(gè)區(qū)域的交集區(qū)域[圖10（c）]。之后，交集區(qū)域的節(jié)點(diǎn)三維坐標(biāo)被獲取，并依據(jù)對(duì)應(yīng)的坐標(biāo)，將空洞區(qū)節(jié)點(diǎn)的波速修改為340 m·s-1。

圖8. 猴子巖水電站地下洞室群布置及ESG微震監(jiān)測系統(tǒng)現(xiàn)場構(gòu)建。

圖9. 2#和3#母線洞之間的微震事件基于均一波速模型的定位結(jié)果。（a）微震事件的三維空間分布；（b）俯視圖；（c）前視圖；（d）左視圖。

圖10. 含空洞波速模型的建立。（a）圈定計(jì)算區(qū)域；（b）生成網(wǎng)格節(jié)點(diǎn)；（c）獲取開挖區(qū)所對(duì)應(yīng)的網(wǎng)格節(jié)點(diǎn)。

4.4.2. 爆破事件定位測試

在這一小節(jié)，基于監(jiān)測的信號(hào)數(shù)據(jù)，使用所提出的方法去定位記錄兩個(gè)鉆孔爆破事件。使用AIC方法提取所記錄的信號(hào)的初至到時(shí)。定位結(jié)果與誤差結(jié)果展示在表3中。兩個(gè)爆破事件的定位結(jié)果被分別降低了大約9 m和3 m，更接近于實(shí)際爆破位置?？梢郧宄乜吹剑捎诳紤]了空洞區(qū)，所提出的方法可以獲得較為準(zhǔn)確的定位結(jié)果。

4.4.3. 微震事件定位

使用本方法對(duì)2013年12月5日至2014年1月16日開挖活動(dòng)中所誘發(fā)的微震事件進(jìn)行定位。微震事件的空間分布如圖11所示。一個(gè)微震事件的射線路徑如圖11所示，展示了在含空洞的模型中P波的傳播路徑。射線在遠(yuǎn)離空洞區(qū)域沿直線傳播，而在接近空洞區(qū)域時(shí)，繞過開挖區(qū)域進(jìn)行傳播。相比于均一波速模型的定位結(jié)果（圖9），并沒有任何的微震事件被定位到空洞區(qū)域內(nèi)部，而主要是分布于主廠房的下游側(cè)邊墻，以及2#和3#母線洞之間的區(qū)域。聚集區(qū)1處的微震事件集中在主廠房下游側(cè)邊墻墻趾附近，標(biāo)高1680～1690 m（圖11），這由開挖活動(dòng)引起的應(yīng)力集中所致。聚集區(qū)2處的微震事件主要分布在母線洞2#和3#之間，標(biāo)高1698～1710 m，表明圍巖內(nèi)部存在開挖誘發(fā)裂隙。通過常規(guī)監(jiān)測技術(shù)（圖12）和現(xiàn)場噴射混凝土破裂（圖13）來驗(yàn)證本文定位方法的定位結(jié)果的合理性。原位多點(diǎn)位移計(jì)的位置如圖12（a）所示，位于2#母線洞上方，標(biāo)高1706.5 m。其孔口絕對(duì)位移過程圖表明，在2014年1月1日至2014年1月9日期間，位移呈上升趨勢[圖12（b）]，其間微震事件數(shù)目呈明顯增加的趨勢。不同區(qū)段的相對(duì)位移圖[圖12（c）]進(jìn)一步表明，多點(diǎn)位移計(jì)的第24 m測點(diǎn)與尾端固定點(diǎn)之間的區(qū)段發(fā)生了明顯的變形，其位置靠近微震事件的聚集區(qū)2（圖11和圖12）。此外，如圖13所示，2#母線洞側(cè)壁的噴射混凝土開裂和剝落進(jìn)一步表明了原位巖體損傷與微震事件聚集區(qū)2之間的關(guān)系，驗(yàn)證了所提出的微震定位方法的可行性。

表3 使用不同的速度模型定位爆破事件的結(jié)果

5. 結(jié)論

準(zhǔn)確的微震事件定位是微震監(jiān)測技術(shù)圈定圍巖損傷區(qū)的關(guān)鍵。現(xiàn)場復(fù)雜情況會(huì)對(duì)最終微震破裂源定位結(jié)果造成消極影響。在本文中，一種基于微震信號(hào)P波理論到時(shí)與實(shí)際到時(shí)之間殘差的微震定位方法被提出，用于含空洞復(fù)雜區(qū)域微震破裂源定位?；诙A差分與窄帶技術(shù)，二階快速行進(jìn)法通過求解程函方程來獲取P波的理論到時(shí)。在獲得理論到時(shí)的基礎(chǔ)上，微震破裂源到微震傳感器的射線路徑可采用線性插值和龍格-庫塔方法求解獲得。利用AIC方法自動(dòng)提取微震信號(hào)的實(shí)際到時(shí)后，通過搜索使得目標(biāo)函數(shù)值達(dá)到最小的網(wǎng)格節(jié)點(diǎn)，實(shí)現(xiàn)微震破裂源定位。

圖11. 位于2#和3#母線洞之間的微震事件，基于含空洞速度模型的定位結(jié)果。（a）微震事件的三維空間分布；（b）俯視圖；（c）前視圖；（4）左視圖。

圖12. 多點(diǎn)位移計(jì)3-8的測量結(jié)果，以及日均微震事件數(shù)目。（a）多點(diǎn)位移計(jì)3-8的位置；（b）3-8的鉆孔位移量與日均微震事件數(shù)；（c）3-8在不同分段的相對(duì)位移量。

圖13. 2#母線洞邊墻處噴射混凝土破裂與剝落。

本文所提出的方法通過數(shù)值實(shí)驗(yàn)得到驗(yàn)證，并應(yīng)用于猴子巖水電站三大洞室的開挖。對(duì)于所記錄的爆破事件的定位結(jié)果表明，與均一波速模型的定位相比，本文所提出的方法可以有效地減小定位誤差。此外，通過安裝多點(diǎn)位移計(jì)獲得的現(xiàn)場監(jiān)測結(jié)果以及2#母線洞側(cè)壁的噴射混凝土裂縫和剝落情況，驗(yàn)證了本方法對(duì)開挖誘發(fā)微震事件的定位結(jié)果的合理性。該方法可用于含空洞復(fù)雜環(huán)境下的微震事件定位，可有效地幫助圈定圍巖內(nèi)部的損傷區(qū)域。

致謝

感謝國家自然科學(xué)基金（52039007）提供的資金支持。

Compliance with ethics guidelines

Ruochen Jiang, Feng Dai, Yi Liu, and Ang Li declare that they have no conflict of interest or financial conflicts to disclose.