面向基础设施的三维探地雷达属性成像及信息可视化研究

朱家松, 雷占占, 罗享寰

(1.深圳大学土木与交通工程学院, 广东 深圳 518060; 2. 深圳大学城市智慧交通与安全运维研究院,广东 深圳 518060; 3. 广东省城市空间信息工程重点实验室, 广东 深圳 518060; 4. 深圳市地铁集团有限公司, 广东 深圳 518026)

0 引言

探地雷达(GPR)是一种广泛应用于城市基础设施内部测量和成像的无损检测方法,具有分辨率高、非破坏性和非接触式连续测量等优点[1]。GPR的成像方式分为一维(A-scan)、二维(B-scan)和三维(C-scan)3种。利用B-scan对城市地下基础设施测量时,需检查一系列相邻的剖面才可以确定地下目标的位置和大小。而C-scan可以从俯视图的角度,以水平切片的方式来呈现测量结果,具有直观且全覆盖的优点,可表现目标物体的空间位置及几何特性。因此,基于C-scan图像的目标提取研究逐渐成为GPR领域的热点[2-3]。

由于三维C-scan图像的空间分辨率可以达到毫米级,这导致包含大粒径颗粒的非均匀介质,如砂土中的孤石、混凝土大孔隙或地下植被根系等非目标反射物会在C-scan上表现为雪花状噪声,从而影响目标物体的准确提取。此外,目前GPR C-scan所展现的测量值仅有雷达波的反射能量,而雷达波包含的相位、色散、波速等信息都未得到充分利用,因此,解译时可依据的信息非常有限。然而,城市基础设施多是钢筋混凝土结构,基础设施检测中不同构筑物往往会产生相似的反射强度,难以将其区分,进一步导致C-scan的多解性问题严重。目前,针对GPR C-scan图像质量低、信息少的问题,国内外已有不少研究。例如: Luo等[4]量化了C-scan空间分辨率与对象参数之间的关系,建立了面向目标对象的标准化C-scan成像流程。近年来,为了从地球物理数据中提取更多的信息来解释和探测目标,以提高数据解译的准确度和效率,研究人员将探测信号的属性提取与分析技术引用到GPR数据解释过程中[5-11],其中基于小波变换和基于复信号技术的瞬时属性提取方法在应用中取得了良好的效果。Moeller等[12]证明了利用频率建立C-scan的可行性,然而对于多种属性融合的C-scan成像效果尚未有量化的比较分析。

C-scan图像质量通常由空间分辨率(像素尺寸)和颜色(像素值)2部分组成。目前,国内外已有许多学者关注C-scan的分辨率,即由步长、天线频率和采样密度等测量参数决定的测量信息密度分布,但是关于GPR C-scan中像素值(信息可视化映射)的研究却鲜有报道。C-scan的解译是以人类视觉感知为基准,并没有关于映射的客观标准规则,解译结果常带有人为主观误差,甚至引发漏判或误判。一般来说,C-scan的像素值由测量信息值映射到色彩区间得到,因而测量信息值的范围和信息映射方法是决定C-scan色彩的2个主要变量[13-17]。Lu等[18]指出,一个合理的色彩是符合人类视觉感知的。因此,优化信息映射的主要目的是提高对比度和消除噪声的影响,突出目标物体的特征,从而提高解译的准确性。

综上,目前三维GPR C-scan图像主要存在可直接读取的信息维度太低以及信息可视化方式主观2个主要问题。因此,本文的研究内容由2部分组成: 1)提出一种融合波属性信息的GPR C-scan图像成像方法,丰富C-scan图像信息; 2)定量分析C-scan不同测量信息可视化的成像效果,建立适合城市典型应用场景的优化成像策略,最终提高C-scan的解读性以及GPR测量的解译精度。

1 属性C-scan成像及色彩优化方法

1.1 属性C-scan成像

