CN110058315A - 一种三维各向异性射频大地电磁自适应有限元正演方法 - Google Patents

Abstract

本发明公开了一种三维各向异性射频大地电磁自适应有限元正演方法，包括：构建表示求解区域的三维地电几何模型，并获取不同区域的控制参数；将三维地电几何模型剖分为若干个互不相交的四面体单元；获取三维各向异性射频大地电磁对应的积分弱形式，并构建稀疏线性方程组，再在地表引入两种不同极化方向的入射场获得两种稀疏线性方程组中的不同右端项，并求解所述两种稀疏线性方程组得到两种极化方向上的每条棱边上的近似电场切向分量；最后基于两种极化方向上每条棱边的近似电场切向分量获取测点处的电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应中一个或多个的参数组合。该方法基于自适应有限元法实现3D RMT各向异性正演模拟。

Description

一种三维各向异性射频大地电磁自适应有限元正演方法

技术领域

本发明属于地球物理电磁勘探方法的正演技术领域，具体涉及一种三维各向异性射频大地电磁自适应有限元正演方法。

背景技术

射频大地电磁法(Radio-magnetotelluric,RMT)是在大地电磁法、可控源音频大地电磁法及探底雷达等地球物理方法的基础上发展起来的一种浅部频率域电磁勘探方法。它一般采用潜艇、无线电台发射的高频电磁波作为信号源，通过在远区采集的电磁场计算视电阻率及相位等参数，以此来研究浅地表电性结构，其探测频率约为10k-300kHz，探测深度一般在100m以内。因RMT法具有较高的探测频率，所以它不仅需要考虑电导率(与传导电流相关参数)的影响，同时还需要考虑介电常数(与位移电流相关参数)的影响。近年来，RMT法已被广泛应用于浅地表工程、水文及环境地质问题中，已成为浅地表无损检测的一种重要方法之一。

为快速、可靠地获得RMT数据的反演结果，我们需要有高效、高精度的正演求解器作为保证。而近年来，人居浅地表的电各向异性情况逐渐被人们所认识，该各向异性不仅包括电导率各向异性，同时还包括介电常数的各向异性。因此，为满足野外实测数据反演的需要，对于RMT问题的正演求解器需要同时考虑电导率和介电常数的各向异性。此外，实际的地质模型是三维的，为了避免1D、2D反演解释带来的虚假异常我们需要借助于3D正演来发展3D解释技术。因此，发展同时考虑电导率和介电常数双参数各向异性条件下的3D RMT高效、高精度正演模拟技术对提高RMT实测数据的反演解释水平具有重要的实际价值。

而现有技术中尽管只考虑电导率各向异性的3D大地电磁各向异性正演模拟技术已趋于成熟，但因该方法探测频率低而不需考虑介电常数的影响以及大地电磁与射频大地电磁法的存在本质差异，因此，并不能适用于3DRMT射频大地电磁正演模拟中。此外，由于射频大地电磁法受介电常数的影响也很严重，因此，我们很有必要提出同时考虑电导率及介电常数双参数各向异性的3DRMT各向异性正演模拟技术。

发明内容

本发明的目的是提供了一种三维各向异性射频大地电磁自适应有限元正演方法，其基于自适应有限元法实现了3D RMT各向异性正演模拟技术，可以在3D任意复杂带地形、双参数各向异性RMT问题上实现高效、高精度正演计算，且同时考虑了电导率和介电常数双参数各向异性情况，为RMT实测数据反演解释奠定了基础。

本发明提供的一种三维各向异性射频大地电磁自适应有限元正演方法，包括如下步骤：

S1：构建表示求解区域的三维地电几何模型，并获取三维地电几何模型中不同区域的控制参数；

所述控制参数为各个区域的电导率主值参数、相对介电常数主值参数、对应的主轴坐标系与观测坐标系之间的旋转角；

S2：将三维地电几何模型剖分为若干个互不相交的四面体单元；

其中，一个四面体单元仅属于三维地电几何模型中的一个区域；

