基于柱坐标系的隧道空间全波场数值模拟与分析

鲁光银，罗帅，朱自强，石克亮，夏成志

基于柱坐标系的隧道空间全波场数值模拟与分析

鲁光银，罗帅，朱自强，石克亮，夏成志

(中南大学 地球科学与信息物理学院，湖南 长沙 410083)

不同于地面地震，隧道空间的特殊结构使得其波场十分复杂，基于实际隧道模型的全波场数值模拟是有效认识隧道空间内波场传播规律的重要手段。关于这方面的研究以往多是基于直角坐标系，而隧道作为类似空气柱的特殊地质模型，采用柱坐标系下的全波场数值模拟对认识隧道空间的波场特征更有意义。基于柱坐标系下的一阶弹性波速度-应力方程，推导任意偶数阶交错网格差分格式，设计了人工截断边界和隧道自由边界，对岩性分界面、断层破碎带、溶洞等隧道典型不良地质模型进行数值模拟与波场特征分析。研究结果表明：基于柱坐标系下的一阶弹性波速度-应力方程，采用交错网格有限差分算法结合FCT技术，可实现隧道空间全波场高精度数值模拟，波场特征符合波的运动学与动力学特征，为隧道空间复杂波场特征的认识和资料解译提供有效的依据。

超前预报；隧道空间；柱坐标系；全波场；数值模拟

随着交通运输的快速发展，越来越多的隧道被建造。然而，隧道作为一项隐蔽性地下工程，在勘察设计阶段很难精确的发现不良地质体。在隧道开挖过程中，有可能遇到断层、溶洞，暗河等不良地质现象，极易诱发突水、涌泥、岩爆等地质灾害，进而造成工期延误，经济损失，甚至人员伤亡等重大安全事故[1−2]。地震法作为隧道超前预报的主要方法之一，包括TSP，TST，HSP，TRT和陆地声纳法等[3−8]，以其长距离、相对高精度而得以广泛应用。但是隧道空间内的波场是十分复杂的，接收记录上往往存在直达波、反射波、面波、纵波、横波转换波[9]，这无疑为数据处理与成果解译增加了难度。基于隧道空间的全波场数值模拟是认识隧道空间内复杂波场传播规律的基础，同时也是观测系统设计，数据采集、处理、解释的重要依据。近年来有不少学者进行了相关研究，其中绝大部分是基于直角坐标系，鲁光银等[10]采用直接法处理隧道自由边界条件，并进行了地质界面模型的正演模拟与偏移成像；刘江平，程飞等[9, 11]采用吸收边界条件处理隧道自由边界，并进行了不同观测系统和不同倾角地质界面对波场特征和偏移成像的影响研究；凌飞等[12]基于黏弹性介质并在隧道周围引入一开挖损伤带对隧道典型不良地质现象进行了全波场数值模拟；宋杰[13]针对隧道工程建设中可能遇到的典型不良地质体，建立了岩性分界面、断层以及岩溶体等多个三维地质模型，进行了三维地震波探测正演模拟，分析总结了隧道典型不良地质体三维波场响应和地震成像特征；王京[14]基于二阶弹性波位移波动方程采用有限元方法对隧道典型不良地质体进行了数值模拟，并重点探讨了有限元方法实现自由边界条件的处理方法；Jetschny等[15]考虑了隧道空间的影响，利用面波和转换横波对隧道前方的地质构造进行预测，并且讨论了这2种波的转换机理；YANG等[16]基于三维弹性波方程交错网格有限差分算法，实现了煤矿巷道模型的三维槽波勘探数值模拟；Harmankaya等[17]实现了利用隧道掘进过程中产生的噪声作为震源进行超前探测的数值模拟，并验证了该方法对散射体的定位能力；LIU等[18]进行了极坐标系下的全波场数值模拟，并处理了复杂的边界条件；查欣杰等[19]通过构建二维含低速异常的隧道介质模型，研究了隧道弹性波场传播规律和异常体边界成像准确性。模拟隧道空间三维波场最经济的方法是把构造模型包括震源近似为轴对称模型。而相较于直角坐标系采用柱坐标系进行数值模拟具有以下2点优势[20]：1) 直角坐标系中的有限差分算法的公式比较简单, 但是处理轴对称地质构造需要比较多的计算成本；2) 对于具有轴对称性的波场问题, 可以转化为二维问题, 大大简化计算。论文基于柱坐标系下的一阶弹性波速度−应力方程，将隧道与围岩边界作自由边界处理，同时引入柱坐标系下的分裂式完全匹配层，采用时间二阶，空间十阶的高阶交错网格有限差分法并结合FCT技术，实现了柱坐标系下隧道空间的高精度全波场数值模拟，通过计算岩性分界面、断层破碎带、溶洞等复杂地质模型，分析了隧道典型不良地质现象的波场特征，为隧道空间复杂波场特征的认识和资料解译提供了依据。

