基于隧道超前预报的GPR正演模拟及应用分析

陈 松，周黎明，罗士新

(1.武汉地质调查中心勘查技术室，武汉 430205;2.长江科学院 水利部岩土力学与工程重点实验室，武汉 430010)

1 研究背景

我国地质预报工作开展相对较晚，在20世纪70年代才开始进行隧道地质预报的应用和研究。20世纪80年代，在南岭隧道、大瑶山隧道的掘进过程中开展了地质预报的试验研究，且大瑶山隧道较成功地引用了地震反射法等［1］。C.W.Zang(2004)等［2］把地质雷达法、TSP法、超前钻探法集合起来，对围岩的完成性、强度以及水系发育的分级等进行了研究。齐传生(2004)等［3］在圆梁山隧道中采用了地质雷达法、TSP202法、HSP法、红外探测法、超前钻等方法综合分析，总结出了各种地质预报的特点、工序。曲海峰(2006)等［4］以TSP法、地质雷达法、地质调查法为手段，研究了综合预报的方法。由上可以看出地质雷达是隧道预报工作者最常选用的一种方法，因此，如何能利用好这一物探方法是预报的前提和基础。

近年来也有不少学者做过探地雷达正演模拟相关的研究工作，如邓世坤等(1993)［5］给出了用射线追踪法合成的二维地电模型的探地雷达图像，用衍射波公式合成了二维管状模型的雷达图像，并与物理模拟作了比较。吴宝杰(2007)［6］、赵宇(2009)［7］详细介绍了地质雷达的原理及二维正演模拟。赵峰(2012)等［8］分析了探地雷达在探测隧道衬砌空洞的有效性以及衬砌空洞的雷达波反射特征。胡俊(2012)等［9］对比分析了探地雷达对多年冻土区厚层地下冰的分布、埋藏深度、赋存情况以及对冻土类型识别的有效性。孙忠辉(2013)等［10］利用软件GprMax模拟地质雷达在初衬中的检测是准确可行的，能够清晰地反映隧道衬砌中钢拱架和钢筋的成像特征。

在隧道环境下，多数正演模拟的对象主要集中在隧道的衬砌上，使用的天线频率高，模拟的厚度小，消耗的时间短，最后的效果比较理想，而将隧道掌子面前方的围岩作为正演模型研究对象的相对较少。因此，本文主要从大尺度的模型出发，通过设计预报的地质模型，使用低频天线来探讨地质雷达在不同模型下的波场响应特征及其对应的解释方法，并通过实测数据的应用效果对比分析，证明地质雷达在隧道环境下进行超前预报的可靠性。

2 地质雷达正演基本原理

地质雷达GPR(Ground Penetrating Radar)是通过发射天线将高频电磁波以脉冲形式从隧道的掌子面传递到隧道前方的围岩中，经过围岩界面的反射返回隧道掌子面，并被雷达的接收天线接收。通过对回波信号进行分析解释，就可以获得隧道掌子面前方的围岩信息。本文中正演的基本原理是基于1966年Kane S.Yee提出的 FDTD(Finite Difference Time Domain)差分算法和 PML(perfectly matched layer)吸收边界条件，其中差分算法是把Maxwell方程组用中心差分近似替代，该方法运算速度快，可以解决各种复杂点的散射［11－12］。

假设介质是各向同性的，在二维空间下，若Z方向上只有电场分量Ez，磁场分量Hz为0，那么麦克斯韦的旋度方程可表示为

式中:H为磁场强度(A/m);E为电场强度(V/m);μ为磁导率(H/m);σ为电导率(s/m);ε为介电常数(F/m);σm为导磁率(Ω/m)。

K.S.Yee根据上面的3个偏微分方程引入FDTD算法，通过在空间建立矩形的差分网格，一个网格节点对应的时变参量为f(x，y，t)，则可以得到FDTD表达式为

式中:i，j为x，y方向的网格号;Δx，Δy为矩形网格分别沿着x，y方向的空间步长;Δt为时间步长。麦克斯韦方程要保持稳定性和收敛性，则需要满足公式

考虑到离散化后的频散现象，如果通过减小网格而抑制频散的话，就会增大相应的计算量，所以，在考虑压制频散与控制计算量的情况下，可以使网格尽量缩小［13－15］。

