三维频率域可控源数值模拟方法

摘要

本发明提供一种三维频率域可控源数值模拟方法，该方法包括：步骤S100：将计算区域定义在电阻率异常体及其邻近包裹层内，并剖分为多个规则的单元体；步骤S200：计算人工源激发条件下的初始电场，采用加权余量法，在计算区域内建立单元体的棱边中点的有限元方程组；步骤S300：使用格林函数定义观测点为边界节点，使得内部节点的二次电场表达边界节点的二次电场；步骤S400：将格林函数表达式代入混合方程组，得到仅关于内部节点二次电场的第一方程组；步骤S500：求解第一方程组，得到单元体所含节点的二次电场值。该方法运用矢量有限元方法与体积分方程法，缩小了计算区域减少了计算量提高了计算精度，形成方程系数矩阵稀疏且条件数好，计算结果精确。

说明书

技术领域

本发明涉及可控源电磁法正演技术领域，具体的涉及一种三维频率域可控源数值模拟方法。

背景技术

可控源电磁法通过人工源激产生的电磁场，与天然源相比具有抗干扰能力强、效率高等特点。可控源电磁法工作区域已经发展到陆地、海洋、航空与地下。可控源电磁法已经被广泛应用于工程勘察、地热资源开发、金属矿开采、石油勘探等领域的各个时期。

可控源电磁法正演可模拟不同地质条件对应的观测数据，分析不同条件下对应数据的形态特征，有助于定性了解勘探中的一些实际问题。随着计算机的发展，通过对观测数据进行三维反演，可实现定量解释地下电性构造，三维反演通过不断使用与并根据所得结果修改预测模型，来寻找与实测数据相对应的构造模型。而对预测模型进行三维正演则是现有方法中判断预测电性构造模型是否正确的唯一途径。三维可控源电磁法正演方法与理论是三维定量反演的基础与核心，该方法高效而精确。

可控源电磁场数值模拟方法包括微分方程法与积分方程法两大类。

微分方程法的关键在于，求解微分方程。求解微分方程的重要手段是有限元法，有限元法尤其适用于复杂电导率模型，特别是矢量有限元法满足电导剖分单元，边界切向电磁场分量连续性，更满足电磁法实际物理意义，矢量有限元从出现起就引起了地球物理学家的极大兴趣。

积分方程法自从电磁学引进以来就被地球物理学家广泛关注，其研究区域只限于电性异常体，积分方程法因计算区域小而具有快速高效的特点。地球物理学家先后使用体积分方程法模拟了均匀半空间中异常体引起的大地电磁响应与层状各向异性介质中异常体的电磁响应。

虽然微分方程法(以下简称有限元方法)与积分方程法被广泛研究与认可，但是两种方法各具有其待改进的地方或缺陷，有限元方法的缺点具体表现在：

(1)矢量有限元方法使用无穷远边界条件，使得其计算区域不但范围大而且由于包含空气层，形成的线性方程组条件数差，在使用迭代法求解时往往效率低，且随频率变低需要迭代次数不断增加，降低了计算效率；

(2)使用矢量有限元方法时，测点电场值往往需通过测点周围网格节点的插值求得，影响精度。甚至在电性分界面附近，使用插值获得垂直电性分界面电磁场时，会获得错误数据；

积分方程法的缺点具体表现在：

(1)使用积分方程法时，当电流源与电磁场测点在同一位置时，将出现奇异性，难以获得高精度计算结果；

(2)积分方程法形成的系数矩阵满秩且稠密，在对大范围区域剖分后，形成的大型稠密矩阵求解难度大，效率低；

(3)积分方程法中使用的脉冲基函数，前提条件为剖分单元体内电场均匀，对高电导率对比度差异模型并不适用。

发明内容

为了解决上述技术问题，本发明提供了一种三维频率域可控源数值模拟方法，该方法结合运用矢量有限元方法与体积分方程法，缩小了计算区域减少了计算量提高了计算精度，形成方程系数矩阵稀疏且条件数好，计算结果精确。

本发明提供一种三维频率域可控源数值模拟方法，包括以下步骤：

