2017, Vol. 32
2. 湖南财政经济学院信息技术与管理学院, 长沙 410205
2. School of Information Technology and Management, Hunan University of Finance and Economics, Changsha 410205, China
大地电磁测深法中有一个重要的研究内容就是电磁场信号的功率谱估计，从测量数据通过估计不同频率的功率谱进而计算其阻抗，就可以推断出不同深度的大地电阻率(Wight et al., 1977；陈乐寿和王光锷，1990).因而正确地估计电场和磁场的功率谱对于改善大地电磁测深的阻抗估计具有重要的意义.目前，最传统也是最主要的功率谱估计方法还是基于Fourier变换的谱估计方法(Kaufman and Keller, 1981；王书明和王家映，2004a)，但随着人们对大地电磁信号和噪声特性的不断深入研究，发现以分析平稳信号为基础的Fourier变换与大地电磁信号的非平稳特性是不相适宜的，必然导致估算带来偏差，尤其是在矿集区和存在近场源干扰的情况下，情况就显得尤为突出(汤井田等，2012；王辉等，2016).在抑制噪声方面，经验模态分解、数学形态滤波以及相结合的远参考方法被成功应用(汤井田等，2012；王大勇等，2015).在谱估计方面，人们也陆续提出了基于小波变换和高阶统计量等的一些现代谱估计方法，他们克服了传统谱估计方法中的一些不足，提高了谱估计的精度(王书明和王家映，2004b；严家斌等，2008).Hilbert-Huang变换作为一种新发展起来的非线性、非平稳信号处理工具，在大地电磁信号处理方面的应用也得到了不同程度的研究，取得了一系列的成果(Huang et al，1998；汤井田等，2008).但其适用范围，尤其是在MT信号的处理中与传统Fourier变换的对比研究，文献中报道的不多.作者试图以仿真和实测大地电磁信号为例，从准确性、稳定性、计算效率等几个方面，在与Fourier变换对比的基础上，研究Hilbert-Huang变换在大地电磁信号处理中的适用性，以供参考和讨论.
1 傅里叶变换与HHT算法简介 1.1 傅里叶变换傅里叶变换的理论可简单描述如下：若f(t)∈L2(R)，则存在Fourier变换对为
|
(1) |
|
(2) |
称f(ω)为f(t)的Fourier变换，且有f(ω)∈L2(R).离散Fourier变换对(DFT)为(Wight et al., 1977)：
|
(3) |
|
(4) |
式中x(n)为时域离散信号，X(k)就称为信号x(n)频域上的Fourier频谱.工程中通常利用傅立叶快速算法(FFT)来实现DFT.
1.2 HHT方法HHT方法处理信号的基本过程分为两部分：经验模态分解(EMD)和Hilbert谱分析(Huang et al，1998, 2003).
(1) EMD分解是一个循环筛分的过程，它将原始信号分解成若干固有模态函数(Intrinsic Mode Function, IMF)分量与一个残余分量之和.IMF满足以下两个条件: (a)极大值点与极小值点的包络线关于时间轴对称; (b)极大点和极小点的个数与过零点个数相等或至多相差1.分解后，原始信号x(t)可表示为所有的IMF及剩余量之和(Huang et al., 1998)，公式为
|
(5) |
(2) Hilbert变换.获得IMF分量以后，就可以对每一阶IMF作Hilbert变换, 得到每一个IMF随时间变化的瞬时频率和瞬时幅值, 综合所有IMF的Hilbert瞬时谱得到原始信号的Hilbert谱(Huang et al., 2003；Battista et al., 2007)，公式为
|
(6) |
进一步对时间积分，可以得到Hilbert边际谱h(ω)为
|
(7) |
式中，T为信号的长度.边际谱是Hilbert谱沿时间轴的积分，它代表了每个频率在统计意义上的全部累加幅值(Cai et al，2009).
2 傅里叶变换与HHT应用的比较研究以仿真的平稳信号、加噪信号和非平稳信号为对象，对傅里叶变换与HHT算法性能进行比较研究.
2.1 对平稳信号的处理设计了一个含有两种成分的电压信号，其表达式为
|
(8) |
合成信号包含一个150 Hz的正弦波平稳信号和一个基频为50 Hz，调制频率为15 Hz的调频调幅信号，信号采样率为1 kHz，时间长度为1 s.
为说明Hilbert谱的特点，图 1a给出了信号的三维Hilbert谱，图 1b为对应的二维Hilbert谱(图 1a在X-Y平面的投影).从时频谱图中可明显看出2个分别呈带状和齿状的能量分布，信号能量随时间和频率的细微变化被刻画的非常精细.能量集中在两个能量带上，集中程度和形成的谱值空间图像正好反应映了调频调幅成分与正弦信号成分的实际变化特征.图 1c为对应的Hilbert谱沿时间轴积分得到的HHT边际谱，图 2a为其放大图，表达了每个频率在整个时间长度内所累积的能量.可见，在50 Hz和150 Hz附近有两个能量很强的谱线，尤其在以50 Hz为中心频率的15 Hz带宽内都有能量分布，这正确反应了调频调幅信号的被调制实质.图 2b为FFT方法计算的功率谱，图中在两段频段内存在能量很强的谱线，150 Hz的谱线精确反应了正弦稳态信号，在50 Hz、35 Hz和65 Hz附近有三根谱线，而在其他频点没有反应能量的谱线存在，这是不符合调频调幅信号的能量分布规律的.
|
图 1 仿真平稳信号的Hilbert谱 (a)三维谱图; (b)二维Hilbert谱; (c)边际谱. Figure 1 Hilbert spectra of the simulated stationary signals (a) Three-dimensional spectra; (b) Two-dimensional Hilbert spectra; (c) Marginal spectra. |
|
图 2 仿真平稳信号的功率谱 (a)HHT边际谱; (b)FFT方法计算的功率谱. Figure 2 Power spectrum of the simulated stationary signal (a) HHT marginal spectrum; (b) Power spectrum calculated by FFT method. |
可见，对于平稳信号，边际谱和傅立叶谱都能非常精确的给出信号的频谱.对于调频调幅信号，Hilbert时频谱及其边际谱都能正确反映其时频特征，而FFT不能正确表达出调频调幅信号的能量分布规律.再者，傅里叶变换原理可以知道，傅里叶谱仅能从频域上表征信号的特征，不能同从时域上得到信号能量随时间的变换情况.而HHT时频谱图对于时间域和频率域的分辨率都很高，能量突变点的定位与检测能力较强，具有很强的非稳态动态变换的时频刻画能力，这在大地电磁信号检测、时频分析等方面有着广阔的应用前景.
2.2 对含白噪声信号的处理大地电磁数据采集的时候，工作环境恶劣，信道有可能受到外界电磁干扰的影响，表现为叠加在原始信号上的随机噪声扰动(白噪声)(Banks，1998；蔡剑华等，2013).用含三个频率成分(分别为2 Hz，20 Hz和100 Hz)的正弦测试信号上叠加-3 dB的随机白噪声，模拟多分量的含噪大地电磁信号，抽样频率为1 kHz，信号长度为1s.时域波形如图 3a中的”original signal”所示.对于平稳信号，最快捷、最常用的谱估计方法还是傅立叶变换，HHT算法具有滤除高频扰动的能力，因而有较强的抗干扰性能，下面分别用FFT和HHT方法对信号进行谱估计.
|
图 3 加噪信号的处理 (a)EMD分解；(b)对应各阶IMF的功率谱；(c)信号的Hilbert谱. Figure 3 Procession of noisy signal (a) EMD decomposition; (b) Power spectrum of each order IMF; (c) Hilbert spectrum of signal. |
对信号进行EMD分解，结果如图 3所示.图 3a为分解得到的各阶IMF分量，从高频到低频，原信号被分解得到四个分量：“imf1”、“imf2”、“imf3”和“imf4”，“imf1”为高频噪声分量，“imf2”-“imf4”是信号的主导成分.图 3b为分解所得各分量对应的频谱.“imf1”对应的频谱图频带很宽，幅值较小，为随机噪声的频谱，“imf2”-“imf4”很精确的反应了2 Hz，20 Hz和100 Hz的频率成份.可见，EMD能很好的将噪声与三个成分分别开来.图 3c为信号的Hilbert时频谱，呈现了噪声和三个分量的能量分布.摈弃“imf1”，将“imf2”-“imf4”叠加重构信号，再计算其边际谱如图 4a所示.为作对比，图 4b给出了FFT方法计算的功率谱.通过比较两种算法得到的功率谱可知：HHT算法本身具有滤除高频分量的能力，而FFT算法对噪声极为敏感，需要增加前置的滤波器，滤除输入数据中的噪声干扰.
|
图 4 加噪信号的功率谱 (a)HHT边际谱; (b)FFT方法计算的功率谱. Figure 4 Power spectrum of noisy signal (a) HHT marginal spectrum; (b) Power spectrum calculated by FFT method. |
大地电磁信号是典型的非线性、非平稳信号，为了检验前述方法用于大地电磁信号谱分析的有效性, 文中设计了4个振幅为1, 频率分别为20 Hz、17.5 Hz、15 Hz和11.25 Hz的合成信号, 此频率序列与我国引进较多的大地电磁测深仪EH-4系统的分析频点是对应的.该仿真信号在时域上分为4段, 4个频率成分各占总数据长度的1/4, 模拟MT信号的非平稳性, 整个时间长度为2s, 采样频率为fc=50 Hz, 采样点数N=100，每段为25个点，信号为
|
(9) |
先用傅立叶方法估算功率谱.由傅里叶变换原理和抽样定理可知，用FFT方法估计该仿真信号的谱，其频率分辨率
|
图 5 非平稳信号的功率谱 (a)HHT边际谱; (b)FFT方法计算的功率谱. Figure 5 Power spectrum of the non-stationary signal (a) HHT marginal spectrum; (b) Power spectrum calculated by FFT method. |
用HHT方法进行对比研究.先用HHT方法得到信号的时频谱图如图 6a所示，从图中可明显看出4个成分的能量分布随时间和频率的动态变化特征，谱图从时域和频域两方面把信号在4个时间段上的频率特征准确的区分出来了.由时频谱图沿时间轴积分得到边际谱如图 6c所示.图 5b是图 6c的放大图.显然，图中4个频率已经分开，已将瞬态频率较准确的放映出来了(图中出现的毛刺与低频信号交接采样点的选取和EMD算法中的端点效应等有关).
|
图 6 非平稳信号的Hilbert谱 (a)三维谱图; (b)二维Hilbert谱; (c)边际谱. Figure 6 Hilbert spectra of the simulated non-stationary signals (a)Three-dimensional spectra; (b) Two-dimensional Hilbert spectra; (c) Marginal spectra. |
以最大幅值对应的频点作为此频率成分的有效频点，表 1给出了傅立叶方法和HHT方法频率成分估算准确度的比较.可以看出，HHT的方法准确度高于傅立叶方法，对于非平稳信号的分析(本例所用的信号严格的说是分段平稳信号)，HHT方法得到的边际谱是有明显优势的，可以减少信号非平稳性带来的估算偏差.对于多成分且非平稳的大地电磁信号，边际谱更符合实际.
|
表 1 仿真非平稳信号FFT和HHT谱估计精度比较 Table 1 Comparison of spectral estimation accuracy of FFT and HHT for the simulated non-stationary signals |
通过以上的比较研究可以得出：HHT的边际谱是一种合理的新的信号谱估计方法，它由Hilbert谱沿时间轴积分得到，代表了整个信号在时间跨度上的幅值累积效应，突破了Fourier变换只对平稳信号有效的不足之处.所以HHT的边际谱更适合对非线性、非平稳的大地电磁信号进行谱估计，且频率分辨率要比傅氏谱的高.
3 实测大地电磁信号分析以安徽某测点测得的大地电磁信号为例，测点位于安徽泥庐枞集区.庐江—枞阳矿集区位于安徽省境内，经过三县两市辖区，涉及约40个乡镇.测区内人口稠密、水系发育、交通网密布、通信电力网发达，另有较多的矿山正在开采，错综复杂的干扰源，为MT数据的采集和处理带来了许多困难，大地电磁信号受人文噪声干扰严重.本次作业，在区内部署了5条综合地球物理测线，测线总长约325 km，设计MT测点655个.文中截取了一屏含4个分量的大地电磁信号，如图 7所示.从信号的时域图可以看出，信号受噪声的干扰很大，几乎淹没了有用信号的信息.下面就以此信号为例，对比分析大地电磁信号的Fourier谱和边际谱.为说明问题，图 8给出了一个磁场分量Hy和一个电场分量Ex的HHT时频谱，从图中可以看出HHT方法对大地电磁信号的动态变化过程刻画得比较清楚，大地电磁数据的突变点、持续时段和频带能量分布均能够清晰地显现，每时段都有各自的频率特性、能量差异，其他方法难于揭示出这些细微性变化.可见，从时频谱图沿时间积分得到的边际谱，是真实的具有此频率成分的各时间点能量的集合，其结果可更接近大地电磁信号的实际频谱.
|
图 7 实测大地电磁信号的时域波形图 Figure 7 Time domain waveform of the measured MT signal |
|
图 8 实测大地电磁信号的HHT时频谱 (a) Hy分量; (b) Ex分量. Figure 8 Hilbert time-frequency spectra of the measured MT signal (a) Hy component; (b) Ex component. |
图 9的左右两组图，给出了该测点大地电磁场信号4个分量的边际谱和傅氏谱.显然，实测信号的边际谱和傅氏谱线的形状和趋势基本一致，但细节上具有明显的区别.Fourier幅值谱表现为MT信号的能量分散在比较宽的频率上，但由边际谱看，MT信号的能量主要集中在低频部分，这更符合大地电磁信号的实际；边际谱的谱线较Fourier谱线平滑，个别频点幅值突出的情况较少，谱线鲁棒性较好，尤其在高频段，边际谱较平稳的特征更为明显；再者，从信号的傅氏谱可以看到电场信号谱线的平稳性较磁场信号谱线平稳性更差，某些频点幅值突出的情况更为严重，这说明了MT信号中电场信号比磁场信号更容易受噪声污染.而在边际谱中这种情况得到了较好的压制，电、磁场信号谱线都较傅氏谱平稳.因此，也进一步证明了边际谱比傅氏谱更具抗干扰能力，更能反应大地电磁信号频率分布的真实情况.
|
图 9 实测MT信号的功率谱 (a)HHT边际谱; (b)FFT方法计算的功率谱. Figure 9 Power spectrum of the measured MT signal (a) HHT marginal spectrum; (b) Power spectrum calculated by FFT method. |
为进一步比较两种方法在谱估计上的计算效率，表 2给出了傅立叶方法和HHT方法估算不同数据长度大地电磁信号的功率谱所消耗的时间对比(计算机配置：CPU：AMD 5200+, 主频：2.7 GHz, 内存：2 G).由表中参数可以看出，相同数据长度下, 基于傅立叶变换的谱估计方法相对HHT方法速度要快，且HHT分析耗时基本与数据长度成线性关系；而对于FFT方法，在数据长度不是很大的情况下，数据点数的多少对速度的影响不是很大.这主要是由于HHT方法需要对信号进行EMD分解，再由得到的各阶IMF分量计算得到Hilbert时频谱进而估算边际谱，过程中需要进行EMD循环分解和Hilbert变换两个步骤耗时比较大.可见，对于大地电磁测深，尤其是现场作业，需要对海量数据进行处理，HHT方法的计算效率将会是制约其在MT领域推广的瓶颈问题.
|
|
表 2 MT信号FFT和HHT谱估计耗时比较 Table 2 Time consuming comparison of FFT and HHT method |
基于傅里叶变换的谱估计方法适合于对平稳信号进行处理，计算速度较快；目前工程应用中的大地电磁测深还是用傅里叶方法计算信号的功率谱，再带入方程进而计算出不同频率的响应参数，而傅里叶方法对所分析的信号必须平稳的要求与大地电磁信号的非平稳特性是不相适宜的，必然导致估算带来偏差.HHT算法具有时频分析和滤除高频分量的能力，对信号能量有很高的时频刻画能力，为MT信号的时频分析与数据筛选提供了条件；基于HHT边际谱的功率谱估计更适合MT信号非平稳特性的实质，从方法上减小了谱估算的误差；但HHT方法计算效率低，是制约其工程应用的瓶颈，寻求快速算法(如曲线拟合、端点效应抑制和筛选停止准则等几个问题的优化处理)是迫切需要解决的问题.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Banks R J. 1998. The effects of non-stationary noise on electromagnetic response estimates[J]. Geophys. J. Int., 135(2): 553–563. DOI:10.1046/j.1365-246X.1998.00661.x |
| [] | Battista B M, Knapp C, McGee T, et al. 2007. Application of the empirical mode decomposition and Hilbert-Huang transform to seismic reflection data[J]. Geophysics, 72(2): H29–H37. DOI:10.1190/1.2437700 |
| [] | Cai J H, Tang J T, Hua X R. 2009. An analysis method for magnetotelluric data based on the Hilbert-Huang Transform[J]. Exploration Geophysics, 40(2): 197–205. DOI:10.1071/EG08124 |
| [] | Cai J H, Wang X C, Hu W W. 2013. A method for MT data denoising based on empirical mode decomposition and wavelet threshold[J]. Oil Geophysical Prospecting, 48(2): 303–307. |
| [] | Chen L S, Wang G E. 1990. Magnetotelluric Method[M]. Beijing: Geological Press: 12-40. |
| [] | Huang N E, Shen Z, Long S R, et al. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationarity time series analysis[J]. Proc. Roy. Soc. A Math. Phys. Eng. Sci., 454(1971): 903–995. DOI:10.1098/rspa.1998.0193 |
| [] | Huang N E, Wu M L C, Long S R, et al. 2003. A confidence limit for the empirical mode decomposition and the Hilbert spectral analysis[J]. Proc. Roy. Soc. A Math. Phys. Eng. Sci., 459(2037): 2317–2345. DOI:10.1098/rspa.2003.1123 |
| [] | Kaufman A A, Keller G V. 1981. The Magnetotelluric Sounding Method[M]. Amsterdam: Elsevier: 595-598. |
| [] | Tang J T, Hua X R, Cao Z M, et al. 2008. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data[J]. Chinese J. Geophys., 51(2): 603–610. |
| [] | Tang J T, Liu Z J, Liu F Y, et al. 2015. The denoising of the audio magnetotelluric data set with strong interferences[J]. Chinese J. Geohpys., 58(12): 4636–4647. DOI:10.6038/cjg20151225 |
| [] | Tang J T, Xu Z M, Xiao X, et al. 2012. Effect rules of strong noise on magnetotelluric (MT) sounding in the Luzong ore cluster area[J]. Chinese J. Geophys., 55(12): 4147–4159. DOI:10.6038/j.issn.0001-5733.2012.12.027 |
| [] | Wang D Y, Zhu W, Fan C S, et al. 2015. Noise processing methods and application study of MT in the ore concentration area[J]. Geophysical and Geochemical Exploration, 39(4): 823–829. |
| [] | Wang H, Cheng J L, Teng X Z, et al. 2016. Source effect on magnetotelluric data due to mining area and its suppression[J]. Progress in Geophysics, 31(3): 1358–1366. DOI:10.6038/pg20160359 |
| [] | Wang S M, Wang J Y. 2004a. Analysis on statistic characteristics of magnetotelluric signal[J]. Acta Seismologica Sinica, 26(6): 669–674. |
| [] | Wang S M, Wang J Y. 2004b. Application of higher-order statistics in magnetotelluric data processing[J]. Chinese J. Geophys., 47(5): 928–934. DOI:10.3321/j.issn:0001-5733.2004.05.027 |
| [] | Wight D E, Bostick F X, Smith H W. 1977. Real time fourier transformation of magnetotelluric Data[R]. EGRL/EERL Technology Report, Texas:University of Texas. |
| [] | Yan J B, Liu G Z, Liu J X. 2008. Application of wavelet transform in processing nature electromagnetic field time series[J]. Geology and Prospecting, 44(3): 75–78. |
| [] | 蔡剑华, 王先春, 胡惟文. 2013. 基于经验模态分解与小波阈值的MT信号去噪方法[J]. 石油地球物理勘探, 48(2): 303–307. |
| [] | 陈乐寿, 王光锷. 1990. 大地电磁测深法[M]. 北京: 地质出版社: 12-40. |
| [] | 汤井田, 化希瑞, 曹哲民, 等. 2008. Hilbert-Huang变换与大地电磁噪声压制[J]. 地球物理学报, 51(2): 603–610. |
| [] | 汤井田, 刘子杰, 刘峰屹, 等. 2015. 音频大地电磁法强干扰压制试验研究[J]. 地球物理学报, 58(12): 4636–4647. DOI:10.6038/cjg20151225 |
| [] | 汤井田, 徐志敏, 肖晓, 等. 2012. 庐枞矿集区大地电磁测深强噪声的影响规律[J]. 地球物理学报, 55(12): 4147–4159. DOI:10.6038/j.issn.0001-5733.2012.12.027 |
| [] | 王大勇, 朱威, 范翠松, 等. 2015. 矿集区大地电磁噪声处理方法及其应用[J]. 物探与化探, 39(4): 823–829. DOI:10.11720/wtyht.2015.4.27 |
| [] | 王辉, 程久龙, 腾星智, 等. 2016. 矿区近场源噪声对大地电磁测深数据的影响及其压制方法[J]. 地球物理学进展, 31(3): 1358–1366. DOI:10.6038/pg20160359 |
| [] | 王书明, 王家映. 2004a. 大地电磁信号统计特征分析[J]. 地震学报, 26(6): 669–674. |
| [] | 王书明, 王家映. 2004b. 高阶统计量在大地电磁测深数据处理中的应用研究[J]. 地球物理学报, 47(5): 928–934. DOI:10.3321/j.issn:0001-5733.2004.05.027 |
| [] | 严家斌, 刘贵忠, 柳建新. 2008. 小波变换在天然电磁场信号时间序列处理中的应用[J]. 地质与勘探, 44(3): 75–78. |