CN116842813B - 一种三维大地电磁正演数值模拟方法 - Google Patents
一种三维大地电磁正演数值模拟方法 Download PDFInfo
- Publication number
- CN116842813B CN116842813B CN202311126708.7A CN202311126708A CN116842813B CN 116842813 B CN116842813 B CN 116842813B CN 202311126708 A CN202311126708 A CN 202311126708A CN 116842813 B CN116842813 B CN 116842813B
- Authority
- CN
- China
- Prior art keywords
- electric field
- representing
- curved surface
- equation
- resistivity
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 35
- 230000005684 electric field Effects 0.000 claims abstract description 69
- 230000010287 polarization Effects 0.000 claims abstract description 26
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 210000004027 cell Anatomy 0.000 claims description 27
- 230000006870 function Effects 0.000 claims description 19
- 230000014509 gene expression Effects 0.000 claims description 18
- 239000011159 matrix material Substances 0.000 claims description 13
- 239000013598 vector Substances 0.000 claims description 12
- 238000013507 mapping Methods 0.000 claims description 9
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 5
- 230000005672 electromagnetic field Effects 0.000 claims description 5
- 230000002159 abnormal effect Effects 0.000 claims description 3
- 238000002224 dissection Methods 0.000 claims description 3
- 230000035699 permeability Effects 0.000 claims description 3
- 238000000926 separation method Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims description 2
- 238000003379 elimination reaction Methods 0.000 claims description 2
- 238000010586 diagram Methods 0.000 description 9
- 238000004590 computer program Methods 0.000 description 6
- 238000011160 research Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000005259 measurement Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000001125 extrusion Methods 0.000 description 1
- 239000003673 groundwater Substances 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
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)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Computer Hardware Design (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Measuring Magnetic Variables (AREA)
Abstract
本发明提出一种三维大地电磁正演数值模拟方法，具体是：依据目标地质体的特征构建任意各向异性电阻率模型；基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度将地质体模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值，依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程；依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；采用正则化约束下的正演方程获得两种极化模式下的电场分量和磁场分量，继而得到视电阻率和阻抗相位。本发明既能够真实模拟三维大地电磁情况，又能保证高精度求解的同时提高电磁正演的计算效率。
Description
技术领域
本发明属于数值模拟技术领域，特别涉及一种三维大地电磁正演数值模拟方法。
背景技术
大地电磁法利用天然源电磁场研究地下地电结构，通过测量不同频率的电磁场达到研究不同深度的地质构造的目的，具有勘探深度较大、分辨率高、较小的高阻屏蔽、工作效率高等优点，在深部地球动力学研究、矿产与油气勘探、地热资源勘探、地下水资源勘探等领域广泛应用。
传统的三维大地电磁法正演主要采用笛卡尔坐标系和平面波场源假设，在进行跨大陆尺度的大地电磁法勘探时，传统的大地电磁法正演受到地球曲率的影响使得正演结果误差较大，而采用球坐标系下的大地电磁法正演可以有效解决地球曲率带来的问题。
对于实际地球介质而言，其经过持续的分异、交代以及挤压等复杂地质作用使得其内部结构、物质运移及赋存状态普遍存在各向异性介质。而目前在进行球坐标系下三维大地电磁法正演时，并没有考虑这种各向异性对大尺度电磁正演模拟的影响，因此，正演结果误差较大，不能真实反映三维大地电磁的情况。为获得大尺度下精确的大地电磁三维正演结果，提高大地电磁反演解释的可靠性，很有必要进行基于球坐标系的大地电磁各向异性正演研究。
综上所述，急需一种操作方便且能够真实反映三维大地电磁的情况的三维大地电磁正演数值模拟方法以解决现有技术中存在的问题。
发明内容
本发明的目的在于提供一种三维大地电磁正演数值模拟方法，具体是：依据目标地质体的特征构建任意各向异性电阻率模型；基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度将地质体模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值，依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程；依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；采用正则化约束下的正演方程获得两种极化模式下的电场分量和磁场分量，继而得到视电阻率和阻抗相位。采用上述方案既能够真实模拟三维大地电磁情况，又能保证高精度求解的同时提高电磁正演的计算效率。具体技术方案如下：
一种三维大地电磁正演数值模拟方法，步骤如下：
步骤一、依据目标地质体的特征构建任意各向异性电阻率模型；
步骤二、基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度/>将任意各向异性电阻率模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值；
步骤三、依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程，如下：
；
其中：为哈密顿算子；/>表示电场强度；/>表示磁导率；/>表示角频率；/>；表示任意各向异性电导率，用3×3的电导率张量表示，/>；
其中，表示/>轴方向的电导率值；/>表示/>轴方向的电导率值；/>表示/>轴方向的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；电导率张量/>是一个对称且正定的矩阵，有/>；
步骤四、在双旋度方程中加入散度校正项，得到正则化约束下的正演方程如下：
；
其中：为正则化因子；
步骤五、依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；
步骤六、基于交错网格有限差分法对双旋度方程离散，并求解该双旋度方程，得到曲面单元上对应的电场分量；
步骤七、根据曲面单元上的电场分量求旋度，计算出不同极化模式下曲面单元上的磁场分量；
步骤八、根据曲面单元上不同极化模式下的电场分量和磁场分量计算出对应测点的视电阻率和阻抗相位。
优选的，步骤五具体如下：
将地球及其上覆空气层划分为球体层状电导率模型，每一层电场的散度为零，则电场的极向-环向分解表达式如下：
；
其中：是东西方向电势，/>为南北方向电势；/>为球体半径，/>为球体半径单位矢量，/>为梯度算子的角分量，/>；/>为拉普拉斯算子的角分量，；
将代入双旋度方程中，经过代数变换后，使得电势满足亥姆霍兹方程如下：
；
其中：代表东西方向电势/>或南北方向电势/>；/>为复波数；
对亥姆霍兹方程通过分离变量变成如下表达式：
；
其中：和/>表示球谐系数；/>表示球谐次数；/>表示球谐阶数；/>为第一类贝塞尔函数，/>为第三类贝塞尔函数，与贝塞尔球谐函数有关；/>是第/>阶球谐函数的归一化线性组合；
空气层通常满足准静态条件，则将贝塞尔函数用常数因子代替转换成下式：
；
电磁场球谐函数分解成东西方向电势和南北方向电势/>，将带入电场的极向－环向分解表达式中，消去东西方向电势/>相关的因式项或消去南北方向电势/>相关的因式项，得到准静态条件下电场的表达式：
；
其中：表示施密特归一化缔合勒让德函数；/>表示/>阶和/>阶对应的贝塞尔函数。
优选的，所述步骤六具体是：
在球坐标系下采用交错网格有限差分法对正则化约束下的正演方程进行离散，其具体形式如下：
；
其中：表示离散旋度算子，表示将单元棱边上的矢量映射到单元的曲面上；/>是/>的转置矩阵，表示将单元面上的矢量映射回到单元棱边上；/>；/>表示加权因子的对角矩阵；/>表示离散梯度算子，表示将单元节点上的量映射到单元的内部边上，即/>，/>表示单元体的曲边边长，/>表示离散梯度算子的拓扑矩阵；离散算子/>则表示将单元内部边上的量映射到内部节点上，通过对/>转置得到；/>表示对角矩阵，为单元棱边的电导率，/>表示离散电场矢量；
直到曲面单元上的电场残差分量满足相对残差小于时，停止迭代计算，得到电场分量。
优选的，所述步骤一具体是：所述目标地质体的特征包括目标地质体的形状、大小和电阻率分布；所述目标地质体为三维异常体。
优选的，所述步骤七中磁场分量如下表达式：
。
优选的，所述步骤八中：根据两种极化模式下的电场分量和磁场分量，计算出两种模式的阻抗：
；
其中：和/>为/>极化计算出/>方向的电场分量和磁场分量，和/>为/>极化计算出/>方向的电场分量和磁场分量；
视电阻率和阻抗相位计算表达式如下：
；
其中，表示/>和/>两种模式下的视电阻率；/>表示/>和/>两种模式下的相位；/>和/>表示虚部和实部；/>表示反三角函数。
应用本发明的技术方案，具有的有益效果是：本发明在球坐标系的大地电磁法中充分考虑了各向异性的影响，各向异性地质模型可以更加真实地描绘地球介质的赋存状态，考虑电各向异性的影响能够提高大地电磁法反演解释的可靠性；本发明基于球坐标系对任意各向异性电阻率模型使用曲面单元进行剖分，采用交错网格有限差分法离散曲面单元上电场满足的双旋度方程，并且采用了适用于全球性、半全球性电性结构的非均一场源替代了传统大地电磁勘探采用的平面波场源，避免了地球曲率对大尺度下的大地电磁法正演的影响；本发明采用了正则化约束下的正演方程，即将散度校正项添加到原始控制方程中来对控制方程进行约束，以保证每次迭代电流的散度为零，此方法避免了额外求解散度方程，可以显著提高球坐标系下大地电磁法正演效率。
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为本实施例中球坐标系网格剖分示意图，其中：(a)为模型剖分示意图，(b)为曲面单元示意图；
图2为本实施例中任意各向异性半空间模型示意图；
图3为本实施例得到的极化方向视电阻率和阻抗相位数值解与解析解对比图，其中：(a)为/>极化方向视电阻率数值解与解析解对比图，(b)为/>极化方向相位数值解与解析解对比图，(c)为/>极化方向视电阻率数值解与解析解相对误差图；(d)为/>极化方向相位数值解与解析解绝对误差图；
图4为本实施例得到的极化方向视电阻率和阻抗相位数值解与解析解对比图，其中：(a)为/>极化方向视电阻率数值解与解析解对比图，(b)为/>极化方向相位数值解与解析解对比图，(c)为/>极化方向视电阻率数值解与解析解相对误差图；(d)为/>极化方向相位数值解与解析解绝对误差图；
图5为本实施例中的计算机设备的示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述来清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
实施例：
一种三维大地电磁正演数值模拟方法，尤其适用于大尺度、大深度的大地电磁勘探，具体包括以下步骤：
步骤一、依据目标地质体的特征构建任意各向异性电阻率模型，所述目标地质体的特征包括目标地质体的形状、大小和电阻率分布；所述目标地质体为三维异常体。
步骤二、基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度/>将任意各向异性电阻率模型以及上覆空气层剖分成若干个曲面单元，详见图1，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值。此处划分为/>个曲面单元，分别为/>、/>、/>三个方向曲面单元的个数。
步骤三、依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程，如下：
；
其中：为哈密顿算子；/>表示电场强度；/>表示磁导率；/>表示角频率；/>；表示任意各向异性电导率，用3×3的电导率张量表示，/>，/>表示/>轴方向的电导率值；/>表示/>轴方向的电导率值；/>表示/>轴方向的电导率值；表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；电导率张量/>是一个对称且正定的矩阵，有。而正定的矩阵都可以通过三次矩阵旋转得到主参考系下的电导率张量，此时仅有对角线元素，其余元素均为零，这样的转换亦便于我们正演模型的生成。即有：
；
其中： ，/>代表/>以及/>其中一个变量，/>分别为三次矩阵旋转所对应的旋转角，即各向异性走向角、倾向角、倾斜角。故有六个量可表示电导率张量。
步骤四、在双旋度方程中加入散度校正项，得到正则化约束下的正演方程如下：
；
其中：为正则化因子。
通过将散度校正项添加到原始控制方程中来对控制方程进行约束，能有效地避免对散度方程额外的求解（使得迭代求解器的收敛性得到最显著的改善），保证每次迭代电流的散度为零。
步骤五、依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；
本实施例中，依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件，如下：
将地球及其上覆空气层划分为球体层状电导率模型，每一层电场的散度为零，则电场的极向-环向分解表达式如下：
；
其中：是东西方向电势，/>为南北方向电势，/>为球体半径，/>为球体半径单位矢量，/>为梯度算子的角分量，/>；/>为拉普拉斯算子的角分量，；
将代入双旋度方程中，经过代数变换后，使得电势满足亥姆霍兹方程如下：
；
其中：代表东西方向电势/>或南北方向电势/>；/>为复波数；
对亥姆霍兹方程通过分离变量变成如下表达式：
；
其中：和/>表示球谐系数；/>表示球谐次数；/>表示球谐阶数；/>为第一类贝塞尔函数，/>为第三类贝塞尔函数，与贝塞尔球谐函数有关；/>是第/>阶球谐函数的归一化线性组合；
空气层通常满足准静态条件，则将贝塞尔函数用常数因子代替转换成下式：
；
电磁场球谐函数分解成东西方向电势和南北方向电势/>，将带入电场的极向－环向分解表达式中，消去东西方向电势/>相关的因式项或消去东西方向电势/>相关的因式项，得到准静态条件下电场的表达式：
；
其中：表示施密特归一化缔合勒让德函数；/>表示/>阶和/>阶对应的贝塞尔函数。
步骤六、基于交错网格有限差分法对双旋度方程离散，并求解该双旋度方程，得到曲面单元上对应的电场分量。
本实施例中，在球坐标系下采用交错网格有限差分法对正则化约束下的正演方程进行离散，其具体形式如下：
；
其中：表示离散旋度算子，表示将单元棱边上的矢量映射到单元的曲面上，/>是/>的伴随矩阵，表示将单元面上的矢量映射回到单元棱边上；/>，/>表示加权因子的对角矩阵，/>表示离散梯度算子，表示将单元节点上的量映射到单元的内部边上，即/>，/>表示单元体的曲边边长，/>表示离散梯度算子的拓扑矩阵；离散算子/>则表示将单元内部边上的量映射到内部节点上，通过对/>转置得到；/>表示对角矩阵，为单元棱边的电导率，/>表示离散电场矢量。
直到曲面单元上的电场残差分量满足相对残差小于时，停止迭代计算，得到电场分量。
步骤七、根据曲面单元上的电场分量求旋度，计算出不同极化模式下曲面单元上的磁场分量如下表达式：
。
步骤八、根据曲面单元上不同极化模式下的电场分量和磁场分量计算出对应测点的视电阻率和阻抗相位，具体是：
根据两种极化模式下的电场分量和磁场分量，计算出两种模式的阻抗：
；
其中：和/>为/>极化计算出/>方向的电场分量和磁场分量，和/>为/>极化计算出/>方向的电场分量和磁场分量；
视电阻率和阻抗相位计算表达式如下：
；
其中，表示/>和/>两种模式下的视电阻率；/>表示/>和/>两种模式下的相位；/>和/>表示虚部和实部；/>表示反三角函数。
应用本实施例的方案，具体是：设计了如图2所示任意各向异性半空间模型，三个主轴方向的电阻率分别为：主轴电阻率ρ1=100Ωm（=1/100），/>主轴电阻率ρ2=500Ωm（=1/500），/>主轴电阻率ρ3=300Ωm（/>=1/300）；主轴坐标系到测量坐标系的旋转角分别为：各向异性的走向角/>=25°，各向异性的倾向角/>=60°，各向异性的倾斜角/>=15°；空气层电阻率/>，模拟区域为/>和/>方向均为-1.6°到1.6°，/>方向为0km到137km，将该模型区域剖分成64×64×64个曲面单元，采取不规则网格剖分。
为了说明对比结果。选取了位于剖面中部的一个测点展示数值误差情况，图3和图4为两种极化模式下的视电阻率和阻抗相位对比曲线图以及误差图，计算周期设置为10s到1000s等对数间隔。由图3和图4可知本发明方法（Spherical）与通用的各向异性解析解（Analytic）在数值模拟结果上基本吻合，均符合绝对误差、相位误差小于5%的要求，从而验证了本发明方法在结果精度上的可靠性。
另外，本实施例还公开了一种计算机设备，该计算机设备可以是服务器，其内部结构图可以如图5所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口和数据库。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的数据库用于存储样本数据。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现上述三维大地电磁正演数值模拟方法。
本领域技术人员可以理解，图5中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
另外，本实施例还公开了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述实施例中三维大地电磁正演数值模拟方法的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（ROM）、可编程ROM（PROM）、电可编程ROM（EPROM）、电可擦除可编程ROM（EEPROM）或闪存。易失性存储器可包括随机存取存储器（RAM）或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM（SRAM）、动态RAM（DRAM）、同步DRAM（SDRAM）、双数据率SDRAM（DDRSDRAM）、增强型SDRAM（ESDRAM）、同步链路（Synchlink） DRAM（SLDRAM）、存储器总线（Rambus）直接RAM（RDRAM）、直接存储器总线动态RAM（DRDRAM）、以及存储器总线动态RAM（RDRAM）等。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。
Claims (6)
1.一种三维大地电磁正演数值模拟方法，其特征在于，步骤如下：
步骤一、依据目标地质体的特征构建任意各向异性电阻率模型；
步骤二、基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度/>和经度/>将任意各向异性电阻率模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值；
步骤三、依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程，如下：
；
其中：为哈密顿算子；/>表示电场强度；/>表示磁导率；/>表示角频率；/>；/>表示任意各向异性电导率，用3×3的电导率张量表示，/>，/>表示/>轴方向的电导率值；/>表示/>轴方向的电导率值；/>表示/>轴方向的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；/>表示当/>方向施加电场，在/>方向会形成电流密度的电导率值；
步骤四、在双旋度方程中加入散度校正项，得到正则化约束下的正演方程如下：
；
其中：为正则化因子；
步骤五、依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；
步骤六、基于交错网格有限差分法对双旋度方程离散，并求解该双旋度方程，得到曲面单元上对应的电场分量；
步骤七、根据曲面单元上的电场分量求旋度，计算出不同极化模式下曲面单元上的磁场分量；
步骤八、根据曲面单元上不同极化模式下的电场分量和磁场分量计算出对应测点的视电阻率和阻抗相位。
2.根据权利要求1所述的三维大地电磁正演数值模拟方法，其特征在于，步骤五具体如下：
将地球及其上覆空气层划分为球体层状电导率模型，每一层电场的散度为零，则电场的极向-环向分解表达式如下：
；
其中：是东西方向电势，/>为南北方向电势；/>为球体半径， />为球体半径单位矢量，为梯度算子的角分量，/>；/>为拉普拉斯算子的角分量，；
将代入双旋度方程中，经过代数变换后，使得电势满足亥姆霍兹方程如下：
；
其中：代表东西方向电势/>或南北方向电势/>；/>为复波数；
对亥姆霍兹方程通过分离变量变成如下表达式：
；
其中：和/>表示球谐系数；/>表示球谐次数；/>表示球谐阶数；/>为第一类贝塞尔函数，/>为第三类贝塞尔函数，与贝塞尔球谐函数有关；/>是第/>阶球谐函数的归一化线性组合；
空气层通常满足准静态条件，则将贝塞尔函数用常数因子代替转换成下式：
；
电磁场球谐函数分解成东西方向电势和南北方向电势/>，将带入电场的极向－环向分解表达式中，消去东西方向电势/>相关的因式项或消去南北方向电势/>相关的因式项，得到准静态条件下电场的表达式，即大地电磁法边界电场值/>：
；
其中：表示施密特归一化缔合勒让德函数；/>表示/>阶和/>阶对应的贝塞尔函数。
3.根据权利要求1-2任意一项所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤六具体是：
在球坐标系下采用交错网格有限差分法对正则化约束下的正演方程进行离散，其具体形式如下：
；
其中：表示离散旋度算子，表示将单元棱边上的矢量映射到单元的曲面上；/>是/>的转置矩阵，表示将单元面上的矢量映射回到单元棱边上；/>；/>表示加权因子的对角矩阵；/>表示离散梯度算子，表示将单元节点上的量映射到单元的内部边上，即，/>表示单元体的曲边边长，/>表示离散梯度算子的拓扑矩阵；离散算子/>则表示将单元内部边上的量映射到内部节点上，通过对/>转置得到；/>表示对角矩阵，/>为单元棱边的电导率，/>表示离散电场矢量；
直到曲面单元上的电场残差分量满足相对残差小于时，停止迭代计算，得到电场分量。
4.根据权利要求3所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤一具体是：所述目标地质体的特征包括目标地质体的形状、大小和电阻率分布；所述目标地质体为三维异常体。
5.根据权利要求3所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤七中磁场分量如下表达式：
。
6.根据权利要求5所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤八中：根据两种极化模式下的电场分量和磁场分量，计算出两种模式的阻抗：
；
其中：和/>为/>极化计算出/>方向的电场分量和磁场分量，和/>为/>极化计算出/>方向的电场分量和磁场分量；
视电阻率和阻抗相位计算表达式如下：
；
其中，表示/>和/>两种模式下的视电阻率；/>表示/>和/>两种模式下的相位；/>和/>表示虚部和实部；/>表示反三角函数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311126708.7A CN116842813B (zh) | 2023-09-04 | 2023-09-04 | 一种三维大地电磁正演数值模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311126708.7A CN116842813B (zh) | 2023-09-04 | 2023-09-04 | 一种三维大地电磁正演数值模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116842813A CN116842813A (zh) | 2023-10-03 |
| CN116842813B true CN116842813B (zh) | 2023-11-14 |
Family
ID=88165628
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311126708.7A Active CN116842813B (zh) | 2023-09-04 | 2023-09-04 | 一种三维大地电磁正演数值模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116842813B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117538945B (zh) * | 2024-01-10 | 2024-03-26 | 中南大学 | 三维大地电磁多分辨率反演方法、装置、设备及介质 |
| CN118013768B (zh) * | 2024-04-10 | 2024-07-12 | 中国科学院地质与地球物理研究所 | 一种行星岩石圈磁场模型系数的确定方法及装置 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110058315A (zh) * | 2019-05-29 | 2019-07-26 | 中南大学 | 一种三维各向异性射频大地电磁自适应有限元正演方法 |
| CN110068873A (zh) * | 2019-05-10 | 2019-07-30 | 成都理工大学 | 一种基于球坐标系的大地电磁三维正演方法 |
| CN111611737A (zh) * | 2020-05-19 | 2020-09-01 | 中南大学 | 一种三维任意各向异性介质的海洋可控源电磁正演方法 |
| CN113553748A (zh) * | 2021-09-22 | 2021-10-26 | 中南大学 | 一种三维大地电磁正演数值模拟方法 |
| CN114036745A (zh) * | 2021-11-08 | 2022-02-11 | 中南大学 | 各向异性介质三维大地电磁正演方法及系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| MX2007015622A (es) * | 2005-06-09 | 2008-02-25 | Exxonmobil Upstream Res Co | Metodo para determinar la anisotropia electrica vertical terrestre en levantamientos electromagneticos marinos. |
-
2023
- 2023-09-04 CN CN202311126708.7A patent/CN116842813B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110068873A (zh) * | 2019-05-10 | 2019-07-30 | 成都理工大学 | 一种基于球坐标系的大地电磁三维正演方法 |
| CN110058315A (zh) * | 2019-05-29 | 2019-07-26 | 中南大学 | 一种三维各向异性射频大地电磁自适应有限元正演方法 |
| CN111611737A (zh) * | 2020-05-19 | 2020-09-01 | 中南大学 | 一种三维任意各向异性介质的海洋可控源电磁正演方法 |
| CN113553748A (zh) * | 2021-09-22 | 2021-10-26 | 中南大学 | 一种三维大地电磁正演数值模拟方法 |
| CN114036745A (zh) * | 2021-11-08 | 2022-02-11 | 中南大学 | 各向异性介质三维大地电磁正演方法及系统 |
Non-Patent Citations (4)
| Title |
|---|
| 《A Scalable Parallel Algorithm for 3-D Magnetotelluric Finite Element Modeling in Anisotropic Media》;Xiaoxiong Zhu;《IEEE Transactions on Geoscience and Remote Sensing》;全文 * |
| 《An Efficient Preconditioner for 3-D Finite Difference Modeling of the Electromagnetic Diffusion Process in the Frequency Domain》;Jian Li，Jianxin Liu等;《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》;全文 * |
| 三维陆地可控源电磁法有限元模型降阶快速正演;张继锋;刘寄仁;冯兵;郑一安;;地球物理学报(第09期);全文 * |
| 面向目标自适应有限元法的带地形三维大地电磁各向异性正演模拟;曹晓月;殷长春;张博;黄鑫;刘云鹤;蔡晶;;地球物理学报(第06期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116842813A (zh) | 2023-10-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116842813B (zh) | 一种三维大地电磁正演数值模拟方法 | |
| Key et al. | Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example | |
| Schillinger et al. | The hp‐d‐adaptive finite cell method for geometrically nonlinear problems of solid mechanics | |
| CN112287534B (zh) | 基于nufft的二维磁异常快速正演模拟方法和装置 | |
| CN106199742B (zh) | 一种频率域航空电磁法2.5维带地形反演方法 | |
| CN110346835B (zh) | 大地电磁的正演方法、正演系统、存储介质及电子设备 | |
| Peixoto et al. | Analysis of grid imprinting on geodesic spherical icosahedral grids | |
| CN110346834B (zh) | 三维频率域可控源电磁的正演方法、系统 | |
| US20200349306A1 (en) | General Scattered Field Simulator | |
| Ren et al. | Recursive analytical formulae of gravitational fields and gradient tensors for polyhedral bodies with polynomial density contrasts of arbitrary non-negative integer orders | |
| CN111103628B (zh) | 大地电磁te极化模式对电场数据二维反演方法和装置 | |
| Barbosa et al. | The direct interpolation boundary element technique applied to three-dimensional scalar free vibration problems | |
| Zhu et al. | 3-D dc resistivity modelling based on spectral element method with unstructured tetrahedral grids | |
| CN113076678B (zh) | 一种频率域二度体重力异常快速数值模拟方法和装置 | |
| CN117538945A (zh) | 三维大地电磁多分辨率反演方法、装置、设备及介质 | |
| CN113792445B (zh) | 一种基于积分方程法的三维大地电磁数值模拟方法 | |
| Dai et al. | The forward modeling of 3D gravity and magnetic potential fields in space-wavenumber domains based on an integral method | |
| Song et al. | Finite element method for modeling 3D resistivity sounding on anisotropic geoelectric media | |
| Spitzer | Electromagnetic modeling using adaptive grids–Error estimation and geometry representation | |
| Jahandari et al. | Forward modeling of direct-current resistivity data on unstructured grids using an adaptive mimetic finite-difference method | |
| CN113627027B (zh) | 非平凡各向异性介质电磁场数值模拟方法及系统 | |
| CN113779818B (zh) | 三维地质体其电磁场数值模拟方法、装置、设备及介质 | |
| CN114970289A (zh) | 三维大地电磁各向异性正演数值模拟方法、设备及介质 | |
| CN114036806A (zh) | 基于热导率各向异性介质的三维地温场数值模拟方法 | |
| CN113673163A (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 |