基于GprMax的地下管线探地雷达图像正演模拟

2022-04-14莫其妙

工程地球物理学报 2022年2期

张迪,莫其妙

(1.南宁市勘测设计院集团有限公司检测中心，广西南宁 530000；2.南宁市勘测设计院集团有限公司规划咨询分院，广西南宁 530000)

1 引言

地下管线属于城市基础设施，是城市赖以生存和发展的物质基础[1]，具有隐蔽性。以南宁市为例,南宁市在进行排水管网普查过程中发现存在部分地下排水管线淤堵、设计图缺失、检查井被覆盖等情况，导致该区域部分管网连接关系不明。若该区域排水管线内部存在错、混接情况，可能会污染城市内湖及周边水域；在该区域进行城市建设很可能会挖断下方不明管线，给施工人员与附近居民带来安全隐患，故查明地下管线意义重大。

近些年，高密度电法[2]、地震映像法[3]、瞬变电磁法[4]等物探方法越来越多被应用在城市勘察中，但多数情况下在沥青、混凝土路面及车流量大的狭窄场地不适用[5]，探地雷达法[6]作为无损探测的新方法为此类不具备开挖、钎探验证区域内的探测提供了技术支持。如陈海涛等[7]使用探地雷达在桥梁预应力管道探测方面取得了良好的效果；吴宝杰等[8]利用探地雷达成功找到了地下废弃化粪池位置；许泽善等[9]、赵明堂[10]使用探地雷达在道路地下病害体方面获得了较好的效果；张汉春等[11]、陈军等[12]、赵永辉等[13]研究了地下管线雷达图像特征，并得到了理想的效果。

通过正演模拟可以高效地提高雷达图像解译经验，为今后辨别雷达时间剖面图提供数据基础。雷勤梅等[14]、张军伟等[15]利用GprMax正演模拟软件进行了探地雷达管线正演模拟；周奇才等[16]利用GprMax正演模拟软件进行了地铁隧道和桩方面的探地雷达正演模拟。本文利用GprMax正演模拟软件结合Python语言，对典型地下管线进行正演模拟并对比验证，最后总结经验，为以后的雷达探测给予指导。

2 探地雷达及正演基本原理

2.1 探地雷达工作原理

探地雷达简称GPR(Ground Penetrating Radar,GPR)，由雷达主机、电子单元、发射和接收天线组成，是利用25～2 300 MHz的无线电磁波对地下结构、埋藏物成像的探测技术。图1为探地雷达的工作原理，工作时发射天线(Tx)向地下发射电磁波，遇到有电性差异的异常体时会产生可被接收天线(Rx)接收的反射波，同时雷达主机记录反射波到达时间及相关信息，进而组成探地雷达时间剖面图。

图1 雷达工作原理Fig.1 GPR schematic diagram

图2为模拟500 MHz天线探测地下埋深0.3 m、尺寸0.2 m管线的波场快照。探地雷达未工作时，雷达发射天线不发射电磁波，如图2(a);探测时探地雷达发射天线向四周发射电磁波，见图2(b)；电磁波探测到地下管线时发生反射，见图2(c)；再由接收天线接收，部分电磁波继续向下传播；随着探测深度与时间推移，电磁波能量不断衰减，见图2(d)。

图2 波场快照Fig.2 Wave field snapshot

不同介质中电磁波传播速度不同，传播速度与介质相对介电常数和电导率有关，传播速度由式(1)表示。由式(2)可得异常体的埋深，带入可得异常体的时深转换公式(3)：

(3)

其中，Vr为电磁波在介质中的速度(m/ns);Hr为在介质中的传播距离(m);εr为介质的相对介电常数;t为电磁波双程走时(ns);μr为介质的相对磁导率，默认为1，c为电磁波在空气中的传播速度(m/ns)，文中取0.3 m/ns。

2.2 分辨率

