电法测井的三维有限差分正演
陈汉波, 吕琴音
桂林理工大学 地球科学学院, 广西 桂林 541004

作者简介: 陈汉波(1990-),男,在读硕士研究生,主要从事电磁场数值模拟与反演研究。E-mail:563179536@qq.com

摘要

研究了电法测井不同测量方式(井—地,地—井,井—井)下点源场井中电法的三维有限差分数值模拟。采用六面体网格剖分方式来对模型进行剖分,运用一维非零元素行压缩存储模式来存储系数矩阵,减少了内存需求和计算量;采用不完全Cholesky共轭梯度(ICCG)方法来求解线性方程组,提高了求解效率;编制了相应的程序实现了井—地、地—井、井—井、倾斜井条件下的电法测井三维有限差分数值模拟。设计的算例结果验证了该算法的正确性和效率性,并且分析了各种情形下的异常特征,为进一步的反演工作打下了基础。

关键词: 电法测井; 有限差分; 不完全Cholesky共轭梯度迭代法; 一维非零元素行压缩
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2015)03-0580-05
Three-dimensional forward of electrical well logging with the finite difference method
CHEN Han-Bo, LYU Qin-Yin
College of Earth Sciences, Guilin University of Technology, Guilin 541004, China
Abstract

The electrical well logging is one of the high efficiency methods to explore deep mineral resources. The study of efficient and high-accuracy simulations is of important theoretical and practical significance. In this paper, 3D forward of electrical well logging was studied with different measurement methods (surface-borehole, borehole-surface, borehole-borehole), the one-dimensional none zero element row compressed storage method was used to store the coefficient matrix and, as a result, the memory requirements and computation were reduced. Using the ICCG method, the authors solved the linear system of equations. The solving efficiency was improved. The corresponding program was developed to achieve different measuring ways (surface-borehole, borehole-surface, borehole-borehole, tilted borehole) to study electrical well logging forward. The results of modeling have verified the correctness of the algorithm and the efficiency. The anomaly characteristics under various conditions are analyzed, which has laid the foundation for further inversion work.

Keyword: electrical well logging; finite difference method; ICCG method; CSR

传统的电法勘探一般都在地表上进行, 而采用井中电法勘探可以加大勘探的深度, 提高分辨率。目前, 该方法主要应用于深部矿产勘查, 在油气资源、水资源以及工程勘察中也有广泛的运用。对采用有限差分法进行井中电法的数值模拟, 已有学者对其展开了研究, 也取得了一些成果[1, 2, 3, 4, 5, 6, 7], 这些研究性工作对于井中电法勘探的运用有着很大的指导意义。但这些研究更多地考虑的是井— 地、地— 井的情况。笔者针对井— 井、倾斜井的测量方式, 对同时处理井— 地、地— 井、井— 井等情形以及如何提高计算效率等问题进行了研究。

1 电法测井三维有限差分正演问题
1.1 基本原理

点源位于地下的A点, 产生的三维电场电位u满足微分方程

其中:σ 为电导率, I为电流强度, (x, y, z)、(x0, y0, z0)分别为观测点和点电流源的坐标, δ 为狄拉克函数。边界条件为

un=0Γs, un+cosθr=0Γ;

其中:Γ s为地面边界, Γ 为地下的无穷边界, n为边界的外法向矢量, r为源点到边界上的节点的矢径, θ 为外法向矢量n和矢径r的夹角。

采用异常场法对井中电法进行数值模拟。假设电流源附近的电导率为σ m, 均匀半空间或者均匀全空间的正常场电位为up, 与供电电源无关的由电阻率分布异常引起的异常场电位为us, 它们之间的关系为

将式(2)、式(3)代入式(1), 得到异常场方程

其边界条件就是:

usn=0,  usn+cosθr=0

1.2 微分方程的离散及网格剖分

利用有限差分法对式(4)进行数值求解时, 首先对其进行体积分, 通过格林函数将体积分转化为面积分, 再利用有限差分法构成差分方程。如节点在地表上, 那么体积

ΔVi, j, k=Δxi+Δxi-1)·(Δyj+Δyj-1)8;

若节点在地下, 则有

Δ Vi, j, k= Δxi+Δxi-1)·(Δyj+Δyj-1)·(Δzk+Δzk-1)8

对式(4)进行体积分, 得到

采用格林定理对式(5)进行体积分; 对求解的区域采用六面体剖分, 然后(如), 可以得到:

其中, Si, j, k为Δ Vi, j, k的表面积。

图1 三维网格剖分及离散示意

将式(6)分解成6个面的面积分, 在网各节点上采用七点差分格式(如图1所示)代替偏导数, 可得到:

c1· us(i-1, j, k)+c2· us(i+1, j, k)+c3· us(i, j-1, k)+c4· us(i, j+1, k)+c5· us(i, j, k+1)+c0· us(i, j, k)= c1* · up(i-1, j, k)+ c2* · up(i+1, j, k)+ c3* · up(i, j-1, k)+ c4* · up(i, j+1, k)+ c5* · up(i, j, k+1)+ c0* · up(i, j, k) (14)

详细的推导过程以及连接系数的表达式见文献[10]。

对于全部的网格节点, 将系数矩阵合成一个总的系数矩阵, 会得到一个线性方程组

Aus=q,

其中:A为耦合系数矩阵, us为待求的异常场的电位, q是源项, 为已知量。

基于有限差分法的井中电法三维正演所形成的系数矩阵为大型的稀疏矩阵, 为了节省存储空间, 所以只存储非零元素。采用一维非零元素行压缩存储模式存储总体系数矩阵。因为系数矩阵中每行最多只有7个非零元素, 因此设计7个一维数组ID0(μ ), …, ID6(μ )来表示C0(μ ), …, C6(μ )中的非零元素的列号, u为行号, μ =1, 2, ……, mi× mj× mk; 节点编号先从x方向开始, 再到y方向, 最后是z方向。采用i, j, k进行节点编号, i=1, 2, …, mi; j=1, 2, …; mj, k=1, 2, …, mk; 这里mimjmk分别表示网格剖分的行号、列号、页号的最大值。节点编号由u=(j-1)× mi+i+(k-1)× mi× mj计算给出。

以上存储过程简单易行, 通过节点的编号就可以确定节点的位置, 这样便实现了CSR存储。通过求解这个方程组, 便可以得到异常场的电位。对于某个节点上的正常场电位具有解析解

u0=I/(kaπσ0r)

那么总的电位u=us, 其中:r为节点到源点的距离; I为电流; 当求解的区域为三维半空间时, ka=2, 当为三维全空间的时候, ka=4。

视电阻率由公式ρ s=k ΔuI计算, K是装置系数, u是总场电位。

2 线性方程组的求解

三维电阻率正演的速率往往取决于线性方程组的求解效率, 所以选择合适的求解方法对于提高电阻率三维正演的速率至关重要。共轭梯度法一直以来都是求解对称正定线性方程组的一种极为有效的方法, 但是它的收敛速率取决于线性方程组的条件数:条件数小则收敛快, 条件数很大时收敛速率就会变得缓慢。为了提高收敛速度, 采用预条件共轭梯度法来求解方程组。笔者采用不完全Cholesky共轭梯度迭代法[9], 即ICCG算法来对方程进行求解。

首先将系数矩阵A分解为

A=LDLT,

其中L为下三角矩阵, D为对角矩阵。然后进行迭代, 满足收敛条件就停止, 不满足则继续循环迭代, 直到满足收敛条件为止。具体计算过程见文献[9]。

3 算例

依据本文的正演算法, 笔者编制了井中电法三维有限差分正演模拟程序, 可以同时处理地— 井、井— 地、井— 井以及斜井等情况下的正演问题, 程序在Intel(R) Core (TM) I3 CPU , 2.53 GHz , 2.00 GB, Windows 7下运行。设计了几个数值算例来验证本文算法的精确性和效率性。

模型1是为了验证本算法的精确性。电阻率ρ =1 Ω · m的均匀介质半空间模型, 点电源位于地表网格中心, 单位点电源的地电场的电位解析解和数值解的相对误差如图2所示。可以看到, 在源附近相对误差是最小的; 远处的测点由于离电源较远, 电位值很小, 导致较大的相对误差, 但是都不超过3% , 并且随着距离的加大, 相对误差没有继续加大。可见, 本算法解点源三维地电场的精度是稳定可靠的。

图2 模型1均匀半无限介质中点电源电位的数值解和解析解的相对误差

模型2如图3所示, 井深400 m, 井口坐标(0, 0, 0)。井的右边有一50 m× 50 m× 50 m的异常体, 电阻率为10 Ω · m, 其中心位于(75, 0, 125)处。背景电阻率ρ 0=200 Ω · m。采用地— 井测量方式, 设供电点分别位于地面上A1=(0, 0, 0), A2=(-110, 0, 0), A3=(110, 0, 0)处。图4为所获视电阻率曲线, 可以看出视电阻率很好地反映出低阻体的存在, 在异常体的上边界视电阻率曲线均达到极小值。

图3 模型2的xOz断面示意(地— 井)

图4 模型2的视电阻率曲线

图5 模型3过xOy断面(井— 地)

模型3如图5所示。井口位于(0, 0, 0)处, ρ 0=100 Ω · m。两个异常体, 低阻异常体ρ 1=10 Ω · m, 大小为30 m× 30 m× 20 m, 位于x:[-40 m, -10 m], y:[-40 m, -10 m], z:[50 m, 70 m]处; 高阻异常体ρ 2=1 000 Ω · m, 大小为30 m× 30 m× 20 m, 在x:[10 m, 40 m], y:[10 m, 40 m], z:[30 m, 50 m]处。