S3：根据麦克斯韦方程组和Garlerkin有限元原理获取三维各向异性射频大地电磁问题对应的积分弱形式，并基于所述积分弱形式、1st order Nédélec基函数、四面体单元每条棱边的近似电场切向分量构建稀疏线性方程组，

S4：在地表引入两种不同极化方向的入射场得到两种所述稀疏线性方程组中的边界条件值(线性方程组右端项)，并求解所述两种含不同边界条件值的稀疏线性方程组获得两种极化方向上分别对应的每条棱边上的近似电场切向分量；

S5：基于两种极化方向上每条棱边的近似电场切向分量获取测点处的目标参量，所述目标参量为：电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应中一个或多个的组合。

本发明提供了一种三维各向异性射频大地电磁自适应有限元正演方法，其同时考虑了电导率及介电常数双参数各向异性对三维射频大地电磁响应的影响，通过构建积分弱形式来得到稀疏线性方程组，其中，积分弱形式是根据麦克斯韦方程组和Garlerkin有限元原理推导获得，其公式涉及到了电导率及介电常数双参数，而稀疏线性方程组是用于求解每条棱边上的近似电场切向分量，进而可以计算出两种不同极化方向的每条棱边上的近似电场切向分量，为计算出测点处电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应目标参数提供了基础。

进一步优选，步骤S3中所述积分弱形式如下：

其中，μ₀＝4π×10^-7H/m，ε₀≈8.854187817×10^-12F/m；

式中，V表示测试函数，表示阻抗率，ω表示角频率，E表示电场，表示导纳率张量，σ、ε_r分别表示各向异性电导率和相对介电常数张量，Ω表示求解区域，表示求解区域的边界，表示边界处的单位外法向量，i为虚数单位，dv、ds分别表示体积微分元和面积微分元(v表示体积、s表示面积)。

进一步优选，所述稀疏线性方程组如下所示：

Ax＝b

其中，

x＝(E₁,E₂,...,E_N)^T

式中，E₁、E₂、E_N分别表示第1条、第2条、第N条棱边的近似电场切向分量，N为所述三维地电几何模型中四面体单元棱边总数，A_ij表示矩阵A中的元素，N_i、N_j表示第i条、第j条棱边的1st order Nédélec基函数。

进一步优选，步骤S4中计算第一个极化方向上每条棱边的近似电场切向分量的过程如下：

S41：基于第一个极化方向入射场获得的边界条件值求解稀疏线性方程组得到每条棱边的近似电场切向分量；

S42：基于步骤S41计算出的每条棱边的近似电场切向分量计算每个四面体单元的后验误差，再基于后验误差计算出每个四面体单元的相对单元误差指示因子；

S43：分别判断每个四面体单元的相对误差指示因子是否大于预设最大误差，若大于，对所述四面体单元进行二次剖分(加密)，返回步骤S3，满足迭代终止条件；

所述迭代终止条件为：迭代次数达到预设最大迭代次数或求解区域的四面体单元总数达到预设最大单元数；

其中，最后得到的所有四面体单元棱边的近似电场切向分量为所需的第一个极化方向上每条棱边的近似电场切向分量。

众所周知，1st order Nédélec矢量有限元的解在很大程度上依赖于正演网格(四面体单元)的密度。随着正演网格密度的增加，解的精度将随之增加，直到收敛于真解。但采用全局加密的方式会造成很多不必要的细化，而为了尽量降低不必要的网格加密带来的时间和内存消耗，本发明基于单元后验误差对正演网格进行自适应地加密。

进一步优选，步骤S4中计算第二个极化方向上每条棱边的近似电场切向分量的过程如下：

采用第一个极化方向上最后得到的所有四面体单元以及第二个极化方向上的边界条件值求解稀疏线性方程组得到第二个极化方向上每条棱边的近似电场切向分量。

由于考虑到同一求解区域在不同极化方向上对于正演网格的差异不大，因此，当第一极化方向上基于单元后验误差对正演网格进行自适应地加密后，第二极化方向上选择与第一极化方向相同的正演网格，提高效率。

