0 引言

大地电磁法(Magnetotelluric，MT)是一类重要的地球物理方法^[1-2]。随着经济的发展和应用环境的多样化，MT方法所面临的噪声问题越来越复杂。如何有效压制各种电磁噪声，获得高质量的阻抗估计，一直是MT方法最重要的任务之一^{[1, 3-5]}。

远参考法是大地电磁法压制噪声最有效的策略之一。Goubau等^[6]首先提出互功率谱估计法，以压制常规单站处理阻抗估计中自功率谱计算所引起的噪声。在此基础上，Gamble等^[7]提出了远参考估计法，该方法要求在距离测点一定范围内观测磁场信号的变化，利用远参考站与测站间信号相干、噪声不相干的性质，以互功率谱代替自功率谱计算，以数据相干性压制噪声的影响。随着人文噪声的影响日益严重，远参考法得到了极大的发展并逐步普及^[8-12]。

远参考技术的发展要求参考数据规模不断地增大；数据的利用从双磁道向电磁多测道、多站叠加及多站处理不断发展。常规参考道数据一般选择两道水平磁场，这是因为磁场传播过程相对更稳定，不同测站间磁场的相关性更好。但当参考站的磁场数据质量不佳时，电场数据也可被用于作为参考道进行互功率谱的计算，即所谓“电参考”^[13]。Epishkin^[14]提出一种利用所有测道的远参考技术，可同时利用电场和磁场分别计算一套权重系数，提高了参考道数据的利用率，进而提高了MT数据处理质量。但该算法仅限于单参考站的处理。Varentsov等^[15]提出利用远参考站的磁场与测站磁场间的相干性分析，计算一系列阈值参数，据此对测站的含噪功率谱数据进行挑选，从而提高阻抗估计的稳健性，即所谓“远参考磁场控制法(RRMC)”。张刚等^[16]将这一思想进一步发展到多点远参考数据，但其实质仍是不同参考站处理结果的叠加，因而对于持续性噪声，仍难以通过功率谱数据的挑选提高数据质量。MT多站叠加技术、“伪远参”及站间依赖关系等方法^[17-19]拓展了远参考法的应用，但这些方法未明确多参考站的应用方式。Larsen等^[20]提出基于远参考站的信噪分量方法，Oettinger等^[21]在此基础上进一步提出一种可利用双远参考数据压制相干噪声的技术，并可同时获得相干噪声的阻抗信息。该方法为双远参考站数据的合理利用提供了可行策略，但当远参考站超过两个时，该方法无法进一步提高数据利用率。

进一步提高参考数据的利用率是远参考技术的发展方向之一。随着阵列观测技术的不断发展^[22]，多参考站数据的获取已具备实际可操作性。多参考站数据不仅可以从测区外直接进行选点观测，还可以从测区内的低噪数据中选择。例如，从测区内低噪数据中选择参考站数据作互参考处理^[23]就是常用的一种选择。一些学者研究了多个远参考站条件下不同参考距离对处理结果的影响^[24-26]，然而算法的实质仍是各参考站数据相互独立地进行处理，未能充分利用多点远参考数据间的相干信息。显而易见，多点远参考站间均可能满足信号相干、噪声不相干的条件。如何更加有效地利用多参考站数据，统一提取多站、多道间的相干信息，增加数据利用率，提高处理质量，仍是重要的研究方向。

阵列电磁数据处理方法的提出与发展为多参考站数据的利用提供了基础^[27-31]。常规阵列处理仅构建单一数据阵列^[28]，天然场源极化参数的提取精度易受影响。本文在常规阵列电磁数据处理方法的基础上，发展了多参考站条件下的阵列电磁数据处理技术，提出分别构建天然场占优的数据阵列和含噪的目标数据阵列，以阵列数据分析方法提取高精度天然电磁场极化参数，继而提高目标数据阵列的天然场响应估计质量。实测数据试验验证了该方法的有效性。

1 阵列电磁数据处理模型

