在野外实际地质情况中，由于岩、矿石的组成、温度、湿度的变化，地球介质的电导率分布往往是连续变化的(徐世浙等，1995；李予国等，1996).阮百尧、徐世浙(1998)给出了矩形网格单元、电导率连续变化2.5维直流电阻率数值模拟方法，从而有效克服了Coggon(1971)、Dey和Morrison(1979)电导率分布均匀的缺陷.为了易于模拟起伏地形和地下复杂地电体的分布，阮百尧(2001)给出了三角网格单元、电导率连续变化2.5维直流电阻率数值模拟方法.本文在此基础上，导出了四边形内三角网格单元剖分，电导率连续变化2.5维直流电阻率数值模拟方法，并对远边界条件做了进一步的简化处理，重新给出了一组Fourier反变换系数，采用分离式Cholesky分解法快速求解线性方程组.通过对一维连续介质解析法理论模型、土槽实验模型、倾斜界面异常体模型的计算与对比分析，本文方法有较高的模拟精度和计算速度，能够有效模拟起伏地形和复杂形态的地电体模型.

1 波数域点源二维直流电场的变分问题

如图 1所示，假设三维地球介质模型的电导率沿y方向(构造方向)上无变化，三维点源直流电场总电位U沿y方向做Fourier变换至波数域，这时可将三维点源直流电场的边值问题转化为2.5维直流电场的边值问题来进行处理.此时，波数域点源二维直流电位 u 所对应的变分问题(徐世浙，1994)为

图 1

Fig. 1

图 1 研究区域Ω有限元网格剖分 Fig. 1 Division of the study region Ω with FEM grids

采用有限单元法求解泛函式(1)的变分问题，具体步骤如下内容.

2.1 网格剖分

如图 1所示，考虑到复杂地球介质体形状的任意性，将研究区域Ω用三角单元网格进行剖分.首先对研究区域进行水平地形情况下的四边形网格剖分(熊彬和阮百尧，2002；潘克家等，2012)，然后在每一个四边形网格内剖分为4个三角单元网格(徐世浙，1994；刘云和王绪本，2012).当模拟起伏地形时，地面网格沿实际地形线进行网格剖分，地下网格采用逐渐放平进行处理，这样得到的起伏网格更加符合复杂地形情况下直流电场的分布规律(吕玉增和阮百尧，2006；吕玉增，2008).纵向起伏网格的坐标按照如下方式进行计算得到：z_i，j=Z_i，j+(D－Z_i，j)E_i/D(i=1，2，…，n;j=1，2，…，m).其中，E_i是地面上第i个测点的地形高程值，D为水平地形情况下纵向网格(z方向)的最大深度值，Z_i，j是水平地形情况下第i个测点、第j个纵向网格节点的z坐标，z_i，j是起伏地形情况下第i个测点、第j个纵向网格节点的z坐标，n为测点个数，m为纵向网格节点个数.采用这样的网格剖分方法，一方面避免了三角单元网格可能出现过于尖锐情况，另一方面则可以利用三角形的斜边去模拟起伏地形线和复杂介质体的倾斜界面(陈小斌，1999；刘云，2012).

图 2

Fig. 2

图 2 三角单元e Fig. 2 Triangle element e

如图 2所示，在三角单元e内，电位 u和电导率σ 设计为线性变化，则在每一个三角单元内有

当三角单元e的斜边12 位于Γ₁边界上时，有σ=N_iσ_i.这里，u_i为三角单元节点处的电位，σ_i为三角单元节点处的电导率，下标(i=1，2，3)为三角单元的设定节点号. N_i= 1(a_ix+b_iz+c_i)/ 2Δ 是关于x和z的线性插值形函数(i=1，2，3).其中，Δ= 1(a₁b₂－a₂b₁)/2 是三角单元的面积.

a₁=z₂－z₃，a₂=z₃－z₁，a₃=z₁－z₂，b₁=x₃－x₂，

