CN113221393B - 一种基于非结构有限元法的三维大地电磁各向异性反演方法 - Google Patents

Abstract

本发明公开了一种基于非结构有限元法的三维大地电磁各向异性反演方法，包括：(1)获取并筛选具有各项异性特征的大地电磁测深数据用于反演；(2)采用非结构有限元法构建电导率张量模型；(3)构建在各向异性介质条件下的大地电磁正则化反演目标函数；(4)对电导率张量模型进行正演得到模型响应对应的预测全阻抗张量，同时伴随正演计算灵敏度矩阵的转置与向量的乘积；(5)计算大地电磁正则化反演目标函数的梯度；(6)采用L‑BFGS算法计算模型参数更新量，依据模型参数更新量更新模型参数；(7)迭代执行步骤(4)～步骤(6)直到到达迭代终止条件为止，获得参数优化的电导率张量模型，实现三维大地电磁各向异性反演。

Description

一种基于非结构有限元法的三维大地电磁各向异性反演方法

技术领域

本发明属于地球物理电磁法反演技术领域，具体涉及一种基于非结构 有限元法的三维大地电磁各向异性反演方法。

背景技术

大地电磁法利用天然的低频信号能够有效探测地球深部的地电结构， 且具有成本低、工作方便、不受高阻层屏蔽以及勘探深度范围广等优势， 已经被广泛应用于深部地质调查研究、矿产等地下资源勘查。

大地电磁法反演的目的是通过地球表面测得的电磁响应来反推地下 空间的真实地电结构分布，并以此辅助地质解释工作者得到更为准确的地 质结构判断。关于大地电磁反演方法，国内外学者做了大量的研究，早期 的一维大地电磁反演方法，具有较快的速度，如一种简易的一维大地电磁 测深反演方法——博斯蒂克法反演及其应用(周虬，1985)以及大地电磁 数据反演：非线性最小二乘法(L.B.Pedersen，李然，1990)，但反演的结 构都是成层状分布的，不够接近真实的地下结构；二维大地电磁反演方法 高效稳定得到广泛应用，如应用非线性共轭梯度法的二维大地电磁反演 (Rodi,2000)，但对于复杂地形和地质结构，二维反演会有较大的误差。 真实的地下电性结构都应该是三维的，因此近几年来，三维大地电磁反演 方法的研究越来越火热。

传统大地电磁法反演技术无论是离散线性求解还是完全非线性求解 大多是基于地下空间的电性结构为均匀分布的各向同性介质，如基于有限 差分法的三维大地电磁反演方法(Egbert和Kelbert,2012)，这与实际的地 下深部地质空间是不吻合的，越来越多的地球物理勘探表明地下深部是各 向异性的，且利用各向同性的反演方法求解带有各向异性信号的大地电磁 测深数据会得到错误的结果，因此需要开发一种能够适应地下各向异性结 构的大地电磁反演方法来发展深部探测技术。

而现有的大地电磁各向异性反演技术只停留在主轴各向异性介质的 反演，对于任意的各向异性介质还无法反演，但已有的大地电磁测深各向 异性数据表明，各向异性是具有方向性，因此还需开发一种能够带电性旋 转角的各向异性三维大地电磁反演方法，来满足实际的地球物理探索工作 需求。

三维大地电磁各向异性反演问题是复杂的多参数反演问题，面临计算 效率低、消耗内存大等难题，迟迟没能应用于实际地球物理探索工作之中。

发明内容

本发明的目的在于提供一种基于非结构有限元法的三维大地电磁各 向异性反演方法，能够实现带角度的各向异性大地电磁测深数据快速三维 反演。

为实现上述发明目的，本发明提供以下技术方案：

一种基于非结构有限元法的三维大地电磁各向异性反演方法，包括以 下步骤：

(1)获取大地电磁测深数据，并依据大地电磁测深数据的各向异性 特征筛选具有各项异性特征的大地电磁测深数据用于反演；

(2)采用非结构有限元法构建电导率张量模型；

(3)根据步骤(1)筛选获得的大地电磁测深数据和电导率张量模型 构建在各向异性介质条件下的大地电磁正则化反演目标函数；

(4)对电导率张量模型进行正演得到模型响应对应的预测全阻抗张 量，同时伴随电导率张量模型的正演过程，以模型响应对模型参数的偏导 数作为灵敏度矩阵，计算灵敏度矩阵的转置与向量的乘积；

