基于格点型有限体积法的直流电测深2.5维正演模拟

欧阳黎明，赵煊煊，童孝忠，谢维，李欢

(1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；2. 湖南省国土空间调查监测所，湖南 长沙，410129；3. 山西省煤炭地质物探测绘院有限公司，山西 晋中，030600)

直流电测深法勘探是根据地壳中不同岩矿石之间的电阻率或电导率差异，通过在地表观测天然或人工建立的电场分布来寻找不同岩体结构、调查地质构造以及解决各类工程地质问题[1]。直流电测深勘探方法在探测几十米至几百米深度的目标体时具有独特的优势，如勘探仪器操作简单、野外勘测成本低廉及定量反演解释效果较好等。为了利用直流电测深法查明地下构造及岩体的空间展布特征，这需要对直流电测深勘探数据进行定性解释与定量解释，其中包括正演数值模拟与反演成像。直流电测深法勘探的正反演解释需要以勘探地质对象的实际电性结构变化为基础。当前，直流电测深法的一维、二维正反演技术是实测资料处理与解释的主要技术[2-3]。

直流电测深正演数值模拟是基于地电模型的电性分布求解场分布，目前，成熟的直流电测深正演数值模拟方法主要有有限差分法、有限单元法和有限体积法。有限差分法(finite difference method)是一种近似的数值模拟方法，其正演计算的基本原理是采用差商来代替微商或偏微商，将偏微分方程定解问题转化为适定的线性方程组来求解，这可很好地近似表征地下介质的电性结构[4-7]。有限单元法(finite element method)是求解地球物理波场变分问题的一种数值模拟近似方法，通过将偏微分方程定解问题转换为变分方程，然后将被积分区域离散化，再按照一定的规律和需要以及异常体的响应精度进行网格离散，最终也能将偏微分方程定解问题转化为适定的线性方程组来求解[8]。有限单元法已广泛应用于直流电测深正演模拟，取得了较好的数值计算效果[9-12]。

有限体积法(finite volume method)也称控制体积法，其基本求解思路是将整个计算区域剖分为有限个连续的单元，并且在每个离散网格节点的周围都按一定的规则形成一个互不重复的控制体，最后将控制方程对每一个控制体进行积分运算，从而将偏微分方程定解问题转化为适定的线性方程组来求解[13]。有限体积法的求解思路简单，通过离散控制方程和定解条件将微分方程转化为代数方程，遵循因变量在控制体积内的守恒原理。无论是常微分方程、偏微分方程、各种类型的二阶线性方程还是高阶线性或非线性方程，均可利用有限体积法将方程求解转化为线性方程组，而后利用计算机求其数值解[14]。近年来，有限体积法在地球物理勘探中逐渐得到应用，众多学者研究了地球物理正演模拟问题[15-17]。通常，有限体积法构造控制体存在格点型(vertex-centered)和格心型(cell-centered)共2 种方式[18]。格心型方式将单一的网格单元作为控制体单元，并假设单元内电位为常数，这需要进一步近似地面-空气界面节点的直流电测深离散电位[19]，必然会给直流电测深的电位模拟带来计算误差。而格点型方式以网格节点为中心形成控制体单元，并将未知电位置于网格离散节点处，从而可容易地计算地表处的电位，减小数值计算误差。

本文从点电源电位所满足的边值问题出发，建立点源2.5维直流电测深的格点型有限体积正演算法，以便为后续的时间域激电和复电阻率法的正反演研究提供数值计算工具。通过试算存在解析解的均匀半空间地电模型和层状介质地电模型，验证格点型有限体积正演算法的准确性和稳定性，并分析直流电测深2.5维正演模拟的误差来源。利用格点型有限体积数值算法计算典型二维地电模型的视电阻率响应，并总结直流电测深的异常响应规律和特点，以便为实测数据的定性解释提供指导。

1 边值问题

直流电测深法勘探通常采用点电源供电，而点电源产生的电位具有三维分布的特征，即u=u(x，y，z)，其满足如下三维变系数泊松方程边值问题[20]：

其中：u为点电源产生的电位；δ(x)为狄拉克函数；ρ为点电源处的电荷密度；ΓS为地表-空气界面，即边值问题的上边界；Γ∞为截断边界，即无穷远边界；r为点电源S到截断边界Γ∞的距离；σ为电导率(电阻率的倒数)，其单位为S/m；cos(r，n)为矢径r与边界外法线方向n的夹角余弦。

对于二维地电模型，若地下介质的电性参数沿走向y不发生变化，即σ=σ(x，z)，为了消除走向坐标y，需要对式(1)中的偏微分方程进行傅里叶变换。由于电位u(x，y，z)是实函数，并且是变量y的 偶 函 数，即u(x，y，z) =u(x，-y，z)，所 以，对u(x，y，z)进行余弦傅里叶变换，且积分区间选择为0至+∞，则余弦傅里叶变换公式为

其中：U(λ，x，z)为空间域电位u(x，y，z)经余弦傅里叶变换后的波数域电位；λ为波数或傅里叶变换变量。

