CN102495402B - 在钢筋强反射回波干扰下的道面灾害目标检测识别方法 - Google Patents

Abstract

一种在钢筋强反射回波干扰下的道面灾害目标识别方法。其主要包括目标检测、识别算法及钢筋回波抑制算法等步骤。本发明提供的在钢筋强反射回波干扰下的道面灾害目标识别方法是在钢筋强反射回波干扰的背景下，首先对探地雷达天线接收到的时空二维回波数据进行预处理，然后进行目标检测，并确定目标的横向位置，再提取每个目标的一道回波数据进行目标识别，从而确定钢筋的横向位置，最后利用波场分离理论将钢筋回波与目标回波分离，并进一步进行抑制，重新检测识别，达到灾害检测识别的目的。本方法具有受噪声、钢筋回波影响小，检测识别率高等优点。

Description

在钢筋强反射回波干扰下的道面灾害目标检测识别方法

技术领域

本发明属于探地雷达信号处理技术领域，特别是涉及一种在钢筋强反射回波干扰下的道面灾害目标检测识别方法。

背景技术

探地雷达(Ground Penetrating Radar，简称GPR)是一种通过向地下发射超高频(10⁶MHz-10⁸MHz)的宽频带窄脉冲电磁波而对地下介质内不可见的目标或界面进行定位的电磁技术。电磁波遇到介质层中电磁特性变化的目标体即发生反射，返回地表时被探地雷达的接收天线接收，并根据其延时、形状及频谱特性等参数解译出目标位置、结构及属性等目标信息。由于探地雷达具有无损、快速、穿透力强及分辨率高等优点，因此逐渐被人们所熟知，被应用于工程勘察，并在许多领域都获得了较成熟的研究成果，如高速公路路面和机场跑道隐性灾害的检测、混凝土建筑结构检测、地下管线和人防工事探测及考古调查等领域。

在探地雷达系统中，为了获得更多反射波的特征，通常使用宽频带接收机记录回波，因此其在接收有用信号的同时，无可避免地会接收到各种干扰和噪声。探地雷达的干扰主要来自以下几个方面：

1、仪器的噪声及系统噪声；

2、来自手机、电台等电磁干扰；

3、汽车、飞机等产生的噪声；

4、来自地下加固钢筋、通信电缆、管线等设施的反射波；

5、来自地下非探测目的体的局部不均匀体的反射波；

一般来说，从频域上即能够分辨出是否存在噪声，并能在频域上抑制掉大部分的噪声能量。后两类干扰是相干干扰，当它们不是探测目的体时，其反射波就变成了干扰波。但对于相干干扰则并非如此，它们在时空域上与感兴趣的目标回波混叠在一起，在频域上也无法分开。甚至某些小目标的回波在能量上远远小于加固钢筋等强反射体的回波，从而导致其回波被湮没在钢筋回波中，因此很难被检测到。总而言之，相干干扰波的存在使得对感兴趣目标的检测识别工作变得十分困难，因此需要尽可能地进行消除。

发明内容

为了解决上述问题，本发明的目的在于提供一种受噪声、钢筋回波影响小，检测识别率高的在钢筋强反射回波干扰下的道面灾害目标识别方法。

为了达到上述目的，本发明提供的在钢筋强反射回波干扰下的道面灾害目标识别方法包括按顺序进行的下列步骤：

(1)对探地雷达时空二维回波数据进行预处理，以去除时空二维回波数据中的直达波分量和噪声；

(2)利用时空二维回波数据进行目标检测，并确定目标的横向位置；

(3)提取每个目标的一道回波数据并进行识别，从而确定钢筋的横向位置；

(4)利用Hyp-curvelet变换将钢筋回波与灾害回波分离，抑制钢筋回波，重构灾害回波，然后重新进行灾害检测识别。

所述的步骤(2)中目标检测的方法是利用Hyp-curvelet变换在一个合适的尺度下对时空二维回波数据进行分解，通过在尺度空间中搜索局部峰值来检测目标及得到目标的横向位置。

所述的步骤(3)中目标识别的方法是首先基于不同种类目标回波时频能量密度分布之间的差异，在时频平面上提取时间方差、频率方差和正负能量时间中心三个特征，并利用该三维特征和最小距离分类器将待识别目标判为最可能的两类目标，然后采用Fisher线性判别分析方法预选出区分该两类目标性能最好的特征子集，最后完成识别，同时得到钢筋的横向位置。

