0 引言

电磁法是一种重要的勘探地球物理方法，主要有被动源电磁法和人工源电磁法(Nabighian, 1987; Chen et al., 2011; Ren et al., 2013; Grayver et al., 2014; 任政勇等, 2017; Chave et al., 2017).相对于被动源电磁法，人工源电磁法具有信号强、分辨率高的特点，场源主要有方波信号和伪随机信号两种，其中以伪随机信号为场源的电磁法可通过一次发送和接收获得多个频率的地电信息，提高了野外工作效率和抗干扰能力(何继善, 2010a, 2010b; Ziolkowski et al., 2011).由于大地是一个开放的环境，在数据采集过程中不可避免地会受到各种噪声的干扰，尤其是人文干扰，幅度大且频率范围宽，常使观测信号畸变，为电磁法工作的开展和推广带来了极大的困扰(Szarka, 1988; 汤井田等, 2017, 2018; 杨洋等, 2018; Yang et al., 2018; Chave, 2017; Li et al., 2019).因此，如何有效压制强电磁干扰，为后续处理解释提供可靠数据仍然是当前电磁勘探中极为重要的研究课题.

伪随机多频电磁信号通过编码一次发送可以同时接收到多个频率信号，同时在人工源电磁法勘探中，通常进行冗余测量，获得各频点的多个频谱数据，具有一定的抗干扰性.然而，测量的信号包含了有效信号和干扰信号.在人工源数据处理中，常将剔除异常值后的频谱数据的均值作为测量结果，从而获取可靠的有效信号.为压制噪声影响，国内外学者提出了多种数据处理方法.稳健阻抗估计是一种主流的提高阻抗估计质量的方法(Egbert and Livelybrooks, 1996; Chave and Thomson, 2003; 周聪等, 2015; Liu et al., 2016).但由于场源不同，人工源电磁法场源极化方向固定，不同方向上的信号能量差异大，在近区平面波假设不再成立，稳健阻抗估计在人工源电磁数据处理中可能导致错误的结果(Streich et al., 2013)；小波变换过分依赖于小波函数的选择，且随着尺度的增大，相应正交基函数的频谱局部性变差，使信号的精细分解受限(Trad and Travassos, 2000; 凌振宝等, 2016)；数字滤波技术需综合考虑电磁测深方法、仪器、有效信号时频特征甚至是工作环境等多种因素，在噪声和有用信号在大范围相互重叠的情况下，去噪效果不佳(Olsen and Hohmann, 1992; 邓靖武等, 2004)；采用人机联作去噪虽能较好改善数据质量，但操作不可避免地引入了主观因素，编辑结果因人而异，且工作量巨大.

张必明等(2015)采用自适应双向均方差阈值法，采用单一阈值控制异常剔除，该方法在少量频点被污染时能有效压制噪声影响.然而，实测的各频点数据受噪声影响程度不一致，采用单一阈值控制可能会出现有效数据被剔除或者噪声数据未被处理的情况，且该方法未考虑曲线整体的圆滑连续性.Mo等(2017)提出了一种基于灰色建模和稳健M估计的人工源电磁数据处理方法，并结合阈值法大幅提高了对异常值占比的容忍度及处理结果的效果，处理后频谱数据曲线趋于光滑.但由于需要多次试验来寻找合适的阈值，数据处理比较耗时.总的来说，现有的各种单一的去噪技术均具有一定的优势，同时也存在一定局限性.

作为平滑滤波的一种，有理函数滤波克服传统滤波算法在保护信号细节时，耗时长、平滑能力低等缺点，在函数近似、插值和图像处理等方面得到广泛应用，并可用于恢复干扰少的频率域电磁场数据曲线的基本形态(Ramponi, 2001; 严家斌和刘贵忠, 2006; 雷达等, 2010).但因其过度依赖于原始曲线的趋势分析，当干扰严重出现连续飞点或某一频段形态失真时，会造成资料的二次污染，且滤波参数的选择不当会使滤波结果脱离真实情况(黄皓平, 1991).

针对上述问题，本文提出了一种基于灰色判别准则和有理函数滤波的伪随机电磁信号数据处理方法.在滤波前，引入灰色判别准则对原始数据进行初步处理，剔除明显的异常值，还原曲线的基本走势.然后对经灰色判别准则处理后得到的曲线进行有理函数滤波，用充分接近真实值的滤波结果约束数据的二次优化处理，从而达到剔除残余噪声的目的，并克服了滤波参数选择不当可能导致的滤波结果脱离真实值的问题.且由于经灰色判别准则处理后大部分噪声干扰消除殆尽，故选择相同滤波参数处理干扰程度不一的数据同样可以有效剔除干扰，且在一定范围内选择任意参数都可大幅改善数据质量.

