水平井随钻电磁波测井实时正反演方法

王磊，刘英明，王才志，范宜仁，巫振观

（1. 中国石油大学（华东）地球科学与技术学院，山东青岛 266580；2. 海洋国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266071；3. 中国石油大学CNPC测井重点实验室，山东青岛 266580；4. 中国石油勘探开发研究院，北京 100083）

0 引言

随钻电磁波测井是水平井地质导向的常用工具[1-2]，提供的相位差和幅度比电阻率曲线常出现分离和“犄角”等复杂响应，导致其定量乃至定性解释困难。同时，随钻时提供的测量曲线数量有限且缺乏方位识别能力，也进一步增加了地层界面预测的难度[3]。

目前，水平井随钻电磁波测井资料的处理主要基于一维（1D）水平层状地层模型[4-5]和梯度算法来进行的，国外学者已经实现了对水平井资料的井斜校正和各向异性地层电阻率反演[5-6]。然而，该资料处理方法需依赖其他测井资料提供的界面位置信息，故仅适用于钻后评价。同时，梯度算法初值依赖性强，导致解释结果可能不唯一。为实时确定井周界面位置，水平井导向或储集层评价中普遍采用交互式解释方法，即不断人工调整界面位置并进行正演，以获得模拟数据与实测数据一致性高的地层模型[7]。通过引入导眼井/邻井测井资料和已知地质信息等约束，交互式解释方法有效地降低了水平井测井解释的不确定性。但是该方法是基于试错法的手动寻优，理论上其工作量大且处理精度有限，如何提升现有解释方法的效率和地层界面的预测精度是目前随钻测井亟待解决的关键问题。

快速 1D正演是水平井随钻电磁波测井资料处理的核心，其难点之一在于核函数的高效积分[8-12]。水平井测井资料处理过程中，积分核函数振荡性强、收敛速度慢，直接积分方法不再适用。目前业界普遍采用快速Hankel变换（FHT）方法正演随钻电磁波测井响应[13]。然而，该方法仍存在诸多局限：发射与接收纵向位置相同时，存在积分奇点[14]；积分精度难以控制；不同接收点的场需进行重复推导。因此，探索新的积分方法，避免阵列线圈响应的重复计算，是提高随钻电磁波测井正演速度和精度的必然要求。

本文首先研究了水平井随钻电磁波测井正演模型构建方法；然后，提出层界面最优匹配技术和索末菲高效直接积分新方法，实现了随钻电磁波测井正演加速；然后，将交互式模型调整策略和梯度寻优算法相结合，形成基于交互式反演的水平井随钻电磁波测井地层界面预测方法；最后，将正反演方法应用于现场资料处理，取得了良好的应用效果。

1 随钻电磁波测井正演模型构建

随钻电磁波测井采用单发双收的轴向天线结构，可同时提供幅度比和相位差视电阻率曲线。该类测井仪器的源距通常小于1.5 m，能够探测数米内的地层。一般而言，井周附近各地层近似认为相互平行[15-17]。此时，正演模型中应考虑井眼、仪器、地层和仪器的相对夹角θ等，即图1a所示的3D正演模型。正演模拟时，该模型可简化为一维。这主要基于两点考虑：①随钻电磁波测井响应通常受井眼影响较小[4]；②井眼影响可经相关校正进行消除。

为验证上述正演模型降维方法的可行性，建立两层地层模型，其中仪器以80°的相对夹角穿过地层。当模型中有/无井眼存在时，基于3D/1D算法计算的视电阻率曲线如图2a、2b所示。模拟时，井眼半径为0.101 6 m，井内分别填充油基钻井液（电阻率为1 000 Ω·m）、淡水钻井液（电阻率为1.0 Ω·m）和盐水钻井液（电阻率为0.1 Ω·m）。由图可知，油基钻井液填充的钻井的随钻电磁波测井响应和无井眼1D测井响应一致，此时降维方法可行。相比而言，使用水基钻井液的钻井的视电阻率值与1D响应差异明显，且钻井液电阻率越低差异越大。此时，建立考虑对应井眼情况的刻度图版，对随钻电磁波测井响应进行井眼校正，校正后的结果见图 2c、2d。可以看出，井眼校正后的测井响应与1D模拟结果基本重合，这表明模型降维正演方法同样适用于水基钻井液填充的钻井。

图2 正演模拟校正前后幅度比电阻率与校正前后相位差电阻率

