基于多源地震干涉法隧道斷層位置預測數(shù)值模擬

范占鋒， 蔡建華， 趙 偉

(1.成都大學建筑與土木工程學院， 成都 610106； 2.中鐵西南科學研究院有限公司， 成都 611731； 3.中鐵十八局集團有限公司， 天津 300222)

隨著中國的交通建設和能源開發(fā)不斷發(fā)展，特別是西部山區(qū)在建或即將修建的高鐵、公路以及水利工程都避繞不開無數(shù)大小斷層破碎帶、高地應力以及大變形區(qū)域[1-4]。在高地應力地區(qū)進行地下工程建設，如何提高隧道超前地質(zhì)預報的準確性是當前工程界面臨的一大關鍵難點。

從20世紀90年代至今，伴隨著中國大型基礎設施的建設，如大型水電站引水隧洞、公路隧道、鐵路隧道等的開挖，隧道超前地質(zhì)預報技術在中國隧道施工中也逐漸趨于完善，預報精度不斷提高，預報方法趨于多樣化。其中彈性波反射法應用最多，其次是電磁波反射法、瞬變電磁法、激發(fā)激化法和核磁共振法等。這些預報方法對不良地質(zhì)體的識別側重點不同。為了提高預報精度，很多專家學者提出了綜合超前地質(zhì)預報方法[5-8]，主要采用長短搭配組合方法。

然而，當隧道賦存地應力時，由于地應力的存在會使得彈性波在節(jié)理巖體中的傳播有一定的影響，進而影響預報精度。趙明階[9]推導了節(jié)理巖體在均勻應力場作用下的彈性波傳播速度和衰減隨應力的變化關系，并通過試驗進行了驗證。劉高等[10]研究了金川礦體的巖體力學參數(shù)及時空演化規(guī)律時指出地應力是影響巖體結構的主要因素，最終可通過彈性波波速來綜合反映。田家勇等[11]基于巖石的聲彈理論研究了不同初始應力條件下沿不同方向的彈性波波速，表明巖石波速和應力之間存在特定的關系。杜存蒼[12]采用數(shù)值模擬軟件FLAC3D分析了籃家?guī)r隧道高地應力段圍巖變形規(guī)律，但未對高地應力段的不良地質(zhì)體如何預報進行解釋說明。近年來，一些學者還嘗試采用新的預報技術，吳豐收等[13]認為目前的隧道超前地質(zhì)預報方法距離短、頻次多，占用施工時間并存在安全隱患，基于互相關地震干涉法理論，提出了將多源地震干涉法(multisource seismic interferometry)應用于隧道超前地質(zhì)預報。汪旭等[14]認為不論是隧道盾構法施工還是鉆爆法施工，其實際震源很復雜并且施工環(huán)境會對信號采集質(zhì)量有影響，從數(shù)值模擬角度探索了將多源地震干涉技術應用于隧道超前地質(zhì)預報。多源地震干涉技術能夠對接收到的透射波地震記錄進行相互干涉運算，能夠從透射波信號中提取隧道前方不良地質(zhì)體的反射波，實現(xiàn)對前方地質(zhì)條件做出判斷的功能，這表明多源地震干涉技術在隧道超前預報領域具有廣泛的應用前景[15-16]。

可見，隧道超前地質(zhì)預報實施效果除了與目前本身技術水平有關外，還與巖體的應力狀態(tài)關系極大。因此，主要基于多源地震干涉法，模擬不同圍巖波速和不同斷層傾角對斷層位置識別的影響，通過數(shù)值模擬建立地應力和斷層傾角與斷層位置識別的修改關系式，進而對白馬隧道高地應力段不良地質(zhì)體的預報結果進行修正。

1 多震源地震波干涉法原理

1.1 多源地震

多源地震的波動方程可表示為[6]

?2u+Sbl(t,X)

(1)

式(1)中：u和v分別為位移和速度；?2為Laplace算子；s(t)δ(X-X0)為震源函數(shù)；X0為震源位置；Sbl(t,X)為多震源波動方程的震源項；Γ(X)為Gamma函數(shù)，表示多震源混合矩陣算子，可表示為

(2)

式(2)中：γ1,γ2,…γn表示第1～n個震源地震記錄；Tn(X)為位置X處震源n的延遲激發(fā)時間；j為虛數(shù)單位；ω為圓頻率。

從式(1)可以看出，多震源激發(fā)環(huán)境下的地震記錄是由單個震源地震記錄的線性疊加所組成。

1.2 地震波干涉法

