Hilbert-Huang变换与大地电磁信号的时频分析

汤井田¹， 蔡剑华^{1, 2}，任政勇¹，化希瑞¹

(1. 中南大学 信息物理工程学院，湖南 长沙，410083；

2. 湖南文理学院 物理与电子科学系，湖南 常德，415000)

摘  要：将Hilbert-Huang变换引入大地电磁信号的时频分析中，介绍HHT(Hilbert-Huang transform)时频分析原理及方法，给出仿真信号的经验模态分解及其时频分布，并对实测大地电磁信号进行HHT时频处理与剖析。研究结果表明：Hilbert能量谱随时频的具体分布具有很强的非稳态动态变换时频刻画能力；时频谱的时间、频率分辨率不受Heisenberg测不准原理的限制，且其时间、频率分辨率都很高，有很好的时频聚集性；HHT方法能用于描述大地电磁信号的非线性时变特征，是大地电磁信号时频分析的有效工具。

关键词：Hilbert-Huang变换；经验模态分解；大地电磁信号；时频分析

中图分类号：P631          文献标识码：A         文章编号：1672-7207(2009)05-1399-07

Hilbert-Huang transform and time-frequency analysis of magnetotelluric signal

TANG Jing-tian¹, CAI Jian-hua^{1, 2}, REN Zheng-yong¹, HUA Xi-rui¹

(1. School of Info-physics and Geomatics Enginerring, Central South University, Changsha 410083, China;

2. Department of Physics and Electrics, Hunan University of Arts and Science, Changde 415000, China)

Abstract: Hilbert-Huang transform was proposed as a method of time-frequency analysis in magnetotelluric (MT) signals processing. Firstly, the time-frequency analysis method of Hilbert-Huang transform and conception of transient frequency were introduced, and the empirical mode decomposition and time-frequency distribution of a simulative signal was carried out. Then the actual MT data series were analyzed and processed based on Hilbert-Huang transform (HHT) spectrum. The results show that Hilbert energy spectrum can clearly show the energy distribution with time and frequency in detail and has strong ability to analyze non-stationary signal. In HHT spectrum, time and frequency resolving power is strong and not restricted by Heisenberg principle. Hilbert-Huang transform can be used to precisely describe the nonlinear and non-stationary characteristics of MT signals, and it is an effective tool of time-frequency analysis in MT signals processing.

Key words: Hilbert-Huang transformation; empirical mode decomposition; magnetotelluric signal; time-frequency analysis

大地电磁(MT)信号具有非线性、非平稳、非最小相位特征，是一种典型的非平稳随机信号^[1-2]。目前，对大地电磁信号进行分析的主要工具有傅里叶变换、现代谱估计和小波变换等。傅里叶变换和现代谱估计建立在稳态信号处理基础之上，它仅能给出信号总体所包含的各种频率成分,不能解决何时出现何频率的问题^[3-4]。小波理论与以往其他研究方法相比，其频谱分辨率和时间定位精度都得到显著提高，但小波变换的有效性依赖于小波函数的选取，有时还会存在随着尺度增大，相应正交基函数的频谱局部性变差的缺陷，使其对信号的更精细分解受到限制。1998年，Huang等^[5-7]提出的Hilbert-Huang变换(Hillbert-Huang transformation, HHT)是最新发展起来的处理非线性、非平稳信号的时频分析方法，其基本的实现步骤分为2步：多分辨经验模态分解EMD (Empirical mode decomposition)和瞬时频率的求解。随后，可以获得信号的时频谱，把信号的幅度(能量)表示成时间和频率的函数。HHT变换是一种对以Fourier变换为基础的改进的信号处理方法，其在电法勘探尤其是在大地电磁信号的处理中应用较少。在此，本文作者把Hilbert- Huang变换引入大地电磁信号的时频分析中，从多方面探讨HHT变换在大地电磁信号分析中的适用性和处理非平稳MT信号的能力。

1  HHT变换原理

1.1  Huang变换

EMD方法即Huang变换，它将信号分解为含有不同时间尺度且满足以下2个定义条件的1组IMF (Intrinsic mode function)：

a. 对于一列数据，极值点和过零点数目必须相等或至多相差1点。

