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

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

【技术实现步骤摘要】
一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法
本专利技术涉及数字信号处理领域，具体涉及到探地雷达的噪声抑制问题，尤其涉及一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法。
技术介绍
探地雷达基于电磁波传播与散射原理，通过向地下发射电磁波信号并接收地下介质不连续处散射的回波实现对地下目标的探测。与电阻率法、低频电磁感应法和地震法等探测方法相比，探地雷达具有探测速度快、探测过程连续、分辨率高、操作方便灵活、费用低、探测范围广(能探测金属和非金属)等优点，在地质、资源、环境、工程和军事等领域得到广泛的应用。在对地下目标进行探测时，受地下杂波和周围介质的干扰等，探地雷达回波信号通常是一种弱信噪比的非平稳信号，因此如从噪声背景中提取目标信息，是探地雷达信号处理的一个重要研究领域。目前常用的探地雷达噪声抑制方法有傅立叶变换、小波变换、S变换和主成分分析法等。傅立叶变换只能反映信号的整体特征,不适用于频率随时间变化的非平稳信号；小波变换中小波基函数、分解层数以及阈值的选择都依赖于主观经验，缺乏自适应性；S变换中时频滤波器的设计较为复杂，也限制了其应用场合。主成分分析法(PrincinalComponetAnalysis，PCA)是一种建立在最小均方误差基础上的线性变换处理方法，其算法核心是用特定的正交矩阵对信号矩阵进行正交变换，得到相互正交的对角主成分矩阵。奇异值分解算法是PCA方法中的常用算法，在该方法中，将奇异值对应的信号分量称为主成分，主成分的选择通常是以主成分对应的奇异值选择为基础的。PCA方法去噪的原理就是选择目标信号的奇异值重构，以去除噪声信号。目前，在传统的探地雷达PCA去噪方法中，一般是对探地雷达的B扫描图像进行奇异值分解，然后选择目标信号奇异值重构。选择奇异值的方法主要有定性的经验方法和特征值能量百分比方法。经验方法根据经验确定目标信号的奇异值，容易受到人的主观性影响；特征值能量百分比方法可确定目标信号奇异值，但该方法比较复杂，阈值稳定性差，受具体探测条件影响大。因此，如何以较小的复杂度准确选择目标信号奇异值，对于提高探地雷达的噪声抑制性能具有重要意义。
技术实现思路
为了解决上述问题，本专利技术提供了一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法，主要包括以下步骤：S101：输入探地雷达B扫描图像X∈RM×N，其中M为道数，N为每道数据的采样点数；S102：根据所述探地雷达B扫描图像的一道数据，构造一个Hankel矩阵；S103：对所述Hankel矩阵进行奇异值分解；S104：求取所述奇异值差分谱；S105：通过计算得到所述奇异值差分谱的均值，所述均值作为阈值；S106：根据所述阈值，确定目标信号奇异值与噪声信号奇异值的分界点，并利用所述目标信号奇异值进行重构，得到去噪后的数据；S107：根据步骤S102～S106，对所述探地雷达B扫描图像中的每一道数据进行处理，得到去噪后的探地雷达B扫描图像X'∈RM×N。进一步地，在步骤S104中，利用公式(1)求取所述奇异值差分谱：qi＝σi-σi+1,i＝1,2,…r-1(1)其中，qi为所述奇异值差分谱，σi和σi+1为所述Hankel矩阵B的奇异值，且B∈Rm×n，r＝min(m,n)，m＝N-n+1，1＜n＜N。进一步地，在步骤S105中，利用公式(2)求取所述奇异值差分谱的均值T，以所述均值T作为奇异值判断的阈值：其中，T为所述奇异值差分谱的均值，即为所述阈值，qi为所述奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N。进一步地，在步骤S106中，利用相邻的三个奇异值差分谱和所述阈值进行比较，得到目标信号奇异值与噪声信号奇异值的分解点k1，如下所示：k1＝i|qi+1＜Tandqi+2＜Tandqi+3＜Ti＝1,2,…,r-3(3)其中，k1为目标信号奇异值与噪声信号奇异值的分解点，qi+1、qi+2和qi+3为相邻的三个奇异值差分谱，r＝min(m,n)，m＝N-n+1，1＜n＜N；进一步地，在步骤S106中，利用所述分解点k1及在所述分解点k1之前的目标信号奇异值进行重构，得到去噪后的目标信号的Hankel矩阵：其中，Bs为去噪后的目标信号的Hankel矩阵，ui∈Rm×1，vi∈Rn×1，σi为所述Hankel矩阵B的奇异值。本专利技术提供的技术方案带来的有益效果是：本专利技术根据目标信号和噪声信号的Hankel矩阵奇异值分布差异，利用奇异值差分谱的均值作为阈值自动确定目标信号奇异值，计算简便，阈值稳定性好，可有效抑制探地信号中的噪声。附图说明下面将结合附图及实施例对本专利技术作进一步说明，附图中：图1是本专利技术实施例中一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法的流程图；图2是本专利技术实施中所述探地雷达的探测模型示意图；图3是本专利技术实施中所述探地雷达得到的B扫描图像示意图；图4是本专利技术实施中含噪图像的示意图；图5是本专利技术实施中第40道数据构造的Hankel矩阵分解奇异值和奇异值差分谱示意图；图6是本专利技术实施中第40道含噪数据和去噪后数据示意图；图7是本专利技术实施中去噪后的结果示意图。具体实施方式为了对本专利技术的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本专利技术的具体实施方式。本专利技术的实施例提供了一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法。请参考图1，图1是本专利技术实施例中一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法的流程图，具体包括以下步骤：S101：输入探地雷达B扫描图像X∈RM×N，其中M为道数，N为每道数据的采样点数；S102：根据所述探地雷达B扫描图像的一道数据，构造一个Hankel矩阵；设所述一道数据为x＝{x1,x2,…,xN}，对x＝{x1,x2,…,xN}构造Hankel矩阵B，如公式(1)所示：其中，1＜n＜N,m＝N-n+1,则B∈Rm×n；依据所述探地雷达回波信号的组成，可以将Hankel矩阵B表示为公式(2)：B＝Bs+Bn(2)其中，Bs表示目标信号构成的Hankel矩阵，Bn表示噪声信号构成的Hankel矩阵；S103：对所述Hankel矩阵进行奇异值分解，如公式(3)所示：B＝USVT(3)其中，B为所述Hankel矩阵，U∈Rm×m和V∈Rn×n分别是由BBT和BTB的特征值向量构成的正交矩阵；m＝N-n+1，1＜n＜N，m＞n时，S＝[diag(σ1,σ2,…,σr),0]，m＜n，S＝[diag(σ1,σ2,…,σr),0]的转置；设U＝[u1,u2,…,um]，V＝[v1,v2,…,vn]，其中ui∈Rm×1，vi∈Rn×1，则所述Hankel矩阵B可写成公式(4)：其中，σ1≥σ2≥…≥σr≥0是所述Hankel矩阵B的奇异值，从大到小排列在矩阵的主对角线上，r＝min(m,n)，m＝N-n+1，1＜n＜N，0为零矩阵；S104：求取所述奇异值差分谱；目标信号间具有较强的相关性，低阶的目标信号奇异值较大，高阶的目标信号奇异值较小；噪声信号的相关性较小，噪声信号奇异值分布均匀且较小。因此，利用奇异值差分谱区分目标信号与噪声信号的变换特性，利用公式(5)求取所述奇异值差分谱：qi＝σi-σi+1,i＝1,2,…r-1(5)其中，本文档来自技高网...