步骤S100：将计算区域定义在电阻率异常体及其邻近包裹层内，并剖分为多个规则的单元体；

步骤S200：计算人工源激发条件下的初始电场，采用加权余量法，在计算区域内建立单元体的棱边中点的有限元方程组，采用节点重排，将有限元方程组的系数分解为关于计算区域的内部节点与边界节点的系数矩阵，

系数矩阵包括第一分块矩阵、第二分块矩阵、第三分块矩阵和第四分块矩阵，第一分块矩阵与第四分块矩阵为对称稀疏矩阵，且第一分块矩阵元素为复数，第四分块矩阵元素包含未知边界结点的二次磁矢量势，第一分块矩阵和第二分块矩阵构成混合方程组，混合方程组以内部节点二次电场与边界节点二次电场为未知数；

步骤S300：使用格林函数定义观测点为边界节点，使得内部节点的二次电场表达边界节点的二次电场；

步骤S400：将格林函数表达式代入混合方程组，得到仅关于内部节点二次电场的第一方程组；

步骤S500：求解第一方程组，得到单元体所含节点的二次电场值；将观测点位置换为观测系统测点位置，重复步骤S300，得到测点的二次电场值。

进一步地，步骤S300还包括以下步骤：采用体积分方程法中的格林函数，使计算区域的内部节点二次电场表达边界节点二次电场。

进一步地，采用格林函数步骤中包括计算异常体内散射电流源产生的二次电场步骤，该步骤中采用矢量有限元线性插值基函数近似单元中心散射电流源。

进一步地，将所得关系式代入上部分方程组中，得到仅关于内部节点二次电场的混合方程组，求解混合方程组，得到内部节点二次电场值。

进一步地，异常体包裹层为剖分步骤中所用网格向x,y,z三个方向分别延伸一个单元的厚度形成。

进一步地，计算初始电场的步骤包括以下步骤：

步骤S210：计算人工源的激发源电流强度J；

步骤S220：计算在层状背景为空气层和/或均匀大地层状条件下，x,y或z方向线载电流发射源在均匀空气和/或大地组成的介质模型中产生的磁矢量势；

步骤S230：通过使用格林函数求得剖分单元体的棱边中点的初始电场值。

本发明的技术效果：

1、本发明提供的三维频率域可控源数值模拟方法，现有的矢量有限元方法需要将计算区域边界扩展至离异常体外很大的范围，才能近似满足第一类边界条件，而本发明提供的方法，将计算区域限定在异常体及其包裹层内，形成的有限元方程组规模小、待求解未知数数量少，提高了工作效率。同时通过将计算区域限定在异常体及其包裹层内，还能避免计算区域内包含空气介质，所得的有限元方程组具有很好的条件数，使用迭代法求解时具有很高的收敛速度，提高了方程求解效率。

2、本发明提供的三维频率域可控源数值模拟方法，直接使用格林函数计算异常体内电流源在测点位置产生的二次电场，避免测点邻近节点时，插值计算带来的误差。提高计算准确性。

3、本发明提供的三维频率域可控源数值模拟方法，在使用格林函数计算散射电流源产生的二次电场时，电流源不再以单元体中心电场结合脉冲基函数表达，而是使用单元体的棱边中心电场结合矢量有限元的函数表达；这种方式更符合物理学常识，单元体内体电流产生二次电场的计算具有更高精度；

4、本发明提供的三维频率域可控源数值模拟方法，将电流源位置与用于计算其产生的二次电场的点的位置分别限定在异常体内与异常体外，避免了格林函数的奇异性，使得格林函数计算精度提高；所得方程组的系数矩阵具有稀疏性质，方便高效求解。

具体请参考根据本发明的三维频率域可控源数值模拟方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。

附图说明

图1是本发明提供的三维频率域可控源数值模拟方法示意图；

图2是本发明优选实施例中计算区域与子单元三维示意图，a)为计算区域三维示意图，其中包括异常体和异常体的包覆层；b)为子单元三维示意图，其中各箭头处的数字和英文表示场所在点位置，数字表示节点编号，图中ABCDEFGH点围成的区域为计算区域，abcdefgh 点围成的区域为异常体区域，计算区域与异常体区域之间的区域为包裹层；

