CN115906559B - 一种基于混合网格的大地电磁自适应有限元正演方法 - Google Patents

Abstract

本发明公开一种基于混合网格的大地电磁自适应有限元正演方法，包括以下步骤：1)获取待计算的大地电磁场区域，并将所述大地电磁场区域划分为近地表区域、过渡区域和剩余区域；2)利用混合网格对近地表区域、过渡区域和剩余区域进行剖分，得到混合网格模型；3)对所述混合网格模型进行优化；4)基于优化后的混合网格模型，完成大地电磁正演计算。本发明可有效提高高频下大地电磁三维正演计算的效率，对三维地下电性结构的大地电磁的响应特征的研究，以及对大地电磁的观测数据的分析和解释有良好的实用价值和应用前景。

Description

一种基于混合网格的大地电磁自适应有限元正演方法

技术领域

本发明涉及地球物理电磁计算领域，具体是一种基于混合网格的大地电磁自适应有限元正演方法。

背景技术

大地电磁测深法是以天然的平面电磁场作为场源，根据电磁场在地下的分布情况随频率变化的原理，利用专门的仪器在地表测量多个频点下的电磁信号，再对这些电磁信号进行处理后，反推出地下电性结构的电阻率和形态。大地电磁具有施工方便、勘探成本低、勘探深度大、不受高阻屏蔽和对低阻分辨率高等优点，在油气勘探、地热资源、地壳和上地幔电性结构研究等领域应用广泛。近几十年来，计算机领域发展的十分迅速，这促进了学者们对大地电磁三维有限元正反演的研究。

目前，基于结构化网格的大地电磁三维有限元正反演已经较为成熟，但这种方法不容易模拟复杂的地质体以及计算效率不高。所以许多学者采用非结构化网格离散计算区域，并在这种网格的基础上进行了大地电磁三维正反演研究，这很好的解决了复杂地质体不易模拟的问题，但依然存在两个问题。第一个问题是在保证计算结果精度的前提下，需要大量的网格单元剖分计算区域，这会产生大量的网格单元数，从而导致计算效率不高。第二个问题是人工剖分网格存在缺陷。大地电磁的频率范围较宽，可从10^-4Hz到10⁴Hz。在不同频率下，地下电磁场的分布也大不相同，很难找到一套网格同时在多个频点下获得高精度的计算结果，并且网格分布最优。针对上述的问题，已经有许多学者们在大地电磁三维有限元计算的过程中采用网格自适应技术来找到不同频点下所需的最优网格，这种方法减少了网格数量和自由度数，计算结果的精度和计算效率都有所提高。

虽然利用网格自适应技术可以在三维大地电磁自适应正演计算中提高结果精度和计算效率，但是在自适应为保证高频下的结果精度，依然会在近地表区域产生大量的四面体单元，这降低了计算效率。

发明内容

本发明的目的是提供一种基于混合网格的大地电磁自适应有限元正演方法，包括以下步骤：

1)获取待计算的大地电磁场区域，并将所述大地电磁场区域划分为近地表区域、过渡区域和剩余区域；

2)利用混合网格对近地表区域、过渡区域和剩余区域进行剖分，得到混合网格模型；3)对所述混合网格模型进行优化；

4)基于优化后的混合网格模型，完成大地电磁正演响应的有限元数值计算。

进一步，所述过渡区域与近地表区域具有相同界面，该界面记为第一界面；

所述过渡区域与剩余区域具有相同界面，该界面记为第二界面。

进一步，所述近地表区域的深度为最高频率下对应趋肤深度δ的3倍，横向尺寸覆盖住整个测点；

所述过渡区域的厚度是近地表区域深度的1～2倍，横向范围包围近地表区域范围，过渡区域水平尺寸大于近地表区域中侧边最长的三棱柱单元侧边的1.62倍；

所述剩余区域为整个计算区域中，去除近地表区域和过渡区域剩余的区域。

进一步，利用混合网格对近地表区域、过渡区域和剩余区域进行剖分的步骤包括：

利用三棱柱网格对近地表区域进行剖分；

利用四面体单元和金字塔单元同时对过渡区域进行剖分；所述四面体单元包括四个三角形面；所述金字塔单元包括一个四边形面和四个三角形面；

