基于GPRMAX的隧道超前地质预报正演模拟与实测数据分析

朱云峰，王齐仁，邓国文

(湖南科技大学　土木工程学院，湘潭　411201)



基于GPRMAX的隧道超前地质预报正演模拟与实测数据分析

朱云峰，王齐仁，邓国文

(湖南科技大学土木工程学院，湘潭411201)

摘要：以电磁波的传播理论为依据，总结了基于GPRMAX2D的地质雷达正演模拟生成的二维剖面图的一般规律，探讨了隧道掘进过程中遇到的典型不良地质体的雷达图像特征，依据正演模拟结果对实测地质雷达数据进行了相应地解译，并用开挖实例验证了解译成果的准确性。

关键词：隧道超前预报；GPRMAX；正演模拟；数据分析

0引言

为确保在隧道掘进过程中人员、设备以及工程的安全，需要查明掌子面前方围岩的地质属性、溶洞、断层破碎带及富水区等一系列不良地质体。虽然在隧道开挖前工程人员就已对隧道所经地区的地质情况有所了解，但是地质勘探所提供的资料往往只在大区域内较准确，由于条件限制并不能对隧道施工过程中所遇到的局部地质构造做出详细地描述[2]，因此，有必要在施工过程中进行超前地质预报工作，以帮助确定围岩级别、预测地质灾害和制定施工方案。探地雷达(GPR)作为一种无损的检测方式[1-4]以其经济、灵活、轻便等优点早已被用于隧道超前地质预报。但该方法是一种间接探测方法，数据解译依赖解释人员的经验，因此，解释成果的准确性也因人而异。由于受到各种条件的制约(如天气、场地、施工时间、现场的机械设备等)，想通过大量的试验来提高雷达数据解释的准确性是不太现实的，而利用计算机对近乎现场的模型做正演模拟则可以解决这一矛盾[5]。

1GPRMAX软件及其理论基础

GPRMAX2D是Dr Antoni Giannopoulos于1996推出的一种基于FDTD(Finite-difference Time-domain)算法和PML(PerfectlyMatched Layer)边界吸收条件的探地雷达正演模拟软件，用于研究探地雷达的成像规律[2-3]。经过十几年的发展，目前GPRMAX不仅可以模拟二维模型还可以模拟三维模型，使目标体的波形剖面图更加直观[6]。

1.1FDTD算法

FDTD算法是K.S.Yee于1996年提出来的。该方法将Maxwell旋度方程组(公式(1))化为一组电场和磁场各分量的偏微分方程，然后利用二阶精度的中心差近似将这一组偏微分算符转换为差分形式，通过解差分方程得到微分方程解的近似值[7]。对不同的介质模型进行电磁波差分法数值模拟，就可以研究其相应的探地雷达时间—深度剖面图。

(1)

1.2PML边界吸收条件

受计算机容量的限制，在用GPRMAX进行模拟时FDTD计算只能在有限的区域进行。因此，为了在有限的地下介质计算区域模拟真实的无限空间的电磁问题，必须在计算区域的截断边界上设置吸收边界条件[8]。

完全匹配层(Perfectly Matched Layer,PML)是由J.P.Berenger首先提出，它在计算区域边界面附近引进虚拟的各向异性有耗媒质，模拟区域内的外行电磁波可以无反射地进入有耗媒质，并在有耗媒质内进行衰减，从而有效吸收模拟区域内出射的外行波。其边界吸收条件公式及推导过程请参考文献[11]。

2岩石的介电常数与反射系数

各种岩石矿物的介电常数都大于“1”，而水的介电常数达81，比各种常见的造岩矿物都大的多。所以介质中水的含量是影响其介电常数的的主要因素之一[9]。而电磁波在岩层中传播遇到介电常数不同的分界面会产生反射和透射，其反射系数如式(2)。

(2)

其中：γ为反射系数；ε为介电常数。

反射系数与界面两侧介质的介电常数的关系如图1所示。电磁波遇到界面发生正反射还是负反射取决于介质的介电性，当岩性不变或变化不大时，含水性的变化就决定了介电常数的不同。因此，可根据电磁波的正负反射来判断目标体是否含水，可以对有地下水引起的地质灾害进行准确预报[9]。

