一种三维大地电磁正演数值模拟方法
阅读说明：本技术 一种三维大地电磁正演数值模拟方法 (Three-dimensional magnetotelluric forward modeling numerical simulation method ) 是由 杨刚强 柳建新 郭荣文 王永斐 于 2021-09-22 设计创作，主要内容包括：本申请公开了一种三维大地电磁正演数值模拟方法,首先根据目标地质体构建电阻率分布模型,然后通过多重网格法粗化电阻率分布模型,并通过交错网格有限差分法离散粗细网格上的双旋度方程,得到系数矩阵,接着对离散后的双旋度方程作4色分块处理,通过二维有限差分法计算出目标地质体模型的边界条件,并计算出细网格上双旋度方程的右端项,最后使用基于4色分块高斯-赛德平滑的多重网格法求解；更换极化模式并重复上述过程,根据不同极化模式下的电场分量和磁场分量,计算出对应测点的视电阻率和阻抗相位。本申请采用线分块高斯-赛德平滑技术,有效去除粗细网格上的高频残差分量,提高了多重网格法的计算效率,达到快速收敛的目的。(The application discloses a three-dimensional magnetotelluric forward modeling numerical simulation method, which comprises the steps of firstly constructing a resistivity distribution model according to a target geologic body, then coarsening the resistivity distribution model by a multi-grid method, dispersing a double-rotation equation on a coarse grid and a fine grid by a staggered grid finite difference method to obtain a coefficient matrix, then carrying out 4-color blocking processing on the dispersed double-rotation equation, calculating the boundary condition of the target geologic body model by a two-dimensional finite difference method, calculating the right-end term of the double-rotation equation on the fine grid, and finally solving by using a multi-grid method based on 4-color blocking Gaussian-Saidel smoothing (ii) a Changing the polarization mode and repeating the above process according to the electric field in different polarization modesMeasuring the magnetic field component, and calculating apparent resistivity and impedance phase of the corresponding measuring point. The method adopts the line blocking Gaussian-Saider smoothing technology, effectively removes the high-frequency residual error component on the coarse and fine grids, improves the calculation efficiency of the multiple grid method, and achieves the purpose of rapid convergence.)
技术领域
本发明涉及地球物理技术领域，具体地，涉及一种三维大地电磁正演数值模拟方法。
背景技术
大地电磁法勘探是一种利用大地中广泛分布的天然变化电磁场（10-4~104Hz）进行深部地质构造研究的方法，因具有探测深度大、工作效率高、支出成本低等优点，被广泛用于工程勘探、金属和油气资源勘查以及深部地球动力学研究等领域，而能高效地开展电磁三维正反演是其得以广泛应用的前提。
在三维电磁法正演数值模拟过程中，需要对涉及磁场旋度和电场旋度的双旋度方程进行求解，但该方程具有丰富的零空间，而且当频率很低时，该方程往往无法模拟出电性变化，因此采用常规的迭代方法如BICGstab求解该方程，其收敛效率低，必须进行散度校正以加快收敛速度。但随着计算频率越低、剖分网格数越多，迭代次数表现出线性快速增长，因此，再用常规迭代方法模拟研究超深部和复杂的地质异常体，很难具有高效的生产效率。
现有文献公开了多重网格法在扩散电磁场中的应用，为多重网格法在大地电磁模拟中的应用奠定了基础，考虑到多重网格法的收敛速度不会随网格节点数的增加表现出线性增加，因此相比于常规的迭代方法，多重网格法更适合用于对大型方程组求解。但其计算效率很大程度上依赖于平滑算法，而对于电磁场的双旋度方程，常用平滑算法的计算效率也较低。
综上所述，研究人员需要对现有模拟方法进行针对性优化，解决随着计算规模增大、计算效率呈非线性下降，计算频率越低、求解问题越病态，以及常规迭代法收敛速度慢等问题，进而满足大地电磁法勘探中对复杂地质体的快速反演成像需求。
发明内容
本申请的目的在于提供三维大地电磁正演数值模拟方法，在模拟深部复杂地质体产生电场和磁场响应时确保结果精度和计算效率。本申请的技术方案如下：
一种三维大地电磁正演数值模拟方法，包括如下步骤：
步骤1）根据目标地质体的形状、大小和电阻率分布特征构建对应的电阻率分布模型；
步骤2）通过多重网格法将所述电阻率分布模型粗化，得到不同尺度下粗细网格的电阻率分布模型，并计算细网格和粗网格空间关系的限制算子和插值算子；
步骤3）通过交错网格有限差分法离散粗细网格上电场满足的双旋度方程，基于线分块高斯-赛德平滑器原理对离散后的双旋度方程作4色分块处理；
步骤4）根据细网格的电阻率分布模型，通过二维有限差分法计算出目标地质体模型的边界条件，并计算出细网格上双旋度方程的右端项；
步骤5）根据细网格上离散后的双旋度方程，通过线分块高斯-赛德平滑器进行前光滑，得到对应的细网格上的电场分量，并计算出细网格上的电场残差分量；
步骤6）根据所述的细网格上的电场残差分量和限制算子，计算粗网格上的电场残差分量；
步骤7）根据粗网格上离散后的双旋度方程以及电场残差分量，求解该双旋度方程，得到粗网格上的电场迭代量；
步骤8）根据粗网格上的电场迭代量和插值算子，计算出细网格上的电场迭代量；
步骤9）根据细网格上的电场分量和电场迭代量，计算出新的电场分量，并使用线分块高斯-赛德平滑器进行后光滑，更新细网格上的电场分量；
步骤10）重复步骤5）~步骤9），直到细网格上的相对残差小于收敛阈值；
步骤11）根据细网格上的电场分量和旋度算子，计算出细网格上的磁场分量；
步骤12）更换极化模式，重复步骤4）~步骤11）；
步骤13）根据细网格上不同极化模式下的电场分量和磁场分量，计算出对应测点的视电阻率和阻抗相位，从而完成一个计算频率的三维大地电磁正演数值模拟。
在一些具体的实施例中，所述步骤1）包括以下具体过程：
确定研究区域、计算频率和观测方式，布置对应的测网并构建三维长方体模型；
将所述三维长方体模型剖分成若干个细长方体单元，分别将沿、、方向剖分的细长方体单元个数表示为、和，每个细长方体单元的长宽高分别为、和；
根据目标地质体的形状、大小和电阻率分布特征对每个细长方体单元的电阻率进行赋值，进而实现构建对应目标地质体的电阻率分布模型。
在一些具体的实施例中，所述步骤2）包括以下具体过程：
将、、三个方向的细长方体单元两两组合形成一个粗长方体单元，即八个细长方体单元组合成一个粗长方体单元，则沿、、方向剖分的粗长方体单元个数为、和；
粗长方体单元的长宽高及电阻率由其所包含的细长方体单元的长宽高及电阻率决定，进而构建出粗网格的电阻率分布模型；
根据细长方体单元和粗长方体单元上的空间关系，分别计算出限制算子和插 值算子；
基于三维大地电磁正演数值模拟的需求，可进一步粗化模型，得到更粗的电阻率 分布模型，并得到对应的限制算子、……和插值算子、……，为模型粗化次数。
在一些具体的实施例中，所述步骤3）中离散双旋度方程包括以下具体过程：
由麦克斯方程组可推出三维大地电磁电场强度满足的双旋度方程：
(1)
式中：为电场强度，为角频率且与计算频率的关系为，为磁导率，为介质电导率且与电阻率的关系为；
采用交错网格有限差分法，可将公式（1）离散为：
(2)
(3)
式中：为长方体单元上棱边中点的电场分量，为由棱边中点到面中心的旋度 算子，为由面中心到棱边中点的旋度算子，为棱边相连长方体单体的电导率体积平 均；
将步骤1）中构建的三维长方体模型的六个面上的电场分量表示为，即为边界电 场分量，而位于三维长方体模型内部的电场分量用表示，即为内部电场分量，则公式（2） 可表示为：
(4)
式中：和分别为由内部电场分量和边界电场分量分离后的系数矩阵，将项替换为，并移动到等式右边作为等式右端项，公式（4）简化为：
(5)
采用上述方法离散粗细网格上电场满足的双旋度方程，得到对应的系数矩阵，，用、和分别表示细网格上的系数矩阵、电场分量和右端项，用、……表 示粗网格上的系数矩阵。
在一些具体的实施例中，所述步骤3）中对离散后的双旋度方程作4色分块处理包括以下具体过程：
基于线分块高斯-赛德平滑器的工作原理，将得到的电场分量分为4类并用不同颜色标识，同种颜色线上所有节点相连的电场分量相互独立并进行同步计算，而不同颜色间通过相邻线上所有节点相连的电场分量进行边界耦合，即
将分解为4个小的系数矩阵、、和，将分解为、、和，将分解为、、和，对于公式（5）的求解转为对以下4个子方程的求解：
(6)。
在一些具体的实施例中，所述步骤4）包括以下具体过程：
选择或极化，将三维长方体模型切分为一系列二维长方体模型；
通过二维有限差分法计算出二维长方体模型的电场分量，进而得到三维长方体模 型剖分后的所有电场分量，将位于三维长方体模型的六个面上的电场分量赋给，并根据计算出双旋度方程的右端项，将位于内部的电场分量赋给，作为初始估计值并用表示。
在一些具体的实施例中，所述步骤5）~9）包括以下具体过程：
将带入公式（5），并用公式（6）对其进行平滑更新，得到对应的细网格上的电场 分量，同时计算出细网格上的电场残差分量：
(7)
配合限制算子，将细网格上的电场残差分量限制为粗网格上的电场残差分 量：
(8)
根据粗网格上离散得到的系数矩阵，使用直接求解器计算出粗网格上的电场 迭代量：
(9)
配合插值算子，将粗网格上的电场迭代量延拓为细网格上的电场迭代量：
(10)
根据细网格上的电场分量和电场迭代量对电场分量进行粗网格校正，进 而得到新的电场分量；
将校正后的新电场分量带入公式（5），并用公式（6）对其进行平滑更新，得到对应 的细网格上新的电场分量。
在一些具体的实施例中，在所述步骤10）中，当细网格上的电场残差分量满足相对 残差＜1×10-10时，停止迭代计算。
在一些具体的实施例中，所述步骤13）包括以下具体过程：
根据两种极化模式下的电场分量和磁场分量，计算出、两种模式的阻抗：
(11)
式中：、、和为极化计算出、方向的电场分量和磁场分量，、、和为极化计算出、方向的电场分量和磁场分量；
根据所得阻抗计算出、两种模式下的视电阻率和阻抗相位：
(12)。
本申请提供的技术方案至少具有如下有益效果：
1、本申请采用了具有高度并行性的线分块高斯-赛德平滑技术，使得本方法中包含了大地电磁数字模拟中电场分量散度为零的信息，因此无需进行散度校正，加快了多重网格的计算效率；此外还可以有效去除粗细网格上的高频残差分量，采用线分块高斯-赛德平滑器对双旋度方程的系数矩阵进行重新排序分解，分解出的4个小的系数矩阵稀疏性变好，使得本方法在低频时也具有较快的收敛速度和计算效率。
2、本申请对于三维地质体模型使用长方体单元进行剖分，对于每个长方体单元赋值的电阻率可以不同，因此可模拟任意电阻率分布的目标地质体。
3、在本申请中：使用交错网格有限差分法离散长方体单元上电场满足的双旋度方程，有利于简单高效地离散粗细网格上的双旋度方程，进而得到对应的系数矩阵；使用二维有限差分数值解作为边界条件，使得数值模拟结果更准确；使用多重网格法求解交错网格有限差分法离散得到双旋度方程，这样可以将细网格上的低频电场残差分量限制为粗网格上的高频残差分量进行剔除，从而达到快速收敛的目的。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本申请提供的三维大地电磁正演数值模拟方法的流程图；
图2为本申请方法中线分块高斯-赛德平滑器的分块原理图，其中，图2（a）为线平滑系统，图2（b）为线分块示意图，1、2、3、4分别代表4种不同的颜色标识；
图3为本申请实施例1中低-高阻组合异常体模型示意图；
图4为本申请实施例1中采用快速多重网格法与国际公开程序ModEM计算出的低-高阻组合异常体在地表测线上的结果，其中，图4（a）为模式下二者所得视电阻率的对比曲线图，图4（b）为模式下二者所得阻抗相位的对比曲线图，图4（c）为模式下二者所得视电阻率的对比曲线图，图4（d）为模式下二者所得阻抗相位的对比曲线图；
图5为本申请实施例1中采用快速多重网格法与传统迭代法BICGstab在有散度校正和无散度校正下的收敛曲线对比图。
具体实施方式
为了便于理解本申请，下面将结合说明书附图和较佳的实施例对本申请中的技术方案作更全面、细致地描述，但本申请的保护范围并不限于以下具体的实施例，基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，均属于本申请保护的范围。
需要特别说明的是，当某一元件被描述为与另一元件存在“固定、固接、连接或连通”关系时，它可以是直接固定、固接、连接或连通在另一元件上，也可以是通过其他中间件间接固定、固接、连接或连通在另一元件上。
除非另有定义，下文中所使用的所有专业术语与本领域技术人员通常理解的含义相同。本文中所使用的专业术语只是为了描述具体实施例的目的，并不是旨在限制本申请的保护范围。
参见图1，一种三维大地电磁正演数值模拟方法，具体过程如下：
第一步、构建三维电阻率分布模型
根据目标地质体的形状、大小和电阻率分布确定研究区域、计算频率、观测方式，布置对应的测网并构建三维长方体模型。
将该三维长方体模型进行网格剖分，剖分成若干个细长方体单元，分别将沿、、方向剖分的细长方体单元个数表示为、和，每个细长方体单元的长宽高分别为、和。
根据目标地质体的形状、大小和电阻率分布特征对每个细长方体单元的电阻率进行赋值，每个细长方体单元的电阻率可不同，因此可模拟任意情况下的目标地质体构建对应的电阻率分布模型，所述空气层中的每个细长方体单元电阻率赋为。
第二步、粗化三维电阻率分布模型
将、、三个方向的细长方体单元两两组合形成一个粗长方体单元，即八个细长方体单元组合成一个粗长方体单元，则沿、、方向剖分的粗长方体单元个数为、和。
粗长方体单元的长宽高由其所包含的细长方体单元的长宽高对应相加得到，粗长方体单元的电阻率由其所包含的细长方体单元的电阻率体积平均获得，这样就构建出粗网格的电阻率分布模型。
根据细长方体单元和粗长方体单元上的空间关系，分别计算出限制算子和插 值算子。
基于三维大地电磁正演数值模拟的需求，可进一步网格粗化，得到更粗的电阻率 分布模型，并得到对应的限制算子、……和插值算子、……，为模型粗化次数。
第三步、离散粗细网格上的双旋度方程
由麦克斯方程组可推出三维大地电磁电场强度满足的双旋度方程：
(1)
式中：为电场强度（V/m），为角频率（Hz）且与计算频率的关系为， 为磁导率（H/m），在不考虑地下岩石磁化的情况下，取，为介质电导率（S/m）且 与电阻率（）的关系为。
采用交错网格有限差分法，可将公式（1）离散为：
(2)
(3)
式中：为长方体单元上棱边中点的电场分量，为由棱边中点到面中心的旋度 算子，为由面中心到棱边中点的旋度算子，为棱边相连长方体单体的电导率体积平 均。
将构建的三维长方体模型的六个面上的电场分量表示为，即为边界电场分量， 而位于三维长方体模型内部的电场分量用表示，即为内部电场分量，则公式（2）可表示 为：
(4)
式中：和分别为由内部电场分量和边界电场分量分离后的系数矩阵，将项替换为，并移动到等式右边作为等式右端项，公式（4）简化为：
(5)
采用上述方法离散粗细网格上电场满足的双旋度方程，得到对应的系数矩阵，用、和分别表示细网格上的系数矩阵、电场分量和右端项，用、…… 表示粗网格上的系数矩阵。
第四步、线分块高斯-赛德平滑器工作原理推演
由麦克斯韦方程组可推出大地电磁中电场强度和磁场强度的关系：
(13)
由公式（13）可知，电场强度可由磁场强度的旋度获得，磁场强度可由电场强度的旋度获得。
离散后的电场分量，沿着一个长方体单元的面作回路积分可得到元面中心点的磁场分量；磁场分量沿所述交错网格有限差分法中双倍体积单元作回路积分可得到长方体单元棱边中点的电场分量。
因此，与方向一条线上所有节点相连的电场分量只与相邻线上节点的电场分量有关，如图2a所示，系统外围的箭头标识表示相邻一条线上的电场分量，系统内部与中心相连的箭头标识表示当前待求的一条线上所有节点相连的电场分量。
对于不相邻的两条线上所有节点相连的电场分量相互独立，因此，可将第三步中的电场分量归为4类，如图2b中的1、2、3、4标识，相同标识表示为同一类且用一种颜色示出。对同种颜色线上所有节点相连的电场分量相互独立，可进行同步计算，而对不同颜色间通过相邻线上的所有节点相连的电场分量进行边界耦合。
那么，可将分解为4个小的系数矩阵、、和，将分解为、 、和，将分解为、、和，1号色对应，2号色对应，3号色对应，4号色 对应。对于公式（5）的求解转为对以下4个子方程的求解：
(6)
对所述的不同颜色间通过相邻线上的所有节点相连的电场分量进行边界耦合，就 是将公式（6）中的等式按顺序进行计算：先从初值分解出，根据公式（6）中的第一个 式子计算出，更新；然后从更新后的中分解出，根据公式（6）中的第二个式子 计算出，更新；接着从更新后的中分解出，根据公式（6）中的第三个式子计算 出，更新；最后从更新后的中分解出，根据公式（6）中的最后一个式子计算出，更新，即完成一次线分块高斯-赛德平滑。
为了防止计算过程中对系数矩阵的LU分解的重复计算，事先将、、和的LU分解结果保存。
第五步、计算离散后的双旋度方程右端项
根据极化方式，将所述的三维长方体模型切为一系列二维长方形模型，其中，极化是沿方向逐层切出平面的二维长方形模型，极化是沿方向逐层切出平面的二维长方形模型。
通过二维有限差分法计算出二维长方体模型的电场分量，进而得到三维长方体模 型剖分后的所有电场分量，将位于三维长方体模型的六个面上的电场分量赋给，并根据计算出双旋度方程的右端项，将位于内部的电场分量赋给，作为初始估计值并用表示。
第六步、求解公式（5）
多重网格法就是将细网格上的低频残差分量限制为粗网格上的高频残差分量进行剔除，线分块高斯-赛德平滑器可实现对网格上的高频残差分量快速有效地剔除。因此采用基于线分块高斯-赛德平滑器的快速多重网格法求解公式（5）。
使用线分块高斯-赛德平滑器进行前光滑，即将第五步中的带入公式（5），并用 公式（6）对其进行平滑更新，得到对应的细网格上的电场分量，同时计算出细网格上的 电场残差分量：
(7)
配合限制算子，将细网格上的电场残差分量限制为粗网格上的电场残差分 量：
(8)
根据粗网格上离散得到的系数矩阵，使用直接求解器计算出粗网格上的电场 迭代量：
(9)
配合插值算子，将粗网格上的电场迭代量延拓为细网格上的电场迭代量：
(10)
根据细网格上的电场分量和电场迭代量对电场分量进行粗网格校正，得到新的电场分量；
使用线分块高斯-赛德平滑器进行后光滑，即将校正后的新电场分量带入公式 （5），并用公式（6）对其进行平滑更新，得到对应的细网格上新的电场分量。
第七步、迭代循环
重复第六步的计算过程，直到细网格上的电场残差分量满足相对残差＜1× 10-10，停止迭代计算，根据公式（13）可知，将计算出的电场分量取旋度，即可获得地质体 模型中的磁场分量。
为了防止直接求解器对系数矩阵的LU分解的重复计算，事先将系数矩阵的LU分解结果保存到计算机中。
第八步、更换极化模式，重复第五步、第六步和第七步。
第九步、计算地面对应测线上测点的视电阻率和阻抗相位
根据两种极化模式下的电场分量和磁场分量，计算出、两种模式的阻抗：
(11)
式中：、、和为极化计算出、方向的电场分量和磁场分量，、、和为极化计算出、方向的电场分量和磁场分量。
根据所得阻抗计算出、两种模式下的视电阻率和阻抗相位：
(12)
下面对本申请提供的三维大地电磁正演数值模拟方法进行检验。
实施例1
为了对本发明提出的方法进行检验，设计了如图3所示的低-高阻组合异常体模型，详情如下：在地下有低阻异常体和高阻异常体沿方向横向排列的组合异常体，其大小均为10km×10km×10km，其埋深均为10km，其间距为10km；三维低-高阻组合异常体模型的空气层电阻率为，大地层中背景电阻率为100，低阻异常体的电阻率为10，高阻异常体的电阻率为1000。取两个异常体间的中点在地面的投影点作为坐标原点O，测线为，方向-25.5km到25.5km范围内均匀间隔设置25个测点。
若采用32×32×32个大小为2.5km×2.5km×2.5km的长方体单元对大地层进行剖分，则模拟区域为：和方向均为-40km到40km，方向为0km到80km；若采用64×64×64个大小为10km×10km×10km的长方体单元对大地层进行剖分，则模拟区域为：和方向均为-32km到32km，方向为0km到64km；若采用128×128×128个大小为10km×10km×10km长方体单元对大地层进行剖分，则模拟区域为：和方向均为-64km到64km，方向为0km到128km。上述三种网格剖分均可有效模拟图3所述的低-高阻组合异常体模型，对于空气层则采用类似的长方体单元剖分。
图4为采用本申请的三维大地电磁正演数值模拟方法与采用国际公开的ModEM程序计算出的低-高阻组合异常体在地表测线上的、两种模式下的视电阻率和阻抗相位对比曲线图，其中剖分网格大小为64×64×64，计算频率为0.1Hz。由图4可知本申请方法与国际公开ModEM程序在数值模拟结果上基本吻合，验证了本申请方法在结果精度上的可靠信。
图5为本申请方法与传统迭代法BICGstab在有散度校正和无散度校正下的收敛曲线对比图，其中剖分网格大小为64×64×64，计算频率为0.1Hz，计算极化模式为，由图5可知本申请方法在几次多重网格迭代内实现快速收敛。
表1
表1为本实施例中基于线分块高斯-赛德平滑器的快速多重网格法在不同计算频率下的迭代次数和计算时间的统计数据，其中剖分网格大小为64×64×64，计算极化模式为由表1可知本申请方法的计算效率不会随着计算频率变低而变慢，对低频时变态的双旋度方程也具有高效的计算速度，在0.1~0.001Hz范围内，本申请方法的计算速度约为加散度校正的传统迭代法的计算速度的6倍。
表2
表2为本实施例中基于线分块高斯-赛德平滑器的快速多重网格法在不同剖分网格大小下的迭代次数和计算时间统计表，其中计算频率为0.1Hz，计算极化模式为。由表2可知本申请方法的迭代次数不会随剖分网格大小变大而变多，在模拟剖分网格较多的问题时，计算优势更加明显。
因此，本申请方法对于开展复杂异常体的大规模三维大地电磁正演数值模拟研究具有重要意义，可以大大的提高计算效率，减少时间带来的成本开销。
以上所述仅为本申请的部分实施例，并非因此限制本申请的专利保护范围，对于本领域的技术人员来说，本申请可以有各种更改和变化。在本申请的精神和原则之内，凡是利用本申请说明书及附图内容所作的任何改进或等同替换，直接或间接运用在其它相关的技术领域，均应包括在本申请的专利保护范围内。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：一种桥梁健康监测方法、系统、计算机及存储介质