b. 在任意点，由局部极大点和极小点构成的2条包络线的平均值为0。每个IMF可以被认为是信号中固有的1个模态函数。

首先，找到信号x(t)的极大值和极小值，通过3次样条拟合，获得信号的上包络线和下包络线，计算上、下包络线在每点的平均值，从而获得1条平均值曲线m₁，用每点的信号x(t)减去对应点的m₁可得到新数据序列h₁：

h₁一般不满足IMF的条件，为此，重复以上过程n次，使所得的平均包络线每点的值趋于0，此时，h_n就是第1阶IMF(C₁)，从而也得到第1阶剩余项R₁。

对R₁重复上述过程，就可以获得第2阶IMF分量。通过EMD对信号一次次筛选，就可以获得多个IMF分量和1个逼近分量R_n，从而信号可表示为：

该分解过程可以理解为多尺度滤波的过程，每个IMF分量都反映了信号的特征尺度，代表非线形信号的内在模态特征。

1.2  Hilbert变换与时频谱

信号经分解后得到多个IMF 的组合。对每个IMF 分量进行Hilbert变换，即可得到每个IMF分量的瞬时频率，综合所有IMF分量的瞬时频谱就可获得Hilbert谱^[5-8]。

获得IMF分量后，就可以对每一阶IMF(设为c(t))进行Hilbert变换。

解析信号的极坐标形式反映了Hilbert变换的物理含义：它是通过一正弦曲线的频率和幅值调制获得信号局部最佳逼近。根据瞬时频率的定义：

对每一阶IMF作Hilbert变换，并求出相应解析函数的幅值谱和瞬时频率，从而原始信号可以表示为：

其数学表达式反映出Hilbert-Huang变换是Fourier变换的一种扩展形式。

式(9)反映了信号的幅值、时间和瞬时频率之间的关系。信号的幅值能表示为时间、瞬时频率的函数H(ω,t)，从而获得信号幅值的时间和频率分布即Hilbert谱，通过对时间的积分可以获得信号的Hilbert边际谱。

边际谱表达了每个频率在全局上的幅度(或能量)，它代表了在统计意义上的全部累加幅度。

2  仿真信号的经验模态分解与时频分析

2.1  计算机仿真信号

仿真信号的EMD分解、固有模态函数IMF的重构和Hilbert变换的时频分析，是检验HHT方法实际应用效果的一种最简单、直观和有效的方法^[9]。参考文献[9]，设计了含有2种成分的已知电压信号，其解析表达式为：

信号由一基频50 Hz、调制频率15 Hz的调频调幅成分以及150 Hz的正弦波叠加而成。按EMD原理，上述信号的一阶IMF应为150 Hz的信号，2阶IMF应为50 Hz的调频调幅信号，更高阶的IMF将是EMD分解的残差，其值越小，表明EMD分解效果越好。

2.2  经验模态分解

进行Hilbert变换的前提是EMD分解，即Huang变换，EMD分解质量是时频谱估计效果的关键^[10-11]。图1(a)所示为仿真信号和6阶IMF时域波形曲线，由7条曲线组成。其中：signal曲线代表式(11)的仿真信号；imf1代表1阶序列(C₁)，是最先被分解出来的高频信号，对应150 Hz的正弦波；imf2代表2阶序列(C₂)，是第2次被分解出来的信号，对应50 Hz的调频调幅信号；imf3，imf4，imf5和residual分别代表3，4和5阶序列(C₃，C₄和C₅)和剩余量R₆。可以看出，它们的变化幅度已是原始幅值的1/200，接近1条直线，说明EMD分解结果与预计结果十分吻合。图1(b)所示为对应信号和6阶IMF频谱，可以看出信号50 Hz的基频及150 Hz的正弦波；imf1和imf2的频谱也恰好只有150 Hz和50 Hz的成分，imf3，imf4，imf5和residual频谱中的振幅陡降到信号振幅的1/200甚至更小，进一步地说明了EMD分解的正确性。

仿真信号由EDM自适应地分解为6阶IMF，在时间域中(图1(a))表现为小尺度到大尺度层层分解，频率域(图1(b))则对应于从高频到低频的滤波过程。因此，IMF序列能更好地反映原始数据固有的物理特性，其分解是内在的和自适应的，每阶IMF 序列都代表了某种特定意义的频带信息，有利于实际应用与解释。