【技术保护点】
1.一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法，其特征在于：包括以下步骤：S101：输入探地雷达B扫描图像X∈RM×N，其中M为道数，N为每道数据的采样点数；S102：根据所述探地雷达B扫描图像的一道数据，构造一个Hankel矩阵；S103：对所述Hankel矩阵进行奇异值分解；S104：求取所述奇异值差分谱；S105：通过计算得到所述奇异值差分谱的均值，所述均值作为阈值；S106：根据所述阈值，确定目标信号奇异值与噪声信号奇异值的分界点，并利用所述目标信号奇异值进行重构，得到去噪后的数据；S107：根据步骤S102～S106，对所述探地雷达B扫描图像中的每一道数据进行处理，得到去噪后的探地雷达B扫描图像X'∈RM×N。

【技术特征摘要】
1.一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法，其特征在于：包括以下步骤：S101：输入探地雷达B扫描图像X∈RM×N，其中M为道数，N为每道数据的采样点数；S102：根据所述探地雷达B扫描图像的一道数据，构造一个Hankel矩阵；S103：对所述Hankel矩阵进行奇异值分解；S104：求取所述奇异值差分谱；S105：通过计算得到所述奇异值差分谱的均值，所述均值作为阈值；S106：根据所述阈值，确定目标信号奇异值与噪声信号奇异值的分界点，并利用所述目标信号奇异值进行重构，得到去噪后的数据；S107：根据步骤S102～S106，对所述探地雷达B扫描图像中的每一道数据进行处理，得到去噪后的探地雷达B扫描图像X'∈RM×N。2.如权利要求1所述的一种基于Hankel矩阵奇异值分解的探地雷达噪声抑制方法，其特征在于：在步骤S104中，利用公式(1)求取所述奇异值差分谱：qi＝σi-σi+1,i＝1,2,…r-1(1)其中，qi为所述奇异值差分谱，σi和σi+1为所述Hankel矩阵B的奇异值，且B∈Rm×n，r＝min(m,n)，m＝N-n+1，1＜n＜N。3.如权利要求2所述的一种基于Hankel矩阵奇异值分解的探...

【专利技术属性】
技术研发人员：薛伟，罗严，朱继超，刘力，戴向阳，
申请(专利权)人：中国地质大学武汉，
类型：发明
国别省市：湖北,42