一种各向异性介质的海洋可控源电磁法有限元正演方法与流程

本发明涉及地球物理勘探领域中地下介质的电导率任意各向异性问题，提出了一种基于二次场矢量位、标量位的海洋可控源电磁三维非结构化有限元数值模拟方法，适用于几何复杂、物性分布的航空电磁、井中电磁、陆面电磁等地球物理勘探方法数值模拟研究。

背景技术：

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

海洋可控源电磁法是利用人工源电磁场研究地球电性结构的一种地球物理勘探方法。近十年来，该方法发展迅速，已经成为实现寻找油气田、研究海洋深部构造等勘探目的一种行之有效的勘探方法。现有的海洋可控源电磁数值模拟方法主要有积分方程法(iem),有限差分法(fdm)，有限单元法(fem)等。

目前为止，对于电导率各向异性介质的海洋可控源电磁数值模拟，大都是仅仅考虑三主轴电导率各向异性，这固然简化问题，然而，由地质运动所造成的地层的隆起、翻转等引起的电导率各向异性往往并不局限于三主轴，甚至可能更为复杂，并且在实际探测中，复杂的海底地形与电导率各向异性的影响对实际数据的解释可能会产生一定的偏差。

技术实现要素：

本发明所要解决的技术问题在于提供一种各向异性介质的海洋可控源电磁法有限元正演方法，实现高效、快速的海洋电磁各向异性介质的数值模拟。

本发明是这样实现的，

一种各向异性介质的海洋可控源电磁法有限元正演方法，该方法包括：

1)对研究区域进行非结构化网格剖分，得到网格单元编号，节点编号和坐标参数；

2)读取设计的地电模型参数，包括背景层参数、频率参数、参考电导率、三个欧拉旋转角度、网格单元编号以及节点坐标；

3)计算背景层相应的一次场电磁各分量；

4)对所有网格单元进行循环计算单元系数矩阵，接着总体合成所有网格单元的系数矩阵，接着根据步骤3)背景层计算出来的一次电磁场计算线性方程组的右端项；

5)加载本质边界条件，求解线性方程组，得到各个节点的二次场矢量位及标量位；

6)对二次场矢量位、标量位进行求导，得到所有节点的电磁场各分量。

进一步地，将研究区域剖分成有限多个四面体单元e。

进一步地，该方法还包括设定一个参考电导率，该参考电导率三个非零对角线元素的电导率，接着设定三个欧拉旋转角度，经过三次欧拉旋转，得到任意方向的电导率张量模型。

进一步地，接着从麦克斯韦方程组出发，引入基于coulomb规范下的磁矢量位a、标量位ψ来表示电场、磁场，在二次场中，总磁矢量位和标量分解成二次场与背景场之和，得到关于电磁场二次位表达式：

其中，为异常电导率，ap为一次场矢量位，as为二次场矢量位，ψp为一次场标量位，ψs为二次场标量位。

进一步地，将电磁场二次位表达式按照xyz轴依次展开，利用伽辽金方法对展开式进行加权，结合矢量恒等式以及散度定理，得到关于二次位的体积分方程组。

进一步地，步骤4)中线性方程组：

ku＝b

k为系数矩阵，u为求解域中各个节点待求的场值，b为右端项。

进一步地，步骤5)中求解线性方程组：ku＝b，加入狄利克雷边界条件：(as,ψs)γ＝0。

进一步地，采用不完全lu分解的预条件因子的idr(s)迭代算法对线性方程组进行求解。

进一步地，步骤6)采用指数加权移动平均最小二乘法求解二次场矢量位、标量位对x，y，z的偏导数。

本发明与现有技术相比，有益效果在于：

本发明提出了一种基于二次场矢量位、标量位的非结构化有限元数值模拟算法来模拟电导率任意各向异性的海洋可控源电磁问题。该方法既可以避免源点的奇异性的影响，又能满足节点有限元对于场的连续性要求，同时采用非结构化网格，有利于构建几何复杂地电模型，更好地模拟真实的地质情况，因此，可以实现高效、快速的海洋电磁各向异性介质的数值模拟。

