一种钢筋回波干扰下的灾害目标成像方法与流程

本发明属于探地雷达信号处理技术领域，特别是涉及一种钢筋回波干扰下的灾害目标成像方法。

背景技术

探地雷达(groundpenetratingradar,gpr)是一种通过利用地下介质的不连续性来探测地下目标的有效工具，其是由置于地表或地上的发射天线向地下发射一超宽带的电磁波，电磁波在地下传播过程中遇到不同电性介质界面时会发生反射和折射，反射回来的电磁波由接收天线接收。根据接收的雷达回波波形、振幅和双程时延等参数便可以推断出地下目标体的空间位置、结构、电磁特性及几何形态，从而达到探测地下结构的目的。由于探地雷达具有探测速度快、探测过程连续、操作方便灵活、分辨率高、不损坏被探测目标等特点，因此使其在国防、公路、铁路、机场建设、探矿、考古等方面得到了广泛应用。

在不影响道路正常运行并且不破坏道面原有结构的基础上，采用高科技的无损检测技术是必然趋势。开展道面无损检测技术研究，将在控制道面施工质量、深入认识道面长期使用性能、改善道面设计、优化道面改建和扩建方案及提高道面养护水平等方面具有重要的意义。用探地雷达进行道面检测操作简单,结果直观精确，且能适应规模大、无破损、速度快、精度高的工程质量检测要求。

在道面中加入钢筋，可以提高道面的耐久性，降低维护费用，与无筋混凝土道面相比具有较好的经济性；由于钢筋的存在，减少了裂缝的数量，并使裂缝不能张开，从而可以增大板的长度和相应地减少接缝的数量，改善了道面的使用品质。

道面初期常见的隐形灾害主要有脱空和裂缝。钢筋的强反射回波会严重影响甚至淹没脱空和裂缝的反射回波，并且钢筋与道面下各层介面的多次反射回波也会严重影响灾害目标的检测，所以有效地抑制钢筋的强反射回波是十分关键的。

探地雷达成像是最直接和最有效的探地雷达目标检测与识别技术，它可以更直观地向非技术人员显示地下目标的相关信息。常用的成像算法是通过对散射场数据进行聚焦成像，其可以得到地下目标的位置及形状信息。

钢筋回波干扰下的灾害目标成像方法既要考虑如何对钢筋回波进行抑制，保留有效的灾害目标回波，又要考虑地下多层介质的特点，研究适用范围广的钢筋回波干扰下的灾害目标成像方法。

技术实现要素：

为了解决上述问题，本发明的目的在于提供一种直观、处理步骤简单且处理速度快及运算量小的钢筋回波干扰下的灾害目标成像方法。

为了达到上述目的，本发明提供的钢筋回波干扰下的灾害目标成像方法包括按顺序进行的下列步骤：

(1)对探地雷达接收的回波信号进行预处理，提取出有效的钢筋及灾害目标反射回波信号；

(2)将上述预处理后的回波信号进行fhyp-curvelet分解，使探地雷达接收的回波信号投影到多个尺度空间，得到每个尺度下的回波信号；

(3)将步骤(2)得到的回波信号中的钢筋回波信号分量进行滤除；

(4)将步骤(3)得到的滤除钢筋回波信号分量的回波信号进行fhyp-curvelet重构，由此将该回波信号变换到时空域，获得时空域回波信号；

(5)将步骤(4)得到的时空域回波信号进行频域—波数域的二维快速傅里叶变换而得到频域—波数域的数据；

(6)在步骤(5)得到的频域—波数域的数据基础上，考虑多层介质波速变化，得到波数域下探地雷达的回波信号；

(7)将步骤(6)得到的波数域下的回波信号在波数域下进行二维快速逆傅里叶变换，最终得到探地雷达成像结果。

在步骤(1)中，所述的对探地雷达接收的回波信号进行预处理，提取出有效的钢筋及灾害目标反射回波信号的方法是：利用非一致性检测选取不含有目标的背景回波，将探地雷达接收的回波信号减去背景回波，即可提取出效的钢筋及灾害目标反射回波信号。

在步骤(2)中，所述的将上述预处理后的回波信号进行fhyp-curvelet分解，使探地雷达接收的回波信号投影到多个尺度空间，得到每个尺度下的回波信号的方法是：综合考虑探地雷达发射信号波形及衰减特性和时空特征后，建立fhyp-curvelet分解的小波函数，利用二维快速傅里叶变换进行双曲线积分，将沿双曲线分布的能量在多个尺度空间下聚集成点。

