CN108680966B - 海洋可控源电磁勘探噪声降噪效果评估方法 - Google Patents

Abstract

本发明提供了海洋可控源电磁勘探噪声降噪效果评估方法，具体涉及海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，通过将MCSEM信号一维正演、MT噪声一维正演、噪声降噪处理、一维OCCAM反演相结合，利用电阻率反演结果与电阻率输入值的均方相对误差作为降噪效果的量化评价指标。既可以利用正演的手段生成MT噪声，也可以直接输入MCSEM发射机不开机时间段的实测数据作为MT噪声，保证了该评估方法具有较强的灵活性。本发明利用所提出的MT噪声降噪效果评估方法，对短时窗结合预白化滤波方法的降噪效果进行了量化评价，并验证了其降噪效果的有效性。本发明提出的降噪效果评估方法应用范围广，填补了噪声降噪效果的量化评估领域的空白。

Description

海洋可控源电磁勘探噪声降噪效果评估方法

技术领域

本发明涉及海洋可控源电磁勘探降噪领域，具体涉及一种海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法。

背景技术

随着全球经济的持续发展，人类对自然资源的需求量不断增加。广阔的海洋蕴含着丰富的资源，具有巨大的开发潜力，因此许多石油公司纷纷把目光转向海洋。而海洋可控源电磁勘探(Marine Controlled-Source Electromagnetic，MCSEM)作为一种海洋油气探测新技术，能有效的识别海底高阻油气储层，明显提高钻井成功率，因此近年来在国际上得到了高度关注。

信号去噪是海洋可控源电磁勘探数据处理中的重要环节，直接关系着资料解释的质量。在MCSEM主要的噪声源中，对于系统内部噪声的处理已相对成熟；对于空气波干扰，在浅水区影响较大，但在深水电磁勘探中影响不大；对于由海水运动造成的电磁噪声，目前也有一些处理方法。但对于大地天然电磁信号(Magneto Telluric，MT)造成干扰的处理，目前国内外针对该问题的研究相对较少。

大地电磁场是由于太阳风与地球磁层、电离层间复杂的相互作用以及雷电活动所形成的一种天然交变电磁场。这些交变电磁场射入地下时，一部分被地下介质吸收衰减，一部分被反射。MT信号的频谱成分较丰富，频带范围相对较宽。频率低于1Hz的大地电磁信号主要来源于太阳活动的粒子流和电辐射，而高于1Hz的大地电磁信号主要来源于人文活动和气象运动。同时MT信号在空间和时间上有着明显的分布规律。高纬度和中纬度的MT信号比较强而低纬度相对较弱；夏季强，冬季弱；白天强，夜间弱。MT信号在海洋中的传播过程与陆地相比是不同的。在海洋环境中，海水是一种高导电体，电磁波因为海水的天然屏蔽作用，其高频成分迅速衰减。因此，海洋环境类似于低通滤波器，允许低频成分顺利通过，屏蔽高频成分。由于海洋的这种特殊环境，使得MT信号的频谱主要集中在低频段，与MCSEM信号的频谱混叠，特别是在天然电磁活动比较强烈时其影响较难处理，成为影响MCSEM探测质量的主要噪声之一。

在MCSEM勘探过程的中、远收发距处，MCSEM信号相对较弱，来自海底地层的有效信号将可能淹没于MT噪声中，从而难以识别和提取。因此，若不进行专门的MT噪声抑制，将可能造成无法有效探测出海底油气储层。因此，针对MT噪声的降噪开展研究是非常有必要的。

Ryhove(2008)提出了一种基于远参考道技术的MT噪声降噪方法，但参考道布设区域与实际勘探区域间的适宜距离不易确定，且去噪过程的计算量较大。刘宁(2015)提出了一种时变双边滤波方法，但此方法滤波参数较多，不易控制，且计算时间相对较长。D.Myer(2011)提出了一种短时窗结合预白化滤波的方法，通过差分运算实现预白化滤波进而有效抑制低频噪声的影响，通过减小叠加处理中的时窗长度降低频谱泄露效应的影响。需指出的是，由于实测数据中信号与噪声相混合，真实信号的情况未知，导致上述几种方法的有效性缺乏降噪效果的量化评价作为佐证。因此，需要提供一种针对现有MT噪声降噪方法降噪效果量化评价的方法。

发明内容

针对现有的大地天然电磁信号噪声降噪方法缺乏降噪效果量化评价的问题，本发明将MCSEM信号一维正演、MT噪声一维正演、噪声降噪处理、一维OCCAM反演相结合，提出了海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法。

本发明采用以下的技术方案：

海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，包括以下步骤：

步骤1：根据需求的不同，判断是否采用海洋可控源电磁勘探发射机不开机时间段的实测数据作为MT噪声，若不采用，则进行步骤2；若采用，则进行步骤6；

步骤2：输入地层层数、深度、电阻率等参数，参数用于MCSEM信号一维正演和MT噪声一维正演；

步骤3：设置用于MT噪声一维正演的海表面磁场强度、接收站个数和坐标、海水层深度和电阻率等输入参数，其中接收站个数和坐标、海水层深度和电阻率等参数也用于MCSEM信号一维正演；

步骤4：利用MT噪声一维正演模型，生成MT噪声数据文件；

步骤5：设置用于MCSEM信号一维正演的发射源长度、初始坐标、航行速度、发射信号周期等其它输入参数，之后，进行步骤10；

