0 引言

勘探地球物理存在众多方法，如磁法、重力、地震、直流电及电磁勘探等(Dai et al., 2016; Ren et al., 2017a, 2017b, 2017c, 2018b; Chou and Chouteau, 2016; Attias et al., 2016), 人们常根据勘探需求选择不同的方法.直流电阻率法(Maillet, 1947)是一种可以有效探测地表浅部的电导率结构的地球物理探测方法，它利用地表观测的电势分布，反演地下结构的电导率分布(Pain et al., 2003; Greenhalgh et al., 2009).随着野外采集仪器的快速发展和数值算法及其理论的快速提高，高精度、高效地求解带地形复杂三维直流电阻率模型是目前研究的热点和难点.考虑地下介质各向异性会带来更多的计算参数(Yin, 2000)，因此目前的正反演算法通常只考虑地下介质的电各向同性或者简化的电导率各向异性(Pain et al., 2003; Linde and Pedersen, 2004).然而即使最后反演得到的模型与野外数据完全吻合，也会隐含不能忽视的错误(Asten, 1974; Yin and Weidelt, 1999; Linde and Pedersen, 2004).当遇到具有层理面的岩石或者裂缝等具有方向依赖性(Greenhalgh, 2009)的物质结构时，直流电阻率模拟必须考虑地下电导率各向异性.

正演是反演的关键，首先，必须开发考虑地下电导率各向异性的高精度直流电阻率正演模拟算法.对于一些简单的各向异性模型，有学者计算了均匀半空间倾斜各向异性介质的解析解(Habberjam, 1975)以及各向异性层状介质的地电响应(Yin and Weidelt, 1999; Yin and Maurer, 2001; Wait, 1990; Li and Uren, 1997).然而对于复杂各向异性电导率模型时，我们只能寻求数值解法.目前主要有以下几种数值解法，分别是积分方程法(Li and Stagnitti, 2004; Li and Uren, 1999, Ren et al., 2014, 任政勇等, 2017d)，有限差分法(Wang and Fang, 2001; Hou et al., 2006)，谱元法(Zhou et al., 2009)和有限元法(Bibby, 1978; Li and Spitzer, 2005; Wang et al., 2013; 严波, 2013).如果考虑带地形复杂几何模型时，基于非结构化网格的有限元法成为最优化的选择，因此本文选择基于非结构化网格的有限元算法来处理任意起伏地形复杂直流电阻率问题(Li and Pek, 2008; Wang et al., 2013).

众所周知，源电极附近快速变化的电位导致源附近电位较大的数值误差.为了提高源附近场的精度，一般采用二次场的求解策略(Lowry et al., 1989).传统的二次场求解方法要求一次场满足背景模型的边值问题(Li and Spitzer, 2002; Blome et al., 2009; 徐世浙等, 1994).然而背景场带起伏地形时，一次场不存在解析解只能通过高阶有限元近似或其他需要付出较大计算成本的数值方法来计算(Rücker et al., 2006).本文基于Penz等(2013)从直流电阻率各向同性模型提出的概念，扩展到各向异性模型，利用均匀半空间各向异性模型的解析解(Li and Uren, 1997; Penz et al., 2013)作为一次场，建立了直流电阻率各向异性问题的“虚拟场边界值问题”，即有效地去除了源奇异性现象，也正确地处理了任意起伏地形情况，从而达到精度的最优化方案.