进一步优选，四面体单元的后验误差以及相对单元误差指示因子按照如下公式计算：

式中，η_e表示一个四面体单元的后验误差，F_i表示四面体单元上的第i个面，表示电流密度在四面体单元的F_i面的不连续性，表示单位外法向量，代表与四面体单元的F_i面相邻的四面体单元的导纳率张量，E_-代表与四面体单元的F_i面相邻的四面体单元上点的电场值，和E₊分别表示当前四面体单元上点的导纳率张量、电场值；

式中，β_e表示当前四面体单元的相对单元误差指示因子，表示所有四面体单元中后验误差的最大值。

进一步优选，步骤S5中所述阻抗张量、磁倾子矢量、视电阻率、相位响应的计算公式如下：

式中，Z_xx、Z_xy、Z_yx、Z_yy分别表示阻抗张量Z的四个分量，分别表示第一个极化方向得到的在x、y轴方向的电场值，分别表示第二个极化方向得到的在x、y方向上的电场值，分别表示第一个极化方向得到的在x、y方向上的磁场值，分别表示第二个极化方向得到的在x、y轴方向的磁场值；

式中，T_x、T_y分别表示在x、y轴方向的磁倾子矢量，分别表示第一个极化方向、第二个极化方向得到的在x轴方向的磁场值；

式中，表示ij方向的视电阻率，ω表示角频率，μ₀＝4π×10^-7H/m；

式中，表示ij方向的相位响应。

进一步优选，三维地电几何模型中任一区域的各向异性电导率和相对介电常数张量是根据该区域的控制参数计算而来，如下：

σ＝R_z(α_S)R_x(α_D)R_z(α_L)σ′R_z(-α_L)R_x(-α_D)R_z(-α_S)

ε_r＝R_z(α_S)R_x(α_D)R_z(α_L)ε′_rR_z(-α_L)R_x(-α_D)R_z(-α_S)

式中，σ、ε_r分别表示各向异性电导率和相对介电常数张量，α_S、α_D、α_L为主轴坐标轴xyz到观测坐标系x″y″′z′三个旋转角，R_z(α)和R_x(α)为基本初等旋转矩阵，σ′和ε′为对角矩阵，如下：

式中，σ_x,σ_y,σ_z为电导率的主值参数，为相对介电常数的主值参数。

有益效果

1、本发明针对实际射频大地电磁探测浅地表存在各向异性问题，基于自适应有限元提出了同时考虑电导率和介电常数双参数各向异性的三维射频大地电磁正演模拟方法。同时考虑了电导率和介电常数双参数各向异性对三维射频大地电磁响应的影响。提出用自适应有限元的算法能进一步提高正演的精度和效率。本发明为快速、高效的三维RMT数据的反演解释提供了核心的正演引擎，进而可以让RMT更好地服务于地质构造、矿产、水文及环境地质勘探等问题中。

2、本发明所述方法还实现了自适应加密过程，1st order Nédélec矢量有限元的解在很大程度上依赖于正演网格(四面体单元)的密度。随着正演网格密度的增加，解的精度将随之增加，直到收敛于真解。但采用全局加密的方式会造成很多不必要的细化，而为了尽量降低不必要的网格加密带来的时间和内存消耗，本发明基于单元后验误差对正演网格进行自适应地加密，既能保证精度，同时降低了不必要的消耗。

附图说明

图1为本发明原理图；

图2为本发明电导率和相对介电常数张量坐标旋转图，(a)、(b)、(c)、(d)为原图以及旋转示意图；

图3为本发明实施例提供的各向异性均匀半空间模型；

图4为本发明实例提供的各向异性半空间模型网格自适应加密过程中的三次网格切片图，(a)、(b)、(c)分别为第1次、第3次和第8次网格切片图；

图5为本发明实例提供的各向异性均匀半空间模型的视电阻率相对误差和相位绝对误差随着网格加密的收敛曲线图，(a)、(b)分别为视电阻率相对误差和相位绝对误差随网格自适应加密的收敛曲线图；