2.3  希尔伯特变换—HHT时频估计

对信号进行自适应EMD分解，得到一系列IMF分量(图1)后，对各阶IMF进行Hilbert变换，综合各阶瞬时频率谱得到信号的HHT时频谱(图2)。

(a) 信号和六阶IMF时域波形；(b)信号和六阶IMF的频谱

图1  仿真信号的EMD分解及各阶IMF的频谱

Fig.1  EMD decomposition and frequency spectrum of signal and IMF components

(a) 仿真信号的二维HHT时频谱；(b) 仿真信号的HHT边际谱；(c) 仿真信号的三维HHT时频谱

图2  仿真信号的HHT时频谱

Fig.2  HHT time-frequency spectrum of signal

从时频谱图(图2(a)和2(c))可看出2个已知成分(分别为带状和波浪状)的能量分布随时间和频率的动态变化特征，据谱值能区分出能量随时间和频率的细微变化。能量的变化不是连续的，而是离散的。带状和波浪状图形表示能量的分布相对集中，集中程度和形成的谱值空间图像正好反映了调频调幅成分与正弦信号成分的实际变化特征^{[9, 12]}。图2(b)所示为对应的HHT边际谱，反映了每个频率在整个时间长度内所累积的能量。可见，在没有能量的区域，谱值一般很小或为0，而不会像其他方法那样，因为加窗造成能量泄露，形成连续的空间时频谱值分布，并因此产生虚假信息，给解释工作带来困难。

已知数据的HHT时频谱图特征说明该方法对于时间域和频率域分辨率都很高，不受测不准原理限制，能量突变的定位与检测能力较强，具有很强的非稳态动态变换的时频刻画能力^[9]。

3  大地电磁信号分析应用实例

3.1  实测大地电磁信号

大地电磁(MT)信号具有非线性、非平稳和非最小相位特征，而Hilbert-Huang变换是最新发展起来的处理非线性、非平稳信号的时频分析方法^[13-14]。以实测EH-4观测信号(在内蒙古某地测得的EH-4仪器低频模式信号)为例，由于数据量非常大，这里只选取其中1个测点4个分量(E_x, E_y, H_x, H_y)中的1个磁场分量(H_x)进行分析。实测大地电磁信号如图3所示(采样频率为12 kHz)，检验HHT方法处理大地电磁资料的有效性和刻画信号非线性的能力。

图3  实测大地电磁信号

Fig.3  Original MT signal

3.2  大地电磁信号的经验模态分解

借助于经验模态分解方法(EMD)，可以获得有限数目的分段固有模态函数(IMF)。由于特征尺度参数是基于实际采集所获得的信号，因此，根据特征尺度参数对数据进行筛分所得到的分段IMF一般都具有明显的物理意义，每个IMF 表征信号在某一特征尺度参数上的模态。图4(a)所示为实测MT信号的EMD分解示意图，该信号依据EMD信号分解技术，自适应地分解为7阶固有模态函数(IMF)，在时间域中表现为小尺度到大尺度的层层滤波。图4(b)所示为信号与各阶IMF的频谱。显然，在频率域中，EMD分解过程表现为从高频到低频的层层滤波。由图4可见，在整个频段内，0~1 kHz成分的振幅较大，在1 kHz以上的振幅较弱，这正是由于数据采集使用的是EH-4仪器的低频模式，高频截止频率为1 kHz，也说明了EMD 分解的正确性。EMD是以信号的极值特征尺度为度量的筛分过程，信号从最小的特征尺度进行筛分，从而获得最短周期的固有模态函数(IMF)，随后，经过一层层筛分，可以获得周期长度逐渐增大的多阶IMF。并且信号经过EMD分解后能够完全重构，没有能量   损失。

(a) MT信号和7阶IMF时域波形；(b)MT信号和7阶IMF的频谱

图4  MT信号的EMD分解及各阶IMF的频谱

Fig.4  EMD decomposition and frequency spectrum of MT signal and IMF components