利用四面体网格对剩余区域进行剖分。

进一步，利用四面体单元和金字塔单元同时对过渡区域进行剖分的步骤包括：

I)将过渡区域划分为边界区域和内部区域；

II)根据近地表区域的剖分方式将边界区域划分为三棱柱单元底边界区域和三棱柱单元侧边界区域；其中，三棱柱单元底边界区域包括若干三棱柱单元底面；所述三棱柱单元侧边界区域包括若干三棱柱单元侧面；

III)利用四面体单元对三棱柱单元底边界区域、内部区域进行剖分；利用金字塔单元对三棱柱单元侧边界区域进行剖分。

进一步，所述金字塔单元用于实现三棱柱网格的四边形面与四面体单元的三角形面的耦合连接。

进一步，对近地表区域、过渡区域和剩余区域进行剖分的工具包括COMSOL中内置的网格剖分工具；

对所述混合网格模型进行优化的工具包括COMSOL中内置的网格细化工具。

进一步，对所述混合网格模型进行优化的方法包括自适应网格细化方法。

进一步，对所述混合网格模型进行优化的方法包括八分法，即将每个网格划分为八个子网格。

进一步，在混合网格更新的过程中，根据计算得到的网格单元误差利用COMSOL软件对指定的网格单元进行细化，步骤包括：

a)将三棱柱单元和金字塔单元设置为零，设置最小体积约束；

b)计算第i个四面体单元的误差η_i，即：

η_i＝η_w,i·η_e,i (1)

其中，特定物理量的后验误差估计量η_e,i、受观测点解精度影响的误差估计器的权重估计η_w,i分别如下所示：

式中，w₊和w_-是内、外影响函数；ε_k+和ε_k-分别是沿四面体单元法向n的内部和外部误差；F_k表示四面体单元的第k个表面；是表面F_k的切向；/>是表面F_k的法向；

其中，内影响函数w₊、外影响函数w_-通过求解EM势方程的对偶形式得到；

EM势方程的对偶形式如下所示：

B^*(w,Φ)＝L(Φ) (4)

式中，V_j为电压；I为电流；Φ为矢量基函数；L(Φ)为平方可积函数；_Ωj为第j个单元积分区域；B*(,)是EM势方程系统B(,)的对偶形式；t为单元总数；v为积分变量；

c)令i＝i+1，并重复步骤b)，直到获取所有四面体单元的误差；

d)按照误差从大到小，选择出m％M个四面体单元，然后对这些单元中体积大于最小体积约束的单元进行细化。m>0，M为四面体单元总数。

值得说明的是，本发明具体涉及利用三棱柱、金字塔和四面体这三种网格剖分计算区域，并用这种混合网格进行了大地电磁三维自适应正演，实现了在保证结果精度的前提下减少网格数量和自由度数，提高了计算效率。

本发明的技术效果是毋庸置疑的，本专利在近地表区域采用三棱柱单元剖分，在其余区域采用四面体单元剖分，并在这混合网格的基础上进行了大地电磁三维自适应正演的研究。

与现有的基于纯四面体网格的大地电磁三维自适应有限元计算技术相比，当自由度数相近时，混合网格比四面体网格更容易获得高精度数值解。特别是对于承载近地表区域的高频数据，混合网格的初始网格无需自适应优化即可获得高精度数值解，所需的自由度数仅为传统四面体自适应网格的几十分之一。对于低频数据，该方法可以保持与传统自适应方法相近的效果。

本发明方法可有效提高高频下大地电磁三维正演计算的效率，对三维地下电性结构的大地电磁的响应特征的研究，以及对大地电磁的观测数据的分析和解释有良好的实用价值和应用前景。

附图说明

图1为混合网格自适应剖分示意图；

图2为四面体单元细化示意图；

图3为三棱柱单元构建示意图；

图4为为五层层状模型示意图；

图5为不同网格在z＝0m处的网格示意图；图5(a)为混合网格在z＝0m处的网格示意图；图5(b)为纯四面体网格在z＝0m处的网格示意图；