本发明主要针对可控源电磁探测地下介质中广泛存在的电导率任意方向各向异性的问题，首先假设一个电导率张量及三个欧拉旋转角度，通过进行三次欧拉旋转，构造了任意方向各向异性的电导率张量模型。并且首次推导了基于二次场磁矢量位、标量位的电导率任意各向异性条件下的可控源电磁非结构化有限元方程，可以有效避免源点奇异性的影响，提高数值解的精度，又能够满足节点有限元对于场的连续性的要求。为了更好模拟复杂地质结构，本发明采用非结构化网格对研究区域进行离散剖分，可构建复杂的地质结构及物性分布，同时避免了对全区域进行同等程度的加密，减少了计算量。针对常规迭代算法求解线性方程组收敛慢的问题，提出了采用不完全lu分解预处理方法与krylov子空间迭代算法中的idr(s)法相结合对线性方程组求解，有效地提高了计算效率。采用指数加权的滑动平均最小二乘法对二次场矢量位、标量位进行求导，相对于常规的求导方法，进一步提高了导数求解的精度。因此，本发明可以实现高精度、快速的电导率呈任意各向异性条件下的海洋可控源电磁法数值模拟。

附图说明

图1是本发明实施例提供的非结构化网格剖分单元的节点编号示意图；

图2是本发明实施例提供的任意各向异性介质电导率张量的构建相应的三次欧拉旋转示意图；其中，图2(a)沿x轴，图2(b)沿y轴，图2(c)沿z轴；

图3是本发明实施例提供的任意各向异性介质海洋可控源电磁法非结构化有限元数值模拟流程图；

图4是本发明实施例建立的水平层状电导率各向异性介质模型；

图5是本发明实施例各向异性层状模型拟解析解与有限元数值解的对比；

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

图7是本发明实施例建立的任意各向异性介质中含有高阻异常体模型非结构化网格剖分示意图；

图8是本发明实施例建立的任意各向异性介质中含有高阻异常体模型电磁场各分量幅值第一平面图，其中图8(a)、图8(b)、图8(c)为参考电导率不旋转时电场振幅平面图；8(d)、8(e)、8(f)、为参考电导率沿着y轴旋转30。时的电场振幅平面等值线图图；8(g)、8(h)、8(l)为参考电导率沿着y轴旋转45。时的电场振幅平面等值线图；8(m)、8(n)、8(o)为参考电导率沿着y轴旋转60。时的电场振幅平面等值线图；

图9是本发明实施例建立的任意各向异性介质中含有高阻异常体模型电磁场各分量幅值第二平面图，图9(a)、图9(b)、图9(c)分别为参考电导率沿着z轴旋转30。时的电场振幅平面等值线图；图9(d)、图9(e)、图9(f)为参考电导率沿着z轴旋转45。时的电场振幅平面等值线图；图9(g)、图9(h)、图9(l)为参考电导率沿着z轴旋转60。时的电场振幅平面等值线图。

具体实施方式

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

参见图1，本发明实施例提供的非结构化网格剖分的四面体单元节点(1、2、3、4)编号示意图。

采用现成的网格剖分软件，应用任意四面体网格单元对研究区域进行离散剖分，在电性分界处及源点区域进行加密，离激发源较远的地方剖分稀疏些，这样便可避免了对全区域进行同等程度的加密，节省了计算内存，减少了计算量。

参见图3，任意电导率各向异性介质可控源电磁有限元模拟方法流程图，包括如下的步骤：

1)对研究区域进行非结构化网格剖分，得到网格单元编号，节点编号和坐标等参数；

2)读取设计的地电模型参数，包括背景层参数、频率参数、参考电导率、三个欧拉旋转角度、网格单元编号、节点坐标等；

3)计算背景层相应的一次场电磁各分量；

4)对所有网格单元进行循环计算系数矩阵，接着总体合成所有单元的单元系数矩阵，同时，根据背景层计算出来的一次电磁场计算线性方程组的右端项。

5)加载本质边界条件，求解线性方程组，得到各个节点的二次场矢量位及标量位；

6)对二次场矢量位、标量位进行求导，得到所有节点的电磁场各分量。