为了达到速度的最优化方案，我们采用了面向目标的自适应有限元算法来求解上述新的“虚拟场边界值问题”.求解上述“虚拟场边界值问题”，可采用基于先验信息的人工有限元网格加密算法(Li and Spitzer, 2002; Rücker et al., 2006)、基于后验误差估计的全局自适应有限元法(Babuška and Rheinboldt, 1979; Zienkiewicz, 2006; 胡恩球等, 1997; 王建华等, 2000)和面向目标的自适应有限元算法.对于复杂的模型，很难确定场变化剧烈程度等先验信息，从而难以设计出能够真实模拟复杂场变化的有限元网格，全局自适应有限元法(Franke et al, 2007; Ren and Tang, 2010; 任政勇和汤井田, 2009; 王飞燕, 2009)采用全局的后验误差估计算法(Ainsworth and Oden, 1997)，即所有的网格单元误差具有相同的权重，来达到降低全局网格上平均数值误差的目的.对于直流电阻率勘探问题来说，我们既需要降低全局网格的数值误差，更需要测量电极处的数值解达到足够的精度.全局自适应有限元法所导致的全局网格再分配方案往往不能有效提供测点附近的数值解精度，并且存在浪费计算量的情况.面向目标的自适应有限元算法首次由Oden和Prudhomme提出(Oden and Prudhomme, 2001)，通过增加测量区域的误差权重，从而达到大幅度提高测点处数值解精度的目的.由于其优越性，面向目标的自适应有限元算法在勘探地球物理领域得到越来越多的应用.Kerry和Weiss(2006)首先利用面向目标的自适应有限元算法解决了2D大地电磁问题.Li and Kerry(2007)利用面向目标的自适应有限元算法解决了2.5D可控源电磁问题.Ren等(2013)利用面向目标的自适应有限元算法解决了3D大地电磁问题.Ren和Tang(2014)利用面向目标的自适应有限元算法解决了3D直流电阻率各向同性问题.对于直流电阻率各向异性问题，目前最新的研究成果为Wang等(2013)利用基于非结构化网格的有限元算法，因此，本文提出的面向目标自适应有限元算法可被看作高精度求解直流电阻率各向异性问题的进一步延伸.

另外，本文的另一目的为测试基于电流密度法向分量的不连续性的后验误差估计算法求解稳态直流电阻率法的能力.后验误差估计算法是指导面向目标的自适应有限元网格自适应修正的关键.在电阻率不连续分界面上，电流密度法向分量呈连续状态，由于网格的不合理性，有限元算法计算出的电流密度法向分量在电阻率不连续分界面上一般呈不连续状态.在前期的大地电磁问题研究中，我们发现基于这一物理现象构建出的后验误差估计算法具有非常优良的性能，能够有效计算出有限元解的数值误差，从而驱动网格的自适应加密.这一研究发现，也被Yin等(2016)成功应用了航空电磁法的高精度计算，另外这一后验误差与其他种类后验误差的对比见我们最新的研究成果(Ren et al., 2018a).

本文从考虑电导率各向异性情况出发，建立了直流电阻率各向异性问题的“虚拟场边值问题”，利用基于电流密度法向分量的不连续性的后验误差估计算法，提出了可以大幅度提高数值解精度、大幅度减小计算消耗的面向目标的自适应有限元算法，利用非结构化网格和高效的并行线性方程组求解器Pardiso，成功地实现了带起伏地形复杂各向异性直流电阻率三维问题高精度、低计算消耗的求解目的.最后，三组电阻率各向异性模型验证了本文算法的正确性和适应性.

1 理论 1.1 电导率张量

对于任意各向异性的地下介质，电导率σ是一个具有6个独立分量的张量，其为实对称正定矩阵(Yin, 2000)，可以通过三次旋转变换转化为正对角矩阵.如图 1所示，x′, y′, z′表示主轴坐标系，对应岩石的自然结构，x, y, z为测量坐标系，表示测点的地理坐标.由测量坐标系变换到主轴坐标系共需要确定三个欧拉角(Yin, 2000；Greenhalgh, 2009；Zhou et al., 2009)，首先固定x轴，以角α顺时针旋转平面(y, z)，变换矩阵为；然后围绕上一步旋转之后的y轴以角β旋转平面(x, z)使得z轴与电导率主轴z′重合，相应的旋转矩阵为；最后固定上一步得到的z轴(即z′轴)以角γ旋转平面(x, y)使得x, y轴分别与x′, y′轴重合，可得变换矩阵.综上可得由主轴坐标系到测量坐标系的总变换旋转矩阵R=(R_zR_yR_x)^-1，即：

(1)

