CN109977585A - 一种高精度大地电磁正演模拟方法 - Google Patents

Abstract

本发明提供一种高精度大地电磁正演模拟方法，包括：建立地球物理模型，采用切比雪夫点剖分地电模型；建立切比雪夫求导矩阵；利用切比雪夫求导矩阵，将大地电磁正演计算的偏微分方程转换成代数方程；边界条件处理；大型非对称稠密复系数线性方程组求解；大地电磁响应计算。为提高大地电磁响应的正演计算精度，本发明提出了一种高精度、高适应性和高效的大地电磁正演模拟方法。本发明可以进行大地电磁高精度正演模拟，进而能大幅度地提高大地电磁反演计算的效率，对于促进大地电磁方法在深部资源勘探界的推广应用具有重要的现实意义。

Description

一种高精度大地电磁正演模拟方法

技术领域

本发明涉及地球物理勘探领域，具体涉及一种用于深部资源勘探的高精度大地电磁正演模拟方法。

背景技术

大地电磁测深是通过改变电磁场频率进行测探的一类电法分支方法，它利用电磁感应的趋肤效应，即高频电磁场穿透浅，低频电磁场穿透深，在场源和接收点间距不变的条件下，改变电磁场的频率来达到测深的目的。在勘探领域，正演是指由源的属性推导出场的分布属性，反演则是由场的分布来推导源的属性。

大地电磁正演问题是建立在谐变场麦克斯韦方程组的基础上，依赖于稳定场问题的亥姆霍兹方程求解。对于一维层状介质模型正演问题，大地电磁场具有解析解。但对于一维连续介质、二维地电模型和三维地电模型的正演问题，大地电磁场的解析求解就变得十分困难。除了少数极简单的构造形态外，一般都不能给出解析解，只有借助于解偏微分方程或积分方程的数值计算方法求出近似解。目前，大地电磁正演模拟的数值方法以有限单元法和有限差分法为主。

有限单元法(Finite Element Method)是一种求偏微分(常微分)方程和方程组定解问题数值解的近似计算方法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。有限单元法的优点是适用于模拟形状不规则的地质体和地形起伏，但它的计算程序编制较复杂，尤其是对三维地电条件的模拟更为复杂。

有限差分法(Finite Difference Method)也是一种求偏微分(常微分)方程和方程组定解问题数值解的近似计算方法，其数值计算的基本原理是采用差商来代替微商，把计算微分问题转变为线性方程组来求解，具有简单、灵活以及通用性强等特点，容易在计算机上编程实现。有限差分法在正演计算过程中简便易算，但是当物性参数复杂分布或场域的几何特征不规则时，其适应能力比较差，且网格间距对正演计算精度影响较大。

因此，有必要设计一种新的计算精度高、计算速度快的大地电磁正演模拟方法。

发明内容

因此，本发明提供一种高精度大地电磁正演模拟方法，包括利用基于切比雪夫点的谱方法进行大地电磁正演计算，提高数值模拟精度，所述方法包括如下步骤：步骤A、建立地球物理模型，采用切比雪夫点剖分地电模型；步骤B、建立切比雪夫求导矩阵；步骤C、利用切比雪夫求导矩阵，将大地电磁正演计算的偏微分方程转换成代数方程；步骤D、边界条件处理；步骤E、大型非对称稠密复系数线性方程组求解；步骤F、大地电磁响应计算。

在一种具体的实施方式中，利用切比雪夫点剖分地电模型，所述大地电磁响应计算方法适合于均匀半空间模型、二层G型地电模型、二维地电模型和二维光滑电阻率模型中的一种或多种；步骤A或B中包括选用多项式插值给剖分的网格节点赋电阻率值或电导率值。

在一种具体的实施方式中，利用切比雪夫求导矩阵将偏微分方程转换成代数方程，步骤C中由偏微分方程转换成代数方程的过程中需要使用模型地电参数，步骤C中转换成代数方程后，即建立起大地电磁正演线性方程组；步骤D中的边界条件包括狄利克雷边界条件、诺伊曼边界条件和洛平边界条件，步骤D中在处理边界条件后还包括导出复系数线性方程组。

