Advances of the inversion of radial slowness profile near borehole with array acoustic logging

Zi WANG; WenZheng YUE; YuMing ZHU; Xin LIU; Bin TIAN

doi:10.6038/pg2024HH0046

Progress in Geophysics >

2024 , Vol. 39 >Issue 5: 2002 - 2012

DOI: https://doi.org/10.6038/pg2024HH0046

Advances of the inversion of radial slowness profile near borehole with array acoustic logging

Zi WANG ^,¹ ,
WenZheng YUE ^,¹^,* ,
YuMing ZHU ¹ ,
Xin LIU ¹ ,
Bin TIAN ²

Expand

¹ State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China
² Testing Company of Shenyang Oil Production Plant of Liaohe Oilfield, CNPC, Shenyang 110316, China

Received date: 2023-04-17

Online published: 2024-12-19

Copyright

Fold

Abstract

The slowness information of near-wellbore formation is of great significance to oil-gas exploration and development, which can reflect the properties of rock. However, the near-wellbore formation presents radial heterogeneity with a variable slowness profile, due to the effects from mud filtrate invasion, wellbore stress, drilling machinery damage and fracturing operation. Accordingly, extracting radial slowness profile from acoustic logging inversion has important application value in identifying invasion depth, evaluating wellbore stability, guiding well completion design, and evaluating fracturing effect. Therefore, we have investigated the current inversion methods constructed by researchers in this paper and divided them into three categories: the inversion with arrival time of waves, the inversion with slowness-time coherence (STC), and the inversion with wave dispersion. After comprehensive research and analysis, we found that the usability and performance of the methods is different due to the difference of target mode waves and their algorithms. Based on the obversions, we further concluded the advances of the methods in extracting the slowness profile near a wellbore from array acoustic logging data, and compared their performances in real applications.

Key words： Radial depth; Slowness profile; Tomographic inversion; Array acoustic logging

Cite this article

Zi WANG , WenZheng YUE , YuMing ZHU , Xin LIU , Bin TIAN . Advances of the inversion of radial slowness profile near borehole with array acoustic logging[J]. Progress in Geophysics, 2024 , 39(5) : 2002 -2012 . DOI: 10.6038/pg2024HH0046

0 引言

由于井眼应力集中、钻井机械损坏、泥浆侵入等外部因素(Baker and Winbow, 1988)的影响，导致井周地层慢度发生径向变化.探测近井地层慢度径向变化具有重要意义(Pistre et al., 2005; 乔文孝等，2011)，因此，国内外学者们针对这个问题开展了诸多研究和探索(王彦仓等，2014; 陈春林等，2023)，在测井尺度上有如下几点：岩石脆性识别(赵龙等，2014; Tang et al., 2016a, b; 钱玉萍等，2022)，井周地层各向异性评价(Plona et al., 2002; 曾富强，2019)，井壁机械破坏评价(Hornby, 1993; Winkler, 2005)，压裂评价(张国栋等，2016; 钱玉萍等，2018; 王文文等，2019)，指导优化完井设计(Sinha et al., 2005, 2006; 赵龙等，2014)，井周应力评价(Sayers et al., 2008; Tang et al., 2016a; 雪宇超，2021)等.国内外学者们在井周轴向、周向地层慢度信息提取方面已经做了大量硏究工作(张海澜等，2004)，而对于径向慢度的信息获取上主要聚焦于基于单极子声源和偶极子声源的慢度径向分布反演.