图6为网格耦合示意图；图6(a)为三棱柱与四面体单元的耦合示意图；图6(b)为三棱柱与金字塔单元的耦合示意图；图6(c)为金字塔与三棱柱单元的耦合示意图；

图7为八分单元的相邻单元的细化示意图；图7(a)-(b)分别是网格四分和二分细化的示意图，其四分面与八分单元共面，二分边与八分单元共边；

图8为1000Hz下自适应响应误差迭代图；图8(a)为1000Hz下自适应响应视电阻率绝对误差迭代图；图8(b)为1000Hz下自适应响应相位绝对误差迭代图；

图9为1Hz下自适应响应误差迭代图；图9(a)为1Hz下自适应响应视电阻率绝对误差迭代图；图9(b)为1Hz下自适应响应相位绝对误差迭代图；

图10为0.04Hz下自适应响应误差迭代图；图10(a)为0.04Hz下自适应响应视电阻率绝对误差迭代图；图10(b)为0.04Hz下自适应响应相位绝对误差迭代图。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

实施例1：

参见图1至图10，一种基于混合网格的大地电磁自适应有限元正演方法，包括以下步骤：

1)获取待计算的大地电磁场区域，并将所述大地电磁场区域划分为近地表区域、过渡区域和剩余区域；

2)利用混合网格对近地表区域、过渡区域和剩余区域进行剖分，得到混合网格模型；

3)对所述混合网格模型进行优化；

4)基于优化后的混合网格模型，完成大地电磁正演响应的有限元数值计算。具体计算参考顾观文,武晔,石砚斌.基于矢量有限元的大地电磁快速三维正演研究[J].物探与化探,2020,44(06):1387-1398。

所述过渡区域与近地表区域具有相同界面，该界面记为第一界面；

所述过渡区域与剩余区域具有相同界面，该界面记为第二界面。

所述近地表区域的深度为最高频率下对应趋肤深度δ的3倍，横向尺寸覆盖住整个测点；趋肤深度是指电磁波的幅值衰减为表面值的0.368时，电磁波所传播的距离。

所述过渡区域的厚度是近地表区域深度的1～2倍，横向范围包围近地表区域范围，过渡区域水平尺寸大于近地表区域中侧边最长的三棱柱单元侧边的1.62倍；

所述剩余区域为整个计算区域中，去除近地表区域和过渡区域剩余的区域。

利用混合网格对近地表区域、过渡区域和剩余区域进行剖分的步骤包括：

利用三棱柱网格对近地表区域进行剖分；

利用四面体单元和金字塔单元同时对过渡区域进行剖分；所述四面体单元包括四个三角形面；所述金字塔单元包括一个四边形面和四个三角形面；

利用四面体网格对剩余区域进行剖分。

利用四面体单元和金字塔单元同时对过渡区域进行剖分的步骤包括：

I)将过渡区域划分为边界区域和内部区域；

II)根据近地表区域的剖分方式将边界区域划分为三棱柱单元底边界区域和三棱柱单元侧边界区域；其中，三棱柱单元底边界区域包括若干三棱柱单元底面；所述三棱柱单元侧边界区域包括若干三棱柱单元侧面；

III)利用四面体单元对三棱柱单元底边界区域、内部区域进行剖分；利用金字塔单元对三棱柱单元侧边界区域进行剖分。

所述金字塔单元用于实现三棱柱网格的四边形面与四面体单元的三角形面的耦合连接。

对近地表区域、过渡区域和剩余区域进行剖分的工具包括COMSOL中内置的网格剖分工具；

对所述混合网格模型进行优化的工具包括COMSOL中内置的网格细化工具。

对所述混合网格模型进行优化的方法包括自适应网格细化方法。

对所述混合网格模型进行优化的方法包括八分法，即将每个网格划分为八个子网格。

在混合网格更新的过程中，根据计算得到的网格单元误差利用COMSOL软件对指定的网格单元进行细化，步骤包括：

a)将三棱柱单元和金字塔单元设置为零，设置最小体积约束；

b)计算第i个四面体单元的误差η_i，即：

η_i＝η_w,i·η_e,i (1)

