基于有限差分法的2.5D频谱激电法正演模拟

Research on 2.5-Dimensional Forward Modeling of Spectral Induced Polarization Based on Finite Difference Method

摘要: 随着地球物理勘探技术的不断进步，电法勘查作为一种对地下介质进行识别的重要手段而广泛应用于资源勘查、环境监测和工程勘察等领域。谱激电法作为电法勘查中的常用方法，由于其在探测地下低阻异常体和复杂地质体方面的独特优势，逐渐成为地球物理勘探中的核心技术之一。在谱激电法的正演模拟中，有限差分法能够通过网格剖分精确地求解地下介质各处的电位值，处理各种介质和几何结构，提供高度精确的模拟结果。该研究旨在结合激发极化法与有限差分法，深入探讨基于有限差分法的谱激电2.5维正演模拟的实现与应用。文章推导了有限差分法的相关公式，得到网格各节点所对应的边界条件与差分方程，并用Fortran90进行编程实现。通过建立不同的地电模型并对其进行正演模拟，通过得到的数据绘制等值线断面图，通过观察均可较好的模拟出所设置的地电模型形态特点。在此基础上修改对应参数，揭示了各电性参数对模拟结果的影响，为实际电法勘探和数据解释提供理论依据和技术支持。

Abstract: With the continuous advancement of geophysical exploration technology, electrical exploration, as an important means of identifying underground media, is widely used in fields such as resource exploration, environmental monitoring, and engineering exploration. Spectral induced polarization method, as a commonly used method in electrical exploration, has gradually become one of the core technologies in geophysical exploration due to its unique advantages in detecting underground low resistance anomalous bodies and complex geological bodies. In the forward simulation of spectral induced polarization method, the finite difference method can accurately solve the potential values of underground media at various locations through grid meshing, handle various media and geometric structures, and provide highly accurate simulation results. This study aims to explore in depth the implementation and application of 2.5-dimensional forward simulation of spectral induced polarization based on the finite difference method, by combining the induced polarization method and the finite difference method. This article derives the relevant formulas of the finite difference method, obtains the boundary conditions and difference equations corresponding to each node of the grid, and implements them using Fortran90 programming. By establishing different geoelectric models and conducting forward simulations, contour section maps can be drawn based on the obtained data, and the morphological characteristics of the geoelectric model can be well simulated through observation. On this basis, the corresponding parameters were modified to reveal the influence of various electrical parameters on the simulation results, providing theoretical basis and technical support for practical electrical exploration and data interpretation.

1. 引言

谱激电法是一种可以通过测量介质在不同频率下的复电阻率响应以表征地下电化学特征，从而对地下异常体、埋藏物或矿体进行识别，从而满足生产需要的一种地球物理方法[1]-[5]。时间域激发极化法属于时间域激电法的衍生技术，可通过时–频转换间接获取频谱特征信息，相比于频谱激电法，前者存在数据质量受限于时间采样密度与噪声水平，且分辨率较低的问题，而后者能获取更为丰富的极化谱信息，更好地识别不同矿物的电化学特性，增强对地下介质成分及形态的识别能力[6] [7]。2.5维正演研究将地下介质视为二维地质体，而场源仍保持三维特性，即实际在空间中呈现三维分布[8]。其具有网格节点数多从而计算效率高、适应实际地质构造、数据解释精度高以及可组合使用提高模型准确率的特点，适用于电法勘查中。

有限差分法作为一种将微分方程转化为差分方程，用以求解各种偏微分方程的数值求解方法之一，广泛应用于流体力学、声波传播、地震模拟及传热分析等方面[9]-[12]。通常在结构化网格上实施，网格节点规律排布，一般为矩形或者六面体(对应二维及三维有限差分法网格剖分)，无需进行复杂的几何计算，且可模拟非均匀介质和复杂地形，应用混合边界条件，可控制误差在1%以内[13]，同时可以与其他数值方法或现在技术进行融合，包括时频转换、多尺度建模等方向[14] [15]。在电法勘探中，使用有限差分法进行计算，可以最大程度上简化计算步骤，减小运算量，压缩计算时间，在保障计算效率的同时拥有较高精度，故本研究将其与频谱激电法相结合，拟对水平地层及地下介质进行模拟分析，从而实现复杂地质条件下极化体响应的高精度模拟，加深对频谱激电法的理解与应用。