图6为本发明实例提供的各向异性梯形山脉模型；

图7为本发明实例提供的各向异性梯形山脉模型的视电阻率和相位曲线图，(a)、(b)分别为视电阻率和的示意图，(c)、(d)分别为相位和的示意图。

具体实施方式

下面将结合实施例对本发明做进一步的说明。

本发明针对现有三维各向异性射频大地电磁模拟技术的缺失，提供了一种高效、高精度地求解三维各向异性射频大地电磁自适应有限元正演方法，本发明能够为快速、可靠的三维射频大地电磁数据解释提供核心的正演引擎。本发明的一种三维各向异性射频大地电磁自适应有限元正演方法是依据下述原理进行的。

1、关于3D各向异性RMT边值问题：

建立三维地电几何模型Ω(表示求解区域)，取时间因子e^-iωt，在求解区域Ω中麦克斯韦方程组的微分形式可表示为：

其中，E、H和J分别为电场强度、磁场强度和电流密度矢量，μ₀、ε₀分别表示真空中的磁导率及介电常数，μ₀＝4π×10^-7H/m和ε₀≈8.854187817×10^-12F/m，表示阻抗率， 表示导纳率，角频率ω＝2πf，σ和ε_r分别为各向异性电导率和相对介电常数张量。

借助于如图2所示的欧拉坐标旋转，任意电导率张量σ或相对介电常数张量ε_r都可以由主轴方向(如图2所示的x、y、z分别为电导率和相对介电常数张量的主轴方向)上的主值σ_x,σ_y,σ_z或和通过主轴坐标轴xyz到观测坐标系x″y″′z′三个旋转角α_S,α_D,α_L唯一确定(如图2所示，由主轴坐标系xyz到观测坐标系x″y″′z′的旋转过程为：首先将主轴坐标系沿着z轴旋转α_S角获得x′y′z坐标系，然后将坐标系x′y′z沿着x′轴旋转α_D角获得x′y″z′坐标系，最后将坐标系x′y″z′沿着z′轴旋转α_L角获得我们的观测坐标系x″y″′z′。注意逆时针旋转时旋转角取正，顺时针旋转时旋转角取负)，任意电导率张量σ或相对介电常数张量ε_r的计算公式如下：

σ＝R_z(α_S)R_x(α_D)R_z(α_L)σ′R_z(-α_L)R_x(-α_D)R_z(-α_S) (5)

ε_r＝R_z(α_S)R_x(α_D)R_z(α_L)ε′_rR_z(-α_L)R_x(-α_D)R_z(-α_S) (6)

其中，σ′和ε′为对角矩阵，R_z(α)和R_x(α)为基本初等旋转矩阵(Pek,2002)，具体表达式如下：

因此，根据地电几何模型Ω不同区域所对应的电导率主值参数σ_x,σ_y,σ_z、相对介电常数主值参数及其各自所对应的旋转角α_S,α_D,α_L可以计算出各个区域的电导率张量σ或相对介电常数张量ε_r。

将方程(1)式带入(2)式中消去磁场强度H可以获得如下三维各向异性RMT问题在求解区域Ω中满足的电场双旋度方程，如下式(11)：

为了保证方程(11)求解的唯一性，根据(1)式，推导出如下电场在求解区域边界上应满足的第二类边界条件：

其中，Ω表示求解区域(含地下和空气)即为初始构建的三维地电模型，表示边界处的单位外法向量，H₀为各向异性均匀半空间或层状介质的解析解，具体的计算方式参考(徐震寰，2015)。方程(11)和(12)式便构成了3D各向异性RMT问题所满足的边值问题，下面将采用Garlerkin有限元法对其进行近似地求解。

2、3D RMT边值问题对应的积分弱形式

通过(11)式可以获得残差矢量R的具体表达式为：

根据加权余量方法，选取测试函数V，令其残差在计算区域Ω内的加权为0：

∫∫∫_ΩV·Rdv＝0 (14)

将残差项(13)式带入上式(14)有：