图3是本发明优选实施例中层状异常体及其包覆层内各层之间磁矢势示意图；

图4是本发明优选实施例中海洋三维标准模型示意图；

图5是本发明优选实施例中海洋三维标准模型计算E_x结果对比图；

图6是本发明优选实施例中海洋三维标准模型计算E_z结果对比图；

图7是本发明优选实施例中陆地三维模型示意图；

图8是本发明优选实施例中对比度10：1陆地模型在1hz与100hz激发源E_sx计算结果对比图；

图9是本发明优选实施例中对比度100：1陆地模型在1hz与100hz激发源E_sx计算结果对比图；

图10是本发明优选实施例中对比度500：1陆地模型在1hz与100hz激发源E_sx计算结果对比图；

图11是本发明优选实施例中混合法系数矩阵散点示意图。

具体实施方式

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

参见图1，本发明提供的三维频率域可控源数值模拟方法，包括以下步骤：

步骤S100：将计算区域定义在电阻率异常体及其邻近包裹层内，并剖分为多个规则的单元体；

步骤S200：计算人工源激发条件下的初始电场，采用加权余量法，在计算区域内建立单元体的棱边中点的有限元方程组，采用节点重排，将有限元方程组的系数分解为关于计算区域的内部节点与边界节点的系数矩阵，

系数矩阵包括第一分块矩阵、第二分块矩阵、第三分块矩阵和第四分块矩阵，第一分块矩阵与第四分块矩阵为对称稀疏矩阵，且第一分块矩阵元素为复数，第四分块矩阵元素包含未知边界结点的二次磁矢量势，第一分块矩阵和第二分块矩阵构成混合方程组，混合方程组以内部节点二次电场与边界节点二次电场为未知数；

步骤S300：使用格林函数定义观测点为边界节点，使得内部节点的二次电场表达边界节点的二次电场；

步骤S400：将格林函数表达式代入混合方程组，得到仅关于内部节点二次电场的第一方程组；

步骤S500：求解第一方程组，得到单元体所含节点的二次电场值；将观测点位置换为观测系统测点位置，重复步骤S300，得到测点的二次电场值。

本发明提供的方法通过加权余量法，形成因考虑边界面积分边界条件而导致具有未知边界节点二次磁矢量势的矢量有限元方程组。通过对节点重排将含有未知边界节点二次磁矢量势的系数转移至右下方分块矩阵，从而便于求解。

优选的，步骤S300还包括以下步骤：采用体积分方程法中的格林函数，使计算区域的内部节点二次电场表达边界节点二次电场。

优选的，异常体包裹层为剖分步骤中所用网格向x,y,z三个方向分别延伸一个单元的厚度形成。

优选的，初始电场为单元体的各棱边中点位置处。

优选的，所得关系式代入上部分方程组中，得到仅关于内部节点二次电场的混合方程组，求解混合方程组，得到内部节点二次电场值。

优选的，采用格林函数步骤中包括计算异常体内散射电流源产生的二次电场步骤，该步骤中采用矢量有限元线性插值基函数近似单元中心散射电流源。散射电流源的位置与二次电场计算点的位置分离，避免了格林函数奇异性问题的出现。

优选的，计算初始电场的步骤包括以下步骤：

步骤S210：计算人工源的激发源电流强度J；

步骤S220：计算在层状背景为空气层和/或均匀大地层状条件下，x，y或z方向线载电流源在均匀空气和/或大地组成的介质模型中产生的磁矢量势；

步骤S230：通过使用格林函数求得剖分单元体的棱边中点的初始电场值。

具体的，本发明提供的三维频率域可控源数值模拟方法，包括以下步骤：

1、地电模型的建立及初始场的计算