对于频域电磁场，假设在I个激励时窗内，系统输入端共有L个电磁场源同时进行激励，输出端共有K个测道同步测量系统对所有输入激励的总响应。第l(l=1, 2, …, L)个场源在第i(i=1, 2, …, I)个时窗内的极化参数记为A_li，L×I个极化参数构成输入阵列A；第k(k=1, 2, …, K)个测道在第i个时窗的观测数据记为X_ki，K×I个观测数据构成输出阵列X。线性时不变系统条件下，系统频率特性响应U_kl不随时间变化，输出与输入满足线性关系^{[27, 30]}

$ {X_{ki}} = \sum\limits_{l = 1}^L {{U_{kl}}{A_{li}} + {R_{ki}}} $

(1)

矩阵形式为

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{UA}} + \mathit{\boldsymbol{R}} $

(2)

式中：R为频域输出端噪声矩阵，可存在于各个测道与时窗中；A为场源的极化信息，或称为极化参数矩阵；K×L维矩阵U反映了地球电磁系统本身的特性，以及激励源与道的空间信息，与激励延时(观测时窗)无关，称为系统响应矩阵。式(1)和式(2)即为阵列电磁勘探模型的频率域输出—输入关系表达式。

式(1)和式(2)亦是大地电磁法阵列数据处理的基本模型。通过合理的实施方案，可获得阵列观测数据。求解上述阵列方程组，并从中提取相应的地电解释参数，可以建立阵列数据处理的基本理论。

2 多参考站阵列数据处理技术

由于天然场源及噪声的随机性和不可预测性，很难准确地知道场源个数，因此采用数据降维类算法(如主成分分析等)求解式(2)是合理的。一般而言，在大地电磁法和音频大地电磁法涉及的观测频段范围内，天然场源信号可被认为是均匀平面电磁波，因其发收距远大于观测尺度，不同方向场源的主要能量可以通过分解、合并等效为两个正交水平场源的作用，或者说天然场源的主要极化可等效至两个正交水平方向上。

根据上述特点，Egbert^[27]设计了基于稳健主成分分析的常规阵列电磁数据处理方案。该方案将远参考与目标测站置于同一阵列，在低噪环境下，可获得较好的处理结果。然而，当环境中存在相干噪声时，含噪数据会严重影响阻抗数据的估计质量。

多点远参考数据为提高天然场源极化参数的估计质量提供了可能。利用多点远参考数据，将远参考测站与目标测站物理分离，分别构建两类阵列数据集，即天然场信号占主导的低噪参考阵列数据集和含噪的目标阵列数据集。先提取高精度的天然场源极化参数，再利用该参数提高目标阵列数据集的处理质量。在此思路下，本文提出一种分步求解方案，即先估计场源极化参数，后求系统响应参数。

2.1 优选参考站数据，分别构建参考数据阵列和目标测站阵列

一般可通过在测区外布设两个甚至多个远参考站获得多点远参考数据。当测区范围较大、噪声分布不均匀时，还可优选测区内的同步低噪测站数据，或在低噪时段中加入参考数据阵列。参考站数据的选择与常规单站远参考、互参考数据的挑选方法相同，以阻抗估计误差、信号相干度和极化方向等多种参数进行信、噪识别，挑选低噪测站和时段，需保证参考站本身的数据品质较高，并且尽可能与待处理的测站噪声不相干。

实际中，根据参考测站的测道组合方式，主要有以下几种参考阵列组合：①经过挑选的磁道组合；②经过挑选的电道组合；③经过挑选的混合测道组合，包括电道、磁道甚至其他可能的梯度测道等。记参考数据矩阵为X^r，对应上述三种组合方式可分别写为

$ {\mathit{\boldsymbol{X}}^{\rm{r}}} = {[H_x^{(1)}\;\;H\;_y^{(1)}\;\; \cdots \;\;H_x^{(J)}\;\;H_y^{(J)}]^{\rm{T}}} $

(3)

