水平井随钻电磁波测井实时正反演方法
王磊1,2,3, 刘英明4, 王才志4, 范宜仁1,2,3, 巫振观1,2,3
1.中国石油大学(华东)地球科学与技术学院,山东青岛 266580
2.海洋国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛 266071
3.中国石油大学CNPC测井重点实验室,山东青岛 266580
4.中国石油勘探开发研究院,北京 100083

第一作者简介:王磊(1989-),男,山东潍坊人,现为中国石油大学(华东)在站博士后,主要从事电法测井正反演研究及随钻地质导向研究。地址:山东省青岛市黄岛区长江西路66号,中国石油大学(华东)地球科学与技术学院,邮政编码:266580。E-mail:upcwanglei199133@163.com

摘要

基于层状地层随钻电磁波测井伪解析公式,提出边界最优匹配技术,来实现对钻遇地层结构的自适应截断;同时提出散射场一次反射/透射波高阶渐近与扣除方法,极大地提高了索末菲积分的收敛速度,克服了传统数字滤波积分方法精度低、适用性弱等问题,实现了对随钻电磁波测井正演的加速。利用邻井或导眼井等先验信息,通过对地层模型进行交互式调整,为地层界面的预测提供最优反演初值。结合梯度寻优算法形成水平井随钻电磁波测井实时交互式反演方法。实际资料处理结果表明,随钻电磁波测井实时交互式反演方法很好地解决了井周地层边界位置确定的难题,为水平井钻井轨迹最优化和油藏解释提供了依据。图10表1参18

关键词: 随钻电磁波测井; 水平井; 实时正演; 交互式反演; 地层界面
中图分类号:P631.8 文献标志码:A
Real-time forward modeling and inversion of logging-while-drilling electromagnetic measurements in horizontal wells
WANG Lei1,2,3, LIU Yingming4, WANG Caizhi4, FAN Yiren1,2,3, WU Zhenguan1,2,3
1. School of Geosciences, China University of Petroleum, Qingdao 266580, China
2. Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China
3. CNPC Key Laboratory for Well Logging, China University of Petroleum, Qingdao 266580, China
4. Research Institute of Petroleum Exploration & Development, PetroChina, Beijing 100083, China
Abstract

Based on the pseudo-analytical equation of electromagnetic log for layered formation, an optimal boundary match method is proposed to adaptively truncate the encountered formation structures. An efficient integral method is put forward to significantly accelerate the convergence of Sommerfeld integral. By asymptotically approximating and subtracting the first reflection/transmission waves from the scattered field, the new Sommerfeld integral method has addressed difficulties encountered by the traditional digital filtering method, such as low computational precision and limited operating range, and realized the acceleration of the computation speed of logging-while-drilling electromagnetic measurements (LWD EM). By making use of the priori information from the offset/pilot wells and interactively adjusting the formation model, the optimum initial guesses of the inversion model is determined in order to predict the nearby formation boundaries. The gradient optimization algorithm is developed and an interactive inversion system for the LWD EM data from the horizontal wells is established. The inverted results of field data demonstrated that the real-time interactive inversion method is capable of providing the accurate boundaries of layers around the wellbore from the LWD EM, and it will benefit the wellbore trajectory optimization and reservoir interpretation.

Keyword: logging-while-drilling electromagnetic measurement; horizontal well; real-time forward modeling; interactive inversion; bed boundary
0 引言

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

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

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

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

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

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

图1 水平井随钻电磁波测井正演模型

为验证上述正演模型降维方法的可行性, 建立两层地层模型, 其中仪器以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)激发的磁场的叠加, 即:

$H\text{=}\int_{0}^{\infty }{h\text{d}{{k}_{\rho }}}\text{=}{{H}_{xx}}\text{si}{{\text{n}}^{2}}\theta +{{H}_{zz}}\text{co}{{\text{s}}^{2}}\theta +\frac{\text{sin}2\theta }{2}\left( {{H}_{xz}}+{{H}_{zx}} \right)$ (1)

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

