董載天 王軍民

摘 ?????要： 由于隧道超前探測的環(huán)境與常規(guī)石油地震勘探的差異，不宜直接將過去的模擬方法拿到隧道中進行應用。結合層狀地層表現(xiàn)出橫向各向同性的特性，推導出在隧道空間適用的二維VTI介質一階qP波方程，分析了其穩(wěn)定條件。建立斷層和軟弱夾層的幾種數(shù)值模型，利用PML邊界解決模型邊界反射，采用交錯網(wǎng)格進行有限差分正演模擬，得到其共炮點道集。結果表明，qP波有限差分法可以快速地模擬不良的地震響應，規(guī)避偽S波和數(shù)值不穩(wěn)定問題，為隧道超前探測資料的解釋提供依據(jù)，提高預測的準確性。

關 ?鍵 ?詞：隧道超前探測;VTI介質;qP波;PML

中圖分類號：TV221.2???????文獻標識碼：?A ??????文章編號： 1671-0460（2020）03-0683-07

Forward Modeling of Tunnel Advanced

Detection qP?Wave Based on VTI Medium

DONG Zai-tian， WANG Jun-min

（Yangtze University， Hubei Wuhan?430100， China）

Abstract： ?Due to the difference between the environment of tunnel pre-detection and conventional petroleum seismic exploration， it is not appropriate to directly apply the past simulation method to the tunnel. Combined with the transverse isotropy?characteristic?of layered stratum， the first-order qP wave equation of the two-dimensional VTI medium suitable for tunnel space was?derived， and the stability conditions were?analyzed. Several numerical models of fault and weak interlayer were?established. The boundary reflection of the model was?solved by PML boundary， and the finite difference forward modeling was?carried out by staggered grid to obtain the common shot gather. The results show that the qP wave finite difference method can quickly simulate the bad seismic response， avoid the pseudo S wave and numerical instability problem， and provide a basis for the interpretation of tunnel advance detection data?to?improve the accuracy of prediction.

Key words： ?tunnel advance detection; ?TI media;??qP wave; ?PML

在我國經(jīng)濟高速發(fā)展的當下，大量公路鐵路工程都有隧道施工。在隧道工程施工過程中，由于掌子面前方地質狀況不詳，面對破碎帶、裂隙和溶洞的支護不當或支護不及時，經(jīng)常出現(xiàn)坍塌、突水突泥、冒頂?shù)戎卮鬄暮?，造成嚴重的人員傷亡和財產(chǎn)損失。顯然在隧道工程中，使用超前探測技術了解掌子面前方地質狀況，對工程施工安全極其重要。超前探測方法從原理上可以劃分為：地震波法、電磁法（瞬變電磁、電阻率、激發(fā)極化、核磁共振）、巖體溫度法等，其中地震波法在遠距離勘探中有較好的分辨率，成為目前隧道超前探測中的常用方法之一。在隧道地震波法探測中，常用的方法有VSP、HSP^[1]、TSP^[2]、TVSP^[3]、陸地聲納^[4]、HSP^[5]、TRT^[6，7]、TST^[8]、ISIS^[9]等方法。由于隧道施工現(xiàn)場的空間狹小且工程工期壓力大，所以隧道超前探測技術必須滿足狹小場地條件下進行大距離、高精度、高效率探測要求^[10]，這就要求在對傳統(tǒng)方法技術（測線與探測對象平行）進行改進（測線與探測對象垂直）的同時，進一步研究地震波在隧道復雜條件下的傳播機理特性^[11]。因此，本文在層狀地層條件下，構建斷層、軟弱夾層等不良地質災害體模型，從一階速度?應力qP波方程出發(fā)，采用交錯網(wǎng)格高階有限差分進行隧道復雜波場的正演模擬計算，以充分認識不同條件下地震波響應特征。

1 ?TI介質的波動方程

