短时傅里叶变换在GPR数据解释中的应用

引用本文

戴前伟, 吴铠均, 张彬. 短时傅里叶变换在GPR数据解释中的应用[J]. 物探与化探, 2016,40(6): 1227-1231.
DAI Qian-Wei, WU Kai-Jun, ZHANG Bin. A study of application of short-time Fourier transform to GPR data interpretation[J]. Geophysical and Geochemical Exploration, 2016,40(6): 1227-1231.
Doi:10.11720/wtyht.2016.6.27 复制到剪切板

Permissions

短时傅里叶变换在GPR数据解释中的应用

戴前伟^1,², 吴铠均¹, 张彬^1,²

1.中南大学地球科学与信息物理学院,湖南长沙 410083

2.中南大学有色金属成矿预测与地质环境监测教育部重点实验室,湖南长沙 410083

通讯作者:吴铠均(1993-),女,中南大学地球科学与信息物理学院研究生在读,从事地球物理勘探工作。E-mail:540416481@qq.com

作者简介: 戴前伟(1968-),男,工学博士、教授,博士生导师,学科专业:地球探测与信息技术。长期从事电磁法勘探理论及应用和工程、环境及灾害地球物理勘探的基础理论及实践的教学和科研工作。

收稿日期: 2016-03-14

基金: 国家自然科学基金项目(41374118)

摘要

小波变换、经验模态分解、奇异值分解是提高探地雷达数据信噪比、突出深部异常的有效手段。笔者提出一种快速的短时傅里叶变换技术,对探地雷达噪声信号进行解译,利用雷达发射信号与噪声信号的频谱差异提高雷达数据的信噪比;并以数值算例进行验证,分别从单道波信号的频率能量差异,主频等值分布差异等方面详细分析了该方法与小波变换、经验模态分解方法的优势。结果表明,短时傅里叶分析技术能实现探地雷达数据的快速解译,为探地雷数据的反演成像提供更精确的约束条件。

关键词: 探地雷达; 短时傅里叶; 频率差异; 主频等值分布; 快速解译

中图分类号:P631 文献标志码:A 文章编号:1000-8918(2016)06-1227-05 doi: 10.11720/wtyht.2016.6.27

A study of application of short-time Fourier transform to GPR data interpretation

DAI Qian-Wei^1,², WU Kai-Jun¹, ZHANG Bin^1,²

1.School of Geosciences and Info-Physics,Central South University,Changsha 410083,China

2.Key Laboratory of Metallogenic Prediction of Nonferrous Metal and Geological Environment Monitoring,Ministry of Education,Central South University,Changsha 410083,China

Abstract

Wavelet transform,empirical mode decomposition and singular value decomposition are effective means to highlight the deep anomaly and improve the signal-to-noise ratio of ground penetrating radar (GPR) data.This paper proposes a fast short-time Fourier transform technique for interpreting ground penetrating radar (GPR) noise signal.The frequency spectrum differences between transmit signal and noise signal are used to improve signal-to-noise ratio of radar data.It is verified by using a numerical example.The method and the wavelet transform as well as the advantage of empirical mode decomposition method are analyzed in detail respectively from single signal frequency energy difference,frequency equivalent distribution difference,and other aspects.The results show that the short time Fourier analysis technology can achieve rapid interpretation of GPR data and,for the ground penetrating radar data inversion imaging,it can provide more precise constraints.

Keyword: ground penetrating radar; short time Fourier transform; difference of frequency; frequency equivalent distribution; fast interpretation

Show Figures

传统信号处理中, 分析和处理信号最常用与最直接的方法是傅里叶变换。一般来说, 以时间为自变量来表示信号, 通过傅里叶变换, 信号也能用频率为自变量来表示, 即频谱。傅里叶变换是相关函数和功率谱之间联系的桥梁。然而, 傅里叶变换作为一种整体变换, 信号不能同时在时域和频率表征出来, 不能看出某种频率分量出现的时间及其变化情况。在实际应用中, 就探地雷达信号而言, 对其进行单一的时域、频域分析已不能满足实际需要, 并且实际探测中, 信号受高频杂波干扰程度较大, 有时甚至淹没目标信号, 因而需要了解信号频谱及其能量大小随时间的变化情况。

