基于黏弹性介质的隧道地震二维正演模拟

凌 飞，肖宏跃，朱夏乐，雷 宛，蒋亚东

(成都理工大学 a.地球物理学院; b.环境与土木工程学院，成都 610059)

基于黏弹性介质的隧道地震二维正演模拟

凌 飞a，肖宏跃a，朱夏乐a，雷 宛b，蒋亚东a

(成都理工大学 a.地球物理学院; b.环境与土木工程学院，成都 610059)

针对隧道中各种不良地质体的超前预报问题，以黏弹性介质为基础，在空气隧道周围添加一定厚度的开挖损伤带，建立符合实际情况的观测系统，通过交错网格有限差分正演，模拟地震波在断层、软弱夹层、溶洞等地质灾害体模型中的传播，根据不同时刻的波场快照来研究相应地质模型地震记录的响应特征。通过对隧道地震二维波场模拟结果进行分析，研究结果表明：对一阶速度-应力黏弹性波动方程，采用交错网格有限差分，并且对隧道添加约束条件，可以快速、有效地实现隧道地震的二维正演模拟，并且波场快照能够清楚地反映地震波在不同地质模型中的传播特征以及波的转换特征，地震记录能够反映不同地质模型的地震响应特征，可以为实际隧道超前预报提供理论依据。

黏弹性介质；隧道地震；交错网格；开挖损伤带；波场模拟

1 研究背景

地下隧道施工过程中，由于地质条件的变化，常易发生突发性的意外情况。据国内隧道施工的不完全统计，施工过程中因塌方、涌水、突泥、岩爆、高瓦斯等地质灾害事故造成的停工时间大约占总工期的30%。通过隧道地质超前预报，及时预报隧道掌子面前方的不良地质体及其性质、位置、规模、产状与成灾可能性，预防涌水、突泥、坍塌等灾害性事故的发生，为正确选择开挖断面、支护设计参数和优化施工方案提供依据。以工程地质分析法为基础，以地球物理探测方法为手段的隧道综合地质超前预报方法已成为隧道工程界的共识，并得到广泛的应用。在隧道超前预报中，常用的地震波法有：负视速度法、水平剖面法(HSP)、隧道地震反射波超前预测法(TSP)、真反射层析成像法(TRT)、隧道地震波层析成像法(TST)、隧道地质超前预报(TGP)等[1]。由于隧道超前地质预报的特点决定了现场工作必须在狭小的隧道施工现场实时完成，因此，超前探测技术必须满足狭小场地条件下进行长距离、高精度、高效率探测要求。大量的工程实践表明，隧道开挖过程中最常见的不良地质体是断层破碎带、岩溶、软弱夹层、采空区等。有不少隧道处在复杂岩性、构造及高地应力地区，尤其是在复杂的岩溶地区，因受岩性、地层组合、构造、地貌等诸多因数的影响控制，岩溶发育复杂多变、大小不一、形态各异，给勘察工作带来更大难度[2]。研究和了解这些地质体的地质地球物理特征将对于准确的探测与预报十分关键。

近年来，许多学者在隧道地震的数值模拟方面做了相关的研究[3-7]。杜立志[3]对不同倾角的地层进行了波场模拟，但并未将隧道作为空气介质考虑；鲁光银等[5]将隧道空间作为自由边界，刘江平等[4]对隧道添加吸收边界条件，两者均对不同的地质模型进行了正演模拟，但均未进行黏弹性介质的数值模拟。

为了模拟结果更接近实际地质情况，本文以开尔芬黏弹性介质为基础，从一阶速度-应力黏弹性波动方程出发，采用交错网格有限差分[8]对含开挖扰动带的复杂隧道介质模型进行二维正演模拟。对模拟结果进行对比分析，为实际隧道超前预报提供理论依据及指导。

2 黏弹性介质波动方程及其有限差分格式

在非均匀各向同性黏弹性介质中，假设不考虑外力项，二维开尔芬(Kelvin)黏弹性波动方程的一阶速度-应力方程形式如下[9]：

(1)

式中：B是密度ρ=ρ(x,z)的倒数，或者叫做疏度；vx(x,z,t)和vz(x,z,t)分别是x和z方向的速度分量；σxx(x,z,t)，τxz(x,z,t)和σzz(x,z,t)是应力分量；λ和μ为均匀弹性各向同性介质的拉梅系数;λ′和μ′为黏性介质的拉梅系数。式(1)中拉梅系数为：

(2)

式中：vp和vs分别为介质的纵波和横波速度；ω为圆频率；Qp，Qs分别为纵波和横波的品质因子。