$H\text{=}\frac{1}{4\text{ }\!\!\pi\!\!\text{ }}\int_{0}^{\infty }{\left\{ \left[ \frac{\text{sin}2\theta }{2}\left( k_{\rho }^{2}{{F}^{\text{TE}, \text{h}}}-\frac{ik_{\rho }^{2}}{{{k}_{n, \text{h, }z}}}\frac{\partial {{F}^{\text{TE, v}}}}{\partial z} \right)-\frac{\text{si}{{\text{n}}^{2}}\theta }{\rho }\left( \frac{\partial {{F}^{\text{TE}, \text{h}}}}{\partial z}+\frac{\omega \mu {{\sigma }_{n\text{, v}}}{{F}^{\text{TM}, \text{h}}}}{{{k}_{n, \text{v}, \ z}}} \right) \right]{{\text{J}}_{1}}({{k}_{\rho }}\rho )+\left( {{k}_{\rho }}{{F}^{\text{TE}, \text{h}}}\text{si}{{\text{n}}^{2}}\theta +\frac{ik_{\rho }^{3}}{{{k}_{n\text{, h, }z}}}{{F}^{\text{TE}, \text{v}}}\text{co}{{\text{s}}^{2}}\theta \right){{\text{J}}_{0}}({{k}_{\rho }}\rho ) \right\}\text{d}{{k}_{\rho }}}$(2)

其中 ${{k}_{n, \text{h}, z}}\text{=}\sqrt{i\omega \mu {{\sigma }_{n, \text{h}}}-k_{\rho }^{2}}$ $k_{n, v, z}=\sqrt{{i \omega \mu \sigma}_{n, v} -k^{2}_{p}}$...

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

