0 引言

电磁法是一种重要的地球物理勘探方法，按照源的种类可分为天然源电磁法(MT、AMT等)和人工源电磁法(CSEM、TEM等)^[1-3]。相比天然源电磁法，人工源电磁法克服了天然场源信号弱且随机性强的缺点^[4-5]，电磁信号具有更高的信噪比和分辨率。可控源以周期性的方波和伪随机信号为主。广域电磁法(Wide Field Electromagnetic Method, WFEM)拥有完备的理论体系及成熟的仪器设备，中国对其具有完全自主知识产权^{[3, 6]}。该方法通过一次激发—接收可获得多个频率的地电信息，提高了野外工作效率和抗干扰能力^[7-8]。与可控源音频大地电磁法(CSAMT)相比，WFEM摒弃了只能在“远区”观测而导致信号微弱的劣势，将“过渡区”与“远区”有机统一，显著扩大了观测适用范围以及探测深度^[9]；与磁偶源频率测深法(MELOS)相比，WFEM保留了计算公式中的高次项，定义了适合全域的广域视电阻率，不需要对测量结果进行“校正”，且只需要观测一个分量就可以获得视电阻率^{[3, 6, 10-11]}。目前，WFEM已经在页岩气勘探^[12-14]、金属矿勘探^[15]、工程物探及城市物探等领域得到越来越广泛的应用。

随着现代工业与科技的发展及人类活动范围的不断扩大，电磁干扰越来越强，噪声压制^[16]一直是困扰广大电磁工作者的关键问题，在一定程度上制约了电磁技术方法的发展。WFEM采用人工场源，信号发射机功率很高(最高可达250kW)，相较于天然场源电磁法，采集到的电磁信号信噪比明显提高。在实际勘探时，WFEM同样受到噪声干扰的问题，去噪仍然是数据处理的关键环节。此外，为了提高WFEM方法的纵向分辨率，须加强WFEM数据去噪方法的研究。因此，如何利用新方法对WFEM数据进行有效去噪是目前亟需解决的关键问题之一。

针对人工源电磁信号受噪声影响的问题，张必明等^[17]提出了自适应双向均方差阈值法，先对原始数据样本进行排序，采用迭代或递归的方式，以中点为界计算前、后两部分数据的均方差；再利用单一的阈值控制异常点，在少量频点信息被噪声淹没时，可有效剔除频谱数据中的粗大误差。但是，仅采用阈值控制可能会造成有效数据的过处理及部分噪声不明显的数据未被处理，同时该方法也没有考虑到整体曲线的圆润性。针对此问题，Mo等^[18]提出了一种灰色系统理论和稳健估计方法，通过灰色建模求解测量数据的标准差，结合阈值法识别、剔除频谱数据中的异常值，将最终保留数据的M估计值视为真实值。陈超健等^[19]提出了一种灰色判别准则和基于有理函数滤波的电磁数据去噪方法，通过灰色判别准则剔除各个频点频谱数据的明显异常值，再利用有理函数滤波充分平滑数据曲线，进而约束数据的二次优化处理，以剔除残余噪声。

上述频率域处理方法主要针对干扰小的电磁数据，可显著提高数据质量。若电磁采集工作受到持续性强人文噪声的影响，大部分频点会发生畸变，依赖于频率域功率谱挑选的去噪方法就会失效。为此，杨洋等^[20]提出频率域基于小波变换和Hilbert解析包络的CSEM噪声评价方法，在不增加野外工作量的同时可提取更丰富的频率信息，针对高阶谐波信号，实现了频率加密，提高了CSEM数据的纵向分辨率。针对时间域CSEM数据处理，Yang等^[21]利用不含噪声或仅含随机噪声的时序片段，求解出噪声频谱，获得了高质量的数据，但该方法的前提是选出一段高质量的时间序列片段，而精准地挑选出高质量数据的位置较难。Li等^[22]提出了基于快速傅里叶变换、互补总体经验模态分解以及与移不变稀疏编码结合的CSEM去噪方法，可精准去除工频干扰，校正基线漂移，并基于稀疏表示和字典学习剔除信号中的强干扰，重构获取高质量的伪随机有用信号，获得的视电阻率曲线和电场曲线更稳定、连续。上述时间域广域电磁信号处理方法均是对采集的电磁数据进行整体处理，虽在一定程度上提高了数据质量，但缺乏对信号和噪声的辨识环节^[23-24]。若能从CSEM信号和噪声类型的特征入手，剖析CSEM信号与噪声之间的定性及定量辨识关系，这会对后续电磁法高精度的信噪分离起到关键作用，并能避免时间域处理技术在噪声压制过程中的过处理，最大限度地保留CSEM有效信号。

