路基下伏地质缺陷弹性波数值模拟

李庆春，吴 华，周学明

1.长安大学地质工程与测绘学院，西安 710054

2.长安大学理学院，西安 710064

3.铁道第三勘察设计院集团有限公司，天津 300251

0 引言

我国西部地区，尤其是西北地区，存在特殊的自然地理环境和恶劣的地质环境，地下地质结构十分复杂，地质缺陷广泛发育，主要有黄土洞穴（湿陷性暗穴、陷穴，地下采空区）、地面塌陷、活断层、地裂缝、滑坡面、破碎带及软弱结构层等。这些隐伏的地质缺陷常常引发大量的地质灾害，给公路及铁路的施工和安全运营造成极大威胁。

在路基下伏地质缺陷探测方面，美国等西方发达国家以物探方法为主，而我国目前则以钻探为主、物探为辅［1－2］。国内近年来在利用地球物理勘探技术查明路基下伏地质缺陷方面做了大量的工作，提出了多种方法：蒋正红等［2］为了能较准确地揭示出路基病害的位置、规模及延展范围等，提出应用地球物理层析成像技术对铁路路基病害进行探测；孙忠弟［3］针对高等级公路下伏空洞勘探，利用电法、微重力、地震、地质雷达等方法进行了较为系统的研究；童立元等［1］从采空区探测、稳定性评价、治理及质量监控技术等方面对高速公路下伏采空区问题进行了综合探讨和分析，指出在采空区探测技术方面，有必要从理论上和实践上来优选最佳的探测技术组合，提高勘测精度；李庆春等［4］利用多道瞬态面波分析、浅层转换波地震技术探测穿过公路的湿陷性黄土暗穴，取得了良好的应用效果；李志辉等［5］采用CT速度成像技术，对岩溶、洞穴、软弱土层等病害存在部位进行定位，取得了良好的应用效果；韩永琦［6］、韩金河［7］采用地震折射层析法分别对青藏铁路路基开裂和石膏岩溶进行了探测，取得了较好的应用效果；罗维斌等［8］利用高密度电法对黄土路基下伏地质缺陷（地裂缝、黄土暗穴、断层破碎带、滑坡以及地下空洞等）进行模拟及实际探测，总结出针对不同缺陷的探测技术组合；张华等［9］将地震层析成像（CT）应用于公路采空区调查。

总之，在公路及铁路建设中，地球物理探测已显示并发挥了重要作用，在路面、路基、桥隧等探测检测方面做过相应工作，使用了包括瑞利面波、声波无损检测、浅层反射地震、隧道地震预测、探地雷达、直流电法、自然电场法等诸多技术。但是针对西部黄土地区路基下伏地质缺陷开展有效的综合地球物理探测，目前尚缺乏实用成熟的整套技术。

典型路基下伏地质缺陷的地球物理数值模拟是节约地球物理探测成本和验证各种地球物理方法探测路基下伏地质缺陷可行性和有效性的重要技术手段，但是，在路基下伏地质缺陷数值模拟方面，前人的研究［10－11］主要集中在缺陷的变形、沉降、稳定性和承载力模拟方面，关于模拟路基下伏地质缺陷对弹性波响应方面的研究为数不多。周学明等［12］采用交错网格高阶有限差分数值模拟技术对一个包含有多种地质缺陷的综合地质模型进行了模拟计算，通过模拟表明工程地震勘探在探测浅层断层、裂缝、溶洞等地质缺陷方面具有有效性。

笔者拟对西部黄土地区路基下伏地质缺陷进行地球物理数值模拟，以验证各种地球物理方法探测路基下伏地质缺陷的可行性和有效性，以期为实际施工中地球物理方法技术组合以及观测系统的优化设计提供技术支持。

1 地震波场记录有限差分数值模拟

数值模拟可模拟研究地震波在地下各种介质中的传播规律，是研究地下构造和和岩性等各种地质因素与地震波运动学和动力学响应特征之间关系的重要技术手段，对于认识和研究地震波在介质中的传播具有重要作用。地震数值模拟已广泛应用在地震勘探的资料采集、处理和解释的各个阶段。笔者根据一阶速度－应力弹性波方程，采用交错网格高阶差分方法，引入高效的完全匹配层（PML）技术解决边界吸收问题，可进行黏弹介质、双相介质和各向异性介质的弹性波数值计算。