探地雷达分辨率分为横向分辨率与垂向分辨率，指在横向与垂向方向上将两个近距离物体辨别开的能力。

横向分辨率通常由式(4)表示：

(4)

垂向分辨率通常由式(5)表示：

(5)

其中，Δh代表横向分辨率(m);Δd代表垂向分辨率(m);λ代表波长(m);d代表目标埋深(m);fc代表天线中心频率(MHz)。可见，一定条件下探地雷达分辨率与天线的中心频率成正相关。

2.3 反射系数

电磁波穿过有相对介电常数差异的介质时会产生反射和透射，两者能量分配主要与电磁波反射系数有关[17]，反射系数可由式(6)表示：

(6)

其中，r表示反射系数;ε1表示第一层介质的相对介电常数;ε2表示第二层介质的相对介电常数。当上层介质相对介电常数大于下层介质相对介电常数时，反射数为正，反射波相位为正，反之为负。这也是工作人员对地下异常体材质、填充物辨别的一种手段[18]。

3 GprMax正演模拟

探地雷达正演模拟是指在计算机上利用模拟软件模拟电磁波在地下介质中的传播过程，通过研究分析在遇到异常体时反射波的相关信息及波场分布特征，进而与实际中探地雷达时间剖面图进行对比验证，提高对探地雷达结果的解释精度。GprMax正演模拟工具由英国爱丁堡大学的Dr Antonis Giannopoulos[19]研发，可用于模拟电磁波传播。通过电磁波与目标物体的相互影响，得到目标物体的探地雷达地质图像[20]。该正演模拟工具基于时域有限差分法(Finite-Difference Time-Domain，简称FDTD)和PML(Perfectly Matched Layer,PML)边界吸收条件，以Yee元胞[21]为空间电磁场离散单元，将Maxwell旋度方程转化为差分方程，进而通过求解差分方程得到微分方程解的近似值[22]。因在时间上迭代求解，有较好的稳定性和收敛性[23]。

在FDTD中，时间增量及空间步长不是相互独立的，为了保证模拟过程中电磁波的关系不被破坏，模拟数值解要稳定，所以时间增量及空间步长需满足稳定性条件(Courant-Freidrichs-Lewy ，简称CFL)：

(7)

其中，Δx、Δy、Δz为空间步长(m);Δt为时间步长(ns);c为电磁波在空气中的传播速度(m/ns)。

3.1 地下管线的二维正演模拟

城市地下管线一般包括明挖与机械顶管两种方式敷设，埋深多在0.3～3.0 m范围内，管径多为0.05～2.5 m，材质有PVC(Polyvinyl Chloride)、砼及铸铁等。为使模拟结果准确且贴近实际，文中使用500 MHz雷达天线，设置四层统一道路结构：首层为空气层，层厚0.15 m；第二层为面层，成分为沥青，层厚0.15 m；第三层为基层，主要成分为混凝土，层厚0.3 m；第四层为路基，主要成为为黏土、砂砾等。正演模型统一为3.5 m×2.5 m、空间网格步长为0.002 5 m×0.002 5 m、时窗为50 ns，使用雷克子波作为激励源，测线为160道的剖面法正演模拟。相关介质物性参数见表1。

表1 模型中介质物性参数Table 1 Physical parameters of model intermediates

3.1.1 层状路面正演模拟

图3为层状路面几何模型与正演模拟图。图3(a)为几何模型图，从上往下依次为蓝色空气层、黑色路面层、灰色基层及红褐色路基层；图3(b)为相应正演模拟图(B-Scan)，图3(c)为第83道单道模拟图(A-Scan)，雷达天线频率为500 MHz。

图3 层状路面层状路面几何模型与正演模拟Fig.3 Geometrical model and forward simulation modeling