电导率任意各向异性介质的可控源电磁边值问题地球物理勘探中，往往是采用低频电磁信号，忽略位移电流,采用正弦谐变时间因子e^-(wt，maxwell′s方程可表示如下形式:

其中：w为角频率，μ0为真空中的磁导率。.s为场源电流分布。为任意各向异性电导率张量，e表示的电场，h表示的磁场，如下所示：

为了计算电导率张量通常事先设定一个参考电导率σ4(其中三个对角线的元素表示三个主轴的电导率)，如下所示：

接着通过三次欧拉旋转，参见图2(a)、图2(b)、图2(c)，便可得到任意各向异性介质的电导率张量：

其中旋转矩阵为：

α为坐标轴沿着x轴旋转角度；β为坐标轴沿着y轴旋转角度；γ为坐标轴沿着z轴旋转角度。，此过程对应步骤3。

引入基于coulomb规范下的磁矢量位a、标量位ψ来表示电场、磁场，如下所示：

将方程(6)、(7)代入方程(1)、(2)得：

在二次场算发中，总磁矢量位和标量为可分解成二次场与背景场之和，如下所示：

a＝ap+as(11)

ψ＝ψp+ψs(12)

其中，为背景电导率，为异常电导率，ap为一次场矢量位，as为二次场矢量位，ψp为一次场标量位，ψs为二次场标量位。

将方程(10)-(12)代入方程(8)(9)可得关于电磁场二次位表达式：

通过观察方程(13)、(14)可知，方程右端项是关于一次位，而不含电流项。因此可避免源点的奇异性的影响。而一次位可通过层状介质或者均匀半空间为背景层求得。考虑如下关系：

ep为一次电场。

方程(13)(14)可改写成如下形式：

这样就避免求解一次场标量位的梯度，提高一次场的精度。

有限单元分析:

将方程(13)(14)按照xyz轴依次展开，可得如下方程：

利用伽辽金方法对方程(15)至(18)进行加权，结合矢量恒等式以及散度定理，最后可得关于二次位的体积分方程组：

其中，n为线性插值函数。

本发明中，将研究区域剖分成有限多个四面体单元e，参见图1，在各个单元中，电导率是固定不变的，单元中的二次矢量位和标量位呈线性变化：

线性插值函数定义如下所示：

其中v^e为单元体积，分别为插值函数的系数。

将(23)(24)方程代入(19)-(22)最后可得离散线性方程组：

ku＝b(25)

其中:

b^e＝(bx^e,by^e,b2^e,bψ^e)^t

u^e＝(asx^e,asy^e,as2^e,ψs^e)^t

为了求解方程(25)，还需要加入相应的边界条件。本发明采用狄利克雷边界条件：(as,ψs)γ＝0(26)。

步骤5中线性方程组的求解：krlov子空间迭代算法因收敛速度快，求解精度高，而且稳定性好等优点而被广泛用于求解大型稀疏线性方程组，特别是当其与预处理技术结合能够有效加快求解线性方程组的速度，为了进一步提高求解的效率，本发明采用不完全lu分解的预条件因子的idr(s)迭代算法对线性方程组进行求解。

步骤6电磁场各分量的求取：

电磁场二次场各分量可通过以下方程进行求解：

指数加权移动平均最小二乘法具体实现过程：

首先设单元中的二次位的线性描述：

fi＝t1xi+t2yi+t3zi+d(29)

这里:xi、yi、zi为所需要求解的二次位的网格节点x、y、z坐标。

显而易见，系数t1、t2、t3即为待求解的二次位的导数，即为电磁场二次场各分量，在此直接给出高斯加权的滑动平均最小二乘法所形成的方程，如下所示：

t＝(x^twx)^-1x^twf(30)，

这里

x＝(x,y,z,m),x＝(x1,x2,x3…xn)^t

y＝(y1,y2,y3…yn)^t,z＝(z1,z2,z3…zn)^t

m＝(1,1,1…1)^t

其中，x为表示网格节点三维空间坐标数组，w为权函数，r表示网格节点到坐标原点的距离。h为r最大值。

权函数

h＝max(r)；

这里β为权函数系数，一般取1。f为待求二次位导数点附近n个节点的二次位值，x、y、z为待求二次位的节点坐标。

为了验证本发明中方法的正确性。设计一个三层水平层状沉积模型，如图4所示。采用海水-沉积层两层介质为背景层。设上半空间中的空气电导率为10^(-12)s·m-1。海水电导率为3.2s·m-1，深度为0.3km。距离海水底面深度1km有一水平方向无限延伸的各向异性高阻体，其厚度为100m，参考电导率为σc＝diag(0.01,0.01,0.025)s·m-¹，向y轴旋转30。沉积层的电导率为1s·m-1。采用沿着x方向的水平电偶极子为激发源，其坐标位于(0,0,970)。发射频率为0.25hz。偶极矩为1。沿着z＝970处进行观测。剖分区域大小(-4km≤x,y≤4km,-1km≤z≤5km)，网格单元总数1467492。在amdathlon(tm)x463quad-coreprocessor2.6ghz，内存12g计算平台上执行。

图5为在计算区域内的有限元数值解与拟解析解的电场x分量幅值对比图，拟解析解的来源为losethlo,ursimb.electromagneticfieldsinplanarlylayeredanisotropicmedia[j]。从图5中可以看出，有限元数值解与拟解析解的数据两者高度吻合，充分验证了本发明关于电导率任意各向异性问题的电磁法数值模拟算法的正确性，可为其他电磁法在各向异性介质方面的研究提供一种新的方法。

图6为建立的任意各向异性介质中含有高阻异常体模型，上半空间空气电导率设为10^(-12)，海水深度为1km，电导率为3.3s·m-1。沉积层中有一高阻体，顶面距离海水底面1km，大小为2kmx2kmx0.1km，中心坐标为(2000,0,2050)，电导率为0.01s·m-1。沉积层的参考电导率设为σc＝diag(2,1,1)。激发源为沿着x方向的水平电偶极子，发射频率为0.25hz，偶极距为100，源点坐标(0,0,950)。接收机位于海水深度z＝950处。非结构化网格剖分如图7所示，剖分区域大小(-4km≤x,y≤4km,-1km≤z≤4km)。令沉积层参考电导率分别沿着y轴、z轴旋转0。、30。、45。、60。后研究电场三分量幅值的异常响应特征。

图8，图9分别为沉积层参考电导率沿着y轴跟z轴旋转不同角度后电场振幅平面分布图。图8中图8(a)、8(b)、8(c)为参考电导率不旋转时电场振幅平面图；8(d)、8(e)、8(f)、为参考电导率沿着y轴旋转30。时的电场振幅平面等值线图图；8(g)、8(h)、8(l)为参考电导率沿着y轴旋转45。时的电场振幅平面等值线图；8(m)、8(n)、8(o)为参考电导率沿着y轴旋转60。时的电场振幅平面等值线图。

图9是本发明实施例建立的任意各向异性介质中含有高阻异常体模型电磁场各分量幅值第二平面图，9(a)、9(b)、9(c)分别为参考电导率沿着z轴旋转30。时的电场振幅平面等值线图；9(d)、9(e)、9(f)为参考电导率沿着z轴旋转45。时的电场振幅平面等值线图；9(g)、9(h)、9(l)为参考电导率沿着z轴旋转60。时的电场振幅平面等值线图。

从图8可知，当参考电导率张量沿着y轴旋转的角度30。、45。时，电场ex分量的幅值在x、y方向都有所延伸，特别在y方向的延展更为明显；而电场ey分量正好相反。相比于ex、ey分量，ez分量幅值变化幅度比较大，由一开始的对称到明显的不对称，且右侧的响应值比左侧的大。当参考电导率张量沿着y轴旋转60。时，电场ex分量幅值分布有所收缩；电场ey分量变化不大；而电场ez幅值分布明显向外扩展，且左侧的幅值变大，但右侧幅值依然比左侧大。通过观察图9，可知，参考电导率沿着z轴旋转30。、45。时，电场ex分量分布发生倾斜向外扩张.ey分量的幅值随着旋转的角度增大而有所增大，且在x、y方向分布也有所延伸；ez分量幅值分布发生明显的倾斜变化。当参考电导率沿着z轴旋转60。时，ex分量分布进一步发生扩展；ey分量的幅值分布有所收缩；ez分量幅值分布变成了斜对称。因此，可知不同方向的电导率各向异性所引起的电场异常响应与三主轴电导率各向异性的响应特征有明显的区别。

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