0 引言

海洋可控源电磁(marine controlled-source electromagnetics，MCSEM)是一种感应类电磁勘探方法，通过测量和处理主动激发源激发的低频电磁场信号，获得海底地层的电阻率分布，从而推断海底储层的含油气性^[1-3]。MCSEM数据的系统处理流程已由国际电磁服务公司OHM和EMGS开发建立，限于商业服务竞争和知识产权保护，目前仍处于垄断和保密阶段^[4]。因此，研发具有自主知识产权的MCSEM数据处理方法并逐渐形成相应的处理系统，对于我国丰富海洋资源的勘探与开发具有重要的实际意义。

国际上，海洋可控源电磁数据预处理技术已有诸多进展。Behrens^[5]开发了基于Matlab环境的海洋电磁数据处理算法，处理流程包括：时频转换、椭圆极化方位校正、导航数据合并和随机噪声估计。Lu等^[6]提出了利用被动频率信号构建的噪声模型消除主动频率信号中噪声的去噪方法。Myer等^[7]应用一阶差分预白技术压制了大地电磁等频谱污染。Mattson等^[8]应用空间平均滤波和奇异值分解技术对拖曳式海洋可控源电磁数据进行去噪分析。Maclennan等^[9]把时间域的等效源方法引入频率域可控源电磁数据去噪之中。

国内，在国家863计划“深水可控源电磁勘探系统开发”课题的推动下，6所大学的不同课题组、2家企业的研究人员共同研发具有自主知识产权的MCSEM勘探系统，目前已取得诸多成果^[10]。中国海洋大学的李予国课题组提出一套MCSEM数据预处理方法，并应用于实际数据处理^[11-12]。中国地质大学(北京)的邓明课题组基于C++11语言开发了一套MCSEM数据可视化预处理软件^[13-14]，并把压缩感知技术引入对MCSEM数据去噪的处理中^[15]。同校的魏文博课题组讨论了不同时域滤波时窗长度对频域信号振幅的影响^[16]，并利用Hilbert-Huang变换对MCSEM数据进行时频分析^[17]。吉林大学的刘财课题组建立了MCSEM数据预处理流程和软件，并改进了压制噪声的平滑滤波和双边滤波方法^[18]。中国石油集团东方地球物理勘探有限公司在GMECS系统中开发和集成了MCSEM预处理模块，成功应用于中国南海北部油气区块^[19]。随着国家863计划“深水可控源电磁勘探系统开发”的结题，我国已建立具有自主知识产权的MCSEM数据处理系统，实现了零的突破^[10]。但与国际横向对比，我国对MCSEM数据处理方法的研究仍处于初步阶段。纵观国内已发表的MCSEM数据处理方法和软件，基本采用经典的海洋可控源电磁数据预处理流程：时频变换、椭圆极化方位校正、导航数据合并，未见关于“飞点”的规则噪声去除方法、时变源矩和时变偏移距的信号畸变修正。

本文从实际勘探的发射信号和接收信号分析开始，基于经典数据预处理流程，提出规则噪声的去除方法、时变源矩的归一化校正和偏移距的航速修正，并进一步利用合成孔径源方法改善数据信噪比。所有处理过程均使用某海域可控源电磁实测数据。

1 发射信号与接收信号分析

为了鉴别接收电磁信号中的有效信号和噪声信号，首先对发射电流和接收电磁信号的组成进行分析。在MCSEM勘探中，通常根据峰值频率对勘探目标的敏感度，选用基频能量高的方波信号^[20]，也称为“Cox”方波。这种双对称方波的能量主要集中于基频和三次谐频，每个周期在切换电流极性时有短暂的时刻没有电流通过，可实现宽频带和可控相位输出。电偶极源激发的电磁场经过海底地层、海水及空气层的作用，返回到多分量电磁接收器。采集到的电磁场信号由来自地层的有效信号、大地电磁场与其他噪声组成。

图 1给出了某海域实测发射机电流波形及其频谱。由图 1a可见，发射机的稳定电流约为800 A，周期为12.5 s。截取100~150 s时窗作傅里叶变换，并输出频段0.02~2.00 Hz于图 1b。图 1b显示，发射机电磁能量主要集中于基频0.08 Hz和三次谐频0.24 Hz，且与稳定电流的幅值比分别为0.823和0.762。

工作期间发射机随勘探船只沿测线拖动，接收器持续记录海底电磁场响应。图 2给出了26个接收器中第一个接收器水平电场E_x的时间序列。图 2a为MCSEM接收器工作时段的接收信号；图 2b为MCSEM发射源工作时段的接收信号，共16:35:47.625时长，其余为大地电磁和背景噪声信号。发射机开始工作时，发射机-接收器间距约为10 900 m，对应于电磁波在海水中的传播时长约为0.32 ms，小于接收机时间采样间隔1/1024≈0.98 ms，所以可以认为发射机工作时段即为MCSEM记录时间序列时段。需要注意的是，由于发射机发射信号与接收器的采样频率不同(发射机为1 024 Hz，而接收器为128 Hz)，需计算两组时间序列的交集时段以确定MCSEM数据的起始与结束时刻。

