CN110058315B - 一种三维各向异性射频大地电磁自适应有限元正演方法 - Google Patents
一种三维各向异性射频大地电磁自适应有限元正演方法 Download PDFInfo
- Publication number
- CN110058315B CN110058315B CN201910457480.7A CN201910457480A CN110058315B CN 110058315 B CN110058315 B CN 110058315B CN 201910457480 A CN201910457480 A CN 201910457480A CN 110058315 B CN110058315 B CN 110058315B
- Authority
- CN
- China
- Prior art keywords
- electric field
- representing
- tetrahedral
- dimensional
- polarization direction
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 55
- 230000003044 adaptive effect Effects 0.000 title claims description 17
- 230000005684 electric field Effects 0.000 claims abstract description 62
- 230000010287 polarization Effects 0.000 claims abstract description 57
- 239000013598 vector Substances 0.000 claims abstract description 22
- 230000004044 response Effects 0.000 claims abstract description 14
- 230000006870 function Effects 0.000 claims description 14
- 230000008569 process Effects 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 239000000126 substance Substances 0.000 claims description 6
- 150000001875 compounds Chemical class 0.000 claims description 4
- 238000012360 testing method Methods 0.000 claims description 4
- 230000035699 permeability Effects 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000001514 detection method Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 238000000354 decomposition reaction Methods 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 241000282414 Homo sapiens Species 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000005672 electromagnetic field Effects 0.000 description 1
- 230000005674 electromagnetic induction Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000009659 non-destructive testing Methods 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 238000010587 phase diagram Methods 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/08—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
- G01V3/083—Controlled source electromagnetic [CSEM] surveying
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/08—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
- G01V3/083—Controlled source electromagnetic [CSEM] surveying
- G01V2003/086—Processing
Landscapes
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Measuring Magnetic Variables (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种三维各向异性射频大地电磁自适应有限元正演方法,包括:构建表示求解区域的三维地电几何模型,并获取不同区域的控制参数;将三维地电几何模型剖分为若干个互不相交的四面体单元;获取三维各向异性射频大地电磁对应的积分弱形式,并构建稀疏线性方程组,再在地表引入两种不同极化方向的入射场获得两种稀疏线性方程组中的不同右端项,并求解所述两种稀疏线性方程组得到两种极化方向上的每条棱边上的近似电场切向分量;最后基于两种极化方向上每条棱边的近似电场切向分量获取测点处的电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应中一个或多个的参数组合。该方法基于自适应有限元法实现3D RMT各向异性正演模拟。
Description
技术领域
本发明属于地球物理电磁勘探方法的正演技术领域,具体涉及一种三维各向异性射频大地电磁自适应有限元正演方法。
背景技术
射频大地电磁法(Radio-magnetotelluric,RMT)是在大地电磁法、可控源音频大地电磁法及探底雷达等地球物理方法的基础上发展起来的一种浅部频率域电磁勘探方法。它一般采用潜艇、无线电台发射的高频电磁波作为信号源,通过在远区采集的电磁场计算视电阻率及相位等参数,以此来研究浅地表电性结构,其探测频率约为10k-300kHz,探测深度一般在100m以内。因RMT法具有较高的探测频率,所以它不仅需要考虑电导率(与传导电流相关参数)的影响,同时还需要考虑介电常数(与位移电流相关参数)的影响。近年来,RMT法已被广泛应用于浅地表工程、水文及环境地质问题中,已成为浅地表无损检测的一种重要方法之一。
为快速、可靠地获得RMT数据的反演结果,我们需要有高效、高精度的正演求解器作为保证。而近年来,人居浅地表的电各向异性情况逐渐被人们所认识,该各向异性不仅包括电导率各向异性,同时还包括介电常数的各向异性。因此,为满足野外实测数据反演的需要,对于RMT问题的正演求解器需要同时考虑电导率和介电常数的各向异性。此外,实际的地质模型是三维的,为了避免1D、2D反演解释带来的虚假异常我们需要借助于3D正演来发展3D解释技术。因此,发展同时考虑电导率和介电常数双参数各向异性条件下的3DRMT高效、高精度正演模拟技术对提高RMT实测数据的反演解释水平具有重要的实际价值。
而现有技术中尽管只考虑电导率各向异性的3D大地电磁各向异性正演模拟技术已趋于成熟,但因该方法探测频率低而不需考虑介电常数的影响以及大地电磁与射频大地电磁法的存在本质差异,因此,并不能适用于3DRMT射频大地电磁正演模拟中。此外,由于射频大地电磁法受介电常数的影响也很严重,因此,我们很有必要提出同时考虑电导率及介电常数双参数各向异性的3DRMT各向异性正演模拟技术。
发明内容
本发明的目的是提供了一种三维各向异性射频大地电磁自适应有限元正演方法,其基于自适应有限元法实现了3D RMT各向异性正演模拟技术,可以在3D任意复杂带地形、双参数各向异性RMT问题上实现高效、高精度正演计算,且同时考虑了电导率和介电常数双参数各向异性情况,为RMT实测数据反演解释奠定了基础。
本发明提供的一种三维各向异性射频大地电磁自适应有限元正演方法,包括如下步骤:
S1:构建表示求解区域的三维地电几何模型,并获取三维地电几何模型中不同区域的控制参数;
所述控制参数为各个区域的电导率主值参数、相对介电常数主值参数、对应的主轴坐标系与观测坐标系之间的旋转角;
S2:将三维地电几何模型剖分为若干个互不相交的四面体单元;
其中,一个四面体单元仅属于三维地电几何模型中的一个区域;
S3:根据麦克斯韦方程组和Garlerkin有限元原理获取三维各向异性射频大地电磁问题对应的积分弱形式,并基于所述积分弱形式、1st order Nédélec基函数、四面体单元每条棱边的近似电场切向分量构建稀疏线性方程组,
S4:在地表引入两种不同极化方向的入射场得到两种所述稀疏线性方程组中的边界条件值(线性方程组右端项),并求解所述两种含不同边界条件值的稀疏线性方程组获得两种极化方向上分别对应的每条棱边的近似电场切向分量;
S5:基于两种极化方向上分别对应的每条棱边的近似电场切向分量获取测点处的目标参量,所述目标参量为:电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应中一个或多个的组合。
本发明提供了一种三维各向异性射频大地电磁自适应有限元正演方法,其同时考虑了电导率及介电常数双参数各向异性对三维射频大地电磁响应的影响,通过构建积分弱形式来得到稀疏线性方程组,其中,积分弱形式是根据麦克斯韦方程组和Garlerkin有限元原理推导获得,其公式涉及到了电导率及介电常数双参数,而稀疏线性方程组是用于求解每条棱边上的近似电场切向分量,进而可以计算出两种不同极化方向的每条棱边上的近似电场切向分量,为计算出测点处电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应目标参数提供了基础。
进一步优选,步骤S3中所述积分弱形式如下:
式中,V表示测试函数,表示阻抗率,ω表示角频率,E表示电场,表示导纳率张量,σ、εr分别表示各向异性电导率和相对介电常数张量,Ω表示求解区域,表示求解区域的边界,表示边界处的单位外法向量,i为虚数单位,dv、ds分别表示体积微分元和面积微分元(v表示体积、s表示面积)。
进一步优选,所述稀疏线性方程组如下所示:
Ax=b
其中,
x=(E1,E2,...,EN)T
式中,E1、E2、EN分别表示第1条、第2条、第N条棱边的近似电场切向分量,N为所述三维地电几何模型中四面体单元棱边总数,Aij表示矩阵A中的元素,Ni、Nj表示第i条、第j条棱边的1st order Nédélec基函数。
进一步优选,步骤S4中计算第一个极化方向上每条棱边的近似电场切向分量的过程如下:
S41:基于第一个极化方向入射场获得的边界条件值求解稀疏线性方程组得到每条棱边的近似电场切向分量;
S42:基于步骤S41计算出的每条棱边的近似电场切向分量计算每个四面体单元的后验误差,再基于后验误差计算出每个四面体单元的相对误差指示因子;
S43:分别判断每个四面体单元的相对误差指示因子是否大于预设最大误差,若大于,对所述四面体单元进行二次剖分(加密),返回步骤S3,满足迭代终止条件;
所述迭代终止条件为:迭代次数达到预设最大迭代次数或求解区域的四面体单元总数达到预设最大单元数;
其中,最后得到的所有四面体单元棱边的近似电场切向分量为所需的第一个极化方向上每条棱边的近似电场切向分量。
众所周知,1st order Nédélec矢量有限元的解在很大程度上依赖于正演网格(四面体单元)的密度。随着正演网格密度的增加,解的精度将随之增加,直到收敛于真解。但采用全局加密的方式会造成很多不必要的细化,而为了尽量降低不必要的网格加密带来的时间和内存消耗,本发明基于单元后验误差对正演网格进行自适应地加密。
进一步优选,步骤S4中计算第二个极化方向上每条棱边的近似电场切向分量的过程如下:
采用第一个极化方向上最后得到的所有四面体单元以及第二个极化方向上的边界条件值求解稀疏线性方程组得到第二个极化方向上每条棱边的近似电场切向分量。
由于考虑到同一求解区域在不同极化方向上对于正演网格的差异不大,因此,当第一极化方向上基于单元后验误差对正演网格进行自适应地加密后,第二极化方向上选择与第一极化方向相同的正演网格,提高效率。
进一步优选,四面体单元的后验误差以及相对误差指示因子按照如下公式计算:
式中,ηe表示一个四面体单元的后验误差,Fi表示四面体单元上的第i个面,表示电流密度在四面体单元的Fi面的不连续性,表示单位外法向量,代表与四面体单元的Fi面相邻的四面体单元的导纳率张量,E-代表与四面体单元的Fi面相邻的四面体单元上点的电场值,和E+分别表示当前四面体单元上点的导纳率张量、电场值;
进一步优选,步骤S5中所述阻抗张量、磁倾子矢量、视电阻率、相位响应的计算公式如下:
式中,Zxx、Zxy、Zyx、Zyy分别表示阻抗张量Z的四个分量,分别表示第一个极化方向得到的在x、y轴方向的电场值,分别表示第二个极化方向得到的在x、y方向上的电场值,分别表示第一个极化方向得到的在x、y方向上的磁场值,分别表示第二个极化方向得到的在x、y轴方向的磁场值;
进一步优选,三维地电几何模型中任一区域的各向异性电导率和相对介电常数张量是根据该区域的控制参数计算而来,如下:
σ=Rz(αS)Rx(αD)Rz(αL)σ′Rz(-αL)Rx(-αD)Rz(-αS)
εr=Rz(αS)Rx(αD)Rz(αL)εr′Rz(-αL)Rx(-αD)Rz(-αS)
式中,σ、εr分别表示各向异性电导率和相对介电常数张量,αS、αD、αL为主轴坐标轴xyz到观测坐标系x″y″′z′三个旋转角,Rz(α)和Rx(α)为基本初等旋转矩阵,σ′和ε′r为对角矩阵,如下:
有益效果
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为本发明实例提供的各向异性梯形山脉模型;
具体实施方式
下面将结合实施例对本发明做进一步的说明。
本发明针对现有三维各向异性射频大地电磁模拟技术的缺失,提供了一种高效、高精度地求解三维各向异性射频大地电磁自适应有限元正演方法,本发明能够为快速、可靠的三维射频大地电磁数据解释提供核心的正演引擎。本发明的一种三维各向异性射频大地电磁自适应有限元正演方法是依据下述原理进行的。
1、关于3D各向异性RMT边值问题:
建立三维地电几何模型Ω(表示求解区域),取时间因子e-iωt,在求解区域Ω中麦克斯韦方程组的微分形式可表示为:
其中,E、H和J分别为电场强度、磁场强度和电流密度矢量,μ0、ε0分别表示真空中的磁导率及介电常数,μ0=4π×10-7H/m和ε0≈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的计算公式如下:
σ=Rz(αS)Rx(αD)Rz(αL)σ′Rz(-αL)Rx(-αD)Rz(-αS) (5)
εr=Rz(αS)Rx(αD)Rz(αL)εr′Rz(-αL)Rx(-αD)Rz(-αS) (6)
其中,σ′和ε′r为对角矩阵,Rz(α)和Rx(α)为基本初等旋转矩阵(Pek,2002),具体表达式如下:
将方程(1)式带入(2)式中消去磁场强度H可以获得如下三维各向异性RMT问题在求解区域Ω中满足的电场双旋度方程,如下式(11):
为了保证方程(11)求解的唯一性,根据(1)式,推导出如下电场在求解区域边界上应满足的第二类边界条件:
其中,Ω表示求解区域(含地下和空气)即为初始构建的三维地电模型,表示边界处的单位外法向量,H0为各向异性均匀半空间或层状介质的解析解,具体的计算方式参考(徐震寰,2015)。方程(11)和(12)式便构成了3D各向异性RMT问题所满足的边值问题,下面将采用Garlerkin有限元法对其进行近似地求解。
2、3D RMT边值问题对应的积分弱形式
通过(11)式可以获得残差矢量R的具体表达式为:
根据加权余量方法,选取测试函数V,令其残差在计算区域Ω内的加权为0:
∫∫∫ΩV·Rdv=0 (14)
将残差项(13)式带入上式(14)有:
采用第一矢量格林定理(Zienkiewiczetal.,2013)并带入边界条件(12)式,我们可以获得3D各向异性RMT问题的积分弱形式:
3、Galerkin有限元分析
采用开源四面体生成器Tetgen(Si,2015)将求解区域Ω剖分成由M个互不相交的四面体单元,构成M个四面体单元边的总数为N(一个四面体单元仅处于求解区域Ω中电导率和相对介电常数的控制参数不变的一个区域,不跨区域)。并基于求解区域Ω中不同区域的模型电导率和相对介电常数的控制参数对每个四面体单元赋上电导率和相对介电常数值。然后对每个单元采用1st order Nédélec基函数可以将计算区域中任意点的电场E表示为:
其中,Nj代表第j条边的1st orderNédélec基函数,(若棱边不属于待求电场点所在的四面体单元,则Nj为零。)Ej代表第j条棱边的近似电场切向分量。基于Galerkin理论,取测试函数V为基函数,即:
V=Ni,i=1,...,N (18)
将(18)式带入方程(16)可获得如下大型稀疏线性方程组:
Ax=b (19)
其中,x=(E1,E2,...,EN)T为所要求解的近似电场切向分量,大型稀疏矩阵A中元素的具体表达式为:
对于上式中的体积分项可以借助于(Jin,2014,P268-269)的公式进行解析计算。
线性方程组(19)式的右端项b中元素的具体表达式为:
对于上式中的面积分我们采用高斯数值积分进行计算。通过调用软件包PARDISO(Schenk,2004)中的LU分解法对线性方程组(19)式进行求解可以获得高精度的电场解。
4、正演网格自适应加密
众所周知,1st orderNédélec矢量有限元的解在很大程度上依赖于正演网格的密度。随着正演网格密度的增加,解的精度将随之增加,直到收敛于真解。但采用全局加密的方式会造成很多不必要的细化,而为了尽量降低不必要的网格加密带来的时间和内存消耗,下面我们将基于单元后验误差对正演网格进行自适应地加密。
因为网格离散化带来的误差将造成电流密度在两个物性相同四面体单元的分界面处不连续,从而违背了真实的物理内部边界条件。因此,我们可以将四面体单元三角面上电流密度不连续性的定量估计作为单元后验误差:
其中,代表与四面体单元的Fi面相邻的四面体单元的导纳率张量,E-代表与四面体单元的Fi面相邻的四面体单元上点的电场值,和E+分别表示当前四面体单元上点的导纳率张量、电场值。其中,是根据公式计算的,四面体单元的各向异性电导率和相对介电常数张量σ、εr是根据公式(5)、(6)计算出的。E+、E-是根据公式计算的,若j边不属于该点所在的四面体单元,此时Nj为零。相对误差指示因子βe:
其中代表所有单元中最大的后验误差。预先给定一个相对单元误差允许的最大值β∈(0,1],然后对满足βe>β的单元进行加密直到满足网格加密迭代结束条:迭代次数达到预先设定的允许的最大迭代次数tmax或单元总数达到预先设定的允许的最大单元数nmax。
5、有限元后处理
为了计算出3D RMT实际应用中所需要的阻抗张量、倾子矢量进而获得视电阻率和相位参数,我们需要在地表引入两种不同极化方向的入射场从而获得两种不同的边界条件值,进而获得线性方程组(19)式的两种不同右端项和然后通过代入两种不同的右端项和于线性方程组(19)式中便可获得两种不同极化方向上任意观测点的电场值E1和E2,再根据电磁感应定律公式(1)便可求出观测点处相对应的磁场H1和H2,进而可以获得阻抗张量响应:
倾子矢量(垂直磁场转换函数)响应:
和视电阻率和相位响应:
基于上述原理,本发明提供的一种三维各向异性射频大地电磁自适应有限元正演方法,包括如下步骤:
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分别为105Ωm,103Ωm和106Ωm,对应的旋转角为αS=30°,αL=60°,αD=90°,相对介电常数为主轴各向异性,主值分别模型的求解区域为Ω=[-200m,200m]3。11个观测点沿y=[-50m,50m]以10m为间距对称地布设于地表,观测频率为100kHz。
图4为网格自适应加密过程中第1次、第3次和第8次网格及其单元相对误差示意图。图5为视电阻率相对误差和相位绝对误差随网格自适应加密的收敛曲线图,其中视电阻率相对误差和相位绝对误差分别为所有测点视电阻率ρxy和ρyx相对误差的平均值以及所有测点相位和绝对误差的平均。该模型的解析解为ρxy=2191.96Ωm,ρyx=6653.93Ωm,从图中可以看出,随着网格的自适应加密,视电阻率相对误差和相位绝对误差在逐渐减小,也就是说我们的有限元解在逐渐逼近真解,从而充分验证了本算法对求解各向异性射频大地电磁问题的有效性。
图6为建立的三维各向异性梯形山脉模型,其电导率为主轴各向异性,各主值分别为σx=103Ωm,σy=104Ωm和σz=105Ωm,相对介电常数表现为主轴各向异性,其主值为梯形山脉的上底为边长为45m的正方形,其中心点坐标为(0,0,-45m),下底为边长为200m的正方形,其中心点位于原点,模型求解区域为Ω=[-1km,1km]3。19个测点沿y=[-150m,150m]对称地分布在梯形山脉两侧,在y=[-150m,-30m]和y=[30m,150m]之间的测点的间距为20m,其余测点之间的间距为10m,测量频率为50kHz。
图7为采用自适应有限元算法获得图5的各向异性梯形山脉模型的视电阻率和相位结果图。从图中可以看出由于地形的存在视电阻率从整体上表现为低阻异常,且主要受x方向的电导率影响,而主要受y方向的电导率影响;从视电阻率图中可大致分辨出山脉地形的边界。
需要强调的是,本发明所述的实例是说明性的,而不是限定性的,因此本发明不限于具体实施方式中所述的实例,凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式,不脱离本发明宗旨和范围的,不论是修改还是替换,同样属于本发明的保护范围。
Claims (8)
1.一种三维各向异性射频大地电磁自适应有限元正演方法,其特征在于:包括如下步骤:
S1:构建表示求解区域的三维地电几何模型,并获取三维地电几何模型中不同区域的控制参数;
所述控制参数为各个区域的电导率主值参数、相对介电常数主值参数、对应的主轴坐标系与观测坐标系之间的旋转角;
S2:将三维地电几何模型剖分为若干个互不相交的四面体单元;
其中,一个四面体单元仅属于三维地电几何模型中的一个区域;
S3:根据麦克斯韦方程组和Garlerkin有限元原理获取三维各向异性射频大地电磁对应的积分弱形式,并基于所述积分弱形式、1st order基函数、四面体单元每条棱边的近似电场切向分量构建稀疏线性方程组;
S4:在地表引入两种不同极化方向的入射场,得到所述稀疏线性方程组中的两种不同右端项,并求解两种含不同边界条件值的稀疏线性方程组,获得两种极化方向上分别对应的每条棱边的近似电场切向分量;
S5:基于两种极化方向上分别对应的每条棱边的近似电场切向分量获取测点处的目标参量,所述目标参量为:电场、磁场、阻抗张量、磁倾子矢量、视电阻率、相位响应中一个或多个的组合。
4.根据权利要求1所述的方法,其特征在于:步骤S4中计算第一个极化方向上每条棱边的近似电场切向分量的过程如下:
S41:基于第一个极化方向入射场获得的边界条件值求解稀疏线性方程组得到每条棱边的近似电场切向分量;
S42:基于步骤S41计算出的每条棱边的近似电场切向分量计算每个四面体单元的后验误差,再基于后验误差计算出每个四面体单元的相对误差指示因子;
S43:分别判断每个四面体单元的相对误差指示因子是否大于预设最大误差,若大于,则对所述四面体单元二次剖分,并返回步骤S3,直至满足迭代终止条件;
所述迭代终止条件为:迭代次数达到预设最大迭代次数或求解区域的四面体单元总数达到预设最大单元数;
其中,最后得到的所有四面体单元棱边的近似电场切向分量为所需的第一个极化方向上每条棱边的近似电场切向分量。
5.根据权利要求4所述的方法,其特征在于:步骤S4中计算第二个极化方向上每条棱边的近似电场切向分量的过程如下:
采用第一个极化方向上最后得到的所有四面体单元以及第二个极化方向上的边界条件值求解稀疏线性方程组得到第二个极化方向上每条棱边的近似电场切向分量。
7.根据权利要求1所述的方法,其特征在于:步骤S5中所述阻抗张量、磁倾子矢量、视电阻率、相位响应的计算公式如下:
式中,Zxx、Zxy、Zyx、Zyy分别表示阻抗张量Z的四个分量,分别表示第一个极化方向得到的在x、y方向上的电场值,分别表示第二个极化方向上得到的在x、y方向上的电场值,分别表示第一个极化方向得到的在x、y轴方向的磁场值,分别表示第二个极化方向得到的在x、y轴方向的磁场值;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910457480.7A CN110058315B (zh) | 2019-05-29 | 2019-05-29 | 一种三维各向异性射频大地电磁自适应有限元正演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910457480.7A CN110058315B (zh) | 2019-05-29 | 2019-05-29 | 一种三维各向异性射频大地电磁自适应有限元正演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110058315A CN110058315A (zh) | 2019-07-26 |
CN110058315B true CN110058315B (zh) | 2020-04-14 |
Family
ID=67325001
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910457480.7A Expired - Fee Related CN110058315B (zh) | 2019-05-29 | 2019-05-29 | 一种三维各向异性射频大地电磁自适应有限元正演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110058315B (zh) |
Families Citing this family (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111611737B (zh) * | 2020-05-19 | 2022-05-20 | 中南大学 | 一种三维任意各向异性介质的海洋可控源电磁正演方法 |
CN111638556B (zh) * | 2020-06-09 | 2022-12-27 | 东华理工大学 | 基于地空分解策略的大地电磁正演方法及装置、存储介质 |
CN111898294B (zh) * | 2020-07-09 | 2023-06-27 | 长安大学 | 一种电偶极源多频三维有限元数值模拟方法 |
AU2021101629A4 (en) * | 2021-02-25 | 2021-05-20 | Institute Of Geology And Geophysics, Chinese Academy Of Sciences | Boundary truncation layer method and apparatus for three-dimensional forward modelling of low-frequency |
CN113221411B (zh) * | 2021-05-08 | 2022-12-09 | 桂林理工大学 | 针对任意形状有损耗介质充电电位数值模拟方法、系统、终端 |
CN113792445B (zh) * | 2021-11-15 | 2022-02-08 | 中南大学 | 一种基于积分方程法的三维大地电磁数值模拟方法 |
CN115906559B (zh) * | 2022-10-31 | 2023-08-15 | 重庆大学 | 一种基于混合网格的大地电磁自适应有限元正演方法 |
CN115935671B (zh) * | 2022-12-20 | 2023-08-18 | 安徽大学 | 一种区域分解电磁仿真方法及系统 |
CN116341332A (zh) * | 2023-03-30 | 2023-06-27 | 重庆大学 | 基于电导率分块连续变化的大地电磁三维有限元正演方法 |
CN116842813B (zh) * | 2023-09-04 | 2023-11-14 | 中南大学 | 一种三维大地电磁正演数值模拟方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102930071A (zh) * | 2012-08-29 | 2013-02-13 | 电子科技大学 | 基于非匹配网格的周期结构的三维电磁场仿真模拟方法 |
CN105204073A (zh) * | 2015-09-18 | 2015-12-30 | 中南大学 | 一种张量视电导率测量方法 |
CN106324689A (zh) * | 2016-06-24 | 2017-01-11 | 杭州迅美科技有限公司 | 一种水平井地层环境下电阻率各向异性识别方法 |
CN106980736A (zh) * | 2017-04-11 | 2017-07-25 | 吉林大学 | 一种各向异性介质的海洋可控源电磁法有限元正演方法 |
CN107748394A (zh) * | 2017-09-30 | 2018-03-02 | 中南大学 | 一种针对rmt数据的双参数反演算法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20090083006A1 (en) * | 2007-09-20 | 2009-03-26 | Randall Mackie | Methods and apparatus for three-dimensional inversion of electromagnetic data |
-
2019
- 2019-05-29 CN CN201910457480.7A patent/CN110058315B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102930071A (zh) * | 2012-08-29 | 2013-02-13 | 电子科技大学 | 基于非匹配网格的周期结构的三维电磁场仿真模拟方法 |
CN105204073A (zh) * | 2015-09-18 | 2015-12-30 | 中南大学 | 一种张量视电导率测量方法 |
CN106324689A (zh) * | 2016-06-24 | 2017-01-11 | 杭州迅美科技有限公司 | 一种水平井地层环境下电阻率各向异性识别方法 |
CN106980736A (zh) * | 2017-04-11 | 2017-07-25 | 吉林大学 | 一种各向异性介质的海洋可控源电磁法有限元正演方法 |
CN107748394A (zh) * | 2017-09-30 | 2018-03-02 | 中南大学 | 一种针对rmt数据的双参数反演算法 |
Non-Patent Citations (4)
Title |
---|
3-D direct current resistivity anisotropic modelling adaptive finite element methods by goal-oriented adaptive finite element methods;Zhengyong Ren et al.;《Geophysical Journal International》;20181231(第212期);第76–87页 * |
A hybrid finite-element and integral-equation method for forward modeling of 3D controlled-source electromagnetic induction;Zhou Feng et al.;《APPLIED GEOPHYSICS》;20181231;第15卷(第3-4期);第536–544页 * |
一种新的三维直流电阻率体积分方程方法;陈煌 等;《地球物理学报》;20171130;第60卷(第11期);第4506-4515页 * |
基于电流密度连续性条件的直流电阻率各向异性问题自适应有限元模拟;任政勇 等;《地球物理学报》;20180131;第61卷(第1期);第331-343页 * |
Also Published As
Publication number | Publication date |
---|---|
CN110058315A (zh) | 2019-07-26 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110058315B (zh) | 一种三维各向异性射频大地电磁自适应有限元正演方法 | |
Key et al. | Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example | |
Li et al. | A finite-element time-domain forward solver for electromagnetic methods with complex-shaped loop sources | |
CN107742015B (zh) | 基于任意偶极-偶极装置的直流激电法三维数值模拟方法 | |
Grayver et al. | 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO 2 storage formation | |
Yin et al. | A goal-oriented adaptive finite-element method for 3D scattered airborne electromagnetic method modeling | |
CN106199742B (zh) | 一种频率域航空电磁法2.5维带地形反演方法 | |
Yin et al. | 3D time-domain airborne EM forward modeling with topography | |
Liang et al. | A new inversion method based on distorted born iterative method for grounded electrical source airborne transient electromagnetics | |
CN112949134A (zh) | 基于非结构有限元方法的地-井瞬变电磁反演方法 | |
CN104656156A (zh) | 音频大地电磁测深三维采集资料的磁参考处理方法 | |
CN110346835A (zh) | 大地电磁的正演方法、正演系统、存储介质及电子设备 | |
CN106443803B (zh) | 基于发射装置实测形态数据的海洋可控源电磁响应计算方法 | |
CN104408021A (zh) | 一种电偶源三维时域有限差分正演成像方法 | |
CN111638556A (zh) | 基于地空分解策略的大地电磁正演方法及装置、存储介质 | |
Wang et al. | 2D joint inversion of CSAMT and magnetic data based on cross-gradient theory | |
CN110068873A (zh) | 一种基于球坐标系的大地电磁三维正演方法 | |
CN115238550A (zh) | 自适应非结构网格的滑坡降雨的地电场数值模拟计算方法 | |
Chen et al. | Three-dimensional inversion of controlled-source audio-frequency magnetotelluric data based on unstructured finite-element method | |
Li et al. | 3-D marine CSEM forward modeling with general anisotropy using an adaptive finite-element method | |
Zhou et al. | A hybrid finite-element and integral-equation method for forward modeling of 3D controlled-source electromagnetic induction | |
CN113297526B (zh) | 一种基于Wenner四极和大地电磁数据的水平分层土壤结构联合反演方法 | |
Dong et al. | An approach to model earth conductivity structures with lateral changes for calculating induced currents and geoelectric fields during geomagnetic disturbances | |
Oh et al. | Interpretation of controlled-source electromagnetic data from iron ores under rough topography | |
CN115238566A (zh) | 一种基于频域电磁响应的人工空洞识别方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200414 |