CN102798898A - 大地电磁场非线性共轭梯度三维反演方法 - Google Patents

Abstract

本发明公开了一种大地电磁场非线性共轭梯度三维反演方法，包括如下步骤：a.建立三维模型，确定模型参数；b.以大地电磁场交错网格有限差分作为正演方法通过模型电阻率获得阻抗响应；c.将正演得到的阻抗与实测阻抗对比，获得阻抗数据偏差；d.判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；e.计算目标函数的梯度，以非线性共轭梯度法作为反演方法获得查找步长和查找方向，从而计算出模型修改量，其中使用阻抗数据误差和当前迭代计算所得的电阻率作为反演方法中的预处理因子；f.修正模型，并返回步骤b。本发明方法具有高效率，高精度的特点。

Description

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

技术领域

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

背景技术

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

目前大地电磁场三维反演技术仍是行业内研究的前沿问题，其中要解决的重点问题是减少反演计算量和内存占用量、提高数据处理的效率。为了提高效率和精度，地球物理学家先后将快速松弛法、共轭梯度法、和非线性共轭梯度法引入大地电磁场的三维反演中。快速松弛法使用降低维度的方法减少计算量，但计算精度有限；共轭梯度法使用矩阵乘积代替了雅克比矩阵的计算，提高了计算效率，但是收敛较慢；非线性共轭梯度法不但取消了雅克比矩阵的计算，而且采用预处理和线性查找的方法提高计算效率，加快反演的收敛速度。因此非线性共轭梯度法是目前解决大地电磁场三维反演问题较好的方法。

所以现有的大地电磁场反演方法中以大地电磁场交错网格有限差分作为正演方法（Smith，1996），以非线性共轭梯度作为反演方法（Newman & Alumbaugh，2000）。但是非线性共轭梯度法中预处理因子的选择十分重要，该参数的选择关系反演效率和精度。现有的非线性共轭梯度法中由于预处理因子的选择要么十分复杂，计算效率较低；要么较为简单，并且与初始模型电阻率相关，降低了算法的有效性。

鉴于上述现有的大地电磁场三维反演方法中存在的问题和缺陷，本发明人依靠多年的工作经验和丰富的专业知识积极加以研究和创新，最终发明了一种大地电磁场非线性共轭梯度三维反演方法，具有高效率，高精度的特点。

发明内容

为了解决现有技术中存在的上述问题，本发明提供了一种大地电磁场非线性共轭梯度三维反演方法，具有高效率，高精度的特点。

为了解决上述技术问题，本发明采用了如下技术方案：

大地电磁场非线性共轭梯度三维反演方法，包括如下步骤：

a.建立三维模型，确定模型参数；

b.以大地电磁场交错网格有限差分作为正演方法通过模型电阻率获得阻抗响应；

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

d.判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；

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

f.修正模型，并返回步骤b。

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

进一步，步骤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_(i)、Δy_(j)、Δz_(k)，其电阻率为ρ(i,j,k)，其中电场取在网格单元边缘的中点，磁场取在网格单元表面的中心。

进一步，步骤b中，通过下述方程组（3）得到阻抗，

${\▿}_{\mathrm{sfd}}\·{\▿}_{\mathrm{sfd}}\×W=0\mathrm{for\; all\; W}$

${\▿}_{\mathrm{sfd}}\×{\▿}_{\mathrm{sfd}}\mathrm{\φ}=0\mathrm{for\; all\φ}$

${\▿}_{\mathrm{sfd}}\×{\▿}_{\mathrm{sfd}}\×W={\▿}_{\mathrm{sfd}}[{\▿}_{\mathrm{sfd}}\·W]-{\▿}_{\mathrm{sfd}}\·\left[{\▿}_{\mathrm{sfd}}W\right]\mathrm{for\; all\; W}---\left(3\right)$

其中W为任意矢量场，φ为任意标量场，

为差值算子，

为散度算子，

为旋度算子，由式3得到三维模型各节点的电场与磁场分量，进而得到阻抗。

进一步，步骤e中，目标函数为：

其中

和Z_n分别表示由野外观测数据得到的阻抗和正演得到的阻抗，第1-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W为对角线元素值为1的斜对角矩阵，

为目标函数，ε为实测数据误差，λ为正则化因子（可根据实测数据情况调节），m为模型电阻率矢量；

目标函数梯度表示为：

其中，

Re表示复数取实部；

