CN105426339A - 一种基于无网格法的线源时域电磁响应数值计算方法 - Google Patents

Abstract

本发明涉及一种基于无网格法的线源时域电磁响应数值计算方法，尤其是可以克服传统数值计算方法中对于网格的依赖，适用于复杂地形下时域电磁探测的数值模拟。本发明基于瞬变电磁法满足的控制方程和定解条件，建立了二维线源边值问题的泛函，利用罚因子法加载本质边界条件，提出旁轴近似方程消除截断边界处的反射波，采用Crack-Nicolson格式进行时间离散，得到递推方程。利用等参单元思想将局部坐标中形状规则的单元离散为节点任意分布的不规则求解对象。采用LU分解方法求解递推方程，最终得到求解区域内各个节点的场值。计算结果表明，该方法形函数光滑性好，模拟精度高，最大误差不超过1×10^-3，实现了电磁法高精度的数值计算。

Description

一种基于无网格法的线源时域电磁响应数值计算方法

技术领域

本发明涉及一种地球物理勘探领域的电磁场数值计算方法，尤其是基于无网格法的线源时域电磁响应数值计算方法。

背景技术

时间域电磁法(Timedomainelectromagneticmethods)或称瞬变电磁法(Transientelectromagneticmethods)，是一种建立在电磁感应原理基础上的人工源电磁探测方法。目前该方法被广泛应用于金属矿勘探、煤矿水文地质调查和工程勘查等领域。

电磁场数值模拟是数据处理和反演解释的基础，在地球物理数值计算中具有基础性、全局性作用。现有瞬变电磁数值计算方法主要有积分方程法(IEM)，有限差分法(FDM)，有限单元法(FEM)等。

CN201410125866.5公开了一种频率域正演方法及装置，该方法通过建立17点格式的差分公式，构造稀疏矩阵，读入子波参数和速度模型，频率循环得到单频波场，对所有频率波场求反傅里叶变换，得到正演结果。该方法增加了差分公式的适用性，通过加载PML(最佳匹配层)边界条件，防止边界反射的波场回折到地表。

CN201110459223.0公开了一种电磁场仿真分析方法，该方法采用棱边单元对三维模型进行网格剖分，引入满足inf-sub条件的标量乘子空间，加入有约束条件的集成后总体矩阵方程求解，得到电磁场分析结果。

美国专利US11756384公开了一种基于三维时域正演建模基础上的反演方法。该方法基于梯度的波形反演理论，通过匹配建模的数据来估计模型参数，在时域建模之后通过离散傅里叶变换(FFT)转换到频率域，完成了三维频域波形反演。

以上所述方法均是基于网格实现的，在计算的过程中都不可避免地涉及到网格的剖分问题，对于单元的形状也有一定的要求，特别是对于物性参数分布复杂和几何特征分布不规则的地电模型适应性差。

发明内容

本发明所要解决的技术问题在于提供一种基于无网格法的线源时域电磁响应数值计算方法，可以克服传统正演方法中对于网格的依赖，适用于复杂地形下时域电磁探测的数值模拟。

本发明是这样实现的，一种基于无网格法的线源时域电磁响应数值计算方法包括：

1)在计算区域内进行节点、电性参数、支持域、背景网格、初始场值设置；

2)在计算区域内对背景网格进行等参数变换，在每个背景单元中分别计算高斯积分点及其权系数；

3)对所有背景网格进行外循环，对所有高斯积分点内循环，搜索背景网格中局部定义域内的有效节点并计算局部定义域内节点处形函数；

4)加载强加边界条件得到线性方程组，加载吸收边界条件得到边界阻尼矩阵，将边界阻尼矩阵加入到线性方程组中；

5)采用LU分解法求解线性方程组，得到各个节点各个时刻的场值。

进一步地，步骤3中局部定义域内利用滑动最小二乘法建立形函数。

进一步地，步骤4中利用罚因子法处理计算模型的强加边界条件，建立和瞬变电磁边值问题的等价泛函数，求取泛函数的变分，得到刚度矩阵K，阻尼矩阵K'和右端项S，形成线性方程组：