为了模拟野外实际工作, 观测系统采用二极装置, 测量方式为井— 地方式。分别在井中A1(0, 0, 20)、A2(0, 0, 25)、A3(0, 0, 30)、A4(0, 0, 35)、A5(0, 0, 40)、A6(0, 0, 45)处固定供电点进行供电, 在地面上测线逐点移动测量电极进行测量, 图6给出了模拟计算的视电阻率等值线。图6a中地表ρ s等值线分布表现为一个高阻异常和一个低阻异常形态, 与异常体在地面的投影位置相对应, 因为高阻异常体更靠近地表和供电点, 所以高阻异常体引起的ρ s异常要明显高于低阻体所引起的ρ s异常。从图6b— 图6f可以看到, 随着供电点深度的加大, 低阻体所引起的ρ s异常明显高于高阻体的, 这主要是由于低阻体更加靠近供电点; 同时, 可以看到高、低阻异常体引起的ρ s异常中心有所偏移, 低阻异常体的异常中心向右上方偏移, 而高阻体的异常中心向左下方偏移。可见, 供电点的深度不同, 在地表上观测到的ρ s异常发生很大的不同。因此, 可以通过改变供电点的深度来判断异常体的分布及延伸范围。

图6 不同供电点的地面ρ s等值线

模型4如图7所示。斜井井口位于(-200, 0, 0), 倾角75° , 井深600 m。直立井井口位于(200, 0, 0), 井深600 m。背景电阻率为200 Ω · m, 异常体电阻率ρ 1=10 Ω · m, 长宽高均为100 m, 其中心位于(0, 0, 300)。

图7 模型4过xOz断面示意(井— 井)

分别在两井中的A1(-300, 0, 250)、A2(200, 0, 250)处供电, 在另一口井中测量, 结果见图8。可以看到:两井的ρ s曲线均在异常体的位置范围达到了极大值; 斜井井轴方向ρ s曲线在异常体的下边界达到了极大值, 直立井井轴方向的ρ s曲线在异常体的上边界达到了极大值。这是因为低阻异常体吸引了电流, 使得钻井所处的周围聚集大量的电流, 测量电极越靠近低阻体, 其相应的电流密度就越大, 最后形成了高阻异常。在实际应用中要注意此现象, 以免做出错误判断。

图8 模型4直立井与倾斜井的ρ s曲线

4 结论

采用有限差分法实现了井中电法三维数值模拟, 充分考虑了地— 井、井— 地、井— 井以及斜井的情形下的正演研究, 根据其原理所编制的程序同时能处理地— 井、井— 地、井— 井、斜井的正演问题。采用一维非零元素行压缩存储模式来对系数矩阵进行存储, 减少了内存需求结合不完全Cholesky(ICCG)迭代法来对有限差分法所形成的线性方程组进行求解, 显著地提高了求解的效率。所选取的模型结果验证了本文的算法的精确性和效率性, 能够很好反映异常场的特征, 为进一步的井中电法反演和解释研究工作奠定了基础。

The authors have declared that no competing interests exist.

参考文献
[1] 高杰, 柯式镇, 魏宝君, . 电法测井数值模拟现状及发展趋势分析[J]. 测井技术, 2010, 34(1): 1-5. [本文引用:1]
[2] 岳建华, 刘志新. 井—地三维电阻率成像技术[J]. 地球物理学进展, 2005, 20(2): 407-411. [本文引用:1]
[3] 徐凯军, 李桐林. 垂直有限线源三维地电场有限差分正演研究[J]. 吉林大学学报: 地球科学版, 2006, 36(1): 137-147. [本文引用:1]
[4] 王志刚, 何展翔, 刘昱. 井地直流电法三维数值模拟及异常规律研究[J]. 工程地球物理学报, 2006, 3(2): 87-92. [本文引用:1]
[5] 刘树才. 矿井直流电法三维正演计算的若干问题[J]. 物探与化探, 2004, 28(2): 170-176. [本文引用:1]
[6] 王磊. 井地电阻率法三维正反演[D]. 北京: 中国地质大学, 2008. [本文引用:1]
[7] 屈有恒. 井地电阻率法及双频激电三维数值模拟与反演研究[D]. 北京: 中国地质大学, 2008. [本文引用:1]
[8] 吴小平, 汪彤彤. 利用共轭梯度算法的电阻率三维有限元正演[J]. 地球物理学报, 2003, 46(3): 428. [本文引用:1]
[9] 吴小平. 利用不完全的Cholesky共轭梯度法求解点源三维地电场[J]. 地球物理学报, 1998, 41(6): 848-855. [本文引用:1]
[10] Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped three-dimensional structures[J]. Geophysics, 1979, 44(4): 753-780. [本文引用:1]