基于单极子声源的井周地层慢度径向分布反演方法发展较早，Hornby(1993)开发了一种基于同步迭代重建技术(Dines and Lytle, 1979)和射线追踪理论的慢度径向剖面反演技术，利用滑行纵波首波到时来计算井眼附近地层的纵波慢度，对到时的准确性有较高要求(Chi et al., 2004).夏克文等(1994)提出了利用反射波到时信息层析成像的方法，将阻尼代数重建技术作为算法.王乃星等(1998)在此基础上研究了全波列测井中反射波的分离和成像问题.李昌彪等(2004)进一步研究了层析成像技术，用离散图像重建技术，将井眼径向层析成像问题转化为线性方程组的求解，计算速度较快，但仍然没有解决反射波初至拾取问题.Zeroug等(2006)提出了另一种基于射线追踪的方法，能够以最小的用户交互实现实时实施反演，但这种方法仅适用于慢度随径向距离单调变化的地层.张国栋等(2016)采用最小二乘法求取满足计算走时和实测走时之间相对误差最小的模型，这种反演方法精度较高但计算效率较低，不适用于井场.Liu等(2021)提出逐步反演方法，使用慢度-时间相关法(Kimball and Marzetta, 1984)拾取径向慢度序列，判断地层是否径向变化，后基于射线理论，但是慢度拾取结果受人为因素影响.

基于偶极子声源的井周地层慢度径向分布反演方法适用性更广，有效解决了单极子声源在慢地层中难以探测到横波信号的问题.Backus和Gilbert(1970)、Burridge和Sinha(1996)提出了将微扰法和BG理论相结合的方法来得到地层横波慢度的径向变化，这种方法牺牲了分辨率(Sinha, 2004)，且井眼附近反演误差较大(Sayers et al., 2008).Tang和Patterson(2010)开发了一种横波慢度径向剖面的约束反演方法，对径向侵入厚度和与原状地层慢度差这两个参数建立理论模型和目标函数，通过求取目标函数的最小值即得到慢度径向剖面，但会导致径向位置产生误差.Yang等(2012)在此基础上研究了套管井中的径向慢度剖面反演方法，函数主要由3个参数控制，以此建立目标函数.该方法同样使用指数函数去近似来表示横波慢度径向分布，结果表明反演结果严重依赖初始参数设置.马明明(2013)通过研究发现噪声和参数误差对微扰法结合BG理论的反演方法径向慢度的影响较小，并分析了如何在反演迭代过程中将计算效率提高(马明明等, 2013).赵龙等(2014)针对使用约束反演方法在实际测井中会遇到高频缺失情况，提出一种改进方法，利用理论频散数据来进行高频约束，通过数值模拟对其正确性进行了验证.叶正伟(2020)考虑了测井仪器对反演结果带来的影响，使用仪器等效理论，消除了仪器对反演结果的影响.Wang和Yue(2024)开发了一种新的反演方法，可与有效反演层厚度小于阵列声波测井仪器接受阵列长度的薄层，提高了慢度剖面反演的轴向分辨率.

本文在分析不同算法的基础上，总结了近年来井周地层慢度径向分布的反演方法，对具有代表性的几种反演方法进行详细说明，介绍其方法原理和反演效果，结合相关实例应用分析，对比分析了各反演方法的优势与适用性，并提出未来研究趋势.

1 基于首波到时信息的径向慢度剖面反演方法

1.1 基于滑行波首波到时的反演方法

Hornby(1993)提出通过折射纵波首波到时来对井壁附近慢度的二维图像进行重建，采用同步迭代重建技术(SIRT)进行反演，并结合射线理论假设以及FDM分析来计算.这种方法开发了一种根据声波传播时间来给出井眼附近和远离井眼地层的速度公式，其基本原理与地震层析成像相同，给出径向慢度与走时的关系式：

其中t_k代表滑行波的走时，g_k是一个非线性的函数，u(x)代表未知待求的慢度函数，x=(xr, xz)表示地层中的径向和垂向位置，T_k(u)表示发射器和接收器间所有的路径，ds是弧长.

因为Hornby方法中T_k(u)对期望结果慢度函数u(x)有复杂的依赖性，所以该反演过程需要线性化.由此，Zeroug等(2006)在Hornby研究的基础上提出了另一种反演方法，将估计的走时与井周慢度变化联系起来，考虑一个径向分层的井筒，s为井眼半径的尺寸，V_i是和H_i是第i层的速度和厚度，v_f是钻井液速度.在多个走时中取其最小值作为首波到时，以此作为反演径向厚度的条件.射线K到达源距为X的到时(TOF)为：