为实现GPR属性提取与属性C-scan成像,本研究基于python开发了GPR成像程序,主要功能包括基于B-scan的属性提取与分析,提取属性为瞬时振幅与瞬时频率、GPR测量值的空间编码以及属性C-scan图像的生成。GPR成像及解译程序的用户界面如图1所示。

(a) 信号坐标编码

(b) 波属性提取

1.1.1 基于复信号技术的属性提取方法

GPR属性提取过程是通过数学统计与计算处理GPR数据的过程。GPR的属性存在许多不同的分类方法,可以从计算的角度把GPR属性提取划分为2种类别: 1)依据单道波进行计算的GPR属性,比如频率、振幅、相位属性; 2)依据多道波进行计算的GPR属性,比如相干性和相似性。在这些理论的基础上,本文采用单道波的复信号分析方法,提取GPR回波在目标物体深度的三瞬属性,即电磁波信号的瞬时相位、瞬时振幅和瞬时频率。其具体方法为: 通过采集GPR的实信号进行希尔伯特变换,将各个点的实信号通过卷积运算得到其实信号的虚部,再将两者构成复信号并以三角函数的形式表达,然后从已分解的三角函数形式中提取各组信号的瞬时反射强度、瞬时频率和瞬时相位[19]。由于瞬时相位主要由应用在GPR二维剖面上的同相轴连续性进行度量,因此,本文采用瞬时频率与瞬时振幅作为三维图像的研究属性。

1.1.2 空间插值成像

提取GPR的回波属性后,将其正交均匀排布,然后数值标准化以消除物理量纲的影响,并通过反距离加权法进行空间插值,以填充测量值采样点间的空隙,形成全覆盖的、高分辨率的C-scan图像。

其中,反距离加权法的假设条件是事物间距越接近,则事物彼此越相似。其假定每个测点都对预测位置有局部影响,且离预测位置越近的测点局部影响越大。在使用反距离加权法插值时,测点的分配权重函数随距离的增大而减小,测点与预测位置间的距离越近,则其被分配到的权重越大[20]。

由于GPR在试验时采集的数据类型为剖面B-scan,所以,需要进行水平方向与垂直方向的重采样,以模拟现实中的地表三维GPR图像。根据Luo等[4]的三维C-scan成像流程,以目标物体直径为C-scan厚度,对属性的剖面图进行垂直方向的重采样,每个划分深度为1组平面C-scan数据集,得到空间分辨率一致且临近像素间无冲突的各属性C-scan图像。

1.2 信息映射优化

在形成全覆盖的C-scan后,需要对C-scan中的测量信息值进行颜色赋值,以具现化测量信息。因C-scan的数据由三维空间信息(坐标编码)和测量值2部分组成,而测量值仅为单一波段的相对反射强度信息,单色的灰度图更符合真实情况。把测量值投影到灰度空间有2个步骤,分别为有效测量值阈值和映射函数。

1.2.1 有效测量值阈值

由于海量非目标杂波会在C-scan图像上表现为雪花状噪声,而细碎的噪声往往集中在测量值分布的高低两端,定义灰度图像的上下边界可有效排除非目标回波[21]。将目标物体和周围材料产生的GPR信号的测量组成为分布直方图来描述其振幅分布,将各介质的测量值归一为0～100,如图2所示。

本研究采用标准差量化描述有效测量值的区间。首先,从直方图中确定平均值;然后,根据标准差寻找与平均值有一定距离的上限和下限。其中,定义期望值为μ=0、标准差σ=0.5条件下的正态分布,则当σ=0.5,1,1.5,2,2.5,3,3.5,4时,其灰度取值范围为全测量值的68.27%～99.99%。上下限之间的测量值为有效测量值,而高于上限和低于下限的测量值则将被去除。

1.2.2 映射函数

确定灰度边界后,将有效测量值映射成灰度值,如下式所示。

S=Tr。

式中:S为转化后的像素值;T为映射函数;r为映射前的测量值。

在8 bit图像中,灰度图像的像素值为0～255,如何分配灰度是信息映射的关键。