其中，特定物理量的后验误差估计量η_e,i、受观测点解精度影响的误差估计器的权重估计η_w,i分别如下所示：

式中，w₊和w_-是内、外影响函数，具体的，分界面两侧的单元用正负区分，单元内部影响函数用w+，单元外部(相邻单元)的影响函数用w-。；ε_k+和ε_k-分别是沿四面体单元法向n的内部和外部误差；F_k表示四面体单元的第k个表面；t_Fk是表面F_k的切向；n_Fk是表面F_k的法向；

其中，内影响函数w₊、外影响函数w_-通过求解EM势方程的对偶形式得到；

EM势方程的对偶形式如下所示：

B^*(w,Φ)＝L(Φ) (4)

式中，V_j为电压；I为电流；Φ为矢量基函数；L(Φ)为平方可积函数；_Ωj为第j个单元积分区域；B*(,)是EM势方程系统B(,)的对偶形式；t为单元总数；v为积分变量；

c)令i＝i+1，并重复步骤b)，直到获取所有四面体单元的误差；

d)按照误差从大到小，选择出m％M个四面体单元，然后对这些单元中体积大于最小体积约束的单元进行细化。m>0，M为四面体单元总数。

实施例2：

一种基于混合网格的大地电磁自适应有限元正演方法，包括以下步骤：

1)获取待计算的大地电磁场区域，并将所述大地电磁场区域划分为近地表区域、过渡区域和剩余区域；

2)利用混合网格对近地表区域、过渡区域和剩余区域进行剖分，得到混合网格模型；3)对所述混合网格模型进行优化；

4)基于优化后的混合网格模型，完成大地电磁正演响应的有限元数值计算。

所述过渡区域与近地表区域具有相同界面，该界面记为第一界面；

所述过渡区域与剩余区域具有相同界面，该界面记为第二界面。

利用混合网格对近地表区域、过渡区域和剩余区域进行剖分的步骤包括：

利用三棱柱网格对近地表区域进行剖分；

三棱柱单元的纵向延伸长度取决于趋肤深度δ，当三棱柱单元的纵向延伸长度小于趋肤深度时，单元内部的电磁场不接近线性，这会使得计算结果不准确。为保证计算结果的准确性，我们通常使棱柱单元网格的厚度不小于趋肤深度。在本发明中设置棱柱单元的厚度(即近地表区域的深度)为最高频率下对应趋肤深度的3倍。横向尺寸覆盖住整个测点部分。

利用四面体单元和金字塔单元同时对过渡区域进行剖分；所述四面体单元包括四个三角形面；所述金字塔单元包括一个四边形面和四个三角形面；

过渡区域分布在近地表区域的侧部和底部。其中底部的过渡区域由四面体单元进行剖分，即近地表区域三棱柱的三角形底面与过渡区域四面体的三角形面进行耦合。侧部的过渡区域由金字塔单元进行剖分，即近地表区域三棱柱的四边形侧面与过渡区域的金字塔的四边形面进行耦合。在本发明中设置过渡区域的厚度为1～2倍的近地表区域的深度，横向范围包围近地表区域范围，根据金字塔单元的黄金比例，金字塔单元底边与塔高比例为0.618，所以设置过渡区域水平尺寸需大于近地表区域中侧边最长的三棱柱单元侧边的1.62倍。

利用四面体网格对剩余区域进行剖分。

模型的剩余区域全部采用四面体单元进行剖分，四面体单元的四个三角形面正好可以耦合过渡区域底部四面体单元的三角形面和过渡区域侧面金字塔单元的三角形面。在本发明中剩余区域的范围没有明确的界定，取决于近地表区域和过渡区域的范围。整个计算区域去除近地表区域和过渡区域剩余的即为近地表区域。

利用四面体单元和金字塔单元同时对过渡区域进行剖分的步骤包括：

1)将过渡区域划分为第一界面区域(边界区域)和其他区域(内部区域)；

2)根据近地表区域的剖分方式将第一界面区域(边界区域)划分为三棱柱单元底面区域(底边界区域)和三棱柱单元侧面区域(侧边界区域)；其中，三棱柱单元底面区域包括若干三棱柱单元底面；所述三棱柱单元侧面区域包括若干三棱柱单元侧面；

