强干扰环境下的大地电磁时间序列处理过程

张赟昀 王培杰 陈小斌1，* 詹 艳 韩 冰 王立凤 赵国泽

1)应急管理部国家自然灾害防治研究院，北京 100085

2)中国地震局地质研究所，地震动力学国家重点实验室，北京 100029

0 引言

极低频电磁台网不仅可监测稳定的高信噪比人工源极低频电磁场信号，也可监测天然场的电磁信号(赵国泽等，2010)。然而，天然场电磁信号包含各种人文噪声，不同台站受到的噪声干扰情况不同，使得不同台站的视电阻率曲线有的频段数据质量较好，有的频段数据质量较差。数据质量差的频段无法为地震短临前兆提供有用信息，影响了极低频电磁台网数据的使用率，因此有必要研究含噪声干扰情况下的数据处理方法。

大地电磁测深法(Magnetotelluric，MT)(Cagniard，1953)是一种天然源频率域的测深方法，探测深度受限于地下电性结构和观测数据有效周期等因素，有效周期越长，探测深度则越大。天然场源信号弱，存在死频带、近场干扰、随机噪声等影响阻抗估计结果的因素，在强干扰区很难获取十分准确的电磁响应，而稳定可靠的阻抗估计是后续反演解释的前提(王家映等，1998)。为得到无偏的大地电磁响应，地球物理学家提出了一系列方法。Sims等(1971)假定天然源大地电磁信号是平稳信号，磁道在一定区域内为平稳分布，且电、磁道噪声服从高斯分布，可采用最小二乘法进行阻抗估计；Goubau等(1978，1984)和Gamble等(1979a，b)提出了远参考道方法(Remote Reference)，通过同时记录测点与参考点的磁场，用互功率谱代替自功率谱以减少相干噪声的影响；Egbert等(1986)和Chave等(1989，2003)提出稳健估计方法(Robust)，通过降低含噪声窗口的傅里叶系数权重压制随机噪声。在时域去噪方面，汤井田等(2008)使用希尔伯特-黄变换(Hilbert-Huang Transformation，HHT)分析大地电磁信号中的噪声分布特征，并进行干扰压制，取得了一定的效果。Jin等(2012)应用HHT处理海洋电磁数据，并与有界影响估计(Bounded-Influence Remote Reference Processing，BIRRP)(Chaveetal.，2003，2004)的处理结果进行了对比，发现HHT得到的结果更稳定。Wang等(2017)提出了STIN(Synthesis of Natural Electric and Magnetic Time-series Using Inter-station Transfer Functions and Time-series from a Neighboring Site)法，通过站间传输函数合成本地道含噪声时段的数据，避免含噪声段数据对阻抗估计的影响。李晋等(2019)结合变分模态分解(Variational Mode Decomposition，VMD)和匹配追踪技术(Matching Pursuit，MP)压制音频大地电磁强干扰。

由于各种新兴的数据处理方法稳定性不足、适用性不强，并没有得到普遍应用，远参考道方法和稳健估计方法仍然是当前生产中应用最广泛的方法。然而在实际工作中我们发现，使用远参考道法和稳健估计法相结合的方式处理在强干扰环境下采集的数据，效果并不理想。其他学者也针对这个问题展开过研究，例如Egbert(1997)利用多台站数据估算背景噪声水平，使用低相干噪声的数据估算阻抗，当至少存在一个台站不受相干噪声影响或为间歇相干噪声源时，可得到比单台站处理更好的阻抗估计结果，但是这种方法要求同时采集多台站的数据。本文希望找到合适的数据处理方案，以提升单台站、单远参考台站数据处理结果的稳定性和可靠性。中国东部是经济发展繁荣的地区，人文干扰强烈，本文对郯庐断裂带中南段强干扰区的大地电磁数据进行处理，选取了其中部分特征明显的测点进行分析，并针对不同的干扰类型提出了一套数据处理方法。

1 数据处理方法

大地电磁法记录相互垂直的天然源电磁场分量，在频率域内MT电磁场与阻抗张量满足如下的线性关系：

(1)

式(1)中，Ex、Ey、Hx、Hy为4个水平电磁场分量，Zxx、Zxy、Zyx、Zyy为4个阻抗张量元素，可简写为

