反射声波测井波场分离方法研究

宫昊, 陈浩, 何晓, 卫建清, 王秀明

(1.中国科学院声学研究所声场与声信息国家重点实验室, 北京 100190;2.清华大学, 北京 100084; 3.中国工程院, 北京 100088)

0 引 言

隐蔽油气藏储层的特点是岩性横向变化快、非均质性较强。裂缝孔洞型及未完全充填的洞穴型缝洞储层产量较高,但在实际测量中发现,井壁缝洞发育并不代表井旁缝洞发育(反之亦然),给井下隐蔽储层的准确识别带来困难[1]。因此,在缝洞型储层的勘探开发中,需要对井壁较远范围内地层展布与裂缝孔洞的分布情况进行描述。此外,在水平井钻井中,需要对井周储层边界进行标定,在测量条件较好情况下,还期望对油气界面(GOC)进行成像[2]。

常规地面地震与井中地震技术(如垂直地震剖面VSP)拥有很好的测量覆盖区域,但是由于分辨率不足而难以描述上述储层的精细构造。因此,这类井周微裂缝、微断层等也被称为亚地震尺度构造[3]。声波测井具有较高的分辨率,但是探测范围局限在井周围1 m左右,亦难以满足探测需求。对于一切利用声波的探测方法,都需要在分辨率与探测范围之间寻求平衡。因此,只能针对不同的探测目标与勘探的不同阶段,采用不同的声测量方法[4-5]。

反射声波测井技术在这一背景下应运而生,它将声源与接收器布置于同一个井中,通过探测从井旁地质界面反射回来的信号,实现对井周围地质构造进行成像[6-8]。常规声波测井与井中地震之间,存在探测盲区(主要针对分辨率以及探测范围2项探测特性),反射声波测井的出现,消除了这一盲区,扩大了地球物理勘查技术的应用效果与应用范围。

利用常规的地震处理方法对反射波数据进行偏移处理,可以得到地质构造的图像。然而,不同于地震测量方法,反射声波测井的声源和接收器均布置在同一井中。因此,由于工作模式与数据采集方式的不同,给数据采集、信号处理以及后期的测井解释工作都带来了一定的困难[9-11]。一个主要的困难是用于偏移成像的反射信号相对于高幅度井中直达波较为微弱,在准确记录反射信号的前提下,只有利用波场分离技术将反射波从全波形中提取出来,才能对反射信号进行成像处理[12-13]。

反射声波测井数据处理思路一般都借鉴了地震资料处理中较为成熟的滤波方法,或者直接使用地震处理软件。Tang等[12]针对反射声波测井模型的特殊性,分析了反射波的传播特征,在参数估计法[14]的基础上提出一种针对反射声波测井数据的波场分离方法,首先利用参数估计法反演出井中不同振型的直达波频谱,对其进行滤除,进而得到幅度较小的反射信号。Hirabayashi等[15]也提出了类似的滤波方法,对表征井中直达波传播的参数进行反演,构建井中直达波形并进行滤除。Wang等[13]将线性预测滤波法与fk滤波法联合应用于反射声波测井中,并且采用倾角叠加法增强反射信号,得到了更为清晰的反射波信号。目前,普遍的思路是通过引入地震中成熟的数据处理流程,针对反射声波测井作业环境与井孔声场的特殊性,对现有滤波技术进行优化,能够取得较好的滤波效果。实际应用中,普遍存在井中幅度较大的直达波难以被完全压制的问题,残余直达波对井外微弱反射波的识别与提取形成干扰。因此,反射声波测井成像结果中普遍存在成像精度低、可靠性较差、多解性强等问题。针对这一问题,本文对广泛应用于偶极反射声波测井中的线性预测滤波法[12]进行优化,在滤波过程中考虑了偶极测井弯曲波频散的影响,提出一种结合了慢度—时间相干法与线性预测滤波法的波场分离方法,以提升滤波效果。通过油田实测数据处理,证明了该方法的有效性。

1 线性预测滤波法

