CN112578471A - 一种探地雷达杂波噪声去除方法 - Google Patents

Abstract

本发明涉及一种探地雷达杂波噪声去除方法，其包括如下步骤：（a）对原始探地雷达的回波信号进行小波变换，得到不同分解次数下的近似值、水平细节、垂直细节与对角线信号；（b）经步骤（a）小波分解后的各部分信号进行K‑SVD算法处理，选取稀疏的部分进行重构；（c）将重构处理后的各部分信号进行小波逆变换，并进行进一步的自适应双边滤波，得到最终去噪声后的图像。本发明通过引入对随机噪声的压制，更加有利于处理雷达信号中的随机噪声干扰，提高雷达图像的信噪比，同时引入自适应双边滤波，有效地保留图像的边缘和细节特征，提高了探地雷达回波信号质量和图像的清晰率和准确率。

Description

一种探地雷达杂波噪声去除方法

技术领域

本发明涉及一种探地雷达杂波噪声去除方法。

背景技术

探地雷达(ground penetrating radar，GPR)是近些年发展最为迅速的地球物理工具之一，它通过向地下待测物体发射高频电磁波，对接收的回波进行相应的数据采样以及相关分析处理来确定目标物分布和性质的一种无损检测的方法。但是在雷达检测的过程中受到地面反射、地下复杂多变的环境和雷达发射和接收天线之间的耦合的影响，会使检测到的雷达数据含有噪声和杂波的干扰，影响对于目标信号的检测和分析，所以对于杂波噪声干扰的抑制对于改善探地雷达回波信号质量和提高图像的准确率十分重要。

主成分分析法是一种建立在最小均方误差基础上的线性变换处理法。探地雷达接收的B-scan信号可以表示为矩阵W_R∈R^M×N，其中M为时间上的采样点数，N为某测线的测量点数目。对W_R进行奇异值分解。

其中U＝{u₁u₂…u_N}∈R^M×N，V＝{v₁v₂…v_N}∈R^M×N为正交矩阵，其列向量分别为对称矩阵

和

的特征向量。D∈R^N×N为对角矩阵，对角成分元素为奇异值，按照从大到小排列，可以表示为D_1,1≥D_2,2≥…≥D_N,N。一般选择y_i＝D_i,iv_i(i＝1,2,…,N)作为“主成分”，则W_R可以认为是主成分y_i与相应特征信号u_i的加权和，即

奇异值D_i,i的大小反映了子信号y_i与信号W_R的相关性。接收信号矩阵可以表示为W_R＝S_R+D_R，其中S_R与D_R分别为目标回波信号矩阵和杂波信号矩阵。杂波对应于强相关性的主成分；而接收信号中随机变化的噪声则对应弱相关的主成分。去除掉能量较大和较小的分量，选取合适的主成分进行数据重构，达到抑制杂波和随机噪声的目的。

随着复信号分析技术的发展，探地雷达信号噪声抑制的方法有空间滤波技术和变换域滤波技术。其中空间滤波技术主要包括均值滤波、中值滤波、Lee滤波等，这些方法相对简单易于实现，但是会造成图像边缘和线性目标的模糊。变换域滤波技术主要包括小波变换、平稳小波、Bandelet变换、Curvelet变换和非下采样Contourlet变换等。变换域滤波相比于经典空间滤波方法，对于图像的边缘及线性目标的保持能力有了很大的提高。但是由于噪声和图像边缘具有相似的频率特性，都是高频信号，因此也会在噪声抑制后的图像在均匀区域和边缘附近常会产生伪吉布斯效应。同时以上技术对于随机噪声的处理效果并不明显。

发明内容

由于地面反射，探地雷达发射和接收天线之间的耦合以及地下复杂多变的情况会严重破坏地下目标检测信号，本发明的目的是提供一种用于去除探地雷达信号中的杂波以及噪声干扰，提取探地雷达有效信号的探地雷达杂波噪声去除方法。

