基于COMSOL的井地电阻率正演研究

冯硕 ，刘得军 ，张颖颖 ，秦民君 ，兀凤娜

（1.中国石油大学（北京）地球物理与信息工程学院，北京 102249；2.中国石油集团测井有限公司生产测井中心，陕西 西安 710201；3.中国石油长庆油田分公司第一采油厂，陕西 延安 716000）

0 引言

目前我国大部分油田已进入高含水开发期，随着开发的不断深入，为实现开发后期油田的稳产挖潜，对剩余油的开发与研究提出了更高的要求［1-15］。在诸多地球物理勘探方法中，井地电阻率测量方法具有对地层影响小、测量成本低等优势，在剩余油分布研究中得到了广泛的应用。

井地电阻率测量方法，是利用钻孔套管、井液作为线电极供电，通过井口周围安置的阵列接收电极测量电位差，然后进行数据处理，同时综合已知地质数据分析结果，得到地下局部电阻率的勘探方法［16］。其主要应用于圈定油气藏边界、油田开发注水检测以及压裂裂缝成像等方面［17］。为了实现井地电阻率方法的精确定量解释，必然需要进行快速高效的正演算法研究。

井地电阻率法的正演仿真，主要是通过已知空间的电阻率分布，利用数学物理方法计算得到地表的电位分布。目前解决井地电阻率测量的正演方法主要有α中心法、有限差分法、有限元法、积分方程法等。α中心法［3］简单快捷但精度不高，现在已经很少使用。有限差分法［2］的数学表达式简单直观，发展较早且应用较广，但边界刻画能力弱，只使用矩形网格，无法对复杂地质构造进行准确地模拟。有限元法［4］将求解区域划分为离散网格，将每个网格看作独立的单元进行求解，这样可以得到若干个有限单元的解，这些解的集合构成整体函数的解。积分方程法［5］对于复杂地电模型中差异边界面无法进行有效离散化，尤其是模型中存在多个非闭合异常体时，因此该方法现在只是处于试验阶段。

水永秀树等［4］利用有限差分法实现了垂直线元三维井地电位正演，随后Tetsuo Aono等［7］讨论了垂直线电流三维成像理论。徐凯军、刘地渊等［8-9］对垂直线元、任意线元三维地电模型进行了正演研究。屈有恒等［10］实现了倾斜线元的三维电场正演。刘海飞等［11］实现了连续介质三维线元井地电位的正演计算。何展翔、王志刚等［12］研究了利用积分方程法对井地电阻法的三维正演模拟。

本文首先利用COMSOL软件进行建模，然后利用COMSOL可灵活处理各类边界条件特点，施加混合边界条件，最后用有限元方法，计算高、低阻异常体，分析模型变化对地表电位分布的影响。

1 理论和方法

1.1 线电源场基本分析方法

在井地电阻率法中，电流主要通过钻孔套管和井液传到地下，在套管直径远小于套管长度的情况下，可以将套管简化为有限长线电流源来处理［13］。

井地电阻率分布规律满足以下拉普拉斯形式的微分方程：

式中：σ为线电源所处地层模型电阻率，Ω·m；Φ为模型中各节点电位，V；f为电流源，A。

线电源的电流强度为I时，f为

式中：δ为狄拉克函数，在源点处为无穷大，其余各点为0；r，rs分别为观测点与线电源的位置半径，m。

对于均匀半空间介质，线电流源产生的电位可以由点电源激发的电位积分得到：

式中：l2，l1分别为线电流源的顶、底埋深，m；zp为测点埋深，m；rp为测点到线电源的径向距离，m。

地层模型中，电阻率为非均匀的情况下，利用解析法无法求出模型中的电位分布，这时需要使用数值算法进行求解。利用有限元进行数值计算时，需要确定边界条件进而求解。

在计算区域Ω中，井地电阻率法的边界条件主要包括2种：

1）地面边界Γs，由于认为空气是绝缘体，选择使用诺依曼边界条件 σ·∂Ф/∂n=0。