反射声波测井数据中由2类波形构成。一种是井中高幅度直达波,从声源沿井孔传至各个接收器。单极反射声波测井中的纵波、横波、斯通利波以及偶极反射声波测井中的弯曲波均属于井中直达波。这类波的走时特征确定,它们的速度(或者慢度)可以根据阵列信号处理得到。井中直达波的走时曲线在接收器阵列上表现为线性,斜率即为井中直达波的速度;而井外反射波,由于难以确定井外反射体的位置,因此,该类反射波的到时与走时特征未知,反射波的走时曲线在接收器阵列上表现为双曲线形态。一般情况下,井中直达波与井外反射波的走时特征拥有明显的差别,线性预测滤波法正是基于这样一种前提下提出的。

线性预测滤波法[12]的思路:利用已知的走时特征,估计得到井中直达波,从测井波形中减去井中直达波,剩余的波形即为反射波。

采用参数估计法来估计得到直达波[14,16]。该方法中,将直达波看作是以已知速度传播的模式,Al(ω)exp(-iωz/vl)(l=1,…,L),其中,ω为角频率,L是井中直达波的振型个数(比如,对单极数据,L=3,包括P波、S波、斯通利波;偶极数据L=1,对应弯曲波),z代表第n个接收器的所处位置,vl是第l种直达波振型对应的传播速度。假设该直达波从第n个接收器传至第m个接收器,在数学上可以把这一过程表示为Al(ω)exp(-iω(m-n)d/vl),其中,d是相邻2个接收器之间的距离。指数m与n之间大小关系的不同,对应该直达波的传播方向不同,m大于n意味着该直达波向前传播;反之,m小于n意味着该直达波向后传播。综合考虑所有的L种直达波,可以表示为

(1)

El=exp(-iωd/vl)

(2)

式中,N是阵列接收器的个数;An(ω)(n=1,…,L)是某一接收器n接收到的井中直达波的L个振型,为待求量;Wn(ω)(n=1,…,N)是测井仪器接收到的实际波形。式(1)的矩阵形式为

E·A=W

(3)

E是1个N×L的复矩阵,其中包含对应于直达波沿井轴传播的指数。对于L

(4)

式中,～表示复共轭;T表示转置。通过方程(4),得出了接收器位置n处直达波各个振型所对应的频谱。将估算得到的各个阵型的频谱进行相加,则得到了对应该接收器位置n处的直达波的完整频谱。

至此,已经通过参数估计的方法得到了第n个接收器对应波形中直达波的频谱,将其从第n个接收器的全波形的频谱中减去,得到的剩余的第n个接收器信号对应的频谱

(5)

将式(5)变换回时间域,得到了对应于第n个接收器的、时域的、去除了直达波的剩余信号。对n=1,…,N重复,则得到完整的阵列反射信号。

该剩余信号中,主要是由井外反射信号组成,还存在一系列数据噪声以及由于井壁不规则所引起的散射波。通常情况下,在经过线性预测法滤波后的数据中,残余的模式波以及其他噪声与反射信号会处于同一数量级,在成像前还需进行倾角叠加[12]、fk滤波[17]、中值滤波等操作,以进一步增强反射信号的强度。

对应偶极子的情况,线性预测法中通常取L=1,即认为只有弯曲波一种振型,vl=1为对应地层的横波速度。不过,值得注意的是,弯曲波是一种频散波,不同频率的弯曲波对应的传播速度是不同的,在低频时弯曲波的传播速度趋向于横波速度,高频时趋向于流体速度[18-19]。因此,偶极反射声波测井中,如果只取L=1,取vl=1为横波速度,只能滤掉对应的弯曲波中以横波速度传播的部分,而以低于横波速度传播的弯曲波部分将留在剩余信号中,对井外反射波的提取形成干扰。

2 基于慢度—时间相干法的线性预测滤波

假设测井仪器有8个接收器(即N=8),以n=8时为例,对式(1)进行进一步讨论,该情况下对应的式(1)为

(6)

式中,A1(ω)、A2(ω)…AL(ω)分别为第8个接收器接收到的井中直达波中不同振型的频谱。如果知道准确的vl=1、vl=2…vl=L,则E1-7(ω)A1(ω)+E2-7(ω)A2(ω)+…+EL-7(ω)AL(ω)意味着第8个接收器的直达波中各个振型在频域上传至第1个接收器并进行叠加。不考虑微弱反射信号的影响,叠加结果理论上应与第1个接收器接收到的实测数据频谱W1(ω)相等。重复该过程,将第8接收器的直达波通过传播矩阵E与各个接收器的实测数据频谱联系在一起。最终,通过实测波形的频谱W,以及由准确传播速度组成的传播矩阵E,共同反演得到准确的第8个接收器接收到的直达波频谱。如果使用了过高和过低的传播速度,则意味着传播矩阵与实际情况不符,进而导致反演结果中估计得到的直达波频谱的不准确。