2 规则噪声压制

对E_x的时间序列信号应用经典预处理方法，得到E_x的振幅随偏移距变化(MVO)曲线，见图 3。在MVO曲线上，我们可以清晰观测到多个明显的异常值——“飞点”。对于明显的异常值，可以采用手动剔除的方法去除^[10]，但隐藏在正常信号之间的异常值无法通过肉眼观察识别，仍保留在采集信号之中。注意，噪声的信号振幅也可能与有效信号相近。

对发射源工作时段的E_x时间序列进行离散短时傅里叶变换，并绘制其功率谱(图 4)。从图 4中可以观测到不同频率下信号的分布，同时也可以观测到多条浅色竖线。我们高亮这些浅色竖线并把E_x的MVO曲线叠加覆盖在功率谱之上，得到图 5。图 5清晰地表明飞点噪声具有独特的功率谱特征。我们抽取图 5中第393、514和1197号时窗的规则噪声功率谱曲线和非规则噪声功率谱曲线，并绘制于图 6中。由图 6a可见，在非激发频段上，规则噪声的功率谱呈现相似的曲线形态，功率谱随频率的增大而减小；反观一般信号(图 6b)，在非激发频段上呈现不规则的变化，无相似性特征。

通过上文的分析我们可以发现，此飞点噪声并非随机噪声，而是具有一定特征的规则噪声。利用其信号的独特性和普遍性特征可以自动压制所有此类规则噪声。我们提出了5步半自动去噪方法：1)在单时窗功率谱曲线上光滑激发频率上的数值；2)人工挑选一条明显的规则噪声，并记录其时窗号；3)计算其余时窗的功率谱曲线与第2)步挑选规则噪声功率谱曲线的相关性；4)筛选出相关系数大于0.85的时窗；5)校正。

MVO曲线去噪效果如图 7所示。通过调整筛选相关系数的阈值可以控制去除规则噪声的程度。另外，可以进一步使用图像识别技术代替第2)步的人工挑选工作，实现全自动去噪。此规则噪声可能是接收机中的硬盘在工作时磁盘旋转引起的噪声^[5]。

3 MCSEM测量数据校正 3.1 时变源矩计算与归一化处理

经时频转换后的MCSEM频域电磁场数据需要依据发射源源矩进行源矩归一化处理。本文所采用电磁测量数据的电偶极源由2条拖缆上的前后两个电极组成，两电极间距约为300 m。对于埋深远大于300 m的地质体而言，该电极组激发的电磁场可以近似认为电偶极子源激发的电磁场。由于拖缆的可延展性和施工环境的复杂性，摆动的拖缆导致电偶极子源的长度随时间变化，所以在作归一化处理时，需要考虑电偶极子源长度的变化，对不同时窗的电磁场数据采用不同的源矩归一化参数。设发射源源矩(M)为

式中：I为电流，文中实测数据使用的电流为稳定电流800 A；L为源长。

图 8为实测电偶极子源源长随时窗变化曲线，平均源长约为304 m，最大与最小源长差约为26 m。图 9为时变源矩与固定源矩归一化的振幅谱，图中虚线框内是局部振幅谱。由图 9可知，20 m左右的源长变化对于E_x振幅谱存在较小影响。

3.2 偏移距船速修正

前文构建MVO曲线时，需要与导航数据合并，即实现对振幅谱与相位谱进行时窗-偏移距转换，其中的偏移距计算至关重要。由于海洋环境的不可控，海上作业难以确保船速均匀，因此，在时窗-偏移距转换时应对时变船速加以修正，以提高MCSEM数据的定位精度。

图 10给出了某工区实测船速随时窗的变化曲线。平均船速1.38 kn，约0.71 m/s，在整个发射源工作时段内，船速在1.06~1.59 kn之间变化。图 11为原始偏移距、修正偏移距和偏移距修正量随时窗变化曲线。由图 11可见，300号和900号时窗附近偏移距受船速影响较大；其中300号时窗存在约100 m的负修正量，而900号时窗存在约190 m的正修正量(正负号表示沿测线正向或负向修正)。图 12为船速修正前后E_x分量的MVO曲线。结合图 11和图 12可知，由于船速不均匀导致近偏移距与中到远偏移距之间存在最大约190 m的偏移据误差。

4 合成孔径源方法随机噪声压制与信号增强

