一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法技术

本发明专利技术提供一种预处理器，该预处理器的算法用于获得系数矩阵A，具体是：将有限差分网格上的节点按照耦合关系的强弱重新排序为四种不同类型的节点，以方便每种类型对应的的节点解除耦合，因为算法中各个类型节点的独立性，该方法是高度并行的一个矢量化算法；获得的系数矩阵的A用于正演计算，能够有效地求解大型线性方程组，该预处理器算法自动满足散度条件，无需额外进行散度校正，能够解决低频情况下迭代求解器收敛缓慢或者不收敛的问题。本发明专利技术还公开了一种基于上述预处理器的三维有限差分电磁正演计算方法，计算结果与参考解的误差很小，表明该计算方法的精度高、可靠性好。

A preprocessor and a three-dimensional finite-difference electromagnetic forward calculation method based on the preprocessor

【技术实现步骤摘要】
一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法
本专利技术涉及地球物理
，具体涉及一种预处理器及采用此处理器的三维有限差分电磁正演计算方法。
技术介绍
大地电磁测深法以天然的平面电磁波为场源，通过在地表观测相互正交的电磁场分量来获取地下电性构造信息。由于它无需人工源、装置轻便、探测深度大以及对低阻体反应灵敏等优点，该方法在金属矿产勘探、石油天然气勘探、地热资源勘查以及岩石圈深部构造的研究方面有着广泛的应用。地球物理正演模拟就是通过解析或者数值的方法，在给定地下介质分布及激发源的情况下，计算相应的地球物理响应的过程。通过地球物理正演模拟，我们可以研究不同地球物理模型下响应的分布规律，从而可以指导实际勘探工作。而无论何种地球物理正演模拟方法，最终都会回归到线性方程组的求解问题，常规求解线性方程组有两种主要的方法：第一种：使用直接求解方法来求解稀疏线性方程组，这种方法比较稳健，但是需要花费大量的计算内存和计算时间。第二种：krylov迭代法，在求解大型稀疏线性方程组中更具吸引力。许多学者利用krylov子空间迭代法对三维电磁正演模拟进行了深入的研究，取得了很大的进展。krylov迭代法的性能基本上由其预处理器决定，常用的预处理器有高斯-赛德尔、超松弛和不完全LU分解预处理器。目前，三维电磁正演模拟主要是直接基于电磁场，存在一个问题：随着频率的减小，不满足散度条件，从而导致离散得到的线性方程组很难收敛。所以，寻找一种精度较高且能够极大地加快线性方程组收敛的正演加速方法是电磁法正演中亟需解决的问题。
技术实现思路
本专利技术的第一目的在于提供一种预处理器，具体技术方案如下：一种预处理器，该预处理器的算法用于获得系数矩阵A，包含以下步骤：S1、读取正演模型，获得各个网格单元的沿x、y和z三个方向的棱边长度以及含所有的网格节点信息的初始矩阵；S2、将所有的网格节点分为四种不同类型(具体是按照耦合关系强弱分)，每一种类型的节点被其他三种类型的节点沿x、y和z三个方向所包围，使得相同类型节点之间相互解除耦合，而同类型节点间互不影响，便于同类型节点的集中求解；S3、将初始矩阵按四种类型分解成四个子矩阵；S4、并行实现相同类型的网格节点的计算，最终获得系数矩阵A；边界条件为最近更新的其它类型的网格节点所附着棱边的电场值。以上技术方案中优选的，并行实现相同类型网格节点的计算(即更新各网格节点)是：由包围该六条棱边的24条边界棱边的电场值求解附着在该节点的六条棱边的电场值，从而形成一个局部求解系统，此局部求解系统的所有同种类型的网格节点形成对应的一个子矩阵，各子矩阵分别进行LU分解，然后分别用GS迭代法求解，更新各相同类型网格节点的值；将更新后的网格节点所附着的电场值继续用于求取其他类型节点所附着的电场值；同种类型的节点所附着的电场值同时进行计算。应用本专利技术的预处理器，将有限差分网格上的节点按照耦合关系强弱重新排序为四种类型，以方便每种类型对应的的节点解耦，该方法是高度并行的一个矢量化算法；算法中各个类型节点的独立性，并行实现也是非常便利；获得系数矩阵的矩阵系数A用于正演计算，能够有效地求解大型线性方程组，该预处理器算法自动满足散度条件▽·σE＝0，无需额外进行散度校正，能够解决低频下迭代求解器不收敛或者收敛缓慢的问题。本专利技术的第二目的在于公开一种三维有限差分电磁正演计算方法，技术方案如下：一种三维有限差分电磁正演计算方法，包括以下步骤：步骤一、读取正演模型，获得各个网格单元网格单元的电阻率值以及频率列表；步骤二、结合麦克斯韦方程组和相关物性方程形成双旋度方程，并进行方程的离散化，得到离散化的方程组；步骤三、根据离散化的方程组结合上述的预处理器所得到的矩阵系数A，获得关于Ex、Ey和Ez的线性方程组Ae＝LUe＝b；第四步、求解线性方程组Ae＝b；第五步、输出计算得到的电场分量Ex、Ey和Ez；第六步、根据电场分量Ex、Ey和Ez得到地下介质的视电阻率ρ和相位p。以上技术方案中优选的，所述正演模型为都柏林模型和卡斯卡迪亚模型。以上技术方案中优选的，麦克斯韦方程组详见表达式1)：其中：▽为求旋度符号，E为电场，ω是角频率，σ为电导率，μ为真空磁导率，i2＝-1，磁场H采用表达式H＝(-iωμ)-1▽×E获得；由麦克斯韦方程组表达式1)，对第一个方程求双旋度，后将第二个方程代入第一个方程中得到双旋度方程表达式2)：▽×▽×E+iωμσE＝Es2)；其中：Es为电场源。以上技术方案中优选的，所述第四步具体是：步骤4.1、根据正演模型和频率列表，计算正演模型边界处的电场值，结合上场源项，得到b；步骤4.2、给定初值ek＝e0，计算残差rk＝b-Aek，给定双共轭梯度法的残差步骤4.3、进行迭代循环，若残差达到要求10-10或者超出最大收敛次数3000次，进入下一步；否则，取k＝k+1，返回步骤4.2。以上技术方案中优选的，所述步骤4.3的迭代采用双共轭梯度迭代算法。以上技术方案中优选的，地下介质的视电阻率ρ通过表达式3)获得，地下介质的相位p通过表达式4)获得：其中：H为磁场，e为待求的电场矢量矩阵。应用本专利技术的计算方法，效果是：视电阻率和相位与参考解的最大误差小，精度高；不需要施加散度校正，自动满足散度条件；所有频率的收敛速率基本保持一致，可靠性好，在频率很低的时候，本专利技术算法的优势更加明显。除了上面所描述的目的、特征和优点之外，本专利技术还有其它的目的、特征和优点。下面将参照图，对本专利技术作进一步详细的说明。附图说明构成本申请的一部分的附图用来提供对本专利技术的进一步理解，本专利技术的示意性实施例及其说明用于解释本专利技术，并不构成对本专利技术的不当限定。在附图中：图1是实施例1中预处理器将所有的网格节点按四种类型进行编号的原理图；图2是实施例1中三维有限差分电磁正演计算方法流程图；图3是在XY和YX两种模式下参考解与实施例1中正演算法的结果对比图；图4是实施例2的卡斯卡迪亚模型电导率y-z方向切片示意图；图5是施加散度校正的超松弛预处理器与实施例2中预处理器正演结果的对比图。具体实施方式以下结合附图对本专利技术的实施例进行详细说明，但是本专利技术可以根据权利要求限定和覆盖的多种不同方式实施。实施例1：一种三维有限差分电磁正演计算方法，采用一种新的预处理器，具体为分块单元高斯—赛德尔预处理器，其算法主要是对系数矩阵中的矩阵系数A的优化处理，具体是：S1、读取正演模型(正演模型采用都柏林模型)，获得含所有的网格节点的初始矩阵；S2、将所有的网格节点按四种类型进行编号(O、G、R、B)，如图1，具体示意了一个包含底面、中间面和顶面的三本文档来自技高网...