$ {\mathit{\boldsymbol{X}}^{\rm{r}}} = {[E_x^{(1)}\;\;E\;_y^{(1)}\;\; \cdots \;\;E_x^{(J)}\;\;E_y^{(J)}]^{\rm{T}}} $

(4)

$ \begin{array}{l} {\mathit{\boldsymbol{X}}^{\rm{r}}} = [E_x^{(1)}\;\;\;E\;_y^{(1)}\;\;H_x^{(1)}\;\;H_y^{(1)}\;\;E\;_x^{(2)}\\ \;\;\;\;\;\;\;\;E\;_y^{(2)}\;\;H_x^{(2)}\;\;H\;_y^{(2)}\;\; \ldots \\ \;\;\;\;\;\;\;\;E\;_x^{(J)}\;\;E_y^{(J)}\;\;H\;_x^{(J)}\;\;\;H_y^{(J)}{]^{\rm{T}}} \end{array} $

(5)

式中：电磁场分量的上角标1、2、…、J表示经过挑选的参考测站序号，J为参考测站总数；E表示电场；H表示磁场。

同样地，可以将待处理目标测站的张量电道、磁道数据写入同一个数据矩阵中，记为Xⁿ。

(2) 利用参考数据阵列，建立阵列方程组

$ {\mathit{\boldsymbol{X}}^{\rm{r}}} = {\mathit{\boldsymbol{U}}^{\rm{r}}}\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{R}}^{\rm{r}}} $

(6)

式中R^r为参考阵列的信号残差。

(3) 根据多个目标测站的数据可建立阵列方程组

$ {\mathit{\boldsymbol{X}}^{\rm{n}}} = {\mathit{\boldsymbol{U}}^{\rm{n}}}\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{R}}^{\rm{n}}} $

(7)

式中Rⁿ为目标阵列的信号残差，高噪条件下，该残差的方差较大。

2.2 估计天然场源极化参数

通过合适的数据选择，可确保参考数据阵列中天然场占主导。如前文所述，天然场源的主要极化可等效至两个正交水平方向上，即L≈2。此时遵循Egbert^[27]和Smirnov等^[29]提出的稳健主成分分析(Robust Principal Component Analysis, RPCA)方法，可提取出数据中的主成分信息

$ [{\mathit{\boldsymbol{U}}^{\rm{r}}}, \mathit{\boldsymbol{\widetilde A}}] = {\rm{RPCA}}({\mathit{\boldsymbol{X}}^{\rm{r}}}) $

(8)

式中：L×I维矩阵$\mathit{\boldsymbol{\widetilde A}} $为天然场源的等效极化参数；RPCA(·)表示获取矩阵主成分的稳健主成分分析函数。

稳健估计式(7)的实现可直接采用现有的RPCA算法，如公开的Libra统计学代码^[31]。但经测试，该算法在数据量较大时，对计算资源要求较高，计算速度慢。因此，本文采用常规阵列处理中所用的RPCA算法^{[27, 29]}，其简要流程如下。

(1) 计算K个测道的不相干噪声方差，构建K×K维对角矩阵C_K。

(2) 阵列观测数据归一化，得到归一化矩阵${\mathit{\boldsymbol{\widetilde X}}^r} = \mathit{\boldsymbol{C}}_K^{ - \frac{1}{2}}{\mathit{\boldsymbol{X}}^{\rm{r}}} $。

(3) 对归一化矩阵作奇异值分解(SVD)，$ {\mathit{\boldsymbol{\widetilde X}}^{\rm{r}}} = {\mathit{\boldsymbol{u}}^{(0)}}{\mathit{\boldsymbol{s}}^{(0)}}{[{\mathit{\boldsymbol{v}}^{(0)}}]^{\rm{T}}}$，其中特征值矩阵s⁽⁰⁾的初始维数取K′=min(K, I); u、v分别为SVD后的左、右奇异矩阵。取I×I维初始时窗权值矩阵W⁽⁰⁾=I(I为单位对角矩阵)。

(4) 设置最大循环次数j^max、迭代终止条件阈值ε，对j(j=0, 1, …, j^max)作循环计算

