CN109283524A - 一种用于提高地质雷达信号分辨率的方法 - Google Patents

Abstract

本公开揭示了一种用于提高地质雷达信号分辨率的信号处理方法，所述方法在连续小波变换域内自适应的拓展地质雷达信号的频率，以得到高分辨率的拓频信号。首先采用基于地质雷达信号的频谱分析，确定各拓频频段和相应的加权系数，然后对地质雷达信号进行连续小波变换，得到相应的时间‑尺度(频率)分布，根据频率和小波尺度的换算关系，对相应尺度范围内的小波变换系数进行加权，最后基于加权的连续小波变换，重构对应的拓频信号，达到提高地质雷达信号分辨率的目的。所述方法能够为隧道衬砌和铁路路基检测提供高分辨率地质雷达数据，从而增加后续检测的准确性。

Description

一种用于提高地质雷达信号分辨率的方法

技术领域

本发明属于地质雷达领域，特别涉及一种用于提高地质雷达信号分辨率的高分辨处理方法。

背景技术

利用地质雷达(Ground-penetrating Radar，GPR)进行探测时，发射天线发射电磁波在介质中传播，经过存在电性差异的介质体或者界面时会产生反射信号，该反射信号被设置于检测表层的接收天线接收。通过分析该接收信号实现对介质内部分布规律的分析。地质雷达探测具有超浅层勘探的独特优势，是一种快速的高分辨率无损检测方法，在隧道以及铁路路基检测中具有广泛的应用。但是，受到介质耗散等因素的影响，电磁波在传播过程中会出现频率降低，幅度衰减等问题，使得接收信号的分辨率降低，从而影响后续的异常检测。目前用于提高地质雷达信号分辨率的方法主要有

现有技术1：谱白化类方法

(1)计算每个A-scan的一维傅里叶变换，得到信号的振幅谱和相位谱；

(2)在信号的有效频宽范围内，设定若干个带通滤波器，将幅度谱分为若干个子带；

(3)计算每个子带内的信号振幅谱加权系数，将每个子带内的信号谱乘以对应的子带加权系数，得到拓频后的频谱；

(4)基于拓频后的频谱，计算其一维傅里叶反变换，得到拓频后的A-scan信号；

该类方法还可以在时间域进行能量均衡。

现有技术1的缺点：

1)该类方法在提高地质雷达数据分辨率的同时，也放大了信号中的随机噪声；

2)人为设定各拓频子带宽度，没有统一的设置规则，易造成不同设置参数的拓频效果差异较大；

现有技术2：反褶积方法

(1)提取地质雷达信号子波表达式；

对每一个反射时间点进行如下操作：

(2)基于子波表达式，构造自相关矩阵；

(3)对自相关矩阵进行奇异值分解，并计算加权系数向量；

(4)计算加权系数向量与观测时间序列的内积得到该时间点的反射强度；

上述为时间域反褶积处理流程，还可以根据卷积定理，将时间域反褶积流程推广到频率域进行反褶积处理。

现有技术2的缺点：

1)地质雷达子波信号难以准确提取，因此影响后续的反射强度计算；

2)通常地质雷达数据中含有较强的随机噪声，这也使得反褶积效果受到影响。

发明内容

基于此，本公开揭示了一种用于提高地质雷达信号分辨率的信号处理方法，其特征在于，该方法包括以下步骤：

S100、对多道A-scan信号进行一维傅里叶变换，并在A-scan信号的有效频带宽度内确定A-scan信号的基础频段和各拓频频段；

S200、计算A-scan信号各拓频频段的加权系数；

S300、对A-scan信号进行连续小波变换，得到对应的时间-尺度域小波变换系数(时频分布)；

S400、根据频率和尺度的换算关系，将不同时间-尺度范围内的小波变换系数乘以对应的加权系数；

S500、计算乘以对应的加权系数后小波变换系数对应的小波反变换，得到该A-scan信号对应的高分辨率信号；

S600、对每个A-scan信号重复步骤S300-S500。

本公开具有如下技术效果：

1)本公开的方法不需要估计地质雷达信号子波；

2)本公开参数设置简单，根据信号频谱分布自适应确定频谱划分和加权系数；