由于电导率各向异性的存在，向地下供入电流之后电流密度的方向不再与电场强度的方向一致，但是在主轴坐标系下电流密度J′和电场强度E′方向相同：J′=σ′·E′.测量坐标系下电流密度J和电场强度E都可以通过主轴坐标系下的电流密度和电场强度变换得到：J=R·J′, E=R·E′.根据电流密度和电场强度的关系：J=σ·E可以得到电导率张量的表达式：

(2)

式中，对角线上的元素表示三个主轴电导率，当三者相同时地下介质呈现各向同性.三个旋转角和电导率主轴的组合可以唯一确定地下各向异性的电导率结构.对于典型的TI(transversely isotropic)介质，我们假定(y, z)平面为岩石层理面，x轴表示垂直层理面的对称轴，则σ₂=σ₃表示沿层理面的纵向电导率，σ₁表示垂直层理面方向的横向电导率(李金铭，2005).当旋转角β=γ=0时，角α反映地层倾向的变化，当α=90°时表示垂直各向异性介质；当旋转角α=β=0时，此时地层为水平各向异性介质，旋转角γ表示地层水平各向异性的方位.

1.2 基于二次虚拟场的电位边值问题

如图 2所示的模型，Γ₀是地表面，Γ_∞是无穷远边界，围绕封闭区域Ω.{S_i}(i=1, 2, …, s)和{M_i}(i=1, 2, …, m)分别表示源电极群和测量电极群，在地表面Γ₀上共布置s个源电极和m个测量电极，源电极S_i的坐标是p_i，点p_o是源电极群{S_i}的中心；p, n分别表示边界上任一点及其对应的外法向向量，σ是地下介质的电导率张量，矢量r_i和r_o分别表示：r_i=p－p_i和r_o=p－p_o.在S_i处注入电流I，其引发电流密度J_i=J(x, y, z)，基于全电流定律(Stratton, 2007)，忽略位移电流，有：

(3)

对方程(3)两边取散度，结合电荷守恒定律的微分形式及稳态场下电场强度和电势的关系E_i=-ΔU_i，我们可以得到在M_i处关于电位U_i的控制方程：

(4)

由于源附近场变化迅速，直接求解总场会产生很大的数值误差，为了避免源的奇异性，我们采用特别的二次场法(Penz et al., 2013)，其中一次场采用各向异性均匀半空间的解析解，我们将总场分解为

(5)

式中U_i^s, U_i^v分别表示一次奇异场和剩余的二次场.在平滑的地表Γ₀，如果S_i的邻域足够小，那么它可视作各向异性的均匀半空间，其电位解析解有以下形式(Li and Uren, 1997)：

(6)

式中B_io=r_i^T·ρ_o·r_i，ρ_o是源电极S_i附近邻域的张量电阻率，如果邻域足够小，那么ρ_o可以假设为一个常张量，在实际计算中，将源电极邻域定义为一系列包含源电极S_i的小区域，常电阻率张量作为该邻域的平均电阻率.此处对于源电极的位置不做任何限制，如果源位于带地形的地表，一次场采用源点处的切向平面往地下延伸构成的均匀半空间作为背景场.将总场减去邻域的一次奇异场，剩余的二次场因为不满足地表的齐次诺依曼条件被命名为二次虚拟场.

与总场类似，我们可以得到源电极S_i产生的奇异场所满足的控制方程：

(7)

将(7)式与(5)式代入方程(4)定义的总场满足的方程，我们可以得到虚拟电势满足的微分方程：

(8)

源电极S_i向地下供电时由于空气无穷大的电阻率，没有电流通过地表流到空气中，所以地表电位的诺依曼条件成立：

(9)

式中n是由地表外法向向量.将(5)式代入(9)式，在地表Γ₀，有

(10)

在无穷远边界Γ_∞上，地形以及不均匀电导体的影响可以忽略，参照各向异性均匀半空间的解析解(Li and Uren, 1997)及类比各向同性的情况(Dey and Morrison, 1979)，假设位于平滑地表的点电源S_i在无穷远边界Γ_∞上产生的总场电位(Li and Uren, 1997)以如下的形式衰减：