步骤6：利用发射机开机时间段的MCSEM实测数据，采用一维OCCAM反演模型，反演得到地层层数、深度、电阻率等参数；

步骤7：读取MCSEM实测数据的相关配置文件，保存用于MT噪声一维正演的接收站个数和坐标、海水层深度和电阻率等输入参数，这些参数也用于MCSEM信号一维正演；

步骤8：选取MCSEM发射机不开机时间段的部分实测数据作为MT噪声数据文件，其数据点数应与MCSEM数据文件的数据点数相同；

步骤9：读取MCSEM实测数据的相关配置文件，保存用于MCSEM正演的发射源长度、初始坐标、航行速度、发射信号周期等输入参数；

步骤10：利用MCSEM信号一维正演模型，生成MCSEM正演数据文件；

步骤11：根据MCSEM正演数据文件，生成无噪数据的MVO曲线；

步骤12：将MCSEM正演数据文件中的电磁场数据与MT噪声数据文件中的电磁场数据相加，生成含噪MCSEM数据文件；

步骤13：根据含噪MCSEM数据文件，生成含噪数据的MVO曲线；

步骤14：利用短时窗结合预白化滤波等MT噪声降噪方法，对含噪MCSEM数据进行MT噪声降噪处理；

步骤15：根据降噪后的MCSEM数据，生成降噪后的MVO曲线；

步骤16：对比降噪前、后的MVO曲线，观察降噪效果，给出定性分析；

步骤17：利用一维OCCAM反演模型，对降噪前、后的MCSEM数据分别进行反演，得到反演后的地层层数、深度、电阻率等结果；

步骤18：分别画出真实电阻率与降噪前、后反演电阻率随深度变化的曲线，观察降噪处理对反演精度的提高，给出定性分析；

步骤19：分别计算降噪前、后反演电阻率与真实电阻率之间的均方相对误差；

步骤20：利用降噪前、后反演电阻率与真实电阻率之间的均方相对误差，对MT噪声的降噪效果给出最终的量化评价。

优选的，MCSEM信号一维正演模型为：

设立N层一维模型，其中，x轴和y轴位于水平面，z轴垂直向下，xyz坐标系满足右手定则，σ_m为第m层的电导率且为各向同性，z_m为第m层层顶的z坐标，h_m为第m层的深度，N为总共的层数；第1层总为空气层，z₁为-∞，将h₁设置为∞，第2层为海水，将z₂设置为0，第N层的深度h_N设置为∞；

假设时间谐变偶极子源的时谐因子为e^-jωt，忽略磁导率μ的变化，另外，角频率ω较低，导致位移电流也忽略不计，则麦克斯韦方程组可表示为：

▽×E＝jωB (1)

▽×B＝μσE+μJ_s (2)

其中，▽为汉密尔顿算子，E为电场强度矢量，B为磁感应强度矢量，σ为电导率，J_s为电偶极源的传导电流密度矢量；磁矢量势A与磁感应强度矢量B之间的关系为：

B＝▽×A (3)

由式(1)-(3)，电场强度矢量E可表示为：

通过汉克尔变换，可计算矢量势A：

其中，λ为信号波长，J₀为第一类零阶贝塞尔函数，为经过傅里叶变换后的矢量势，r为水平距离矢量，r为r的模；在MCSEM勘探中，经常使用水平电偶极子作为发射源；若电偶极源沿y轴方向放置，则矢量势具有如下形式：

其中，是满足的中间变量，且矢量势和是柱对称；当将电偶极子放置于第i层、深度坐标为z_s的位置处时，则第m层、深度坐标为z的位置处，其对应的矢量势分量和分别为：

其中，δ_mi为克罗内克函数，a_m、b_m、c_m、d_m是相关的系数项；

对于发射源所在层以上的各层，其横电模式的反射系数和横磁模式的反射系数分别定义为：

对于发射源所在层以下的各层，其横电模式的反射系数和横磁模式的反射系数分别定义为：

根据电磁场边界条件可得：

其中，h_m＝z_m+1-z_m，和对应γ_m+1、h_m+1，和对应γ_m-1、h_m-1；另外，和的表达式分别如下：

其中，和对应γ_m+1、σ_m+1，和对应γ_m-1、σ_m-1；

显然，顶层和底层处的反射系数为0，即有：利用电磁场的边界条件，进行递推即可求出a_m，b_m，c_m，d_m的值；然后，将a_m，b_m，c_m，d_m的值代入式(7)和式(8)，可计算出矢量势分量；最后，将式(5)代入式(3)和式(4)，即可分别计算出磁感应强度矢量B和电场强度矢量E；以上是以电偶极源沿y轴方向放置为例，最终得出的相关关系式。若电偶极源沿x轴方向放置，则推导过程类似，可得到一系列类似的关系式；当电偶极源在水平面沿任意方向放置时，则可通过叠加原理进行计算。

优选地，MT噪声一维正演为：

MT噪声一维正演所采用的坐标系与MCSEM信号一维正演模型所使用的坐标系保持一致；接下来，为了推导过程中符号表示的简便，同时又不失一般性，假设MCSEM发射机的摆放使得电场强度矢量的方向位于x轴，磁场强度矢量的方向位于y轴；上述假设下，电场强度矢量将只有x分量，磁场强度矢量将只有y分量；