(3)

(4)

(5)

(6)

(7)

(8)

(9)

3 隧道约束条件及边界条件

由于超前预报是基于空气隧道介质模型来进行，所以正演模拟时，需要对隧道添加约束条件。首先是对隧道中空气密度的处理。当数值计算至空气隧道壁时，因为剪切模量突然下降至0，并且压力会下降一个数量级，这将导致数值模拟发生错误。为了避免这种错误的发生，需要将隧道中空气的密度设置成一个较大值，本文将其密度设置成150 kg/m3。

其次是在隧道周围添加开挖损伤带(EDZ)。由于隧道在开挖过程中，会对隧道壁周围的岩体造成损伤。通常情况下，越靠近隧道壁的围岩，其损伤程度越大，而检波器则需要安置在未受损伤的围岩中。围岩的损伤程度取决于许多因素，如开挖方法、岩体力学特性、初始地应力场、裂隙分布等[10]。本文将隧道损伤带的宽度设置为0.5 m。

最后，对隧道里的数值模拟采用时间二阶、空间四阶的弹性波有限差分公式[5]进行计算，并且进行吸收边界处理[4]。

对于模型边界条件的处理，本文采用了PML完全匹配层吸收边界[11]，可以较好地吸收传播到模型边界的地震波。

4 观测系统及地质模型的建立

为了使模拟结果更接近实际情况，本文建立了与实际TSP隧道超前预报相应的观测系统以及各种地质模型。正演模拟采用的观测系统如图1所示，隧道长50 m，宽8 m，震源位于隧道壁一侧，沿隧道壁一共安置了24个检波点，最小偏移距为10 m，道间距为1.5 m，震源埋深2 m，检波点埋深1.5 m。

图1 观测系统示意图

本文建立了3类地质模型，即断层、软弱夹层以及溶洞模型，如图2所示。通过建立不同类型的模型来了解并研究地震波在不同地质体中的传播规律，以及地震记录的响应特征；建立不同倾角(模型分界面与隧道中心轴线的夹角)的断层和软弱夹层模型来研究随界面倾角变化波场以及地震记录的变化规律。本文建立了界面倾角为90°和60°的2种地质模型，模型的分界面距离掌子面80 m，软弱夹层的宽度为10 m。通过建立不同物理特性参数的溶洞模型来研究不同充填物质对地震记录的影响。溶洞位于隧道掌子面正前方，溶洞中心距离掌子面83 m，其半径为3 m。模型具体参数见表1。

图2 地质模型

模型材料Vp/(m·s-1)Vs/(m·s-1)ρ/(kg·m-3)QpQs隧道3400．1150围岩30001732．020005040EDZ22001270．016004030断层25001443．018004535软弱夹层20001154．014003022溶洞11500866．010003020溶洞23400．1150105

图3 直立断层模型波场快照

5 波场数值模拟

通过添加隧道约束条件，运用黏弹性介质中的差分格式，针对建立的各种不良地质体实现二维波场模拟。本文模型大小为250 m×100 m，网格间隔Δx=Δz=0.5 m，采样间隔Δt为0.05 ms，地震记录长度为175 ms，震源采用200 Hz的Ricker子波。

5.1 断层模型模拟结果及分析

对直立断层模型不同时刻的波场快照(图3所示)分析，当t=35 ms时，可以观察到由震源激发产生的直达纵波(P)、横波(S)、沿隧道壁传播的瑞利面波(R)和绕射波。瑞利面波在此刻到达隧道掌子面，其部分能量转换成横波(RS)继续向前传播，也有部分能量的瑞利面波会沿隧道壁往后传播。t=55 ms时，直达纵波(P)在断层分界面处产生了反射波和透射波，此时可以清楚地看到反射纵波(PP)和透射纵波(TP)。t=75 ms时，由于纵横波存在波速差异，导致纵横波之间的距离逐渐增加，此时可以明显看到反射转换横波(PS)和透射转换横波(TPS)。t=95 ms时，横波(S)已到达断层分界面并产生了反射横波(SS)。

在分析地震记录响应特征的过程中，先确定不同波的同相轴，再比较相同道数记录的峰值大小，从而判断不同波之间的能量强弱关系。图4和图5分别是在弹性和黏弹性介质中对直立断层模型做正演模拟得到的地震记录。与弹性介质相比，黏弹性介质中的反射波能量明显减弱，原因是黏弹性介质考虑了地层对地震波能量的吸收作用，使模拟更接近实际地质情况。

图4 直立断层模型地震记录(弹性介质)

图5 直立断层模型地震记录(黏弹性介质)