另一种阵列信号处理的方法,慢度—时间相干法(STC),是使用二维网格(一维为时间,另一维为慢度)搜索法,找出阵列波形的相关函数极值位置对应的波的到时与慢度[20]。是目前广泛使用的较为成熟的计算慢度的方法。

该方法在计算相关函数时,会将每一接收器的波形数据在时域上往回传播,移动到第一接收器的位置。然后将N个数据点(第一个换能器接收的数据加上N-1个经时移处理的波形数据)相加求和。该时域的计算过程,实际上与线性预测法中E·A的频域处理拥有同样的原理。

对比线性预测滤波法与慢度—时间相干法,线性预测滤波法是根据实测数据W与传播矩阵E,估算井中直达波A,且传播矩阵E的准确与否直接影响最终滤波效果。慢度—时间相干法(STC)是根据实测数据W,估算井中直达波各个振型的速度,实际上就是估算传播矩阵E。线性预测滤波法是频域的处理,慢度—时间相干法是时域的处理,但都基于同样的理论背景。

因此,本文提出一种数据处理流程:在处理偶极子反射声波测井数据时,先不限定井中直达波只有弯曲波一个振型,而是利用STC方法,从实测数据W中求得实际的慢度范围(由于频散,弯曲波对应的慢度应该是一个范围),即求得实测数据W对应的准确传播矩阵E;在线性预测滤波的过程中,不仅考虑以横波速度传播的弯曲波,而且考虑以整个慢度范围传播的各个振型,利用实测数据W与通过STC方法求得的传播矩阵E,估算井中直达波A。

3 实测数据处理

为进一步阐明上述方法的有效性,采用Xmac-F1正交偶极子阵列声波测井仪器,对该仪器在某油田测量得到的偶极反射波测井数据进行处理。该仪器源距为3.200 4 m,有8个四分量接收器(各个分量相互正交),分别记录XX、YY、XY、YX的分量信号,8个接收器之间间距为0.152 4 m。数据接收的深度范围为×720～×460 m(用×720与×460来代替测井仪器真实所处的深度),共有1 707个深度点。波形的采样间隔为36 μs,共有400个时间点,对应记录长度为14.4 ms。

首先取测井仪器在某一深度(×582 m)测得的全波波形进行处理。对该全波波形进行带通滤波后如图1(a)所示。图1(a)中有幅度较强的弯曲波,经过速度分析,该深度点的横波速度约为3 000 m/s,根据该速度绘出时距曲线,与弯曲波到时符合较好。

对该波形进行STC处理,得到如图1(b)所示的慢度—时间相关图。从图1(b)可以看出,大约有2个极值区域,一个对应S1=330 μs/m(对应速度3 000 m/s),另一个对应大约S2=430 μs/m(对应速度2 300 m/s)。S1对应横波速度,S2则是由于频散造成的弯曲波中低于横波速度的部分。

图1 ×582 m深度全波波形预处理结果

记录速度的文件显示,该深度点提取的横波慢度为337.211 μs/m,对应速度约为2 970 m/s,即对应图1(b)中的S1。实际处理中,是在滤波过程中从速度文件中读取该深度点对应的横波速度,构成式(3)中的传播矩阵E,进行线性预测滤波。

图2 传统线性预测滤波法处理结果

按常规偶极反射声波测井的处理方法,取L=1(即假设井中直达波中只有一种振型,以横波速度传播的弯曲波),速度取横波速度v=3 000 m/s,对波形进行线性预测滤波。滤波结果见图2(a),从波形图中可以看出,相对于图1(a),直达波波形有所压制,但滤波效果差。采用横波速度滤波,红色时距曲线标定的以横波速度传播的部分被滤掉了,但是在红色曲线之后到达的波形,未被有效压制。

图2(b)的相干图也印证了该现象,在相干图中,慢度S1对应的极值区域被有效压制了,而慢度S2对应的极值区域对应未被滤掉的低速传播的弯曲波部分。