E=ZH

(2)

阻抗的特征通常用视电阻率ρ和相位φ来表征。式(3)、(4)中，i，j∈{x，y}：

(3)

(4)

ERH=ZΗRH

(5)

其最小二乘解形式如式(6)，式中，上标H为共轭转置符号：

Z=ERH(HRH)-1

(6)

展开为

(7)

式(7)中，星号代表共轭，尖括号代表平均值，即功率谱的平均值。

在大地电磁测深中，实测数据通常包括信号和噪声：

Em=Es+En

Hm=Hs+Hn

Rm=Rs+Rn

(8)

式(8)中，下标m表示实测数据(measured)，s表示信号(signal)，n表示噪声(noise)。实测数据的张量阻抗最小二乘解为

(9)

由于水平电场变化剧烈，阻抗求解不稳定，通常用磁场作为参考道：

(10)

互相关可突出2个信号中相干性强的部分，压制相干性弱的部分。若2个信号中的噪声不相干，则互相关可以压制噪声；若2个信号的噪声相干性强，互相关反而会突出噪声信号。因此，最小二乘法阻抗估计的关键在于选取与Es和Hs相干性强、与En和Hn相干性弱的参考道。当观测点附近存在强人工源噪声信号时，En和Hn具有强相干性，采用最小二乘法得到的阻抗估计会严重失真。Goubau等(1978，1984)和Gamble等(1979a，b)提出远参考法，以远参考点的磁场作为参考道R参与最小二乘计算，可有效压制相干噪声。

稳健估计(Robust)技术是一种通过自适应加权来降低飞点对阻抗估计影响的方法。Robust将最小二乘法的计算结果作为初始阻抗，计算系统观测值与预测值之间的误差，根据误差设置权函数对观测值加权，并利用加权最小二乘法计算新的阻抗。重复上述步骤直至计算的阻抗不再变化。Larsen等(1996)发现当相干噪声较大时，基于Robust技术的单点最小二乘法阻抗估计会给出错误的结果，而基于Robust技术和远参考的阻抗估计结果则接近真实阻抗。但Robust技术的前提是预测值接近真实值，这要求数据中干扰较弱的时间长度远大于干扰强的时间长度。而强人文干扰的地区显然不满足这种条件，Robust处理强干扰持续时间长的数据往往会放大干扰信号而压制有效信号。

虽然近场效应具有明显的特征，但对于近场干扰不太强烈的测点，其曲线具有与近场类似的特征，视电阻率上翘但角度不到45°，相位下降但不到0°，很难与正常的测点曲线区分。若近场影响频段较小，很容易被误认为是高阻异常。Rhoplus(Parkeretal.，1996)是一种直接反演视电阻率和相位的一维反演方法，该方法考虑视电阻率和相位数据的关联性和数据的连续性，可用于大地电磁视电阻率和相位数据的一致性检测，以判断曲线是否合理并矫正畸变数据(谭捍东等，2004；Sprattetal.，2005；周聪等，2015)。本文中所有数据的处理结果都已导入大地电磁资料处理与解释可视化集成系统MT-Pioneer(陈小斌等，2004)中，剔除飞点后与Rhoplus曲线拟合进行验证，并以Rhoplus曲线为标准反复进行数据编辑。

由于强人文干扰地区的相干噪声存在的时间长，对全时段采集的数据做Robust处理难以得到准确的阻抗，因此有必要对数据进行分段处理，从中挑选出受干扰较小时段的数据参与叠加。若分段数较少，很可能每段中都包含强噪声，不能挑选出不含噪声时段的数据；若分段数太多，则低噪声时段的数据有可能被截断，也不能获取较好的阻抗估计结果。本文尝试采用不同的分段数进行处理，分析分段数对处理结果的影响。利用Robust方法处理强干扰数据的效果并不稳定，本文分析了Robust方法在强干扰区数据处理中的应用效果。在强干扰区数据质量差的时段占比很大，需要在进行功率谱叠加前予以剔除，本文中所有功率谱都在叠加前经过人工挑选。

2 实测数据采集与处理

2.1 工区概况及数据采集