所述的步骤(4)中钢筋回波抑制的方法是首先利用Hyp-curvelet变换将探地雷达时空二维回波信号在多个尺度下进行投影，目标回波的能量将聚集且与其他目标分离，进而结合步骤(3)中得出的钢筋的横向位置，在尺度空间中消除钢筋回波分量，最后将数据反变换回时空域，得到抑制了钢筋回波的数据。

本发明提供的在钢筋强反射回波干扰下的道面灾害目标识别方法是在钢筋强反射回波干扰的背景下，首先对探地雷达天线接收到的时空二维回波数据进行预处理，然后进行目标检测，并确定目标的横向位置，再提取每个目标的一道回波数据进行目标识别，从而确定钢筋的横向位置，最后利用波场分离理论将钢筋回波与目标回波分离，并进一步进行抑制，重新检测识别，达到灾害检测识别的目的。本方法具有受噪声、钢筋回波影响小，检测识别率高等优点。

附图说明

图1为本发明提供的在钢筋强反射回波干扰下的道面灾害目标检测识别方法流程图。

图2(a)和图2(b)分别为本发明所采用的一段探地雷达的时空二维回波数据的几何模型及其B-scan横向堆积图。

图3为使用Hyp-curvelet变换在一个合适的尺度下分解的结果图。

图4为待识别目标数据的时频特征在特征空间的分布图。

图5为钢筋回波抑制后重构的B-scan数据图。

图6(a)和图6(b)分别为使用重构数据重新检测和识别的结果图。

具体实施方式

下面参照附图和具体实施例对本发明提供的在钢筋强反射回波干扰下的道面灾害目标检测识别方法进行详细说明。

为了方便灾害目标识别性能分析且不失一般性，本实施例在机场跑道两层结构中使用发射脉冲波形为Ricker波、中心频率为1GHz、收发分置的探地雷达系统进行灾害检测识别。在其它道面结构下使用其它探地雷达系统，本发明的方法同样适用。

如图1所示，本发明提供的在钢筋强反射回波干扰下的道面灾害目标检测识别方法包括按顺序进行的下列步骤：

(1)对探地雷达(GPR)时空二维回波(B-scan)数据进行预处理，以去除时空二维回波数据中的直达波分量和噪声，以便后续处理；

(2)利用时空二维回波数据进行目标检测，并确定目标的横向位置；

探地雷达沿着测试线扫描时，天线与目标的相对距离不断变化而导致目标回波信号在时空域呈现近似双曲线的形状，我们认为可通过在时空域中检测双曲线来实现目标检测。

Hyp-curvelet变换实质上是双曲线Radon变换切片上的一维小波变换，它延续了双曲线Radon变换在描述双曲线奇性上的优越性能，其定义为：

W_γ(b)＝ψ_a*R_hypf(τ，b)

      ＝<f，ψ_γ>＝∫∫f(t，x)ψ_γdtdx                    (1)

其中，“*”表示卷积运算；

为双曲线Radon变换结果；f为GPR B-scan数据；综合考虑GPR信号的时空域特征和衰减特性后，Hyp-curvelet小波ψ_γ可表示为：

${\mathrm{\&psi;}}_{\mathrm{\&gamma;}}(t,x)=a^{-\frac{1}{2}}\mathrm{\&psi;}\left(\frac{g(t,x)-b}{a}\right)\exp [-a_0{\left(\frac{x-b_x}{a_x}\right)}^2]$

$={\left(a_x{a}_t\right)}^{-\frac{1}{2}}\{1-{[\frac{t-b_t}{a_t}-\sqrt{{\left(\frac{x-b_x}{a_x}\right)}^2+{\mathrm{\&tau;}}_0^2}]}^2\}\exp \{-\frac{1}{2}{[\frac{t-b_t}{a_t}-\sqrt{{\left(\frac{x-b_x}{a_x}\right)}^2+{\mathrm{\&tau;}}_0^2}]}^2-a_0{\left(\frac{x-b_x}{a_x}\right)}^2\}---\left(2\right)$

式中a₀为空间方向上的衰减因子；τ₀是一个正常数；a_t和a_x分别为时间和空间上的尺度参数；b_t和b_x分别为时间和空间上的位移参数。

综上所述，Hyp-curvelet变换在合适的尺度下都能将目标回波信号的能量聚集在双曲线的顶点处，可利用Hyp-curvelet变换将探地雷达接收的时空二维数据投影到一个尺度空间中，并在这个尺度空间搜索局部峰值来进行目标检测。我们将框架γ＝{a_t，a_x，b_t，b_x}离散化为：

γ∈Γ＝{2^J，a_x0，m，n}，J＝3，m＝1，2，…，N_t，n＝1，2，…，N_x    (3)

a_x0可由下式计算得出：