图6是倾角为60°断层模型的地震记录，观察黏弹性介质中2种倾角断层模型的地震记录发现，t≤60 ms时，地震记录主要由以下几种波构成：直达纵波(P)、由隧道壁以及开挖损伤带产生的反射波以及面波(R)，这几种波的能量都很强。t>60 ms时，地震记录主要反映来自断层分界面的反射波。对比图5中直立断层模型的x和z分量地震记录，发现x分量记录中反射纵波(PP)能量强于z分量，而z分量记录中反射转换横波(PS)和反射横波(SS)的能量均强于x分量。对比60°倾角断层模型的x和z分量地震记录，发现z分量记录中反射转换横波(PS)和反射横波(SS)的能量依然强于x分量。分别对比2种断层模型的x分量和z分量记录，较直立断层模型，60°倾角断层模型的x分量反射纵波(PP)能量增强，但是反射横波(SS)能量减弱；其z分量记录反射纵波(PP)和反射转换横波(PS)能量均增强，但是反射横波(SS)能量减弱。

图6 倾角60°断层模型地震记录

5.2 软弱夹层模型模拟结果及分析

不同倾角的软弱夹层模型地震记录如图7和图8所示，x分量和z分量记录均可以看出地震波到达软弱夹层2个分界面产生的反射波，x分量记录的反射纵波能量要强于z分量记录，而z分量记录的反射横波能量要强于x分量记录。此外，2个分量记录中，2条反射纵波同相轴的波形反向(由PP1和PP2表示)，这也验证了在坚硬围岩中间存在着一软弱夹层。对比2个软弱夹层模型的x分量和z分量记录，对于倾角为60°的软弱夹层模型而言，其x分量和z分量记录中，反射纵波和反射转换纵波能量增强，但是反射横波及反射转换横波能量减弱。

图7 直立软弱夹层模型地震记录

图8 倾角60°软弱夹层模型地震记录

图9 溶洞模型1(水)地震记录

图10 溶洞模型2(空气)地震记录

5.3 溶洞模型模拟结果及分析

溶洞模型1全部充填水，溶洞模型2全部充填空气，空气溶洞的物理特性参数与隧道一样。图9和图10为2种不同填充物质溶洞的x和z分量地震记录。对于2种溶洞模型的x分量记录而言，反射纵波(PP)能量明显强于反射横波(SS)能量，富水溶洞的反射纵波(PP)和反射横波(SS)能量均弱于空气溶洞，并且富水溶洞其反射波的叠加现象比空气溶洞明显。对于z分量记录而言，反射横波(SS)能量突出，反射纵波(PP)能量减弱。富水溶洞的反射纵波(PP)和反射横波(SS)能量依然弱于空气溶洞，但观察初至的反射横波(SS)之后的反射波，富水溶洞的反射波能量强于空气溶洞。

6 结 论

(1) 对一阶速度-应力黏弹性波动方程，采用交错网格有限差分，模型外部边界采用PML完全匹配层吸收边界条件，对隧道添加约束条件，能够快速、准确以及有效地实现隧道地震的二维正演模拟；波场快照能够清楚地反映地震波在不同地质模型中的传播特征以及波的转换特征，地震记录能够反映不同地质模型的地震响应特征，可以为实际隧道超前预报提供理论依据。

(2) 由于在黏弹性介质中地层对弹性波的吸收衰减作用，使地震波能量减弱，更符合实际地质的情况。通过对隧道地震模拟结果进行分析，可以更好地了解地震波在地层中的传播和衰减规律。

(3) 对本文建立的坐标系而言，x分量记录中，反射纵波能量突出，有利于对P波的提取与研究；z分量记录中，反射横波能量突出，有利于对S波的提取与研究。

(4) 对于不同倾角的地质模型，当地质界面倾角减小时，反射纵波能量增强，反射横波能量减弱；对于软弱夹层，夹层第1个分界面的反射波表现为正视速度特征，而第2个界面的反射波表现为负视速度特征；对于溶洞，富水溶洞反射波能量弱于空气溶洞，其反射横波叠加现象比空气溶洞明显。

[1] 赵永贵. 国内外隧道超前预报技术评析与推介[J]. 地球物理学进展，2007，22(4)：1344-1352.(ZHAO Yong-gui. Analysis and Recommendation of Tunnel Prediction Techniques at Home and Abroad[J]. Progress in Geophysics，2007，22(4)：1344-1352. (in Chinese))