图 1a为阵列声波测井仪器在径向变化地层中滑行波路径示意图，该类型反演方法对滑行波首波到时准确性有一定要求，反演过程中需要对初始径向慢度模型进行迭代，计算量较大.在定性判断径向慢度变化趋势上，该方法效果良好，图 1b是该反演方法在非常规储层压裂评价中的应用，非压裂段的速度剖面具有很好的重复性，压裂层段(红线区域)在压裂前、后径向速度剖面差异明显，判断压裂效果良好.

显示原图|下载原图ZIP|生成PPT

图1 基于滑行波首波到时的反演方法(钱玉萍等, 2018)

(a)非均质地层中滑行波路径示意图；(b)该方法在实际资料中评价压裂效果的应用.

Figure 1 The inversion with arrival time of slide waves (Qian et al., 2018)

(a)Schematics of slide wave ray paths in a heterogeneous formation; (b)Application of the method in evaluating fracturing effect in actual data.

1.2 基于单极声源反射波首波到时的反演方法

夏克文等(1994)提出利用阵列声波反射波首波到时进行径向慢度剖面成像，利用反射波初至来反演慢度径向剖面基本思想是：把地层划分成若干个网格，每个网格具有独立的慢度值，假设地层在接收器和发射器之间的这段距离内为均质地层，得到反射波的传播路径，计算射线在各网格内的传播距离，得到反射波走时，对比计算得到的走时和仪器接受的反射波初至，当两者的差别小于某个阈值，认为各单元格慢度能够表征该地层的实际慢度，反之则需要不断改变各单元格慢度，直到计算得到的走时与拾取的反射波初至之间的误差满足既定阈值要求.

图 2a为阵列声波测井仪器利用反射波首波到时反演径向速度的示意图.测井仪器声系由发射器和接收器阵列组成.在井眼位置建立坐标系，z轴方向向上，x轴为井眼的径向.反射波在通过侵入带地层时，不同接受阵列获得的走时不同，对初始空间网格进行赋值后不断迭代修改速度模型，最终获得径向速度剖面，如图 2b所示.

显示原图|下载原图ZIP|生成PPT

图2 基于反射波首波到时的反演方法

(a)阵列声波测井反射波路径；(b)成像区域网格化图.

Fig 2 The inversion with the arrival time of reflection wave

(a) Reflection wave path of array acoustic logging; (b)Imaging area grid.

成像的目的是利用反射初至信息确定井眼径向慢度图像f=(x, z)，接收器阵列得到的反射初至波走时数据满足方程：

式中L为反射初至波传播路径，ζ、θ分别为与射线的位置和波的入射角有关的量.阵列声波测井不同于全方位扫描的医学CT，因为这种测井方式是在有限角度下进行，射线对目的层覆盖不均匀, 所测的数据是不充分的.

采用离散图象重建技术，将所测区域目标体模型剖分成n个互不重叠的网格，如图 2b所示，将式(3)作线性化处理得到：

式中d_i代表反射波初至走时的观测值，a_ij代表第i条射线在第j个网格中的长度，f_i为所求慢度值.式(4)可以用矩阵形式d=Af(A∈R^m×n, d∈R^m)表示.解决井眼径向层析成像问题转化成求解f=(f₁, f₂, …f_n)^T.常用的求解方法有奇异值分解法、阻尼最小二乘法、代数重建技术、联合迭代技术、共轭梯度算法、最小二乘法、QR分解等(许令周和范宜仁，2011)，本文不作过多介绍.

该类型方法能够有效识别径向速度异常，且反演精度高.图 3a为层状异常识别速度模型，模型大小是3×3的九个网格，网格颜色代表模型纵波速度，单位为m/s.经过1805次迭代，得到反演结果如图 3b所示.结果表明反演速度值与理论模型误差较小，误差分析如图 3c所示，最小相对误差为0.175%, 最大相对误差为2.35%, 反演效果令人满意.该方法在处理实际资料中存在两个问题，其一是提取实际资料中与井壁平行的反射波存在局限，其二是速度模型较大时迭代速度较慢，不适用井场环境.

