海洋大地电磁场非线性共轭梯度三维并行反演方法

摘要

本发明公开了一种海洋大地电磁场非线性共轭梯度三维并行反演方法，将固定电阻率值的海水层加入到初始模型，推导带有海水层模型各节点的大地电磁场场值的插值表达式和计算过程，以实现正演计算；推导带有海水层模型各节点的模型修改量(计算过程中主要包括目标函数、目标函数梯度、查找方向和查找步长等变量参数)的计算表达式和求解过程，以实现反演计算。本发明适用于海洋或其他包含水域的大地电磁场非线性共轭梯度三维反演。

说明书

技术领域

本发明涉及大地电磁测深方法的处理技术领域，具体地说是一种海洋大地电磁场非线性共轭梯度三维并行反演方法。

背景技术

大地电磁测深方法广泛应用于深部地质探测和矿产、水、石油、地热等资源勘查。而目前二维反演是其数据处理的主流方法，三维反演方法的研究较少，国外科学家分别提出了快速松弛法、共轭梯度法、非线性共轭梯度法等方法，国内已经实现了快速松弛法和共轭梯度法。

目前，对海洋资源的开发已经成为重要的发展方向。现有的大地电磁三维反演均为陆地资料处理方法，只能反演地表(空气与地面界面)测量的大地电磁场数据，不能处理海底等水底大地测量的数据。

发明内容

有鉴于此，本发明实施例提供一种海洋大地电磁场非线性共轭梯度三维并行反演方法，主要目的是提供一种适用于包含水域数据的大地电磁场反演的方法。

为达到上述目的，本发明主要提供如下技术方案：

一方面，本发明实施例提供了一种海洋大地电磁场非线性共轭梯度三维反演方法，包括如下步骤：

a.在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成N_x、N_y、N_z个小的长方体网格单元，间距为Δx_(i)(i＝1,...,N_x)、Δy_(j)(j＝1,...,N_y)、Δz_(k)(k＝1,...,N_z)，其中模型空间包括空气层、海水层和地下半空间；海底节点(i，j，k)的长方体网格单元的长度、宽度和高度分别为x_seafloor-i、y_seafloor-j、z_seafloor-k，电阻率为ρ_{seafloor-i,j,k}，其中电场取在各网格单元边缘的中点，用和表>

b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其计算，实现正演；

c.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；

d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；

e.计算目标函数的梯度，以非线性共轭梯度法作为反演方法获得查找步长和查找方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得的电阻率计算反演方法中的预处理因子；

f.使用模型修改量修正模型，并返回步骤b。

作为优选，其中根据计算机的CPU数量自动将每个频点的数据分配给所有的CPU进行正演以及目标函数梯度的并行计算。

作为优选，节点的场值表达式，E_seafloor-x如下所示：

其中seafloor表示海底及地层网格节点，sea表示海水层网格节点。

作为优选，使用电场值E，通过麦克斯韦方程可以求出各网格节点的磁场值H，便可以得到模型的阻抗响应Z，

作为优选，反演目标函数为：

其中，和Z_seafloor-n专指海底观测数据和正演数据，ε_seafloor-n为观测数据误差，第1-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W＝W_air+W_sea+W_seafloor为预调节矩阵，m＝m_air+m_sea+m_seafloor为网格剖分模型，air表示空气层部分，sea表示海水层部分，seafloor表示海底及地下半空间部分，

目标函数梯度表示为：

三维阻抗灵敏度张量表示为：

其中，K＝K_air+K_sea+K_seafloor为正演系数矩阵，E＝E_air+E_sea+E_seafloor为场值并分为两种源条件下的场值，

为正演模型有限差分网格至接收点插入两种极化源电磁场向量的线性组合；

得到：

p_(i)＝-C_(i)^-1g_(i)+β_(i)p_(i-1)>

其中γ为模型电阻率和数据误差相关的变量，从而得到更新的模型：

m_{seafloor-(i+1)}＝m_seafloor-(i)+α_(i)p_(i)>

与现有技术相比，本发明的有益效果在于：