由图3(b)、图3(c)知，在2.5 ns、4.0 ns、5.8 ns、11.5 ns时场强发生了较大变化，表明此处物质相对介电常数和电导率差异大。顶部反射最强的为直达波，往下依次为面层、基层与路基层反射波。因路面层较薄，雷达子波又有一定的宽度和衰减震相，导致路面反射波与直达波出现了相互叠加现象，基层与路基层之间因相对介电常数差异出现了能量较弱的反射波。根据公式(3)可得面层厚0.135 m，基层厚0.285 m，与模型参数设置的数值基本吻合。

3.1.2 不同形状地下管线正演模拟

图4为探地雷达对圆形管道与矩形暗涵的响应特征，设置圆形混凝土管线：管顶距路面0.9 m，半径0.25 m，壁厚0.02 m，内部充气；矩形混凝土暗涵：管顶距路面0.9 m，规格0.6 m×0.5 m，壁厚0.02 m，内部充气。图4(b)、图4(c)分别为圆形管道与矩形暗涵的B-Scan，图4(d)、图4(e)分别为其第83道的A-Scan,雷达天线频率为500 MHz。

图4 圆形管道与矩形暗涵几何模型与正演模拟Fig.4 Geometrical model and forward simulation modeling

19 ns时两者B-Scan均出现开口向下的管顶双曲线型反射波，反射波能量强，两侧绕射波弱。矩形暗涵上顶板在水平方向有一定延伸，对应矩形暗涵的上顶板。因暗涵两侧壁有一定角度，故反射波相较圆形管道较杂乱，两者均有顶底反射及多次波出现。涵洞两侧绕射波之间的距离即为涵洞的宽，“X”型反射波组是矩形暗涵的代表剖面图。图4中涵洞左翼对应第66道的位置，右翼对应第95道的位置，可推断暗涵宽度为0.58 m，与模型中涵洞宽0.6 m接近。根据探地雷达工作原理，可由顶底反射波推算异常体尺寸。由圆形管道及矩形暗涵A-Scan，两者均在20.5 ns及23.5 ns处出现明显顶底反射波，其后均伴有多次波。顶底反射波时差为3 ns，得出传播距离为0.45 m，加上0.02 m壁厚，得圆形管道直径、暗涵高为0.47 m，与模型尺寸基本接近。圆形管道与矩形暗涵顶反射与路基层反射双程走时差为10.5 ns，电磁波在路基层速度为0.086 m/ns，可得电磁波传播距离为0.476 m，得管顶距路面层为0.876 m，与设置管顶埋深0.9 m接近。

3.1.3 不同埋深的管线正演模拟

高频雷达天线分辨率高，但在介质中衰减快，多进行浅部、小尺度目标体探测，如沥青层厚度、钢筋混凝土结构探测等。低频雷达天线分辨率低，探测深度大，多用来进行城市道路路基病害、隧道超前地质预报及地下大管径管线探测。为研究不同埋深目标体对低频和高频天线的响应，现进行如下模拟：选用250 MHz、500 MHz及900 MHz天线分别对埋深为0.65 m、1.05 m、1.45 m、1.85 m，半径为0.15 m金属管进行正演模拟。图5为几何模型与正演模拟图。图5(b)、图5(c)、图5(d)分别对应250 MHz、500 MHz、900 MHz天线B-Scan。

图5 不同埋深管线几何模型与正演模拟Fig.5 Geometrical model and forward modeling

由图5可见，250 MHz天线分辨率最差，500 MHz天线分辨率较好，900 MHz天线分辨率最好。但900 MHz天线的反射波能量弱，在实际运用中能量衰减快。500 MHz天线与其它两个频率天线相比无论是分辨率还是反射波能量都较均衡，故在文中主要选择500 MHz天线进行正演模拟。不同频率天线随着管线埋深不断增大，双曲线能量也不断减弱，双曲线曲率也不断变小。故探测深部目标体时，在保证分辨率前提下尽可能使用低频天线。

3.1.4 不同材质管线正演模拟