b₂=x₁－x₃，b₃=x₂－x₁，c₁=x₂z₃－x₃z₂，

c₂=x₃z₁－x₁z₃，c₃=x₁z₂－x₂z₁，

2.3 单元分析

将泛函式(1)在整个研究区域Ω中离散化，积分区域可以表示为所有三角单元e的线性组合，即

单元分析1 u 对x求偏导数，则有

其中，(∂N /∂x)^T=(∂N₁/∂x ∂N₂ /∂x ∂N₃/∂x)= 1(a₁ a₂ a₃)/2Δ，u ^T_e=(u₁ u₂ u₃).

显然，(∂u/∂x)²= u ^T_e(∂N/∂x)(∂N/∂x)^T u _e.同理可得，(∂ u /∂z)²= u ^T_e(∂N/∂z)(∂N/∂z)^T u _e.

因此，对于式(2)中第一个积分项(徐世浙，1994；阮百尧，2001；刘云，2012)，有

单元分析2 对于式(2)中第二个积分项(徐世浙，1994；阮百尧，2001；刘云，2012)，有

单元分析3 对于式(2)中的第三个积分式(徐世浙，1994；阮百尧，2001)，根据δ函数的性质可知，只有点电源A(x_A，z_A)对该积分式才有贡献.因此有

单元分析4 对于式(2)中的第四个积分项只对远边界Γ_∞进行线积分(阮百尧，2001；徐世浙，1994).假设，当三角单元的边落在Γ₁边界上时，在波数域k中，kK₁(kr_A)cos(r_A，n)/K₀(kr_A)为一个常数值C，可提至积分式之外.因此，线积分

2.4 系数矩阵总体合成及求变分

对每个单元系数矩阵 k _e按照总体节点号进行扩展，得到，并相加得到总体系数矩阵 K，即

其中，P =(0 … 1/2 I … 0)^T是与点电源有关的列向量.系数矩阵 K 的数学表现形式为 K =(i，j=1，…，ND)，NE为总体单元数，ND为总体节点号，e为单元号，i、j为总体节点号变量.显然，K 是一个对称、正定的大型稀疏矩阵，可按照变带宽方式进行存储(刘德贵等，1980；阮百尧和熊彬，2000).令式(7)的泛函F(u)的变分为零(徐世浙，1994；刘云，2012)，便可得到线性方程组

2.5 计算加速的处理技术

由于远边界Γ_∞离点电源较远，即r_A→∞.此时有，K₁(kr_A)/K₀(kr_A)→1，u→0.因此，波数域点源二维直流电位u所满足的远边界条件(黄俊革，2003；黄俊革等，2003；强建科，2006)为 ∂u/∂n =－ kK₁(kr_A)cos(r_A，n)u/K₀(kr_A)→0.这时，可将远边界条件从第三类边界条件简化成为第二类齐次边界条件∂u/∂n=0，简化后的第二类齐次边界条件属于有限单元法的自然边界条件，在变分方程的推导过程中可以自动满足(徐世浙，1994).对照式(7)中 K系数矩阵，这时可以简化为 K =.此时，得到的 K 矩阵中的元素值只与地电模型参数和波数有关系，而与点电源没有关系.参考附录程序所示，若在求解方程式(8)的过程中，可以先使用Cholesky分解法将 K系数矩阵进行分解，并保留分解结果.随着供电电极的每一次移动，右侧列向量P中点电源的存储位置也相应发生变化，这时只需要将每一次的不同列向量P 通过顺代和回代就可以实现快速求解(阮百尧，2001；宋滔，2013).

图 3

Fig. 3

图 3 四边形内三角网格剖分方式 Fig. 3 The triangle division method in the quadrangle