MT信号具有非平稳特征，数据采集过程可能出现某一时段信号强度大，信噪比较高，而在另一时段信号强度小，信噪比也较小的现象，这是由于所采取的方法是多次采集。但是，若所采集的大部分信号强度低，则会严重降低数据质量^[11]。若能够在时间序列中筛选信号强度大的时段资料进行数据处理，而舍弃信号强度小的时段资料，则必然会提高数据质量。基于傅里叶变换的传统方法难以完成上述过程，因为傅里叶变换是一种全局变换，要么完全在时域，要么完全在频域，因此，无法同时表述信号的时-频分布特性。而时频特性对非平稳信号分析来说是非常重要的。Hilbert-Huang变换则能够很好地解决这一问题^[8-14]。

3.3  大地电磁信号的HHT时频分析

图5所示是对图4(a)所示频谱经过Hilbert变换得到的时频谱。其中：图5(a)所示为采用二维平面等值线图表示能量变换特征的时频演化过程，图中各点表示能量，颜色越亮，表示能量越高，反之，能量越低；图5(b)所示为MT信号的HHT边际谱，对应图5(a)中每个频率在整个时间长度内所累积的能量，显然，能量集中在0~1 kHz频段内；图5(c)所示为用Hilbert变换计算的三维时频谱。可见，时频谱具有很好的时频聚集性，瞬时能量谱直观地表明能量主要不均匀地分布在几个时间区间内，能很好地区分能量随时间和频率的细微变化。显然，在0~1 kHz频段内，整个时域能量较强，且基本上均匀分布，在1 kHz以上能量较弱，且分布不均匀，呈毛绒状。

(a) MT信号的二维HHT时频谱；(b) MT信号的HHT边际谱；(c) MT信号的三维HHT时频谱

图5  MT 信号的HHT时频谱

Fig.5  MT signal spectrum plotted in time-frequency based on HHT

从上述大地电磁信号的HHT时频谱演化特征可以看出，大地电磁信号非线性动态变化特征可被HHT方法较好地刻画，体现在：

a. HHT方法对大地电磁信号的动态变化过程刻画得比较清楚，反映了大地电磁信号的非线性、非平稳性，每时段都有各自的频率特性和能量差异，其他方法难于揭示这些细微性变化。这种变换特征为对地下介质的研究提供了一种新的思路与方法。

b. 时频谱阵图上的能量主要集中在0~1 kHz内，这与实际EH-4低频模式的频率设置范围为0.001~1 kHz相符。

c. 在IMF序列曲线上和时频谱图像中， 大地电磁数据的突变点、持续时段和频带能量分布均能够清晰地显现，表明HHT方法具有完全局部时频特性。

4  结  论

a. EMD可自适应地从原始信号中按不同的时间尺度分解出反映原始信号本身固有特性的IMF。IMF分量通常具有实际的物理含义, 每阶IMF序列都代表了某种特定意义的频带信息，为高精度的大地电磁信号的时频分析提供了保证，也给大地电磁信号的信号检测、噪声压制提供了新的途径。

b. 大地电磁信号经过EMD分解获得的IMF序列具有窄带特性，可用Hilbert变换计算得到三维时频谱值，从而可用于大地电磁的功率谱估计和对阻抗进行进一步估算等。

c. HHT时频谱在大地电磁信号的分析中，能够反映大地电磁信号随时间和频率动态变化在不同阶段的主要特征，其时频局部定位能力较强, 且时频分辨率很高, 能够显示能量突变点信息，这为大地电磁信号的识别、数据筛选和分段平稳的研究提供了有利条件。

d. HHT方法不需要选择基函数，且完全取消了窗函数的作用，其结果不受核函数影响与时频测不准原理限制，具有完全的局部时频特性，可用于准确描述大地电磁信号的时变特征，有很强的时频刻画能力。HHT方法可用于大地电磁信号检测、噪声压制、功率谱计算、阻抗估计等。

参考文献：

[1] 王书明, 王家映. 大地电磁信号统计特征分析[J]. 地震学报, 2004, 26(6): 669-674.
WANG Shu-ming, WANG Jia-ying. Analysis on statistic characteristics of magnetotelluric signal[J]. Acta Seismologica Sinica, 2004, 26(6): 669-674.