由此可以推测,只有将所有的极值区域从慢度—时间相干图中滤掉,才能对应地从时域波形中将频散的弯曲波完全压制。

因此,对应图1(b)中整个极值区域的上下边界,上边界大约对应慢度为250 μs/m(速度4 000 m/s),下边界对应慢度为500 μs/m(速度2 000 m/s),将其均分为5等分,即取L=5(假设井中直达波有5种不同速度传播的振型),其中v1=4 000 m/s,v2=3 500 m/s,v3=3 000 m/s,v4=2 500 m/s,v5=2 000 m/s。构建传播矩阵E,进行线性滤波。滤波结果如图3(a)中蓝色波形所示,相对于图1(a)的原始波形与图2(a)中只采用横波速度滤波的波形,弯曲波得到了较大的压制。

此外,将采用L=4(其中v1=4 000 m/s,v2=3 300 m/s,v3=2 700 m/s,v4=2 000 m/s)滤波的结果用红色波形表示,采用L=3(其中v1=4 000 m/s,v2=3 000 m/s,v3=2 000 m/s)滤波的结果用黑色波形表示,同时画在图3(a)中。

图3 结合了慢度—时间相干法与线性预测滤波法的波场分离方法处理结果

观察L=5的蓝色波形、L=4的红色波形和L=3的黑色波形,滤波之后井中弯曲波均受到较大压制,且残余弯曲波幅度相近。这说明线性滤波中速度的选取以及直达波振型个数L的选取对滤波效果的影响不大,该方法具有较高的稳定性。

同样,图3(b)的时间—慢度相干图中,所有的极值区域均从慢度—时间相干图中被滤掉,对应的时域波形中频散的弯曲波亦被有效压制。

线性滤波之后[见图3(a)],仍然有残余噪声存在,其中包括由于井壁不规则所产生的散射波、实际波形采集时不可避免的数据噪声。在偏移成像之前,还需要对线性滤波后的波形沿时距曲线进行叠加增强处理以及fk滤波,以进一步提升反射信号信噪比。

对整个深度范围,×720～×460 m段重复上述数据处理流程(见图4)。图4(a)为YY分量经过带通滤波后的波形,图中弯曲波幅度较大,井外反射波微弱,难以识别。图4(b)与4(c)为经过线性预测滤波、沿时距曲线叠加增强、fk滤波后的波形图。2幅图的区别:图4(b)中,仅利用横波速度作为滤波速度,取L=1(其中v1=3 000 m/s),进行线性预测滤波;图4(c)中,根据慢度—时间相干图提供的速度范围,采用L=4(其中v1=4 000 m/s,v2=3 300 m/s,v3=2 700 m/s,v4=2 000 m/s)进行线性预测滤波。

对比图4(b)和4(c)可知,考虑了频散效应的线性预测滤波法,相对于传统的滤波方法,井中弯曲波被进一步的压制,井外反射信号更加清晰。尤其是对于微弱的反射信号2,图4(b)中由于与残余弯曲波幅度差别较大而显得较微弱;而图4(c)中,则在×660深度附近清晰可见下行反射信号。

图4 实际测井资料处理结果

4 结 论

(1) 线性预测滤波法是偶极反射声波测井中常用的滤波方法,主要根据井中直达波与井外反射波的走时特征不同进行滤波。偶极反射声波测井中,弯曲波为主要的井中直达波,其到时与地层横波速度相近,实际处理中常常设定弯曲波速度即为横波速度,进行线性滤波。由于没有考虑到弯曲波的频散特征,在滤波中仅仅将以横波速度传播的部分滤掉,而以低于横波速度传播的部分残留在处理后的波形中,对反射波的提取形成干扰。

(2) 慢度—时间相干法是在时域对阵列信号进行处理的方法,可以识别出阵列信号中存在的振型数以及各个振型的传播速度,即某一振型中不同部分的传播速度。

