基于构造系统函数的大地电磁时间序列模拟方法

张宝强　裴建新*②　王　启

(①中国海洋大学海洋地球科学学院，山东青岛 266100； ②海底科学与探测技术教育部重点实验室，山东青岛 266100； ③广州海洋地质调查局，广东广州 510760)

1　引言

大地电磁测深(MT)法是一种在地质探测和矿产普查中用途非常广泛的地球物理方法，通过分析视电阻率曲线和相位曲线的特征获取地下结构电性分布。MT观测数据的阻抗计算以及对阻抗特性产生严重影响的噪声处理是研究的重要内容之一[1-4]。近些年，国内外学者在MT阻抗估算和噪声处理研究上做出了很多努力，频率域模拟电磁场数据已用于验证和评估MT阻抗估计方法[5-7]和提取被污染信号的预处理技术[8,9]。与此同时，时间域数据处理方法的发展[10-13]也推动了MT数据处理方法的进步。如今，国内外学者使用不同的方法模拟MT时间序列； Larsen等[6]使用实测磁场时间序列、通过失真函数修改估计的一维传递函数，在测试区域中生成电场时间序列； Chen[14]在测试经验模态分解算法时，采用的模拟数据是由实测磁场和一维地球模型的阻抗合成的频率域电场数据。鉴于实测MT数据均含有噪声，且通常地下电性结构是未知的，所以实测数据在测试时间域阻抗估计方法和信号处理技术可靠性方面具有一定的局限性。Lamarque[9]和Goubau等[5]分别使用简单的高斯随机序列和均匀分布在(-1,1)的随机序列作为无噪声的MT时间序列及磁场分量的实部和虚部。在此基础上，有学者利用二维阻抗张量生成了电场数据； Yee[7]和Filipo等[8]利用数字递归滤波器设计了随频率增大而衰减的指数函数描绘天然磁场的主要频谱特征，将滤波器系数和伪随机序列进行卷积计算获得磁场时间序列，并用模拟的磁场时间序列和模型阻抗的脉冲响应进行卷积计算获得无噪的电场时间序列；Loddo等[15]提出了一种合成宽频带电磁场时间序列的方法，研究了场的极化特征，其思想是将构造的系统函数作为频率域磁场，建立层状模型阻抗的频率域积分解析式，使用傅里叶积分变换算法获得磁场和阻抗的时间域脉冲响应，磁场脉冲响应与一个随机数序列做褶积运算生成合成磁场时间序列，然后与阻抗脉冲响应做褶积运算获得合成电场时间序列； 俄罗斯学者Varentsov等[16]使用合成数据库(COMDAT Project)比较不同的MT数据处理方法，其中，由选定的随机函数乘以简单的场得到磁场频谱，并结合任意一维标量阻抗计算电场频谱； 严家斌[17]、蔡剑华等[18]、陈知富等[19]在研究中使用随机数作为MT各分量的时间序列；凌振宝等[20]将多个单频率的正弦波线性叠加模拟电磁场数据； 汤井田等[2]使用实测磁场数据和简单二维电阻率模型阻抗计算合成电场数据。

本文提出一种MT时间序列模拟方法，不仅可以避免直接使用实测数据作为合成数据而引入干扰信息的不确定性，还可以改善仅使用随机数序列及正弦波叠加序列作为合成数据而难以体现MT基本频谱特性的问题，该方法支持不同的采样率，得到的模拟数据可用于进一步的数据处理方法研究。

2　方法原理

Loddo等[15]用构造的系统函数作为磁场单边谱刻画MT频谱。在此基础上，本文提出了基于构造系统函数、通过模型阻抗实现包含地下地电结构信息的电磁场数据的模拟方法。

2.1　电磁场频谱的求取

在确定模拟时间序列的采样频率Fs和序列长度N之后，输入系统函数；为使模拟的电磁场频谱和后续模拟的时间序列中带有一定限度的波动和跳跃、且随机数的DFT序列的量级在整个频段保持在1附近而不改变系统函数刻画的磁场频谱的量级[7,8,15,16,21,22]，需引入一均值为0、方差为1、长度为N的随机序列，对其进行归一化FFT计算[22]，然后在频率域与系统函数乘积，该计算结果就作为模拟的频率域磁场。