图1　电磁波的正负反射Fig.1　The positive and negative reflection of electromagnetic waves

当ε1<ε2时，γ<0电磁波从高阻介质过渡到低阻介质产生正反射，即反射波与入射波相位相反；当ε1>ε2时，γ>0电磁波从低阻介质过渡到高阻介质产生负反射 ，即反射波和入射波相位相同。图1为四层介质电磁波的正负反射示意图，各层介质的介电常数关系为：ε1<ε2>ε3<ε4。

3典型不良地质体模型及正演模拟

针对在隧道掘进过程中常遇到的一些典型的不良地质体，对以下三种异常做正演模拟：①断层破碎带(填充物为碎石、水)；②溶洞(无填充物)；③富水区域(填充物为水)。建立的模型尺寸都为1.5 m×1.05 m，模拟网格步长Δx=Δy=0.002 5 m，模拟天线采用900 MHz，天线间距为0.005 m，天线步长为0.01 m，电磁波类型为雷克子波(ricker)，模型中各种介质所需要的物理参数见表1。

表1　模型中所需介质的物理参数Tab.1 The physical parameters of the medium in the model

图2为完整围岩模型及正演模拟图，从图2(b)可以看出，电磁波只有在掌子面才有较强烈的反射，同相轴较连续也可以反推出掌子面较平整。因为围岩的介电常数未发生变化因此掌子面前方几乎没有电磁波发生反射。

图3为断层破碎带模型，图3中大矩形空白部分的填充物为水，圆、三角形、小矩形为碎石和淤泥。从图3(b)可以看出,断层的电磁波反射呈现出杂乱无章，这是因为其介质的分布是没有规律的。分析图3(c)可以看出，电磁波在到达大矩形(填充物为水)时，发生了正反射。这是因为岩石的介电常数比水的介电常数小的多，此时的反射系数γ<0，因此，由前面的理论部分也可知电磁波在从岩石过渡到水中时会发生正反射。

图4为溶洞模型及正演模拟图，其中溶洞设置为近圆形，填充物为空气。在GPRMax中空气(free_space)的参数是默认的，不需要特意设定，因此，在写代码时直接调用free_space即可[10]。从图4(b)中可以看出，图形整体呈现出多次反射的双曲线形，电磁波在遇到溶洞顶部时反射最为强烈，在下方18 ns以后还有一次电磁波反射较为强烈，通过电磁波在空气中传播的速度与其双程走时计算得出的距离与模型中溶洞的直径对比可以得出，其应该是电磁波在遇到溶洞底部时的反射。从图4(c)中可看出，电磁波在由岩石过渡到空气中时会发生负反射，这是因为岩石的介电常数大于空气的介电常数。

图5为富水区(矩形)模型与正演模拟图，模型中矩形的填充物全是水。从图5(b)中可以看出，矩形物体的波形与圆形物体的波形有明显的区别，圆形物体的波形在顶端是拱形的，而矩形物体的波形顶端较平并逐渐向左右两边呈曲线延伸。分析图5(c)时可以发现其与图3有着相似的规律，电磁波在分界面也都发生了正反射。为了进一步验证这一规律的正确性，下面的模型运用控制变量法将其与图4进行比较。

图2　完整围岩模型及正演模拟图Fig.2　A complete model of surrounding rock and the forward simulation diagram(a)模型；(b)堆积波形；(c)单道波形

图3　断层破碎带模型及正演模拟图Fig.3　Fault fracture zone model and forward modeling(a)模型；(b)堆积波形；(c)单道波形

图4　溶洞模型及正演模拟图Fig.4　Cave models and forward modeling(a)模型；(b)堆积波形；(c)单道波形

图6的模型是围岩中存在一个富水的圆形溶洞，其基本位置、大小与图4中的溶洞(填充物为空气)是一致的。从图6(b)与图4(b)对比可以看出，含水的溶洞波形图没有含空气的溶洞波形图连贯，并且图6中电磁波在图的底部有强烈反射，两层反射之间的距离比图4大，这是因为电磁波在水中的传播速度比在空气传播慢的缘故。分析单道波形图可知电磁波在遇到溶洞顶部时发生了正反射，这与图4中单道波形图的负反射是相反的。这一结论可以为隧道中潜在的由地下水引起的地质灾害，提供可靠地预判。

