基于无单元Galerkin法电流场数值模拟

Numerical Simulation of Current Field Based on Element-Free Galerkin Method

摘要: 直流电场勘探中，一般采用点电源通过傅氏变换来进行电流场二维正演，相对于利用线源进行二维电流场正演来说，其计算较为复杂；同时，无单元法解决了传统有限元方法的缺点，无单元Galerkin法只对节点进行处理，不受网格的影响，前期处理十分简单，具备精确的计算值、能够求解二次连续的优势。本文通过推导线源直流电的微分方程，采用强加边界条件即第一类边界条件，基于无单元Galerkin法建立了电位场目标函数。然后，由经典的均匀地下半空间模型，通过无单元Galerkin法正演模拟计算与理论值计算结果进行比较分析，验证了该方法的有效性及准确性，最后，在均匀地下半空间模型的基础上，建立了高阻和低阻异常体，通过正演计算，得到了其电位分布与电荷密度分布结果。结果表明该方法适用于对地下构造的正演模拟分析，其能够准确的找出地下含异常的位置，为实际中直流电场反演地下介质信息奠定了基础。

Abstract: In direct current field exploration, two-dimensional current field forward modeling is generally carried out by point power source through Fourier transform. Compared with two-dimensional current field forward modeling by line source, its calculation is more complicated. At the same time, the element-free method solves the shortcomings of the traditional finite element method. The element-free Galerkin method only deals with nodes and is not affected by the grid. The preliminary processing is very simple, and it has the advantages of accurate calculation values and quadratic continuity. In this paper, the objective function of potential field is established based on the element-free Galerkin method by deriving the differential equation of linear source direct current and applying the imposed boundary condition (the first kind of boundary condition). Then, based on the classical underground half-space model, the effectiveness and accuracy of this method are verified by comparing the forward simulation results of element-free Galerkin method with the theoretical calculation results. Finally, on the basis of the uniform underground half-space model, the high and low resistance abnormal bodies are established. Through forward calculation, the potential distribution and charge density distribution were obtained. The results show that this method is suitable for forward modeling and analysis of underground structure, and it can accurately find the location of underground anomaly, which lays a foundation for inversion of underground media information by DIRECT current electric field in practice.

1. 概述

基于偏微分方程的直流(DC)电阻率法数值模拟是计算地球物理中的一个热点问题。较为成熟的数值方法含有限差分法(FDM)和有限元法 [1] [2] [3]。关于有限差分法，推导过程通常很简单，计算成本也很低。然而，有限差分法对复杂模型的适应性差限制了它的应用，因为它的网格是规则分布的。近年来，有限元法因其良好的适应性和精度成为电阻率正演模拟最常用的数值方法 [4] [5] [6]。然而，在有限元计算中，网格离散化是必不可少的，这导致了一些缺点：1) 虽然区域划分的思想非常巧妙，但是分析区域的网格划分可能非常费力和耗时，例如，当网格的数量很大并且计算区域特别由非结构化网格离散化时，通常很难提高局部精度。2) 对于不同类型的单元网格，由于有限元法中形函数的约束条件复杂，形函数的确定是一项困难的工作。

为了避免这些问题，一种无单元Galerkin方法(EFGM) [7] [8] 被引入到直流电场正演模拟中。无单元Galerkin方法没有有限元法的典型缺点。无单元Galerkin方法可以通过任意细化局部区域的节点来提高局部精度。其在地球物理正演研究 [9] [10] [11] [12] [13] 中得到了广泛的运用，最近，麻昌英 [14] 利用径向点插值(RPIM)形函数实现了对直流电场的正演研究，取得了很好的效果。据知，这项工作是第一次尝试应用于基于无单元Galerkin方法构造形函数的DC电阻率建模问题。本文采用无单元Galerkin法对线源直流电场进行二维正演研究，探究在强加边界条件下的稳态电流场响应特征。

2. 线源直流电场数学描述

在二维的直流电阻率法研究中，令z的方向为某空间内的垂直深度，x为垂直于z方向的水平距离，y方向为垂直于空间xz平面的无限延伸走向。因此二维的直流电场中电位满足的微分方程如下表示：

$\frac{\partial }{\partial x}\left(\frac{1}{\rho }\frac{\partial U}{\partial x}\right)+\frac{\partial }{\partial z}\left(\frac{1}{\rho }\frac{\partial U}{\partial z}\right)=\frac{\partial q}{\partial t}$ (1)

为了通过电流来表示直流电场中的电位分布，对上式中等式的右端进行处理，把电荷密度转换成了与电流有关的等式。得到如下方程：

$\nabla \left(\sigma \left(x,z\right)\cdot \nabla U\left(x,z\right)\right)=f_s$ (2)

