极点特征聚类的公路隐藏裂缝自动识别算法

周丽军

(山西交通科学研究院集团有限公司 智能装备研究院，太原 030032)

目前对公路下隐藏病害的无损探测主要依赖于探地雷达(GPR)技术，即向地下空间发射高频电磁波，利用地下物体与环境之间的介电差异对接收回波进行处理，提取地下介质分布情况[1-2]。应用最广泛的GPR探测技术主要是基于B-Scan扫描，通常得到的扫描图像是基于同一场景的多方位多角度测量，并根据周围介电环境重构目标像从而得到目标的相对位置，但是并不能通过信号幅值或双曲线峰值区分不同材质的目标[3]。而且仅通过B-scan图像很难对目标进行精确的解析，识别难度大，主观因素高[4]。

近年来，有学者将人工智能技术用在探地雷达公路病害检测中。文献[5]建立人工神经网络多层模型，用于道路病害的自动诊断。文献[6]利用深度字典学习方法处理探地雷达回波信号，区分多个不同形状的掩埋目标。文献[7]将向量量化神经网络方法用于铁路路基的病害识别中。文献[8]利用聚类的方法进行雷达图像异常边缘的自动识别。这些基于B-scan图像的方法还主要集中在物理模型中，多数是识别雷达图像的双曲线特征，对目标的介电属性识别方面的进展还十分有限[9]。文献[10]通过不同核匹配追踪算法能有效识别水害和空洞，可识别城市道路土基病害。

随着时域脉冲雷达的发展，基于雷达回波瞬态响应的目标识别方法得到了较大的发展，相比于成像技术，其优势在于仅需要目标时域测量的单个瞬态响应。根据SEM技术，目标依赖于其瞬态响应回波晚时部分的复自然谐振(CNR)，即极点特征，其分布规律只与目标属性有关，而与入射波的入射方式等无关，此性能能够很好的用于目标识别领域[11]。如文献[12]利用极点差异识别不同形状或不同方位的目标。近年来，这些技术也已经应用于GPR探测的地下目标识别。利用瞬态谐振现象估计与识别地下目标，往往需要得到目标的冲激响应。文献[13]从时域积分方程的角度研究了目标的脉冲响应以及其所包含的内部谐振电流，并指出目标的瞬态散射特性可以完全由其脉冲响应表示。文献[11]利用解卷积方法获取目标的冲激响应，并对冲激响应进行奇异值分解，获取目标晚时响应的CNR，研究半空间介电目标的谐振行为。文献[14]与文献[15]基于瞬态响应的内、外部谐振分别研究了CNR与目标几何属性的关系。但是上述文献尚缺乏谐振与目标介电特征的深入研究。

针对上述问题提出一直基于极点特征聚类的地质雷达公路隐藏裂缝自动识别方法。通过对回波信号进行奇异值分解，提取体现目标信号的大奇异值，重构回波信号，由此不仅能够对回波信号进行降维，还去除了噪声干扰。再利用解卷积过程计算目标的冲激响应，利用奇点展开法(SEM)提取其晚时响应部分的极点信息，构建极点特征空间，并用模糊C均值聚类方法对极点进行分类，根据已知环境样本对极点进行初始聚类中心优化并取其均值，再利用此均值初始聚类中心对待识别目标回波的极点特征空间进行聚类，输出分类结果。

1 算法描述

1.1 奇异值分解的降维方法

由于采集的数据量大，并且伴有很多噪声，对回波数据进行奇异值分解，获得前M个主奇异值，用此重构原始信号。方法如下：

对噪声数据y(kΔt)，其中Δt是采样间隔，k为采样点数，构造Hankel矩阵Y：

Y=

(1)

式中:N为数据采样点数，L为函数束参数。

对Y进行奇异值分解以获得特征矢量及特征值：

Y=U∑VH

(2)

式中：矩阵U、V分别为Y的左奇异矩阵和右奇异矩阵，由矩阵YYH、YHY的特征矢量组成，上标H表示共轭转置。对角矩阵∑由Y的特征值组成:

(3)

由于σc(c=1，2，…，N-L)按照从大到小顺序排列，比较每个σc和最大奇异值的比值直至满足σc/σmax≫10-p(p代表精度，默认值为采样数据的小数位数)，即可确定M的取值。取∑的前M列∑′，重新构造数据，显然重构数据维数远低于原始数据。并且由此得到的数据去除了噪声干扰，突出了目标信号信息。

