CN108983158B - 一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法 - Google Patents

Abstract

本发明公开了一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法。首先，利用探地雷达B扫描图像的每一道数据，构造Hankel矩阵；其次，对所述Hankel矩阵进行奇异值分解；然后，求取奇异值差分谱，并以差分谱的均值作为阈值判断目标信号奇异值和噪声信号奇异值的分界点；最后，利用目标信号奇异值进行重构，得到去噪后的数据。本发明的有益效果是：本发明根据目标信号和噪声信号的Hankel矩阵奇异值分布差异，利用奇异值差分谱的均值作为阈值自动确定目标信号奇异值，计算简便，阈值稳定性好，可有效抑制探地信号中的噪声。

Description

一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法

技术领域

本发明涉及数字信号处理领域，具体涉及到探地雷达的噪声抑制问题，尤其涉及一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法。

背景技术

探地雷达基于电磁波传播与散射原理，通过向地下发射电磁波信号并接收地下介质不连续处散射的回波实现对地下目标的探测。与电阻率法、低频电磁感应法和地震法等探测方法相比，探地雷达具有探测速度快、探测过程连续、分辨率高、操作方便灵活、费用低、探测范围广(能探测金属和非金属)等优点，在地质、资源、环境、工程和军事等领域得到广泛的应用。

在对地下目标进行探测时，受地下杂波和周围介质的干扰等，探地雷达回波信号通常是一种弱信噪比的非平稳信号，因此如从噪声背景中提取目标信息，是探地雷达信号处理的一个重要研究领域。

目前常用的探地雷达噪声抑制方法有傅立叶变换、小波变换、S变换和主成分分析法等。傅立叶变换只能反映信号的整体特征,不适用于频率随时间变化的非平稳信号；小波变换中小波基函数、分解层数以及阈值的选择都依赖于主观经验，缺乏自适应性；S变换中时频滤波器的设计较为复杂，也限制了其应用场合。主成分分析法(Princinal ComponetAnalysis，PCA)是一种建立在最小均方误差基础上的线性变换处理方法，其算法核心是用特定的正交矩阵对信号矩阵进行正交变换，得到相互正交的对角主成分矩阵。奇异值分解算法是PCA方法中的常用算法，在该方法中，将奇异值对应的信号分量称为主成分，主成分的选择通常是以主成分对应的奇异值选择为基础的。PCA方法去噪的原理就是选择目标信号的奇异值重构，以去除噪声信号。

目前，在传统的探地雷达PCA去噪方法中，一般是对探地雷达的B扫描图像进行奇异值分解，然后选择目标信号奇异值重构。选择奇异值的方法主要有定性的经验方法和特征值能量百分比方法。经验方法根据经验确定目标信号的奇异值，容易受到人的主观性影响；特征值能量百分比方法可确定目标信号奇异值，但该方法比较复杂，阈值稳定性差，受具体探测条件影响大。因此，如何以较小的复杂度准确选择目标信号奇异值，对于提高探地雷达的噪声抑制性能具有重要意义。

发明内容

为了解决上述问题，本发明提供了一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法，主要包括以下步骤：

S101：输入探地雷达B扫描图像X∈R^M×N，其中M为道数，N为每道数据的采样点数；

S102：根据所述探地雷达B扫描图像的一道数据，构造一个Hankel矩阵；

S103：对所述Hankel矩阵进行奇异值分解；

S104：求取所述奇异值差分谱；

S105：通过计算得到所述奇异值差分谱的均值，所述均值作为阈值；

S106：根据所述阈值，确定目标信号奇异值与噪声信号奇异值的分界点，并利用所述目标信号奇异值进行重构，得到去噪后的数据；

S107：根据步骤S102～S106，对所述探地雷达B扫描图像中的每一道数据进行处理，得到去噪后的探地雷达B扫描图像X'∈R^M×N。

进一步地，在步骤S104中，利用公式(1)求取所述奇异值差分谱：

q_i＝σ_i-σ_i+1,i＝1,2,…r-1 (1)

其中，q_i为所述奇异值差分谱，σ_i和σ_i+1为所述Hankel矩阵B的奇异值，且B∈R^m×n，r＝min(m,n)，m＝N-n+1，1＜n＜N。

进一步地，在步骤S105中，利用公式(2)求取所述奇异值差分谱的均值T，以所述均值T作为奇异值判断的阈值：

其中，T为所述奇异值差分谱的均值，即为所述阈值，q_i为所述奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N。

进一步地，在步骤S106中，利用相邻的三个奇异值差分谱和所述阈值进行比较，得到目标信号奇异值与噪声信号奇异值的分解点k₁，如下所示：

k₁＝i|q_i+1＜T and q_i+2＜T and q_i+3＜T i＝1,2,…,r-3 (3)