$ \begin{array}{l} {\mathit{\boldsymbol{Y}}^{(i)}} = {[Y_1^{(i)}, \cdots , Y_I^{(i)}\left] = \right[{\mathit{\boldsymbol{s}}^{(j)}}]^{ - 1}}{[{\mathit{\boldsymbol{u}}^{(j)}}]^*}^{\rm{r}}{\mathit{\boldsymbol{\widetilde X}}^r}\\ r_i^{(j)} = Y_i^{(j)}\frac{I}{{K\prime }}\\ w_i^{(j)} = \left\{ \begin{array}{l} 1\;\;\;\;\;\;\;\;\;\;\;\;\;r_i^{(j)}{r_0}\\ \frac{{{r_0}}}{{r_i^{(j)}}}\;\;\;\;\;\;\;r_i^{(j)} \ge {r_0} \end{array} \right.\\ {\mathit{\boldsymbol{W}}^{(j)}} = {\rm{diag}}(w_1^{(j)}, \cdots , w_I^{(i)}) \end{array} $

作${\rm{SVD}},{\mathit{\boldsymbol{Y}}^{(j)}}{\mathit{\boldsymbol{W}}^{(j)}} = {\mathit{\boldsymbol{u}}^{(j + 1)}}{\mathit{\boldsymbol{s}}^{(j + 1)}}{[{\mathit{\boldsymbol{v}}^{(j + 1)}}]^{\rm{T}}} $

如果$\left\| {{\mathit{\boldsymbol{u}}^{(j + 1)}}{\mathit{\boldsymbol{s}}^{(j + 1)}}{\mathit{\boldsymbol{s}}^{(j + 1)}}{{[{\mathit{\boldsymbol{u}}^{(j + 1)}}]}^*} - \mathit{\boldsymbol{I}}} \right\|＜\varepsilon $, 停止迭代；其中：Y、r、w、W分别为预测矩阵、误差、权系数及权系数矩阵；diag(·)表示求对角矩阵。

(5) 利用最终权系数矩阵W，对$ {\mathit{\boldsymbol{\widetilde X}}^{\rm{r}}}{\bf{W}}$作SVD，${\mathit{\boldsymbol{\widetilde X}}^{\rm{r}}}\mathit{\boldsymbol{W}} = \mathit{\boldsymbol{U\bar S}}{\mathit{\boldsymbol{V}}^{\rm{T}}},\mathit{\boldsymbol{U}}、\mathit{\boldsymbol{\bar S}}、\mathit{\boldsymbol{V}}$分别为SVD所得的各参数矩阵。

(6) 计算场源极化参数$\mathit{\boldsymbol{\widetilde A}}:\mathit{\boldsymbol{\widetilde A}} = \mathit{\boldsymbol{\overline S}} {\mathit{\boldsymbol{V}}^{\rm{T}}} $。

2.3 估计目标测站的阻抗张量

利用最小二乘法，对目标数据矩阵Xⁿ及式(7)进行带参考信号的稳健回归分析，可求得系统响应

$ {\mathit{\boldsymbol{U}}^{\rm{n}}} = ({\mathit{\boldsymbol{X}}^{\rm{n}}}\mathit{\boldsymbol{W}}\prime \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{\widetilde A}}^{\rm{T}}}){(\mathit{\boldsymbol{\widetilde A}}W\prime \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{\widetilde A}}^{\rm{T}}})^{ - 1}} $

(9)

式中W′为对式(5)进行分析时的稳健估计权值矩阵。

获得Uⁿ的估计后，再利用阵列响应参数的提取方法计算阻抗及其视电阻率和相位^[30]。阻抗的计算公式为

$ {\mathit{\boldsymbol{Z}}^{{\rm{MT}}}} = [{\mathit{\boldsymbol{U}}^{\rm{E}}}{({\mathit{\boldsymbol{U}}^{\rm{H}}})^{\rm{T}}}]{[{\mathit{\boldsymbol{U}}^{\rm{H}}}{({\mathit{\boldsymbol{U}}^{\rm{H}}})^{\rm{T}}}]^{ - 1}} $

