公路隧道岩溶病害探测数值模拟研究*

钟安然，王应伟，陈海涛

(1.昭通市交通建设工程质量安全监督局，云南 昭通 657000； 2.昭通市宜昭高速公路项目指挥部，云南 昭通 657000；3.中铁大桥科学研究院有限公司，湖北 武汉 430034)

0 引言

岩溶灾害是目前国内外隧道建设面临的重难点问题。在岩溶区施工的隧道常发生突泥、突水、塌方等现象，造成人员伤亡、施工设备损坏，严重影响工期。因此，在隧道施工过程中及时进行掌子面前方及底板岩溶探测，可起到指导隧道安全施工的重要作用。目前用于隧道岩溶探测的方法主要有隧道超前地质预报法、探地雷达法、水平声波剖面法、超前钻孔法等[1]，不同方法具有各自适应性和缺点，探地雷达法因具有对施工影响小、预报效率高等特点，在隧道岩溶探测中广泛应用[2]，但在探地雷达工程检测和资料解释中，受现场检测条件、环境噪声干扰及技术人员经验制约，常导致在具体异常解释推断方面出现分歧，且缺乏足够的数值模拟和实测典型数据支持[3]。因此，笔者将数值模拟与工程实践相结合，为探地雷达技术在隧道岩溶病害探测中的应用进行研究。

1 基于GprMax软件的正演模拟基本原理

宏观上，所有的电磁现象均可通过Maxwell方程组进行描述，早在1966年，Yee[4]通过将Maxwell旋度方程引入空间离散方式，将其由微分方程转化为差分方程，并模拟出理想导体在时间域的电磁响应，经完善后发展形成时域有限差分(finite difference time domain，FDTD)模拟电磁场传播的数值计算方法。

假定在空间中一定区域内没有电磁场场源，且区域内介质是各向同性的，那么Maxwell方程组中的2个旋度方程在该特定区域内可写成[5-6]：

(1)

(2)

式中：E为电场强度(V/m)；H为磁场强度(A/m)；ε为介电常数(F/m)，σ为介质电导率(S/m)；μ为磁导率(H/m)；ρ为磁损耗磁阻率(Ω/m)；t为时间(s)。

在隧道岩溶探测中，模拟的对象是二维剖面，且探地雷达探测中主要利用TM(横磁模式)电磁波[3,6]，则Maxwell方程组中的2个旋度方程在直角坐标系中的分量形式可简化为[5]：

(3)

(4)

(5)

式中：Hx，Hy分别为磁场x，y向分量；Ez为电场z向分量。

由式(3)～式(5)可知，TM电磁波仅有Hx，Hy，Ez分量。运用Yee氏网格模型，利用中心差分代替对时间、空间坐标的微分，将连续变量离散化，即可推导出探地雷达正演模拟方程[7]。在二维情况下，Yee氏差分网格如图1所示。

图1 二维Yee氏差分网格示意

2 二维时域有限差分基本参数选取原则

2.1 空间步长的确定

在时域有限差分网格中，数值模拟的传播速度将随频率改变，导致非物理因素引起的脉冲波形畸变、人为的各向异性及虚假的折射现象，即出现数值色散现象，造成数值不稳定[7]。

二维空间中TM电磁波数值色散方程为[8]：

(6)

式中：kx，ky分别为波矢量沿x，y向的分量；ω为角频率；v为被模拟均匀介质中的光速；Δt为时间步长；Δx，Δy分别为x，y向空间步长。

当Δt，Δx，Δy均趋于0时，色散可减小至任意程度，但由于计算机内存空间及运算速度等因素限制，时间步长和空间步长不可能无限小，因此需选择合适的时间步长和空间步长。通常在时域有限差分中，网格空间步长Δl最大值不超过电磁波波长λ的1/10。因数值色散引起的误差是可接受的，通常取Δl=0.1λ，有效数字保留至小数点后三位即可，本研究案例据此取Δl=0.030m。

2.2 时间步长的确定

早在1975年，Taflove[9]探讨了Yee氏差分算法的稳定性问题，并给出了时间步长Δt的限定条件。对于二维TM电磁波正演问题，其对应的稳定条件为[10]：

(7)

式中：c为真空中的光速。

考虑到探地雷达正演计算的简便性，可令Δl= Δx=Δy，继而有：

(8)

2.3 吸收边界厚度的确定