1.1 波动方程

根据弹性介质位移、应力和应变之间的关系，可以推导出二维均匀各向同性介质弹性波一阶速度－应力方程［13－16］：

式中：vx、vz分别为质点位移速度的水平和垂直分量；σx、σz分别为x和z方向的正应力；τzx为切应力分量；ρ为密度；t为时间；λ和μ为拉梅系数。

1.2 差分格式

1976 年，Madariaga［17］提出了一种较为先进的波动方程交错网格有限差分解法，并首次利用交错网格来模拟弹性介质内圆形扩展破裂产生的波动。该方法与常规差分格式相比，具有局部精度高、收敛速度快等优点，且计算量和占用内存与常规方法相当。交错网格差分格式不仅要求空间网格交错，而且要求在时间上也要交错。交错网格中速度和应力节点如图1所示。

图1 交错网格中速度和应力节点示意图Fig.1 Sketch map of velocity and stress node in staggered grid

时间上的高阶差分会极大地增加计算量，在地震勘探所需的时间采样上时间频散并不明显，增大时间步长可能会导致算法不稳定。而空间频散是数值频散的主要部分，对空间导数需采用高阶差分格式来有效抑制数值频散。考虑到运算效率和精度，笔者选择了时间二阶、空间十阶的差分格式［12］：

其中：

将式（2）、（3）代入式（1）即可得到弹性波方程的高阶差分格式。

1.3 边界条件

1）自由边界条件

在水平自由表面上（z＝0），应力分量必须满足σzz＝0，τzx＝0。由于采用的是交错网格，所有的速度及应力分量不可能同时位于自由表面，设定vx，σxx，σzz位于自由表面。当采用高阶差分格式时，还要求τzx，σzz在高阶差分算子上一定是关于z＝0反对称的，即虚像原理，并在自由界面上设置N层真空层，即纵、横波速度分别为0（vP＝0，vS＝0）。则在计算速度分量时满足：

其中：i表示x方向上的离散值；k表示时间上的离散值。

2）完全匹配层吸收边界条件

采用完全匹配层吸收边界条件吸收来自人工边界的反射，如图2所示。

将vx分解为x方向部分和z方向部分，则相应的匹配层方程［12］可写为

其中，d（x）、d（z）采用 Coffino等［18］导出的如下衰减模型：

图2 完全匹配层吸收边界示意图Fig.2 Sketch map of perfect matched layer absorbing boundary

式中：vPmax为最大纵波速度；δ为匹配层宽度；R为理想的边界反射系数（这里取为10－6）；系数a＝0.25，b＝0.75。当d（x），d（z）不等于0时表示衰减；d（x），d（z）为0时表示不衰减。

利用完全匹配层作为吸收边界的基本做法是在研究区域的四周引入完全匹配层，左右边界和上下边界采用完全匹配层吸收边界。左右匹配层在x方向衰减，在z方向不衰减，即d（x）≠0、d（z）＝0。上下匹配层z方向衰减，x方向不衰减，即d（x）＝0，d（z）≠0。角点区域x方向和z方向都衰减，即d（x）≠0，d（z）≠0。模拟瑞利波场时，研究区域的上边界取为自由边界。

3）稳定性条件

在模型计算中，空间网格尺寸和时间间隔的选择至关重要，既要保证计算精度，还得兼顾计算效率。对于空间十阶差分格式，保证整个计算区域能稳定计算的高阶交错网格差分显式的稳定性条件为

特别地，对于Δx＝Δz的情况，稳定性条件为

2 典型地质缺陷地球物理数值模拟及分析

为考查不同地质缺陷的地震波场响应特征，检验不同地震方法的探测能力，设计了一个综合地质地球物理模型，如图3所示。该模型包含凹陷、断层、地裂缝、地下低速体、软弱夹层及滑坡等典型地质缺陷。

模型长400m，深80m，共分为7层。在地面坐标50m处有一个隐伏地裂缝，裂缝宽0.25m，从地下2m开裂到地下40m。在地面坐标182m、地下20m处有一个宽和高均为4m的规则地下低速体。另外，在地面坐标224m、地下45m处有一个宽7 m、高10m的不规则低速体。利用弹性波场有限差分模拟程序，计算了多分量地震记录和面波地震记录各150炮，之后分别进行了P－P、P－SV波成像和高阶面波反演横波速度成像。