$K^{\′}\stackrel{\·}{E}+KE=S$

进一步地，步骤4中根据旁轴近似方程加载吸收边界条件，将左行波作用于右行入射的波，右行波作用于左行入射的波，下行波作用于上行入射的波，相互抵消后消除截断边界处的反射波，得到边界阻尼矩阵：

$F=\frac{1}{v}\×({\∫}_{{\mathrm{\Γ}}_l}{N}_1{N}_2\left(n_x+n_y\right)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_d}{N}_1{N}_2\left(n_x+n_y\right)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_r}{N}_1{N}_2(n_x+n_y\left)d\mathrm{\Γ}\right)$

N₁和N₂均为无网格法边界处两个节点的形函数，n_x和n_y为边界外法线的方向余弦，v为地下介质中的电磁波传播速度。

进一步地，步骤5中线性方程组的时间变量的离散采用Crack-Nicolson格式，令K′+F＝K”，最终形成的递推方程为：

(2K”+ΔtK_t+Δt)E_t+Δt＝ΔtS_t+Δt+S_tΔt+2K”E_t-ΔtK_tE_t

其中E_t为当前时刻的场值，E_t+Δt为待求的下一时刻的场值，S为线性方程组的右端项，Δt为离散的时间步长，K为刚度矩阵。

本发明与现有技术相比，有益效果在于：本发明针对传统的网格化数值计算方式的不足，提出了无网格法模拟的方案。以二维线源瞬变电磁法为例，利用变分原理详细推导了等价线性方程组的形成过程，采用Crack-Nicolson方案离散时间变量，形成最终的递推方程，求取了线源各个节点各个时刻的电磁响应，完整地展现了无网格法电磁数值模拟的全部过程。由于无网格法摆脱了传统方法中对于网格的依赖，理论上节点可以任意布置，为含有起伏地形等复杂情况下电磁法的数值模拟提供新的参考。在实践中同时发现无网格法形函数光滑性好，模拟精度高，为实现电磁法高精度反演打下了坚实的基础。

附图说明

图1是无网格法高斯点、场节点、影响域，局部定义域与背景单元示意图；

图2是二维线源时域电磁探测响应计算模拟示意图；

图3是截断边界处的旁轴近似方程示意图；

图4是任意四边形下的等参数变换示意图；

图5是二维线源时域电磁响应计算方法流程图；

图6是本发明一个实施例节点处的无网格法数值解与理论解对比图；

图7是本发明一个实施例一条剖面中无网格法数值解与理论解对比图。

具体实施方式

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

实施例

参见图5结合图1所示，一种基于无网格法的线源时域电磁响应数值计算方法，包括：

1)设置模型参数，包括节点设置、电性参数设置、背景网格设置、初始场值设置；

如图1所示的无网格法高斯点、场节点、影响域Ω_I，局部定义域Ω_x与背景单元示意图，包括背景积分网络1中的求解域Ω，在一个背景积分网络1的求解域Ω内包括用黑点表示的节点3，以及用空心圆表示的高斯积分点或计算点2，影响域Ω_I是指一个节点的影响区域，影响域Ω_I的中心为场节点，而局部定义域Ω_x则为计算点进行最小二乘拟合所选择的区域，中心常常为一个计算点。在实际计算中常用节点的影响域取代计算点的局部定义域，即在构造一个计算点的无网格形函数时，若计算点位于场节点的影响域中，则该场节点属于构造该形函数的场节点，也即如果场节点的影响域包含该计算点，则该场节点参与构造该计算点的形函数。

2)在计算区域内对背景网格进行等参数变换，在每个背景单元中分别计算高斯积分点及其权系数；

3)对所有背景网格进行循环；

4)对该背景网格的所有高斯积分点循环；

判断高斯积分点是否位于求解域Ω外，如是，则忽略该高斯积分点，处理下个高斯点；

如高斯积分点位于求解域Ω内，计算此计算点的形函数，带入求积公式中计算积分，并将其组装到线性方程组的刚度矩阵和阻尼矩阵中；