$a_{x0}=\frac{\mathrm{vdt}}{2\mathrm{dx}}=\frac{\mathrm{cdt}}{2\sqrt{\mathrm{\&epsiv;}}\mathrm{dx}}---\left(4\right)$

其中v为电磁波在介质中传播的速度；dt为雷达时间采样间隔；dx为雷达空间采样间隔；c为光速，ε为介质的相对介电常数。

(3)提取每个目标的一道回波数据(A-scan)并进行识别，从而确定钢筋的横向位置；

通过步骤(2)检测到的并不一定就是我们感兴趣的目标，因此需要进一步识别，尤其是识别钢筋回波。此时我们根据步骤(2)得到的目标横向位置提取目标的一道回波数据s(t)进行具体分析。不同目标回波的时频分布是不一样的，对于大小不同的目标，其回波能量在时间上分布不同，如裂缝主要是纵向分布，其上、下表面反射回波分离且能量相差不大，因此其回波能量在时间上分布范围较大；脱空主要是水平向分布，其上、下表面反射回波混叠，因此其回波能量在时间上分布范围较小；而钢筋则因为电磁波对其的穿透能力较弱，其反射回波的能量主要聚集于上表面，且多次波的能量远小于一次回波的能量，因此其回波能量在时间上分布的范围也较小。同时目标电磁特性之间的差异会导致其反射回波的频谱密度不同，因此频率上的特征主要表征的是目标在材质上的差别。于是我们在时频平面上提取三个特征，即时间方差、频率方差和正负时间中心之差，分别定义为：

${\mathrm{\&sigma;}}_t^2=\frac{1}{E}\underset{t,f}{\&Integral;\&Integral;}{(t-t_0)}^2{\mathrm{WD}}^2(t,f)\mathrm{dtdf}---\left(5\right)$

${\mathrm{\&sigma;}}_f^2=\frac{1}{E}\underset{t,f}{\&Integral;\&Integral;}{(f-f_0)}^2{\mathrm{WD}}^2(t,f)\mathrm{dtdf}---\left(6\right)$

$\mathrm{dt}=t_{0+}-t_{0-}=\frac{1}{E_{+}}\underset{t,f}{\&Integral;\&Integral;}t{\mathrm{WD}}_{+}^2(t,f)\mathrm{dtdf}-\frac{1}{E_{-}}\underset{t,f}{\&Integral;\&Integral;}t{\mathrm{WD}}_{-}^2(t,f)\mathrm{dtdf}---\left(7\right)$

式中WD为s(t)的Wigner变换；

E为回波总能量，表示为：

t₀为回波能量的时间中心，表示为：

f₀为回波能量的频率中心，表示为：