目前,GPR成像领域广泛应用线性变换[13],因为它涉及较少的人为操作且提供了最完整的反射强度值。然而,针对不同的测量环境以及特定的测量物体,简单的线性转换未必是最合适的方式。例如: 当成像目标是一个不均匀的界面时,使用线性变换会生成一个噪声较多且很模糊的C-scan图像,而模糊的C-scan图像难以描绘测量区域的细节信息。特别是当相邻区域具有相似的介电特性时,低对比度的颜色会导致图像结构信息的损失。除了线性映射以外,广泛使用的4种不同的灰度映射函数分别为对数函数、平方根函数、平方函数和指数函数。其中: 对数函数和平方根函数关注较暗区域,可将像素值过大的区域做适当灰度减弱;而平方函数和指数函数关注较亮区域,即将像素值过小的区域做适当灰度增大。本文采用定量分析的方法详细评估每种函数在基础设施成像中的适用度。

2 基于实验室环境的成像试验

2.1 试验数据

以典型的GPR应用案例——混凝土墙体检测为例,混凝土墙GPR测线设计图如图3所示。

图3 混凝土墙GPR测线设计图

在混凝土墙上的正交网格进行GPR检测,并收集墙体的GPR回波数据。GPR测线由17条平行于X轴和16条平行于Y轴的测线组成,测线间隔为100 mm。为了更容易地识别管道的影响,在设计网格时将网格横轴纵轴交点中心设置在L形PVC管道的渗流点处。数据采集使用天线中心频率为2 GHz的GPR设备,设置采集样本的时间窗为6 s,道间距为3 mm,每个信号采集512个样本。

在进行三维C-scan成像前,利用Reflexw软件对采集的原始GPR回波数据进行基本的信号处理,其步骤包括去直流漂移、增益调整、零时校正、带通滤波以及滑动平均。

信号处理前后的GPR B-scan样本如图4所示。虽然墙体内钢筋与PVC管道呈双曲线状反射,且清晰可见,但两者的图谱特征太过相似,即使两者引发的反射信号相位相反,但仍难以区分。

(a) 信号处理前

(b) 信号处理后

2.2 试验结果

2.2.1 属性灵敏度分析

在利用波属性成像前,需确认雷达波属性是否对不同介质响应不一,且需明确振幅或频率对介质的灵敏性,以区分属性在融合成像上的权重。在试验区中选择20个点位,涵盖混凝土、钢筋以及PVC管道,并提取这些点位的GPR回波信号。对已经过信号处理的回波信号进行复信号分析,最终得到关于混凝土、钢筋以及PVC管道的三瞬属性。从中选取瞬时振幅与瞬时频率,并提取相应的属性数据,将所提取的瞬时属性的变化趋势绘制成如图5所示的曲线。

(a) 瞬时振幅

(b) 瞬时频率

在图5(a)中,钢筋与PVC管道的曲线起伏变化较大,混凝土的曲线相对平稳,而且钢筋在每一探测位置的GPR回波峰值振幅均比PVC管道大,因此,可以利用该振幅区分钢筋与PVC管道。在图5(b)中,钢筋曲线的中心频率变化较大,而PVC管道与混凝土的频率相对平稳,且每个点位钢筋的中心频率都比管道的频率大。由频率的灵敏度可见,钢筋与PVC管道的响应频率曲线存在明显不同,可依据该属性区别钢筋与PVC管道。

2.2.2 属性C-scan成像

采用全测量值线性变换的钢筋混凝土墙体C-scan与理想C-scan对比,结果如图6所示。

(a) 瞬时振幅

(b) 瞬时频率

(c) 精调最优

首先,采用自建程序确定试验区内钢筋网和PVC管道的深度,以15～30 mm为1个切片深度,根据采样频率确定C-scan的合适深度为150～300 mm。然后,根据当前业界普遍采用的全测量值线性投影的方法形成属性C-scan原始图像,如图6(a)所示。灵敏度分析表明,PVC管道与混凝土的介电常数相接近,其振幅反射灰度值也相似,因此,振幅灰度图无法判断PVC管道的位置,但频率灰度图可以看清PVC管道的部分形态。基于对试验场景的了解,人为精调得到接近真实情况的理想C-scan(如图6(c)所示),可清晰界定同深度的钢筋及PVC管道。