网格步长要小于波传播时最小波长 λmin的1/10，则有

3 隧道模型下的正演模拟

3.1 吸收边界

受到计算量的影响，用FDTD计算求解时，一般只求取空间有限区域的迭代解，在实际的计算中，通常用有限计算空间来代替无限的空间。实际应用表明，PML在大的入射角上吸收的效果很好，在一般的计算中该算法会显得很复杂，占用很大的内存［16－19］。

3.2 模型参数设置

隧道模型参照隧道开挖面，设计了2个侧壁和1个掌子面。下面将给定模拟的相关参数。我们选取的天线频率f0为50 MHz。根据电磁波在介质中传播的最小波速Vmin=c/=0.6 × 108，则有模拟的最高频率 fm=3f0=150 MHz，λmin=0.4 m。根据式(3)和式(4)，有稳定性条件 Δx≤0.04m，Δy≤0.04m，Δt≤30 ps，给定的时间窗口为400 ns。根据模拟的效果，本节选取的网格步长为Δx=Δy=0.01 m。天线间距为0.22 m，移动步长为0.20 m。天线探测方向为从掌子面的左端测向右端。模型中的几种介质电性参数见表1。

表1 试算模型电性参数Table 1 Dielectric parameters for model calculation

本文采用的是GprMax2D软件进行正演分析，给出了9种地质预报的模型，具体参数见表2。

表2 试算模型属性Table 2 Properties of model for the calculation

通过试算并根据隧道超前预报的场地条件，本文设计的隧道正演模型总长为30 m，总宽度为20 m，左右侧壁宽度均为5 m，掌子面处于已开挖的5 m处，开挖空间长5 m，宽10 m。模型Ⅰ至模型Ⅵ为设计的简单模型，分别为距掌子面前方10 m处发育有水平状的破碎带、距掌子面5 m处垂直掌子面的破碎带、左倾30°的破碎带、右倾30°的破碎带、中心距掌子面10 m处发育有半径为1 m的空溶洞及充水溶洞;模型Ⅶ至模型Ⅸ为花岗岩地层、页岩地层、灰岩地层与水系(水系带中部为梯形石块)、溶洞(分别充满空气、水、破碎块)、透镜体的组合模型。以上设计的模型在隧道开挖中会经常遇见，具有代表性。隧道围岩介质中花岗岩的相对介电常数为5，电导率为0.001(s/m);页岩的相对介电常数为7，电导率为0.01(s/m);灰岩的相对介电常数为8，电导率为1×10－9(s/m);破碎带的相对介电常数为25，电导率为0.003(s/m);透镜体参照破碎带的电性参数。隧道模型见图1。

3.3 正演模拟计算及分析

把模型参数文件导入到GprMax2D软件进行计算，可以得到正演模拟数据;再把得到的数据文件导入Fortran程序，经过增益恢复、能量均衡等初步处理生成*.grd格式数据，利用Surfer软件成图即完成整个正演的工作，成果图如图2所示。

从图2正演成果分析可知:地质雷达识别异常体的顶界面位置非常准确，如模型Ⅰ顶界面的埋深10 m、模型Ⅱ的顶界面埋深5 m、模型Ⅴ及模型Ⅵ溶洞的顶界面埋深9 m、模型Ⅷ组合溶洞埋深7 m等位置均与正演结果对应，但是由于电磁波受到破碎带、水等介质的干扰及吸收影响，地质雷达对异常体下界面位置判断相对有出入;同时地质雷达对断层破碎带的倾角判断也比较敏感，如模型Ⅰ、模型Ⅲ、模型Ⅳ正演结果分别对应了水平层状倾向、左倾30°、右倾30°的角度变化，模型 IX中倾斜断层面(图1中标注为9)与正演的同相轴(图2中标注为9)位置也比较吻合，由于受到上层透镜体的反射吸收影响，断层面反射波同相轴的能量较弱;地质雷达对圆形溶洞的形状预测比较准确，如模型Ⅴ、模型Ⅵ、模型Ⅷ正演剖面中反射波的弧形顶部对应了溶洞的顶部形状，但是对于模型Ⅷ中同时出现3个相邻的溶洞识别效果一般;地质雷达对于岩层的界面位置识别也较为准确，如模型Ⅶ的页岩与灰岩界面3、模型Ⅷ的花岗岩与页岩界面5及页岩与灰岩界面6对应了正演剖面中的位置(图2中分别标注为3，5，6);由于模型吸收边界没有起到对隧道侧壁边界反射波完全吸收的作用，边界拐点的绕射和散射干扰等对最后的正演结果有一定的影响。因此，地质雷达在隧道环境下的正演结果与理论模型对应得比较一致，通过以上分析电磁波的波场特征可以为解释实测的雷达数据作理论依据，对提高解释的精度有很大帮助。