其中，K^-1为刚度矩阵的逆，E为场值，E₁和E₂分别为不同场源模式下计算得到的电场场值，g的各分量为正演模型有限差分网格至接收机在位置j的地点插入两种极化源电磁场向量的线性组合，得到：令 $g_T^1=-2\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left({\mathrm{\ΔZ}}_n\right)g_n^T{}_{*1},$ $g_T^2=-2\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left({\mathrm{\ΔZ}}_n\right)g_n^T{}_{*2},$ ¹v^T=¹g^TK^-1、²v^T=²g^TK^-1，经过两次正演计算：K¹v=¹g和K²v=²g便可得到v的解，从而得到目标函数中数据偏差项对模型参数的导数；

非线性共轭梯度的查找方向和查找步长由单减或沿着查找方向计算的线性查找确定，采用如下方程组（10）：

m₀=given

m₀为初始模型，m_k为第k次迭代的模型电阻率参数，p为查找方向，α为查找步长；

每一次迭代中，查找方向随预处理因子和目标函数梯度变化采用如下方程组（11）：

引用Polak-Ribiere方法：

p_k为第k次迭代的查找方向，C为预处理因子，

是一个最速下降方向，最小化

在m_k上的方向导数，这里C=(γI+λL^TL)^-1，其中γ是与实测数据误差和当前迭代模型电阻率相关的矢量，I为主对角元素值为1的主对角阵，L为与模型网格节点插值元素有关的预调节矩阵，通常，对h求解方程组将预调节矩阵作用于

便可有由方程组（11）和式（12）求得查找方向；

定义单减方程最小化为φ_k，并且高斯-牛顿方程近似为φ_k：

第l次线性查找产生模型修改量：

m_k,l=m_k+α_k,lp_k,l=0,1,2,…（14）

其中，α_k,0=0，这样最小值就表示成：

其中，H为目标函数的海森矩阵；

通过迭代求解，最终求解最佳拟合模型。

进一步，目标函数值包括由野外采集数据得到的电性阻抗与模型阻抗响应的偏差和当前迭代电阻率与初始模型电阻率的差，目标函数值的标准偏差达到数据误差水平即为目标函数值足够小。

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

本发明的大地电磁场非线性共轭梯度三维反演方法在已有的非线性共轭梯度算法、三维大地电磁场模拟方法、以及三维反演预处理算法基础上，改进了预处理方法，使用阻抗数据误差和当前迭代计算所得的模型电阻率作为反演方法中的预处理因子，减少了每次迭代中8次拟正演计算，很大程度上提高了计算效率，并且保证了计算精度。

本发明方法能够将野外采集到的数据通过反演计算转换成地下三维电性结构信息，反映更为真实的地下以及深部电性结构信息，根据测量数据能够反映出地壳浅表、地壳浅层、地壳深层、壳幔边界、上地幔的电阻率信息。能够为探索成矿区带不同深度尺度上地球物理特征与成矿之间的关系、形成矿产资源立体探测的技术解决方案和深部资源勘查提供依据技术支持，指明找矿方向，查明成矿系统的结构，重塑深部成矿动力学过程，为建立陆内成矿理论提供更可靠的信息。

附图说明

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

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

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

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

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

具体实施方式

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

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

a.建立三维模型，确定模型参数；

b.以大地电磁场交错网格有限差分作为正演方法通过模型电阻率获得阻抗响应；

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

d.判断目标函数值（包括由野外采集数据得到的电性阻抗与模型阻抗响应的偏差和当前迭代电阻率与初始模型电阻率的差）是否足够小（一般情况下目标函数值的标准偏差达到数据误差水平即为满足要求），当目标函数值足够小时结束迭代，否则进入下一步；

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

f.修正模型，并返回步骤b。

作为优选，根据Newman在2000年发表的文章，设随时间变化的谐波函数为e^iωt，

ω是角频率，电场强度E满足矢量方程：

$\▿\×{\mathrm{\μ}}_0/\mathrm{\μ}\▿\×E+\mathrm{i\ω}{\mathrm{\μ}}_0\mathrm{\σE}=0---\left(1\right)$

$H=\▿\×E/(-\mathrm{i\ω\μ})---\left(2\right)$

其中，H为磁场强度，σ为电导率，μ为磁导率，μ₀为真空磁导率，为旋度算子。

三维模型建立如下：图2为本发明的三维模型的交错采样网格示意图，如图2所示，在笛卡尔坐标系下沿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)。

图3为图2中的单个长方体网格单元示意图，图3所示的是编号为(i，j，k)的长方体网格单元，其长度、宽度和高度分别为Δx_(i)、Δy_(j)、Δz_(k)，电阻率为ρ(i,j,k)。其中电场取在各网格单元边缘的中点，磁场取在各网格单元表面的中心。按照此设置将Maxwell方程组转换成为：