2.2.3 成像优化结果评价与试验结果讨论

为量化评价不同的测量值可视化方式对于C-scan图像的优化效果,本文选取峰值信噪比(peak signal-to-noise ratio, PSNR)[22]与结构相似性(structural similarity, SSIM)[23]作为成像优化指标。PSNR越大,表示图像质量越高;而SSIM的取值范围在0～1,当SSIM为1时,表示试验C-scan与理想C-scan的结构完全一样。因为SSIM的计算不受色彩饱和度的影响,仅依据图像中物体几何形状的完整度,适合用于本研究的目标物体区分评价。

对于测量信息可视化的2个变量——测量值阈值和映射函数,本文采用控制变量法研究其合适取值。首先,映射函数固定为常用的线性映射不变,测量信息为振幅,改变测量值阈值σ,形成不同测量值阈值的振幅C-scan图像,与理想C-scan对比计算SSIM(如图7所示)。比较SSIM数值后得知,阈值σ=1.0时得到最高SSIM,为0.90。然后,固定映射函数为线性映射不变,测量信息为频率,改变测量值阈值σ,同样当阈值σ=1.0时得到最高SSIM,为0.85。

同理可得,本文分别以振幅和频率为对象,研究映射函数的选取。固定阈值σ=1.0,改变不同映射函数,形成映射函数不同的C-scan图像,与理想C-scan对比计算SSIM。映射函数为平方函数时,SSIM最高,即成像最接近理想C-scan。试验结果表明,可视化设置为有效阈值σ=1.0和平方映射函数时,振幅和频率图像的PSNR和SSIM均达到最高值,分别为32.74/0.90和29.96/0.85。

为探索融合属性是否可提高C-scan图像的可解读性,将振幅和属性的图像按一定权重进行叠加,并根据上述试验结果,测试测量值阈值属性和映射函数的成像效果。结果表明,当有效测量值阈值σ=1.0,且灰度映射函数为平方函数时,融合属性的SSIM最高,为0.90,高于振幅和频率的最大值,且变化趋势与振幅和频率相符。融合属性C-scan的PSNR高于单属性频率图像,与单属性振幅图像接近(见图8)。

(a) 瞬时振幅优化

(b) 瞬时频率优化

(c) 属性融合优化

2.2.3.1 属性区别信息映射

上述试验结果表明在选定同样的有效测量值阈值(相同的σ值)与灰度映射方式时,瞬时振幅属性C-scan的SSIM与PSNR明显高于瞬时频率属性C-scan,而且多属性融合后的SSIM结果明显高于单属性图像,因为同时具备振幅和频率的图像所能读取的信息更加丰富,所以与真实值更加接近。

2.2.3.2 有效测量值阈值σ

在属性与灰度映射方式相同时,C-scan的成像效果仅由有效测量值阈值决定。上述试验结果表明,噪声主要集中在分布直方图的两端,应适当截选。分析上述结果可知,当σ=1.0时,即包含75.72%原始数据范围时,优化效果最好。此时由于除掉了大部分的非目标回波,使得目标对象更加突出,因此,PSNR与SSIM均达到了最大值。当有效测量值范围逐步增大时(σ增大),优化效果逐渐减弱。当σ大于1.5时,试验结果趋于稳定。

2.2.3.3 信息映射方式

在属性与有效测量值阈值固定时,C-scan的成像效果仅由灰度映射函数决定。试验结果表明,在该组数据中平方函数和线性函数明显优于其他信息映射函数,主要是由于钢筋和混凝土的介电常数差异大, GPR在界面的反射强度大。平方函数与指数函数能较好地提高部分低反射区域的强度,使得目标对象更加清晰。同时,由于介质较为均匀,去噪需求较低,而线性函数一般较符合现实情况,可保留较多的测量值。对数函数与平方根函数不适用于低反射强度区域,因此,无论是PSNR还是SSIM均低于常规图像数值。