对式(1)中的偏微分方程和边界条件同时进行傅里叶余弦变换，可得波数域电位U(λ，x，z)满足如下带参数(波数λ)的二维变系数亥姆霍兹方程边值问题[21]：

式中：Q为稳恒电流密度；K1和K0分别为第二类的1阶、0阶修正贝塞尔函数。

2 有限体积正演算法

2.1 控制方程离散

采用矩形网格单元将直流电测深2.5维正演的二维计算区域进行离散，如图1所示。令

图1 二维地电模型离散化Fig.1 Discretization for two-dimensional geo-electric model

二维离散区域的内部节点示意图如图2 所示。取二维地电模型离散化后网格内部的任意一个节点，按照有限体积法的基本思想，在控制容积中对式(3)中的偏微分方程两端进行积分：

图2 二维离散区域的内部节点示意图Fig.2 Internal nodes of two-dimensional discretized region

2.2 边界条件处理

上边界节点满足第二类边界条件方程，其控制单元形成的边界为线段a—b—IV—V—VI—VII(见图3)。因此，在该线段上的积分为

图3 剖分网格上边界节点示意图Fig.3 Diagram of nodes located on the top edge of mesh

而

将式(10)和式(11)代入式(5)，可得上边界节点(不含2个角点)满足的代数方程为

其中：

下边界节点满足第三类边界条件，其控制单元形成的边界为线段I—II—III—b—a—VIII，如图4所示，则在该线段上的积分为

图4 剖分网格的下边界节点示意图Fig.4 Diagram of nodes located on the bottom edge of mesh

且

将式(13)和式(14)代入式(5)，可得下边界节点(不含两个角点)满足的代数方程为

其中：

左右边界节点均满足第三类边界条件方程，因此，对左右边界节点满足的方程，采用上面相同的方式进行处理。

对于离散区域的角点处理，以上边界右端角点为例，给出其处理办法。上边界右角点控制单元形成的闭曲线为a—d—VI—VII(见图5)。因此，在该闭曲线上的积分为

图5 剖分网格的上边界右端角点示意图Fig.5 Diagram of right corner node on the bottom of mesh

且

将式(16)和式(17)代入式(5)，可得上边界右角点满足的代数方程为

其中：

2.3 线性方程组求解

将计算区域的所有内部节点和边界节点进行格心型有限体积法处理，最终得出波数域中的线性方程组为

通过求解该适定线性方程组(19)，即得波数域中所有离散网格节点的电位U。对于多个供电电源，需要对方程(19)进行多次求解，然后对波数域电位进行傅里叶逆变换，便可得到点源激励下二维地电断面的空间域电位，并根据相应的观测装置计算视电阻率。

采用有限体积法解直流电测深正演问题，最终将形成大型、稀疏、对称的线性方程组，常采用直接法或迭代法求解。若供电点较多，则采用直接法计算，计算效率较高，但占用内存较大；若供电点较少，则采用迭代法计算，计算效率较高，并且占用内存较少[20-22]。鉴于此，根据地电模型和计算配置情况选择合理的求解方法，本文采用不完全LU 分解(即上三角与下三角分解)预处理的稳定双共轭梯度算法(ILU-BICGSTAB 迭代法)[23]，实现直流电测深2.5维正演模拟的有限体积线性方程组求解。

2.4 傅里叶逆变换及视电阻率计算

当求出若干波数域中点源二维直流电测深的电位U(λ，x，z)时，采用傅里叶逆变换就可以求出三维点源空间的电位u(x，y，z)。在剖面(y=0)上，傅里叶逆变换公式为

利用数值积分，将式(20)写成[24]

其中：N为单位长度内波周的个数，为了获得高精度的数值解，在实际计算过程中取N=10；，为剖面上的点至点电源处的距离；λi为利用最优化原理得到的离散波数[25]；gi为相应的傅里叶逆变换系数。离散波数λi及相应的傅里叶逆变换系数gi见表1。

表1 离散波数λi及相应的傅里叶逆变换系数giTable 1 Discrete wavenumbers λi and corresponding coefficients of inverse Fourier transform gi

直流电测深2.5 维正演模拟时需计算视电阻率，可以根据电位的模拟结果换算出不同装置的视电阻率，其计算公式为[26]

其中：K为测量装置系数；Δu为空间域电位差；I为点源供电电流强度，一般可取为1 A。

图6 所示为格点型有限体积法正演模拟2.5 维直流电测深响应的基本流程图，通过正演最终获得三维源、二维地电模型的视电阻率。

图6 格点型有限体积正演算法流程Fig.6 Flow-process diagram of vertex-centered finitevolume forward algorithm

3 正演算法测试与模型试算

采用数值计算的MATLAB 软件，编写直流电测深2.5维正演模拟的格点型有限体积法程序。本文的数值测试平台为CPU i5-12400，RAM 16G；正演程序编制系统为MATLAB R2021b，该系统可方便调用0阶和1阶贝塞尔函数。

3.1 算法测试

3.1.1 均匀半空间模型

