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

1.本发明属于数值模拟技术领域，特别涉及一种三维大地电磁正演数值模拟方法。


背景技术：

2.大地电磁法利用天然源电磁场研究地下地电结构，通过测量不同频率的电磁场达到研究不同深度的地质构造的目的，具有勘探深度较大、分辨率高、较小的高阻屏蔽、工作效率高等优点，在深部地球动力学研究、矿产与油气勘探、地热资源勘探、地下水资源勘探等领域广泛应用。
3.传统的三维大地电磁法正演主要采用笛卡尔坐标系和平面波场源假设，在进行跨大陆尺度的大地电磁法勘探时，传统的大地电磁法正演受到地球曲率的影响使得正演结果误差较大，而采用球坐标系下的大地电磁法正演可以有效解决地球曲率带来的问题。
4.对于实际地球介质而言，其经过持续的分异、交代以及挤压等复杂地质作用使得其内部结构、物质运移及赋存状态普遍存在各向异性介质。而目前在进行球坐标系下三维大地电磁法正演时，并没有考虑这种各向异性对大尺度电磁正演模拟的影响，因此，正演结果误差较大，不能真实反映三维大地电磁的情况。为获得大尺度下精确的大地电磁三维正演结果，提高大地电磁反演解释的可靠性，很有必要进行基于球坐标系的大地电磁各向异性正演研究。
5.综上所述，急需一种操作方便且能够真实反映三维大地电磁的情况的三维大地电磁正演数值模拟方法以解决现有技术中存在的问题。


技术实现要素：

