基于EMD分解的探地雷达信号瞬时频率分析

1、第37卷 第4期 2009年8月煤田地质与勘探COAL GEOLOGY & EXPLORATIONVol. 37 No.4 Aug. 2009文章编号: 1001-1986(2009)04-0064-04基于EMD分解的探地雷达信号瞬时频率分析杨秋芬(西安文理学院物理系，陕西 西安 710065)摘要: 介绍了Hilbert-Huang变换中经验模式分解(EMD)的基本原理；讨论了实际探地雷达信号处理中EMD分解的终止条件；给出了利用内蕴模式函数(IMF)计算信号瞬时频率的计算公式。实际探地雷达剖面的HHT(Hilbert-Huang Transform)分析表明，由IMF得到的瞬时频率剖面对

2、埋地目标IMF 关 键 词：HHT；内蕴模式函数；瞬时频率中图分类号：P631 文献标识码：A DOI: 10.3969/j.issn.1001-1986.2009.04.017The instantaneous frequency analysis of GPR data usingempirical mode decompositionYANG Qiufen(Departement of physic, Xian University of Arts and Science, Xian 710065, China)Abstract: The foundational theory abou

3、t HHT(Hilbert-Huang Transform)based on EMD (empirical mode decom-position)is introduced in brief. The stop condition of EMD on the rear GPR data is discussed. The formula of cal-culating instantaneous frequency using IMF (intrinsic mode function)is presented. The identification of real buried object

4、 results shows that the IMF method is superior to the ordinary method. The multi-resolution of IMF is ana-lyzed.Key words: HHT; intrinsic mode function; instantaneous frequencyN. E. Huang.等1提出的Hilbert-Huang变换(Hilbert-Huang Transform, 简称HHT)是一种适用于非平稳、非线性信号分析的自适应时频分析方法。HHT首先把信号通过经验模式分解(Empirical Mode

5、Decomposition, 简称EMD)成若干个具有窄带特性的内蕴模式函数(Instrinsic Mode Function，简称IMF)。IMF具有时间上的局域化特征，属于窄带信号，满足2个条件：其一，对整个时间序列来说，极值点数与过零点的个数相同或其差值为1；其二，在任意一点处，极大值包络与极小值包络的均值为零。然后利用Hilbert变换构造解析信号，得到瞬时频率和振幅。在算法研究方面，Huang等2提出了归一化HHT，对信号进行归一化处理；P. Flandrin等3通过对高斯白噪声的分析，发现EMD方法表现为时域的二进滤波；Wu等4也发现了类似的规律，并应用于基于噪声统计规律的去噪；S

6、. Kizhner等5对EMD首先把快变分量提取出来的原因进行了理论证明；C. Dameral等6则提出了二维EMD快速算法；S. R. Qin等7为消除3次样条求EMD包络引起的“过冲”，提出了用分段幂级数函数求包络的方法；Q. Chen等8则采用三次B样条代替3次样条作为插值函数，提出基于B样条的EMD，具有更好的特性；G. Rilling等9通过对2个正弦信号叠加而成的合成信号的分析，发现了2个分量的频率、幅度的相对值与频率分辨率之间的关系，并给出了理论证明；Tanaka等10把EMD方法推广到了复数域；另外，邓拥军11、黄大吉12、张郁山13等还提出了镜像延拓、神经网络延拓、边筛选边延

7、拓等多种边界处理方法。HHT方法已在各领域得到了应用，如故障诊断、医学信号分析、地震信号分析等1417。瞬时参数是探地雷达信号解释中的有力工具之一，对于地层的构造形态及地下介质的物性认识都非常重要18。传统的瞬时参数大部分都基于Hilbert变换，利用Hilbert变换计算信号瞬时参数，要求信号必须是窄带的。探地雷达信号属于宽带信号，直接采用Hilbert计算其瞬时参数可能带来错误的物理解释，如负频率等。另外，通过EMD分解，得到的内蕴模式信号之间也具有类似小波分析等时频分析方法的多分辨特性。本文给出了探地雷达信号EMD分解算法，给收稿日期: 2009-01-19作者简介:杨秋芬(1971),