地震波干涉法是通過對記錄到的地震信號進行干涉得到新的地震信號。新地震信號不僅包含了原始地震信號的特性，也包含了原始信號所不具有的某些重要特征[17]。目前，該技術還處于理論階段，距離實際應用還有一定差距，原因在于實際工程中震源復雜，產(chǎn)生振動的頻率和能量等難以識別[18]。同時，實際工程中隧道所處的地質(zhì)條件復雜對成像也會造成一定的困難。采用數(shù)值模擬手段，對影響成像結果的主要因素進行分析研究，有利于該技術的實際應用。

Claerbout最早通過對水平層狀介質(zhì)模型研究，利用自相關運算將透射波記錄和反射記錄建立了有效聯(lián)系，即透射波地震記錄的自相關等價于自激自收模式記錄，Wapenaar等[19]使用格林定理嚴格地證明了Claerbout的研究，為地震干涉法奠定了堅實的數(shù)學和物理基礎。將透射波記錄和反射波記錄建立有效聯(lián)系，形成對隧道掌子面前方異常體的目標成像，互相關干涉法原理可表示為

(3)

式(3)中：T(xA,xi, -t)和T(xB,xi,t)分別為從地下震源xi傳到A檢波點和B檢波點的透射地震記錄；xA和xB分別為A檢波點和B檢波點的位置；t為時間；δ為Dirac函數(shù)；R(xA,xB,t)為以B為震源，A為檢波點的反射地震記錄；*表示卷積；相關型地震波干涉法提取的地震信號既包含因果部分R(xA,xB,t)，也包含非因果部分和零時刻的脈沖響應R(xA,xB, -t)。

2 數(shù)值計算模型及驗證

2.1 計算模型

在深埋隧道建設中，地應力使斷層破碎帶或節(jié)理裂隙密集發(fā)育帶巖體被壓密、節(jié)理閉合，地震波在此類介質(zhì)中傳播速度增加[20-21]。據(jù)此，建立尺寸為長×寬=400 m×400 m的二維平面模型，如圖1所示，隧道剖面寬10 m，信號激發(fā)端和接收端隧道洞段分別開挖50 m。在隧道中軸線上，沿水平方向150 m處設置1條斷層，斷層距離信號接收端掌子面為100 m，貫穿整個模型，斷層與掌子面平行，如圖1(a)所示。同時設置幾種傾斜斷層，如圖1(b)所示，具體計算參數(shù)如表1所示。

圖1 發(fā)育1條斷層的二維平面數(shù)值計算模型Fig.1 Two-dimensional plane numerical calculation model for developing a fault

表1中，假設隧道圍巖波速Cs按照200 m/s的增量從4 000 m/s增長到4 600 m/s，表示地應力使巖體波速增加。斷層傾角θ按照15°增量從0°變化到45°，斷層波速設為1 800 m/s。利用有限差分交錯網(wǎng)格法，通過MATLAB編程，分析圍巖波速和斷層傾角兩因素對斷層位置識別的影響。

表1 數(shù)值計算模型參數(shù)Table 1 Parameters of the numerical simulation model

2.2 算法驗證

首先驗證所用數(shù)值模擬算法是否正確。假設激發(fā)震源函數(shù)采用典型的雷克(Ricker)子波，主頻為35 Hz，震源激發(fā)掌子面位于模型左側，共布設10個震源，震源間距均為1 m；透射波信號接收點位于右側，共布設21個檢波器，檢波器間距均為0.5 m。斷層傾角θ=0°和圍巖波速Cs=4 000 m/s時的數(shù)值模型速度分布如圖2所示。由圖2可知，圍巖波速分布均勻，在此基礎上進行正演計算。

圖2 斷層傾角0°和圍巖波速4 000 m/s時的數(shù)值模型 速度分布Fig.2 Numerical model velocity distribution when θ=0° and Cs = 4 000 m/s

斷層傾角θ=0°時正演計算的不同時間波場傳播分布如圖3所示。由圖3可知，地震波傳播時間從5 ms增加到16 ms時，由于圍巖為各向同性均質(zhì)體，觀察到半橢圓形的帶狀區(qū)域逐漸增大，沒有發(fā)生波的反射，即透射波未到達圍巖與斷層交匯處。當傳播時間增加到30 ms時，在設置的斷層位置處(150 m)發(fā)生了明顯的反射現(xiàn)象，透射波與反射波之間形成了“透鏡”型交匯區(qū)域。

圖3 斷層傾角為0°時5～30 ms的波場傳播分布Fig.3 The wave field propagation distribution from 5 ms to 30 ms when the fault dip angle is 0°

通過對原始透射波信號進行濾波、反褶積及互相關干涉處理后，可得到斷層傾角θ=0°和圍巖波速Cs=4 000 m/s時的地震干涉結果如圖4所示。根據(jù)楊為民等[6]的結論，該成像結果等價于自激自收剖面。因此，成像結果中的時間皆為雙程旅行時間，即在0.05 s附近區(qū)域存在異常界面。根據(jù)圍巖波速4 000 m/s，可計算出單程旅行時異常界面(斷層)到激發(fā)掌子面的距離為100 m，加上已開挖洞段50 m，實際預報斷層位置位于150 m，與模型中設置的斷層位置完全吻合，驗證了數(shù)值模擬算法的正確性。

