CN107121705A - 一种基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法 - Google Patents

Abstract

本发明公开了一种基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法，包括以下步骤：将探地雷达原始回波信号与相同长度的零均值白噪声信号进行随机拟合，得到两道信号；运用独立成分分析算法对这两道信号进行分离处理，输出峰度值较大的信号记为x′(m)，峰度值较小的信号记为n′(m)；对x′(m)进行自动反相校正，再运用完全总体经验模态算法对反相校正后的信号进行分解得到P个信号分量，计算各信号分量的峰度值；再计算信号n′(m)的峰度值作为阈值；最后将峰度值大于阈值的所有信号分量累加，作为原始回波信号x(m)去噪后的信号。本方法提高了计算效率和去噪效果。

Description

一种基于自动反相校正和峰度值比较的探地雷达回波信号去 噪算法

技术领域

本发明属于探地雷达探测与应用技术领域，具体涉及探地雷达回波信号去噪处理。

背景技术

探地雷达(Ground Penetrating Rada，GPR)通过发射天线向地下发射宽带电磁波，接收天线接收散射回波，通过对散射回波进行处理，对地下未知区域进行无损探测和参数反演。电磁波在地下介质中传播时，遇到存在电性差异的分界面时发生散射，根据接收到的电磁散射回波来反演地下未知区域的异常体位置、形态、埋藏深度等参数。GPR作为探测地下目标的一种重要工具，在人道主义救援中探测地下生命迹象，在维和行动中定位清除地雷、未爆炸武器，在安全保卫中探测判断地下异常物，在施工中定位地下管道、线缆等方面正得到越来越多的研究与使用。

GPR回波信号由直达波、界面反射波、目标散射波、随机噪声等成分构成，在频率域和时间域互有重叠，难以区分。GPR发送的电磁波在地下传播的过程中，受到地下介质结构复杂多变、仪器参数或者噪声的多种因素的影响，会发生衰减、频散和其他干扰，这在很大程度上限制了GPR的探测分辨率和数据解释效果。因此，为了获得质量良好的回波信号，需要抑制噪声和杂波，提取所需的目标回波信号，最大限度地还原地下的异常体的散射回波。GPR信号处理水平对定位和识别目标起着决定性作用，因此开展GPR信号处理方法的研究工作有着重要的研究价值与现实意义。

国内外开展了大量的GPR信号去噪处理方面的研究工作。文献1“Jing Li,CaiLiu,Zhaofa Zeng,Lingna Chen,GPR Signal Denoising and Target Extraction Withthe CEEMD Method,IEEE Geoscience and Remote Sensing Letters,2015,12(8):1615-1619.”提出了一种基于完全总体经验模态分解(Complete Ensemble Empirical ModeDecomposition,CEEMD)的GPR信号处理方法，通过对GPR回波信号进行CEEMD分解，人工方式进行噪声辨识和去除，进而实现信号去噪。运用Hilbert-Huang变换，证明了该算法比经验模态分解和总体经验模态分解的谱分辨率更高；文献2“谢朋成,探地雷达探测浅地层目标的算法研究,东北林业大学,2016”提出了一种基于独立成分分析(Independent ComponentAnalysis，ICA)的GPR探测浅地层目标的算法研究。文献3“陈玲娜,基于CEEMD和PCA的探地雷达数据处理研究与应用,吉林大学,2016”对复杂GPR信号进行分频处理，通过选择固定频段的本征模态函数(Intrinsic Mode Function,IMF)进行重构从而去除高频杂波和低频隐失波的影响，再通过从原始信号中去掉高阶主成分信息，提取目标信号。上述的处理方法中，基于ICA的去噪算法存在信号相位不定性的问题，基于CEEMD的去噪算法则需要通过人工判读的方式对各IMF分量进行判别，再进行信号重构，效率低。

因此，有必要设计一种能够自动判定相位并自动进行IMF分量提取的GPR信号去噪算法。

发明内容

本发明所要解决的技术问题是，针对ICA算法的信号相位不定性和CEEMD算法中问题，提供一种基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法，提高了计算效率和去噪效果。

本发明的技术方案如下：

一种基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法，包括以下步骤：

步骤1：首先将探地雷达单道含噪原始回波信号x(m),m＝1,2,...,M与相同长度的零均值白噪声信号n(m),m＝1,2,...,M进行随机拟合，得到两道信号；然后根据独立成分分析算法的分离特性，运用独立成分分析算法对这两道信号进行分离处理，输出长度为M的峰度值不相等的两道信号；最后分别计算两道信号的峰度值，峰度值较大的信号记为x′(m),m＝1,2,...,M，峰度值较小的信号记为n′(m),m＝1,2,...,M；