3)本公开的连续小波域信号自适应拓频方法可以有效的提高地质雷达剖面的时间分辨率以及增强深层弱反射的强度；

4)通过对频带的有效划分和加权处理，本公开所述方法对随机噪声有很好的抑制作用。

附图说明：

图1为本公开一个实施例中的GPR数据自适应频谱拓展流程；

图2为本公开一个实施例中的频谱划分示意；

图3(A)-图3(F)为本公开一个实施例的处理结果：图3(A)为隧道衬砌地质雷达数据B-scan剖面；图3(B)为连续小波域自适应频谱拓展结果；图3(C)为抽取的水平距离x＝1.44m的A-scan信号及其拓频结果；图3(D)为x＝1.44m处A-scan信号及其拓频结果的傅里叶频谱分析；图3(E)为抽取的水平距离x＝4.83m的A-scan信号及其拓频结果；图3(F)为x＝4.83m处A-scan信号及其拓频结果的傅里叶频谱分析；

图4(A)-图4(F)为本公开一个实施例中的处理结果：图4(A)铁路路基地质雷达数据B-scan剖面；图4(B)连续小波域自适应频谱拓展结果；图4(C)抽取的水平距离x＝2m的A-scan信号及其拓频结果；图4(D)x＝2m处A-scan信号及其拓频结果的傅里叶频谱分析；图4(E) 抽取的水平距离x＝4m的A-scan信号及其拓频结果；图4(F)x＝4m 处A-scan信号及其拓频结果的傅里叶频谱分析。

具体实施方式：

下面结合附图1至附图4(F)和具体实施方式对本发明做进一步详细的说明。

本发明提供一种用于隧道衬砌和铁路路基地质雷达数据的频谱拓展方法，包括对地质雷达数据的频谱分析、小波域自适应拓频处理等。

在一个实施例中，本公开揭示了一种用于提高地质雷达信号分辨率的信号处理方法，所述方法包括以下步骤：

S100、对多道A-scan信号进行一维傅里叶变换，并在A-scan信号的有效频带宽度内确定A-scan信号的基础频段和各拓频频段；

S200、计算A-scan信号各拓频频段的加权系数；

S300、对A-scan信号进行连续小波变换，得到对应的时间-尺度域小波变换系数(时频分布)；

S400、根据频率和尺度的换算关系，对不同时间-尺度范围内的小波变换系数乘以对应的加权系数；

S500、计算乘以对应的加权系数后小波变换系数对应的小波反变换，得到A-scan信号对应的高分辨率信号；

S600、对每个A-scan信号重复步骤S300-S500。

在本实施例中，所述方法基于连续小波变换的多分辨率思想，在连续小波变换域实现对信号的分频拓展处理，以达到提高分辨率的目的。所述方法不需要估计子波信号，对频带能够实现有效划分，而且具有很好的抗噪性能，能够为后续的地质雷达数据处理步骤提供高保真的高分辨率数据。

本实施例的连续小波域信号自适应频率拓展方法可以有效的提高地质雷达数据的时间分辨率以及弱反射的幅度，有利于后期的异常检测。相比较于现有技术1，本实施例提供的方法能够在提高数据分辨率的同时，不增加噪声能量，具有较高的信噪比，且根据倍频程规则实现频带的稳定划分；相比较于现有技术2，本公开提供的方法不需要估计子波信号，仅需要基于傅里叶振幅谱对信号的小波变换结果进行有效的划分即可。

在一个实施例中，所述步骤S100具体包括如下步骤：

S101、对每一道A-scan信号进行一维傅里叶变换，然后计算对应傅里叶振幅谱的平均振幅谱；

S102、基于对应的平均振幅谱选定A-scan信号的基准频率点，作为标准，根据基准频率点确定基础频率段，然后在频率轴上确定各拓展频段；

所述拓展频段包括：高频段的一次、二次和高次谐波段，和低频段的一次次谐波、二次次谐波以及高次次谐波频率段。

在本实施例中，对若干个A-scan进行一维傅里叶变换，将选取的若干A-scan的傅里叶变换的模值对应点相加，然后除以A-scan个数，得到傅里叶振幅谱的平均振幅谱；基于对应的平均振幅谱选定信号的基准频率点，作为标准，根据基准频率点确定基础频率段，然后在频率轴上确定高频段的一次、二次和高次谐波段，和低频段的一次次谐波、二次次谐波以及高次次谐波频率段。

