海洋大地电磁场非线性共轭梯度三维并行反演方法
【专利摘要】本发明公开了一种海洋大地电磁场非线性共轭梯度三维并行反演方法，将固定电阻率值的海水层加入到初始模型，推导带有海水层模型各节点的大地电磁场场值的插值表达式和计算过程，以实现正演计算；推导带有海水层模型各节点的模型修改量(计算过程中主要包括目标函数、目标函数梯度、查找方向和查找步长等变量参数)的计算表达式和求解过程，以实现反演计算。本发明适用于海洋或其他包含水域的大地电磁场非线性共轭梯度三维反演。
【专利说明】
海洋大地电磁场非线性共辆梯度H维并行反演方法
技术领域
[0001] 本发明设及大地电磁测深方法的处理技术领域，具体地说是一种海洋大地电磁场 非线性共辆梯度=维并行反演方法。
【背景技术】
[0002] 大地电磁测深方法广泛应用于深部地质探测和矿产、水、石油、地热等资源勘查。 而目前二维反演是其数据处理的主流方法，=维反演方法的研究较少，国外科学家分别提 出了快速松弛法、共辆梯度法、非线性共辆梯度法等方法，国内已经实现了快速松弛法和共 辆梯度法。
[0003] 目前，对海洋资源的开发已经成为重要的发展方向。现有的大地电磁=维反演均 为陆地资料处理方法，只能反演地表(空气与地面界面)测量的大地电磁场数据，不能处理 海底等水底大地测量的数据。
【发明内容】
[0004] 有鉴于此，本发明实施例提供一种海洋大地电磁场非线性共辆梯度=维并行反演 方法，主要目的是提供一种适用于包含水域数据的大地电磁场反演的方法。
[0005] 为达到上述目的，本发明主要提供如下技术方案：
[0006] -方面，本发明实施例提供了一种海洋大地电磁场非线性共辆梯度=维反演方 法，包括如下步骤：
[0007] a.在笛卡尔坐标系下沿x、y、zS个坐标轴将模型空间划分成Nx、Ny、Nz个小的长方 体网格单元，间距为 Ax(i)(i = l, . . .，Nx)、Ay(j)(j = l, . . . ,Ny)、Az(k)化=1，. . .，Nz)，其中 模型空间包括空气层、海水层和地下半空间；海底节点（i，j，k)的长方体网格单元的长度、 宽度和高度分别为^63|1。。1-：1、7363|1。。^心63|1。。1-1<，电阻率为口363|1。。1-：1,川，其中电场取在各网 格单兀边缘的中点，用革廊。心咖。r-化和怎!脚。表不；磁场取在各网格单 兀表面的中屯、，用巧w?/foor-、,一…：、巧W如。:和马础。《-年1任,.-1。,4表不；
[000引b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其 计算，实现正演；
[0009] C.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；
[0010] d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值 足够小时结束迭代，否则进入下一步；
[0011] e.计算目标函数的梯度，W非线性共辆梯度法作为反演方法获得查找步长和查找 方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得 的电阻率计算反演方法中的预处理因子；
[0012] f.使用模型修改量修正模型，并返回步骤b。
[0013] 作为优选，其中根据计算机的CPU数量自动将每个频点的数据分配给所有的CPU进 行正演W及目标函数梯度的并行计算。
[0014] 作为优选，节点的场值表达式，EseaflDDr-X如下所示：
[0015]
[0016] 其中seafloor表示海底及地层网格节点，sea表示海水层网格节点。
[0017] 作为优选，使用电场值E，通过麦克斯韦方程可W求出各网格节点的磁场值H，便可 W得到模型的阻抗响应Z，
[001 引 （5)
[0019]
[0020]
[0021] (7)
[0022] 其中，式娜。。?"^0和Zseafloor-Ii 专指海底观测数据和正演数据，Eseaf iDDr-n为观测数据误 差，第I-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W = Wair+Wsea+ Wseaflwr为预调节矩阵，m = mair+msea+mseaflQQr为网格剖分模型，air表示空气层部分，sea表示 海水层部分，seafloor表示海底及地下半空间部分，
[0023] 目标函数梯度表示为：
[0027]=维阻抗灵敏度张量表示为：[002引
[0024]
[0025] 巧>
[0026]
(11)
[00巧]其中，K = Kair+Ksea+Kseafloor为正演系数矩阵，E = Eair+Esea+Eseafloor为场值并分为两 种源条件下的场值，
[0030]
为正演模型有限差分网 格至接收点插入两种极化源电磁场向量的线性组合；
[003。得到；
[0032]
(12)
[(
[(
[(
[(
[(
[0038] 其中丫为模型电阻率和数据误差相关的变量，从而得到更新的模型：
[0039] Hlseaf loor-(i+：L)_Hlseafloor-W+口 (i)P(i) ( 18 ) 〇
[0040] 与现有技术相比，本发明的有益效果在于：
[0041] 本发明的海洋大地电磁场非线性共辆梯度=维反演方法能够将野外采集到的海 底(或湖、河等水底)大地电磁场测深数据通过反演计算转换成地下=维电性结构信息，反 映更为真实的地下电性结构信息，根据海底测量的数据反映出海底测区地壳浅表、地壳浅 层、地壳深层、壳慢边界、上地慢的电阻率信息。该反演方法在保证计算精度的条件下，计算 效率很高，适用于目前国内外十分主流的Windows操作系统和Linux、化ix操作系统，并且可 W在普通个人电脑上推广使用。该新技术能够为探索海洋成矿区带不同深度尺度上地球物 理特征与成矿之间的关系、形成海洋矿产资源立体探测的技术解决方案和深部资源勘查 提供依据技术支持，指明海洋找矿方向，重塑深部成矿动力学过程，为建立海洋成矿理论提 供更可靠的信息。
【附图说明】
[0042] 图1为本发明的海洋大地电磁场非线性共辆梯度=维反演方法的流程框图；
[0043] 图2为本发明的S维模型的示意图；
[0044] 图3为图2中的单个长方体网格单元示意图；
[0045] 图4为本发明方法的应用例1的模型示意图；
[0046] 图5为本发明方法的应用例2的模型示意图；
[0047] 图6为本发明方法的应用例2的S维反演模型的二维切片图。
【具体实施方式】
[004引下面结合附图和具体实施例对本发明作进一步详细描述，但不作为对本发明的限 定。
[0049] 本发明中的海洋大地电磁场是指与陆地大地电磁场相对应的，包括水域数据的大 地电磁场，如海底、湖底、河底等的大地电磁场。本发明实施例中未尽之处，均可从现有技术 中获得。
[0050] 图1为本发明的海洋大地电磁场非线性共辆梯度=维反演方法的流程框图。如图1 所示，海洋大地电磁场非线性共辆梯度=维反演方法，包括如下步骤：
[0051] a.在笛卡尔坐标系下沿x、y、zS个坐标轴将模型空间划分成Nx、Ny、Nz个小的长方 体网格单元，间距为 Ax(i)(i = l, . . .，Nx)、Ay(j)(j = l, . . . ,Ny)、Az(k)化=1，. . .，Nz)，其中 模型空间包括空气层、海水层和地下半空间；海底节点（i，j，k)的长方体网格单元的长度、 宽度和局度分力U为Xseaf loor-i、yseaf loor-j、Zseaf loor-k ,电阻率为Pseaf loor-i, j , k ,其中电场取在各网 格单兀边缘的中点，用&蝴。.《'-、,_1..2.,..4.、&?卿。|--1',.,_1;.4和&蝴。。》-_-,,,_1:表不；磁场取在各网格单 兀表面的中心、，)? 如。9^-我针轉。化-,和和震'S碑扭佛'-年.,也y化4表/J、;
[0052] b.通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其 计算，实现正演；
[0053] C.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；
[0054] d.通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值 足够小时结束迭代，否则进入下一步；
[0055] e .计算目标函数的梯度，W非线性共辆梯度法作为反演方法获得查找步长和查找 方向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得 的电阻率计算反演方法中的预处理因子；
[0056] f.使用模型修改量修正模型，并返回步骤b。
[0057] 本发明方法将固定电阻率值(海水的电阻率值是根据实际情况设计的，在正、反演 过程中保持不变)的海水层加入到初始模型，推导带有海水层模型各节点的场值插值表达 式和计算过程，W实现正演计算;推导带有海水层模型各节点的模型修改量(计算过程中主 要包括目标函数、目标函数梯度、查找方向和查找步长等变量参数)的计算表达式和求解过 程，W实现反演计算;添加基于频点的并行计算结构，W提高效率，最终达到快速获得准确 地下电性结构分布的目的。本发明方法能够将野外采集到的海底(或湖、河等水底)大地电 磁测深数据通过反演计算转换成地下=维电性结构信息，反映更为真实的地下深部电性结 构信息，根据测量数据能够反映出海底地壳浅表、地壳浅层、地壳深层、壳慢边界、上地慢的 电阻率信息。能够为探索海洋成矿区带不同深度尺度上地球物理特征与成矿之间的关系、 形成海洋矿产资源立体探测的技术解决方案和海底深部资源勘查提供依据技术支持，指明 找矿方向，查明成矿系统的结构，重塑深部成矿动力学过程，为建立海洋地区成矿理论提供 更可靠的信息。
[005引本发明实施例中建立的S维模型参见图2和图3。在笛卡尔坐标系下沿x、y、zS个 坐标轴将模型空间划分成Nx、Ny、Nz个小的长方体网格单元，间距为Ax(I)Q = I, ...，Nx)、 Ay(j)(j = l,. ..，Ny)、Az(k)化=1，. ..，Nz)，其中模型空间包括空气层、海水层和地下半空 间；海底点（i , j , k)的长方体网格单兀的长度、宽度和局度分别为Xseaf loor-i、Yseaf loor-j、 ZseaflDDr-k，电阻率为PseaflDDr-i, j,k，其中电场取在各网格单元边缘的中点，用 、每邮。。和A姐脚。!'-Z,.w4.,s表不；磁场取在各网格单兀表面的中心，用 巧《1脚。''-、,.,1。_,:、巧《/脚。1'-.1',_,,:.,,;_,,:和馬>却。。,'-心;,,_,。表不。图3为图2中的单个长方体网格单兀 示意图，图3所示的是编号为(i，j，k)的长方体网格单元。
[0化9]海底大地电磁场=维正演：
[0060]设随时间变化的谐波函数为eTi"t，麦克斯韦方程式组用微分形式表示为：
(1)
[0061 ]电场强度E满足矢量方程：
[0062] (2)
[0063]
[0064] 其中：
，O是角频率，O为电导率，y为磁导率，趴为真空磁导率。
[0065] 通过麦克斯韦方程组和边界条件的离散化，可W得到该节点的场值表达式， Eseaf Ioor-X如下所不：
[0066]
(4)
[0067] 其中seafloor表示海底及地层网格节点，sea表示海水层网格节点；
[0068] 使用电场值E，通过麦克斯韦方程可W求出各网格节点的磁场值H，
[0072] 海底大地电磁非线性共辆梯度=维反演：[0073] 海底=维反演的目标函数表示为：
[0069] 俩可Pi得到標巧的阳杭晌巧Z。
[0070； (5)
[0071;
[0074]
(7)
[007 引其中.Zw<?floor-B^CZseafloor-n 在运里专指海底观测数据和正演数据，EseaflDDr-n为观测数 据误差，第I-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W = Wair+ Wsea+Wseafloor为预调节矩阵，m = mair+msea+mseafloor为网格剖分模型，air表示空气层部分，sea 表示海水层部分，seafloor表示海底及地下半空间部分。
[0076]目标函数梯度表示为：
[0080] =维阻抗灵敏度张量表示为：[0081]
[0077；
[0078; C9)：
[0079；
(11)
[0082] 其中，K = Kair+Ksea+Kseaf Ioor为正演系数矩阵，E = Eair+Esea+Eseaf Ioor为场值并分为两 种源条件下的场值，
[0083]
弓正演模型有限差分网 格至接收点插入两种极化源电磁场向量的线性组合。
[0084] 得到；
[0085]
C12)
[0091] 其中丫为模型电阻率和数据误差相关的变量，从而得到更新的模型：
[0092] niseafloor-(i+：L)二mseafloor-(:L)+口(i)P(i) (18)
[0093] 重复上述迭代过程直至梦达到最小，便可得到最优模型。
[0094] 在OPENMP并行技术基础上，实现了该反演方法的并行运算算法，将不同频点的数 据分派给不同的CPU核屯、分别计算目标函数的梯度，大幅提高了计算效率
[0095]
(19)
[0096] 本发明实施例中，
[0097] 第一步，接收反演参数、初始模型参数和数据及误差等参数；
[0098] 第二步，根据应用计算机的CPU数量自动将每个频点的数据分配给所有的CPU进行 正演和目标函数梯度的并行计算；
[0099 ]第=步，并行计算结束，合并数据，计算查找方向和查找步长；
[0100] 第四步，计算模型改变量并对模型参数进行更新；
[0101] 第五步，当目标函数I巧y I足够小时（|巧。| < 6件），结束迭代，否则令i = i + l，进入第 二步。
[0102] 本发明采用将各频点数据分散于不同线程进行并行运算不会影响计算精度。本发 明根据计算机线程数量将各频点W及网格节点数据分散传输于不同的线程，由于该并行反 演方法是高效低损耗的方法，对内存和CPU的要求较低，因此本发明方法能够适用于目前主 流的PC机上，基本上不受硬件条件制约，能够在多种操作系统上推广使用。
[0103] 并行的非线性共辆梯度反演算法中的并行计算部分主要是通过类似正演模拟的 计算求解不同频点对应的目标函数梯度，因此对于各个计算机线程，分别计算对应频点的 查找步长和查找方向，并行计算结束后进行参数规整求和即可。
[0104] 下面设置不同的海洋地下空间模型，使用交错网格有限差分法计算正演响应并模 拟反演使用的实测数据，建立加密的反演初始模型，应用本发明方法进行反演。
[0105] 模型一
[0106] 图4为本发明方法的应用例1的模型示意图，如图4所示，（网格为5*5*6,测点9个， 频点为〇.〇13、0.13、13、1〇3、10〇3)反演结果如下：
[0107]
[010引数据频率为0.01s、0.1 s、Is、10s、100s，9个测点位于网格中屯、，迭代82次后满足终 止条件，RMS为0.056，主频为2.93細Z的8线程计算机运行时间为42s，（a)为y、z方向的电阻 率断面，中屯、600m*1000m范围内电阻率为9Q -m,其上方网格电阻率为95Q -m,下方网格 电阻率为75 Q .m,左右网格电阻率为91 Q *m;(b)为x、z方向的电阻率断面，中屯、600m* 1000m范围内电阻率为9Q .m,其上方网格电阻率为99Q .m,下方网格电阻率为78Q .m, 左右网格电阻率为88 Q ? m。可W看出本发明方法的反演收敛速度较快，计算效率和计算精 度均有较好表现。
[0109] 模型二
[0110] 图5为本发明方法的应用例2的模型示意图，如图5所示，模型包括海水层和地下半 空间，海水层电阻率0.3Q -m，厚度1.8km，地下半空间电阻率IOOQ -m，包含2个低阻异常 体和2个高阻异常体，异常体电阻率分别为IOQ -m和1000 Q -m，异常体顶深(异常体顶界 到海底的距离）1.8虹1，长10km，宽10km、厚3虹1、间距4虹1，数据频率范围O.l-lOOOs，共40个频 点，网格数为30*30*40,400个测点位于网格中屯、。在本应用例中迭代60次后满足终止条件， RMS为0.2,主频为2.63G化的4核计算机运行时间为2小时，由结果可知（图6所示），两个低阻 异常体和两个高阻异常体均能够较好的分辨出来，异常体中屯、电阻率分别为9 Q ? m和300 Q ? m。
[0111] W上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围 由权利要求书限定。本领域技术人员可W在本发明的实质和保护范围内，对本发明做出各 种修改或等同替换，运种修改或等同替换也应视为落在本发明的保护范围内。
【主权项】
1. 海洋大地电磁场非线性共辄梯度三维并行反演方法，其特征在于，包括如下步骤： a.在笛卡尔坐标系下沿x、y、z三个坐标轴将模型空间划分成Nx、Ny、Nz个小的长方体网 格单元，间距为 Ax(i)(i = l,…，NX)、Ay(j)(j = l,…，Ny)、Az(k)(k=l,…，NZ)，其中模型 空间包括空气层、海水层和地下半空间；海底节点（i，j，k)的长方体网格单元的长度、宽度 和高度分力U为Xseaf loor-i、yseaf loor-j、Zseaf loor-k，电阻率为Pseaf loor-i, j , k，其中电场取在各网格单元边缘的中 ' ' 羡示；磁场取在各网格单元表 面的中心J 康示； b .通过麦克斯韦方程组和边界条件的离散化得到模型各节点的场值表达式及其计算， 实现正演； C .将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差； d .通过阻抗数据偏差计算目标函数，判断目标函数值是否足够小，当目标函数值足够 小时结束迭代，否则进入下一步； e .计算目标函数的梯度，以非线性共辄梯度法作为反演方法获得查找步长和查找方 向，从而计算出模型各节点的模型修改量，其中使用阻抗数据误差和当前迭代计算所得的 电阻率计算反演方法中的预处理因子； f.使用模型修改量修正模型，并返回步骤b。2. 根据权利要求1所述的海洋大地电磁场非线性共辄梯度三维并行反演方法，其特征 在于，其中根据计算机的CPU数量自动将每个频点的数据分配给所有的CPU进行正演以及目 标函数梯度的并行计算。3. 根据权利要求1所述的海洋大地电磁场非线性共辄梯度三维并行反演方法，其特征 在于，节点的场值表达式，E_f lciOT-x如下所示：其中seafloor表示海底及地层网格节点，sea表示海水层网格节点。4. 根据权利要求1所述的海洋大地电磁场非线性共辄梯度三维并行反演方法，其特征 在于，使用电场值E，通过麦克斯韦方程可以求出各网格节点的磁场值H，便可以得到模型的 阻抗响应Z:(5)5. 根据权利要求1所述的海洋大地电磁场非线性共辄梯度三维并行反演方法，其特征 在于，反演目标函数为：其中，Cm,和z_fl_-n专指海底观测数据和正演数据，观测数据误差，第 I-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W = Wair+Wsea+ Wseafi。。!·为预调节矩阵，m=mair+msea+mseafi。。!·为网格剖分模型，air表示空气层部分，sea表示 海水层部分，seafloor表示海底及地下半空间部分， 目标函数梯度表示为：其中，K = Kair+Ksea+Kseafl。。!·为正演系数矩阵，E = Eair+Esea+Eseafl。。!·为场值并分为两种源 条件下的场值，3正演模型有限差分网格至 接收点插入两种极化源电磁场向量的线性组合； 得到：其中γ为模型电阻率和数据误差相关的变量，从而得到更新的模型： IIlseafloor-(i+l) - Hlseafloor-(i)~^Q(i)P(i) ( 18) 〇
【文档编号】G01V3/38GK106019394SQ201610268705
【公开日】2016年10月12日
【申请日】2016年4月27日
【发明人】张昆, 严加永, 蔡德超, 董浩
【申请人】中国地质科学院矿产资源研究所