在一种具体的实施方式中，切比雪夫谱方法形成大型非对称稠密复系数线性方程组的求解，步骤E中采用DILU-BICGSTAB算法迭代求解方程组，步骤E中求解得到各节点的电磁场值，步骤F中计算得到视电阻率和阻抗相位。

本发明至少具有如下有益效果：

为提高大地电磁响应的正演计算精度，本发明提出了一种高精度、高适应性和高效的大地电磁正演模拟方法。本发明可以进行大地电磁高精度正演模拟，进而能大幅度地提高大地电磁反演计算的效率，对于促进大地电磁方法在深部资源勘探界的推广应用具有重要的现实意义。

具体例如，本发明步骤E中，在DILU-BICGSTAB求解过程中，由于预条件矩阵M与系数矩阵A具有相同的稀疏性，在迭代过程中只是系数矩阵非零元素的运算，计算量小，计算速度快。

附图说明

图1为大地电磁正演模型的切比雪夫点离散图，其中图1(a)为TM极化模式的模型剖分示意图，图1(b)为TE极化模式的模型剖分示意图。

图2为切比雪夫谱方法计算二维大地电磁响应导出的线性方程组系数矩阵非零元素分布图，对应本发明步骤E。

图3为实施例1均匀半空间地电模型及切比雪夫点网格剖分图；图3(a)为电导率0.1S/m的均匀半空间模型，图3(b)为网格剖分示意图。

图4为实施例1均匀半空间地电模型的切比雪夫谱方法正演模拟结果。

图5为实施例2一维地电模型及网格剖分图；图5(a)为层状介质地电模型，图5(b)为切比雪夫点剖分示意图。

图6为实施例2一维地电模型的切比雪夫谱方法正演模拟结果；图6(a)为大地电磁视电阻率-频率曲线，图6(b)为大地电磁阻抗相位-频率曲线。

图7为实施例3二维地电模型的正演模拟结果对比图；图7(a)为COMMEMI 2D-1二维地电模型，图7(b)为切比雪夫谱方法与有限差分法模拟的视电阻率对比，图7(c)为切比雪夫谱方法与有限差分法模拟的阻抗相位对比。

图8为实施例4二维光滑电阻率模型。

图9为实施例4二维光滑电阻率模型的正演模拟结果对比图；图9(a)为切比雪夫谱方法计算的视电阻率-频率拟断面图，图9(b)为切比雪夫谱方法计算的阻抗相位-频率拟断面图，图9(c)为有限单元法计算的视电阻率-频率拟断面图，图9(d)为有限单元法计算的阻抗相位-频率拟断面图。

图10为本发明流程图。

具体实施方式

本发明解决的技术问题是，针对现有正演模拟技术的不足，提供了一种新的高精度大地电磁正演模拟方法。

本发明的技术方案：

(1)建立地球物理模型，采用切比雪夫(Chebyshev)点剖分地电模型：

在等距点插值过程中，会出现龙格(Runge)现象，即插值函数在区间的边界处出现震荡。为了消除龙格现象，引入切比雪夫点来代替等距点。在区间[-1,1]内的切比雪夫点定义为：

这些切比雪夫点可以理解为半个单位圆上等距的点在横轴上的投影。图1为大地电磁正演模型的切比雪夫点离散图，其中图1(a)为横磁(TM)极化模式的模型剖分示意图，图1(b)为横电(TE)极化模式的模型剖分示意图。

(2)建立切比雪夫求导矩阵：

在计算区间内确定N+1个切比雪夫点x₀,x₁,…,x_N，在这些点上对函数进行多项式插值，得到最高次幂小于或等于N的插值函数u(x)，然后求插值函数在切比雪夫点处的导数u′(x₀),u′(x₁),…,u′(x_N)。由于上述过程是线性的，可以写成矩阵形式，总结其中的规律并构造切比雪夫求导矩阵，将求导运算转化为矩阵运算。任意N阶切比雪夫求导矩阵D_N中的每个元素的表达式为：

且

