水泥路面脱空区域GPR信号的时频域分析

2024-04-01余秋琴罗婷倚李有鑫

现代雷达 2024年1期

余秋琴,罗婷倚,朱欣,李有鑫,张军*,刘斌

(1. 广西北投公路建设投资集团有限公司, 广西南宁 530028)(2. 公路养护装备国家工程实验室, 陕西西安 710064)(3. 道路施工技术与装备教育部重点实验室, 陕西西安 710064)

0 引言

截至2021年底,全国公路总里程达519.8万千米,由于养护不及时导致大量路面病害得不到及时处置,迫切需要实现公路路面病害的精准识别,为路面预养护提供技术指导。探地雷达(GPR)是公路病害有效的无损检测方法[1-2],然而地下结构复杂,路面病害类型众多,导致GPR信号解译困难[3],急需研究GPR信号的有效解译方法。

实际道路交通情况复杂,很难开展现场实验验证。因此,采用正演模拟来掌握各种病害特征或地下目标的成像分布规律成为当前研究热点[4-8]。GprMax采用时域有限差分算法(FDTD)软件[9],可有效模拟电磁波在不同病害特征的响应及其图谱,广泛应用于公路、隧道的脱空病害成像特征分析[10-11],其结果有助于对实际探测空洞识别和雷达图像的解释。探地雷达信号特征提取方法主要分为时域,频域和时频域等方法。如采用多个时频特征表征病害特征,并结合机器学习方法,可实现病害的智能识别[12]。也有研究者通过信号变换来提高信噪比,从而提高GPR解译精度[13-15];文献[16]采用连续小波变换(CWT)解决复信号分析易受噪声干扰的问题。目前CWT的研究工作主要集中在提高信噪比等方面,而没有用于提取病害特征。

针对GPR信号解译困难的问题,以水泥路面为对象,采用正演模拟和反演实验相结合的方式,研究不同脱空形态和填充物对GPR信号响应的影响。在此基础上,为提取脱空的特征,采用编程软件设计了GPR信号处理模块,用于处理正演模拟和室内模型数据,并结合连续小波变换和3D可视化方法实现了脱空特征提取,本文的研究方法可为表达水泥路面脱空区域的几何特征提供依据。

1 GPR信号特征提取方法

1.1 GPR信号特征提取理论

探地雷达的回波信号是瞬变非平稳信号,传统傅里叶变换由于基函数的特性,无法表达信号的局部时频域性质,进而导致有些信号特征不明显。小波变换是一种窗口大小固定而形状可改变的时频局部化分析方法,用于对信号详细的时频分析或准确的信号瞬态定位,特别是对于瞬时频率增长迅速的信号[17]。因此,联合时频分析优点,针对GPR信号提出使用连续小波变化进行特征提取,可以凸显雷达信号的时频特征,连续小波变化表达式为

(1)

式中:x(t)为GPR原信号;ψa,b(t)为母波;a和b分别是母波的尺度和平移参数。

根据小波变换内积定理,小波变换的幅值平方在尺度-位移平面的加权求和等于信号在时域的总能量,如式(2)。因此,小波变换的幅平方可看作是信号能量的时频分布的一种形式[18]。

(2)

(3)

式中:ψ(ω)为ψ(t)的傅里叶变换。

考虑到在时间与频率的局部化之间有很好的平衡,以及提取更多的特征信息,母小波采用非正交复Morlet (Gabor)小波,该小波函数是一个零均值并且由复余弦调制的高斯函数,它的时频方差相等,没有尺度函数。母小波的频域解析函数定义见式(4),相应的时域解析式见式(5)。

(4)

(5)

式中:U(ω)是Heaviside 阶跃函数;ω为角频率。

1.2 GPR信号特征提取方法

脱空信号的GPR后处理方法如图1所示,包括抽取平均道、静矫正、增益、背景去除和F-K偏移等。对于正演模拟数据,由于信号缺少干扰,不做抽取平均道工作。

图1 脱空信号的基本处理模块Fig.1 Basic processing protocol for void signal

根据图1的要求,在编程软件中开发了GPR的后处理模块,数据处理流程如图2所示,首先读取GPR数据,接着对数据进行预处理,最后将GPR数据进行可视化显示,并给用户提供操作主面板。设计的处理模块可对正演模拟数据和常用的商用雷达数据进行处理,从而方便了数据格式转换,为病害分析提供便利。

图2 GPR处理模块的工作流程Fig.2 Flow chart of GPR processing module in softuare

2 路面脱空模型设计

