基于高阶统计量的大地电磁数据处理与仿真

蔡剑华^{1, 2}，胡惟文¹，任政勇²，雷立云¹

(1. 湖南文理学院 信息研究所，湖南 常德，415000；

2. 中南大学 信息物理工程学院，湖南 长沙，410083)

摘  要：针对大地电磁(MT)数据具有非线性、非平稳的特征而不能满足传统功率谱方法要求的问题，从高阶统计量的理论入手，提出由信号的高阶谱恢复功率谱，再由恢复的功率谱估算MT响应函数的高阶统计量方法，建立通过高阶谱恢复功率谱的数学模型，并对合成数据和实测MT数据进行仿真分析。研究结果表明：利用高阶累积量对高斯噪声不敏感的特点，可实现在高斯白噪声和有色噪声污染下的瞬态信号频率与功率谱的正确估计；高阶统计方法消除了传统功率谱方法中大地电磁响应函数出现个别频点分散和形态异常的现象，使地质参数的估计精度和资料的可解释性得到明显提高。

关键词：高阶统计量；大地电磁信号；数据处理；噪声抑制

中图分类号：P631          文献标志码：A         文章编号：1672-7207(2010)04-1556-05

Magnetotelluric data processing and simulation based on higher-order statistics

CAI Jian-hua^{1, 2}, HU Wei-wen¹, REN Zheng-yong², LEI Li-yun¹

(1. Information Institute, Hunan University of Arts and Science, Changde 415000, China;

2. School of Info-physics and Geomatics Engineering, Central South University, Changsha 410083, China)

Abstract: Based on the fact that magnetotelluric data (MT) series are non-linear and non-stationary, and donot meet the basic requirements of the traditional power spectrum, the MT data processing method based on the high-order statistics was proposed. This method used higher-order spectra to reconstruct power spectra and estimated MT responses with the reconstructed power spectra. A model was established using higher-order spectra to reconstruct power spectra. Meanwhile, simulations were conducted using the man-made data and recorded MT data. The results show that the presented method is superior to the traditional power spectral method in suppressing Gaussian colored noise and can extract more useful information. Compared with the traditional method, it eliminates scattering at individual frequencies with shape distortion in MT response curves. Therefore, the accuracy of estimation and the interpretability of data are significantly improved.

Key words: higher-order statistics; magnetotelluric signal; data processing; noise suppression

在大地电磁测深法中，以往基于二阶统计量(如二阶矩阵、相关函数和功率密度函数等)的各种信号处理方法，对信号和地质模型作了许多假设和要求，如假设大地电磁信号为高斯信号，噪声为高斯白噪声，大地系统为线性最小相位系统等^[1-2]，但实际情况往往并非如此。许多学者已经证明大地电磁信号具有非线性、非平稳和非最小相位特征，是一种典型的非平稳随机信号^[3-5]，由于传统功率谱方法本身内在的计算方式不能有效抑制高斯噪声，致使MT 解释结果较差。高阶统计量是近20多年来发展起来的一种新的信号分析和处理理论，它与传统功率谱法相比具有显著的优点：(1) 高阶累积量对高斯有色噪声恒为0 dB；(2) 高阶累积量含有系统的相位信息^[6-7]，这就降低了对信号和模型的要求，如信号可以是非高斯性，地质模型可以是非最小相位系统 。高阶统计信号处理已经渗透到信号处理的各个应用领域^[5-9]，而国内把高阶统计量应用于MT资料处理的研究较少。本文作者针对大地电磁测深法，把高阶统计量方法引入MT资料处理中，通过数值模型和仿真试验，探讨高阶统计量在MT资料处理中的适用性和抑制噪声的能力。

1  高阶统计量与高阶谱

高阶统计量是描述随机过程高阶(二阶以上)统计特性的一种数学工具，包括高阶累积量和高阶矩。与自相关函数的傅里叶变换定义为函数的功率谱类似，高阶累积量的多维傅里叶变换定义为高阶谱(或称多谱)^[10-12]。

定义1  设{x(n)}为零均值的k阶平稳随机过程，则该过程的k阶矩m_kx(τ₁, …, τ_k-1) 定义为：

m_kx(τ₁,…, _k-1)=mon{x(n), x(n+τ₁), …, x(n+τ_k-1)}   (1)