本文选取中国东部郯庐断裂带(合肥—宿迁段)的实测数据进行处理分析。测区位于安徽省东北部，该区大部分为平原，经济发展较好。区内村庄城镇密布，电网纵横，有风力发电站、太阳能发电站、电气化铁路以及大型工厂，人文干扰强烈。以往采集的数据表明，在该区域很难采集到质量好的数据。特别是在0.1Hz附近的死频带内，天然场源信号弱，近场信号强，导致在该频带内难以得到稳定的阻抗数据。

图1 大地电磁测点及远参考点分布图Fig.1 Distribution of MT station and remote reference station.

2.2 远参考处理

Shalivahan等(2002)分别使用3个远参考点对同一测点进行了处理，发现当远参考点距离观测点较近时(80km、115km)，远参考技术难以改善死频带的处理结果，当远参考点距离观测点足够远时(215km)，才能在死频带得到稳定平滑的阻抗估计。为改善数据质量，本文设置了远参考点，采用远参考技术对工区内的数据进行处理，并根据远参考技术的处理效果对远参考点的位置进行了调整。最初的远参考点(远参一)位于工区南面安徽省与江西省交界处的山区内(29°55′46″N，117°23′38″E)(图1)，距离工区内测点的平均距离约为350km。收取第1轮远参考数据后，使用该远参考点对工区内的测点进行远参考处理，发现该点自身的数据质量不佳，存在不连续的干扰信号，从图2 中可以看出该远参考点的中低频段与Rhoplus曲线不能很好地拟合。810测点参考后数据质量稍有改善，高频部分与Rhoplus曲线拟合得更好，中频近场形态明显被压制，但中低频部分仍有飞点存在，拟合效果不佳，且与远参考点拟合不好的频带一致。随后，我们在山东省临沂市蒙阴县(35°31′57″N，118°1′8″E)设置了第2个远参考点(远参二)，该远参考点处于工区北部，距离工区内测点的平均距离约为260km。远参二的低频部分非常连续，但0.1～10s处的视电阻率快速上升又急剧下降，虽然曲线平滑，但与Rhoplus曲线拟合不上，相位曲线更是在1～10s处下掉，因此推测该频带可能存在稳定干扰。使用该点对工区内测点进行远参考处理，发现工区内的数据有所改善，可能的原因是高频电磁信号衰减快，与测区干扰相干度低，或远参考站的干扰频带与测点的干扰频带不一致。但仍有部分测点在参考后出现异常，如图2 中305测点参考后，在0.1～10s范围内视电阻率能完全拟合，但相位拟合不好，推测参考站的本地噪声也会影响远参考处理的效果。远参三设置在内蒙古鄂托克前旗(107°49′40″N，38°1′27″E)，电磁环境非常好，几乎没有人文噪声，且距离工区的平均距离约为1050km，与工区之间存在相干噪声的可能性很低，远参考处理能很好地改善工区的数据质量。图2 中远参三的视电阻率和相位均与Rhoplus曲线拟合得很好，仅在低频接近10000s处的yx向视电阻率相位未拟合好。509测点未进行远参考处理的数据在10s处不连续，处理后数据十分连续，与Rhoplus曲线的拟合情况非常好，且和远参考点类似，仅在接近10000s处的yx向视电阻率、相位拟合不佳。

图2 不同远参考站的数据处理效果图Fig.2 Data processing effect of different remote reference stations.红点、蓝点分别为xy、yx向视电阻率和阻抗相位，橙线、绿线分别为xy、yx向Rhoplus反演曲线，下同

由于 3 个远参考站没有同时记录，仅挑选具有代表性的测点进行比较。不难发现，所有测点的数据经过远参考处理后均有改善。一般情况下，只要远参考站本身的数据干扰较小，经过远参考处理后的数据都会变好。因此，只要条件允许，都应该布设远参考站。下文中的处理过程均是在远参考处理的基础上进行的。为简化叙述，不再对此作单独说明。

2.3 Robust处理

工区范围较大，测点间距大，不同测点受到干扰的情况不同。为测试不同干扰情况下Robust技术的处理效果，将测点使用Robust技术处理的结果与不使用Robust(简称非Robust)技术处理的结果进行对比。

图3 709测点的视电阻率和相位曲线Fig.3 Apparent resistivity and phase curves of MT station 709.