其中，k₁为目标信号奇异值与噪声信号奇异值的分解点，q_i+1、q_i+2和q_i+3为相邻的三个奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N；

进一步地，在步骤S106中，利用所述分解点k₁及在所述分解点k₁之前的目标信号奇异值进行重构，得到去噪后的目标信号的Hankel矩阵：

其中，B_s为去噪后的目标信号的Hankel矩阵，u_i∈R^m×1，v_i∈R^n×1，σ_i为所述Hankel矩阵B的奇异值。

本发明提供的技术方案带来的有益效果是：本发明根据目标信号和噪声信号的Hankel矩阵奇异值分布差异，利用奇异值差分谱的均值作为阈值自动确定目标信号奇异值，计算简便，阈值稳定性好，可有效抑制探地信号中的噪声。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法的流程图；

图2是本发明实施中所述探地雷达的探测模型示意图；

图3是本发明实施中所述探地雷达得到的B扫描图像示意图；

图4是本发明实施中含噪图像的示意图；

图5是本发明实施中第40道数据构造的Hankel矩阵分解奇异值和奇异值差分谱示意图；

图6是本发明实施中第40道含噪数据和去噪后数据示意图；

图7是本发明实施中去噪后的结果示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法。

请参考图1，图1是本发明实施例中一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法的流程图，具体包括以下步骤：

S101：输入探地雷达B扫描图像X∈R^M×N，其中M为道数，N为每道数据的采样点数；

S102：根据所述探地雷达B扫描图像的一道数据，构造一个Hankel矩阵；设所述一道数据为x＝{x₁,x₂,…,x_N}，对x＝{x₁,x₂,…,x_N}构造Hankel矩阵B，如公式(1)所示：

其中，1＜n＜N,m＝N-n+1,则B∈R^m×n；

依据所述探地雷达回波信号的组成，可以将Hankel矩阵B表示为公式(2)：

B＝B_s+B_n (2)

其中，B_s表示目标信号构成的Hankel矩阵，B_n表示噪声信号构成的Hankel矩阵；

S103：对所述Hankel矩阵进行奇异值分解，如公式(3)所示：

B＝USV^T (3)

其中，B为所述Hankel矩阵，U∈R^m×m和V∈R^n×n分别是由BB^T和B^TB的特征值向量构成的正交矩阵；m＝N-n+1，1＜n＜N，m＞n时，S＝[diag(σ₁,σ₂,…,σ_r),0]，m＜n，S＝[diag(σ₁,σ₂,…,σ_r),0]的转置；

设U＝[u₁,u₂,…,u_m]，V＝[v₁,v₂,…,v_n]，其中u_i∈R^m×1，v_i∈R^n×1，则所述Hankel矩阵B可写成公式(4)：

其中，σ₁≥σ₂≥…≥σ_r≥0是所述Hankel矩阵B的奇异值，从大到小排列在矩阵的主对角线上，r＝min(m,n)，m＝N-n+1，1＜n＜N，0为零矩阵；

S104：求取所述奇异值差分谱；目标信号间具有较强的相关性，低阶的目标信号奇异值较大，高阶的目标信号奇异值较小；噪声信号的相关性较小，噪声信号奇异值分布均匀且较小。因此，利用奇异值差分谱区分目标信号与噪声信号的变换特性，利用公式(5)求取所述奇异值差分谱：

q_i＝σ_i-σ_i+1,i＝1,2,…r-1 (5)

其中，q_i为所述奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N，σ_i和σ_i+1为所述Hankel矩阵B的奇异值，且B∈R^m×n；

S105：通过计算得到所述奇异值差分谱的均值，所述均值作为阈值；利用公式(6)求取所述奇异值差分谱的均值T，以所述均值T作为奇异值判断的阈值：

其中，T为所述奇异值差分谱的均值，即为所述阈值，q_i为所述奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N；

S106：根据所述阈值，确定目标信号奇异值与噪声信号奇异值的分界点，并利用所述目标信号奇异值进行重构，得到去噪后的数据；利用相邻的三个奇异值差分谱和所述阈值进行比较，得到目标信号奇异值与噪声信号奇异值的分解点k₁，如公式(7)所示：

k₁＝i|q_i+1＜T and q_i+2＜T and q_i+3＜T i＝1,2,…,r-3 (7)

其中，k₁为目标信号奇异值与噪声信号奇异值的分解点，q_i+1、q_i+2和q_i+3为相邻的三个奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N；

利用所述分解点k₁及在所述分解点k₁之前的目标信号奇异值进行重构，得到去噪后的目标信号的Hankel矩阵：

其中，B_s为去噪后的目标信号的Hankel矩阵，u_i∈R^m×1，v_i∈R^n×1，σ_i为所述Hankel矩阵B的奇异值；