1 方法原理 1.1 灰色判别准则

灰色判别准则实际是利用灰色系统理论中的累加生成法，通过累加后得到累加曲线，从曲线变化规律获取信息并用于离群值判别，原理如下(施冬等, 1998; 朱坚民等, 2010; Kayacan et al., 2010; Ayvaz et al., 2014)：

对于某测点的某一频点的伪随机电磁频谱数据，可以表示为序列x=(x(k)|k=1, 2, …, L)，其中，x_k表示第k个频谱数据，L表示该频点频谱数据的个数.当实测数据受到噪声干扰时，电磁频谱数据序列存在无规律的、随机的或有明显的跳动.将序列从小到大排列成序列x⁽⁰⁾，有：

(1)

其中，x(k)≤x(k+1)(k=1, 2, …, L－1).

将排序后的电磁频谱数据序列x⁽⁰⁾经一次累加生成后，可以获得新序列x⁽¹⁾：

(2)

如图 1中实线2所示.其平均累加序列(如图 1中虚线1所示)可以表示为：

(3)

其中，表示频谱序列x所有元素的平均值.为减少异常值对判别结果的影响，采用稳健估计值代替传统的算数平均值.那么，实际测量累加序列x⁽¹⁾与平均值累加序列的差值(如图 1中虚线3所示)可以表示为：

(4)

其中：表示实测累加序列x⁽¹⁾中与平均值最大差值.通常情况下，h=3.75(Wang et al., 2002).p的数值取决于频点频谱数据的个数L的奇偶性，可以表示为：

(5)

根据测量值不确定度的灰分析，Δx系列随着测量次数k的增大，由0开始逐渐增加至最大值，随后由于随机误差的叠加相消效应，又逐渐减小到0.若测量中不出现异常值，在随机误差的作用下，测量累加曲线变化并约束在一个三角区域内，如图 1所示.

如果，x⁽¹⁾(1)≤Δx(1)，那么认为相应的最小测量值x(1)为受到噪声干扰的异常值；如果，x⁽¹⁾(k－1)≤Δx(k－1)，那么认为相应的最大测量值x(k)为受到噪声干扰的异常值，然后相应的剔除对应位置上的异常值，从而达到剔除明显异常值的目的.将处理后得到的序列命名为y=(y_j|j=1, 2, …, M)，M为保留频谱数据的个数.

1.2 有理函数滤波

将经过灰色判别准则剔除明显异常值后的伪随机电磁频谱数据进行有理函数滤波，进一步消除残余噪声的影响.有理函数通常定义为两个多项式的比值，只需较少参数就能表示信号的变换关系(Ramponi, 2001; 雷达等, 2010).对于序列y中的连续三个伪随机电磁频谱数据(y_j－1, y_j, y_j+1)，其中，1≤j－1＜j＜j+1≤M，其对称的线性低通滤波系统可以表示为：

(6)

式中，0 < λ < 1/3.上述滤波器的响应函数随着频率的增加而单调减小，具有低通滤波的特性，其截止频率为λ的反函数.因此，在使用该方法进行滤波时，具有压制高频信号的特点，从而使滤波后的信号模糊化.为了方便数据处理，将公式(6)写成如下形式：

(7)

式中，截止频率N=1/λ，I表示迭代的次数.

为了避免这类线性滤波器在进行数据处理时导致信号模糊化的问题，需要控制改正项的变化，并对信息的重要性进行估计，使其对重要的高频信息的影响变小.本文通过估计梯度值的大小来判定信息的重要性，采用中心差分代替微分，改进公式(7)中的有理函数滤波模型，改进后的有理函数滤波模型如下：

(8)

式中，K称为正参数，为一正数，即K>0.从而使得改进后的利用中心差分进行估计的滤波模型具有带通滤波响应特征.从式(8)可知，当K=0时，滤波模型退化为线性滤波方程，K→∞时滤波器将不被激活.

1.3 处理流程

综上所述，基于灰色判别准则和有理数函数滤波的数据处理算法步骤如下：

(1) 基于灰色判别准则剔除各个频点频谱数据中的明显异常值，取各个频点存留频谱数据的均值，得到序列X₁(i), i=1, 2, …, n，n为频点个数；

(2) 对X₁(i)进行有理函数滤波，得到滤波曲线X₂(i)；