(5)依据灵敏度矩阵的转置与向量的乘积、预测全阻抗张量与大地 电磁测深数据对应的实测全阻抗张量获得的数据拟合差，来计算大地电磁 正则化反演目标函数的梯度；

(6)根据大地电磁正则化反演目标函数的梯度，并采用L-BFGS算 法计算模型参数更新量，依据模型参数更新量更新模型参数；

(7)迭代执行步骤(4)～步骤(6)直到到达迭代终止条件为止，获 得参数优化的电导率张量模型，实现三维大地电磁各向异性反演。

与现有技术相比，本发明具有的有益效果至少包括：

本发明提供的基于非结构有限元法的三维大地电磁各向异性反演方 法中，反演采用的数值模拟算法为基于三维非结构有限元法开发的全张量 电导率条件下正演算法，能够模拟任意电导率结构下的电磁场信号，其具 有较快速度下以较高精度模拟复杂地质地形结构并节省内存的优点；同时 反演方法基于有限内存拟牛顿(L-BFGS)算法进行求解，并提出了有效的正 则化约束方法及灵敏度求解方法，能够在节省大量内存的情况下快速有效 地计算复杂的各向异性多参数反演问题。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对 实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地， 下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员 来讲，在不付出创造性劳动前提下，还可以根据这些附图获得其他附图。

图1是本发明实施例提供的基于非结构有限元法的三维大地电磁各向 异性反演方法的流程图；

图2是本发明实施例提供的大地电磁测区示意图，其中，(a)为测区俯 视图(Z＝0m)，图中黑点为大地电磁测点位置，虚线框所示为地下异常体的 埋藏位置和横向形态大小图(b)为模型横切图(Y＝0m)，图中实线框所示为 异常体的埋藏位置和垂向形态大小；

图3是本发明实施例提供的各向异性识别结果图，其中，(a)、(c)、(e) 为不同频率下异常体中心点上方测点

示意图，(b)、(d)、(f)为不同频率 下该测区左上角边缘测点

示意图；

图4是本发明实施例提供的各向异性反演迭代参数图，其中，(a)为反 演过程中RMS变化图；(b)为反演过程中目标函数值变化图；(c)为反演过 程中线性搜索步长示意图；(d)为反演过程中正则化因子示意图；

图5是本发明实施例提供的大地电磁各向异性反演三维电阻率示意图， 色标为对数分布，其中，(a)、(c)、(e)分别为主轴电阻率ρ_x、ρ_y、ρ_z的横 向剖面图(Z＝4.5Km)；(b)、(d)、(f)分别为主轴电阻率ρ_x、ρ_y、ρ_z的垂向 剖面图(Y＝0Km)；

图6是本发明实施例提供的大地电磁各向异性反演三维电性旋转角度 示意图，色标为对数分布，其中，(a)、(c)、(e)分别为α、β、γ的横向剖 面图(Z＝4.5Km)；(b)、(d)、(f)分别为α、β、γ的垂向剖面图(Y＝0Km)；

图7是本发明实施例提供的大地电磁各向同性反演三维电阻率示意图， 色标为对数分布，其中，图(a)为电阻率ρ的横向剖面图(Z＝4.5Km)；图(a) 为电阻率ρ的垂向剖面图(Y＝0Km)；

图8是本发明实施例提供的各向同性反演迭代参数图，其中，(a)为反 演过程中RMS变化图；(b)为反演过程中目标函数值变化图；(c)为反演过 程中线性搜索步长示意图；(d)为反演过程中正则化因子示意图；

图9是本发明实施例提供的正演模型图；

图10是本发明实施例提供的精度验证图，其中，(a)、(b)、(c)分别为 非对角阻抗实部拟合图、非对角阻抗虚部拟合图、非对角阻抗实、虚部相 对误差示意图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及 实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实 施方式仅仅用以解释本发明，并不限定本发明的保护范围。

图1是本发明实施例提供的基于非结构有限元法的三维大地电磁各向 异性反演方法的流程图。如图1所示，实施例提供的三维大地电磁各向异 性反演方法包括以下步骤：

S1，获取大地电磁测深数据，并依据大地电磁测深数据的各向异性特 征筛选具有各项异性特征的大地电磁测深数据用于反演。

实施例中，采用极坐标下的视电阻率响应特征，作为大地电磁测深数 据的各向异性特征。设在方位角

时，极轴

对应笛卡尔坐标系下的x 轴方向，方位角方向