1.2 目标冲激响应极点提取与极点特征空间构成

入射波从目标体表面返回到接收机而不再作用到目标上时，目标体上残留的被激发(或诱导)的电磁波会产生与目标介电性能相关的谐振。理论上，此谐振电磁波在方位上独立，仅依赖于目标的形状与介电特性，被称为晚时响应。为了得到目标的谐振信息，需要将目标的瞬态响应从系统响应中提取出来。为了减少数据量，去除噪声，采用各回波信号经过奇异值分解后的重构信号。

首先测量背景信号，即没有目标时的回波响信号yb(t)，再测量有目标时的目标回波响应yt(t)，设目标的冲激响应为h(t)，则卷积过程可以写成：

yb(t)*h(t)=yt(t)

(4)

于是可以通过解卷积过程得到目标的冲激响应h(t).

此时目标的冲激响应h(t)可以写成如下形式：

h(t)=hop(t)+u(t-t0)∑Riesit

(5)

其中hop(t)表示冲激响应的早时响应，上式右边第二部分则是冲激响应的晚时响应，与目标的固有属性有关，其中t0为晚时响应的开始时刻，其计算方法可参考文献[16]。式中Ri表示第i个谐振态的复幅值，si=αi+j2πfi表示第i个谐振态的极点，其中αi与fi分别为衰减因子和谐振频率。

根据获得的衰减因子αi与谐振频率fi构成极点特征空间如下：

[A，F]={(α1，f1)，(α2，f2)，…，(αM，fM)}

(6)

其中M为重构数据获得的可用极点个数。将前M个谐振态的极点分布在式(6)的二维极点特征空间，能得到不同目标冲激响应的特征分布，从而通过极点特征进一步识别不同目标。

1.3 模糊C均值分类算法

模糊C均值(FCM)聚类算法是一种无监督方法，它是通过对目标函数进行迭代优化获得各个样本点对所属的聚类中心的隶属度，从而对样本进行划分分类[17]。

(7)

(8)

由于模糊C均值聚类算法依赖于聚类中心的初始化值，本文根据多组已知的背景、空气裂缝、充水裂缝的极点特征空间分布调整优化初始聚类中心，使其趋于稳定。

1.4 算法流程

第一步：降维。为了减少庞大的数据量，利用奇异值分解方法对对探地雷达的背景回波、空气裂缝回波、充水裂缝回波信号进行数据降维，同时也能对采集信号进行降噪。

第二步：构建特征空间。对回波信号提取冲激响应，计算极点特征，并构建以衰减因子与谐振频率为极点的特征空间模型，分别计算出背景、空气裂缝、冲水裂缝的极点特征空间。

第三步：对背景、空气裂缝、冲水裂缝的极点特征空间中的衰减因子-谐振频率对进行模糊C均值聚类，根据已知的极点特征空间与聚类结果调整聚类算法的初始聚类中心。

第四步：对投影到极点特征空间的待识别信号极点根据调整优化后的初始聚类中心进行聚类，输出分类结果。

2 实验结果分析

本文使用GSSI公司的SIR-20系列天线对地下空气裂缝，充水裂缝进行实验分析，为了模拟空气裂缝与充水裂缝，分别取空PVC管与充满水的PVC管进行试验，管长1 m，直径4 cm.天线的中心频率为1 GHz，维度为49.5 cm×21 cm×55.6 cm.由于瞬态响应是基于平面波信号源，为了将发射信号近似成平面波，预先对天线放置高度进行经验值估算，如图1所示，在天线架高77.9 cm时，可以近似形成面积范围139 cm×111 cm的平面波区域，即目标在此区域内其入射波可近似为平面波，后向接收回波为瞬态响应回波。实测场景如图2所示。

图1 地质雷达入射波近似平面波示意图Fig.1 Gram of approximated plane wave from GPR

图2 实测场景Fig.2 Real scene

分别获取30组背景雷达回波信号、50组空气裂缝雷达回波信号、50组充水裂缝雷达回波信号作为极点特征空间训练样本，再各取5组空气裂缝与充水裂缝回波作为测试样本。

