1 引言

以CSAMT为代表的可控源电磁法是通过观测和研究人工供电产生电磁场的频率域测深方法(汤井田和何继善, 2005).在不考虑地下介质电性各向异性的情形下, 野外数据采集、数据处理、反演成像解释方面的研究日趋成熟(汤井田等, 2015).该方法作为一种重要的地球物理勘探方法, 具有抗干扰能力强，信噪比高等特点.在寻找深部隐伏矿床, 油气资源探测, 地热以及水文、工程地质勘查等方面都取得了比较好的应用效果(何继善, 1990; 罗延钟和万乐, 1996; 石昆法, 1999; 底青云等, 2002; 柳建新等, 2009; 何继善, 2010; 王显祥等, 2014).

2.5维电磁问题是二维问题和三维问题的一种折中, 它兼顾了物理的真实性和计算成本.在各向同性介质数值模拟方面, Stoyer和Greenfield(1976)、Doherty(1988)分别将有限差分法和积分方程法应用于磁偶极源的2.5维可控源电磁法数值模拟中; Coggon(1971)首次将有限元方法应用于2.5维电磁问题中, 取得了较好的效果.随后大量学者(Lee and Morrison, 1985; Unsworth, 1993; Mitsuhata, 2000; 底青云等, 2004; Li and Key, 2007; 戴世坤等, 2013)利用有限单元法实现了可控源电磁法2.5维数值模拟.其中, Mitsuhata(2000)引入伪δ函数避免场源的奇异性, 采用等参单元模拟复杂地形, 提高了数值模拟精度.底青云等(2004)实现了基于二次场的可控源电磁法2.5D正演, 并通过与野外实测剖面数据进行比较表明其算法的可靠性.Li和Key(2007)、Li和Constable(2007)采用非结构化网格(Key and Weiss, 2006)剖分, 将自适应有限单元法引入海洋可控源电磁法2D数值模拟中, 提高了复杂模型的数值模拟精度, 并研究了水深对电磁响应的影响规律.戴世坤等(2013)基于双二次插值的有限单元法对可控源电磁法2.5维问题进行了模拟, 并在此基础上实现了对电磁场值直接反演的算法.在实际的地球物理勘探中, 地下介质异常复杂, 尤其对于沉积地层, 宏观上表现出来的电学性质往往是各向异性的, 而各向异性对可控源电磁响应的影响很大.为了建立更加符合地下实际电性分布的地质模型, 实现精细化反演成像, 近年来, 不少学者针对电导率各向异性介质进行了2.5D数值模拟研究.Kong等(2008)在各向同性数值模拟的基础上, 对水平各向异性介质海洋可控源电磁法2.5D数值模拟进行了研究.Li和Dai(2011)采用自适应有限元法实现了基于二次场的2.5D海洋可控源电磁法正演, 并研究了各向异性参数对海洋可控源电磁响应的影响规律.

目前, 对于2.5维电磁问题, 为了降低场源奇异性影响和截断边界对计算区域的影响, 可控源电磁法2.5D数值模拟一般将总场分解为一次场和二次场分别进行求解(Kong et al., 2008; Li and Dai, 2011), 在水平地形条件下可以达到很高的数值模拟精度.但是, 由于带地形的一次场没有解析解, 传统的地形校正方法精度较低, 该方法很难满足任意复杂地形条件下的陆地可控源电磁法各向异性介质数值模拟精度要求.因此, 研究基于总场的高效、高精度可控源电磁法2.5维各向异性介质正演具有重要意义, 而边界条件和场源奇异性是制约总场算法的两个重要因素.本文首先从总场算法(Mitsuhata, 2000)出发, 详细推导了考虑电导率和介电常数各向异性、适应复杂地形的全张量2.5维微分控制方程, 避免二次场算法中带地形一次场不易求取的问题.采用结构化网格加密收缩方法(Zhang et al., 2016)有效地降低了场源奇异性的影响，引入行波分解边界条件降低了截断边界对计算区域的影响, 并利用有限单元法进行傅里叶变换, 实现了适应复杂地形、任意各向异性介质的可控源电磁法2.5维有限元数值模拟, 为下一步高效、高精度、精细化各向异性介质反演成像奠定了基础.其次, 通过与二次场算法对比主剖面电磁场振幅响应, 验证本文算法的准确性和可靠性, 并研究了任意各向异性倾角对可控源电磁响应的影响规律.最后, 对山丘条件下的典型高、低阻模型进行了模拟, 分析了各向异性介质电磁场响应规律.