对应笛卡尔坐标系下的y轴方向，将大地电磁测深 数据对应的实测全阻抗张量由笛卡尔坐标系转换到极坐标系下：

依据极坐标系的实测全阻抗张量

和

计算视电阻率：

其中，笛卡尔坐标系下的大地电磁测深数据对应的实测全阻抗张量表 示为

然后将视电阻率的状态分布呈现的极性特征分布作为各项异性特征 分布，筛选具有各项异性特征的大地电磁测深数据用于反演。

在极坐标系下，原点到各方位点的矢量长度代表视电阻率的大小，对 于各向同性介质视电阻率为近似正圆形，对于各向异性介质视电阻率的分 布呈现极性特征分布，

和

极轴指示的各向异性主轴方向呈共轭关系。 通过极坐标系下的视电阻率识别，能够得到测区各测点局部范围内的大致 各向异性角度以及主轴电导率信息，进而获得各项异性特征分布。在实施 例中，可以将各向异性角度用作电导率张量模型的初始模型参数，有利于 反演的准确性。

S2，采用非结构有限元法构建电导率张量模型。

构建电导率张量模型时，在笛卡尔坐标系下将初始模型空间按照不同 区域对应不同划分细粒度的原则将电导率张量模型剖分成非结构四面体 单元，扩边区域以及空气层采用大划分细粒度的大单元，地下空间以及测 点区域采用小划分细粒度的小单元，以此来提高效率和精度。模型文件包 括节点文件、单元文件以及单元属性文件，模型参数通过属性编号与四面 体单元文件联系在一起，每个单元具有6个模型参数，即模型参数 m＝(m₁m₂ m₃ m₄ m₅ m₆)^T＝(σ_x σ_y σ_z α β γ)^T，其中，σ_x、σ_y、σ_z分别为三个主轴电导率，α、β、γ为三个主轴方向的电性旋转角。

S3，根据S1筛选获得的大地电磁测深数据和电导率张量模型构建在 各向异性介质条件下的大地电磁正则化反演目标函数。

实施例中，同时结合大地电磁测深数据的各项异性介质条件约束、相 邻单元的模型参数条件约束和主轴电导率参数条件约束构建大地电磁正 则化反演目标函数为：

其中，d表示大地电磁测深数据向量，采用全阻抗张量形式，m为电 导率张量模型参数向量，包含六个用于表征任意电各向异性张量的各向异 性参数，即m＝(m₁ m₂ m₃ m₄ m₅m₆)^T＝(σ_x σ_y σ_z α β γ)^T，λ_m和λ_c来 为正则化因子；

φ_d为数据拟合项，用于表示大地电磁测深数据的各项异性介质条件约 束，具体形式为：

其中，f(m)表示模型参数m对应的电导率张量模型通过正演得到的 模型响应，C_d表示大地电磁测深数据的协方差矩阵，T为转置算符；实施 例提供的三维大地电磁各向异性反演方法的数值模拟算法基于三维非结 构有限元法开发的全张量电导率来正演，能够模拟任意电导率结构下的电 磁场信号；

为了有效的约束多参数反演，模型约束项设置为两类模型约束项φ_m和 φ_c的线性组合；即φ_m为模型约束项，用于表示相邻单元的模型参数条件约 束，具体形式为：

其中，m₀表示模型的初始参数，一般电导率信息取均匀半空间分布， 而角度信息如果识别出来可以通过插值的方法给一个预估的空间分布，使 反演能够快速收敛避免线性求解方法陷入局部最小，C_m表示模型的粗糙的 算子，选用L₂范数来约束模型的单元与相邻单元间的光滑程度，其有限元 离散形式如下：

其中，M为四面体单元总个数，i、j表示四面体单元索引，V_i表示第 i个四面体单元体积，N(i)表示与第i个四面体单元相邻的相邻单元的总个 数，

表示处于相邻关系的第i个和第j个四面体单元间的第k个模型 参数之差，每个单元有6个模型参数，d_ij为相邻单元的质心距离；w_j表示 相邻单元的体积加权，

V_all为第j个单元的所有相邻单元的体积之 和；

φ_c为模型约束项，用于表示相邻单元的主轴电导率参数条件约束，具 体形式为：

φ_c＝(m₁ m₂ m₃)^TW_c(m₁ m₂ m₃)

