基于多源地震干涉法隧道断层位置预测数值模拟

范占锋， 蔡建华， 赵 伟

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

随着中国的交通建设和能源开发不断发展，特别是西部山区在建或即将修建的高铁、公路以及水利工程都避绕不开无数大小断层破碎带、高地应力以及大变形区域[1-4]。在高地应力地区进行地下工程建设，如何提高隧道超前地质预报的准确性是当前工程界面临的一大关键难点。

从20世纪90年代至今，伴随着中国大型基础设施的建设，如大型水电站引水隧洞、公路隧道、铁路隧道等的开挖，隧道超前地质预报技术在中国隧道施工中也逐渐趋于完善，预报精度不断提高，预报方法趋于多样化。其中弹性波反射法应用最多，其次是电磁波反射法、瞬变电磁法、激发激化法和核磁共振法等。这些预报方法对不良地质体的识别侧重点不同。为了提高预报精度，很多专家学者提出了综合超前地质预报方法[5-8]，主要采用长短搭配组合方法。

然而，当隧道赋存地应力时，由于地应力的存在会使得弹性波在节理岩体中的传播有一定的影响，进而影响预报精度。赵明阶[9]推导了节理岩体在均匀应力场作用下的弹性波传播速度和衰减随应力的变化关系，并通过试验进行了验证。刘高等[10]研究了金川矿体的岩体力学参数及时空演化规律时指出地应力是影响岩体结构的主要因素，最终可通过弹性波波速来综合反映。田家勇等[11]基于岩石的声弹理论研究了不同初始应力条件下沿不同方向的弹性波波速，表明岩石波速和应力之间存在特定的关系。杜存苍[12]采用数值模拟软件FLAC3D分析了篮家岩隧道高地应力段围岩变形规律，但未对高地应力段的不良地质体如何预报进行解释说明。近年来，一些学者还尝试采用新的预报技术，吴丰收等[13]认为目前的隧道超前地质预报方法距离短、频次多，占用施工时间并存在安全隐患，基于互相关地震干涉法理论，提出了将多源地震干涉法(multisource seismic interferometry)应用于隧道超前地质预报。汪旭等[14]认为不论是隧道盾构法施工还是钻爆法施工，其实际震源很复杂并且施工环境会对信号采集质量有影响，从数值模拟角度探索了将多源地震干涉技术应用于隧道超前地质预报。多源地震干涉技术能够对接收到的透射波地震记录进行相互干涉运算，能够从透射波信号中提取隧道前方不良地质体的反射波，实现对前方地质条件做出判断的功能，这表明多源地震干涉技术在隧道超前预报领域具有广泛的应用前景[15-16]。

可见，隧道超前地质预报实施效果除了与目前本身技术水平有关外，还与岩体的应力状态关系极大。因此，主要基于多源地震干涉法，模拟不同围岩波速和不同断层倾角对断层位置识别的影响，通过数值模拟建立地应力和断层倾角与断层位置识别的修改关系式，进而对白马隧道高地应力段不良地质体的预报结果进行修正。

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

1.1 多源地震

多源地震的波动方程可表示为[6]

∇2u+Sbl(t,X)

(1)

式(1)中：u和v分别为位移和速度；∇2为Laplace算子；s(t)δ(X-X0)为震源函数；X0为震源位置；Sbl(t,X)为多震源波动方程的震源项；Γ(X)为Gamma函数，表示多震源混合矩阵算子，可表示为

(2)

式(2)中：γ1,γ2,…γn表示第1～n个震源地震记录；Tn(X)为位置X处震源n的延迟激发时间；j为虚数单位；ω为圆频率。

从式(1)可以看出，多震源激发环境下的地震记录是由单个震源地震记录的线性叠加所组成。

1.2 地震波干涉法

地震波干涉法是通过对记录到的地震信号进行干涉得到新的地震信号。新地震信号不仅包含了原始地震信号的特性，也包含了原始信号所不具有的某些重要特征[17]。目前，该技术还处于理论阶段，距离实际应用还有一定差距，原因在于实际工程中震源复杂，产生振动的频率和能量等难以识别[18]。同时，实际工程中隧道所处的地质条件复杂对成像也会造成一定的困难。采用数值模拟手段，对影响成像结果的主要因素进行分析研究，有利于该技术的实际应用。