时频分析的主要任务是描述信号频谱含量随时间的变化。基本思想是要描述信号在不同时间和频率内的能量密度和强度, 即要设计出时间和频率的联合函数。该联合函数即为时频分布, 利用它来分析信号, 可以在每一时刻看出信号在瞬时频率附近的能量聚集情况, 物理意义明确, 对整个信号采用单一分辨率进行研究, 可以反映信号的整体时频趋势。目前常用的时频分析方法有短时傅里叶、小波变换等。短时傅里叶变换的优点在于概念直接、简单, 实现容易, 成为分析非平稳信号的有力工具。

1 基于短时傅里叶变换的时频分析技术

1.1 短时傅里叶变换的基本原理

短时傅里叶变换(STFT)可以视为加窗的傅里叶变换, 表达式为:

其中:x(m)为反射信号, n为时间变量。w(n)为窗函数, 作用是截取出x(n)在n时刻附近的一段信号, 并对其进行傅里叶变换, 当n发生变化时, 窗函数也随之移动, 由此得到信号的频谱随时间n的变化规律。此时, 傅里叶变换成为一个二维域的(n, ω )函数(图1)。

	Figure Option View Download New Window
	图1 窗函数的移动

1.2 窗函数的选择

由式(1)可知, 短时傅里叶变换使用固定的窗函数w, 一旦确定, 其形状就不能发生改变, 因此也确定了短时傅里叶变换的分辨率。

为了改善频谱泄露情况, 笔者采用海明窗函数。窗函数的长度也同时影响着时间分辨率与频率分辨率, 窗函数越短, 时间分辨率越高, 而频率分辨率越低; 窗函数越长则频率分辨率越高, 而时间分辨率降低。因而, 两者不能同时任意减小, 只能相对折衷, 为了使得窗函数时频聚集性良好, 其宽度应与信号的局域平稳长度相适应, 笔者取窗长为雷达子波长度的整数倍, 点数为121。

1.3 时频分析的方法

时频分析方法可以把GPR高频电磁波的频谱与其到达的时间联系起来, 分析不同频率的电磁波沿水平或垂直方向的变化规律。首先选取合适的时窗对一道数据进行短时傅里叶变换, 时窗的长度要同时兼顾短时傅里叶变换的时间分辨率与频率分辨率, 每一个时窗的时频分析结果中, 其傅里叶变换振幅谱的极值频率可以看作该时间段内的主频率值。随着时窗的移动, 对一道数据中不同延迟时间进行逐段分析, 可以得到这一道的不同延迟时间的主频信息, 再把多道数据的时频分析结果横向连接起来, 根据探地雷达主频范围及变化规律, 就可以区分出高频杂波和雷达波, 从而辨别目标体埋藏情况。

时频分析主频的计算方法是对采集的每一道数据进行短时傅里叶变换, 提取主频, 按一定时间步长计算时窗中心时间的主频。根据各采集道的计算结果, 以测线长度为横坐标(x轴), 时窗中心时间为纵坐标(y轴)绘制时频域等值色图, 在图上可以直观地表现出异常频谱沿测线方向和沿时间的分布, 而根据地震记录的时间特征和频率特征推断出异常体的位置。

2 探地雷达数字模拟算例

2.1 模拟信号获取及加噪