2. 有限差分法理论

2.1. 激发极化法

激发极化法的机理目前有几种公认的假说[16]，对于电子导体，包括大多数金属矿、石墨矿及其矿化岩石等的激发极化机理问题，被认为是电子导体与其周围溶液的界面上发生电位的结果。在向地下供入稳定电流的情况下，岩石、矿石与其中所含水溶液进行复杂的化学反应，从而产生了测量电极间的电位差随时间变大，一段时间后数值趋于饱和，而停止供电后测得的电位差在极短时间内极速下降，产生陡崖式落差，随后继续随时间延长逐渐减小为零，恢复初始状态的现象，这种现象被称为激发极化效应[17]。在此基础上引入交变电流激发的激电效应，后被称为“频率域”激电效应。

将稳定的直流电更改为超低频信号发生器，通过保持交变电流$\tilde{I}$ 的幅值不变，依次改变其频率f，由此即可得到测量电极之间交变电位差$\Delta \tilde{U}$ 与频率f之间的关联所在。一般频率的取值范围在$n\times {10}^{-2}$ 到$n\times {10}^2$ Hz间，此时电场随频谱的变化原因为岩矿石的激发极化，而不是介电极化和电磁耦合。由此可以得到，在一定电流密度条件下，$\Delta \tilde{U}$ 与$\tilde{I}$ 呈现线性关系，则可得到交流电阻率为

$\tilde{\rho }=K\cdot \Delta \tilde{U}/\tilde{I}$

式中，K为装置系数。由于交变电位差$\Delta \tilde{U}$ 与频率f之间存在联系，则交流电阻率是也随频率f的变化而变化，即为关于频率f的复变函数，因此交流电阻率$\tilde{\rho }$ 也被称为复电阻率，即$\rho \left(i\omega \right)$ 。

为了方便与无激电效应时介质的真电阻率进行分辨，将发生体极化效应的电阻率成为——“等效电阻率”[2]，其随频率或充电时间的变化而变化，在充电时间趋近于0或频率趋近于无穷时，总电位差趋于无激电效应的一次电位差$U_1$ ，此时的等效电阻率就等于介质的真电阻率$\rho$ ；反之当在充电时间趋近于无穷或频率趋近于0时，总电位差则趋于饱和电位，此时的等效电阻率被称为“极限等效电阻率”，即电阻率${\rho }^{*}$ 。其与真电阻率之间存在的关系为

${\rho }^{*}=\frac{\rho }{1-\eta }$

描述频率域激电效应的参数有三种，包括复电阻率频谱、频散谱以及相位[18]。本研究中的主要研究对象之一——频谱激电法，则是通过改变频率大小从而实现对视复现祖率的实分量、虚分量或振幅和相位的频率观测，进而对地下地质情况进行分析研究。交变电流场中的激电效应与稳定电流场中的激电效应有对应关系，则可利用高频$f_D$ 和低频$f_G$ 两个频率下的总电位差幅值$\left|\text{Δ}\tilde{U}\left(f_D\right)\right|$和$\left|\text{Δ}\tilde{U}\left(f_G\right)\right|$对频散率进行计算，得到的数据也叫百分频率效应或视幅频率[17]，如下式所示，

$P\left(f_D,f_G\right)=\frac{\left|\text{Δ}\tilde{U}\left(f_D\right)\right|-\left|\text{Δ}\tilde{U}\left(f_G\right)\right|}{\left|\text{Δ}\tilde{U}\left(f_G\right)\right|}\times 100\text{\%}$

以此来表示频率域激电效应的强弱。

2.2. Cole-Cole模型

W H Pelton et al(1978)依据对岩、矿石基本结构的认识及对大量岩、矿石标本与露头的测量结果，于1978年提出激发极化效应引起的复电阻率频谱可用以下Cole-Cole模型[19]表示

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

式中，${\rho }_0$ 为频率为零时的电阻率；m为极化率(或充电率)，相当于时间域的极限极化率，这两个参数用以表征介质导电性和激电效应强弱的参数(即强度参数)，且${\rho }_0=\rho /\left(1-m\right)$ ；$\omega$ 为角频率；$\tau$ 为时间常数，具有时间量纲，决定频率特性曲线最大斜率对应的频率；c为频率相关系数，决定频率特性曲线的变化快慢，无量纲。下图1为基于Cole-Cole模型绘制的各分量频谱曲线图。