(3) 计算X₁(i), X₁(i+1), X₁(i+2)构成三条折线段的斜率k₁, k₂, k₃.若存在k₁ < k₃ < k₂或k₂ < k₃ < k₁，则X₁(i+1)为局部极值点，并对相应频点的频谱数据进行处理；

(4) 采用公式(9)，将从小到大排序后的频谱数据分为两段[1, m₁]、[m₂, M]，基于灰色建模计算标准差σ₁、σ₂.若σ₁>σ₂，剔除第一个数据，反之剔除最后一个数据.更新数据，M=M－1，重新分段，直到数据个数小于设定的值.取每次处理后的稳健均值，计算与前一个频点的斜率k₁，保留与滤波曲线对应斜率最为接近的值；

(9)

式中，M为单个频点频谱数据个数；INT()为取整函数，k为奇偶因子.N为偶数时，k=0，N为奇数时，k=1.

(5) 处理完所有频点后，为消除残余异常对处理结果的影响，转步骤(3)，再处理一次转步骤(6)；

(6) 结束算法.

算法在处理过程中只剔除端点的数据点，删除后更新数据，在接下来的每次迭代中，算法不断趋向样本数据中分布最集中的一部分，使得算法具有自适应优化的特点；通过不依赖数据分布规律和数量的灰色建模方法计算标准差，通过稳健估计算法分配权重的方式来削弱离群值对估计结果的影响，减小甚至避免后续处理出现异常误判或有效值误删的情况，故算法适用离群值比例大的小样本数据；经灰色判别准则处理，大致还原真实信号的形态走势，大大降低异常对滤波的影响，以保证滤波曲线能高度逼近真实信号；用滤波曲线约束数据自动处理，避免人为控制，极大提高了处理效率和数据的可靠性.另外算法实现简单，设备配置要求低，适用于各种计算环境.

2 数据仿真实验 2.1 仿真结果评价标准

野外噪声干扰错综复杂.分析大量伪随机多频电磁信号，比较常见的强噪声干扰有方波噪声、工频噪声、充放电噪声、脉冲噪声等，这些噪声都是由于复杂的外界因素和人文因素引起的，其能量强、频带广(汤井田等, 2018; 张必明等, 2015; 李晋等, 2017; Liu et al., 2017, 2019).为测试算法效果并分析噪声的影响规律，在以下仿真实验中，向实测无明显干扰的伪随机多频电磁时域数据的各个频组中逐点添加同种噪声.各种类型噪声强度为原始数据的30~60倍，数据长度约为原始数据的1/5.将加噪前的时域数据分别分成若干段，每段时域数据长度为最低频率周期的整数倍，利用数字相干检测技术提取每段数据主频点的振幅.对每个频点的振幅求均值，可得到加噪前频谱数据曲线z(i), i=1, 2, …, n.用本文方法对加噪后的数据进行处理，得到处理后的频谱数据曲线s(i), i=1, 2, …, n.通过计算z(i)与s(i)的曲线相似度NCC、信噪比SNR和误差ERR来评价处理效果(Cai et al., 2009; 汤井田等, 2017; Li et al., 2017)，具体表达式如下所示：

(10)

(11)

(12)

NCC为1时表示两个信号完全一样，为0表示正交，为－1表示两个信号的绝对值相同而符号相反.曲线相似度NCC越接近1，信噪比SNR越大，误差ERR越接近0，表示处理效果越好.

2.2 不同滤波参数的仿真结果

频谱曲线滤波是算法的重要一步，选择合适的滤波参数可高精度还原理想曲线，提高数据质量.为分析滤波参数对数据处理的影响，对图 2所示受不同干扰影响的加噪数据进行处理.图 2a和图 2c分别为仿真所添加的工频干扰和脉冲干扰的频谱曲线，加入工频干扰后的曲线在中频段出现不同程度的凸起，其中48 Hz处的频谱约增大一个数量级，脉冲干扰影响范围在100 Hz以下，能量在10~0.1 Hz频段达到最大值，电位差频谱曲线的100 Hz以后频段受到不同程度的影响，0~1 Hz频段最为严重.采用本文方法进行去噪处理，处理结果如图 2b和图 2d所示，频谱曲线变得光滑，逼近原始数据.