设计探地雷达应用几何模型(图2), 模型长 2 m, 高0.5 m, 坐标(0.88, 0.25)处有一圆形空洞, 半径为2.5 cm。利用GPRMAX V2.0进行探地雷达探测正演数值模拟, 天线发射频率为900 MHz, 获取100道数据, 绘制雷达波反射剖面(图3)。再用GPRMAX模拟出一个相同情况下不含空洞的数据模型, 用含目标体和直达波的数据与只含直达波的数据相减, 由此去除直达波, 做出去除直达波后该组数据的剖面(图4)。

	Figure Option View Download New Window
	图2 介质模型

	Figure Option View Download New Window
	图3 雷达反射剖面

	Figure Option View Download New Window
	图4 去除直达波后雷达反射剖面

采用MATLAB产生高斯白噪声信号, 并将噪声信号幅度扩大46倍后与去除直达波的信号叠加, 从而对信号进行加噪。从加噪后的信号剖面(图5)可以看出, 有效信号已经完全淹没在噪声信号中, 无法判断地下空洞埋藏情况。

	Figure Option View Download New Window
	图5 加噪后不含直达波的雷达反射剖面

2.2 信号短时傅里叶时频分析

首先选取含噪雷达信号(图5)中的第一道数据做短时傅里叶变换, 从雷达反射剖面(图3)中可以推断, 第一道数据未能检测到异常, 因而第一道的时频图(图6)中, 只有能量不高的杂波信号频率, 杂乱无序地分布在时频图各区域。而在第44道数据时频图(图7)中, 虽然也有杂乱分布的杂波信号频率, 但仍能分辨出能量较高的目标体频率位置, 并且能与雷达波主频范围及其发生的时间相对应, 与不加噪声的第44道信号时频图(图8)相一致。

	Figure Option View Download New Window
	图6 第1道含噪雷达数据时频

	Figure Option View Download New Window
	图7 第44道含噪雷达数据时频

	Figure Option View Download New Window
	图8 第44道未含加噪雷达数据时频

对100道数据分别进行短时傅里叶分析, 并把每一道中心时间的主频横向连接, 作出整个测线剖面的主频等值图(图9)。由于高斯噪声信号随机性较大, 在短时傅里叶的时窗还未移动到目标信号位置时, 中心时间的主频不一致, 在主频等值图中十分凌乱, 时窗移动到目标信号位置后, 由于目标体反射的电磁波能量较规律, 目标信号范围内可以得到相对较一致的主频值, 并在横向范围得到体现。大约在0.9 m处, 5 ns时检测到雷达发射频率范围的信号, 由此推断, 该位置存在地质异常。

	Figure Option View Download New Window
	图9 主频等值图

3 探地雷达实测信号的分析

图10为探地雷达实验记录实测数据剖面, 发射频率为900 MHz, 401道, 探测目标为一圆形空洞, 横坐标为雷达水平移动距离, 纵坐标为时间。从第150道数据时频来看, 频率能量大的区域集中在8 ns左右, 900 MHz附近, 该道数据位置对应剖面图横坐标(图10)大约0.74 m处, 检测到异常的时间在剖面(图11)中也显示为8 ns附近; 同理, 第200道数据对应剖面(图10)横坐标1 m处, 其时频(图12)能量较大区域的时间6 ns与检测到异常的时间也相符。同时, 从主频等值图(图13)中, 也能从频率差异分辨出异常的存在。

	Figure Option View Download New Window
	图10 探地雷达实测数据剖面

	Figure Option View Download New Window
	图11 实测数据第150道时频

	Figure Option View Download New Window
	图12 实测数据第200道时频

	Figure Option View Download New Window
	图13 探地雷达实测数据主频等值图

4 结论

1)用短时傅里叶对雷达信号进行分析, 结合了时间域与频率域, 弥补了传统傅里叶方法一维方向的不足, 更直观地体现频率随时间的变化情况。

2)对某一道数据进行短时傅里叶分析, 可以从其时频图中看出该道反射波不同频率实幅差异, 再结合雷达发射频率, 推断目标体存在情况。