图6为不同质管线的几何与正演模拟图，在路基层中设置三根管径0.25 m、埋深0.9 m、壁厚0.05 m不同材质的管线。图6(a)从左到右依次为混凝土管、PVC管、铸铁管，图6(b)、图6(c)、图6(d)为使用500 MHz天线时相应B-Scan。

图6 不同材质管线几何模型与正演模拟Fig.6 Geometrical model and forward simulation

由正演模拟图可知三种材质管线B-Scan均为开口向下双曲线，目标与背景介质介电常数相差越大，反射波幅值就越大，故管顶反射波能量依次增强。金属管线与背景介质介电常数差异大，电磁波在金属管线管顶发生全反射，所以金属管线管顶反射波能量最强。

背景介质相对介电常数大于PVC与混凝土的相对介电常数，小于金属的相对介电常数，故PVC与混凝土管顶反射波同相轴相位与直达波相位一致，金属管管顶反射波相位与直达波相位相反，与正演模拟图一致。

3.1.5 不同管径大小正演模拟

为研究不同管径大小管线的雷达剖面图，现设计三根平行埋设的混凝土管，管径分别为0.03 m、0.1 m、0.3 m，管顶距地面均为0.9 m，雷达天线频率为500 MHz。

从正演模拟图7中可知，随着管径不断增加，管顶反射波能量不断增大，双曲线两翼延伸越长，双曲线曲率也越大，故管径越大越容易识别。管径为0.03 m与0.1 m的管线因管径较小，顶、底反射叠加在一起，未观察到管底反射。A处异常为上部管线绕射波相互叠加形成，在管线探测工作时要注意区分。

3.1.6 分辨率

雷达天线频率越高分辨率越高，但能量衰减快、探测深度小。1.5 GHz、2.0 GHz等高频天线常用来进行路面层厚、钢筋混凝土结构检测；250 MHz、500 MHz等中频天线常用来进行管线探测；100 MHz、50 MHz 等低频天线常用来做地质勘察、超前地质预报。探地雷达分为横向分辨率和纵向分辨率，为了研究不同频率天线的分辨率，依次对100 MHz、250 MHz、500 MHz和900 MHz天线进行正演模拟。

图7 不同管径管线几何模型与正演模拟Fig.7 Geometrical model and forward simulation

1)横向分辨率

近间距管线探测一直是雷达探测一个难点，为了获得100 MHz、250 MHz、500 MHz、900 MHz天线横向分辨变率大小，设计近间距平行管线：设置六根埋深、半径均为0.1 m的内部充气金属管，管线间水平距离分别为0.25 m、0.3 m、0.4 m、0.5 m、0.6 m。

图8(b)、图8(c)、图8(d)、图8(e)分别为100 MHz、250 MHz、500 MHz、900 MHz天线的B-Scan，当相邻管线间距足够大时管线较易区分；当相邻管线间距较小时，两者双曲线相互干扰，甚至重叠，并伴有绕射波出现，易对探测结果进行误判，对断定管线数目带来难度；反射波到时依次减小，但管顶反射波在水平上的位置不变，遵循管间中心定律。100 MHz天线仅能辨别有异常体，但无法区分管线数目；250 MHz天线勉强可分辨出间距为0.5～0.6 m的管线；500 MHz天线可分辨出间距为0.4～0.6 m的管线；900 MHz可分辨出0.3～0.6 m的管线。

图8 近间距管线几何模型与正演模拟Fig.8 Geometrical model and forward simulation

2)纵向分辨率

图9(a)为管顶距地面0.85 m，与下部管线间距分别为0.05 m、0.2 m、0.5 m，半径为0.1 m的金属管线组，分别使用100 MHz、250 MHz、500 MHz、900 MHz天线研究其纵向分辨率。如图9(b)、图9(c)、图9(d)、图9(e)分别为相应的B-Scan。