(11)

式中B_i=r_i^T·ρ·r_i，C是一个常数.结合总场的分解方程(5)和一次场的奇异解析式(6)，我们得到位于边界Γ_∞上关于虚拟场的混合边界条件：

(12)

综上所述，关于源电极S_i，对于二次虚拟场的边值问题如下所示：

(13)

基于虚拟场的微分方程，由伽里金有限单元法(Zienkiewicz and Taylor, 2000)，可得：

(14)

式中T是一个测试函数，且T属于希尔伯特函数空间H¹(Ω)(Brenner and Scott, 2007)，该空间中函数及其梯度都是有限的，应用格林恒等式再代入边界条件(13)式，将虚拟场在各个单元上应用线性形函数插值(Jin, 2014)，采用形函数作为测试函数T，且当点p位于无穷远边界Γ_∞上，即距点源S_i足够远时，r_i≈r_o，这样便可忽略单个源的特殊性，从而形成适用于多源系统的线性方程组(Ren and Tang, 2014)：

(15)

式中K是所有源共用的系数矩阵，F_i是只与源S_i有关的右端项，他们都包含各个单元的面积分和体积分，我们使用开源代码(Si, 2015)将模型域剖分成Delaunay四面体网格，采用高斯积分规则(Zhu et al., 2005)计算各个积分项，最后使用基于Intel MKL的精确高效的直接LU求解器PARDISO(Schenk and Gärtner, 2004)求解大型稀疏矩阵线性方程组.K和F_i具体表达式如下：

(16)

式中B_i=r_i^T·ρ·r_i & B_o=r_o^T·ρ·r_o.

1.3 面向目标的后验误差估计技术

有限元初始网格并不能保证其数值解达到足够的精度，我们需要估计单元基本误差并对初始网格进行修正即自适应计算.目前自适应有限元方法中主要存在两类后验误差估计技术被广泛采用，分别是基于梯度恢复技术的后验误差估计方法和基于残差的后验误差估计方法(Ainsworth and Oden, 1997).残差型的后验误差估计技术利用局部区域的残值(Ainsworth and Oden, 1993)来估计后验误差，这里采用由于数值误差引起的边界条件的不连续性来估计单元的基本误差，这些数值误差来源于模型误差，网格离散误差，线性形函数插值误差等因素.在计算区域Ω上我们假定其中一个网格单元为Ω_k，它的一个邻元为Ω_k^ne，我们将Ω_k和Ω_k^ne上的电导率和电势分别用σ, U_i和σ^ne, U_i^ne标记.在他们的共面上，其上电流密度的不连续性可以表示如下(Stratton, 2007)：

(17)

式中[·]表示L2范数操作，n为由单元Ω_k指向单元Ω_k^ne的单位法向量，i表示源S_i，F_j表示单元Ω_k的第j个邻面，据此我们可以计算网格单元的基本残差E_k：

(18)

式中s, m分别指代源的数量、单元Ω_k邻面的个数.

式(18)对单元相对误差的计算属于全局型的误差估计，即所有网格单元的基本数值误差具有相同的误差权重，然而在测量剖面处通常需要获得更高的数值精度，即需要增加测量区域的误差权重.为了解决这一问题，应用面向目标的概念(Oden and Prudhomme, 2001)，可由此导出面向目标的后验误差.

大型线性方程组(15)可以写成双线性形式：

(19)

式中U_i^v, T∈H¹(Ω)，H¹(Ω)表示希尔伯特函数空间(Brenner and Scott, 2007)，b(T, U_i^v)表示未知项，S(T)表示源项，两者表达式如下：

(20)

通过构造虚拟场变分公式(19)的对偶问题可以得到影响函数W(Ren and Tang, 2014)的分布.应用对偶加权误差估计法(Oden and Prudhomme, 2001)，为了构造式(19)的对偶变分，我们定义误差项E(e^U)(Ren and Tang, 2014)：

(21)