2 随钻电磁波测井快速正演方法

2.1 1D层状各向异性介质偶极子源电磁场递推算法

在图1b所示的1D介质中，源激发的电磁场可通过递推算法获得。对随钻电磁波测井，接收点的磁场H可表示为水平磁偶极子（HMD）和垂向磁偶极子（VMD）激发的磁场的叠加，即：

在 1D介质中，VMD只激发横电（TE）波，而HMD既产生TE波又激发横磁（TM）波。同时，磁偶极子源产生的两种波相互分离。假定轴向发射线圈处为地层坐标系原点，当发射和接收线圈分别位于m和n层（n＞m）时，（1）式可展开为：

为简化起见，本文仅给出 TM波的传播系数表达式，TE波传播系数的推导过程可见文献[9]。对水平磁偶极子激发的TM波，FTM,h表达式如下：

边界处，反射系数和透射系数的公式可写为：

2.2 积分核函数特性

为计算接收线圈的磁场，还需实现（2）式所示的索末菲积分。仪器相对倾角分别为0°，30°，60°，80°和 90°时，图 3给出了两层模型中核函数随积分变量kρ的变化关系。其中，仪器源距和频率分别为0.5 m和2 MHz，且发射线圈固定于地层界面处。可以看出，随着倾角的增加，积分核函数的振荡性加剧、收敛性降低。当仪器与地层相互平行时，核函数不收敛，导致总场的计算困难。

总场的振荡性是由散射场和一次场共同引起，且后者具有解析表达式。图3b给出了谱域散射场核函数（hs）的收敛性。与总场相比，散射场的振荡性变弱，且在水平井时仍收敛。利用相同的地层模型，图 4进一步给出了水平井散射场积分谱与仪器距地层界面距离（D）的关系图。随D的不断增大，核函数急剧衰减、收敛性迅速增强。一般而言，相对倾角小于 60°或D大于0.1 m时，积分上限截断至100，即可保证磁场的计算误差在0.1%以内。

图3 不同井斜角度条件下核函数随积分变量的变化关系

图4 散射场的收敛性与仪器距地层界面距离的关系

2.3 散射场直接积分新方法

理论上，散射场可分解为多次反射/透射波的叠加，且多次波的贡献随反射/透射次数的增加而迅速减小。利用该性质，笔者提出一种基于一次反射/透射波扣除的积分新方法。该方法分为 3步：①谱域散射场中扣除一次反射/透射场的贡献；②被扣除项解析计算；③对剩余谱域核函数直接积分。散射场Hs来自两部分，即0阶和1阶贝塞尔函数积分。将一次反射/透射场分离，变为：

分离后的积分项变多，且一次反射/透射场的积分项无解析解。问题变为如何实现一次反射/透射波的高效计算。

当kρ趋近于无穷大时，kn,z≈kρ，且广义反射系数约等于狭义反射系数。此时，对TE和TM波的狭义反射系数进行泰勒展开，即：

基于（13）式和（14）式，可进一步将（11）式变为：

上式中，第2项为一次反射/透射波的逼近。此时，只需对第 1项进行直接积分即可。利用同样的思路，可处理一阶贝塞尔函数积分，此处不再赘述。

为验证新方法的有效性，建立图5a所示的地层模型。当仪器与地层的夹角为90°和89°时，图5b对比了新方法处理前后散射场的收敛性。由图5b可知，将一次透射/反射波高阶逼近且扣除后，剩余散射场的核函数随kρ的增加而急剧衰减。采用新方法后，积分核函数收敛性可提高3个数量级以上。此时，积分区间取[0，50]，即可满足工程精度要求。表1进一步对比了新的直接积分方法和传统FHT方法的计算效率及精度。当计算误差为万分之一时，FHT方法至少需要260个滤波点[18]，而新方法只需40个高斯勒让德积分（GLQ）点。本文给出的计算方法精度更高，且速度是传统FHT方法的6倍以上。同时，新方法的速度基本不受接收线圈数量的影响，更适用于阵列化随钻电磁波测井仪器的模拟。

图5 地层模型（a）及其新方法处理前后散射场收敛性（b）

表1 索末菲积分方法对比（θ=89°）

2.4 模型边界自适应匹配方法