Claerbout最早通过对水平层状介质模型研究，利用自相关运算将透射波记录和反射记录建立了有效联系，即透射波地震记录的自相关等价于自激自收模式记录，Wapenaar等[19]使用格林定理严格地证明了Claerbout的研究，为地震干涉法奠定了坚实的数学和物理基础。将透射波记录和反射波记录建立有效联系，形成对隧道掌子面前方异常体的目标成像，互相关干涉法原理可表示为

(3)

式(3)中：T(xA,xi, -t)和T(xB,xi,t)分别为从地下震源xi传到A检波点和B检波点的透射地震记录；xA和xB分别为A检波点和B检波点的位置；t为时间；δ为Dirac函数；R(xA,xB,t)为以B为震源，A为检波点的反射地震记录；*表示卷积；相关型地震波干涉法提取的地震信号既包含因果部分R(xA,xB,t)，也包含非因果部分和零时刻的脉冲响应R(xA,xB, -t)。

2 数值计算模型及验证

2.1 计算模型

在深埋隧道建设中，地应力使断层破碎带或节理裂隙密集发育带岩体被压密、节理闭合，地震波在此类介质中传播速度增加[20-21]。据此，建立尺寸为长×宽=400 m×400 m的二维平面模型，如图1所示，隧道剖面宽10 m，信号激发端和接收端隧道洞段分别开挖50 m。在隧道中轴线上，沿水平方向150 m处设置1条断层，断层距离信号接收端掌子面为100 m，贯穿整个模型，断层与掌子面平行，如图1(a)所示。同时设置几种倾斜断层，如图1(b)所示，具体计算参数如表1所示。

图1 发育1条断层的二维平面数值计算模型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。利用有限差分交错网格法，通过MATLAB编程，分析围岩波速和断层倾角两因素对断层位置识别的影响。

表1 数值计算模型参数Table 1 Parameters of the numerical simulation model

2.2 算法验证

首先验证所用数值模拟算法是否正确。假设激发震源函数采用典型的雷克(Ricker)子波，主频为35 Hz，震源激发掌子面位于模型左侧，共布设10个震源，震源间距均为1 m；透射波信号接收点位于右侧，共布设21个检波器，检波器间距均为0.5 m。断层倾角θ=0°和围岩波速Cs=4 000 m/s时的数值模型速度分布如图2所示。由图2可知，围岩波速分布均匀，在此基础上进行正演计算。

图2 断层倾角0°和围岩波速4 000 m/s时的数值模型 速度分布Fig.2 Numerical model velocity distribution when θ=0° and Cs = 4 000 m/s

断层倾角θ=0°时正演计算的不同时间波场传播分布如图3所示。由图3可知，地震波传播时间从5 ms增加到16 ms时，由于围岩为各向同性均质体，观察到半椭圆形的带状区域逐渐增大，没有发生波的反射，即透射波未到达围岩与断层交汇处。当传播时间增加到30 ms时，在设置的断层位置处(150 m)发生了明显的反射现象，透射波与反射波之间形成了“透镜”型交汇区域。

图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所示。根据杨为民等[6]的结论，该成像结果等价于自激自收剖面。因此，成像结果中的时间皆为双程旅行时间，即在0.05 s附近区域存在异常界面。根据围岩波速4 000 m/s，可计算出单程旅行时异常界面(断层)到激发掌子面的距离为100 m，加上已开挖洞段50 m，实际预报断层位置位于150 m，与模型中设置的断层位置完全吻合，验证了数值模拟算法的正确性。

图4 断层倾角0°和围岩波速4 000 m/s时的地震干涉结果Fig.4 The result of the seismic wave interference when Cs = 4 000 m/s and θ = 0°

3 不同数值模拟工况计算

在验证所用算法正确性的基础上，根据表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所示。可观察到在断层位置发生了明显的透、反射现象，透射波传播形态基本相似，在反射波和透射波之间都形成了一个“透镜”型区域，反射波与透射波传播方向在一条直线上。所不同的是，从图5可以看出，随着围岩波速从4 000 m/s增加到4 600 m/s，“透镜”型区域面积在不断增大，其原因在于围岩波速增加使得从断层位置处反射回来的波速传播加快所造成。