5)搜索局部定义域内的有效节点；

6)计算局部定义域内节点处形函数，加载强加边界条件得到线性方程组，加载吸收边界条件得到边界阻尼矩阵，将边界阻尼矩阵加入到线性方程组中，利用Crack-Nicolson方法离散时间变量，得到最终的递推方程；

7)采用LU分解法求解递推方程，得到各个节点各个时刻的场值。

步骤6中，局部支持域Ω_x内滑动最小二乘法建立形函数：

如图1所示，在计算点的局部定义域Ω_x内，构造场函数的近似形式

$u^h\left(x\right)=\underset{j=1}{\overset{m}{\Σ}}{p}_j\left(x\right)a_j\left(x\right)=P^T\left(x\right)a\left(x\right)---\left(1\right)$

p^T(x)代表基函数，m为基函数的维数，a(x)为m维系数矢量，基函数p^T(x)通常为Pascal三角形所决定的单项式以保证其最小完备性。

在局部范围内构造加权离散L2范数J得：

$J=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\ω}}_i\left(x\right){\[u^h\left(x\right)-u_i\]}^2=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\ω}}_i\left(x\right){\[\underset{j=1}{\overset{m}{\Σ}}{p}_j\left(x_i\right)a_j\left(x\right)-u_i\]}^2$

ω_i(x)＝w(x-x_i)称为权函数，n为权函数非零域内的节点个数。以J取最小值得关于a(x)的支配方程：

$\∂J/\∂a=0$

展开可得如下线性方程：

A(x)a(x)＝B(x)U

其中：

将a(x)代入场函数的近似形式(1)，可得：

$u^h\left(x\right)=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{m}{\Σ}}{P}_j\left(x\right){\left(A^{-1}\right(x\left)B\right(x\left)\right)}_{ji}{u}_i=\underset{i=1}{\overset{n}{\Σ}}{N}_i\left(x\right)u_i---\left(2\right)$

N_i(x)即为节点i处的形函数，它是坐标(x,y)的函数。写成矩阵形式即为：

N(x)＝P^T(x)A^-1(x)B(x)(3)

步骤6中，二维线源瞬变电磁无网格法离散：

如图2所示，在二维空间中建立x-z坐标系，z轴以向下为正，二维线源瞬变电磁法的边值问题可以如下描述：

