2021, Vol. 64
直流电阻率法是电法勘探领域一个重要的分支已在环境和工程领域获得广泛应用.该方法利用接地电极向地下供电、建立地下稳定电流场，通过观测测量电极间的电位差以达到推断和解释地下目标体的目的.在野外勘探中我们常常采用剖面或者测深装置进行高密度电阻率测量.为此，需要在地面布设数十至数百个测量电极，通过多芯电缆连接到观测设备，从而实现多种电极装置的分布式测量(Loke and Barker, 1996).由于环境工程领域勘探目标大多为三维地质体，采用这种阵列测量方法可以准确地重构出三维地质体的空间分布，因此高密度电阻率法可获得丰富的地质信息(Loke et al., 2013).
随着高密度电阻率法在环境、工程、地下水资源等勘探领域的广泛应用，三维反演和解释技术研究受到广泛关注(Udphuay et al., 2011; Kneisel et al., 2014; Neyamadpour, 2019; Nthaba et al., 2020).同时，随着电阻率法向三维精细化探测方向发展，能够提取更多的细节信息，为数据解释打下良好的基础.然而，对精细结构的刻画将导致庞大的三维数据体，对这些三维数据体进行反演解释需要高效和高精度的三维正演模拟技术(Günther et al., 2006; Rücker et al., 2006).常用的电阻率三维模拟算法包括积分方程法(IE)、有限差分法(FD)和有限元法(FE)(Wang and Signorelli, 2004; Yoon et al., 2016; Ren et al., 2018a).积分方程法仅对异常体进行剖分，计算效率高，但不适合复杂异常体的数值计算(Sheng et al., 2006; Ueda and Zhdanov, 2006).有限差分算法简单、求解速度快，但由于要求计算域被剖分成规则单元，对复杂模型的适用性较差(Jaysaval et al., 2014; Davydycheva et al., 2003).有限元法由于网格剖分灵活、计算精度高，是目前最主流的方法(Song et al., 2017；谷宇等，2020).传统有限元法仅在单元端点或棱边上产生未知数，计算精度对网格的依赖性较为严重，为达到较高的精度通常需要增加单元内未知数的个数或者网格数量(Ren et al., 2013).对于规则的八叉树网格，通过引入悬挂点对局部区域加密，在控制网格数量的同时提高精度(Grayver and Bürg, 2014).非结构网格及自适应网格加密技术目前已发展得较为成熟，它能够根据后验误差不断迭代得到最优的计算网格(Key and Ovall, 2011).然而，其结果是未知数大幅增加，并且网格的迭代过程也会带来时间损耗.为此，人们尝试通过增加插值基函数阶数来提高数值精度，目前在电法勘探领域大多只用到二阶插值基函数(Chen et al., 2018).
本文研究非结构谱元法(Spectral Element Method, SEM)进行电阻率三维各向异性正演模拟问题.谱元法属于高阶有限元算法的一种.传统高阶有限元法采用等间距插值，由于龙格现象的存在，数值解在靠近端点位置会产生震荡，因此随着基函数阶数的升高，精度不一定会得到提升(Runge, 1901).谱元法采用非线性正交基函数描述场的扩散，在端点附近插值点更加密集，因此能够有效提高数值精度(Maday et al., 1993).谱元法结合了有限元法与谱方法的双重优点.它不仅具有有限元法的网格灵活性，还具有谱方法良好的数值精度和指数收敛性，可以通过细化网格或提高插值基函数阶数达到提升计算精度的目的.Komatitsch等(2002)将谱元法用于全球地震波场的精确模拟，Huang等(2019)验证了谱元法在航空电磁正演模拟中的有效性.然而，所有这些基于谱元法的三维正演均采用六面体网格，这种几何结构清晰的网格限制了谱元法对于地形及复杂异常体的刻画能力，给地球物理探测用于解决复杂工程问题带来困难(Yin et al., 2016b; Wang et al., 2013).相比之下，基于非结构四面体网格的谱元法可以很好地解决这些问题.它除了继承传统谱元法通过细化网格或提高插值函数阶数达到很高的计算精度外，还可利用非结构网格模拟起伏地形和地下复杂构造.Liu等(2014)分别利用三角形网格有限元法和谱元法进行二维时域弹性波场模拟，对比结果验证了谱元法能够达到更高的计算精度.Zhu等(2020)将四面体网格谱元法引入电法勘探领域，在与自适应加密的有限元算法比较中发现谱元法只需求解较少的未知数，就能够达到更高的精度要求.由于四面体网格坐标之间存在的复杂耦合关系，无法使用六面体网格对应的基函数和插值节点(Cohen, 2002)，因此，本文我们将利用Proriol-Koornwinder-Dubiner(PKD)正交多项式(Proriol, 1957; Dubiner, 1991; Koornwinder, 1975) 构建基函数.另外，我们还使用Warp & Blend插值节点集，保证良好的插值属性，改善算法的稳定性(Warburton, 2006).
各向异性是环境与工程问题中经常遇到的地质现象.从微观上，各向异性是由于温度压力等因素造成矿物结构发生变化，形成了页岩和板岩等层理发育的岩石；从宏观上，由于构造运动造成地下介质变质分层或者产生裂隙，这些各向同性的薄层相互叠加或者定向排列的构造裂隙充水后，导致电阻率呈现各向异性特征(Herwanger et al., 2004; Aissaoui et al., 2019；蔡义宇等，2020). 在这些层理或裂隙发育地区利用传统的各向同性模型进行数据解释会产生错误结果(Yin and Weidelt, 1999).因此，利用各向异性模型进行数据反演解释非常必要.虽然针对由介质不均匀性造成的各向异性本质上可以利用三维精细模型进行模拟，然而受计算条件限制，利用不均匀介质模拟这些精细结构难度很大.因此，当这些介质的层理或裂隙尺寸比采用的电法勘探技术的有效尺度小很多时，我们可以利用各向异性描述这种不均匀性，从而大大减少计算量.需要注意的是，各向异性与电性不均匀性在数学描述上存在本质差异.不均匀性指示了岩矿石电阻率随位置发生变化，与方向无关，因此描述各向同性不均匀电导率模型仅需要一个随位置变化的标量函数.然而，各向异性是同一位置上电导率随电流流动方向发生变化，此时电流密度与不同方向的电场存在耦合关系，因此描述各向异性电导率是一个3×3的张量(Yin, 2000).在直流电阻率法各向异性数值模拟方面，Zhou等(2009)使用高斯正交格网(GQG)方法处理复杂各向异性介质模型，能够达到与谱方法近似的收敛速度.Li和Spitzer(2005)结合有限元正演模拟和Cholesky预处理共轭梯度法，研究了各向异性模型的电阻率分布特征.Ren等(2018b)针对各向异性电阻率模型研究了不同的非结构网格自适应加密策略，取得了很好的应用效果.本文针对直流电阻率法验证高阶谱元法模拟各向异性模型的有效性.
在本文后续讨论中，我们首先对各向异性模型进行数学描述，并以此为基础建立任意各向异性介质电阻率法三维正演理论，进而我们利用正交多项式构建谱元法基函数，并结合插值节点与数值积分节点完成有限元线性方程中系数矩阵的建立，最后利用直接求解器实现方程的快速求解.在对本文谱元法进行精度验证后，我们系统分析各向异性对电阻率响应的影响特征.最后，我们重点探讨各项异性特征识别方法以及地形效应和地形校正技术.
1 各向异性模型针对任意各向异性模型，我们用一个对称正定的电导率张量对其进行描述(Yin, 2000)，即：
|
(1) |
其中，σxy=σyx, σxz=σzx, σyz=σzx.电导率张量和电阻率张量之间满足σ=ρ-1.另外，电导率张量σ可由主轴电导率张量经三次欧拉旋转获得，即：
|
(2) |
其中，主轴电导率张量σp为
|
(3) |
而σx、σy、σz为主轴电导率.
式(2)中的总旋转矩阵D可表示为
|
(4) |
其中：
|
(5) |
分别为绕x、y、z轴的旋转矩阵，χ1、χ2和χ3为对应的欧拉旋转角.图 1给出任意各向异性介质的主轴电导率张量旋转示意图.其中，实线和虚线分别表示沿不同坐标轴欧拉旋转前后的介质模型，而阴影面表示同一个面旋转前后的位置.图 1a中(x, y, z)为测量坐标系, 主轴电导率绕x轴旋转χ1之后，坐标系变为(x1, y1, z1)；参见图 1b，(x1, y1, z1)再绕y轴旋转χ2，主轴坐标系变为(x2, y2, z2)；同理，(x2, y2, z2)绕z轴旋转χ3之后，可得最终的主轴坐标系(x3, y3, z3)，如图 1c所示.由此完成了沿三个方向的欧拉旋转，此时三个主轴电导率不再沿测量坐标系方向，电导率张量变为3×3的对称正定矩阵.
|
图 1 任意各向异性旋转示意图(参考Yin et al., 2016a) (a) 绕x轴旋转；(b) 绕y轴旋转；(c) 绕z轴旋转. Fig. 1 Sketch of arbitrary anisotropic rotation(refer to Yin et al., 2016a) (a) Rotation around x axis; (b) Rotation around y axis; (c) Rotation around z axis. |
由上述讨论可知，我们只需要三个主轴电导率σx、σy、σz和三个欧拉角χ1、χ2、χ3即可对地下各向异性介质进行定量描述.参见图 2，按照主轴电导率及定向不同可以将地下介质划分为垂直方向的横向各向同性(Vertically Transverse Isotropy，VTI)、水平方向的横向各向同性(Horizontally Transverse Isotropy，HTI)、倾斜各向异性(Dipping anisotropy)、三轴各向异性(Triple-axis anisotropy)及任意各向异性(Arbitrary anisotropy)(刘云鹤等，2018).VTI介质(图 2a)层理面水平，介质沿层理面的电导率不随方向变化，其电导率张量为一个对角阵，且σx=σy≠σz；而HTI介质(图 2b)的层理面直立，可由VTI介质绕一个水平轴旋转90°得到，其电导率沿直立的层理面同样不随方向变化，电导率张量也为一个对角阵，且有σy=σz≠σx或者σx=σz≠σy；三轴各向异性(图 2c)对应的电导率张量也是对角阵，但是其三个元素不相同σx≠σz≠σy.倾斜各向异性介质(图 2d)是由VTI介质绕水平轴旋转一个任意角度，其电导率张量包含5个分量，其中旋转轴对应的电导率不随旋转发生变化；任意各向异性(图 2e)由VTI介质或三轴各向异性绕多个坐标轴旋转得到，其电导率张量为3×3的满秩矩阵.考虑到地球物理中各向异性大多是由层理或定向排列的裂隙引起的，人们通常引入沿层理方向上电导率为σL及垂直层理方向上为σT(满足σT＜σL)，并定义平均电导率
|
图 2 各向异性模型及电导率张量 (a) VTI介质；(b) HTI介质；(c) 三轴各向异性介质；(d) 倾斜各向异性介质；(e) 任意各向异性介质. Fig. 2 Anisotropic model and conductivity tensor (a) VTI medium; (b) HTI medium; (c) Triple-axis anisotropic medium; (d) Dipping anisotropic medium; (e) Arbitrarily anisotropic medium. |
为求解三维直流电阻率正演问题，我们首先建立如下坐标系：x、y轴位于水平地表，z轴垂直向下.假设电流强度为I的电流源位于r0=(x0, y0, z0)，则各向异性介质中任意位置r=(x, y, z)处的电位U满足如下定解问题(Dey and Marrison, 1979)：
|
(6) |
其中，σ为大地电导率张量，δ为狄拉克函数，n是与边界相关的法向量，而：
|
(7) |
对于电势U，我们引入n阶正交基函数进行谱展开，即：
|
(8) |
则根据伽辽金有限元理论，并选取相同的n阶谱元基函数与权函数φi(i=1, …, Nt)，Nt=(n+1)(n+2)(n+3)/6为每个四面体单元的节点数，我们可将公式(6)转化为积分弱形式后得到如下线性方程组：
|
(9) |
其中，左端项包括全局电势x及系数矩阵K，右端项b为源项.在每个单元上我们有：
|
(10) |
|
(11) |
其中，Ωe表示四面体单元，Γe表示单元边界.公式(10)中的第二项仅在包含外边界的四面体单元中出现.
在每个单元内，φ=[φ1, …, φNt]是一组正交完备的多项式基，可表示为
|
(12) |
其中，r1, r2, …rNt表示插值节点，V表示对称正定的范德蒙德矩阵，而内部元素ψj(ri)为[-1, 1]正交区间内的Proriol-Koornwinder-Dubiner(PKD)多项式，即：
|
(13) |
其中
在正演模拟中，合理的插值节点分布是避免数值解震荡的关键问题所在.前人研究发现Gauss-Lobatto-Legendre(GLL)点集是最优的六面体网格插值节点(Pasquetti and Rapetti, 2010, 图 3a)，然而对于四面体网格，简单地将GLL节点映射到四面体单元中会在顶点位置产生大量重合的节点(图 3b)，因此我们需要在四面体中构建合理的插值点集.本文选用具有良好插值属性的Warp & Blend点集，其主要思想是将等距节点拉伸到最优位置，使得节点分布具有向端点聚拢的特性(Warburton, 2006).图 3c展示了4阶谱插值节点分布.
|
图 3 不同网格四阶插值节点分布 (a) 传统谱元法六面体网格节点(每个单元有125个GLL节点)；(b) 将六面体GLL节点直接映射到四面体网格上(顶点位置上产生大量重合的节点)；(c) 四面体网格W&B节点(每个单元仅有25个节点，数量远小于六面体). Fig. 3 Distribution of fourth-order interpolation nodes on different grids (a) Traditional spectral element method with hexahedral mesh (each unit has 125 GLL nodes); (b) Hexahedral GLL nodes are directly mapped to the tetrahedral mesh (a large number of overlapping nodes are generated at vertex positions); (c) W&B nodes in tetrahedral mesh (each element has only 25 nodes, much less than that of hexahedron). |
在得到插值节点和谱插值基函数后，我们可将公式(10)中每个单元的系数矩阵计算出来.需要强调的是，由于基函数梯度的三重积分项没有解析形式，故采用高斯数值积分计算，即：
|
(14) |
式中，我们已将任意电导率张量展开.同时，为了对基函数进行积分，我们将物理域Ωe转换到参考坐标域Ωs.Nf为参考域内高斯数值积分节点数，wi为权重，J3D为雅克比矩阵，表示了物理坐标系与参考坐标系之间的转换关系，可表示为
|
(15) |
综上所述，本文针对任意各向异性介质的电导率张量模型，从直流电阻率法满足的泊松方程出发，采用灵活的四面体网格对复杂模型进行精细剖分，结合高精度的谱插值基函数与插值节点建立了伽辽金法有限元控制方程，使用多波前直接求解器MUMPS进行线性方程求解(Amestoy et al., 2000)，并利用公式(8)中的插值方法计算任意接收点处的电势.
3 算例分析 3.1 精度验证我们首先使用一个3000 m×3000 m×3000 m的层状各向异性模型验证本文算法的精度，并以Wait(1990)中给出的解析解作为参考，比较使用不同阶次谱元插值基函数时精度变化特征.我们采用二极装置，假设发射源位于坐标原点，沿y轴布设30个间距为2 m的测点，初始模型共有105012个网格.为计算方便，在后续算例中我们采用电阻率张量代替电导率张量.设第一层介质层厚为15 m，主轴电阻率为ρx1=200 Ωm, ρy1=ρz1=50 Ωm，旋转角χ11=χ21=χ31=0°; 第二层主轴电阻率为ρx2=20 Ωm, ρy2=ρz2=5 Ωm，旋转角χ12=χ22=χ32=0°，两层均为沿yoz层理面的HTI介质.
图 4a、b分别给出了不同阶数插值基函数计算得到的视电阻率和相对误差曲线.其中，图 4a中靠近发射源附近的视电阻率趋近于第一层平均电阻率ρM1=100 Ωm，距离较远位置则趋近第二层的平均电阻率ρM2=10 Ωm.可以看出，传统有限元法(n=1)视电阻率数值明显小于解析解.然而，当我们使用高阶算法(n=2，n=3，n=4)，视电阻率曲线均能与解析解曲线很好地吻合.为了更加直观地对比不同阶数插值基函数模拟精度的差异，我们在图 4b展示了相对误差曲线.阶数分别为n=1，n=2，n=3，n=4的谱元法模拟结果与解析解的最大误差分别为24.99%，5.08%，0.75%，0.10%，平均误差为13.39%，1.65%，0.24%，0.05%. 由此可见，随插值函数阶数升高相对误差明显地下降，说明谱元法能够大幅提高直流各向异性正演模拟精度.考虑到高阶谱元均能获得很高的计算精度，后面算例中均采用三阶插值基函数.
|
图 4 测线沿y方向的视电阻率及与解析解的相对误差(图中的n代表谱插值基函数阶数) Fig. 4 Apparent resistivities and relative errors of analytical solution along survey line in y direction (n represents the order of the spectral interpolation basis function) |
本文利用图 5中的模型来研究直流电法中各向异性的特征识别方法.该模型为一个埋藏于均匀半空间的边长为100 m×100 m×100 m的立方体，顶部埋深为5 m.模型尺寸为6000 m×6000 m×3000 m. 围岩主轴电阻率为ρy1=400 Ωm, ρx1=ρz1=100 Ωm，异常体主轴电阻率为ρy2=40 Ωm, ρx2=ρz2=10 Ωm.
|
图 5 立方体异常各向异性模型及测量装置示意图 Fig. 5 Sketch of anisotropic cube model and survey configuration |
我们采用图 5中的二级装置研究两组不同收发距(AM1：r=40 m和AM2：r=150 m)条件下的视电阻变化特征.为获得视电阻率极性图，我们进行36组测量：发射源点位置A固定，接收点绕A点依次顺时针旋转10°，完成一个闭合环路观测.结合二极装置极距与探测深度的关系，我们知道当极距r=40 m时，能够探测到异常体的电阻率分布，而当r=150 m时，能够探测到围岩的电阻率.图 6a给出当异常体电阻率绕z轴旋转0°、30°、45°、60°、90°、135°时的视电阻率极性图(图中从原点到极性图端点的长度表征了视电阻率大小，而其方向表征了二极装置的定向)，图 6b给出围岩电阻率绕x轴旋转0°、30°、45°、60°、90°、135°的视电阻率极性图.由图可以看出，当r=40m时，随着沿z轴方向旋转角度的变化，视电阻率曲线形态基本不变，均为轴对称图形，并且有两个正交的对称轴.然而，曲线绕中心发生了旋转，其中长轴方向角与异常体的旋转角吻合，即视电阻率极性图中长轴方向与走向方向一致.当接收半径r=150 m时，随着沿x轴方向旋转角度的变化视电阻率曲线呈现明显的不对称.长轴的方向基本不变，但是短轴的一侧沿着y轴正方向逐渐被拉伸.当旋转角达到90°时曲线为正圆形.当旋转角为45°和135°时视电阻率极性图关于x轴对称，但两者短轴的拉伸方向发生了翻转.另外，从图可以看出短轴的大头指示了地层倾向，这与围岩电阻率的变化规律相反，呈现电各向异性反常现象.由此可以得出结论，利用视电阻率极性图我们可以确定各向异性主轴方向与倾向，为未来各向异性电阻率数据三维反演提供重要的旋转角度参数信息，减少多解性.
|
图 6 各向异性围岩和异常体主轴电阻率发生旋转时视电阻率极性图 (a) 围岩旋转角χ1=0°, 0°, 0°，异常体旋转角χ2=0°, 0°, 0°/30°/45°/60°/90°/135°，接收半径r=40 m; (b) 围岩旋转角χ1=0°/30°/45°/60°/90°/135°, 0°, 0°, 异常体旋转角χ2=0°, 0°, 0°，接收半径r=150 m. Fig. 6 Polar plots of apparent resistivity when principal resistivity axes of anomaly body and anisotropic host rock rotate (a) Rotation angle of host rock χ1=0°, 0°, 0°, rotation angle of abnormal body χ2=0°, 0°, 0°/30°/45°/60°/90°/135°, receiving radius is r=40 m; (b) Rotation angle of host rock χ1=0°/30°/45°/60°/90°/135°, 0°, 0°, rotation angle of abnormal body χ2=0°, 0°, 0°, receiving radius is r=150 m. |
下面进一步研究多个各向异性参数变化时视电阻率极性图分布特征.图 7中我们给出两组极距的模拟结果.对比图 7a、d中小极距极性图发现，异常体绕x轴旋转45°和135°时，两者曲线关于x轴对称，大头分别指向y轴负方向和正方向.大极距视电阻率极性图长轴对应的方位角与背景电阻率绕z轴旋转角基本相同(分别为135°和45°)，两支曲线关于y轴对称.由于受异常体各向异性的影响，大极距曲线在立方体倾向方向也发生了一定的拉伸.图 7b中围岩与异常体分别绕x轴旋转30°和60°，大极距与小极距对应视电阻率极性图的大头方向相反，分别指向y轴正、负方向，而图 7e中除绕x方向旋转以外，两者均绕z轴再旋转45°，因此两支曲线长轴也发生了旋转，旋转角接近45°，并且曲线大头的方向也发生了改变，大致指示走向与倾向的变化.同样地，我们将模型主轴电导率先绕z轴旋转30°和60°(图 7c)，再绕x轴旋转45°(图 7f)后计算视电阻率极性图.从图 7c可以看出两支曲线分别关于30°和60°呈现轴对称，而从图 7f可以看出两只曲线均发生了拉伸，但拉伸方向相反.由此，我们可以得出结论：当地下介质各向异性属性较为复杂时，通过分析视电阻率极性图的分布形态，我们仍可确定各向异性主轴的走向和倾向.这些参数的确定对未来电阻率数据三维反演中设定合理反演参数、减少各向异性数据反演多解性至关重要.
|
图 7 各向异性围岩和异常体主轴电阻率发生旋转时视电阻率极性图 (a) 围岩旋转角χ1=0°, 0°, 135°, 异常体旋转角χ2=45°, 0°, 0°；(b) 围岩旋转角χ1=30°, 0°, 0°, 异常体旋转角χ2=60°, 0°, 0°；(c)围岩旋转角χ1=0°, 0°, 30°, 异常体旋转角χ2=0°, 0°, 60°；(d) 围岩旋转角χ1=0°, 0°, 45°, 异常体旋转角χ2=135°, 0°, 0°；(e)围岩旋转角χ1= 30°, 0°, 45°, 异常体旋转角χ2=60°, 0°, 45°；(f)围岩旋转角χ1=45°, 0°, 30°, 异常体旋转角χ2=45°, 0°, 60°. Fig. 7 Polar plots of apparent resistivities when the principal resistivity axes of anomaly body and anisotropic host rock rotate (a) Rotation angle of host rock χ1=0°, 0°, 135°, rotation angle of abnormal body χ2=45°, 0°, 0°; (b) Rotation angle of host rock χ1=30°, 0°, 0°, rotation angle of abnormal body χ2=60°, 0°, 0°; (c) Rotation angle of host rock χ1=0°, 0°, 30°, rotation angle of abnormal body χ2=0°, 0°, 60°; (d) Rotation angle of host rock χ1=0°, 0°, 45°, rotation angle of abnormal body χ2=135°, 0°, 0°; (e) Rotation angle of host rock χ1=30°, 0°, 45°, rotation angle of abnormal body χ2=60°, 0°, 45°; (f) Rotation angle of host rock χ1=45°, 0°, 30°, rotation angle of abnormal body χ2=45°, 0°, 60°. |
为检验本文算法的计算效率，我们还对比了针对不同网格的谱元和有限元算法的计算结果.为此，我们将图 5中的立方体主轴电阻率绕x轴旋转90°，并在图 8中给出了三种不同情况的视电阻率极性图.图 8a为多次自适应网格加密的有限元计算结果(Ren and Tang, 2010，开源代码https://software.seg.org/2010/0002/index.html).其中，Ap_0为初始网格极性图，Ap_1、Ap_2、Ap_3为网格自适应加密1、2及3次的极性图.可以看出随着网格加密极性图越来越光滑，视电阻率计算精度越来越高.图 8b、c分别给出粗、细网格上不同阶数谱元法的极性图.有图可以看出，当谱元基函数阶数为n=2时，粗、细网格上的视电阻率极性图均很光滑.进一步对比图 8b、c可以发现，在细网格上n=3与n=4的视电阻率计算结果已很好地收敛.
|
图 8 不同网格与不同算法视电阻率极性图 (a) 自适应加密网格的有限元算法加密不同次数；(b) 粗网格谱元法不同阶数；(c) 细网格谱元法不同阶数. Fig. 8 Polar plots of apparent resistivities with different grids and different algorithms (a) Adaptive finite element algorithm with different grid refinement; (b) Coarse grid spectral element method with different orders; (c) Fine grid spectral element method with different orders. |
如前所述，加密网格与提高阶数均能提升谱元法数值模拟精度.本节我们对不同计算方法之间的计算效率进行对比.以细网格4阶谱元模拟结果作为参考, 我们计算针对不同方法及不同网格模拟结果的相对误差.表 1给出了对比结果.图 8a、b、c分别对应表中网格加密的有限元法、粗网格谱元法和细网格谱元法.由表可以看出，通过简单加密网格有限元法相对误差逐渐减小，但多次网格加密之后的相对误差仍然较大.相比之下，谱元法能够在保持初始网格不变的条件下，通过提高阶数使得模拟精度得到明显的提升.自适应加密一次网格(AP_1)与2阶插值(n=2)未知数数量相近，但谱元法的相对误差远小于有限元法，在计算时间上也略少于后者.对比网格自适应加密三次(AP_3)有限元与4阶谱元(n=4)的计算结果可知，谱元法只需求解较少的未知数就能达到更高的精度.进一步对比粗、细网格下不同阶次谱元法的网格数、自由度、求解时间、计算误差等参数发现，采用更加精细的网格比粗网格能够获得更高的计算精度.由此可以得出结论：基于谱元法我们可以通过采用双重策略获得高精度的数值计算结果.
|
表 1 自适应加密有限元法与谱元法参数对比 Table 1 Forward modelling data comparison between adaptive finite element method and spectral element method |
为检验四面体网格谱元法对于复杂各向异性模型的适应性，我们设计了一个如图 9所示的高度差为20 m，直径为100 m的三维山脊地形模型.山脊顶部中心点的坐标为(0, 0, -20).在地形下方10 m处埋藏一个半径为15 m的球体，模型区域为3000 m× 3000 m×3000 m，共有165793个四面体网格.测线沿x轴布设，间距为2m.由于在源点附近、电阻率分界面及地表面等区域会产生急剧变化的场，为保证算法稳定，我们对这些区域进行了局部网格加密(图 9b).
|
图 9 三维山脊地形模型及其剖分网格 Fig. 9 3D ridge terrain model with a sphere embedded and grid |
下面首先研究各向异性模型的电流分布情况.将源点置于地表面(-25, 0, -10)处.围岩设为倾斜各向异性介质，主轴电阻率为ρx1=200 Ωm, ρy1=ρz1= 50 Ωm，各向同性球体电阻率为ρx2=ρy2=ρz2=10 Ωm.根据公式：
|
(16) |
|
(17) |
即可计算出任意点的电流密度矢量J.由此我们绘制了如图 10所示的xz面上电流密度等值线切片图.其中，图 10a—d分别对应于围岩主轴电导率绕y轴旋转0°、45°、90°、135°的结果.图中球形异常体的轮廓用白线圈出.观察不同旋转轴对应的电流场分布，我们可以看出电流向层理面倾斜的方向集聚，这是由电流向低阻方向聚焦产生的通道效应.随着旋转角的增大，电流等值线由垂向逐渐发生倾斜，直至旋转角为90°时电流等值线变为近乎水平.旋转角为135°时电流倾斜方向与45°正好相反.由于围岩各向异性的影响，低阻异常体内部的电流密度也发生了畸变.
|
图 10 围岩为倾斜各向异性介质时xz平面上电流密度分布 主轴电阻率绕x轴的旋转角(a)χ1=0°, 0°, 0°；(b) χ1=0°, 45°, 0°；(c) χ1=0°, 90°, 0°；(d) χ1=0°, 135°, 0°. Fig. 10 Current density distribution on the xz- plane when host rock is dipping anisotropic The rotation angle of theprincipal resistivity axis around the x axis are(a)χ1=0°, 0°, 0°; (b) χ1=0°, 45°, 0°; (c) χ1=0°, 90°, 0°; (d) χ1=0°, 135°, 0°. |
为进一步研究地形效应及其校正问题，我们针对联合剖面装置计算了四组不同的模型：(1)围岩为HTI介质，主轴电阻率分别为ρx1=200 Ωm, ρy1=ρz1= 50 Ωm，球体电阻率ρx2=ρy2=ρz2=10 Ωm；(2)异常体与围岩均为HTI介质，围岩主轴电阻率ρx1=200 Ωm, ρy1=ρz1=50 Ωm，而球体主轴电阻率ρx2=20 Ωm, ρy2=ρz2=5 Ωm；(3) 异常体为HTI介质，主轴电阻率ρx2=20 Ωm, ρy2=ρz2=5 Ωm，围岩电阻率ρx1=ρy1=ρz1=100 Ωm；(4)异常体与围岩均为各向同性介质：围岩电阻率ρx1=ρy1=ρz1=100 Ωm，球体电阻率ρx2=ρy2=ρz2=10 Ωm.
针对联合剖面装置形式，我们计算了上述四个模型的视电阻率剖面曲线，并与纯地形(不存在异常体)的视电阻率曲线进行对比，如图 11a—d所示.每组联合剖面测量产生A、B两支曲线，每对曲线均有三个公共交点，其中山脊顶部对应高阻交点，在两侧山脚处出现低阻交点.由于地形效应远大于异常体响应，为消除地形效应我们利用如下地形校正公式(Fox et al., 1980)：
|
(18) |
|
图 11 各向异性山脊模型联合剖面装置视电阻率曲线(模型参见图 8) (a) 围岩为HTI介质；(b) 围岩与球体异常均为HTI介质；(c) 球体异常体为HTI介质；(d) 围岩与球体异常均为各向同性介质；(e)—(h) 地形校正后视电阻率曲线. Fig. 11 Apparent resistivity of the composite profile device for the anisotropic ridge model shown in Fig. 8 (a) Host rock is HTI; (b) Host rock and sphere are both HTI; (c) Sphere is HTI; (d) Host rock and sphere are isotropic; (e)—(h) Apparent resistivity curves after terrain correction. |
其中，ρs为实测的电阻率，ρsL为纯地形电阻率，而ρsG为校正后的视电阻率.对于各向异性介质背景电阻率取为平均电阻率
利用式(18)，我们对上述四个模型的视电阻率进行地形校正，结果如图 11e—h所示.从图可以看出，四组校正后的数据均很好地展示球形良导体上的联合剖面曲线.曲线左右支交点均位于球体中心点位置(x=0).各向异性围岩在地形校正后的视电阻率值远大于各向异性异常体的响应.对比图 11e、f和图 11g、h，我们还发现各向异性的球体使得视电阻率的极小值更小.必须指出，由于地形和各向异性之间的复杂耦合关系，导致地下电流场分布异常复杂(参见图 10)，因此简单应用地形校正无法完全消除地形效应的影响.然而，从图 11给出的算例可以看出，即使大地存在起伏地形和复杂各向异性，视电阻率响应经地形校正后仍能很好地确定异常体位置与导电性特征.
4 结论本文将四面体网格的谱元法成功应用于各向异性介质直流电阻率法正演模拟之中，并通过层状模型的解析解验证了算法的正确性.相比于传统有限元算法，谱元法在提高数值模拟精度上具有独特优势.它可以灵活利用加密网格或提高谱插值基函数阶数，改善计算精度.通过对各向异性模型模拟及响应特征分析发现，利用视电阻率极性图的长短轴走向以及极性图不对称性，可以很好地确定地下介质各向异性主轴电阻率及其走向和倾向.针对复杂地形条件下各向异性大地模型，采用传统的地形校正方法仍然可以有效地消除地形效应，获得地下导电体的可靠信息.这将有助于对各向异性环境下地下三维电性结构进行精确反演和解释.
致谢 感谢中南大学邱乐稳博士、陈煌博士以及吉林大学殷长春教授研究团队其他成员对本文研究提供的帮助，感谢两位审稿人提出的宝贵意见.
Aissaoui R, Bounif A, Zeyen H, et al. 2019. Evaluation of resistivity anisotropy parameters in the eastern Mitidja basin, Algeria, using azimuthal electrical resistivity tomography. Near Surface Geophysics, 17(4): 359-378. |
Amestoy P R, Duff I S, L'Excellent J Y. 2000. Multifrontal parallel distributed symmetric and unsymmetric solvers. Computer Methods in Applied Mechanics and Engineering, 184(2-4): 501-520. DOI:10.1016/S0045-7825(99)00242-X |
Cai Y Y, Xiao T J, Song T. 2020. Influence of two-dimensional magnetotelluric anistropic parameters on apparent resistivities. Progress in Geophysics (in Chinese), 35(1): 86-93. DOI:10.6038/pg2020CC0478 |
Chen H B, Li T L, Shi H Y, et al. 2018. An edge-based finite element method for 3D marine controlled-source electromagnetic forward modeling with a new type of second-order tetrahedral edge element. International Journal of Applied Electromagnetics and Mechanics, 57(2): 217-233. DOI:10.3233/JAE-170153 |
Cohen G C. 2002. Higher-order numerical methods for transient wave equations. //Scientific Computation. New York: Springer-Verlag.
|
Davydycheva S, Druskin V, Habashy T. 2003. An efficient finite-difference scheme for electromagnetic logging in 3D anisotropic inhomogeneous media. Geophysics, 68(5): 1525-1536. DOI:10.1190/1.1620626 |
Dey A, Morrison H F. 1979. Resistivity modeling for arbitrarily shaped three-dimensional structures. Geophysics, 44(4): 753-780. DOI:10.1190/1.1440975 |
Dubiner M. 1991. Spectral methods on triangles and other domains. Journal of Scientific Computing, 6(4): 345-390. DOI:10.1007/BF01060030 |
Fox R C, Hohmann G W, Killpack T J, et al. 1980. Topographic effects in resistivity and induced-polarization surveys. Geophysics, 45(1): 75-93. DOI:10.1190/1.1441041 |
Grayver A V, Bürg M. 2014. Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method. Geophysical Journal International, 198(1): 110-125. DOI:10.1093/gji/ggu119 |
Gu Y, Yin C C, Qiu C K, et al. 2020. Three-dimensional anisotropic forward modeling for marine CSEM method in time-domain based on rational Krylov subspace method. Progress in Geophysics (in Chinese), 35(6): 2332-2350. DOI:10.6038/pg2020EE0295 |
Günther T, Rücker C, Spitzer K. 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography-Ⅱ. Inversion. Geophysical Journal of the Royal Astronomical Society, 166(2): 506-517. DOI:10.1111/j.1365-246X.2006.03011.x |
Herwanger J V, Pain C C, Binley A, et al. 2004. Anisotropic resistivity tomography. Geophysical Journal International, 158(2): 409-425. DOI:10.1111/j.1365-246X.2004.02314.x |
Huang X, Yin C C, Farquharson C G, et al. 2019. Spectral-element method with arbitrary hexahedron meshes for time-domain 3D airborne electromagnetic forward modeling. Geophysics, 84(1): E37-E46. DOI:10.1190/geo2018-0231.1 |
Jaysaval P, Shantsev D, De La Kethulle de Ryhove Sébastien. 2014. Fast multimodel finite-difference controlled-source electromagnetic simulations based on a Schur complement approach. Geophysics, 79(6): E315-E327. DOI:10.1190/geo2014-0043.1 |
Karniadakis G, Sherwin S. 2005. Spectral/hp element methods for computational fluid dynamics. Oxford: Oxford University Press.
|
Key K, Ovall J. 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling. Geophysical Journal International, 186(1): 137-154. DOI:10.1111/j.1365-246X.2011.05025.x |
Kneisel C, Emmert A, Kästl J. 2014. Application of 3D electrical resistivity imaging for mapping frozen ground conditions exemplified by three case studies. Geomorphology, 210: 71-82. DOI:10.1016/j.geomorph.2013.12.022 |
Komatitsch D, Ritsema J, Tromp J. 2002. The spectral-element method, Beowulf computing, and global seismology. Science, 298(5599): 1737-1742. DOI:10.1126/science.1076024 |
Koornwinder T. 1975. Two-variable analogues of the classical orthogonal polynomials. //Askey R ed. Theory and Application of Special Functions. New York: Academic Press, 435-495.
|
Li Y G, Spitzer K. 2005. Finite element resistivity modelling for three-dimensional structures with arbitrary anisotropy. Physics of the Earth and Planetary Interiors, 150(1-3): 15-27. DOI:10.1016/j.pepi.2004.08.014 |
Liu Y H, Yin C C, Cai J, et al. 2018. Review on research of electrical anisotropy in electromagnetic prospecting. Chinese Journal of Geophysics (in Chinese), 61(8): 3468-3487. DOI:10.6038/cjg2018L0004 |
Liu Y S, Teng J W, Lan H Q, et al. 2014. A comparative study of finite element and spectral element methods in seismic wavefield modeling. Geophysics, 79(2): T91-T104. DOI:10.1190/geo2013-0018.1 |
Loke M H, Barker R D. 1996. Practical techniques for 3D resistivity surveys and data inversion. Geophysical Prospecting, 44(3): 499-523. DOI:10.1111/j.1365-2478.1996.tb00162.x |
Loke M H, Chambers J E, Rucker D F, et al. 2013. Recent developments in the direct-current geoelectrical imaging method. Journal of Applied Geophysics, 95: 135-156. DOI:10.1016/j.jappgeo.2013.02.017 |
Maday Y, Meiron D, Patera A T, et al. 1993. Analysis of iterative methods for the steady and unsteady stokes problem: application to spectral element discretizations. SIAM Journal on Scientific Computing, 14(2): 310-337. DOI:10.1137/0914020 |
Neyamadpour A. 2019. 3D electrical resistivity tomography as an aid in investigating gravimetric water content and shear strength parameters. Environmental Earth Sciences, 78(19): 583. DOI:10.1007/s12665-019-8603-7 |
Nthaba B, Shemang E, Hengari A, et al. 2020. Characterizing coal seams hosted in Mmamabula Coalfield, central Botswana using pseudo-3D electrical resistivity imaging technique. Journal of African Earth Sciences, 167: 103866. DOI:10.1016/j.jafrearsci.2020.103866 |
Pasquetti R, Rapetti F. 2010. Spectral element methods on unstructured meshes: which interpolation points?. Numerical Algorithms, 55(2-3): 349-366. DOI:10.1007/s11075-010-9390-0 |
Proriol J. 1957. Sur une famille de polynomes a deux variables orthogonaux dans un triangle. Comptes Rendus Hebdomadaires des Seances de l'academie des Sciences, 245: 2459-2461. |
Ren Z Y, Chen H, Tang J T, et al. 2018a. A volume-surface integral approach for direct current resistivity problems with topography. Geophysics, 83(5): E293-E302. DOI:10.1190/geo2017-0577.1 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling. Geophysical Journal International, 194(2): 700-718. DOI:10.1093/gji/ggt154 |
Ren Z Y, Qiu L W, Tang J T, et al. 2018b. 3-D direct current resistivity anisotropic modelling by goal-oriented adaptive finite element methods. Geophysical Journal International, 212(1): 76-87. DOI:10.1093/gji/ggx256 |
Ren Z Y, Tang J. 2010. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method. Geophysics, 75(1): H7-H17. DOI:10.1190/1.3298690 |
Runge C. 1901. Vber empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten. Zeitschrift für Mathematik und Physik, 46: 224-243. |
Rücker C, Günther T, Spitzer K. 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography-Ⅰ. modelling. Geophysical Journal International, 166(2): 495-505. DOI:10.1111/j.1365-246X.2006.03010.x |
Sheng F, Gao G Z, Verdin C T. 2006. Efficient 3D electromagnetic modelling in the presence of anisotropic conductive media, using integral equations. Exploration Geophysics, 37(3): 239-244. DOI:10.1071/EG06239 |
Song T, Liu Y, Wang Y. 2017. Finite element method for modeling 3d resistivity sounding on anisotropic geoelectric media. MathematicalProblems in Engineering, 2017: 8027616. DOI:10.1155/2017/8027616 |
Udphuay S, Günther T, Everett M E, et al. 2011. Three-dimensional resistivity tomography in extreme coastal terrain amidst dense cultural signals: application to cliff stability assessment at the historic D-Day site. Geophysical Journal International, 185(1): 201-220. DOI:10.1111/j.1365-246X.2010.04915.x |
Ueda T, Zhdanov M S. 2006. Fast numerical modeling of multitransmitter electromagnetic data using multigrid quasi-linear approximation. IEEE Transactions on Geoscience and Remote Sensing, 44(6): 1428-1434. DOI:10.1109/TGRS.2006.864386 |
Wait J R. 1990. Current flow into a three-dimensionally anisotropic conductor. Radio Science, 25(5): 689-694. DOI:10.1029/RS025i005p00689 |
Wang T, Signorelli J. 2004. Finite-difference modeling of electromagnetic tool response for logging while drilling. Geophysics, 69(1): 152-160. DOI:10.1190/1.1649383 |
Wang W, Wu X P, Spitzer K. 2013. Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids. Geophysical Journal International, 193(2): 734-746. DOI:10.1093/gji/ggs124 |
Warburton T. 2006. An explicit construction of interpolation nodes on the simplex. Journal of Engineering Mathematics, 56(3): 247-262. |
Yin C C. 2000. Geoelectrical inversion for a one-dimensional anisotropic model and inherent non-uniqueness. Geophysical Journal International, 140(1): 11-23. DOI:10.1046/j.1365-246x.2000.00974.x |
Yin C C, Qi Y F, Liu Y H. 2016a. 3D time-domain airborne EM modeling for an arbitrarily anisotropic earth. Journal of Applied Geophysics, 131: 163-178. DOI:10.1016/j.jappgeo.2016.05.013 |
Yin C C, Weidelt P. 1999. Geoelectrical fields in a layered earth with arbitrary anisotropy. Geophysics, 64(2): 426-434. DOI:10.1190/1.1444547 |
Yin C C, Zhang B, Liu Y H, et al. 2016b. A goal-oriented adaptive finite-element method for 3d scattered airborne electromagnetic method modeling. Geophysics, 81(5): E337-E346. DOI:10.1190/geo2015-0580.1 |
Yoon D, Zhdanov M S, Mattsson J, et al. 2016. A hybrid finite-difference and integral-equation method for modeling and inversion of marine controlled-source electromagnetic data. Geophysics, 81(5): E323-E336. DOI:10.1190/geo2015-0513.1 |
Zhou B, Greenhalgh M, Greenhalgh S A. 2009. 2.5-D/3-D resistivity modelling in anisotropic media using gaussian quadrature grids. Geophysical Journal of the Royal Astronomical Society, 176(1): 63-80. DOI:10.1111/j.1365-246X.2008.03950.x |
Zhu J, Yin C C, Liu Y S, et al. 2020. 3-D dc resistivity modelling based on spectral element method with unstructured tetrahedral grids. Geophysical Journal International, 220(3): 1748-1761. DOI:10.1093/gji/ggz534 |
蔡义宇, 肖调杰, 宋滔. 2020. 二维大地电磁各向异性参数对视电阻率的影响研究. 地球物理学进展, 35(1): 86-93. DOI:10.6038/pg2020CC0478 |
谷宇, 殷长春, 邱长凯, 等. 2020. 基于有理Krylov子空间法的可控源海洋电磁三维各向异性正演研究. 地球物理学进展, 35(6): 2332-2350. DOI:10.6038/pg2020EE0295 |
刘云鹤, 殷长春, 蔡晶, 等. 2018. 电磁勘探中各向异性研究现状和展望. 地球物理学报, 61(8): 3468-3487. DOI:10.6038/cjg2018L0004 |