时间域激发极化法三维有限元正演

骆红梅

（福建省地质测绘院，福建 福州 350011）

时间域激发极化法三维正演国内外有不少的研究，大多都是在直流电正演的基础上引入等效电阻率或COLECOLE 模型来做等效计算，然后采用有限元技术来实现激发极化场的正演[1]。但是三维有限元的计算速度和庞大的内存需求一直是困扰激发极化法三维反演的主要问题[2]。国内外近几年针对以上问题发表了不少文章，主要是采用定半带宽存储稀疏矩阵，直接采用LDLT 分解法求解方程组；采用一维非零元素压缩存储模式，然后运用预条件共轭梯度法（PCG）求解方程组。

1 正演方法及验证

1.1 正演方法

本文采用三维有限元数值模拟的方法，通过详细推导三维地电场变分问题所满足的变分方程，在直流电法三维有限元正演算法[3]的基础上加入COLE-COLE 模型，利用数字滤波求解γ 函数，可实现时间域激发极化法三维正演计算。在正演计算时，采用六面体矩形网格对变分方程进行离散，发射源附近适当加密；有限元形成的大型稀疏对称方程组采用CSR 存储格式进行存储，大大节约了内存空间。考虑到现在的计算机大多都是多核多线程的，为了使计算机资源得到最大限度的利用，使求解方程组的速度得到提高，所以采用并行求解技术。而MKL 库中有许多现成的并行求解算法可以利用，其中并行求解器PARDISO，不仅采用CSR 存储格式，而且为多线程多核心并行计算，可以大大提高时域激发极化发三维正演的计算效率和速度。综上所述步骤就可以实现时域激发极化发有限元三维正演。

1.2 精度验证

首先计算两正负点源在均匀半空间中产生的场值，并与解析解做对比，来验证程序的正确性。图1 为计算点到源点距离随相对误差分布图。从图可以看出随着计算点与源点距离的减小误差逐渐增大，在距离源点4m 处，相对误差小于1.6%。所以当计算点距离源点大于4m 的地方，计算结果是正确，精度是可靠的。

图1 测点离源的距离的误差分布图

1.3 高低阻异常体模型试算

采用中梯装置，建立如图2 的高、低阻异常模型。高阻异常体参数：围岩电阻率10Ω·m、极化率0.01、时间常数3.0、频率相关系数0.1；异常体电阻率分别为20Ω·m、50Ω·m、100Ω·m、500Ω·m，极化率0.3，频率相关系数0.3，时间常数5.0。

图2 模型示意图

从图3 高阻异常视极化电阻率分布图可以看出，在视极化电阻率图上，高阻异常体在它的正上方会产生一个视极化电阻率高阻异常，且随着异常体与围岩电阻率差异的增大，它的视极化电阻率异常幅值越大。

图3 高阻异常视极化电阻率分布图

低阻异常体与围岩参数：围岩电阻率200Ω·m、极化率0.01、时间常数3.0、频率相关系数0.1；异常体电阻率10Ω·m、50Ω·m、100Ω·m、150Ω·m，极化率0.3、频率相关系数0.1、时间常数分别为3.0。从图4 低阻异常视极化电阻率分布图可以看出，低阻高极化异常体在它的正上方产生低的视极化电阻率异常，且随着异常体与围岩的电阻率差异越大，它的视极化电阻率异常越明显，与高阻异常有相同的规律。

高、低阻模型的时间域有限元三维模型的响应特征说明本文采用的正演方法是合理的，有效的。

图4 低阻异常视极化电阻率分布图

1.4 速度对比

为了对比程序速度上的优势，特意编写了用定半带宽存储，直接采用LDTD 分解法求解方程组的程序和用CSR 存储，用预条件共轭梯度法（PCG）解方程的程序来作对比。在同一台计算机（双核四线程，内存2G）上，对三个程序的计算速度进行对比，见表1。

表1 不同计算方法CPU 计算时间对比表

可以看出本文所采用的CSR 存储格式和Pardiso 解方程的方法比其他两种方法的计算速度高。相比定半带宽存储模式，内存需求少；相比PCG 法，计算速度优势明显，且随着剖分节点数和计算机核数的增加，计算速度优势越明显。

2 结论

通过详细推导三维地电场变分问题所满足的变分方程，并对方程进行有限元离散分析，得到要求解的大型稀疏对称方程组。采用一维非零元素压缩存储的CSR 模式和PARDISO 并行求解器求解方程。设计了均匀半空间模型，用解析解与数值解相互拟合，证明了该方法是正确的、精度是可靠的。采用中梯装置试算了高、低阻异常模型，计算了异常体的视极化电阻率，证明该方法是合理的、有效的。最后将三种计算方法的速度进行了比较，在保证精度和较小的存储需求的前提下，本文的方法显著地提高了计算速度，为解决激发极化法三维反演问题提供了一定的基础。