$\left\{\begin{array}{c}\begin{array}{cc}{\▿}^{2}E(x,z,t)-\mathrm{\μ}\mathrm{\σ}\frac{\∂E(x,z,t)}{\∂t}=0& t>0\end{array}\\ E{|}_{t=t_0}=E(x,z)\\ \begin{array}{cc}E{|}_{\mathrm{\Γ}1}=s\left(t\right)& t>0\end{array}\end{array}\right.$

E(x,z,t)为待求的场量，Γ₁为上边界，t₀为初始时刻。μ是地下媒质的磁导率，σ是媒质的电导率，s(t)为上边界的场值。利用罚因子法处理上边界(强加边界条件)，建立和上述边值问题的等价泛函：

$\begin{array}{c}I\left(u\right)=\frac{1}{2}\underset{\mathrm{\Ω}}{\∫\∫}\left({\▿}^{T}E\right(x,z,t\left)\right)(\▿E(x,z,t\left)\right)d\mathrm{\Ω}+\mathrm{\μ}\mathrm{\σ}\underset{\mathrm{\Ω}}{\∫\∫}E(x,z,t)\frac{\∂E(x,z,t)}{\∂t}d\mathrm{\Ω}\\ +\frac{a}{2}\underset{\mathrm{\Γ}1}{\∫}{\left(E\right(x,z,t)-s(t\left)\right)}^{2}d\mathrm{\Γ}\end{array}---\left(4\right)$

a为罚因子，为一固定的常数。将E(x,z,t)＝N(x,z)E(t)代入(4)中得：

$\begin{array}{c}I=\underset{\mathrm{\Ω}}{\∫\∫}\[\frac{1}{2}\[{\left(\frac{\∂N(x,z)E\left(t\right)}{\∂x}\right)}^T\left(\frac{\∂N(x,z)E\left(t\right)}{\∂x}\right)+{\left(\frac{\∂N(x,z)E\left(t\right)}{\∂z}\right)}^T\left(\frac{\∂N(x,z)E\left(t\right)}{\∂z}\right)\]\]d\mathrm{\Ω}\\ +\underset{\mathrm{\Ω}}{\∫\∫}\mathrm{\μ}\mathrm{\σ}{\left(N\right(x,z\left)E\right(t\left)\right)}^{T}N(x,z)\stackrel{.}{E\left(t\right)}d\mathrm{\Ω}\\ +\frac{a}{2}\underset{\mathrm{\Γ}1}{\∫\∫}{\left(N\right(x,z\left)E\right(t)-s(t\left)\right)}^{2}d\mathrm{\Γ}\end{array}$

其中，N(x,z)指的是形函数。

将最后一项展开，可得：

$\begin{array}{c}I=\frac{1}{2}{E}^T\left(t\right)\underset{\mathrm{\Ω}}{\∫\∫}\[{\left(\frac{\∂N(x,z)E\left(t\right)}{\∂x}\right)}^T\left(\frac{\∂N(x,z)E\left(t\right)}{\∂x}\right)+{\left(\frac{\∂N(x,z)E\left(t\right)}{\∂z}\right)}^T\left(\frac{\∂N(x,z)E\left(t\right)}{\∂z}\right)\]d\mathrm{\Ω}E\left(t\right)\\ +{\mathrm{\μ\σE}}^T\left(t\right)\underset{\mathrm{\Ω}}{\∫\∫}{N}^T(x,z)N(x,z)d\mathrm{\Ω}\stackrel{.}{E\left(t\right)}\\ +\frac{a}{2}{E}^T\left(t\right)\underset{\mathrm{\Γ}1}{\∫}{N}^T(x,z)N(x,z)d\mathrm{\Γ}E\left(t\right)-a\underset{\mathrm{\Γ}1}{\∫}N(x,z)s\left(t\right)d\mathrm{\Γ}E\left(t\right)+\frac{a}{2}\underset{\mathrm{\Γ}1}{\∫}{s}^2\left(t\right)d\mathrm{\Γ}\end{array}$

求取泛函I对E^T(t)的偏导数，可得：

$\begin{array}{c}\frac{\∂I}{\∂E^T\left(t\right)}=\underset{\mathrm{\Ω}}{\∫\∫}\[{\left(\frac{\∂N(x,z)}{\∂x}\right)}^T\left(\frac{\∂N(x,z)}{\∂x}\right)+{\left(\frac{\∂N(x,z)}{\∂z}\right)}^T\left(\frac{\∂N(x,z)}{\∂z}\right)\]d\mathrm{\Ω}E\left(t\right)\\ +\mathrm{\μ}\mathrm{\σ}\underset{\mathrm{\Ω}}{\∫\∫}{N}^T(x,z)N(x,z)d\mathrm{\Ω}\stackrel{.}{E\left(t\right)}\\ +a\underset{\mathrm{\Γ}1}{\∫}{N}^T(x,z)N(x,z)d\mathrm{\Γ}E\left(t\right)-a\underset{\mathrm{\Γ}1}{\∫}N(x,z)s\left(t\right)d\mathrm{\Γ}\end{array}$

令 $\∂I/\∂E^T=0,$ 可得：

$\begin{array}{c}\underset{\mathrm{\Ω}}{\∫\∫}\[{\left(\frac{\∂N(x,z)}{\∂x}\right)}^T\left(\frac{\∂N(x,z)}{\∂x}\right)+{\left(\frac{\∂N(x,z)}{\∂z}\right)}^T\left(\frac{\∂N(x,z)}{\∂z}\right)\]d\mathrm{\Ω}E\left(t\right)+a\underset{\mathrm{\Γ}1}{\∫}{N}^T(x,z)N(x,z)d\mathrm{\Γ}E\left(t\right)+\\ \mathrm{\μ}\mathrm{\σ}\underset{\mathrm{\Ω}}{\∫\∫}{N}^T(x,z)N(x,z)d\mathrm{\Ω}\stackrel{.}{E\left(t\right)}=a\underset{\mathrm{\Γ}1}{\∫}N(x,z)s\left(t\right)d\mathrm{\Γ}\end{array}$

可以简化为如下线性方程组：

$K^{\′}\stackrel{\·}{E}+KE=S---\left(5\right)$

其中：

$K=\underset{\mathrm{\Ω}}{\∫\∫}(\frac{\∂N^T}{\∂x}\frac{\∂N}{\∂x}+\frac{\∂N^T}{\∂z}\frac{\∂N}{\∂z})d\mathrm{\Ω}+a\underset{\mathrm{\Γ}}{\∫}{N}^{T}Nd\mathrm{\Γ}$

$\begin{array}{cc}{K}^{\′}=\mathrm{\μ}\mathrm{\σ}\underset{\mathrm{\Ω}}{\∫\∫}{N}^{T}Nd\mathrm{\Ω}& S=a\underset{\mathrm{\Γ}1}{\∫}Ns\left(t\right)d\mathrm{\Γ}\end{array}$

其中，K称为刚度矩阵，K'为阻尼矩阵，S为右端项。

步骤6中，吸收边界条件的加载：

如图3所示，在均匀各向同性介质中，根据StanfordUniversityClaerbout(1985)推导出来的旁轴近似方程，下行波可以表示为：

$\frac{\∂E}{\∂t}=v\frac{\∂E}{\∂y}$

左行波方程为：

$\frac{\∂E}{\∂t}=-v\frac{\∂E}{\∂x}$

右行波方程为：

$\frac{\∂E}{\∂t}=v\frac{\∂E}{\∂x}$

由法向导数的定义知：

$\frac{\∂}{\∂n}=\frac{\∂}{\∂x}{n}_x+\frac{\∂}{\∂y}{n}_y$

其中，n_x和n_y为边界外法线的方向余弦，v为地下介质中的电磁波传播速度。

将左行波作用于右行入射的波，右行波作用于左行入射的波，下行波作用于上行入射的波，相互抵消后实现了消除截断边界处的反射波。利用高斯定理可得：

${\∫}_{\mathrm{\Ω}}\frac{\∂}{\∂x}\left(\frac{\∂E^T}{\∂x}N\right)dxdy={\∫}_{{\mathrm{\Γ}}_l}\frac{\∂E^T}{\∂n}Nd\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_d}\frac{\∂E^T}{\∂n}Nd\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_r}\frac{\∂E^T}{\∂n}Nd\mathrm{\Γ}---\left(6\right)$