模拟计算时，模型物性参数纵横波速度之比相对于实际浅层地质情况略微偏小，但不影响各种物探方法成像效果的验证。

2.1 反射地震波场模拟及成像

波场模拟时采用高阶交错网格有限差分和PML吸收边界技术，模型长×宽为400m×80m，取网格间距Δx＝Δz＝0.25m，时间步长Δt＝0.1 ms，取主频为80Hz的Ricker子波作为震源函数。合成的单炮记录为纵波源激发，单边接收，每炮48道，采样率为0.5ms，采样长度为512ms，每道1 024个采样点，道间距为2m，炮间距为2m，最小偏移距为4m。

图4为数值模拟产生的典型单炮地震记录（前两炮），由图4可知，来自各层的反射可以清晰识别。

图5和图6分别是采用常规反射叠加处理技术产生的P－P、P－SV波成像剖面。由图5和图6可知，由于模型复杂，在水平叠加剖面上，还存在绕射现象，经叠后偏移处理，成像剖面的效果大大改善。再经时深转换，反射层位与真实模型层位之间的最大误差不超过2m。在进行P－SV波偏移时，采用等效偏移距偏移（EOM）方法，其中P－SV波的偏移速度模型是通过转换波高精度非双曲时距曲线“双曲化”处理后，再按常规速度分析方法进行速度扫描分析得到的，具体方法原理见文献［19］。

图3 典型地质缺陷综合地质地球物理模型Fig.3 A comprehensive geological and geophysical model with typical geological defects

图4 前两炮典型单炮反射地震记录Fig.4 Typical reflection seismic records of the first two shots

真实模型中的凹陷、断层、浅层低速体、软弱夹层及滑坡等主要地质缺陷在反射地震成像剖面上均能清晰成像，且对应关系良好。真实模型中的地裂缝由于尺度较小且近垂直展布，不能形成有效反射，在成像剖面中反映不明显。同时，较深层的地下低速体在成像剖面中反映也不明显，其原因可能是该低速体没有有效的水平反射界面、埋深较深、侧面倾角较陡，因此，在地表无法接收来自低速体的有效反射信息所致。

图5 模型地震记录P－P波成像剖面Fig.5 Section from P－P wave imaging of the model

2.2 瑞利波场模拟及成像

弹性波模拟时，上、下、左、右4个边界均采用吸收边界，因此模拟记录中不出现面波（被边界吸收）。面波模拟时，仍采用高阶交错网格有限差分技术，取上边界为自由边界，其他边界为PLM吸收边界。模型长×宽为400m×80m，取网格间距Δx＝Δz＝0.25m，时间步长Δt＝0.05ms，取主频为25Hz的Ricker子波作为震源函数。合成的单炮记录为纵波源激发，单边接收，每炮48道，采样率为0.5ms，采样长度为512ms，每道1 024个采样点，道间距为2m，炮间距为2m，最小偏移距为4m。

图7为数值模拟产生的典型单炮面波地震记录：x分量中波场成分简单，面波发育；z分量中波场成分相对复杂，但面波能量仍占主要地位。

图8是z分量记录的频散曲线提取结果与理论频散曲线对比图。由图8可知，由波场记录提取得到的频散曲线与理论频散曲线基本上是一致的，说明了面波记录的正确性。二者存在微小误差的原因是波场模拟时所建模型在横向方向界面是变化的，而理论频散曲线的计算是基于水平层状介质基础上得到的。另外，由图8可知，由模拟记录提取得到的频散曲线最低频率成分为4Hz左右，对应的瑞利波相速度约为400m/s；因此瑞利波的最大波长（波长等于相速度/频率）为100m左右。根据半波长原理，可估算出本文所给模型瑞利波的最大探测深度为40～50m。

图6 模型地震记录P－SV波成像剖面Fig.6 Section from P－SV wave imaging of the model

图7 第一炮典型单炮面波地震记录Fig.7 Typical surface wave records of the first shot

图8 z分量记录频散曲线的提取结果与理论频散曲线对比图Fig.8 Extracted dispersion from z component and its theoretical dispersion curves