采用第一矢量格林定理(Zienkiewicz et al.,2013)并带入边界条件(12)式，我们可以获得3D各向异性RMT问题的积分弱形式：

3、Galerkin有限元分析

采用开源四面体生成器Tetgen(Si,2015)将求解区域Ω剖分成由M个互不相交的四面体单元，构成M个四面体单元边的总数为N(一个四面体单元仅处于求解区域Ω中电导率和相对介电常数的控制参数不变的一个区域，不跨区域)。并基于求解区域Ω中不同区域的模型电导率和相对介电常数的控制参数对每个四面体单元赋上电导率和相对介电常数值。然后对每个单元采用1st order Nédélec基函数可以将计算区域中任意点的电场E表示为：

其中，N_j代表第j条边的1st orderNédélec基函数，(若棱边不属于待求电场点所在的四面体单元，则N_j为零。)E_j代表第j条棱边的近似电场切向分量。基于Galerkin理论，取测试函数V为基函数，即：

V＝N_i,i＝1,...,N (18)

将(18)式带入方程(16)可获得如下大型稀疏线性方程组：

Ax＝b (19)

其中，x＝(E₁,E₂,...,E_N)^T为所要求解的近似电场切向分量，大型稀疏矩阵A中元素的具体表达式为：

对于上式中的体积分项可以借助于(Jin,2014,P268-269)的公式进行解析计算。

线性方程组(19)式的右端项b中元素的具体表达式为：

对于上式中的面积分我们采用高斯数值积分进行计算。通过调用软件包PARDISO(Schenk,2004)中的LU分解法对线性方程组(19)式进行求解可以获得高精度的电场解。

4、正演网格自适应加密

众所周知，1st orderNédélec矢量有限元的解在很大程度上依赖于正演网格的密度。随着正演网格密度的增加，解的精度将随之增加，直到收敛于真解。但采用全局加密的方式会造成很多不必要的细化，而为了尽量降低不必要的网格加密带来的时间和内存消耗，下面我们将基于单元后验误差对正演网格进行自适应地加密。

因为网格离散化带来的误差将造成电流密度在两个物性相同四面体单元的分界面处不连续，从而违背了真实的物理内部边界条件。因此，我们可以将四面体单元三角面上电流密度不连续性的定量估计作为单元后验误差：

其中项代表电流密度在四面体单元第F_i(i＝1,...,4)面的不连续性，其具体表达式为：

其中，代表与四面体单元的F_i面相邻的四面体单元的导纳率张量，E_-代表与四面体单元的F_i面相邻的四面体单元上点的电场值，和E₊分别表示当前四面体单元上点的导纳率张量、电场值。其中，是根据公式计算的，四面体单元的各向异性电导率和相对介电常数张量σ、ε_r是根据公式(5)、(6)计算出的。E₊、E_-是根据公式计算的，若j边不属于该点所在的四面体单元，此时N_j为零。相对单元误差指示因子β_e：

其中代表所有单元中最大的后验误差。预先给定一个相对单元误差允许的最大值β∈(0,1]，然后对满足β_e＞β的单元进行加密直到满足网格加密迭代结束条：迭代次数达到预先设定的允许的最大迭代次数t_max或单元总数达到预先设定的允许的最大单元数n_max。

5、有限元后处理

为了计算出3D RMT实际应用中所需要的阻抗张量、倾子矢量进而获得视电阻率和相位参数，我们需要在地表引入两种不同极化方向的入射场从而获得两种不同的边界条件值，进而获得线性方程组(19)式的两种不同右端项和然后通过代入两种不同的右端项和于线性方程组(19)式中便可获得两种不同极化方向上任意观测点的电场值E¹和E²，再根据电磁感应定律公式(1)便可求出观测点处相对应的磁场H¹和H²，进而可以获得阻抗张量响应：

倾子矢量(垂直磁场转换函数)响应：

和视电阻率和相位响应：

基于上述原理，本发明提供的一种三维各向异性射频大地电磁自适应有限元正演方法，包括如下步骤：

