物探与化探 2018 , 42 (1): 220-224 https://doi.org/10.11720/wtyht.2018.1.28

Orginal Article

探地雷达信号消噪中的时频谱分解重排算法

吴楠, 吴舰, 吴志坚

贵州师范大学机械与电气工程学院,贵州贵安新区 550025

Research on denoising of ground penetrating radar signals using the time-frequency spectral decomposition reassignment algorithm

WU Nan, WU Jan, WU Zhi-Jian

School of Mechanical and Electrical Engineering, Guizhou Normal University, Gui'an 550025, China

中图分类号: P631

文献标识码: A

文章编号: 1000-8918(2018)01-0220-05

责任编辑: WU NanWU JanWU Zhi-Jian

收稿日期: 2017-03-29

修回日期: 2017-07-3

网络出版日期: 2018-01-20

基金资助: 贵州省科技厅科技支撑计划项目(黔科合SY[2010]3023号)

作者简介:

作者简介: 吴楠(1975-),男,硕士,讲师,主要从事网络技术与控制系统、优化算法研究工作。Email:wn_211@sina.com

展开

摘要

为了改善现有的时频域算法对于探地雷达信号消噪效果不佳的现状,本文利用时频谱分解重排算法实现了探地雷达信号的优化消噪处理。由于采用窗口函数对信号进行划分,传统的时频变换方法得到的结果是模糊化的信号时频分布,这使得在时频消噪过程中,如果阈值取不好,信号的有用成份极容易受到损伤。而采用时频谱分解重排算法可以对信号的时频变换结果进行重排,使得信号能量更加集中在真实频率点附近,从而获得更高的分辨率,而且可以显著降低消噪过程中对有用信息损伤的几率。仿真数据和实际资料对比结果表明,时频重排谱分解算法比传统消噪方法具有更好的消噪效果。

关键词： 时频谱分解重排算法 ; 探地雷达 ; 短时傅里叶变换 ; 消噪 ; 瞬时频率

Abstract

In order to improve the denoising effect of the current time frequency (TF) method on ground penetrating radar signals, the authors accomplished the denoising in the TF domain by using the spectral decomposition reassignment algorithm. The application of a window function in traditional TF transforms to dividing the signal causes a blurred TF distribution, i.e., TF energies of the signal are distributed around the true instantaneous frequencies, which makes it easier to damage useful information of the signal during denosing if the threshold is inappropriate. The spectral decomposition reassignment algorithm is capable of making energy of the signal more concentrated around the true instantaneous frequencies by relocating TF result of the signal, thus obtaining higher TF resolution and reducing the opportunity of damaging useful information of the signal during denosing. Processing results of both synthetic and real data suggest that the denoising effect of the reassignment method is better than that of traditional TF including Gabor transform methods.

Keywords： spectral decomposition reassignment algorithm ; ground penetrating radar ; short time Fourier transform ; denoising ; instantaneous frequency

PDF (3155KB) 元数据多维度评价相关文章收藏文章

本文引用格式导出 EndNote Ris Bibtex

吴楠, 吴舰, 吴志坚. 探地雷达信号消噪中的时频谱分解重排算法[J]. 物探与化探, 2018, 42(1): 220-224 https://doi.org/10.11720/wtyht.2018.1.28

WU Nan, WU Jan, WU Zhi-Jian. Research on denoising of ground penetrating radar signals using the time-frequency spectral decomposition reassignment algorithm[J]. , 2018, 42(1): 220-224 https://doi.org/10.11720/wtyht.2018.1.28

0 引言

探地雷达(ground penetrating radar, GPR)已被广泛应用于公路检测、污染监测、管线追踪、考古调查、军事排雷等多个领域^[1-5]。随着对探测目标分辨率要求的不断提高,雷达使用天线的频率也在不断提高,例如在公路检测中已有使用2 GHz天线进行探测的例子^[6]。由于GPR数据在实际采集过程中不可避免地会受到各种类型的随机噪声干扰,为了降低这些噪声对GPR数据分辨率的影响,在对采集到的数据进行应用之前,必须要对原始数据进行消噪处理。目前对信号消噪的方法很多,如基于小波变换的方法、反演的方法等,其中应用较多的是频率域方法^[7]和时频域方法^[8-9];而对于非平稳信号来说,效果最好的是时频域方法。此类方法的基本思路是先将时间域信号进行时频变换,然后根据设定的阈值对某些特定时间和频率点的数据进行处理,最后再将处理后的时频域数据反变换为时间域信号。这样可以根据时间和频率两种量度对信号的含噪声区域进行定位和处理,因此时频域消噪方法会比频率域消噪方法获得更好的效果。但是,传统的时频变换方法使用的都是窗口函数,这使得信号的时频能量在时频图上会分布在信号的瞬时频率周围,在进行消噪时,如果阈值取得不合适,会导致部分有效信号被当作噪声除去。

