基于多尺度有限元的复杂模型航空电磁三维快速正演方法

技术领域

本发明涉及地球物理电磁计算技术领域，具体涉及一种基于多尺度有限元的复杂模型航空电磁三维快速正演方法。

背景技术

目前在三维航空电磁领域，实现三维正演模拟的方法主要包括积分方程法、有限差分法与有限元法等。相较于其他正演模拟方法，有限元法在复杂介质的三维电磁数据模拟中具有较大优势。国内外众多学者对其进行了研究，Pridmore(1978，1979，1981)等人、Unsworth(1993)等人、李长伟等人(2010)和蔡红柱等人(2015)均成功将有限元法用于三维电磁数据的模拟。

正演是反演核心技术部分，正演算法中，有限差分算法比较简单、便于编程，在当前三维正演中使用最多，其次是积分方程法和有限体积法，最后是有限元法，虽然目前有限元法已经被广泛用于各方法正演计算，但由于航空电磁系统响应模拟和灵敏度矩阵计算必须求解超大规模线性方程组，计算效率较低，因此无法应用于大型模型航空电磁数据三维精细反演。

发明内容

为了解决现有技术中心计算效率较低，无法应用于大型模型航空电磁数据三维精细反演的问题，本发明提供一种基于多尺度有限元的复杂模型航空电磁三维快速正演方法，该方案包括：设计航空电磁三维正演模拟的复杂模型，将计算区域剖分为六面体网格得到三维正演初始模型；将每个网格单元作为粗网格，进行剖分得到多个细网格；获取三维正演初始模型中的局部网格；计算局部网格中细网格的质量矩阵、刚度矩阵和源项，并合成整体矩阵与右端向量；构建多尺度有限元插值算子对整体矩阵和右端向量进行降维处理；构建多尺度有限元线性方程组进行求解，得到粗网格棱边处电场；根据粗网格棱边处电场获取接收点处磁场值。本发明降低了正演方程的阶数，进而提高了航空电磁三维数值模拟的计算效率。

本发明采用如下技术方案：基于多尺度有限元的复杂模型航空电磁三维快速正演方法，包括：

设计航空电磁三维正演模拟的复杂模型，根据所述复杂模型将计算区域剖分为六面体网格，得到三维正演初始模型；

将三维正演初始模型中的每个网格单元作为粗网格，对每个网格单元进行剖分，得到每个网格单元对应的多个细网格；

利用模型网络提取方法获取三维正演初始模型中的核心网格，对核心网格进行延展得到局部网格；所述局部网格包含粗网格和细网格；

利用坐标变换计算局部网格中每个细网格的质量矩阵、刚度矩阵，并结合高斯积分和局部网格中细网格的位置关系合成有限元方法中细网格的整体矩阵；

根据麦克斯韦方程组构建多尺度有限元插值算子；利用所述插值算子对所述细网格的整体矩阵进行降维处理，得到细网格所在粗网格的整体矩阵；

根据粗网格的整体矩阵构建多尺度有限元线性方程组，求解所述多尺度有限元线性方程组，得到粗网格棱边处电场；

根据粗网格棱边处电场获取细网格棱边处的电场，并根据细网格棱边处电场获取接收点处磁场值。

进一步的，所述复杂模型包括：起伏地表、异常地质体、收发装置类型、发射频率、飞行高度以及网格之间的对应关系；

进一步的，对核心网格进行延展得到局部网格的方法为：

采用狄利克雷边界条件对网格进行扩边处理，以收发装置为中心进行延展生成局部网格。

进一步的，利用坐标变换计算局部网格中每个细网格的质量矩阵、刚度矩阵，表达式为：

其中，

进一步的，根据麦克斯韦方程组构建多尺度有限元插值算子的方法为：

将插值算子带入到电场双旋度方程中，令双旋度方程右端量为零，令粗网格边界处插值算子值等于粗网格相应位置矢量基函数值，求解该方程得到多尺度有限元插值算子：

其中，v

进一步的，利用坐标变换与高斯积分计算局部网格中每个细网格的质量矩阵、刚度矩阵和源项时，频率域电场满足以下双旋度方程：

其中，E

利用伽辽金法与第一格林公式得出上述双旋度方程的弱形式为：

为确保得到电场的唯一解，在计算区域的外边界上添加Dirichlet边界条件：

E

进一步的，构建多尺度有限元线性方程组表示为：

K

其中，K

进一步的，根据粗网格棱边处电场获取细网格棱边处的电场的方法为：

获取每个粗网格单元插值算子为h

本发明的有益效果是：根据本发明提出的技术手段，将多尺度有限元用于麦克斯韦方程组的求解，通过基于粗网格对航空电磁响应进行数值模拟，很大程度地降低了正演方程的阶数，进而提高了航空电磁三维数值模拟的计算效率，通过坐标变换的方法引入六面体网格，具有与四面体网格一样拟合不规则物性分界面能力，能够实现对任意起伏地表复杂模型的电磁响应进行快速高精度模拟。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例的一种基于多尺度有限元的复杂模型航空电磁三维快速正演方法流程示意图；

