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

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

【技术实现步骤摘要】

本专利技术涉及大地电磁测深方法的处理
，具体地说是一种海洋大地电磁场非线性共轭梯度三维并行反演方法。
技术介绍
大地电磁测深方法广泛应用于深部地质探测和矿产、水、石油、地热等资源勘查。而目前二维反演是其数据处理的主流方法，三维反演方法的研究较少，国外科学家分别提出了快速松弛法、共轭梯度法、非线性共轭梯度法等方法，国内已经实现了快速松弛法和共轭梯度法。目前，对海洋资源的开发已经成为重要的发展方向。现有的大地电磁三维反演均为陆地资料处理方法，只能反演地表(空气与地面界面)测量的大地电磁场数据，不能处理海底等水底大地测量的数据。
技术实现思路
有鉴于此，本专利技术实施例提供一种海洋大地电磁场非线性共轭梯度三维并行反演方法，主要目的是提供一种适用于包含水域数据的大地电磁场反演的方法。为达到上述目的，本专利技术主要提供如下技术方案：一方面，本专利技术实施例提供了一种海洋大地电磁场非线性共轭梯度三维反演方法，包括如下步骤：a.在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成Nx、Ny、Nz个小的长方体网格单元，间距为Δx(i)(i＝1,...,Nx)、Δy(j)(j＝1,...,Ny)、Δz(k)(k＝1,...,Nz)，其中模型空间包括空气层、海水层和地下半空间；海底节点(i，j，k)的长方体网格单元的长度、宽度和高度分别为xseafloor-i、yseafloor-j、zseafloor-k，电阻率为ρseafloor-i,j,k，其中电场取在各网格单元边缘的中点，用和表 示；磁场取在各网格单元表面的中心，用和 表示；b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其计算，实现正演；c.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；e.计算目标函数的梯度，以非线性共轭梯度法作为反演方法获得查找步长和查找方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得的电阻率计算反演方法中的预处理因子；f.使用模型修改量修正模型，并返回步骤b。作为优选，其中根据计算机的CPU数量自动将每个频点的数据分配给所有的CPU进行正演以及目标函数梯度的并行计算。作为优选，节点的场值表达式，Eseafloor-x如下所示：其中seafloor表示海底及地层网格节点，sea表示海水层网格节点。作为优选，使用电场值E，通过麦克斯韦方程可以求出各网格节点的磁场值H，便可以得到模型的阻抗响应Z，作为优选，反演目标函数为：其中，和Zseafloor-n专指海底观测数据和正演数据，εseafloor-n为观测数据误差，第1-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W＝Wair+Wsea+Wseafloor为预调节矩阵，m＝mair+msea+mseafloor为网格剖分模型，air表示空气层部分，sea表示海水层部分，seafloor表示海底及地下半空间部分，目标函数梯度表示为：三维阻抗灵敏度张量表示为：其中，K＝Kair+Ksea+Kseafloor为正演系数矩阵，E＝Eair+Esea+Eseafloor为场值并分为两种源条件下的场值，为正演模型有限差分网格至接收点插入两种极化源电磁场向量的线性组合；得到：p(i)＝-C(i)-1g(i)+β(i)p(i-1) (14)其中γ为模型电阻率和数据误差相关的变量，从而得到更新的模型：mseafloor-(i+1)＝mseafloor-(i)+α(i)p(i) (18)。与现有技术相比，本专利技术的有益效果在于：本专利技术的海洋大地电磁场非线性共轭梯度三维反演方法能够将野外采集到的海底(或湖、河等水底)大地电磁场测深数据通过反演计算转换成地下三维电性结构信息，反映更为真实的地下电性结构信息，根据海底测量的数据反映出海底测区地壳浅表、地壳浅层、地壳深层、壳幔边界、上地幔的电阻率信息。该反演方法在保证计算精度的条件下，计算效率很高，适用于目前国内外十分主流的Windows操作系统和Linux、Unix操作系统，并且可以在普通个人电脑上推广使用。该新技术能够为探索海洋成矿区带不同深度尺度上地球物理特征与成矿之间的关系、形成海洋矿产资源立体探测的技术解决方案和深部资源 勘查提供依据技术支持，指明海洋找矿方向，重塑深部成矿动力学过程，为建立海洋成矿理论提供更可靠的信息。附图说明图1为本专利技术的海洋大地电磁场非线性共轭梯度三维反演方法的流程框图；图2为本专利技术的三维模型的示意图；图3为图2中的单个长方体网格单元示意图；图4为本专利技术方法的应用例1的模型示意图；图5为本专利技术方法的应用例2的模型示意图；图6为本专利技术方法的应用例2的三维反演模型的二维切片图。具体实施方式下面结合附图和具体实施例对本专利技术作进一步详细描述，但不作为对本专利技术的限定。本专利技术中的海洋大地电磁场是指与陆地大地电磁场相对应的，包括水域数据的大地电磁场，如海底、湖底、河底等的大地电磁场。本专利技术实施例中未尽之处，均可从现有技术中获得。图1为本专利技术的海洋大地电磁场非线性共轭梯度三维反演方法的流程框图。如图1所示，海洋大地电磁场非线性共轭梯度三维反演方法，包括如下步骤：a.在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成Nx、Ny、Nz个小的长方体网格单元，间距为Δx(i)(i＝1,...,Nx)、Δy(j)(j＝1,...,Ny)、Δz(k)(k＝1,...,Nz)，其中模型空间包括空气层、海水层和地下半空间；海底节点(i，j，k)的长方体网格单元的长度、宽度和高度分别为xseafloor-i、yseafloor-j、zseafloor-k，电阻率为ρseafloor-i,j,k，其中电场取在各网格单元边缘的中点，用和表示；磁场取在各网格单元表面的中心，用和 表示；b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其计算，实现正演；c.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；e.计算目标函数的梯度，以非线性共轭梯度法作为反演方法获得查找步长和查找方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得的电阻率计算反演方法中的预处理因子；f.使用模型修改量修正模型，并返回步骤b。本专利技术方法将固定电阻率值(海水的电阻率值是根据实际情况设计的，在正、反演过程中保持不变)的海水层加入到初始模型，推导带有海水层模型各节点的场值插值表达式和计算过程，以实现正演计算；推导带有海水层模型各节点的模型修改量(计算过程中主要包括目标函数、目标函数梯度、查找方向和查找步长等变量参数)的计算表达式和求解过程，以实现反演计算；添加基于频点的并行计算结构，以提高效率，最终达到快速获得准确地下电性结构分布的目的。本专利技术方法能够将野外采集到的海底(或湖、河等水底)大地电磁测深数据通过反演计算转换成地下三维电性结构信息，反映更为真实的地下深部电性结构信息，根据测量数据能够反映出海底地壳浅表、地壳浅层、地壳深层、壳幔边界、上地幔的电阻率信息。能够为本文档来自技高网...