圖4 斷層傾角0°和圍巖波速4 000 m/s時的地震干涉結果Fig.4 The result of the seismic wave interference when Cs = 4 000 m/s and θ = 0°

3 不同數(shù)值模擬工況計算

在驗證所用算法正確性的基礎上，根據(jù)表1中斷層傾角θ和圍巖波速Cs，分別考慮4種圍巖波速和4個斷層傾角的計算模型，(首先設置斷層傾角θ為0°，計算圍巖波速Cs為4 000、4 200、4 400、4 600 m/s的成像結果，依次類推，計算斷層傾角為0°、15°、30°和45°的情況，并假設斷層波速始終為1 800 m/s)，其目的是主要考慮不同圍巖波速及斷層與掌子面的夾角關系對最終成像結果的影響，進而分析地應力對預報斷層位置的影響。

3.1 工況1：斷層傾角為0°

首先模擬斷層與掌子面平行，即斷層傾角θ=0°時的情形。圍巖波速Cs=4 000 m/s時的波場速度分布如圖2所示。其他3種圍巖波速的波場速度分布與之類似。不同圍巖波速條件下傳播30 ms后的波場傳播如圖5所示。可觀察到在斷層位置發(fā)生了明顯的透、反射現(xiàn)象，透射波傳播形態(tài)基本相似，在反射波和透射波之間都形成了一個“透鏡”型區(qū)域，反射波與透射波傳播方向在一條直線上。所不同的是，從圖5可以看出，隨著圍巖波速從4 000 m/s增加到4 600 m/s，“透鏡”型區(qū)域面積在不斷增大，其原因在于圍巖波速增加使得從斷層位置處反射回來的波速傳播加快所造成。

通過對原始透射波信號進行濾波、反褶積及互相關干涉處理后，斷層傾角θ=0°時不同圍巖波速的地震干涉結果如圖6所示。可以看出，4種圍巖波速對應的異常信號出現(xiàn)時刻分別為0.05、0.047、0.044、0.042 s，即從斷層處反射回來的時間在逐漸縮短。根據(jù)2節(jié)計算過程，可求出每種圍巖波速的斷層位置距離激發(fā)掌子面位置分別為100、98.7、96.8、96.6 m。加上已開挖洞段的50 m，可知預報斷層位置分別為150、148.7、146.8、146.6 m。預報斷層位置與模型中的斷層位置相比分別向掌子面方向移動了0、1.3、3.2、3.4 m，即預報斷層位置逐漸向掌子面靠近。此現(xiàn)象說明在進行隧道超前地質(zhì)預報中，如遇到賦存地應力隧道時，為了準確推測斷層破碎帶等不良地質(zhì)體位置，應對物探測試結果進行一定的修正，這對地震波法隧道超前地質(zhì)預報有一定的參考意義。

圖6 斷層傾角0°、不同圍巖波速時的地震干涉結果Fig.6 The results of seismic wave interference when θ = 0° and wave velocity of surrounding rock is different

3.2 工況2：斷層傾角為15°

對原始透射波信號進行濾波、反褶積及互相關干涉處理后，斷層傾角θ=15°時不同圍巖波速的地震干涉結果如圖7所示。可以看出，4種圍巖波速依次對應的異常信號出現(xiàn)時刻分別為0.049、0.046、0.043、0.041 s。根據(jù)斷層傾角為0°時的計算過程，可得出每種圍巖波速的斷層位置與激發(fā)掌子面的距離分別為98.0、96.6、94.6、94.3 m。加上已開挖洞段的50 m，可得到預報的斷層位置分別為148.0、146.6、144.6、144.3 m，這與模型中的斷層位置相比分別向掌子面移動了2.0、3.4、5.4、5.7 m。與斷層傾角為0°的情形相比(圖6)，傾角為15°的斷層向掌子面移動距離更多。

圖7 斷層傾角15°、不同圍巖波速時的地震干涉結果Fig.7 The results of seismic wave interference when θ=15° and wave velocity of surrounding rock is different

3.3 工況3和工況4：斷層傾角為30°和45°

根據(jù)斷層傾角0°時的計算方法，得出工況3和工況4預報斷層位置如表2所示。從表2可知，工況3和工況4向隧道掌子面移動的最大距離分別為14.9 m和31.8 m，這表明斷層傾角和圍巖波速對預報位置的影響非常顯著。因此，在地應力隧道進行超前地質(zhì)預報時圍巖波速是不可忽視的重要因素之一。

表2 工況3和工況4斷層預報位置Table 2 Predicting location of fault of case 3 and case 4