模糊C均值(Fuzzy C-means, FCM)聚类属于无监督机器学习方法，主要通过优化目标函数计算每个样本点对所有样本类中心的隶属度，从而决定样本点的类属以自动进行数据分类，主要应用于图像分割、地震属性分析等^[25-26]。本文主要以广域电磁法7频波2频组(7-2频组)和7频波3频组(7-3频组)的数据为例进行分析与处理，其中7-2频组的频率为1、2、4、8、16、32、64Hz，7-3频组的频率为0.75、1.5、3、6、12、24、48Hz。由于广域电磁法伪随机信号属于周期信号，具备一定的特征及规律性，因而根据时域统计特征能对广域电磁7频波数据进行信号与噪声的定量识别。本文通过对比电场曲线和视电阻率曲线，分析、评价方法的处理效果。通过提取时域统计特征并结合FCM聚类方法进行广域电磁伪随机信号处理，提出了一种特征分析与聚类识别相结合的人工源电磁伪随机信号处理方法。首先，建立伪随机信号和典型噪声的样本库，分析样本库信号的时、频域特征；然后，利用时域统计特征结合FCM聚类进行信、噪识别，重构有用信号，并利用数字相干技术^[27]提取有效频点的频谱；最后，通过本文方法对合成数据及实测数据进行处理，验证方法的可靠性与有效性。

1 伪随机信号的时域特征

分析采集的广域电磁7频波时域信号可知，时间域序列中的伪随机信号是有用信号^[3]，但时域信号难免会受到各种异常波形的影响，导致信号发生跳变或衰减，从而在频谱上出现频点信息的淹没或丢失。因此，识别异常波形并压制其影响能有效提升数据质量。

本文通过提取简单的时域统计特征参数进行信、噪的定量识别，通过探讨量纲特征参数和无量纲特征参数识别广域电磁信号与噪声。信号中的无量纲特征参数主要包括峰值因子、脉冲因子、裕度因子、波形因子、峭度和偏度，量纲特征参数主要有最大值、最小值、峰峰值、平均值、方根幅值、方差、标准差和有效值等。本文主要研究最大值、峰峰值、峰值因子、脉冲因子和裕度因子等特征参数，这些参数有助于快速、直观地表征广域电磁信号与噪声之间的关系，为后续聚类算法提供有效的特征参数。本文重点介绍无量纲特征参数。

1.1 无量纲特征参数 1.1.1 峰值因子

峰值因子F_FZ是信号峰值X_p与有效值X_rms的比值，表征峰值在波形中的极端程度。其计算公式为

${F_{{\rm{FZ}}}} = \frac{{{X_{\rm{p}}}}}{{{X_{{\rm{rms}}}}}}$

(1)

1.1.2 脉冲因子

脉冲因子F_MC是检测信号中是否存在冲击的统计指标，其计算公式为

${F_{{\rm{MC}}}} = \frac{{{X_{\rm{p}}}}}{{|\bar X|}}$

(2)

式中|X|表示信号X的平均值的绝对值。

由式(1)和式(2)可知，F_MC和F_FZ的区别在分母，由于同一组数据的绝对值的平均值小于有效值，因而F_MC必定大于F_FZ，根据这两个特征参数可检测信号中的突变成分和冲击因素。

1.1.3 裕度因子

裕度因子F_YD的物理意义与峰值因子和脉冲因子相似，常用于检测设备与信号中的损伤和变化程度。其定义为

${F_{{\rm{YD}}}} = \frac{{{X_{\rm{p}}}}}{{{X_{\rm{r}}}}}$

(3)

式中X_r表示信号幅值的方根。

1.2 样本库特征分析

为了剖析实测广域电磁7频波数据中的伪随机信号与异常干扰波形之间的定量辨识关系，基于实测WFEM数据的信号和噪声特点，利用Matlab构建一个由模拟典型干扰和模拟伪随机信号形成的样本库，如图 1所示。从样本库中任意选取三类信号的时域波形及其频谱，其中样本库包含了30个伪随机信号、30个脉冲干扰信号和30个衰减干扰信号，每个样本信号的采样长度均为1200，采样率为300Hz。

通过观测噪声样本库(图 1)可知，时域脉冲干扰和时域衰减干扰导致原始伪随机信号(图 1a)发生异常突变(图 1b、图 1c)，幅值增大，相应的频谱及主频信息也出现不同程度的失真。由图 1a可见，未受干扰的伪随机信号呈周期性，幅值稳定，频谱也相对稳定，其频点信息被完整地保留。分析图 1b、图 1c可见，若时域信号中存在着干扰或异常波形，信号在时、频域出现严重畸变，无法反映原始伪随机7频波信号的固有特征。