在一个实施例中，步骤S200所述的计算各拓频频段加权系数具体为：

S201、计算信号基础频段内频谱的平均能量密度P为

其中，F(ω)为A-scan的傅里叶变换，ω₀表示基准频率点，ω⁰表示基础频段的下界；

S202、按照步骤S201中的公式，依次计算高频段的一次谐波段、二次谐波段、低频段的一次次谐波段和二次次谐波段的平均能量密度 P₁，P₂，P¹，P²；

S203、确定一次谐波段、二次谐波段、低频段一次次谐波段和二次次谐波段的加权值分别为

在一个实施例中，计算步骤S300所述的连续小波变换的关系式如下：

式中，t表示时间，f(t)表示A-scan信号；τ为时移，τ属于实数集，小波尺度因子a＞0；为小波函数的复共轭，对母小波函数ψ(t)经过伸缩a和时移τ得到

在一个实施例中，步骤S400所述的频率和尺度的换算关系如下：

其中，a为小波尺度因子，ω为信号频率，ω_b为母小波函数的频率。

在一个实施例中，步骤S500所述的连续小波反变换的关系式如下：

式中，为频谱带宽已扩展的A-scan信号时间-尺度分布，

为重构的高分辨率A-scan信号。

在一个实施例中，所述A-scan信号的基准频率点选取为A-scan 信号平均振幅谱的最大值对应的频率值。

在一个实施例中，向A-scan信号的基准频率点的低频段方向幅度下降约30％时对应的频率定义为基础频率段。

在一个实施例中，在频率轴上小于基础频段的频率定义为低频，高于基础频段的频率定义为高频，然后向高频段随着频率升高依次确定一次谐波段(ω₀，ω₁)，二次谐波段(ω₁，ω₂)以及高次谐波段，向低频段随着频率降低依次确定一次次谐波段(ω¹，ω⁰)，二次次谐波段(ω²，ω¹)。

在一个实施例中，所述L个A-scan信号平均振幅谱的计算公式为：

其中：F_l(ω)，l＝1，2，…，L表示L个A-scan信号的傅里叶变换。

在一个实施例中，如图1所示，本发明所提供的地质雷达数据处理方法，具体通过如下步骤实施：

准备探地雷达数据，进行如下处理：

步骤1：对若干A-scan信号进行一维傅里叶变换，并确定各拓频频段，主要包括：

(1)对若干个A-scanf(t)进行一维傅里叶变换，然后计算对应傅里叶振幅谱的平均振幅谱，一维傅里叶变换由以下关系式决定：

选取L个A-scan信号，分别计算其傅里叶变换F_l(ω)，l＝1，2，…，L，平均振幅谱由以下关系式决定：

(2)参照图2，基于对应的平均振幅谱选定信号的基准频率点ω₀，作为标准，ω₀通常选取为平均振幅谱最大值对应的频率值，之后向低频段方向幅度下降约30％时对应的频率定义为ω⁰，确定基础频率段 (ω⁰，ω₀)，在频率轴上小于基础频段的频率定义为低频，高于基础频段的频率定义为高频，然后向高频段随着频率升高依次确定一次谐波段(ω₀，ω₁)，二次谐波段(ω₁，ω₂)以及高次谐波段，向低频段随着频率降低依次确定一次次谐波段(ω¹，ω⁰)，二次次谐波段(ω²，ω¹)等；其中，待扩展的高频段的一次和二次谐波频率段频宽分别为基础频段的一倍和二倍大小，更高次谐波以此类推；待扩展的低频段的一次次谐波和二次次谐波频宽分别为基础频段的一半和四分之一大小；

步骤2：计算各拓频频段的加权系数，其确定方法为：

(1)计算信号基础频段内频谱的平均能量密度P为

其中，F(ω)为A-scan的傅里叶变换；

(2)按照(1)中的公式，依次计算一次谐波段、二次谐波段、低频段一次次谐波段和二次次谐波段的平均能量密度P₁，P₂，P¹，P²；

(3)最后确定待扩展的各频段加权值分别为