2）无限远处边界面Γ∞，为保持电位的物理特性，消除电性突变，选择使用混合边界条件∂Ф/∂n-λФ=0，Ф为均匀介质情况下利用解析公式得到的电位值。

井地电阻率模型中电位满足：

1.2 COMSOL正演准确性验证

COMSOL是一款以MATLAB为内核的仿真软件［14］，相比大多数仿真软件，利用COMSOL进行电法仿真的主要优势在于能够自由定义各类边界条件，得到更为准确的正演结果；正确利用COMSOL中的physics选项可以有效简化对边界条件的处理。

地表边界面Γs，可通过选择COMSOL中的Electric insulation边界进行设定。对无限远处边界面Γ∞，需要选择 Distributed resistance 边界条件 n·J=σ（V-Vref）/d,并将混合边界条件变形为n·J=-λσФ，然后把Vref选项赋值0，d选项赋值-1/λ，即完成混合边界条件的处理。

采用验证模型对该方法的精确性进行验证。验证模型空间为1 000 m×1 000 m×1 000 m的正方体，电阻率为10 Ω·m，线电流长度为100 m，供电电流为1 A，线电源通过地表供电，供电点坐标为（0，0，0）。分别对无限远边界施加第1类、第2类及混合边界条件，经过正演计算，得到 5，10，15，20，30，40 m 处的电位值（见表1）。由表1可知，在观测范围区域内，最大绝对差值为0.001 1 V，最大相对误差小于0.25%，这样的计算精度是令人满意的［15］。

表1 不同边界条件电位值计算结果

2 算例分析

正演计算井地电阻率模型为一个正方体，为1 000 m×1 000 m×1 000 m，背景电阻率为 10 Ω·m。 异常模型的高阻体电阻率为100 Ω·m，大小为200 m×200 m×100 m，中心点位置在（110，110，-80）。 线电流长度为200 m，供电电流为5 A，线电源通过地表供电，供电点坐标为（0，0，0）。

通过COMSOL计算，得到模型中每个有限元节点位置、地表面上的电位值（见图1）。

图1 高阻异常体模型电位分布

由图1可以看出，在含有高阻异常体的情况下（异常体位置如虚线所示），0.062 63 V等势线出现明显等位圆压缩现象，并在压缩曲线上有一个小突起，其方向指向高阻异常体方向。电位分布情况与井地模型中高阻物质对电流排斥作用相符合。

对模型中异常体的电阻率进行修改，改为1 Ω·m，电阻率与背景电阻率的比例保持不变，电位分布情况见图2。

图2 低阻异常体模型电位分布

由图2可以看出，在有低阻异常体的情况下，0.087 25 V等势线明显向低阻体位置拉伸。电位分布情况与井地模型中低阻物质对电流吸引作用相符。

模型中2个异常体，电阻率分别为1，100 Ω·m，大小均为200 m×200 m×100 m，其中心点位置分别在（-110，-110，-80）、（110，110，-80），电位分布见图 3。

图3 异常体水平排列时模型电位分布

由图3可以看出，当高、低阻异常体位于同一水平面时，0.074 31 V等势线在高阻异常电阻区域出现压缩现象，并依然保留一个沿高阻异常体对角线方向的小突起，在低阻异常电阻区域明显被拉大。仿真结果表明，井地电阻率法能够大致分辨出同一高度面上，不同电阻率异常体的位置。

模型中设置2个不同高度面异常体，电阻率分别为 1，100 Ω·m，大小均为 200 m×200 m×100 m。将 2个异常体上下排列，异常体中心点位置分别放置在（-110，-110，-180）、（-110，-110，-80），电位分布见图 4。

由图4可以看出，模型中同时存在高阻异常体与低阻异常体，且当低阻体置于高阻体上方时，0.06300 V等势线表现出低阻异常体拉伸现象；与此同时，在异常体存在方向上，等势线也存在一个小突起，说明电位分布图依然保留有高阻体存在时的特征。仿真结果表明，当高阻异常体在下，低阻异常体在上时，根据地表电位能够大致判断地下电阻分布情况。