图 2为样本库信号的特征参数统计结果。分析图 2可知，含干扰信号的特征参数明显大于伪随机信号的特征参数，伪随机信号的特征参数值稳定且幅值小。基于这些时域统计特征参数(最大值、峰峰值、F_FZ、F_MC和F_YD)，可较好地识别伪随机信号和含干扰的信号。因此，结合这些特征参数可提高后续聚类分析的信噪区分度和信噪识别效果。

2 FCM聚类分析

FCM聚类是一种典型的无监督机器学习方法^[28]，该方法按照“物以类聚，人以群分”的思想将样本点按照某种规律进行划分，这些规律是通过样本点的某些特征确定的，无需事先给定或约束。FCM聚类算法通过定义下面的最小化目标函数得到最优解^[29]

${W_m} = \sum\limits_{i = 1}^c {\sum\limits_{j = 1}^N {u_{ij}^m} } {\left\| {{x_i} - {c_j}} \right\|^2}\quad m \geqslant 1$

(4)

式中：m表示聚类的簇数；N表示样本数；C表示聚类中心数；u_ij表示第j个待分类对象属于第i个聚类中心的隶属度；x_i表示第i个样本；c_j表示第j个聚类中心；‖·‖表示数据相似性的度量。

通过选用0~1的随机数确定初始化隶属度矩阵U_C×N，使其元素u_ij满足以下约束条件

$\sum\limits_{i = 1}^c {{u_{ij}}} = 1\quad \forall j = 1, 2, \cdots , N$

(5)

针对上述约束问题，利用拉格朗日乘数法对目标函数进行求导，可得

${u_{ij}} = \frac{1}{{\sum\limits_{b = 1}^c {{{\left( {\frac{{\left\| {{x_i} - {c_j}} \right\|}}{{\left\| {{x_i} - {c_b}} \right\|}}} \right)}^{\frac{2}{{m - 1}}}}} }}$

(6)

计算每组的聚类中心V_i

${V_i} = \frac{{\sum\limits_{j = 1}^N {u_{ij}^m} {x_j}}}{{\sum\limits_{j = 1}^N {u_{ij}^m} }}$

(7)

使得目标函数最小，即欧式距离最短，相似度最高。迭代更新u_ij和V_i，直到前后两次隶属度的最大变化值低于预设的误差阈值，即目标函数达到最小值。

图 3所示为对图 1样本库中的90个样本提取的峰峰值、峰值因子、脉冲因子和裕度因子分别进行K均值聚类和FCM聚类的结果。

分析图 3可知：虽然K均值聚类分析能将样本库中的90个样本分为两类，但对于海量的电磁数据特征值较乱时，K均值聚类无法选取合适的聚类中心，导致聚类效果不够理想；FCM聚类分析通过计算每个样本点到聚类中心的欧式距离，自动将两类含干扰的信号与伪随机信号划分为不同的类别，聚类效果明显。因此，FCM聚类有利于后续的干扰信号剔除，保留伪随机信号。

3 仿真实验分析

为了验证方法的识别效果，对模拟样本库信号中的典型干扰类型进行仿真实验。图 4所示为模拟图 1所示样本库中的7-2频组数据的信、噪识别结果及其频谱，以及作为对比的小波去噪结果及其频谱。

从图 4可知，含噪信号(图 4a)经时域特征提取，根据聚类分析结果能准确地识别信号中的异常干扰，保留伪随机信号(图 4b)。结合频谱分析(图 4b)可知，主频信息受干扰影响，主频值不稳定，且异常程度不同。对比时频分析小波去噪结果(图 4c)可知，针对时域中的噪声进行去噪处理的同时，伪随机信号也被滤除，过处理导致无法还原伪随机信号的形态；同时，对应的主频信息因残余少量噪声导致主频幅值降低。经本文方法处理后(图 4d)，典型的干扰信号被有效识别，并被准确去除，噪声频谱也同时被滤除，伪随机信号的有效波形及其频谱原始特征被有效重构、还原。

为了定量分析本文方法的有效性，统计了模拟7-2频组信号经不同方法去噪处理后的电场幅值及相对误差(表 1)。

分析表 1可知，由于噪声添加在不同的时刻，导致时间域和频率域信号出现突变或混乱现象，且含噪信号在1~16Hz的5个频点的电场幅值均大于真实值，最大达1.301mV，对应频率为16Hz，相对误差为34.95%。因而，噪声在时域会严重影响伪随机信号的特征，在频域会影响不同的主频信号及其谐波。根据小波去噪方法处理结果(图 4c)可见，7-2频点的主频幅值整体更低，相对误差平均达14.10%，小波去噪方法虽能压制异常波形，但严重过处理导致了该方法失效。因而，经本文方法处理后的电场幅值更接近于真实电场值，相对误差也较低，最大误差仅0.52%。

