CN106199532A - 基于混合傅立叶‑小波分析的探地雷达信号降噪方法 - Google Patents

Abstract

本发明提供了一种基于混合傅立叶‑小波分析的探地雷达信号降噪方法，通过对探地雷达回波信号进行多尺度小波分解，并利用阈值函数对小波系数进行更新，继承了傅立叶变换频域高分辨率的优点，并克服了傅立叶变换不适应于非平稳信号，不能提取局部时间段的频率特征信息的不足，充分利用了小波变换局部化与自动变焦的优点，弥补了小波变换对于低频噪声去除效果不佳的缺陷；充分考虑探地雷达的具体特性，给出了使用过程中每一步骤的具体公式，对于实际应用中不同的探地雷达，可以直接代入其采样率，中心频率、带宽等参数进行去噪处理，操作简便。

Description

基于混合傅立叶-小波分析的探地雷达信号降噪方法

技术领域

本发明涉及雷达技术领域，更具体地涉及一种基于混合傅立叶-小波分析的探地雷达信号降噪方法。

背景技术

探地雷达作为一种新型的地球物理探测仪器，具有探测效率高、操作简单、采样迅速、无损伤探测、精度较高、探测分辨率高、抗干扰能力较强等优点，已成为浅层地下探测的有效手段。在探地雷达探测中，为了得到更多的反射波特征，通常利用宽频带进行记录，因此在记录各种有效波的同时，不可避免地记录了各种干扰噪声。通常，没有经过处理的剖面无法进行地质解释。因而，必须对接收的探地雷达回波信号进行去噪处理，以改善数据质量，为进一步的地质解释提供清晰可靠的探地雷达信号与剖面图。

傅立叶分析是传统的信号去噪与分析方法，它将信号变换到频率域进行处理，对于平稳白噪声去噪效果良好。但傅立叶变换是一种全局变换，在实际应用中，由于探地雷达回波信号是非平稳信号，含有许多突变与尖峰，此时傅立叶分析不能给出信号在某个时间点上的变化情况，信号在时间域的任何一个突变，都会影响信号的整个频谱，即无法表述探地雷达信号的时频局域特性。

针对傅立叶变换不能提取局部时间段的频率特征信息，Gabor提出了短时傅立叶变换，该方法基本思想是用一个随时间平移的窗函数将原来的非平稳信号分成若干平稳或近似平稳段，然后逐段确定其频谱。然而该方法有两个缺陷：一是窗口不能随频率变化，二是缺乏离散正交基。

小波分析是一种有效的时频分析方法，由Morlet提出，不仅继承和发展了短时傅立叶变换的局部化思想，而且克服了它的缺点。该方法能同时在时域中对信号进行分析，并且在频域内分辨率高时，时域内分辨率响应较低；在频域内分辨率低时，时域内分辨率响应较高；即具有自动变焦功能，因此能有效地区分信号中的突变部分和噪声，从而实现信号的去噪。但在探地雷达实际信号处理中，小波分析只对于信号中较高频分量去除效果较好，无法有效去除信号中的低频噪声。

在实现本发明的过程中，申请人发现上述现有技术存在如下技术缺陷：常用的傅里叶变换无法描述信号时频局部特性，而由于地下介质造成的衰减，探地雷达信号是典型的非平稳信号，因此傅立叶变换不能提取局部时间段的频率特征信息，去噪效果不佳；传统的离散小波变换阈值去噪方法虽然具备较好的时频特性，但大多数情况下只适用于去除信号中的高频分量，而对于信号中固有的低频噪声抑制效果较差；上述两种方法在操作时需要根据经验选取合适参数以达到良好效果，在实际工程应用中不够简便。

发明内容

(一)要解决的技术问题

为了解决现有技术问题，本发明提供了一种基于混合傅立叶-小波分析的探地雷达信号降噪方法。

(二)技术方案

本发明提供了一种基于混合傅立叶-小波分析的探地雷达信号降噪方法，包括：步骤A：输入探地雷达回波信号，并对所述探地雷达回波信号进行保守频域滤波；步骤B：选取小波基并确定小波分解级数；步骤C：对所述探地雷达回波信号进行小波多尺度分解，得到各级分解尺度的小波系数；步骤D：计算各级分解尺度的阈值并选择阈值函数，得到更新后的各级分解尺度的小波系数；步骤E：基于所述更新后的各级分解尺度的小波系数，对探地雷达回波信号进行重构；以及步骤F：对重构后的探地雷达回波信号进行频域修正，输出去噪后的探地雷达回波信号。

(三)有益效果

从上述技术方案可以看出，本发明的基于混合傅立叶-小波分析的探地雷达信号降噪方法具有以下有益效果：

(1)本发明对探地雷达回波信号进行多尺度小波分解，并利用阈值函数对小波系数进行更新，继承了傅立叶变换频域高分辨率的优点，并克服了傅立叶变换不适应于非平稳信号，不能提取局部时间段的频率特征信息的不足，充分利用了小波变换局部化与自动变焦的优点，弥补了小波变换对于低频噪声去除效果不佳的缺陷；