${{F}^{\text{TM}, \text{h}}}=\left\{ \begin{array}{* {35}{l}} {{A}_{m}}\left( {{\text{e}}^{i{{\lambda }_{m}}{{k}_{m, \text{v}, z}}z}}+\tilde{R}_{m, m+1}^{\text{TM}, \text{h}}{{\text{e}}^{i{{\lambda }_{m}}{{k}_{m, \text{v}, z}}\left( 2{{z}_{m}}-z \right)}} \right) & \left( m=n \right) \\ {{A}_{n}}\left( {{\text{e}}^{i{{\lambda }_{n}}{{k}_{n, \text{v, }z}}\left( z-{{z}_{n-1}} \right)}}+\tilde{R}_{n, n+1}^{\text{TE}, \text{h}}{{\text{e}}^{i{{\lambda }_{n}}{{k}_{n, \text{v}, z}}\left( 2{{z}_{n}}-z-{{z}_{n-1}} \right)}} \right) & \left( m\ne n \right) \\ \end{array} \right.$ (3)

${{A}_{n}}=\left\{ \begin{matrix} \frac{1+\tilde{R}_{m, m-1}^{\text{TM}, \text{h}}{{\text{e}}^{-2i{{\lambda }_{m}}{{k}_{m, \text{v}, z}}{{z}_{m-1}}}}}{1-\tilde{R}_{m, m-1}^{\text{TM}, \text{h}}\tilde{R}_{m, m+1}^{\text{TM}, \text{h}}{{\text{e}}^{2i{{\lambda }_{m}}{{k}_{m\text{, v, }z}}\left( {{z}_{m}}-{{z}_{m-1}} \right)}}}\ \ \ \left( m=n \right)\text{ } \\ \frac{T_{n, n+1}^{\text{TM}, \text{h}}{{\text{e}}^{i{{\lambda }_{n}}{{k}_{n\text{, v, }z}}{{z}_{n}}}}A_{n-1}^{\text{TM}, \text{h}}}{1-\tilde{R}_{n+1, n+2}^{\text{TE}, \text{h}}{{\text{e}}^{2i{{k}_{n+1, \text{h, }z}}\left( {{z}_{n+1}}-{{z}_{n}} \right)}}R_{n+1, n}^{\text{TM}, \text{h}}}\ \ \ \left( m\ne n \right) \\ \end{matrix} \right.$ (4)

边界处, 反射系数和透射系数的公式可写为:

$\left\{ \begin{align} & \tilde{R}_{n, n-1}^{\text{TM}, \text{h}}=\frac{R_{n, n-1}^{\text{TM, h}}+\tilde{R}_{n-1, n-2}^{\text{TM}, \text{h}}{{\text{e}}^{2i{{k}_{n-1, \text{h}, z}}\left( {{z}_{n-1}}-{{z}_{n-2}} \right)}}}{1+R_{n, n-1}^{\text{TM, h}}\tilde{R}_{n-1, n-2}^{\text{TM, h}}{{\text{e}}^{2i{{k}_{n-1, \text{h}, z}}\left( {{z}_{n-1}}-{{z}_{n-2}} \right)}}} \\ & \tilde{R}_{n, n\text{+1}}^{\text{TM}, \text{h}}=\frac{R_{n, n\text{+1}}^{\text{TM}, \text{h}}+\tilde{R}_{n\text{+1}, n\text{+2}}^{\text{TM}, \text{h}}{{\text{e}}^{2i{{k}_{n+1, \text{h}, z}}\left( {{z}_{n+1}}-{{z}_{n}} \right)}}}{1+R_{n, n\text{+1}}^{\text{TM}, \text{h}}\tilde{R}_{n\text{+1}, n\text{+2}}^{\text{TM}, \text{h}}{{\text{e}}^{2i{{k}_{n+1, \text{h}, z}}\left( {{z}_{n+1}}-{{z}_{n}} \right)}}} \\ \end{align} \right.$(5)

$\left\{ \begin{align} & T_{n, n+1}^{\text{TM}, \text{h}}=\frac{2{{\lambda }_{n}}{{\sigma }_{n, \text{v}}}{{k}_{n+1, \text{v}, z}}}{{{\lambda }_{n+1}}{{\sigma }_{n+1, \text{v}}}{{k}_{n, \text{v}, z}}+{{\lambda }_{n}}{{\sigma }_{n, \text{v}}}{{k}_{n+1, \text{v}, z}}} \\ & R_{n, n+1}^{\text{TM}, \text{h}}=1-T_{n, n+1}^{\text{TM}, \text{h}} \\ & R_{n+1, n}^{\text{TM}, \text{h}}=-R_{n, n+1}^{\text{TM}, \text{h}} \\ \end{align} \right.$(6)

将式中的TM, h变为TE, h或TE, v, 同时将λ nkn, v, z替换为kn, h, z, 即可获得TE波的广义反射系数$\tilde{R}_{n, n+1}^{\text{TE}, \text{h}}$和 $\widetilde{R}^{TE, v}_{n, n+1}$.。狭义反射系数$R_{n, n+1}^{\text{TE}, \text{h}}$和 $R^{TE, v}_{n, n+1}$..表达式如下:

$R_{n, n+1}^{\text{TE}, \text{h}}\text{=}R_{n, n+1}^{\text{TE}, \text{v}}\text{=}\frac{{{k}_{n, \text{h}, z}}-{{k}_{n+1, \text{h}, z}}}{{{k}_{n, \text{h}, z}}-{{k}_{n+1, \text{h}, z}}}$ (7)

2.2 积分核函数特性

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

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

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

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

2.3 散射场直接积分新方法

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

${{H}^{\text{s}}}\text{=}H_{0}^{\text{s}}+H_{1}^{\text{s}}$ (8)

$H_{0}^{\text{s}}=\int_{0}^{\infty }{\left[ {{\xi }_{1}}\left( {{k}_{n, z}} \right){{\text{e}}^{i{{k}_{n, z}}z}}+{{\xi }_{2}}\left( {{k}_{n, z}} \right){{\text{e}}^{-i{{k}_{n, z}}z}} \right]{{\text{J}}_{0}}\left( {{k}_{\rho }}\rho \right)\text{d}{{k}_{\rho }}}$ (9)

$H_{1}^{\text{s}}\text{=}\int_{0}^{\infty }{\left[ {{\zeta }_{1}}\left( {{k}_{n, z}} \right){{\text{e}}^{i{{k}_{n, z}}z}}+{{\zeta }_{2}}\left( {{k}_{n, z}} \right){{\text{e}}^{-i{{k}_{n, z}}z}} \right]{{\text{J}}_{1}}\left( {{k}_{\rho }}\rho \right)\text{d}{{k}_{\rho }}}$ (10)

将一次反射/透射场分离, $H_{0}^{s}$变为:

$H_{0}^{\text{s}}=\int_{0}^{\infty }{\left[ {{\xi }_{1}}\left( {{k}_{n, z}} \right){{\text{e}}^{i{{k}_{n, z}}z}}+{{\xi }_{2}}\left( {{k}_{n, z}} \right){{\text{e}}^{-i{{k}_{n, z}}z}}-\xi _{1}^{1}\left( {{k}_{n, z}} \right){{\text{e}}^{i{{k}_{n, z}}{{z}_{m}}}}-\xi _{2}^{1}\left( {{k}_{n, z}} \right){{\text{e}}^{-i{{k}_{n, z}}{{z}_{m+1}}}} \right]{{\text{J}}_{0}}\left( {{k}_{\rho }}\rho \right)\text{d}{{k}_{\rho }}+}\int_{0}^{\infty }{\left[ \xi _{1}^{1}\left( {{k}_{n, z}} \right){{\text{e}}^{i{{k}_{n, z}}{{z}_{m}}}}+\xi _{2}^{1}\left( {{k}_{n, z}} \right){{\text{e}}^{-i{{k}_{n, z}}{{z}_{m+1}}}} \right]{{\text{J}}_{0}}\left( {{k}_{\rho }}\rho \right)\text{d}{{k}_{\rho }}}$ (11)

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

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

$\left\{ \begin{matrix} {{\left. R_{n, n+1}^{\text{TE, h}} \right|}_{{{k}_{\rho }}\to \infty }}\text{=}\frac{i\omega \mu }{k_{\rho }^{2}}\frac{{{\sigma }_{n+1, \text{h}}}-{{\sigma }_{n, \text{h}}}}{4}+O\left( k_{\rho }^{-4} \right) \\ {{\left. R_{n, n+1}^{\text{TM, h}} \right|}_{{{k}_{\rho }}\to \infty }}\text{=}\frac{{{\lambda }_{n+1}}{{\sigma }_{n+1, \text{v}}}-{{\lambda }_{n}}{{\sigma }_{n, \text{v}}}}{{{\lambda }_{n+1}}{{\sigma }_{n+1, \text{v}}}+{{\lambda }_{n}}{{\sigma }_{n, \text{v}}}}+O\left( k_{\rho }^{-2} \right) \\ \end{matrix} \right.$ (12)

忽略上式的高阶项, $\xi _{1}^{1}$和$\xi _{2}^{1}$通常可构建为:

$\left\{ \begin{align} & {{\left. \xi _{1}^{1}\left( {{k}_{n, z}} \right) \right|}_{{{k}_{\rho }}\to \infty }}\text{=}\frac{{{k}_{\rho }}}{{{k}_{n, z}}}{{{\tilde{\xi }}}_{1}} \\ & {{\left. \xi _{2}^{1}\left( {{k}_{n, z}} \right) \right|}_{{{k}_{\rho }}\to \infty }}\text{=}\frac{{{k}_{\rho }}}{{{k}_{n, z}}}{{{\tilde{\xi }}}_{2}} \\ \end{align} \right.$ (13)

上式中, ${{\tilde{\xi }}_{1}}$和${{\tilde{\xi }}_{2}}$与积分变量无关。(13)式的索末

菲积分可利用如下恒等式解析计算:

$\int_{0}^{\infty }{\frac{{{\text{e}}^{i{{k}_{n, z}}z}}}{{{k}_{n, z}}}{{k}_{\rho }}{{\text{J}}_{0}}\left( {{k}_{\rho }}\rho \right)\text{d}{{k}_{\rho }}}\text{=}\frac{{{\text{e}}^{i{{k}_{n}}r}}}{r}$ (14)

其中 ${{k}_{n}}=\sqrt{i\omega \mu {{\sigma }_{n, \text{h}}}}$

基于(13)式和(14)式, 可进一步将(11)式变为:

$H_{0}^{\text{s}}=\int_{0}^{\infty }{\left[ {{\xi }_{1}}\left( {{k}_{n, z}} \right){{\text{e}}^{i{{k}_{n, z}}z}}+{{\xi }_{2}}\left( {{k}_{n, z}} \right){{\text{e}}^{-i{{k}_{n, z}}z}}-\frac{{{k}_{\rho }}}{{{k}_{n, z}}}\left( {{{\tilde{\xi }}}_{1}}{{\text{e}}^{i{{k}_{n, z}}{{z}_{m}}}}+{{{\tilde{\xi }}}_{2}}{{\text{e}}^{-i{{k}_{n, z}}{{z}_{m+1}}}} \right) \right]{{\text{J}}_{0}}\left( {{k}_{\rho }}\rho \right)\text{d}{{k}_{\rho }}}+\int_{0}^{\infty }{\left( {{{\tilde{\xi }}}_{1}}{{\text{e}}^{i{{k}_{n, z}}{{z}_{m}}}}+{{{\tilde{\xi }}}_{2}}{{\text{e}}^{-i{{k}_{n, z}}{{z}_{m+1}}}} \right)\frac{{{k}_{\rho }}}{{{k}_{n, z}}}{{\text{J}}_{0}}\left( {{k}_{\rho }}\rho \right)\text{d}{{k}_{\rho }}}$(15)

上式中, 第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以内的地层界面。由于不同电阻率和频率条件下仪器探测范围不一, 该边界截断方法可能会导致正演精度不够或计算耗时等问题。为此, 本节提出了一种模型边界自适应匹配方法。

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

表示:

$\delta \text{=}\sqrt{{2}/{\left( \omega \mu {{\sigma }_{n, \text{h}}} \right)}\; }$ (16)

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

$s\approx \prod\nolimits_{n=m+1}^{l}{\frac{2{{\sigma }_{n}}}{{{\sigma }_{n-1}}+{{\sigma }_{n}}}}\exp \left[ -{{z}_{m}}/{{\delta }_{m}}-\sum\nolimits_{n=m+1}^{l}{\left( {{d}_{n}}/{{\delta }_{n}} \right)} \right]$ (17)

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的最小二乘问题, 反演目标函数满足下式:

$C\left( \mathbf{m} \right)=\frac{1}{2}\left[ \left\| \mathbf{S}\left( \mathbf{m} \right)-{{R}_{\text{a}}} \right\|_{2}^{2}+\chi \left\| \mathbf{m}-{{\mathbf{m}}_{\text{ref}}} \right\|_{2}^{2} \right]$ (18)

更新方式如下:

${{\chi }_{t}}\text{=}{{\chi }_{0}}\left\| \mathbf{S}\left( {{\mathbf{m}}_{t-1}} \right)-{{R}_{\text{a}}} \right\|_{2}^{2}$ (19)

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

${{\mathbf{m}}_{t+1}}={{\mathbf{m}}_{t}}-\frac{{{\mathbf{J}}^{\text{T}}}({{\mathbf{m}}_{t}})\left[ \mathbf{S}({{\mathbf{m}}_{t}})-{{R}_{\text{a}}} \right]+{{\chi }_{t}}({{\mathbf{m}}_{t}}-{{\mathbf{m}}_{\text{ref}}})}{{{\mathbf{J}}^{\text{T}}}({{\mathbf{m}}_{t}})\mathbf{J}({{\mathbf{m}}_{t}})+{{\chi }_{t}}\mathbf{I}}$ (20)

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

4 应用实例

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

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

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

图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; $H_{0}^{\text{s}}$— — 与0阶贝塞尔函数相关的散射磁场强度, A/m; $H_{1}^{\text{s}}$— — 与1阶贝塞尔函数相关的散射磁场强度, A/m; I— — 单位矩阵; J— — 雅克比矩阵; J0— — 0阶贝塞尔函数; J1— — 1阶贝塞尔函数; k— — 波数, 无因次; kρ — — 积分变量; kn, h, z, kn, v, z— — 第n层地层水平和垂直方向的波数; m, n, l— — 地层序号; m— — 待反演的参数矢量; N— — 模型中地层总数; mref— — 已知参考模型; Ο — — 取高阶无穷小函数; Ra— — 实测视电阻率数据, Ω · m; RP33H— — 频率2 MHz、源距0.838 2 m(33 in)的相位差电阻率, Ω · m; $R_{n, n+1}^{\text{TE, h}}$, $R_{n, n+1}^{\text{TE, v}}$— — zn处水平和垂向源激发的上行TE波的狭义反射系数; $R_{n, n+1}^{\text{TM, h}}$, $R_{n+1, n}^{\text{TM, h}}$— — zn处上行和下行TM波的狭义反射系数; $\tilde{R}_{n, n+1}^{\text{TM, h}}$, $\tilde{R}_{n+1, n}^{\text{TM}, \text{h}}$— — zn处上行和下行TM波的广义反射系数; $\tilde{R}_{n, n+1}^{\text{TE, h}}$, $\tilde{R}_{n, n+1}^{\text{TE}, \text{v}}$— — zn处水平和垂向源激发的上行TE波的广义反射系数; r— — 接收点距发射点的距离, m; s— — 源激发的电磁场在界面处的衰减率; S(m)— — 正演响应, Ω · m; t— — 迭代次数; $T_{n, n+1}^{\text{TM, h}}$— — zn处上行TM波的狭义透射系数; z— — 垂向坐标, m; zn— — 第n个地层界面的垂向坐标, m; σ n, h, σ n, v— — 第n层地层水平和垂向电导率, S/m; ρ — — 径向距离, m; χ — — 正则化系数; χ 0— — 初始正则化系数; δ — — 趋肤深度, m; θ — — 仪器轴与地层界面法向的相对夹角, (° ); ω — — 角频率, rad/s; μ — — 自由空间的磁导率, H/m; λ — — 地层电阻率各向异性系数; ξ 1, ζ 1— — $H_{0}^{s}$和$H_{1}^{s}$积分中上行波的幅度项; ξ 2, ζ 2— — $H_{0}^{s}$和$H_{1}^{s}$积分中下行波的幅度项; $\xi _{1}^{1}$, $\xi _{2}^{1}$— — 一次反射/透射波的幅度; ${{\tilde{\xi }}_{1}}$— — 谱域上行波扣除项系数; ${{\tilde{\xi }}_{2}}$— — 谱域下行波扣除项系数; $\left\| \cdot \right\|_{2}^{2}$— — 求向量的2范数。

(编辑 黄昌武)

参考文献
[1] 高杰, 辛秀艳, 陈文辉, . 随钻电磁波电阻率测井之电阻率转化方法与研究[J]. 测井技术, 2008, 32(6): 503-507.
GAO Jie, XIN Xiuyan, CHEN Wenhui, et al. Resistivity derivation in electromagnetic wave propagation resistivity logging while drilling[J]. Well Logging Technology, 2008, 32(6): 503-507. [本文引用:1]
[2] HU Song, LI Jun, GUO Hongbo, et al. Analysis and application of the response characteristics of DLL and LWD resistivity in horizontal well[J]. Applied Geophysics, 2017, 14(3): 351-362. [本文引用:1]
[3] WANG Lei, FAN Yiren. Fast inversion of logging while drilling azimuthal resistivity measurements for geosteering and formation evaluation[J]. Journal of Petroleum Science and Engineering, 2019, 176: 342-351. [本文引用:1]
[4] 范宜仁, 胡旭飞, 邓少贵, . 倾斜层理地层随钻电磁波测井响应特征[J]. 石油勘探与开发, 2019, 46(4): 675-683.
FAN Yiren, HU Xufei, DENG Shaogui, et al. Logging while drilling electromagnetic wave responses in inclined bedding formation[J]. Petroleum Exploration and Development, 2019, 46(4): 675-683. [本文引用:2]
[5] LI Qiming, LIU Chengbing, MAESO C, et al. Automated interpretation for LWD propagation resistivity tools through integrated model selection[J]. Petrophysics, 2004, 45(1): 14-26. [本文引用:2]
[6] WANG Hanming, GARROW P, SILVA P. Correcting LWD propagation resistivity in horizontal sinusoidal wells for better water saturation estimate: Case study in a deep water, channelized turbidite reservoir[R]. New Orland s, Louisiana, USA: SPE Annual Technical Conference and Exhibition, 2013. [本文引用:1]
[7] PARDO D, TORRES-VERDÍN C. Fast 1D inversion of logging-while-drilling resistivity measurements for improved estimation of formation resistivity in high-angle and horizontal wells[J]. Geophysics, 2015, 80(2): E111-E124. [本文引用:1]
[8] ZHOU Qiang. Uncertainty in geosteering and interpretation of horizontal wells: The necessity for constraints and geometric models[J]. The Leading Edge, 2015, 34(5): 492-499. [本文引用:1]
[9] HONG Decheng, XIAO Jiaqi, ZHANG Guoyan, et al. Characteristics of the sum of cross-components of triaxial induction logging tool in layered anisotropic formation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(6): 3107-3115. [本文引用:3]
[10] GAO Jie, XU Chenhao, XIAO Jiaqi. Forward modeling of multi-component induction logging tools in layered anisotropic dipping formations[J]. Journal of Geophysics and Engineering, 2013, 10(5): 54-67. [本文引用:1]
[11] MICHALSKI K A, ZHENG Dalian. Electromagnetic scattering and radiation by surfaces of arbitrary shape in layered media(I): Theory[J]. IEEE Transactions on Antennas and Propagation, 1990, 38(3): 335-344. [本文引用:1]
[12] 姚东华, 汪宏年, 杨守文, . 用传播矩阵法研究层状正交各向异性地层中多分量感应测井响应[J]. 地球物理学报, 2010, 53(12): 3026-3037.
YAO Donghua, WANG Hongnian, YANG Shouwen, et al. Study on the responses of multi-component induction logging tool in layered orthorhombic anisotropic formations by using propagator matrix method[J]. Chinese Journal of Geophysics, 2010, 53(12): 3026-3037. [本文引用:1]
[13] ANDERSON W L. Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering[J]. Geophysics, 1979, 44(7): 1287-1305. [本文引用:1]
[14] LI Dawei, WILTON D R, JACKON D R, et al. Accelerated computation of triaxial induction tool response for arbitrary deviated wells in planar-stratified transversely isotropic formations[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(6): 902-906. [本文引用:1]
[15] 王磊, 范宜仁, 袁超, . 随钻方位电磁波测井反演模型选取及适用性[J]. 石油勘探与开发, 2018, 45(5): 914-922.
WANG Lei, FAN Yiren, YUAN Chao, et al. Selection criteria and feasibility of the inversion model for azimuthal electromagnetic logging while drilling(LWD)[J]. Petroleum Exploration and Development, 2018, 45(5): 914-922. [本文引用:2]
[16] WU Zhenguan, FAN Yiren, WANG Jiawei, et al. Application of 2. 5-D finite difference method in Logging-While-Drilling electromagnetic measurements for complex scenarios[J]. IEEE Geoscience and Remote Sensing Letters, 2020, 17(4): 577-581. [本文引用:1]
[17] WU Zhenguan, DENG Shaogui, HE Xuquan, et al. Numerical simulation and dimension reduction analysis of electromagnetic logging while drilling of horizontal wells in complex structures[J]. Petroleum Science, 2020, 17(3): 645-657. [本文引用:1]
[18] GUPTASARMA D, SINGH B. New digital linear filters for Hankel J0 and J1 transforms[J]. Geophysical Prospecting, 1997, 45: 745-762. [本文引用:1]