通过对原始透射波信号进行滤波、反褶积及互相关干涉处理后，断层倾角θ=0°时不同围岩波速的地震干涉结果如图6所示。可以看出，4种围岩波速对应的异常信号出现时刻分别为0.05、0.047、0.044、0.042 s，即从断层处反射回来的时间在逐渐缩短。根据2节计算过程，可求出每种围岩波速的断层位置距离激发掌子面位置分别为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，即预报断层位置逐渐向掌子面靠近。此现象说明在进行隧道超前地质预报中，如遇到赋存地应力隧道时，为了准确推测断层破碎带等不良地质体位置，应对物探测试结果进行一定的修正，这对地震波法隧道超前地质预报有一定的参考意义。

图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种围岩波速依次对应的异常信号出现时刻分别为0.049、0.046、0.043、0.041 s。根据断层倾角为0°时的计算过程，可得出每种围岩波速的断层位置与激发掌子面的距离分别为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°

根据断层倾角0°时的计算方法，得出工况3和工况4预报断层位置如表2所示。从表2可知，工况3和工况4向隧道掌子面移动的最大距离分别为14.9 m和31.8 m，这表明断层倾角和围岩波速对预报位置的影响非常显著。因此，在地应力隧道进行超前地质预报时围岩波速是不可忽视的重要因素之一。

表2 工况3和工况4断层预报位置Table 2 Predicting location of fault of case 3 and case 4

不考虑其他影响因素的前提下，对上述数值模拟结果中的修正值进行二元二次函数回归分析，得到初步的断层位置修正数值计算表达式为

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进行预报。测试掌子面岩性主要为灰黑色薄-中层状炭质板岩，如图8所示，岩体较破碎，岩层产状255∠75°，属软岩。主要发育一组节理，节理产状160∠40°，延伸长度0.5～3.5 m，间距0.4～0.8 m，裂隙微张，波浪状形态，面粗糙，无充填，掌子面干燥。

图8 白马隧道YK39+041掌子面照片Fig.8 The photo of Baima tunnel face (YK39+041)

掌子面前方围岩XOY方向和YOZ方向的纵波速度分布如图9所示。从物探测试图推测前方80 m(YK39+041～YK38+961)范围内，围岩基本以薄-中层板岩夹含炭质板岩为主，岩质软，岩体较破碎，节理裂隙较发育。其中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段为实际预报异常段落。推测在上述段落内岩体破碎，岩性以炭质板岩发育为主，节理裂隙密集发育，开挖后拱部易坍塌掉块，局部炭质富集，可能有有毒气排出，围岩含水；建议围岩级别为Ⅴ级，施工时应加强超前支护。

隧道开挖验证，当掌子面开挖至YK39+027时，掌子面发育一断层，断层走向与隧道轴线呈大角度相交。掌子面岩性主要为板岩，层间结合差，修正的围岩基本质量指标[BQ]=179.5～252.5，节理裂隙密集发育。通过钻孔取芯，发现YK39+002.5～YK38+992.5段处于断层影响带范围内，岩性主要为薄-中层板岩，局部夹炭质板岩，其中在掌子面YK38+995右上方存在股状水，瓦斯浓度低于预警值。隧道YK38+974～YK38+960.5段岩性主要为薄-中层板岩，岩层产状为250～265∠70°。节理裂隙发育，节理产状160～175∠45°，延伸长度0.5～4.0 m，间距0.4～1.2 m，裂隙微张，平直形态，面粗糙，石英充填，掌子面湿润，隧道开挖情况与修正后的预报结果基本一致。

5 结论

采用多源地震波干涉法和工程应用的方法分析了地应力对预报断层位置的影响。得出如下主要结论。

(1)多源地震波干涉法数值模拟结果表明：地应力增大将使围岩波速增大，随着围岩波速和断层倾角的增大，断层位置逐渐向掌子面方向移动，且断层倾角增大使得断层位置移动的距离较围岩波速移动的距离更多。

(2)根据数值模拟结果中的修正值，通过二元二次函数回归分析初步得到赋存地应力隧道中预报断层位置的修正公式。

(3)采用水平声波剖面法对白马隧道高地应力段掌子面前方围岩破碎段进行了预报，实际开挖显示围岩破碎洞段较物探预报位置向掌子面移动一定距离，进一步说明为提高隧道超前地质预报精度，需要对物探结果进行适当修正。

所采用的多源地震干涉法主要是对掌子面前方断层位置的预报，对于非断层或小的节理裂隙密集带等不良地质体的判断，该方法的适用性还需要进一步的研究。