模型的地层界面数是制约随钻电磁波测井正演速度的关键之一，且两者基本成线性反比关系[15]。地质导向模型通常包含众多地层，故正演中必须进行适当截取。目前多采用固定阈值的方法截断模型边界，如正演中只考虑距仪器10 m以内的地层界面。由于不同电阻率和频率条件下仪器探测范围不一，该边界截断方法可能会导致正演精度不够或计算耗时等问题。为此，本节提出了一种模型边界自适应匹配方法。

在均匀地层中，电磁场的衰减通常用趋肤深度δ表示：

1D介质中，源激发的电磁场在界面处的衰减率s是趋肤深度、层厚和透射系数的函数。假定源在m层，仪器上部l层界面处s可近似为：

当s小于1%时，电磁波基本衰减完毕，更远处地层对测井响应影响极小，故正演中取s=1%为模型截断边界。为验证上述方法的可行性，建立图6a所示各向异性俄克拉荷马模型[9]。仪器自上而下穿过地层时，采用不同边界截断方法计算的幅度比电阻率响应和相对误差见图6。可以看出，当截断距离固定为8 m时，传统方法的模拟结果与精确解一致。将截断距离减小至4 m时，传统方法的模拟精度急剧下降。特别是在横向深度25～170 m处，传统方法的相对误差常高于5%。相比而言，采用自适应边界匹配方法的计算结果相对误差可控制在0.5%以内，这表明新方法准确、可靠。

图6 不同截断方法幅度比电阻率响应与计算精度对比（仪器源距和频率分别为1.016 m和400 kHz）

为进一步说明新方法的优势，图 7对比了两种边界截断方法所需的正演模型层数。与传统方法相比，新方法能够对钻遇地层进行自适应截断，正演模型层数急剧减少。采用新方法后，随钻电磁波测井模拟速度可提升一倍以上。采用Intel（R）i7-8700 CPU计算时，本文给出的正演方法每秒可计算约16 000个测井点，能够满足水平井随钻电磁波测井实时交互反演的需要。

图7 不同边界截断方法的正演模型层数对比

3 水平井随钻电磁波测井交互式反演方法

3.1 随钻电磁波测井交互反演方法

为确定井周地层界面位置，本节给出了一种水平井随钻电磁波测井交互式反演方法，具体流程见图8。该方法的核心为邻井建模、反演模型人工调整和地层界面梯度寻优。邻井建模是指基于导眼井/邻井信息，提取地层电阻率与岩性序列，以建立参考导向/解释模型。将参考模型与邻近测井点的先验约束（已知的反演结果）相结合，可以确定初始反演模型。然后，利用梯度算法和1D快速正演程序，对地层界面寻优。若反演结果不满足精度误差，则不断手工调整地层界面位置进行梯度寻优，直至模拟与实测结果相吻合。

图8 水平井随钻电磁波测井反演流程图

3.2 边界反演算法与初值选取

随钻电磁波测井反演可转换为求实测资料Ra与模拟响应S的最小二乘问题，反演目标函数满足下式：

更新方式如下：

为获取目标函数的最小值，采用正则化 Gauss-Newton方法求取。在第t次迭代时，仪器附近的地层界面位置可用下式确定：

采用正则化梯度方法时，反演结果的精度常取决于反演初值的精度。为准确预测地层界面位置，本文采用交互式多初值反演策略，该策略实现方式如下：①利用邻近测井点界面信息作初值；②基于初始模型反演结果，结合 Ciflog软件，可视化手工调整地层并建立新的初始反演模型；③重复步骤②，直至反演结果满足精度误差。交互式调整反演初值时，应当遵循以下原则：①视电阻率值远高于地层序列电阻率值时，则将井眼距地层上或下界面的距离变小；②仪器在高电阻率地层时，若重构响应小于实测值，则适当缩小仪器距地层界面的距离d，反之则增大d。一般而言，对模型进行3～5次交互式调整和梯度迭代，即可获取准确的地层界面位置。利用1D快速正演算法，该交互式反演方法每秒可处理几十至几百个测井点，满足了随钻电磁波测井资料实时处理的需要。

4 应用实例

将交互式正反演方法应用于吉林油田某水平井，该井靶层为粉砂岩且层内夹薄泥岩层。该井水平段地层呈水平展布，地层结构沿横向变化较为缓慢，井眼与地层界面的相对夹角为60°～133°。根据邻井资料提取地层序列，并建立初始地质模型（见图9）。在横向距离214.6 m处，自然伽马曲线急剧下降，而视电阻率值迅速增大至 20 Ω·m，故判断仪器由此入靶。入靶后，通过电阻率曲线难以直接判断仪器进出层情况。此时，参考地震剖面地层走向，然后根据自然伽马曲线确定井轨迹进/出地层界面的位置。对比第1道内的两条曲线可知，初步调整后的地质模型的自然伽马值与实测结果基本吻合。然而，基于该模型正演后的随钻电磁波测井响应与实测值差距较大，无法满足精细解释的需要。