[2] 宋先海,顾汉明,肖柏勋. 我国隧道地质超前预报技术述评[J]. 地球物理学进展，2006，21(2)：605-613. (SONG Xian-hai, GU Han-ming, XIAO Bo-xun. Overview of Tunnel Geological Advanced Prediction in China[J]. Progress in Geophysics，2006，21(2)：605-613.(in Chinese))

[3] 杜立志. 隧道施工地质地震波法超前探测技术研究[D]. 长春：吉林大学，2008. (DU Li-zhi. Study of Tunnel Seismic Geological Prediction Technology in the Process of Construction[D]. Changchun：Jilin University，2008. (in Chinese))

[4] 刘江平，程 飞，范承余，等. 基于隧道空间全波场二维数值模拟与特征分析[J]. 岩土工程学报，2012，34(9)：1705-1711.(LIU Jiang-ping, CHENG Fei, FAN Cheng-yu,etal. Two-dimensional Numerical Simulation of Tunnel-based Seismic Full-wave Fields[J]. Chinese Journal of Geotechnical Engineering，2012，34(9)：1705-1711.(in Chinese))

[5] 鲁光银，熊 瑛，朱自强. 隧道反射波超前探测有限差分正演模拟与偏移处理[J]. 中南大学学报(自然科学版)，2011，42(1)：136-141. (LU Guang-yin, XIONG Yin, ZHU Zi-qiang. Detection Simulation ahead of Tunnel Face and Reverse-time Migration with Reflection Wave Method[J]. Journal of Central South University (Science and Technology)，2011，42(1)：136-141. (in Chinese))

[6] KNEIB G, LEYKAM A. Finite-difference Modelling for Tunnel Seismology[J]. Near Surface Geophysics，2004，2(2)：71-93.

[7] BOHLEN T, LORANG U, RABBEL W,etal. Rayleigh-to-shear Wave Conversion at the Tunnel Face: From 3D-FD Modeling to Ahead-of-drill Exploration[J]. Geophysics，2007，72(6):67-79.

[8] ROBERTSSON J O A, BLANCH J O, SYMES W W. Viscoelastic Finite-difference Modeling[J]. Geophysics, 1994, 59(9):1444-1456.

[9] 奚 先,姚 姚. 二维粘弹性随机介质中的波场特征分析[J]. 地球物理学进展，2004，19(3)：608-615. (XI Xian，YAO Yao. The Analysis of the Wave Field Characteristics in 2-D Viscoelastic Random Medium[J]. Progress in Geophysics，2004，19(3)：608-615. (in Chinese))

[10]吉小明. 隧道开挖的围岩损伤扰动带分析[J]. 岩石力学与工程学报，2005，24(10)：1697-1702. (JI Xiao-ming. Study on Mechanical and Hydraulic Behavior of Tunnel Surrounding Rock Masses in Excavation-Disturbed Zone[J]. Chinese Journal of Rock Mechanics and Engineering，2005，24(10)：1697-1702. (in Chinese))

[11]COLLINO F, TSOGKA C. Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media[J]. Geophysics，2001，66(1):294-307.

(编辑：黄 玲)

Two-dimensional Forward Seismic Modeling for TunnelBased on Viscoelastic Medium

LING Fei1，XIAO Hong-yue1，ZHU Xia-le1，LEI Wan2，JIANG Ya-dong1

(1.College of Geophysics，Chengdu University of Technology，Chengdu 610059，China； 2.College of Environment and Civil Engineering, Chengdu University of Technology, Chengdu 610059,China)

In view of the problems in advanced detection of defective geological structures in tunnel, a geological tunnel model was established based on viscoelastic medium to simulate the seismic wave propagation in fault, weak interlayer, karst cave and other geological hazard body. The tunnel is surrounded by the excavation damaged zone. A practical observation system and the staggered grid finite-difference method were employed for the simulation. According to the wave field snapshots at different instants, the response of seismic record acquired from the corresponding forward modelling was researched. Result shows that the two-dimensional forward simulation based on viscoelastic medium can be achieved by using a staggered grid finite-difference algorithm and adding tunnel constraints. The wave field snapshots could clearly reflect the seismic wave propagation and transformation in different geological models. The seismic record reflects the seismic response characteristics of different models, and offers a theoretical basis for the seismic prediction of actual tunnels.

viscoelastic medium; tunnel seismic; staggered grid; excavation damaged zone; wave field simulation

2014-07-07；

2014-08-18

凌 飞(1989-)，男，四川成都人，硕士，主要研究方向为工程与环境地球物理勘探，(电话) 028-84073510(电子信箱)825729075@qq.com。

10.3969/j.issn.1001-5485.2015.05.023

2015,32(05):121-126

P631.4

A

1001-5485(2015)05-0121-06