式中e_i^U指代源S_i产生的虚拟场有限元数值解的误差，E(e^U)是虚拟电位总的数值误差项的函数，它衡量的是在子域集合Ω_M内所有源产生的虚拟电位误差的平均值，在实际计算中我们将各个测点的邻域统称为子域集Ω_M.将方程(19)中源项S(T)用误差项E(e^U)替换，并基于偶函数b(, )的自我伴随特征(Ren et al., 2013)，我们可以得到下面的变分形式：

(22)

式中W∈H¹(Ω)，W被称为影响函数(Oden and Prudhomme, 2001)，它对产生面向目标的概念起着关键作用，(22)式可以写成矩阵形式：

(23)

所以影响函数W离散解与U_i^v的有限元数值解具有相同的系数矩阵K，它的计算基于和U_i^v相同的网格.

经过换算，影响函数的变分公式(22)对应的微分形式为

(24)

由上式我们可以看出影响函数W具有满足地表齐次诺依曼边界条件的电位分布.对于单个源S_i，如果考虑所有测量电极处的虚拟电位误差项，则是在所有测点所在的邻域中分别注入单位幅度电流I，在整个域Ω所产生的电位分布.将影响函数W分成数值解W^FEM和误差项e^W，基于方程(22)，我们可计算测量剖面总的数值误差Θ为

(25)

基于伽里金正交性(Brenner and Scott, 2007)，忽略双线性形式b(, )(式(20))中的面积分项，有

(26)

式中Θ_k用于衡量每个单元的基本误差，利用柯西不等式(Brenner and Scott, 2007; Ren et al., 2013)，上式变为

(27)

式中Θ_k^U, Θ_k^W分别反映了虚拟场的误差和影响函数的误差，由式(24)我们意识到在测点区域附近影响函数的变化十分迅速，所以相对其他区域，测量剖面由于影响函数产生的误差Θ_k^W也会更大，这表明在考虑了全局数值误差Θ_k^U的基础上，使测量区域的单元基本误差在所有单元中占有更高的权重，从而有目的地加密测量电极处的网格，在后面的模型计算中我们将会佐证这一点.由于梯度误差Δe_i^U和Δe_i^W不容易获得，我们采用电流密度的法向分量残差的L2范数来替代相应的能量范数，即使用式(18)来分别估计虚拟场和影响函数场的电流密度法向分量残差，这种根据场的变化规律设置的后验误差估计将能更好地模拟物性界面.

为了实现特定单元的网格加密，我们定义一个单元相对误差指标(Ren and Tang, 2014)：

(28)

式中Θ_max是所有单元误差中的最大值，标记单元相对误差β_k大于设定的阈值β_s的网格单元，在下一次网格再生中进行加密，可以通过设置更小的阈值β_s来增加每次网格修正的单元数以缩减算法的收敛时间.流程图 3给出了自适应算法(Ren and Tang, 2014)的具体实现.

本文算法的终止条件包含总的单元数值误差阈值、最大允许迭代次数、最多未知数个数，我们通常使用后两者作为程序的终止条件.如果循环终止条件达到要求，获得基于最终网格的虚拟场，加上一次奇异场电位便可得到测量电极处的总场，然后可用于计算不同排列装置的标量视电阻率或张量视电阻率(Bibby, 1986；Bibby and Hohmann, 1993).

2 算例 2.1 代码验证

本文算例的运行平台为：Intel(R) Core(TM) i5-4590 CPU @3.30GHz(4 cores)和8GB RAM.我们首先采用各向异性的两层模型(图 4)来验证代码的精度，两层介质具有相同的水平各向异性，源电极位于坐标原点，沿x轴正方向以间隔为0.5 m的间距布设20个测量电极，σ₁表示第一层介质电导率张量的三个主轴电导率及三个旋转欧拉角：σ₁:ρ₁=100 Ωm, ρ₂=ρ₃=25 Ωm, α=0°, β=0°, γ=90°，σ₂表示第二层介质的有关参数：σ₂:ρ₁=10 Ωm, ρ₂=ρ₃=2.5 Ωm，α=0°, β=0°, γ=90°.采用pole-pole装置，参照解析解(Wait, 1990)，我们计算得到了测量剖面的视电阻率和相对误差曲线图，如图 5所示，本文算法计算的结果与解析解有极好的一致性，计算节点数为28373，最大误差小于0.556%.