S1：构建表示求解区域的三维地电几何模型，并获取三维地电几何模型中不同区域的控制参数；

其中，控制参数为各个区域的电导率主值参数σ_x,σ_y,σ_z、相对介电常数主值参数对应的主轴坐标系与观测坐标系之间的旋转角α_S,α_D,α_L，因此在获知各个区域的控制参数后，根据上述公式(5)、(6)可以计算出各个区域的电导率张量σ或相对介电常数张量ε_r，同理，可以计算出任意一个四面体单元的电导率张量σ或相对介电常数张量ε_r。

S2：将三维地电几何模型剖分为若干个互不相交的四面体单元；

其中，采用开源的四面体生成器Tetgen将三维地电模型剖分成多个互不相交的四面体单元。由于一个四面体单元仅属于三维地电几何模型中一个区域，因此，在已知区域内的控制参数，根据公式(5)、(6)可以计算出任意一个四面体单元的电导率张量σ或相对介电常数张量ε_r。

S3：根据麦克斯韦方程组和Garlerkin有限元原理获取三维各向异性射频大地电磁对应的积分弱形式，并基于所述积分弱形式、1st order Nédélec基函数、四面体单元每条棱边的近似电场切向分量构建稀疏线性方程组。

从上述原理性内容描述可知，针对一个求解区域Ω可以得到式(19)形式的稀疏线性方程。

S4：在地表引入两种不同极化方向的入射场得到两种所述稀疏线性方程组中的边界条件值，并求解所述两种不同极化方向入射场分别对应的所述稀疏线性方程组得到两种极化方向上分别对应的每条棱边上的近似电场切向分量。

本发明实施例中优选先处理第一个极化方向(采用自适应加密算法)，再采取相同的方式处理第二个极化方向或者采用第一个极化方向的网格来处理第二个极化方向的数据。其中，处理第一个极化方向的数据的过程如下：

S41：基于第一个极化方向对应的边界条件值求解稀疏线性方程组得到每条棱边的近似电场切向分量；此时，稀疏线性方程组中的求解区域为本发明步骤S1中设定的三维地电几何模型，即本发明首先将求解区域进行初次剖分为若干个四面体单元，采用LU分解法对此时稀疏线性方程组进行求解，获得当前加密网格中每个棱边上的近似电场切向分量。

S42：基于步骤S41计算出的每条棱边的近似电场切向分量计算每个四面体单元的后验误差，再基于后验误差计算出每个四面体单元的相对单元误差指示因子；

S43：分别判断每个四面体单元的相对误差指示因子是否大于预设最大误差，若大于，对所述四面体单元作进一步剖分(加密)，返回步骤S3，满足迭代终止条件；

所述迭代终止条件为：迭代次数达到预设最大迭代次数或求解区域的四面体单元总数达到预设最大单元数。利用每个四面体单元的相对单元误差指示因子来鉴别四面体单元是否需要进一步剖分，若需要进一步剖分，则将该四面体单元剖分成两个四面体单元，再代入稀疏方程中求解当前更小的四面体单元上棱边的近似电场切向分量。

其中，最后得到的所有四面体单元棱边的近似电场切向分量为最终所需的第一个极化方向上每条棱边的近似电场切向分量。

本发明实施例优选在处理第二极化方向的数据时，基于第一极化方向最后的加密网格进行剖分求解区域求解稀疏方程，进而提高计算效率；其他可行的实施例中，在处理第二极化方向的数据时，还可以选择和第一极化方向上相同的方式，即初始剖分四面体单元后，再判断各个四面体单元是否需要进一步剖分，进而得到最终的加密网格以及每条棱边上的近似电场切向分量。

S5：基于两种极化方向上每条棱边的近似电场切向分量获取观测点(观测点坐标已知)处目标参量，所述目标参量为：电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应中一个或多个的组合。