1 波动方程及其有限差分离散

考虑到隧道空间的近似轴对称结构以及三维数值模拟的计算成本，本文将隧道空间作轴对称处理。对于具有轴对称的完全弹性介质，介质参数仅是和的函数，三维模型可以简化为二维。此时柱坐标系下的一阶弹性波速度−应力方程为：

式中：为介质密度；和为拉梅常数；v和v分别为轴向(方向)和径向(方向)的速度；和分别为轴向和径向的正应力；为切应力。

对上述方程采用交错网格[21]进行离散，得到任意偶数阶差分格式：

2 稳定性条件与边界条件

正演模拟过程中为保证离散数值计算的稳定性需要设置稳定性条件，与直角坐标系不同, 柱坐标系要求时间步长更小,稳定性条件为[22]：

式中：min为最小波长；min(Δ,Δ)为空间最小网格步长；max为最大波速。

若考虑为深埋隧道，轴对称柱坐标系下隧道空间全波场数值模拟需要处理2种边界条件，即计算区域外的人工截断边界和隧道自由边界条件(图1)。对于人工截断边界采用LIU[23]提出的柱坐标系下的分裂式完全匹配层处理；而隧道自由边界条件则采用横向各项同性介质替换法处理[24]，即用横向各向同性介质近似代替自由边界，直接令自由边界处的正应力为0，切应力为0的条件则通过对自由边界上物性参数的设定，在波动方程交错网格有限差分迭代求解中实现。二维情况下的处理方式为：

式中：τxx和τrr分别为隧道掌子面和侧壁的正应力；ρ0，λ0和μ0为自由边界上的密度和拉梅常数。ρ，λ和μ为弹性介质的密度和拉梅常数。

3 观测系统与地质模型

为了获取较好的波场记录，采用如图2所示的观测系统。隧道长70 m，宽10 m，与实际隧道相仿。炮点置于隧道掌子面与右边墙的交点处，埋深1 m；在隧道左右边墙各布置了47检波器，偏移距15 m，道间距0.5 m，埋深1.5 m。

图2 观测系统示意图

为验证本文算法的正确性，以及分析隧道典型不良地质现象的波场特征设计了如图3所示的3个地质模型。模型Ⅰ为垂直岩性分界面，界面距掌子面40 m；模型Ⅱ为断层破碎带，破碎带宽10 m，倾角45°；模型Ⅲ为溶洞，溶洞半径5 m，各模型的弹性参数见表1。

(a) ModelⅠ；(b) ModelⅡ；(c) Model Ⅲ

表1 模型弹性参数表

4 数值模拟结果及分析

针对上述理论模型采用高阶差分格式进行数值模拟，其中模型大小为，200 m×200 m，Δ=Δ=0.25 m，Δ=0.025 ms，震源选择中心频率为100 Hz的雷克子波。

4.1 岩性分界面模型

如图4为垂直界面模型径向分量20，35和50 ms的波场快照。由图4可知，在=20 ms时，由于没有遇到岩性分界面，在波场快照中可以看到直达纵波(P波)、直达横波(S波)、隧道自由表面上产生的面波(R波)；同时可以看到，面波延隧道侧壁传播并在掌子面与侧壁交点处转换为绕射横波(RS波)[25]，绕射横波传至侧壁时又重新转换成面波；在=35 ms时，在岩性分界面处处产生的反射纵波(PP波)和透射纵波(TP波)，以及由P波产生的反射转换横波(PS波)以及透射转换横波(TPS波)；在=50 ms时SS波出现，但由横波产生的转换纵波(SP波)能量较弱，不太明显。从数值模拟结果可以看出，波场传播特征符合波的运动学和动力学规律，证明了本文算法的正确性。