步骤2：判断x′(m)与x(m)是否反相；如果反相，则置相位判别因子α＝-1，否则置α＝1；得到信号α·x′(m)(即若反相则进行自动反相校正，如果不反相，则保持不变)；

步骤3：首先运用完全总体经验模态算法对α·x′(m)进行分解，得到P个信号分量，记为y_p(m),p＝1,...,P；m＝1,...,M；然后计算各信号分量的峰度值，分别记为k_p,p＝1,...,P；再计算信号n′(m)的峰度值，记为k；最后以k作为阈值，将峰度值k_p大于k的所有信号分量累加，作为原始回波信号x(m)去噪后的信号，记为

进一步地，所述步骤1中，零均值白噪声信号n(m)的幅度和方差任意设置。

进一步地，所述步骤1中，随机拟合的方法为：首先将x(m)和n(m)组合成2×M的矩阵，记为K，该矩阵的第一行为x(m)，第二行为n(m)，M表示x(m)的长度；然后生成一个2×2的矩阵记为L；最后将L左乘K，得到矩阵U，其第一行数据记为U₁，第二行数据记为U₂；U₁和U₂即为随机拟合得到的两道信号。

进一步地，所述步骤1和步骤3中，记某道信号为r(m),(m＝1,2…,M)，其峰度值计算公式如下：

其中，代表r(m),(m＝1,2…,M)的均值。

进一步地，所述步骤2中，根据max(|x′(m)-x(m)|)>max(|x(m)|),m＝0,...,M-1是否成立，判断x′(m)与x(m)是否反相；如果成立，说明x′(m)与x(m)反相；否则说明x′(m)与x(m)不反相；式中，|·|表示对一维数组x′(m)-x(m)中各元素取绝对值，max(|x(m)|)表示取一维数组|x(m)|中各元素的最大值。

有益效果：

本发明提出了一种基于自动反相校正和峰度值比较的GPR回波信号的去噪算法，针对独立成分分析算法分离后信号相位不定性问题，设计了相位判别因子，实现了独立成分分析算法分解后的信号相位自动判别和校正。针对完全总体经验模态分解算法分解后，各IMF分量需要人工判别的问题，设计了基于峰度值比较的IMF分量自动筛选，阈值选择为ICA算法分离出来的噪声信号的峰度值。本发明避免了了独立成分分析算法分解后的相位不定性，且在总体经验模态分解算法分解后无需传统的人工方式进行各IMF分量的筛选，提高了计算效率和去噪效果。

附图说明

图1示出了本方法的流程图。

图2示出了GPR正演模型图。

图3示出了图2正演模拟得到的GPR无噪回波信号图。

图4示出了人为加入噪声后模拟得到的含噪GPR信号图和产生的等长度的随机噪声信号；图4(a)为含噪GPR信号，图4(b)为根据图4(a)含噪GPR信号随机产生的等长度的随机噪声信号。

图5示出了图4所示的两个信号随机拟合后得到的两道信号，5(a)和(b)分别为随机拟合后得到的两道信号。

图6示出了图5所示的两道信号经过ICA算法后分离出的两道信号。

图7示出了图6中的第一道信号经过CEEMD分解后的各IMF分量波形图，为图7(a)为IMF5～IMF8，7(b)为IMF9～IMF12，7(c)为IMF13～IMF15。

图8示出了各IMF分量的峰度值计算结果和ICA分解得到的噪声信号的峰度值阈值。

图9示出了图7中的各IMF分量采用峰度值阈值比较后重构的信号。

图10示出了本方法和常规算法的去噪误差和信噪比的变化曲线对比图

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明。