6.本发明的目的在于提供一种三维大地电磁正演数值模拟方法，具体是：依据目标地质体的特征构建任意各向异性电阻率模型；基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度将地质体模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值，依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程；依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；采用正则化约束下的正演方程获得两种极化模式下的电场分量和磁场分量，继而得到视电阻率和阻抗相位。采用上述方案既能够真实模拟三维大地电磁情况，又能保证高精度求解的同时提高电磁正演的计算效率。具体技术方案如下：一种三维大地电磁正演数值模拟方法，步骤如下：步骤一、依据目标地质体的特征构建任意各向异性电阻率模型；步骤二、基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度将任意各向异性电阻率模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值；步骤三、依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足
的双旋度方程，如下：；其中：为哈密顿算子；表示电场强度；表示磁导率；表示角频率；；表示任意各向异性电导率，用3
×
3的电导率张量表示，；其中，表示轴方向的电导率值；表示轴方向的电导率值；表示轴方向的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；电导率张量是一个对称且正定的矩阵，有；步骤四、在双旋度方程中加入散度校正项，得到正则化约束下的正演方程如下：；其中：为正则化因子；步骤五、依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；步骤六、基于交错网格有限差分法对双旋度方程离散，并求解该双旋度方程，得到曲面单元上对应的电场分量；步骤七、根据曲面单元上的电场分量求旋度，计算出不同极化模式下曲面单元上的磁场分量；步骤八、根据曲面单元上不同极化模式下的电场分量和磁场分量计算出对应测点的视电阻率和阻抗相位。
7.优选的，步骤五具体如下：将地球及其上覆空气层划分为球体层状电导率模型，每一层电场的散度为零，则电场的极向-环向分解表达式如下：；其中：是东西方向电势，为南北方向电势；为球体半径，为球体半径单位矢量，为梯度算子的角分量，；为拉普拉斯算子的角分量，
；将代入双旋度方程中，经过代数变换后，使得电势满足亥姆霍兹方程如下：；其中：代表东西方向电势或南北方向电势；为复波数；对亥姆霍兹方程通过分离变量变成如下表达式：；其中：和表示球谐系数；表示球谐次数；表示球谐阶数；为第一类贝塞尔函数，为第三类贝塞尔函数，与贝塞尔球谐函数有关；是第阶球谐函数的归一化线性组合；空气层通常满足准静态条件，则将贝塞尔函数用常数因子代替转换成下式：；电磁场球谐函数分解成东西方向电势和南北方向电势，将带入电场的极向－环向分解表达式中，消去东西方向电势相关的因式项或消去南北方向电势相关的因式项，得到准静态条件下电场的表达式：；其中：表示施密特归一化缔合勒让德函数；表示阶和阶对应的贝塞尔函数。
8.优选的，所述步骤六具体是：在球坐标系下采用交错网格有限差分法对正则化约束下的正演方程进行离散，其具体形式如下：；其中：表示离散旋度算子，表示将单元棱边上的矢量映射到单元的曲面上；是的转置矩阵，表示将单元面上的矢量映射回到单元棱边上；；表示加权因子的对角矩阵；表示离散梯度算子，表示将单元节点上的量映射到单元的内部边上，即，表示单元体的曲边边长，表示离散梯度算子的拓扑矩阵；离散算子
则表示将单元内部边上的量映射到内部节点上，通过对转置得到；表示对角矩阵，为单元棱边的电导率，表示离散电场矢量；直到曲面单元上的电场残差分量满足相对残差小于时，停止迭代计算，得到电场分量。
9.优选的，所述步骤一具体是：所述目标地质体的特征包括目标地质体的形状、大小和电阻率分布；所述目标地质体为三维异常体。
10.优选的，所述步骤七中磁场分量如下表达式：。
11.优选的，所述步骤八中：根据两种极化模式下的电场分量和磁场分量，计算出两种模式的阻抗：；其中：和为极化计算出方向的电场分量和磁场分量，和为极化计算出方向的电场分量和磁场分量；视电阻率和阻抗相位计算表达式如下：；其中，表示和两种模式下的视电阻率；表示和两种模式下的相位；和表示虚部和实部；表示反三角函数。
12.应用本发明的技术方案，具有的有益效果是：本发明在球坐标系的大地电磁法中充分考虑了各向异性的影响，各向异性地质模型可以更加真实地描绘地球介质的赋存状态，考虑电各向异性的影响能够提高大地电磁法反演解释的可靠性；本发明基于球坐标系对任意各向异性电阻率模型使用曲面单元进行剖分，采用交错网格有限差分法离散曲面单元上电场满足的双旋度方程，并且采用了适用于全球性、半全球性电性结构的非均一场源替代了传统大地电磁勘探采用的平面波场源，避免了地球曲率对大尺度下的大地电磁法正演的影响；本发明采用了正则化约束下的正演方程，即将散度校正项添加到原始控制方程中来对控制方程进行约束，以保证每次迭代电流的散度为零，此方法避免了额外求解散度方程，可以显著提高球坐标系下大地电磁法正演效率。
13.除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
14.构成本技术的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：图1为本实施例中球坐标系网格剖分示意图，其中：(a)为模型剖分示意图，(b)为曲面单元示意图；图2为本实施例中任意各向异性半空间模型示意图；图3为本实施例得到的极化方向视电阻率和阻抗相位数值解与解析解对比图，其中：(a)为极化方向视电阻率数值解与解析解对比图，(b)为极化方向相位数值解与解析解对比图，(c)为极化方向视电阻率数值解与解析解相对误差图；(d)为极化方向相位数值解与解析解绝对误差图；图4为本实施例得到的极化方向视电阻率和阻抗相位数值解与解析解对比图，其中：(a)为极化方向视电阻率数值解与解析解对比图，(b)为极化方向相位数值解与解析解对比图，(c)为极化方向视电阻率数值解与解析解相对误差图；(d)为极化方向相位数值解与解析解绝对误差图；图5为本实施例中的计算机设备的示意图。
具体实施方式
15.为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述来清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
实施例：
16.一种三维大地电磁正演数值模拟方法，尤其适用于大尺度、大深度的大地电磁勘探，具体包括以下步骤：步骤一、依据目标地质体的特征构建任意各向异性电阻率模型，所述目标地质体的特征包括目标地质体的形状、大小和电阻率分布；所述目标地质体为三维异常体。
17.步骤二、基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度将任意各向异性电阻率模型以及上覆空气层剖分成若干个曲面单元，详见图1，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值。此处划分为个曲面单元，分别为、、三个方向曲面单元的个数。
18.步骤三、依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程，如下：；其中：为哈密顿算子；表示电场强度；表示磁导率；表示角频率；；
表示任意各向异性电导率，用3
×
3的电导率张量表示，，表示轴方向的电导率值；表示轴方向的电导率值；表示轴方向的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；电导率张量是一个对称且正定的矩阵，有。而正定的矩阵都可以通过三次矩阵旋转得到主参考系下的电导率张量，此时仅有对角线元素，其余元素均为零，这样的转换亦便于我们正演模型的生成。即有：；其中： ，代表以及其中一个变量，分别为三次矩阵旋转所对应的旋转角，即各向异性走向角、倾向角、倾斜角。故有六个量可表示电导率张量。
19.步骤四、在双旋度方程中加入散度校正项，得到正则化约束下的正演方程如下：；其中：为正则化因子。
20.通过将散度校正项添加到原始控制方程中来对控制方程进行约束，能有效地避免对散度方程额外的求解（使得迭代求解器的收敛性得到最显著的改善），保证每次迭代电流的散度为零。
21.步骤五、依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；本实施例中，依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件，如下：将地球及其上覆空气层划分为球体层状电导率模型，每一层电场的散度为零，则电场的极向-环向分解表达式如下：
；其中：是东西方向电势，为南北方向电势，为球体半径，为球体半径单位矢量，为梯度算子的角分量，；为拉普拉斯算子的角分量，；将代入双旋度方程中，经过代数变换后，使得电势满足亥姆霍兹方程如下：；其中：代表东西方向电势或南北方向电势；为复波数；对亥姆霍兹方程通过分离变量变成如下表达式：；其中：和表示球谐系数；表示球谐次数；表示球谐阶数；为第一类贝塞尔函数，为第三类贝塞尔函数，与贝塞尔球谐函数有关；是第阶球谐函数的归一化线性组合；空气层通常满足准静态条件，则将贝塞尔函数用常数因子代替转换成下式：；电磁场球谐函数分解成东西方向电势和南北方向电势，将带入电场的极向－环向分解表达式中，消去东西方向电势相关的因式项或消去东西方向电势相关的因式项，得到准静态条件下电场的表达式：；其中：表示施密特归一化缔合勒让德函数；表示阶和阶对应的贝塞尔函数。
22.步骤六、基于交错网格有限差分法对双旋度方程离散，并求解该双旋度方程，得到曲面单元上对应的电场分量。
23.本实施例中，在球坐标系下采用交错网格有限差分法对正则化约束下的正演方程进行离散，其具体形式如下：；
其中：表示离散旋度算子，表示将单元棱边上的矢量映射到单元的曲面上，是的伴随矩阵，表示将单元面上的矢量映射回到单元棱边上；，表示加权因子的对角矩阵，表示离散梯度算子，表示将单元节点上的量映射到单元的内部边上，即，表示单元体的曲边边长，表示离散梯度算子的拓扑矩阵；离散算子则表示将单元内部边上的量映射到内部节点上，通过对转置得到；表示对角矩阵，为单元棱边的电导率，表示离散电场矢量。
24.直到曲面单元上的电场残差分量满足相对残差小于时，停止迭代计算，得到电场分量。
25.步骤七、根据曲面单元上的电场分量求旋度，计算出不同极化模式下曲面单元上的磁场分量如下表达式：。
26.步骤八、根据曲面单元上不同极化模式下的电场分量和磁场分量计算出对应测点的视电阻率和阻抗相位，具体是：根据两种极化模式下的电场分量和磁场分量，计算出两种模式的阻抗：；其中：和为极化计算出方向的电场分量和磁场分量，和为极化计算出方向的电场分量和磁场分量；视电阻率和阻抗相位计算表达式如下：；其中，表示和两种模式下的视电阻率；表示和两种模式下的相位；和表示虚部和实部；表示反三角函数。
27.应用本实施例的方案，具体是：设计了如图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个曲面单元，采取不规则网格剖分。
28.为了说明对比结果。选取了位于剖面中部的一个测点展示数值误差情况，图3和图4为两种极化模式下的视电阻率和阻抗相位对比曲线图以及误差图，计算周期设置为10s到1000s等对数间隔。由图3和图4可知本发明方法（spherical）与通用的各向异性解析解（analytic）在数值模拟结果上基本吻合，均符合绝对误差、相位误差小于5%的要求，从而验证了本发明方法在结果精度上的可靠性。
29.另外，本实施例还公开了一种计算机设备，该计算机设备可以是服务器，其内部结构图可以如图5所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口和数据库。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的数据库用于存储样本数据。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现上述三维大地电磁正演数值模拟方法。
30.本领域技术人员可以理解，图5中示出的结构，仅仅是与本技术方案相关的部分结构的框图，并不构成对本技术方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
31.另外，本实施例还公开了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述实施例中三维大地电磁正演数值模拟方法的步骤。
32.本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本技术所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（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）等。
33.以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
34.以上所述实施例仅表达了本技术的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。