${\▿}_{\mathrm{sfd}}\·{\▿}_{\mathrm{sfd}}\×W=0\mathrm{for\; all\; W}$

${\▿}_{\mathrm{sfd}}\×{\▿}_{\mathrm{sfd}}\mathrm{\φ}=0\mathrm{for\; all\φ}$

${\▿}_{\mathrm{sfd}}\×{\▿}_{\mathrm{sfd}}\×W={\▿}_{\mathrm{sfd}}[{\▿}_{\mathrm{sfd}}\·W]-{\▿}_{\mathrm{sfd}}\·\left[{\▿}_{\mathrm{sfd}}W\right]\mathrm{for\; all\; W}(\mathrm{Smith},1996)---\left(3\right)$

其中W为任意矢量场，φ为任意标量场，为差值算子，

为散度算子，为旋度算子，由式3得到三维模型各节点的电场与磁场分量，进而得到阻抗。

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

目标函数表示为（Newman & Alumbaugh，2000）：

其中

和Z_n分别表示由野外观测数据得到的阻抗和正演得到的阻抗，第1-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W为对角线元素值为1的斜对角矩阵，

为目标函数，ε为实测数据误差，γ为正则化因子（可根据实测数据情况调节），m为模型电阻率矢量，阻抗表达式如下：

$Z=\left(\begin{array}{cc}{Z}_{\mathrm{xx}}& Z_{\mathrm{xy}}\\ Z_{\mathrm{yx}}& Z_{\text{yy}}\end{array}\right)---\left(5\right)$

目标函数梯度表示为：

其中，

Re表示复数取实部。

根据Rodi在1976年给出的表达式，三维阻抗的灵敏度张量表示为如下方程组：

${\∂Z}_{\mathrm{xx}}/{\∂m}_k=-g_{\mathrm{xx}}^T{}_1{K}^{-1}(\∂K/{\∂m}_k{E}_1)-g_{\mathrm{xx}}^T{}_2{K}^{-1}(\∂K/{\∂m}_k{E}_2)$

${\∂Z}_{\mathrm{xy}}/{\∂m}_k=-g_{\mathrm{xy}}^T{}_1{K}^{-1}(\∂K/{\∂m}_k{E}_1)-g_{\mathrm{xy}}^T{}_2{K}^{-1}(\∂K/{\∂m}_k{E}_2)$

${\∂Z}_{\mathrm{yx}}/{\∂m}_k=-g_{\mathrm{yx}}^T{}_1{K}^{-1}(\∂K/{\∂m}_k{E}_1)-g_{\mathrm{yx}}^T{}_2{K}^{-1}(\∂K/{\∂m}_k{E}_2)$

${\∂Z}_{\mathrm{yy}}/{\∂m}_k=-g_{\mathrm{yy}}^T{}_1{K}^{-1}(\∂K/{\∂m}_k{E}_1)-g_{\mathrm{yy}}^T{}_2{K}^{-1}(\∂K/{\∂m}_k{E}_2)---\left(7\right)$

其中，K^-1为刚度矩阵的逆，E为场值，E₁和E₂分别为不同场源模式下计算得到的电场场值，g的各分量为正演模型有限差分网格至接收机在位置j的地点插入两种极化源电磁场向量的线性组合。得到：

令 $g_T^1=-2\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left({\mathrm{\ΔZ}}_n\right)g_n^T{}_{*1},$ $g_T^2=-2\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left({\mathrm{\ΔZ}}_n\right)g_n^T{}_{*2},$ ¹v^T=¹g^TK^-1、²v^T=²g^TK^-1。

经过两次正演计算：K¹v=¹g和K²v=²g便可得到v的解，从而得到目标函数中数据偏差项对模型参数的导数。

非线性共轭梯度的查找方向和查找步长由单减或沿着查找方向计算的线性查找确定，采用方程组如下：

m₀=given

m₀为初始模型，m_k为第k次迭代的模型电阻率参数，p为查找方向，α为查找步长

每一次迭代中，查找方向随预处理因子和目标函数梯度变化，采用方程组如下：

引用Polak-Ribiere方法：

p_k为第k次迭代的查找方向，C为预处理因子，

是一个最速下降方向，最小化

在m_k上的方向导数。

定义单减方程最小化为φ_k，并且高斯-牛顿方程近似为φ_k：

第l次线性查找产生模型修改量：

m_k,l=m_k+α_k,lp_k,l=0,1,2,…（14）

其中，α_k,0=0，这样最小值就表示成：

其中，H为目标函数的海森矩阵。

通过迭代求解，最终求解最佳拟合模型。