为了使模拟数据包含地下的地电结构信息，建立一维水平层状电阻率模型，使用正频率序列计算模型阻抗，利用频率域电磁场和阻抗关系式计算模拟的频率域电场。

2.2　电磁场时间序列的计算

傅里叶分析建立了函数时间域与频率域之间的联系，由模拟的频率域磁场和电场分别做IFFT计算得到磁场和电场时间序列；根据实信号DFT的正负频率复共轭对称性质[21-23]，在引入随机数序列之前对系统函数做复共轭对称处理，扩展到负频率使其成为双边对称谱。复共轭对称处理公式为

Re[f(-ω)]=Re[f(ω)]

Im[f(-ω)]=-Im[f(ω)]

(1)

式中f(ω)代表正频率ω的频域磁场h(ω)或阻抗序列z(ω)。同样地，将正演的阻抗序列做复共轭对称处理扩展为双边对称谱，最后将电场和磁场的双边对称谱利用IFFT算法即可获得模拟的实数信号。

基于以上思路，模拟MT时间序列的基本流程如图1所示。

图1　模拟MT时间序列的基本流程

3　一维磁场数据模拟

3.1　频率域磁场数据模拟

模拟MT数据之前需要分析天然电磁场的基本频谱特性。天然电磁场的产生有多种原因，高于6Hz的电磁场源主要来自于高能电磁现象，如闪电、磁暴等气象活动，低于6Hz的电磁场来源于地球磁场和地球电离层。不同频带的天然电磁场的基本特征已有大量学者做了详细阐述[24-26]。磁场水平分量的平均振幅值通常随着周期的增长而增大，在104～10-4s内可以达到10nT； 此外，在1～6Hz及1k～2kHz之间分别存在着一个极小值，这两个极小值代表了低频和高频天然电磁场场源机理的不同。图2所示为一个通过线性时不变系统模拟的、经平滑处理的水平磁场分量的平均频谱，可描述天然磁场的主要频谱特点。

图2　线性时不变系统模拟的水平磁场频谱[15]

一般地磁测量、大地电磁、音频大地电磁方法仅研究有限频带范围内的电磁场，因此本研究也仅在有限的频带范围内模拟MT时间序列数据。Yee[7]和Filipo等[8]提出，可以直接以Z域中极点和零点的方式编写系统的离散传递函数来近似天然场的频率特性; Oppenheim等[27]认为还有一种替代方案，即通过双线性变换将传递函数从连续域变换到离散域。在模拟MT时间序列时，这些方法都是适用的。通过使用连续域中的极点和零点形式，一个系统函数可以表示为[28]

(2)

式中：s=σ+iω是拉普拉斯复数变量，σ是系统阻尼因子； 复数s01,s02,…,s0m分别代表系统函数的m个零点；sp1,sp2,…,spn分别代表系统函数的n个极点；K是比例系数。

函数的解析性质完全由极点和零点在s复数平面上的分布决定。为了使系统稳定，极点必须在s复平面的左半部分，并且任何落在复轴iω上的极点必须是单极点；此外，当变量s虚部为0(即为实数)时，对应的函数值也是实数，要求极点和零点必须是实数或者复共轭对。单极点和单零点、复共轭对极点和复共轭对零点的分子和分母可以分解为[15]

H(s)=

(3)

式中：i、j、l、k分别表示单零点、单极点、复共轭零点和复共轭极点；K1为比例系数；ζ1和χ1代表每个临界频率的阻尼系数。

对于图2所示的天然磁场的频谱，坡度向下变化(能量降低)可以由式(3)的极点刻画，而向上变化(能量升高)可以由式(3)的零点刻画，这样四个极点值和四个零点值就基本表达了天然磁场频谱的总体变化趋势，单极点和单零点都选择实数，如表1所示。对于复共轭对，选择系数ζ1=1、χ1=1，就相当于二重极点和二重零点都选为实数。为了使系统稳定，所有极点频率都取负值[15]，为了使磁场频谱曲线低频处的渐近线为10nT，选择比例系数K1=10。

表1　模拟磁场使用的极点和零点的频率绝对值和特性

满足以上条件后，式(3)转化为

(4)

式中：zi(i=1,2,3,4)为零点；pi(i=1,2,3,4)为极点。

3.2　时间域磁场数据模拟

