2016, Vol. 59
2. 中国科学院地质与地球物理研究所, 北京 100029
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
过套管电阻率测井是油气井在开发过程中油气藏动态监测和剩余油地层重新评价的重要测井方法之一(Askey et al.，2002; Kaufman，1990; Kaufman and Wightman，1993; 尹军强等，1998；高杰等，2008; 刘福平等，2007; Vail，1993; 谢树棋等，1999)，但自从20世纪30年代过套管电阻率测井被提出后的几十年间一直未取得突破性进展(Vail，1989，1991; Vail et al.，1993，1995; Schenkel and Morrision，1994，1990;Xie et al.，1999; 刘福平等，2013)．直到20世纪90年代初Kaufman发表了基于传输线方程的套管井电阻率测井近似理论模型和测量理论后(Benimeli et al.，2002; Kaufman，1990; Kaufman and Wightman，1993; 刘福平等，2013)，在过套管电阻率测井理论上和在测井仪器的研制上才取得了实质性进展(Aulia et al.，2001; Zhou et al.，2002；高杰等，2008; 刘福平等，2007)．然而到目前为止，尚未涉及到过套管电阻率测井的参数反演问题.由于金属套管的电阻率与地层的电阻率存在巨大的差别(有的相差1010)，致使有限元、有限差分法离散方程的矩阵系数大小存在巨大差别，用这两种方法实现套管电阻率测井的正演计算仍然存在困难(目前还没见到成功的算例)，传输线方法已成为过套管电阻率测井成熟的正演计算方法(高杰等，2008; Fondyga et al.，2004；刘福平等，2013)，因此过套管电阻率测井参数反演的正演计算利用传输线方法目前应为一种好的选择(目前过套管电阻率测井的地层参数反演及Jacobi矩阵的计算尚未见到算法报道).要实现地层参数的反演关键是如何实现Jacobi矩阵的计算，反演计算时间的长短和精度取决于Jacobi矩阵的计算方法(魏宝君和张庚骥，2002; 丁继才等，2007; 张美根等，2003; 刑光龙等，2007；刘福平和杨长春，2003; Liu et al.，2012)，因此针对层状地层开展Jacobi矩阵算法研究具有重要意义.本文研究了套管壁电势、电流对地层横向电阻导数的Jacobi矩阵计算方法，导出了过套管电阻率测井反演地层参数Jacobi矩阵的解析表示，并探讨了Jacobi矩阵的特性，为进一步实现电阻率测井资料反演和测井数据解释都将具有实际意义.
2 场分布的正演计算设金属套管单位长度的电阻为Rc，金属套管单位长度所对应的地层的横向电阻为T(漏电电阻)，则传输线方程可写为(高杰等，2008; 刘福平等，2007)
|
(1) |
其中U，I分别为套管井的电势和流过套管的电流，
|
(2) |
其中Ap*、Bp*为待定系数，ξp=Tpαp.由第p、p+1层边界电流和电势连续条件得系数关系
|
(3) |
其中
|
(4) |
考虑n+1层地层后有
|
(5) |
由于在n+1层中z可以取无限远，为保证电流I(z)有限，应取An+1=0，设电源在坐标原点，则I1(0)=I0，其中I0为由电源流出而流向上半个空间的电流.利用电流源条件和方程(5)可得Bn+1、A1、B1，再利用公式(3)可得任意一层中的解.
3 过套管电阻率测井的Jacobi矩阵沿套管将求解模型剖分成一系列小单元，设每一个小单元模型参数为常数，设参数化后模型参数向量为m=(m1，m2，…，mNz)，Nz为参数总数，第i个观测数据Ui(如电势)，则有Ui=Ui(m)，i=1，2，…，Mq，Mq为观测数据的个数，设模型初值为m0=(m10，m20，…，mNz0)，在m0点将Ui展开成Taylor级数并取一级近似有
|
(6) |
矩阵形式为
|
(7) |
其中
|
(8) |
|
(9) |
(8)式即为过套管电阻率测井观测数据对地层模型参数的导数(称Jacobi矩阵或敏感系数).
4 Jacobi矩阵的计算方法为了利用解析结果计算Jacobi矩阵，设观测数据也有Nz个(这样可以得到计算Jacobi矩阵的解析结果，利用这个解析解可以计算任意一个观测数据的Jacobi矩阵)，与网格个数相同.将方程(1)对Tj求导得第p层中电势、电流对Tj的导数方程
|
(10) |
其中
|
(11) |
在纵向多层状地层的界面电流、电势是连续的(方程同时对地层参数求导等式仍成立)，因此由于某一网格ΔTj的微小变化，在相邻两层中所引起的变化ΔU、ΔI(电荷守恒)也应该相等，所以Gpj、Dpj在地层界面是连续的(Zhang et al.，2003)，即
|
(12) |
由于套管壁上电压与电流的关系(高杰等，2008; 刘福平等，2007)为
|
(13) |
|
(14) |
其中
|
(15) |
由(13)(14)及(15)式得
|
(16) |
实际上当观测数据少于网格个数时，因为p只是取整数，当p≠j时(第p、j层不在同一地层)，方程(10)为齐次方程，其通解为
|
(17) |
由(13)式得
|
(18) |
其中
当p=j时(即第p、j层为同一地层)，观测点和方程(10)为非齐次方程，其特解
|
(19) |
非齐次方程解为
|
(20) |
同理
|
(21) |
其中Ip、Up为方程(1)的解.若纵向有Nz层地层，利用Gpj、Dpj在地层界面连续可得两相邻地层均为齐次方程时的系数递推关系
|
(22) |
其中
|
(23) |
在方程(10)中每选一个j则有Nz个微分方程，只有当p=j时，微分方程为非齐次，在该层中由上界面连续边界条件(12)及(20)、(21)式得
|
(24) |
|
(25) |
写成矩阵形式有
|
(26) |
其中
|
(27) |
在第p层中由下界面连续边界条件(12)及(20)、(21)式得
|
(28) |
其中
|
(29) |
由(22)式得
|
(30) |
|
(31) |
将(30)、(31)式代入(26)、(28)式得
|
(32) |
由电流源(16)式得
|
(33) |
在(32)、(33)式中只有A1、B1、An+1、Bn+1为未知常数，考虑到在z→∞时Gpj有限，则An+1=0，利用(32)、(33)式可计算A1、B1、Bn+1，再利用(26)(28)及(22)式可确定全部待定系数，从而实现Jacobi矩阵的解析计算.
5 过套管电阻率测井的地层电阻率反演方法设反演目标函数为
|
(34) |
其中Ud为过套管电阻率测井电势实际测量结果，U为地层模型电势响应(由地层模型计算的电势分布).设f(m)=U-Ud，则f(m)为地层模型电势响应与过套管电阻率测井电势实际测量结果所构成的参差函数，若设mk为m的第k次反演迭代近似.使反演目标函数极小得迭代公式
|
(35) |
|
(36) |
其中Gk=ΔF(m)，λk为第k次迭代计算的步长.在式(35)中有时会出现矩阵GkTGk的奇异或接近奇异的情形，为提高计算的速度和稳定性，把正定对角矩阵加到GkTGk上去，使矩阵变成条件数较好的对称正定矩阵(Marquardt方法)，反演迭代公式修改为
|
(37) |
其中I为n阶单位矩阵，αk为一个正实常数，根据计算精度要求通过试算适当选取.
6 数值算例及分析 6.1 Jacobi矩阵特性分析下面取四层状地层模型为算例，电极供电电流I0=6 A，电源位于坐标原点，套管电阻率为ρc=2.0×10-7 Ωm，套管半径a=0.1 m，套管厚度Δa=0.01 m，四层状地层界面位于dp=15，18，22 m，取观测数据点的坐标为z0=10，17，20，24 m(表 1及图 1均使用上述统一参数).表 1的3种情况的差别仅是地层的电阻率不同，表的第2，3，4，5列分别为在z0=10，17，20，24 m测量点电势对第1，2，3，4层地层横向电阻的导数，表 1中Case 1四层状地层的电阻率比较接近(ρb=1.1，1.0，1.1，1.0 Ωm)，表 1中Case 3四层状地层的电阻率差距较大(ρb=20，1，20，1 Ωm).表 1显示地层电阻率对电势敏感系数(电势对地层横向电阻的导数)影响较大.图 1为电势对地层横向电阻的导数随z坐标的变化曲线，由于金属套管的高导电性的影响(电势沿套管的变化缓慢，接近线性)，使得电势敏感系数随z坐标的变化基本呈线性关系，曲线1、2、3、4分别为测量点位置变化时(z坐标是变化的)，测量电势对第1，2，3，4层地层横向电阻的导数，图 1中曲线1，2是重合的.在该算例中第4层是一个半无界的空间层，电势敏感系数与其他层的差别较大.
|
表 1 电势对地层横向电阻的导数 Table 1 Derivations of potential with respect to transverse resistance of formation |
|
图 1 电势对地层横向电阻的导数曲线 Fig. 1 Curves of derivations of potential with respect to transverse resistance of formation |
我们还计算了Jacobi矩阵(电势敏感系数)方程的条件数，对应表 1的3种情况的Jacobi矩阵方程条件数分别是Con=2.47×104、4.73×104、6.63×105，该表说明Jacobi矩阵方程病态较重，且随层状地层电阻率差别的增大方程条件数变差，因此利用电势的过套管电阻率测井参数反演，选择什么样的反演方法十分重要，这仍是需要进一步研究和探索的新问题.为提高计算的稳定性，我们在公式(37)中把正定对角矩阵加到GkTGk上去，试图改变反演迭代矩阵的条件数.设E=GkTGk+αkI，针对表 1的3种情况的计算条件，其对应反演迭代矩阵E的条件数分别为
|
其中关于反演公式(37)中αk，经我们反复试验发现取GkTGk每一行代数和的模作为αk，取得了较好的效果(该算例αk是这样选取的)，这样αk在反演过程中还实现了动态变化.
6.2 数值反演算例为考察Jacobi矩阵的可靠性及反演方法的可行性，下面以三层状和五层状地层实现了过套管地层电阻率反演，除地层参数外，其余计算条件与前面图 1算例相同.图 2为三层状地层模型，其初值为15 Ωm，反演的最大相对误差为1.3%，反演结果与地层模型取得很好的逼近.图 3为五层状地层模型，其初值分为2段，第1，2层为8 Ωm，第3，4，5层为13 Ωm，反演结果见表 2，其最大相对误差为2.2%，反演结果与地层模型也取得很好的逼近，验证了方法的可靠性与可行性.算例说明该方法是可以用于过套管电阻率的反演.
|
图 2 地层电阻率反演结果(三层状地层) Fig. 2 Inversion result of formation resistance (3-layer model) |
|
图 3 地层电阻率反演结果(五层状地层) Fig. 3 Inversion result of formation resistance (5-layer model) |
|
|
表 2 层间电阻率反演结果 (五层状地层模型,反演迭代次数为k=8) Table 2 Inversion results of interlayer resistance (5-layer formation model, with iteration number k=8) |
本文针对层状地层，导出了过套管电阻率测井套管壁电势、电流对地层横向电阻导数(Jacobi系数矩阵)的解析表示，实现了地层参数反演中Jacobi系数矩阵的计算，选择Marquardt反演方法实现了过套管测井的地层电阻率反演.由于Jacobi系数矩阵计算的解析表示，可获得较快的计算速度.本文通过算例对Jacobi系数矩阵的特性进行了探讨，结果表明由于金属套管的高导电性影响使得电势敏感系数随z坐标的变化接近线性关系，且变化缓慢，地层电阻率的差异对电势敏感系数有较大影响；本算法是基于传输线方程实现的，同样克服了由于金属套管的电阻率与地层电阻率存在巨大差别给有限元、有限差分法模拟过金属套管问题所造成的困难.算例还考察了Jacobi矩阵特性，发现由Jacobi矩阵组成的反演方程属病态方程，针对反演迭代矩阵通过试验发现，利用每次反演迭代矩阵每行代数和的模作为Marquardt反演方法中每行的迭代校正因子，较好地改善了反演迭代系数矩阵条件数.本文提出了地层参数反演中Jacobi系数矩阵的计算方法及地层电阻率反演算法，并用解析的形式实现了Jacobi系数矩阵的计算，不仅计算精度高而且计算速度快(推导过程没有新的近似)，反演算例表明，反演结果能较好地逼近地层模型，为进一步的过套管电阻率测井数据处理提供了理论依据和方法.