在层状背景条件介质中，电导率异常体的几何形状及电导率分布，使用规则六面体剖分异常体后，将剖分网格向x,y,z三个方向分别延伸一个单元的厚度，形成异常体包裹层。如当异常体为立方体时，层状背景使其内部分层。具体参见图3～5所示。在此实施例中，如图2a) 所示，计算区域为立方体，以异常体为内核，异常体的各面均垂直自身各面向外延伸一个单位的厚度，所形成的立方体为异常体的包覆层。

2、人工源条件下初始电场的计算

确定电偶极激发源的坐标r′、长度dl、电流I、频率f，电磁场激发源的坐标、及在均匀半空间或层状介质条件中的激发电流源J具有的时谐因子e^iwt，在可控源电磁法频率范围内可以忽略介质内位移电流，频率域电磁场表达为：以下各类不同的电磁场，均用表示，但根据具体的表达式进行区分。

ω＝2πf,

J＝Idl，

其中，σ为介质电导率，ω为激发电流源频率，μ为介质磁导率。此处介质磁导率统一取值为空气磁导率。E与H分别是介质中的电场与磁场，将介质内的电磁、场磁分别分为初始场与二次场：

E＝E_p+E_s,

H＝H_p+H_s,

其中，E、H的下标p和s分别代表电流源在介质中引起的E、H初始电磁场部分和电导率异常体存在引起E、H的二次电磁场部分。初始电磁场满足：

其中，σ_*是均匀半空间或层状介质背景电导率。

参见Tai,C.T.,1971,Dyadic Green's functions in electromagnetic theory:International Textbook Co.，激发电流在空间任意位置产生的电场有如下表达式：

E_P(r)＝g(r,r′)j(r′),

其中，r为初始电场观测位置，g(r,r′)是电场张量的格林函数，V是r′为中心的离散单元体积。

电场张量的格林函数g(r,r′)，表达式：

其中，A为单位长度电偶极源产生的磁矢量势；

层状异性介质条件下不同方向单位电偶极源产生的磁矢量势，参见Xiong,Z.H.,Y.Z.Luo, S.Wang,and G.Wu,1986,Induced-polarization and electromagneticmodeling for three-dimensional body buried in a two-layer anisotropic earth:Geophysics,51,no.2,2235-2246：中公开的推导方法得到。

3、格林函数计算观测点二次电场

在异常体离散单元内存在散射电流源j_s，其在任意位置处产生的二次电场与磁场可表示为：

j_s(r_i)＝(σ(r_i)-σ_*)[E_P(r_i)+E_S(r_i)]，

其中，r_i为第i个离散单元中心点位置，第i离散单元内散射电流在观测点r处产生的二次电场E_si(r)表达式为：

当异常体被剖分为N个离散单元时，在r处的二次电场表达为N个散射电流源在r处的累加：

在层状介质中x方向水平电流源产生的散射电场：

j_x(r_i)＝(σ(r_i)-σ_*)E_x(r_i)

E的上标代表了电偶极源的方向，下标S表示二次场、x与z表示场方向；y方向水平电流源产生同种效果的y电场与z电场；z方向电流源产生的散射电场

j_z(r_i)＝(σ(r_i)-σ_*)E_z(r_i)

以矩阵形式表达为

4、加权余量法及节点重排

介质中二次电场满足以下微分方程：

根据，Jin,J.,2002.Finite element method in electromagnetics,J.W.Wileyand Sons.

将离散子单元坐标为的点的电场定义为单元棱边中心点电场插值表达：

为单元插值形函数,其中取值范围都在(-1,1)；ei、ej、ek分别代表子单元中x、y、 z三个方向。E_pj为离散单元棱边中点初始电场，可以通过已知激发电偶极源计算求得；E_sj为离散单元棱边中点二次电场，为待求取未知数；在单个离散单元内二次电场满足以下微分方程进行变分：

使用矢量变换：

令并使用矢量散度定理得

根据电偶极源计算磁矢量势公式：

有

在单元之间界面有边界条件：

A_ea＝A_eb

A_ea，A_eb，为界面两边无穷趋近于界面的A_e，n为垂直界面方向。在整个区域内，由于单元之间面法向方向相反，将只存在计算区域表面积分项FS