为验证算法的正确性和收敛性性，建立如图3所示的实施例，该实施例为均匀各向异性半空间模型，各向异性电导率张量的主值σ_x,σ_y,σ_z分别为10⁵Ωm,10³Ωm和10⁶Ωm，对应的旋转角为α_S＝30°,α_L＝60°,α_D＝90°，相对介电常数为主轴各向异性，主值分别模型的求解区域为Ω＝[-200m,200m]³。11个观测点沿y＝[-50m,50m]以10m为间距对称地布设于地表，观测频率为100kHz。

图4为网格自适应加密过程中第1次、第3次和第8次网格及其单元相对误差示意图。图5为视电阻率相对误差和相位绝对误差随网格自适应加密的收敛曲线图，其中视电阻率相对误差和相位绝对误差分别为所有测点视电阻率ρ_xy和ρ_yx相对误差的平均值以及所有测点相位和绝对误差的平均。该模型的解析解为ρ_xy＝2191.96Ωm,ρ_yx＝6653.93Ωm,从图中可以看出，随着网格的自适应加密，视电阻率相对误差和相位绝对误差在逐渐减小，也就是说我们的有限元解在逐渐逼近真解，从而充分验证了本算法对求解各向异性射频大地电磁问题的有效性。

图6为建立的三维各向异性梯形山脉模型，其电导率为主轴各向异性，各主值分别为σ_x＝10³Ωm,σ_y＝10⁴Ωm和σ_z＝10⁵Ωm，相对介电常数表现为主轴各向异性，其主值为梯形山脉的上底为边长为45m的正方形，其中心点坐标为(0,0,-45m),下底为边长为200m的正方形，其中心点位于原点，模型求解区域为Ω＝[-1km,1km]³。19个测点沿y＝[-150m,150m]对称地分布在梯形山脉两侧，在y＝[-150m,-30m]和y＝[30m,150m]之间的测点的间距为20m，其余测点之间的间距为10m，测量频率为50k Hz。

图7为采用自适应有限元算法获得图5的各向异性梯形山脉模型的视电阻率和相位结果图。从图中可以看出由于地形的存在视电阻率从整体上表现为低阻异常，且主要受x方向的电导率影响，而主要受y方向的电导率影响；从视电阻率图中可大致分辨出山脉地形的边界。

需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。

Claims (8)

1.一种三维各向异性射频大地电磁自适应有限元正演方法，其特征在于：包括如下步骤：

S1：构建表示求解区域的三维地电几何模型，并获取三维地电几何模型中不同区域的控制参数；

所述控制参数为各个区域的电导率主值参数、相对介电常数主值参数、对应的主轴坐标系与观测坐标系之间的旋转角；

S2：将三维地电几何模型剖分为若干个互不相交的四面体单元；

其中，一个四面体单元仅属于三维地电几何模型中的一个区域；

S3：根据麦克斯韦方程组和Garlerkin有限元原理获取三维各向异性射频大地电磁对应的积分弱形式，并基于所述积分弱形式、1st order Nédélec基函数、四面体单元每条棱边的近似电场切向分量构建稀疏线性方程组，

S4：在地表引入两种不同极化方向的入射场得到所述稀疏线性方程组中的两种不同右端项，并求解所述两种含不同边界条件值的稀疏线性方程组获得两种极化方向上分别对应的每条棱边上的近似电场切向分量；

S5：基于两种极化方向上每条棱边的近似电场切向分量获取测点处的目标参量，所述目标参量为：电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应中一个或多个的组合。

2.根据权利要求1所述的方法，其特征在于：步骤S3中所述积分弱形式如下：

其中，

式中，V表示测试函数，表示阻抗率，ω表示角频率，E表示电场，表示导纳率张量，σ、ε_r分别表示各向异性电导率和相对介电常数张量，Ω表示求解区域，表示求解区域的边界，表示边界处的单位外法向量，i为虚数单位，μ₀、ε₀分别表示真空中的磁导率及介电常数，dv、ds分别表示体积微分元和面积微分元，v表示体积、s表示面积。

3.根据权利要求2所述的方法，其特征在于：所述稀疏线性方程组如下所示：

Ax＝b