不考慮其他影響因素的前提下，對上述數(shù)值模擬結果中的修正值進行二元二次函數(shù)回歸分析，得到初步的斷層位置修正數(shù)值計算表達式為

y=-0.110 5+0.016 1α2+(-0.000 1)Δvα+0.007 7Δv+(-0.049)α

(4)

式(4)中：y為修正值；Δv為速度增加值；α為斷層傾角。

4 工程驗證

采用水平聲波剖面法(horizontal sonic profiling，HSP)對白馬隧道右洞YK39+041～YK38+961進行預報。測試掌子面巖性主要為灰黑色薄-中層狀炭質(zhì)板巖，如圖8所示，巖體較破碎，巖層產(chǎn)狀255∠75°，屬軟巖。主要發(fā)育一組節(jié)理，節(jié)理產(chǎn)狀160∠40°，延伸長度0.5～3.5 m，間距0.4～0.8 m，裂隙微張，波浪狀形態(tài)，面粗糙，無充填，掌子面干燥。

圖8 白馬隧道YK39+041掌子面照片F(xiàn)ig.8 The photo of Baima tunnel face (YK39+041)

掌子面前方圍巖XOY方向和YOZ方向的縱波速度分布如圖9所示。從物探測試圖推測前方80 m(YK39+041～YK38+961)范圍內(nèi)，圍巖基本以薄-中層板巖夾含炭質(zhì)板巖為主，巖質(zhì)軟，巖體較破碎，節(jié)理裂隙較發(fā)育。其中YK39+035～YK39+020段、YK39+010～YK39+000段及YK38+981～YK38+968段存在較強反射，且波速明顯降低。

圖9 YK39+041掌子面縱波速度測試成果Fig.9 Results of vertical wave velocity in YK39+041 tunnel face (YK39+041)

由于該段埋深約348 m，最大水平主應力高達14.92 MPa。不考慮傾角等其他因素的影響下，將波速增加值代入式(4)，得到斷層位置的修正值為7.5 m，故須對異常位置進行修正。最終確定YK39+027.5～YK39+012.5、YK39+002.5～YK38+992.5、YK38+973.5～YK38+960.5段為實際預報異常段落。推測在上述段落內(nèi)巖體破碎，巖性以炭質(zhì)板巖發(fā)育為主，節(jié)理裂隙密集發(fā)育，開挖后拱部易坍塌掉塊，局部炭質(zhì)富集，可能有有毒氣排出，圍巖含水；建議圍巖級別為Ⅴ級，施工時應加強超前支護。

隧道開挖驗證，當掌子面開挖至YK39+027時，掌子面發(fā)育一斷層，斷層走向與隧道軸線呈大角度相交。掌子面巖性主要為板巖，層間結合差，修正的圍巖基本質(zhì)量指標[BQ]=179.5～252.5，節(jié)理裂隙密集發(fā)育。通過鉆孔取芯，發(fā)現(xiàn)YK39+002.5～YK38+992.5段處于斷層影響帶范圍內(nèi)，巖性主要為薄-中層板巖，局部夾炭質(zhì)板巖，其中在掌子面YK38+995右上方存在股狀水，瓦斯?jié)舛鹊陀陬A警值。隧道YK38+974～YK38+960.5段巖性主要為薄-中層板巖，巖層產(chǎn)狀為250～265∠70°。節(jié)理裂隙發(fā)育，節(jié)理產(chǎn)狀160～175∠45°，延伸長度0.5～4.0 m，間距0.4～1.2 m，裂隙微張，平直形態(tài)，面粗糙，石英充填，掌子面濕潤，隧道開挖情況與修正后的預報結果基本一致。

5 結論

采用多源地震波干涉法和工程應用的方法分析了地應力對預報斷層位置的影響。得出如下主要結論。

(1)多源地震波干涉法數(shù)值模擬結果表明：地應力增大將使圍巖波速增大，隨著圍巖波速和斷層傾角的增大，斷層位置逐漸向掌子面方向移動，且斷層傾角增大使得斷層位置移動的距離較圍巖波速移動的距離更多。

(2)根據(jù)數(shù)值模擬結果中的修正值，通過二元二次函數(shù)回歸分析初步得到賦存地應力隧道中預報斷層位置的修正公式。

(3)采用水平聲波剖面法對白馬隧道高地應力段掌子面前方圍巖破碎段進行了預報，實際開挖顯示圍巖破碎洞段較物探預報位置向掌子面移動一定距離，進一步說明為提高隧道超前地質(zhì)預報精度，需要對物探結果進行適當修正。

所采用的多源地震干涉法主要是對掌子面前方斷層位置的預報，對于非斷層或小的節(jié)理裂隙密集帶等不良地質(zhì)體的判斷，該方法的適用性還需要進一步的研究。