在内部子单元中，求上述一阶变分驻点等价于求下列方程组：

此处，i对应子单元电场编号，x方向电场i＝1，2，3，4；y方向电场i＝5，6，7，8；y方向电场i＝9，10，11，12；对x方向棱边有：

上式可改写为矩阵形式方程组：

对变分求y、z二次电场分量求一阶导数可得：

组成单元三分量方程组：

K_e＝K_eX+K_eY+K_eZ，K_e＝KS_eX+KS_eY+KS_eZ

在子单元中：

其中，a,b,c分别为子单元的x，y，z三个方向边长，按matlab符号积分工具很容易求得子单元刚度矩阵。组装子单元刚度矩阵形成计算区域总体刚度矩阵后，将计算区域边界棱边中心点与区域内部单元棱边中心点，按先内部点后边界点顺序移项排列，得到计算区域内部电场与边界电场分离表达式：

其中K_II与K_BB为对称系数矩阵，K_II是联系内部节点二次电场E_Is之间的关系数，K_BB是联系边界节点二次电场E_Bs之间的关系数；K_IB与K_BI是联系内部节点二次电场与边界节点二次电场的关系数；KS_II、KS_BB、KS_IB、KS_BI分别与K_II、K_BB、K_IB、K_BI有相同形式，区别在于前者联系计算区域初始电场E_IP与E_BP。取方程组的上部分得到混合方程组

K_IIE_Is+K_IBE_Bs＝KS_IIE_IP+KS_IBE_BP

5、建立第一方程组

对于混合方程组

K_IIE_Is+K_IBE_Bs＝KS_IIE_IP+KS_IBE_BP

其中，边界节点二次场可以按格林函数表示为：

如图3所示，计算区域由内部异常体离散单元及单层包裹层组成；散射电流源j_s位于异常体离散单元中心位置处，第i个单元中心处散射电流源j_s(r_i)使用单元棱边中心点电场值结合矢量有限元线性插值函数表达

形成关于仅内部节点二次电场的第一方程组：

K_IIE_Is+K_IBG_BIσ_IE_Is＝S_I-K_IBG_BIσ_IE_IP

当计算区域边界有m个节点，内部有n个节点时：

在混合方法方程组中，K_II为对称稀疏带状矩阵，K_IB为稀疏矩阵，得到已内部结点二次电场为未知数的第一方程组：

K′E_Is＝S′

其中K′为部分稀疏，部分稠密矩阵，通常情况下n远大于m，在n与m相差越大时，K′的稀疏性越好。

6、内部节点及测点二次电场值的计算

求解方程组可得到内部节点处二次场，再次使用格林函数，将观测点定义为观测系统观测点计算测点二次场。

以下结合具体实例对本发明提供的方法在不同发射装置及模型的应用进行详细说明。

本发明计算了三维标准海洋模型，计算区域的确定于剖分参见图2，计算结果与已发表标准数据做了对比。如图4所示，背景介质为海水层与海底层，海水电导率为3.33s/m,海底电导率为1s/m；异常体为水平饼状电阻模型，厚度为100m，顶面距海底1000m，电导率为0.01s/m，具有5km直径，模型中心点位置位于x轴正下方。人工发射源为x方向1hz单位电偶极源，位于模型左边缘上方离海底100m高度，使用同线观测系统，测点沿海底面布置。

如图5所示，本发明方法计算所得x方向电场在幅值与相位都与公认的标准结果具有很好的拟合度。

如图6所示，本发明方法计算z方向幅值与公认的标准结果具有很好的一致性。对海洋三维标准模型的数值模拟结果与已公认结果的对比验证了本发明对海洋频率域可控源电磁法数值模拟的准确性。

为了说明本发明方法对陆地频率域可控源电磁法数值模拟的准确性及对大对比度电导率模型的适用性，本发明算法对不同电导率对比度的陆地三维模型进行了高频100hz与低频1hz 激发条件下数值模拟。