3)利用四面体单元对三棱柱单元底面区域(底边界区域)、其他区域(内部区域)进行剖分；利用金字塔单元对三棱柱单元侧面区域(侧边界区域)进行剖分；

所述金字塔单元用于实现三棱柱网格的四边形面与四面体单元的三角形面的耦合连接。

对近地表区域、过渡区域和剩余区域进行剖分的工具包括COMSOL中内置的网格剖分工具；

对所述混合网格模型进行优化的工具包括COMSOL中内置的网格细化工具。

对所述混合网格模型进行优化的方法包括自适应网格细化方法。

对所述混合网格模型进行优化的方法包括八分法，即将每个网格划分为八个子网格。

在混合网格更新的过程中，根据计算得到的网格单元误差利用COMSOL软件对指定的网格单元进行细化。

由于混合网格中的三棱柱单元和金字塔单元不进行自适应细化，因此在计算混合网格单元误差时，只计算四面体单元的误差。在自适应过程中将三棱柱单元和金字塔单元设置为零。计算网格误差时，设置自适应的最大迭代次数，最大单元数，最小体积约束和误差阈值，从误差大的单元开始选择，选择出单元总数的1％个单元，然后指定对这些单元中体积大于最小体积约束的单元进行细化。当满足自适应的停止迭代要求时则退出迭代过程。

其中，网格单元误差的计算公式如下：

将网格自适应中第i个单元的总体后验误差估计η_i设置为面向目标的形式：

η_i＝η_w,i·η_e,i

式中，η_e,i是特定物理量的后验误差估计量，它可以指示物理量的正常不连续性或切线不连续性，例如元素界面的电流密度和磁电流密度。η_w,i是受观测点解精度影响的误差估计器的权重估计，用于在较小范围内定位目标区域。

对于四面体，后验误差估计器可以求四个面上的误差之和，可以用以下形式表示：

其中，在四面体单元的第k个表面F_k上，ε_k+和ε_k-分别是沿四面体单元法向n的内部和外部误差，误差的大小由L2范数测量。误差ε可以是电流密度J不连续性和磁电流密度B不连续性。

当必须测量物理量的切向不连续性时，假设t是元素表面F_k的切向，则后验误差估计值可以表示为：

加权后验误差估计器的权重η_w,i可以表示为：

其中w₊和w_-是元素内外的影响函数。EM势方程可以表示为 其中Φ代表矢量基函数，L是平方可积函数，可以表示为：

通过求解EM势方程的对偶弱形式可以得到影响函数，对偶形式如下：

B^*(w,Φ)＝L(Φ)

其中B*(,)是EM势方程系统的B(,)的对偶形式。当Φ可以作为电磁势有限元方程中的基函数时，由于自混合性质，B*(,)等于B(,)。

实施例3：

一种基于混合网格的大地电磁自适应有限元正演方法，步骤如下：

(1)混合网格的剖分

层状模型如图4所示。第一层到第四层的层厚均为4km，第一层到第五层的电阻率分别为100Ω·m、10Ω·m、50Ω·m、10Ω·m和100Ω·m。计算区域大小为60km×60km×200km。这里选择1000Hz、1Hz和0.04Hz三个频率点进行计算。

以传统的基于四面体网格的自适应技术作为参考项，评价在自适应算法中使用混合网格的计算效果。图5是混合网格和纯四面体网格初始网格在地面处的网格示意图，两种网格的初始网格的自由度数接近，分别为74154和73690。这里混合网格中的三棱柱单元水平跨度为600m，总厚度为500m，划分为8层，首层厚度约为16m，三棱柱单元厚度的增长率约为1.37。纯四面体网格水平跨度为200m，单元增长率为1.2。

(2)混合网格的不同单元之间的耦合

图6是不同网格之间耦合的示意图，其中的阴影面代表不同网格的耦合面。图6a、b和c分别是三棱柱单元和四面体单元的具体耦合示意图、三棱柱单元和金字塔单元的具体耦合示意图以及金字塔单元和四面体单元的具体耦合示意图。根据这些网格耦合示意图可知本专利利用金字塔单元帮助三棱柱单元和四面体单元的具体耦合方式为：首先，金字塔单元的矩形与三棱柱单元的侧面直接耦合，然后四面体单元与金字塔单元的三角形面进行耦合。

