一种各向异性介质大地电磁无网格数值模拟方法

摘要

本发明针对地下空间介质中广泛存在的各向异性问题，提出了一种基于离散节点构造形函数的一种大地电磁无网格数值模拟方法。本发明假设各向异性介质的电性主轴一个垂直于层面方面，另外一个平行于层面方向，构造了各向异性介质模型的电导率张量模型。从各向异性大地电磁边值问题出发，推导了对应无网格径向基点插值法的等价线性方程组，形函数采用复合2次径向基函数（MQ-RBF）构造，通过Krylov子空间的

说明书

技术领域

本发明涉及地球物理勘探领域中地下介质的各向异性问题，提出了一种基于离散节点构造形函数的无网格数值模拟方法，适用于大地电磁等地球物理勘探方法。

背景技术

地球内部的结构和属性是地球物理学研究的核心内容。21世纪的地球科学是介质，构造和深层动力学过程的横向不均匀性和各向异性的时代。随着现代观测技术的不断进步以及认识水平的不断提高，各向异性问题逐渐引起了人们的广泛重视，成为地球物理学研究的热点。

大地电磁测深法(MagnetotelluricSounding，MT)是利用天然交变电磁场研究地球电性结构的一种地球物理勘探方法。该方法具有对地下高导构造反映敏感的特点，被广泛应用于地壳、上地幔的研究中。现有的大地电磁正演计算方法主要有积分方程法(IEM)，有限差分法(FDM)，有限单元法(FEM)等。

CN201410370018.0公开了一种时移大地电磁信号采集和反演方法，该方法消除了时移数据中的系统误差和正演过程中的计算误差，减小了系统误差对于所得的不同时刻模型差异的影响，使得反演所得的模型变化更接近真实情况。

CN201010597160.0公开了一种各向异性三维叠前时间偏移方法，该方法考虑地球介质速度各向异性对地震波传播的走时和幅值的影响，能在偏移过程中自主决定偏移速度和各向异性参数，该方法主要用在于地震勘探中反射地震资料处理中。

以上方法都是基于网格实现的，在计算实现的过程中存在着生成网格成本高，场值变化剧烈的地方精度低，自适应分析困难等缺点。本发明针对以上的不足，提出用无网格径向基点插值法(RadicalPointInterpolationMethod,RPIM)模拟大地电磁的各向异性问题，可以实现复杂模型的正演计算。

发明内容

本发明所要解决的问题在于提供一种针对地下空间各向异性介质的，不依赖于网格构造形函数，能够进行复杂介质大地电磁正演数值模拟方法。

本发明是这样实现的，一种各向异性介质大地电磁无网格数值模拟方法，包括如下的步骤：

1)从各向异性大地电磁边值问题出发，构造等价泛函，推导对应无网格径向基点插值法的等价线性方程组；

2)读取当前模型参数，包括频率参数，节点坐标，背景单元，支持域，形状参数，极化模式等；

3)对当前极化模式进行判断，若为TE极化模式，则计算区域包含空气层，若为TM极化模式，则不含空气层；

4)对所有背景网格进行循环，对背景网格的所有高斯积分点循环，搜索该高斯积分点支持域内的有效节点，计算支持域内节点处的形函数，求取系数矩阵和右端项；

5)加载本质边界条件，求解线性方程组，得到各个节点的场值；

6)由视电阻率计算公式，代入场值求取地面处各个方向的视电阻率和相位。

进一步地，步骤1中，各向异性大地电磁边值问题，假设一个电性主轴垂直于层面方面，另外一个主轴平行于层面方向，构造各向异性介质模型的电导率张量模型。

进一步地，步骤4中，计算支持域内节点处的形函数，形函数采用复合2次径向基函数(MQ-RBF)构造。

进一步地，步骤5中，线性方程组的求解采用Krylov子空间的正则化拟残量极小方法(Quasi-minimalResidualmethod，QMR)，实现了大型稀疏线性方程组高效，精确的求解。

本发明与现有技术相比，有益效果在于：本发明针对实际大地电磁探测地下介质中广泛存在的各向异性问题，假设一个电性主轴垂直于层面方面，另外一个主轴平行于层面方向，构造了各向异性介质模型的电导率张量模型。针对传统的正演模拟方法依赖于网格，物性参数复杂分布适应性差的缺点，提出用无网格径向基点插值法模拟大地电磁的各向异性问题，可以有效地克服单纯多项式基点插值方法中存在的奇异性问题，形函数光滑稳定，实现了电磁法高精度，自适应的数值模拟。