在步骤(3)中，所述的将步骤(2)得到的回波信号中的钢筋回波信号分量进行滤除的方法是：根据钢筋与灾害目标的密度不同且与道面地下媒质的密度也有差异的特点，在多个尺度空间下将钢筋回波分量滤除，保留灾害目标回波分量。

在步骤(4)中，将步骤(3)得到的滤除钢筋回波信号分量的回波信号进行fhyp-curvelet重构，由此将该回波信号变换到时空域，获得时空域回波信号的方法是：基于框架理论，利用二维快速逆傅里叶变换将fhyp-curvelet分解后的某一尺度下的信号转换到时空域，最后将多个尺度下的重构信号进行叠加。

在步骤(5)中，所述的将步骤(4)得到的时空域回波信号进行频域—波数域的二维快速傅里叶变换而得到频域—波数域的数据的方法是：将信号进行频域—波数域的二维快速傅里叶变换的方法是：将步骤(4)得到的时空域信号做水平位置和时间的二维快速傅里叶变换，得到水平波数域—频域—深度所表示的信号形式。

在步骤(6)中，所述的在步骤(5)得到的频域—波数域的数据基础上，考虑多层介质波速变化，得到波数域下探地雷达的回波信号的方法是：考虑多层介质波速变化的特点，建立深度波数域、水平波数域、频域及变化的波速之间的关系，将信号变换为考虑了波速变化的波数域下的雷达信号。

在步骤(7)中，所述的将步骤(6)得到的波数域下的回波信号进行波数域下的二维快速逆傅里叶变换，最终得到探地雷达成像结果的方法是：将步骤(6)得到的波数域下的雷达信号进行二维快速逆傅里叶变换，最终得到雷达成像结果。

本发明提供的钢筋回波干扰下的灾害目标成像方法能在灾害目标回波信号被钢筋强反射多次回波信号淹没的情况下，有效抑制钢筋回波对灾害目标的影响，且考虑电磁波在多层介质下传播波速变化的特点，整个方法利用快速傅里叶变换进行不同域下的信号变换，处理速度快，计算量小。

附图说明

图1为本发明提供的钢筋回波干扰下的灾害目标成像方法流程图。

图2为本发明中探地雷达实际接收到的二维b-scan回波信号图。

图3为本发明中利用非一致性检测滤除背景后的目标反射回波信号图。

图4(a)—(d)为本发明中fhyp-curvelet分解在四个尺度空间中的投影结果。

图5(a)—(d)为本发明中四个尺度空间下钢筋滤波后的结果。

图6为本发明中fhyp-curvelet重构后的空时域灾害目标回波信号。

图7为本发明中的灾害目标成像结果。

图8为未进行钢筋抑制的成像结果。

具体实施方式

为能进一步了解本发明的发明内容、特点及功效，下面参照附图和具体实施例对本发明提供的钢筋回波干扰下的灾害目标成像方法进行详细说明。

如图1所示，本发明提供的钢筋回波干扰下的灾害目标成像方法包括按顺序进行的下列步骤：

(1)对探地雷达接收的回波信号进行预处理，提取出有效的钢筋及灾害目标反射回波信号；

根据不同的需要，探地雷达系统可采用不同的扫描方式采集数据，主要有三种不同的扫描方式，分别称之为a-scan、b-scan、c-scan，扫描得到的分别是一维的a-scan回波信号、二维的b-scan回波信号、三维的c-scan回波信号。本发明主要对b-scan回波信号进行处理。

首先去除探地雷达接收的b-scan回波信号中收发天线间的直接耦合波：将雷达天线对准开阔无物体的空间，或置于暗室内，录取对空信号，此信号可以近似为直接耦合信号，再将探地雷达接收的b-scan回波信号减去此直接耦合信号即有效去除了直接耦合波，获得b-scan回波信号。

然后去除道面下各层媒质表面反射回波：先对上述已去除直接耦合波的b-scan回波信号在时间维进行非一致性检测，检测到的非一致性部分的回波信号即为有目标反射信号，将其滤除，一致性部分的回波信号即为各层媒质表面反射回波信号，选取此部分a-scan回波信号作为背景回波，将每一道探地雷达接收的a-scan回波信号减去背景回波即可以有效去除各层媒质表面反射回波信号，由此提取出有效的钢筋及灾害目标反射回波信号；