其余更高次的谐波和次谐波以此类推；

对B-scan剖面上的每一个A-scan重复以下步骤：

步骤3：对A-scan f(t)进行连续小波变换，得到对应的时间-尺度域分布剖面；

式中，t表示时间，f(t)表示A-scan信号；τ为时移，τ属于实数集，尺度因子a＞0；为小波函数的复共轭，是对母小波函数经过伸缩 a和时移τ得到的，基本频率ω_b＝6.0，为虚数单位。

步骤4：根据频率和尺度的换算关系，对各个尺度范围内的小波变换系数乘以对应的加权系数，具体包括：

(1)根据频率ω和小波尺度a的基本换算关系确定各次谐波频段对应的尺度范围。具体的，一次谐波频段对应的尺度范围为 (a₁，a₀)，则有

(2)计算对应加权后的时间-尺度域分布，对于一次谐波段内的连续小波变换分布，则有：

其余各次谐波和次谐波以此类推；

步骤5：基于时频域扩展后的连续小波域时间-尺度域分布，计算其对应的连续小波反变换，重构相应的高分辨率信号，连续小波反变换由以下公式决定；

式中，为重构的高分辨率A-scan信号。

本发明具有如下有益效果：

1)本发明的方法不需要估计地质雷达信号子波；

2)本发明参数设置简单，根据信号频谱分布自适应确定频谱划分和加权系数；

3)本发明的连续小波域信号自适应拓频方法可以有效的提高地质雷达剖面的时间分辨率以及增强深层弱反射的强度；

4)通过对频带的有效划分和加权处理，本发明的方法对随机噪声有很好的抑制作用。

下面利用本发明提供的数据处理方法对(1)实测的隧道衬砌地质雷达B-scan数据和(2)实测的铁路路基地质雷达B-scan数据进行了处理；处理流程包括步骤1-步骤5，并进行了单个A-scan的波形图及频谱对比。

图3(A)-图3(F)为实测的隧道衬砌地质雷达B-scan剖面及处理结果：

图3(A)为实测的B-scan数据，发射天线的主频为900MHz，隧道管片厚度为35厘米，其中分布两层钢筋，主筋半径为0.01米，配筋半径为0.005米。从图3(A)上可以较为明显地看到浅层的强能量层状杂波以及管片中钢筋反射的双曲线波形。而由于实际介质的频散作用，管片与后方围岩的界面反射在图3(A)上较弱，很难用于对隧道衬砌异常进行有效检测。图3B为高分辨拓频处理结果，可以看出剖面的时间方向分辨率得到了有效提高，而且深层的反射波能量得到了有效增强。图3(C)-图3(F)为抽取的A-scan信号，其中黑色虚线对应未做拓频处理的A-scan信号，黑色实线为拓频后的 A-scan信号。图3(C)为水平距离x＝1.44m的A-scan信号，对比可以看出，经过连续小波域自适应拓频后，多层反射波的振幅得到了有效增强，突出了更多的反射界面。图3(D)为对应图3(C)的频谱分析，可以看到，经过拓频处理后，在有效频带内信号的低频、高频成分得到了增强，随机噪声占据的高频频段并没有得到显著提高，这验证了拓频结果的高可信度。图3(E)为水平距离x＝4.83m处未拓频和已拓频的A-scan，图3(E)为对应的傅里叶振幅谱，可以看出，低频、高频成分得到了提升，因此增加了有效频宽。

图4(A)-图4(F)为实测的铁路路基地质雷达B-scan剖面及处理结果：

图4(A)为实测的用于铁路路基检测的B-scan数据，其中发射天线的主频为900MHz，地下分布有网状的钢筋结构。图中可以看到钢筋网反射的双曲线结构以及层状结构。但是，由于介质衰减、频散等因素的影响，双曲线结构模糊，且与水平界面反射的信号相互干扰。利用本发明专利的方法，对图4(A)B-scan剖面上的每个A-scan 进行拓频处理，提高其时间分辨率。图4(B)为对应的拓频处理剖面，可以看到钢筋反射的双曲线结构得到了有效增强，10ns左右的反射界面得以显现。图4(C)-图4(F)为抽取的A-scan信号，其中黑色虚线对应未做拓频处理的A-scan信号，黑色实线为拓频后的 A-scan信号。对比可以看出，基于小波域的自适应拓频技术有效的拓展了地质雷达信号的带宽，且将噪声限制在一定水平。这为后面的路基病害检测提供了有效的处理结果。