2 各向异性介质2.5D问题

在准静态条件下, 假设时谐因子取e^iωt, 其中为角频率, t为时间.则带源的频率域麦克斯韦方程为(Mitsuhata, 2000)：

(1)

(2)

式中, E为电场强度, H为磁场强度, M_s为外加源磁化强度, J_s为外加源电流密度, 为阻抗率, μ为磁导率, 为导纳率张量.

在导纳率张量中, σ为任意各向异性的电导率张量, ε为任意各向异性的介电常数张量.具体形式为

(3)

其中x、y、z表示直角坐标系三个坐标方向.这两个张量均可以通过旋转建立与主坐标轴的关系.对于电导率张量, 用三个主坐标轴电导率σ_x、σ_y、σ_z和各向异性走向角(α_s)、各向异性倾角(α_d)、各向异性偏角(α_l)表示.即：

(4)

其中：

设x轴水平向右为正, z轴垂直向上为正, y轴为构造走向.式(1)—(2) 为三维问题, 沿构造走向y方向做一维傅里叶变换，公式为

(5)

其中(x, y, z)为空间坐标, k_y为波数, 为波数域电磁场值, F为空间域电磁场值.可得二维问题公式为：

(6)

(7)

(8)

(9)

(10)

(11)

联立式(6)—(11), 可以得到关于两个平行于构造走向方向的电磁场分量和的耦合微分控制方程为：

(12)

(13)

其中是2×2的矩阵, 这四个矩阵都是与波数k_y以及电导率张量单元σ_ij(i, j=x, y, z)相关的, 具体形式为：

其中：

其他系数是与波数和电导率张量有关的函数, 具体表达式为:

该微分方程考虑了电导率和介电常数的各向异性, 忽略变化很小的磁导率各向异性, 既可在低频段作为可控源电磁法数值模拟的控制方程，也可在高频段作为探地雷达数值模拟的控制方程, 具有普适性.对于可控源电磁法所涉及的频率, 地下介质σ≫ωε，因而可以忽略介电常数的各向异性.本文只考虑电导率的倾斜各项异性, 即那么, 两个平行于构造走向方向的电磁场分量耦合微分方程式(12)—(13) 可以简化为：

(14)

(15)

电场分量和磁场分量可以简化为：

(16)

(17)

(18)

(19)

这样, 就得到了波数域以和为变量的耦合控制方程式(14) 和(15), 其余的辅助场分量由式(16)—(19) 求得.

3 有限单元法

本文采用三角形单元进行剖分, 三角形单元中采用线性插值方式(徐世浙, 1994).运用Galerkin法(徐世浙, 1994)将式(14)—(15) 满足的微分控制方程推导为有限元方程(附录A)，公式为：

(20)

(21)

式(21) 即为和满足的有限元方程.其中, N_i^e为第e个单元第i个节点对应的线性形函数.通过单元分析和总体合成(Li and Key, 2007; Kong et al., 2008; Li and Dai, 2011)得到一个大型、稀疏的复线性方程组.公式为

(22)

其中k_y为波数, A(k_y)是大型稀疏复矩阵, x(k_y)是波数域的解向量, S(k_y)是源向量.每一个波数都需要计算一次稀疏矩阵, 求解一次方程组, 因此线性方程组的求解速度和精度直接影响整个数值模拟速度和精度.本文在电脑配置为CPU-Inter Core i7，主频2.80 GHz，内存8.00 GB上运行程序, 均采用PARDISO_64直接求解器(张钱江等, 2016)求解方程.其他波数域电磁场分量通过式(16)—(19) 求取, 其中, 采用五点差分法计算.

4 提高计算精度的措施

在基于总场计算的可控源电磁法2.5维各向异性介质正演中, 场源奇异性、边界条件和反傅里叶变换是制约总场算法数值精度的几个重要因素.本文引入网格加密-收缩算法, 通过均匀加密源线网格节点以降低场源奇异性的影响范围，通过收缩空气层和地下深层网格节点以提升数值计算效率; 引入行波分解边界条件降低截断边界的影响; 提出任意采样反傅里叶变换方法快速、高精度地计算出空间域电磁场分量.