图2为本发明实施例的一种六面体网格示意图；

图3为本发明实施例的一种坐标变换示意图；

图4为本发明实施例的一种算法验证精度结果对比示意图；

图5为本发明实施例的一种多尺度有限元与传统有限元计算时间对比示意图；

图6为本发明实施例的一种多尺度有限元航空电磁含异常体正演模拟结果对比示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，给出了本发明实施例的一种基于多尺度有限元的复杂模型航空电磁三维快速正演方法结构示意图，包括：

设计航空电磁三维正演模拟的复杂模型，根据所述复杂模型将计算区域剖分为六面体网格，得到三维正演初始模型；

在一个具体实施例中，复杂模型包括起伏地表、异常地质体、收发装置类型、发射频率、飞行高度以及网格之间的对应关系，且由于测点分布在整个测区范围内，必须依据完整的信息对模型进行三维精确正演；进而采用规则网格对计算区域进行剖分(该种网格应用更加方便，网格与棱边的位置关系明确，且容易构造粗尺度网格与细尺度网格之间的插值算子，非常适合多尺度有限元方法运用)。

将三维正演初始模型中的每个网格单元作为粗网格，对每个网格单元进行剖分，得到每个网格单元对应的多个细网格；

在一个具体实施例中，利用尺度更小的结构网格对计算区域进一步剖分，视为多尺度有限元中的细网格，之后指定每个粗网格中包含的细网格编号，由多个细网格直接组成对应的粗网格，能够清晰表明粗网格与细网格之间的对应关系，方便后续插值算子的构建，如图2所示，较粗线段为粗尺度网格棱边，每个粗网格内包含27个细网格；需要说明的是，在设计正演模型时，可以任意调整指定粗网格内细网格个数，来达到提高网格质量，降低计算时间等需求。

利用模型网络提取方法获取三维正演初始模型中的核心网格，对核心网格进行延展得到局部网格；所述局部网格包含粗网格和细网格；

在一个具体实施例中已经构建了粗细两套网格，但由于航空电磁计算区域巨大，若在整个计算区域内对网格进行精细剖分，将消耗大量计算资源，因此本实施例采用模型网格提取方法生成局部网格核心部分，之后为了减少边界效应的影响，本实施例中采用Dirichlet狄利克雷边界条件，对网格进行扩边处理，以收发装置为中心进行延展生成局部网格，能够在保证计算精度的前提下提高计算速度，降低计算内存需求。

在一个具体实施例中，生成局部网格后，根据起伏地表信息对两套网格进行拉伸，使形变后的六面体网格能够精确模拟起伏地形与各种不规则地质体，并将异常地质体与背景场的电性信息映射到对应网格中，通过对局部网格进行纵向拉伸来实现拟合起伏地表的目的；由于局部网格以收发装置为中心，因此当收发装置位置发生变化时需重新对网格进行拉伸，但由于网格拉伸时并不改变各个网格大小的比例，因此在构建插值算子时并不需要重新计算插值算子的边界条件。

利用坐标变换计算局部网格中每个细网格的质量矩阵、刚度矩阵，并结合高斯积分和局部网格中细网格的位置关系合成有限元方法中细网格的整体矩阵；

如图3所示，为本发明实施例的一种坐标变换示意图，当运用四面体单元与规则单元对网格进行剖分时，可解析计算质量矩阵与刚度矩阵，但对于六面体单元并没有解析解，只能运用坐标变换的方法进行计算，利用坐标变换与高斯积分计算局部网格中每个细网格的质量矩阵、刚度矩阵，表达式为：

其中，

根据六面体单元之间的节点与棱边关系构建矢量基函数之后，可用如下高斯积分公式进行数值计算：

其中，

利用坐标变换计算局部网格中每个细网格的质量矩阵、刚度矩阵时，频率域电场满足以下双旋度方程：

其中，E

利用伽辽金法与第一格林公式得出上述双旋度方程的弱形式为：

为确保得到电场的唯一解，在计算区域的外边界上添加Dirichlet边界条件：

E

在获得各个细网格单元的整体矩阵后，利用网格单元的位置关系合成三维正演大型线性方程组：K

根据麦克斯韦方程组构建多尺度有限元插值算子；利用所述插值算子对所述细网格的整体矩阵进行降维处理，得到细网格所在粗网格的整体矩阵；

根据麦克斯韦方程组构建多尺度有限元插值算子的方法为：

将插值算子带入到电场双旋度方程中，且由于插值算子仅作用于粗网格内部，令双旋度方程右端量为零，为保证插值算子解的唯一性，令粗网格边界处插值算子值等于粗网格相应位置矢量基函数值，求解该方程即可得到多尺度有限元插值算子，具体为：

其中，v

利用所述插值算子对所述整体矩阵进行降维处理，具体为：

K

其中K

根据粗网格的整体矩阵构建多尺度有限元线性方程组，求解所述多尺度有限元线性方程组，得到粗网格棱边处电场；构建多尺度有限元线性方程组表示为：K