3)对含噪雷达信号的分析, 从横向连接的主频等值图中, 可以把噪声信号与目标信号的频率区分开, 省去滤波这一步骤, 把被淹没的目标信号从复杂的噪声信号中识别出, 在实际应用中易适用于初步地质情况的推断。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	李华, 鲁光银, 何现启, 等. 探地雷达的发展历程及其前景探讨[J]. 地球物理学展, 2010, 25(4): 1492-1502. [本文引用:1]
[2]	刘喜武, 张宁, 勾永峰, 等. 地震勘探信号时频分析方法对比与应用分析[J]. 地球物理学进展, 2008, 23(3): 743-753. [本文引用:1]
[3]	徐浩, 刘江平, 范承余, 等. 隧道衬砌病害的探地雷达波场模拟与特征分析[J]. 中南大学学报: 自然科学版, 2013, 44(11): 4581-4587. [本文引用:1]
[4]	吴宝杰, 杨桦, 张伟光. 探地雷达数据的S变换时频分析[J]. 上海地质, 2008, 1(3): 13-15+19. [本文引用:1]
[5]	戴前伟, 王洪华, 冯德山, 等. 基于双二次插值的探地雷达有限元数值模拟[J]. 地球物理学进展, 2012, 27(2): 736-743. [本文引用:1]
[6]	刘秀娟, 邓世坤, 徐保林. 探地雷达信号小波变换去噪[J]. 工程勘察, 2006, 1(10): 66-71. [本文引用:1]
[7]	李大心. 探地雷达方法及应用[M]. 北京: 地质出版社, 2000. [本文引用:1]
[8]	肖兵, 鲍光淑, 赵秋梅, 等. 探地雷达复信号分析及改进[J]. 中南工业大学学报, 1997, 28(1): 4-6. [本文引用:1]
[9]	王超, 沈斐敏. 一维HHT变换在探地雷达数据处理中的应用[J]. 工程地质学报, 2015, 23(2): 328-334. [本文引用:1]
[10]	单娜琳, 程志平, 丁彦礼. 地震映像数据的时频分析方法及应用[J]. 地球物理学进展, 2007, 20(6): 1740-1745. [本文引用:1]
[11]	庞存锁, 刘磊, 单涛. 基于短时分数阶傅里叶变换的时频分析方法[J]. 电子学报, 2014, 42(2): 347-352. [本文引用:1]
[12]	冯德山. 基于小波多分辨探地雷达正演及偏移处理研究[D]. 长沙: 中南大学, 2006. [本文引用:1]
[13]	冯德山, 戴前伟. 探地雷达小波域三维波动方程偏移[J]. 地球物理学报, 2008, 51(2): 566-574. [本文引用:1]
[14]	郭立, 崔喜红, 陈晋. 基于GprMax正演模拟的探地雷达根系探测敏感因素分析[J]. 地球物理学进展, 2012, 27(4): 1754-1763. [本文引用:1]
[15]	刘康. 基于MATLAB的探地雷达数据处理研究及软件开发[D]. 北京: 中国地质大学(北京), 2011. [本文引用:1]
[16]	闫永斌. 探地雷达图像识别与处理[D]. 北京: 北京化工大学, 2008. [本文引用:1]
[17]	张虎. 探地雷达数值模拟与成像方法研究[D]. 武汉: 长江大学, 2013. [本文引用:1]
[18]	Xing H, Li J, Chen X, et al. GPR reflection position identification by STFT[C]//Tenth International Conference on Ground Penetrating Radar, 2004. [本文引用:1]
[19]	Lopera O, Milisavljevi N, Daniels D, et al. Time frequency domain signature analysis of gpr data for land mine identification[C]//The 4^th International Workshop on Advanced Ground Penetrating Radar, 2007. [本文引用:1]
[20]	Zhou H L, Mao T, Chen X L. Time-frequency representations for classification of ground penetrating radar echo signal[C]// Proceedings of 2005 International Symposium on Intelligent Signal Processing and Communication Systems, 2005. [本文引用:1]