其中，

x＝(E₁,E₂,...,E_N)^T

式中，E₁、E₂、E_N分别表示第1条、第2条、第N条棱边的近似电场切向分量，N为所述三维地电几何模型中四面体单元棱边总数，A_ij表示矩阵A中的元素，N_i、N_j表示第i条、第j条棱边的1st order Nédélec基函数。

4.根据权利要求1所述的方法，其特征在于：步骤S4中计算第一个极化方向上每条棱边的近似电场切向分量的过程如下：

S41：基于第一个极化方向入射场获得的边界条件值求解稀疏线性方程组得到每条棱边的近似电场切向分量；

S42：基于步骤S41计算出的每条棱边的近似电场切向分量计算每个四面体单元的后验误差，再基于后验误差计算出每个四面体单元的相对单元误差指示因子；

S43：分别判断每个四面体单元的相对误差指示因子是否大于预设最大误差，若大于，对所述四面体单元二次剖分，返回步骤S3，直至满足迭代终止条件；

所述迭代终止条件为：迭代次数达到预设最大迭代次数或求解区域的四面体单元总数达到预设最大单元数；

其中，最后得到的所有四面体单元棱边的近似电场切向分量为所需的第一个极化方向上每条棱边的近似电场切向分量。

5.根据权利要求4所述的方法，其特征在于：步骤S4中计算第二个极化方向上每条棱边的近似电场切向分量的过程如下：

采用第一个极化方向上最后得到的所有四面体单元以及第二个极化方向上的边界条件值求解稀疏线性方程组得到第二个极化方向上每条棱边的近似电场切向分量。

6.根据权利要求4所述的方法，其特征在于：四面体单元的后验误差以及相对单元误差指示因子按照如下公式计算：

式中，η_e表示一个四面体单元的后验误差，F_i表示四面体单元上的第i个面，表示电流密度在四面体单元的F_i面的不连续性，表示单位外法向量，代表与四面体单元的F_i面相邻的四面体单元的导纳率张量，E_-代表与四面体单元的F_i面相邻的四面体单元上点的电场值，和E₊分别表示当前四面体单元上点的导纳率张量、电场值；

式中，β_e表示当前四面体单元的相对单元误差指示因子，表示所有四面体单元中后验误差的最大值。

7.根据权利要求1所述的方法，其特征在于：步骤S5中所述阻抗张量、磁倾子矢量、视电阻率、相位响应的计算公式如下：

式中，Z_xx、Z_xy、Z_yx、Z_yy分别表示阻抗张量Z的四个分量，分别表示第一个极化方向得到的在x、y方向上的电场值，分别表示第二个极化方向上得到的在x、y方向上的电场值，分别表示第一个极化方向得到的在x、y轴方向的磁场值，分别表示第二个极化方向得到的在x、y轴方向的磁场值；

式中，T_x、T_y分别表示在x、y轴方向的磁倾子矢量，分别表示第一个极化方向、第二个极化方向得到的在x轴方向的磁场值；

式中，表示ij方向的视电阻率，ω表示角频率，μ₀表示真空中的磁导率；

式中，表示ij方向的相位响应。

8.根据权利要求1所述的方法，其特征在于：三维地电几何模型中任一区域的各向异性电导率和相对介电常数张量是根据所述区域的控制参数计算而来，如下：

σ＝R_z(α_S)R_x(α_D)R_z(α_L)σ′R_z(-α_L)R_x(-α_D)R_z(-α_S)

ε_r＝R_z(α_S)R_x(α_D)R_z(α_L)ε′_rR_z(-α_L)R_x(-α_D)R_z(-α_S)

式中，σ、ε_r分别表示各向异性电导率和相对介电常数张量，α_S、α_D、α_L为主轴坐标轴xyz到观测坐标系x″y″′z′三个旋转角，R_z(α)和R_x(α)为基本初等旋转矩阵，σ′和ε′为对角矩阵，如下：

式中，σ_x,σ_y,σ_z为电导率的主值参数，为相对介电常数的主值参数。