图9是由z分量面波记录提取的频散曲线反演得到的横波速度断面图。由图9可知，模型中的凹陷、地裂缝、断层、浅层低速体、软弱夹层及滑坡等主要地质缺陷均有反映，尤其是反射剖面无法成像的地裂缝在面波反演速度断面中有明显的反映。与反射地震成像相比，面波成像在浅层结构探测中具有更好的纵、横向分辨能力，分辨率可达1m左右，但面波的反演深度有限。因此，对于不同深度、不同规模的地质体，应采用不同的地球物理探测方法，如利用多道瞬态面波解决浅层问题，用反射地震探测更深的目标，用多波地震解决裂缝等精细构造问题，以达到对地质缺陷进行立体探测与定位的效果。

3 结论与讨论

1）典型路基下伏地质缺陷的数值模拟，是考查不同地质缺陷地震波场响应特征、检验不同地震方法的探测能力、验证其可行性和有效性的重要技术手段。因此，在正式布置勘探之前，对地质缺陷进行波场数值模拟是非常必要的，可以最大限度地降低探测成本。

2）常规反射地震P波和SV波成像技术对真实模型中的凹陷、断层、浅层低速体、软弱夹层及滑坡等主要地质缺陷均能清晰成像，且对应关系良好，最大误差不超过2m。真实模型中的地裂缝由于尺度较小且近垂直展布，不能形成有效反射，在成像剖面中反映不明显。

3）与反射地震成像相比，面波成像在浅层结构探测中具有更好的纵、横向分辨能力，但探测深度有限。因此，对于不同深度、不同规模的地质体，应采用不同的地球物理探测方法，如：利用多道瞬态面波解决浅层问题，用反射地震探测更深的目标，用多波地震解决裂缝等精细构造问题，以达到对地质缺陷进行立体探测与定位的效果。

4）数值模拟结果验证了地球物理方法探测路基下伏地质缺陷的可行性，表明不同的地球物理探测方法对地质缺陷的探测效果和适用范围不同。试验结果对实际施工中地球物理方法技术组合以及观测系统的优化设计具有重要的指导意义。

图9 面波反演横波速度断面图Fig.9 Shear wave velocity section from surface wave inversion

（References）：

［1］童立元，刘松玉，邱钰，等.高速公路下伏采空区问题国内外研究现状及进展［J］.岩石力学与工程学报，2004，23（7）：1198－1202.Tong Liyuan，Liu Songyu，Qiu Yu，et al.Current Research State of Problems Associated with Mined－Out Regions Under Expressway and Future Development［J］.Chinese Journal of Rock Mechanics and Engineering，2004，23（7）：1198－1202.

［2］蒋正红，蒋伟宏.铁路路基病害地球物理层析成像方法简介［J］.路基工程，1996（3）：25－26.Jiang Zhenghong，Jiang Weihong.Introduction to Geophysical Tomography Methods on Rail－Bed Disasters［J］.Subgrade Engineering，1996（3）：25－26.

［3］孙忠弟.高等级公路下伏空洞勘探、危害程度评价及处治研究报告集［M］.北京：科学出版社，1999.Sun Zhongdi.Reports of Cavity Exploration Under High－Grade Highway， Hazard Evaluation and Treatment［M］.Beijing：Science Press，1999.

［4］李庆春，刘永华，冯兵，等.湿陷性黄土暗穴的综合地球物理探测试验研究［J］.公路交通科技，2005，22（6）：138－142.Li Qingchun，Liu Yonghua，Feng Bing，et al.Integrated Trial Studies on Geophysical Exploration for Hidden Loess Cavities［J］.Journal of Highway and Transportation Research and Development，2005，22（6）：138－142.

［5］李志辉，张立，刘争平.CT技术在铁路路基病害勘探中的应用［J］.物探化探计算技术，2006，28（2）：137－141.Li Zhihui，Zhang Li， Liu Zhengping.The CT Technique Application to the Survey of Rail－Bed Disasters［J］.Computing Techniques for Geophysical and Geochemical Exploration，2006，28（2）：137－141.

［6］韩永琦.地震折射层析法在青藏铁路路基开裂调查中的应用［J］.工程地球物理学报，2009，6（1）：1－5.Han Yongqi.The Application of Investigating the Roadbed Cracks of the Qinghai－Tibet Railway Using the Seismic Refraction Tomography Method［J］.Chinese Journal of Engineering Geophysics，2009，6（1）：1－5.