附图说明

图1是本发明实施例提供的无网格法场节点、高斯点、求解域及其边界，支持域与背景网格示意图；

图2是本发明实施例提供的各向异性介质的电性主轴与地下空间坐标系示意图；

图3是本发明实施例提供的各向异性介质大地电磁无网格径向基点插值法数值模拟流程图；

图4是本发明实施例建立的层状各向异性介质模型，其中旋转角度为30°；

图5是本发明实施例各项异性层状模型拟解析解与无网格解的对比；

图6是本发明实施例建立的均匀介质中含有各向异性异常体模型；

图7是本发明实施例建立的均匀介质中含有各向异性异常体模型各个方向视电阻率计算结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

参见图1，本发明实施例提供的无网格法计算点1、求解域2，背景单元3，场节点4，支持域5与求解域边界6示意图。

求解域是整个边值问题所要求解的区域，支持域是构造无网格法形函数所选择的区域，计算点是支持域的中心。在计算中，计算点常常选择为高斯积分点。背景单元是便于进行域积分而划分的区域，背景单元并不依赖节点而存在。

参见图3，各向异性介质大地电磁无网格数值模拟方法流程图，包括如下的步骤：

1)从各向异性大地电磁边值问题出发，构造等价泛函，推导对应无网格径向基点插值法的等价线性方程组；

2)读取当前模型参数，包括频率参数，节点坐标，背景单元，支持域，形状参数，极化模式等；

3)对当前极化模式进行判断，若为TE极化模式，则计算区域包含空气层，若为TM极化模式，则不含空气层；

4)对所有背景网格进行循环，对背景网格的所有高斯积分点循环，搜索该高斯积分点支持域内的有效节点，计算支持域内节点处的形函数，求取系数矩阵和右端项；

5)加载本质边界条件，求解线性方程组，得到各个节点的场值；

6)由视电阻率计算公式，代入场值求取地面处各个方向的视电阻率和相位。

当步骤4中计算系数矩阵和右端项结束后检测高斯积分点循环是否结束，结束后进一步检测背景网格是否循环完毕，高斯积分点循环未结束时，返回步骤4对高斯积分点循环，背景网格循环完毕后进入下一步，当背景网格未循环完毕时，则返回到步骤4对背景网络进行循环。

步骤1中各向异性大地电磁边值问题

在介质中，频率域电磁场方程是：

在地下介质空间中，对大地电磁探测所涉及的频率而言，电导率σ>>ωε，因而也可不考虑ε的各向异性问题。电导率σ在各向同性介质中是标量，在各向异性介质中是张量。本发明只考虑σ的各向异性问题。

如图2所示，建立如图所示的坐标系，设平行层面的电导率为σ_//，垂直层面的电导率为σ_⊥,在层面坐标系x'y'z'中x'平行层面,y'垂直层面,z'平行走向，电导率张量为