【技术保护点】
1.一种预处理器，其特征在于，该预处理器的算法用于获得系数矩阵A，包含以下步骤：/nS1、读取正演模型，获得各个网格单元的沿x、y和z三个方向的棱边长度以及含所有的网格节点信息的初始矩阵；/nS2、将所有的网格节点分为四种不同类型，每一种类型的节点被其它三种类型的节点沿x、y和z三个方向所包围；/nS3、将初始矩阵按四种类型分解成四个子矩阵；/nS4、并行实现相同类型的网格节点的计算，最终获得系数矩阵A；边界条件为最近更新的其它类型的网格节点所附着棱边的电场值。/n

【技术特征摘要】
1.一种预处理器，其特征在于，该预处理器的算法用于获得系数矩阵A，包含以下步骤：
S1、读取正演模型，获得各个网格单元的沿x、y和z三个方向的棱边长度以及含所有的网格节点信息的初始矩阵；
S2、将所有的网格节点分为四种不同类型，每一种类型的节点被其它三种类型的节点沿x、y和z三个方向所包围；
S3、将初始矩阵按四种类型分解成四个子矩阵；
S4、并行实现相同类型的网格节点的计算，最终获得系数矩阵A；边界条件为最近更新的其它类型的网格节点所附着棱边的电场值。


2.根据权利要求1所述的预处理器，其特征在于，并行实现相同类型网格节点的计算具体是：
由包围该六条棱边的24条边界棱边的电场值求解附着在该节点的六条棱边的电场值，从而形成一个局部求解系统，此局部求解系统的所有同种类型的网格节点形成对应的一个子矩阵，各子矩阵分别进行LU分解，然后分别用常规的GS迭代法求解，更新各相同类型网格节点的值；将更新后的网格节点所附着的电场值继续用于求取其他类型节点所附着的电场值；
同种类型的节点所附着的电场值可同时进行计算。


3.一种三维有限差分电磁正演计算方法，其特征在于，包括以下步骤：
步骤一、读取正演模型，获得各个网格单元的电阻率值以及频率列表；
步骤二、结合麦克斯韦方程组和相关物性方程形成双旋度方程，并进行方程的离散化，得到离散化的方程组；
步骤三、根据离散化的方程组结合权利要求1-2任意一项所述的预处理器所得到的矩阵系数A，获得关于Ex、Ey和Ez的线性方程组Ae＝LUe＝b；
第四步、求解线性方程组Ae＝b；
第五步、输出计算得到的电场分量Ex、E...

【专利技术属性】
技术研发人员：李健，柳建新，郭荣文，刘嵘，
申请(专利权)人：中南大学，
类型：发明
国别省市：湖南;43