图5为垂直界面模型模拟记录。由图5可知，径向分量中的SS波、PS波能量强于轴向分量，而轴向分量中的PP波能量则强于径向分量；轴向分量中炮检同侧接收的PP波能量弱于炮检异侧，而径向分量中炮检同侧接收的PS波、SS波弱炮检异侧；界面一次反射纵波到达掌子面时发生反射，该反射波向前传至界面处再次反射形成多次反射波(MV波)，当其到达检波器的时间和界面一次反射波相近时就会与界面一次反射波相互叠加，若该多次反射波能量较强，就会严重影响有效波的识别，这在轴向分量的记录中体现的十分明显；同时，由单炮记录可估算出P波和R波的视速度分别约为1 879 m/s和1 000 m/s，而PP波和SS波的旅行时分别约为48.5 ms和83.5 ms这与模型设计相符，进一步证明了数值模拟的正确性。

(a) t=20 ms；(b) t=35 ms；(c) t=50 ms

(a) X分量炮检同侧接收；(b) X分量炮检异侧接收；(c) R分量炮检同侧接收；(d) R分量炮检异侧接收

4.2 断层破碎带模型

图6为断层破碎带模型模拟记录。如图6所示，与岩性分界面模型相比断层破碎带模型仍然具有轴向分量PP波能量较强，径向分量PS波、SS波能量较强的特征；炮检同侧接收的记录中顶、底界面的反射纵波(PP1，PP2)能量较强且同相轴较为清晰，可较好的区分破碎带顶、底界面，但在一定程度上受波场叠加干扰，而炮检波异侧接收记录波场叠加干扰严重，影响了有效波的识别；同时顶、底界面的反射横波(SS1，SS2)受波场叠加干扰较小，利于顶、底界面的的识别与横波成像；不同于单反射界面模型，断层破碎带模型在其顶底界面均会产生反射及透射，并且在顶底界面之间还会发生多次反射，因此模拟记录上波场比较复杂。

4.3 溶洞模型

图7为溶洞模型模拟记录。如图7所示，对于上述溶洞作模型，溶洞顶部的反射纵波(PP1)在轴向分量及径向分量的炮检同侧的接收记录中有比较清晰的体现，但溶洞底部的反射纵波(PP2)由于受波场叠加影响，对溶洞底部的识别会有一定程度的干扰，径向分量炮检异侧接收的反射纵波能量十分微弱，在记录中几乎无法体现；溶洞顶、底的反射横波(SS1，SS2)在径向分量的记录中有较好的体现，利于横波成像；由于溶洞内部的低速填充物对地震波能量具有聚集作用[26]，所以溶洞底部的反射波PP2能量较溶洞顶部反射波PP1能量更强；不同于界面反射，溶洞作为局部地质体，产生的绕射波会严重干扰有效波的识别。

(a) X分量炮检同侧接收；(b) X分量炮检异侧接收；(c) R分量炮检同侧接收；(d) R分量炮检异侧接收

(a) X分量炮检同侧接收；(b) X分量炮检异侧接收；(c) R分量炮检同侧接收；(d) R分量炮检异侧接收

5 结论

1) 考虑隧道空间的情况下，基于柱坐标系下的一阶弹性波速度−应力方程，采用高阶交错网格算法结合FCT技术，实现了柱坐标系下隧道空间典型不良地质模型的高精度全波场数值模拟，系统分析了其传播机理与波场特征，为隧道空间复杂波场的认识以及数据处理和解译提供了有效的依据。

2) 岩性分界面模型，当界面倾角为90°时，掌子面及地质界面间的多次反射波对有效波的识别会产生干扰，在数据处理和解译时应当加以区分。

3) 断层破碎带模型，炮检同侧接收的记录中，顶、底界面的反射纵波能量较强且同相轴较为连续，可较好的区分顶、底界面，但多次波以及波场叠加干扰会影响有效波的识别。

4) 溶洞顶底边界的反射纵波以及反射横波在记录中均有体现，但是溶洞底部反射纵波受波场叠加的影响，会对溶洞顶底边界的识别造成一定程度的影响，而反射横波能够有效的区分溶洞顶底边界，利于横波成像。

[1] 张杨, 杨君, 周黎明, 等. TSP在隧道工程施工中的常见干扰和对岩体裂隙水及软弱夹层等的预报研究[J]. 地球物理学进展, 2018, 33(2): 892−899. ZHANG Yang, YANG Jun, ZHOU Liming, et al. Common interference and the prediction of rock fissure water and weak interlayer in tunnel construction using TSP[J]. Progress in Geophysics, 2018, 33(2): 892−899.