在0.1Hz附近的死频带内，若没有近场干扰，信号整体比较弱，天然场的信号强度和环境噪声的信号强度处于同一数量级。使用非Robust技术处理时，对每段数据的谱以相同的权重进行叠加，微弱的天然场信号会淹没在噪声中。如果使用Robust技术处理，对相干性较强的谱以较大的权重进行叠加，能够突出具有一定相干度的天然场信号，从而获得更准确的阻抗估计。以709测点为例(图3)，从未经编辑的曲线中可以看出，709测点的视电阻率和阻抗相位在高频部分和低频部分均展示出基本形态。对比Robust处理和非Robust处理的未经编辑的曲线可以发现，经Robust处理的曲线形态更好，0.1Hz附近死频带内的频点更接近Rhoplus曲线，而非Robust处理的曲线中飞点现象严重。对二者的功率谱进行精细挑选后得到最终结果，可以发现Robust的处理结果和Rhoplus曲线能很好地拟合，而非Robust的处理结果在死频带仍存在飞点，不能与Rhoplus曲线拟合。

图4 209测点的视电阻率和相位曲线Fig.4 Apparent resistivity and phase curves of MT station 209.

而对于持续存在近场干扰源的点，干扰信号的电场和磁场相干度更高，Robust会将近场信号识别为期望的信号，处理过程中会压制天然场信号，放大近场干扰信号，使阻抗估计偏向于近场信号的阻抗。工区内靠近大型城市、大型工厂、电气化铁路等近场源的测点曲线都呈现视电阻率45°上升、相位趋于0°的近场特征。以209测点为例(图4)，该点SW向1.5km处有村庄，W向1.5km处有工厂，存在潜在的干扰源，使用Robust处理后得到的曲线存在明显的近场特征，视电阻率曲线在1Hz附近接近45°上升，yx方向的相位曲线接近0°。由于分段进行Robust处理后每段的视电阻率和相位均呈现近场形态，数据不具备可编辑的空间，编辑后近场形态改善不明显，且中频仍存在飞点，因此对该点采用Robust处理只能得到近似于近场的阻抗。采用非Robust方法进行处理，即分组功率谱叠加时不对每段数据赋予权重，可得到少量无近场特征或近场特征弱的谱，经过人工编辑，滤除以近场为主的功率谱，仅保留少量符合大地电磁曲线特征的谱，叠加后可以发现视电阻率和相位都能与Rhoplus曲线较好地拟合，近场形态改善明显。

工区内数据处理的难点在于死频带和近场干扰这2种情形，同时也不乏二者同时存在的测点，特别是近场频带和死频带有重合部分的测点，其数据处理十分困难。这种情况下的数据中包含近场信号、有效信号和随机噪声，其中近场信号最强，有效信号和随机噪声较弱。采用Robust处理会放大相干度高的近场信号和有效信号，压制相干度低的随机噪声，但近场信号较强，因此阻抗估计结果仍以近场阻抗为主；不采用Robust处理则会让死频带内的天然场信号和随机噪声淹没在近场信号中，仍然无法得到准确的阻抗估计。针对这种数据，我们把2种方法处理得到的谱放在一起编辑，兼顾二者的优点。以测点307为例(图5)，分别用2种方法进行处理后，均不能得到理想的阻抗估计，最终编辑得到的曲线中均存在飞点，而将2种处理方法得到的功率谱放在一起进行编辑后，曲线在0.1Hz处有明显改善，且低频曲线更为平滑。其可能的原因有2点：1)参与挑选和叠加的功率谱增加了一倍，组合得到合理结果的概率增大，对结果的改善存在偶然性；2)不同时间段的近场信号、有效信号和随机噪声强弱存在变化，Robust和非Robust处理得到的谱都参与编辑，近场干扰强时选择非Robust处理得到的谱可以降低近场的影响，而在近场干扰弱时，选择Robust处理得到的谱有利于从背景噪声中挑出有效信号。值得注意的是，307测点只是个例，该方法的通用性未经考证，仅供参考。

图5 307测点稳健估计(Robust)处理、非稳健估计(非Robust)处理及二者叠加处理得到的视电阻率和相位曲线Fig.5 Apparent resistivity and phase curve of MT station 307 obtained by robust processing，non-robust processing and superposition processing.

2.4 分段数