其中， $\sigma \left(x,z\right)$ 为电导率， $f_s$ 为源项。

对于单极线源与偶极线流源来说，源项 $f_s$ 的表示方法有所差别。表达式如下：

$f_s=\{\begin{array}{l}-I_A\delta \left(x-x_A\right)\delta \left(z-z_A\right)\\ -I_A\delta \left(x-x_A\right)\delta \left(z-z_A\right)-I_B\delta \left(x-x_B\right)\delta \left(z-z_B\right)\end{array}$ (3)

其中，A和B为供电位置，I表示供电的电流大小。

式(2)即为稳定的二维电流场的电流与电位之间的微分方程，通过给定的电流值与供电位置，我们可以求解出空间内相应电位分布。同时，我们还需要对边界条件进行处理，由电流的连续性可以得到，在边界的左右面有：

${\sigma }_1\frac{\partial U}{\partial n}={\sigma }_2\frac{\partial U}{\partial n}$ (4)

其中 $n$ 为法线向量， ${\sigma }_1$ 和 ${\sigma }_2$ 为分界面的内电导率与外电导率。

基于直流电阻率电法中的三类边界条件的选择，本文采用第一类边界条件进行处理，在求解范围内再延拓一定边界，延拓边界就是强行给定的计算求解域外的网格边界，意味着为无穷远处的地下空间，当电流传导到无穷远处时，其电位可令其为 $U=0$。这种边界条件的设定就叫做强加边界条件，在直流电场中被划为第一类边界条件。

针对二维地电模型来说，当使用的收发装置类型不同时，产生的混合边界条件将会有所区别，经过转换，可以的到不同边界条件下的变分式：

单极线电流源：

$\begin{array}{l}F\left(U\right)={\displaystyle {\int }_{\Omega }\left[\frac{\sigma }{2}\left[\frac{\partial }{\partial x}\left(\frac{\partial U}{\partial x}\right)+\frac{\partial }{\partial z}\left(\frac{\partial U}{\partial z}\right)\right]-I_A\delta \left(x-x_A\right)\delta \left(z-z_A\right)U\right]\text{d}\Omega }\\ +\frac{1}{2}{\displaystyle {\int }_{{\Gamma }_{\infty }}\frac{\sigma }{r_A\ln \left(1/r_A\right)}}\cos \left(r_A,n\right)U^2\text{d}\Gamma \\ \delta F\left(U\right)=0\end{array}$ (5)

偶极线电流源：

$\begin{array}{l}F\left(U\right)={\displaystyle {\int }_{\Omega }\frac{\sigma }{2}}\left[\frac{\partial }{\partial x}\left(\frac{\partial U}{\partial x}\right)+\frac{\partial }{\partial z}\left(\frac{\partial U}{\partial z}\right)\right]\text{d}\Omega \\ -{\displaystyle {\int }_{\Omega }\left[I_A\delta \left(x-x_A\right)\delta \left(z-z_A\right)-I_B\delta \left(x-x_B\right)\delta \left(z-z_B\right)U\right]\text{d}\Omega }\\ +\frac{1}{2}{\displaystyle {\int }_{{\Gamma }_{\infty }}\frac{\sigma \left[\left(1/r_A\right)\cos \left(r_A,n\right)-\left(1/r_B\right)\cos \left(r_B,n\right)\right]}{\ln \left(1/r_A\right)-\ln \left(1/r_B\right)}{U}^2\text{d}\Gamma }\\ \delta F\left(U\right)=0\end{array}$ (6)

其中Ω为需要求解的空间，Г为左右边界和底边界， $r$ 为测点与电源之间的矢量距离， $\cos \left(r,n\right)$ 为边界的法相方向同 $r$ 之间的夹角余弦值。本文采用单极线电流源方法进行计算。

通过无单元Galerkin法对直流电场的电位进行形函数的构造，则其电位的场函数可以用n个离散的节点表示为：

$U\left(x\right)=n_1{u}_1+\cdots +n_i{u}_i=N^{\text{T}}U$ (7)

用式(7)表示式(5)则有：

$\begin{array}{c}F\left(U\right)={\displaystyle {\int }_{\Omega }\left[\frac{\sigma }{2}\left[\frac{\partial }{\partial x}\left(\frac{\partial N^{\text{T}}U}{\partial x}\right)+\frac{\partial }{\partial z}\left(\frac{\partial N^{\text{T}}U}{\partial z}\right)\right]-I_A\delta \left(x-x_A\right)\delta \left(z-z_A\right)N^{\text{T}}U\right]\text{d}\Omega }\\ +\frac{1}{2}{\displaystyle {\int }_{{\Gamma }_{\infty }}\frac{\sigma }{r_A\ln \left(1/r_A\right)}}\cos \left(r_A,n\right)U^2\text{d}\Gamma \end{array}$ (8)

通过求解上式的变分问题，将对每个节点进行微分，得到如下表达式：

$\begin{array}{c}\frac{\partial F}{\partial U_i}={\displaystyle \underset{j=1}{\overset{n}{\sum }}{U}_j\sigma {\displaystyle {\iint }_e\left[\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}\right]\text{d}x\text{d}z}}-{\displaystyle {\iint }_{e}f}{N}_i\text{d}x\text{d}z\\ +U_i{\displaystyle {\int }_{\Gamma }\frac{\sigma }{r_A\ln \left(1/r_A\right)}\cos \left(r_A,n\right)N_i{N}_j^{\text{T}}\text{d}\Gamma }\end{array}$ (9)