[2] 周轮, 李术才, 许振浩, 等. 隧道综合超前地质预报技术及其工程应用[J]. 山东大学学报(工学版), 2017, 47(2): 55−62. ZHOU Lun, LI Shucai, XU Zhenhao, et al. Integrated advanced geological prediction technology of tunnel and its engineering application[J]. Journal of Shandong University (Engineering Science), 2017, 47(2): 55−62.

[3] 李术才, 刘斌, 孙怀凤, 等. 隧道施工超前地质预报研究现状及发展趋势[J]. 岩石力学与工程学报, 2014, 33(6): 1090–1113. LI Shucai, LIU Bin, SUN Huaifeng, et al. State of art and trends of advanced geological prediction in tunnel construction[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(6): 1090−1113.

[4] 赵永贵. 国内外隧道超前预报技术评析与推介[J]. 地球物理学进展, 2007, 22(4): 1344−1352. ZHAO Yonggui. Analysis and recommendation of tunnel prediction techniques at home and abroad[J]. Progress in Geophysics, 2007, 22(4): 1344−1352.

[5] ZHAO Yonggui, JIANG Hui, ZHAO Xiaopeng. Tunnel seismic tomography method for geological prediction and its application[J]. Applied Geophysics, 2006, 3(2): 69− 74.

[6] Dickmann T, Sander B. Drivage concurrent tunnel seismic prediction[J]. Felsbau-Rock and Soil Engineering, 1996(14): 406−411.

[7] 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.

[8] Neil D W, Haramy K Y, Hanson D H, et al. Tomography to evaluate site condition during tunneling[J]. Geotechnical Speccial Publicaion, 1999(90): 71−281.

[9] 刘江平, 程飞, 范承余, 等. 基于隧道空间全波场二维数值模拟与特征分析[J]. 岩土工程学报, 2012, 34(9): 1705−1711. LIU Jiangping, CHENG Fei, FAN Chengyu, et al. Two- dimensional numerical simulation of tunnel-based seismic full-wave fields[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(9): 1705−1711.

[10] 鲁光银, 熊瑛, 朱自强. 隧道反射波超前探测有限差分正演模拟与偏移处理[J]. 中南大学学报(自然科学版), 2011, 42(1): 136−141. LU Guangyin, XIONG Ying, ZHU Ziqiang. 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.

[11] CHENG Fei, LIU Jiangping, QU Niannian, et al. Two-dimensional pre-stack reverse time imaging based on tunnel space[J]. Journal of Applied Geophysics, 2014 (104): 106−113. DOI:10.1016/j.jappgeo.2014.02. 013.

[12] 凌飞, 肖宏跃, 朱夏乐, 等. 基于黏弹性介质的隧道地震二维正演模拟[J]. 长江科学院院报, 2015, 32(5): 121−126. LING Fei, XIAO Hongyue, ZHU Xiale, et al. Two- dimensional forward seismic modeling for tunnel based on viscoelastic medium[J]. Journal of Yangtze River Scientific Research Institute, 2015, 32(5): 121−126.

[13] 宋杰. 隧道施工不良地质三维地震波超前探测方法及其工程应用[D]. 济南: 山东大学, 2016. SONG Jie. The three-dimensional seismic ahead prospecting method and its application for adverse geology in tunnel construction[D]. Jinan: Shandong University, 2016.

[14] 王京. 隧道空间有限元地震波场模拟与逆时偏移成像[D]. 武汉: 中国地质大学, 2017. WANG Jing. Finite element seismic wave field simulation and reverse time migration based on tunnel space[D]. Wuhan: China University of Geosciences, 2017.

[15] Jetschny S, Bohlen T, Kurzmann A. Seismic prediction of geological structures ahead of the tunnel using tunnel surface waves[J]. Geophysical Prospecting, 2011, 59(5): 934−946.

[16] YANG Sitong, WEI Jiuchuan, CHENG Jiulong, et al. Numerical simulations of full-wave fields and analysis of channel wave characteristics in 3-D coal mine roadway models[J]. Applied Geophysics, 2016, 13(4): 621−630.

[17] Harmankaya U, Kaslilar A, Wapenaar K, et al. Locating scatterers while drilling using seismic noise due to tunnel boring machine[J]. Journal of Applied Geophysics, 2018, 152: 86−99.

[18] LIU Yutao, LIU Jiangping, CHENG Fei, et al. Full wave field finite difference modeling of tunnel space in polar coordinates[C]// International Geophysical Conference, Qingdao, China, 17-20 April 2017, Qingdao, China. Society of Exploration Geophysicists and Chinese Petroleum Society, 2017.

[19] 查欣洁, 高星, 王伟, 等. 隧道工程勘察中的超前预报成像方法研究[J]. 地球物理学报, 2018, 61(3): 1151− 1156. ZHA Xinjie, GAG Xing, WANG Wei, et al. Advanced prediction migration method research in tunnel engineering investigation[J]. Chinese J Geophys, 2018, 61(3): 1151−1156.

[20] Hiroshi Takenaka, Hiroki Tanaka, Taro Okamoto, et al. Quasi-cylindrical 2.5D wave modeling for large-scale seismic surveys[J]. Geophysical Research Letters, 2003, 30(21):2086.DOI:10.1029/2003GL018068.

[21] Graves R W. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite difference[J]. Bulletin of the Seismological Society of America, 1996, 86(4): 1091−1106.

[22] LIU Qinghuo, Sinha B K. A 3D cylindrical PML/FDTD method for elastic waves in fluid-filled pressurized boreholes in triaxially stressed formations[J]. Geophysics, 2003, 68(5): 1731−1743.

[23] LIU Q H. Perfectly matched layers for elastic waves in cylindrical and spherical coordinates[J]. The Journal of the Acoustical Society of America, 1999, 105(4): 2075− 2084.

[24] Mittet R. Free-surface boundary conditions for elastic staggered-grid modeling schemes[J]. Geophysics, 2002, 67(5): 1616−1623.

[25] Bohlen T, Lorang U, Rabbel W, et al. Rayleigh-to-shear wave conversion at the tunnel face: From 3D-FD modeling to ahead-of-drill exploration[J]. Geophysics, 2007, 72(6): T67−T79.

[26] 董良国, 黄超, 刘玉柱, 等. 溶洞地震反射波特征数值模拟研究[J]. 石油物探, 2010, 49(2): 121−124, 15. DONG Liangguo, HUANG Chao, LIU Yuzhu, et al. Numerical simulation of seismic wave propagation in cave carbonate reservoir[J]. Geophysical Prospecting for Petroleum, 2010, 49(2): 121−124, 15.

Full wave filed numerical simulation and analysis of tunnel space based on cylindrical coordinate system

LU Guangyin, LUO Shuai, ZHU Ziqiang, SHI Keliang, XIA Chengzhi

(School of Geosciences and Info-physic, Central South University, Changsha 410083, China)

Differing from ground seismic exploration, the special structure of tunnel space makes its wave field very complex. Full-wave field numerical simulation based on actual tunnel model is an important method to effectively understand the propagation law of wave field in tunnel space. Previous studies on this field are mostly based on Cartesian coordinate system. As a special geological model similar to air column, numerical simulation of full wave field in cylindrical coordinate system is more meaningful for understanding the characteristics of wave field in tunnel space. An arbitrary even-order staggered grid difference scheme was deduced, artificial truncated boundary and tunnel free boundary were designed, and numerical simulation and wave field characteristics analysis were performed for typical problematic geological models of tunnels such as lithologic interface, fault fracture zone and karst cave based on the first-order elastic wave velocity-stress equation in cylindrical coordinates. The results show that the staggered grid finite difference method based on the first-order elastic wave velocity-stress equation in cylindrical coordinates combined with FCT technology can achieve high-precision numerical simulation of the full-wave field in tunnel space. The wave field characteristics are in accordance with the kinematic and dynamic characteristics of the wave, which provides an effective basis for understanding the complex wave field characteristics in tunnel space and interpreting the data.

advance prediction; tunnel space; cylindrical coordinates; full wave field; numerical simulation

P631.4

A

1672 − 7029(2020)02 − 0388 − 08

10.19713/j.cnki.43−1423/u.T20190267

2019−04−08

国家自然科学基金资助项目(41974148)；湖南省安全生产监督管理局资助项目(201907)

鲁光银(1976−)，男，湖北宜昌人，教授，博士，从事工程地球物理勘探研究；E−mail：13975894898@139.com

(编辑 蒋学东)