显示原图|下载原图ZIP|生成PPT

图3 基于反射波首波到时的概率加权ART反演方法(许令周和范宜仁，2011)

(a)理论速度模型; (b)实际反演结果; (c)相对误差分析.

Fig 3 Probability weighted ART inversion method with arrival time of reflection wave (Xu and Fan, 2011)

(a)Theoretical velocity model; (b) Inversion results; (c) Relative error analysis.

2 基于慢度-时间相关(STC)的逐步反演方法

Liu等(2021)提出一种基于分步反演的地层慢度径向分布成像方法，该方法利用源距在一定范围内增加时，探测深度也会对应增大这一特点，建立了接收器组数和径向位置的联系(耿祎晗等，2017).不同于同时反演慢度和径向位置两个参数的传统方法，该研究的分步反演过程分为两个步骤进行：首先对声波阵列的接收器组合依次采用慢度-时间相关(STC)法(Kimball and Marzetta, 1984)提取不同组的慢度序列，然后基于射线理论(Baker, 1984; Baker and Winbow, 1988; 耿祎晗等, 2017)获得对应速度的层厚.

慢度-时间相关法是从阵列声波测井资料获得慢度的有效方法，其相关系数公式为：

式中ρ(s, T)为慢度-时间相关法的相似系数，s为慢度值，T代表接收器接收信号的时间，T_w是时间窗的长度，X_m(t)代表第m个接收器在时刻t所接收到的声波信号，d为接收器之间的距离，N为接收阵列总数.实际计算时第一步需要选取n∈[2, N－1]个接收阵列为一组，获取慢度序列，判断是否为均质地层.如果为侵入地层，此时需要进行第二步根据射线理论确定相应的层厚.滑行波传播路径T_pn和层厚H_i的关系表达式为：

根据式(6)可得，对于固定源距z，第n层的临界厚度H_n为：

当源距处于z₁和z₂(z₁＜z₂)时H_n1＞H_n, H_n2＜H_n，当n＞2时用两层层厚的平均值即

作为临界厚度，当n=1时，H₁为：

使用该方法对不同层厚(0.2 m、0.4 m、0.6 m和0.8 m)的非均质地层(地层纵波速度从4000~5000 m/s)模拟数据进行验证.分步反演的第一步是建立速度和组数的关系，图 4a是使用2个接收器作为一组，拾取不同组慢度值的结果，拾取方法不需要首波到时，受噪声影响较小.分步反演的第二步是建立速度和径向位置的关系，图 4b中红色线段是通过式(7)和式(8)得到的速度和径向位置的关系，蓝色线段是理论速度模型，对比发现该方法反演准确率较高，且反演速度较快，但是在拾取慢度时，结果一定程度上受到人为因素影响(窗长、接收器组合)，且需要长源距仪器提供足够多的数据(Chi et al., 2004).

显示原图|下载原图ZIP|生成PPT

**图4 不同层厚的非均质地层模型反演结果(Liu et al., 2021)**

(a)拾取的纵波速度序列；(b)纵波速度径向剖面.

Fig 4 Inversion results of a heterogeneous formation model with different thickness of velocity variation region(Liu et al., 2021)

(a)The sequence of extracted P-wave velocities; (b)Radial profiles of P-wave velocity.

3 基于频散的径向慢度剖面反演方法

3.1 基于偶极声源的微扰法结合BG理论反演方法

Burridge和Sinha(1996)通过结合微扰法和BG理论，提出了一种反演地层横波慢度径向分布的方法，通过微扰法得知地层剪切模量的变化对弯曲波相速度的影响最大，反演过程可以将地层相速度的变化看成完全由剪切模量的扰动引起.此方法的基本思想是：选定一个已知横波速度v_s的参考模型，计算其理论频散曲线v_refer，实际测得的弯曲波频散曲线v_real与参考模型的差异δv=v_real－v_refer，此差异视为两者横波慢度在径向上的差异引起.建立δv和δv_s(r)的关系，则实际横波径向慢度分布即v_s(r)=v_s+δv_s(r).反演主要公式为：