图1 隧道模型示意图Fig.1 Schematic diagram of tunnel models

图2 正演成果Fig.2 Results of forward simulation

图3为从模型Ⅰ至模型Ⅸ正演结果中提取异常体中心位置的单道数据，通过对比分析总结有:电磁波遇到断层破碎带及溶洞等界面的反射，电磁波的相位会发生变化，如图3中红色标注的为地质体引起的相位变化，模型Ⅰ—Ⅳ、模型Ⅵ、模型Ⅶ电磁波由介电常数小的介质进入介电常数大的介质，波相位由正变负，而模型Ⅴ电磁波从围岩进入溶洞的空气介质，介电常数由大变小，对应的波相位由负变正，因此，根据波相位的正负相对变化也可以预测掌子面前方围岩的变化情况。同时透过破碎带、透镜体区域后的电磁波衰减很快，波的振幅变小，随着深度的增加能量已经很弱。在后期资料解释中，我们要结合单道数据综合判断异常体引起的波形特征，提高解释的精度。

对比图2和图3可以发现，地质雷达在理论的隧道模型下预测掌子面前方的地质异常体比较准确，波形图上能够定性地反映地质体的形态，并能定量地给出埋深及发育的规模、倾向倾角等重要的信息，同时也可以根据单道数据的波相位变化分析反射层两端岩体的介电常数变化，进而推断掌子面前方岩性的变化。因此，通过正演模拟，从理论上验证了地质雷达在隧道环境下做超前预报的可靠性，同时也为实际应用提供了分析波形的依据。

通过大尺度的模型计算，可以很明显地看到计算量和计算时间在大大增加，一个剖面的计算时间大约是23 h，这也是不同于对衬砌模拟(计算1个剖面大约数分钟的时间)的地方。同时模型尺度变大，电磁波会随着深度的增加迅速衰减，50 MHz的雷达天线探测的有效深度在30 m以内。

图3 单道数据振幅图Fig.3 Amplitude of single-channel data

4 GPR实测资料特征分析

隧道Ⅰ:此隧道区域地质露头多为牛屋组(Ptn)中厚层状风化岩层的岩体，其节理裂隙比较发育。布置的雷达测线在掌子面上的长度为13 m，雷达的天线中心频率是50 MHz，天线间距1.5 m，天线移动步长0.25 m，时间窗口长度1 200 ns。

隧道Ⅱ:此隧道区地貌单元属构造剥蚀低山—丘陵地貌，隧道走向呈现近北东—南西向，与山脊线走向近于平行。该区大地构造属于昆仑—秦岭纬向构造体系东秦岭南亚带，该构造体系是区内成生发展最早的构造体系，为区内主要构造格架，其他后期构造在此构造之上叠加、穿切、横跨，构造形态复杂，褶皱发育，断裂交织，具有多期活动特征，是个相对复杂的构造地区。布置的雷达测线长度为5.5 m，天线的中心频率为100 MHz，天线间距1 m，天线移动步长0.25 m，时间窗口长360 ns。

从图4(a)雷达剖面分析，可以发现在深度位置17～20 m处出现较为明显的连续的强能量同相轴，且同相轴的展布范围较大，与本文正演的模型I的波形特征比较接近，推断此处围岩的节理裂隙较发育，有破碎带的存在，且破碎带呈现平行掌子面的发育方向。后期隧道开挖过程中，掌子面前方17 m后岩层比较破碎，掉块严重，说明岩体节理裂隙发育，跟探测的结果比较一致，达到了超前预报的效果。同时也可以看到，雷达剖面深度在30 m以后，同相轴的能量已经很弱，无法提供有效的信息。