本发明的海洋大地电磁场非线性共轭梯度三维反演方法能够将野外采集到的海底(或湖、河等水底)大地电磁场测深数据通过反演计算转换成地下三维电性结构信息，反映更为真实的地下电性结构信息，根据海底测量的数据反映出海底测区地壳浅表、地壳浅层、地壳深层、壳幔边界、上地幔的电阻率信息。该反演方法在保证计算精度的条件下，计算效率很高，适用于目前国内外十分主流的Windows操作系统和Linux、Unix操作系统，并且可以在普通个人电脑上推广使用。该新技术能够为探索海洋成矿区带不同深度尺度上地球物理特征与成矿之间的关系、形成海洋矿产资源立体探测的技术解决方案和深部资源 勘查提供依据技术支持，指明海洋找矿方向，重塑深部成矿动力学过程，为建立海洋成矿理论提供更可靠的信息。

附图说明

图1为本发明的海洋大地电磁场非线性共轭梯度三维反演方法的流程框图；

图2为本发明的三维模型的示意图；

图3为图2中的单个长方体网格单元示意图；

图4为本发明方法的应用例1的模型示意图；

图5为本发明方法的应用例2的模型示意图；

图6为本发明方法的应用例2的三维反演模型的二维切片图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细描述，但不作为对本发明的限定。

本发明中的海洋大地电磁场是指与陆地大地电磁场相对应的，包括水域数据的大地电磁场，如海底、湖底、河底等的大地电磁场。本发明实施例中未尽之处，均可从现有技术中获得。

图1为本发明的海洋大地电磁场非线性共轭梯度三维反演方法的流程框图。如图1所示，海洋大地电磁场非线性共轭梯度三维反演方法，包括如下步骤：

a.在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成N_x、N_y、N_z个小的长方体网格单元，间距为Δx_(i)(i＝1,...,N_x)、Δy_(j)(j＝1,...,N_y)、Δz_(k)(k＝1,...,N_z)，其中模型空间包括空气层、海水层和地下半空间；海底节点(i，j，k)的长方体网格单元的长度、宽度和高度分别为x_seafloor-i、y_seafloor-j、z_seafloor-k，电阻率为ρ_{seafloor-i,j,k}，其中电场取在各网格单元边缘的中点，用和表示；磁场取在各网格单元表面的中心，用和>

b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其计算，实现正演；

c.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；

d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；

e.计算目标函数的梯度，以非线性共轭梯度法作为反演方法获得查找步长和查找方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得的电阻率计算反演方法中的预处理因子；

f.使用模型修改量修正模型，并返回步骤b。

本发明方法将固定电阻率值(海水的电阻率值是根据实际情况设计的，在正、反演过程中保持不变)的海水层加入到初始模型，推导带有海水层模型各节点的场值插值表达式和计算过程，以实现正演计算；推导带有海水层模型各节点的模型修改量(计算过程中主要包括目标函数、目标函数梯度、查找方向和查找步长等变量参数)的计算表达式和求解过程，以实现反演计算；添加基于频点的并行计算结构，以提高效率，最终达到快速获得准确地下电性结构分布的目的。本发明方法能够将野外采集到的海底(或湖、河等水底)大地电磁测深数据通过反演计算转换成地下三维电性结构信息，反映更为真实的地下深部电性结构信息，根据测量数据能够反映出海底地壳浅表、地壳浅层、地壳深层、壳幔边界、上地幔的电阻率信息。能够为探索海洋成矿区带不同深度尺度上地球物理特征与成矿之间的关系、形成海洋矿产资源立体探测的技术解决方案和海底深部资源勘查提供依据技术支持，指明找矿方向，查明成矿系统的结构，重塑深部成矿动力学过程，为建立海洋地区成矿理论提供更可靠的信息。

本发明实施例中建立的三维模型参见图2和图3。在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成N_x、N_y、N_z个小的长>(i)(i＝1,...,N_x)、Δy_(j)(j＝1,...,N_y)、Δz_(k)(k＝1,...,N_z)，其中模型空间包括空气层、海水层和地下半空间；海底节点(i，j，k)的长方体网格单元的长度、宽度和高度分别为x_seafloor-i、y_seafloor-j、z_seafloor-k，电阻率为ρ_{seafloor-i,j,k}，其中电场取在各网格单元边缘的中点，用>

海底大地电磁场三维正演：

设随时间变化的谐波函数为e^-iωt，麦克斯韦方程式组用微分形式表示为：

电场强度E满足矢量方程：

其中，ω是角频率，σ为电导率，μ为磁导率，μ₀为真空磁导率。

通过麦克斯韦方程组和边界条件的离散化，可以得到该节点的场值表达式，E_seafloor-x如下所示：

其中seafloor表示海底及地层网格节点，sea表示海水层网格节点；

使用电场值E，通过麦克斯韦方程可以求出各网格节点的磁场值H，

便可以得到模型的阻抗响应Z。