其中μ为参考模型剪切模量，δv_i/v_i为实际模型与参考模型弯曲波相速度相对差，r表示径向位置，δμ(r′)/μ为r′处横波慢度的相对变化.G_i(r)为待求解的核函数.式(9)将已知δμ(r′)/μ求解δv_i/v_i的正演问题和已知δv_i/v_i求解δμ(r′)/μ的反演问题统一，即可得到径向慢度分布.

图 5a是通过加权频谱相干法(唐晓明和郑传汉，2004)径向非均质性地层(横波速度1000~1200 m/s)的偶极横波频散曲线.高频段弯曲波因衰减原因信号弱，频散曲线出现异常点.图 5b是根据式(9)求解的径向速度剖面，因为偶极横波的高频部分携带近井地层的信息，高频缺失会导致该方法在近井眼处(径向深度0.2 m内)误差较大.

显示原图|下载原图ZIP|生成PPT

图5 BG理论法反演结果(叶正伟，2020)

(a)频散拾取结果；(b)径向速度剖面.

Fig 5 Inversion result of BG theory inversion method(Ye, 2020)

(a)Frequency dispersion result; (b)Radial velocity profile.

3.2 基于偶极声源的高频约束反演方法

Tang和Patterson(2010)为了解决近井眼误差较大的现象，提出了一种高频约束的反演方法，通过对井眼附近的高频波进行约束，使得近井眼处反演精度提高.他利用频散数据在不同频率的探测深度不同构造反演目标函数，用于确定径向速度变化量和径向变化地层厚度，反演目标函数为：

其中Δr和ΔV表示径向变化地层的厚度和速度变化量，V_m是以Δr和ΔV为函数模型的频散曲线，V_d是实际资料的频散曲线，Ω为给定的频率范围(2~10 kHz)，E是各个频率上残差的总和.式(10)中的第二项用于求解非唯一性问题的约束函数，其中V_h是利用井眼附近的地层速度V₁得到的均质频散曲线，Ω′是Ω的高频子集(9~10 kHz)，λ为权重因子通常取4.地层横波速度发生径向变化，通常为单调递增或者单调递减.根据这种现象建立地层径向变化模型，即：

其中，V₀是原状地层的横波速度，r₀为井眼半径.将式(10)所求得的Δr和ΔV带入式(11)即可得到横波径向慢度剖面，如图 6所示.对比图 6和图 5b的反演结果不难发现BG理论在近井眼处存在横波速度先降低再升高的现象(黑创，2013)，而高频约束法则很好的反应了地层横波速度不断增大的特点，表明高频约束方法有效地解决了近井眼处反演准确率不足的问题.

显示原图|下载原图ZIP|生成PPT

图6 高频约束法径向速度剖面

Fig 6 Radial velocity profile of high-frequency constrained method

该方法可用于估计岩石脆性和可破碎性(Tang and Patterson, 2016)，图 7为高频约束反演方法在致密气层的综合应用实例(赵龙等，2014)，图中有上、下两个储集层分别是在4740~4758 m和4776~4808 m，孔隙度曲线表明储层为低孔隙致密储层.图中第4道是理论计算的岩石脆性曲线，其中下部储层4776~4808 m脆性较大.第5道给出利用偶极横波反演的横波速度径向变化剖面，从图中可看出, 在4740~4758 m井附近地层径向横波速度变化较小，而在4776~4808 m附近地层横波速度径向变化特征较明显，径向变化深度约25 cm，对比发现地层径向横波速度变化与岩石的脆性有很好的对应关系，说明横波径向慢度剖面可以有效指示岩石的脆性，有效指导压裂施工.4776~4808 m储集层压裂前后各向异性明显增大，压裂效果理想，说明横波径向慢度剖面可以对压裂效果进行评价.