最后需要说明的是，以上模型和实际资料算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。

Claims (10)

1.一种用于提高地质雷达信号分辨率的信号处理方法，其特征在于，所述方法包括以下步骤：

S100、对多道A-scan信号进行一维傅里叶变换，并在A-scan信号的有效频带宽度内确定A-scan信号的基础频段和各拓频频段；

S200、计算A-scan信号各拓频频段的加权系数；

S300、对A-scan信号进行连续小波变换，得到对应的时间-尺度域小波变换系数；

S400、根据频率和尺度的换算关系，对不同时间-尺度范围内的小波变换系数乘以对应的加权系数；

S500、计算所述乘以对应的加权系数后小波变换系数对应的连续小波反变换，得到A-scan信号对应的高分辨率信号；

S600、对每个A-scan信号重复步骤S300-S500。

2.根据权利要求1所述的方法，其特征在于，优选的，所述步骤S100包括如下步骤：

S101、对多道A-scan信号进行一维傅里叶变换，然后计算对应傅里叶振幅谱的平均值，即为A-scan信号的平均振幅谱；

S102、基于对应的平均振幅谱选定A-scan信号的基准频率点，根据基准频率点确定基础频率段，然后在频率轴上确定各拓频频段；

所述拓频频段包括：高频段的一次、二次和高次谐波段，低频段的一次次谐波、二次次谐波以及高次次谐波频率段。

3.根据权利要求2所述的方法，其特征在于，步骤S200所述计算A-scan信号各拓频频段加权系数包括以下步骤：

S201、计算A-scan信号基础频段内频谱的平均能量密度P，

其中，|F(ω)|为A-scan的平均振幅谱，ω₀表示基准频率点，ω⁰表示基础频段的下界；

S202、按照步骤S201中的公式，依次计算A-scan信号高频段的一次谐波段、二次谐波段，低频段的一次次谐波段和二次次谐波段的平均能量密度P₁，P₂，P¹，P²；

S203、确定A-scan信号高频段的一次谐波段、二次谐波段、低频段的一次次谐波段和二次次谐波段的加权值分别为

4.根据权利要求1所述的方法，其特征在于，计算步骤S300所述的连续小波变换的关系式如下：

式中，W(a，τ)表示小波变换系数，t表示时间，f(t)表示A-scan信号；τ为时移，τ属于实数集，小波尺度因子a＞0；为小波函数的复共轭，对母小波函数ψ(t)经过伸缩a和时移τ得到

5.根据权利要求1所述的方法，其特征在于，步骤S400所述的频率和小波尺度因子的换算关系如下：

其中，a为小波尺度因子，ω为信号频率，ω_b为母小波函数的频率。

6.根据权利要求1所述的方法，其特征在于，步骤S500所述的连续小波反变换的关系式如下：

式中，为频谱带宽已扩展的A-scan信号时间-尺度分布， 为重构的高分辨率A-scan信号。

7.根据权利要求2所述的方法，其特征在于：所述A-scan信号的基准频率点ω₀选取为A-scan信号平均振幅谱的最大值对应的频率值。

8.根据权利要求7所述的方法，其特征在于：在A-scan信号的平均振幅谱上，向A-scan信号的基准频率点的低频段方向幅度下降约30％时对应的频率定义为基础频率段的下界ω⁰。

9.根据权利要求8所述的方法，其特征在于：在频率轴上小于基础频段的频率定义为低频，高于基础频段的频率定义为高频，然后向高频段随着频率升高依次确定一次谐波段(ω₀，ω₁)，二次谐波段(ω₁，ω₂)以及高次谐波段，向低频段随着频率降低依次确定一次次谐波段(ω¹，ω⁰)，二次次谐波段(ω²，ω¹)。

10.根据权利要求2所述的方法，其特征在于：所述L个A-scan信号平均振幅谱的计算公式为：

其中：F_l(ω)，l＝1，2，…，L表示L个A-scan信号的傅里叶变换。