如图 3所示，本文所采用的四边形内三角网格剖分方式，比较单一的四边形网格剖分来说，其中增加了一个中间节点5.但是不难发现，节点5只与周围的四个网格节点有关系，而与其它的网格节点没有任何直接联系.可以在单元分析中，运用消元法将这个节点消去(徐世浙，1994).因此，节点5是一个虚节点，并不包含在总体节点数中参与计算.等到解方程结束后，也可以根据周围四个角点的已知电位值和网格单元系数之间的关系，将节点5的电位值直接计算出来.这样，一方面既不增加节点的总数，节省了计算量，另一方面在每个三角单元内仍具有各自的不同物性.具体做法如下(徐世浙，1994；刘云，2012)：

节点5的消除，将k_i，j(i，j=1，2，3，4，5)替换为k′ _i，j，则有

k′ _i，j=k_i，j－ k_i，5k_5，j/k_5，5(i，j=1，2，3，4)，

3 Fourier反变换及视电阻率的计算

3.1 Fourier反变换

在以上内容中，可以求出若干波数域k_i中点源二维直流电场的电位u(k_i)值，然后采用Fourier反变换，就可以求解出三维点源空间的总电位U.其中，Fourier反变换的计算公式为

表 1(Table 1)

表 1 Fourier反变换的ki和gi Table 1 ESR dating results of Ge signal in volcanic quartz from Xiaoshan site in Cangzhou，Hebei Province 1 2 3 4 5 ki 3.9928588E-03 3.6820836E-02 0.1434751 0.4675921 1.554464 gi 8.4590055E-03 3.6928654E-02 0.1119829 0.3464288 1.275571

表 1 Fourier反变换的ki和gi Table 1 ESR dating results of Ge signal in volcanic quartz from Xiaoshan site in Cangzhou，Hebei Province

二极测量装置的视电阻率表示为(阮百尧，2001)ρ_a(A，M)=K(A，M)U(A，M)/I，三极测量装置的视电阻率表示为ρ_a(A，M，N)=(K(A，M，N)U(A，M)－U(A，N))/I，四极测量装置的视电阻率表示ρ_a(A，B，M，N)=(K(A，B，M，N)U(A，M)－U(A，N)－U(B，M)+U(B，N))/I .其中，K为装置系数(程志平，2007)，A、B为供电电极，M、N为测量电极.供电电流强度I为常数值，一般可设为1A.

4 模型计算

为验证本文的数值模拟方法，首先与一维连续介质模型的解析法，以及与土槽实验地形模型做计算与对比分析，之后用该算法对复杂界面异常体模型做计算分析.

在对以下所有模型的数值模拟中，研究区域的网格剖分方法具体如下：在横向网格剖分中，目标研究区域采用等间隔网格，网格系数为1.左、右两边扩展的稀疏网格数各为15个，网格系数依次分别为1，1，1，1，2，2，2，2，4，4，4，8，16，32，64.横向网格的实际大小等于点距乘以网格系数.在纵向网格剖分中，网格系数为20个，依次分别为1，1，1，1，1，1，1，1，2，2，2，2，2，4，4，4，8，16，32，64.纵向网格的实际大小等于点距乘以网格系数.值得注意的是，在4.1节中的Fourier反变换系数的选择，就是按照这一组纵向网格系数进行的优化选取得到的.当模拟起伏地形时，起伏地形网格则按照第3.1节方法进行剖分.

4.1 一维连续介质模型

如图 4所示，为一维连续介质地电模型，Mallick和Jain(1979)曾用解析法计算过这个模型.顶层电阻率ρ₁为50 Ωm，层厚度H₁为2 m.底 层电阻率ρ₃为1000 Ωm.中间层为连续介质层，层厚度H₂为10 m，中间层连续介质电阻率ρ₂按线性变化规律从ρ₁变化到ρ₃.测量装置为对称四极装置.

图 4

Fig. 4

图 4 一维连续介质地电模型 Fig. 4 The 1-D geoelectric model of continuous medium