除采用不同的Robust处理策略外，还可通过增加分段数来改善数据质量。死频带和近场频带都位于中频带，中频带的数据具有足够长的采样时间，增大分段数会减小分段处理中每段的时间，得到更多组功率谱。当某时间段内的干扰较小或近场干扰暂停时，该段数据可以得到相对准确的阻抗估计，人工编辑时将其保留。相对于每段时间长、段数少的处理方式，这种每段时间短、段数多的处理方式能保留更多的有效功率谱。以503测点和504测点为例(图6，7)，截取4000min的数据，分别采用5段、10段、20段、40段的分段数，每段对应的时长分别为800min、400min、200min、100min，经过精细编辑后得到如图6、7所示的响应曲线。不难发现，随着分段数的逐渐增加，响应曲线与Rhoplus曲线的拟合越来越好。增大分段数无疑会增加人工编辑的工作量，但在电磁环境复杂的情况下，为改善数据质量，增加工作量是值得的。

图6 503测点不同分段数对应的视电阻率和相位曲线Fig.6 Apparent resistivity and phase curves of MT station 503 corresponding to different segment numbers.从左至右对应的每段时长分别为800min、400min、200min和100min

图7 504测点不同分段数对应的视电阻率和相位曲线Fig.7 Apparent resistivity and phase curves of MT station 504 corresponding to different segment numbers.从左至右对应的每段时长分别为800min、400min、200min和100min

图9 305测点昼夜不同时长的视电阻率及相位曲线对比图Fig.9 Comparison of apparent resistivity and phase curves of MT station 305 with different time and length.

2.5 分时段处理

人文干扰的强度以天为周期，白天人类活动多时干扰强，夜晚人类活动少时干扰弱。使用305测点进行昼夜数据的对比研究，分别选取10：00—16：00和23：00—5：00作为昼夜数据对比时段。图8 对比了白天和夜晚未编辑与编辑后的数据，从图中可明显看出夜晚数据的质量优于白天。白天的死频带数据经编辑后仍存在大量飞点，与Rhoplus曲线拟合较差。而夜晚的死频带数据情形稍好，低频部分能得到连续光滑的曲线。为了突出夜晚数据质量较高的特征，选取了10h白天的数据，与5h夜晚的数据进行对比，从图9 中未经过人工选谱的原始数据可以看出，虽然白天的采集时间是晚上的2倍，但视电阻率曲线不连续、相位曲线在10～0.1Hz杂乱无章。经人工选谱后白天和夜晚的数据均得到改善，但仍是夜晚的数据比白天好。

对305测点白天(8：00—20：00)和夜晚(20：00—8：00)的数据分别截取12h、24h、36h、48h、60h进行处理，其中第1段为第1晚的数据，后面每段为前面一段加上与之连续的白天或夜晚12h的数据。从图10 中可以看出，随着采集时间的增加，视电阻率曲线的质量逐渐变好，60h数据的视电阻率曲线与Rhoplus曲线拟合得很好，相位曲线与Rhoplus曲线的拟合尚可。从左往右观察相位曲线可以发现，相位曲线的质量并不是逐步改善的。从图10 第2张和第4张小图中可以看出，加入1个白天的数据后，相位反而变差了。由此可知，由于白天的数据干扰较强，对曲线的整体改善有限，因此在数据编辑时可以有所取舍，若白天数据质量较差，可在处理时截去整段数据中白天时段的数据。同样地，建议在进行野外数据采集时，应尽量多采集夜晚的数据，即在下午布设仪器，根据周围的干扰情况确定采集天数后，在上午结束采集，这样可保证夜晚的数据必然存在且比白天多。

图10 305测点不同时长数据对应的视电阻率及相位曲线图Fig.10 Apparent resistivity and phase curves of MT station 309 with different duration.

2.6 采集时长

