基于EMD和KL变换的时空联合探地雷达杂波抑制

图1 探地雷达接收信号模型

Fig.1 Ground penetrating radar receiving signal model

1.1 时间维处理

探地雷达杂波中对观测影响最大的是噪声,噪声强度过大不仅掩盖隐藏目标体,而且也会降低后续成像的质量。EMD分解具有快速自适应的优点,已经在探地雷达中取得了一定的应用,但由于噪声强度的差异,单纯的EMD分解很难达到理想去噪效果。EMD分解的IMF分量阶数越高,噪声强度越弱,噪声主导的IMF分量并非完全是噪声,可能包含有微弱的有用信号。噪声强度越强,噪声IMF分量越多,则包含的有用信号越多,如果直接去除噪声IMF分量,必定会造成有用信号丢失,如果去除不足,达不到很好的去噪效果。时间维的处理主要针对强随机噪声,利用EMD方法分离探地雷达每道A-scan回波信号中的高频噪声成分,通过能量的方法区分含噪分量和信号分量,对含噪分量引入阈值处理方法克服直接去除的缺点,达到抑制噪声同时保留有用信号细节的目的,最大限度保留信号中的有用成分。

1.1.1 经验模态分解

经验模态分解是Hilbert-Huang变换的重要组成部分,其关键问题是如何把一个非线性非平稳序列分解为有限的IMF分量和一个趋势项。原始信号x(t)经过EMD分解^[4]表示为:

(1)x

(t)

\overset{n}{\sum_{i = 1}}

c_i

(t)

+r_n

(t)

式中:c_i(t)为n个IMF分量,r_n $(t)$ 为残余分量。EMD分解的优点:直观、基于经验资料以及适应性强。其分解快速而且有效,同时它又是基于信号的局部变化特性的,所以适合于探地雷达信号的分析。

1.1.2 噪声与信号分界处选择

信号经过EMD分解为有限个频率由高到低的IMF分量,其中低阶的IMF对应于信号的高频部分,一般为尖锐部分或噪声;高阶IMF对应于信号的低频部分,基本不含有噪声。Wu等^[12]通过大量实验研究指出:IMF分量的能量密度与其平均周期的乘积为一个常数,能量密度与平均周期成反比,即阶数越高能量越弱。基于文献[12]的思想,定义第i阶IMF分量c_i $(t)$ 的能量为

(2)

E_{i}^{}

\frac{1}{N} \overset{N}{\sum_{t = 1}} {[c_{i}^{} (t)]}^{2}

。

当阶数i=k能量达到最低点时,即认为此处IMF分量为噪声成分占主导模态与有用信号占主导模态的分界点,从第i=k+1阶IMF分量开始均为有用信号占主导模态部分。

1.1.3 阈值函数确定

对含噪的前K阶IMF分量采用下列软阈值函数进行处理:

(3)$IMF_{i}^{\prime}=\begin{cases}IMF_{i}(t)-Thr_{i}&\text{若 }IMF_{i}(t)\geq Thr_{i}\\0&\text{若}|IMF_{i}(t)|<Thr_{i}\\IMF_{i}(t)+Thr_{i}&\text{若 }IMF_{i}(t)\leq-Thr_{i}&&\end{cases}$ 。

根据Donoho^[13]提出的噪声阈值选取方法,相应的阈值采用:

(4)Thr_i=σ_i

\sqrt{2 l n (N)}

式中:N为各IMF分量的长度,σ_i为IMF分量中所含噪声分量的方差,基于文献[13]的思想,定义σ_i为

(5)σ_i=MAD_i/0.6745,

式中:MAD_i代表第i阶IMF分量的绝对中值偏差,定义为:

(6)MAD_i=
Median

\{a b s [I M F_{i} (t) - M e d i a n (I M F_{i} (τ))]\}

。