WD₊(t，f)为WD(t，f)的正能量部分，即： ${\mathrm{WD}}_{+}(t,f)=\left\{\begin{array}{cc}\widetilde{\mathrm{WD}}(t,f)& \widetilde{\mathrm{WD}}(t,f)>0\\ 0,& \widetilde{\mathrm{WD}}(t,f)<0\end{array}\right.;$

WD_-(t，f)为WD(t，f)的负能量部分；

E₊和E_-分别为正、负总能量。

特征的区分能力主要体现在不同种类目标之间的类间距是否足够大，而同种类目标之间的类内距是否足够小。我们定义第i类目标与第j类目标之间的距离为：

$D_{i,j}=\sqrt{({\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_i-{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j){\left(\stackrel{\&OverBar;}{\mathrm{\&Lambda;}}\right)}^{-1}{({\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_i-{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j)}^T,}i=\mathrm{1,2},...,p---\left(8\right)$

其中p为总目标的类数；为第i类目标特征向量的均值，表示此类目标在特征空间的中心；

为第i类目标特征的协方差矩阵，表示此类目标在特征空间的大小；为目标特征协方差矩阵的均值，表示所有类目标在特征空间中大小的均值。

在机场跑道的典型结构模型下对五类目标(裂缝，气脱空，水脱空，钢筋和小气孔)进行仿真实验，并建立了训练样本空间。选取特征向量δ＝[lgσ_t，lgσ_f，sgn(dt)lg|dt|]，各类目标之间的类间距及其类内距如表1所示。

表1三维特征空间中目标之间的类间距及其类内距

	裂缝	气脱空	水脱空	钢筋	小气孔
						裂缝	2.8473
气脱空	15.5708	0.6464
						水脱空	16.5786	5.6816	1.9762
钢筋	13.463	9.6556	6.2738	0.8863
						小气孔	10.9595	4.7206	6.7765	7.7653	1.8594

从表1的内容可以看出，各类目标的类内距远小于类间距，因此所提取的三个特征具有良好的分类性能。其中气脱空与钢筋的类内距较小，因此其聚类效果最好。而气脱空与小气孔由于属性相同，且其纵向尺寸相近，因此其类间距相对较小。同时钢筋与水脱空对于背景媒质来说均属于光密媒质，且能量都聚集在小段时间范围中，因此其类间距相对较小。

Fisher线性判别分析是一种通过类间方差与类内方差来寻找一个有效的分类方向的方法，对于具体的两类问题Fisher判别函数定义如下：

$J\left({\mathrm{\&delta;}}_k\right)=\frac{{|{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_{i,k}-{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_{j,k}|}^2}{{\mathrm{\&Lambda;}}_{i,k}+{\mathrm{\&Lambda;}}_{j,k}},k=1,...,K---\left(9\right)$

其中

和 ${\mathrm{\&Lambda;}}_{i,k}=\frac{1}{M-1}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{\left({\mathrm{\&delta;}}_{i,k}\right[m]-{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_{i,k})\&CenterDot;\left({\mathrm{\&delta;}}_{i,k}\right[m]-{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_{i,k})}^T$ 分别为第i类目标的第k个特征的均值与方差，M为样本数，K为特征数，J(δ_k)称为特征δ_k的Fisher比，特征在训练样本集上的Fisher比越大说明该特征的类别区分度越好，即包含越多的鉴别信息。针对具体的两类问题，可根据Fisher判别函数筛选出类别区分度好的特征子集δ_sub。我们定义平均Fisher比为：

$\stackrel{\&OverBar;}{J}=\frac{1}{m}\underset{k}{\mathrm{\&Sigma;}}J\left({\mathrm{\&delta;}}_k\right)---\left(10\right)$

如果

则将δ_k纳入此两类识别问题的特征子集δ_sub。

由于所提特征具有良好的区分能力，因此可采用最简单的最小距离分类器来进行判别。但是待识别目标的特征分量正处在两类目标之间时，直接使用最小距离分类器进行判决容易发生误判，此时需要对识别后的目标进行再识别。使用Mahalanobis Distance作为待识别的目标与各类目标在特征空间中的距离，即：

$S(X,Y_i)={\left[{(X-Y_i)}^T{\mathrm{\&Lambda;}}_i^{-1}(X-Y_i)\right]}^{\frac{1}{2}},i=\mathrm{1,2},...,p---\left(11\right)$

其中X为待识别目标的特征向量，Y_i为第i类目标的特征向量，p为总目标类别数。我们定义判别函数为：

$\left\{\begin{array}{c}{k}_1=\underset{i}{\arg \min }S(X,Y_i)\\ k_2=\underset{i\&NotEqual;k_2}{\arg \min }S(X,Y_i)\end{array}\right.---\left(12\right)$

即认为待识别的目标可能是与其最接近的一类或次接近的一类。

然后利用由Fisher线性判别方法预选出在上述两类目标上鉴别性能最强的特征子集δ_sub并进行再次判决，即：

$k=\underset{i=k_1,k_2}{\arg \min }S(\stackrel{\&OverBar;}{X},{\stackrel{\&OverBar;}{Y}}_i)---\left(13\right)$

和

分别为待识别目标特性向量X和第i类目标特征向量Y_i在特征子集δ_sub上所对应的特征向量。

(4)利用Hyp-curvelet变换将钢筋回波与灾害回波分离，抑制钢筋回波，重构灾害回波，然后重新进行灾害检测识别；

钢筋回波与灾害目标回波在时空域的特征相似，均呈现近似双曲线特征。同时当目标相距不够远时，其回波在时空域上会发生混叠，因而不能直接在时空域消除钢筋回波。我们认为钢筋的干扰问题用以下方法来解决，即首先通过Hyp-curvelet变换将探地雷达时空二维回波数据投影到多个尺度空间中，滤除所有尺度下的钢筋回波分量后再将数据反变换回时空域，可得到只含有灾害目标回波的数据，最后重新检测识别。

Hyp-curvelet不是正交基而是一个框架(ψ_γ)_γ∈Γ，如式(3.7)所示，观测信号不仅能够通过框架系数<f(t，x)，ψ_γ>_γ∈Γ来刻画，而且还能从这个系数中重构出来，即：

其中，

为(ψ_γ)_γ∈Γ的对偶框架：

其中F为(ψ_γ)_γ∈Γ的框架算子

是H到 $l^2\left(J\right)=\{c={\left(c_j\right)}_{j\&Element;J};{\left|\right|c\left|\right|}^2=\underset{j\&Element;J}{\mathrm{\&Sigma;}}{\left|c_j\right|}^2<\&infin;\}$ 的线性算子。

对于Lf＝w的求解问题，我们可采用共轭梯度法进行计算，这是由于共轭梯度法具有指数级的收敛速度，且不需要知道框架的边界。框架的上下边界越接近，迭代收敛速度越快。框架边界不仅与小波的性能有关，而且与尺度参数及位移参数的采样间隔有关。因此在Hyp-curvelet变换中参数设置尤其重要，这影响着框架是否是紧框架，从而影响重构的结果精度以及迭代时的收敛速度。事实上，在尺度空间中抑制钢筋回波的同时无可避免地会损失一些有用信息及残留一些钢筋回波分量。由于框架元素之间是线性相关的，因此Hyp-curvelet变换能带来大量的冗余信息。尺度空间中的相关信息可互相补偿，因此部分重构也能有效地逼近灾害目标回波，从而达到了抑制钢筋回波的目的。综合考虑后，我们将框架γ＝{a_t，a_x，b_t，b_x}离散为：

$\mathrm{\&gamma;}\&Element;\mathrm{\&Gamma;}=\{2^J,a_{x_0},m,n\},J=\mathrm{0,1},...,N_J,m=\mathrm{1,2},...,N_t,n=\mathrm{1,2},...,N_x---\left(17\right)$

即时间尺度参数按幂级数离散，空间尺度参数固定为一个与电磁波在介质中传播速度有关的常数，位移参数的离散采用均匀采样，按GPR数据实际采样间隔采样。

本实施例中使用高保真的仿真软件GPRmax生成的时空二维回波信号，其中参数设置如下：收发天线间隔为10cm，离地面4cm，发射脉冲波形为Ricker波，中心频率为1GHz，时间窗为12ns，时间与空间采样点数分别为2036和250。地下场景为机场跑道两层结构，其中第一层为34cm的混凝土面层，其介电常数与电导率分别为9和0.02；第二层为20cm的碎石基层，其介电常数与电导率分别为12和0.025。该段B-scan数据中包含有4根钢筋和2条裂缝，其几何模型如图2(a)所示，其B-scan横向堆积图如图2(b)所示，信噪比为10dB。

图3为使用Hyp-curvelet变换在一个尺度下分解的结果图。可以看到双曲线的能量被累积起来而形成了局部峰值，而噪声的能量聚集效果不明显，从而达到了目标检测的目的。

图4为待识别目标数据的时频特征在特征空间的分布图。可以看到4个目标均与钢筋样本接近。本发明第一次识别时将目标判为最可能的两类目标，该4个目标均被判为钢筋与裂缝，利用Fisher线性分析方法预选出的区分钢筋与裂缝性能最好的特征子集进行分类后，所有目标均被判为钢筋。

图5为钢筋回波抑制后重构的B-scan数据图。可以看到随机噪声、钢筋回波及其多次波大部分被抑制了，而裂缝回波虽然稍有损伤，但微小的差别不会影响目标检测与识别。图5中给出的结果说明本发明能有效地抑制噪声和钢筋回波。

图6为使用重构数据重新检测和识别的结果图。如6(a)所示，裂缝上表面的回波和裂缝与第二层界面之间的反射回波均形成了局部峰值。图6(b)为待识别目标数据的特征在特征空间的分布，可以看到2个目标均与裂缝样本接近，同时目标被正确识别为裂缝。

Claims (1)

1.一种在钢筋强反射回波干扰下的道面灾害目标识别方法，所述的在钢筋强反射回波干扰下的道面灾害目标识别方法包括按顺序进行的下列步骤：

(1)对探地雷达时空二维回波数据进行预处理，以去除时空二维回波数据中的直达波分量和噪声；

(2)利用时空二维回波数据进行目标检测，并确定目标的横向位置；

(3)提取每个目标的一道回波数据并进行识别，从而确定钢筋的横向位置；

(4)利用Hyp-curvelet变换将钢筋回波与灾害回波分离，抑制钢筋回波，重构灾害回波，然后重新进行灾害检测识别；其特征在于：

所述的步骤(3)中目标识别的方法是首先基于不同种类目标回波时频能量密度分布之间的差异，在时频平面上提取时间方差、频率方差和正负能量时间中心三个特征，并利用该三维特征和最小距离分类器将待识别目标判为最可能的两类目标，然后采用Fisher线性判别分析方法预选出区分该两类目标性能最好的特征子集，最后完成识别，同时得到钢筋的横向位置。

Images