有限元数值模拟中，电极数为200个，最大移动间隔为99，取测线中心处第100号点测点为研究对象.则解析法与数值模拟的结果对比如图 5所示.横坐标为极距(AB/2，单位为m)，纵坐标为视电阻率(ρ_a，单位为Ωm).其中黑色虚线曲线(Analytical)为解析法的计算结果(Mallick and Jain， 1979)，蓝色叉点曲线(Numerical 1)为常规方法的带第三类边界条件数值模拟结果，红色圈点曲线(Numerical 2)为简化后的带第二类齐次边界条件数值模拟结果.

图 5

Fig. 5

图 5 有限元数值模拟与解析法结果比较 Fig. 5 Comparison between analytical and numerical solutions by FEM

由图可见，曲线首支除有2个测点由于受电源点场源影响，有较大的相对误差(最大相对误差为7.1％)，其余各测点拟合都很好，两条数值模拟曲线和解析法曲线的均方根误差小于0.2％.曲线Numerical 1计算耗时为12.7 s，曲线Numerical 2计算耗时为2.5 s. 计算结果表明：一方面，本文方法能够有效地模拟连续介质模型.另一方面，在有较高模拟精度的前提下，采用了远边界条件的简化处理技术和消去网格虚节点的处理技术后，大大提高了计算速度.

4.2 地形模型

如图 6下部所示为二维山脊地形模型，周熙襄和钟本善(1986)对这个模型做过土槽物理模拟实验.山坡的倾角为45°，电阻率为1 Ωm的均匀介质，山顶高0.3 m，剖面长2 m，测量装置为二极测深装置，点距为0.1 m，极距(AO)为0.1 m，电极总数为20个.

图 6

Fig. 6

图 6 二维山脊地形模型与数值模拟结果 Fig. 6 The 2-D chine terrain model and result of the numerical modeling

土槽物理模拟实验和本文数值模拟结果如图 6上部曲线图所示.其中，横坐标为测点号，纵坐标为视电阻率(ρ_a，单位为Ωm)，黑色虚线曲线(Experiment)为土槽物理模拟实验结果(周熙襄和钟本善，1986)，红色圈点曲线(Numerical)为本文数值模拟结果.由于受到二维物理模型几何尺寸的局限性，山峰顶点处测点的相对误差较大，为14.2％，其余各测点拟合都很好，数值模拟曲线和土槽物理模拟曲线的均方根误差为1.5％.计算结果表明，本文采用的三角网格模拟起伏地形、电导率连续变化的数值模拟技术更加符合野外实际地质情况.

4.3 复杂界面异常体模型

如图 7所示，模型设计为均匀介质半空间中放置一“W”字形低阻异常体模型.背景电阻率为100 Ωm，异常体电阻率为20 Ωm.异常体上顶面距地面20 m，下底面距地面60 m.异常体上顶面外宽180 m，下底面外宽90 m.直流电阻率法的测量装置为偶极—偶极测深装置，相邻电极间隔为10 m，发射偶极大小为1，接收偶极大小为1，初始隔离系数为1，测道数为30，排列数为60.

图 7

Fig. 7

图 7 复杂界面异常体模型 Fig. 7 The complex interface abnormity body model

通过计算，测线上电极总数为92个，则有限元数值模拟网格剖分为：网格数为121×20(横向网格数为121，左、右稀疏网格数各为15；纵向网格数为20).大型稀疏矩阵的压缩存储方式按变带宽方式进行存储(刘德贵等，1980；阮百尧和熊彬，2000).参考附录程序所示，采用分离式Cholesky分解法求解线性方程组(阮百尧和熊彬，2000；宋滔，2013).

如图 8所示，为本文数值模拟方法得到的偶极—偶极测量装置视电阻率拟断面图，横坐标为测点位置(Distance，单位为m)，纵坐标为视深度(Depth，单位为m).计算耗时为0.5 s.从图可见，“W”字形异常体的3个上顶面处所对应的3个“八”字低阻异常形态清晰可见，表明本文采用三角网格模拟复杂异常体的方法是有效的.

图 8

Fig. 8

图 8 偶极—偶极测量装置视电阻率拟断面图 Fig. 8 Pseudo section of apparent resistivity on dipo1e-dipole survey array