图7所示，一个顶面埋深为100m的200×200×100m尺度异常体处于电导率为0.01s/m的均匀半空间围岩介质中。激发水平电偶极源为x轴方向，长度为100m、电流为10A，源中心处于坐标原点处；异常体中心坐标为(0，5100，150)m，测线首尾端点坐标为(0,5000,0)m 与(400,5000,0)m，测点间距离为20m。在使用标准结果对比验证本发明方法精度的同时，本发明方法又与传统积分方程法做了精度对比。为了保持与积分方法在存在算法更多的相同条件使对比更具有说服力，本发明方法与传统积分方程法对异常体进行了同种方式同样大小网格的剖分；为了得到不同条件下的对比数据，激发电偶极源频率分别使用了1hz与100hz，设置了异常体电导率与背景电导率不同对比度。

在电导对比度为10：1，x方向二次电场相应对比如图8所示，对于高频于低频激发条件下，本发明方法与传统积分方程法计算结果的实部与虚部与均与标准结果有较好的拟合，但本发明方法计算结果与标准结果拟合更好。

在电导对比度为100：1，x方向二次电场相应对比如图9所示，对于高频于低频激发条件下，本发明方法计算结果实部与虚部都与标准结果有较好拟合，而积分方程法计算结果与标准结果开始有明显的偏差。

在电导对比度为500：1，x方向二次电场相应对比如图10所示，结合方法计算结实部结果在高低频条件下都与标准结果有最好很好的拟合。

为了体现本发明提供方法的效率优势，对陆地三维模型使用了粗细网格剖分，并与同等网格条件下的积分方程法做了效率比较，所得结果列于表1中。

表1结合方法与积分方程法的计算效率参数表

张量格林函数数量方程自由度计算时间(1001Hz)积分方程法(粗网格)250,00015008583s混合法(粗网格)624,0001925132126s积分方程法(细网格)16,000,00012,00053625301s混合法(细网格)16,192,00013,65036053570s

由表1可见，采用本发明提供的方法计算时间、计算量远低于各类现有方法。

图11给出了粗网格条件下本发明算法形成的方程矩阵散点图，证实了本方法形成方程组系数矩阵具有稀疏性。

本发明算法计算陆地三维模型电磁响应具体实施如下(频率为1hz，异常体电导率与背景电导率对比度为10：1)：

步骤一：以边长为20m的相同大小立方体作为单元体，将良导异常体剖分为10×10×5个该单元体立方体；分别向x，y，z坐标正负方向延伸一个边长，使得异常体外形成一个单元体厚度的包裹层，如此计算区域尺度为240×240×140m，并被12×12×7个立方体剖分；

步骤二：首先计算激发源电流强度J，其中dl＝100m，I＝10A，f＝1hz；然后计算在层状背景为空气层与均匀大地层状条件下，x方向水平单位电偶极源在空气与均匀大地组成的介质模型中产生的磁矢量势，进而通过引入格林函数计算得到各离散单元棱边中心位置初始电场；

步骤三：使用加权余量法推导得到关于所有单元棱边中点的二次电场及边界单元棱边中点的二次磁矢量势的矢量有限元方程组；单元刚度矩阵计算中，子单元尺度a＝b＝c＝20m,异常体内子单元σ＝0.1s/m，σ_*＝0.01s/m，包裹层内子单元σ＝σ_*＝0.01s/m；

步骤四：对矢量有限元方程组进行移项，使得系数矩阵形成左上与右下分块矩阵为对称稀疏矩阵，且包含未知边界结点二次磁矢量势部分全部在右下方矩阵元素；

步骤五：使用重排后的方程组上半部分，形成以内部节点二次电场与边界结点二次电场为未知数的上部分方程组；

步骤六：使用格林函数，计算10×10×5个异常体单元内散射电流源，对计算区域边界结点形成的二次电场，得到使用异常体内部节点表达计算区域边界节点的表达式将表达式带入上半部分方程组得到关于计算区域内部节点的二次电场方程组；

步骤七：使用稳定双共轭梯度法求解方程组，可以得到异常体内节点二次电场值，再次使用格林函数，将观测点位置为观测系统测点位置得到测点二次电场。

本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。

通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。