时域有限差分法需在电磁场全部空间建立Yee氏网格计算空间，但由于计算机内存空间有限，不可能利用计算机直接在无限大的网格空间中计算电磁场，因此在计算过程中须在某处将网格空间截断，使之成为有限空间[11-12]。为使电磁波不在差分网格截断处发生明显反射，须设置吸收边界条件，将传播至边界处的电磁波吸收，进而确保计算机模拟的有限空间与自然界无穷大空间之间的差异达到最小[13]。

在PML边界条件中，层厚通常为3～9个网格单元厚度，本文在进行探地雷达时域有限差分正演模拟时，采用PML吸收边界条件，边界层选取10个网格单元。

3 隧道岩溶数值模拟分析

岩溶发育种类繁多，形态各异，成因不一，影响岩溶发育的主要因素包括岩性、岩层产状与构造、流水及地下水作用。岩溶发育过程可简化为：溶蚀裂隙→溶孔→溶洞，溶洞内充填的物质主要为空气、水、黏土(夹杂碎石)等。

3.1 单个溶洞模型的建立

为验证数值模拟方法的正确性，并了解灰岩中典型岩溶病害波场特征，对隧道底板中单个圆形溶洞内全部充填空气、全部充填土、全部充填水、上半部分充填空气下半部分充填水、上半部分充填空气下半部分充填土、上半部分充填水下半部分充填土的模型进行了正演模拟。各模型参数如表1所示，相对介电常数参考实际围岩情况而定。模型1～6模拟的地质模型如图2所示。

表1 单个溶洞模型参数

图2 模型1～6模拟的地质模型

3.2 单个溶洞正演模拟与波场特征分析

根据实际检测使用的参数，模拟时采用的中心频率f=100MHz，采样时间长度为150ns，道间距为0.3m。模型尺寸为7m×15m(深度×宽度)，模拟网格Δx=Δy=0.02m，总道数为45。PML吸收边界采用10个网格单元厚度，激励源选用雷克子波。模型1～6正演雷达剖面如图3所示。

图3 模型1～6正演雷达剖面

模型1雷达响应特征明显，正演雷达剖面主要表现为：洞顶(剖面长度7m、双程走时49ns处)附近呈明显的弧形正反射同相轴；洞底(剖面长度7m，双程走时52.5ns处)为较弱的负反射同相轴，且反射能量低于洞顶；洞顶、底界面波形未完全分离。模型2雷达响应特征明显，正演雷达剖面主要表现为：洞顶(剖面长度7m、双程走时49ns处)附近表现为能量较弱的弧形负反射同相轴；洞底(剖面长度7m、双程走时70ns处)为较弱的正反射同相轴；洞顶、底界面波形完全分离。

模型3雷达响应特征明显，正演雷达剖面主要表现为：洞顶界面(剖面长度7m、双程走时49ns处)雷达波反射同相轴能量较空洞和土洞强，呈弧形强负反射同相轴，且洞底界面(剖面长度7m、双程走时91ns处)反射波更清晰，并与顶界面完全分离。这是因为模型中设定的介质电导率小，雷达波在纯水中的衰减较慢，而水与围岩的介电常数差异大，洞顶界面处雷达波由围岩进入水中反射系数为负，极性发生反转，导致顶界面处呈强负反射。

由于雷达波在水中的波速最小，雷达波穿越溶洞的双程走时最大，故充水溶洞顶、底界面分离程度最大，然后为充土溶洞，最后为充气溶洞。这说明当溶洞内充填单一介质时，根据雷达剖面反射特征可识别出溶洞顶部位置，根据充填介质的不同，可较好地识别出溶洞底部位置，以便对溶洞尺寸进行分析。

由图3可知，洞顶(剖面长度7m、双程走时49ns处)界面反射均呈弧形同相轴形态，洞内2种充填介质的分界面(模型4剖面长度7m、双程走时52.5ns处，模型5剖面长度7m、双程走时52.5ns处，模型6剖面长度7m、双程走时91ns处)较易识别，但洞底反射界面较难识别。由此可推断，当洞内充填介质数量≥2时，可确定洞顶界面，但不易确定洞底界面，故无法准确对溶洞实际尺寸进行评估。

综上所述，模型1～6正演雷达波形特征如表2所示。

表2 模型1～6正演雷达波形特征

3.3 岩溶裂隙模型正演模拟与波场特征分析