考虑水泥路面不同脱空形状及含水的状态,构建了两个脱空模型,模型尺寸为1.6 m×0.45 m,脱空区域位置主要分布在深度0.2 m附近,已知800 MHz电磁波在空气中的波长是0.375 m。如图3所示,模型A考虑了不同脱空尺寸和脱空形状(正方形、圆形和梯形)的影响,脱空区域内的介质均为空气。其中 A-1和A-2边长为10 cm和8 cm的正方形脱空,A-3和A-4直径为10 cm和6 cm的圆形脱空,A-5是脱空梯形,上下边长分别为10 cm和20 cm。模型B考虑脱空区域的含水状态(全充气、全充水、空气和含水和潮湿物料)的影响。其中B-1和B-2为20 cm×15 cm矩形脱空,分别考虑脱空区域为空气和水的两种介质;B-3为20 cm×25 cm矩形,考虑脱空区域为水和空气的混合体,水的体积占比为总体积的2/5;B-4为直径10 cm的拼接圆,介质为高含水区域,介电常数为25,对应含水率10%[19],电导率为0.102 S/m。此外,混凝土的介电常数为7,电导率为0.005 S/m,水的介电常数为81。

图3 脱空模型示意图Fig.3 Void models of cement pavement

参照水泥路面脱空检测情况[20-21],仿真模型天线中心频率采用800 MHz,单道A-Scan数据采样时窗为12 ns,详细仿真参数见表1。

表1 FDTD正演模拟参数表Tab.1 FDTD forward simulation parameters

3 结果分析

3.1 不同形状对脱空特征的影响

将正演模型图3创建的两个模型导入仿真软件进行仿真,获得的仿真数据经过静矫正、增益、背景去除,滑动平均等处理后结果如图4所示。两模型总长为1.6 m,步长0.005 m,总道数为(1.6-0.01)/0.005=318(道),单道A-Scan有512个采样点。

图4 处理后的B-ScanFig.4 B-Scan image after processing

进一步对脱空区域的特征进行分析,图4a)为充气脱空区域的图谱特征,圆形脱空区域呈现抛物线特征;矩形为平直的“白黑白”的亮斑,平直程度与矩形的长度相关,而两端形成短直边,总体呈现拱桥形状;梯形区域因两次反射在图谱上表现为梯形形状,且带有矩形的特征。图4b)为空气、充水或高含水率下的脱空图谱特征。对比B-1和B-2区域的特征,由于二者介电常数差异,图谱表现为“白黑白”和“黑白黑”的亮斑,波相出现明显的反相特征,为此,可通过波相特性判别脱空区域是否含水。B2和B4由于介电常数高,GPR信号在其底部出现多次反射,如B-4潮湿混凝土不仅在其病害区域出现双曲线特征,还在底部存在二次反射。此现象表明,如果脱空区域含水,会导致电磁波反相,从而根据反射特征可判断脱空区域是否含水。

3.2 脱空信号的时频特征分析

为掌握脱空区域的信号特征,将图4b)结果经过F-K处理后,再进行FFT频谱计算,结果如图5所示。

图5 模型B的时频域分析Fig.5 Time-frequency domain analysis of model B

如图5a)所示,脱空区域的GPR信号幅值要远高于正常信号幅值。GPR信号在B-2(纯水)与B-1(纯空气)、B-3(水和空气混合)处的反射回波在第1 ns时极性相反,因为B-2的第一反射面反射系数R小于零,而B-1、B-3的反射系数是大于零的,导致波相相反,与前面的B-Scan图谱特征一致。信号在B-3区域的1 ns、2.5 ns、4 ns等三处有反射回波,对应了反射面为混凝土-空气-水-混凝土;信号在B-4(含水区域)与B-2的7 ns和10 ns存在二次反射,极性与B-2号1 ns处相同,原因是信号B-4和B-2的部分反射回波因在地表面和介质面之间经过多次反射造成时间延迟差异。因此B-2和B-4表明存在高介电常数的物质会导致GPR信号多次反射,这个特征可有利于脱空区域判断是否含水。从图5b)可以得到各介质回波的频谱主峰能量大依次为:B-3>B-2>B-1>B-4>正常,各谱峰能量集中600 MHz～1 500 MHz且主峰位置存在差异,而且在高频处,各信号能量差异较明显;但仅依靠时域和频域信号仍难以掌握脱空区域的信号特征。

3.3 脱空信号的时频域分析

上述分析表明脱空与正常路面在时域和频域的信号存在差异,但是难以在时间轴上对病害区域进行定位,为此,将上述模型的数据进行进一步连续小波变换分析,对比不同脱空状态下的时频域差异,结果如图6所示。

图6 不同含水脱空信号的连续小波变换Fig.6 CWT result for different water content signal

通过统一颜色域及坐标的取值范围,从图6中可以明显看出不同脱空状态信号时频差异明显,含水、脱空或潮湿区域的能量比正常GPR信号能量大。图6a)为纯空气脱空,能量集中区域的时间约为0.53 ns～2.57 ns,对应的深度为0.03 m～0.145 7 m,图6b)为纯含水区域,能量集中在0.64 ns～2.4 ns ,按照波速0.113 4 m/ns计算,对应深度0.036 m～0.136 m;而实际模型B-1上、下反射面距地表面为0.05 m～0.20 m,因此能量区域可反应脱空病害的深度信息。此外,图6a)与图6c)处脱空含水相比,在9 ns处有差异,介质为脱空含水的高频对应亮度要比充气脱空区域的亮度和区域要更亮和更宽广,在时间域跨度也更宽。图6b)与图6d)相比,在第7 ns和第9 ns处有差异,图6b)在第2 ns左右显得更亮,这些特性表明,CWT可以在时频域上将脱空区域通过高能量区突出病害区域。