其中，Γ_l、Γ_d、Γ_r分别表示左边界、底边界和右边界。这里没有上边界是因为上边界为自由边界条件。N是指形函数。

将左行波，右行波和下行波方程分别代入到(6)式，可以得到：

$\begin{array}{c}{\∫}_{{\mathrm{\Γ}}_l}\frac{\∂E^T}{\∂n}Nd\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_d}\frac{\∂E^T}{\∂n}Nd\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_r}\frac{\∂E^T}{\∂n}Nd\mathrm{\Γ}\\ =\frac{1}{v}\stackrel{\·}{E}\×\left({\∫}_{{\mathrm{\Γ}}_l}{N}_1{N}_2\right(n_x+n_y)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_d}{N}_1{N}_2(n_x+n_y)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_r}{N}_1{N}_2(n_x+n_y)d\mathrm{\Γ})\end{array}$

则边界阻尼矩阵为：

$F=\frac{1}{v}\×({\∫}_{{\mathrm{\Γ}}_l}{N}_1{N}_2\left(n_x+n_y\right)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_d}{N}_1{N}_2\left(n_x+n_y\right)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_r}{N}_1{N}_2(n_x+n_y\left)d\mathrm{\Γ}\right)---\left(7\right)$

将边界阻尼矩阵F加入到线性方程组(5)中，便得到了完整的二维线源瞬变电磁无网格法数值模拟的等价方程：