重排是一种对时频谱进行有效提高的方法,在1976年由Kodera,De Villedary和Gendrin提出^[10],该方法通过将信号分散的时频能量进行重排,使得能量能够更加集中到信号的瞬时频率附近,从而得到分辨率更高的时频结果。目前对于重排算法的研究,大多数集中在对声音信号的处理上^[11-12],用在物理探测信号处理上很少。为了克服降噪对GRP数据分辨率产生的影响,文中用时频谱分解重排算法实现了探地雷达信号降噪处理,取得了较好的效果,最后通过仿真和实际测试进行了验证。

1 重排谱分解算法

时频变换方法是重排谱分解算法的基本前提,一个信号的时频谱可以用短时傅里叶变换(short time fourier transform, STFT)模的平方表示:

$\begin{matrix} T (t, w) = | X (t, w) |^{2} 。 (1) \end{matrix}$

其中:T(t,w)是时频谱;X(t,w)是信号s(t)的短时傅里叶变换,其表达式为

$\begin{matrix} X (t, ω) = \int s (τ) h (τ - t) e^{- jω (τ - t)} dτ, (2) \end{matrix}$

这里窗函数h(t)采用高斯窗,它可以表示为

$\begin{matrix} h (t) = 1 / (\sqrt[]{2 πσ}) e^{- \frac{t^{2}}{2 σ^{2}}}; (3) \end{matrix}$

参数变量σ决定了时间分辨率Δt=σ/ $\begin{matrix} \sqrt[]{2} \end{matrix}$ ,频率分辨率由它的倒数决定:Δw=1/( $\begin{matrix} \sqrt[]{2} \end{matrix}$ σ)。

对于短时傅里叶变换时频谱而言,起初瞬时频率和群时延被用于再现一个信号的频谱特性,但是当信号中包含多个频率分量时,一维瞬时频率只能给出这些频率的平均值,而要对这两个量进行计算还要解决解缠问题。

为了解决上述问题,提出了时频谱分解重排算法。重排的思想就是通过移动能量点,得到一个最接近频谱位置的满意点,而这个点往往在时频谱内被高斯窗口所隐藏,所以需要在Gabor变换结果的基础上进行谱分解的重新排列。这首先需要计算作为目标的时间和频率坐标,也就是说,要先算出将原时频谱上的能量集中到新时频图中的时间—频率点。从式(2)中可以看出,短时傅里叶变换可以被认为是局部信号的频率组成,如果式(2)中的窗函数取高斯窗,此时的短时傅立叶变换就为Gabor变换,其Gabor谱的极坐标形式为

$\begin{matrix} X (t, ω) = A (t, ω) \exp [jψ (t, ω)], (4) \end{matrix}$

其中:ψ(t,ω)为X(t,ω)的相位,信号目标点的重排函数时间坐标和频率坐标可以通过对ψ(t,ω)求偏导数获得,分别记为 $\begin{matrix} \dot{t} \end{matrix}$ (t,ω)和 $\begin{matrix} \dot{ω} \end{matrix}$ (t,ω):

$\begin{matrix} \begin{matrix} \dot{t} (t, ω) = t - dψ (t, ω) / dω, (5) \\ \dot{ω} (t, ω) = dψ (t, ω) / dt 。 (6) \end{matrix} \end{matrix}$

这时重排函数 $\begin{matrix} \dot{t} \end{matrix}$ (t,ω)和 $\begin{matrix} \dot{ω} \end{matrix}$ (t,ω)可以被认为是在二维时频谱域中的瞬时频率和群时延,重排的过程就是移动在原谱中的每个点(t,w)到重排坐标( $\begin{matrix} \dot{t} \end{matrix}$ , $\begin{matrix} \dot{w} \end{matrix}$ ),新时频谱中的重排目标点坐标可以通过式(5)、式(6)对原谱坐标进行移动得到。需要强调的是,这里重排的都为Gabor谱中的非零点,对于每个需要重排的点(t,w)都可以由式(5)、式(6)得到唯一对应的点( $\begin{matrix} \dot{t} \end{matrix}$ , $\begin{matrix} \dot{w} \end{matrix}$ ),而被高斯窗所掩盖的谱都有被点(t,w)位置决定的唯一中心,但是,对于时频谱内间隔很小的两个点,高斯窗口将产生重叠,虽然这会影响重排精度,但不会影响重排谱分解对数据的重构能力。