4.1 场源奇异性处理

可控源电磁法2.5D数值模拟中的场源奇异性问题, 通常采用的解决办法主要是一次场、二次场分别求解(Li and Dai, 2011)和伪δ函数近似场源(Mitsuhata, 2000).但一次场、二次场分别求解的办法的不足是:在地表起伏或者存在复杂地质体时, 很难找到精确的一次场; 伪δ函数近似场源的办法很难处理多源条件下的场源奇异性问题.本文采用结构化网格加密收缩(Zhang et al., 2016)的办法:一方面可以通过局部加密降低单场源奇异性的影响范围, 在多场源条件下, 可以实现在浅层水平方向均匀加密网格单元, 降低场源奇异性的影响范围; 另一方面在深层通过网格隔层收缩的方式将多个单元合成一个单元, 以达到减小网格节点规模, 提升数值计算效率的目的.图 1为网格加密收缩示意及效果图, 图 1a为常规网格, 图 1b为网格加密, 图 1c为网格收缩, 图 1d为全空间模型(电阻率为100 Ωm)的Ey响应.其中, 图 1d中，发射源为沿地表x方向的水平长导线源, 长度为70 m, 发射电流幅值为50 A, 发射频率为0.1 Hz.从图 1d和图 5可以看出网格加密有效地降低了场源奇异性的影响, 网格隔层收缩数值稳定性较好, 与解析解的吻合度较高.

4.2 吸收边界条件加载

目前, 可控源电磁法2.5维数值模拟中主要采用Dirichlet边界条件，并通过扩边处理降低截断边界带来的影响(Mitsuhata, 2000).但是, 该边界条件在靠近边界区域仍存在很大的误差.为了进一步降低截断边界的影响, 薛东川和戴世坤(2008)借鉴地震波模拟中的吸收边界条件(Reynolds, 1978), 将波数域电磁场控制方程分解为两个传播方向相反的单程波方程, 以沿边界外法向衰减的单程波方程作为该边界上的吸收边界条件, 称为行波分解法吸收边界条件.基于二次插值的吸收边界条件具体构造方法及具体形式见文献(薛东川和戴世坤, 2008).由于本文采用线性插值, 故忽略二阶及以上积分项, 得到基于线性插值的边界积分表达式为：

(23)

(24)

(25)

(26)

为了验证该边界条件的有效性, 与常规的齐次边界条件和解析解进行对比.设计一水平各向同性全空间模型:水平电阻率为10 Ωm, 垂向电阻率为10 Ωm.发射源沿x方向布设在坐标原点, 发射电流幅值为20 A, 发射频率为10 Hz.图 2为300 m旁侧剖面上电场E_y和磁场H_y不同边界条件的数值解和解析解对比图, 图 2a为E_y分量, 图 2b为H_y分量.其中, 红色实线为齐次边界条件; 蓝色实线为本文吸收边界条件; 绿色点线为解析解.从图中可以看出行波分解法吸收边界条件能够有效压制边界反射, 显著提高计算精度, 且效果明显优于常规的齐次边界条件, 与解析解的吻合度较高.

4.3 任意采样反傅里叶变换

计算得到波数域电磁场分量后, 通常采用算法相对比较复杂的三次样条插值算法进行傅里叶反变换求取空间域电磁场分量(Mitsuhata, 2000).为了更加精确、高效地求取反傅里叶变换积分, 引入有限单元法中的长度坐标(徐世浙, 1994), 如图 3所示.积分公式为

(27)

令：

(28)

本文基于二次插值进行有限单元积分, 则设G是波数单元中的二次插值函数(式中λ₁、λ₂、λ₃为任意常数)，公式为

(29)

将式(27) 和式(28) 代入式(29), 可得：

(30)

其中e为第e个单元, Ne是单元个数, k_y^j、k_y^p、k_y^m为任意相邻的三个波数值.

其中G_j、G_p、G_m分别是一个波数单元内相邻三个波数对应的波数域电磁场G(x, y, z, ω, k_y^j)、G(x, y, z, ω, k_y^p)、G(x, y, z, ω, k_y^m)的简记形式.