图5　富水模型及正演模拟图(矩形)Fig.5　Rich water model and forward modeling diagram(rectangular)(a)模型；(b)堆积波形；(c)单道波形

图6　富水模型及正演模拟图(圆形)Fig.6　Rich water model and forward modeling diagram(round)(a)模型；(b)堆积波形；(c)单道波形

4实测数据分析

数据是在福建省南龙铁路Ⅱ标段某隧道内采集的，该段隧道分为正洞和斜井两部分，斜井与正洞交接处的桩号为DK30+348。地质雷达使用的是中国电波传播研究所的LTD2100，天线采用的是100 MHz。为了更好地进行对比，实测数据在后期处理过程中增益指数都设定为5.5。探测数据如图7所示，其中，a剖面为正洞中的数据，b、c、d剖面为斜井中的数据。

图7中a剖面的掌子面桩号为DK30+642。该段围岩条件较好，从图7(a)中可以看出，该次预测范围内的波形图几乎没有强烈的电磁波反射，围岩的局部范围可能存在空隙和碎石，另外，掌子面的同相轴较连续说明掌子面是比较完整的，其受爆破的扰动不是很大，围岩级别应该为Ⅱ级。

图7(b)剖面掌子面的桩号为X2DK0+605，从图7(b)中我们可以看到一条明显的断层或破碎带，在破碎带的下方离掌子面大概9.0 m左右有段弧形的反射波，雷达电磁波在此反映较强烈，根据正演模拟结果推断其应该是个溶洞，洞顶较明显但洞底不是很明显。

图7中c剖面掌子面的桩号为X2DK0+435，可以看出，在距离掌子面大概6 m位置有一条从右向左，从外向内的电磁波反射信号带。因为这一信号并非像图7(b)那样多次反射，因此,推断这可能是条裂隙。在12 m左右有一条明显的弧形反射信号，在其下方3 m左右又出现一条弧形反射信号，根据正演模拟结果推测，其应该是一个近乎圆形的溶洞，上方的弧形为洞顶，下方的弧形为洞底，直径大概在3 m。

图7(d)剖面掌子面的桩号为X2DK0+256.4。从图7(d)中可以看到，两片杂乱无章的强反射区域，将其放大并仔细观察可以看出，许多波形图首先发生的是正反射。参考正演模拟得到的结果，推测该区域围岩较破碎，裂隙发育，含水丰富，为富水破碎带。

将分析结果与之后挖掘到该位置的实际情况进行对比(图8)，揭示的实测数据分析与真实地质异常体基本吻合。这说明基于GPRMAX的隧道超前预报正演模拟基础上做的剖面异常推断,是较准确可行的。

图7　雷达实测堆积波形图Fig.7　Radar measured waveform figure(a)完整围岩；(b)断层破碎带；(c)溶洞；(d)局部富水

图8　开挖结果验证图Fig.8　Certified figure of excavation results(a)隧道中裂隙；(b)隧道中溶洞；(c)隧道中富水区

5结论

1) GPRMAX2D结合Matlab语言用于隧道超前预报的正演模拟是可行的，模拟得出的波形图是比较准确的，这对提高解释人员的解释精度有一定的帮助。并且模拟程序较简单适合初学地质雷达超前预报正演模拟的人员使用。

2) 从单道波形图中电磁波在分界面发生的是正反射还是负反射可以判断地质异常体介电常数的变化规律，若发生正反射则该区域含水量可能增多。

3) 用GPRMAX2D进行隧道中典型地质异常体的超前预报正演模拟所生成的波形图存在一定的相似性，因此，需要结合加深炮孔、水平钻探等其他超前预报手段以提高隧道超前预报的准确性。

参考文献：

[1]曾昭发，刘四新，王者江，等.探地雷达方法原理及应用[M].北京：科学出版社，2006.

ZENG Z F,LIU S X,WANG Z J,et al.The principle and application of GPR method[M].Beijing:Science press,2006.(In Chinese)

[2]周奇才，李炳杰，郑宇轩，等.基于GPRMAX2D的探地雷达正演模拟[J].工程地球物理学报,2008,5(4)：396-399.