对于每个期望点的值都是来自于X(t,ω)对这个点进行重排的所有值的和,在得到新谱目标点坐标后,下一步就可以将原时频谱中的能量重新集中到目标点上,所有期望点的边际就构成了重排谱:

$\begin{matrix} X_{R} (t, ω) = \iint X (t, ν) δ [t - \dot{t} (τ, ν), ω - \dot{ω} (τ, ν)] dτdν 。 (7) \end{matrix}$

为了设计一个能够恢复信号的重构过程,这个复杂的频谱就不能只是对幅值进行重排,而必须对幅值和相位都进行重排。与能量重排不同的是,这里移动的是式(2)中X(t,w)谱的幅值而不是方程(1)中T(t,w)的幅值。可是如果对于式(5)和(6)直接进行偏导数求解会遭受相位解缠的问题。为解决这个问题,令U(t,w)=X(t,w)e^-ⁱ^wt=|F(t,w)|eⁱ⁽^φ-wt⁾,然后引入Gabor变换的一般定义:

$\begin{matrix} U (t, w) = \int s (τ) h (τ - t) e^{- iwτ} dτ 。 (8) \end{matrix}$

重排函数 $\begin{matrix} \dot{t} \end{matrix}$ (t,ω)和 $\begin{matrix} \dot{ω} \end{matrix}$ (t,ω)就可以通过计算3个Gabor变换得到:

$\begin{matrix} \begin{matrix} \dot{w} (t, w) = w - Im \frac{U^{D h} (t, w)}{U (t, w)}, (9) \\ \dot{t} (t, w) = t + Re \frac{U^{T h} (t, w)}{U (t, w)} 。 (10) \end{matrix} \end{matrix}$

式中,U^Dh(t,w)和U^Th(t,w)是通过窗函数 Dh(t)=dh(t)/dt和Th(t)=th(t)分别得到的Gabor谱,这样就避免了相位解缠问题。同时,通过简单的推导可得到Dh(t)=- $\begin{matrix} \frac{1}{σ^{2}} \end{matrix}$ Th(t)和U^Dh(t,w)=- $\begin{matrix} \frac{1}{σ^{2}} \end{matrix}$ U^Th(t,w),这时只要计算出U(t,w)和U^Th(t,w)就可以估算出重排函数。

为了展示时频重排算法对于信号时频谱分辨率的提高,设计了一个包含两部分如图1a所示的仿真信号s(t):

$\begin{matrix} \begin{matrix} s (t) = s_{1} (t) + s_{2} (t), \\ s_{1} (t) = \cos (4 πt), s_{1} (t) = \sin (2 πt 1.5); \end{matrix} \end{matrix}$

分别使用短时傅里叶变换和时频谱分解重排算法对信号s(t)进行时频变换,结果见图1b和图1c。对比二图可以看到,原先的时频变换结果中,信号的能量分布在瞬时频率周围,因此幅度谱表现为2条有宽度的能量带;而时频谱分解重排算法将信号的时频变换结果进行重新排列,将分散在信号真实的瞬时频率附近的能量集中到瞬时频率处,使得结果变为2条线,从而取得了更高的时频分辨率。由于信号的时频能量更加集中,使得有用信号和噪声的区别更加明显,而且采用时频谱分解重排算法进行消噪处理,可以尽可能地避免损伤有用信号。

显示原图| 下载原图ZIP| 生成PPT

图1 仿真信号短时傅立叶变换和时频谱分解重排方法处理结果对比

2 仿真结果

为了进行仿真测试,建立了一个地下包含一个楔形介质的模型,如图2a所示。仿真中使用的天线主频率为900 MHz,每一道数据采样18 ns,采样频率为1.7 GHz,共采集90道数据。在生成的数据中加入信噪比为10 dB的高斯白噪声以及一个全时间段都存在的、频率为700 MHz的单频率噪声,仿真信号时域剖面如图2b所示。

显示原图| 下载原图ZIP| 生成PPT

图2 仿真模型(a)及仿真信号波形(b)