(10)

式中U^E和U^H分别为电场和磁场张量矩阵。上式省略了测站的角标编号。

分析以上过程可知，与常规方法不同，本文将多参考站数据集中于统一的阵列数据中进行分析处理，而不是简单的多点参考数据叠加。由于以上阵列数据分析的引入，提高了多参考站数据的利用率和天然场信息提取精度。基于噪声的统计特征，理论上参与计算的数据越多，噪声压制效果越好。

3 应用实例 3.1 实测多参考站观测试验

2016年10月，中南大学在中国东部地区开展了多点远参考观测及处理试验，测点分布如图 1所示。共包含目标测站1个(S1)，位于安徽省霍山地区，该处背景噪声较高；同步在背景噪声相对较低的地区布设参3个参考站，分别位于安徽省滁州市定远县(R1)、安徽省黄山市(R2)、江西省景德镇市(R3)，与S1的距离分别约为160、200、280km；各测站的张量布设方向均为正南或正北，以北向为x方向。R1、R2参考站观测水平4分量电磁场(E_x、E_y、H_x、H_y)，R3仅观测水平2分量磁场(H_x、H_y)。数据采集仪器采用加拿大凤凰公司生产的V5-2000型大地电磁仪，采样率为15Hz，同步采集时长约18h。

3.2 观测数据分析

阵列数据处理的理论依据之一是认为不同空间测站的观测数据中包含大量相干信息，并可从中提取出相应的场源极化参数。为此，分析了实例中同步阵列观测的时间域、时频谱数据，如图 2所示。不难看出，S1、R1、R2和R3尽管相距数百千米，所处的噪声环境各不相同，但各测站的同步电磁场间时间序列的总体轮廓及细节变化均具有相似性，功率谱随时间的变化趋势也一致，只是强度存在差异。这说明电磁数据处理可考虑空间多测站间的同步信息，进行多站联合处理。

为进一步分析各测站的含噪情况，进行了不同频率的极化方向对比分析，如图 3所示。可以看出，S1测站1Hz和0.1Hz的极化方向呈集中趋势，并且在整个观测时段趋势一致，说明该测站处包含稳定的人文噪声场，且影响贯穿整个观测时段。由于噪声超过了稳健估计的“崩溃点”^[32]，常规稳健估计方法很难得到可靠的处理结果。并且，因噪声包含在所有观测时段中，数据挑选类方法往往也难以凑效。各参考测站的极化方向分布相对均匀，仅R2测站略呈集中趋势。考虑到0.1~1Hz处于MT的“死频带”范围^[33]，结合后文的视电阻率、相位数据分析可知，各参考测站在0.1~1Hz范围内未包含明显的人文噪声，但天然场信号强度较低，分散的极化方向可能为天然场和随机噪声的混合表现。

3.3 处理结果对比

经过时频转换，先对张量观测的S1、R1和R2测站进行了常规Robust处理和常规阵列处理，获得了视电阻率和相位结果(图 4)。常规阵列处理中R3的两道磁场也参与了计算，只是因该点未观测电场，因此未获得视电阻率、相位数据。可以看出，S1测站在0.1~1Hz的MT“死频带”范围内出现了明显的畸变，视电阻率曲线在低于0.1Hz频段明显脱节，yx模式的数据表现尤为明显。而参考站R1和R2在0.1Hz附近频段也出现了一定程度的畸变，这正是MT“死频带”畸变的显著表现。由于“死频带”内天然场信号强度极低，甚至可能低于仪器的本底噪声，因此即使在人文噪声较低的测区，也难以获得高质量的天然场响应。这也是目前MT数据处理中最棘手的问题之一。经过常规阵列处理后，利用多测道数据压制了随机噪声的干扰，S1测站在0.1~1Hz内的形态得到明显改善，但yx模式视电阻率曲线在0.1Hz前后频段仍表现出一定的脱节，显示出相干噪声的影响仍然存在。这说明在阵列规模不大的条件下，常规阵列处理可较好地压制随机噪声的影响，但对于相干噪声，其处理效果不容乐观。