为了验证有限体积正演算法和正演程序的正确性，选取1个均匀半空间地电模型，其电导率为0.01 S/m。采用格点型有限体积法模拟2 个异性点电源产生的电位，取计算区域的长度为200 m、宽度为100 m，且将求解区域离散为100×50 的网格。正负电流源布置于(-20 m，0 m)和(20 m，0 m)处，供电电流为I=1 A，将2.5 维数值计算结果与均匀半空间电位解析解[27]进行对比，如图7所示。从正演模拟结果可知：利用有限体积法得到的电位数值解与解析解拟合得较好，从而验证了正演算法和MATLAB 正演程序的正确性；另外，在供电点源附近的模拟结果存在误差，这可以通过加密点源附近网格来提高模拟精度，也可利用边界/奇异校正处理方法来提高数值计算精度[25]。

图7 均匀半空间模型的点电源电位数值解与解析解的对比Fig.7 Comparison of analytical and numerical solutions of the point source potentials in the homogeneous halfspace model

3.1.2 层状介质模型

为了进一步验证算法的准确性，构建1个水平层状地层模型，其模型参数为σ1=0.005 S/m，σ2=0.05 S/m和h1=4m，如图8所示。分别采用直流电测深一维正演算法[28]和本文建立的2.5维有限体积正演算法计算该2层水平地电模型的视电阻率。选用的测量装置为温纳测深装置，供电电极距设计为1.5～75.0 m。

图8 2层水平层状模型示意图Fig.8 Schematic diagram of two-layered horizontal model

一维正演和2.5维正演计算的视电阻率曲线对比见图9。从图9 可知视电阻率数值解与解析解较吻合，这说明本文建立的2.5维正演算法模拟的视电阻率精度满足实用要求。在点源附近和供电极距大于70 m时的误差较大，但其相对误差均在2%以内。直流电测深2.5维有限体积正演的误差主要来自3个方面：1) 供电点源附近奇异性问题带来的误差；2) 截断边界条件和有限体积离散带来的误差；3) 离散傅里叶逆变换数值积分带来的误差。

图9 2层地电模型数值解与解析解对比Fig.9 Comparison of analytical and numerical solutions for two-layered geo-electrical model

3.2 模型试算分析

在直流电测深2.5维正演模拟过程中，无穷远截断边界条件的近似处理会产生数值计算误差。当勘测点距离截断边界越近时，误差必然会越大，只有足够远的截断边界方能满足直流电测深勘探的实际要求精度。但足够远的边界处理方式又势必增加网格单元的数量，增加求解有限体积线性方程组所需的时间。为了实现无穷远边界的仿真模拟，将直流电测深计算区域剖分为目标区域和网格外延区域，进而使第三类齐次截断边界不受异常体的影响。目标区域为地质构造体的赋存部分，也是直流电测深勘探数据的地表采集区域，其剖分方式为均匀网格；而网格外延区域用于实现无穷远边界仿真，其剖分方式为非均匀网格，网格剖分的步长按大于1的倍数等比递增，在保证正演计算精度的情况下，减少网格剖分单元数，进而提高正演计算效率。

为了分析典型模型的直流电测深的视电阻率响应，构建1个二维简单异常体地电模型。在电导率为0.01 S/m的均匀半空间中，存在1个电导率为0.1 S/m 的矩形异常体(长4 m、宽2 m)，其顶部距离地表2 m，具体的模型参数如图10所示。

图10 二维地电模型示意图Fig.10 Schematic diagram of a two-dimensional geoelectrical model

研究区域x和z方向的剖分范围分别为[-35，35] m 和[-30，0] m，采用1 m 的等间隔单元进行剖分，合计剖分71×31=2 201个节点。外延区域的最小剖分单元边长为2 m，且外延节点数为10 个。选用的测量装置为温纳测深装置，电极距设置为1 m。图11(a)所示为电阻率ρs有限体积数值解等值线图。根据封闭的等值线图，可以定性判别出低阻异常体(或高导异常体)的电阻率分布特征，并能准确判别出异常体的空间展布特征。图11(b)所示为Res2dmod 软件[29]计算所得电阻率等值线图，这与有限体积法的计算结果一致，且两者的电阻率绝对误差最大值仅为3.4 Ω·m(见图11(c))。

图11 低电阻率异常体的正演模拟结果Fig.11 Forward modeling results of low-resistivity anomalous body

4 结论

1) 从点电源电位满足的变系数泊松方程边值问题出发，建立了直流电测深2.5维正演模拟的格点型有限体积算法。格点型有限体积法利用散度定理可直接将变系数亥姆霍兹方程边值问题转换成适定线性方程组来求解，原理简单，易于编程实现。

2) 通过模拟均匀半空间模型中2个异性点电源产生的电位和层状介质模型的直流电测深视电阻率响应，其结果验证了格点型有限体积正演算法和MATLAB 正演程序的准确性，这为下一步开展直流电测深2.5维反演成像提供了数值计算工具。

3) 采用格点型有限体积法对典型二维地电模型进行了试算，将试算结果与有限差分数值计算结果进行对比分析，进一步验证了正演算法的有效性。本文建立的格点型有限体积正演算法不仅完善了目前2.5 维直流电阻率测深数值模拟方法，提高了数据处理和解释精度，而且该方法具有通用性，对类似地球物理正演问题同样具有实用价值。