对角线上的元素是由对角线以外的元素来计算，即

(3)利用切比雪夫求导矩阵，将大地电磁正演计算的偏微分方程转换成代数方程:

为了实现切比雪夫谱方法计算二维地电模型的大地电磁响应，首先需要空间变量离散化为向量y＝(y₀,y₁,…,y_M)^T和z＝(z₀,z₁,…,z_M)^T，相应地，ρ(y,z)被离散化为下列矩阵形式

为了计算TM极化模式亥姆霍兹方程中偏微分项，需要将u按坐标纵轴z方向排序为一维数组，即

同时，将电阻率参数ρ(电导率的倒数)按坐标纵轴z方向排序为一维数组，即

对于根据切比雪夫求导矩阵可得：

同理可得：

于是，大地电磁正演的偏微分方程可以按切比雪夫求导矩阵写成代数方程形式：

(4)第一类边界、第二类边界条件和第三类边界条件处理：

正演计算的上边界为第一类边界，即狄利克雷边界条件，采用强加边界条件的形式：

下边界条件为第三类边界，即洛平边界条件，写成代数方程形式：

左、右边界为第二类边界，即诺伊曼边界条件，写成代数方程形式：

(5)大型非对称稠密复系数线性方程组求解：

采用切比雪夫谱方法计算大地电磁响应，最终导出的线性方程组的系数矩阵具有稠密形式，且非对称非稀疏，如图2所示。构造DILU-BICGSTAB算法求解线性方程组：

1)给定系数矩阵A，向量b、初始值x₀最大迭代次数k_max，最大容许误差ε_max以及预优矩阵M，计算r₀＝b-Ax₀，并令k＝1,

2)如果k≤k_max且ε≥ε_max转3)，否则停止，输出x_k；

3)如果ρ_k-1＝0或者ω_k-1＝0，算法终止，输出失败信息，否则转4)

4)如果k＝1，则p₁＝r₀，否则计算β_k-1＝ρ_k-1α_k-1/ρ_k-2ω_k-1,p_k＝r_k+1+β_k-1(p_k-1-ω_{k- 1}V_k-1)；

5)求解方程计算s_k＝r_k-1-α_kV_k；

6)ε＝||s_k||,如果ε＞ε_max，否则停止输出x_k；

7)求解ω_k＝tTs/tTt，r_k＝s-ω_kt，ε＝||r_k||，令k＝k+1,转2)。

在DILU-BICGSTAB求解过程中，由于预条件矩阵M与系数矩阵A具有相同的稀疏性，在迭代过程中只是系数矩阵非零元素的运算，计算量小，计算速度快。

(6)大地电磁响应计算：

当计算出节点的电磁场值后，再利用数值方法求出场值沿垂向的偏导数，这样便可计算视电阻和阻抗相位。

对于TE极化模式：

对于TM极化模式：

本发明可以进行大地电磁高精度正演模拟，大幅度地提高了大地电磁反演计算的效率，对于促进大地电磁方法在深部资源勘探界的推广应用具有重要的现实意义。

以下结合附图和具体实施方式对本发明作进一步的说明。

本发明涉及的大地电磁高精度正演计算方法包括以下步骤：

步骤1、建立地质模型，离散计算区域：确定测点的个数与坐标、电阻率异常体位置坐标文件，采用切比雪夫点剖分计算区域，并选用多项式插值给剖分的网格节点赋电阻率值或电导率值。

步骤2、构造切比雪夫求导矩阵，并将大地电磁正演计算的亥姆霍兹方程转换为代数方程，并加入狄利克雷(Dirichlet)边界条件、诺伊曼(Neumann)边界条件和洛平(Robin)边界条件，导出复系数线性方程组。

步骤3、求解线性方程组得到各节点的电磁场值，再利用数值方法求出场值沿垂向的偏导数，进而可计算视电阻率和阻抗相位。

以下为本发明计算大地电磁响应的实例。

实施例1