归一化的DFT序列做乘积，作为模拟的频率域磁场。这样可以使得磁场曲线能够平稳变化，更接近实际情况。得到频率域磁场数据之后，由IFFT便可获得磁场时间序列。

为了验证本文方法，模拟一段采样频率为2Hz、采样点数为10000的电磁场时间序列数据，其频率范围为0.0002～1Hz，参与计算的正频率点数为5000个，磁场振幅谱由式(4)构造，如图3所示。图4为模拟的磁场频谱和时间序列，时间序列记录长度为5000s，总体上模拟的磁场信号随时间变化平稳，具体数据统计参量：均值为-9.7868×10-4nT，均方差为1.5348×10-5nT，最小值为-9.7×10-3nT，最大值为1.05×10-2nT，能量为0.1630nT2。

图3　系统函数构造的磁场振幅谱

图4　模拟的磁场频谱和时间序列，频率范围为0.0002～1Hz，时间长度为5000s

4　一维电场数据模拟

4.1　一维层状模型阻抗计算

电场可利用关系式E(ω)=Z(ω)·H(ω)计算，所以必须获得地下介质的阻抗信息。本文引入一维水平层状模型，如图5所示，设大地由n层水平层状介质所组成，各层的电阻率为ρ1,ρ2,…,ρn, 厚度为h1,h2,…,hn-1。可利用同一层中两个不同深度处阻抗值之间的关系及界面上阻抗连续条件推导地面阻抗[29]。

图5　n层介质模型

由于层状介质中各层的电阻率ρ和传播系数k不同，任一层介质中的波动方程为

(5)

假定第m层的顶面深度为zm，底面深度为zm+1，则厚度hm=zm+1-zm。第m层的特征阻抗为Z0m=-iωμ/km，于是可得阻抗的递推公式

(6)

需注意，在计算模型阻抗的过程中，参与计算的频率点与计算磁场单边谱的频率点需保持一致，频率范围为(0，Fs/2)，为使模拟的电场频谱为对称谱，需将阻抗扩展到负频率部分使其对称。图6为设计的一个一维水平四层模型，使用相同的频率点计算该模型的大地阻抗的实部和虚部(图7)。

4.2　时间域电场数据模拟

获得模拟磁场数据和模型的阻抗之后，在频率域根据电磁场关系式E(ω)=Z(ω)·H(ω)计算电场的频谱序列，对电场对称频谱做离散反傅里叶变换，就可得到模拟的电场时间序列。下列几点需注意： ①模拟电场的零频率值无法根据关系式得到，本文取模拟的磁场零频率对应的DFT值作为电场零频率的DFT值； ②需考虑采样长度N的情况，尤其是当N为偶数时，奈奎斯特频率处的DFT值必须是实数。图8是由模拟磁场和阻抗计算的模拟电场频谱和时间序列，其时间长度为5000s，频率范围为0.0002～1Hz。总体来看模拟的电场信号也是随时间平稳变化的，如图8d所示，其数据统计参量： 均值为-1.2292V/m，均方差为7.2409×10-5V/m，最小值为-1.2504V/m，最大值为-1.2055V/m，能量为1.5111×104(V/m)2。

图6　一维水平层状模型

图7　模型阻抗的实部和虚部

通过模拟的电磁场时间序列，可进一步计算得到模型的视电阻率和相位曲线，结果如图9所示。

4.3　模拟数据与实测数据对比

为了检验模拟效果，选取一段某海域实测的海洋MT时间序列数据，为TE模式的水平电场x分量和磁场y分量，采样频率为10Hz，时间长度为15min，计算其振幅谱，然后用本文提出的方法模拟电磁场数据，如图10所示。对比实测数据与模拟数据振幅谱可以看到，模拟的电磁场振幅谱基本体现了大地电磁场在低、中频段随着频率的增大，能量逐渐衰减的总体趋势，即信号频谱呈下降趋势，而中频区段之后频谱缓慢爬升的一般特征[29,30]。

图8　模拟的电场频谱和时间序列

图9　模拟数据计算的视电阻率(左)和相位(右)曲线

图10　实测海洋MT数据(左)与模拟MT数据(右)的时间序列和振幅谱

5　二维水平非均匀介质电磁场数据合成