(2)将上述预处理后的回波信号进行fhyp-curvelet分解，使探地雷达接收的回波信号投影到多个尺度空间，得到每个尺度下的回波信号；

奥地利数学家radon于1917年提出了radon变换，其基本思想是对被积函数在给定路径上进行积分，根据被积函数的积分路径可以分为线性radon变换和非线性radon变换。由于探地雷达接收的灾害目标和钢筋的回波信号在时空域呈现近似双曲线特征，这里给出双曲线radon变换的定义：

式中f(t,x)为gpr接收的回波信号，τ为时间截距，q为双曲线曲率参数。双曲线radon变换的实质是将gpr接收的回波信号沿着双曲线进行积分，并将沿双曲线分布的能量在radon域中聚集成一个点。

fhyp-curvelet分解实质上是双曲线radon变换切片上的一维小波变换，fhyp-curvelet分解是利用快速傅里叶变换将探地雷达接收的回波信号投影到多个尺度空间，那么尺度a下的fhyp-curvelet分解可以表示为：

其中wa(bt,bx)为分解系数，ψa(t,x)为hyp-curvelet小波，φa(u,v)为hyp-curvelet小波ψa(t,x)的二维离散傅里叶变换，if(u,v)为gpr接收的回波信号f(t,x)的二维逆离散傅里叶变换，at和bt分别为时间上的尺度参数及位移参数，ax和bx分别为空间上的尺度参数及位移参数。

如式(2)所示，将gpr接收的回波信号f(t,x)投影到k个尺度空间中可分为以下几步：

1)对gpr接收的回波信号f(t,x)进行二维逆fft变换得到if(u,v)，并令a＝0；

2)在尺度a下，对hyp-curvelet小波ψa(t,x)进行二维fft变换得到φa(u,v)，并求得y(u,v)＝φa(u,v)if(u,v)，对y(u,v)进行二维fft变换得到每个尺度下的分解系数wa(bt,bx)；

3)若所有尺度都进行过计算，则结束分解；反之，令a＝a+1，转到步骤2)。

根据上式可以将gpr接收的回波信号f(t,x)投影到k个需要的尺度空间，需要经过2k+1次二维fft或ifft变换和k次n×m维矩阵点乘。

(3)将步骤(2)得到的回波信号中的钢筋回波信号分量进行滤除；

根据钢筋与灾害目标的密度不同且与道面地下媒质的密度也有差异的特点，我们可以根据光疏媒质与光密媒质的信号传输特性区分钢筋回波分量与灾害目标回波分量。

钢筋的密度比机场道面介质的密度大，属于光密媒质，灾害目标填充的多为空气，属于光疏媒质。探地雷达接收的光疏媒质与光密媒质的回波信号初始相位相反，且钢筋回波分量强度远远大于灾害目标，所以可以根据这两点识别出钢筋回波信号分量，将钢筋回波信号分量在多个尺度空间下滤除。

(4)将步骤(3)得到的滤除钢筋回波信号分量的回波信号进行fhyp-curvelet重构，由此将该回波信号变换到时空域，获得时空域回波信号；

将存在钢筋干扰的gpr接收的回波信号投影变换到一个新域，在该域下的多个尺度空间中，不同目标的回波信号是相互分离的，将钢筋回波信号分量滤除后，需要把含有灾害目标的有效回波信号变换到时空域，进行fhyp-curvelet重构。

根据框架理论，fhyp-curvelet利用快速傅里叶变换对尺度a下的信号进行重构，可以表示为：

其中为的二维离散傅里叶变换，wa(u,v)为分解系数wa(bt,bx)的二维离散傅里叶变换。iza(t,x)为的二维离散逆傅里叶变换。

利用在k个尺度下的分解系数重构f可分为以下几步：

1)分别对尺度a下的分解系数wa和进行二维离散傅里叶变换得到wa和

2)计算对za进行二维离散逆傅里叶变换得到iza(t,x)；

3)将k个尺度下得到的iza(t,x)进行累加；

根据式(3)，可以将滤除钢筋回波信号分量后的k个尺度空间的回波信号进行重构，得到时空域回波信号f(t,x)，同样需要经过2k+1次二维fft或ifft变换和k次n×m维矩阵点乘。

