CN111638556B - 基于地空分解策略的大地电磁正演方法及装置、存储介质 - Google Patents

基于地空分解策略的大地电磁正演方法及装置、存储介质 Download PDF

Info

Publication number
CN111638556B
CN111638556B CN202010515462.2A CN202010515462A CN111638556B CN 111638556 B CN111638556 B CN 111638556B CN 202010515462 A CN202010515462 A CN 202010515462A CN 111638556 B CN111638556 B CN 111638556B
Authority
CN
China
Prior art keywords
magnetic field
node
phi
scalar
psi
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.)
Active
Application number
CN202010515462.2A
Other languages
English (en)
Other versions
CN111638556A (zh
Inventor
周峰
任政勇
陈煌
邓居智
张志勇
汤井田
李广
蓝泽鸾
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
East China Institute of Technology
Original Assignee
East China Institute of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by East China Institute of Technology filed Critical East China Institute of Technology
Priority to CN202010515462.2A priority Critical patent/CN111638556B/zh
Publication of CN111638556A publication Critical patent/CN111638556A/zh
Application granted granted Critical
Publication of CN111638556B publication Critical patent/CN111638556B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/081—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 the magnetic field is produced by the objects or geological structures

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Electromagnetism (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Measuring Magnetic Variables (AREA)

Abstract

本发明公开了基于地空分解策略的大地电磁正演方法及装置、存储介质。该方法包括：包括：依据数字地形数据(DEM)来构建求解区域三维MT的地电几何模型；将地电几何模型离散为一系列四面体单元；获取基于地空分解策略的A‑Φ‑Ψ耦合势三维MT正演满足的积分弱解形式并对积分弱分解形式进行基于节点离散的转换得到稀疏线性方程组；加载入射场置于稀疏线性方程组的右端项并求解所述稀疏线性方程得到求解区域中四面体单元各个节点的矢量位A、标量位Φ，标量位Ψ；根据矢量位A、标量位Φ，标量位Ψ计算出观测点的参数。本发明不仅对开展大规模三维大地电磁法反演研究具有推动作用，对于促进学科发展，提升电磁法资料解释的学术水平具有重要现实意义与价值。

Description

基于地空分解策略的大地电磁正演方法及装置、存储介质
技术领域
本发明属于应用地球物理技术领域，具体涉及一种基于地空分解策略的大地电磁正演方法及装置、存储介质。
背景技术
大地电磁法(MT)是一种利用天然交变电磁场源来探测地球内部电性结构的一种地球物理勘探方法，被广泛应用于矿产资源勘查、油气与天然气资源勘探、地热资源勘探以及深部地质构造的探测等领域的研究，并取得了良好的应用效果。
目前，任意复杂的野外勘探环境会直接影响大地电磁法资料解释的准确度，如由平缓地区向复杂山区过渡的勘探区域；其次，大规模的三维数据集涌现使得现有的MT资料解释技术面临着内存消耗大等难题。除此之外，基于现有的三维MT正演求解公式在求解大尺度问题时，正演求解公式离散后得到的系统矩阵面临着条件数差的难题，严重影响了线性方程组求解效率。针对上述问题，要实现应用于深部矿产资源勘探、油气与天然气资源的勘探的大规模、任意复杂地形下的三维MT数据集的资料解释技术，快速、高精度的正演技术是其前提条件。因此，有必要寻找一种高精度、高效三维大地电磁法正演方法。
发明内容
本发明的目的是提供一种基于地空分解策略的大地电磁正演方法，该技术是基于地空分解策略A-Φ-Ψ耦合势实现了三维MT自适应有限元正演，可以实现高效、高精度正演计算，该技术克服传统MT求解系统面临的条件数差、求解效率低以及精度差等难题，具有高效率、高精度以及高稳定性的能力，为实现大尺度三维大地电磁法实测数据反演解释提供了核心动力。
一方面，本发明提供的一种基于地空分解策略的大地电磁正演方法，包括如下步骤：
步骤1：依据数字地形数据(DEM)来构建求解区域三维MT的地电几何模型；
步骤2：将所述地电几何模型离散为一系列不重叠的四面体单元；
步骤3：获取基于地空分解策略的A-Φ-Ψ耦合势三维MT正演满足的积分弱解形式，并采用节点基函数和四面体单元的节点对所述积分弱分解形式进行基于节点离散的转换得到稀疏线性方程组，其中，A表示磁场矢量位、Φ表示电场标量位，Ψ表示磁场标量位，节点为离散四面体单元之间的连接枢纽；
步骤4：加载入射场置于稀疏线性方程组的右端项并求解所述稀疏线性方程得到求解区域中四面体单元各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ；其中，空气区域内节点的磁场矢量位A、电场标量位Φ为0，地下区域内节点的磁场标量位Ψ为0；
步骤5：根据所述磁场矢量位A、电场标量位Φ，磁场标量位Ψ计算出观测点的电场、磁场一个或多个参数的组合。
进一步优选，步骤3中所述积分弱解形式如下：
Figure BDA0002529922970000021
式中，V表示测试函数，E表示电场，B为磁场,Ω表示求解区域，
Figure BDA00025299229700000211
表示求解区域的边界，Ω0为空气空间区域,Ω1为地下空间区域,Γ0表示空气层的上顶面，Γ1表示地下区域的下底面，Γ为地下区域和空气区域的分界面，n表示边界外法向量，dv、ds分别表示体积微分元和面积微分元，v表示体积、s表示面积，
Figure BDA0002529922970000022
表示阻抗率，
Figure BDA0002529922970000023
表示导纳率，μ为磁导率。
进一步优选，步骤3中所述稀疏线性方程组如下：
Figure BDA0002529922970000024
式中，
Figure BDA0002529922970000025
Af,Am均为稀疏的子矩阵，
Figure BDA0002529922970000026
均为地空分界面的面积分项形成的子矩阵,
Figure BDA0002529922970000027
均表示零矩阵，B0,B1,B2,B3,B4均为稀疏线性方程组的右端项，X0、X1、X2、X3、X4分别为稀疏线性方程组中的未知数矩阵，且由求解区域内所有节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ构成。
进一步优选，所述稀疏线性方程组中的参数表示如下：
Figure BDA0002529922970000028
Figure BDA0002529922970000029
Figure BDA00025299229700000210
Figure BDA0002529922970000031
Figure BDA0002529922970000032
Figure BDA0002529922970000033
Figure BDA0002529922970000034
Figure BDA0002529922970000035
Figure BDA0002529922970000036
Figure BDA0002529922970000037
Figure BDA0002529922970000038
Figure BDA0002529922970000039
Figure BDA00025299229700000310
Figure BDA00025299229700000311
Figure BDA00025299229700000312
Figure BDA00025299229700000313
Figure BDA00025299229700000314
Figure BDA00025299229700000315
Figure BDA00025299229700000316
Figure BDA00025299229700000317
Figure BDA00025299229700000318
Figure BDA00025299229700000319
式中，i,j＝1～n；m,k＝1～n1，
Figure BDA00025299229700000320
Figure BDA00025299229700000321
表示笛卡尔坐标系三个方向，
Figure BDA00025299229700000322
分别表示对x求导，
Figure BDA00025299229700000323
分别表示对y求导，
Figure BDA00025299229700000324
分别表示对z求导，[·]x表示x方向上的分量，
Figure BDA00025299229700000325
表示地下区域第一个节点，第二个节点、第n节点在x方向的磁场矢量位Ax；
Figure BDA00025299229700000326
表示地下区域第一个节点，第二个节点、第n节点在y方向的磁场矢量位Ay；
Figure BDA0002529922970000041
表示地下区域第一个节点，第二个节点、第n节点在z方向的磁场矢量位Az，x、y、z方向的磁场矢量位Ax、Ay、Az构成磁场矢量位A；Φ1、Φ2、Φn表示地下区域第一个节点，第二个节点、第n节点的电场标量位；Ψ1、Ψ2、
Figure BDA0002529922970000042
表示空气区域第一个节点，第二个节点、第n1节点的磁场标量位，T表示矩阵转置，n为地下区域的节点总数，n1为空气区域的节点总数，节点对应为四面体单元的顶点；
Li,Lj分别表示地下区域内节点i，节点j对应的节点基函数，分别表示空气区域中第k个、第m个节点对应的节点基函数，
Figure BDA0002529922970000043
表示子矩阵
Figure BDA0002529922970000044
中节点i和节点j对应的元素，
Figure BDA0002529922970000045
表示子矩阵B0中节点i对应的元素，其他子矩阵的定义类似，在此不多赘述。
Figure BDA0002529922970000046
表示阻抗率，
Figure BDA0002529922970000047
表示导纳率，Ω0为空气空间区域，Ω1为地下空间区域，Γ0表示空气层的上顶面，Γ1表示地下区域的下底面，Γ为地下区域和空气区域的分界面，n表示边界外法向量，dv、ds分别表示体积微分元和面积微分元，v表示体积、s表示面积，μ为磁导率。
Figure BDA0002529922970000048
分别表示极化场源在均匀半空间或层状介质地电模型下激发后产生的背景矢量位A的三个分量，Eb,Bb分别表示极化场源在均匀半空间或层状介质地电模型下激发后产生的背景电场和磁场。因此，本发明加载入射场利用上述公式可以计算出系数线性方程组的右端项。
进一步优选，步骤4中加载入射场置于稀疏线性方程组的右端项并求解所述稀疏线性方程得到求解区域中四面体单元各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ的过程如下：
①：构建基于地空分解策略的A-Φ-Ψ耦合势三维MT对偶问题，并求解对偶问题得到对偶场WA-Wφ-WΨ；以及基于加载的入射场计算出所述稀疏线性方程组的右端项，并求解稀疏线性方程组得到求解区域内当前四面体单元离散对应的节点的磁场矢量位A、电场标量位Φ和磁场标量位Ψ；
其中，将对偶问题转换为线性方程组，并由三分量单位源积分表达式求解得到线性方程组的右端项，再求解线性方程组得到对偶场WA-WΦ-Wψ，WA表示磁场矢量位A的对偶场，WΦ表示为电场标量位Φ的对偶场，Wψ为磁场标量位Ψ的对偶场；
②：基于步骤①计算出的每个节点上的磁场矢量位A、电场标量位Φ和磁场标量位Ψ以及磁场矢量位A的对偶场WA、电场标量位Φ的对偶场WΦ、磁场标量位Ψ的对偶场Wψ计算出大地电磁法问题的对偶问题分别对应的单元误差分布，再基于单元误差分布计算出四面体单元的相对误差单元指示因子；
③：分别判断每个四面体单元的相对误差指示因子是否大于预设最大误差，若大于预设最大误差，对所述四面体单元进行细化离散，再返回步骤3，直至满足迭代终止条件；
所述迭代终止条件为：迭代次数达到预设最大迭代次数或求解区域的四面体单元总数达到预设最大单元数；
其中，最后所得到的所有四面体单元节点上的磁场矢量位A、电场标量位Φ和磁场标量位Ψ。
其中，预设最大误差、预设最大迭代次数以及预设最大单元数可以是实际需求实验得到的经验值。
进一步优选，大地电磁法问题的对偶问题分别对应的单元误差分布如下所示：
Figure BDA0002529922970000051
Figure BDA0002529922970000052
式中，e、w分别表示大地电磁法问题、对偶问题分别对应的单元误差分布，n表示边界外法向量，I为虚数单元，ω为角频率，ds表示面积微分元，s表示面积，μ为磁导率，
Figure BDA0002529922970000053
为四面体的表面，符号[·]表示场在面
Figure BDA0002529922970000054
的跳跃值；
所述相对误差单元指示因子如下：
Figure BDA0002529922970000055
其中，γ＝|ew|
式中，γe为四面体单元的相对误差单元指示因子，γ为四面体单元的加权误差，
Figure BDA0002529922970000056
为所有四面体单元中加权误差的最大值。
进一步优选，步骤4中加载入射场置于稀疏线性方程组的右端项分为两个极化方向分别进行入射场加载得到两个极化方向对应的各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ，所述方法还包括：
基于两个极化方向对应的各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ得到阻抗张量、视电阻率、相位响应中一个或多个的组合。
其中，在第一个极化方向按照上述①-③得到各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ后。第二个极化方向同样可以按照上述①-③得到各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ后。此外，第二个极化方向，还可以第一个极化方向自适应加密后的正演网格直接进行计算，这是考虑到同一求解区域在不同极化方向上正演网格的差异不大。
进一步优选，对所述稀疏线性方程进行求解之前还包括获取所述地电几何模型中不同区域的物性参数；
其中，所述物性参数包括各个区域的电阻率参数、介电常数参数以及磁导率参数。
另一方面，本发明提供的一种上述方法的装置，包括：
模型构建模块：用于构建求解区域三维MT的地电几何模型；
四面体单元离散模块：用于将地电几何模型离散为一系列不重叠的四面体单元；
运算模块：用于基于地空分解策略的A-Φ-Ψ耦合势进行正演计算，最终得到观测点的目标参数。
此外，本发明提供一种存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述方法的步骤。
有益效果
1、本发明提供的一种基于地空分解策略的大地电磁正演方法是基于地空分解策略的A-Φ-Ψ耦合势三维MT有限元算法，该算法所在的节点有限元空间进行，相对于传统的电场方程的矢量有限元来说，不存在空解的现象。同时节点有限元形成的矩阵条件数相比于矢量有限元来说，矩阵条件数的值更加均衡，方程组求解更为稳定。本发明提出的基于地空分解策略的A-Φ-Ψ耦合势三维MT有限元算法能够开展迭代求解，而电场方程的矢量有限元不适合迭代计算，只能适合直接求解器。迭代求解器与直接求解器最为明显的优势在于计算机内存需求小，同时不需要系统矩阵进行求逆，节约时间，提高正演求解效率以及节约了内存消耗。因此，本发明提供的基于地空分解策略的A-Φ-Ψ耦合势正演求解系统能够适合大规模三维大地电磁法正演的快速求解，对三维大地电磁法反演研究具有推动作用，且对于促进学科发展，提升电磁法资料解释的学术水平也具有重要现实意义与价值。
2、本发明依据三维MT有限元数值精度主要依赖于网格剖分，传统技术对于简单地电模型可以满足设计要求，对复杂起伏地形的地电模型来说，网格的设计会消耗大量时间，往往得不到理想的效果，为此本发明优选方案中依据数字地形数据(DEM)来构建三维MT所需要的地电几何模型，采用非结构化四面体网格剖分器实现求解区域的离散，得到一系列的不重叠的非结构化四面体网格。
3、本发明的优选方案中构建基于地空分解策略的A-Φ-Ψ耦合势的对偶问题，依据单元加权误差来实现后验误差估计，从而可以实现基于地空分解策略的A-Φ-Ψ耦合势自适应MT正演算法，该自适应加密算法为了考虑目标区或测点处场值的计算精度，能够在正演计算过程中不断调整求解区域的网格剖分密度，对非目标区或测点处进行粗化处理，对目标区和测点区进行不断细化处理，该技术过程弥补了常规有限元不能有效提高测点处场值的求解精度的缺陷，同时也可避免过度网格加密带来计算资源的浪费。
附图说明
图1为MT勘探示意图。其中，Γ0表示空气层的上顶面，Γ1表示地下区域的下底面，Γ为起伏地表，n表示外法向，Ω0表示为空气空间，Ω1表示为地下空间。
图2为基于地空分解策略的A-Φ-Ψ耦合势自适应三维大地电磁法正演流程图；
图3为几何地电模型示意图，其中，(a)图为x-y平面的示意图，(b)图为x-z剖面的示意图；
图4为三维地形模型网格示意图；
图5为频率为2Hz，地形模型得到的ρxx和ρxy视电阻率等值线图，其中，(a)图为ρxx视电阻率等值线图，(b)图为ρxy视电阻率等值线图；
图6为频率为2Hz，地形模型得到的ρyx和ρyy视电阻率等值线图，其中，(a)图为ρyx视电阻率等值线图，(b)图为ρyy视电阻率等值线图。
具体实施方式
下面将结合实施例对本发明做进一步的说明。
本发明针对传统MT求解系统面临的条件数差、求解效率低以及精度差等难题，提供了一种基于地空分解策略的大地电磁正演方法。该技术是基于地空分解策略A-Φ-Ψ耦合势实现了3D MT自适应有限元正演模拟技术，可以在任意三维复杂地形MT问题上实现高效、高精度正演计算，为实现大尺度三维大地电磁法实测数据反演解释提供了核心动力。其是依据如下原理进行的：
1、关于3D各向同性MT边值问题：
从大地电磁法满足的Maxwell方程出发，导出磁场矢量位A和电场标量位Φ耦合的矢量Helmholtz总场方程如下：
Figure BDA0002529922970000081
其中，
Figure BDA0002529922970000082
Figure BDA0002529922970000083
表示导纳率，
Figure BDA0002529922970000084
表示阻抗率，i为虚数单位，A表示磁场矢量位，Φ和Ψ分别为电场标量位、磁场标量位，σ为电导率(电阻率的倒数)，ω＝2πf为角频率，f为频率，ε≈8.854187817×10-12F/m，ε为介电常数，μ0真空磁导率。
众所周知，在MT频段，可假设不存在位移电流，在地球外部(即空气空间)，磁场B表示为：
Figure BDA0002529922970000085
由公式(2)可知，磁场标量位Ψ满足下面的方程，
Figure BDA0002529922970000086
磁场标量位Ψ只在空气区域进行求解，根据磁场连续条件，磁场矢量位A和磁场标量位Ψ在空气界面上满足下面的方程，
Figure BDA0002529922970000087
同时，在空气中，磁场矢量位A和电场
Figure BDA0002529922970000088
法向不连续，
Figure BDA0002529922970000089
(J为电流密度)以及磁场矢量位A在地空边界消失，因此，
Figure BDA00025299229700000810
其中，Γ表示地空界面，如图1所示。
基于地-空区域分解的A-Φ-Ψ的3D MT问题满足的微分方程组为：
Figure BDA00025299229700000811
为了确保微分方程组解的唯一性，其需要结合下面的边界条件,
Figure BDA00025299229700000812
其中，Γ为地下区域和空气区域的分界面，即边界，n表示边界外法向量。
2、3D各向同性MT边值问题对应的积分弱解形式。
通过(6)式可以获得残差矢量RΨ,RA,RΦ的具体表达式为：
Figure BDA0002529922970000091
根据加权余量方法，选取测试函数V，令其残差在计算区域Ω(求解区域)内的加权为0：
Figure BDA0002529922970000092
将残差项(8)带入到(9)，可得基于地-空区域分解的A-Φ-Ψ的三维MT满足的微分方程转变成积分弱解形式，如下：
Figure BDA0002529922970000093
Figure BDA0002529922970000094
Figure BDA0002529922970000095
式中，V表示测试函数，E表示电场，B为磁场,Ω表示求解区域，
Figure BDA0002529922970000096
表示求解区域的边界，Γ0表示空气层的上顶面，Γ1表示地下区域的下底面，dv、ds分别表示体积微分元和面积微分元(v表示体积、s表示面积)。
对上述积分弱解形式进行转换，得到如下所示的稀疏线性方程组：
Figure BDA0002529922970000097
式中，
Figure BDA0002529922970000098
Af,Am均为稀疏的子矩阵，
Figure BDA0002529922970000099
均为地空分界面的面积分项形成的子矩阵，
Figure BDA00025299229700000910
均表示零矩阵，B0,B1,B2,B3,B4均为稀疏线性方程组的右端项。
Figure BDA00025299229700000911
Figure BDA00025299229700000912
Figure BDA0002529922970000101
X3＝(Φ1,Φ2,…Φn)T (17)
Figure BDA0002529922970000102
Figure BDA0002529922970000103
Figure BDA0002529922970000104
Figure BDA0002529922970000105
Figure BDA0002529922970000106
Figure BDA0002529922970000107
Figure BDA0002529922970000108
Figure BDA0002529922970000109
Figure BDA00025299229700001010
Figure BDA00025299229700001011
Figure BDA00025299229700001012
Figure BDA00025299229700001013
Figure BDA00025299229700001014
式中，i,j＝1～n；m,k＝1～n1，
Figure BDA00025299229700001015
Figure BDA00025299229700001016
表示笛卡尔坐标系三个方向，
Figure BDA00025299229700001017
分别表示对x求导，
Figure BDA00025299229700001018
分别表示对y求导，
Figure BDA00025299229700001019
分别表示对z求导，
Figure BDA00025299229700001020
表示地下区域第一个节点，第二个节点、第n节点在x方向的磁场矢量位Ax；
Figure BDA0002529922970000111
表示地下区域第一个节点，第二个节点、第n节点在y方向的磁场矢量位Ay；
Figure BDA0002529922970000112
表示地下区域第一个节点，第二个节点、第n节点在z方向的磁场矢量位Az，x、y、z方向的磁场矢量位Ax、Ay、Az构成磁场矢量位A；Φ1、Φ2、Φn表示地下区域第一个节点，第二个节点、第n节点的电场标量位；Ψ1、Ψ2、
Figure BDA0002529922970000113
表示空气区域第一个节点，第二个节点、第n1节点的磁场标量位，T表示矩阵转置，n为地下区域的节点总数，n1为空气区域的节点总数，节点对应为四面体单元的顶点；
Li,Lj,Lk,Lm表示节点基函数，E表示电场，B为磁场,
Figure BDA0002529922970000114
表示阻抗率，
Figure BDA0002529922970000115
表示导纳率，Ω0为空气空间区域,Ω1为地下空间区域,Γ0表示空气层的上顶面，Γ1表示地下区域的下底面，Γ为地下区域和空气区域的分界面，n表示边界外法向量，dv、ds分别表示体积微分元和面积微分元，v表示体积、s表示面积，μ为磁导率。该混合求解系统的未知数为4n+n1，n为地下区域的节点总数，n1为空气区域的节点总数。为了求解出各个子矩阵的积分稀疏，需要将测试函数用1阶的节点函数代替，即V＝L，L为节点基函数。
其中，针对四面体单元中任意非节点位置，其矢量位A、标量位Φ以及标量位Ψ的计算如下：
Figure BDA0002529922970000116
Figure BDA0002529922970000117
Figure BDA0002529922970000118
其中，Lj表示四面体单元第j个节点对应的节点基函数。若推至n个节点求和，则公式(35)-(37)中存在：若节点不位于该四面体单元上，则此时Lj为0。
最后，采用GMRES迭代求解器对方程组(13)进行计算，得到各个节点的矢量位A以及标量位Φ和Ψ。其中，需要说明的是，上述计算得到的矢量位A以及标量位Φ是针对地下区域的，空气区域中节点的矢量位A以及标量位Φ是为0，上述计算得到的标量位Ψ是针对空气区域，地下区域中节点的标量位Ψ为0。因此，通过上述计算可以获知地下区域和空气区域中所有节点对应的矢量位A以及标量位Φ和Ψ。
3、网格自适应加密
为了实现面向目标自适应有限元技术，需要构建基于地空分解策略的A-Φ-Ψ耦合势3D MT对偶问题。首先在测点区域引入简单的“辅助人工源”，“辅助人工源”不是实际的物理问题，是人为构建的数学问题，其宗旨是用来加密测点区域。采用Ωj表示观测点的区域，则通过在Ωj区域注入三分量的单位源I＝(1,1,1)，来构建大地电磁法问题对应的统一对偶问题：
Figure BDA0002529922970000121
其中WA表示磁场矢量位A的对偶场，WΦ表示为电场标量位Φ的对偶场，WΨ为空气中标量位Ψ的对偶场，δj为delta函数，t为三分量的单位源的个数，观察方程(38)的结构形式，大地电磁法满足A-Φ-Ψ耦合势的总场方程与对偶场WA-Wφ-WΨ具有相同的结构，从而具有相同有限元线性系统矩阵。因此，方程(38)对应的对偶问题线性方程组为：
KY＝C (39)
其中，K是对偶问题形成的刚度矩阵，其与方程(13)稀疏矩阵是一致的，Y为在单元节点上的对偶场，
Figure BDA0002529922970000122
由三分量单位源积分表达式求得。
e为四面体单元上的后验误差估计值，
Figure BDA0002529922970000123
为四面体单元的表面，符号[·]表示场在面
Figure BDA0002529922970000124
的跳跃值。大地电磁法问题和对偶问题分别对应的单元误差分布，e和w：
Figure BDA0002529922970000125
采用加权的思想，最终获得的四面体单元的加权误差γ：
γ＝|ew| (41)
最后，可以估计出当前单元的相对误差单元指示因子
Figure BDA0002529922970000131
式中，
Figure BDA0002529922970000132
表示所有四面体单元中加权误差的最大值。
4、观测点处的场值计算。
为了计算出3D MT实际应用，可以利用本领域的理论公式计算出电场、磁场，如电场公式：
Figure BDA0002529922970000133
此外，本发明为了得到视电阻率和相位等参数，我们需要加载不同极化方向的场源，得到不同极化方向的视电阻率和相位，具体计算观测点处视电阻率和相位的表达式如下：
Figure BDA0002529922970000134
式中，Zxx、Zxy、Zyx、Zyy分别表示阻抗张量Z的四个分量，
Figure BDA0002529922970000135
分别表示第一个极化方向场源加载得到的x、y方向的电场值，
Figure BDA0002529922970000136
分别表示第二个极化方向场源加载得到的x、y方向的电场值，
Figure BDA0002529922970000137
分别表示第一个极化方向场源加载得到的在x、y方向的磁场值，
Figure BDA0002529922970000138
分别表示表示第二个极化方向场源加载得到的x、y方向的磁场值。视电阻率的表达式如下，
Figure BDA0002529922970000139
式中，
Figure BDA00025299229700001310
表示ij方向的视电阻率，ω表示角频率，μ0＝4π×10-7H/m；
Figure BDA00025299229700001311
式中，
Figure BDA00025299229700001312
表示ij方向的相位响应。
基于上述原理，本发明实施例提供的一种基于地空分解策略的大地电磁正演方法，包括如下步骤：
步骤1：依据DEM数据构建求解区域三维MT的地电几何模型，并建立不同网格区域的物性参数；
物性参数为网格区域的电阻率参数、介电常数参数以及磁导率参数，DEM数据是测地卫星、InSAR数据、地形图来获取实测区域Digital Evaluation Map(DEM)数据。
步骤2：采用Tetgen程序将所述地电几何模型离散为一系列不重叠。其中，每一个四面体单元隶属于三维MT地电几何模型中的一个区域。
步骤3：获取基于地空分解策略的A-Φ-Ψ耦合势三维MT正演满足的积分弱解形式，并采用节点基函数和四面体单元的节点对所述积分弱分解形式进行基于节点离散的转换得到稀疏线性方程组，本实施例中，节点为四面体单元的顶点，其他可行的实施例中，可以选用棱边，再基于本领域的理论对本实施例中提供的公式进行相应变化；
步骤4：加载两个不同极化方向的入射场置于稀疏线性方程组的右端项并求解所述稀疏线性方程得到求解区域四面体单元各个节点在两个不同极化方向对应的磁场矢量位A、电场标量位Φ，磁场标量位Ψ；
步骤5：根据所述磁场矢量位A、电场标量位Φ，磁场标量位Ψ计算出观测点的电场、磁场、阻抗张量、视电阻率、相位响应一个或多个参数的组合。
此外，本发明还提供一种实现上述方法的装置，其包括：模型构建模块、四面体单元离散模块以及运算模块，其中，模型构建模块用于构建求解区域三维MT的地电几何模型；四面体单元离散模块用于将地电几何模型离散为一系列不重叠的四面体单元；运算模块用于基于地空分解策略的A-Φ-Ψ耦合势进行正演计算，最终得到观测点的目标参数。具体请参照方法相关内容，本发明对此不进行具体的限定。
本发明还提供一种计算机设备，其包括存储器和处理器，存储器上存储了可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述一种基于地空分解策略的大地电磁正演方法。本发明还提供一种存储介质，其上存储有计算机程序，计算机程序被处理器执行时所述一种基于地空分解策略的大地电磁正演方法的步骤。为了进一步验证所开发系统求解地形问题的精度和收敛性，建立如图2梯形山的地形模型，在坐标轴中心设置一个梯形的地形模型，梯形的高度为450m，其顶部的宽度大小为450m×450m，底部宽度大小为2000m×2000m，整个求解区域的大小为20km×20km×20km，纯地形模型的电导率为0.01s·m-1，空气的电导率为1e-8s·m-1，沿着y轴布设20条观测剖面，每条观测剖面的范围为y∈[-1500m,1500m]，激发频率分别为2Hz，具体情况下如图3的x-y和x-z剖面所示，图4展示了网格剖分示意图。从图5可知，视电阻率ρxx数值非常小，高值区表示山脊所对应的位置；而视电阻率ρxy在山脊的位置表现出低阻现象。视电阻率ρxy在x方向分界面明显，y方向的分界面在该视电阻率等值线平面图较难与实际的地形起伏有一一对应关系。从图6可知，视电阻率ρyy数值非常小，高值区表示山脊所对应的位置；视电阻率ρyx等值线图在y方向分界面明显。
需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。

Claims (10)

1.一种基于地空分解策略的大地电磁正演方法，其特征在于：包括如下步骤：
步骤1：构建求解区域三维MT的地电几何模型；
步骤2：将所述地电几何模型离散为一系列不重叠的四面体单元；
步骤3：获取基于地空分解策略的A-Φ-Ψ耦合势三维MT正演满足的积分弱解形式，并采用节点基函数和四面体单元的节点对所述积分弱解形式进行基于节点离散的转换得到稀疏线性方程组，其中，A表示磁场矢量位、Φ表示电场标量位，Ψ表示磁场标量位，节点为离散四面体单元之间的连接枢纽；
步骤4：加载入射场置于稀疏线性方程组的右端项并求解所述稀疏线性方程得到求解区域中四面体单元各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ，其中，空气区域内节点的磁场矢量位A、电场标量位Φ为0，地下区域内节点的磁场标量位Ψ为0；
步骤5：根据所述磁场矢量位A、电场标量位Φ，磁场标量位Ψ计算出观测点的电场、磁场一个或多个参数的组合。
2.根据权利要求1所述的方法，其特征在于：步骤3中所述积分弱解形式如下：
Figure FDA0003936360290000011
Figure FDA0003936360290000012
Figure FDA0003936360290000013
式中，V表示测试函数，E表示电场，B为磁场,Ω表示求解区域，
Figure FDA0003936360290000014
表示求解区域的边界，Ω0为空气空间区域，Ω1为地下空间区域，Γ0表示空气层的上顶面，Γ1表示地下区域的下底面，Γ为地下区域和空气区域的分界面，n表示边界外法向量，dv、ds分别表示体积微分元和面积微分元，v表示体积、s表示面积，
Figure FDA0003936360290000015
表示阻抗率，
Figure FDA0003936360290000016
表示导纳率，μ为磁导率。
3.根据权利要求1所述的方法，其特征在于：步骤3中所述稀疏线性方程组如下：
Figure FDA0003936360290000017
式中，
Figure FDA0003936360290000018
均为稀疏的子矩阵，
Figure FDA0003936360290000021
均为地空分界面的面积分项形成的子矩阵,
Figure FDA0003936360290000022
均表示零矩阵，B0,B1,B2,B3,B4均为稀疏线性方程组的右端项，X0、X1、X2、X3、X4分别为稀疏线性方程组中的未知数矩阵，且由求解区域内所有节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ构成。
4.根据权利要求3所述的方法，其特征在于：所述稀疏线性方程组中的参数表示如下：
Figure FDA0003936360290000023
Figure FDA0003936360290000024
Figure FDA0003936360290000025
X3＝(Φ1,Φ2,…Φn)T，
Figure FDA0003936360290000026
Figure FDA0003936360290000027
Figure FDA0003936360290000028
Figure FDA0003936360290000029
Figure FDA00039363602900000210
Figure FDA00039363602900000211
Figure FDA00039363602900000212
Figure FDA00039363602900000213
Figure FDA00039363602900000214
Figure FDA00039363602900000215
Figure FDA00039363602900000216
Figure FDA00039363602900000217
Figure FDA00039363602900000218
Figure FDA00039363602900000219
Figure FDA00039363602900000220
Figure FDA00039363602900000221
Figure FDA0003936360290000031
Figure FDA0003936360290000032
Figure FDA0003936360290000033
式中，i,j＝1～n；m,k＝1～n1，
Figure FDA0003936360290000034
Figure FDA0003936360290000035
表示笛卡尔坐标系三个方向，
Figure FDA0003936360290000036
分别表示对x求导，
Figure FDA0003936360290000037
分别表示对y求导，
Figure FDA0003936360290000038
分别表示对z求导，[·]x表示x方向上的分量，
Figure FDA0003936360290000039
表示地下区域第一个节点，第二个节点、第n节点在x方向的磁场矢量位Ax；
Figure FDA00039363602900000310
表示地下区域第一个节点，第二个节点、第n节点在y方向的磁场矢量位Ay；
Figure FDA00039363602900000311
表示地下区域第一个节点，第二个节点、第n节点在z方向的磁场矢量位Az，x、y、z方向的磁场矢量位Ax、Ay、Az构成磁场矢量位A；Φ1、Φ2、Φn表示地下区域中第一个节点，第二个节点、第n节点的电场标量位；Ψ1、Ψ2、Ψn1表示空气区域中第一个节点，第二个节点、第n1节点的磁场标量位，T表示矩阵转置，n为地下区域的节点总数，n1为空气区域的节点总数，节点对应为四面体单元的顶点；
Li,Lj分别表示地下区域内节点i，节点j对应的节点基函数，Lk,Lm分别表示空气区域中第k个、第m个节点对应的节点基函数，
Figure FDA00039363602900000312
表示子矩阵
Figure FDA00039363602900000313
中节点i和节点j对应的元素，
Figure FDA00039363602900000314
表示子矩阵B0中节点i对应的元素，
Figure FDA00039363602900000315
分别表示极化场源在均匀半空间或层状介质地电模型下激发后产生的背景磁场矢量位的三个分量，Eb,Bb分别表示极化场源在均匀半空间或层状介质地电模型下激发后产生的背景电场和磁场，
Figure FDA00039363602900000316
表示阻抗率，
Figure FDA00039363602900000317
表示导纳率，Ω0为空气空间区域，Ω1为地下空间区域，Γ0表示空气层的上顶面，Γ1表示地下区域的下底面，Γ为地下区域和空气区域的分界面，n表示边界外法向量，dv、ds分别表示体积微分元和面积微分元，v表示体积、s表示面积，μ为磁导率。
5.根据权利要求1所述的方法，其特征在于：步骤4中加载入射场置于稀疏线性方程组的右端项并求解所述稀疏线性方程得到求解区域中四面体单元各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ的过程如下：
①：构建基于地空分解策略的A-Φ-Ψ耦合势三维MT对偶问题，并求解对偶问题得到对偶场WA-Wφ-WΨ；以及基于加载的入射场计算出所述稀疏线性方程组的右端项，并求解稀疏线性方程组得到求解区域内当前四面体单元离散对应的节点的磁场矢量位A、电场标量位Φ和磁场标量位Ψ；
其中，WA表示磁场矢量位A的对偶场，WΦ表示为电场标量位Φ的对偶场，Wψ为磁场标量位Ψ的对偶场；
②：基于步骤①计算出的每个节点上的磁场矢量位A、电场标量位Φ和磁场标量位Ψ以及磁场矢量位A的对偶场WA、电场标量位Φ的对偶场WΦ、磁场标量位Ψ的对偶场Wψ计算出大地电磁法问题的对偶问题分别对应的单元误差分布，再基于单元误差分布计算出四面体单元的相对误差单元指示因子；
③：分别判断每个四面体单元的相对误差指示因子是否大于预设最大误差，若大于预设最大误差，对所述四面体单元进行细化离散，再返回步骤3，直至满足迭代终止条件；
所述迭代终止条件为：迭代次数达到预设最大迭代次数或求解区域的四面体单元总数达到预设最大单元数；
其中，最后所得到的所有四面体单元节点上的磁场矢量位A、电场标量位Φ和磁场标量位Ψ。
6.根据权利要求5所述的方法，其特征在于：大地电磁法问题的对偶问题分别对应的单元误差分布如下所示：
Figure FDA0003936360290000041
Figure FDA0003936360290000042
式中，e、w分别表示大地电磁法问题、对偶问题分别对应的单元误差分布，n表示边界外法向量，I为虚数单元，ω为角频率，ds表示面积微分元，s表示面积，μ为磁导率，
Figure FDA0003936360290000043
为四面体单元的表面，符号[·]表示场在面
Figure FDA0003936360290000044
的跳跃值，
Figure FDA0003936360290000045
表示导纳率；
所述相对误差单元指示因子如下：
Figure FDA0003936360290000046
其中，γ＝|ew|
式中，γe为四面体单元的相对误差单元指示因子，γ为四面体单元的加权误差，
Figure FDA0003936360290000047
为所有四面体单元中加权误差的最大值。
7.根据权利要求1所述的方法，其特征在于：步骤4中加载入射场置于稀疏线性方程组的右端项分为两个极化方向分别进行入射场加载得到两个极化方向对应的各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ，所述方法还包括：
基于两个极化方向对应的各个节点的磁场矢量位A、电场标量位Φ，磁场标量位Ψ得到阻抗张量、视电阻率、相位响应中一个或多个的组合。
8.根据权利要求1所述的方法，其特征在于：对所述稀疏线性方程进行求解之前还包括获取所述地电几何模型中不同区域的物性参数；
其中，所述物性参数包括各个区域的电阻率参数、介电常数参数以及磁导率参数。
9.一种基于权利要求1-8任一项所述方法的装置，其特征在于：包括：
模型构建模块：用于构建求解区域三维MT的地电几何模型；
四面体单元离散模块：用于将地电几何模型离散为一系列不重叠的四面体单元；
运算模块：用于基于地空分解策略的A-Φ-Ψ耦合势进行正演计算，最终得到观测点的目标参数。
10.一种存储介质，其上存储有计算机程序，其特征在于：所述计算机程序被处理器执行时实现权利要求1至8任一项所述方法的步骤。
CN202010515462.2A 2020-06-09 2020-06-09 基于地空分解策略的大地电磁正演方法及装置、存储介质 Active CN111638556B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010515462.2A CN111638556B (zh) 2020-06-09 2020-06-09 基于地空分解策略的大地电磁正演方法及装置、存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010515462.2A CN111638556B (zh) 2020-06-09 2020-06-09 基于地空分解策略的大地电磁正演方法及装置、存储介质

Publications (2)

Publication Number Publication Date
CN111638556A CN111638556A (zh) 2020-09-08
CN111638556B true CN111638556B (zh) 2022-12-27

Family

ID=72329957

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010515462.2A Active CN111638556B (zh) 2020-06-09 2020-06-09 基于地空分解策略的大地电磁正演方法及装置、存储介质

Country Status (1)

Country Link
CN (1) CN111638556B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113761762B (zh) * 2021-08-03 2023-10-20 西北核技术研究所 用于电场/温度有限元数值解的后验误差估计方法
CN113792445B (zh) * 2021-11-15 2022-02-08 中南大学 一种基于积分方程法的三维大地电磁数值模拟方法
CN114970289B (zh) * 2022-07-25 2022-10-25 中南大学 三维大地电磁各向异性正演数值模拟方法、设备及介质
CN115906559B (zh) * 2022-10-31 2023-08-15 重庆大学 一种基于混合网格的大地电磁自适应有限元正演方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6901333B2 (en) * 2003-10-27 2005-05-31 Fugro N.V. Method and device for the generation and application of anisotropic elastic parameters
CA2661505C (en) * 2006-08-24 2018-01-02 Exxonmobil Upstream Research Company Electromagnetic data processing system
US20090083006A1 (en) * 2007-09-20 2009-03-26 Randall Mackie Methods and apparatus for three-dimensional inversion of electromagnetic data
KR101403296B1 (ko) * 2013-12-09 2014-06-03 한국지질자원연구원 3차원 항공 자력 탐사 시스템 및 이를 이용한 3차원 항공 자력 탐사 방법
CN104122585B (zh) * 2014-08-08 2017-07-21 中国石油大学(华东) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
CN106980736B (zh) * 2017-04-11 2019-07-19 吉林大学 一种各向异性介质的海洋可控源电磁法有限元正演方法
CN110058315B (zh) * 2019-05-29 2020-04-14 中南大学 一种三维各向异性射频大地电磁自适应有限元正演方法
CN110346835B (zh) * 2019-07-22 2020-11-17 中国科学院地球化学研究所 大地电磁的正演方法、正演系统、存储介质及电子设备
CN110989002B (zh) * 2019-11-07 2021-01-05 吉林大学 地空时域电磁系统浅部低阻异常体数据解释方法

Also Published As

Publication number Publication date
CN111638556A (zh) 2020-09-08

Similar Documents

Publication Publication Date Title
CN111638556B (zh) 基于地空分解策略的大地电磁正演方法及装置、存储介质
Zhdanov et al. Electromagnetic inversion using quasi-linear approximation
Li et al. A finite-element time-domain forward solver for electromagnetic methods with complex-shaped loop sources
Oldenburg et al. Three dimensional inversion of multisource time domain electromagnetic data
Grayver et al. 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO 2 storage formation
Streich 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy
Abubakar et al. 2.5 D forward and inverse modeling for interpreting low-frequency electromagnetic measurements
Newman et al. Solution accelerators for large-scale three-dimensional electromagnetic inverse problems
CN106980736B (zh) 一种各向异性介质的海洋可控源电磁法有限元正演方法
Hou et al. Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials
Liu et al. Adaptive finite element modelling of three-dimensional magnetotelluric fields in general anisotropic media
Zhdanov et al. Large-scale 3D inversion of marine magnetotelluric data: Case study from the Gemini prospect, Gulf of Mexico
Chakravarthi et al. 3D gravity inversion of basement relief—A depth-dependent density approach
US20090083006A1 (en) Methods and apparatus for three-dimensional inversion of electromagnetic data
SHEN Modeling of 3‐D Electromagnetic Responses in Frequency domain by Using the Staggeredgrid Finite Difference Method
Jaysaval et al. Fast multimodel finite-difference controlled-source electromagnetic simulations based on a Schur complement approach
Zhang et al. 3D inversion of time-domain electromagnetic data using finite elements and a triple mesh formulation
Kong et al. Three‐Dimensional Inversion of Magnetotelluric Data for a Resistivity Model With Arbitrary Anisotropy
Huang et al. Spectral-element method with arbitrary hexahedron meshes for time-domain 3D airborne electromagnetic forward modeling
CN111611737B (zh) 一种三维任意各向异性介质的海洋可控源电磁正演方法
Vachiratienchai et al. An efficient inversion for two-dimensional direct current resistivity surveys based on the hybrid finite difference–finite element method
Wang et al. 2D joint inversion of CSAMT and magnetic data based on cross-gradient theory
Yang et al. Survey decomposition: A scalable framework for 3D controlled-source electromagnetic inversion
Akça et al. Extraction of structure-based geoelectric models by hybrid genetic algorithms
Cai et al. Effective 3D-transient electromagnetic inversion using finite-element method with a parallel direct solver

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