根据图4(b)地质雷达剖面图，可以比较清楚的看到在掌子面右上方向前50～80 ns的区段内(深度2.8～4.0 m，参考波速0.1 m/ns)，同相轴呈弯曲抛物线状态，与本文模型Ⅵ正演的波形特征非常接近，而且信号振幅较强，与周围波形存在明显差异，该异常区由岩溶洞穴带形成，结合隧道地质调查，初步判断此区段为岩溶洞穴。经过后期隧道开挖，验证了此位置预测的洞穴的存在，达到了预报的效果。

图4 隧道雷达剖面Fig.4 GPR image for tunnels

5 结论

地质雷达在短距离内识别隧道掌子面前方的断层破碎带、溶洞等地质体上高效、快捷、准确，本文不仅从理论模型上进行了验证，采集的实测数据也说明了其可靠性。根据正演及应用实例分析可以总结以下几点:

(1)利用GprMax2D软件与Fortran语言可以很好地实现隧道超前预报模型下的数值模拟，效果理想。

(2)分析正演成果图可知，电磁波在识别破碎带、溶洞的顶界面位置、发育方向及形状上准确，波相位会随着破碎带两端介质电性差异产生变化，综合这2点可以很好地获取掌子面前方反射波带回来的围岩信息，为探明隧道前方围岩情况提供很好的理论支撑。

(3)采集的2段隧道超前预报数据分别验证了水平发育方向破碎带、溶洞下的波场特征，反射波形特征给解释人员提供了很重要的参考资料，可指导施工的顺利进行。

［1］徐济川，黄少霞.大瑶山隧道的突泥涌水机制［J］.铁道工程学报，1996，(2):83－89.(XU Ji-chuan，HUANG Shao-xia.Mechanism of Mud and Water Burst in Dayao Mountain Tunnel［J］.Journal of Railway Engineering Society，1996，(2):83－89.(in Chinese))

［2］ZANG C W，HUANG H W，ZHANG Z X.Forecasting the Strata Condition of a Long Road Tunnel by Using Fuzzy Synthetic Judgment［J］.International Journal of Rock Mechanics Sciences，2004，41(3):406－407.

［3］齐传生，王洪勇.圆梁山隧道综合超前地质预报技术［J］.铁道勘察，2004，(5):52－56.(QI Chuan-sheng，WANG Hong-yong.Technique of Comprehensive Geological Prediction in Yuanliangshan Tunnel［J］.Railway Investigation and Surveying，2004，(5):52－56.(in Chinese))

［4］曲海峰，刘志刚，朱合华.隧道信息化施工中综合超前地质预报技术［J］.岩石力学与工程学报，2006，25(6):1246－ 1250.(QU Hai-feng，LIU Zhi-gang，ZHU He-hua.Technique of Synthetic Geologic Prediction Ahead in Tunnel Informational Construction［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25(6):1246－1250.(in Chinese))

［5］邓世坤，王惠濂.探地雷达图像的正演合成与偏移处理［J］.地球物理学报，1993，36(4):528－535.(DENG Shi-kun，WANG Hui-lian.Forward Synthesizing and Migration Processing for GPR Image［J］.Chinese Journal of Geophysics，1993，36(4):528－535.(in Chinese))

［6］吴宝杰.探地雷达二维正演模拟及其工程实例［D］.杭州:浙江大学，2007.(WU Bao-jie.Two-dimensional Forward Simulation of Ground Penetrating Radar and Engineering Examples［D］.Hangzhou:Zhejiang University，2007.(in Chinese))

［7］赵 宇.GPR及TSP在隧道超前地质预报中的解译标志研究［D］.成都:成都理工大学，2009.(ZHAO Yu.Study on the Interpretation of GPR and TSPin Tunnel Geological Forecast［D］.Chengdu:Chengdu University of Technology，2009.(in Chinese))