【技术保护点】
海洋大地电磁场非线性共轭梯度三维并行反演方法，其特征在于，包括如下步骤：a.在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成Nx、Ny、Nz个小的长方体网格单元，间距为Δx(i)(i＝1,...,Nx)、Δy(j)(j＝1,...,Ny)、Δz(k)(k＝1,...,Nz)，其中模型空间包括空气层、海水层和地下半空间；海底节点(i，j，k)的长方体网格单元的长度、宽度和高度分别为xseafloor‑i、yseafloor‑j、zseafloor‑k，电阻率为ρseafloor‑i,j,k，其中电场取在各网格单元边缘的中点，用和表示；磁场取在各网格单元表面的中心，用和表示；b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其计算，实现正演；c.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；e.计算目标函数的梯度，以非线性共轭梯度法作为反演方法获得查找步长和查找方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得的电阻率计算反演方法中的预处理因子；f.使用模型修改量修正模型，并返回步骤b。...

【技术特征摘要】
1.海洋大地电磁场非线性共轭梯度三维并行反演方法，其特征在于，包括如下步骤：a.在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成Nx、Ny、Nz个小的长方体网格单元，间距为Δx(i)(i＝1,...,Nx)、Δy(j)(j＝1,...,Ny)、Δz(k)(k＝1,...,Nz)，其中模型空间包括空气层、海水层和地下半空间；海底节点(i，j，k)的长方体网格单元的长度、宽度和高度分别为xseafloor-i、yseafloor-j、zseafloor-k，电阻率为ρseafloor-i,j,k，其中电场取在各网格单元边缘的中点，用和表示；磁场取在各网格单元表面的中心，用和表示；b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其计算，实现正演；c.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；e.计算目标函数的梯度，以非线性共轭梯度法作为反演方法获得查找步长和查找方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得的电阻率计算反演方法中的预处理因子；f.使用模型修改量修正模型，并返回步骤b。2.根据权利要求1所述的海洋大地电磁场非线性共轭梯度三维并行反演方法，其特征在于，其中根据计算机的CPU数量自动将每个频点的数据分配给所有的CPU进行正演以及目标函数梯度的并行计算。3.根据权利要求1所述的海洋大地电磁场非线性共轭梯度三维并行反演方法，其特征在于，节点的场值表达式，Eseafloor-x如下所示： iωμ 0 E s e a f l o o r - x i - 1 / 2 , j , k x s e a f l o o r - i - 1 / 2 y s e a f l o o r - j z s e a f l o o r - k / ρ s e a f l o o r - i - 1 / 2 , j , k = ( E s e a f l o o r - y i , j + 1 / 2 , k - E s e a f l o o r - y i - 1 , j + 1 / 2 , k ) z s e a f l o o r - k - ( E s e a f l o o r - x i - 1 / 2 , j + 1 , k - E s e a f l o o r - x i - 1 / 2 , j , k ) z s e a f l o o r - k x s e a f l o o r - i - 1 / 2 / y s e a f l o o r - j + 1 / 2 - ( E s e a f l o o r - y i , j - 1 / 2 , k - E s e a f l o o r - y i - 1 , j - 1 / 2 , k ) z s e a f l o o r - k - ( E s e a f l o o r - x i - 1 / 2 , j , k - E s e a f l o o r - x i - 1 / 2 , j - 1 , k ) z s e a f l o o r - k x s e a f l o o r - i - 1 / 2 / y s e a f l o o r - j - 1 / 2 + ( E s e a f l o o r - z i , j , k + 1 / 2 - E s e a f l o o r - z i - 1 , j , k + 1 / 2 ) y j - ( E s e a - x i - 1 / 2 , j + 1 , k + 1 - E s e a f l o o r - x i - 1 / 2 , j , k ) x s e a f l o o r - i - 1 / 2 y j / z s e a - k + 1 / 2 - ( E s e a - z i , j , k - 1 / 2 - E s e a - z i - 1 , j , k - 1 / 2 ) y j + ( E s e a f l o o r - x i - 1 / 2 , j , k - E s e a - x i - 1 / 2 , j , k - 1 ) y s e a f l o o r - j x s e a f l o o r - i - 1 / 2 / z s e a - k - 1 / 2 ]]>其中seafloor表示海底及地层网格节点，sea表示海水层网格节点。4.根据权利要求1所述的海洋大地电磁场非线性共轭梯度三维并行反演方法，其特征在于，使用电场值E，通过麦克斯韦方程可以求出各网格节点的磁场值H，便可以得到模型的阻抗响应Z： E s e a f l o o r - x ( i , j , k ) T E s e a f l o o r - y ( i , j , k ) T = Z s e a f l o o r - x x ( i , j , k ) Z s e a f l o o r - x y ( i , j , k ) Z s e a f l o o r - y x ( i , j , k ) Z s e a f l o o r - y y ( i , j , k ) H s e a f l o o r - x ( i , j , k ...

【专利技术属性】
技术研发人员：张昆，严加永，蔡德超，董浩，
申请(专利权)人：中国地质科学院矿产资源研究所，
类型：发明
国别省市：北京;11