本文给出四边形内三角单元网格剖分、电导率连续变化2.5维直流电阻率有限元数值模拟方法，比较常规的四边形或三角单元剖分、电导率分布均匀的模拟方法来说：(1)对起伏网格采用了逐渐放平的处理，三角单元内电位值、电导率连续变化的数值模拟技术，一方面更加符合复杂地形情况下直流电场的分布规律，另一方面又能很好模拟野外实际地形和复杂地电体模型，提高了模拟精度.(2)对远边界条件采用第二类齐次边界条件，以及采用分离式Cholesky分解法快速求解线性方程组，使计算速度提高了5倍.通过对一维连续介质模型的模拟计算，以及对土槽地形模型的模拟计算，验证了数值模拟算法的可靠性.对复杂异常体模型进行模拟，结果表明本方法是有效的.本文方法考虑的是2.5维模型情况，三维电导率连续变化的数值模拟方法以及快速计算问题将是本文的后续研究工作.

附录 分离式Cholesky分解法求解线性方程组Fortran语言程序
1. 子程序 SUBROUTINE CHOLESKY_DEC(N,NP,ID,A)
(1)哑元说明
N 整变量，输入参数，方程组的阶数或解向量个数.
NP 整变量，输入参数，为下三角矩阵变带宽内的元素个数.
A 输入、输出参数，NP个元素的一维实数组.输入时存放下三角矩阵变带宽内的各个元素，输出时存放Cholesky分解后的各个元素.
ID 输入参数,N个元素的一维整数组.存放下三角矩阵中的对角线元素在A中的位置.
(2)程序
SUBROUTINE CHOLESKY_DEC(N,NP,ID,A)
INTEGER N,NP
INTEGER ID(N)
REAL A(NP)
I0=1
DO 1 I=1,N
IM=ID(I)
IF(I==1) GOTO 4
JM=I-IM+I0-1
IF(JM .LT. 1) THEN
J0=1
ELSE
J0=ID(JM)+1
END IF
DO 2 JT=I0,IM-1
J=I-IM+JT
JM=ID(J)
DO 3 KT=I0,JT-1
K=I-IM+KT
JK=JM-(J-K)
IF(JK .LT. J0) GOTO 3
A(JT)=A(JT)-A(KT)*A(JK)
3 CONTINUE
J0=JM+1
A(JT)=A(JT)/A(JM)
2 CONTINUE
DO 5 K=I0,IM-1
A(IM)=A(IM)-A(K)*A(K)
5 CONTINUE
4 A(IM)=SQRT(A(IM))
I0=IM+1
1 CONTINUE
RETURN
END
2. 子程序 SUBROUTINE SOLVE(N,NP,ID,A,B)
(1)哑元说明
N,NP,A,ID 同子程序CHOLESKY_DEC.
B 输入、输出参数， N个元素的一维实数组.输入时存放方程组右端的N维列向量，输出时存放解向量.
(2)程序
SUBROUTINE SOLVE(N,NP,ID,A,B)
INTEGER N,NP
INTEGER ID(N)
REAL R
REAL (NP),B(N)
R=DOT_PRODUCT(B,B)
IF(R==0.)THEN
B=0.
RETURN
END IF
I0=1
DO 1 I=1,N
IM=ID(I)
J=1
DO 2 KT=IM-1,I0,-1
B(I)=B(I)-A(KT)*B(I-J)
J=J+1
2 CONTINUE
B(I)=B(I)/A(IM)
I0=IM+1
1 CONTINUE
IN=ID(N)
B(N)=B(N)/A(IN)
DO 3 I=N,2,-1
NN=ID(I)-ID(I-1)-1
DO 4 J=1,NN
IJ=I-J
B(IJ)=B(IJ)-A(ID(I)-J )*B(I)
4 CONTINUE
B(I-1)=B(I-1)/A(ID(I-1))
3 CONTINUE
RETURN
END