其中，m₁ m₂ m₃分别表示所有单元对应的主轴电导率σ_x、σ_y、σ_z构 成的向量；

由于地球物理反演问题是不适定的，通过控制正则化因子λ_m和λ_c来使 得反演能够得到符合实际物理意义的解，对于识别出只含各向同性信号的 数据集可以使λ_c适当的增大使反演趋于各向同性，而对于一般的各向异性 信号响应则使λ_c＜＜λ_m，得到各向异性条件下的解。所有的正则化因子在反 演里都是按冷却原理给定的，一开始给出较大的初始值，当反演迭代过程 中相邻迭代间的RMS减小量小于给定的阈值时，按比例系数k更新正则化因子。

S4，对电导率张量模型进行正演得到模型响应对应的预测全阻抗张量。

实施例中，在对电导率张量模型进行正演时，电导率张量由各项异性 的主轴电导率和三个主轴方向的电性旋转角构建得到，在此基础上，采用 伽辽金有限单元法从各向异性介质下大地电磁场满足的电场双旋度方程 出发推导控制方程，选用矢量插值基函数对推导控制方程进行单元分析得 到单元矩阵形式，再利用单元间编号关系进行压缩存储组装成总体矩阵； 具体过程为：

频率域三维大地电磁电场双旋度方程：

其中，

为 旋度，i为虚数单位，ω为角频率，μ₀为真空中的磁导率；

各向异性三维大地电磁的单元矩阵形式：B^eE^e+iωμ₀A^eE^e＝C^e；

质量矩阵：

刚度矩阵：

s为单元棱边局部编号，t为第s条棱边对应棱边的局部编号，N_s和N_t分别表示四面体单元内棱边s和棱边t的电场矢量插值基函数，上标e表 示四面体单元。

对于大地电磁场来说，单元内无源项C^e＝0；任意电导率张量：

其中，R_x(α)R_y(β)R_z(γ)为三个主轴电导率方向的欧拉旋转矩阵(右手 螺旋)：

则总体矩阵形式：KE＝0,K＝B+iωμ₀A；

然后，对总体矩阵施加狄里克雷边界条件得到电场正演线性方程组， 并对电场正演线性方程组求解得到每个单元棱边电场，实施例中，可以采 用通过MUMPS直接求解器求解每个单元棱边电场，再利用矢量插值函数 插值出电场及磁场，由于采用TE、TM这两类极化方式即可计算模型的预 测全阻抗张量。具体过程为：

各向异性三维大地电磁正演线性方程组矩阵形式：KE＝b；其中，b为 右端项，包含独立极化模式下的边界棱边电场值。

磁矢量插值关系：

全阻抗张量：

本实施例中，选用全阻抗张量响应是因为该响应相对于其他电磁响应 对方向性最为敏感。

S5，同时伴随电导率张量模型的正演过程，以模型响应对模型参数的 偏导数作为灵敏度矩阵，计算灵敏度矩阵的转置与向量的乘积。

实施例中，采用伴随正演方式计算灵敏度矩阵的转置与向量的乘积。 该方式能够节省大量时间和内存并且避免了显示求取灵敏度矩阵的困难， 且这种方法并不减少反演的分辨率。

模型响应对模型参数的偏导数作为灵敏度矩阵J_k：

电导率张量模型的电场正演线性方程组求解的是棱边电场，故电场正 演线性方程组对第k个模型参数m_k求导得到

而阻抗响应通过矢量 插值函数L计算得到，则灵敏度矩阵写作如下形式：

其中，K＝B+iωμ₀A，B和A分别表示质量矩阵和刚度矩阵，单元 总体矩阵对模型参数的导数，其实先是电导率张量对各个模型参数的导数 再通过下式的数值积分方法积分得到：

系数a＝iωμ₀，i为虚数单位，ω为角频率，μ₀为真空中的磁导率，V^e表示单元体积，

表示由各项异性的主轴电导率和三个主轴方向的电性旋 转角构建得到的电导率张量，s为单元棱边局部编号，t为第s条棱边对应 棱边的局部编号，N_s和N_t分别表示四面体单元内棱边s和棱边t的电场矢 量插值基函数，s＝1,2,…,6；t＝1,2,…,6，地下均匀介质中电导率张量正 定对称以符合能量耗散为正的物理规律，故：

其中，σ_xx、σ_yy、σ_zz、σ_xy、σ_xz、σ_yz分别表示电导率张量的各个组成 元素，由于电导率张量是通过主轴电导率经欧拉旋转得到，故张量元素是 主轴电导率与各向异性角度的耦合关系；