假设射入海表面的平面波以e^-jωt规律进行变化，另外，设入射到海表面的磁场强度为H，即z＝0时，H_y(0)＝H，这里的常数H实际上表征了天然电磁场活动的强烈程度，另外，为了体现MT噪声随时间、空间随机变化的特点，H应满足实测MT资料的统计特性；根据边界条件以及电场和磁场的关系，第m层、深度坐标为z处的电磁场分量为：

其中A_m，B_m为待定系数，为波数；

利用数学递推方法和边界条件，可得：

其中，R_N＝1，和R_m的表达式为：

将式(19)、式(20)代入式(17)、式(18)，便可计算出MT信号的电磁场表达式。

优选地，一维OCCAM反演模型为：

OCCAM反演的实质是将无约束目标函数最小化，定义目标函数φ为：

φ＝||Rm||²+μ^-1(||Wd-WF(m)²-X_* ²) (23)

其中，R为一阶差分算子矩阵，m为模型矢量，μ^-1为拉格朗日乘子，W为权重矩阵，d为观测数据矢量，F(m)为模型矢量m对应的正演数据矢量，X_* ²为反演的目标拟合差；R和W的表达式分别为：

其中，σ₁,σ₂,…,σ_M为反演数据的标准差，M是观测数据的个数；

求取目标函数φ的梯度，使φ达到最小值，并令▽_mφ＝0，此时模型矢量m应满足：

μ^-1(WJ)^TWJm-μ^-1(WJ)^TWd+R^TRm＝0 (26)

其中，J＝▽_mF为雅可比矩阵，其元素可表示为：

其中，i＝1,...,M，j＝1,...,l，l为模型层数；

在求解式(26)时，需将F进行线性泰勒展开，即：

F(m+Δ)≈F(m)+JΔ (28)

其中，Δ表示迭代过程中模型的增量；

给定初始模型矢量m₁，并将式(28)代入式(23)，使φ最小化，得到k次迭代后的模型矢量：

m_k+1(u)＝[μR^TR+(WJ_k)^TWJ_k]^-1(WJ_k)^TW(d-F[m_k]+J_km_k) (29)

通过线性搜索μ，寻找F[m_k+1]与实测数据的最佳拟合差：

X_k+1(μ)＝||W(d-F[m_k+1(μ)])|| (30)

初始迭代阶段的拟合差可能大于目标拟合差，经过式(29)不断迭代得到满足精度的模型；在拟合精度满足要求的情况下，选择最大的μ对应的模型，因为此时对应的模型最光滑；在迭代运算过程中，如果均方根拟合差x_rms达到1时，则OCCAM反演结果收敛；其定义为：

其中，s_i为第i列数据的标准差。

本发明具有的有益效果是：

本发明提供的海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，将MCSEM信号一维正演、MT噪声一维正演、噪声降噪处理、一维OCCAM反演相结合，利用电阻率反演结果与电阻率输入值的均方相对误差作为降噪效果的量化评价指标。既可以利用正演的手段生成MT噪声，也可以直接输入MCSEM发射机不开机时间段的实测数据作为MT噪声，保证了该评估方法具有较强的灵活性。本发明利用所提出的MT噪声降噪效果评估方法，对短时窗结合预白化滤波方法的降噪效果进行了量化评价，并验证了其降噪效果的有效性。本发明提出的降噪效果评估方法应用范围广，也可推广至其它类型噪声降噪效果的量化评估，填补了噪声降噪效果的量化评估领域的空白。

附图说明

图1为本发明的流程图。

图2为N层一维模型的示意图。

图3为MT噪声一维正演的流程图。

图4为MCSEM数据的正演流程图。

图5为短时窗结合预白化滤波的降噪处理流程图。

图6为降噪前、后OCCAM反演流程图。

图7为实施例1的空气层、海水层、海底地层参数图。

图8为无噪MCSEM信号电通道MVO曲线图。

图9为无噪MCSEM信号磁通道的MVO曲线图。

图10为含噪MCSEM信号电通道的MVO曲线图。

图11为含噪MCSEM信号磁通道的MVO曲线图。

图12为降噪前、后含噪MCSEM信号电通道的MVO曲线图与无噪MCSEM信号电通道MVO曲线图。

图13为降噪前、后含噪MCSEM信号磁通道的MVO曲线图与无噪MCSEM信号磁通道MVO曲线图。

图14为MCSEM信号电通道降噪前真实电阻率与反演电阻率的变化曲线。

图15为MCSEM信号电通道降噪后真实电阻率与反演电阻率的变化曲线。

图16为实施例2中的空气层、海水层、海底地层参数图。

图17为无噪MCSEM信号电通道的MVO曲线图。

图18为无噪MCSEM信号磁通道的MVO曲线图。

图19为含噪MCSEM信号电通道的MVO曲线图。

图20为含噪MCSEM信号磁通道的MVO曲线图。

图21为降噪前、后含噪MCSEM信号电通道的MVO曲线图与无噪MCSEM信号电通道MVO曲线图。

图22为降噪前、后含噪MCSEM信号磁通道的MVO曲线图与无噪MCSEM信号磁通道MVO曲线图。

图23为MCSEM信号电通道降噪前真实电阻率与反演电阻率的变化曲线。

图24为MCSEM信号电通道降噪后真实电阻率与反演电阻率的变化曲线。

具体实施方式