(3)自适应网格细化

本专利采用的是常规的有限元算法，并不支持网格出现悬挂节点。对需误差大的网格单元八分后，为了避免八分单元和与八分单元相邻的单元之间出现悬挂节点，这些相邻的单元也需要细化。八分单元的相邻单元可分为共单点单元(后续提到的共点单元指共单点单元)、共边单元(共用两个点)以及共面单元(共用三个点)。其中，共点单元不需要随着八分单元细化，而共边单元和共面单元需要随着八分单元细化。共面单元与八分单元的共面处被分为了四个三角形(图7a)，所以共面单元随着八分单元四分。这里的共面单元四分后，共面单元的相邻单元也需要细分，但是随着共面单元细化的相邻单元刚好是八分单元的共边单元，因此共面单元的相邻单元的细化情况此处不再额外介绍。共边单元与八分单元的共用边被等分为两段(图7b)，所以共边单元随着八分单元二分。

(4)混合网格的构建与更新

三棱柱网格在COMSOL中分为两步构建。第一步先利用COMSOL中的“自由三角形网格”工具将地表面用自由三角形单元离散。第二步利用“扫掠”工具从地表开始，向地下扫掠，从而构建出三棱柱单元。

四面体单元可以利用COMSOL中的“自由四面体网格”和“大小”工具构建。其中“大小”工具用于控制四面体单元的最大单元大小、最大单元增长率以及最小单元大小。最大单元增长率是指四面体单元的增长率不超过设定值。金字塔单元只出现在四面体单元与三棱柱单元耦合处，由“自由四面体网格”网格剖分工具插入。

在自适应过程中，利用COMSOL中内置的网格细化工具根据网格误差更新网格。每次选择网格单元数的1％且误差大的单元进行细化。

实验效果

图8是1000Hz下自适应过程中的响应误差图。从图中可以看到初始混合网格的视电阻率和相位的绝对误差分别小于0.01％和0.1％。四面体网格经过5次自适应迭代后相位的绝对误差小于0.1％，经过7次自适应迭代后视电阻率的绝对误差与混合网格接近，小于0.02％。然而，此时四面体网格的DoFs为2201100，约为混合网格初始网格自由度的29倍。1000Hz的电磁响应信号主要集中在浅层近地表。初始的混合网格无需进行自适应即可以获得高精度的数值解，其自由度数远小于四面体网格生成的自由度数。这表明，在高频情况下，三棱柱单元大大减少了获得高精度数值解所需的自由度数，提高了计算效率。

随着频率的降低，感应电磁场携带的地层信息由浅到深变化。当频率为1Hz时，混合网格在自由度数方面相对于传统四面体网格的优势不能直接体现在初始混合网格中(图9)，需要在自适应优化的迭代过程中体现出来。在相位误差随自适应迭代变化的曲线中，当自由度数相近时，混合网格的误差小于纯四面体网格的误差。在自适应迭代的视电阻率误差曲线中，混合网格与纯四面体网格的误差下降曲线相近。

当频率为0.04Hz时，7次自适应迭代后两种网格的视电阻率和相位的绝对误差都分别小于0.0004％和0.0003％。虽然混合网格的相位误差曲线下降速度比纯四面体网格快(图10a)，但是两种网格的视电阻率误差曲线下降速度接近(图10b)。这表明混合网格对计算效率的提升随着频率的降低而减弱。

总的来说，高频下混合网格能在保持高精度结果的前提下有效提高计算效率，低频下能保持与纯四面体网格相近的计算效率。

实施例4：

一种基于混合网格的大地电磁自适应有限元正演方法，包括以下步骤：

1)获取待计算的大地电磁场区域，并将所述大地电磁场区域划分为近地表区域、过渡区域和剩余区域；

2)利用混合网格对近地表区域、过渡区域和剩余区域进行剖分，得到混合网格模型；

3)对所述混合网格模型进行优化；

4)基于优化后的混合网格模型，完成大地电磁正演响应的有限元数值计算。