记

则J^T＝G^T(K^-1)^TL^T；

则灵敏度矩阵的转置与向量的乘积记作J^Tv，令u＝(K^-1)^TL^Tv；

由于K为对称矩阵，K^T＝K，故解伴随正演方程组Ku＝L^Tv便可求得u， 伴随正演方程组与正演方程只有右端项不同，而大地电磁需要求解两类极 化方式(TE、TM)下的电场解，在求解时，储存两次电场正演线性方程组 所得系数矩阵K，回代入方程组Ku＝L^Tv作两次伴随正演可最终求得 J^Tv＝G^Tu。

S6，依据灵敏度矩阵的转置与向量的乘积、预测全阻抗张量与大地电 磁测深数据对应的实测全阻抗张量获得的数据拟合差，来计算大地电磁正 则化反演目标函数的梯度。

对大地电磁正则化反演目标函数对模型参数求一次导数得到梯度：

其中，数据拟合项的梯度

令v＝res即数据拟合差，模型 参数约束项的梯度为：

其中，

其中，M为四面体单元总个数，i、j表示四面体单元索引，V_i表示第 i个四面体单元体积，N(i)表示与第i个四面体单元相邻的相邻单元的总个 数，

表示处于相邻关系的第i个和第j个四面体单元间的第k个模型 参数之差，w_j表示相邻单元的体积加权，d_ij为相邻单元的质心距离；

分别求出对应梯度组装成总的大地电磁正则化反演目标函数的梯度 矩阵，设定导数控制精度ε＞0，一般为较小的实数，如10^-6；设RMS终止 值r₀≥1。

S7，根据大地电磁正则化反演目标函数的梯度，并采用L-BFGS算法 计算模型参数更新量，依据模型参数更新量更新模型参数。

大地电磁反演问题是非线性的，经泰勒展开离散化后可采用线性求解 算法求解，求解的过程实质上是最小化目标函数的过程。L-BFGS(有限内 存拟牛顿法)算法，通过迭代来获得正定近似海森矩阵，保证反演的下降 方向，并只储存有限步数的列向量对来构建海森矩阵的逆矩阵，而不储存 海森矩阵的逆矩阵，节省了大量的内存且避免了复杂计算，是一种高效的 求解最小化目标函数解的算法。其具体流程如下：

1)设最大迭代次数N，k迭代次数的索引，为给定初始海森矩阵的逆 矩阵D₀＝I，I是单位矩阵，给定步数阈值δ，通常为5～10之间的正整数；

2)计算p_k＝-D_kg_k，其中g_k为目标函数梯度，海森矩阵逆矩阵的迭代式 如下：

每当步数超过δ，就遗忘前面多出的向量对模型位移量s和梯度差y。

3)设1＞β₂＞β₁＞0，在给定的步长区间[α_min,α_max]内，线性搜索步长α_k， 满足Wolfe条件:

4)模型更新量Δm_k＝α_kp_k。

在获得模型更新量时，即模型参数更新量更新模型参数。

为加快收敛速度，减少虚假异常，对反演中的模型参数在对数域施加 约束，具体地，依据实际的地质资料和其他资料来提取物理性质特征，依 据物理性质特征设置模型参数的上界a_k和下界b_k；

依据模型参数的上界a_k和下界b_k计算约束项x_k和实际模型改变量 Δm_k'：

x_k＝ln(m_k-a_k)-ln(b_k-m_k)

依据实际模型改变量Δm_k'更新模型参数；

同时依据约束项x_k调整灵敏度矩阵为：

依据新灵敏度矩阵J_k'，迭代执行S4～S7。

S8，迭代执行S4～S7直到到达迭代终止条件为止，获得参数优化的 电导率张量模型，实现三维大地电磁各向异性反演。

实施例中，迭代终止条件包括：(1)预测全阻抗张量与实测全阻抗张 量的均方根误差RMS小于误差阈值，即RMS≤r₀；(2)大地电磁正则化反 演目标函数收敛，即||g_k||≤ε；(3)达到预设最大迭代次数N，只要满足任 何一个迭代终止条件就迭代终止，得到符合地下真实地电结构的反演解。

其中，RMS的计算公式为：

其中，e_i为预估噪声水平，n为观测数据个数。