由式（15）可见，α的每一次计算都用到海森矩阵，但是海森海森为目标函数的二次倒数，其计算量十分巨大，为了避免大量计算，提高计算效率，需要使用预处理方法代替海森矩阵的计算。

在非线性共轭梯度法的共轭梯度迭代中使用了预处理因子C。该因子可能会对非线性共轭梯度方法的效率有一些影响。对其选取主要考虑两个方面：使用预处理矩阵的计算量大小以及使梯度向量成为有效的搜索方向的能力。

这里C=(γI+λL^TL)^-1，其中γ是与实测数据误差和当前迭代模型电阻率相关的矢量，I为主对角元素值为1的主对角阵，L为与模型网格节点插值元素有关的预调节矩阵。通常，对h求解方程组

将预调节矩阵作用于

其原理是存在一个有效算子，在一定意义下起到近似海森矩阵

逆的作用。求解上述方程组的计算量远远小于一次正演计算，但是它取代了8次拟正演计算，因此对于非线性共轭梯度反演减少了大量计算量。

并行计算

本发明方法通过使用阻抗数据误差和当前迭代计算所得的电阻率作为反演方法中的预处理因子在每一次迭代中已经减少了8次拟正演计算，但是三维正演的计算效率仍是缩短消耗时间的瓶颈。阻抗、查找步长和查找方向三个参数需要进行正演或者拟正演计算，而就目前的正演模拟现状而言，三维交错有限差分方法是相对高效的方法之一，但其串行计算速度受数据频点数量限制。另一个正演模型的难点是存储大规模网格节点数据，对于三维反演问题，其模型网格数量可能达到百万量级，而进行正演模拟时，刚度矩阵便是一个元素数量为万亿量级的矩阵，这对内存的消耗十分巨大。

为此，首先对所有稀疏矩阵进行一维定位存储，减少内存消耗，并采用并行计算反演方法。该方法会根据计算机线程数量自动分配每个线程的工作量，将各频点数据分散给不同的线程，进行并行运算。显著减少了解算时间。

交错网格有限差分正演模拟方法在计算刚度矩阵、场值和阻抗时需要所有网格节点参数参与计算，因此，如果将模型分散为小块体分布于各计算线程，会破坏正演模拟的边界条件，从而影响计算精度。但是测点不同频点的数据彼此毫不相关，所以本发明采用将各频点数据分散于不同线程进行并行运算不会影响计算精度。

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

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

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

模型一

图4为本发明方法的应用例1的模型示意图，如图4所示，设置为一个异常体模型，低阻异常体电阻率为10Ω·m、宽为1000m、高为600m，顶部埋深400m，地下背景介质电阻率100Ω·m，模型网格设置如下：

数据频率为0.01s、0.1s、1s、10s、100s，9个测点位于网格中心，迭代94次后满足终止条件，RMS为0.0328，主频为2.93GHz的8线程计算机运行时间为66s，（a）为y、z方向的电阻率断面，中心600m*1000m范围内电阻率为10Ω·m，其上方网格电阻率为102Ω·m，下方网格电阻率为77Ω·m，左右网格电阻率为89Ω·m；（b）为x、z方向的电阻率断面，其上方网格电阻率为102Ω·m，下方网格电阻率为77Ω·m，左右网格电阻率为90Ω·m。可以看出本发明方法的反演收敛速度较高，进行94次迭代后RMS达到最小值0.0328，计算效率和计算精度均有较好表现。

模型二

图5为本发明方法的应用例2的模型示意图，如图5所示，设置为四个异常体模型，100Ω·m的半空间背景中包含2个低阻异常体和2个高阻异常体，异常体电阻率分别为10Ω·m和1000Ω·m，异常体顶深2km，长2km，宽2km、深3km、间距2km，数据频率范围0.01-100s，共17个频点，网格数为30*30*30，64个测点位于网格中心。在本应用例中迭代86次后满足终止条件，RMS为0.045，主频为2.93GHz的8线程计算机运行时间为73274s，由结果可知，两个低阻异常体和两个高祖异常体均能够较好的分辨出来，异常体中心电阻率分别为20Ω·m和174Ω·m。

使用本发明方法以日本KAYABE地区实测数据为例。

数据频率范围64-2Hz，共6个频点，网格数为23*23*30，150个测点位于网格中心。使用本发明方法的反演结果显示，全张量阻抗反演最终RMS为2.09，测区地下500m深度范围内主要体现为地阻异常。测区地下浅表电阻率分布较为凌乱，南部出现小面积高阻异常，而50m-200m深处，低阻异常由东向西收拢，主要表现为西南和东南两个低阻异常体，至300m以下，西南部以及中部的低阻异常体逐渐变小，并且与东南部低阻异常体分离。反演在第70次迭代时满足拟合差与模型偏差总和最小而结束反演，反演收敛速度较快，基本上在7次反演之后便能够达到较小的RMS值，计算效率和计算精度较高。本发明方法能够准确的反映出地下介质的电性结构，对异常体的边界、埋深以及延伸范围都有较好的反映。

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