$(K^{\′}+F)\stackrel{\·}{E}+KE=S$

在时域内进行时间离散：

$\mathrm{\θ}{\left(\frac{\∂E}{\∂t}\right)}_{t+\mathrm{\Δ}t}+(1-\mathrm{\θ}){\left(\frac{\∂E}{\∂t}\right)}_t=\frac{E_{t+\mathrm{\Δ}t}-E_t}{\mathrm{\Δ}t}$

令K′+F＝K”，当时，称为C-N(Crank-Nicolson)格式，最终的递推方程为：

(2K”+ΔtK_t+Δt)E_t+Δt＝ΔtS_t+Δt+S_tΔt+2K”E_t-ΔtK_tE_t(8)

其中E_t为当前时刻的场值，E_t+Δt为待求的下一时刻的场值，S为线性方程组的右端项，Δt为离散的时间步长，K为刚度矩阵。

初始场值由边值条件给定，由方程(8)实现递推从而得到各个时刻各个节点处的场值。

步骤2中，任意背景网格积分的等参数变换：

对于任意设置的背景网格积分区域如任意四边形,可以利用等参单元方法确定高斯积分点和权重。

如图4所示，对于四边形的积分区域，建立(ξ,η)→(x,y)的坐标映射，得到一阶偏导数的雅克比矩阵为：

$J\≡\frac{\∂(x,y)}{\∂(\mathrm{\ξ},\mathrm{\η})}=\left[\begin{array}{cc}\underset{i=1}{\overset{m}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\ξ}}{x}_i& \underset{i=1}{\overset{m}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\ξ}}{y}_i\\ \underset{i=1}{\overset{m}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\η}}{x}_i& \underset{i=1}{\overset{m}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\η}}{y}_i\end{array}\right]=\left[\begin{array}{cccc}\frac{\∂N_1}{\∂\mathrm{\ξ}}& \frac{\∂N_2}{\∂\mathrm{\ξ}}& ...& \frac{\∂N_m}{\∂\mathrm{\ξ}}\\ \frac{\∂N_1}{\∂\mathrm{\η}}& \frac{\∂N_2}{\∂\mathrm{\η}}& ...& \frac{\∂N_m}{\∂\mathrm{\η}}\end{array}\right]\left[\begin{array}{cc}{x}_1& x_1\\ x_2& y_2\\ .& .\\ .& .\\ .& .\\ x_m& y_m\end{array}\right]$

$\[J(\mathrm{\ξ},\mathrm{\η})\]=\left[\begin{array}{cc}\frac{\∂x}{\∂\mathrm{\ξ}}& \frac{\∂y}{\∂\mathrm{\ξ}}\\ \frac{\∂x}{\∂\mathrm{\η}}& \frac{\∂y}{\∂\mathrm{\η}}\end{array}\right]=\left[\begin{array}{cc}\underset{i=1}{\overset{4}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\ξ}}{x}_i& \underset{i=1}{\overset{4}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\ξ}}{y}_i\\ \underset{i=1}{\overset{4}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\η}}{x}_i& \underset{i=1}{\overset{4}{\Σ}}\frac{\∂N_i}{\∂\mathrm{\η}}{y}_i\end{array}\right]$

其中， $N_i=\frac{1}{4}(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1+{\mathrm{\η}}_i\mathrm{\η}),(i=1,2,3,4)$

对于二维单元积分的高斯积分有：

$\begin{array}{c}{\∫}_{S_e}gdS={\∫}_{-1}^{+1}{\∫}_{-1}^{+1}g\left(x\right(\mathrm{\ξ},\mathrm{\η}),y(\mathrm{\ξ},\mathrm{\η}\left)\right)\left|J(\mathrm{\ξ},\mathrm{\η})\right|d\mathrm{\ξ}d\mathrm{\η}\\ ={\∫}_{-1}^1\underset{j=1}{\overset{n}{\Σ}}{H}_{j}g\left(x\right({\mathrm{\ξ}}_j,\mathrm{\η}),y({\mathrm{\ξ}}_j,\mathrm{\η}\left)\right)\left|J(\mathrm{\ξ},\mathrm{\η})\right|d\mathrm{\η}\\ =\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{H}_i{H}_{j}g\left(x\right({\mathrm{\ξ}}_j,{\mathrm{\η}}_i),y({\mathrm{\ξ}}_j,{\mathrm{\η}}_i\left)\right)\left|J({\mathrm{\ξ}}_j,{\mathrm{\η}}_i)\right|\end{array}$

