CN113553748A - 一种三维大地电磁正演数值模拟方法 - Google Patents

一种三维大地电磁正演数值模拟方法 Download PDF

Info

Publication number
CN113553748A
CN113553748A CN202111102810.4A CN202111102810A CN113553748A CN 113553748 A CN113553748 A CN 113553748A CN 202111102810 A CN202111102810 A CN 202111102810A CN 113553748 A CN113553748 A CN 113553748A
Authority
CN
China
Prior art keywords
electric field
grid
dimensional
equation
field component
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.)
Granted
Application number
CN202111102810.4A
Other languages
English (en)
Other versions
CN113553748B (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.)
Central South University
Original Assignee
Central South University
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 Central South University filed Critical Central South University
Priority to CN202111102810.4A priority Critical patent/CN113553748B/zh
Publication of CN113553748A publication Critical patent/CN113553748A/zh
Application granted granted Critical
Publication of CN113553748B publication Critical patent/CN113553748B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • 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/38—Processing data, e.g. for analysis, for interpretation, for correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Geology (AREA)
  • Geophysics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Environmental & Geological Engineering (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请公开了一种三维大地电磁正演数值模拟方法，首先根据目标地质体构建电阻率分布模型，然后通过多重网格法粗化电阻率分布模型，并通过交错网格有限差分法离散粗细网格上的双旋度方程，得到系数矩阵，接着对离散后的双旋度方程作4色分块处理，通过二维有限差分法计算出目标地质体模型的边界条件，并计算出细网格上双旋度方程的右端项，最后使用基于4色分块高斯‑赛德平滑的多重网格法求解
Figure 162914DEST_PATH_IMAGE001
；更换极化模式并重复上述过程，根据不同极化模式下的电场分量和磁场分量，计算出对应测点的视电阻率和阻抗相位。本申请采用线分块高斯‑赛德平滑技术，有效去除粗细网格上的高频残差分量，提高了多重网格法的计算效率，达到快速收敛的目的。

Description

一种三维大地电磁正演数值模拟方法
技术领域
本发明涉及地球物理技术领域，具体地，涉及一种三维大地电磁正演数值模拟方法。
背景技术
大地电磁法勘探是一种利用大地中广泛分布的天然变化电磁场（10-4~104Hz）进行深部地质构造研究的方法，因具有探测深度大、工作效率高、支出成本低等优点，被广泛用于工程勘探、金属和油气资源勘查以及深部地球动力学研究等领域，而能高效地开展电磁三维正反演是其得以广泛应用的前提。
在三维电磁法正演数值模拟过程中，需要对涉及磁场旋度和电场旋度的双旋度方程进行求解，但该方程具有丰富的零空间，而且当频率很低时，该方程往往无法模拟出电性变化，因此采用常规的迭代方法如BICGstab求解该方程，其收敛效率低，必须进行散度校正以加快收敛速度。但随着计算频率越低、剖分网格数越多，迭代次数表现出线性快速增长，因此，再用常规迭代方法模拟研究超深部和复杂的地质异常体，很难具有高效的生产效率。
现有文献公开了多重网格法在扩散电磁场中的应用，为多重网格法在大地电磁模拟中的应用奠定了基础，考虑到多重网格法的收敛速度不会随网格节点数的增加表现出线性增加，因此相比于常规的迭代方法，多重网格法更适合用于对大型方程组求解。但其计算效率很大程度上依赖于平滑算法，而对于电磁场的双旋度方程，常用平滑算法的计算效率也较低。
综上所述，研究人员需要对现有模拟方法进行针对性优化，解决随着计算规模增大、计算效率呈非线性下降，计算频率越低、求解问题越病态，以及常规迭代法收敛速度慢等问题，进而满足大地电磁法勘探中对复杂地质体的快速反演成像需求。
发明内容
本申请的目的在于提供三维大地电磁正演数值模拟方法，在模拟深部复杂地质体产生电场和磁场响应时确保结果精度和计算效率。本申请的技术方案如下：
一种三维大地电磁正演数值模拟方法，包括如下步骤：
步骤1）根据目标地质体的形状、大小和电阻率分布特征构建对应的电阻率分布模型；
步骤2）通过多重网格法将所述电阻率分布模型粗化，得到不同尺度下粗细网格的电阻率分布模型，并计算细网格和粗网格空间关系的限制算子和插值算子；
步骤3）通过交错网格有限差分法离散粗细网格上电场满足的双旋度方程，基于线分块高斯-赛德平滑器原理对离散后的双旋度方程作4色分块处理；
步骤4）根据细网格的电阻率分布模型，通过二维有限差分法计算出目标地质体模型的边界条件，并计算出细网格上双旋度方程的右端项；
步骤5）根据细网格上离散后的双旋度方程，通过线分块高斯-赛德平滑器进行前光滑，得到对应的细网格上的电场分量，并计算出细网格上的电场残差分量；
步骤6）根据所述的细网格上的电场残差分量和限制算子，计算粗网格上的电场残差分量；
步骤7）根据粗网格上离散后的双旋度方程以及电场残差分量，求解该双旋度方程，得到粗网格上的电场迭代量；
步骤8）根据粗网格上的电场迭代量和插值算子，计算出细网格上的电场迭代量；
步骤9）根据细网格上的电场分量和电场迭代量，计算出新的电场分量，并使用线分块高斯-赛德平滑器进行后光滑，更新细网格上的电场分量；
步骤10）重复步骤5）~步骤9），直到细网格上的相对残差小于收敛阈值；
步骤11）根据细网格上的电场分量和旋度算子，计算出细网格上的磁场分量；
步骤12）更换极化模式，重复步骤4）~步骤11）；
步骤13）根据细网格上不同极化模式下的电场分量和磁场分量，计算出对应测点的视电阻率和阻抗相位，从而完成一个计算频率的三维大地电磁正演数值模拟。
在一些具体的实施例中，所述步骤1）包括以下具体过程：
确定研究区域、计算频率和观测方式，布置对应的测网并构建三维长方体模型；
将所述三维长方体模型剖分成若干个细长方体单元，分别将沿
Figure DEST_PATH_IMAGE001
、
Figure 835060DEST_PATH_IMAGE002
、
Figure DEST_PATH_IMAGE003
方向剖分的细长方体单元个数表示为
Figure 486621DEST_PATH_IMAGE004
、
Figure DEST_PATH_IMAGE005
和
Figure 453309DEST_PATH_IMAGE006
，每个细长方体单元的长宽高分别为
Figure DEST_PATH_IMAGE007
、
Figure 848519DEST_PATH_IMAGE008
和
Figure DEST_PATH_IMAGE009
；
根据目标地质体的形状、大小和电阻率分布特征对每个细长方体单元的电阻率进行赋值，进而实现构建对应目标地质体的电阻率分布模型。
在一些具体的实施例中，所述步骤2）包括以下具体过程：
将
Figure 600574DEST_PATH_IMAGE001
、
Figure 423036DEST_PATH_IMAGE010
、
Figure 879950DEST_PATH_IMAGE003
三个方向的细长方体单元两两组合形成一个粗长方体单元，即八个细长方体单元组合成一个粗长方体单元，则沿
Figure 78850DEST_PATH_IMAGE001
、
Figure 747729DEST_PATH_IMAGE010
、
Figure 678776DEST_PATH_IMAGE003
方向剖分的粗长方体单元个数为
Figure DEST_PATH_IMAGE011
、
Figure 433105DEST_PATH_IMAGE012
和
Figure DEST_PATH_IMAGE013
；
粗长方体单元的长宽高及电阻率由其所包含的细长方体单元的长宽高及电阻率决定，进而构建出粗网格的电阻率分布模型；
根据细长方体单元和粗长方体单元上的空间关系，分别计算出限制算子
Figure 622647DEST_PATH_IMAGE014
和插 值算子
Figure DEST_PATH_IMAGE015
；
基于三维大地电磁正演数值模拟的需求，可进一步粗化模型，得到更粗的电阻率 分布模型，并得到对应的限制算子
Figure 83715DEST_PATH_IMAGE016
、
Figure DEST_PATH_IMAGE017
……
Figure 247980DEST_PATH_IMAGE018
和插值算子
Figure 489606DEST_PATH_IMAGE019
、
Figure 217259DEST_PATH_IMAGE020
……
Figure 595151DEST_PATH_IMAGE021
，
Figure 195897DEST_PATH_IMAGE022
为模型粗化次数。
在一些具体的实施例中，所述步骤3）中离散双旋度方程包括以下具体过程：
由麦克斯方程组可推出三维大地电磁电场强度满足的双旋度方程：
Figure 659239DEST_PATH_IMAGE023
(1)
式中：
Figure 941316DEST_PATH_IMAGE024
为电场强度，
Figure 173714DEST_PATH_IMAGE025
为角频率且与计算频率
Figure 945361DEST_PATH_IMAGE026
的关系为
Figure 895999DEST_PATH_IMAGE027
，
Figure 965456DEST_PATH_IMAGE028
为磁导率，
Figure 317939DEST_PATH_IMAGE029
为介质电导率且与电阻率
Figure 260488DEST_PATH_IMAGE030
的关系为
Figure 698422DEST_PATH_IMAGE031
；
采用交错网格有限差分法，可将公式（1）离散为：
Figure 322302DEST_PATH_IMAGE032
(2)
Figure 529292DEST_PATH_IMAGE033
(3)
式中：
Figure 642742DEST_PATH_IMAGE034
为长方体单元上棱边中点的电场分量，
Figure 567972DEST_PATH_IMAGE035
为由棱边中点到面中心的旋度 算子，
Figure 792280DEST_PATH_IMAGE036
为由面中心到棱边中点的旋度算子，
Figure 309237DEST_PATH_IMAGE037
为棱边相连长方体单体的电导率体积平 均；
将步骤1）中构建的三维长方体模型的六个面上的电场分量表示为
Figure 328008DEST_PATH_IMAGE038
，即为边界电 场分量，而位于三维长方体模型内部的电场分量用
Figure 740535DEST_PATH_IMAGE039
表示，即为内部电场分量，则公式（2） 可表示为：
Figure 768534DEST_PATH_IMAGE040
(4)
式中：
Figure 887800DEST_PATH_IMAGE041
和
Figure 343052DEST_PATH_IMAGE042
分别为
Figure 977295DEST_PATH_IMAGE043
由内部电场分量和边界电场分量分离后的系数矩阵，将
Figure 543406DEST_PATH_IMAGE044
项替换为
Figure 32025DEST_PATH_IMAGE045
，并移动到等式右边作为等式右端项，公式（4）简化为：
Figure 658178DEST_PATH_IMAGE046
(5)
采用上述方法离散粗细网格上电场满足的双旋度方程，得到对应的系数矩阵，，用
Figure 779718DEST_PATH_IMAGE047
、
Figure 149520DEST_PATH_IMAGE048
和
Figure 305695DEST_PATH_IMAGE049
分别表示细网格上的系数矩阵、电场分量和右端项，用
Figure 40432DEST_PATH_IMAGE050
、
Figure 649268DEST_PATH_IMAGE051
……
Figure 557181DEST_PATH_IMAGE052
表 示粗网格上的系数矩阵。
在一些具体的实施例中，所述步骤3）中对离散后的双旋度方程作4色分块处理包括以下具体过程：
基于线分块高斯-赛德平滑器的工作原理，将得到的电场分量分为4类并用不同颜色标识，同种颜色线上所有节点相连的电场分量相互独立并进行同步计算，而不同颜色间通过相邻线上所有节点相连的电场分量进行边界耦合，即
将
Figure 567863DEST_PATH_IMAGE053
分解为4个小的系数矩阵
Figure 722769DEST_PATH_IMAGE054
、
Figure 818901DEST_PATH_IMAGE055
、
Figure 530505DEST_PATH_IMAGE056
和
Figure 395693DEST_PATH_IMAGE057
，将
Figure 472234DEST_PATH_IMAGE058
分解为
Figure 55662DEST_PATH_IMAGE059
、
Figure 570957DEST_PATH_IMAGE060
、
Figure 290651DEST_PATH_IMAGE061
和
Figure 787360DEST_PATH_IMAGE062
，将
Figure 858084DEST_PATH_IMAGE045
分解为
Figure 911491DEST_PATH_IMAGE063
、
Figure 751271DEST_PATH_IMAGE064
、
Figure 231931DEST_PATH_IMAGE065
和
Figure 727634DEST_PATH_IMAGE066
，对于公式（5）的求解转为对以下4个子方程的求解：
Figure 319153DEST_PATH_IMAGE067
(6)。
在一些具体的实施例中，所述步骤4）包括以下具体过程：
选择
Figure 13439DEST_PATH_IMAGE068
或
Figure 930580DEST_PATH_IMAGE069
极化，将三维长方体模型切分为一系列二维长方体模型；
通过二维有限差分法计算出二维长方体模型的电场分量，进而得到三维长方体模 型剖分后的所有电场分量，将位于三维长方体模型的六个面上的电场分量赋给
Figure 900197DEST_PATH_IMAGE038
，并根据
Figure 295407DEST_PATH_IMAGE070
计算出双旋度方程的右端项，将位于内部的电场分量赋给
Figure 109779DEST_PATH_IMAGE071
，作为初始估计值并用
Figure 932241DEST_PATH_IMAGE072
表示。
在一些具体的实施例中，所述步骤5）~9）包括以下具体过程：
将
Figure 136958DEST_PATH_IMAGE072
带入公式（5），并用公式（6）对其进行平滑更新，得到对应的细网格上的电场 分量
Figure 335858DEST_PATH_IMAGE073
，同时计算出细网格上的电场残差分量：
Figure 4737DEST_PATH_IMAGE074
(7)
配合限制算子
Figure 998100DEST_PATH_IMAGE014
，将细网格上的电场残差分量限制为粗网格上的电场残差分 量：
Figure 939380DEST_PATH_IMAGE075
(8)
根据粗网格上离散得到的系数矩阵
Figure 676392DEST_PATH_IMAGE076
，使用直接求解器计算出粗网格上的电场 迭代量：
Figure 199778DEST_PATH_IMAGE077
(9)
配合插值算子
Figure 629622DEST_PATH_IMAGE078
，将粗网格上的电场迭代量延拓为细网格上的电场迭代量：
Figure 871247DEST_PATH_IMAGE079
(10)
根据细网格上的电场分量
Figure 84054DEST_PATH_IMAGE080
和电场迭代量
Figure 727525DEST_PATH_IMAGE081
对电场分量进行粗网格校正，进 而得到新的电场分量；
将校正后的新电场分量带入公式（5），并用公式（6）对其进行平滑更新，得到对应 的细网格上新的电场分量
Figure 62691DEST_PATH_IMAGE080
。
在一些具体的实施例中，在所述步骤10）中，当细网格上的电场残差分量满足相对 残差
Figure 791613DEST_PATH_IMAGE082
＜1×10-10时，停止迭代计算。
在一些具体的实施例中，所述步骤13）包括以下具体过程：
根据两种极化模式下的电场分量和磁场分量，计算出
Figure 57378DEST_PATH_IMAGE083
、
Figure 555355DEST_PATH_IMAGE084
两种模式的阻抗：
Figure 327002DEST_PATH_IMAGE085
(11)
式中：
Figure 277641DEST_PATH_IMAGE086
、
Figure 97829DEST_PATH_IMAGE087
、
Figure 450313DEST_PATH_IMAGE088
和
Figure 392861DEST_PATH_IMAGE089
为
Figure 830796DEST_PATH_IMAGE068
极化计算出
Figure 438364DEST_PATH_IMAGE090
、
Figure 910933DEST_PATH_IMAGE002
方向的电场分量和磁场分量，
Figure 24383DEST_PATH_IMAGE091
、
Figure 949614DEST_PATH_IMAGE092
、
Figure 173922DEST_PATH_IMAGE093
和
Figure 173102DEST_PATH_IMAGE094
为
Figure 457452DEST_PATH_IMAGE069
极化计算出
Figure 869979DEST_PATH_IMAGE090
、
Figure 897978DEST_PATH_IMAGE010
方向的电场分量和磁场分量；
根据所得阻抗计算出
Figure 281160DEST_PATH_IMAGE083
、
Figure 470833DEST_PATH_IMAGE084
两种模式下的视电阻率和阻抗相位：
Figure 370656DEST_PATH_IMAGE095
Figure 936766DEST_PATH_IMAGE096
(12)。
本申请提供的技术方案至少具有如下有益效果：
1、本申请采用了具有高度并行性的线分块高斯-赛德平滑技术，使得本方法中包含了大地电磁数字模拟中电场分量散度为零的信息，因此无需进行散度校正，加快了多重网格的计算效率；此外还可以有效去除粗细网格上的高频残差分量，采用线分块高斯-赛德平滑器对双旋度方程的系数矩阵进行重新排序分解，分解出的4个小的系数矩阵稀疏性变好，使得本方法在低频时也具有较快的收敛速度和计算效率。
2、本申请对于三维地质体模型使用长方体单元进行剖分，对于每个长方体单元赋值的电阻率可以不同，因此可模拟任意电阻率分布的目标地质体。
3、在本申请中：使用交错网格有限差分法离散长方体单元上电场满足的双旋度方程，有利于简单高效地离散粗细网格上的双旋度方程，进而得到对应的系数矩阵；使用二维有限差分数值解作为边界条件，使得数值模拟结果更准确；使用多重网格法求解交错网格有限差分法离散得到双旋度方程，这样可以将细网格上的低频电场残差分量限制为粗网格上的高频残差分量进行剔除，从而达到快速收敛的目的。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本申请提供的三维大地电磁正演数值模拟方法的流程图；
图2为本申请方法中线分块高斯-赛德平滑器的分块原理图，其中，图2（a）为线平滑系统，图2（b）为线分块示意图，1、2、3、4分别代表4种不同的颜色标识；
图3为本申请实施例1中低-高阻组合异常体模型示意图；
图4为本申请实施例1中采用快速多重网格法与国际公开程序ModEM计算出的低-高阻组合异常体在地表测线上的结果，其中，图4（a）为
Figure 176118DEST_PATH_IMAGE083
模式下二者所得视电阻率的对比曲线图，图4（b）为
Figure 536692DEST_PATH_IMAGE083
模式下二者所得阻抗相位的对比曲线图，图4（c）为
Figure 923811DEST_PATH_IMAGE084
模式下二者所得视电阻率的对比曲线图，图4（d）为
Figure 293612DEST_PATH_IMAGE084
模式下二者所得阻抗相位的对比曲线图；
图5为本申请实施例1中采用快速多重网格法与传统迭代法BICGstab在有散度校正和无散度校正下的收敛曲线对比图。
具体实施方式
为了便于理解本申请，下面将结合说明书附图和较佳的实施例对本申请中的技术方案作更全面、细致地描述，但本申请的保护范围并不限于以下具体的实施例，基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，均属于本申请保护的范围。
需要特别说明的是，当某一元件被描述为与另一元件存在“固定、固接、连接或连通”关系时，它可以是直接固定、固接、连接或连通在另一元件上，也可以是通过其他中间件间接固定、固接、连接或连通在另一元件上。
除非另有定义，下文中所使用的所有专业术语与本领域技术人员通常理解的含义相同。本文中所使用的专业术语只是为了描述具体实施例的目的，并不是旨在限制本申请的保护范围。
参见图1，一种三维大地电磁正演数值模拟方法，具体过程如下：
第一步、构建三维电阻率分布模型
根据目标地质体的形状、大小和电阻率分布确定研究区域、计算频率、观测方式，布置对应的测网并构建三维长方体模型。
将该三维长方体模型进行网格剖分，剖分成若干个细长方体单元，分别将沿
Figure 371159DEST_PATH_IMAGE090
、
Figure 168213DEST_PATH_IMAGE002
、
Figure 777049DEST_PATH_IMAGE003
方向剖分的细长方体单元个数表示为
Figure 950541DEST_PATH_IMAGE004
、
Figure 898906DEST_PATH_IMAGE005
和
Figure DEST_PATH_IMAGE097
，每个细长方体单元的长宽高分别为
Figure 866862DEST_PATH_IMAGE007
、
Figure 962994DEST_PATH_IMAGE008
和
Figure 674598DEST_PATH_IMAGE009
。
根据目标地质体的形状、大小和电阻率分布特征对每个细长方体单元的电阻率进行赋值，每个细长方体单元的电阻率可不同，因此可模拟任意情况下的目标地质体构建对应的电阻率分布模型，所述空气层中的每个细长方体单元电阻率赋为
Figure 726736DEST_PATH_IMAGE098
。
第二步、粗化三维电阻率分布模型
将
Figure 865594DEST_PATH_IMAGE001
、
Figure 449022DEST_PATH_IMAGE010
、
Figure 698738DEST_PATH_IMAGE003
三个方向的细长方体单元两两组合形成一个粗长方体单元，即八个细长方体单元组合成一个粗长方体单元，则沿
Figure 621694DEST_PATH_IMAGE001
、
Figure 931453DEST_PATH_IMAGE010
、
Figure 2177DEST_PATH_IMAGE003
方向剖分的粗长方体单元个数为
Figure DEST_PATH_IMAGE099
、
Figure 242534DEST_PATH_IMAGE012
和
Figure 82314DEST_PATH_IMAGE100
。
粗长方体单元的长宽高由其所包含的细长方体单元的长宽高对应相加得到，粗长方体单元的电阻率由其所包含的细长方体单元的电阻率体积平均获得，这样就构建出粗网格的电阻率分布模型。
根据细长方体单元和粗长方体单元上的空间关系，分别计算出限制算子
Figure 235078DEST_PATH_IMAGE014
和插 值算子
Figure 527519DEST_PATH_IMAGE015
。
基于三维大地电磁正演数值模拟的需求，可进一步网格粗化，得到更粗的电阻率 分布模型，并得到对应的限制算子
Figure 384617DEST_PATH_IMAGE016
、
Figure 78903DEST_PATH_IMAGE017
……
Figure 123607DEST_PATH_IMAGE018
和插值算子
Figure 168924DEST_PATH_IMAGE019
、
Figure 564133DEST_PATH_IMAGE020
……
Figure 112926DEST_PATH_IMAGE021
，
Figure 200968DEST_PATH_IMAGE022
为模型粗化次数。
第三步、离散粗细网格上的双旋度方程
由麦克斯方程组可推出三维大地电磁电场强度满足的双旋度方程：
Figure 654952DEST_PATH_IMAGE023
(1)
式中：
Figure 853852DEST_PATH_IMAGE024
为电场强度（V/m），
Figure 257151DEST_PATH_IMAGE025
为角频率（Hz）且与计算频率
Figure 516094DEST_PATH_IMAGE026
的关系为
Figure 208107DEST_PATH_IMAGE027
，
Figure 945119DEST_PATH_IMAGE028
为磁导率（H/m），在不考虑地下岩石磁化的情况下，取
Figure DEST_PATH_IMAGE101
，
Figure 468504DEST_PATH_IMAGE029
为介质电导率（S/m）且 与电阻率
Figure 85299DEST_PATH_IMAGE030
（
Figure 61345DEST_PATH_IMAGE102
）的关系为
Figure 602048DEST_PATH_IMAGE031
。
采用交错网格有限差分法，可将公式（1）离散为：
Figure 245519DEST_PATH_IMAGE032
(2)
Figure 315106DEST_PATH_IMAGE033
(3)
式中：
Figure 44028DEST_PATH_IMAGE034
为长方体单元上棱边中点的电场分量，
Figure 122842DEST_PATH_IMAGE035
为由棱边中点到面中心的旋度 算子，
Figure 11033DEST_PATH_IMAGE036
为由面中心到棱边中点的旋度算子，
Figure 517100DEST_PATH_IMAGE037
为棱边相连长方体单体的电导率体积平 均。
将构建的三维长方体模型的六个面上的电场分量表示为
Figure 733318DEST_PATH_IMAGE038
，即为边界电场分量， 而位于三维长方体模型内部的电场分量用
Figure 615823DEST_PATH_IMAGE039
表示，即为内部电场分量，则公式（2）可表示 为：
Figure 905990DEST_PATH_IMAGE040
(4)
式中：
Figure 848539DEST_PATH_IMAGE041
和
Figure 286473DEST_PATH_IMAGE042
分别为
Figure 707090DEST_PATH_IMAGE043
由内部电场分量和边界电场分量分离后的系数矩阵，将
Figure 914081DEST_PATH_IMAGE044
项替换为
Figure 761951DEST_PATH_IMAGE045
，并移动到等式右边作为等式右端项，公式（4）简化为：
Figure 687182DEST_PATH_IMAGE046
(5)
采用上述方法离散粗细网格上电场满足的双旋度方程，得到对应的系数矩阵，用
Figure 911489DEST_PATH_IMAGE047
、
Figure 238566DEST_PATH_IMAGE048
和
Figure 916059DEST_PATH_IMAGE049
分别表示细网格上的系数矩阵、电场分量和右端项，用
Figure DEST_PATH_IMAGE103
、
Figure 328586DEST_PATH_IMAGE104
……
Figure DEST_PATH_IMAGE105
表示粗网格上的系数矩阵。
第四步、线分块高斯-赛德平滑器工作原理推演
由麦克斯韦方程组可推出大地电磁中电场强度
Figure 825426DEST_PATH_IMAGE106
和磁场强度
Figure DEST_PATH_IMAGE107
的关系：
Figure 7009DEST_PATH_IMAGE108
(13)
由公式（13）可知，电场强度可由磁场强度的旋度获得，磁场强度可由电场强度的旋度获得。
离散后的电场分量，沿着一个长方体单元的面作回路积分可得到元面中心点的磁场分量；磁场分量沿所述交错网格有限差分法中双倍体积单元作回路积分可得到长方体单元棱边中点的电场分量。
因此，与
Figure 852474DEST_PATH_IMAGE001
方向一条线上所有节点相连的电场分量只与相邻线上节点的电场分量有关，如图2a所示，系统外围的箭头标识表示相邻一条线上的电场分量，系统内部与中心相连的箭头标识表示当前待求的一条线上所有节点相连的电场分量。
对于不相邻的两条线上所有节点相连的电场分量相互独立，因此，可将第三步中的电场分量归为4类，如图2b中的1、2、3、4标识，相同标识表示为同一类且用一种颜色示出。对同种颜色线上所有节点相连的电场分量相互独立，可进行同步计算，而对不同颜色间通过相邻线上的所有节点相连的电场分量进行边界耦合。
那么，可将
Figure 752297DEST_PATH_IMAGE053
分解为4个小的系数矩阵
Figure 318408DEST_PATH_IMAGE054
、
Figure 541447DEST_PATH_IMAGE055
、
Figure 167601DEST_PATH_IMAGE056
和
Figure 554720DEST_PATH_IMAGE057
，将
Figure 393363DEST_PATH_IMAGE058
分解为
Figure 549538DEST_PATH_IMAGE059
、
Figure 346592DEST_PATH_IMAGE060
、
Figure 345641DEST_PATH_IMAGE061
和
Figure 519134DEST_PATH_IMAGE062
，将
Figure 529815DEST_PATH_IMAGE045
分解为
Figure 497771DEST_PATH_IMAGE063
、
Figure 328324DEST_PATH_IMAGE064
、
Figure 774348DEST_PATH_IMAGE065
和
Figure 905116DEST_PATH_IMAGE066
，1号色对应
Figure DEST_PATH_IMAGE109
，2号色对应
Figure 981656DEST_PATH_IMAGE110
，3号色对应
Figure 565084DEST_PATH_IMAGE111
，4号色 对应
Figure 207942DEST_PATH_IMAGE112
。对于公式（5）的求解转为对以下4个子方程的求解：
Figure 193216DEST_PATH_IMAGE067
(6)
对所述的不同颜色间通过相邻线上的所有节点相连的电场分量进行边界耦合，就 是将公式（6）中的等式按顺序进行计算：先从初值
Figure 502975DEST_PATH_IMAGE072
分解出
Figure 573699DEST_PATH_IMAGE063
，根据公式（6）中的第一个 式子计算出
Figure 361526DEST_PATH_IMAGE059
，更新
Figure 201306DEST_PATH_IMAGE072
；然后从更新后的
Figure 681966DEST_PATH_IMAGE072
中分解出
Figure 974407DEST_PATH_IMAGE064
，根据公式（6）中的第二个式子 计算出
Figure 831505DEST_PATH_IMAGE060
，更新
Figure 916004DEST_PATH_IMAGE072
；接着从更新后的
Figure 567566DEST_PATH_IMAGE072
中分解出
Figure 612882DEST_PATH_IMAGE065
，根据公式（6）中的第三个式子计算 出
Figure 8091DEST_PATH_IMAGE061
，更新
Figure 291305DEST_PATH_IMAGE072
；最后从更新后的
Figure 379347DEST_PATH_IMAGE072
中分解出
Figure 646380DEST_PATH_IMAGE066
，根据公式（6）中的最后一个式子计算出
Figure 579701DEST_PATH_IMAGE062
，更新
Figure 435530DEST_PATH_IMAGE072
，即完成一次线分块高斯-赛德平滑。
为了防止计算过程中对系数矩阵的LU分解的重复计算，事先将
Figure 694473DEST_PATH_IMAGE113
、
Figure 183224DEST_PATH_IMAGE114
、
Figure 185815DEST_PATH_IMAGE115
和
Figure 949678DEST_PATH_IMAGE116
的LU分解结果保存。
第五步、计算离散后的双旋度方程右端项
根据极化方式，将所述的三维长方体模型切为一系列二维长方形模型，其中，
Figure 51627DEST_PATH_IMAGE068
极化是沿
Figure 683465DEST_PATH_IMAGE001
方向逐层切出
Figure 958589DEST_PATH_IMAGE117
平面的二维长方形模型，
Figure 461114DEST_PATH_IMAGE069
极化是沿
Figure 61860DEST_PATH_IMAGE002
方向逐层切出
Figure 712153DEST_PATH_IMAGE118
平面的二维长方形模型。
通过二维有限差分法计算出二维长方体模型的电场分量，进而得到三维长方体模 型剖分后的所有电场分量，将位于三维长方体模型的六个面上的电场分量赋给
Figure 56546DEST_PATH_IMAGE038
，并根据
Figure 682087DEST_PATH_IMAGE070
计算出双旋度方程的右端项，将位于内部的电场分量赋给
Figure 188155DEST_PATH_IMAGE071
，作为初始估计值并用
Figure 404373DEST_PATH_IMAGE072
表示。
第六步、求解公式（5）
多重网格法就是将细网格上的低频残差分量限制为粗网格上的高频残差分量进行剔除，线分块高斯-赛德平滑器可实现对网格上的高频残差分量快速有效地剔除。因此采用基于线分块高斯-赛德平滑器的快速多重网格法求解公式（5）。
使用线分块高斯-赛德平滑器进行前光滑，即将第五步中的
Figure 224561DEST_PATH_IMAGE072
带入公式（5），并用 公式（6）对其进行平滑更新，得到对应的细网格上的电场分量
Figure 701679DEST_PATH_IMAGE073
，同时计算出细网格上的 电场残差分量：
Figure 378648DEST_PATH_IMAGE074
(7)
配合限制算子
Figure 816583DEST_PATH_IMAGE014
，将细网格上的电场残差分量限制为粗网格上的电场残差分 量：
Figure 689729DEST_PATH_IMAGE075
(8)
根据粗网格上离散得到的系数矩阵
Figure 21354DEST_PATH_IMAGE076
，使用直接求解器计算出粗网格上的电场 迭代量：
Figure 72486DEST_PATH_IMAGE077
(9)
配合插值算子
Figure 187597DEST_PATH_IMAGE078
，将粗网格上的电场迭代量延拓为细网格上的电场迭代量：
Figure 411905DEST_PATH_IMAGE079
(10)
根据细网格上的电场分量
Figure 676665DEST_PATH_IMAGE080
和电场迭代量
Figure 961015DEST_PATH_IMAGE081
对电场分量进行粗网格校正
Figure 107963DEST_PATH_IMAGE119
，得到新的电场分量；
使用线分块高斯-赛德平滑器进行后光滑，即将校正后的新电场分量带入公式 （5），并用公式（6）对其进行平滑更新，得到对应的细网格上新的电场分量
Figure 322912DEST_PATH_IMAGE080
。
第七步、迭代循环
重复第六步的计算过程，直到细网格上的电场残差分量满足相对残差
Figure 504495DEST_PATH_IMAGE082
＜1× 10-10，停止迭代计算，根据公式（13）可知，将计算出的电场分量
Figure 959747DEST_PATH_IMAGE120
取旋度，即可获得地质体 模型中的磁场分量
Figure 859570DEST_PATH_IMAGE121
。
为了防止直接求解器对系数矩阵的LU分解的重复计算，事先将系数矩阵的LU分解结果保存到计算机中。
第八步、更换极化模式，重复第五步、第六步和第七步。
第九步、计算地面对应测线上测点的视电阻率和阻抗相位
根据两种极化模式下的电场分量和磁场分量，计算出
Figure 425681DEST_PATH_IMAGE083
、
Figure 399453DEST_PATH_IMAGE084
两种模式的阻抗：
Figure 25606DEST_PATH_IMAGE122
(11)
式中：
Figure 147146DEST_PATH_IMAGE123
、
Figure 516948DEST_PATH_IMAGE124
、
Figure 860073DEST_PATH_IMAGE125
和
Figure 657128DEST_PATH_IMAGE126
为
Figure 265964DEST_PATH_IMAGE068
极化计算出
Figure 173877DEST_PATH_IMAGE001
、
Figure 309192DEST_PATH_IMAGE002
方向的电场分量和磁场分量，
Figure 277148DEST_PATH_IMAGE127
、
Figure 373280DEST_PATH_IMAGE128
、
Figure 84884DEST_PATH_IMAGE129
和
Figure 153334DEST_PATH_IMAGE130
为
Figure 292191DEST_PATH_IMAGE069
极化计算出
Figure 610040DEST_PATH_IMAGE001
、
Figure 125335DEST_PATH_IMAGE010
方向的电场分量和磁场分量。
根据所得阻抗计算出
Figure 300489DEST_PATH_IMAGE083
、
Figure 610248DEST_PATH_IMAGE084
两种模式下的视电阻率和阻抗相位：
Figure DEST_PATH_IMAGE131
Figure 415393DEST_PATH_IMAGE132
(12)
下面对本申请提供的三维大地电磁正演数值模拟方法进行检验。
实施例1
为了对本发明提出的方法进行检验，设计了如图3所示的低-高阻组合异常体模型，详情如下：在地下有低阻异常体和高阻异常体沿
Figure 406482DEST_PATH_IMAGE090
方向横向排列的组合异常体，其大小均为10km×10km×10km，其埋深均为10km，其间距为10km；三维低-高阻组合异常体模型的空气层电阻率为
Figure 246262DEST_PATH_IMAGE133
，大地层中背景电阻率为100
Figure 726922DEST_PATH_IMAGE134
，低阻异常体的电阻率为10
Figure 284943DEST_PATH_IMAGE102
，高阻异常体的电阻率为1000
Figure 63412DEST_PATH_IMAGE102
。取两个异常体间的中点在地面的投影点作为坐标原点O，测线为
Figure DEST_PATH_IMAGE135
，
Figure 23277DEST_PATH_IMAGE090
方向-25.5km到25.5km范围内均匀间隔设置25个测点。
若采用32×32×32个大小为2.5km×2.5km×2.5km的长方体单元对大地层进行剖分，则模拟区域为：
Figure 612522DEST_PATH_IMAGE090
和
Figure 392259DEST_PATH_IMAGE010
方向均为-40km到40km，
Figure 787468DEST_PATH_IMAGE003
方向为0km到80km；若采用64×64×64个大小为10km×10km×10km的长方体单元对大地层进行剖分，则模拟区域为：
Figure 788791DEST_PATH_IMAGE090
和
Figure 611254DEST_PATH_IMAGE002
方向均为-32km到32km，
Figure 878287DEST_PATH_IMAGE003
方向为0km到64km；若采用128×128×128个大小为10km×10km×10km长方体单元对大地层进行剖分，则模拟区域为：
Figure 77187DEST_PATH_IMAGE090
和
Figure 683749DEST_PATH_IMAGE010
方向均为-64km到64km，
Figure 677113DEST_PATH_IMAGE003
方向为0km到128km。上述三种网格剖分均可有效模拟图3所述的低-高阻组合异常体模型，对于空气层则采用类似的长方体单元剖分。
图4为采用本申请的三维大地电磁正演数值模拟方法与采用国际公开的ModEM程序计算出的低-高阻组合异常体在地表测线上的
Figure 431442DEST_PATH_IMAGE083
、
Figure 620984DEST_PATH_IMAGE084
两种模式下的视电阻率和阻抗相位对比曲线图，其中剖分网格大小为64×64×64，计算频率为0.1Hz。由图4可知本申请方法与国际公开ModEM程序在数值模拟结果上基本吻合，验证了本申请方法在结果精度上的可靠信。
图5为本申请方法与传统迭代法BICGstab在有散度校正和无散度校正下的收敛曲线对比图，其中剖分网格大小为64×64×64，计算频率为0.1Hz，计算极化模式为
Figure 144369DEST_PATH_IMAGE068
，由图5可知本申请方法在几次多重网格迭代内实现快速收敛。
表1
Figure 308634DEST_PATH_IMAGE136
表1为本实施例中基于线分块高斯-赛德平滑器的快速多重网格法在不同计算频率下的迭代次数和计算时间的统计数据，其中剖分网格大小为64×64×64，计算极化模式为
Figure 487943DEST_PATH_IMAGE068
由表1可知本申请方法的计算效率不会随着计算频率变低而变慢，对低频时变态的双旋度方程也具有高效的计算速度，在0.1~0.001Hz范围内，本申请方法的计算速度约为加散度校正的传统迭代法的计算速度的6倍。
表2
Figure DEST_PATH_IMAGE137
表2为本实施例中基于线分块高斯-赛德平滑器的快速多重网格法在不同剖分网格大小下的迭代次数和计算时间统计表，其中计算频率为0.1Hz，计算极化模式为
Figure 28646DEST_PATH_IMAGE068
。由表2可知本申请方法的迭代次数不会随剖分网格大小变大而变多，在模拟剖分网格较多的问题时，计算优势更加明显。
因此，本申请方法对于开展复杂异常体的大规模三维大地电磁正演数值模拟研究具有重要意义，可以大大的提高计算效率，减少时间带来的成本开销。
以上所述仅为本申请的部分实施例，并非因此限制本申请的专利保护范围，对于本领域的技术人员来说，本申请可以有各种更改和变化。在本申请的精神和原则之内，凡是利用本申请说明书及附图内容所作的任何改进或等同替换，直接或间接运用在其它相关的技术领域，均应包括在本申请的专利保护范围内。

Claims (9)

1.一种三维大地电磁正演数值模拟方法，其特征在于，包括如下步骤：
步骤1）根据目标地质体的形状、大小和电阻率分布特征构建对应的电阻率分布模型；
步骤2）通过多重网格法将所述电阻率分布模型粗化，得到不同尺度下粗细网格的电阻率分布模型，并计算细网格和粗网格空间关系的限制算子和插值算子；
步骤3）通过交错网格有限差分法离散粗细网格上电场满足的双旋度方程，基于线分块高斯-赛德平滑器原理对离散后的双旋度方程作4色分块处理；
步骤4）根据细网格的电阻率分布模型，通过二维有限差分法计算出目标地质体模型的边界条件，并计算出细网格上双旋度方程的右端项；
步骤5）根据细网格上离散后的双旋度方程，通过线分块高斯-赛德平滑器进行前光滑，得到对应的细网格上的电场分量，并计算出细网格上的电场残差分量；
步骤6）根据所述的细网格上的电场残差分量和限制算子，计算粗网格上的电场残差分量；
步骤7）根据粗网格上离散后的双旋度方程以及电场残差分量，求解该双旋度方程，得到粗网格上的电场迭代量；
步骤8）根据粗网格上的电场迭代量和插值算子，计算出细网格上的电场迭代量；
步骤9）根据细网格上的电场分量和电场迭代量，计算出新的电场分量，并使用线分块高斯-赛德平滑器进行后光滑，更新细网格上的电场分量；
步骤10）重复步骤5）~步骤9），直到细网格上的相对残差小于收敛阈值；
步骤11）根据细网格上的电场分量和旋度算子，计算出细网格上的磁场分量；
步骤12）更换极化模式，重复步骤4）~步骤11）；
步骤13）根据细网格上不同极化模式下的电场分量和磁场分量，计算出对应测点的视电阻率和阻抗相位，从而完成一个计算频率的三维大地电磁正演数值模拟。
2.根据权利要求1所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤1）包括以下具体过程：
确定研究区域、计算频率和观测方式，布置对应的测网并构建三维长方体模型；
将所述三维长方体模型剖分成若干个细长方体单元，分别将沿
Figure 75085DEST_PATH_IMAGE001
、
Figure 179176DEST_PATH_IMAGE002
、
Figure 958913DEST_PATH_IMAGE003
方向剖分的细长方体单元个数表示为
Figure 354122DEST_PATH_IMAGE004
、
Figure 168495DEST_PATH_IMAGE005
和
Figure 928640DEST_PATH_IMAGE006
，每个细长方体单元的长宽高分别为
Figure 195674DEST_PATH_IMAGE007
、
Figure 394574DEST_PATH_IMAGE008
和
Figure 63452DEST_PATH_IMAGE009
；
根据目标地质体的形状、大小和电阻率分布特征对每个细长方体单元的电阻率进行赋值，进而实现构建对应目标地质体的电阻率分布模型。
3.根据权利要求2所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤2）包括以下具体过程：
将
Figure 243767DEST_PATH_IMAGE001
、
Figure 998096DEST_PATH_IMAGE010
、
Figure 735108DEST_PATH_IMAGE003
三个方向的细长方体单元两两组合形成一个粗长方体单元，即八个细长方体单元组合成一个粗长方体单元，则沿
Figure 258493DEST_PATH_IMAGE001
、
Figure 688338DEST_PATH_IMAGE010
、
Figure 867646DEST_PATH_IMAGE003
方向剖分的粗长方体单元个数为
Figure 142770DEST_PATH_IMAGE011
、
Figure 786241DEST_PATH_IMAGE012
和
Figure 121407DEST_PATH_IMAGE013
；
粗长方体单元的长宽高及电阻率由其所包含的细长方体单元的长宽高及电阻率决定，进而构建出粗网格的电阻率分布模型；
根据细长方体单元和粗长方体单元上的空间关系，分别计算出限制算子
Figure 37279DEST_PATH_IMAGE014
和插值 算子
Figure 116094DEST_PATH_IMAGE015
；
基于三维大地电磁正演数值模拟的需求，可进一步粗化模型，得到更粗的电阻率分布 模型，并得到对应的限制算子
Figure 614071DEST_PATH_IMAGE016
、
Figure 385718DEST_PATH_IMAGE017
……
Figure 274040DEST_PATH_IMAGE018
和插值算子
Figure 156545DEST_PATH_IMAGE019
、
Figure 509029DEST_PATH_IMAGE020
……
Figure 451577DEST_PATH_IMAGE021
，
Figure 76463DEST_PATH_IMAGE022
为模型粗化次数。
4.根据权利要求3所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤3）中离散双旋度方程包括以下具体过程：
由麦克斯方程组可推出三维大地电磁电场强度满足的双旋度方程：
Figure 497080DEST_PATH_IMAGE023
(1)
式中：
Figure 969649DEST_PATH_IMAGE024
为电场强度，
Figure 83099DEST_PATH_IMAGE025
为角频率且与计算频率
Figure 8329DEST_PATH_IMAGE026
的关系为
Figure 170321DEST_PATH_IMAGE027
，
Figure 231817DEST_PATH_IMAGE028
为磁导 率，
Figure 516168DEST_PATH_IMAGE029
为介质电导率且与电阻率
Figure 928695DEST_PATH_IMAGE030
的关系为
Figure 143645DEST_PATH_IMAGE031
；
采用交错网格有限差分法，可将公式（1）离散为：
Figure 325227DEST_PATH_IMAGE032
(2)
Figure 514900DEST_PATH_IMAGE033
(3)
式中：
Figure 414723DEST_PATH_IMAGE034
为长方体单元上棱边中点的电场分量，
Figure 980834DEST_PATH_IMAGE035
为由棱边中点到面中心的旋度算 子，
Figure 220185DEST_PATH_IMAGE036
为由面中心到棱边中点的旋度算子，
Figure 580759DEST_PATH_IMAGE037
为棱边相连长方体单体的电导率体积平 均；
将步骤1）中构建的三维长方体模型的六个面上的电场分量表示为
Figure 967878DEST_PATH_IMAGE038
，即为边界电场 分量，而位于三维长方体模型内部的电场分量用
Figure 337680DEST_PATH_IMAGE039
表示，即为内部电场分量，则公式（2） 可表示为：
Figure 415226DEST_PATH_IMAGE040
(4)
式中：
Figure 212281DEST_PATH_IMAGE041
和
Figure 821117DEST_PATH_IMAGE042
分别为
Figure 994609DEST_PATH_IMAGE043
由内部电场分量和边界电场分量分离后的系数矩阵，将
Figure 942973DEST_PATH_IMAGE044
项替换为
Figure 910929DEST_PATH_IMAGE045
，并移动到等式右边作为等式右端项，公式（4）简化为：
Figure 7061DEST_PATH_IMAGE046
(5)
采用上述方法离散粗细网格上电场满足的双旋度方程，得到对应的系数矩阵，用
Figure 718665DEST_PATH_IMAGE047
、
Figure 770804DEST_PATH_IMAGE048
和
Figure 909661DEST_PATH_IMAGE049
分别表示细网格上的系数矩阵、电场分量和右端项，用
Figure 493089DEST_PATH_IMAGE050
、
Figure 742805DEST_PATH_IMAGE051
……
Figure 728079DEST_PATH_IMAGE052
表示粗网格上的系数矩阵。
5.根据权利要求4所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤3）中对离散后的双旋度方程作4色分块处理包括以下具体过程：
基于线分块高斯-赛德平滑器的工作原理，将得到的电场分量分为4类并用不同颜色标识，同种颜色线上所有节点相连的电场分量相互独立并进行同步计算，而不同颜色间通过相邻线上所有节点相连的电场分量进行边界耦合，即
将
Figure 975520DEST_PATH_IMAGE053
分解为4个小的系数矩阵
Figure 46244DEST_PATH_IMAGE054
、
Figure 99651DEST_PATH_IMAGE055
、
Figure 939431DEST_PATH_IMAGE056
和
Figure 607042DEST_PATH_IMAGE057
，将
Figure 899483DEST_PATH_IMAGE058
分解为
Figure 756580DEST_PATH_IMAGE059
、
Figure 450867DEST_PATH_IMAGE060
、
Figure 40111DEST_PATH_IMAGE061
和
Figure 85428DEST_PATH_IMAGE062
，将
Figure 480637DEST_PATH_IMAGE045
分解为
Figure 29430DEST_PATH_IMAGE063
、
Figure 117472DEST_PATH_IMAGE064
、
Figure 571456DEST_PATH_IMAGE065
和
Figure 770356DEST_PATH_IMAGE066
，对于公式（5）的求解转为对 以下4个子方程的求解：
Figure 173655DEST_PATH_IMAGE067
(6)。
6.根据权利要求5所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤4）包括以下具体过程：
选择
Figure 432598DEST_PATH_IMAGE068
或
Figure 124611DEST_PATH_IMAGE069
极化，将三维长方体模型切分为一系列二维长方体模型；
通过二维有限差分法计算出二维长方体模型的电场分量，进而得到三维长方体模型剖 分后的所有电场分量，将位于三维长方体模型的六个面上的电场分量赋给
Figure 861623DEST_PATH_IMAGE038
，并根据
Figure 385008DEST_PATH_IMAGE070
计算出双旋度方程的右端项，将位于内部的电场分量赋给
Figure 814852DEST_PATH_IMAGE071
，作为初始估计 值并用
Figure 977849DEST_PATH_IMAGE072
表示。
7.根据权利要求6所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤5）~9）包括以下具体过程：
将
Figure 518552DEST_PATH_IMAGE073
带入公式（5），并用公式（6）对其进行平滑更新，得到对应的细网格上的电场分 量
Figure 162023DEST_PATH_IMAGE074
，同时计算出细网格上的电场残差分量：
Figure 497189DEST_PATH_IMAGE075
(7)
配合限制算子
Figure 226111DEST_PATH_IMAGE014
，将细网格上的电场残差分量限制为粗网格上的电场残差分量：
Figure 242608DEST_PATH_IMAGE076
(8)
根据粗网格上离散得到的系数矩阵
Figure 740586DEST_PATH_IMAGE077
，使用直接求解器计算出粗网格上的电场迭代 量：
Figure 246653DEST_PATH_IMAGE078
(9)
配合插值算子
Figure 462871DEST_PATH_IMAGE079
，将粗网格上的电场迭代量延拓为细网格上的电场迭代量：
Figure 532327DEST_PATH_IMAGE080
(10)
根据细网格上的电场分量
Figure 884811DEST_PATH_IMAGE081
和电场迭代量
Figure 827359DEST_PATH_IMAGE082
对电场分量进行粗网格校正，进 而得到新的电场分量；
将校正后的新电场分量带入公式（5），并用公式（6）对其进行平滑更新，得到对应的细 网格上新的电场分量
Figure 265294DEST_PATH_IMAGE081
。
8.根据权利要求7所述的三维大地电磁正演数值模拟方法，其特征在于，在所述步骤 10）中，当细网格上的电场残差分量满足相对残差
Figure 623594DEST_PATH_IMAGE083
＜1×10-10时，停止迭代计算。
9.根据权利要求8所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤13）包括以下具体过程：
根据两种极化模式下的电场分量和磁场分量，计算出
Figure 96164DEST_PATH_IMAGE084
、
Figure 944034DEST_PATH_IMAGE085
两种模式的阻抗：
Figure 869265DEST_PATH_IMAGE086
(11)
式中：
Figure 280523DEST_PATH_IMAGE087
、
Figure 607599DEST_PATH_IMAGE088
、
Figure 891950DEST_PATH_IMAGE089
和
Figure 304477DEST_PATH_IMAGE090
为
Figure 66897DEST_PATH_IMAGE068
极化计算出
Figure 186162DEST_PATH_IMAGE001
、
Figure 828365DEST_PATH_IMAGE002
方向的电场分量和磁场分量，
Figure 728188DEST_PATH_IMAGE091
、
Figure 294299DEST_PATH_IMAGE092
、
Figure 268071DEST_PATH_IMAGE093
和
Figure 894224DEST_PATH_IMAGE094
为
Figure 281343DEST_PATH_IMAGE069
极化计算出
Figure 385566DEST_PATH_IMAGE001
、
Figure 728691DEST_PATH_IMAGE010
方向的电场分量和磁场分量；
根据所得阻抗计算出
Figure 525746DEST_PATH_IMAGE084
、
Figure 134582DEST_PATH_IMAGE085
两种模式下的视电阻率和阻抗相位：
Figure 308074DEST_PATH_IMAGE095
Figure 256438DEST_PATH_IMAGE096
(12)。
CN202111102810.4A 2021-09-22 2021-09-22 一种三维大地电磁正演数值模拟方法 Active CN113553748B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111102810.4A CN113553748B (zh) 2021-09-22 2021-09-22 一种三维大地电磁正演数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111102810.4A CN113553748B (zh) 2021-09-22 2021-09-22 一种三维大地电磁正演数值模拟方法

Publications (2)

Publication Number Publication Date
CN113553748A true CN113553748A (zh) 2021-10-26
CN113553748B CN113553748B (zh) 2021-11-30

Family

ID=78106417

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111102810.4A Active CN113553748B (zh) 2021-09-22 2021-09-22 一种三维大地电磁正演数值模拟方法

Country Status (1)

Country Link
CN (1) CN113553748B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113779818A (zh) * 2021-11-15 2021-12-10 中南大学 三维地质体其电磁场数值模拟方法、装置、设备及介质
CN113792445A (zh) * 2021-11-15 2021-12-14 中南大学 一种基于积分方程法的三维大地电磁数值模拟方法
CN114048823A (zh) * 2021-11-25 2022-02-15 成都理工大学 基于全卷积网络电阻率反演模型建立方法
CN114065511A (zh) * 2021-11-15 2022-02-18 中南大学 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
CN114297905A (zh) * 2022-03-10 2022-04-08 中南大学 一种二维大地电磁场的快速数值模拟方法
CN114912310A (zh) * 2022-04-11 2022-08-16 中南大学 一种基于正则化修正方程的三维大地电磁数值模拟方法
CN116842813A (zh) * 2023-09-04 2023-10-03 中南大学 一种三维大地电磁正演数值模拟方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180138330A1 (en) * 2009-08-18 2018-05-17 U.S. Army Research Laboratory Attn: Rdrl-Loc-I Photo detection device using resonance and related method
CN108229082A (zh) * 2018-04-12 2018-06-29 吉林大学 一种基于数据空间快速计算的联合反演方法
CN110068873A (zh) * 2019-05-10 2019-07-30 成都理工大学 一种基于球坐标系的大地电磁三维正演方法
CN110826283A (zh) * 2019-11-15 2020-02-21 中南大学 一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法
CN111323830A (zh) * 2020-01-14 2020-06-23 东华理工大学 一种基于大地电磁和直流电阻率数据的联合反演方法
CN112327374A (zh) * 2020-10-15 2021-02-05 广州市市政工程设计研究总院有限公司 Gpu探地雷达复杂介质dgtd正演方法
CN112666612A (zh) * 2020-11-02 2021-04-16 中国铁路设计集团有限公司 基于禁忌搜索的大地电磁二维反演方法
CN113221393A (zh) * 2021-01-29 2021-08-06 吉林大学 一种基于非结构有限元法的三维大地电磁各向异性反演方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180138330A1 (en) * 2009-08-18 2018-05-17 U.S. Army Research Laboratory Attn: Rdrl-Loc-I Photo detection device using resonance and related method
CN108229082A (zh) * 2018-04-12 2018-06-29 吉林大学 一种基于数据空间快速计算的联合反演方法
CN110068873A (zh) * 2019-05-10 2019-07-30 成都理工大学 一种基于球坐标系的大地电磁三维正演方法
CN110826283A (zh) * 2019-11-15 2020-02-21 中南大学 一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法
CN111323830A (zh) * 2020-01-14 2020-06-23 东华理工大学 一种基于大地电磁和直流电阻率数据的联合反演方法
CN112327374A (zh) * 2020-10-15 2021-02-05 广州市市政工程设计研究总院有限公司 Gpu探地雷达复杂介质dgtd正演方法
CN112666612A (zh) * 2020-11-02 2021-04-16 中国铁路设计集团有限公司 基于禁忌搜索的大地电磁二维反演方法
CN113221393A (zh) * 2021-01-29 2021-08-06 吉林大学 一种基于非结构有限元法的三维大地电磁各向异性反演方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
JIAN LI ET AL.: "An Efficient Preconditioner for 3-D Finite Difference Modeling of the Electromagnetic Diffusion Process in the Frequency Domain", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
李健，等: "一种频率域三维有限差分电磁正演的高效预处理器", 《中国地球科学联合学术年会 2019》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113779818A (zh) * 2021-11-15 2021-12-10 中南大学 三维地质体其电磁场数值模拟方法、装置、设备及介质
CN113792445A (zh) * 2021-11-15 2021-12-14 中南大学 一种基于积分方程法的三维大地电磁数值模拟方法
CN113779818B (zh) * 2021-11-15 2022-02-08 中南大学 三维地质体其电磁场数值模拟方法、装置、设备及介质
CN114065511A (zh) * 2021-11-15 2022-02-18 中南大学 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
CN114065511B (zh) * 2021-11-15 2024-08-13 中南大学 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
CN114048823A (zh) * 2021-11-25 2022-02-15 成都理工大学 基于全卷积网络电阻率反演模型建立方法
CN114297905A (zh) * 2022-03-10 2022-04-08 中南大学 一种二维大地电磁场的快速数值模拟方法
CN114297905B (zh) * 2022-03-10 2022-06-03 中南大学 一种二维大地电磁场的快速数值模拟方法
CN114912310A (zh) * 2022-04-11 2022-08-16 中南大学 一种基于正则化修正方程的三维大地电磁数值模拟方法
CN114912310B (zh) * 2022-04-11 2024-04-12 中南大学 一种基于正则化修正方程的三维大地电磁数值模拟方法
CN116842813A (zh) * 2023-09-04 2023-10-03 中南大学 一种三维大地电磁正演数值模拟方法
CN116842813B (zh) * 2023-09-04 2023-11-14 中南大学 一种三维大地电磁正演数值模拟方法

Also Published As

Publication number Publication date
CN113553748B (zh) 2021-11-30

Similar Documents

Publication Publication Date Title
CN113553748B (zh) 一种三维大地电磁正演数值模拟方法
WO2022227206A1 (zh) 一种基于全卷积神经网络的大地电磁反演方法
Spitzer A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods
CN106199742B (zh) 一种频率域航空电磁法2.5维带地形反演方法
CN113408167B (zh) 基于场路耦合的直流偏磁计算方法
CN112949134B (zh) 基于非结构有限元方法的地-井瞬变电磁反演方法
CN109977585B (zh) 一种高精度大地电磁正演模拟方法
CN102798898A (zh) 大地电磁场非线性共轭梯度三维反演方法
CN108509693A (zh) 三维频率域可控源数值模拟方法
CN113569447B (zh) 一种基于舒尔补方法的瞬变电磁三维快速正演方法
CN115238550B (zh) 自适应非结构网格的滑坡降雨的地电场数值模拟计算方法
CN114065511B (zh) 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
CN117538945B (zh) 三维大地电磁多分辨率反演方法、装置、设备及介质
CN113255230A (zh) 基于mq径向基函数的重力模型正演方法及系统
CN116522713A (zh) 基于分块坐标下降算法的航空电磁数据分区三维反演方法
CN110346834A (zh) 三维频率域可控源电磁的正演方法、系统
CN113779818B (zh) 三维地质体其电磁场数值模拟方法、装置、设备及介质
CN108873084B (zh) 一种基于单位分解积分的直流电阻率无单元正演方法
CN114970289B (zh) 三维大地电磁各向异性正演数值模拟方法、设备及介质
CN110826283A (zh) 一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法
CN114065577A (zh) 一种直流电阻率小波伽辽金三维正演方法
CN114297905B (zh) 一种二维大地电磁场的快速数值模拟方法
Günther et al. Electrical Resistivity Tomography (ERT) in geophysical applications-state of the art and future challenges
Tsourlos et al. Efficient 2D inversion of long ERT sections
CN115563791A (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