Figure 1. Spectral curves of each component of the Cole-Cole model

图1. Cole-Cole模型各分量的频谱曲线图

3. 有限差分法原理

地球物理正演的数值计算方法多种多类，常用的主要是：有限差分法、有限单元法、边界单元法和积分方程法这几种，也可通过将这几种方法联合起来进行数值模拟计算的混合算法。本次研究中采用的是有限差分法方法，本节主要内容为直流电阻率法2.5维有限差分法理论。

3.1. 基本微分方程

从稳定电流场中电位满足的基本方程出发，可得到点源场满足的微分方程

$\frac{\partial }{\partial x}\left(\sigma \frac{\partial U}{\partial x}\right)+\frac{\partial }{\partial y}\left(\sigma \frac{\partial U}{\partial y}\right)+\frac{\partial }{\partial z}\left(\sigma \frac{\partial U}{\partial z}\right)=f$ (1)

对于二维地点断面，选取y轴平行于地质体走向，则地电参数如介质的电导率和电位只随x，z方向发生变化，在y方向没有变化，此时电位满足二维偏微分方程

$\frac{\partial }{\partial x}\left(\sigma \frac{\partial U}{\partial x}\right)+\frac{\partial }{\partial z}\left(\sigma \frac{\partial U}{\partial z}\right)=f$ (2)

对等式两端进行傅里叶变换，得到

$\frac{\partial }{\partial x}\left(\sigma \frac{\partial V}{\partial x}\right)+\frac{\partial }{\partial z}\left(\sigma \frac{\partial V}{\partial z}\right)-{\lambda }^2\sigma V=f$ (3)

式中

$\sigma =\sigma \left(x,z\right),V=V\left(x,\lambda ,z\right),f=-\frac{I}{2}\delta \left(x-x_A\right)\delta \left(z-z_A\right)$

即将三维偏微分方程变为二维偏微分方程，若给定若干个$\lambda$ 值求解方程(3)，计算出$V\left(x,\lambda ,z\right)$ 后，即可进行傅里叶反变换从而计算出电位值。

3.2. 网格剖分

取用平行于x轴和z轴的直线将求解区域离散化，分割成许多小的矩形，这里的每个小矩形被称为单元，如图2。为确保精度，网格范围应该尽可能取大，划分尽可能多的单元，同时也应考虑计算难度与计算速度，根据实际研究对象与选取的模型决定网格的疏密，此处选择使用中部加密，两侧稀疏，由中部向两侧逐渐放大的网格。网格中的交点称为节点，给其设定x方向的节点编号为$i=1,2,\cdots ,n$ ，z方向的节点编号为$j=1,2,\cdots ,m$ 。记节点$\left(i,j\right)$ 到节点$\left(i+1,j\right)$ 的距离为$\Delta x_i$ ，节点$\left(i,j\right)$ 到节点$\left(i,j+1\right)$ 的距离为$\Delta z_j$ ，同样节点$\left(i,j\right)$ 代表一个它附近的网格区域$\Delta A_{ij}$ 。

对式(3)沿$\Delta A_{ij}$ 积分得到

$\begin{array}{l}-{\displaystyle {\iint }_{\Delta A_{ij}}\nabla \cdot \left(\sigma \left(x,z\right)\nabla V\left(x,k,z\right)\right)\text{d}x\text{d}z}+{\displaystyle {\iint }_{\Delta A_{ij}}{\lambda }^2\sigma \left(x,z\right)V\left(x,k,z\right)\text{d}x\text{d}{z}^{\omega }}\\ ={\displaystyle {\iint }_{\Delta A_{ij}}\frac{I}{2}\delta \left(x-x_0\right)\delta \left(z-z_0\right)\text{d}x\text{d}z}\end{array}$ (4)

代入格林公式并利用中心差分近似$\frac{\partial V}{\partial n}$ 可以得到差分方程