在我國，層狀地層分布廣泛，約占陸地面積的77.3%^[12]，在隧道施工過程中大多要面對這樣的地質環(huán)境，而地下層狀在橫向各向同性地層往往表現(xiàn)出橫向各向同性的物理響應特征。介于筆者使用波動方程法進行地震數(shù)值模擬，所以推導出橫向各向同性（TransverseIsotropy）介質的波動方程尤為重要。

1.1 ?波動方程導出

波動方程的建立涉及應力、應變和位移這三個物理量，需要使用彈性體的本構方程、運動微分方程和柯西方程。由本構方程，可知彈性體應力與應變關系表達式：

????????（1）

式中：?— 應力張量;

e?— 應變張量;

C?— 彈性常數(shù)。

彈性常數(shù)共有81個，由柯西方程可知，應變張量有對稱性：

???????????????（2）

其中：u、v —位移分量;

x，?y?—笛卡爾坐標。

而應力張量也有此種對稱性，即。因此彈性常數(shù)減少為36個，同時彈性常數(shù)張量相對于主對角線對稱，有。使用單下標替代雙下標11→1，22→2，33→3，23（32）→4，13（31）→5，12（21）→6，則可的表達式：

??????（3）

上式可簡記為：

?????????????（4）

從彈性體的角度，描述完應變與應力關系、應變與位移關系后，通過牛頓第二定律可建立運動平衡微分方程：

??????????????（5）

式中：—?彈性體密度，kg/m;

U?— 位移張量;

— 應力張量;

t — 時間，s;

F —體力張量;

L —偏導算子。

偏導算子L為：

???????（6）

利用Cauchy方程可將（5）式表達為：

?????????（7）

1.2 ?TI介質中qP波波動方程

地表巖土物理性質極其復雜不便于數(shù)值模擬研究，地球物理學家為了能借助波動物理研究其波傳播特性，根據(jù)各項異性彈性體介質的基本對稱性進行了分類。其中三方、六方各項異性介質常叫做橫向各向同性（Transverse Isotropy，簡稱TI）介質，是目前實際應用研究中使用最廣泛的模型。TI介質又根據(jù)對稱軸的水平和垂直細分為，VTI介質和HTI介質（圖1）。

兩介質都需要5個獨立的彈性參數(shù)來描述，其彈性矩陣分別為：

分別將兩式帶入公式（7），可得VTI介質三維波動方程：

由于VTI介質與地下實際介質的響應特質較為一致且研究最為廣泛，所以筆者主要使用VTI介質的彈性波波動方程，因此此處不再給出HTI介質的波動方程。

在利用此方程進行TI介質的彈性波正演模擬中，波動方程的包含三個分量，占用計算機大量內存，同時由于縱波速度快于橫波，而網(wǎng)格的空間步長和采樣時間步長主要受最波的小速度影響，為穩(wěn)定計算，只能減小空間和采樣的間隔，使得計算的時間點數(shù)與空間點數(shù)大幅提高，進一步減慢了計算速度，不利于計算機模擬。且目前實際勘探仍以縱波為主，多分量數(shù)據(jù)中進行橫波分離仍然是個難點問題。

因此筆者借用Alkhalifah提出的聲學假設條件，直接將橫波速度設為零，從而在復雜的 VTI 介質彈性波方程中分裂出qP波方程^[14]。該方程不產(chǎn)生橫波，雖然現(xiàn)實中上不可實現(xiàn)，但在運動學上能對彈性波方程進行有效地近似，便于研究TI介質中qP波的傳播規(guī)律。

根據(jù)Alkhalifah的假設導出的三維 VTI 介質中的qP波方程為：

???（11）

式中：F?— 為波場函數(shù);

t?— 為時間，s;

x、y、z?— 為三維空間坐標;

v 、v_v?— qP波的動校正速度和垂向速度，m/s;

η?— 各向異性參數(shù)。

已知：e、d 為 Thomsen 參數(shù)^[15]。

由于隧道超前探測觀測方向為掌子面正前方，而常規(guī)地震勘探觀測方向為地面垂向（如圖2），所以筆者后面推導VTI介質中的波動方程時只保留xy分量。