均匀半空间模型，其电导率值为0.1S/m，如图3(a)所示。正演模拟过程中，取计算区域的长度为8km、宽度为4km，横向剖分单元数Ny＝40，而纵向剖分单元数分别取Nz＝40、Nz＝30、Nz＝20和Nz＝10，其网格剖分示意图见图3(b)。将切比雪夫谱方法数值计算结果与均匀半空间模型磁场响应的解析表达式计算结果做对比，如图4所示，即便剖分网格数目较小时数值近似解与解析解仍然吻合很好，说明本发明基于切比雪夫谱算法的大地电磁正演计算方法的正确性。

实施例2

二层G型地电模型，其模型参数为ρ₁＝10ohm-m，ρ₂＝100ohm-m和h₁＝1000m，如图5所示。切比雪夫谱方法正演计算过程中，横向网格单元取Ny＝40，而纵向网格单元分别取Nz＝100、Nz＝50和Nz＝20。图6给出了切比雪夫谱方法计算横磁(TM)极化模式下G型地电模型所得的视电阻率和阻抗相位曲线，与理论值曲线吻合的较好，这进一步说明本发明正演算法的准确性。但随着纵向单元剖分间距的增加，切比雪夫谱方法的计算精度会下降，主要体现为高频段的视电阻率值和相位值误差增大。通过模拟对比分析，建议取近地表第一个单元的间距Δz＜δ_min/20(这里的δ_min为最高频率值对应的趋肤深度)。

实施例3

二维地电模型如图7(a)所示，在电导率为100ohm-m的围岩中，存在电导率为0.5ohm-m的高导异常体，异常体距离顶部250m。数值模拟过程中，横向网格单元数与纵向网格单元数均取为50。将切比雪夫谱方法数值计算结果与有限差分数值结果做对比，如图7(b-c)所示，本发明的计算结果与美国COMMEMI组织提供的参考解吻合很好，且优于有限差分数值解。

实施例4

二维光滑电阻率模型如图8所示，其原始模型地电参数为：在均匀半空间模型中(电阻率值设置为50ohm-m)，存在3个电性低阻异常体(电阻率值分别设置为10ohm-m、5ohm-m和5ohm-m)以及1个电性高阻异常体(电阻率值设置为100ohm-m)。数值模拟过程中，横向网格单元数取为40，纵向网格单元数均取为50。将切比雪夫谱方法数值计算结果与有限元数值结果做对比，如图9所示，本发明的计算结果与有限元数值计算精度相近，两者均具有高精度的计算结果。

此外，图10为本发明流程图。图10的第四步“利用切比雪夫求导矩阵实现偏微分方程转代数方程，建立大地电磁正演线性方程组”之前还包括“建立切比雪夫求导矩阵”。

以上内容是结合具体的优选实施方式对本发明作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演和替换，都应当视为属于本发明的保护范围。

Claims (10)

1.一种高精度大地电磁正演模拟方法，其特征在于，包括利用基于切比雪夫点的谱方法进行大地电磁正演计算，提高数值模拟精度，所述方法包括如下步骤：

步骤A、建立地球物理模型，采用切比雪夫点剖分地电模型；

步骤B、建立切比雪夫求导矩阵；

步骤C、利用切比雪夫求导矩阵，将大地电磁正演计算的偏微分方程转换成代数方程；

步骤D、边界条件处理；

步骤E、大型非对称稠密复系数线性方程组求解；

步骤F、大地电磁响应计算。

2.根据权利要求1所述的方法，其特征在于，利用切比雪夫点剖分地电模型，所述大地电磁响应计算方法适合于均匀半空间模型、二层G型地电模型、二维地电模型和二维光滑电阻率模型中的一种或多种；步骤A或B中包括选用多项式插值给剖分的网格节点赋电阻率值或电导率值。

3.根据权利要求1所述的方法，其特征在于，利用切比雪夫求导矩阵将偏微分方程转换成代数方程，步骤C中由偏微分方程转换成代数方程的过程中需要使用模型地电参数，步骤C中转换成代数方程后，即建立起大地电磁正演线性方程组；步骤D中的边界条件包括狄利克雷边界条件、诺伊曼边界条件和洛平边界条件，步骤D中在处理边界条件后还包括导出复系数线性方程组。