(3) 处理偶极子反射声波测井数据时,先不限定井中直达波仅有弯曲波一个振型,而是利用慢度—时间相干法,从实测数据中求得实际的慢度范围(由于频散,弯曲波对应的慢度应该是一个范围),进而求得实测数据对应的传播矩阵。在线性预测滤波的过程中,不仅考虑以横波速度传播的弯曲波,而应该考虑以整个慢度范围传播的各个振型。通过对实测数据的处理发现,考虑了频散效应、结合了慢度—时间相干法的线性预测滤波法,较常规滤波方法而言,井中弯曲波被进一步压制,井外微弱的反射波被有效提取。并且,振型个数以及对应速度的选取,对滤波效果的影响较小。这意味着该方法具有较好的稳定性,适于用在实际处理中。

参考文献:

[1] 李宁, 肖承文, 伍丽红, 等. 复杂碳酸盐岩储层测井评价: 中国的创新与发展 [J]. 测井技术, 2014, 38(1): 1-10.

[2] ESMERSOY C, CHANG C, KANE M, et al. Acoustic Imaging of Reservoir Structure from a Horizontal Well [J]. The Leading Edge, 1998, 17(7): 940-946.

[3] COATES R, KANE M, CHANG C, et al. Single-well Sonic Imaging: High-definition Reservoir Cross-sections from Horizontal Wells [C]∥The SPE/CIM International Conference on Horizontal Well Technology, F, 2000. Society of Petroleum Engineers.

[4] 薛梅, 楚泽涵, 边环玲, 等. 远探测声波反射波测井方法研究及声系设计方案 [J]. 测井技术, 2002, 26(1): 35-39.

[5] 车小花, 乔文孝, 阎相祯. 相控线阵声波辐射技术在反射声波测井中的应用探讨 [J]. 测井技术, 2004, 28(2): 108-111.

[6] HORNBY B E. Imaging of Near-borehole Structure Using Full-waveform Sonic Data [J]. Geophysics, 1989, 54(6): 747-757.

[7] TANG X M. Imaging Near-borehole Structure Using Directional Acoustic-wave Measurement [J]. Geophysics, 2004, 69(6): 1378-1386.

[8] TANG X M, PATTERSON D J. Single-well S-wave Imaging Using Multicomponent Dipole Acoustic-log Data [J]. Geophysics, 2009, 74(6): WCA211-WCA223.

[9] 李超, 岳文正, 金行林, 等. 声反射成像测井数据处理研究进展 [J]. 测井技术, 2013, 37(1): 13-20.

[10] GONG H, CHEN H, HE X, et al. Eliminating the Azimuth Ambiguity in Single-well Imaging Using 3C Sonic Data [J]. Geophysics, 2015, 80(1): A13-A17.

[11] ZHANG Y D, HU H. A Technique to Eliminate the Azimuth Ambiguity in Single-well Imaging [J]. Geophysics, 2014, 79(6): D409-D416.

[12] TANG X, ZHENG Y, PATTERSON D. Processing Array Acoustic-logging Data to Image Near-borehole Geologic Structures [J]. Geophysics, 2007, 72(2): E87-E97.

[13] BING W, GUO T, HUA W, et al. Extracting Near-borehole P and S Reflections from Array Sonic Logging Data [J]. Journal of Geophysics and Engineering, 2011, 8(2): 308.

[14] TANG X M. Predictive Processing of Array Acoustic Waveform Data [J]. Geophysics, 1997, 62(6): 1710-1714.

[15] HIRABAYASHI N, LEANEY W S, HALDORSEN J B. Wavefield Separation for Borehole Acoustic Reflection Surveys Using Parametric Inversion [C]∥ The 2008 SEG Annual Meeting, F, 2008. Society of Exploration Geophysicists.

[16] 唐晓明, 郑传汉, 赵晓敏. 定量测井声学 [M]. 北京: 石油工业出版社, 2004.

[17] LI Y, ZHOU R, TANG X, et al. Single-well Imaging with Acoustic Reflection Survey at Mounds, Oklahoma, USA [C]∥The 64th EAGE Conference & Exhibition, F, 2002.

[18] 张海澜, 王秀明, 张碧星. 井孔中的声场和波 [M]. 北京: 科学出版社, 2004.

[19] KURKJIAN A L, CHANG S K. Acoustic Multipole Sources in Fluid-filled Boreholes [J]. Geophysics, 1986, 51(1): 148-163.

[20] KIMBALL C V, MARZETTA T L. Semblance Processing of Borehole Acoustic Array Data [J]. Geophysics, 1984, 49(3): 274-281.