在实现一维MT时间序列模拟的基础上，进一步合成得到二维水平非均匀介质电磁场数据。假设地下介质是非均匀的，并且测量轴与介质的电性主轴重合，则任何一个椭圆偏振波都可以看成是两个线性偏振波的合成，这两个线性偏振波就如同分别在两个不同电阻率的均匀介质中传播。那么对应的场可以分解为分别平行和垂直于构造走向[29]的两个部分。在这种情况下，每个电场分量和与之垂直的磁场分量有关，相应的电磁场阻抗关系为

(7)

式中:ZTE=Zxy;ZTM=Zyx;Zxy≠Zyx。

汤井田等[2]和Loddo等[15]对给定的两个地层模型进行一维正演，得到两组随频率变化的阻抗值，将两个阻抗值分别作为一个矩阵的非对角元素，合成的矩阵即可用来模拟地下二维结构。利用该方法，本文在一维模拟的基础上通过分量合成的方法获得二维水平非均匀介质电磁场数据，建立两个不同的一维水平层状电阻率模型(表2)，将获得的两个阻抗响应分别作为Zxy和Zyx。在模拟两个磁场分量的时候，均使用式(4)构造磁场单边频谱，但是两个分量采用的随机序列是不同的，以此来获得两个不同的磁场分量Hx和Hy，并由式(7)计算得到电场分量Ex和Ey。

表2　两个一维水平层状模型的电性参数

设定采样频率为10Hz，时间长度为1h，采样点数为36000，计算得到该模型的电磁场水平分量如图11所示。在上述水平层状电阻率模型下模拟的水平电场和磁场是变化缓慢的平稳信号。此外，相比磁场，电场变化相对更加平稳，体现了大地电磁信号随时间变化的基本特征，由模拟时间序列计算TE模式和TM模式的模型视电阻率和相位曲线(图12)，可见视电阻率和相位曲线光滑连续，反映了地下模型介质的电阻率和相位特征。由于模拟的每一个采样点的电场和磁场都满足阻抗关系式 ，所以模拟信号不包含任何噪声，这是因为任何包含在电场或磁场中的相关或不相关噪声都不会满足该阻抗关系式。

图11　模拟二维MT水平分量时间序列

图12　模拟数据计算的视电阻率(上)和相位(下)曲线

6　结束语

本文实现了一种模拟MT时间序列的方法，该方法不仅在模拟电磁场时结合了天然电磁场的频谱特性，而且模拟的电磁场数据包含了地下介质的视电阻率和相位信息，可以任意选取采样频率和序列长度，适用于一维电阻率模型和二维水平非均匀介质的电阻率模型。由该方法获得的模拟数据与实测的海洋MT数据比较发现，模拟的电磁场振幅谱基本体现了大地电磁场在低、中频段随着频率的增加，能量衰减的总体趋势。该方法模拟的的数据并未考虑噪声的影响，因此下一步的研究应包含噪声和更复杂的电阻率模型，使其适用范围更广。

[1]张全胜,王家映.大地电磁测深资料的去噪方法.石油地球物理勘探,2004，39(增刊1):17-23.

Zhang Quansheng,Wang Jiaying.An method of noise elimination for magnetotelluric sounding data.OGP,2004，39(S1):17-23.

[2]汤井田,张弛,肖晓等.大地电磁阻抗估计方法对比.中国有色金属学报,2013，23(9):2351-2358.

Tang Jingtian,Zhang Chi,Xiao Xiao et al.Comparison of methods for magnetotelluric impedance estimation.The Chinese Journal of Nonferrous Metals,2013，23(9):2351-2358.

[3]蔡剑华,王先春,胡惟文.基于经验模态分解与小波阈值的MT信号去噪方法.石油地球物理勘探,2013,48(2):303-307.

Cai Jianhua,Wang Xianchun,Hu Weiwen.De-noising of MT signal based on empirical mode decomposition and wavelet threshold method.OGP,2013,48(2):303-307.

[4]张翔,严良俊,苏朱刘等.山区MT资料处理新技术及效果.石油地球物理勘探,2007,42(4):454-456.

Zhang Xiang,Yan Liangjun,Su Zhuliu et al.New technology and effects of mountainous MT data processing.OGP,2007,42(4):454-456.

[5]Goubau W M,Gamble T D,Clarke J.Magnetotelluric data analysis:removal of bias.Geophysics，1978,43(6):1157-1166.

