0 引言

激发极化法是电法勘探的经典方法之一，可分为时间域激电法（TDIP）和频率域激电法（FDIP），在金属矿产、非金属矿产、油气等资源、能源勘探中获得广泛而有效的应用（傅良魁，1983；李金铭，2005；李帝铨等，2007；郑冰，2015；向葵等，2016；赵荣春等，2016；孙仁斌等，2017；赵后越，2017；景强，2018）。近年来，激发极化法还不断拓展到水文地质、环境地质监测与应用等领域（Revil et al，2010；Okay et al，2013；Schwartz et al，2014；Boaga，2017）。在金属矿勘探中，激发极化法因其独特的激发极化特性，比直流电阻率方法更有效，成为金属矿勘探中不可或缺的方法。由某些矿物引起的激电异常通过电阻率难以分辨，而激发极化法与电阻率的结构特点无关，利用该方法则可观测到较强的充电率特征。在实际矿产资源勘探中起伏地形难以避免，与其他电磁勘探方法相比，利用激发极化方法观测到的异常理论上基本不受地形影响，这是其众所周知的另一大优势，当然这只是在地下电阻率、极化率各向同性情况下导出的结论（傅良魁，1983；李金铭，2005）。

激发极化法与电磁法有着密切的联系，其数值模拟主要基于“等效电阻率法”的时间域激发极化进行（Pelton et al，1978）。基于稳定电流场方程，根据极化率和电阻率已知的公式可求得等效电阻率，利用等效电阻率去替换初始电阻率，可求得极化总场，一次场可由解析公式求得，然后求得时间域极化响应。频率域激发极化的现有模型如Cole-Cole模型和Dias模型（Dias，1968）。Cole-Cole模型和Dias模型描述了复电阻率，该复电阻率取决于频率，通过引入其他参数，如时间常数、频率相关系数以及充电率，可以得出有关极化材料详细特征的结论，可由含电容器的等效电路表示。

在现有理论中，只基于各向同性介质对复杂地形影响下的激发极化法进行了研究。然而，地下介质的各向异性是客观存在的，一些常见岩石具有明显的各向异性。Linde等（2004）在石灰岩地区测量的各向异性系数高达3.7和4.5。Schmutz等（2000）发现，直流电法和瞬变电磁法在资料反演解释中必须引入较大的各向异性系数。Herwanger等（2004）指出，如果采用各向同性井间成像，将产生无法解释的层状条带图案。Li等（2013）对任意各向异性介质中的可控源电磁法响应进行了自适应二维有限元分析，结果表明，如果忽略各向异性，资料解释会导致储层位置的错误估计。因此，对于激发极化法而言，各向异性的研究具有重要意义。

目前所知，仅Liu等（2017）基于各向异性介质进行了激发极化法的数值模拟，结果表明，在不考虑电磁效应的情况下，实际解释中必须考虑各向异性。然而，其只计算了简单各向异性介质的响应，未考虑地形的影响，难以满足激发极化法勘探的需要。目前，各向异性和地形对激发极化资料解释的综合影响尚未见报道。本研究旨在发展一种适用于带地形、非结构有限元三维各向异性激发极化的数值模拟方法。

本文主要介绍了激发极化的基本方程、边值问题、TDIP和FDIP响应，讨论了各向异性和地形条件下的激发极化响应，并进行规律总结。

1 激发极化模型正演 1.1 基本方程和边值问题

在地球导电介质中，电势的基本方程和边值问题需在地空表面满足Neumann边界条件，在无穷远边界满足混合边界条件。公式如下

$\nabla \cdot(\sigma(\boldsymbol{r}) \nabla u(\boldsymbol{r}))=-I \delta\left(\boldsymbol{r}-\boldsymbol{r}_0\right)$

(1)

$\left.(\partial u(\boldsymbol{r}) / \partial \boldsymbol{n})\right|_{\varGamma_0}=0$

(2)

$\left.\left[\partial u(\boldsymbol{r}) / \partial \boldsymbol{n}+\cos \left(\boldsymbol{r}-\boldsymbol{r}_0, \boldsymbol{n}\right) \partial u(\boldsymbol{r}) / \partial \boldsymbol{r}\right]\right|_{\varGamma_{\mathrm{R}}}=0$