3 基于道路现场的成像试验

3.1 试验数据背景

道路GPR C-scan成像试验以香港石门一段长为20 m的道路为测量对象。该道路模型分为2个部分,左半部分为标准的沥青道路结构,右半部分为标准的钢筋混凝土道路结构,其面层包含0.1 m×0.2 m网状钢筋,路基下为不均匀砂土。道路结构中埋设2种地下管线设施,分别为PVC管道与球墨铸铁管道。石门道路地下结构及GPR测线图如图9所示,正交网格为GPR测线路径。GPR数据由IDS双频系统(200/600 MHz)按照网格测线采集,并经过必须的信号处理。

S为管道直径; D为管道埋深。

3.2 属性成像与优化结果

3.2.1 属性融合成像结果

经相关预处理并通过上述属性提取与三维成像流程后,完成该案例的属性C-scan图像。

由于该试验场地预埋材料复杂,砂土、混凝土与PVC管道均具有接近的介电常数,且3种材质的反射振幅值接近,导致振幅灰度图像较为模糊,无法有效识别对应位置的介质。经对比,实地试验中有效测量值阈值为1.0,映射函数为平方函数时,C-scan成像的SSIM最高,即效果最好。其优化后的C-scan图像如图10所示,2根管线均清晰可见,其PSNR与SSIM均为最高,分别为30.78和0.87。

(a) 瞬时振幅

(b) 瞬时频率

(c) 属性融合

3.2.2 属性成像及映射优化的参数评价

石门道路试验场地变参数成像效果如图11所示。

上述试验结果表明,在选定同样的有效测量值阈值(σ值)与灰度映射方式时,瞬时振幅C-scan的SSIM与PSNR明显低于瞬时频率C-scan,且多属性融合后的图像SSIM结果明显高于单属性图像。这主要是因为同时具备振幅和频率的图像所能读取的信息更加丰富,能表现出更多的测量区域介质特征。

在确定相同的属性与灰度映射方式时,当测量值阈值σ=1.0,即包含75.72%原始数据范围时,成像效果最好。试验的地下环境较复杂,噪声较多,需要更窄的有效测量值,此时由于除掉了大部分的背景噪声,使得目标检测物体更加突出,因此,PSNR与SSIM达到最大值。当所选取范围逐步增大时,优化效果逐渐减弱,试验结果变化趋于稳定。

在确定相同的属性与有效测量值阈值时,实地试验中平方函数和线性函数明显优于其他灰度映射函数,主要是由于钢筋混凝土、管道和砂土的GPR反射强度差异不显著,而平方函数与指数函数能较好地提高部分低反射强度区域,使得优化效果更好。综上所述,城市基础设施的结构相对复杂,介质种类较多,但单一介质内空间变异性不大,平方函数的映射方式可强化局部特征。

4 结论与建议

本研究提出一种基于波属性的GPR C-scan成像方法,通过组合振幅和频率提高C-scan信息维度,研究标准化的C-scan像素值映射方法,建立适合城市典型应用场景的面向对象C-scan成像策略。

通过开发三维GPR C-scan成像及解译程序,在实验室和实地道路构建典型GPR应用场景,定量分析三维GPR测量值可视化的2个主要参数,即测量值阈值与灰度映射函数的取值对C-scan成像的影响,并以PSNR与SSIM作为图像质量的评价指标。混凝土墙体和实地道路的试验结果一致,对于钢筋混凝土结构,适当去除高频低频噪声,保留70%有效测量值,信息映射方式为平方函数时,组合属性的三维GPR图像成像效果最好,可有效突出目标物体。对于环境不太复杂的城市基础设施,需要较窄的有效测量值和更强化的灰度映射方式。因此,面向不同的测量场景和需求,建议建立规范的测量值可视化流程,可在保证C-scan目标物体区分精度的同时,减少人为操作误差,进一步提高GPR测量的精度和准确度。