在地面坐标系中，取走向为z轴(与z'平行)，x轴与z轴垂直，保持水平，y轴垂直向上。在坐标系xyz中，电导率σ张量为

σ＝Aσ'A^T(3)

其中A为坐标变换张量

其中α是x'轴与地下空间坐标系x轴的夹角。

将(1)式按分量展开，并考虑到将电导率张量代入，引入二维算子经化简可得各向异性大地电磁边值问题如下：

其中u指的是场函数E_z或者H_z，Ω指的是求解域，AB是求解域的上边界，AC和BD指的是求解域左右边界，CD指的是求解域下边界，式中：ω代表圆频率，μ代表介质的磁导率，在TE极化模式下，u＝E_z，λ＝σ_//-iωε，在TM极化模式下，u＝H_z，λ＝iωμ，τ的表达式如下：

步骤4中支持域内径向基点插值法建立形函数

如图1所示，在求解域Ω内，径向基点插值函数(RPIM)可表示为

式中R_i(X)为径向基函数(RBF)，n为RBFs的个数，p_j(X)为空间坐标X^T＝(x,y)中的单项式，m为多项式基函数的个数。a_i和b_j为待定常数。本发明选用的径向基函数为复合2次径向基函数(MQ-RBF)，其表达式为：

R_i(x,y)＝(r_i²+(a_cd_c)²)^q(a_c0)(9)

为确定式(8)中的a_i和b_j，需形成计算点X的支持域，其中包括n个场节点。使式(8)满足计算点X周围的n个节点值以确定系数a_i和b_j，这将产生n个线性方程，一个节点对应一个方程。这些方程可表示为矩阵形式

U_s＝R₀a+P_mb(10)

式中的节点场值向量U_s为

U_s＝{u₁u₂Lu_n}^T(11)

RBFs的力矩矩阵为

多项式力矩矩阵为

RBFs的系数向量为

a^T＝{a₁a₂La_n}(14)

多项式系数向量为

b^T＝{b₁b₂Lb_m}(15)

在式(12)中，R_i(r_k)的r_k的表达式为

然而式(8)中有n+m个变量，使用下面m个约束条件添加m个方程

联立式(10)和式(17)可得到如下矩阵方程

式中

求解式(18)可得到

可将式(8)重写为

利用式(21)可得到

式中的RPIM形函数可表示为

最终RPIM形函数Φ(x)可表示为

Φ^T(x)＝{φ₁(x)φ₂(x)Lφ_n(x)}(25)

式(23)可重写为

步骤5中加载本质边界条件

建立各向异性大地电磁边值问题(5)的等价泛函：

求取泛函F(u)对场量u的偏导数，将u(x)＝Φ^T(x)U带入，可得：

其中：

K₂＝(k_ij),k_ij＝∫_ΩλΦ_iΦ_jdΩ

K₃＝(k_ij),k_ij＝∫_CDτkΦ_iΦ_jdΩ

令泛函F(u)对场量u的偏导数为0，从而得到如下矩阵表达式

Ku＝b(29)

其中：K＝K₁-K₂+K₃，u＝{u₁u₂Lu_n}^T，b＝[00...0]。

上边界AB是本质边界，其场值始终是1，为了加载本质边界，将系数矩阵K和右端项b修正为如下形式：

α是一个固定的较大的数，一般可取10⁶～10⁸之间。这样做的优点是使系数矩阵只在两处发生了变化，本质边界条件处理较为简单。

步骤5中线性方程组的求解：

Krlov子空间方法被认为是一种有效的求解大型稀疏线性方程组方法。

基于Krylov子空间的迭代方法收敛速度快，求解精度高，而且稳定性好。Krlov子空间方法常常与预处理技术结合求解线性方程组，本发明采用矩阵分解成一个正交矩阵与一个上三角矩阵的QR预处理方法和正则化拟残量极小方法(QMR)进行线性方程组的求解

步骤6视电阻率和相位的求取：

TE极化模式和TM极化模式下的波阻抗可以通过如下公式计算得到：

TE极化模式:

TM极化模式:

各个方向的阻抗：Z_xy＝Z_TEsin²a+Z_TMcos²a(34)

Z_yx＝-(Z_TEcos²a+Z_TMsin²a)

由卡尼亚视电阻率计算公式，视电阻的表达式为：

将各个方向的阻抗计算公式(34)、(35)代入公式(36)即可得到对应方向上的视电阻率，式(36)的计算结果是复数，对复数求模即为视电阻率，对该复数求幅角即为相位。

为验证算法的正确性，建立如图4所示的实施例，该实施例为层状模型，求解域横向规模为20km，纵向长度为10km。支持域无量纲尺寸α＝1.2，权函数采用四次样条函数，第一层和第三层为各向同性介质，电阻率为100Ω·m，第二层为各向异性介质，平行层面电阻率为10Ω·m，垂直层面电阻率为1000Ω·m。第一层和第二层层厚分别1km和2km。横纵向节点均采用均匀分布，横向间距为200m，纵向间距为100m，形状参数a_c＝1.3，q＝0.5。

图5为在计算区域内地面中心点处的无网格径向基点插值法数值解与拟解析解的测深曲线对比图，拟解析解的来源为Magnetotelluricimpedancesandparametricsensitivitiesfor1-Danisotropiclayeredmedia(Computers&Geosciences，JosefPekaetal，2002,28(8):939-950)。从图中可以看出，在各个频率下无网格法数值解与拟解析解的数据两者高度吻合，充分验证了无网格法求解各向异性大地电磁问题的有效性，为各向异性问题的电磁法正演数值模拟提供了新的思路和方法。

图6为建立的均匀介质中含有的各向异性异常体模型，采用和实施例图4相同的无网格参数。图7为无网格法求解图6的含有各向异性异常体模型四个方向视电阻率r_xx，r_xy，r_yx和r_yy计算结果图，r_xx和r_yy两值相当，但和r_xy，r_yx相差较大,不能较好地反应地层的电性结构。从图中可以看出，xy方向视电阻率和yx方向视电阻率变化趋势相同，较好地反映了地下空间存在的各向异性异常特征。

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