(5)将步骤(4)得到的时空域回波信号进行频域—波数域的二维快速傅里叶变换而得到频域—波数域的数据；

假设在介质中的传播速度是真实速度的一半，电磁波产生于目标散射点，向探地雷达天线传播，由探地雷达的记录剖面f(x,z＝0,t)延拓得到剖面f(x,z,t＝0)，从而确定出目标的位置，其中x为时空域中雷达数据的水平位置，z为深度，t为传播时间。

取f(x,z＝0,t)关于x和t的二维傅立叶变换，得到：

由式(4)可以得到探地雷达接收的回波信号在频域—波数域的数据，其中kx为时空域中探地雷达数据的水平位置x在波数域下的傅里叶变换。

(6)在步骤(5)得到的频域—波数域的数据基础上，考虑多层介质波速变化，得到波数域下探地雷达的回波信号；

在深度z处的波场可以表示为：

其中kz为深度z在波数域下的傅里叶变换。

由于电磁波传播速度采用实际速度v的一半，因此：

因而在单一介质中实际速度v是常数。

假设剖面f(x,z,t)是f(kx,z,ω)关于kx和ω的二维逆傅立叶变换，则：

所以:

上述方法没有考虑地下多层媒质的波速变化情况，现将f(x,z,t＝0)表示为另外一种形式，以推导波速变化的情况：

将式带入式(9)得:

该式中kx的积分是一般意义的逆傅里叶变换，但是频率ω的积分不是，由式(10)可得：

对f(kx,z,t＝0)做z的傅里叶变换，则可得f(kx,kz)：

将式(11)带入式(12)得：

其中

令

则a(kx,kz,ω)是a(kx,z,ω)对于z的傅里叶变换，将式(15)表示为波速变化的形式：

由式(16)推得的波数域下的回波信号f(kx,kz)包含了波速变化的因素：

(7)将步骤(6)得到的波数域下的回波信号在波数域下进行二维快速逆傅里叶变换，最终得到探地雷达成像结果。

将波数域下的回波信号f(kx,kz)进行二维快速逆傅里叶变换，得到探地雷达的成像公式可以表示为：

即可以适用于实际中电磁波波速变化的情况。

本发明的效果可以通过以下仿真结果进一步说明。

仿真数据描述：实验中所有数据均由gprmax2.0仿真软件产生，本实验采用双基地探地雷达，天线间距为0.1m，距地面高度为0.04m，发射信号采用ricker波脉冲形式，中心频率为1000mhz，采集数据时间窗为12ns，共采集250道数据。仿真模型为一般机场场道的面层和上基层，层厚从上到下依次为0.34m,0.20m，介电常数依次为9,12。本实验中所采用的机场道面仿真数据中含有四根钢筋，一条裂缝，钢筋均匀分布在面层中距地表面15cm处，每条钢筋间距为50cm，裂缝长23cm，宽7mm。

图2为本发明中探地雷达实际接收到的二维b-scan回波信号图。如图2所示，由于地表反射回波信号很强，淹没了地下目标的回波信号。图3为利用非一致性检测滤除背景后的目标反射回波信号图。如图3所示，背景的两层媒质表面反射回波信号被有效滤除，凸显出地下目标的反射回波信号，但是由于道面浅层钢筋反射回波信号很强，严重影响了灾害目标的回波信号，也无法有效地检测出灾害目标。

图4为fhyp-curvelet分解在四个尺度空间中的投影结果。如图4所示，在四个不同的尺度空间中，仍无法观测到灾害目标的回波分量，但是钢筋回波分量彼此分离且聚集成点，在该四个尺度空间中对钢筋数据进行滤波。图5为四个尺度空间中钢筋滤波后的结果。如图5所示，在四个尺度空间中仍残留了部分钢筋回波的分量，但是可以较明显地观测到灾害目标的回波分量。图6为fhyp-curvelet重构后的空时域灾害目标回波信号。如图6所示，灾害目标回波被重构的较为理想，钢筋回波分量稍有残留。图7为本发明中的灾害目标成像结果，由图7所示，在空域中，灾害目标的双曲线回波被很好地聚焦，可以较好地显示出灾害目标的空间位置。图8为本发明中未进行钢筋抑制的成像结果。如图8所示，成像结果能明显观测到钢筋的空间位置，而无法观测到灾害目标的空间位置。