Claims (6)

1.大地电磁场非线性共轭梯度三维反演方法，其特征在于，包括如下步骤：

a.建立三维模型，确定模型参数；

b.以大地电磁场交错网格有限差分作为正演方法通过模型电阻率获得阻抗响应；

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

d.判断目标函数值是否足够小，当目标函数值足够小时结束迭代，否则进入下一步；

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

f.修正模型，并返回步骤b。

2.根据权利要求1所述的大地电磁场非线性共轭梯度三维反演方法，其特征在于，其中根据计算机的CPU数量自动将每个频点的数据分配给所需的CPU进行正演以及目标函数梯度的并行计算。

3.根据权利要求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_(i)、Δy_(j)、Δz_(k)，其电阻率为ρ(i,j,k)，其中电场取在网格单元边缘的中点，磁场取在网格单元表面的中心。

4.根据权利要求1所述的大地电磁场非线性共轭梯度三维反演方法，其特征在于，步骤b中，通过下述方程组（3）得到阻抗，

${\▿}_{\mathrm{sfd}}\·{\▿}_{\mathrm{sfd}}\×W=0\mathrm{for\; all\; W}$

${\▿}_{\mathrm{sfd}}\×{\▿}_{\mathrm{sfd}}\mathrm{\φ}=0\mathrm{for\; all\φ}$

${\▿}_{\mathrm{sfd}}\×{\▿}_{\mathrm{sfd}}\×W={\▿}_{\mathrm{sfd}}[{\▿}_{\mathrm{sfd}}\·W]-{\▿}_{\mathrm{sfd}}\·\left[{\▿}_{\mathrm{sfd}}W\right]\mathrm{for\; all\; W}---\left(3\right)$

其中W为任意矢量场，φ为任意标量场，

为差值算子，

为散度算子，

为旋度算子，由式3得到三维模型各节点的电场与磁场分量，进而得到阻抗。

5.根据权利要求1所述的大地电磁场非线性共轭梯度三维反演方法，其特征在于，步骤e中，目标函数为：

其中

和Z_n分别表示由野外观测数据得到的阻抗和正演得到的阻抗，第1-N个数据为阻抗各分量的实部，第N+1-2N个数据为阻抗各分量的虚部，W为对角线元素值为1的斜对角矩阵，

为目标函数，ε为实测数据误差，λ为正则化因子，m为模型电阻率矢量；

目标函数梯度表示为：

其中，

Re表示复数取实部；

其中，K^-1为刚度矩阵的逆，E为场值，E₁和E₂分别为不同场源模式下计算得到的电场场值，g的各分量为正演模型有限差分网格至接收机在位置j的地点插入两种极化源电磁场向量的线性组合,从而得到目标函数中数据偏差项对模型参数的导数；

非线性共轭梯度的查找方向和查找步长由单减或沿着查找方向计算的线性查找确定，如下方程组（10）：

m₀=given

m₀为初始模型，m_k为第k次迭代的模型电阻率参数，p为查找方向，α为查找步长；

每一次迭代中，查找方向随预处理因子和目标函数梯度变化如下述方程组（11）：

引用Polak-Ribiere方法：

p_k为第k次迭代的查找方向，C为预处理因子，

是一个最速下降方向，最小化

在m_k上的方向导数，这里C=(γI+λL^TL)^-1，其中γ是与实测数据误差和当前迭代模型电阻率相关的矢量，I为主对角元素值为1的主对角阵，L为与模型网格节点插值元素有关的预调节矩阵，通常，对h求解方程组

将预调节矩阵作用于便可有由方程组（11）和式（12）求得查找方向；

定义单减方程最小化为φ_k，并且高斯-牛顿方程近似为φ_k：

第l次线性查找产生模型修改量：

m_k,l=m_k+α_k,lp_k,l=0,1,2,…（14）

其中，α_k,0=0，这样最小值就表示成：

其中，H为目标函数的海森矩阵；

通过迭代求解，最终求解最佳拟合模型。

6.根据权利要求1所述的大地电磁场非线性共轭梯度三维反演方法，其特征在于，目标函数值包括由野外采集数据得到的电性阻抗与模型阻抗响应的偏差和当前迭代电阻率与初始模型电阻率的差，目标函数值的标准偏差达到数据误差水平即为目标函数值足够小。

Images