(2)本发明的方法充分考虑探地雷达的具体特性，给出了使用过程中每一步骤的具体公式，对于实际应用中不同的探地雷达，可以直接代入其采样率，中心频率、带宽等参数进行去噪处理，操作简便。

附图说明

图1为本发明实施例的原始探地雷达单道信号；

图2为本发明实施例的各分解尺度的小波系数；

图3为本发明实施例的软阈值函数；

图4为本发明实施例的去噪后探地雷达单道信号；

图5为本发明实施例的基于混合傅立叶-小波分析的探地雷达信号降噪方法的流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

本发明实施例的基于混合傅立叶-小波分析的探地雷达信号降噪方法，参照图5，该探地雷达信号降噪方法包括：

步骤A：输入探地雷达回波信号，并对探地雷达回波信号进行保守频域滤波。

步骤A包括：输入探地雷达回波信号，将探地雷达回波信号通过零相位带通滤波器，零相位带通滤波器的上截止频率f_H和下截止频率f_L的计算公式为：

$\begin{array}{c}{f}_H=f+0.75B\\ f_L=f-0.75B\end{array}---\left(1\right)$

其中，f为探地雷达的中心频率，B为探地雷达的带宽。

步骤B：选取小波基并确定小波分解级数。

步骤B包括：选取小波基，确定小波分解级数，小波分解级数n的计算公式为：

其中，表示向上取整，f_s为探地雷达的采样率，f为探地雷达的中心频率，B为探地雷达的带宽。

在本发明中，可以选取Haar小波基、Morlet小波基、Meyer小波基、db5小波基等作为小波基函数，优选db5小波基。

步骤C：对探地雷达回波信号进行小波多尺度分解，得到各级分解尺度的小波系数。

优选地，利用Mallat算法对探地雷达回波信号进行小波多尺度分解；其中，所述各级分解尺度的小波系数包括n级分解尺度的高频系数W₁ ^H、W₂ ^H、...、W_n ^H和最高分解尺度的低频系数W^L，其中，n为步骤B中确定的小波分解级数。

步骤D：计算各级分解尺度的阈值并选择阈值函数，得到更新后的各级分解尺度的小波系数。

步骤D包括：阈值T的计算公式为：