在迭代的过程中，若没有满足迭代终止条件，则判断相邻两次迭代的 预测全阻抗张量与实测全阻抗张量的均方根误差RMS的减小量是否小于 减小量阈值χ，χ＞0，一般取2％，若小于减小量阈值则更新正则化因子， 继续迭代计算。

数值模拟方法精度验证：

正演的精度影响着反演的准确性和收敛速度，因此对本发明采用的三 维各向异性数值模拟算法进行了精度验证，选用Pek和Santos(2002)发表 的一维各向异性模型正演得到响应与一维半解析解对比，如图9所示，该 模型第一层和基底为各向同性介质层，中间两层为方位各向异性介质层， 探测的频带从0.0001Hz到100Hz取31个对数等间隔的频率分布。结果如 图10所示，其中非对角阻抗的平均相对误差小于0.7％，实、虚部最大相 对误差均小于4％，满足实际工作的需求，适用于各向异性反演。

实验例中各向异性反演：

为验证本发明的三维大地电磁各向异性反演方法的有效性，设计如下 理论模型算例，大地电磁测区示意图如图2所示，测区内布置11条测线， 每条测线有11个测点，点、线距为2.5km，各向异性异常体埋在地下100 Ω·m的均匀半空间内，其位置为图中红线所示范围，大小为10km×10km ×5km，埋深2km，异常体主轴电阻率为500/10/50Ω·m,方位电性角度为 45°。探测频率为0.01Hz到100Hz取9个对数等间隔频率。正演模型的 剖分采用非结构网格按照测点附近细网格剖分，由测区往扩边区域逐渐粗 化，对异常体局部加密，总的正演网格单元数量为272373个。将正演得 到的全阻抗张量响应加入5％的高斯噪声来模拟实测数据。先对合成数据 进行了各向异性识别，由于

与

为共轭关系，一般作

进行分析，如 图3所示，低频部分，在异常体上方测点的视电阻率极性图表现为强各向 异性特征，X方向主轴电导率相对于Y方向为高阻，方位电性角度在45° 左右，而测区边缘测点的视电阻率极性图呈近似正圆形，为各向同性特征， 对应围岩为各向同性介质特征；高频部分，全区测点的视电阻率极性图都 呈近似正圆形，对应了上覆层为各向同性介质的特征。反演的目标区域大 小为40km×40km×20km，初始模型为均匀半空间，网格的剖分除无异常 体外与正演一致，反演目标区域对应网格单元数为163140个，由于识别 出各向异性特征，设置初始正则化因子λ_c＜＜λ_m。对合成数据采用本发明的 三维大地电磁各向异性反演方法迭代12次收敛，RMS下降到0.996165， 耗时1小时37分钟。反演迭代参数图如图4所示，从图4中可以看出， RMS与目标函数在反演初期下降快，到后期变缓，RMS最终下降到 0.996165，说明了反演的稳健性和有效性，正则化因子改变缓慢，线性搜 索步长在大多数迭代过程中都为1，只做一次正演就能满足Wolfe条件， 体现了本反演方法具有较高效率和稳定性。总之，本发明的三维大地电磁 各向异性反演方法，具有计算速度快、稳定下降收敛的优点。各向异性反 演结果如图5和图6所示，本发明的三维大地电磁各向异性反演方法成功 恢复了地下异常体的结构和位置，ρ_x电阻率相对围岩为高阻特征而ρ_y方向 电阻率相对围岩为低阻特征，由于大地电磁平面波特性三维各向异性ρ_z为 不可恢复量(Yin,2003)，反演同时获得了一定的角度信息。三维大地电磁 各向异性反演得到的异常体位置对应真实模型较为准确，异常体大小形态 也较为清晰。通过理论算例证明了本发明的三维大地电磁各向异性反演方 法的有效性和稳定性，且较高的效率、较少的网格和内存使得本发明的基 于非结构有限元法三维大地电磁各向异性反演方法能够适用于普通PC。

实验例中各向同性反演：

把各向异性合成数据利用三维大地电磁各向同性反演方法进行反演， 图8所示为反演过程中反演参数的变化，反演迭代了37次，耗时近3小 时48分钟，得到的结果采用与三维大地电磁各向异性反演方法(本发明) 结果相同的色标进行作图，如图7所示，在异常体区域的电阻率为低阻， 而两边出现了虚假的高阻体异常，异常体的形态与实际相差甚远。说明， 当地下介质是明显的各向异性时，采用各向同性的反演方法，不能恢复出 地下的真实电性分布。进而证明本发明的三维大地电磁各向异性反演方法 的优越性。