图4 异常体上下排列时模型电位分布

将异常体调换位置，当高阻异常体置于低阻异常体上方时，电位分布见图5。

图5 高阻异常体置于低阻异常体上方模型电位分布

由图5可知，高、低电阻体安置位置的互换，地表电位分布已经几乎无法判断是否有异常体的存在，等势线只存在微弱的不对称现象，无法分析地下异常体电阻率分布情况。仿真结果表明，对井地电阻率法的研究应当考虑上、下不同电阻率分布的地层的影响，在进行反演研究时也需要参照其他前期地质资料，对主要反演区域的上、下地层进行考虑。

3 结论

1）仿真结果表明，井地电阻率法对低阻物质较敏感，并不适用于石油领域勘探初期。注水油田进入高含水开发期，地层将出现明显的电阻率差异，井地电阻率法可以进行剩余油气分布测量、注水前沿检测等项目。

2）实验表明，异常体纵向分布情况不同，就会得到不同的地表电位分布，在进行反演研究时也要参照前期地质资料。

3）用具体仿真实例进行分析，验证了在使用混合边界条件的情况下，运用COMSOL软件能够得到更加精确的正演结果。

［1］付德奎，冯振雨，曲金明，等.剩余油分布现状及展望［J］.断块油气田，2007，14（2）：39-41.

［2］杨承林，于春娥，石淑灿.应用数值模拟法研究复杂断块剩余油分布［J］.断块油气田，2007，14（2）：36-38.

［3］Petrick Wm R，Sill Wm R，Wards S H.Three-dimensional resistivity inversionusingalphacenters［J］.Geophysics，1981，46（8）：1148-1162.

［4］水永秀树，青野哲雄，牛岛惠辅.流电电位法（充电法）三维模拟［J］.金旭，译.物理探查，1991，44（4）：215-226.

［5］Yutaka Sasaki.3-D resistivity inversion using the finite-element method［J］.Geophysics，1995，59（1）：1839-1848.

［6］底青云，王妙月.积分法三维电阻成像［J］.地球物理学报，2001，44（6）：843-851.

［7］Aono T，Mizunaga H，Ushijima K.Direct imaging fractures by a 4-D geoelectrical technique at Sumikawa Geothermal Field［A］.Japan：Proc.of 27th Sanford Geothermal Workshop，2002.

［8］徐凯军，李桐林.垂直有限线源三维地电场有限差分正演研究［J］.吉林大学学报：地球科学版，2006，36（1）：137-141.

［9］刘地渊，徐凯军，赵广茂，等.任意形状线电流源三维地电场研究［J］.地球物理学进展，2006，21（2）：395-399.

［10］屈有恒，张贵宾，晋风明.倾斜线源的三维电场数值模拟研究［J］.物探化探计算技术，2007，29（5）：431-435.

［11］刘海飞，陈德鹏，戴前伟，等.连续电性介质线源井-地电位三维有限元数值模拟［J］.桂林理工大学学报，2011，31（1）：28-37.

［12］王志刚，何展翔，刘昱.井地直流电法三维数值模拟及异常体规则研究［J］.工程地球物理学报，2006，3（2）：87-92.

［13］Kaufman A A，Wightman W E.A transmission-line model for electrical logging through casing［J］.Geophysics，1993，58（12）：1739-1747.

［14］王小龙，冯宏，田华光，等.基于COMSOL MULRIPHYSICS的直流电法正演模拟［J］.煤田地质与勘探，2011，39（5）：76-80.

［15］汤井田，肖晓，杜华坤，等.ANSYS在直流电法正演中的应用［J］.地球物理学进展，2006，21（3）：987-992.

［16］高军强，李书忠，李油建.激电法找油的原理及应用［J］.断块油气田，2003，10（1）：29-30.

［17］陈亮，牛艳平，张一伟.高含水期剩余油分布的多学科研究［J］.断块油气田，1999，6（5）：29-33.