本发明采用如下技术方案：

一种探地雷达杂波噪声去除方法，其包括如下步骤：

(a)对原始探地雷达的回波信号进行小波变换，得到不同分解次数下的近似值、水平细节、垂直细节与对角线信号；

(b)经步骤(a)小波分解后的各部分信号进行K-SVD算法处理，选取稀疏的部分进行重构；

(c)将重构处理后的各部分信号进行小波逆变换，并进行进一步的自适应双边滤波，得到最终去噪声后的图像。

其中，步骤(a)中，采用三次小波分解。

其中，步骤(a)中，所述小波变换为二维小波变换。

其中，所述步骤(b)具体包括如下步骤：

(1)将经步骤(a)小波分解后的各部分信号重叠分块，并矢量化表示；

(2)用DCT字典作为初始字典，在此基础上更新；

(3)用初始字典将步骤(1)的分块进行OMP稀疏编码，由正交匹配追踪算法求解稀疏系数；

(4)从字典中挑选与之最相关的原子，更新残差，迭代进行，直至更新所有原子；

(5)使用最后一次迭代得到的更新稀疏系数和更新字典重建纯净图像X分块；

(6)将各分块平均处理得到重构后的各个部分信号。

其中，更新字典D和更新稀疏系数a_ij通过求解如下最优化问题得到：

其中，参数λ为正则化参数，实值μ_ij表示第i行第j块对应的参数；

代表无噪图像中第i行的第j个小图像块在字典D下的稀疏系数。

其中，采用交替方向乘子法，通过变量的交替更新和迭代，求解更新字典D和更新稀疏系数a_ij以及纯净图像X。

其中，步骤(c)中，所述小波逆变换为二维小波逆变换。

其中，步骤(c)中，使用双边滤波器进行自适应双边滤波。

其中，所述自适应双边滤波通过如下公式进行：

其中，

为输出图像，f为输入图像，Ω是以像素点x为中心的领域窗口，ω为滤波核；

双边滤波的滤波核由两部分乘积组成：空间核φ与值域核ψ。

其中，滤波核采用高斯函数形式所示：

ω(x，y)＝φ(x，y)ψ(x，y)

其中，σ_d为空域高斯函数的标准差，σ_r为值域高斯函数的标准差。

本发明的有益效果在于：本发明采用基于改进后的小波变换下的K-SVD算法进行去除雷达信号噪声杂波信号，引入随机信号噪声的压制，主要目的是更有利于处理雷达信号中的随机噪声干扰，提高雷达信号的信噪比，适应较为复杂的随机噪声性质，尽可能的保持信号的特征，更加有效地改善探地雷达回波信号质量和提高图像准确率。同时采用双边滤波，平滑滤波的同时能大量保留图像的边缘和细节特征。

附图说明

图1本发明方法的流程示意图。

图2本发明中K-SVD算法的流程示意图。

图3原始含噪图像。

图4经均值滤波法处理后的图像。

图5经中值滤波法处理后的图像。

图6经本发明的方法处理后的图像。

具体实施方式

下面将结合本申请附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本申请及其应用或使用的任何限制。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

如图1所示，本发明是一种基于改进后小波变换的K-SVD算法去除探地雷达杂波噪声的方法。利用小波变换的方法在时频域内滤除随机噪声，通过尺度伸缩可以改变时间和频率的分辨率，实现在不同尺度下的时频域处理。经二维小波变换得到各个部分进行K-SVD算法去噪，选取稀疏的部分进行重构，对重构各个部分信号进行二维小波反变换，实现二维小波域的K-SVD算法去噪。

具体的，其包括如下步骤。

一、对原始探地雷达的回波信号进行小波变换，得到不同分解次数下的近似值、水平细节、垂直细节与对角线信号。