由图9可知，100 MHz未能分辨出下部管线，仅对上部管线有反应，纵向分辨率差；250 MHz与500 MHz天线可分辨0.2 m及0.5 m垂直间距管线，但500 MHz分辨率更高；900 MHz天线均可分辨出上部和下部管线，且反射波双曲线清晰，但能量相对较弱。

图9 垂向叠置管线几何模型与正演模拟Fig.9 Geometrical model and forward simulation

在一定条件下天线频率越高，分辨率也越高。在满足探测深度条件下，尽可能选择高频天线探测。此外，在判别近间距管线时要慎重，避免因所选择天线分辨率不够而引起误判。

3.1.7 不同空间位置管线正演模拟

城市高速发展不可避免地增加了地下管线的密度，相邻管线及管线相对位置的不同都会给管线探测带来难度。下面模拟常见不同管线空间位置雷达剖面图，见图10。其中图10(a)、图10(c)分别为“上大下小”、“上小下大”模型图，图10(b)、图10(d)为使用500 MHz天线时相应正演模拟图。

图10 几何模型与正演模拟Fig.10 Geometrical model and forward simulation

当管径上大下小时，可见上部大管径管线顶、底反射，下部小管径反射不明显，勉强可见双曲线两翼，且能量较弱。当上部为小管径管线时，可清晰辨别出管顶反射，下部大管径管线反射波两翼清晰，因小管径管线在其正上方，电磁波衰减，故下方管线管顶处反射波能量弱。因此在进行垂直埋设管线探测时，要特别注意“上大下小”情况，必要时加密测线，避免漏测。

图11(a)、图11(c)分别为“正三角”、“倒三角”模型图，图11(b)、图11(d)为相应正演模拟图。在“正三角”时，由相应正演模拟图可分辨出管线个数和相对空间位置，底部绕射波也较明显。“倒三角”模拟图中，上部两管线反射波双曲线清晰可见，中部为两管线反射波交叉叠加产生，下部为底部管线顶反射。故在“倒三角”管线探测时要认真研读雷达剖面图，必要时可复测，对比研究是近间距管线交叉产生的反射波还是下部管线管顶反射波。

图11 几何模型与正演模拟Fig.11 Geometrical model and forward simulation

综上所述，在对不同空间位置管线探测时，对雷达剖面图要认真研究，必要时可通过复测或加密测线，避免漏测或因雷达剖面图信息匮乏发生误判。

3.1.8 不同填充介质正演模拟

在管线普查工作中常有排水管线积水、淤堵，为了研究该类填充介质管线的雷达剖面图，现设计管径0.25 m、壁厚0.02 m的混凝土管线，内部分别充气、满水，如图12(a)所示。

图12(b)、图12(c)分别为使用500 MHz天线时充气、满水正演模拟图，充气管线管顶反射波相位与直达波相位一致，因水的相对介电常数大，满水线管顶反射波相位与直达波相位相反，且伴有多次波，充气管线仅管顶反射波明显。电磁波在空气中的传播速度最快，在水中的传播速度最慢，故充水管线的顶底反射到达时差最大。在一定条件下，可参考反射波相位、能量及顶、底到达时差推断管内填充物质。

图12 几何模型与正演模拟Fig.12 Geometrical model and forward simulation

3.1.9 管顶回填不密实

填土不密实区主要集中在路基层，主要是由明挖敷设管线时因施工质量问题导致管顶回填区域不密实，电磁波穿过该区域会引发“假信号”，这给管线探测带来一定难度，并给今后的路面变形、塌陷等病害埋下隐患。图13(b)～图13(d)分别为使用500 MHz天线时对应充气、给水及金属管线雷达剖面图，三者路基层至管顶范围的不密实区均有能量较弱的反射信号，类似圆形管道反射波形态，但同相轴错断、杂乱、不连续。一般回填不密实区域下方管线均可被识别，但对充气非金属管线影响大，管顶反射波能量被削弱，对给水、金属管线影响不大。

图13 几何模型与正演模拟Fig.13 Geometrical model and forward simulation