因此，由式（11）可得二維VTI介質中qP波方程：

直接求解上式需要消耗大量內存和時間，所以利用，可將上式轉化為：

其中為應力。借鑒Hestholm推導二維VTI介質中qP波方程的一階應力-速度方程形式的方法，可得：

式中：— 彈性體密度，kg/m;

v_x、v?_y— x、y方向速度，m/s;

?— 計算過程中的中間變量。

2 ?交錯網(wǎng)格

在實際計算過程中必須對介質模型進行離散，將介質模型網(wǎng)格化。實際應用過程中，交錯網(wǎng)格相比規(guī)則網(wǎng)格有更好的差分算子局部特性，更高的精度和更好的收斂性。所以筆者采用交錯網(wǎng)格（如圖3）將放置在網(wǎng)格節(jié)點。

利用Taylor級數(shù)展開法，用網(wǎng)格節(jié)點上數(shù)值的差商近似代替前式中的導數(shù)。推導出了時間域 2 階、空間域2N階精度的交錯網(wǎng)格離散差分格式。

式中：?— 空間間隔，m;

— 時間間隔，s;

i、j、k— 空間和時間節(jié)點;

U、V— 速度分量v_x和v_y的離散值，m/s;

P、Q、R— 應力和的離散值。

3 ?有限差分模擬中的重要問題

3.1 ?PML邊界條件

針對隧道超前探測的彈性波數(shù)值模擬而設計的數(shù)據(jù)模型，僅考慮深埋隧道時，隧道空間為自由邊界，計算區(qū)域外部邊界則都是人工截斷邊界。人工邊界會使向外傳播的彈性波完全反射回去，從而使正演模擬計算結果與實際情況大不符。為了盡量降低人工截斷邊界所產(chǎn)生的影響，引入Berenger針對電磁波傳播提出了完美匹配層（Perfect Matched Layer，簡稱 PML）吸收邊界條件，削弱人工截斷邊界的反射，改善正演模擬效果。

根據(jù)PML的推導思路，可給出式（14）的吸收邊界方程：

式中：d（x）、d（y） — 阻尼因子;

R?— 反射系數(shù);

?— 完美匹配層吸收邊界的厚度，m;

h?— 完美匹配層吸收邊界的深度，m;

v?— 匹配層中的彈性波波速，m/s。

3.2??穩(wěn)定性

波動方程的差分格式按照時間進行遞推，前一時刻的波函數(shù)數(shù)值解的舍入誤差必然會對下一時刻的波場產(chǎn)生影響。這就需要對誤差傳播情況進行估計，使其不會隨著時間推移而快速增長，以致影響數(shù)值模擬的效果進而破壞計算，使模擬無法繼續(xù)進行。筆者此處借鑒Virieux在1986年提出的推導方法，得到二維VTI介質中二階時間差分精度、2N階空間差分精度的穩(wěn)定條件：

??????（17）

當中表示模型中最大速度的平方。當模型橫縱空間步長一致時，亦可將方程改寫為：

?????????????（18）

3.3??震源加載

此部分還有一個重要問題就是加載震源。本文為保證成像的分辨率，選擇使用零相位雷克子波作為震源。其時間域表達形式為：

????????（19）

式中：t₀?—初始時刻，ms;

f?—主頻， Hz。

4 ?數(shù)值模擬實例

建立隧道空間的模型（如圖4），正演模型的尺寸為300 m×100 m，模型四周使用PML吸收邊界。

隧道尺寸為60 m×8 m，空腔中充填空氣，波速340?m/s;?圍巖波速3 600 m/s。采樣間隔為0.2 ms，采樣時長為0.1 s。

觀測方式如所示，掌子面中心單震源，左右兩側墻壁各布置4個檢波器，間隔2 m，水平最小炮間距10 m。震源與檢波器皆置于墻內，埋深2 m。