小波变换(wavelet transform，WT)是一种新的变换分析方法，不仅具有短时傅立叶变换局部化的思想，同时又克服了窗口大小不随频率变化等缺点。小波变换的方法是在时频域内滤除随机噪声，小波通过尺度伸缩可改变时间和频率的分辨率，实现在不同尺度下的时频域处理。

因为探地雷达回波信号是非平稳信号，采用小波变换来处理探地雷达回波数据是比较合适的。二维小波变换将原始图像分解为近似解、水平细节、垂直细节与对角线细节部分，对近似解进行多次分解，对图像剖分更加细化。

根据雷达数据的实际情况考虑，本发明算法采用三次小波分解。

对二维小波变换得到各个部分进行K-SVD算法去噪，选取稀疏的部分进行重构，对重构各个部分信号进行二维小波反变换，实现二维小波域的K-SVD算法去噪。

二、经小波分解后的各部分信号进行K-SVD算法处理，选取稀疏的部分进行重构。

K-SVD(K-Singular Value Decomposition，K-SVD)，K次迭代的奇异值分解算法，是一种字典训练算法，依据误差最小原则，对误差项进行奇异值分解，选择使误差最小的分解项作为更新的字典原子和对应的原子系数，通过迭代优化输入数据在当前字典的表示和更新字典中的元素从而得到优化的解。

K-SVD算法用于解决图像去噪问题时，含噪雷达图像Y表示为：

Y＝X+B 式(1)

式(1)中，纯净图像X被标准差为σ均值为0的高斯噪声B污染。X、B和Y代表的都是大小

的方形图像。向量化后都是长度为N的列向量。K-SVD去噪算法是假设X是可稀疏的，而随机噪声B是不可稀疏的，利用学习到的字典对信号稀疏表示进行去噪。

图像越大对应的字典就越大，为了降低计算的复杂程度，字典学习通常在大小为

的数据块上进行。在去噪的过程中首先将图像按照一定步长重叠分为若干小块。无噪图像的第i行第j个块可以表示为x_ij＝R_ijX，x_ij∈Rⁿ，R_ij代表分块算子。在训练好的字典D∈R^n×p下，可以得到长度为p的稀疏系数。将图像块排列为列向量，每个干净的图像块可以由稀疏系数a_ij和字典D表示为：

x_ij＝Dα_ij 式(2)

则去噪问题转化为求解字典D和稀疏系数的问题。K-SVD算法中字典D和稀疏系数通过求解如下最优化问题得到：

式(3)中，参数λ为正则化参数，实值μ_ij表示第i行第j块对应的参数。0范数表明稀疏系数向量中非零值的个数。0范数值越小，稀疏性越好。

代表无噪图像中第i行的第j个小图像块在字典D下的稀疏系数。

式(3)中，第一项为保真项

保证去噪后的图像要和含噪图像相似，不能丢失过多细节。第二项∑_ij||α_ij||₀要保证所得稀疏系数尽可能稀疏。最后一项

是训练好的字典和稀疏系数组合，和对应的图像块之间差越小越好，因为差越小对应去噪后残留噪声越少。

K-SVD去噪就是要从式(3)中求解字典、稀疏系数和纯净图像。对于一个式子中同时含有三个变量的最优化问题，很难直接求解，所以采用交替方向乘子法(ADMM)，通过变量的交替更新和迭代来求解这三个子问题。

通过变量的交替更新和迭代，将式(3)求解将其分为三个子问题：

第一个问题是稀疏编码：给定字典和X，求解稀疏系数相当于求解如下代价函数的最小值：

对每个小块都进行这步操作，相当于一个滑动窗走遍了含噪图像的所有像素，这个稀疏编码问题可以由正交匹配追踪算法求解。

第二个问题是更新字典：固定稀疏系数和X更新字典：

第三个问题是计算无噪图像：交替进行以上两步计算直到分块算子遍布所有小块，并且在每个小块的计算过程中字典中原子全部参与了更新之后，便可以求出纯净图像X：