根据常见的岩溶裂隙发育特征，在模型0的基础上添加1组全部充填土的岩溶裂隙模型，即模型7(见图4)。模型内设定1条斜向岩溶裂隙，裂隙宽度为1.0～1.5m。模型尺寸为7m×15m(深度×宽度)，模拟网格Δx=Δy=0.02m。模拟时采用的中心频率f=100MHz，采样时间长度为150ns，道间距为0.3m，总道数为45。PML吸收边界采用10个网格单元厚度，激励源选用雷克子波。

图4 模型7

模型7正演雷达剖面如图5所示。由图5可知，剖面内有1条明显的倾斜同相轴，同相轴呈强反射特征，整体连续，表明模拟结果与设定模型基本吻合，但正演雷达剖面中的裂隙倾角小于模型中设定的倾角，这是由于受裂隙倾角较大的影响，雷达记录的异常反射信息并非来自该点的正下方，从而产生偏移现象。为得到正确的解释结果，须进行偏移归位处理，常见的偏移处理方法包括Kirchhoff偏移、时域有限差分偏移、F-K偏移等[14]。

图5 模型7正演雷达剖面

4 隧道岩溶探测实例与分析

某在建高速公路隧道隧址区岩溶发育强烈，施工过程中偶有底板坍塌导致施工机械被困的安全事故发生，为探明隧道底板岩溶发育情况，采用PULSE EKKO PRO型探地雷达，配置100MHz天线进行探地雷达岩溶探测，测线纵向布置于隧道底板中线，长60m，时窗设定为350ns，点测模式，点距0.3m。

该段底板岩溶探测探地雷达剖面如图6所示。由图6可知，剖面内有1条明显的斜向同相轴，同相轴呈强反射特征，整体连续，同相轴下方伴有多个具有明显强反射特征的弧形绕射，初步推断本次探测范围内同时发育溶洞与裂隙，且岩溶裂隙向隧道底板延伸，在雷达剖面长度15.5m处，岩溶裂隙距隧道底板最小距离约为2m。

图6 底板岩溶探测探地雷达剖面

为验证正演模拟在实际探测中的有效性，根据实测雷达剖面初步推断结果，进行地质模型正演模拟，建立模型8，在其内设定1条斜向岩溶裂隙及14处埋深不一、半径各异的溶洞，模型总长60m，深16m，选取探地雷达天线的中心频率为100MHz，天线距为0.6m，道间距取0.3m，时窗设定为350ns，共模拟采集了200道数据，每道数据各4 947个采样点。考虑到实际探测过程中，天线不可能完全紧贴隧道底板，因此，模型中将天线置于远离地面0.02m高的位置处。由于该模型较大，为尽量节省计算时间，并保证模拟效果，选取的空间步长Δx=Δy=0.03m，模型8正演雷达剖面如图7所示。

图7 模型8正演雷达剖面

对比实测雷达剖面与模拟雷达剖面，发现如不考虑水平干扰波，模拟雷达剖面与实测雷达剖面特征基本相似。为验证探测结果，在实测雷达剖面长度15.5m处进行了钻孔及井中电视验证。井中电视结果如图8所示，由图8可知，从深度2.1m起，有明显溶洞出现，说明探地雷达探测结果准确，正演解释推断合理。

图8 井中电视结果

5 结语

1)通过分析GprMax正演模拟理论基础，依据数值色散方程给出了时域有限差分法空间步长确定原则。依据解的稳定性条件，确定了基本参数时间步长。根据PML边界条件，给出了吸收边界层厚选取的建议值。

2)采用GprMax软件对隧道底板下方不同类型溶洞模型进行了正演模拟。模拟结果表明，当隧道底板下方基岩中存在溶洞时，因围岩与溶洞充填物交界面两侧存在较大的电性差异，从而易形成强烈的反射波和绕射波，这种回波在时间剖面上表现为不同极性的弧形同相轴，弧形同相轴顶部为溶洞距隧道底板的最小距离；当溶洞内部为单一介质充填时，可通过识别洞顶、底界面位置评估溶洞尺寸；当溶洞内部充填的介质数量≥2时，可确定洞顶界面，但不易确定洞底界面，无法准确对溶洞实际尺寸进行评估；当基岩中存在岩溶裂隙发育带时，会产生相应的回波，常呈现为非水平的较连续同相轴特征；当岩溶裂隙倾角较大时，雷达剖面中显示的异常信息与实际情况存在偏移现象，需进行偏移归位处理。

3)探地雷达隧道岩溶探测工程实例解释结果得到了钻孔电视的验证，表明探地雷达正演模拟有助于探测成果图像中异常的识别，尤其是在灰岩地区，利用探地雷达进行隧道岩溶探测是有效的。