以上所述的具体实施方式对本发明的技术方案和有益效果进行了详 细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制 本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等， 均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于非结构有限元法的三维大地电磁各向异性反演方法，包括以下步骤：

(1)获取大地电磁测深数据，并依据大地电磁测深数据的各向异性特征筛选具有各项异性特征的大地电磁测深数据用于反演；

(2)采用非结构有限元法构建电导率张量模型；

(3)根据步骤(1)筛选获得的大地电磁测深数据和电导率张量模型构建在各向异性介质条件下的大地电磁正则化反演目标函数；

(4)对电导率张量模型进行正演得到模型响应对应的预测全阻抗张量，同时伴随电导率张量模型的正演过程，以模型响应对模型参数的偏导数作为灵敏度矩阵，计算灵敏度矩阵的转置与向量的乘积；

(5)依据灵敏度矩阵的转置与向量的乘积、预测全阻抗张量与大地电磁测深数据对应的实测全阻抗张量获得的数据拟合差，来计算大地电磁正则化反演目标函数的梯度；

(6)根据大地电磁正则化反演目标函数的梯度，并采用L-BFGS算法计算模型参数更新量，依据模型参数更新量更新模型参数；

(7)迭代执行步骤(4)～步骤(6)直到到达迭代终止条件为止，获得参数优化的电导率张量模型，实现三维大地电磁各向异性反演；

其中，步骤(3)中，同时结合大地电磁测深数据的各项异性介质条件约束、相邻单元的模型参数条件约束和主轴电导率参数条件约束构建大地电磁正则化反演目标函数为：

φ(d,m)＝φ_d+λ_mφ_m+λ_cφ_c

其中，d表示大地电磁测深数据向量，采用全阻抗张量形式，m为电导率张量模型参数向量，包含六个用于表征任意电各向异性张量的各向异性参数，即m＝(m₁ m₂ m₃ m₄ m₅ m₆)^T＝(σ_x σ_y σ_z α β γ)^T，其中，σ_x、σ_y、σ_z分别为三个主轴电导率，α、β、γ为三个主轴方向的电性旋转角，λ_m和λ_c来为正则化因子；

φ_d为数据拟合项，用于表示大地电磁测深数据的各项异性介质条件约束，具体形式为：

其中，f(m)表示模型参数m对应的电导率张量模型通过正演得到的模型响应，C_d表示大地电磁测深数据的协方差矩阵，T为转置算符；

φ_m为模型约束项，用于表示相邻单元的模型参数条件约束，具体形式为：

其中，m₀表示模型的初始参数，C_m表示模型的粗糙的算子，选用L₂范数来约束模型的单元与相邻单元间的光滑程度；

φ_c为模型约束项，用于表示单元的主轴电导率参数条件约束，具体形式为：

φ_c＝(m₁ m₂ m₃)^TW_c(m₁ m₂ m₃)

其中，m₁ m₂ m₃分别表示所有单元对应的主轴电导率σ_x、σ_y、σ_z构成的向量。

2.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，步骤(1)中，将大地电磁测深数据对应的实测全阻抗张量由笛卡尔坐标系转换到极坐标系下后，依据极坐标系的实测全阻抗张量计算视电阻率，然后将视电阻率的状态分布呈现的极性特征分布作为各项异性特征分布，筛选具有各项异性特征的大地电磁测深数据用于反演。

3.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，步骤(2)中，构建电导率张量模型时，按照不同区域对应不同划分细粒度的原则将电导率张量模型剖分成非结构四面体单元，扩边区域以及空气层采用大划分细粒度的大单元，地下空间以及测点区域采用小划分细粒度的小单元，每个单元具有6个模型参数，即模型参数m＝(m₁ m₂ m₃ m₄ m₅ m₆)^T＝(σ_x σ_y σ_z α β γ)^T，其中，σ_x、σ_y、σ_z分别为三个主轴电导率，α、β、γ为三个主轴方向的电性旋转角。

4.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，步骤(4)中，在对电导率张量模型进行正演时，电导率张量由各项异性的主轴电导率和三个主轴方向的电性旋转角构建得到，在此基础上，采用伽辽金有限单元法从各向异性介质下大地电磁场满足的电场双旋度方程出发推导控制方程，选用矢量插值基函数对推导控制方程进行单元分析得到单元矩阵形式，再利用单元间编号关系进行压缩存储组装成总体矩阵；