技术特征：
1.一种三维大地电磁正演数值模拟方法，其特征在于，步骤如下：步骤一、依据目标地质体的特征构建任意各向异性电阻率模型；步骤二、基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度将任意各向异性电阻率模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值；步骤三、依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程，如下：；其中：为哈密顿算子；表示电场强度；表示磁导率；表示角频率；；表示任意各向异性电导率，用3
×
3的电导率张量表示，，表示轴方向的电导率值；表示轴方向的电导率值；表示轴方向的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；表示当方向施加电场，在方向会形成电流密度的电导率值；步骤四、在双旋度方程中加入散度校正项，得到正则化约束下的正演方程如下：；其中：为正则化因子；步骤五、依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；步骤六、基于交错网格有限差分法对双旋度方程离散，并求解该双旋度方程，得到曲面单元上对应的电场分量；步骤七、根据曲面单元上的电场分量求旋度，计算出不同极化模式下曲面单元上的磁场分量；步骤八、根据曲面单元上不同极化模式下的电场分量和磁场分量计算出对应测点的视电阻率和阻抗相位。2.根据权利要求1所述的三维大地电磁正演数值模拟方法，其特征在于，步骤五具体如下：将地球及其上覆空气层划分为球体层状电导率模型，每一层电场的散度为零，则电场的极向-环向分解表达式如下：
；其中：是东西方向电势，为南北方向电势；为球体半径， 为球体半径单位矢量，为梯度算子的角分量，；为拉普拉斯算子的角分量，；将代入双旋度方程中，经过代数变换后，使得电势满足亥姆霍兹方程如下：；其中：代表东西方向电势或南北方向电势；为复波数；对亥姆霍兹方程通过分离变量变成如下表达式：；其中：和表示球谐系数；表示球谐次数；表示球谐阶数；为第一类贝塞尔函数，为第三类贝塞尔函数，与贝塞尔球谐函数有关；是第阶球谐函数的归一化线性组合；空气层通常满足准静态条件，则将贝塞尔函数用常数因子代替转换成下式：；电磁场球谐函数分解成东西方向电势和南北方向电势，将带入电场的极向－环向分解表达式中，消去东西方向电势相关的因式项或消去南北方向电势相关的因式项，得到准静态条件下电场的表达式，即大地电磁法边界电场值：；其中：表示施密特归一化缔合勒让德函数；表示阶和阶对应的贝塞尔函数。3.根据权利要求1-2任意一项所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤六具体是：在球坐标系下采用交错网格有限差分法对正则化约束下的正演方程进行离散，其具体形式如下：
；其中：表示离散旋度算子，表示将单元棱边上的矢量映射到单元的曲面上；是的转置矩阵，表示将单元面上的矢量映射回到单元棱边上；；表示加权因子的对角矩阵；表示离散梯度算子，表示将单元节点上的量映射到单元的内部边上，即，表示单元体的曲边边长，表示离散梯度算子的拓扑矩阵；离散算子则表示将单元内部边上的量映射到内部节点上，通过对转置得到；表示对角矩阵，为单元棱边的电导率，表示离散电场矢量；直到曲面单元上的电场残差分量满足相对残差小于时，停止迭代计算，得到电场分量。4.根据权利要求3所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤一具体是：所述目标地质体的特征包括目标地质体的形状、大小和电阻率分布；所述目标地质体为三维异常体。5.根据权利要求3所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤七中磁场分量如下表达式：。6.根据权利要求5所述的三维大地电磁正演数值模拟方法，其特征在于，所述步骤八中：根据两种极化模式下的电场分量和磁场分量，计算出两种模式的阻抗：；其中：和为极化计算出方向的电场分量和磁场分量，和为极化计算出方向的电场分量和磁场分量；视电阻率和阻抗相位计算表达式如下：；其中，表示和两种模式下的视电阻率；表示和两种模式下的相位；和表示虚部和实部；表示反三角函数。

技术总结
本发明提出一种三维大地电磁正演数值模拟方法，具体是：依据目标地质体的特征构建任意各向异性电阻率模型；基于球坐标系对任意各向异性电阻率模型进行网格剖分，沿深度、纬度和经度将地质体模型以及上覆空气层剖分成若干个曲面单元，并根据任意各向异性电阻率分布给每个曲面单元电阻率赋值，依据频率参数和任意各向异性电阻率模型，构建曲面单元上电场所满足的双旋度方程；依据球谐函数近似大地电磁场源，实现球谐函数场源下的大地电磁法边界条件；采用正则化约束下的正演方程获得两种极化模式下的电场分量和磁场分量，继而得到视电阻率和阻抗相位。本发明既能够真实模拟三维大地电磁情况，又能保证高精度求解的同时提高电磁正演的计算效率。正演的计算效率。正演的计算效率。


技术研发人员：侯生 郭荣文 李健 柳建新
受保护的技术使用者：中南大学
技术研发日：2023.09.04
技术公布日：2023/10/8