4.根据权利要求1所述的方法，其特征在于，切比雪夫谱方法形成大型非对称稠密复系数线性方程组的求解，步骤E中采用DILU-BICGSTAB算法迭代求解方程组，步骤E中求解得到各节点的电磁场值，步骤F中计算得到视电阻率和阻抗相位。

5.根据权利要求1所述的方法，其特征在于，为了消除等距点插值过程中的龙格现象，步骤A中具体地：引入切比雪夫点来代替等距点，在区间[-1,1]内的切比雪夫点定义为：

这些切比雪夫点可以理解为半个单位圆上等距的点在横轴上的投影。

6.根据权利要求5所述的方法，其特征在于，建立切比雪夫求导矩阵，步骤B中具体地：

在计算区间内确定N+1个切比雪夫点x₀,x₁,…,x_N，在这些点上对函数进行多项式插值，得到最高次幂小于或等于N的插值函数u(x)，然后求插值函数在切比雪夫点处的导数u′(x₀),u′(x₁),…,u′(x_N)；由于上述过程是线性的，可以写成矩阵形式，总结其中的规律并构造切比雪夫求导矩阵，将求导运算转化为矩阵运算；任意N阶切比雪夫求导矩阵D_N中的每个元素的表达式为：

且

对角线上的元素是由对角线以外的元素来计算，即

7.根据权利要求6所述的方法，其特征在于，利用切比雪夫求导矩阵、将大地电磁正演计算的偏微分方程转换成代数方程，步骤C中具体地：

为了实现切比雪夫谱方法计算二维地电模型的大地电磁响应，首先需要空间变量离散化为向量y＝(y₀,y₁,…,y_M)^T和z＝(z₀,z₁,…,z_M)^T，相应地，ρ(y,z)被离散化为下列矩阵形式

为了计算TM极化模式亥姆霍兹方程中偏微分项，需要将u按坐标纵轴z方向排序为一维数组，即

同时，将电阻率参数ρ(电导率的倒数)按坐标纵轴z方向排序为一维数组，即

对于根据切比雪夫求导矩阵可得：

同理可得：

于是，大地电磁正演的偏微分方程可以按切比雪夫求导矩阵写成代数方程形式：

8.根据权利要求7所述的方法，其特征在于，偏微分方程三类边界条件的切比雪夫谱方法处理，步骤D中具体地：

正演计算的上边界为第一类边界，即狄利克雷边界条件，采用强加边界条件的形式：

下边界条件为第三类边界，即洛平边界条件，写成代数方程形式：

左、右边界为第二类边界，即诺伊曼边界条件，写成代数方程形式：

9.根据权利要求8所述的方法，其特征在于，求解大型非对称稠密复系数线性方程组，步骤E中具体地：

构造DILU-BICGSTAB算法求解线性方程组：

1)给定系数矩阵A，向量b、初始值x₀最大迭代次数k_max，最大容许误差ε_max以及预优矩阵M，计算r₀＝b-Ax₀，并令k＝1，

2)如果k≤k_max且ε≥ε_max转3)，否则停止，输出x_k；

3)如果ρ_k-1＝0或者ω_k-1＝0，算法终止，输出失败信息，否则转4)；

4)如果k＝1，则p₁＝r₀，否则计算β_k-1＝ρ_k-1α_k-1/ρ_k-2ω_k-1,p_k＝r_k+1+β_k-1(p_k-1-ω_k-1V_k-1)；

5)求解方程计算s_k＝r_k-1-α_kV_k；

6)ε＝||s_k||,如果ε＞ε_max，否则停止输出x_k；

7)求解ω_k＝t^Ts/t^Tt，r_k＝s-ω_kt，ε＝||r_k||，令k＝k+1,转2)。

10.根据权利要求9所述的方法，其特征在于，利用数值法近似计算大地电磁响应，步骤F中具体地：

当计算出节点的电磁场值后，再利用数值方法求出场值沿垂向的偏导数，这样便可计算视电阻和阻抗相位；

对于TE极化模式：

对于TM极化模式：