在电磁响应修正的基础上，为了同时压制随机噪声和增强有效信号，我们应用合成孔径源方法。合成孔径源(synthetic aperture source，SAS)方法最初由Fan等^[21]引入海洋可控源电磁领域，用于增强高阻油气储层的可探测性。本文应用这一方法实现数个单源响应的优化叠加，从而构建一个加长稳健的合成源响应^[22-23]，增强来自海底储层的有效信号，以提高海底电磁响应的信噪比。

4.1 方法原理

对于来自N个激发源的复电磁场数据，合成孔径源信号S定义为^[23-24]

式中：F为来自单源s_n接收位置r的电磁场信号；为电磁场的波数，ω为角频率，μ为磁导率，σ为电导率；Δr_n= | r_n-r₁ |为第n个源与第1个源之间的距离；c₁和c₂分别为相移因子和振幅加权因子。为分析合成孔径源的应用效果，我们将标准化幅度作为可探测性D的度量：

式中，S^b为背景模型下的合成孔径源信号。

现在我们应用粒子群算法^[25-26]优化c₁和c₂。在粒子群算法中，视向量x为一个在高维度坐标系中可自由移动的粒子。由于通常有多个粒子(例如10~100)同时以特定交互规则搜索目标，所以称其为粒子群。每个粒子代表最优化问题的一个解，对于合成孔径源问题即为权重因子向量(c₁, c₂)；最优化问题的目标函数则是可探测性D。通过单个粒子简单的运动和粒子间特定的交互规则，寻找可探测性的最大值。

每个粒子初始位置设置为搜索空间中的随机值。当粒子群开始搜索时，所有粒子同时向下一个位置移动：

式中：x_i^k为第i个粒子在k次迭代时的位置向量；v_i^k为第i个粒子在k次迭代时的速度向量；p_i为第i个粒子目前寻找到的最佳位置向量；p_g为所有粒子目前寻找到的最佳位置向量；α和β为正常数；U为单位随机数生成函数。

由式(5)可见，在每次迭代中每个粒子的速度向量由三部分组成：第一部分是上一次迭代的速度，我们设置初始速度向量为v_i⁰=(0.01, 0.01)；第二部分是第i个粒子当前位置与其历史最佳位置的距离(自我认知部分)；第三部分是第i个粒子当前位置与所有粒子的历史最佳位置的距离(社会认知部分)。我们设置α=0.5，β=1.5^[27]。

当粒子群开始运动后，每个粒子对比其当前的D和其历史最优D值。如果目前D即为历史最优值，则p_i更新为x_i。在所有粒子完成移动后，如果p_i是历史最优值，则p_g更新为p_i。

下面，给出粒子群最优化算法步骤：1)初始化参数；2)计算每个粒子的目标函数；3)顺序更新每个粒子的最佳位置向量和所有粒子的最佳位置向量；4)所有粒子移动至下一位置；5)循环步骤2)—4)直至粒子移动向量小于某一阈值。

4.2 数值模拟和实测数据应用

鉴于文中的测量数据均为二维测线，下文分别基于单激发源测线的模拟数据和实测数据，应用合成孔径源方法，分析其增强有效信号的效果。模拟数据来自于三维崎岖海底储层模型。表 1列出了去除海底崎岖的背景模型参数。

100 Ω·m的异常体规模设置为4 km×4 km×100 m，其顶界面位于海底以下2 km处，异常体轮廓和崎岖地形轮廓如图 13所示。接收器以均匀间距放置于海底，x方向间距为500 m，y方向间距为2 km，如图 13所示。

我们计算合成孔径源(10个单源)的响应，然后比较单源与合成孔径源可探测性的差异。图 14给出了合成孔径源与单源的可探测性分布。对于海水厚度为300 m的情形，合成孔径源方法提高了可探测性，在优化相移因子和振幅加权因子后，最大可探测性提高了5倍。

对于实测数据，利用测线上前6个接收器响应进行合成孔径源处理，并与单源(第一个接收器)作对比。图 15为单源与合成孔径源电场E_x的MVO曲线。从图 15看到，合成孔径源长大约6 km，其响应幅值与单源相比在合成源右侧偏移距大幅度增高，有效信号能量增强，进而提高信噪比。

5 结论

本文对经典海洋可控源电磁数据预处理流程提出了预处理信号的进一步修正与增强方法。利用电磁信号的功率谱信息识别并去除规则噪声，分时窗矫正船速和发射源长度以及应用合成孔径源方法实现弱信号增强。将文中建立的处理方法应用于数值模拟数据及某工区的实测数据，得到如下认识：

1) 信号功率谱特征是识别规则噪声的有效属性，使用人工或图像识别技术可以半自动或自动压制规则噪声。

2) 对频率域电磁响应的时变源矩归一化和船速偏移距修正，提高了电磁响应对地下目标电性参数变化的分辨准度。

3) 合成孔径源方法的应用，在不增加采集成本的条件下，进一步增强了有效信号的强度，提高了信噪比。