S107：根据步骤S102～S106，对所述探地雷达B扫描图像中的每一道数据进行处理，得到去噪后的探地雷达B扫描图像X'∈R^M×N。

在本发明实施例中，采用时域有限差分(Finite Difference Time Domain,FDTD)方法得到所述探地雷达的数据，请参考图2，图2是本发明实施中所述探地雷达的探测模型示意图，模型参数设置如下：

(1)模型中介质为混凝土，其相对介电常数为6.0，电导率为0.0001S/m，探地雷达中心频率为900MHZ，发射天线的发射脉冲为Ricker子波；

(2)仿真区域宽度为3m，深度为2m，目标为3个半径为0.2m的理想球导体，球心距地表的埋深约0.6m，且3个目标球体的球心所在水平位置分别为0.9m、1.5m和2.1m；

(3)道间距Dx为0.035m，一共包含80道数据，每道数据有2036个采样点。

请参考图3，图3是本发明实施中所述探地雷达得到的B扫描图像示意图，一共包含80道数据，道间距Dx为0.035m，每道数据有2036个采样点。

请参考图4，图4是本发明实施中含噪图像的示意图，对图3所示的B扫描图像添加白噪声，得到含噪图像，所述含噪图像的信噪比为0.9846dB，一共包含80道数据，道间距Dx为0.035m，每道数据有2036个采样点。

请参考图5，图5是本发明实施中第40道数据构造的Hankel矩阵分解奇异值和奇异值差分谱示意图，选取第40道数据，构造Hankel矩阵，对所述Hankel矩阵进行奇异值分解，得到的奇异值如图(a)所示，奇异值差分谱如图(b)所示。

请参考图6，图6是本发明实施中第40道含噪数据和去噪后数据示意图，通过计算得到的奇异值差分谱均值为58.7382，以此为阈值得到目标信号奇异值和噪声信号奇异值的分界点为18。对前18个奇异值进行重构可得目标信号Hankel矩阵，得到第40道含噪数据(a)和第40道去噪后数据(b)，其中第40道含噪数据信噪比为0.9597dB，第40道去噪后数据信噪比为10.0801dB。

请参考图7，图7是本发明实施中去噪后的结果示意图，对B扫描图像的所有道数据进行去噪处理，得到去噪后的结果，其信噪比为9.1896dB。

本发明的有益效果是：本发明根据目标信号和噪声信号的Hankel矩阵奇异值分布差异，利用奇异值差分谱的均值作为阈值自动确定目标信号奇异值，计算简便，阈值稳定性好，可有效抑制探地信号中的噪声。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法，其特征在于：包括以下步骤：

S101：输入探地雷达B扫描图像X∈R^M×N，其中M为道数，N为每道数据的采样点数；

S102：根据所述探地雷达B扫描图像的一道数据，构造一个Hankel矩阵；

S103：对所述Hankel矩阵进行奇异值分解；

S104：求取所述奇异值差分谱；利用公式(1)求取所述奇异值差分谱：

q_i＝σ_i-σ_i+1,i＝1,2,…r-1 (1)

其中，q_i为所述奇异值差分谱，σ_i和σ_i+1为所述Hankel矩阵B的奇异值，且B∈R^m×n，r＝min(m,n)，m＝N-n+1，1＜n＜N；

S105：通过计算得到所述奇异值差分谱的均值，所述均值作为阈值；利用公式(2)求取所述奇异值差分谱的均值T，以所述均值T作为奇异值判断的阈值：

其中，T为所述奇异值差分谱的均值，即为所述阈值，q_i为所述奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N；

S106：根据所述阈值，确定目标信号奇异值与噪声信号奇异值的分界点，并利用所述目标信号奇异值进行重构，得到去噪后的数据；利用相邻的三个奇异值差分谱和所述阈值进行比较，得到目标信号奇异值与噪声信号奇异值的分解点k₁，如下所示：

k₁＝i|q_i+1＜T and q_i+2＜T and q_i+3＜T i＝1,2,…,r-3 (3)

其中，k₁为目标信号奇异值与噪声信号奇异值的分解点，q_i+1、q_i+2和q_i+3为相邻的三个奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N；

利用所述分解点k₁及在所述分解点k₁之前的目标信号奇异值进行重构，得到去噪后的目标信号的Hankel矩阵：

其中，B_s为去噪后的目标信号的Hankel矩阵，u_i∈R^m×1，v_i∈R^n×1，σ_i为所述Hankel矩阵B的奇异值；

S107：根据步骤S102～S106，对所述探地雷达B扫描图像中的每一道数据进行处理，得到去噪后的探地雷达B扫描图像X′∈R^M×N。

Images