进一步利用随钻电磁波测井反演方法，结合具体仪器模式（斯伦贝谢公司的MCR仪器），以准确获取井周地层界面位置。经逐点交互式反演后，最终确定的水平井解释模型见图10。与图9相比，综合解释后的模型更为平滑，且实测相位电阻率曲线与模拟重构响应一致性高，证明了解释模型的准确性和可靠性。解释结果表明本井储集层钻遇率高，故对全井段进行压裂。该水平井放喷初期日产油86.5 m3，为纯油层。

图9 基于随钻自然伽马测井资料的水平井测井解释结果（不同颜色代表不同地层）

图10 水平井随钻电磁波测井交互式反演结果（不同颜色代表不同地层）

5 结论

针对水平井中地层界面位置预测的难题，本文给出了一种随钻电磁波测井实时正反演方法，并将其应用于实际资料处理，得到以下认识。

利用散射场一次反射/透射波的逼近与扣除方法，索末菲积分的收敛性可提高 3个数量级。相比于传统方法，新积分方法适用于任意倾角情况，其精度高、速度快，且更适用于阵列电测井仪器响应的正演。

模型边界自适应截断方法大大减少了正演模型的层数，平衡了传统方法计算速度与精度之间的矛盾。将新的积分方法和边界处理方法相结合后，本文研发的随钻电磁波测井正演算法每秒可计算16 000个测井点以上，满足了实时正反演的需求。

水平井随钻电磁波测井交互式反演方法充分利用了邻井信息和典型电磁波测井响应规律等先验信息和人为经验，能够实时、准确提取井周地层界面位置，减少了地层解释的多解性和不确定性。

符号注释：

An——第n个地层中波场的幅度，无因次；C(m)——反演目标函数；d——地层层厚，m；D——仪器距离地层界面的垂直距离，m；FTE,h，FTE,v——水平和垂直方向的源在接收点所在层激发 TE波的传播系数；FTM,h——水平方向的源在接收点所在层激发TM波的传播系数；h——磁场的积分核函数，A/m；hs——散射磁场的积分核函数，A/m；H——磁场强度，A/m；Hpq——p方向的源激发的q方向的磁场分量，A/m；Hs——磁场的散射场强度，A/m；——与 0阶贝塞尔函数相关的散射磁场强度，A/m；——与1阶贝塞尔函数相关的散射磁场强度，A/m；I——单位矩阵；J——雅克比矩阵；J0——0阶贝塞尔函数；J1——1阶贝塞尔函数；k——波数，无因次；kρ——积分变量；kn,h,z，kn,v,z——第n层地层水平和垂直方向的波数；m，n，l——地层序号；m——待反演的参数矢量；N——模型中地层总数；mref——已知参考模型；O——取高阶无穷小函数；Ra——实测视电阻率数据，Ω·m；RP33H——频率2 MHz、源距0.838 2 m（33 in）的相位差电阻率，Ω·m；，——zn处水平和垂向源激发的上行TE波的狭义反射系数；，——zn处上行和下行TM波的狭义反射系数；，——zn处上行和下行 TM波的广义反射系数；，——zn处水平和垂向源激发的上行TE波的广义反射系数；r——接收点距发射点的距离，m；s——源激发的电磁场在界面处的衰减率；S(m)——正演响应，Ω·m；t——迭代次数；——zn处上行 TM波的狭义透射系数；z——垂向坐标，m；zn——第n个地层界面的垂向坐标，m；σn,h，σn,v——第n层地层水平和垂向电导率，S/m；ρ——径向距离，m；χ——正则化系数；χ0——初始正则化系数；δ——趋肤深度，m；θ——仪器轴与地层界面法向的相对夹角，（°）；ω——角频率，rad/s；μ——自由空间的磁导率，H/m；λ——地层电阻率各向异性系数；ξ1，ζ1——和积分中上行波的幅度项；ξ2，ζ2——和积分中下行波的幅度项；，——一次反射/透射波的幅度；——谱域上行波扣除项系数；——谱域下行波扣除项系数；——求向量的2范数。