从图 3和图 4中可以看出，被噪声污染的频谱曲线比原始曲线的相似度NCC低，信噪比SNR为负数，误差大.保持其他两个参数不变，依次改变式(8)中迭代次数I、正参数K和截止频率N.两图曲线变化趋势一致：(1)固定K=0.1，N=100，随着I值的增加，曲线相似度NCC和信噪比SNR值变大，误差ERR下降，当I达到某一数值时，各项评价指标均达到最优.I值的持续增大导致仿真效果缓慢变差，最后趋于稳定；(2) N值的选取对仿真效果的影响与I值类似，固定K=0.1，I=70，当N等于某一个值时，处理效果最好；(3)固定I=70，N=100，K值取较小时，结果最佳，随着K的增大，NCC和SNR迅速下降，ERR急速上升，随后缓慢变化.其中，对于干扰程度较大的信号，选择更大的迭代次数可以更接近原始信号.

综上分析，过大的I、K和过小的N会导致滤波曲线与待处理数据曲线相差甚远，虽然曲线平滑但会丢失大量特征信息；反之，滤波曲线虽足够接近待处理数据，但因保留过多原始信息不能有效剔除离群值.故需选择在滤波曲线充分逼近原始信号同时能保持圆滑、连续走势的滤波参数.此外，由于经灰色判别准则处理后大部分噪声干扰消除殆尽，故选择相同滤波参数处理干扰程度不一的数据同样可以有效剔除干扰，且在一定范围内选择任意参数都可大幅改善数据质量，说明本文算法具有较强的适应性，为数据的批量处理提供一种可能的途径.

采用灰色判别准则剔除明显异常后，根据滤波参数的作用规律，选择滤波参数I=30，N=70，K=0.1和I=100，N=70，K=0.1分别对受工频、脉冲干扰影响的电位差频谱曲线滤波，得到参考曲线.图 4b、d中，经本文方法处理后，噪声干扰被有效消除，能量高达原始信号几个数量级的异常得以剔除，高精度还原出加噪前的曲线.表 1列出处理前后的定量评价标准.本文方法处理后的伪随机多频电磁频谱数据曲线与原始曲线相似度NCC分别提高到0.9999、0.9995，信噪比SNR上升约100 dB，误差显著下降.

2.3 强干扰信号处理

在上述高质量伪随机电磁数据中分别加入充放电噪声和方波噪声，两类噪声的频谱曲线如图 5a、c所示.两种噪声干扰范围主要集中在中低频段.其中，充放电噪声的主要频率范围在100 Hz以下频段，能量主要集中在10~0.1 Hz，相应频段的电位差频谱曲线有明显的突变，上升近3个数量级；方波噪声的频谱能量影响范围主要在10 Hz以后的频段，在电位差频谱曲线中表现为10 Hz以后的曲线迅速抬升，最大可达原始数据的3个数量级.上述噪声都造成伪随机多频电磁数据质量严重下降，中低频段出现不同程度的飞点，曲线形态不明确；噪声幅度大于原始信号的若干倍甚至几个数量级，信噪比SNR降至负数，曲线相似度NCC低，误差极大.

采用本文方法进行去噪处理，并将之与传统的采用自适应双向均方差阈值法(张必明等, 2015)进行去噪处理的结果进行对比，如图 5b、d所示.两类方法均能基本消除噪声造成的异常，处理数据的NCC、SNR以及ERR三项定量评价标准都明显改善，如表 2所示.经自适应双向均方差法处理的频谱曲线存在部分频点数据被过度剔除导致曲线凹陷，或异常未剔除干净的现象.由于该方法采用单个阈值控制处理全部频点频谱数据，而各频点受干扰程度不一，导致有效数据的误删或异常得不到有效剔除，且控制阈值的选择在一定程度上引入了人为控制的因素.采用本文方法处理后的频谱曲线的中低频段突跳频点得到有效恢复，电位差频谱曲线干净圆滑，连续性提高，与原始频谱曲线基本重合.处理后的频谱曲线与原始曲线的相似度高达0.9997，信噪比分别提升至62.573 dB和59.759 dB，误差值显著下降.

3 实测数据处理

为验证本文方法对伪随机电磁实测数据的处理效果，选择重庆某地区的伪随机多频电磁法勘探数据进行去噪处理.该工区内测量环境复杂，干扰众多.部分测点受附近矿山、高压电线、公路等的干扰，噪声源能量大，甚至完全淹没有效信号.现选择存在明显强干扰影响的3线进行处理，并从两个方面对数据质量综合评价：一是基于统计学方法，利用公式(13)计算各个频点数据的相对均方误差大小，一般而言相对均方误差在10%及以下表明数据可靠，相对均方误差在5%及以下表明数据质量良好(张必明等, 2015)；二是在对数坐标系中观测频谱曲线的形态.由于伪随机信号各主频点之间的相关性高，未受或受干扰程度小的曲线形态平滑，表明数据质量好.δ_r的计算公式如下：