2.2 带地形各向异性模型适应性

下面我们将使用山脉峡谷模型证明本算法对于带地形各向异性模型的有效性，如图 6所示，该模型为均匀半空间中嵌入一个带地形的倾斜各向异性(Li and Spitzer, 2005)的异常体，其中山脉和峡谷的最大地形都为10 m且关于顶点对称，我们将源电极布于A处，沿x轴正方向以间距10 m布设19个测量电极经过山顶和谷底，异常体的电阻率为倾斜各向异性，电导率相关参数见图 6，均匀半空间的背景电阻率为100 Ωm.

图 7测试了面向目标自适应算法的收敛性，计算了其三次网格的视电阻率和相对误差曲线，因为该模型没有解析解，所以采用全局网格同等加密的方式加密至网格节点数达到11037605的数值解作为拟解析解(Ren et al., 2013).图 7可以看出视电阻率曲线向准解析解迭代收敛.表 1是三次网格水平的单元参数及统计误差，经过四次迭代，视电阻率的最大相对误差从8.72%极大地减少到了0.82%，表明本文提出的算法可以有效处理地形，非结构化网格可以很好地自适应拟合地形的变化.

图 8(a)、(b)是初始网格和第五次网格的相对单元误差β_k的网格等值线图，具体视角是做y=50 m的切片后将(y, z)平面旋转一定角度.图 8a中源电极附近单元相对误差值较大，表明该区域需要加密，所以在图 8b中得到了细密的剖分，从而消除了源附近的奇异性.图 8b中山脉和峡谷处以及电导率的分界面上都具有稠密的网格密度，而在远离异常体的区域则网格较为稀疏.山脉峡谷地形起伏大需要极小的网格单元才可以精确地拟合，电导率分界面上由于急剧变化的场导致电流密度法向分量的数值误差相对更大，因此这些区域会得到重点加密，其他区域较为稀疏的网格密度下线性插值已经可以模拟出场的变化趋势.此外，图中测点周围得到了细密的网格剖分，这将极大地改善数值解.因为网格单元加密到一定的程度时，测量剖面的精度严重依赖于局部的网格密度(Brenner and Scott, 2007)，所以面向目标的算法可以通过加密测点局部区域从而实现高精度的收敛性.

图 8c是由基于第五次网格计算得到的影响函数W的网格密度分布，由于其服从地表齐次诺依曼边界条件的电位分布(见方程(24))，所以在图中测量剖面上测点处出现红色的亮点，这表明了由对偶变分构造的虚拟源的存在，因而可以产生在测点附近着重加密的面向目标效应.

2.3 埋藏多个各向异性的立方块

最后我们简要分析一下地下埋藏多个各向异性异常体的情况.如图 9所示，均匀半空间中嵌入四个埋深和尺寸相同的立方块A、B、C、D，其中心埋深和尺寸均为3 m，为了排除异常体的特征随源取向的影响，采用两对偶极源S1-S2, S3-S4用于张量视电阻率的测量(Bibby, 1986; Bibby and Hohmann, 1993).如图 10所示，首先通过偶极子S1-S2测量测点P处的电场分量E₁₁, E₁₂，然后测量偶极子S3-S4产生的电场分量E₂₁, E₂₂.根据电场与电流密度的关系：E=ρ·J，即：

(29)

式中J₁₁, J₁₂和J₂₁, J₂₂分别是各向同性的均匀半空间中偶极子S1-S2和S3-S4产生的电流密度分布.由式(29)我们可以得到P₂旋转不变量(Bibby and Hohmann, 1993; Li and Spitzer, 2005; Wang et al., 2013)的表达式：

(30)

P₂不变量具有与标量视电阻率相同的量纲，各向同性均匀半空间下两者等效(Wang et al., 2013).