通過正演模擬可得到地震響應記錄（圖6），為共炮點道集（CSP）?？梢钥闯鲋边_波能量很強，無邊界反射和偽S波干擾。根據(jù)費馬原理計算的初至時間在4 ms左右，基本與記錄吻合。

4.1 ?不良地質體模型及其地震響應

隧道施工過程中會遇到很多危害施工安全的不良地質體，常見的有：斷層破碎帶、軟弱夾層等。由于二維平面不存在傾角，所以后續(xù)模型中需要關注不良地質體的走向。

4.1.1 ?斷層模型及響應

斷層破碎帶是由于斷層上下盤搓動引起斷面附近的巖體發(fā)生破碎，形成破碎帶。其穩(wěn)定性較差，容易引發(fā)隧道塌方;在地下水富集區(qū)域，破碎帶變成了一條良好水通道，極易引發(fā)突水突泥災害。由于破碎帶與斷層的伴生特性且破碎帶本身理化性質過于復雜，筆者此處將斷層破碎帶的模型設置為斷面兩側阻抗不同，破碎帶并不在模型當中直接體現(xiàn)。數(shù)值模型構建時，首先選擇走向與隧道方向垂直，距掌子面80和120 m的斷層，斷層右盤波速為4 000 m/s（圖7）。

對比兩模型的CSP道集（如圖8）可以看出，當斷面與掌子面距離減小時，阻抗界面的反射波逐漸靠近強振幅的直達波，初至有可能被遮蔽，此時僅靠時域記錄是難以判斷前方斷面距離的;當兩者間距拉大，反射波整體向后時移，初至明顯但振幅減弱，現(xiàn)場施工噪音較強時，會難以見到完整波形。

4.1.2??軟弱夾層及響應

軟弱夾層成因復雜且尚無統(tǒng)一的分類原則，此處只針對沉積型軟弱夾層進行分析。此種軟弱夾層往往形成與河流沉積環(huán)境，層中粉砂巖、泥巖含量偏高，強度低質地較軟弱^[17]。當構造作用使夾層被剪切、破碎時，無論區(qū)域是否富水，都會產(chǎn)生低速低密度的軟弱夾層。軟弱夾層的數(shù)值模型參考斷層模型，將夾層設置在掌子面前方80 m處且走向與隧道方向垂直，此時一次反射波不會被直達波遮蔽。分別構建10、20 m厚的夾層，波速為800 m/s（如圖9）。

此時獲得CSP道集（如圖10）。由于此時夾層與圍巖阻抗差異過大，且檢波器內置在墻壁中，導致阻抗界面反射波振幅較強。對于薄夾層來說，子波在夾層中旅行時較短，夾層右側依舊為高阻抗差界面，所以10 m薄夾層的地震反射記錄上出現(xiàn)多次波。而層厚為20 m的夾層，子波旅行時被拉長此時反射波沒有在記錄時間內返回到檢波器，故道集上無明顯反射波。

將軟弱夾層與掌子面之間距離拉大至120 m（模型如圖11）。

從兩模型CSP道集（如圖12）可以看出，此時兩種模型的地震響應頗為一致，顯然夾層右界面的反射波由于旅行時較長沒有在記錄時長內返回到檢波器位置。

構建一個右盤為低速巖體的斷層（如圖13），此時獲得其CSP道集（如圖14），其與掌子面距離為120 m的軟弱夾層的地震響應基本一致?？梢姽潭ú蓸訒r長有時無法分辨的不同的不良地質體，此時需要增加采樣時長或者利用其他地質勘探資料來確定不良地質體形態(tài)。

5??結 論

隧道超前探測尚處于發(fā)展階段，由于其探測空間的復雜以及其與傳統(tǒng)地震勘探的差異，使得隧道超前探測在實際應用過程中存在一些問題。筆者從彈性體基本理論出發(fā)，逐步推導出適用于隧道超前探測二維正演模擬公式，并模擬了斷層、軟弱夾層的地震響應。從中得到一些結論：