[2] 王书明, 王家映. 关于大地电磁信号非最小相位性的讨论[J]. 地球物理学进展, 2004, 19(2): 216-221.
WANG Shu-ming, WANG Jia-ying. Discussion on the non- minimum phase of magnetotelluric signals[J]. Progress in Geophysics, 2004, 19(2): 216-221.

[3] Cohen L. Time-frequency analysis[M]. New Jersey: Prentice Hall, 1995.

[4] 谢 中, 程迎军, 徐清燕. 电解气泡析出时电位波动的频谱分析[J]. 中国有色金属学报, 2003, 13(4): 1011-1016.
XIE Zhong, CHENG Ying-jun, XU Qing-yan. Spectral analysis of potential fluctuations on electrolytic gas evolution[J]. The Chinese Journal of Nonferrous Metals, 2003, 13(4): 1011-1016.

[5] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-station time series analysis[C]//Proceeding of the Royal Society of London, 1998: 903-995.

[6] Yue H Y, Guo H D. A SAR interferogram filter based on the empirical mode decomposition method[C]//Proceedings of Geosoience and Remote Sensing Symposium. IGARSS O1, 2001: 206-2063.

[7] 李夕兵, 张义平, 左宇军, 等. 岩石爆破振动信号的EMD滤波与消噪[J]. 中南大学学报: 自然科学版, 2006, 37(1): 150-154.
LI Xi-bing, ZHANG Yi-ping, ZUO Yu-jun, et al. Filtering and de-noising of rock blasting vibration signal with EMD[J]. Journal of Central South University: Science and Technology, 2006, 37(1): 150-154.

[8] 汤井田, 化希瑞, 曹哲民, 等. Hilbert-Huang变换与大地电磁噪声压制[J]. 地球物理学报, 2008, 51(2): 603-610.

TANG Jing-tian, HUA Xi-rui, CAO Zhe-ming, et al. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data[J]. Chinese J Geophys, 2008, 51(2): 603-610.

[9] 武安绪, 吴培稚, 兰从欣, 等. Hilbert-Huang 变换与地震信号的时频分析[J]. 中国地震, 2005, 21(2): 207-216.
WU An-xu, WU Pei-zhi, LAN Cong-xin, et al. Hilbert-Huang transform and time-frequency analysis of seismic signal[J]. Earthquake Research in China, 2005, 21(2): 207-216.

[10] Livanos G, Ranganathan N J. Heart sound analysis using the S-transform[J]. IEEE Comp Cardiology, 2000, 27(3): 587-590.

[11] 林 君, 项葵葵, 朱宝龙, 等. MT信号现场处理的实现技术研究[J]. 数据采集与处理, 1997, 12(1): 52-55.
LIN Jun, XIANG Kui-kui, ZHU Bao-long, et al. Study on implementation of MT signal processing in sit[J]. Journal of Data Acquisition & Processing, 1997, 12(1): 52-55.

[12] Bradley M B, Camelia K. Application of the empirical mode decomposition and Hilbert-Huang transform to reflection seismic data[J]. Geophysics, 2007, 72(3): H29-H37.

[13] 张海勇. 一种新的非平稳信号分析方法—局域波分析[J]. 电子与信息学报, 2003, 25(10): 1327-1334.
ZHANG Hai-yong. A new method for non-stationary signal analysis—Local wave analysis[J]. Journal of Electronics and Information Technology, 2003, 25(10): 1327-1334.

[14] 盖 强, 张海勇, 徐晓刚. Hilbert-Huang变换的自适应频率多分辨分析研究[J]. 电子学报, 2005, 33(3): 563-566.
GAI Qiang, ZHANG Hai-yong, XU Xiao-gang. Study of adaptive frequency multi-resolution analysis of the Hilbert-Huang transform[J]. Journal of Electronics, 2005, 33(3): 563-566.

收稿日期：2008-10-28；修回日期：2009-01-07

基金项目：国家高技术研究发展计划(“863”计划)项目(2006AA06Z105)

通信作者：蔡剑华(1979-)，男，湖南郴州人，博士研究生，讲师，从事电法勘探及信号处理研究；电话：0736-7186121；E-mail: cjh1021cjh@163.com