阈值的选取是去噪效果的关键,对处理效果起决定作用的主要是噪声分量的方差σ_i。考虑到EMD分解后各个IMF分量的含噪是越来越小的,若对各个IMF分量采用固定的阈值会造成过处理的发生,导致有用信号的丢失。所以将选取阈值的大小与EMD分解的阶数进行关联,由于阶数越大噪声越弱即二者成反比,所以采用新的阈值公式:

(7)Thr_i=σ_i

\sqrt{2 l n (N)}

/i 。

新的公式中,对于越高阶的IMF分量,相应的阈值越小,与阶数越高噪声越小的特点相吻合,在保证噪声去除的同时更好地实现了有用信号的保留。将噪声主导前k阶IMF通过阈值去噪后得到的结果与有用信号占主导模态的IMF分量及残余分量进行叠加,去噪后的重构信号为:

(8)y(t)=

\overset{k}{\sum_{i = 1}}

F_{i}^{'}

\overset{n}{\sum_{j = k + 1}}

IMF_j+r_n(t) 。

1.2 空间维处理

探地雷达数据采集一般采用剖面法,发射天线和接收天线以固定间距沿着侧线同步移动,可以认为相邻的测点间的目标回波信号具有较强的相关性;随机干扰是由各种不可知因素产生的,没有统一的规律,在整个记录上随机出现,因此认为随机干扰之间不具有相关性。空间维处理采用KL变换从数据中提取相关信号,消除时间维处理残留的随机噪声和干扰。探地雷达回波数据的相关性由KL变换后的协方差矩阵的特征值来度量,选取特征值大的成分重建的信号相关性强,选取特征值小的成分重建的信号相关性弱。在噪声强度较大的情况下,KL变换难以取得较好的压制杂波效果,但经过时间轴降噪处理后的残留噪声干扰强度较弱,满足KL变换的条件,选取大的特征值成分重建信号,去除随机的杂波干扰。KL变换原理文献[14]中有详细描述,这里不再赘述。

2 实验与结果

2.1 仿真验证

单道探地雷达数据进行EMD分解。图2a为GPRmax产生的探地雷达典型回波数据,主要组成为100采样点处的直达波,强度最大,350采样点附近目标物体,设置的为金属铁块的回波(图中框出的部分)。加入噪声后如图2b,由于目标回波信号强度弱,噪声在很多位置处的强度已经与信号相似。

图2

图2 原始探地雷达信号(a)及加噪处理后信号(b)

Fig.2 Original ground penetrating radar signal(a) and signal after noise processing(b)

对加噪后的信号进行EMD分解,产生10个IMF分量,残余的2个分量省略如图3。根据图中各个本征模态函数的曲线可以看出:从IMF₁到IMF₈,噪声强度是依次降低的,其中IMF₁和IMF₂基本为加入的白噪声,说明EMD分解具有分离高频噪声的功能,去除掉IMF₁和IMF₂就可达到一定的去噪目的。IMF₃主要为直达波分量。由IMF₄可见,在350采样点目标体回波处波形明显,说明EMD分解特定的IMF分量具有突出异常信号的能力。由仿真可见EMD非常适合探地雷达数据处理,可以在不同分量上对信号进行解释。

图3

图3 EMD分解后的IMF₁~IMF₈

Fig.3 IMF₁~IMF₈ after EMD decomposition

为了更精确还原原始回波数据,对噪声分量进行阈值处理操作。计算各IMF相应的能量,最低点位于6,即前六阶IMF分量为噪声主导分量,对前六阶IMF分量进行阈值处理,之后与后四阶IMF分量进行重组,处理结果如图4a,噪声基本被滤除, 350采样点处信号波形还原准确。与经典的小波去噪图4b进行对比,小波基函数采用db6小波,四级软阈值去噪, EMD阈值重组后的信号波形更加平稳,对噪声去除更加彻底,对有用信号保留完整,与原始的信号更接近,去噪效果更加理想。

图4

图4 EMD改进阈值去噪(a)与小波去噪(b)对比

Fig.4 Comparison between EMD improved threshold denoising(a) and wavelet denoising(b)

信噪比和均方误差通常是衡量去噪效果的定性指标,其公式如下:

(9)信噪比:SNR=10lg

(\frac{\overset{N}{\sum_{n = 1}} x^{2} (n)}{\overset{N}{\sum_{n = 1}} {[\hat{x} (n) - x (n)]}^{2}})

(10)均方根误差:RMSE=

{[\frac{1}{N} \overset{N}{\sum_{n = 1}} {[\hat{x} (n) - x (n)]}^{2}]}^{\frac{1}{2}}

式中:x $(n)$ 为原始信号, $\hat{x} (n)$ 为估计信号。

计算加噪后信号、小波阈值去噪后信号、EMD阈值去噪后信号的信噪比及均方根误差如表1所示。从表中数据可见EMD阈值处理后的信噪比更高,均方根误差更小,定量的描述也可看出EMD方法性能更优越。由于噪声强度较大,如500采样点(图4a中小框标注)处还存在类似信号的回波,单纯依靠时间维处理很难将其彻底去除,因此引入空间维的KL变换对残余的这些干扰进行压制。

表1 信噪比与均方根误差对比

Table 1 Comparison of SNR and RMSE

	SNR	RMSE
加噪信号	10.058	0.0099
小波阈值去噪	18.9591	0.0036
EMD阈值去噪	20.1949	0.0031

新窗口打开| 下载CSV

2.2 实测数据验证

对于实际搭建的探地雷达系统采集的数据进行去噪处理。实验主要有2个,实验一选取单道时间轴数据,数据噪声强烈,主要为了验证EMD阈值去噪的效果;实验二选取多道数据,数据时间维噪声相对较弱但空间维杂波干扰较强烈,主要为了验证时空联合杂波抑制效果。

2.2.1 实验一

自行搭建的探地雷达系统如图5a,主要由脉冲源、天线、安捷伦数据采集卡和低噪放等组成。脉冲源为二阶高斯脉冲,6dB带宽320 MHz~1.56 GHz;发射接收天线型号ATB1G-12-003,均为极化喇叭天线;采集卡型号U1065A,采样频率8 GHz;低噪放型号ADL5523。实验的对象为水泥板,收发分置天线间距35.5 cm,离地20 cm,天线方向垂直朝下。脉冲源产生脉冲信号经过发射天线发送,照射到水泥板后回波被接收天线接收,经过采集卡保存在电脑端,每个测点数据都经过多次采集之后取平均,减少人为误差的影响,去除一部分随机噪声,最后经过归零处理,结果如图5b。

图5

图5 实测场景及数据处理

a—实测场景;b—实测采集回波数据;c—EMD改进阈值去噪;d—小波去噪

Fig.5 Actual measurement scenario and data processing

a—actual measurement scenario;b—measured echo data;c—EMD improved threshold denoising;d—wavelet denoising

回波数据虽然经过统计处理,但噪声依旧非常明显,除直达波由于强度比较大显示的较明显外,其他位置处的回波信息基本淹没在噪声中。对数据进行EMD阈值处理结果如图5c所示。处理后的结果可见,随机噪声被压制,波形平稳,直达波分量处的强度值基本没有改变。在不同位置多次采集数据处理后发现,水泥板的回波与直达波基本叠加在一起,400采样点处的异常并非是地下目标的回波,而是信号源产生。这是信号源存在的问题,但可以作为衡量去噪效果的标准。与小波去噪图5d对比,EMD阈值去噪效果更明显,对400采样点处信号源的异常还原更准确和清晰,能有效地保留信号中的突变成分,具有更好的细节区分能力。

2.2.2 实验二

主要实验仪器:脉冲源,低噪放,安捷伦示波器,喇叭天线,数据经过安捷伦示波器进行显示并采集保存。示波器型号ium 80000,采样频率40 GHz,其他设备与试验一相同。实验场景:沙坑中埋有一个拱形的铁片,埋深5 cm,铁片圆面的半径26.5 cm,弦长40 cm,收发分置天线间距35.5 cm,离地10 cm。实验中天线每次移动2 cm,总计采集41次,每一次采集多个周期取平均减少随机误差影响。采集的原始数据组成B-scan如图6a。原始信号随机噪声较弱但干扰厉害,直达波、杂波突出,导致铁皮的回波信号不明显,基本无法识别。使用减平均法去除直达波,结果如图6b。去除直达波后,一部分水平杂波被去除,主要为地表直达波和周围墙体可能产生的水平杂波干扰,圆形铁皮的形状已经基本能分辨出来,铁皮回波主要集中在5~9 ns之间,但是剩余干扰和噪声依旧很强。对每道信号进行EMD阈值去噪重组,结果如图6c,整个剖面部分噪声被去除,剩余的杂波干扰也得到了很大的抑制,例如水平位置0.2~0.4处抑制非常明显,对目标铁皮回波信号基本没有损伤,目标更加突出。对于处理后还存在的干扰主要原因可能是沙坑没有仔细抹平,一些地方突出比较明显,这些干扰不满足噪声去除条件,但是其具有非相关性,在整个雷达剖面中不是连续出现。经过降噪处理后的信号最后进行KL变换处理,选取前2个大的特征值重建信号去除残余的干扰,结果如图6d。由图可见,剩余的杂波干扰基本被压制,整个B-scan剖面中目标体的回波清晰。

图6

图6 原始数据及数据处理结果

a—原始数据;b—去除直达波数据;c—EMD改进阈值去噪;d—时空二维杂波压制

Fig.6 Raw data and data processing results

a—original data;b—remove direct wave data;c—EMD improved threshold denoising;d—time-space two-dimensional clutter suppression

根据文献[1]图像熵标准衡量杂波抑制效果,熵越小,图像越清晰,表示含目标可能性越大。图像熵公式如式(11),m为测点数,n为时间采样点数。

(11)Q=

\frac{{[\overset{m}{\sum_{i = 1}} \overset{n}{\sum_{j = 1}} x_{i j}^{2}]}^{2}}{\overset{m}{\sum_{i = 1}} \overset{n}{\sum_{j = 1}} x_{i j}^{4}}

。

计算去直达波后的回波、EMD阈值处理、KL变换处理后的Q值如表2所示,说明了时间轴和空间轴对杂波的压制都取得了作用,经过时空二维处理后的图像熵值只有原始的一半,说明了时空处理方法的有效性。

表2 图像熵值对比

Table 2 Comparison of image entropy

	去直达波	EMD阈值	KL变换
Q	2090	1553	1099

新窗口打开| 下载CSV

3 结论

随着探地雷达应用环境越来越复杂,单一杂波抑制方式难以满足地下解释和后续处理要求。本文提出的时空二维联合去杂波方法互为补充,EMD阈值处理能有效去除强噪声,克服KL变换的缺点;KL变换能够进一步去除EMD阈值处理残余的非相关干扰,获得清晰的回波数据。EMD阈值去噪克服了EMD去噪对强噪声不理想的特点,根据IMF含噪特点引入随阶数变化的阈值函数,其具有自适应特点,不需要像小波变换那样根据情况选择小波基。两组实测实验采集场景和方式不同,得到的数据随机噪声强度和杂波类型不同,验证了EMD阈值去噪以及时空二维联合去杂波是有效的。通过杂波压制,雷达剖面更加清晰,有利于进一步的成像识别等处理。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

申家全, 闫怀志, 胡昌振.

基于图像熵的探地雷达杂波抑制效果评价

[J]. 电波科学学报, 2011, 26(2):267-271,413.

Shen

J Q

, Yan

H Z

, Hu

C Z

Evaluation method based on image entropy for clutter suppression effect of ground-penetrating radar

[J]. Chinese Journal of Radio Science, 2011, 26(2):267-271,413.

[2]

石显新, 杨秋芬, 侯彦威.

二维物理小波标架去除探地雷达信号随机噪声

[J]. 电波科学学报, 2012, 27(6):1186-1192,1268.

Shi

X X

, Yang

Q F

, Hou

Y W

Eliminating random noise interference of ground penetrating radar signal by resorting to 2-D physical wavelet fame

[J]. Chinese Journal of Radio Science, 2012, 27(6):1186-1192,1268.

[3]

吴宝杰, 杨桦, 张伟光.

探地雷达数据的S变换时频分析

[J]. 上海地质, 2008, 107(3):13-15,19.

B J

, Yang

, Zhang

W G

S transform time-frequency analysis of GPR data

[J]. Shanghai Geology, 2008, 107(3):13-15,19.

[4]

Huang

N E

, Wu

M L C

, Long

S R

et al.

A confidence limit for the empirical mode decomposition and Hilbert spectral analysis

[J]. Proceedings of the Royal Society of London Series A:Mathematical,Physical and Engineering Sciences,2003, 459(2037):2317-2345.

[5]

Flandrin

, Rilling

, Goncalves

Empirical mode decomposition as a filter bank

[J]. IEEE Signal Processing Letters, 2004, 11(2):112-114.

DOI:10.1109/LSP.2003.821662 URL [本文引用: 1]

[6]

易洪春, 刘树才, 贺克升,

等.

基于 EMD 去噪方法研究及其在地质勘探中的应用

[J]. 物探与化探, 2013, 37(3):533-537.

H C

, Liu

S C

, He

K S

et al.

Data denoising based on EMD and its ap-plication to geological exploration

[J]. Geophysical and Geochemical Explora-tion, 2013, 37(3):533-537.

[7]

冯德山, 戴前伟, 余凯.

基于经验模态分解的低信噪比探地雷达数据处理

[J]. 中南大学学报:自然科学版, 2012, 43(2):596-604.

Feng

D S

, Dai

Q W

, Yu

GPR signal processing under low SNR based on empirical mode decomposition

[J]. Journal of Central South University:Science and Technology, 2012, 43(2):596-604.

[8]

Tesfamariam, Gebremichael

, Mali

et al.

Clutter reduction techniques for GPR based buried landmine detection

[C]// Thuckalay:2011 International Conference on Signal Processing,Communication,Computing and Networking Technologies.IEEE, 2011:182-186.

[9]

Khan

U S

, Al-Nuaimy

Background removal from GPR data using Eigenvalues

[C]// Lecce: Proceedings of the XIII Internarional Conference on Ground Penetrating Radar.IEEE, 2010:1-5.

[10]

申家全, 闫怀志, 胡昌振.

基于主成分自动选择准则的探地雷达杂波抑制

[J]. 电波科学学报, 2010, 25(1):83-87.

Shen

J Q

, Yan

H Z

, Hu

C Z

Auto-selected rule on principal component analysis in ground penetrating radar signal denoising

[J]. Chinese Journal of Radio Science, 2010, 25(1):83-87.

[11]

Solimene

, Cuccaro

, Dell’Aversano

et al.

Ground clutter removal in GPR surveys

[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(3):792-798.

DOI:10.1109/JSTARS.4609443 URL [本文引用: 1]

[12]

Z H

, Huang

N E

A study of the characteristics of white noise using the empirical mode decomposition method

[J]. Proceedings of the Royal Society of London Series A:Mathematical,Physical and Engineering Sciences,2004, 460(2046):1597-1611.

[本文引用: 3]

[13]

Donoho

D L

De-noising by soft-thresholding

[J]. IEEE Transactions on Information Theory, 1995, 41(3):613-627.

DOI:10.1109/18.382009 URL [本文引用: 3]

[14]

杨天春, 张辉, 李好.

探地雷达信号的 KL 变换处理分析

[J]. 湖南科技大学学报:自然科学版, 2013, 28(2):66-69.

Yang

T C

, Zhang

, Li

Signal processing of ground penetrating radar by KL transform

[J]. Journal of Hunan University of Science & Technology:Natural Science Edition, 2013, 28(2):66-69.