K-SVD具有良好的图像去噪性能。该降噪算法学习出让含噪图像尽量稀疏的字典，使得在该字典下图像是可稀疏的，随机噪声是不可稀疏的。

如图2所示，基于K-SVD的噪声压制方法如下：

首先，将含噪雷达图像重叠分块，并矢量化表示。

其次，用DCT字典作为初始字典，在此基础上更新。

去噪时，用该字典将含噪图像小块进行OMP稀释编码，由正交匹配追踪算法求解稀疏系数。每次从字典中挑选与之最相关的原子，更新残差，迭代进行。

更新原子：根据待处理信号和稀疏系数通过SVD分解更新字典中原子。字典逐列更新，每次更新其中一个原子，找到所有与之对应的稀疏系数，进而得到误差矩阵。将误差矩阵进行SVD分解。最后更新字典中原子的值为左奇异矩阵的第一个列向量。更新稀疏系数中与该原子对应的值为右奇异矩阵的第一个列向量乘以第一个奇异值。

稀疏编码和更新原子的过程为块内去噪阶段，噪声在这个阶段被滤除。有数据被滤除是因为SVD分解时只保留了最大奇异值对应的部分。块内去噪在重建图像之前可以根据雷达实际图像混乱程度重复执行多次，在一定范围内(字典D收敛为止)，即在字典收敛的范围内进行不断迭代更新，直至收敛为止，重复执行次数越多去噪效果越好，稀疏系数越稀疏。重复次数过多时由于稀疏系数中非零值过于稀少，导致重建的图像细节丢失，信号中细节和噪声一起被滤除。

字典中原子全部参与了更新以后，用最后一次迭代得到的稀疏系数和字典重建小块。

最后，将各小块平均处理得到去噪后结果。

三、将重构处理后的各部分信号进行小波逆变换，并进行进一步的自适应双边滤波，得到最终去噪声后的图像。

双边滤波器(Bilateral Filter)是一种非线性滤波器，该算法结合空间信息和亮度相似性对图像进行滤波处理，在平滑滤波的同时能大量保留图像的边缘和细节特征。和传统的影响平滑化算法不同，双边滤波器除了使用像素之间几何上的靠近程度之外，还考虑了像素之间广度和色彩的差异，使得双边滤波器能够有效的将影像上的噪声去除，更好地保存了影像的边缘信息。定义如式(7)所示。

式(7)中，

为输出图像，f为输入图像，Ω是以像素点x为中心的领域窗口，ω为滤波核。双边滤波的滤波核由两部分乘积组成：空间核φ与值域核ψ。两个滤波核通常采用高斯函数形式如式(8)、式(9)、式(10)所示。

ω(x，y)＝φ(x，y)ψ(x，y) 式(8)

其中，σ_d为空域高斯函数的标准差，σ_r为值域高斯函数的标准差。空域滤波系数是由像素间的空间距离决定，距离越小，系数越大。值域滤波系数由像素间的相似度决定，像素值越接近系数越大。

在灰度变化平缓区域，值域滤波系数接近1，此时空域滤波起到主要作用，双边滤波器退化为传统的高斯低通滤波器，对图像进行平滑操作。在图像变化剧烈的部分(图像边缘)，像素间差距较大，值域滤波起主要作用，因而能有效保持边缘信息。

四、仿真验证。

通过均值滤波、中值滤波和本发明的Wavelet-KSVD-BF算法(块内去噪重复28次)对含有噪声的探地雷达图像信号(图3)进行仿真对比，具体结果如图4～图6所示。