(3)

式中，σ表示电导率张量，δ为狄拉克三角函数，∇为拉普拉斯算子，r是观测位置，I为电流强度，Γ₀是地空边界，Γ_R是无穷远边界，r₀为电流源位置，n表示法向分量。

1.2 空间离散化

利用Galerkin方法对基本方程和边值问题进行简化（Ciarlet，2002），公式如下

$F(u)=\int_{\varOmega}\left[\frac{1}{2} \sigma(\nabla u)^2-2 I \delta\left(\boldsymbol{r}-\boldsymbol{r}_0\right) u\right] \mathrm{d} \varOmega+\frac{1}{2} \int_{\varGamma_R} \sigma \frac{\cos (\boldsymbol{r}, \boldsymbol{n})}{\boldsymbol{r}} u^2 \mathrm{~d} \varGamma$

(4)

式中，Ω为目标区域。通过模型离散化，电势可以用插值函数u = U_e^T N近似表示，公式如下

$F(u)=1 / 2 \sum _{\mathrm{e}}^{N_{\mathrm{e}}} \boldsymbol{U}_{\mathrm{e}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{e}}^1 \boldsymbol{U}_{\mathrm{e}}-I u\left(\boldsymbol{r}_0\right)+1 / 2 \sum _{\varGamma_R} \boldsymbol{U}_{\mathrm{e}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{e}}^2 \boldsymbol{U}_{\mathrm{e}}$

(5)

式中，U_e为单元节点电势向量，N为形函数，N_e为单元总数目，和分别为体单元矩阵和面单元矩阵。将所有局部矩阵组合成一个全局系统方程，公式如下

$\boldsymbol{A} \boldsymbol{U}(\boldsymbol{r})=\boldsymbol{B}$

(6)

式中，$\boldsymbol{A}=\sum\left(\boldsymbol{K}_{\mathrm{e}}^1+\boldsymbol{K}_{\mathrm{e}}^2\right)$，为一个全局对称的系统矩阵；B表示电流源项；U(r)为待求数值的未知电势向量。

1.3 时间域激发极化响应

极化率是时间域激发极化法的核心参数，极化率的定义（Soueid Ahmed et al，2018）为

$\eta=1-\frac{\sigma_0}{\sigma_{\infty}}$

(7)

式中，η∈[0, 1]；σ₀为直流电法中介质的原始电导率，与初始电阻率ρ₀成反比，即ρ₀=1/σ₀。因此，极化率可表示为

$\eta=1-\frac{\rho_0}{\rho_{\infty}}$

(8)

由于文中研究的是三维各向异性介质，故等效电阻率ρ_∞(r)、初始电阻率ρ₀(r)和极化率η(r)均为三阶张量，则

$\boldsymbol{\rho}_{\infty}(\boldsymbol{r})=\boldsymbol{\rho}_0(\boldsymbol{r}) /(1-\boldsymbol{\eta}(\boldsymbol{r}))$

(9)

$\boldsymbol{\eta}(\boldsymbol{r})=\left(\begin{array}{ccc}\eta_x & 0 & 0 \\ 0 & \eta_y & 0 \\ 0 & 0 & \eta_z\end{array}\right)$

(10)

时域激发极化模型需要计算2次，利用等效电阻率法可实现时域激发极化数值模拟（Seigel，1959）。视极化率η_s(r)可表示为二次电势U_s(r)与总电势U(r)之比，公式如下

$\eta_{\mathrm{s}}(\boldsymbol{r})=U_{\mathrm{s}}(\boldsymbol{r}) / U(\boldsymbol{r})=\left(U(\boldsymbol{r})-U_{\mathrm{p}}(\boldsymbol{r})\right) / U(\boldsymbol{r})$

(11)

因此，视极化率可以通过一次电势和总电势来求解。

1.4 频率域极化响应

在FDIP中，当介质为极化时，可以用Cole-Cole模型在多频率下的复电阻率ρ^*(ω)来描述介质的各个部分（Pelton et al，1978），公式如下