下面结合附图和具体实施例对本发明的具体实施方式做进一步说明：

结合图1至图24，本发明将MCSEM信号一维正演、MT噪声一维正演、噪声降噪处理、一维OCCAM反演相结合，提出了海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，为提高方法的灵活性，既可以利用正演的手段生成MT噪声，也可以直接输入MCSEM发射机不开机时间段的实测数据作为MT噪声。

海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，包括以下步骤：

步骤1：根据需求的不同，判断是否采用海洋可控源电磁勘探发射机不开机时间段的实测数据作为MT噪声，若不采用，则进行步骤2；若采用，则进行步骤6；

步骤2：输入地层层数、深度、电阻率等参数，参数用于MCSEM信号一维正演和MT噪声一维正演；

步骤3：设置用于MT噪声一维正演的海表面磁场强度、接收站个数和坐标、海水层深度和电阻率等输入参数，其中接收站个数和坐标、海水层深度和电阻率等参数也用于MCSEM信号一维正演；

步骤4：利用MT噪声一维正演模型，生成MT噪声数据文件；

步骤5：设置用于MCSEM信号一维正演的发射源长度、初始坐标、航行速度、发射信号周期等其它输入参数，之后，进行步骤10；

步骤6：利用发射机开机时间段的MCSEM实测数据，采用一维OCCAM反演模型，反演得到地层层数、深度、电阻率等参数；

步骤7：读取MCSEM实测数据的相关配置文件，保存用于MT噪声一维正演的接收站个数和坐标、海水层深度和电阻率等输入参数，这些参数也用于MCSEM信号一维正演；

步骤8：选取MCSEM发射机不开机时间段的部分实测数据作为MT噪声数据文件，其数据点数应与MCSEM数据文件的数据点数相同；

步骤9：读取MCSEM实测数据的相关配置文件，保存用于MCSEM正演的发射源长度、初始坐标、航行速度、发射信号周期等输入参数；

步骤10：利用MCSEM信号一维正演模型，生成MCSEM正演数据文件；

步骤11：根据MCSEM正演数据文件，生成无噪数据的MVO曲线；

步骤12：将MCSEM正演数据文件中的电磁场数据与MT噪声数据文件中的电磁场数据相加，生成含噪MCSEM数据文件；

步骤13：根据含噪MCSEM数据文件，生成含噪数据的MVO曲线；

步骤14：利用短时窗结合预白化滤波等MT噪声降噪方法，对含噪MCSEM数据进行MT噪声降噪处理；

步骤15：根据降噪后的MCSEM数据，生成降噪后的MVO曲线；

步骤16：对比降噪前、后的MVO曲线，观察降噪效果，给出定性分析；

步骤17：利用一维OCCAM反演模型，对降噪前、后的MCSEM数据分别进行反演，得到反演后的地层层数、深度、电阻率等结果；

步骤18：分别画出真实电阻率与降噪前、后反演电阻率随深度变化的曲线，观察降噪处理对反演精度的提高，给出定性分析；

步骤19：分别计算降噪前、后反演电阻率与真实电阻率之间的均方相对误差；

步骤20：利用降噪前、后反演电阻率与真实电阻率之间的均方相对误差，对MT噪声的降噪效果给出最终的量化评价。

在上述步骤中，MCSEM信号一维正演模型为：

设立N层一维模型，如图2所示，其中，x轴和y轴位于水平面，z轴垂直向下，xyz坐标系满足右手定则，σ_m为第m层的电导率且为各向同性，z_m为第m层层顶的z坐标，h_m为第m层的深度，N为总共的层数。第1层总为空气层，z₁为-∞，将h₁设置为∞，第2层为海水，将z₂设置为0，第N层的深度h_N设置为∞。尽管在实际应用中，MCSEM的发射机和接收机都放置于海水层，但为了模型的一般性，发射机和接收机可设置在任意层中。根据上述模型，利用电磁场的相关理论，可推导出任意层任意深度处的电磁场表达式。

假设时间谐变偶极子源的时谐因子为e^-jωt，为简化起见，忽略磁导率μ的变化，另外，角频率ω较低，导致位移电流也忽略不计，则麦克斯韦方程组可表示为：

▽×E＝jωB (1)

▽×B＝μσE+μJ_s (2)

其中，▽为汉密尔顿算子，E为电场强度矢量，B为磁感应强度矢量，σ为电导率，J_s为电偶极源的传导电流密度矢量；磁矢量势A与磁感应强度矢量B之间的关系为：

B＝▽×A (3)

由式(1)-(3)，电场强度矢量E可表示为：

通过汉克尔变换，可计算矢量势A：

其中，λ为信号波长，J₀为第一类零阶贝塞尔函数，为经过傅里叶变换后的矢量势，r为水平距离矢量，r为r的模；在MCSEM勘探中，经常使用水平电偶极子作为发射源。若电偶极源沿y轴方向放置，则矢量势具有如下形式：

其中，是满足的中间变量，且矢量势和是柱对称；当将电偶极子放置于第i层、深度坐标为z_s的位置处时，则第m层、深度坐标为z的位置处，其对应的矢量势分量和分别为：

其中，δ_mi为克罗内克函数，a_m、b_m、c_m、d_m是相关的系数项；

对于发射源所在层以上的各层，其横电模式的反射系数和横磁模式的反射系数分别定义为：

对于发射源所在层以下的各层，其横电模式的反射系数和横磁模式的反射系数分别定义为：