对第15道数据进行短时傅里叶变换(STFT)和时频谱分解重排算法处理,结果见图3a、图3b。从图3a中可以明显看到加入的700 MHz的干扰,由于这个干扰在整个18 ns内都存在,且在5~10 ns的范围内噪声的能量和信号的能量产生了重合(图3a中箭头所示),虽然雷达发射波的主频率为900 MHz,但是由于STFT在时频分解时会将信号能量扩散到瞬时频率附近,造成了此处反射波的能量在100~2 500 MHz之间都有分布,在这种情况下进行时频消噪,就很难去除重合部分的噪声。相比之下,时频谱分解重排算法所得到的时频结果(图3b中白圈所示)就不存在这个问题。由于信号时频谱的分辨率得到了提高噪声和信号的时频能量没有混叠在一起,这使得可以在保留有用信号的基础上,将噪声尽量消除。基于短时傅里叶变换和时频谱分解重排算法两种方法的时频消噪结果分别为图3c、图3d。从图中的6~8 ns以及13~18 ns的范围可以明显看到,采用时频谱分解重排算法处理的消噪结果更好,噪声去除的更干净,且有用信号在消噪过程中受到的损伤更小。

显示原图| 下载原图ZIP| 生成PPT

图3 仿真信号消噪结果对比

3 实际资料

实测GPR数据为混凝土公路的检测资料。该公路面层为混凝土路面,下方为碎石基层。检测所使用的探地雷达的天线为900 MHz,图4是部分实际数据的剖面,在时频域对数据进行消噪处理,STFT的消噪结果如图5a,可以看到,在4~6 ns以及8~10 ns之间,干扰信号并没有能完全消除。而图5b显示,采用重排算法的消噪时,4~6 ns的区域,干扰信号基本上被完全消除了;在8 ~10 ns之间的区域,时频重排算法的结果也更好。图6为第300道数据的STFT和重排算法处理结果。

显示原图| 下载原图ZIP| 生成PPT

图4 某混凝土公路实际GPR检测剖面

显示原图| 下载原图ZIP| 生成PPT

图5 GPR数据消噪结果比较

显示原图| 下载原图ZIP| 生成PPT

图6 第20到数据2种消噪方法对比

从实际资料处理的结果来看,时频谱分解重排算法方法确实可以取得比STFT效果更好的时频消噪结果,这主要是因为时频谱分解重排算法的时频分辨率比STFT更高,使得在时频谱图上,干扰信号和有用信号可以被有效地区分开来。

4 结论

传统时频变换中窗函数的引入,导致探地雷达信号的时频能量分散在真实瞬时频率周围。进行时频消噪处理时,不合适的阈值取值会导致信号部分成份的损失。重排谱分解算法可以对信号时频变换结果进行重新排列,将信号的时频能量尽量集中到真实瞬时频率附近,从而减小能量分散程度,并得到分辨率更高的时频变换结果。能量的集中使得有效能量变强,降低了由于阈值取值不合适而引起的能量损失。文中的仿真和实际资料处理结果都表明,采用时频谱分解重排算法消噪能够取得比传统时频

消噪更好的结果。

(本文编辑:沈效群)

The authors have declared that no competing interests exist.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]	杨峰,杜翠,梁胤程,等. 基于探地雷达波速层析的金属矿区勘探仿真 [J].地球物理学进展,2014,29(3):1331-1336. [本文引用: 1]
[2]	方云,翟国林,乔梁. 探地雷达探测技术在奉先寺保护工程中的应用 [J].物探与化探,2014,38(4):815-819.
[3]	江红南. 探地雷达在干旱区盐渍化土壤层定量探测中的应用 [J].物探与化探,2014,38(4):800-803.
[4]	Shao W B, Bouzerdoum A, Phung S L, et al. Automatic classification of ground-penetrating radar signals for railway-ballast assessment [J]. IEEE Trans. Geosci. Remote Sens3, 2011, 49(10): 3961-3971.
[5]	Takahashi K, Preetz H, Igel J. Soil properties and performance of landmine detection by metal detector and ground penetrating radar-Soil characterization and its verification by a field test [J]. Journal of Applied Geophysics, 2011, 73: 368-377. [本文引用: 1]
[6]	江玉乐, 张楠. 探地雷达在隧道工程检测中的应用 [J].勘察科学技术,2008(1):58-61. [本文引用: 1]
[7]	宋瑞霞,李亚楠,张巧霞,等. 一种空域和频域相结合的图像消噪方法 [J].计算机工程与应用,2011,47(34):184-186. [本文引用: 1]
[8]	章浙涛,朱建军,匡翠林,等. 不同噪声下几种滤波方法的比较研究及改进 [J].大地测量与地球动力学,2014,34(1):127-130. [本文引用: 1]
[9]	李明峰,欧江霞,王永明,等. 基于改进小波阈值去噪法的变形预测研究 [J].大地测量与地球动力学,2014,34(3):68-71. [本文引用: 1]
[10]	Kodera K,Villedary C D,Gendrin R. A new method for the numerical analysis of non-stationary signals [J].Physics of the Earth & Planetary Interiors,1976,12(2):142-150. [本文引用: 1]
[11]	Fulop SA,Fitz K. Using the reassigned spectrogram to obtain a voiceprint [J].Journal of the Acoustical Society of America,2006,119(5):3337. [本文引用: 1]
[12]	Fulop S A. Accuracy of formant measurement for synthesized vowels using the reassigned spectrogram and comparison with linear prediction [J].Journal of the Acoustical Society of America,2010,127(4):2114-2117. [本文引用: 1]