基于本文所述的多参考站阵列数据处理方法，对测站S1进行了处理。结合前述算法可知，多参考站阵列处理与常规阵列处理的区别主要在于天然场源极化参数的估计策略不同。图 5是两种算法估计得到的极化参数。由于极化参数为复数，图中显示的是其绝对值。可以看出，尽管两种算法得到的天然场源极化参数变化规律基本一致，但细节仍存在诸多不同，说明了本文方法有别于常规阵列处理，取得的处理结果也不同。进一步地，图 6给出了多参考站阵列处理的阻抗视电阻率计算结果。对比可见，常规阵列处理的S1点(图 4d上)yx模式视电阻率曲线的相对高(0.2~2Hz)、低(0.01~0.6Hz)两个频段出现脱节。而多参考站阵列处理(图 6g)曲线在全频段连续，显然更合理。由大地电磁响应曲线连续的特征^[1]可判断，图 4d上中0.2~2Hz频段的曲线为畸变形态，包含相了干噪声的影响。这说明多参考站阵列处理可取得优于常规阵列处理的结果。

图 6中还给出了利用常规远参考方法处理所得的结果，这里采用了不同的参考道组合。对比利用3个远参考站磁道分别对S1进行常规远参考处理的结果(图 6g~图 6i)与图 4所示的S1常规处理结果，可以发现3个远参考站数据对S1在0.1~1Hz的畸变均有一定程度的改善，人文噪声的畸变影响得到了压制。相对而言，R1的参考效果相对最差，视电阻率在多个频点仍有较强的畸变；R2的参考效果略好，但yx模式的误差棒最大。R3的参考效果最优，仅在0.8Hz、0.2Hz左右几个频点处出现了跳变。这也说明，实际中，通过优选远参考可以得到更好的处理效果。对比单参考处理与两两组合的多参考处理，可以发现，R1、R2组合多参考处理结果(图 6d)优于他们各自的单参考处理(图 6a、图 6b))，前者不仅减少了畸变的频点，也压制了误差棒；R1、R3的多参考处理结果(图 6e)和R2、R3的多参考处理结果(图 6f)分别优于R1的单参考处理(图 6a)和R2的单参考处理(图 6b)，前者均有效压制了0.8Hz处R3单参考处理结果里的跳变点(图 6c)。整体上，3个参考站的磁道组合多参考结果(图 6g)仅在0.2Hz频点处出现了跳变，取得了最佳处理效果。

图 6h和图 6i分别是3个参考站的电道组合及混合测道多参考处理结果，分别对应式(4)、式(5)两种参考数据组合情况。可以看出，远参考电道数据的参考效果不如磁道数据，混合测道的处理结果优于电道组合，但不如磁道组合。这与常规远参考处理的测道选择认识一致^{[6, 14, 34]}。实际中，可通过噪声分析，选择含噪较低的测道组合进行处理，并经过对比评价，优选最终的处理方案。特别需注意的是，电道参考可能会较大程度地改变曲线形态，实际应用时需进一步论证其可靠性。

4 结论

(1) 多参考站阵列数据处理技术有别于常规远参考、互参考及多点远参考叠加处理方法，该方法将多参考站、多参考道数据整合于统一的参考阵列数据集，利用阵列处理方法提取天然场源极化参数，可有效提高多参考站条件下的数据利用率。

(2) 实测含“死频带”畸变的MT数据处理结果表明，多参考站阵列数据处理技术可获得优于常规单站处理、常规阵列处理以及常规单参考站处理的结果。当存在多参考站条件时，可采用多种测站及测道组合进行处理MT数据，并对比优选处理结果。

需指出的是，强干扰条件下，“死频带”畸变MT数据的处理仍是十分棘手的问题。本文提供了一种改善质量的可行策略，但提高远参考站本身的质量并尽可能提高采集数据的信噪比仍是改善目标测站数据处理质量的最关键因素。