根据粗网格棱边处电场获取细网格棱边处的电场，并根据细网格棱边处电场获取接收点处磁场值；根据粗网格棱边处电场获取细网格棱边处的电场的方法为：获取每个粗网格单元插值算子为h

如图4所示，本实施例给出了一种利用Sasaki(2003)验证本发明算法的精度结果对比示意图，其中，图4(a)为山峰模型示意图；图4(b)为计算区域细网格示意图；图4(c)为计算区域粗网格示意图；图4(d)为传统有限元相对误差结果图；图4(e)为多尺度有限元相对误差结果图；设定发射线圈与接收机距地表高度始终为30m，收发装置为水平共面装置，发射频率16kHz，收发距10m，山峰模型与计算结果参考Sasaki和Nakazato(2003)使用有限差分法的计算结果，具体模型参数可参考图4(a)的山峰模型示意图；从图中可看出本发明正演模拟结果与Sasaki正演结果拟合较好，最大相对误差不超过4％，传统有限元方法使用图4(b)中的网格剖分，而本发明多尺度有限元方法的计算精度与传统有限元方法类似，且具有较高的计算精度。

图5为本实施例给出的一种多尺度有限元与传统有限元计算时间对比示意图，以验证本发明多尺度有限元高效的计算能力，在同一台计算机上利用相同的计算资源分别进行传统有限元与多尺度有限元正演计算，多尺度有限元方法采用的细尺度网格也为传统有限元方法运用的正演模型网格，每个粗尺度网格由2×2×2个细尺度网格组成；根据图上信息可以看出传统有限元程序计算一个测点在一个频率下的电磁响应所需时间为多尺度有限元方法的11倍左右，求解矩阵所需时间为多尺度有限元方法的60倍左右，且计算时间并不随频率变化而变化。采用该网格剖分方法，成功将整体矩阵的维度从938762×938762降低到120816×120816，在精度满足要求的情况下，大幅降低计算时间，节省计算内存。且由于求解矩阵所需时间随着矩阵维度呈指数级上升，随着网格数量的增加，本发明中多尺度有限元方法的优势将进一步扩大。

本实施例进一步还提供了使用多尺度有限元方法模拟山峰地形下存在低阻或高阻异常体模型的航空电磁响应计算结果，并通过将计算结果与无异常体的正演结果进行对比，分析山峰地形下的异常体对航空电磁响应的影响特征，如图6所示，为多尺度有限元航空电磁含异常体正演模拟结果图，其中，图6(a)为山峰模型示意图；图6(b)为计算区域细网格示意图；图6(c)为100Hz响应对比结果图；图6(d)为10kHz响应对比结果图；图6(e)为30kHz响应对比结果图，山峰的底部与顶部宽度分别为200m与40m，高度为40m，在山峰下方20m深处有长宽均为40m，高30m的低阻体或高阻体，低阻体、高阻体和围岩的电导率分别为0.1S/m、0.002S/m和0.01S/m，本实施例分别模拟了100Hz，10kHz，30kHz三个频率，采用图6(b)所示的网格剖分方法，其余参数设置均与精度验证时相同。

根据图6正演结果的对比可得出：所有频率下的模拟结果都受地形起伏影响，使目标异常被地形异常淹没；将起伏地形下方存在的低阻体或高阻体正演结果与无异常体的正演结果进行对比，可明显地看出，低阻体对计算结果的影响更大，更容易被分辨出来；装置的发射频率越高，起伏地表对电磁响应的影响越大，验证了高频信号中的信息主要来自浅地表介质，可以证明地表起伏对电磁响应的影响很大，能够轻易覆盖掉地下异常体响应；因此在进行反演的过程中应充分考虑地表起伏对实测数据造成的影响，本发明将可形变六面体能够成功应用于多尺度有限元方法进行航空电磁正演模拟，同时具有四面体网格对地形的精确模拟能力与矩形块网格最少棱边数的优势，能够快速进行正演模拟并有效分辨起伏地形。

由图6的模拟结果可以看出在起伏地表趋势变化处，响应曲线会出现极值，较容易区分出地表变化位置；起伏地表对电磁响应的影响很大，在处理航空电磁数据时应充分考虑地形对实测数据的影响；装置的发射频率越高，起伏地表对电磁响应的影响越大，验证了高频信号中的信息主要来自浅地表介质；可以证明可形变六面体能够成功应用于多尺度有限元方法，具有与四面体网格同样的拟合不规则物性分界面能力，能够进行快速三维航空电磁正演模拟并有效分辨起伏地形。

根据本发明提出的技术手段，将多尺度有限元用于麦克斯韦方程组的求解，通过基于粗网格对航空电磁响应进行数值模拟，很大程度地降低了正演方程的阶数，进而提高了航空电磁三维数值模拟的计算效率，通过坐标变换的方法引入六面体网格，具有与四面体网格一样拟合不规则物性分界面能力，能够实现对任意起伏地表复杂模型的电磁响应进行快速高精度模拟。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。