根据电磁场边界条件可得：

其中，h_m＝z_m+1-z_m，和对应γ_m+1、h_m+1，和对应γ_m-1、h_m-1；另外，和的表达式分别如下：

其中，和对应γ_m+1、σ_m+1，和对应γ_m-1、σ_m-1；

显然，顶层和底层处的反射系数为0，即有：利用电磁场的边界条件，进行递推即可求出a_m，b_m，c_m，d_m的值；然后，将a_m，b_m，c_m，d_m的值代入式(7)和式(8)，可计算出矢量势分量；最后，将式(5)代入式(3)和式(4)，即可分别计算出磁感应强度矢量B和电场强度矢量E。以上是以电偶极源沿y轴方向放置为例，最终得出的相关关系式。若电偶极源沿x轴方向放置，则推导过程类似，可得到一系列类似的关系式；当电偶极源在水平面沿任意方向放置时，则可通过叠加原理进行计算。

在上述步骤中，MT噪声一维正演为：

根据吉洪诺夫-卡尼亚大地电磁理论，假设天然电磁场源为垂直入射到海表面的平面电磁波，利用图2所示的一维模型，可推导得出MT信号的电磁场表达式。理想情况下，MCSEM接收机中电通道、磁通道传感器的摆放应与发射机的方位保持对应。在实际工程应用中，接收机一般从海面被投放后自由沉底，在其着底后传感器的方位存在随机性。为简便起见，本文不考虑这种随机方位的影响。因此，所采用的坐标系与研究MCSEM正演时所使用的坐标系保持一致。接下来，为了推导过程中符号表示的简便，同时又不失一般性，假设MCSEM发射机的摆放使得电场强度矢量的方向位于x轴，磁场强度矢量的方向位于y轴。上述假设下，电场强度矢量将只有x分量，磁场强度矢量将只有y分量。

假设射入海表面的平面波以e^-jωt规律进行变化，另外，设入射到海表面的磁场强度为H，即z＝0时，H_y(0)＝H，这里的常数H实际上表征了天然电磁场活动的强烈程度，另外，为了体现MT噪声随时间、空间随机变化的特点，H应满足实测MT资料的统计特性；根据边界条件以及电场和磁场的关系，第m层、深度坐标为z处的电磁场分量为：

其中Am，Bm为待定系数，为波数；

利用数学递推方法和边界条件，可得：

其中，RN＝1，和R_m的表达式为：

将式(19)、式(20)代入式(17)、式(18)，便可计算出MT信号的电磁场表达式。

上述步骤中，一维OCCAM反演模型为：

OCCAM反演是一种带平滑约束且稳定收敛的最小二乘方法。该方法在搜索与正演数据最大程度拟合的模型时，要求模型光滑，并在可以拟合的多个模型中选择最光滑的模型。因该方法受初始模型影响较小，在地球物理学中有着广泛的应用。

OCCAM反演的实质是将无约束目标函数最小化，定义目标函数φ为：

φ＝||Rm||²+μ^-1(||Wd-WF(m)²-X_* ²) (23)

其中，R为一阶差分算子矩阵，m为模型矢量，μ^-1为拉格朗日乘子，W为权重矩阵，d为观测数据矢量，F(m)为模型矢量m对应的正演数据矢量，X_* ²为反演的目标拟合差；R和W的表达式分别为：

其中，σ₁,σ₂,…,σ_M为反演数据的标准差，M是观测数据的个数；

求取目标函数φ的梯度，使φ达到最小值，并令▽_mφ＝0，此时模型矢量m应满足：

μ^-1(WJ)^TWJm-μ^-1(WJ)^TWd+R^TRm＝0 (26)

其中，J＝▽_mF为雅可比矩阵，其元素可表示为：

其中，i＝1,...,M，j＝1,...,l，l为模型层数；

在求解式(26)时，需将F进行线性泰勒展开，即：

F(m+Δ)≈F(m)+JΔ (28)

其中，Δ表示迭代过程中模型的增量；

给定初始模型矢量m₁，并将式(28)代入式(23)，使φ最小化，得到k次迭代后的模型矢量：

m_k+1(u)＝[μR^TR+(WJ_k)^TWJ_k]^-1(WJ_k)^TW(d-F[m_k]+J_km_k) (29)

通过线性搜索μ，寻找F[m_k+1]与实测数据的最佳拟合差：

X_k+1(μ)＝||W(d-F[m_k+1(μ)])|| (30)

初始迭代阶段的拟合差可能大于目标拟合差，经过式(29)不断迭代得到满足精度的模型；在拟合精度满足要求的情况下，选择最大的μ对应的模型，因为此时对应的模型最光滑；在迭代运算过程中，如果均方根拟合差x_rms达到1时，则OCCAM反演结果收敛；其定义为：

其中，s_i为第i列数据的标准差。

本发明所提出的MT噪声降噪效果评估方法，既支持通过正演的手段生成MT噪声，也可以输入实测数据作为MT噪声。因此，设计了两个实验，分别针对仿真MT噪声和实测MT噪声开展降噪效果评估。

如图3所示，给出了生成仿真MT噪声的流程图。图3的流程是方法步骤中的步骤2至步骤4部分的细化，利用该图所示流程可生成通过正演手段形成的MT噪声数据文件。