特别地, 当y=0, 即主剖面的空间域电磁场积分表达式可简化为

(31)

其中l=k_y^m－k_y^j是单元的长度.

5 算例分析 5.1 精度验证

为了验证本文算法的正确性和可靠性, 将该算法的数值模拟结果与水平各向异性层状介质模型的解析解进行对比.模型设计参考已有文献(Li and Dai, 2011), 如图 4所示.空气电阻率为1×10⁸ Ωm; 第一层为海水层, 海水电阻率为0.3 Ωm, 海水厚度为1020 m; 第二层为各向异性介质, 主坐标轴电阻率分别为ρ_x=ρ_y=1.0 Ωm, ρ_z=4 Ωm, 厚度为990 m; 第三层为高阻油气薄层, 电阻率为50 Ωm, 厚度为100 m, 第四层为各向同性介质, 电阻率为1 Ωm.发射源沿着x方向布设在距离海底30 m的海水中, 发射频率为0.25 Hz, 发射电流幅值为1 A, 发射线源长度为70 m.接收装置布设在y=100 m的海底-10~10 km范围内.

图 5给出了水平各向异性层状介质模型电磁场响应的数值解和解析解.其中, 实线为电磁场振幅解析解, *为电磁场振幅数值解.从图中可以看出, 本文算法计算得到的电磁场与解析解的吻合度较高, 相对误差均在1.5%以下，表明本文算法正确、且精度较高.

5.2 高阻油气藏模型

为了进一步验证本文算法的可靠性和有效性, 研究任意各向异性倾角电磁场分量的振幅响应规律.本文设计海洋高阻油气藏模型(Li and Dai, 2011), 如图 6f所示.其中, 空气电阻率为1×10⁸ Ωm, 海水电阻率为0.3 Ωm, 海水厚度为1 km, 海底地层为各向异性介质, x主坐标轴电阻率为1 Ωm, y主坐标轴电阻率为1 Ωm, z主坐标轴电阻率为4 Ωm, 各向异性倾角α_d分别为0°、30°、45°、60°和90°.根据各向异性倾角的不同, 可以将其分为三种类型:当α_d=0°时, 称为垂直轴对称的横向各向异性(Transverse Isotropy Vertical, 简记为TIV), 当α_d=30°、45°、60°时, 称为倾斜轴对称的横向各向异性(Transverse Isotropy Dip, 简记为TID), 当α_d=90°时, 称为水平轴对称的横向各向异性(Transverse Isotropy Horizontal, 简记为TIH).在海底以下1 km处赋存一个各向同性的高阻油气藏, 电阻率为50 Ωm, 宽为6 km, 高为100 m.发射频率为0.25 Hz, 发射源沿x方向布设在距离海底50 m的海水中, 发射电流幅值为1 A, 发射线源长度为70 m.接收装置沿x方向均匀布设在海底-10~10 km范围内.

首先通过与二次场算法(Li and Dai, 2011)对比主剖面E_x振幅响应, 验证本文算法的可靠性, 如图 6a所示.从图中可以看出随着收发距的增大, 振幅曲线逐渐分离; 由于垂向电阻率最大, 在横向各向异性(TIV)情况下, 电场E_x幅值最大, 随着各向异性倾角的逐渐增大振幅衰减变快.与文献Dai和Li(2011)的研究规律一致, 进一步表明本文算法的可靠性和有效性.在此基础上, 研究任意各向异性倾角对其他电磁场分量的影响规律, y=100 m旁侧剖面电磁场分量的振幅响应, 如图 6b—e所示.从图中可以看出各向异性倾角对异常体正上方的E_y分量影响较大; 而对于其他分量而言, 各向异性倾角的影响随着收发距的增大而增大, 尤其对于E_x分量和H_y分量表现最为明显.且存在各向异性时, 振幅响应曲线表现为相对增强, 这是由于电阻率增大, 与各向同性情况下相比, 电磁场值衰减变慢.

5.3 带地形各向异性模型