基于探地雷达波速层析的金属矿区勘探仿真

2014

... 探地雷达(ground penetrating radar, GPR)已被广泛应用于公路检测、污染监测、管线追踪、考古调查、军事排雷等多个领域^[1-5].随着对探测目标分辨率要求的不断提高,雷达使用天线的频率也在不断提高,例如在公路检测中已有使用2 GHz天线进行探测的例子^[6].由于GPR数据在实际采集过程中不可避免地会受到各种类型的随机噪声干扰,为了降低这些噪声对GPR数据分辨率的影响,在对采集到的数据进行应用之前,必须要对原始数据进行消噪处理.目前对信号消噪的方法很多,如基于小波变换的方法、反演的方法等,其中应用较多的是频率域方法^[7]和时频域方法^[8-9];而对于非平稳信号来说,效果最好的是时频域方法.此类方法的基本思路是先将时间域信号进行时频变换,然后根据设定的阈值对某些特定时间和频率点的数据进行处理,最后再将处理后的时频域数据反变换为时间域信号.这样可以根据时间和频率两种量度对信号的含噪声区域进行定位和处理,因此时频域消噪方法会比频率域消噪方法获得更好的效果.但是,传统的时频变换方法使用的都是窗口函数,这使得信号的时频能量在时频图上会分布在信号的瞬时频率周围,在进行消噪时,如果阈值取得不合适,会导致部分有效信号被当作噪声除去. ...

探地雷达探测技术在奉先寺保护工程中的应用

2014

探地雷达在干旱区盐渍化土壤层定量探测中的应用

2014

Automatic classification of ground-penetrating radar signals for railway-ballast assessment

2011

Soil properties and performance of landmine detection by metal detector and ground penetrating radar-Soil characterization and its verification by a field test

2011

探地雷达在隧道工程检测中的应用

2008

一种空域和频域相结合的图像消噪方法

2011

不同噪声下几种滤波方法的比较研究及改进

2014

基于改进小波阈值去噪法的变形预测研究

2014

A new method for the numerical analysis of non-stationary signals

1976

... 重排是一种对时频谱进行有效提高的方法,在1976年由Kodera,De Villedary和Gendrin提出^[10],该方法通过将信号分散的时频能量进行重排,使得能量能够更加集中到信号的瞬时频率附近,从而得到分辨率更高的时频结果.目前对于重排算法的研究,大多数集中在对声音信号的处理上^[11-12],用在物理探测信号处理上很少.为了克服降噪对GRP数据分辨率产生的影响,文中用时频谱分解重排算法实现了探地雷达信号降噪处理,取得了较好的效果,最后通过仿真和实际测试进行了验证. ...

Using the reassigned spectrogram to obtain a voiceprint

2006

Accuracy of formant measurement for synthesized vowels using the reassigned spectrogram and comparison with linear prediction

2010

〈

〉

探地雷达信号消噪中的时频谱分解重排算法

Research on denoising of ground penetrating radar signals using the time-frequency spectral decomposition reassignment algorithm

0 引言

1 重排谱分解算法

2 仿真结果

3 实际资料

4 结论

参考文献 文献选项 原文顺序 文献年度倒序 文中引用次数倒序 被引期刊影响因子

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子