地表起伏对地下管线GPR探测的影响.pdf

第 47 卷第 4 期物 探 与 化 探Vol.47,No.4 2023 年 8 月GEOPHYSICAL&GEOCHEMICAL EXPLORATION Aug.,2023doi:10.11720/wtyht.2023.1516曾波,刘硕,杨军,等.地表起伏对地下管线 GPR 探测的影响J.物探与化探,2023,47(4):1064-1070.http:/doi.org/10.11720/wtyht.2023.1516Zeng B,Liu S,Yang J,et al.Influence of surface undulations on GPR-based underground pipeline detectionJ.Geophysical and Geochemical Explora-tion,2023,47(4):1064-1070.http:/doi.org/10.11720/wtyht.2023.1516地表起伏对地下管线 GPR 探测的影响曾波1,刘硕2,杨军3,冯德山2,袁忠明3,柳杰3,王珣2(1.广州市城市规划勘测设计研究院,广东 广州 510060;2.中南大学 地球科学与信息物理学院,湖南 长沙 410083;3.广州市市政工程设计研究总院有限公司,广东 广州 510060)摘 要:地下管线是城市的重要设施,承担着能源输送、信息传递等功能,为城市生活提供了便利和保障。探地雷达(ground penetrating radar,GPR)作为一种高分辨率、高精度、非开挖、非破坏性的探测技术,在管线测量中具有巨大的优势。然而地表地形复杂,多有起伏,对于 GPR 探测地下管线有很大的影响。因此,本文采用有限单元法对地下管线探测进行了数值模拟,该方法可以与非结构网格结合,更好地拟合地表起伏地形;此外,介绍了如何进行高度校正,将所得剖面数据与地形相吻合,更容易分析异常体特征。最后开展了两个数值实验,分析了起伏地表对于不同埋深、不同间距、不同材质及不同填充物管线探测的影响,为 GPR 数据解释提供理论基础。实验结果表明,因地表起伏原因,波形和反射波能量将发生畸变,并不能作为判断管线信息的唯一依据。因此,需进行高度校正,利用双曲线的顶点来判断管线的埋深、材质等信息。关键词:探地雷达;管线探测;地表起伏;有限单元法;高度校正中图分类号:P631 文献标识码:A 文章编号:1000-8918(2023)04-1064-07收稿日期:2022-10-24;修回日期:2023-06-09基金项目:国家自然科学基金项目(42074161、42104143);湖南省自然科学基金项目(2021JJ30806、2022JJ40584);湖南省研究生科研创新项目(CX20220188)第一作者:曾波(1977-),男,高级工程师,主要从事地下管线探测和地下空间测绘工作。Email:75813052 通讯作者:冯德山(1978-),男,博士,教授,博导,从事地球物理正反演与数据处理研究工作。Email:fengdeshan 0 引言随着城市的发展和进步,地下管线发挥着越来越重要的作用。城市地下管线种类繁多,有电线电缆、用于水和油气运输的管道以及复杂的地下排水系统1。然而,由于城市建设初期地下空间信息的不完善,在城市的持续建设过程中,地下管线往往会面临老化、改道等问题。这些问题不仅会影响到正常的城市生活,而且还会造成火灾、断水、停电等事故2。为了避免此类事故的发生,需要对地下管线进行全面的调查、记录和系统管理。管线探测仪是目前城市地下管道检测中较为常用的设备,但此类设备只能检测金属管道3。而探地雷达(ground penetrating radar,GPR)不仅能探测金属和非金属管道,还具有高效、显示实时、使用灵活、非破坏性等特点,是管线检测的首选设备4-5。利用探地雷达进行地下管线探测时,因地下介质的非均匀性、管线错综复杂以及人为干扰等因素,会导致电磁波传播过程十分复杂,雷达数据难以解释6。因此,需要开展探地雷达正演模拟,了解地下 GPR 波的传播特征,为数据的解释奠定理论基础。梁小强等7应用时域有限差分算法(finite difference time domain,FDTD)8-9开展了不同情况下探地雷达探测效果的分析,了解并掌握了管线探测与各种影响参数之间的关系;王书等10将 FDTD算法应用到了地下水污染调查中,说明 GPR 可划分地下水的含水层界面;Gao 等11采用 FDTD 算法对探地雷达非平稳信号进行了分析,讨论了不同特征的地下管线的雷达信号及泄漏油的变化情况,提高了地下管线和污染土壤的雷达解释精度。但 FDTD算法通常用于直交网格,不能与非结构化网格结合,4 期曾波等:地表起伏对地下管线 GPR 探测的影响对于边界几何不规则的地质体,易产生虚假反射12-13。雷建伟等14采用共形辛 Euler 算法高效、精细地识别了地下非金属管道,有效降低了由于阶梯近似造成的虚假绕射波,但是该算法还是基于结构化网格,对于拟合非规则边界有一定的缺陷。时域 有 限 单 元 法(finite element time domain,FETD)15-16能够与非结构化网格结合,较好地拟合地质体的边界,在地球物理领域数值模拟得到了广泛应用17-18。因此,文中将采用有限单元法对不同情况下的管线进行正演模拟,分析雷达剖面的图像特征。此外,近地表地形复杂,有一定起伏。实际探测管线时,各种反射波、绕射波、干扰波混淆在一起,无法区分与解释19-20。侯爵等21开展了起伏地表下的地震偏移成像数值实验,表明其方法的灵活性。刘四新等22利用 FDTD 算法对典型的起伏地形模型进行数值模拟,并重点介绍了机载探地雷达的发展应用。但将 FETD 算法应用到起伏地形下管线的GPR 探测中,目前还未见到相关报道。因此,文中对比了水平地表和起伏地表对于管线探测的影响,分析了不同间距、不同深度、不同材质、不同填充物管线的雷达剖面特征,对实际探测 GPR 数据的解释工作提供了理论依据。1 方法原理方法原理分两部分来阐述,第一部分为有限单元法推导 Maxwell 方程组,得到其离散格式,用来模拟 GPR 模型;第二部分为高度校正,实际工作时,一般会在地表进行探测,因此要做高度校正,将雷达剖面向实际地形契合,方便 GPR 数据的处理和解释。1.1 有限单元法求解 GPR 波动方程由电磁波传播理论可知23,雷达波在有耗媒介中的传播遵循 Maxwell 方程组。时间域二维 TM 模式下电场标量波动方程可以表示为-12Ezx2-12Ezy2+Ezt+2Ezt2+Jzt=0,(1)其中,为磁导率(H/m);x,y 为电磁波传播方向;Ez为 z 方向上的电场分量(V/m);为电导率(S/m);t 为时间(s);Jz为 z 方向上的外加电流密度(A/m2)。采用 FETD 算法求解上式。将整个模拟区域 采用三角单元进行离散,假设每个单元内部的电性参数相同,则单元内部线性插值函数为16Eex,y()=3j=1NejEj,(2)其中,Nej为形函数,只与该三角单元的顶点坐标、面积和 x、y 有关;Ej为 3 个顶点上的电场数值。将式(2)代入到式(1),乘以形函数,并在单元内部积分,则可以写为ee+KeEe+CeEe+Fe=0,(3)式中:Me为质量矩阵;Ke为刚度矩阵;Ce为阻尼矩阵;Ee为时间的一次导数;e为时间的二次导数;Fe为源项。它们具体的表达式为e=eeNTNdxdy,(4)Ke=Kex+Key=1eNTxNx+NTyNy()dxdy,(5)Ce=eeNTNdxdy,(6)Fe=eeNTJztdxdy。(7)将式(3)扩展到全局,并进行单元合成,可以得到整个求解域的 GPR 有限元微分方程组:+KE+CE+F=0,(8)式中的一阶及二阶时间导数采用中心差分近似:En=12tEn+1-En-1,(9)n=1t()2En+1-2En+En-1,(10)其中,En为 n 时刻的电场值;En-1为 n-1 时刻的电场值;En+1为 n+1 时刻的电场值。通过这种时间离散方案,即可得到 FETD 算法求解 Maxwell 方程组的递推公式,模拟 GPR 电磁场在时域上的进程,其迭代公式为Mt()2+C2t()En+1=2Mt()2-K()En+C2t-Mt()2()En-1-fn,(11)因上式方程左边每一个时间步的求解不包括刚度矩阵 K,因此是显式算法。其数值稳定性条件为t 2 M-1K(),(12)式中:M-1K()表示矩阵 M-1K 的谱半径,也为其最大特征值。1.2 高度校正在起伏地表做管线探测时,得到的雷达剖面中异常体位置与实际情况有一定的差距,因此需要做高度校正,将数据与实际地形相契合,这样有利于5601物 探 与 化 探47 卷 GPR 数据的处理与解释。所得到的雷达剖面,横向为水平距离,纵向为时间,因此需要得到高度差所引起的“滞后时间”值,具体公式为t=2Hv,(13)式中:t 为时间;H 为高度差,即某一标准水平线与地表测点之间的垂直距离;v 为电磁波传播速度,与背景介质有关:v=cr r,(14)其中,c 为真空中的光速,为 3108m/s;r为背景介质的相对介电常数;r为背景介质的相对磁导率,一般为 1。2 数值案例以下数值案例中的模型均用 Comsol 软件进行Delaunay 三角剖分,在异常变化较大的区域会对网格进行加密,更加贴合其边界,提高计算精度。2.1 地形对不同埋深、不同间距管线探测的影响设计一个不同埋深、不同间距的金属管线模型,如图 1 所示。模拟区域大小为 2.5 m1.2 m。其背景相对介电常数为 6,电导率为 0.01 S/m,相对磁导率为 1。模拟区域中有 4 个金属管(从左到右分别为 A,B,C,D 管),如图中红色圆形所示,其相对介电常数为 1,电导率为 107 S/m,半径为 0.12 m,从左到右圆心依次为(0.5 m,1 m)、(0.94 m,1 m)、(1.5 m,0.5 m)和(1.84 m,0.5 m),间距分别为 0.2 m 和0.1 m。图 1a 为水平地表模型的网格示意,总共 70 424 个单元,35 563 个节点,00.2 m 为空气层;图1b 为起伏地表模型的网格示意,整个模拟区域的单元数为 66 414,节点数为 33 566,与水平地表模型的网格数大约相同。地表平均高度为 0.2 m,上下起伏约为 0.2 m。两者均在地表采用自激自收方式采集数据,激励源为中心频率 600 MHz 的 Ricker 子波,蓝色三角表示激发点(接收点)。采样间隔为0.01 m,总共接收 250 道数据;时间采样间隔为 0.008 ns,总记录时长为 30 ns。图 2a、b 为水平地表和起伏地表的正演剖面。从中可以看到,金属管线作为导体,将高频电磁波全反射,反射系数为-1,因此其相位与直达波相反。在图 2a 中,因为是水平地表,其双曲线的位置对应a水平地表;b起伏地表aflat surface;bragged surface图 1 不同埋深、不同间距的模型网格示意Fig.1 Grid diagram of the model with different burial depths and different spacingsa水平地表;b起伏地表aflat surface;bragged surface图 2 正演剖面Fig.2 Forward profiles6601 4 期曾波等:地表起伏对地下管线 GPR 探测的影响管线的埋深,左边两个金属管线(A 和 B)的反射波振幅明显小于右边两个管线(C 和 D)。此外,因为C 管线和 D 管线的间距较近,所以在其剖面下方,会有一道很明显的反射曲线,这是电磁波在两个管线之间来回传播的反射波;在 A 管和 B 管的剖面下侧也可以看到,只不过因其间距偏大一些,与金属管的反射双曲线离得较远。在图 2b 中,可以明显看到 4个金属管的强反射波,由于起伏地表的原因,不能对应管线的埋深,C 管的反射波也因起伏地表而变形,非对称的双曲线。图 3 为高度校正后的剖面。为了与起伏地表对比,同样将水平地表的剖面图做了统一的高度校正,即图 1 模型中 0 m 处的位置。由图可见,做完高度校正的两个剖面,金属管的双曲线位置大约在同一高度,但振幅却不相同。而且,B 金属管的双曲线并不完整,因起伏地表和另外金属管的作用下,其反射波被其他杂波所湮没。因此,起伏地表所得到的剖面,并不能根据其振幅和双曲线的位置判断其埋深,需要做高度校正,且剖面较水平地表的剖面较为混乱。a水平地表;b起伏地表aflat surface;bragged surface图 3 高度校正后的正演剖面Fig.3 Forward profiles with height correction2.2 地形起伏对不同材质、不同填充物管线探测的影响图 4 为不同材质、不同填充物的模型示意。模拟区域大小为 10 m4 m。其背景相对介电常数为9,电导率为 0.001 S/m,相对磁导率为 1。模拟区域中有 3 个管线,最左边的红色管线为金属管(A管),其相对介电常数为 1,电导率为 107 S/m,半径为 0.3 m,圆心为(2.5 m,2.5 m);中间的白色管线为充满空气的 PVC 管(B 管),其相对介电常数为1,电导率为 0,半径为 0.3 m,圆心为(5 m,2.5 m);最右边的蓝色管线为充满水的 PVC 管(C 管),其相对介电常数为 80,电导率为 0.002 S/m,半径为 0.3 m,圆心为(7.5 m,2.5 m)。图 4a 为水平地表模型的网格示意,总共 65 378 个单元,33 048 个节点,00.8 m 为空气层;图 4b 为起伏地表模型的网格示意,整个模拟区域的单元数为 66 003,节点数为33 364,与水平地表模型的网格数大约相同。地表平均高度为 0.8 m,上下起伏约 0.8 m。两者均在地表采用自激自收方式采集数据,激励源为中心频率100 MHz 的 Ricker 子波,蓝色三角表示激发点(接收点)。采样间隔为 0.037 m,总共接收 270 道数据;时间采样间隔为 0.03 ns,总记录时长为 150 ns。a水平地表;b起伏地表aflat surface;bragged surface图 4 不同材质、不同填充物的模型网格示意Fig.4 Grid diagram of the model with different materials and fillers7601物 探 与 化 探47 卷 图 5a、b 为水平地表和起伏地表的正演剖面。由图可见,A 金属管的反射波最强,且相位与直达波相反;B 空气管的反射波较弱,相位与直达波相同;C水管由于与背景介质的相对介电常数相差较大,因此形成了多次反射波,其相位与直达波相反。图 5a中,3 个管线的反射曲线大约在同一时刻,因此在水平地表探测管线时,可以根据反射曲线顶点的位置判断其埋深,根据其相位和反射振幅、是否有多次波来判断管线材质、填充物等信息。在图 5b 中,因为在起伏地表上探测,3 个管线的反射曲线已不在同一高度,但是同样看到,A 金属管的反射能量最强,相位与直达波相反;B 空气管的反射能量最弱,相位与直达波相同。图 6 为高度校正后的剖面。为了与起伏地表对比,同样将水平地表的剖面图做了统一的高度校正,即图 4 模型中 0 m 处的位置。由图可见,直达波的形态与起伏地表相契合,管线的双曲线大约校正在了同一高度。对比图 6a、b,起伏地表下的水管反射波能量要强于水平地表下的反射能量,且两翼并不完整,这是由于地表在水管正上方凹陷,垂直距离近,且在上坡地段探测时,能量被下凹界面反射,到达水管的能量较弱,因此双曲线两翼能量较弱。a水平地表;b起伏地表aflat surface;bragged surface图 5 正演剖面Fig.5 Forward profilesa水平地表;b起伏地表aflat surface;bragged surface图 6 高度校正后的正演剖面Fig.6 Forward profiles with height correction3 结论1)采用有限单元法对 GPR 管线探测进行了正演模拟,该算法可以与非结构网格结合,能更好地贴合实际地表。推导了 FETD 算法求解 GPR 波动方程的公式,并介绍了高度校正的具体实现过程,有利于雷达数据的处理与解释。2)通过两个数值案例,对不同间距、不同埋深、不同材质、不同填充物的管线进行正演模拟,分析了起伏地表对探测的影响。实验结果表明:起伏的地表会导致探测到的反射波变得畸形,电磁波能量不8601 4 期曾波等:地表起伏对地下管线 GPR 探测的影响均匀,并不是对称的双曲线;此外,不能简单地通过管线在探测结果中的上下位置来判断埋深,需要进一步进行高度校正。高度校正过后,可以根据双曲线顶点来推断管线的埋深,也可根据振幅、相位、是否有多次波来大致推断管线的材质和填充物。需要注意的是,因为地表起伏的原因,波形较为复杂、混乱,其反射波的能量不能作为唯一的判断依据,因此应结合其他信息来进行 GPR 数据解释。参考文献(References):1赵欣,王希良,刘珍岩,等.复杂条件下的地下管线探测模拟J.物探与化探,2014,38(6):1307-1312.Zhao X,Wang X L,Liu Z Y,et al.Simulation of underground pipelines under complicated condition J.Geophysical and Geo-chemical Exploration,2014,38(6):1307-1312.2 Park B,Kim J,Lee J,et al.Underground object classification for urban roads using instantaneous phase analysis of Ground Pene-trating Radar(GPR)data J.Remote Sensing,2018,10(9):1-24.3 姚显春,闫茂,吕高,等.地质雷达探测地下管线分类判别方法研究J.地球物理学进展,2018,33(4):1740-1747.Yao X C,Yan M,Lyu G,et al.Research on underground pipeline classification and discrimination method based on geological radar detection J.Progress in Geophysics,2018,33(4):1740-1747.4 韩佳明,仲鑫,景帅,等.探地雷达在黄土地区城市地质管线探测中的应用J.物探与化探,2020,44(6):1476-1481.Han J M,Zhong X,Jing S,et al.The application of geological radar to urban pipeline detection in the loess area J.Geophysical and Geochemical Exploration,2020,44(6):1476-1481.5 曾昭发,刘四新,冯晅,等.探地雷达原理与应用M.北京:电子工业出版社,2005.Zeng Z F,Liu S X,Feng X,et al.Theory and application of ground penetrating radarM.Beijing:Electronics Industry Press,2005.6 肖敏,陈昌彦,贾辉,等.金属管线对探地雷达探测道路地下病害的干扰J.物探与化探,2016,40(5):1046-1050.Xiao M,Chen C Y,Jia H,et al.The study of the interference region around metal pipeline in underground disease detection of urban road J.Geophysical and Geochemical Exploration,2016,40(5):1046-1050.7 梁小强,杨道学,张可能,等.FDTD 数值模拟在 GPR 管线探测中的应用J.地球物理学进展,2017,32(4):1803-1807.Liang X Q,Yang D X,Zhang K N,et al.Application of FDTD nu-merical simulation of Ground Penetrating Radar in pipeline detec-tion J.Progress in Geophysics,2017,32(4):1803-1807.8 Yee K.Numerical solution of initial boundary value problems invol-ving maxwells equations in isotropic media J.IEEE Transac-tions on Antennas&Propagation,1966,14(3):302-307.9 葛德彪,闫玉波.电磁波时域有限差分法M.西安:西安电子科技大学出版社,2005.Ge D B,Yan Y B.Finite-difference time-domain method for elec-tromagnetic wavesM.Xian:Xidian University Press,2005.10 王书,闫天龙.地下水污染调查中探地雷达有限差分数值模拟J.物探与化探,2016,40(5):1051-1054.Wang S,Yan T L.The extraction of sounding curves from the data of high-density resistivity method for intepretation J.Geophysical and Geochemical Exploration,2016,40(5):1051-1054.11 Gao L,Song H,Liu H,et al.Model test study on oil leakage and underground pipelines using ground penetrating radar J.Rus-sian Journal of Nondestructive Testing,2020,56(5):435-444.12 李静,曾昭发,吴丰收,等.探地雷达三维高阶时域有限差分法模拟研究J.地球物理学报,2010,53(4):974-981.Li J,Zeng Z F,Wu F S,et al.Three dimensional high-order FDTD simulation for GPR J.Chinese Journal of Geophysics,2010,53(4):974-981.13 李静,刘津杰,曾昭发,等.基于变换光学有限差分探地雷达数值模拟研究J.地球物理学报,2016,59(6):2280-2289.Li J,Liu J J,Zeng Z F,et al.Study of GPR simulation based on the transformation optics FDTDJ.Chinese Journal of Geophysics,2016,59(6):2280-2289.14 雷建伟,方宏远,李银萍,等.基于共形辛 Euler 算法的非金属地下管道精细化高效探地雷达正演模型J.地球物理学报,2020,63(8):3192-3204.Lei J W,Fang H Y,Li Y P,et al.GPR forward model of under-ground non-metallic pipeline based on parallel conformal symplec-tic Euler algorithm J.Chinese Journal of Geophysics,2020,63(8):3192-3204.15 徐世浙.地球物理中的有限单元法M.北京:科学出版社,1994.Xu S Z.The finite element method in geophysicsM.Beijing:Sci-ence Press,1994.16 Jin J M.The finite element method in electromagnetics M.John Wiley&Sons,2014.17 冯德山,王珣.基于卷积完全匹配层的非规则网格时域有限元探地雷达数值模拟J.地球物理学报,2017,60(1):413-423.Feng D S,Wang X.Convolution perfectly matched layer for the finte-element time-domain method modeling of Ground Penetrating Radar J.Chinese Journal of Geophysics,2017,60(1):413-423.18 Zhang Z,Wang H H,Wang M L,et al.Non-Split PML boundary condition for finite element time-domain modeling of ground pene-trating radar J.Journal of Applied Mathematics and Physics,2019,7(5):1077-1096.19 王洪华,王敏玲,张智,等.基于 Pade 逼近的 Cole-Cole 频散介质 GPR 有限元正演J.地球物理学报,2018,61(10):4136-4147.Wang H H,Wang M L,Zhang Z,et al.Simulation of GPR in Cole-Cole dispersive media by finite element method based on Pade ap-proximation J.Chinese Journal of Geophysics,2018,61(10):4136-4147.20 王洪华,吕玉增,王敏玲,等.基于 PML 边界条件的二阶电磁波动方程 GPR 时域有限元模拟J.地球物理学报,2019,62(5):1929-1941.Wang H H,Lyu Y Z,Wang M L,et al.A perfectly matched layer for second order electromagnetic wave simulation of GPR by finite 9601物 探 与 化 探47 卷 element time domain method J.Chinese Journal of Geophysics,2019,62(5):1929-1941.21 侯爵,刘有山,兰海强,等.基于起伏地形平化策略的弹性波逆时偏移成像方法J.地球物理学报,2018,61(4):1434-1446.Hou J,Liu Y S,Lan H Q,et al.Elastic reverse time migration using a topograpghy flattening scheme J.Chinese Journal of Geophys-ics,2018,61(4):1434-1446.22 刘四新,冯彦谦,傅磊,等.机载探地雷达的进展以及数值模拟J.地球物理学进展,2012,27(2):727-735.Liu S X,Feng Y Q,Fu L,et al.Advances and numerical simulation of airborne ground penetrating radar J.Progress in Geophysics,2012,27(2):727-735.23 Taflove A,Brodwin M E.Numerical solution of steady-state electro-magnetic scattering problems using the time-dependent maxwells equations J.IEEE Transactions on Microwave Theory and Tech-niques,1975,23(8):623-630.Influence of surface undulations on GPR-based underground pipeline detectionZENG Bo1,LIU Shuo2,YANG Jun3,FENG De-Shan2,YUAN Zhong-Ming3,LIU Jie3,WANG Xun2(1.Guangzhou Urban Planning and Design Survey Research Institute,Guangzhou 510060,China;2.School of Geosciences and Info-Physics,Central South University,Changsha 410083,China;3.Guangzhou Municipal Engineering Design&Research Institute Co.Ltd.,Guangzhou 510060,China)Abstract:As important urban facilities,underground pipelines perform the functions of energy transfer and information transmission,providing convenience and guarantee for urban life.Ground-penetrating radar(GPR),as a high-resolution,high-precision,trenchless,and non-destructive detection technique,has great advantages in pipeline surveys.However,undulating surfaces with complex terrain greatly influence GPR-based detection of underground pipelines.Therefore,this study conducted numerical simulations of the under-ground pipeline detection using the finite element method,which can be combined with an unstructured grid to fit the undulating sur-faces effectively.Furthermore,this study introduced the height correction method to match the obtained geologic sections with terrain,making it easier to analyze the anomaly characteristics.Finally,through numerical experiments,this study analyzed the influence of un-dulating surfaces on the detection of pipelines with different burial depths,spacings,materials,and fillers,providing a theoretical basis for GPR data interpretation.The experimental results show that waveforms and reflected wave energy,subjected to distortion due to sur-face undulations,cannot be used as the sloebasis for judging pipeline information.Therefore,height correction is required,and the verte-xes of hyperbolas can be used to judge the burial depths and materials of pipelines.Key words:ground penetrating radar;pipeline detection;surface undulating;finite element method;height correction(本文编辑:叶佩)0701