(13)

式中，δ为均方误差.

以3线的3115号测点为例讨论处理效果，原始数据和采用本文方法进行处理的结果如图 6所示.该测点在0.75~3072 Hz较为平滑，相对均方误差小，此频段数据质量较好；4096 Hz、0.5 Hz及以下频段曲线跳变剧烈，电位差频谱曲线连续性差.对0.5 Hz及以下各频点原始数据分布进行对比分析，如图 7a，数据中存在明显的、能量强的飞点，且随着频率的降低，异常幅度及受干扰范围都有所增加.基于灰色判别准则处理后，各尖峰形态基本消除，误差棒显著减小，曲线走势恢复正常.图 7b中，噪声造成的异常均有效剔除，各频点数据分布集中.但由于灰色判别准则仅对单个频点的数据处理，未考虑整个频谱曲线.故在灰色判别准则处理的基础上，选择合适的滤波参数进行有理函数滤波，将能正确反映曲线形态的滤波曲线作为参考曲线，约束频谱数据的二次优化处理，经优化处理后的曲线相对误差变小，曲线光滑度以及连续性有所提高，数值回归正常，如图 6c所示.

图 8为该测线中四个测点处理前后的电位差频谱曲线图.每个测点都受到不同程度的干扰，主要体现在4096~2048 Hz、48 Hz以及1 Hz以后频段曲线发生畸变，连续性较差.尤其是3170号测点，整个低频段跳变剧烈，呈现“锯齿”状，其能量较周围频点高达一个数量级.选择滤波参数I=60, K=0.1, N=50，经本文方法处理后，曲线形态明显改善，噪声导致的严重跳变、锯齿得到有效剔除，有用信号得以保留.原本偏倚的数据、脱节的“飞点”得以校正，整个频段曲线光滑，无毛刺.

图 9给出了该工区3线连续22个测点的原始视电阻率和采用本文方法处理后的视电阻率曲线图.从图中可以看出，经过本文方法处理后，各测点处的视电阻率曲线变得平稳，相邻测点的视电阻率曲线走势比较一致：1024~8192 Hz频段，视电阻率呈平缓的下降趋势；随着频率的降低，视电阻率升高并在128 Hz附近到达极值点，随后逐渐下降；1~0.1 Hz频段的曲线变化平稳，小于1 Hz频段视电阻率递增.为了对比处理前后的数据质量，计算3线各测点处各频点电位差频谱数据的相对均方差，并取平均值，分布情况如图 10所示.图 10a中，约半数测点的平均相对均方误差分布在6%~10%，其他测点的平均相对均方误差分布于10%~20%区域，表明数据受到噪声的干扰.经本文方法处理后，数据质量得到明显改善，各测点的平均相对均方误差显著下降，大部分测点平均相对均方误差 < 5%，数据质量良好，少数测点分布在5%~10%.上述实验表明，本文方法可用于在频率域对受到噪声干扰的数据进行处理，在频谱数据个数少(范围为15~24)、干扰幅值大(高达有效信号幅值的几个数量级)的情况下，仍可有效压制甚至去除强干扰，保留了原始观测数据的整体特征和其他频段的高质量数据，提高数据信噪比，恢复原有曲线形态.

4 结论

针对伪随机电磁法存在的强干扰噪声，本文提出了一种基于灰色判别准则与有理函数滤波的数据处理方法，通过系统仿真实验及实测数据处理得出以下结论：

(1) 本文方法通过灰色判别准则剔除显著异常值，结合有理函数滤波进行二次优化处理，处理过程中引入稳健估计、灰色建模等方法，在保留高质量频点数据的同时，可实现对频谱数据中异常的有效识别和剔除.

(2) 在仿真实验中，选取合适的滤波参数可高度还原加噪前数据曲线.噪声强度对滤波参数选择的影响较小，且在合适范围内选任意滤波参数均可明显提升数据质量，表明本文方法具有可行性和较强的适应性，为数据批量预处理提供一种可能.

(3) 经实测数据预处理对比分析，本文方法可显著改善强干扰背景下的电磁数据质量，采用本文方法处理后，视电阻率曲线变得光滑平稳，相对误差减小，整体连续性提高.

此外，本文提出的方法能在高效处理数据的同时，保证数据的质量和曲线整体的光滑度，为提高海量伪随机电磁数据的质量和处理效率提供了一种可能性.