$\begin{array}{c}{\displaystyle {\oint }_{L_{ij}}{\sigma }_{ij}\frac{\partial V_{ij}}{\partial n}\text{d}l}=\frac{\Delta x_i{\sigma }_{i,j-1}}{2}\left(\frac{V_{i,j-1}-V_{i,j}}{\Delta z_{j-1}}\right)+\frac{\Delta z_{j-1}{\sigma }_{i,j-1}}{2}\left(\frac{V_{i+1,j}-V_{i,j}}{\Delta x_i}\right)+\frac{\Delta z_j{\sigma }_{i,j}}{2}\left(\frac{V_{i+1,j}-V_{i,j}}{\Delta x_i}\right)\\ +\frac{\Delta x_i{\sigma }_{i,j}}{2}\left(\frac{V_{i,j+1}-V_{i,j}}{\Delta z_j}\right)+\frac{\Delta x_{i-1}{\sigma }_{i-1,j}}{2}\left(\frac{V_{i,j+1}-V_{i,j}}{\Delta z_j}\right)+\frac{\Delta z_j{\sigma }_{i-1,j}}{2}\left(\frac{V_{i-1,j}-V_{i,j}}{\Delta x_{i-1}}\right)\\ +\frac{\Delta z_{j-1}{\sigma }_{i-1,j-1}}{2}\left(\frac{V_{i-1,j}-V_{i,j}}{\Delta x_{i-1}}\right)+\frac{\Delta x_{i-1}{\sigma }_{i-1,j-1}}{2}\left(\frac{V_{i,j-1}-V_{i,j}}{\Delta z_{j-1}}\right)\end{array}$ (5)

对式(4)左端第二项有

$\begin{array}{l}{\displaystyle {\iint }_{\Delta A_{ij}}{\lambda }^2{\sigma }_{i,j}{V}_{i,j}\text{d}{x}_i\text{d}{z}_j}=G{V}_{i,j}\\ ={\lambda }^2{V}_{i,j}\left(\frac{{\sigma }_{i-1,j-1}\Delta x_{i-1}\Delta z_{j-1}}{4}+\frac{{\sigma }_{i,j-1}\Delta x_i\Delta z_{j-1}}{4}+\frac{{\sigma }_{i,j}\Delta x_i\Delta z_j}{4}+\frac{{\sigma }_{i-1,j}\Delta x_{i-1}\Delta z_j}{4}\right)\end{array}$ (6)

式中$G$ 即为${\lambda }^2$ 与式中括号内多项式的乘积。

最后整理式(4)可以得到

$C_L^{ij}{V}_{i-1,j}+C_R^{ij}{V}_{i+1,j}+C_T^{ij}{V}_{i,j-1}+C_B^{ij}{V}_{i,j+1}+C_P^{ij}{V}_{i,j}=\frac{I}{2}{\displaystyle {\iint }_{\Delta A_{ij}}\delta \left(x-x_0\right)\delta \left(z-z_0\right)\text{d}x\text{d}z}$ (7)

其中

$c_{ij}^L=-\frac{\text{Δ}{z}_{j-1}{\sigma }_{i-1,j-1}+\text{Δ}{z}_j{\sigma }_{i-1,j}}{2\text{Δ}{x}_{i-1}}$

$c_{ij}^R=-\frac{\text{Δ}{z}_{j-1}{\sigma }_{i,j-1}+\text{Δ}{z}_j{\sigma }_{i,j}}{2\text{Δ}{x}_i}$

$c_{ij}^T=-\frac{\text{Δ}{x}_{i-1}{\sigma }_{i-1,j-1}+\text{Δ}{x}_i{\sigma }_{i,j-1}}{2\text{Δ}{z}_{j-1}}$

$c_{ij}^B=-\frac{\text{Δ}{x}_{i-1}{\sigma }_{i-1,j}+\text{Δ}{x}_i{\sigma }_{i,j}}{2\text{Δ}{z}_j}$

$G={\lambda }^2\left(\frac{{\sigma }_{i,j-1}\text{Δ}{x}_i\text{Δ}{z}_{j-1}}{4}+\frac{{\sigma }_{i,j}\text{Δ}{x}_i\text{Δ}{z}_j}{4}+\frac{{\sigma }_{i-1,j}\text{Δ}{x}_{i-1}\text{Δ}{z}_j}{4}+\frac{{\sigma }_{i-1,j-1}\text{Δ}{x}_{i-1}\text{Δ}{z}_{j-1}}{4}\right)$