［7］韩金河.地震折射层析法在青藏铁路石膏岩溶勘察中的应用［J］.铁道勘察，2009（4）：18－20.Han Jinhe.Application of Seismic Refraction Tomography Method in Geologic Investigation of Anhydrite Karst Region on Qinghai－Tibet Railway［J］.Railway Investigation and Surveying，2009（4）：18－20.

［8］罗维斌，李庆春.电法在公路下伏地质缺陷精细探测中的应用［J］.工程勘察，2010（6）：78－85.Luo Weibin，Li Qingchun.Electrical Methods Applied for Road Failures Fine Detection［J］.Journal of Geotechnical Investigation ＆Surveying，2010（6）：78－85.

［9］张华，陈小宏，刘松.地震勘探技术在公路采空区调查中的应用［J］.地下空间与工程学报，2011，7（2）：371－374.Zhang Hua，Chen Xiaohong，Liu Song.Application of Seismic Prospecting Technology in Investigation on the Goaf Under Highway［J］.Chinese Journal of Underground Space and Engineering，2011，7（2）：371－374.

［10］王秋生.公路路基下伏空洞的数值分析［J］.岩土力学，2005，26（增刊）：215－218.Wang Qiusheng.Numerical Aanalysis of an Excavation Under Roadbed［J］.Rock and Soil Mechanics，2005，26（Sup.）：215－218.

［11］邓杨.隐伏溶洞山区公路路堤稳定性数值模拟分析与处治［J］.北方交通，2013（6）：23－26.Deng Yang.Numerical Simulation Analysis and Treatment on the Embankment Stability of Mountain Road Hidden Karst Caves［J］.Northern Communications，2013（6）：23－26.

［12］周学明，李庆春，李志华，等.地震波交错网格高阶差分数值模拟研究［J］.铁道工程学报，2011（8）：1－12.Zhou Xueming，Li Qingchun，Li Zhihua，et al.Research on Staggered－Grid High－Order Finite－Difference Numerical Simulation of Seismic Wave［J］.Journal of Railway Engineering Society，2011（8）：1－12.

［13］Bohlen T，Saenger E H.Accuracy of Heterogeneous Staggered－Grid Finite－Difference Modeling of Rayleigh Waves［J］.Geophysics，2006，71（4）：109－115.

［14］周竹生，刘喜亮，熊孝雨.弹性介质中瑞雷面波有限差分法正演模拟［J］.地球物理学报，2007，50（2）：567－573.Zhou Zhusheng，Liu Xiliang，Xiong Xiaoyu.Finite－Difference Modelling of Rayleigh Surface Wave in Elastic Media［J］.Chinese Journal of Geophysics，2007，50（2）：567－573.

［15］周晓华，陈祖斌，曾晓献，等.交错网格有限差分法模拟微动信号［J］.吉林大学学报：地球科学版，2012，42（3）：852－857.Zhou Xiaohua，Chen Zubin，Zeng Xiaoxian，et al.Simulation of Microtremor Using Staggered－Grid Finite Difference Method［J］.Journal of Jilin University：Earth Science Edition，2012，42（3）：852－857.

［16］杨庆节，刘财，耿美霞，等.交错网格任意阶导数有限差分格式及差分系数推导［J］.吉林大学学报：地球科学版，2014，44（1）：375－385.Yang Qingjie，Liu Cai，Geng Meixia，et al.Staggered Grid Finite Difference Scheme and Coefficients Deduction of Any Number of Derivatives［J］.Journal of Jilin University：Earth Science Edition，2014，44（1）：375－385.

［17］Madariaga R.Dynamics of an Expanding Circular Fault［J］.Bulletin of the Seismological Society of America，1976，66（3）：639－666.

［18］Collino F，Tsogka C.Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media［J］.Geophysics，2001（1）：294－307.

［19］沈鸿雁，李庆春.转换波高精度非双曲时距曲线“双曲化”处理［J］.石油地球物理勘探，2007，42（3）：251－255.Shen Hongyan， Li Qingchun.High－Precision Hyperbolic Simplified Processing for Non－Hyperbolic Time－Distance Curve of Converted Wave［J］.Oil Geophysical Prospecting，2007，42（3）：251－255.