8、 女, 陕西乾县人, 助教, 硕士, 从事探地雷达信号处理研究工作.第4期杨秋芬: 基于EMD分解的探地雷达信号瞬时频率分析 65 出了利用内蕴模式函数计算瞬时频率的方法。通过对实际探的雷达资料的分析，比较了基于HHT的瞬时频率与直接采用Hilbert变换计算瞬时频率的 剖面能量时，EMD分解就可以停止，则确定了其分解级数。 1.3 基于EMD分解探地雷达信号的瞬时频率 差异。 分析通过EMD分解，得到具有窄带特性的内蕴模1 探地雷达信号HHT分析原理 式函数剖面imfk(x, t)；沿空间方向逐点对内蕴模式函数imfk(x, t)进行一维Hilbert变换，得到其对应的1.1 EMD算法简介

9、1虚部剖面Himfk(x, t)，其中H表示一维Hilbert对任意一个一维信号s(t)，EMD方法可表示为：变换因子。这样，得到了与原始探地雷达剖面对应 s(t)=imfk(t)+r(t)， (1)k得解析信号剖面Dk(x, t)，则其中 imfk(t)代表信号的第k个IMF分量；r(t)为残留Dk(x,t)=imfk(x,t)+iHimfk(x,t)， (3) 分量。具体算法大量文献都有介绍，这里不再赘述。 式中 i表示虚数单位。另外，忽略了残量剖面。由1.2 探地雷达信号的EMD分解 式(3)可以分别得到解析信号剖面Dk(x, t)的瞬时幅度若对探地雷达测量的二维测线进行分析，设接ak(

10、x, t)和瞬时相位k(x, t)及瞬时频率剖面k(x, t)： 收到信号表示为d(x, t)，其中x表示测点、t表示时 ak(x,t)=(imfk(x,t)2+(Himfk(x,t)21/2 ； (4) 间。忽略随机噪声的影响，信号中包含直达波及目Himfk(x,t)(x,t)arctan( ； =(5) k标散射信号，则任意测点x处探地雷达信号的EMDimfk(x,t)分解可表示为： Himfk(x,t)dk(x,t)d=(x,t)arctan(2) k d(x,t)=imfk(x,t)+r(x,t)， 。(6) dtdtimf(x,t)kk式中 k表示第k级IMF分量。对二维测线沿空间2

11、 实际资料算例 方向逐点进行EMD分解，可以得到整条测线对应的各级IMF分量剖面。在对实际探地雷达信号进行图1a为某沙坑中实测探地雷达资料，横向70道信号，采样间隔1 cm，纵向512个采样点，时窗EMD分解时，当残留分量剖面的能量远远小于原始图1 探地雷达剖面的EMD分解a 原始探地雷达剖面；b 第一级IMF分量剖面；c 第二级IMF分量剖面；d 第三级IMF分量剖面 Fig. 1 EMD of GPR section 66 煤田地质与勘探 第37卷20 ns。埋地目标为长52.5 cm、宽40.5 cm的铜板，埋深为22.5 cm，雷达信号中心频率1 GHz。测线方向垂直于宽度方向。图1a

12、中能清晰观察到直达波信号及目标散射信号。由于直达波干扰，散射信号的起始位置难以判定。图1b、图1c和图1d分别表示图即埋地目标位置。当2列波在时间上进行叠加，在叠加位置处，接收信号产生较大奇性，表示为高频分量。原始剖面信号由于是宽带信号，信号突变位置能量较弱，其瞬时频率指示不明显。经过EMD分解，IMF信号具有窄带特性，突出了信号突变位置，则对应的瞬时频率可清晰表示目标信号，如图3、图4中箭头1a经EMD分解得到的前三级IMF分量剖面，可以看出，与原始探地雷达信号相比，IMF分量剖面由于具有窄带特性，其信号能量幅度变化比较平缓。另外，第一级频率最高，随IMF分级的增加，其对应的频率逐渐降低，也

13、就是其具有多分辨率特性。图2a是由图1a得到的瞬时频率剖面，图2b、图2c和图2d分别对应图1b、图1c和图1d的瞬时频率剖面。与直接由宽带信号得到的瞬时频率剖面相比，由具有窄带特性的IMF剖面得到的瞬时频率剖面可以清晰的解释埋地目标位置。图2中箭头A所示为发射及接收天线直接耦合波及通过地表传播的电磁波叠加位置；箭头B表示散射信号到达位置，A、B所示。图3为图1a、图1b、图1c和图1d中对图4为图3对应的瞬时频率信号。应的第35道信号；另外，由各级IMF剖面得到的瞬时频率剖面的多分辨率特性清晰，可以联合判断目标位置及空间展布。3 结 论通过对实测探地雷达信号的HHT分析可以看出，由内蕴模式函

14、数计算得到的瞬时频率能清晰地识别埋地目标，并减小直达波对浅层埋地目标的影响。通过识别层状目标，表明HHT分析对于公路质量检测中地层位置及层厚的估计有较好应用前景。图2 探地雷达信号内蕴模式函数瞬时频率分析Fig. 2 Instantanecus frequency of intrinsic mode function of GRP singalsa 原始探地雷达信号的瞬时频率剖面；b 第一级IMF分量的瞬时频率剖面；c 第二级IMF分量的瞬时频率剖面；d 第三级IMF分量的瞬时频率剖面第4期杨秋芬: 基于EMD分解的探地雷达信号瞬时频率分析 67 参考文献1 HUANG N E, SHEN Z

15、, LONG RS, et al. The empirical mode de-composition and the Hilbert spectrum for nonlinear and non- station-ary time series analysisJ. Proc.R.Soc.Lond. A, 1998, 454：903995. 2 HUANG N E. Normalized Hilbert transform and instantaneous fre-quencyJ. NASA Patent Pending, GSC, 2003, 14：6731.3 FLANDRIN P,

16、RILLING G, GONCALVES P. Empirical mode de-composition as a filter bankJ. Signal Processing Letters, IEEE, 2004, 11 (2)：112114.4 WU Z, HUANG N E. A study of the characteristics of white noiseusing the empirical mode decomposition methodJ. ProcRSocLon-don A, 2004, 460：15971611.5 KIZHNER S, BLANK K, FL

17、ATLY T, et al. On certain theoreticaldevelopments underlying the Hilbert Huang transformC. Paris：Aerospace Conference, IEEE, 2006.6 DAMERVAL C, MEIGNEN S, PERRIER V. A fast algorithm for图3 单道探地雷达信号的EMD分解Fig. 3 EMD of GPR single trace signalsa 图1a第35道信号；b 对应图3a的第一级IMF分量剖面；c 对应图3a的第二级IMF 分量剖面；d对应图3a的第

18、三级IMF分量剖面bidimensional EMDJ. Signal Processing Letters, IEEE, 2005, 12 (10): 701-704.7 QIN S R, ZHONG Y M. A new envelope algorithm of Hilbert-HuangTransformJ. Mechanical System and Signal Processing, 2006, 20: 1941-19528 QIUHUI C, HUANG N E, RIEMENSCHNEIDER S, et al. AB-spline approach for empirica

19、l mode decompositionsJ. Advances in Computational Mathematics, 2006, 24：171195.9 RILLING G, FLANDRIN P. One or two frequencies? The empiricalmode decomposition answersJ. Transactions on Signal Processing, IEEE, 2007, 56(1): 85-9510 TANAKA T, MANDIC D P. Complex empirical mode decomposi-tionJ. Signal Processing Letters, IEEE, 2007, 14 (2): 101-104. 11 邓拥军, 王伟, 钱成春, 等. EMD方法及Hilbert变换中边界问题的处理J. 科学通报, 2001, 46(03): 257-263.12 黄大吉, 赵进平, 苏纪兰. 希尔伯特-黄变换的端点延拓J. 海洋学报(中文版), 2003, 25(01): 1-11.13 张郁山, 梁建文, 胡聿贤. 应用自回归模型处理EMD方法中的边界问题J. 自然科学进展, 2003, 13(10): 105