为提高MT噪声的逼真程度，还可以直接输入MCSEM发射机不开机时间段的实测数据作为MT噪声。图4给出了MCSEM正演流程图，其是方法步骤中的步骤6至步骤10部分的细化，利用图4所示流程可生成MCSEM正演数据文件。值得注意的是，当采用实测数据生成MT噪声时，MCSEM信号正演中所需要的地层电阻率、层数、各层深度等参数不能随意设置，而是应根据MCSEM发射机开机时间段实测数据的OCCAM反演结果进行设置。分别得到MT噪声数据文件和MCSEM正演数据文件以后，如图1所示，可分别生成无噪MCSEM信号和含噪MCSEM信号的MVO曲线，从而可定性观察MT噪声对MVO曲线的影响。

得到了含噪MCSEM数据以后，按图5所示流程进行短时窗结合预白化滤波的降噪处理。图5的流程是方法步骤中的步骤14和步骤15部分的细化，利用该图所示流程可得到MT降噪处理以后的MCSEM数据文件，并生成降噪后的MVO曲线，从而可定性观察降噪处理后MVO曲线的改善。

为定量分析降噪效果，针对含噪MCSEM数据、降噪后MCSEM数据，分别按图6所示流程进行OCCAM反演。

图6流程是方法步骤中的步骤17部分的细化，利用该图所示流程可得到降噪前、后的OCCAM反演结果。接下来，如图1流程图所示，分别画出真实电阻率与降噪前、后反演电阻率随深度变化的曲线，分别计算降噪前、后反演电阻率与真实电阻率之间的均方相对误差，从而给出MT噪声降噪效果的量化评价。

实施例1

设置空气层、海水层、海底地层的电阻率、深度如图7所示。设定发射信号波形为周期为10s的方波，电偶源长度为250m，发射机位于海底上方40m处，平台航行速度1m/s，发射机初始坐标(0m,0m,960m)，接收机数量1个、位置位于海底、位置坐标为(0m,0m,1000m)。

将上述参数带入本发明的评估方法步骤内，分别得到MT噪声数据文件和MCSEM数据文件，相加后得到含噪MCSEM数据。图8为无噪MCSEM信号电通道MVO曲线图，图9为无噪MCSEM信号磁通道的MVO曲线图，图10为含噪MCSEM信号电通道的MVO曲线图，图11为含噪MCSEM信号磁通道的MVO曲线图。将图8与图10对比，图9与图11对比，可发现，MT噪声在中、远收发距处对MCSEM信号产生了较强影响，使MVO曲线产生较大震荡，降低了MCSEM数据的信噪比，将严重影响数据解释质量。

接下来，利用短时窗结合预白化滤波的方法，对含有MT噪声的MCSEM数据进行了降噪处理。图12为降噪前、后含噪MCSEM信号电通道的MVO曲线图与无噪MCSEM信号电通道MVO曲线图，图13为降噪前、后含噪MCSEM信号磁通道的MVO曲线图与无噪MCSEM信号磁通道MVO曲线图。从图12和图13可以看出，中、远收发距处的MT噪声水平得到降低，降噪后的MVO曲线相对更接近无噪数据的MVO曲线。

最后，将降噪前与降噪后的MCSEM数据分别进行了一维OCCAM反演。图14为MCSEM信号电通道降噪前真实电阻率与反演电阻率的变化曲线，图15为MCSEM信号电通道降噪后真实电阻率与反演电阻率的变化曲线。对比图14和图15发现，经降噪处理后，所得到的反演电阻率曲线明显更接近于真实电阻率的变化。

实施例2

为提高MT噪声的逼真程度，直接输入发射机不开机时间段的实测数据作为MT噪声。选择了赤道几内亚某海域2009年12月的某次实验，利用其MCSEM发射机开机时间段的实测数据开展了一维OCCAM反演，所反演得到的信息如图16所示。根据该次实验相关配置文件中的参数，设定发射信号波形为周期为10s的方波，电偶源长度为300m，发射机位于海底上方50m处，平台航行速度1m/s，发射机初始坐标(0m,0m,850m)，选取了其中1个位置坐标为(0m,0m,900m)的接收站点、将其MCSEM发射机不开机时间段的部分实测数据作为MT噪声。

利用上述输入参数，正演得到MCSEM数据文件，与从实测数据提取得到的MT噪声文件相加后得到含噪MCSEM数据。图17为无噪MCSEM信号电通道的MVO曲线图，图18为无噪MCSEM信号磁通道的MVO曲线图。图19为含噪MCSEM信号电通道的MVO曲线图，图20为含噪MCSEM信号磁通道的MVO曲线图。

将图17与图19对比，图18与图20对比，MT噪声在中、远收发距处对MCSEM信号产生了较强影响，使MVO曲线产生较大震荡，降低了MCSEM数据的信噪比，将严重影响数据解释质量。

接下来，利用短时窗结合预白化滤波的方法，对含有MT噪声的MCSEM数据进行了降噪处理。图21为降噪前、后含噪MCSEM信号电通道的MVO曲线图与无噪MCSEM信号电通道MVO曲线图，图22为降噪前、后含噪MCSEM信号磁通道的MVO曲线图与无噪MCSEM信号磁通道MVO曲线图。从图21和图22可以看出，中、远收发距处的MT噪声水平得到降低，降噪后的MVO曲线相对更接近无噪数据的MVO曲线。