为了探究采集时间对数据的影响，我们首先使用合成数据研究采集时长与有效周期的关系。选取一个电阻率为 100Ω·m的均匀半空间模型，假设电场信号E(ω)稳定为1mV/km，根据H(ω)=E(ω)/Z(ω)=1/Z(ω)，设计不同的采样时长，用傅里叶逆变换生成电场和磁场的时间序列，然后对时间序列进行处理得到视电阻率和阻抗相位。如图11 所示，由于频谱泄漏和栅栏效应的影响，时长为8s的数据经傅里叶变换后其电磁场的振幅谱并不能完全准确地还原至所能处理的最小频率0.125Hz处。实际上，振幅谱在约2Hz处开始偏离理论值，在约1Hz处达到最低值，相应的阻抗在2Hz处开始发生偏移，1Hz后的阻抗已远远偏离理论值。采用不同长度的合成数据进行计算都能得到一致的结果。因此我们认为，大地电磁数据理论上的最大有效周期是数据时长的1/8。

图11 采样率为4096Hz、时长为8s的合成数据的处理结果Fig.11 Synthetic data processing results with a sampling rate of 4096Hz and a length of 8s.a 振幅谱；b 视电阻率；c 相位

为了验证以上的理论推测，我们选取极低频电磁台网中丰宁台的数据分析。该台建有一个连续观测的大地电磁观测站，可以采集长时间段的大地电磁数据。同时，该观测站的仪器布设相比较常规的野外仪器布设，其电极和磁探头埋设更深，数据采集更稳定。图12a 中24h的数据表明该台站的数据非常稳定，仅通过一个谱就可以得到连续的阻抗。截取不同时长的数据进行处理，得到有效周期与最短数据时长的对应关系如图12b 所示，最短数据时长约为周期的8～10倍。

图12 丰宁台实测数据的处理结果Fig.12 Measured data processing result of Fengning station.a 24h数据处理得到的视电阻率和相位曲线；b 获取有效阻抗所需的最短数据时长

从图10 可见，305测点的中频存在干扰，随着数据时长的增加，数据逐渐改善。有效谱的个数与频率成正比，在高频部分可得到的有效谱个数远多于低频部分。随着观测时长的增加，高频部分的视电阻率相较于低频部分改善得更快。因此，对于中高频干扰，延长观测时间即可改善数据质量。

对比丰宁台和305测点的数据可以发现，不同噪声条件下的数据质量相差很大。对于低噪声干扰地区只需记录10倍有效周期时长的数据，得到一个有效谱即可；而高噪声干扰地区可能需要记录100倍有效周期甚至更长时间的数据才能得到可靠结果。因此，在野外观测时，应根据实际情况合理调整观测时长，才能在满足数据质量要求的前提下以较高的效率完成数据采集。

3 结论

本文结合实测数据分析了数据处理时存在的问题，并针对不同问题提出了数据处理技术方案：

(1)设置远离观测点的、低噪声干扰的远参考点是处理强干扰环境下大地电磁数据的必要条件。对于死频带内信号弱的测点，建议采用Robust技术进行处理；对于近场干扰强且持续时间长的测点，建议采用非Robust技术进行处理；对于2种情况都存在的数据，可以尝试同时使用Robust和非Robust方法得到的功率谱参与挑选。功率谱的人工挑选对于结果的影响十分重要，以上3种处理方法得到的功率谱都必须经过人工精细挑选后叠加，才能得到最终的阻抗估计。

(2)延长采集时间增加有效谱个数，同时适当增加分段数，减少编辑前的功率谱叠加，并对大量功率谱进行人工挑选可有效改善阻抗估计结果。干扰信号以天为周期，白天干扰强，夜晚干扰弱，在数据采集时应尽量保证夜晚数据采集时长，同时在数据处理时有所取舍，可以舍去部分白天的时间序列，也可以在功率谱挑选时去除白天含干扰的功率谱。

(3)理想情况下，最短数据采集时长为最长有效周期的8倍。因此，低噪声干扰区记录10倍有效周期时长的数据就可以获得比较好的阻抗估计结果，而高噪声干扰地区可能需要记录100倍有效周期甚至更长时间的数据才能得到可靠的结果。在野外观测需要根据实际噪声环境合理设计观测时间。

本文的研究区内存在高压线、发电站、大型工厂等诸多干扰，但并未单独就某种干扰进行具体分析。本文所述的Robust、分段数、分时段处理等方法适用于干扰持续存在但部分时段干扰较弱的情况。总而言之，数据受干扰的情况不尽相同，数据处理的技术手段也不应一成不变。在更好的数据处理技术尚未成熟时，应灵活应用已有的技术对大地电磁数据进行处理。