首先对回波数据进行降维，既减少了数据量，又进一步去除了噪声干扰。在奇异值分解的降维方法中，选取合适的特征值重构原始信号，这里精度取0.01，使得重构信号既能完整的表征目标信息，又能去除冗余与噪声，50组空气裂缝雷达回波样本与50组充水裂缝雷达回波样本的主要特征值个数及重构信号与原始信号的误差如图3所示。

图3 重构信号与原始信号的误差Fig.3 Error between reconstructed signal and original signal

从图中可以看出，空气裂缝与充水裂缝的重构信号与原始信号误差都在0.01以内，在此误差范围内，充水裂缝重构信号所需的主特征值数量要大于空气裂缝重构信号的主特征值数量。而相比于分解得到的1 024个特征值，已经极大减少了计算维数。

从图4可以看出，比较空气裂缝与充水裂缝的冲激响应函数，在早时响应(<3.2 ns)两者的冲激响应函数相似，而在5.38 ns处，充水裂缝的冲激响应回波明显平缓于空气裂缝，即在此处充水裂缝回波体现的频率低于空气裂缝。

图4 空气裂缝与充水裂缝的冲激响应函数Fig.4 Impulse response functions of air crack and water-filled crack

根据SEM算法对冲激响应计算复自然谐振信息，即极点特征，得到极点的衰减因子与谐振频率，如图5所示。

图5 探地雷达回波的极点特征分布Fig.5 Pole feature distribution of GPR echoes

从图中可以看出，对沙地背景进行探测，获得的极点如“□”所示，主要分布在0.2 GHz，0.8～1 GHz，1.3 GHz附近，并且衰减因子的绝对值比较小，靠近0.空气裂缝的极点分布如“○”所示，主要分布在1 GHz以上，即高频分量较多。而充水裂缝的极点分布如“*”所示，其中一部分与空气裂缝中的高频分量类似，而另一部分则出现在0.3～1 GHz之间，相比于背景，出现了较多的高频成分，相比于空气裂缝，出现了较多的低频成分。这一部分极点正是区分空气裂缝与充水裂缝的关键特征信息。

对上述极点特征空间进行模糊C均值聚类，得到的聚类结果如图6所示。

从图6中可以看出，利用模糊C均值聚类方法将极点特征空间分成4类，同时根据背景、空气裂缝、充水裂缝的已知经验信息调整初始聚类中心，每组样本调整优化后的初始聚类中心与其均值如图箭头处所示，对测试样本的极点特征利用图6的初始聚类中心均值作为其初始聚类中心，获得的分类结果可以按照如下规则进行判断：

图6 极点特征聚类Fig.6 Clustering of pole feature

对于极点p，如果满足：

p∈①，p∈②，且p∉③，则属于空气裂缝；

p∈③，则属于充水裂缝。

分别取5组空气裂缝与充水裂缝雷达回波信号的测试样本，按同上方法计算其极点，得到的极点在极点特征空间中的分布如图7所示。

图7 测试样本的极点在极点特征空间中的分布Fig.7 Distribution of poles from test samples in pole feature space

图7中“+”的极点主要分布在①与②区域，并且没有极点分布在③区域，根据上述规则可判定为空气裂缝，而“△”的极点大多分布在①与③区域，根据上述规则，凡是出现在③区域的极点，可判定为充水裂缝的极点。结合给定的样本参数，此判定均正确。

3 小结

针对地质雷达无损检测中回波数据量大，人工解读困难等问题，本文提出基于极点特征聚类的公路隐藏裂缝自动识别方法。首先对回波信号进行奇异值分解，提取体现目标信号的大奇异值，重构回波信号，由此不仅对回波信号进行降维，还去除了噪声干扰。在精度范围允许值内，分析了不同介电目标所需的大奇异值数量。再通过解卷积过程获得目标的冲激响应，利用奇点展开法提取其晚时响应部分的极点信息，构建极点特征空间，并对极点特征进行模糊C均值聚类，为了降低初始聚类中心对聚类结果的影响，对已知环境下的学习样本进行聚类，并根据聚类结果调整优化初始聚类中心，获取调整后的均值，将待识别目标回波的极点投影到极点特征空间按照均值初始化聚类中心进行分类。通过100组学习样本与10组测试样本的实验，验证了算法的有效性。