$C_{ij}^P=-\left(C_{ij}^L+C_{ij}^R+C_{ij}^T+C_{ij}^B-G\right)$

式(7)即为网格内部节点$\left(i,j\right)$ 上的差分方程，由此得到，节点上的电位V只与其相邻节点上的V有关，相应连接系数C只与网格大小与岩石的电导率分布有关。

Figure 2. Rectangular grid partition of the solution domain

图2. 求解区的矩形网格剖分

3.3. 边界条件

已知在直角坐标系三维空间均匀介质中的点源电位随距离成反比，即

$U\left(x,y,z\right)=\frac{A}{r}$

式中$r={\left(x^2+y^2+z^2\right)}^{\frac{1}{2}}$ 为点源到测点的距离，A为常数因子，对上式两端进行傅里叶变换，得到

$\Phi \left(x,k,z\right)=A{K}_0\left(kr\right)$ (8)

式中$K_0$ 为零阶修正贝塞尔函数，点源到测点的距离$r={\left(x^2+z^2\right)}^{\frac{1}{2}}$ ，而

$\frac{\partial \Phi }{\partial n}=-Ak{K}_1\left(kr\right)\cos \theta$ (9)

式中$\theta$ 为矢径r和边界处法线n之间的夹角，$K_1$ 为一阶修正贝塞尔函数。则式(6)变为

$\frac{\partial \Phi }{\partial n}+k\frac{K_1\left(kr\right)}{K_0\left(kr\right)}\Phi \cos \theta =0$ (10)

对不同边界上的节点应采用与其对应的边界条件，式(7)将得到进一步简化。将边界节点分为地面节点、地面上左右角点、底边节点、底边上左右角点、左边界节点与右边界节点六部分，根据不同边界条件计算得到对应差分方程，其基本形式均符合式(7)。

3.4. 总体系数矩阵特点

网格上所有节点都满足式(7)差分方程，当$i=1,2,\cdots ,n$ ，以及$j=1,2,\cdots ,m$ 时，总共可以得到$m\times n$ 个线性方程，将该线性方程组写成矩阵形式，如以下形式

$CV=B$

式中

$C=\left(\begin{array}{cccc}{C}_{11}& C_{12}& \cdots & C_{1,mn}\\ C_{21}& C_{22}& \cdots & C_{2,mn}\\ ⋮& ⋮& \ddots & ⋮\\ C_{mn,1}& C_{mn,2}& \cdots & C_{mn,mn}\end{array}\right)$

$V=\left(\begin{array}{c}{V}_1\\ V_2\\ ⋮\\ V_{mn}\end{array}\right)$

$B=\left(\begin{array}{c}0\\ ⋮\\ 0\\ B_p\\ 0\\ ⋮\\ 0\end{array}\right)$

此处C为总体系数矩阵，是网格几何性质与物性分布的函数，V是待求各网格节点上的变换电位值，B为与供电电流有关的右端项，仅在供电节点p上有值，即

$B_p=\frac{I}{2}$

除供电节点处外的其余各元素均为0。

此系数矩阵为对称正定、稀疏带状和对角占优的矩阵，故可采用带状矩阵的$L{L}^{\text{T}}$ 分解法或$LD{L}^{\text{T}}$ 分解法进行求解，得到V后经过傅里叶反变换获得空间域的电位U。

4. 算法验证与实例

4.1. 算法验证

在本研究中首先对单套网格的直流电有限差分正演进行算法验证，采用对称四极装置在一套网格上进行计算，获得单点测深数据，并将得到的有限差分法计算结果与解析解进行对比，分析产生的误差，从而确定精度。单套网格中优先对地电模型进行建立，采用二层水平层状介质的地电模型，靠近地面的一层电阻率赋值为10 Ω⋅m，第二层为100 Ω⋅m，建立网格大小为252 × 100，共计25,553个节点数，对应建立的地电模型大小为300 m × 150 m。得到的数据对比结果如下图3所示。

经过图形与具体数据误差分析计算可以得到，直流电有限差分法计算结果拟合较好，误差在0.1%~0.6%之间，说明该算法精度可以得到有效保证。若需要进一步缩小误差，则需要对网格进行进一步调整，即在异常体周围进行网格的局部加密，但随之而来也需要考虑到计算效率的问题，因而需要更多的尝试与实验。