$\rho_j^*(\omega)=\rho_{j 0}\left\{1-m_j\left[1-\frac{1}{1+\left(i \omega \tau_j\right)^{c_j}}\right]\right\}$

(12)

式中，ρ_j0为第j种介质的电阻率；m_j为充电率；τ_j表示体积极化效应的时间常数，通常在10^-2—10² s之间变化，是描述岩石结构的敏感指标；c_j为频率相关系数，数值区间为0.5—1.0。

当介质表现为各向异性时，Cole-Cole模型参数也可以表示为张量，因此可用于计算各向异性介质在不同频率下的复电阻率响应ρ^*(ω)，公式如下

$\boldsymbol{\rho}^*(\omega)=\boldsymbol{\rho}_0\left\{1-\boldsymbol{m}\left[1-\frac{1}{1+(i \omega \boldsymbol{\tau})^c}\right]\right\}$

(13)

FDIP中的频率效应F_e被定义（Yang et al，2008）为

$F_{\mathrm{e}}=\left(\rho_{\mathrm{L}}-\rho_{\mathrm{H}}\right) / \rho_{\mathrm{H}}$

(14)

式中，ρ_L为低频下确定的视电阻率，ρ_H为高频下确定的视电阻率。

2 算法验证

有学者提出一个验证各向异性介质数值的近似结果，将$\rho_n=\sqrt{\rho_{\mathrm{L}} \rho_{\mathrm{T}}}$以及λz替换ρ和z，并给出两层各向异性模型平均视电阻率的解析公式（Pekşen et al，2014），公式如下

$\begin{aligned} \rho_{n \mathrm{a}}= & \rho_{n 1}\left\lfloor 1+2 p \sum\nolimits_{m=1}^{\infty} K_n^M\left(\frac{1}{\left(r_1^2+4 m^2 \lambda^2 z^2\right)^{1 / 2}}-\right.\right. \\ & \left.\left.\left.\frac{1}{\left(r_2^2+4 m^2 \lambda^2 z^2\right)^{1 / 2}}-\frac{1}{\left(r_3^2+4 m^2 \lambda^2 z^2\right)^{1 / 2}}+\frac{1}{\left(r_4^2+4 m^2 \lambda^2 z^2\right)^{1 / 2}}\right)\right]\right)\end{aligned}$

(15)

式中， $K_n=\left(\rho_{n 2}-\rho_{n 1}\right) /\left(\rho_{n 2}+\rho_{n 1}\right)$ ；$p=\left\{\left[1 / r_1-1 / r_2\right]-\left[1 / r_3-1 / r_4\right]\right\}^{-1}$；ρ_n1和ρ_n2分别表示第一层和第二层的平均电阻率；λ为各向异性系数；z为第一层厚度；r_i（i = 1, 2, 3, 4）为电流与电位电极之间的距离。

为了验证算法的正确性，使用一个两层各向异性模型，参数见图 1所示。该两层水平层状各向异性（VTI）模型具有一个等效的各向同性层状模型，其具有解析解（Pekşen et al，2014）。在图 1所示等效的各向同性层状模型中，第一层厚度λh，各向同性介质平均电阻率为100 Ω·m；第二层为均匀半空间，各向同性介质平均电阻率为10 Ω·m。使用对称四极装置进行测量，测线最大距离为180 m，电位电极在地表上间隔2 m分布。将两层各向异性模型的数值模拟结果与两层各向同性模型的解析解结果进行对比，最大相对误差小于0.5%（图 2），从而验证了算法的准确性。

3 模型算例 3.1 时间域激发极化模型

图 3所示为在平坦地形模型下方有一个深度为150 m的正方体异常，研究各向同性（λ = 1）、水平层状各向异性VTI（λ = 2, 4）和垂直层状各向异性TTI（λ = 2, 4）条件下，该异常体的时间域激发极化响应结果。

设模型尺寸为5 000 m×5 000 m×5 000 m，立方体异常位于x—y平面原点下方，埋深150 m，大小为100 m×100 m×100 m。该模型的背景围岩（图 3中A）电阻率为100 Ω·m，各向同性极化率为η₁ = η₂ = η₃ = 0.01 mV/V。$\lambda=\sqrt{\rho_{\mathrm{T}} / \rho_{\mathrm{L}}}$，λ表示各向异性系数，ρ_L表示横向或平行层理方向电阻率，ρ_T表示纵向或垂直层理方向电阻率，异常体参数见表 1。