即：

$\frac{\partial F}{\partial U}=KU-F$ (10)

其中， $K$ 为节点刚度矩阵，形式为：

$K_{ij}={\displaystyle {\iint }_{\Omega }\sigma \left[\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}\right]\text{d}x\text{d}z}+{\displaystyle {\int }_{\Gamma }\frac{\sigma }{r_A\ln \left(1/r_A\right)}\cos \left(r_A,n\right)N^{\text{T}}N}\text{d}\Gamma$ (11)

$F$ 为包括源项的矩阵：

$F={\displaystyle {\iint }_{e}f}{N}_i\text{d}x\text{d}z$ (12)

通过对电位场微分方程的变分求泛函，再代入无单元Galerkin法的形函数的构造，最终得到整个空间域内各节点的总体刚度矩阵，即为需要解决的单极线电流源直流电场的无单元模型。式(4~10)可以表示为：

$KU=F$ (13)

其中， $K$ 为各节点刚度矩阵之和， $U$ 为节点电势， $F$ 为右端项。对于第一类边界条件来说，由于进行了网格延拓，即外边界为无穷远处，在外边界处，其各处计算点的电势为0，因此右端项 $F$ 只表示为含有源项的矩阵。因此 $F$ 的表达式又可以表示为：

$F={\displaystyle {\iint }_{e}f}{N}_i\text{d}x\text{d}z=\{\begin{array}{l}I在节点\left(x_0,z_0\right)上\\ 0不在节点\left(x_0,z_0\right)上\end{array}$ (14)

3. 线源直流电场正演模拟

3.1. 均匀半空间地电模型

为了验证无单元Galerkin法在电阻率直流电场中的有效性及正确性，同时使直流电场中的电位分布有更加清晰的认识，本小节运用无单元Galerkin法对经典的半空间地电模型进行了正演数值模拟。

本部分电阻率设为100 Ω∙m，供电电极A坐标为(0,0)，MN为接收点。其电流供电为1 A，半空间求解区域为长50 m，地下深度为25 m，并在半空间的左右边界及其底边界增加了延拓边界令其外边界电位值为0，网格尺寸设置为50*25，节点间距为1 m。网格布置如图1所示。通过推导的电流场无单元Galerkin形式的目标函数，得到图2电位等值线图，在均匀的电阻率值下，半空间内的电位等势线是一个围绕单极线源点中心向外均匀扩散半圆形状态。

Figure 1. Grid layout of uniform half-space model

图1. 均匀半空间模型网格布置图

Figure 2. Potential contour map of homogeneous half-space model

图2. 均匀半空间模型电位等值线图

从图2可以看出，在均匀半空间内电位离源点越远其值越小，并且围绕源点以圆域的方式逐渐减弱。越靠近源点其电位值越大，通过与其他文献中有限元法的电位模拟电位等值线图比较，其电位分布符合实际情况，证明了本方法的有效性及其正确性。本文还通过与理论值的计算结果进行验证分析，进一步证明了该方法的实用性。

均匀半空间内单级线源的计算公式为：

$U=\frac{I\rho }{\pi }\ln \left(\frac{1}{\left|r\right|}\right)+c$ (15)

式中c为常数，r为电位场节点与源点的距离。

通过无单元Galerkin法的计算与本身理论的计算，得到表1理论值与数值解的计算对比结果，从对比结果中能够直观的看出，其绝对误差相差精度为10⁻² V，当距离源点为1 m的时候，理论值与数值解的误差很大，为0.09 V。当距离源点为10 m的时候，理论值与数值解的误差减小到0.01 V，总结其中的规律可知，两者之间的计算误差十分的小，且距离源点越远的范围内其误差越来越小，其拟合效果更好。通过计算结果的比较，无单元Galerkin法能够准确的求解出地下空间的电位分布值，能够用来进行对直流电场的正演模拟研究。

Table 1. Comparison between theoretical potential value and numerical solution of homogeneous half-space model

表1. 均匀半空间模型电位理论值与数值解对比

图3为电流密度的走势图，从图中可以看出，当源点电位为正电流时，越靠近源点时，其电流密度越大，越远离源点时，其电流密度越小，当其扩散到边界时，其电荷密度接近为0。并相对于正电荷的源点向外扩散，且其电荷分布均匀，符合均匀半空间场的电荷分布情况。其与电位等值线相互吻合，其走向垂直于电位等值线面，同时其电荷密度呈左右对称分布，两边的扩散形式相同。

Figure 3. Current density diagram of uniform half-space model

图3. 均匀半空间模型电流密度图

3.2. 含异常体的均匀半空间地电模型

为了使无单元Galerkin法电流场模拟合理的运用到对有异常体的地下构造模型研究中，并且能够准确的找出异常体位置，分析其异常体对电位的传导与电荷分布的影响，为电流场反演推导地下物性结构奠定基础。本部分设计了两个异常体模型，两个异常体分别在相同位置具有高于或低于背景电阻率的电阻率值。然后通过对其电位分布与电荷分布的正演模拟，为电流场反演奠定坚实的基础。在本模型中，背景电阻率值与网格的布置将沿用上一小节的设置。不同的是，在平行于x轴走向的−10 m至−5 m，垂直于地下深度的z轴−5 m至−10 m处布置了一个5*5 m的正方体异常块。异常体将设置成高阻或低阻，其中高阻异常体为1000 Ω∙m，低阻异常体为10 Ω∙m。其异常体位置模型布置如图4所示。其中，电极电源位置为(0,0)和网格尺寸大小为50*25、电流值为1 A。

Figure 4. Abnormal body position in half space

图4. 半空间异常体位置图

图5、图6分别为高阻异常体、低阻异常体的电位分布图。从图5可以看出，电阻率值高的异常体其电位等势线越密集，电位差变小，且在矩形异常体位置处的电位等值线出现弯曲，向靠近源点处弯曲，其电位的大小急剧减小，说明高电阻率的异常体将减弱电场的传导，加剧电位的衰减；反之，从图6可以看出，电阻率值低于背景场值的异常体处，电位等势线稀疏，电位差变大，且在矩形异常体位置处的电位等值线出现弯曲，向远离源点处弯曲，其电位的大小缓慢变换，说明低电阻率的异常体将增强电场的传导，减缓电位的衰减。

Figure 5. Isogram of potential of abnormal body with high resistance in half space

图5. 半空间高阻异常体电位等值线图

Figure 6. Isogram of half-space low resistance abnormal body potential

图6. 半空间低阻异常体电位等值线图

图7、图8分别为其电流密度分布情况。从图7可以看出，当异常体为高阻时，电流密度变小，排斥电荷，在异常体位置处的电荷密度远小于其他部分电荷，通过电荷密度分布可以看出，当电阻率异常体越高其电荷聚集越少，且根据电荷分布走向可以看出，高阻率异常体将影响电荷的扩散方向，在异常体周边的电荷将绕过异常体继续向外出扩散，并在其周边的电荷将比其它均匀空间内的电荷聚集得多；从图8可以看出，异常体为低阻时，电流密度变大，吸引电荷，在异常体位置电荷密度远大于其他部分电荷，通过电荷密度分布可以看出，当电阻率异常体越小其电荷聚集越多，且根据电荷分布走向可以看出，低阻率异常体将影响电荷的扩散方向，在异常体周边的电荷将穿过异常体继续向外出扩散，并在异常体内部的电荷将比其它均匀空间内的电荷聚集得多。更加准确的显示出了异常体的位置及异常体的电阻率值与背景场的关系，能够很好的反应地下异常体情况。因此本方法适用于地下空间电流场模型的正演研究，能够为生产实际提供有力的论证。

Figure 7. Current density diagram of high resistance abnormal body in half space

图7. 半空间高阻异常体电流密度图

Figure 8. Current density diagram of low-resistance abnormal body in half space

图8. 半空间低阻异常体电流密度图

4. 结论

本文采用新兴的数值方法(无单元Galerkin法)对单极线线源二维直流电场的公式推导，然后对不同的均匀半空间模型进行正演模拟，分析其电位分布及电荷密度，得出了以下结论：

1) 采用无单元Galerkin法对线源直流电法进行正演模拟，只采用节点信息，摆脱了网格的束缚，因其解高次连续的特点，提高了计算精度。同时，采用第一类边界条件，提高了计算效率。

2) 通过对均匀地下半空间的地电模型进行正演研究并与理论值进行对比分析，验证了该方法的有效性及正确性。

3) 建立了含有高阻、低阻异常体的均匀半空间地电模型，通过分析其电位分布及电荷密度，能够准确的反应出地下异常体位置，对地下构造分析有很好的效果，为生产实际提供了理论依据。