GPR正演模拟如图2所示，地下介质有三层，第一二层的厚度均为15cm，第一二三层的相对介电常数分别为20，10和15，一个半径为4cm的金属管埋在第二层中，金属管的中心距第二层上表面6cm，发射和接收天线位于金属管中心的正上方，距地表的高度为5cm。发射信号为中心频率为900MHz的Ricker子波，时窗取为40ns，采用时域有限差分法仿真得到接收天线所接收的地下区域的散射回波，如图3所示，时长为40ns，采样点数为6784，该信号为无噪的回波信号s(m)。为模拟生成含噪的GPR信号，人为地加入了随机白噪声，白噪声的幅值和方差分别为4.3328和1.0513，加入噪声后的GPR回波信号如图4(a)所示，SNR＝19。下面对该信号运用本发明的去噪算法进行去噪处理。根据图4(a)含噪GPR信号随机产生的等长度的随机噪声信号如图4(b)所示，再将图4(a)含噪GPR信号和图4(b)随机噪声信号进行随机拟合，即与一个2×2的矩阵相乘，得到两道混合信号，随机拟合结果如图5(a)和(b)所示。图5(a)和(b)输入ICA算法进行处理，得到两道输出信号，该两道信号的峰度值为81.6651和2.9763，峰度值高的输出信号记为x′(m)，如图6(a)所示，峰度值低的输出信号记为n′(m)，如图6(b)所示，对图6(a)信号进行判断是否信号反相，如果反相则对信号进行自动相位校正，如果未反相则信号不变。此处，经过公式判断，未发生反相，则α＝1。再对α·x′(m)进行CEEMD算法的处理，得到图7所示的15个IMF分量。对图7中的每一个IMF分量计算峰度值，峰度值分别为[2.1453,2.1983,2.7546,2.8556,75.0391,41.8097,25.3151,16.4506,9.3599,7.3725,3.2763,3.7535,3.9701,1.6372,1.7456]，如图8所示。如前所述，ICA算法分离出来的噪声信号n′(m)，如图6(b)所示，其峰度值为2.9763。将图7中的每一个IMF分量的峰度值与图6(b)噪声信号的峰度值进行比较，将峰度值低于2.9763的IMF分量剔除，即剔除IMF1，IMF2，IMF3，IMF4，IMF14，IMF15分量，保留IMF5-IMF13分量。将IMF5-IMF13分量进行累加得到的重构信号作为去噪后的GPR回波信号z(m)，如图9所示。将该信号与图3所示的无噪GPR回波信号相对比，并运用公式计算误差，得到均方误差为0.001085。

为定量评估该方法的去噪算法与常规CEEMD的去噪算法的性能，生成了不同SNR的含噪GPR回波信号，SNR的设置范围为[0,20]，点数为21点。每一SNR情况下，根据图3的无噪的GPR回波，生成对应幅度的随机噪声信号，并将二者相加，作为含噪的GPR回波。分别运用本专利申请的算法和常规CEEMD算法对该回波进行去噪处理，得到各自的去噪回波。再分别与图3所示的无噪GPR回波进行对比分析，使用公式计算均方误差值。遍历整个SNR设置区间，得到各SNR情况下的均方误差曲线，如图10所示。

从图10中可见，本专利申请的去噪算法和常规CEEMD去噪算法随着SNR的增大，均方误差都随之降低。但在同一SNR情况下，本专利申请的去噪算法的均方误差更低，去噪效果更好。而且，本专利申请的去噪算法无需常规CEEMD去噪算法的人工方式剔除噪声分量的步骤，全程自动处理，效率更高。

Claims (4)

1.一种基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法，其特征在于，包括以下步骤：

步骤1：首先将探地雷达单道含噪原始回波信号x(m),m＝1,2,...,M与相同长度的零均值白噪声信号n(m),m＝1,2,...,M进行随机拟合，得到两道信号；然后运用独立成分分析算法对这两道信号进行分离处理，输出长度为M的峰度值不相等的两道信号；最后分别计算两道信号的峰度值，峰度值较大的信号记为x′(m),m＝1,2,...,M，峰度值较小的信号记为n′(m),m＝1,2,...,M；

步骤2：判断x′(m)与x(m)是否反相；如果反相，则置相位判别因子α＝-1，否则置α＝1；得到信号α·x′(m)；

步骤3：首先运用完全总体经验模态算法对α·x′(m)进行分解，得到P个信号分量，记为y_p(m),p＝1,...,P；m＝1,...,M；然后计算各信号分量的峰度值，分别记为k_p,p＝1,...,P；再计算信号n′(m)的峰度值，记为k；最后以k作为阈值，将峰度值k_p大于k的所有信号分量累加，作为原始回波信号x(m)去噪后的信号，记为

2.根据权利要求1所述的基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法，其特征在于，所述步骤1中，零均值白噪声信号n(m)的幅度和方差任意设置。

3.根据权利要求1所述的基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法，其特征在于，所述步骤1中，随机拟合的方法为：首先将x(m)和n(m)组合成2×M的矩阵，记为K，该矩阵的第一行为x(m)，第二行为n(m)，M表示x(m)的长度；然后生成一个2×2的矩阵记为L；最后将L左乘K，得到矩阵U，其第一行数据记为U₁，第二行数据记为U₂；U₁和U₂即为随机拟合得到的两道信号。

4.根据权利要求1所述的基于自动反相校正和峰度值比较的探地雷达回波信号去噪算法，其特征在于，所述步骤2中，根据max(|x′(m)-x(m)|)>max(|x(m)|),m＝0,...,M-1是否成立，判断x′(m)与x(m)是否反相；如果成立，说明x′(m)与x(m)反相；否则说明x′(m)与x(m)不反相；式中，|·|表示对一维数组x′(m)-x(m)中各元素取绝对值，max(|x(m)|)表示取一维数组|x(m)|中各元素的最大值。