3.1.10 被掩埋井正演模拟

因后期道路维修等原因，市政道路检查井部分井盖被沥青、混凝土等材料覆盖，这对以后的管道维护与普查带来了巨大不便。为了解决该问题，现模拟被掩埋检查井模型，见图14。

图14 几何模型与正演模拟Fig.14 Geometrical model and forward simulation

设置检查井井盖为金属材质，井筒、井室及管道均为混凝土材质，使用500 MHz天线。图14(b)为检查井井盖正演模拟图，电磁波在金属界面发生了全反射，顶部能量极强，易被识别。图14(c)为测线方向与管线敷设一致时的正演模拟图，可清晰辨别管线顶底反射，类似于层状介质雷达剖面图。图14(d)为雷达测线与管线敷设方向垂直时的正演模拟图，为典型圆形管道雷达剖面图。寻找掩埋检查井一般先确定井盖位置，然后在井盖四周布置网格测线，根据雷达剖面图确定雷达测线是与管线敷设方向平行(图14c)还是垂直(图14d)，进而给管线走向提供科学依据。

4 城市工程勘察实践

为了验证正演数值模拟的可靠性，本文结合南宁市排水管网普查工作中的已知管线为工程实例，使用MALA探地雷达进行管线探测，提高对不同形状、材质等雷达剖面图的分辨能力，为以后的探测工作提供强有力的技术指导。

图15 PE给水管剖面Fig.15 PE water supply pipe profile

实例1为某污水泵站出水管，管径0.5 m，材质为PE(Polyethylene)的充水管，设计埋深为2.1 m。图15为使用250 MHz屏蔽天线探测并进行简单滤波处理的成果图，在水平距离3.5 m、深度2.2 m处出现双曲线状反射波，符合圆形管线反射波特征。管顶反射波相位与直达波相位相反，这是因为基层介电常数比水的介电常数小，电磁波在由基层进入水中时电磁波相位反向。由于电磁波在水中衰减较快，地下介质复杂且250 MHz天线在市政道路中普遍仅能探测2.0 m左右，故仅可见管顶反射。此外，在管顶反射波上方两侧均有较规则异常信号，分析应为敷设管道时所开挖的沟槽两壁引起的绕射波。

实例2为钢材质压力管，设计埋深2.0 m，管径0.6 m。图16为使用250 MHz屏蔽天线探测的雷达剖面图，管线位于水平距离1.96 m、埋深约2.05 m处，反射波能量强；在3.7 m处出现了二次反射。管顶上方区域同相轴错断，推断为管顶回填不密实，管顶上方约0.6 m处疑为钢管混凝土保护层反射波，波形特征整体符合金属管线反射波形态。

实例3为圆形管道与矩形暗涵的雷达剖面图，图17、图18分别为使用250 MHz天线与500 MHz天线探测同一条测线的成果图。图17中250 MHz天线雷达剖面图中电磁波能量强，探测深度大，但相较于图18的500 MHz天线雷达剖面图，图17的分辨率低、电磁波能量强，且图18可反映水平距离8.0～12.0 m处矩形暗涵顶板的钢筋反射波。与文中高频天线分辨率高,探测深度小的结论一致。

由图17、图18可知：①处反射波为双曲线型，符合圆形管道反射波形态，推断为圆形混凝土管线，管顶距地面约0.86 m；②处顶部呈水平板状反射波形态对应暗涵顶部盖板，通过暗涵顶部盖板钢筋分布与多次波“X”字型形态可推断为地下暗涵，暗涵宽约4.6 m，通过顶部反射波到时推断暗涵顶部距地面约1.1 m；③处为金属消防管，在管顶发生全反射，反射波能量强，推断埋深约0.67 m，探测结果与设计图基本一致。此外，图17、图18中底部电磁波衰减较快，表现为低频信号、多次波丰富，疑管道内部部分充水，尤其圆形管道反射波杂乱疑管内淤堵含泥。