2019, Vol. 62
声波测井的主要任务之一是为了获得地层的纵横波速度，以满足地震勘探、岩石力学分析等其他地质工程方面的应用.现代多极子阵列声波测井仪器可以同时采集井孔中的单极和偶极阵列波形数据，但在提取地层波速时必须考虑井孔波导频散效应对数据的影响.近年来，随着深水资源的勘探与开发，需要在浅层欠压实的超慢地层(本文将横波速度小于1000 m·s-1的地层称为超慢地层)中测井.这时需要同时考虑多极子阵列波形数据中的偶极弯曲波和单极纵波的频散效应.慢地层中的井孔纵波，由于向地层辐射横波，会造成井孔方向传播的衰减，并伴随有较强的频散效应，称为泄漏纵波(唐晓明和郑传汉，2004).尤其是在重泥浆、大井径井眼中，泄漏纵波的频散效应更为显著(Hornby and Pasternack, 2000).
为了从阵列波形数据中准确获得地层的速度，国内外学者展开了大量的研究，以校正数据中的频散效应.这些方法主要分为两类：数据驱动方法(Huang and Yin, 2005; Zheng et al., 2006; Tang et al., 2010; 孔凡童等, 2017)和模型驱动方法(Tang et al., 1995; Kimball, 1998; Geerits and Tang, 2003; Lee et al., 2016; Su et al., 2016).数据驱动方法需要从阵列声波数据中提取频散数据，通过寻找低频(模式波截止频率附近)“极限”，获得地层的纵横波慢度，但受仪器测量频带的限制，低频数据往往信噪比较低.模型驱动方法一般采用井孔波导理论计算的频散曲线拟合测量数据中提取的频散曲线，对数据有效频带之外的低频段具有预测和外推的优势(Lee et al., 2016).其他一些模型驱动方法(Tang et al., 1995; Kimball, 1998)直接采用理论计算的频散曲线匹配频域数据的相位，避免了从数据中提取频散曲线的过程，但这些早期的方法没有考虑测井仪器的影响.Lee等(2016)指出了现场数据处理中测井仪器对偶极弯曲波频散效应的重要影响.
以上这些模型驱动方法，往往只考虑偶极弯曲波的频散效应，忽略了泄漏纵波的频散效应，而在纵波慢度不准确时，反演得到的横波慢度也存在一定的不确定性.Su等(2016)提出了一种随钻多极子声波数据的联合反演方法，可以同时校正随钻纵波和横波数据的频散效应，从而得到准确的纵横波慢度，但相对于相位匹配方法，该方法需要从波形数据中提取频散数据，反演结果受所提取频散数据可靠性的限制.
本文根据仪器等效理论和相位匹配方法，提出了一种对多极子阵列声波测井数据的联合反演处理方法.该方法同时考虑了单极和偶极数据的频散效应以及测井仪器在其中的影响，采用仪器等效模型计算的多极子模式波频散曲线，同时匹配单极和偶极波形数据的相位，以此来构建反演的目标函数.这种新方法可以直接处理多极子阵列声波数据，且不需提取频散数据，不仅能够同时获得地层的纵横波慢度，而且可以减少纵横波速度反演的不确定性，为多极子阵列声波测井数据处理提供了一种高效、便捷的处理方法.
1 多极子模式波灵敏度分析多极子阵列声波测井数据的频散效应是井孔、地层、仪器等参数共同作用的结果.基于井孔波导理论的模型驱动方法的关键是含仪器井孔波导多极子模式波理论频散曲线的计算.为此人们提出了一些描述声波仪器的理论模型(White，1983；Chen and Willen, 1984；唐晓明和郑传汉，2004；Sinha et al., 2009；Su et al., 2011)，其中比较切实可行的是唐晓明和郑传汉(2004)、Su等(2011)提出的仪器等效理论模型，不需要考虑仪器内部结构，仅使用两个参数(仪器外径和仪器等效模量)就能模拟测井仪器对多极子波声场响应和频散特征的影响(Su et al., 2011).
上述仪器等效理论模型提供了描述测井仪器对井孔声场响应影响便捷而有效的方法.均匀无限大地层，充液井中，仪器居中放置(图 1)时，其激发的井中多极子模式波的频散方程为：
|
图 1 多极子阵列声波测井模型图 Fig. 1 Model for logging in an earth formation with a multipole acoustic tool |
|
(1) |
其中，k为轴向波数；ω为圆频率；F表示各向同性地层的弹性参数，包括纵、横波速度和密度；B表示井孔流体的声学参数，包括流体速度、密度和井径.声波测井仪器的响应利用仪器等效半径a和等效模量MT表示.Det表示矩阵Dij(i, j=1, 2, 3, 4；表达式见附录)的行列式.对于给定的频率ω，通过求解频散方程(式(1))，可以得到井中传播多极子(n=0，单极子；n=1，偶极子)模式波的波数k.相慢度Sphase(ω)、相速度Vphase(ω)可以通过(2)式计算得到.
|
(2) |
通过求解的多极子模式波频散曲线，就可以分析该模式波对于模型参数的灵敏度.为了说明地层、井孔、仪器等参数对于多极子模式波频散效应的影响，首先需要对欠压实、超慢地层中的泄漏纵波、弯曲波进行灵敏度分析(模型参数见表 1).模型对于某一参数的灵敏度实际上就是归一化后的波的相速度对于该参数的偏微分(唐晓明和郑传汉，2004)：
|
表 1 计算多极子模式波的地层及流体参数 Table 1 Parameters of mud and formation for calculating multipole acoustic wav |
|
(3) |
其中，P表示模型的某一参数.
图 2a给出了采用仪器等效理论模型计算的泄漏纵波和弯曲波的频散曲线.图 2b给出了泄漏纵波对于模型参数的灵敏度，描述了地层纵波速度(VPfm)、横波速度(VSfm)、井内流体速度(VPfd)和仪器等效半径(Rtool)对于单极泄漏纵波传播的影响.在低频范围(截止频率附近)内，泄漏纵波主要受地层纵波速度控制；随着频率的增加，其他模型参数的影响逐渐增大；但在高频范围，主要受井内流体速度的控制.图 2c给出了弯曲波对于模型参数的灵敏度.尽管弯曲波的传播主要受地层横波速度的控制，但是其他模型参数的影响也不可忽略.通过上述分析，可以看出，在超慢地层中，泄漏纵波和弯曲波的传播均受到地层—井孔—仪器模型及其参数的影响.因此，可以利用泄漏纵波和弯曲波对模型参数的灵敏度反演该参数，特别是地层的纵、横波速度.
|
图 2 超慢地层理论频散曲线及灵敏度分析 (a)泄漏纵波和弯曲波理论频散曲线; (b)泄漏纵波灵敏度分析; (c)弯曲波灵敏度分析. Fig. 2 Theoretical dispersion curves for a very slow formation and sensitivity analysis (a) Theoretical dispersion curves of leaky-P and flexural wave; (b) Sensitivity analysis for leaky-P wave; (c) Sensitivity analysis for flexural wave. |
通过对泄漏纵波和弯曲波的灵敏度分析可以得知，在超慢地层中，泄漏纵波对地层的纵波和横波速度都具有相当的灵敏度.尽管纵波速度对弯曲波影响较小，但如果该参数未知的情况下，反演的地层横波速度仍有一定的不确定性.本文综合考虑泄漏纵波和弯曲波的频散效应，将两者的灵敏度特性相结合，提出一种针对多极子阵列声波数据的联合反演处理方法.该方法结合Tang等(1995)的相位匹配法与Kimball(1998)的频散慢度相干法(DSTC)，构建反演的目标函数，具体表示为单极纵波相位匹配相干函数与偶极弯曲波相位匹配相干函数的加权和.
|
(4) |
(4) 式中，上(下)角标lp和fl分别表示泄漏纵波和弯曲波；Wi(f)表示第i个接收器的频谱；f为频率；δ是接收器间距；α(f)表示采用仪器等效理论模型计算的相慢度频散曲线；p和s是待反演的地层纵、横波慢度，用来计算α(f)；M和N，分别表示单极和偶极接收器的个数；Ω表示反演时选择的频段；λ(依据单极和偶极数据质量的好坏，λ取值介于0~1之间，一般为0.5)是加权因子，用于调节纵波和弯曲波在目标函数中的权重.反演目标函数的求解可以采用多维下降单纯型法(Press et al., 1992)，不需计算目标函数对于模型参数的导数，从而使极值的求取高效、便捷.值得注意的是，弯曲波是频散波，单极波形数据是否频散取决于地层的类型(快速地层、慢速地层、超慢地层).(4)式中的相干函数(lp或者fl)，如果不对频率进行积分，也可以用来提取阵列波形中波速随频率变化的频散数据(唐晓明和郑传汉, 2004).将利用某一深度点的现场数据(单极、偶极模式仪器工作主频分别为8 kHz、3 kHz)分析数据中的频散效应以及联合反演方法对该效应的有效处理.
图 3给出了该深度点单极、偶极波形数据(图 3a、3c)以及常规STC方法处理结果(图 3b、3d).从STC相干图中，求取极大值，分别得到地层的纵、横波慢度：481 μs·m-1、1073 μs·m-1.从波形图及处理结果可以看出，该地层为软地层，泄漏纵波、弯曲波的频散效应都比较显著.值得注意的是，泄漏纵波的STC处理结果(图 3b)中有两个峰值，左下方的峰值对应低频数据，右上方的峰值对应高频数据.这是因为，低频时泄漏纵波的相速度趋于地层纵波慢度，高频时趋于井内流体的慢度(见图 2a).因此，接下来的反演过程中，可以选取高频泄漏纵波的慢度作为迭代初值，通过最小二乘方法拟合泄露纵波和弯曲波的频散曲线，标定流体的慢度.为了进一步说明频散效应的影响，从实测波形数据中提取了频散数据(图 3e黑色圈线)，并与STC处理结果计算的理论频散曲线(图 3e黑色虚线)进行了对比，显然，STC处理方法得到的纵横波慢度高于实际地层.将利用联合反演方法来处理图 3中的数据.
|
图 3 多极子阵列波形数据、STC处理结果和现场频散数据与处理结果计算频散曲线 (a)泄漏纵波波形数据; (b)泄漏纵波STC处理结果; (c)弯曲波波形数据; (d)弯曲波STC处理结果; (e)现场频散数据(散点)与STC(黑色虚线)、联合反演(黑色实线)结果计算理论频散曲线. Fig. 3 Multipole acoustic wave data, STC processing results and field dispersion data and theoretical dispersion curves calculated with processing results (a) Leaky-P wave data; (b) STC processing result for leaky-P wave; (c) Flexural wave data (d) STC processing result for flexural wave; (e) Field dispersion data (markers) and theoretical dispersion curves calculated with STC method (black dashed line) and inversion result (black line). |
为了分析联合反演方法的特征与优势，我们在图 4a和4b中分别给出了单独采用单极数据和偶极数据反演纵横波慢度的变密度相干图(这可以通过在(4)式中分别令λ=0和1计算得到).相干值越高，说明相位匹配越好.图 4a中显示了一个与纵、横波慢度轴斜交的高值条带区间，这说明泄漏纵波的频散同时受到了纵波和横波慢度的影响(参见图 2b).然而，图 4a中的高值区间狭长，说明单独采用纵波数据反演地层的纵、横波慢度，有较大的不确定性.图 4b中弯曲波的高值区间几近平行于横轴(纵波慢度轴)，说明弯曲波的频散主要受地层横波慢度的影响，地层纵波慢度对其也有微弱影响，但这种影响不可忽略.这也说明单独利用弯曲波数据反演地层纵、横波速度也是比较困难的.
|
图 4 单独反演、联合反演方法处理结果 (a)泄漏纵波单独反演目标函数相干图; (b)偶极弯曲波单独反演目标函数相干图; (c)联合反演目标函数相干图. Fig. 4 The individual inversion and the joint inversion processing result (a) VD (variable density) display of the individual inversion semblance function for leaky-P data; (b) VD display of the individualinversion semblance function for flexural data; (c) VD display of the joint inversion semblance function. |
利用加权因子λ(令式(4)中λ=0.5)将上述单极和偶极数据单独反演的结果相结合，可以得到联合反演方法的处理结果，变密度相干图(图 4c)中显示了具有明显峰值的极大值区间，该峰值接近于1(图 4c中箭头所示)，指示了反演的结果：纵波慢度为465 μs·m-1，横波慢度为995 μs·m-1.为了进一步验证反演结果的准确性，利用反演结果计算泄漏纵波和弯曲波的相慢度频散曲线(图 3e中黑色实线)，并与实际波形数据中提取了频散数据(实际反演过程不需要)进行了对比.理论计算频散曲线与实际数据提取的频散曲线(图 3e中黑色圈线)拟合良好，说明了反演结果的准确性.这一反演实例说明了联合反演方法不仅可以同时反演得到精确的纵、横波慢度，而且不需要从波形数据中提取频散曲线，提高了处理效率.
3 现场数据处理多极子阵列声波数据主要是用来求取地层的纵横波速度剖面.利用联合反演方法来处理现场实测波形数据，包括中速地层和超慢地层数据，以验证该方法的有效性和准确性.在反演之前，应先分析波形数据的频谱特征，选取信噪比较高的频带，以保证反演结果的准确性.
3.1 中速地层处理实例利用联合反演处理方法，对一个井段长度为100 m的中速地层的多极子阵列波形数据进行处理，图 5中第1道分别给出了自然伽马测井曲线(红色)和井径曲线(蓝色).单极、偶极波形数据以变密度的形式在第2道和第4道给出.第3道为单极数据的STC处理结果，蓝色线表示地层纵波慢度，黑色的线表示地层的横波慢度.值得注意的是，该地层为中速地层，单极纵波频散不明显，不需要频散校正，采用常规STC方法就能够得到准确的地层纵、横波慢度，可以用来验证本文联合反演方法的有效性.第5道为偶极数据的STC处理结果，黑色实线表示偶极数据得到的地层横波慢度.第6道和第8道分别给出了利用本文处理方法对单极、偶极数据联合反演得到的地层纵、横波慢度曲线(黑色)，及其反演目标函数的变密度相干图.这一图像是通过在反演目标函数的峰值点(即反演结果)，分别沿着纵波和横波慢度轴取值，然后将这些取值(为一曲线)沿着深度绘图形成.纵、横波慢度曲线很好的追踪了该相干图的极大值，可以作为质量监控的一种方法.为了进一步验证联合反演处理方法的准确性，在第7道给出了联合反演处理方法得到的纵波慢度曲线(红色)与单极数据STC处理结果(蓝色)的对比图，两者之间的差异以百分比的形式显示在第7道的左侧灰色填充部分.从第7道中明显看出，反演得到的结果与STC处理的结果基本重合，误差在1%左右.类似地，在第9道给出了单极数据STC处理得到的横波慢度曲线(蓝绿色)与联合反演得到的横波慢度曲线(红色).单极STC结果基本覆盖反演结果说明了联合反演处理方法确实得到的真实的地层慢度.联合反演处理结果与偶极数据STC处理结果的差异以灰色填充曲线展示在图 5第9道左侧，频散校正量平均约在9%左右.
|
图 5 中速地层多极子阵列声波数据处理实例 Fig. 5 Example of multipole acoustic log data processing for a moderately fast formation |
图 6给出了应用联合反演方法得到的超慢地层多极子阵列波形数据的处理结果，深度间隔为335 m.图 6a中的第1道给出了伽马曲线(红色)和井径曲线(黑色).单极波形和偶极弯曲波数据分别在第2道和第3道以变密度的形式给出.第4道和第6道分别给出了反演目标函数的变密度相干图.第4道中的蓝色曲线表示反演得到的纵波慢度随深度变化的曲线，第六道中的黑色曲线表示反演得到的横波慢度曲线，这两条曲线很好地追踪了相干性图中的极值条带.为了验证该方法对该地层频散数据处理的有效性，在第五道和第七道分别给出了反演得到的纵、横波慢度曲线与STC处理结果的对比图，阴影部分表示反演结果与STC处理结果的差异(以百分比的形式给出)，它们表示了STC处理结果需要进行的频散校正量.因为该地层是欠压实、超慢地层，单极泄漏纵波和偶极弯曲波的频散效应都十分明显.对于泄漏纵波数据的平均频散校正量约为4%~5%，弯曲波的约为9%~10%.为了进一步验证处理结果的准确性，在该井段的上部(X705)和下部(X855m)分别取了两个深度点的处理结果，采用仪器等效理论模型，计算了理论的频散曲线.从图 6b中明显可以看出，利用反演结果计算的理论频散曲线(红色实线)与实测数据提取的频散曲线(黑色圈线)拟合非常好，特别是在低频范围，拟合曲线变得平缓，对应着地层的纵横波慢度.理论模拟与频散数据的对比证明了该处理方法的有效性.
|
图 6 超慢地层多极子阵列声波数据处理实例及验证 (a)超慢地层多极子阵列波形数据处理实例; (b)两深度点处理结果验证. Fig. 6 Example of multipole acoustic log data processing for an unconsolidated slow formation and validation of the processing result (a) Example of multipole acoustic log data processing for an unconsolidated slow formation; (b) Validation of the processing result from two depth point. |
针对多极子阵列波形数据中普遍存在的频散效应，本文提出了一种可以同时获得地层纵、横波慢度的联合反演处理方法.该反演方法中的正演计算采用等效仪器理论得到频散曲线，以此来匹配实测纵波和弯曲横波频域数据中的相位，不仅减少了纵横波慢度反演的不确定性，而且避免了从实测数据中提取频散曲线的过程，简化了处理流程，提高了从频散数据提取地层波速的精度.现场数据的处理实例表明，该联合反演方法有效处理了阵列波形数据中的频散效应，得到了准确的地层纵横波速度，为现代多极子阵列声波测井数据提供了一个行之有效的处理方法.
附录A基于仪器等效理论，多极子模式波频散方程中，矩阵D元素的表达式如下：
|
(A1) |
其中，ρf为井中流体的密度；
|
(A2) |
对于一个未刻槽钢结构测井仪器，等效仪器模量MT≈E/(1-υ)，其E为杨氏模量，υ为泊松比(详见Tang和Cheng，2004)).
Chen S T, Willen D E. 1984. Shear Wave Logging In Slow Formations. Society of Petrophysicists and Well-Log Analysts.
|
Geerits T W, Tang X M. 2003. Centroid phase slowness as a tool for dispersion correction of dipole acoustic logging data. Geophysics, 68(1): 101-107. DOI:10.1190/1.1543197 |
Hornby B E, Pasternack E S. 2000. Analysis of full-waveform sonic data acquired in unconsolidated gas sands. Petrophysics, 41(5): 363-374. |
Huang X J, Yin H Z. 2005. A data-driven approach to extract shear and compressional slowness from dispersive waveform data.//75th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 384-387, doi: 10.1190/1.2144349.
|
Kimball C V. 1998. Shear slowness measurement by dispersive processing of the borehole flexural mode. Geophysics, 63(2): 337-344. DOI:10.1190/1.1444333 |
Kong F T, Zhuang C X, Su Y D, et al. 2017. A data-driven dispersive processing of LWD quadrupole wave data. Chinese Journal of Geophysics (in Chinese) (in Chinese), 60(12): 4938-4944. DOI:10.6038/cjg20171232 |
Lee S Q, Tang X M, Su Y D, et al. 2016. Model-based dispersive processing of borehole dipole wave data using an equivalent-tool theory. Geophysics, 81(1): D35-D43. DOI:10.1190/geo2015-0210.1 |
Press W H, Flannery B P, Teukolsky S A, et al. 1992. Numerical Recipes in C: the art of Scientific Computing. Cambridge: Press Syndicate of the University of Cambridge.
|
Sinha B K, Ikegami T, Johnson D L, et al. 2009. Use of an effective tool model in sonic logging data processing: U. S. Patent, 7529152.
|
Su Y D, Tang X M, Hei C, et al. 2011. An equivalent-tool theory for acoustic logging and applications. Applied Geophysics, 8(1): 69-78. DOI:10.1007/s11770-011-0272-6 |
Su Y D, Jiang C, Zhuang C X, et al. 2016. Joint inversion of logging-while-drilling multipole acoustic data to determine P-and S-wave velocities in unconsolidated slow formations. Geophysics, 81(5): D553-D560. DOI:10.1190/geo2016-0133.1 |
Tang X M, Reiter E C, Burns D R. 1995. A dispersive-wave processing technique for estimating formation shear velocity from dipole and Stoneley waveforms. Geophysics, 60(1): 19-28. DOI:10.1190/1.1443747 |
Tang X M, Zheng C H. 2004. Quantitative Borehole Acoustic Methods (in Chinese). Zhao X M trans. Beijing: Petroleum Industry Press.
|
Tang X M, Li C, Patterson D. 2010. A curve-fitting technique for determining dispersion characteristics of guided elastic waves. Geophysics, 75(3): E153-E160. DOI:10.1190/1.3420736 |
White J E. 1983. Underground Sound, Application of Seismic Waves. Amsterdam: Elsevier Science Publishing Company Inc.
|
Zheng Y B, Huang X J, Tang X M, et al. 2006. Application of a new data-driven processing method to LWD and wireline dispersive compressional and shear waveform data.//76th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, doi: 10.2118/103328-MS.
|
孔凡童, 庄春喜, 苏远大, 等. 2017. 随钻四极横波测量的数据驱动频散处理方法. 地球物理学报, 60(12): 4938-4944. DOI:10.6038/cjg20171232 |
唐晓明, 郑传汉. 2004.定量测井声学.赵晓敏译.北京: 石油工业出版社.
|