4 实测数据分析 4.1 实测信/噪识别与处理

实测数据来自四川盆地某页岩气广域电磁法勘探项目。以某一段含典型干扰的实测CSEM信号(采样率1200Hz)为例(图 5a)，采用本文方法进行特征提取与信噪识别，结果见图 5b和图 5c。

分析图 5可知，实测信号中包含了脉冲干扰和衰减干扰，通过本文方法提取时域统计特征参数，结合FCM聚类，处理结果有效地识别了噪声(图 5b)和伪随机信号(图 5c)得到有效识别并得以保留。

CSEM的发射信号和接收信号都是伪随机信号，若接收到的时域信号为明显的伪随机信号且无任何异常波形的干扰，则频谱也会稳定，无主频点出现紊乱现象。对实测的无强干扰波形的7-2频组信号(图 6a、图 6b)，对其添加脉冲干扰和衰减噪声，得到加噪数据(图 6c、图 6d)，经小波去噪和本文方法处理，提取的噪声及其频谱见图 6e和图 6f，重构数据及频谱见图 6g和图 6h，经小波去噪后的数据及频谱见图 6i和图 6j。对应的广域电磁法E_x分量频域电场曲线见图 7。

从图 6a的实测数据(详见前4s的放大图)可知，数据整体上属于周期性伪随机信号，电场曲线相对平稳且无异常跳变。因此，该数据不存在强干扰的影响，数据质量较高^[20]。通过在此信号中添加典型干扰，信号发生突变，其频谱也相应地发生变化，体现在电场曲线上主频点的跳变及频谱混乱(图 6b和图 6c)。分析小波方法处理结果(图 6i、图 6j)可知，由于小波基函数、分解层数和阈值较难选取，小波去噪效果也不甚理想。由于时间域信号过处理现象和大量有用信息的丢失，频率域中噪声残余使得电场曲线不够稳定、不光滑。经过本文方法处理后(图 6g和图 6h)，基于电场幅值和频谱信息可有效识别干扰信号段并将其剔除，重构信号基本消除了人为添加的典型干扰，成功地还原出原始伪随机信号。进一步地采用数字相干提取技术(图 7)，对比加噪前、后的WFEM电场曲线可知，噪声干扰是主频点波动的主要原因，将其从数据中有效识别并剔除后，利用去噪后的数据可以基本恢复原始WFEM电场曲线的趋势及形态特征。

4.2 实测点分析

广域电磁7频波数据的频段范围为0.0117~8192Hz，实测点仅涉及7-2、7-3频组数据，因而本文分析的电场曲线和视电阻率曲线仅涉及这两个频组信号。这里以测点S1和S2为例，实测数据经本文方法处理后的结果见图 8和图 9。

分析图 8、图 9可知，测点S1和测点S2观测电场曲线和视电阻率曲线在频率低于24Hz时均呈锯齿状，发生跳变，说明观测数据的时间域波形中包含了脉冲干扰和衰减干扰；采用本文方法提取特征参数并进行FCM聚类，处理后的数据可有效识别噪声干扰，并进行剔除，利用数字相干技术提取的归一化电场曲线和计算得出的视电阻率曲线平稳、连续，可见本文方法基本消除了噪声干扰引起的频点跳变现象，重构了高质量的伪随机信号。

5 结论

针对CSEM数据受噪声干扰的问题，提出了基于特征提取和聚类识别的广域电磁数据伪随机信号处理方法。基于仿真实验及实测数据处理结果得到如下结论。

(1) 所提方法针对时间域数据进行信、噪辨识处理，能快速、准确地识别信号中的干扰或异常波形，重构并还原出伪随机信号，提升了数据质量，弥补了传统CSEM数据仅在频率域内开展去噪处理的短板。

(2) 模拟实验和实测数据处理结果表明，经本文方法处理得到了稳定、连续的归一化电场曲线和视电阻率曲线，能有效应用于电磁数据等周期信号的特征提取及信号与噪声辨识处理，为广域电磁法伪随机信号数据处理提供了技术支持。

本文方法简单且易于操作，能快速而有效地识别噪声，可提高广域电磁数据质量，为后续反演解释提供可靠数据。但本文方法仅在时间域进行信号与噪声识别处理，未分析如何去噪，因此，如何去除识别的噪声干扰是后续的研究重点。