而k阶累积量c_kx(τ₁, …, τ_k-1)定义为：

c_kx(τ₁, …, τ_k-1)=cum{x(n), x(n+τ₁), …, x(n+τ_k-1)}  (2)

式中：mon(…)，cum(…)为求k元实变向量的矩和累积量。对于一个零均值的平稳随机过程{x(n)}，其高阶累计量也可定义为：

c_kx(τ₁, …, τ_k-1)=E{x(n), x(n+τ₁), …, x(n+τ_k-1)} -E{g(n),

g(n+τ₁), …, g(n+τ_k-1)}            (3)

其中：k≥3；{g(n)}是一个与{x(n)}具有相同功率谱的高斯过程。

定义2  设高阶累积量c_kx(τ₁, …, τ_k-1)是绝对可和的，若

＜，       (4)

则k阶累积谱定义为k阶累积量的k-1维傅里叶变换，即

   (5)

式中：为频率。习惯上，高阶累积谱简称作高阶谱或多谱，最常用的高阶谱是三阶谱(也称双谱)：

   (6)

和四阶谱(也称三谱)：

  (7)

2  数值模型及估计算法

一个线性的非高斯信号的高阶谱总可以写成^[13]：

S_kx(ω₁, …, ω_k-1)=

γ_kuH(ω₁)…H(ω_k-1)×H^*(ω₁+…+ω_k-1)       (8)

其中：γ_ku是1个与输入信号u 对应的标量常数；H(ω)是线性时不变系统的频率传递函数。若非高斯信号x(n)是由非高斯白噪声u(n)激励一个线性时不变系统所产生，则上式对所有阶数均成立，并且有

P_x(ω₁) =γ_2u|H(ω₁)|²                (9)

B_x(ω₁, ω₂)=γ_3u H(ω₁)H(ω₂)H^*(ω₁+ω₂)     (10)

T_x(ω₁, ω₂, ω₃)=γ_4u H(ω₁)H(ω₂)H(ω₃)×

H^*(ω₁+ω₂+ω₃)                (11)

式中：P_x, B_x和T_x分别是信号的功率谱，双谱和三谱。综合考虑以上3式，一个线性过程的功率谱可以从它的双谱和三谱重构出来，其结果最多相差1个常数(假定H(0)≠0)，即有

B_x(ω, 0)=(γ_3u /γ_2u)P_x(ω)H(0)           (12)

T_x(ω, 0, 0)=(γ_{4 u}/γ_{3 u})P_x(ω)H²(0)          (13)

式中：P_x为功率谱。当加性噪声为高斯有色噪声    时，利用双谱或三谱重构功率谱比直接估计功率谱具有更好的性能。因为双谱和三谱对高斯有色噪声有抑制作用，而功率谱对任何有色噪声都是敏感的。其中常数γ_2u, γ_3u和γ_4u可由信号的累积量和系统的冲激响应来确定^[14-15]。

考虑到若信号是对称分布的，则其双谱很小或近似等于0，所以，本文作者采用三谱的一维切片重构功率谱：

P_x(ω, 0)=(γ₂ /γ_4u)T_x(ω, 0, 0)H²(0)      (14)

式中：γ_2u=cum(0)和γ_4u=cum(0, 0, 0)分别表示0延迟的二阶累积量和四阶累积量；H₂(0)是1个常量；三谱 T_x(ω, 0, 0)由四阶累积量的一维切片c_4x(τ, 0, 0)的傅里叶变换得到。采用互三谱的一维切片重构互功率谱：

P_xy(ω)=(γ_{xy2u /}γ_xyxx4u)T_xyxx (ω, 0, 0)H²(0)     (15)

式中：γ_xy2u=cum(0)和γ_xyxx4u=cum(0, 0, 0)分别表示0延迟的二阶互累积量和四阶互累积量；H²(0)是1个常量；γ_xyxx4u(0, 0, 0)由互四阶累积量的一维切片c_4xyxx(τ, 0, 0)的傅里叶变换得到。至此，利用信号的高阶谱重构其功率谱，然后，利用已有的方法估算视电阻率和相位等MT响应参数，该算法的流程如图1所示。先由大地电磁信号的时间序列计算得到双谱和多谱，再由它们重构得到功率谱，最后，代入已有的计算公式计算阻抗张量，进一步估算视电阻率和相位。