设计如图 7所示的模型, 图 7a为低阻体; 图 7b为高阻体.空气电阻率为1×10⁸ Ωm, 山丘高50 m, 顶部延伸3 km, 底部延伸3.1 km.地下围岩为各向同性介质, 距水平地面为1 km的位置有一各向异性异常体, 距发射源的水平距离分别为1 km, 大小均为3000 m×100 m.分别通过低阻各向异性异常体模型和高阻各向异性异常体模型针对各向异性对可控源电磁响应的影响.图 7a为低阻体模型:围岩电阻率为100 Ωm, 异常体的水平电阻率ρ_h和垂向电阻率ρ_v分别为(ρ_h=10 Ωm, ρ_v=10 Ωm)、(ρ_h=10 Ωm, ρ_v=1 Ωm)、(ρ_h=20 Ωm, ρ_v=10 Ωm).图 7b为高阻体模型:围岩电阻率为5 Ωm, 异常体的水平电阻率ρ_h和垂向电阻率ρ_v分别为(ρ_h=50 Ωm, ρ_v=50 Ωm)、(ρ_h=50 Ωm, ρ_v=100Ωm)、(ρ_h=25 Ωm, ρ_v=50 Ωm).发射源沿x方向布设在坐标原点, 发射电流幅值为50 A, 发射频率为1 Hz.

图 8为山丘地形条件下低阻体模型沿测线方向的电场E_x归一化振幅响应, 图 8a振幅响应, 图 8b归一化振幅响应, 图 8c相位响应.电场E_x归一化振幅为存在各向异性异常体与不存在各向异性异常体时的电场振幅之比.图中蓝线为各向同性异常体的电场响应, 红线为异常体的水平电阻率不变, 垂向电阻率变化的电场响应, 绿线为异常体垂向电阻率不变, 水平电阻率变化的电场响应.从图可知, 低阻异常体的垂向电阻率的变化对E_x的影响很小, 而水平电阻率的变化对E_x影响明显.这是由于在低阻异常体区域, E_x主要沿垂向方向, 受水平电阻率的影响.

图 9为山丘地形条件下高阻体模型沿测线方向的电场E_x归一化振幅响应, 图 9a振幅响应, 图 9b归一化振幅响应, 图 9c相位响应.电场E_x归一化振幅为存在各向异性异常体与不存在各向异性异常体时的电场振幅之比.图中蓝线为各向同性异常体的电场响应, 红线为异常体的水平电阻率不变, 垂向电阻率变化的电场响应, 绿线为异常体垂向电阻率不变, 水平电阻率变化的电场响应.从图可知, 高阻异常体的垂向电阻率的变化对E_x的影响很大, 而水平电阻率的变化对E_x几乎没有影响.这是由于在高阻异常体区域, E_x主要沿水平方向, 受垂向电阻率的影响.

6 结论

(1) 从带源的麦克斯韦方程组出发, 基于总场推导了考虑电导率和介电常数各向异性、适应复杂地形的全张量2.5维微分控制方程, 并采用Galerkin法推导了有限元方程, 通过单元分析、总体合成、直接解法求解方程组, 反傅里叶变换等一系列步骤实现了各向异性介质可控源电磁法2.5维数值模拟, 并通过理论模型验证了其正确性, 为下一步开展任意复杂地形、任意各向异性介质2.5维反演和精细化处理解释提供有利条件.

(2) 在各向异性介质可控源电磁法2.5维数值模拟中, 通过网格加密收缩方法、行波分解法边界条件的加载以及基于形函数的任意采样傅里叶反变换的应用, 有效地保证了数值模拟的精度和速度.

(3) 通过分析典型高阻油气藏模型的各向异性介质可控源电磁响应, 与已有算法的研究规律一致, 进一步表明本文算法的可靠性和稳定性.在此基础上, 研究了各向异性倾角对不同电磁场分量的影响规律.最后, 针对山丘条件下的典型高、低阻模型, 分析了垂向电阻率和水平电阻率对电磁场响应的影响规律.

附录A Galerkin法推导有限元方程

式(14)—(15) 余量分别为:

(A1)

(A2)

令余量在整个积分区域上的加权积分为零, 即:

(A3)

(A4)

其中, N_i^e为第e个单元第i个节点的基函数, Ω为积分区域, De为第e个单元, Ne为剖分单元总个数.

已知格林公式为

(A5)

(A6)

其中n_x和n_z为边界外法线的方向余弦, l为计算区域外边界.

将式(A5)—(A6) 代入式(A3)—(A4), 可得：

(A7)

(A8)