实施例5：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，所述过渡区域与近地表区域具有相同界面，该界面记为第一界面；

所述过渡区域与剩余区域具有相同界面，该界面记为第二界面。

实施例6：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，所述近地表区域的深度为最高频率下对应趋肤深度δ的3倍，横向尺寸覆盖住整个测点；

所述过渡区域的厚度是近地表区域深度的1～2倍，横向范围包围近地表区域范围，过渡区域水平尺寸大于近地表区域中侧边最长的三棱柱单元侧边的1.62倍；

所述剩余区域为整个计算区域中，去除近地表区域和过渡区域剩余的区域。

实施例7：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，利用混合网格对近地表区域、过渡区域和剩余区域进行剖分的步骤包括：

利用三棱柱网格对近地表区域进行剖分；

利用四面体单元和金字塔单元同时对过渡区域进行剖分；所述四面体单元包括四个三角形面；所述金字塔单元包括一个四边形面和四个三角形面；

利用四面体网格对剩余区域进行剖分。

实施例8：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，利用四面体单元和金字塔单元同时对过渡区域进行剖分的步骤包括：

1)将过渡区域划分为边界区域和内部区域；

2)根据近地表区域的剖分方式将边界区域划分为三棱柱单元底边界区域和三棱柱单元侧边界区域；其中，三棱柱单元底边界区域包括若干三棱柱单元底面；所述三棱柱单元侧边界区域包括若干三棱柱单元侧面；

3)利用四面体单元对三棱柱单元底边界区域、内部区域进行剖分；利用金字塔单元对三棱柱单元侧边界区域进行剖分。

实施例9：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，所述金字塔单元用于实现三棱柱网格的四边形面与四面体单元的三角形面的耦合连接。

实施例10：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，对近地表区域、过渡区域和剩余区域进行剖分的工具包括COMSOL中内置的网格剖分工具；

对所述混合网格模型进行优化的工具包括COMSOL中内置的网格细化工具。

实施例11：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，对所述混合网格模型进行优化的方法包括自适应网格细化方法。

实施例12：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，对所述混合网格模型进行优化的方法包括八分法，即将每个网格划分为八个子网格。

实施例13：

一种基于混合网格的大地电磁自适应有限元正演方法，主要内容见实施例4，其中，在混合网格更新的过程中，根据计算得到的网格单元误差利用COMSOL软件对指定的网格单元进行细化，步骤包括：

1)将三棱柱单元和金字塔单元设置为零，设置最小体积约束；

2)计算第i个四面体单元的误差η_i，即：

η_i＝η_w,i·η_e,i (1)

其中，特定物理量的后验误差估计量η_e,i、受观测点解精度影响的误差估计器的权重估计η_w,i分别如下所示：

式中，w₊和w_-是内、外影响函数；ε_k+和ε_k-分别是沿四面体单元法向n的内部和外部误差；F_k表示四面体单元的第k个表面；是表面F_k的切向；/>是表面F_k的法向；

其中，内影响函数w₊、外影响函数w_-通过求解EM势方程的对偶形式得到；

EM势方程的对偶形式如下所示：

B^*(w,Φ)＝L(Φ) (4)

式中，V_j为电压；I为电流；Φ为矢量基函数；L(Φ)为平方可积函数；_Ωj为第j个单元积分区域；B*(,)是EM势方程系统B(,)的对偶形式；t为单元总数；v为积分变量；

3)令i＝i+1，并重复步骤2)，直到获取所有四面体单元的误差；

4)按照误差从大到小，选择出m％M个四面体单元，然后对这些单元中体积大于最小体积约束的单元进行细化。常数m>0，M为四面体单元总数。

Claims (7)

1.一种基于混合网格的大地电磁自适应有限元正演方法，其特征在于，包括以下步骤：

1)获取待计算的大地电磁场区域，并将所述大地电磁场区域划分为近地表区域、过渡区域和剩余区域；

2)利用混合网格对近地表区域、过渡区域和剩余区域进行剖分，得到混合网格模型；

3)对所述混合网格模型进行优化；

4)基于优化后的混合网格模型，完成大地电磁正演响应的有限元数值计算；