ZHOU Q C,LI B J,ZHENG Y X，et al.Forward simulation of GPR image based on GPRMax2D[J].Journal of engineering geophysics,2008,5(4):396-399.(In Chinese)

[3]丁亮,陈璐，孙烨，等.基于GPRMAX2D的混凝土无损检测方法[J].黑龙江大学自然科学学报,2013,30(6)：819-825.

DING L，CHEN L，SUN Y，et al.Concrete nondestructive testing methods on GPRMAX2D[J].Heilongjiang university journal of natural science,2013,30(6):819-825.(In Chinese)

[4]李敏瑞，马景辉.GPRMAX2D在巷道质量检测中的应用[J].煤矿安全,2013,44(9)：166-167.

LI R M,MA J H.The application of GPRMAX2D in the Quality of Tunnels[J].Coal mine safety,2013,44(9):166-167.(In Chinese)

[5]孙忠辉，刘金坤，张新平，等.基于GPRMAX2D的隧道衬砌地质雷达检测正演模拟与实测数据分析[J].工程地球物理学报,2013,10(5)：730-735.(In Chinese)

SUN Z H,LIU J K,ZHANG X P，et al.The Tunnel Lining Detection Forward Numeral Simulation and Measured Data Analysis Based on GPRMAX[J].Journal of engineering geophysics,2013,10(5):730-735.(In Chinese)

[6]宋申宇，于会山.基于GPRMAX的探地雷达图像正演模拟[J].科技信息,2010,10(7)：35-36.

SONG S Y,YU H S.Forward Simulation of GPR Image Based on GprMax[J].Science and technology information,2010,10(7):35-36.(In Chinese)

[7]郭立，崔喜红，陈晋.基于正演模拟的探地雷达根系探测敏感因素分析[J].地球物理学进展,2012,27(4)：1755-1760.

GUO L,CUI X H,CHEN J.Sensitive factors analysis in using GPR for detecting plant roots based on forward modeling[J].The earth physics are reviewed,2012,27(4):1755-1760.(In Chinese)

[8]杨天春，李好，张辉.探地雷达数值模拟中广义完全匹配层吸收边界条件的模拟分析[J].计算物理,2014,31(2)：200-201.

YANG T C,LI H,ZHANG H.The Simulation analysis of the absorbing boundary’s condition of the generalized PML On GPR numerical simulation[J].Computational physics,2014,31(2):200-201.(In Chinese)

[9]王齐仁.隧道地质灾害超前探测方法研究[D].长沙:中南大学,2008.

WANG Q R.Tunnel geological hazard research advanced detection methods[D].Changsha:Central south university,2008.(In Chinese)

[10]Dr Antonis Giannopoulos GprMax2D/3D User’s Manual[EB/OL].[2015].http://download.csdn.net/detail/xajh 20102010/5109754.

[11]葛德彪，阎玉波.电磁场时域有限差分法[M].西安：西安电子科技大学出版社，2002.

GE D B,YAN Y B.Electromagnetic field finite difference time domain method[M].Xi’an:Xian university of electronic science and technology press,2002.(In Chinese)

The forward modeling and measured data analysis of the tunnel advanced prediction based on GPRMAX

ZHU Yun-feng，WANG Qi-ren，DENG Guo-wen

(Hunan university of science and technology institute of civil engineering，xiangtan411201，China)

Abstract:Based on theory of propagation of electromagnetic waves,we summed the universal law of a two-dimensional section of geological radar which is generated by the forward modeling based on GPRMAX.And the typical unfavorable geologic bodies are discussed during the process of tunnel excavation.We make the corresponding judgment about the measured geological radar data with reference to the results of the forward modeling,and check the result of the judgment by a living example.

Key words:the tunnel advanced prediction;GPRMAX;forward modeling;the data analysis

收稿日期：2015-01-30改回日期：2015-03-31

基金项目:湖南科技大学研究生创新基金(S130009)

作者简介：朱云峰(1989—)，男，硕士，主要从事工程物探，E-mail:feng_zy1989@163.com。

文章编号：1001-1749(2016)02-0185-06

中图分类号：P 631.3

文献标志码：A

DOI：10.3969/j.issn.1001-1749.2016.02.07