在地表采用对称四极装置测视极化率，电位电极间距5 m，电流源电极最大距离为1 000 m，最小距离为10 m，电流强度为1 A。以各向同性异常为例，设初始电阻率为10 Ω·m，极化率为0.2 mV/V，计算得到异常体等效电阻率为12.5 Ω·m，即可用等效电阻率来计算三维有限元建模中的总电势U(r)。

图 4(a)—(c)显示了由对称四极装置获得的视极化率剖面图，在平坦地形中，当异常的电阻率和极化率分别表现为各向同性和各向异性时，在时间域激发极化响应中的视极化率剖面表现不同。当异常为各向同性时（λ= 1），观测的异常体视极化率约1%，与背景围岩相近。而当VTI异常为各向异性时（λ= 2），观测的最大视极化率约为1.6%；当VTI异常表现出较强的各向异性时（λ = 4），观测的最大视极化率大于2.2%。与各向同性结果[图 4(a)]相比，各向异性结果[图 4(c)]具有更清晰的极值中心，异常位置定位更准确。由图 4(d)—(e)可见，将η₂和η₃的值进行交换后，对于任意各向异性TTI介质，模拟结果均无法定位异常位置。由图 4(f)可见，在各向异性VTI（λ = 4）条件下，视极化率对异常体的反应最强；TTI介质的视极化率在背景极化率值（1%）附近波动，说明在实际激发极化勘探中，在地表观测不到实际存在于地下的TTI各向异性异常；同时，VTI介质的各向异性越强，异常定位越容易，而TTI介质在平坦地形下也很难被检测到。

3.2 频率域激发极化模型

为了探讨地形对频率域激发极化模型的影响，设计如图 5的示例模型。图 5中平坦地形、台地地形、地堑地形的球体异常模型被用来研究FDIP响应结果。模型尺寸为1 000 m×1 000 m×500 m。在地表采用中间梯度装置观测，电位电极分布范围为200 m×480 m，测点间隔10 m，电流电极间距800 m。球体异常半径10 m，位于x—y平面中心，异常顶部距离地表为30 m。背景围岩参数如下：电阻率ρ_x = ρ_y = ρ_z = 100 Ω·m，充电率为m_x = m_y = m_z = 0.05，时间常数为τ_x = τ_y = τ_z = 1 s，频率相关系数为c = 0.25；球体异常参数如下：电阻率ρ_x = ρ_y = ρ_z = 10 Ω·m，充电率为m_x = m_y = m_z = 0.05，时间常数为τ_x = τ_y = τ_z = 1 s，频率相关系数为c = 0.25，2个计算频率分别为2 Hz和2/13 Hz。平坦地形、台地地形和地堑地形的频率效应等值线见图 6。

由图 6可见，可在地面观测到球体异常的频率效应响应，以及异常体正上方地表中心位置出现椭球形的频率效应F_e最大区域，该区域与图中黑色虚线表示的异常的真实位置较为一致，根据该频率效应值可以判断异常体为高极化异常；台地地形的最大频率效应最弱（红线），地堑地形的最大频率效应最强，平坦地形最大频率效应则介于二者之间[图 6(d)]。因此，台地地形对真实异常的响应存在减弱效果，地堑地形对真实异常的响应存在增强效果，但3种地形对异常真实位置的反映均较为准确。

4 结论

在金属矿产勘查过程中，地下极化介质分布具有任意性和各向异性特征。采用非结构有限元法建立了考虑各向异性的三维激电模型。在此基础上，首次研究了各向异性对激电介质的综合影响。在TDIP结果中，真实介质中各向异性的变化对观测到的地表激电数据有较大影响。VTI异常在平坦地形下容易被准确定位，各向异性越大，信号越强，定位越准确。而TTI介质即使在平坦地形下也很难被探测到。在实际勘探中，在激发极化数据反演工作中应考虑介质的各向异性特征。在FDIP结果中，台地地形对异常结果响应存在减弱作用，地堑地形对异常响应则存在增强效果。