然后，对总体矩阵施加狄里克雷边界条件得到电场正演线性方程组，并对电场正演线性方程组求解得到每个单元棱边电场，再利用矢量插值函数插值出电场及磁场，由于采用TE、TM这两类极化方式即可计算模型的预测全阻抗张量。

5.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，步骤(4)中，计算灵敏度矩阵的转置与向量的乘积的过程为：

模型响应对模型参数的偏导数作为灵敏度矩阵J_k：

f(m)表示模型参数m对应的电导率张量模型通过正演得到的模型响应，电导率张量模型的电场正演线性方程组求解的是棱边电场，故电场正演线性方程组对第k个模型参数m_k求导得到

而阻抗响应通过矢量插值函数L计算得到，则灵敏度矩阵写作如下形式：

其中，K＝B+iωμ₀A，B和A分别表示质量矩阵和刚度矩阵，单元总体矩阵对模型参数的导数，其实先是电导率张量对各个模型参数的导数再通过下式的数值积分方法积分得到：

系数a＝iωμ₀，i为虚数单位，ω为角频率，μ₀为真空中的磁导率，V^e表示单元体积，

表示由各项异性的主轴电导率和三个主轴方向的电性旋转角构建得到的电导率张量，s为单元棱边局部编号，t为第s条棱边对应棱边的局部编号，N_s和N_t分别表示四面体单元内棱边s和棱边t的电场矢量插值基函数，s＝1,2,…,6；t＝1,2,…,6，地下均匀介质中电导率张量正定对称以符合能量耗散为正的物理规律，故：

其中，σ_xx、σ_yy、σ_zz、σ_xy、σ_xz、σ_yz分别表示电导率张量的各个组成元素，由于电导率张量是通过主轴电导率经欧拉旋转得到，故张量元素是主轴电导率与各向异性角度的耦合关系；

记

则J^T＝G^T(K^-1)^TL^T；

则灵敏度矩阵的转置与向量的乘积记作J^Tv，令u＝(K^-1)^TL^Tv；

由于K为对称矩阵，K^T＝K，故解伴随正演方程组Ku＝L^Tv便可求得u，在求解时，储存两次电场正演线性方程组所得系数矩阵K，回代入方程组Ku＝L^Tv作两次伴随正演可最终求得J^Tv＝G^Tu。

6.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，步骤(5)中，对大地电磁正则化反演目标函数对模型参数求一次导数得到梯度：

▽φ(d,m)＝▽φ_d+λ_m▽φ_m+λ_c▽φ_c

其中，数据拟合项的梯度▽φ_d＝-2J^Tv，令v＝res即数据拟合差，模型参数约束项的梯度为：

其中

其中，M为四面体单元总个数，i、j表示四面体单元索引，V_i表示第i个四面体单元体积，N(i)表示与第i个四面体单元相邻的相邻单元的总个数，

表示处于相邻关系的第i个和第j个四面体单元间的第k个模型参数之差，w_j表示相邻单元的体积加权，d_ij为相邻单元的质心距离；

分别求出对应梯度组装成总的大地电磁正则化反演目标函数的梯度矩阵。

7.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，还包括：依据实际的地质资料来提取物理性质特征，依据物理性质特征设置模型参数的上界a_k和下界b_k；

依据第k个模型参数m_k的上界a_k和下界b_k计算约束项x_k和实际模型改变量Δm_k'：

x_k＝ln(m_k-a_k)-ln(b_k-m_k)

依据实际模型改变量Δm_k'更新模型参数；

同时依据约束项x_k调整灵敏度矩阵为：

依据新灵敏度矩阵J_k'，迭代执行步骤(4)～步骤(6)。

8.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，在迭代的过程中，若没有满足迭代终止条件，则判断相邻两次迭代的预测全阻抗张量与实测全阻抗张量的均方根误差的减小量是否小于减小量阈值，若小于减小量阈值则更新正则化因子，继续迭代计算。

9.如权利要求1所述的基于非结构有限元法的三维大地电磁各向异性反演方法，其特征在于，所述迭代终止条件包括：(1)预测全阻抗张量与实测全阻抗张量的均方根误差小于误差阈值；(2)大地电磁正则化反演目标函数收敛；(3)达到预设最大迭代次数，只要满足任何一个迭代终止条件就迭代终止。

Images