[6]Larsen J C,Mackie R L,Manzella A et al.Robust smooth magnetotelluric transfer functions.Geophysical Journal International,2007,124(3):801-819.

[7]Yee E.The reconstruction of the magnetotelluric impedance tensor:An adaptive parametric time-domain approach.Geophysics,1988,53(3):1080-1087.

[8]Filipo W A S,Hohmann G W.Computer simulation of low-frequency electromagnetic data acquisition.Geophysics,1982,48(9):1219-1232.

[9]Lamarque G.Improvement of MT data processing using stationary and coherence tests.Geophysical Prospecting,1999,47(6):819-840.

[10]Ernst T.A comparison of two methods of the transfer function calculation using the least-square criterion in time and frequency domain.Publications of Institute of Geophysics Polish Academy of Sciences,1981,143:13-24.

[11]Kunetz G.Processing and interpretation of magnetotelluric soundings.Geophysics,1972,37(6):1005-1021.

[12]Mcmechan G A,Barrodale I.Processing electromagnetic data in the time domain.Geophysical Journal International,1985,81(1):277-293.

[13]Spagnolini U.Time-domain estimation of MT impe-dance tensor.Geophysics,1994,59(5):712-721.

[14]Chen J.Using empirical mode decomposition to process magnetotelluric data.Christian-Albrechts-Universität,2013.

[15]Loddo M,Schiavone D,Siniscalchi A.Generation of synthetic wide-band electromagnetic time series.Geophysics,2009,45(2):289-301.

[16]Varentsov I M,Sokolova E Y.Generation of synthetic magnetotelluric data.Izvestiya Physics of the Solid

Earth,1995,30(6):554-562.

[17]严家斌.大地电磁信号处理理论及方法研究[学位论文].湖南长沙：中南大学,2003.

Yan Jiabin.The Study on Theory and Method of Magnetotelluric Signal Processing[D].Central South University,Changsha,Hu’nan,2003.

[18]蔡剑华,李晋.基于频率域小波去噪的大地电磁信号工频干扰处理.地质与勘探,2015,51(2):353-359.

Cai Jianhua,Li Jin.Suppression of power line interfe-rence on MT signals based on the frequency domain wavelet method.Geology and Exploration,2015,51(2):353-359.

[19]陈知富,邓居智,陈辉等.基于小波变换的大地电磁信号去噪研究.工程地球物理学报,2012,9(6):732-737.

Chen Zhifu,Deng Juzhi,Chen Hui et al.De-noising research of magnetotelluric signal based on wavelet transform.Chinese Journal of Engineering Geophy-sics,2012,9(6)：732-737.

[20]凌振宝,王沛元,万云霞等.强人文干扰环境的电磁数据小波去噪方法研究.地球物理学报,2016,59(9):3436-3447.

Ling Zhenbao,Wang Peiyuan,Wan Yunxia et al.A combined wavelet transform algorithm used for de-noising magnetotellurics data in the strong human noise.Chinese Journal of Geophysics,2016,59(9)：3436-3447.

[21]Bracewell R,Kahn P B.The Fourier Transform and Its Applications.McGraw-Hill,2000.

[22]Rao K R,Kim D N,Hwang J J.Fast Fourier Transform - Algorithms and Applications.Springer Netherlands,2010.

[23]程乾生.数字信号处理.北京：北京大学出版社,2010.

Chen Qiansheng.Digital Signal Processing.Peking University Press,Beijing,2010.

[24]Labson V F,Becker A,Morrison H F et al.Geophysical exploration with audiofrequency natural magnetic fields.Geophysics,1985,50(4):656-664.

[25]Serson P H.Instrumentation for induction studies on land.Physics of the Earth & Planetary Interiors,1973,7(3):313-322.

[26]Filloux J H.Techniques and instrumentation for study of natural electromagnetic induction at sea.Physics of the Earth & Planetary Interiors,1973,7(3):323-338.

[27]Oppenheim A V,Schafer R W.Digital Signal Processing.Prentice-Hall,New Jersey，1975.

[28]Seshu S,Balabanian N.Linear Network Analysis.Wiley-Blackwell,New Jersey,1959.

[29]李金铭.地电场与电法勘探.北京：地质出版社,2005.

[30]刘国栋,陈乐寿.大地电磁测深研究.北京：地震出版社,1984.