Figure 3. Finite difference simulation results of a single grid horizontal layered medium

图3. 单套网格水平层状介质有限差分模拟结果

4.2. 算法实例

在直流电有限差分法正演模拟的基础上，引入Cole-Cole模型变为频率域谱激电正演模拟。本研究在频率域谱激电算法中选取模型为不极化围岩中存在一低阻极化体，见图4，该地电模型中的各介质激电参数可见下表1。

Figure 4. Schematic diagram of low impedance polarizer model

图4. 低阻极化体模型示意图

Table 1. Induced polarization parameters of various media in the low impedance polarized body geoelectric model

表1. 低阻极化体地电模型各介质激电参数

	电阻率$\rho$ (Ω⋅m)	极化率m	时间常数$\tau$ (s)	频率相关系数c
围岩	100	0	0	0
极化体	10	0.6	0.1；1；10	0.25

对该模型中的极化体时间常数和埋深参数进行更改，从而观测其视复电阻率参数和频散率数值，并对得到的结果进行绘图，得到等值线断面图。

4.2.1. 不同埋深参数情况下的正演模拟

在保证围岩及极化体电阻率、时间常数、频率相关系数都不变的情况下，对极化体埋深进行改变，同时网格数据不变。依次对极化体埋深H取10 m、20 m、30 m，经过有限差分法正演模拟得到三组不同的视复电阻率数据，并根据这些数据对实部、虚部、视频谱的等值线断面图进行绘制。下图5为围岩电阻率为100 Ω∙m，极化率为0，极化体电阻率为10 Ω∙m，极化率为0.6，时间常数为1 s，频率相关系数为0.25，频率为1 Hz时的不同埋深参数情况下视复电阻率实部与虚部的等值线断面图。其中横坐标x为测点位置，纵坐标为供电极距，等值线数据分别为视复电阻率实部与虚部数值。

(a) (b)

(c) (d)

(e) (f)

Figure 5. Cross section of the contour lines of the real and imaginary parts of the apparent complex resistivity under different burial depths. (a) (b) are burial at a depth of 10 m respectively; (c) (d) are burial at a depth of 20 m respectively; (e) (f) are burial at a depth of 30 m respectively

图5. 不同埋深条件下视复电阻率实部与虚部等值线断面图。(a) (b)分别为埋深为10 m；(c) (d)分别为埋深为20 m；(e) (f)分别为埋深为30 m

通过对比图5，可以较为明显地观察到各等值线断面图均呈现左右两端高值、中心低值的对称分布，低阻异常值向下延伸，且供电极距越小时，上层出现复电阻率较大值，趋于围岩电阻率值，较好的可以模拟出低阻极化体的模型。随着埋深的增加，其呈现出的视复电阻率数值中低值区域的出现随着供电极距的增加而加深，且浅部随着埋深的增加呈现更趋于围岩电阻率的趋势，即得到结论为对应埋深越深，视复电阻低率值出现的区域越向下，反映出的是低阻极化体埋深位置的变动。

下图6为围岩电阻率为100 Ω∙m，极化率为0，极化体电阻率为10 Ω∙m，极化率为0.6，时间常数为1 s，频率相关系数为0.25，频率为1/13 Hz时的不同埋深参数情况下视复电阻率实部与虚部的等值线断面图。其中横坐标x为测点位置，纵坐标为供电极距，等值线数据为视频散率数值。

(a)

(b) (c)

Figure 6. Cross section of video dispersion contour lines under different frequency conditions (a) for a burial depth of 10 m; (b) for a burial depth of 20 m; (c) for a burial depth of 30 m

图6. 不同埋深条件下视频散率等值线断面图。(a) 为埋深为10 m；(b) 为埋深为20 m；(c) 为埋深为30 m

图6展示了在不同埋深条件下视频散率等值线的分布情况，图(a)、(b)和(c)分别对应埋深为10 m、20 m和30 m时的视频散率等值线剖面图。从图中可以看出，随着埋深的增加，等值线的分布逐渐趋于对称和平滑，最大值区域的位置逐渐向中心靠拢，频散率数值的集中程度有所提升，且呈现不断向下延伸的趋势。图(a)中最大视频散率位于靠近中心偏上的区域，图(b)中等值线呈现更对称的形态，最大值位于图像中心区域，图(c)中整体散率数值有所下降，最大值也较前两图略小。三幅图的对比体现了埋深对视频散率的显著影响，随着埋深的增加，高视频散率区域的强度逐渐减弱，极化体的视频散率向更深处扩展。