过渡区域分布在近地表区域的侧部和底部；

整个计算区域去除近地表区域和过渡区域剩余的即为近地表区域；

利用混合网格对近地表区域、过渡区域和剩余区域进行剖分的步骤包括：

利用三棱柱网格对近地表区域进行剖分；

利用四面体单元和金字塔单元同时对过渡区域进行剖分；所述四面体单元包括四个三角形面；所述金字塔单元包括一个四边形面和四个三角形面；

利用四面体网格对剩余区域进行剖分；

利用四面体单元和金字塔单元同时对过渡区域进行剖分的步骤包括：

a)将过渡区域划分为边界区域和内部区域；

b)根据近地表区域的剖分方式将边界区域划分为三棱柱单元底边界区域和三棱柱单元侧边界区域；其中，三棱柱单元底边界区域包括若干三棱柱单元底面；所述三棱柱单元侧边界区域包括若干三棱柱单元侧面；

c)利用四面体单元对三棱柱单元底边界区域、内部区域进行剖分；利用金字塔单元对三棱柱单元侧边界区域进行剖分；

所述金字塔单元用于实现三棱柱网格的四边形面与四面体单元的三角形面的耦合连接。

2.根据权利要求1所述的一种基于混合网格的大地电磁自适应有限元正演方法，其特征在于，所述过渡区域与近地表区域具有相同界面，该界面记为第一界面；

所述过渡区域与剩余区域具有相同界面，该界面记为第二界面。

3.根据权利要求1所述的一种基于混合网格的大地电磁自适应有限元正演方法，其特征在于，所述近地表区域的深度为最高频率下对应趋肤深度δ的3倍，横向尺寸覆盖住整个测点；

所述过渡区域的厚度是近地表区域深度的1～2倍，横向范围包围近地表区域范围，过渡区域水平尺寸大于近地表区域中侧边最长的三棱柱单元侧边的1.62倍；

所述剩余区域为整个计算区域中，去除近地表区域和过渡区域剩余的区域。

4.根据权利要求1所述的一种基于混合网格的大地电磁自适应有限元正演方法，其特征在于，对近地表区域、过渡区域和剩余区域进行剖分的工具包括COMSOL中内置的网格剖分工具；

对所述混合网格模型进行优化的工具包括COMSOL中内置的网格细化工具。

5.根据权利要求1所述的一种基于混合网格的大地电磁自适应有限元正演方法，其特征在于，对所述混合网格模型进行优化的方法包括自适应网格细化方法。

6.根据权利要求1所述的一种基于混合网格的大地电磁自适应有限元正演方法，其特征在于，对所述混合网格模型进行优化的方法包括八分法，即将每个网格划分为八个子网格。

7.根据权利要求1所述的一种基于混合网格的大地电磁自适应有限元正演方法，其特征在于，在混合网格更新的过程中，根据计算得到的网格单元误差利用COMSOL软件对指定的网格单元进行细化，步骤包括：

1)将三棱柱单元和金字塔单元设置为零，设置最小体积约束；

2)计算第i个四面体单元的误差η_i，即：

η_i＝η_w,i·η_e,i (1)

其中，特定物理量的后验误差估计量η_e,i、受观测点解精度影响的误差估计器的权重估计η_w,i分别如下所示：

式中，w₊和w_-是内、外影响函数；ε_k+和ε_k-分别是沿四面体单元法向n的内部和外部误差；F_k表示四面体单元的第k个表面；是表面F_k的切向；/>是表面F_k的法向；

其中，内影响函数w₊、外影响函数w_-通过求解EM势方程的对偶形式得到；

EM势方程的对偶形式如下所示：

式中，V_j为电压；I为电流；Φ为矢量基函数；L(Φ)为平方可积函数；_Ωj为第j个单元积分区域；B*(,)是EM势方程系统B(,)的对偶形式；t为单元总数；v为积分变量；

3)令i＝i+1，并重复步骤2)，直到获取所有四面体单元的误差；

4)按照误差从大到小，选择出m％M个四面体单元，然后对这些单元中体积大于最小体积约束的单元进行细化；常数m>0，M为四面体单元总数。