海底大地电磁非线性共轭梯度三维反演：

海底三维反演的目标函数表示为：

其中和Z_seafloor-n在这里专指海底观测数据和正演数据，ε_seafloor-n为观测数据误差，第1-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W＝W_air+W_sea+W_seafloor为预调节矩阵，m＝m_air+m_sea+m_seafloor为网格剖分模型，air表示空气层部分，sea表示海水层部分，seafloor表示海底及地下半空间部分。

目标函数梯度表示为：

三维阻抗灵敏度张量表示为：

其中，K＝K_air+K_sea+K_seafloor为正演系数矩阵，E＝E_air+E_sea+E_seafloor为场值并分为两种源条件下的场值，

为正演模型有限差分网格至接收点插入两种极化源电磁场向量的线性组合。

得到：

p_(i)＝-C_(i)^-1g_(i)+β_(i)p_(i-1)>

其中γ为模型电阻率和数据误差相关的变量，从而得到更新的模型：

m_{seafloor-(i+1)}＝m_seafloor-(i)+α_(i)p_(i)>

重复上述迭代过程直至达到最小，便可得到最优模型。

在OPENMP并行技术基础上，实现了该反演方法的并行运算算法，将不同频点的数据分派给不同的CPU核心分别计算目标函数的梯度，大幅提高了计算效率

本发明实施例中，

第一步，接收反演参数、初始模型参数和数据及误差等参数；

第二步，根据应用计算机的cpu数量自动将每个频点的数据分配给所有的cpu进行正演和目标函数梯度的并行计算；

第三步，并行计算结束，合并数据，计算查找方向和查找步长；

第四步，计算模型改变量并对模型参数进行更新；

第五步，当目标函数足够小时结束迭代，否则令i＝i+1，进入第二步。

本发明采用将各频点数据分散于不同线程进行并行运算不会影响计算精度。本发明根据计算机线程数量将各频点以及网格节点数据分散传输于不同的线程，由于该并行反演方法是高效低损耗的方法，对内存和CPU的要求较低，因此本发明方法能够适用于目前主流的PC机上，基本上不受硬件条件制约，能够在多种操作系统上推广使用。

并行的非线性共轭梯度反演算法中的并行计算部分主要是通过类似正演模拟的计算求解不同频点对应的目标函数梯度，因此对于各个计算机线程，分别计算对应频点的查找步长和查找方向，并行计算结束后进行参数规整求和即可。

下面设置不同的海洋地下空间模型，使用交错网格有限差分法计算正演响应并模拟反演使用的实测数据，建立加密的反演初始模型，应用本发明方法进行反演。

模型一

图4为本发明方法的应用例1的模型示意图，如图4所示，(网格为5*5*6，测点9个，频点为0.01s、0.1s、1s、10s、100s)反演结果如下：

数据频率为0.01s、0.1s、1s、10s、100s，9个测点位于网格中心，迭代82次后满足终止条件，RMS为0.056，主频为2.93GHz的8线程计算机运行时间为42s，(a)为y、z方向的电阻率断面，中心 600m*1000m范围内电阻率为9Ω·m，其上方网格电阻率为95Ω·m，下方网格电阻率为75Ω·m，左右网格电阻率为91Ω·m；(b)为x、z方向的电阻率断面，中心600m*1000m范围内电阻率为9Ω·m，其上方网格电阻率为99Ω·m，下方网格电阻率为78Ω·m，左右网格电阻率为88Ω·m。可以看出本发明方法的反演收敛速度较快，计算效率和计算精度均有较好表现。

模型二

图5为本发明方法的应用例2的模型示意图，如图5所示，模型包括海水层和地下半空间，海水层电阻率0.3Ω·m，厚度1.8km，地下半空间电阻率100Ω·m，包含2个低阻异常体和2个高阻异常体，异常体电阻率分别为10Ω·m和1000Ω·m，异常体顶深(异常体顶界到海底的距离)1.8km，长10km，宽10km、厚3km、间距4km，数据频率范围0.1-1000s，共40个频点，网格数为30*30*40，400个测点位于网格中心。在本应用例中迭代60次后满足终止条件，RMS为0.2，主频为2.63GHz的4核计算机运行时间为2小时，由结果可知(图6所示)，两个低阻异常体和两个高阻异常体均能够较好的分辨出来，异常体中心电阻率分别为9Ω·m和300Ω·m。

以上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或等同替换，这种修改或等同替换也应视为落在本发明的保护范围内。