4.2.2. 不同频率情况下的正演模拟

频率作为研究激电效应的一个重要参数，对于同一个极化地电模型采用不同的频率进行模拟会产生较大的差别。在此基础上沿用低阻极化体模型，分别对不同埋深条件下的极化体模型采用不同的频率，分别取频率f为1 Hz，1/13 Hz，下图7为频率为1/13 Hz时视复电阻率等值线断面图。通过与频率为1 Hz时的图5进行对比，可以得到不同频率条件下的正演模拟情况。

从图7中也可以明显地观察到各等值线断面图均呈现左右两端高值、中心低值的对称分布，模拟出低阻极化体的模型程度较好，且随着埋深的增加，呈现出的视复电阻率数值中低值的出现随着供电极距的增加而加深，即得到结论为对应埋深越深，视复电阻低率值出现的越深，进一步验证图5的结论。再通过结合图5和图7，可以得到在不同频率条件下，复电阻率实部随着频率的增大，对低阻部分的模拟情况更好，而虚部随频率增大并无明显变化。

(a) (b)

(c) (d)

(e) (f)

Figure 7. Cross section of the contour lines of the real and imaginary parts of the apparent complex resistivity under different frequency conditions (a) (d) are buried at a depth of 10 m; (b) and (e) are buried at a depth of 20 m; (c) and (f) are buried at a depth of 30 m

图7. 不同频率条件下视复电阻率实部与虚部等值线断面图。(a) (d)为埋深为10 m；(b) (e)为埋深为20 m；(c) (f)为埋深为30 m

4.2.3. 不同时间常数情况下的正演模拟

时间常数$\tau$ 反映了极化体内部极化过程的特征，在影响复电阻率的频散特性中起着关键作用。下面对不同时间常数情况下的低阻极化体模型进行正演模拟，极化体依旧沿用之前的参数，下图8为围岩电阻率为100 Ω∙m，极化率为0，极化体埋深为10 m，电阻率为10 Ω∙m，极化率为0.6，时间常数分别为0.1s、1 s、10 s，频率相关系数为0.25，频率为1 Hz时的视复电阻率实部与虚部等值线断面图。其中横坐标x为测点位置，纵坐标为供电极距，等值线数据分别为视复电阻率实部与虚部数值。

观察到图(a)和图(d)展示了较短时间常数($\tau$ 为0.1 s时)下的视复电阻率断面图，表现出较为低阻异常值集中在顶部，两侧逐渐扩大至围岩电阻率的形态；图(b)和图(e)中，随着时间常数的延长至1 s，视复电阻率实部随供电极距的增大其低值出现区域收窄，整体形态无较大变化；图(c)和图(f)表示时间常数为10 s时，对低阻极化体的模拟更为明显，数值也更趋于低阻极化体的电阻率参数。而纵观三个不同时间常数的视复电阻率虚部等值线断面图，可以看到整体形态保持不变，随着供电极距的增大，低值出现区域不断收窄。

(a) (b)

(c) (d)

(e) (f)

Figure 8. Cross section of the contour lines of the real and imaginary parts of the apparent complex resistivity at a frequency of 1Hz under different time constants (a) (d) have a time constant of 0.1 s; (b) and (e) have a time constant of 1 s; (c) and (f) have a time constant of 10 s

图8. 不同时间常数情况下频率为1 Hz的视复电阻率实部与虚部等值线断面图。(a) (d)为时间常数为0.1 s；(b) (e)为时间常数为1 s；(c) (f)为时间常数为10 s

下图9为围岩电阻率为100 Ω∙m，极化率为0，极化体埋深为10 m，电阻率为10 Ω∙m，极化率为0.6，时间常数分别为0.1 s、1 s、10 s，频率相关系数为0.25，频率为1/13 Hz时的视复电阻率实部与虚部等值线断面图。其中横坐标x为测点位置，纵坐标为供电极距，等值线数据分别为视复电阻率实部与虚部数值。