进一步可分别计算每个切片(每道GPR信号的CWT变换结果)最大的幅值,结果如表2,可以看出介质为正常混凝土的幅值最小,介质为含水脱空的幅值最大,比正常混凝土的幅值大很多;介质为潮湿混凝土的幅值也比正常差异明显。从图6c)与图6e)可以明显看出在相同频率分量处,脱空病害要比正常区域的更亮,即病害区域的能量较大,且能量集中区域与脱空深度位置相关。

表2 不同序号CWT最大值Tab.2 Maximum CWT values for different targets

3.4 脱空病害三维可视化

CWT可有效突出病害区域,进一步将模型B中所有GPR数据进行CWT变换,再重构为四维CWT信息图谱。为了实现绘图展示,用色彩亮度表示为第四维度,以方便读者理解,图形仍是3D图。图7为模型B的CWT 3D切片图,坐标轴分别对应时间、频率和深度的数组大小,非具体物理值。按照频率轴和水平距离(道数)进行切片,得到8个CWT切片。按照测量道数的水平距离切片,可以看出脱空以及含水脱空区域与正常GPR信号的图谱差异大,能量集中区域与病害区域一一对应。而按频率轴进行的横向切片,可进一步突出脱空特征。CWT的3D切片图相比图6,可进一步表明脱空的特征及其图谱特征值差异,表明CWT的3D图是表征脱空的病害的又一方法。

图7 模型CWT 3D可视化Fig.7 CWT for 3D visualization of model

根据表2求取的正常混凝土的最大幅值作为背景阈值,去除图7的3D模型背景,可得到图8结果,背景已被去除,可进一步突显出整个模型B中不同脱空病害的空间区域分布。这为图6不同脱空信号的时频差异性提供很好的解释,而这种方式也为构建脱空3D区域可视化提供参考。

图8 背景去除的3D CWT特征Fig.8 CWT features after background removal

考虑到图5b)不同病害频谱特征差异,截取频率分量为1 388.64 MHz的横截面,并对其数据3D可视化,如图9a)所示。图中可明显分辨不同病害特征的时频及能量差异性分布,表明采用CWT可有效揭示不同含水脱空的时频域特征。此外,根据每道A-Scan数据的CWT最大亮度值对应的频率分量,提取时间-亮度信号向量(主频率信号),将所有道的数据作类似处理,然后重构B-Scan能量谱,如图9b)所示,能量区域差异不仅反应不同病害的差异性,也反应了不同介质的分布状况。图中B-3在前100采样点与B-1颜色值相近,而在采样点100-200与B-2很相近,因此可确定B-3介质为B-1和B-2的混合,该方法为分析病害特征状态提供方法参考。

图9 能量谱分析Fig.9 Energy spectrum analysis

4 工程应用

为了验证CWT可有效提取脱空特征的正确性,构建了室内物理模型如图10a)所示,模型大小为2.07 m×0.4 m,采样频率是40 GHz,时窗为11.675 ns,总道数为207道。考虑了不同脱空尺寸的影响,并设计间隔为15 cm的10个脱空圆孔,直径依次为100 m、90 m～30 m,25 m,16 mm,各圆心位置位于20 cm中心线上。模型材料为水泥混凝土,并采用US 1G天线进行测试,使用所设计的软件后处理模块对测量数据进行基本处理和F-K偏移处理,及CWT处理;图10b)为室内物理模型3D CWT切片图,可以有效揭示脱空区域的时频域特征。随着圆直径逐渐变小,每道截面的亮斑区域不断缩小,其亮度也随之降低。

图10 实验模型结果Fig.10 Laboratory model CWT result

为了表明CWT的有效性,按照病害的深度将采样时间为2.34 ns结果进行截图,病害在深度方向的结果如图10c)所示,可以得出前7个孔洞尺寸较大的脱空能有效观测,并且能量集中程度反映了孔洞尺寸大小;且随着空洞直径逐渐减小,各空洞的频率分布带也在减小,但主要的能量频带相对集中。为表征特征频率对病害结果的影响,将脱空的敏感频率为1 302.1 MHz结果进行切片, 横截面映射到B-Scan能量谱如图10d)所示,可以得出,除第一个空洞幅值很大,当圆型直径尺寸的变化量为3 cm(空洞2和空洞5的直径差),脱空区域幅值会突变比较明显,椭圆亮斑也会变得越来越窄暗。第一空洞幅值特别大是因为现场采用雷达测试在开始处出现干扰,而空洞10因直径太小与背景融合无法分辨。此外,采用CWT明显提升高频信号的时间分辨率,通过矫正时间,可以确定前9个病害位置。验证了采用CWT可有效表征脱空区域,为水泥路面脱空识别提供参考。