图1  高阶统计量法流程图

Fig.1  Analysis process diagram of higher-order statistics method

3  数值仿真与分析

为突出本文作者所提出的方法的应用效果，将仿真结果与利用传统功率谱方法得到的结果进行比较，传统方法采用周期图谱估算方法。

3.1  合成数据试验

首先，用合成数据进行试验，检验由高阶谱重构功率谱模型的有效性。合成数据实验信号包含频率分别为20, 70, 200, 500, 700, 1 000和3 000 Hz的正弦信号，加有独立同分布噪声、非线性漂移和直流分量。

图2(a)和2(b)所示分别为基于周期图和高阶统计量的功率谱。对比图2(a)和2(b)可以看出：用高阶统计量法估计出了淹没在独立同分布噪声中的所有频率成分，而在基于传统周期图法的频谱图中相应频率的幅度较小，尤其在频率为1 kHz和3 kHz时。显然，高阶统计量能够有效地识别出在高斯强干扰下的非高斯微弱信号，提高信噪比。在对信号进行分析时，不但可以反映信号的整体与局部变化，而且对信号有更高的分辨率。可见：高阶统计量可用在强噪声背景下提取MT数据的有用信号，为有效地识别地质体提供一种新方法。

图2  合成数据的功率谱

Fig.2  Power spectrum of simulated signal

3.2  实测MT资料分析

3.2.1  功率谱重构

图3和图4所示分析数据采用内蒙毛登206线50号点EH-4采集数据。大地电磁测深利用4个场资料进行谱估计，分别是x和y方向的2个电场分量(E_x, E_y)和2个磁场分量(H_x, H_y)。图3~4所示显示了H_y的自谱和它与H_x的互谱，其他分量的功率谱与此类似。分别用周期图法和高阶统计量法估计低频段数据的自功率谱和互功率谱。自功率谱周期图法估计采用分段平均的办法，分段作快速傅里叶变换，取各段的平均值作为估计的结果。高阶统计量估计采用本文提出的信号的高阶谱重构其功率谱的模型，由三谱的一维切片重构功率谱。从图3可以看出：相比传统功率谱相比，基于高阶统计量重构的功率谱幅值较大，消除了传统谱中50 Hz等频率成分的奇异值，更多有用成分被识别，说明高阶累积法对噪声被较好地压制，可以检测和识别微弱信号，提高信噪比。同时，从图4可知：EH-4数据在中频段(1~3 kHz)信噪比是很低的，采用高阶累积量估计，信噪比有较大提高。

图3  实测信号H_y的自功率谱

Fig.3  Self-power spectrum of practical signal H_y

3.2.2  视电阻率估算

分别用恢复的功率谱和传统功率谱估算视电阻率。仿真结果如图5所示。从图5可知：由传统功率谱得到的大地电磁响应函数出现了个别频点分散、形态异常等现象，在进行反演解释时，这使得许多地质特征难以有效提取出来；而高阶统计量法的数据比传统谱法的数据更集中，曲线更平稳。就整个曲线形态来看，估算的参量比用常规方法处理的结果更理想，资料的可解释性得到明显提高。

在大地电磁测深中，由信号的功率谱估算视电阻率时，通常用到4 种稳健的计算方法^[2]，并把这4 种方法计算结果的平均值作为最终结果。为了定量评估本文提出的方法的效果，引入全信息矢量相干度C_P：

              (16)

式中：R 表示4 种算法结果的平均值；D 表示所有4种算法结果与R之间差的模的平均值。若噪声影响越小，则C_P越接近1。经计算，用本文提出的方法估算的视电阻率ρ_xx，ρ_xy，ρ_yx和ρ_yy的C_P分别为0.923，0.986，0.954和0.937；用传统方法得到的相应的C_P分别为-1.280，0.690，-1.420和0.920。由此可知，本文作者提出的方法的抗干扰能力强于传统功率谱方法的抗干扰能力，因此，通过该方法处理的视电阻率曲线更加合理。

图4  实测信号H_y与H_x的互功率谱

Fig.4  Cross-power spectrum of practical signal of H_y and H_x

图5  估算的视电阻率分量比较