观察图(a)和图(d)展示了较短时间常数($\tau$ 为0.1 s时)下的视复电阻率断面图，表现出较为低阻集中在顶部，两侧逐渐扩大至围岩电阻率的形态；图(b)和图(e)中，随着时间常数的延长至1 s，视复电阻率实部随供电极距的增大其低值出现区域扩大，整体形态无较大变化；图(c)和图(f)表示时间常数为10 s时，对低阻极化体的模拟相比较下弱化，数值也更趋于高阻围岩的电阻率参数。而纵观三个不同时间常数的视复电阻率虚部等值线断面图，可以看到整体形态保持不变，随着供电极距的增大，低值出现区域不断收窄，这一点同图7得到的结论一致。通过对比图8和图9，观察到时间常数相同时，高频与低频所对应的视复电阻率实部特征相反，即随着时间常数由0.1 s到1 s再到10 s，低频时实部数据随着供电极距越大，范围扩大且对低阻极化体的模拟程度越不明显，高频时实部数据则随着供电极距的增大而范围缩小，且从图像来看对低阻极化体的模拟越明显。

(a) (b)

(c) (d)

(e) (f)

Figure 9. Cross section of the real and imaginary parts of the apparent complex resistivity at a frequency of 1/13 Hz under different time constants. (a) (d) have a time constant of 0.1 s; (b) and (e) have a time constant of 1 s; (c) and (f) have a time constant of 10 s

图9. 不同时间常数情况下频率为1/13 Hz的视复电阻率实部与虚部等值线断面图。(a) (d)为时间常数为0.1s；(b) (e)为时间常数为1 s；(c) (f)为时间常数为10 s

下图10为围岩电阻率为100 Ω∙m，极化率为0，极化体埋深为10 m，电阻率为10 Ω∙m，极化率为0.6，时间常数分别为0.1 s、1 s、10 s，频率相关系数为0.25，频率分别为1 Hz与1/13 Hz时的不同时间常数条件下的视频散率。其中横坐标x为测点位置，纵坐标为供电极距，等值线数据为视频散率数值。

(a)

(b) (c)

Figure 10. Cross section of video dispersion contour lines under different time constant conditions. (a) (d) have a time constant of 0.1 s; (b) and (e) have a time constant of 1 s; (c) and (f) have a time constant of 10 s

图10. 不同时间常数条件下的视频散率等值线断面图。(a) 为时间常数为0.1 s；(b) 为时间常数为1 s；(c) 为时间常数为10 s

在图(a)中时间常数为0.1秒时，视频散率值主要集中在较浅的深度区域，且异常区的响应比较集中和明显。图(b)中，随着时间常数的延长至1秒，散率异常区的形态发生了变化，开始出现更大的视频散率数值，深层介质的响应开始变得更加明显。到了图(c)中的时间常数为10秒时，视频散率的异常区域进一步扩展，并且在深度方向上呈现出更加广泛的衰减。同时通过对比这三个断面图，可以看到围岩的视频散率数值逐渐向上层延伸。长时间的激发使得地下更深层次的介质也产生了响应，但这种响应由于深度较大而逐渐衰减。图9展示了在对低阻极化体模型进行正演模拟时，随着时间常数的增加，对极化体的形态模拟得越好。

5. 结论

1) 在埋深方面，模拟结果表明，随着埋深的增加，视复电阻率的低值区域向更深的深度延伸，且数值逐渐向围岩的电阻率趋近。这一趋势与直流电与时间域激电正演模拟中的结果一致，进一步验证了埋深对模拟结果的显著影响。

2) 当频率发生变化时，发现复电阻率的实部随着频率的增大而更加精准地模拟了低阻部分的形态，而虚部则几乎不受频率变化的影响，保持了相对稳定的数值。这表明在频率较高时，复电阻率的实部对于低阻异常体的模拟效果更为敏感的特性。

3) 在时间常数的变化上，高阻与低阻正演模拟得到了相反的规律，即当高阻情况下，时间常数的增大有助于改善对低阻极化体的模拟效果，而在低阻情况下，时间常数减小则能够更好地模拟低阻极化体。这一现象表明，时间常数对模拟结果的影响取决于极化体的电性特征，而虚部的变化则表现得较为平稳，在大多数情况下存在向中部收拢的趋势。

基金项目

广西自然科学基金(2024GXNSFBA010257)；广西人才基地专项基金(桂科AD23026231)资助。