图 11中展示了P₂不变量的地表等值线图，黑色的虚线框是立方块在地表的投影，背景均匀半空间的电阻率为100 Ωm，图中四个立方块异常体的电阻率相关参数分别为：A:10、B:1000、C:ρ₁/ρ₂/ρ₃=10/100/10, α/β/γ=0°/0°/0°、D:ρ₁/ρ₂/ρ₃=10/100/10，α/β/γ=0°/0°/90°.立方块A，B分别是各向同性的低阻和高阻异常体，在图中得到了较为准确的响应，立方块的中心位置和大概范围都得到了很好的刻画，C, D表示方位角γ相差90°的水平各向异性立方块，在图中立方块C和D都得到了一定程度的对称拉伸，两者拉伸的方向与主轴电阻率的大小及方位角γ有关，立方块C拉伸较大的方向为其垂直层理面的方向，相对应的横向电阻率相对于高阻背景电性差异较小，从而使得该方向电流密度较大，而立方块D为C旋转90°的方向，故其拉伸方向也相差90°，这也是各向异性悖论的体现(李金铭，2005).此时立方块C，D的中心位置仍然得到了准确的反映，但在两个低阻的各向异性立方块之间存在相对高阻的假异常.

图 12则是四个不同倾斜各向异性立方块的P₂不变量响应，背景均匀半空间的电阻率为5 Ωm，四个立方块具有相同的主轴电阻率：ρ₁/ρ₂/ρ₃=10/100/10，倾角不同：A:α/β/γ=0°/0°/0°，B:α/β/γ=30°/0°/0°，C:α/β/γ=60°/0°/0°，D:α/β/γ=90°/0°/0°.总的来看，随着倾角增大，各个立方块的等值线拉伸幅度逐渐减少，且逐渐呈现为更加低阻的异常.立方块A倾角为0°，此时对应水平各向异性，在图中表现为对称的拉伸椭圆，由于立方块A的纵向电阻率与低阻背景的电性差异相对较小使得该方向的电流密度分布更为密集，所以其拉伸方向为平行层理面的方向.随着倾角增大至30°，立方块B呈现不对称的拉伸，拉伸方向与A相同，此时异常的中心偏移了立方块的投影中心.倾角继续增大至60°，立方块C也呈现不对称的拉伸，但拉伸幅度有所减少，当倾角达到90°时，D立方块的P₂等值线表现为无拉伸的圆形异常，且其中心与立方块的投影中心吻合，此时倾斜各向异性已经转化为垂直各向异性，在地表水平方向不能检测到各向异性.

由上述复杂模型可知，地下电阻率各向异性介质在直流电阻率法野外观测数据中具有不可忽略的反映，不同的各向异性排列方式会产生不同幅度的地表响应.但是，目前仍缺乏一种从观测数据中判断地下各向异性的有效方案，这一问题需要进一步深入研究.

3 结论

本文提出的基于非结构化网格的面向目标的自适应有限元算法，可以处理带起伏地形的任意复杂三维直流电阻率问题.不同于传统的二次场法，本文的一次场采用各向异性均匀半空间的解析解，减去这个简单的一次奇异背景场，可导出满足带任意地形的二次虚拟场的边值问题，在两层模型中验证了算法的高精度以及在山脉峡谷模型中证明了算法对地形的有效性.

通过估计电流密度法向分量的不连续性作为后验误差，可以很好地指导非结构化网格自动加密以适应物性界面场的急剧变化.结合面向目标的概念，使算法在考虑全局数值误差的基础上着重加密测点区域，可以实现测量剖面更高效地达到收敛精度，极大地节约了计算成本.

在对于产生面向目标特性起到关键作用的影响函数的网格密度分布图中(图 8c)，当测点离源电极的距离逐渐加大时，由对偶变分(Oden and Prudhomme, 2001)构造的虚拟源的特性逐渐减弱，即在远离源的测点处的局部加密能力得到削弱.我们可以考虑关于距离的加权因子，以提升远点的面向目标特性，为大范围的电法勘探提供高精度的数值解，这需要进一步的研究与测试.

致谢

感谢Mark Blome博士构建了最初的山脉峡谷模型.