Fig.5  Contrast of apparent resistivities between classical method and herein proposed method

4  结论

(1) 大地电磁信号具有非线性、非平稳和非最小相位特征，是一种典型的非平稳随机信号，不满足以傅里叶变化为基础的传统方法的要求。高阶累积量具有对高斯有色噪声恒为0 dB的特点且含有系统的相位信息, 降低了资料处理中对MT信号和地质模型的要求。

(2) 本文所提出的由信号的高阶谱恢复功率谱的模型是有效的, 由恢复的功率谱估算MT响应函数，在抑制高斯有色噪声和提取信号中有用信息方面优于传统功率谱方法，消除了传统功率谱方法中大地电磁响应函数出现个别频点分散、形态异常的现象, 给实际应用与地质解释工作带来了方便。

参考文献：

[1] 陈乐寿, 王光锷. 大地电磁测深法[M]. 北京: 地质出版社, 1990: 12-40.
CHEN Le-shou, WANG Guang-e. Magnetotelluric method[M]. Beijing: Geogloical Press, 1990: 12-40.

[2] 石应骏, 刘国栋, 吴广耀, 等. 大地电磁测深法[M]. 北京: 地震出版社, 1984: 22-36.
SHI Ying-jun, LIU Guo-dong, WU Guang-yao, et al. Magnetotelluric method[M]. Beijing: Seismological Press, 1984: 22-36.

[3] 王书明, 王家映. 大地电磁信号统计特征分析[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.

[4] 王书明, 王家映. 关于大地电磁信号非最小相位性的讨论[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.

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

[6] 唐斌, 尹成. 基于高阶统计的非最小相位地震子波恢复[J]. 地球物理学报, 2001, 44(3): 404-410.
TANG Bin, YIN Cheng. Non-minimum phase seismic wavelet reconstruction based on higher order statistics[J]. Chinese J Geophys, 2001, 44(3): 404-410.

[7] Lazear G D. Mixed-phase wavelet estimation using fourth-order cumulants[J]. Geophysics, 1993, 58(7): 1042-1051.

[8] Velis D R, Ulrych T J. Simulated annealing wavelet estimation via fourth-order cumulant matching[J]. Geophysics, 1996, 61(6): 1939-1948.

[9] Nikias C L, Petropulu A P. Higher-order spectral analysis: A nonlinear signal processing framework[M]. New Jersey: Prentice-Hall, 1993: 48-76.

[10] 李夕兵, 张义平, 左宇军, 等. 岩石爆破振动信号的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.

[11] 李宏伟, 程乾生. 高阶统计量与随机信号分析[M]. 武汉: 中国地质大学出版社, 2002: 1-27.
LI Hong-wei, CHENG Qian-sheng. Analysis of higher-order statistics and random signals[M]. Wuhan: China University of Geosciences Press, 2002: 1-27.

[12] 中国生, 徐国元, 江文武. 基于小波变换的爆破地震信号去噪的应用[J]. 中南大学学报: 自然科学版, 2006, 37(1): 155-159.
ZHONG Guo-sheng, XU Guo-yuan, JIANG Wen-wu. Application of de-noising in blasting seismic signals based on wavelet transform[J]. Journal of Central South University: Science and Technology 2006, 37(1): 155-159.

[13] Gardner W A, Spooner C M. The cumulant theory of cyclostationary time-series. Part Ⅰ: Foundation[J]. IEEE Trans Signal Processing, 1994, 42(12): 3387-3408.

[14] Spooner C M, Gardner W A. The cumulant theory of stationary time-series. Part Ⅱ: Development and applications[J]. IEEE Trans Signal Processing, 1994, 42(12): 3409-3429.

[15] 王书明, 王家映. 高阶统计量在大地电磁测深数据处理中的应用研究[J]. 地球物理学报, 2004, 47(4): 928-934.
WANG Shu-ming, WANG Jia-ying. Application of higher-order statistics in magnetotelluric data processing[J]. Chinese J Geophys, 2004, 47(4): 928-934.

收稿日期：2009-08-13；修回日期：2009-10-28

基金项目：国家高技术研究发展计划(“863”计划)项目(2006AA06Z105)；湖南省“十一五”重点建设学科资助项目(2006年)

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

(编辑 刘华森)