显示原图|下载原图ZIP|生成PPT

图7 致密地层横波速度剖面与多种测井结果对比示意图(赵龙等，2014)

Figure 7 Comparison of radial shear-wave velocity profile with several logging results (Zhao et al., 2014)

4 结论及展望

论文阐述了目前主流慢度径向分布反演方法，分析不同反演方法的基本原理、求解思路、过程以及得到的结果，评价了每种方法的优势与不足.综合分析不难发现：

(1) 基于单极子声源的反演方法趋于借助源距和探测距离正相关的特性，建立慢度与径向位置的联系.基于偶极子声源反演方法趋于借助频散曲线低频部分对应原状地层慢度，高频部分对应近井地层慢度的特性，建立慢度与径向位置的联系.

(2) 慢度径向分布反演的难点不在于慢度拾取，而在于将拾取的慢度值和径向位置对应.基于单极声源的反演方法在径向位置确定上通常通过射线追踪作为反演方法，基于偶极子声源则是通过实际频散曲线和理论频散曲线的差异与径向慢度变化建立联系，确定径向位置.

(3) 基于反射波到时层析成像方法有着完善的反演理论，但是反演要求高质量反射波.目前反射波提取方法常用于远探测地质体反演，平行于井轴的反射波难以被提取，所以该方法难以反演实际测井资料.

(4) 目前反演方法为了保证反演稳定性和提升反演效率，或多或少牺牲了部分径向分辨率，使得反演结果是一个平均值，能准确反应慢度变化的趋势，但是特定径向位置的慢度会有一些误差.

随着我国瞄准深地前沿领域，向地球深部进军这一战略问题，地层径向慢度信息在井周应力评价、井眼破损评价、岩性识别等方面发挥着重要作用，笔者提出几点展望：

(1) 目前慢度径向分布反演以二维为主，随着测井仪器不断发展，接收阵列能够提供具有方位信息的数据，开发三维慢度径向分布反演方法可以更全面的评价井周的地层特性.

(2) 源距和探测深度正相关，通常测井仪器的源距径向最大探测深度在不超过1 m，通过开发长源距测井仪器或组合不同源距的测井资料，获得更长源距的数据，可以增加径向探测深度以获得更丰富的地层信息.

(3) 逐步反演法给慢度径向分布反演提供了一个新的思路，即慢度和径向位置不需要同时反演.基于这个想法，单极子声源的反演方法可以尝试将源距和探测径向深度做标定，找到源距与径向位置的关系；偶极子声源的反演方法可以在主频一定的情况下，尝试找到频散曲线中频率和径向位置的联系.这样可以避免迭代，以最小的用户交互快速反演出径向剖面.

感谢各位审稿人对论文提出的建设性建议和意见.

References

Publishing order | Descend order by publishing year | Descend order by cited within

Backus

, Gilbert

. Uniqueness in the inversion of inaccurate gross earth data. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1970, 266(1173 123 192.

模态框（Modal）标题

Abstract

Cite this article

0 引言

1 基于首波到时信息的径向慢度剖面反演方法

1.1 基于滑行波首波到时的反演方法

图1 基于滑行波首波到时的反演方法(钱玉萍等, 2018)

1.2 基于单极声源反射波首波到时的反演方法

图2 基于反射波首波到时的反演方法

图3 基于反射波首波到时的概率加权ART反演方法(许令周和范宜仁，2011)

2 基于慢度-时间相关(STC)的逐步反演方法

图4 不同层厚的非均质地层模型反演结果(Liu et al., 2021)

3 基于频散的径向慢度剖面反演方法

3.1 基于偶极声源的微扰法结合BG理论反演方法

图5 BG理论法反演结果(叶正伟，2020)

3.2 基于偶极声源的高频约束反演方法

图6 高频约束法径向速度剖面

图7 致密地层横波速度剖面与多种测井结果对比示意图(赵龙等，2014)

4 结论及展望

References

**图4 不同层厚的非均质地层模型反演结果(Liu et al., 2021)**