最后，将降噪前与降噪后的MCSEM数据分别进行了一维OCCAM反演。图23为MCSEM信号电通道降噪前真实电阻率与反演电阻率的变化曲线，图24为MCSEM信号电通道降噪后真实电阻率与反演电阻率的变化曲线。对比图23和图24，经降噪处理后，所得到的反演电阻率曲线明显更接近于真实电阻率的变化。为给出针对实测数据的量化降噪效果，计算了反演电阻率与真实电阻率之间的均方相对误差。其计算公式为：

其中，σ_rel表示均方相对误差，L为测量点数，ρ_real,i为第i次测量的真实电阻率，ρ_inv,i为第i次测量的反演电阻率。针对所选取的实测数据，表1给出了降噪前、后的反演电阻率与真实电阻率之间的均方相对误差。从该表中可看出，降噪后反演数据的均方相对误差为231％，降噪前反演数据的均方相对误差为357％。相比于不进行降噪处理的直接反演结果，均方相对误差减小了126％。由此可见，所采用的短时窗结合预白化滤波的MT噪声降噪方法确实取得了比较明显的降噪效果，因而也有力地证明了该降噪处理方法的有效性。

表1实施例2中反演电阻率与真实电阻率之间的均方相对误差

本发明提出的海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，利用电阻率反演结果与电阻率输入值的均方相对误差作为降噪效果的量化评价指标。为提高灵活性，本发明既可以利用正演的手段生成MT噪声，也可以直接输入MCSEM发射机不开机时间段的实测数据作为MT噪声。利用该评估方法对实测数据进行处理，验证了短时窗结合预白化滤波降噪方法的有效性。本发明既可应用于MT噪声降噪方法有效性的量化评估，也可推广至其它类型噪声降噪效果的量化评估。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (4)

1.海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，其特征在于，包括以下步骤：

步骤1：根据需求的不同，判断是否采用海洋可控源电磁勘探发射机不开机时间段的实测数据作为MT噪声，若不采用，则进行步骤2；若采用，则进行步骤6；

步骤2：输入地层层数、深度、电阻率参数，参数用于MCSEM信号一维正演和MT噪声一维正演；

步骤3：设置用于MT噪声一维正演的海表面磁场强度、接收站个数和坐标、海水层深度和电阻率输入参数，其中接收站个数和坐标、海水层深度和电阻率参数也用于MCSEM信号一维正演；

步骤4：利用MT噪声一维正演模型，生成MT噪声数据文件；

步骤5：设置用于MCSEM信号一维正演的发射源长度、初始坐标、航行速度、发射信号周期输入参数，之后，进行步骤10；

步骤6：利用发射机开机时间段的MCSEM实测数据，采用一维OCCAM反演模型，反演得到地层层数、深度、电阻率参数；

步骤7：读取MCSEM实测数据的相关配置文件，保存用于MT噪声一维正演的接收站个数和坐标、海水层深度和电阻率参数，这些参数也用于MCSEM信号一维正演；

步骤8：选取MCSEM发射机不开机时间段的部分实测数据作为MT噪声数据文件，其数据点数应与MCSEM数据文件的数据点数相同；

步骤9：读取MCSEM实测数据的相关配置文件，保存用于MCSEM正演的发射源长度、初始坐标、航行速度、发射信号周期参数；

步骤10：利用MCSEM信号一维正演模型，生成MCSEM正演数据文件；

步骤11：根据MCSEM正演数据文件，生成无噪数据的MVO曲线；

步骤12：将MCSEM正演数据文件中的电磁场数据与MT噪声数据文件中的电磁场数据相加，生成含噪MCSEM数据文件；

步骤13：根据含噪MCSEM数据文件，生成含噪数据的MVO曲线；

步骤14：利用短时窗结合预白化滤波MT噪声降噪方法，对含噪MCSEM数据进行MT噪声降噪处理；

步骤15：根据降噪后的MCSEM数据，生成降噪后的MVO曲线；

步骤16：对比降噪前、后的MVO曲线，观察降噪效果，给出定性分析；

步骤17：利用一维OCCAM反演模型，对降噪前、后的MCSEM数据分别进行反演，得到反演后的地层层数、深度、电阻率结果；

步骤18：分别画出真实电阻率与降噪前、后反演电阻率随深度变化的曲线，观察降噪处理对反演精度的提高，给出定性分析；

步骤19：分别计算降噪前、后反演电阻率与真实电阻率之间的均方相对误差；

步骤20：利用降噪前、后反演电阻率与真实电阻率之间的均方相对误差，对MT噪声的降噪效果给出最终的量化评价。

2.根据权利要求1所述的海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，其特征在于，MCSEM信号一维正演模型为：

设立N层一维模型，其中，x轴和y轴位于水平面，z轴垂直向下，xyz坐标系满足右手定则，σ_m为第m层的电导率且为各向同性，z_m为第m层层顶的z坐标，h_m为第m层的深度，N为总共的层数；第1层总为空气层，z₁为-∞，将h₁设置为∞，第2层为海水，将z₂设置为0，第N层的深度h_N设置为∞；

假设时间谐变偶极子源的时谐因子为e^-jωt，忽略磁导率μ的变化，另外，角频率ω较低，导致位移电流也忽略不计，则麦克斯韦方程组可表示为：

▽×E＝jωB (1)

▽×B＝μσE+μJ_s (2)