$\left\{\begin{array}{c}T={\mathrm{\σ}}_n\sqrt{2log\left(N\right)}\\ {\mathrm{\σ}}_n=median\left({W_1}^H\right)/0.6475\end{array}\right.---\left(3\right)$

其中，N为各级分解尺度的小波系数的个数，σ_n为噪声标准差估计值，median(·)表示求取中位数，W₁ ^H为第一级分级尺度的高频系数，利用公式(3)，分别将n级分解尺度的高频系数W₁ ^H、W₂ ^H、...、W_n ^H和最高分解尺度的低频系数W^L的小波系数个数N代入公式(3)，得到各级分解尺度的阈值。

阈值函数T_soft(y)的公式为：

$T_{soft}\left(y\right)=\left\{\begin{array}{c}0,\left|y\right|<T\\ \mathrm{sgn}\left(y\right)\left(\left|y\right|-T\right),\left|y\right|\&GreaterEqual;T\end{array}\right.---\left(4\right)$

其中，T为阈值，y为各级分解尺度的小波系数值，T_soft(y)为更新后的各级分解尺度的小波系数值，各级分解尺度的阈值代入公式(3)，分别由n级分解尺度的高频系数W₁ ^H、W₂ ^H、...、W_n ^H的小波系数值得到更新后的n级分解尺度的高频系数，由最高分解尺度的低频系数W^L的小波系数值得到更新后的最高分解尺度的低频系数。

步骤E：基于更新后的各级分解尺度的小波系数，对探地雷达回波信号进行重构。

优选地，利用Mallat算法对探地雷达回波信号进行重构。

由此可见，本发明继承了傅立叶变换频域高分辨率的优点，并克服了傅立叶变换不适应于非平稳信号，不能提取局部时间段的频率特征信息的不足，充分利用了小波变换局部化与自动变焦的优点，弥补了小波变换对于低频噪声去除效果不佳的缺陷。

步骤F：对重构后的探地雷达回波信号进行频域修正，输出去噪后的探地雷达回波信号。

步骤F包括：使重构后的探地雷达回波信号通过零相位带通滤波器，得到去噪后的探地雷达回波信号并输出，零相位带通滤波器的上截止频率f_H和下截止频率f_L的计算公式为：

$\begin{array}{c}{f}_H=f+0.75B\\ f_L=f-0.75B\end{array}---\left(5\right)$

其中，f为探地雷达的中心频率，B为探地雷达的带宽。

由此可见，本发明的方法充分考虑探地雷达的具体特性，给出了使用过程中每一步骤的具体公式，对于实际应用中不同的探地雷达，可以直接代入其采样率，中心频率、带宽等参数进行去噪处理，操作简便。

以下通过具体示例对本发明的方法进行说明：

(1)输入探地雷达回波信号，如图1所示，其采样率为1GHz，中心频率为25MHz，带宽为25MHz。对其进行零相位带通滤波，根据公式(1)计算得到带通滤波器的上、下截止频率为44MHz、6MHz。

(2)选取db5小波基作为小波变换的小波基，根据公式(2)确定分解级数为6级。

(3)使用Mallat算法对探地雷达回波信号进行小波多尺度分解，得到各分解尺度的高频系数W₁ ^H、W₂ ^H、...、W₆ ^H与最高分解尺度的低频系数W^L，参见图2。

(4)根据公式(3)计算各个分解尺度的阈值，并选取软阈值函数作为阈值函数，软阈值函数图参见图3，根据公式(4)计算更新后的小波系数。

(5)使用Mallat算法对更新后的小波系数进行重构。

(6)对重构后的信号进行零相位带通滤波，带通滤波器的上、下截止频率与第一步相同，输出去噪后的探地雷达信号，本例得到的去噪后探地雷达单道信号参见图4，可以看到，探地雷达信号噪声明显得到了有效去除。

至此，已经结合附图对本发明实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明的一种基于层次化记忆网络的问答方法有了清楚的认识。

需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件的定义并不仅限于实施例中提到的各种方式，本领域普通技术人员可对其进行简单地更改或替换，例如：

(1)实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来限制本发明的保护范围；

(2)上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于混合傅立叶-小波分析的探地雷达信号降噪方法，其特征在于，包括：

步骤A：输入探地雷达回波信号，并对所述探地雷达回波信号进行保守频域滤波；

步骤B：选取小波基并确定小波分解级数；

步骤C：对所述探地雷达回波信号进行小波多尺度分解，得到各级分解尺度的小波系数；

步骤D：计算各级分解尺度的阈值并选择阈值函数，得到更新后的各级分解尺度的小波系数；

步骤E：基于所述更新后的各级分解尺度的小波系数，对探地雷达回波信号进行重构；以及

步骤F：对重构后的探地雷达回波信号进行频域修正，输出去噪后的探地雷达回波信号。

2.如权利要求1所述的探地雷达信号降噪方法，其特征在于，所述步骤A包括：

输入探地雷达回波信号，将所述探地雷达回波信号通过零相位带通滤波器，所述零相位带通滤波器的上截止频率f_H和下截止频率f_L的计算公式为：

$\begin{array}{c}{f}_H=f+0.75B\\ f_L=f-0.75B\end{array}---\left(1\right)$

其中，f为探地雷达的中心频率，B为探地雷达的带宽。

3.如权利要求1所述的探地雷达信号降噪方法，其特征在于，所述步骤B的小波分解级数的计算公式为：

其中，n为小波分解级数，表示向上取整，f_s为探地雷达的采样率，f为探地雷达的中心频率，B为探地雷达的带宽。

4.如权利要求1所述的探地雷达信号降噪方法，其特征在于，所述步骤B的小波基为db5小波基。

5.如权利要求1所述的探地雷达信号降噪方法，其特征在于，在所述步骤C中，利用Mallat算法对探地雷达回波信号进行小波多尺度分解。

6.如权利要求3所述的探地雷达信号降噪方法，其特征在于，在所述步骤C中，所述各级分解尺度的小波系数包括n级分解尺度的高频系数W₁ ^H、W₂ ^H、...、W_n ^H和最高分解尺度的低频系数W^L。

7.如权利要求6所述的探地雷达信号降噪方法，其特征在于，所述步骤D的各级分解尺度的阈值的计算公式为：

$\left\{\begin{array}{c}T={\mathrm{\&sigma;}}_n\sqrt{2log\left(N\right)}\\ {\mathrm{\&sigma;}}_n=median\left({W_1}^H\right)/0.6475\end{array}\right.---\left(3\right)$

其中，T为阈值，N为各级分解尺度的小波系数的个数，σ_n为噪声标准差估计值，median(·)表示求取中位数，W₁ ^H为第一级分级尺度的高频系数。

8.如权利要求7所述的探地雷达信号降噪方法，其特征在于，所述步骤D的阈值函数的计算公式为：

$T_{soft}\left(y\right)=\left\{\begin{array}{c}0,\left|y\right|<T\\ \mathrm{sgn}\left(y\right)\left(\left|y\right|-T\right),\left|y\right|\&GreaterEqual;T\end{array}\right.---\left(4\right)$

其中，T为阈值，y为各级分解尺度的小波系数值，T_soft(y)为更新后的各级分解尺度的小波系数值。

9.如权利要求5所述的探地雷达信号降噪方法，其特征在于，在所述步骤E中，利用Mallat算法对探地雷达回波信号进行重构。

10.如权利要求1所述的探地雷达信号降噪方法，其特征在于，所述步骤F包括：使重构后的探地雷达回波信号通过零相位带通滤波器，得到去噪后的探地雷达回波信号并输出，所述零相位带通滤波器的上截止频率f_H和下截止频率f_L的计算公式为：

$\begin{array}{c}{f}_H=f+0.75B\\ f_L=f-0.75B\end{array}---\left(5\right)$

其中，f为探地雷达的中心频率，B为探地雷达的带宽。