通过分别将探地雷达原始图像数据在均值滤波法、中值滤波法以及本发明中的Wavelet-KSVD-BF算法下进行仿真对比，可以明显发现均值滤波和中值滤波还是可以去除地表直达波和部分噪声，但是效果却不太好，图像数据还存在着噪声的干扰，图像和目标信号较为模糊；相比之下Wavelet-KSVD-BF算法基本去除掉了雷达数据中的杂波噪声的干扰，并且基本保持了目标信号的特征，达到了较好的去噪滤波效果。证明了本算法的具有良好的去噪滤波效果，可以有效提取探地雷达的有效信号。

本发明算法通过引入对随机噪声的压制，更加有利于处理雷达信号中的随机噪声干扰，提高雷达图像的信噪比，同时引入自适应双边滤波，有效地保留图像的边缘和细节特征，提高了探地雷达回波信号质量和图像的清晰率和准确率。

以上实施方式仅为本发明的优选实施例，而并非本发明可行实施的穷举。对于本领域一般技术人员而言，在不背离本发明原理和精神的前提下对其所做出的任何显而易见的改动，都应当被认为包含在本发明的权利要求保护范围之内。

Claims (10)

1.一种探地雷达杂波噪声去除方法，其特征在于，其包括如下步骤：

(a)对原始探地雷达的回波信号进行小波分解，得到不同分解次数下的近似值、水平细节、垂直细节与对角线信号；

(b)经步骤(a)小波分解后的各部分信号进行K-SVD算法处理，选取稀疏的部分进行重构；

(c)将重构处理后的各部分信号进行小波逆变换，并进行进一步的自适应双边滤波，得到最终去噪声后的图像。

2.根据权利要求1所述的一种探地雷达杂波噪声去除方法，其特征在于，步骤(a)中，采用三次小波分解。

3.根据权利要求2所述的一种探地雷达杂波噪声去除方法，其特征在于，步骤(a)中，所述小波变换为二维小波变换。

4.根据权利要求3所述的一种探地雷达杂波噪声去除方法，其特征在于，所述步骤(b)具体包括如下步骤：

(1)将经步骤(a)小波分解后的各部分信号重叠分块，并矢量化表示；

(2)用DCT字典作为初始字典，在此基础上更新；

(3)用初始字典将步骤(1)的分块进行OMP稀疏编码，由正交匹配追踪算法求解稀疏系数；

(4)从字典中挑选与之最相关的原子，更新残差，迭代进行，直至更新所有原子；

(5)使用最后一次迭代得到的更新稀疏系数和更新字典重建纯净图像X分块；

(6)将各分块平均处理得到重构后的各个部分信号。

5.根据权利要求4所述的一种探地雷达杂波噪声去除方法，其特征在于，更新字典D和更新稀疏系数a_ij通过求解如下最优化问题得到：

其中，参数λ为正则化参数，实值μ_ij表示第i行第j块对应的参数；

代表无噪图像中第i行的第j个小图像块在字典D下的稀疏系数。

6.根据权利要求5所述的一种探地雷达杂波噪声去除方法，其特征在于，采用交替方向乘子法，通过变量的交替更新和迭代，求解更新字典D和更新稀疏系数a_ij以及纯净图像X。

7.根据权利要求1所述的一种探地雷达杂波噪声去除方法，其特征在于，步骤(c)中，所述小波逆变换为二维小波逆变换。

8.根据权利要求1所述的一种探地雷达杂波噪声去除方法，其特征在于，步骤(c)中，使用双边滤波器进行自适应双边滤波。

9.根据权利要求8所述的一种探地雷达杂波噪声去除方法，其特征在于，所述自适应双边滤波通过如下公式进行：

其中，

为输出图像，f为输入图像，Ω是以像素点x为中心的领域窗口，ω为滤波核；

双边滤波的滤波核由两部分乘积组成：空间核φ与值域核ψ。

10.根据权利要求8所述的一种探地雷达杂波噪声去除方法，其特征在于，滤波核采用高斯函数形式所示：

ω(x,y)＝φ(x,y)ψ(x,y)

其中，σ_d为空域高斯函数的标准差，σ_r为值域高斯函数的标准差。

Images