［8］赵 峰，周 斌，武永胜.探地雷达在隧道衬砌空洞检测的正演模拟应用研究［J］.铁道建筑，2012，(8):99－103.(ZHAO Feng，ZHOU Bin，WU Yong-sheng.Study on Application of Forward modeling of Inspecting Void in Tunnel Lining by Ground Probing Radar［J］.Railway Engineering，2012，(8):99－ 103.(in Chinese))

［9］胡 俊，俞祁浩，游艳辉，等.探地雷达在多年冻土区正演模型研究及应用［J］.物探与化探，2012，36(3):457－460.(HU Jun，YU Qi-hao，YOU Yan-hui，et al.Study and Application of Ground Penetration Radar to the Forward Modeling of Permafrost Area［J］.Geophysical＆Geochemical Exploration，2012，36(3):457－460.(in Chinese))

［10］孙忠辉，刘金坤，张新平，等.基于GprMax的隧道衬砌地质雷达检测正演模拟与实测数据分析［J］.工程地球物理学报，2013，10(5):730－735.(SUN Zhonghui，LIU Jin-kun，ZHANG Xin-ping，et al.The Tunnel Lining Detection Forward Numeral Simulation and Measured Data Analysis Based on GprMax［J］.Chinese Journal of Engineering Geophysics，2013，10(5):730－735.(in Chinese))

［11］周奇才，李炳杰，郑宇轩，等.基于GPRMax2D的探地雷达图像正演模拟［J］.工程地球物理学报，2008，5(4):396－ 399.(ZHOU Qi-cai，LI Bing-jie，ZHENG Yu-xuan，et al.Forward Simulation of GPR Image Based on GPRMax2D［J］.Chinese Journal of Engineering Geophysics，2008，5(4):396－399.(in Chinese))

［12］YEEK S.Numerical Solution of Initial Boundary Value Problem Involving Maxell’s Equations in Isotropic Media［J］.IEEE Transactions on Antennas and Propagation，1966，14(3):302－307.

［13］喻振华，冯德山，戴前伟，等.复杂地电模型的探地雷达时域有限差分正演［J］.物探化探计算技术，2005，27(4):279－283.(YU Zhen-hua，FENG De-shan，DAI Qian-wei，et al.FDTD Method Forward Simulation of the Complex Geoelectricity GPR Model［J］.Computing Techniques for Geophysical and Geochemical Exploration，2005，27(4):279－283.(in Chinese))

［14］冯德山，戴前伟，何继善，等.探地雷达GPR正演模拟的时域有限差分实现［J］.地球物理学进展，2006，21(2):630－636.(FENG De-shan，DAI Qian-wei，HE Ji-shan，et al.Finite Difference Time Domain Method of GPR Forward Simulation［J］.Progress in Geophysics，2006，21(2):630－636.(in Chinese))

［15］何兵寿，张会星.地质雷达正演中的频散压制和吸收边界改进方法［J］.地质与勘探，2000，36(3):59－63.(HE Bing-shou，ZHANG Hui-xing.The Suppression of Numerical Dispersion and Improvement of Absorbing Boundary Conditions in Forward Modeling of GPR［J］.Geology and Prospecting，2000，36(3):59－ 63.(in Chinese))

［16］BERENGER J P.A Perfectly Match Layer for the Absorption of Electromagnetic Waves［J］.Journal of Computation Physics，1994，114(2):185－200.

［17］HASTINGS F D，SCHNEIDER J B，BROSCHAT S L.Application of the Perfectly Matched Layer(PML)Absorbing Boundary Condition to Elastic Wave Propagation［J］.The Journal of the Acoustical Society of America，1996，100(5):3061－3069.

［18］凌 飞，肖宏跃，朱夏乐，等.基于黏弹性介质的隧道地震二维正演模拟［J］.长江科学院院报，2015，32(5):121－ 126.(LING Fei，XIAO Hong-yue，ZHU Xia-le，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.(in Chinese))

［19］吕乔森，徐 颖，胡志强.提高隧道地震波地质预报(TSP)精度方法研究［J］.长江科学院院报，2012，29(12):41－ 45.(LV Qiao-sen，XU Ying，HU Zhiqiang.Methods of Improving Tunnel Seismic Prediction Accuracy［J］.Journal of Yangtze River Scientific Research Institute，2012，29(12):41－45.(in Chinese))