1. 引言
近年来随着城市基础设施建设快速发展,城区改造、道路开挖、管道改建等地下工程建设中,很容易导致各类管道破裂的现象,造成财产损失 [1]。城市建设施工过程中精确定位地下管线分布是避免工程事故的必要前提。探地雷达是一种地球物理检测方法,它是一种电磁法探测技术 [2],可用于埋藏于地下的各种管道的探测,如金属管线、水泥管线、PVC管线、地下电缆、输水管线、空洞、混凝土中的钢筋网,还可划分回填层及计算各层层厚等 [3] [4]。
常规的管线探测通常利用管线探测仪完成探测任务,但由于其只针对金属管线,且探测深度浅,在应用范围方面存在不足 [5]。而对于埋深较大的地下管线探测问题一直难以解决,主要原因在于:第一,埋藏深度过大,各种无损探测方法分辨率低难以捕捉到管线的特征;第二,对于分辨率较高的高频探测方法衰减快又难以达到4 m左右的深度;第三,管线埋设一般在场地狭小障碍较多的复杂地方,由于其探测范围狭小,不利于使用其它大规模的物探手段。探地雷达应用的范围极广,在城市地下管线探测中,不仅能够用来探测非金属管线,而且能够弥补管线仪探测深度浅等缺点 [6]。而且探地雷达能够利用屏蔽天线减少电磁波在地面上部反射造成的干扰;同时通过选择低频天线(16~100 MHz)增加有效探测范围(探测深度0~50 m);在后续处理中利用滤波等方法再压制外部的各种频率的电磁干扰信号,从而解决埋深较大的管线探测问题。在很多管线探测工程实例中,经常遇到埋藏大于2.5 m的目标体,且存在各种电磁干扰(如高压输电线的电晕、无线电发射设备、交流电网强大的电磁场以及各种电力设备等)。
本文结合北京华能热电厂深部地下管线探测工程,探测管线埋深大约为3~5 m。由于电厂内电磁信号干扰严重,且地表为绿化草地覆盖,给本次探测带来了很大的难题,因此可先通过数值模拟研究目标体反射波场特征,提高实测数据信号分析能力,便可增加处理结果的可靠性。
2. 探地雷达方法原理
2.1. 探地雷达国内外研究现状
电磁波能量在地下介质中会产生衰减,上世纪50年代Bailey将探地雷达应用与介质层相对单一的冰层探测。直至上世纪70年代,探地雷达应用于月球表面探测,该技术得到了更多的发展,使得应用范围更加广泛,如岩层结构探测、地下水位探测、工程地质探测等。在全球探地雷达快速发展的同时,国内针对该研究也得到了长足的发展。自90年代以来,通过自主研制和国外仪器的引进等方式,探地雷达在工程地质探测及科学研究方面得到了普遍的使用,如公路路基检测、地质灾害预测、矿产工程探测、隧道超前预报、考古及文物保护、隧道衬砌病害探测、地下管线探测等。
2.2. 工作原理
如图1所示,探地雷达(也叫地质雷达或透地雷达)法是一种浅层地球物理勘探方法,它是由激发天线T将电磁波以脉冲形式传入地下,经地下电磁差异界面反射后传回地面,被接收天线R所接收,是一种反射法的探测技术。电磁波旅行时为:
(1)
当地下介质中的介电常数已知时,便可确定波速v,可根据测得的时间t值,通过上式便可求得反射目标体的深度 [7] [8] [9]。
探地雷达是以自激自收方式记录反射波的波形,所以其图像为每一道的反射波的波形连续组合起来的剖面。波形的正负峰分别以不同颜色填充,主要以灰度图显示,这样在电磁差异界面会形成明显的同相轴用来表示目标体的某一截面。在地下管线探测中常沿着管线横向截面拖动天线进行探测,通过数据处理获得探地雷达图像,并据此判断地下管线的位置和大小 [10]。
Figure 1. Working principle diagram of GPR
图1. 探地雷达工作原理图
探地雷达在野外实际应用中它的探测深度主要取决于天线频率和地下介质的介电常数,介电常数的影响主要是因为介电常数越大对电磁波的吸收就越强。而天线频率与分辨率成正比,所以在能够达到探测深度时尽量选取高频天线。而它的探测效果与野外探地雷达工作时的参数选择,工区的现场环境干扰,以及后期的数据处理有关 [11]。
2.3. 正演模拟原理
正演模拟采用时域有限差分法,该方法将麦克斯韦方程组进行差分离散,从而在时域进行电磁波的模拟,本次模拟100 MHz天线对大埋深为5 m,半径为1 m的管线进行正演,分析其不同排列方式,不同材质的管线反射波特征规律,从而提高实测数据信号分析能力,对工区实测数据中的管线反射波特征起到指导作用。
2.3.1. 时域有限差分FTDT理论
麦克斯韦方程组统一了电场和磁场的基本规律,微分形式包含法拉第电磁感应定律和全电流定律两个旋度方程,以及静电场高斯定律、磁通连续性原理两个散度方程 [12]。Yee提出的时域有限差分法,是基于麦克斯韦旋度方程进行差分离散。麦克斯韦旋度方程 [13] [14]:
(2)
E——电场强度,单位V/m;H——磁场强度,单位A/m;J——电流密度,单位A/m2;反映电损耗。D——电通量密度,单位C/m;B——磁通密度,单位Wb/m;Jm——磁流密度,单位V/m2;反映磁耗。
麦克斯韦方程组主要研究地下介质的磁场和电场特性。本构关系描述了物质的基本粒子在外部电磁场作用下的整体宏观特性。在各向同性线性介质中的本构关系为:
——磁导率,单位H/m;真空中 = 4 × 10−7 H/m;
——电导率,单位S/m;真空中 = 0。
——导磁率,单位Ω/m;真空中 = 0;
——介电常数,单位F/m;真空中 = 8.85 × 10−12 F/m
在直角坐标系中,以E, H为变量展开麦克斯韦旋度方程分别为:
(3)
令
代表E, H在直角坐标系中的某个分量,时间和空间域中离散符号取为:
(4)
Δx——x方向空间步长;Δt——时间步长;j——y方向空间步长个数;Δy——y方向空间步长;n——时间步长个数;k——z方向空间步长个数;Δy——y方向空间步长;i——x方向空间步长个数。
对
关于的时间t和空间x,y,z方向的一阶偏导数取中心差分近似为:
(5)
模拟地质雷达的电磁波在地下介质层中传播时,Yee提出一种计算电磁波数值的时域有限差分法,根据麦克斯韦方程组的差分近似形式,通过计算每个单位网格中的电场和磁场分量来计算获得整个仿真空间的电磁场值 [18],如图2所示。
Yee元胞中E分量与H分量在空间上都是垂直分布。E与H相互环绕。相邻E(H)相隔一个空间(时间)步长,相邻E, H相隔半个空间(时间)步长。面上的磁场分量近似相等,Yee元胞棱边的电场分量近似相等。x,y,z方向上的时间步长是相等的。
2.3.2. 激励源
数值模拟采用的激励源为100 MHz雷克子波,波形及其频谱如图3所示。雷克子波形状简单、零相位、只有一个正峰、两侧各有一个旁瓣、延续时间很短、收敛快。适用于雷达正演模型 [16]。
2.3.3. 数值色散
时域有限差分法对麦克斯韦方程进行数值计算时,会引起色散,即速度将随频率改变从而引起脉冲波型畸变等一些虚假现象。当离散间隔∆x, ∆y, ∆z满足以下条件能够减小这种色散 [17]:
是介质中雷达波最下波长值,空间离散长度取满足上式的最小值。
2.3.4. 解的稳定性
FTDT方法是以差分方程的解代替麦克斯韦旋度方程组的解。在时域有限差分法中空间步长(∆x, ∆y, ∆z)和时间步长∆t满足一定的关系条件,离散后的解是收敛且稳定的 [18]。
时间步长稳定性条件:
二维空间时间步长与空间步长稳定性条件:
三维空间时间步长与空间步长稳定性条件:
2.3.5. 吸收边界条件
PML吸收边界条件,是在模拟地下地质体时,由于地下实际为无限大空间电磁波会无限传波,但仿真的地下空间是有限的。当电磁波传播到边界面时会发生波的反射,为了模拟无限大空间我们要将传播到边界面的电磁波衰减掉。PML边界能够大大的吸收电磁波且不发生反射 [19]。
2.4. 二维地下管线正演模拟
地下管线探测中常沿着管线横向截面拖动天线进行探测,在对地下管线截面模拟中(图4),雷达正演
模拟的天线中心频率为100 MHz,周围土壤介质的介电常数为9,待模拟的管线截面为混凝土材质假定介电常数为5。管线半径为1 m,坐标(5 m, 5 m)处。其合成的雷达剖面图如图5所示,图中可见圆在雷达正演模拟中呈双曲线,深度4 m处为弧形顶部,且为两层平行弧线。上层相位呈现蓝黄蓝,下层为黄蓝黄,这是由于电磁波从高介电常数介质进入低介电常数介质时反射系数为正,而从低介电常数介质到高介电常数中反射系数为负,相位发生反转所致 [20]。电磁波在介质中传播的波场快照如图6所示,当电磁波在传播过程中遇到地下管线时,一部分电磁波发生波的反射,一部分透射过管线在管线下界面亦发生反射和透射。
30 ns 70 ns 80 ns
Figure 6. Wave field propagation snapshot
图6. 波场传播快照
2.5. 多排管线模拟
图7为一排管线截面排布,模拟的天线频率为100 MHz。周围土壤介质的介电常数为9,模拟3个相同介电常数均为5的并排管线,管线半径为1 m,管线埋深5 m,管线间隔6 m。合成的雷达剖面图如图8所示。从B-scan图中我们可以看出介电常数不变的并排管线反射波特征表现为相位一致,管线顶底界面的反射波时差为40 ns,计算可得管线直径约为2 m与正演模型一致。
图9为一排管线截面排布,模拟的天线频率为100 MHz。周围土壤介质的介电常数为9,模拟3个并排管线介电常数依次为(5, 9, 20),管线半径为1 m,管线埋深5 m,管线间隔6 m。合成的雷达剖面图如图10所示。从图中我们可以看出随着反射系数绝对值依次增大,反射波越明显,管线下界面的反射波位置随着介电常数变大随之加深,这是因为电磁波在不同介质中的传播速度不同,介电常数越大传播速度越慢所致。
第一根管线模拟介电常数为5小于周围介质,电磁波由高介电介质传播到低介电介质相位不变;第二根管线介电常数与周围介质一致,无明显电磁波反射特征;第三根管线介电常数为20高于周围介质,由图10可以看出:反射波相位与第一根管线相位相反、发生反转。因此根据反射波相位特征判断地下存在不同材质的管线及分布。
图11为纵向排布的圆柱状管线截面,模拟的天线频率为100 MHz。周围土壤介质的介电常数为9,管线截面介质介电常数5。管线半径为1 m,其横坐标都在10 m处,纵坐标以此相隔3 m纵向排开。其合成的雷达剖面图如图12所示,从图中我们可以看出随着管线埋藏深度依次增大,反射波越弱。所以当地下埋藏多层管线时,由于上部管线会反射大部分的能量,透射后再反射能量越来越小,即其同相轴越来越模糊、开口角度越来越大。
图13为纵向排布的圆柱状管线截面,模拟的天线频率为100 MHz。周围土壤介质的介电常数为9,
管线截面介质介电常数从上到下依次为(5, 20)。管线半径为1 m,其横坐标都在10 m处,纵坐标以此相隔3 m纵向排开。其合成的雷达剖面图如图14所示,从图中我们可以看出随着管线埋藏深度依次增大,反射波越弱,即其同相轴越来越模糊、开口角度越来越大。上下管线反射波相位相反,所以也可以根据反射波相位特征判断地下存在不同材质的管线及分布。
3. 工程实例
以下通过对北京市朝阳区华能北京热电厂地下管线探测工程实例阐述探地雷达在地下深埋管线探测方面的应用。根据搜集、现场勘察的本场区地质资料可知所存在的难点问题在于,工区位于厂区内部的草坪之上对电磁波具有一定的吸收作用,外界电磁干扰严重主要干扰源有(高压输电线、无线电发射设备、交流电网、路灯杆以及各种电力设备等),管线的埋藏深度大约在3~5 m左右。综合场地地质条件及表1中物性参数可见,场区内第四系粘土层、砂砾质粘土层与砂岩、花岗岩之间均存在较明显的电性差异。如果场区内存在有地下管线,由于水泥,金属物性参数与第四系地层之间存在着更为明显的差异,这为用反映电性差异的地质雷达法来探查地下管线提供了有利条件 [21]。
Table 1. Physical parameters of some geotechnical media
表1. 部分岩土介质的物性参数
地质雷达测线布置主要在地下管线附近及道路共布设16条测线,在图15中用绿色标注线号。现场用皮尺依据标志物的位置对测线进行测量,并对物探测线控制点进行标注,测线尽量远离外界干扰源。如遇障碍物无法穿越测线位置时发生偏移的在原始记录表中做好记录并在图15平面布置图上进行说明。
地质雷达法是依据地下介质界面对电磁波的反射同相轴的连续性进行判断,数据利用美国的Radan7软件进行数据处理成图 [22]。本次工作中所选工作参数为:天线主频:100 MHz (屏蔽天线);时窗:140 ns、320 ns;采样点数:512;测量方式:距离模式;叠加次数:自动。
探地雷达法资料处理的目的是压制随机和规则的干扰波,提高信噪比与分辨率,在探测剖面上易于识别目标体。地质雷达法是依据地下介质界面对雷达波的反射同相轴的连续性进行判断,采用美国Radan7软件进行数据处理成图。其处理流程如下:数据编辑→距离归一化→调整地面时间零点→时间深度转换→增益调整→滤波→偏移→输出时间剖面图→根据速度分析进行时深转换。
图16为100 MHz天线探测数据中的一段(005测线),从黑框处可以很明显地发现在测线里程2 m,埋深4 m处有雷达波同相轴出现与正演模型一致的弧形现象,此类反射信号为地下管线的雷达响应图像。而在里程7 m,深3 m处的范围很广的弧形特征为一路灯杆的干扰特征 [23] [24]。
通过对采集的雷达数据后处理再结合地质资料,获得探测区地下管线分布(图17),图中共探测了4条主管线,主要呈东西走向L型分布。东西向四条主要管线基本与东西向主路平行,其方向在西侧向北延伸。探测详细结果如图所示,本次探测在后期开挖中得到了很好的验证。
Figure 17. Detection result of Beijing Huaneng thermal power plant
图17. 华能北京热电厂地下管线探测结果
4. 结论
通过数值模拟对电磁波反射信号特征识别与探地雷达在强干扰源背景下探测深部管线,得到如下结论:
1) 低频探地雷达探测深度大,能够探测到4 m埋深的管线,其效果优于管线探测仪;
2) 在存在外界强干扰的电厂区内,通过对地面上的电磁波屏蔽以及干扰波主频的滤除,能有效定位目标管线的位置;
3) 同一管线随着介电常数依次增大,反射波越明显,但下界面的反射波位置依次加深。这是由于电磁波在不同介质中的传播速度不同,介电常数越大传播速度越慢所致,通过反射波相位特征可判断管线材质;
4) 当地下埋藏多层管线时,由于上部管线会反射大部分的能量,透射后再反射能量越来越小,即其同相轴越来越模糊、开口角度越来越大。
基金项目
国家自然科学基金(41641040);国家自然科学基金(41790443)。