式中，H_i,H_j，(ξ_j,η_i)分别为一维情况下的高斯积分权系数和积分点，|J(ξ_j,η_i)|为对应的雅克比行列式值。

如图2所示的实施例，无限长导线源置于点(0,0)处，线源电流方向平行于地质体走向，影响域无量纲尺寸α＝1，权函数采用三次样条函数，节点横向不均匀分布，源附近较密，最小间距1m，纵向节点均匀分布，间距为10m，均匀半空间电阻率为100Ω·m。

图6为在计算区域内抽取的任意一点的数值解与理论解的对比图，图7为选取的地下任一深度内的剖面曲线，两者高度吻合，最大误差不超过1×10^-3，充分验证了无网格法高精度模拟的特点，为电磁法的数值模拟提供了新的思路和方法。

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

Claims (4)

1.一种基于无网格法的线源时域电磁响应数值计算方法，其特征在于：包括如下的步骤：

1)在计算区域内进行节点、电性参数、支持域、背景网格、初始场值设置；

2)在计算区域内对背景网格进行等参数变换，在每个背景单元中分别计算高斯积分点及其权系数；

3)对所有背景网格进行外循环，对所有高斯积分点内循环，搜索背景网格中局部定义域内的有效节点并计算局部定义域内节点处形函数；

4)加载强加边界条件得到线性方程组，加载吸收边界条件得到边界阻尼矩阵，将边界阻尼矩阵加入到线性方程组中；

5)采用LU分解法求解线性方程组，得到各个节点各个时刻的场值。

2.按照权利要求1所述的一种基于无网格法的线源时域电磁响应数值计算方法，其特征在于：

步骤4中，从线源瞬变电磁满足的控制方程和边界条件出发，建立和时域电磁边值问题的等价泛函数，利用罚因子法处理计算模型的强加边界条件，求取泛函数的变分，得到刚度矩阵K，阻尼矩阵K'和右端项S，E为待求的场值，形成线性方程组：

$K^{\′}\stackrel{\·}{E}+KE=S$

3.按照权利要求2所述的一种基于无网格法的线源时域电磁响应数值计算方法，其特征在于：

步骤4中，吸收边界条件采用旁轴近似方程加载，边界行波相互抵消后消除了反射波，得到边界阻尼矩阵：

$F=\frac{1}{v}\×({\∫}_{{\mathrm{\Γ}}_l}{N}_1{N}_2\left(n_x+n_y\right)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_d}{N}_1{N}_2\left(n_x+n_y\right)d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_r}{N}_1{N}_2(n_x+n_y\left)d\mathrm{\Γ}\right)$

N₁和N₂均为无网格法边界处两个节点的形函数，n_x和n_y为边界外法线的方向余弦，v为地下介质中的电磁波传播速度。

将边界阻尼矩阵F加入到线性方程组中，得到了无网格法线源时域电磁响应数值计算的等价方程组：

$(K^{\′}+F)\stackrel{\·}{E}+KE=S$

4.按照权利要求3所述的一种基于无网格法的线源时域电磁响应数值计算方法，其特征在于：

步骤5中，线性方程组的时间变量的离散采用Crack-Nicolson格式，令K′+F＝K”，最终形成的递推方程为：

(2K”+ΔtK_t+Δt)E_t+Δt＝ΔtS_t+Δt+S_tΔt+2K”E_t-ΔtK_tE_t

其中E_t为当前时刻的场值，E_t+Δt为待求的下一时刻的场值，S为线性方程组的右端项，Δt为离散的时间步长，K为刚度矩阵。