其中，▽为汉密尔顿算子，E为电场强度矢量，B为磁感应强度矢量，σ为电导率，J_s为电偶极源的传导电流密度矢量；磁矢量势A与磁感应强度矢量B之间的关系为：

B＝▽×A (3)

由式(1)-(3)，电场强度矢量E可表示为：

通过汉克尔变换，可计算矢量势A：

其中，λ为信号波长，J₀为第一类零阶贝塞尔函数，为经过傅里叶变换后的矢量势，r为水平距离矢量，r为r的模；在MCSEM勘探中，经常使用水平电偶极子作为发射源；若电偶极源沿y轴方向放置，则矢量势具有如下形式：

其中，是满足的中间变量，且矢量势和是柱对称；当将电偶极子放置于第i层、深度坐标为z_s的位置处时，则第m层、深度坐标为z的位置处，其对应的矢量势分量和分别为：

其中，δ_mi为克罗内克函数，a_m、b_m、c_m、d_m是相关的系数项；

对于发射源所在层以上的各层，其横电模式的反射系数和横磁模式的反射系数分别定义为：

对于发射源所在层以下的各层，其横电模式的反射系数和横磁模式的反射系数分别定义为：

根据电磁场边界条件可得：

其中，h_m＝z_m+1-z_m，和对应γ_m+1、h_m+1，和对应γ_m-1、h_m-1；另外，和的表达式分别如下：

其中，和对应γ_m+1、σ_m+1，和对应γ_m-1、σ_m-1；

显然，顶层和底层处的反射系数为0，即有：利用电磁场的边界条件，进行递推即可求出a_m，b_m，c_m，d_m的值；然后，将a_m，b_m，c_m，d_m的值代入式(7)和式(8)，可计算出矢量势分量；最后，将式(5)代入式(3)和式(4)，即可分别计算出磁感应强度矢量B和电场强度矢量E；以上是以电偶极源沿y轴方向放置为例，最终得出的相关关系式；若电偶极源沿x轴方向放置，则推导过程类似，可得到一系列类似的关系式；当电偶极源在水平面沿任意方向放置时，则可通过叠加原理进行计算。

3.根据权利要求2所述的海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，其特征在于，MT噪声一维正演为：

MT噪声一维正演所采用的坐标系与MCSEM信号一维正演模型所使用的坐标系保持一致；接下来，为了推导过程中符号表示的简便，同时又不失一般性，假设MCSEM发射机的摆放使得电场强度矢量的方向位于x轴，磁场强度矢量的方向位于y轴；上述假设下，电场强度矢量将只有x分量，磁场强度矢量将只有y分量；

假设射入海表面的平面波以e^-jωt规律进行变化，另外，设入射到海表面的磁场强度为H，即z＝0时，H_y(0)＝H，这里的常数H实际上表征了天然电磁场活动的强烈程度，另外，为了体现MT噪声随时间、空间随机变化的特点，H应满足实测MT资料的统计特性；根据边界条件以及电场和磁场的关系，第m层、深度坐标为z处的电磁场分量为：

其中A_m，B_m为待定系数，为波数；

利用数学递推方法和边界条件，可得：

其中，R_N＝1，和R_m的表达式为：

将式(19)、式(20)代入式(17)、式(18)，便可计算出MT信号的电磁场表达式。

4.根据权利要求1所述的海洋可控源电磁勘探中大地天然电磁信号噪声降噪效果评估方法，其特征在于，一维OCCAM反演模型为：

OCCAM反演的实质是将无约束目标函数最小化，定义目标函数φ为：

φ＝||Rm||²+μ^-1(||Wd-WF(m)||²-X_* ²) (23)

其中，R为一阶差分算子矩阵，m为模型矢量，μ^-1为拉格朗日乘子，W为权重矩阵，d为观测数据矢量，F(m)为模型矢量m对应的正演数据矢量，X_* ²为反演的目标拟合差；R和W的表达式分别为：

其中，σ₁,σ₂,···,σ_M为反演数据的标准差，M是观测数据的个数；

求取目标函数φ的梯度，使φ达到最小值，并令▽_mφ＝0，此时模型矢量m应满足：

μ^-1(WJ)^TWJm-μ^-1(WJ)^TWd+R^TRm＝0 (26)

其中，J＝▽_mF为雅可比矩阵，其元素可表示为：

其中，i＝1,...,M，j＝1,...,l，l为模型层数；

在求解式(26)时，需将F进行线性泰勒展开，即：

F(m+Δ)≈F(m)+JΔ (28)

其中，Δ表示迭代过程中模型的增量；

给定初始模型矢量m₁，并将式(28)代入式(23)，使φ最小化，得到k次迭代后的模型矢量：

m_k+1(u)＝[μR^TR+(WJ_k)^TWJ_k]^-1(WJ_k)^TW(d-F[m_k]+J_km_k) (29)

通过线性搜索μ，寻找F[m_k+1]与实测数据的最佳拟合差：

X_k+1(μ)＝||W(d-F[m_k+1(μ)])|| (30)

初始迭代阶段的拟合差可能大于目标拟合差，经过式(29)不断迭代得到满足精度的模型；在拟合精度满足要求的情况下，选择最大的μ对应的模型，因为此时对应的模型最光滑；在迭代运算过程中，如果均方根拟合差x_rms达到1时，则OCCAM反演结果收敛；其定义为：

其中，s_i为第i列数据的标准差。