（1）利用三維VTI介質波動方程，根據(jù)隧道超前探測的觀測方式，推出了適用于隧道勘探正演模擬的二維VTI介質qP方程，后導出其1階應力-速度方程，進行時間2階空間2N階的有限差分正演模擬。其計算迅速，規(guī)避了偽S波干擾和計算不穩(wěn)定性，便于快速得到不同地質模型的地震響應特征。

（2）利用PML吸收邊界，處理了隧道后方人工邊界引起的反射。采用交錯網(wǎng)格和，使方程能較好地模擬斷層和軟弱夾層的地震響應，減少頻散。

（3）二維平面模擬中，采樣時長與斷層、軟弱夾層的物性參數(shù)組合不當，會導致兩種不同的不良地質體產(chǎn)生相同的地震響應。在數(shù)值模擬時，可改變模型參數(shù)，使其顯現(xiàn)不同的地震響應。而在現(xiàn)場實際情況來說，需要現(xiàn)場工程師結合其他地質勘探資料綜合分析。

參考文獻：

[1]柳楊春. HSP地質超前預報技術及其應用[J]. 西部探礦工程， 1997 （05）： 40-42.

[2]劉志剛， 劉秀峰. TSP（隧道地震勘探）在隧道隧洞超前預報中的應用與發(fā)展[J]. 巖石力學與工程學報， 2003， 22?（8）：?1399-1402.

[3]曾昭璜. 隧道地震反射法超前預報[J]. 地球物理學報， 1994 （2）.

[4]鐘世航. 陸地聲納法的原理及其在鐵路地質勘測和隧道施工中的應用[J]. 中國鐵道科學， 1995 （4）： 48-55.

[5]Inazaki T， Isahai H， Kawamura S， et al. Stepwise application of horizontal seismic profiling for tunnel prediction ahead of the face[J]. The Leading Edge， 1999， 18（12）： 1429?1431.

[6]Otto R， Button E， Bretterebner H， Schwab P.The application of TRT-trae reflection tomography-at the Unterwald Tunnel[J]. Geophysics，?2002， 20 （2）： 51-56.

[7]Yamamoto T，?Shirasagi S，?Yokota Y，?Koizumi Y.?Imaging geological conditions ahead of a tunnel face using three-dimensional seismic reflector tracing system[J]. International Journal of the JCRM， 2011， 6（1）： 23-31.

[8]趙永貴. TST隧道超前預報技術及應用[C].中國地球物理學會第22屆年會論文集，2006.

[9]Borm G， Giese R， Klose C， et al. ISIS integrated seismic imaging system for the geologic prediction ahead in underground construction[C].65th Annual Conference and Exhibition， EAGE， Extended Abstracts. Norway， 2003： 887?899.

[10]宋先海， 顧漢明， 肖柏勛. 我國隧道地質超前預報技術述評[J]. 地球物理學進展， 2006 （02）： 605-613.

[11]張平松，吳健生.中國隧道及井巷地震波法超前探測技術研究分析[J]. 地球科學進展， 2006 （10）： 1033-1038.

[12]張學民. 巖石材料各向異性特征及其對隧道圍巖穩(wěn)定性影響研究[D]. 中南大學， 2007.

[13]程玖兵， 康瑋， 王騰飛. 各向異性介質qP波傳播描述Ⅰ： 偽純模式波動方程[J]. 地球物理學報， 2013， 56 （10）： 3474-3486.

[14]Alkhalifah， Tariq. An acoustic wave equation for anisotropic media[J]. Geophysics， 2000， 65（4）： 44-48.

[15]Thomsen， Leon. Weak Elastic Anisotropy[J]. Geophysics， 1986， 51：1954-1966.

[16]Clayton R， Engquist B. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bull. seismol. soc. amer， 1977， 67?（6）：?1529-1540.

[17]胡濤，任光明，聶德鑫， 等. 沉積型軟弱夾層成因及強度特征[J]. 中國地質災害與防治學報，2004， 15 （1）： 124-128.