基于有限差分方法的溶洞地层井孔声波数值模拟

徐方慧,王祝文,武焕平

(吉林大学地球探测科学与技术学院,吉林长春130026)

碳酸盐岩和火成岩油气藏中溶洞发育,并且溶洞比裂缝具有更强的非均质性,使得利用常规测井资料评价溶洞型储层非常困难[1-3]。地层溶洞的存在使声波测井的波形发生变化[1,4],因此掌握溶洞对井孔声波的影响规律将对溶洞的检测和评价起到重要作用。目前国内外关于裂缝对井中声波传播影响方面的研究成果见诸于许多文献[5-8],但对溶洞性储层声波测井响应的数值模拟相对不多。BURNS等[9]指出,地震波在地层传播时遇到裂缝、溶洞等非均匀体会产生散射现象。

偶极子声源向井外辐射的能量以SH横波为主[10-11],其中S-S散射波为主要成分[12],黑创等[13]利用偶极子源研究了非均匀地层偶极声波测井的散射效应,并分析了散射波的衰减特性。李丹等[14]针对井旁不同尺度溶洞的偶极反射声波测井响应进行了数值模拟,得到了偶极横波的回波信号,探究了地层横波波长与井旁溶洞尺度之间的关系。赵军等[15]认为溶洞的充填度与纵横波时差之间具有很好的对应关系。根据充填情况的不同,溶洞可以分为未充填溶洞、砂泥充填溶洞、角砾充填溶洞和方解石充填溶洞4种类型[16]。溶洞的存在使地层的纵波慢度增大,密度测井值降低,同时深、浅双侧向电阻率值普遍较低[17]。利用声波测井和电成像测井能全面评价地层溶洞,斯通利波受到溶洞的影响能量会减弱[1,4]。

前人研究成果表明,利用声波测井评价储层中的溶洞是可行的。但其成果多基于实际测井数据,定性给出了溶洞对纵波等常规测井曲线的影响,并没有给出溶洞的尺寸、位置和数量对井孔纵波、横波以及散射波的影响规律。我们利用单极子声源的弹性波动方程有限差分方法对储层中不同半径、位置和数量的含水溶洞的全波列测井响应进行了数值模拟,分析了溶洞对纵、横波首波的影响,同时研究了由溶洞形成的下行散射波的特征。

在利用弹性波动方程有限差分程序模拟井孔内、外波场的过程中,矩阵计算占据绝大部分计算量;与此同时需要划分细小的网格来满足小溶洞的计算需求,大量增加了数值模拟的时间成本。因此,我们采用统一计算设备架构(Compute Unified Device Architecture,CUDA)平台的GPU并行计算[18-20]技术来提高计算效率。将程序的CPU串行计算模式转换为GPU并行计算模式,以减少计算时间。

1 方法简介

1.1 弹性波有限差分方程

为了研究方便,数值模拟只考虑二维情况下x和z方向上的波动方程。选用二维直角坐标系,采用交错网格有限差分法[21-23]模拟弹性波在发育有溶洞的均匀各向同性弹性介质模型中的传播。介质的物性参数和应力、速度分量在交错网格中的位置如图1所示。

图1 场量和介质参数在交错网格中的位置示意

图1中空心圆表示剪切应力;实心圆表示正应力和地层参数;空心方形表示z方向上的速度分量vz;实心方形表示x方向上的速度分量vx。括号中的参数i和j分别表示x和z方向上的网格点数。在x、z方向上网格步长分别为Δx和Δz,Δx=Δz。二维直角坐标系下,给出时间上二阶、空间上四阶的弹性波在均匀介质中传播的速度-应力有限差分方程为:

(1)

(2)

(3)

(4)

(5)

σif(i,j)=0.5[f(i+1,j)+f(i,j)]

(6)

σjf(i,j)=0.5[f(i,j)+f(i,j+1)]

(7)

(8)

为了减小在计算过程中的数值频散,纵波速度和声源频率需要满足以下不等式:

(9)

式中:Δmax为模型中最大网格步长;vmin为模型中最小的纵波速度;fmax为声源的最高频率。除此之外弹性波有限差分方程需满足稳定性条件:

(10)

式中:vmax为模型中最大的纵波速度。选取中心频率为10kHz的雷克子波作为声源。

溶洞一般指地层中直径大于2mm的近似球体的空隙空间。溶洞不一定与井孔相交而可能通过散射波影响井孔声波。为了满足溶洞的计算需求,网格差分步长Δx=Δz=2mm。根据稳定性条件可知,Δt=0.2μs。地层中裂缝的宽度一般小于1mm,需要采用变网格差分方法进行模拟。

1.2 模型边界处理

考虑到模型边界处可能会形成人工虚假反射波,故在模型中加入了基于阻抗匹配的完全匹配层分裂吸收边界(PML)[24]。PML的吸收效果主要受吸收层数和衰减因子的影响。本文采用的是COLLINO等[25]提出的一种较为常用的衰减因子:

(11)

(12)

式中:x是计算点到内层边界的距离;H表示吸收边界的宽度;d0是与Rc有关的衰减参数,Rc为理论反射系数,通常取0.001。

1.3 CUDA并行实现

为了满足溶洞的模拟精度,文中选择的网格差分步长较小,因而大大增加了波场模拟的时间成本,我们建立了弹性波动方程有限差分的并行计算模型,算法实现流程如图2所示。由弹性波动方程有限差分公式可知,当前速度场每个网格的变化量由前一时间节点应力场相关网格的数值计算得出,而当前应力场变化量则由计算后的速度场相关网格的数值计算得出。可见,在弹性波动方程有限差分计算模型中,待计算场量仅与已获得的计算结果相关,而待计算场量各个网格之间在当前计算步骤中没有关联。因此,并行计算方式可以简单、高效地实现波动方程的有限差分,而无需考虑计算过程中不同计算结果间的同步问题。在计算模型初始化时,通过CPU分别为速度场vxx、vzz、vxz和应力场τxx、τzz、τxz以及介质参数在GPU内存中分配二维的存储空间,对应力场、速度场和声源设计各自的核函数,通过CPU控制GPU分别完成速度场和应力场的分步计算。利用CUDA并行计算来提高模拟程序的运算速度。模拟实验发现,对于本文的均匀介质模型,常规程序的计算时间为8808.459s,并行程序计算时间为584.776s。CUDA并行计算有效提高了计算效率。

图2 CUDA并行计算流程

1.4 误差分析

为了验证本文有限差分方法计算结果的正确性,利用CUDA并行计算有限差分方法和实轴积分法(RAI)模拟了均匀弹性地层中声波的传播[26]。地层参数如表1所示。声源的中心频率为10kHz。图3对比了两种方法的计算结果。由图3可以看出,两种方法的波形匹配较好,说明了二维直角坐标系下的有限差分方法是有效的。这两种方法波形的差异可能是数据离散和归一化引起的。

图3 有限差分法和实轴积分法计算结果对比

2 数值模拟分析

2.1 建立计算模型

本文设计了如图4所示的试算模型。模型尺寸为4m×2m。井旁有一圆形溶洞,井孔、溶洞中充满了水,井孔外为均匀弹性介质,溶洞半径为R,溶洞距井壁的距离为L,井孔的半径a=10cm。声源和接收器均放置在井孔中。声源距溶洞中心的垂直距离s为2m。各介质的声学参数如表1所示。

图4 溶洞地层的井孔模型示意

表1 地层和井孔参数

2.2 溶洞地层的波场快照

图5a和图5b分别给出了溶洞地层中弹性波在0.90ms和1.44ms的波场快照。单极子声源向地层辐射纵波和横波。当溶洞直径小于0.25倍横波波长时,溶洞可近似为点散射体,点散射体只能形成一次散射波。直径大于0.25倍横波波长的溶洞自身能形成多次散射波,散射波的叠加会使散射波前面畸变[27]。声波沿着地层传播到模型的边界PML区域时没有形成人工边界的反射,说明加入的PML吸收边界是有效的。在0.90ms时,只有纵波与溶洞相遇。由于纵波产生的散射波较弱,在图5a中不能被观察到,只能看到声源激发的原始波场。在1.44ms时,地层横波与溶洞相遇,除了声源激发的原始波场外,在图5b中能清楚看到横波产生的散射波—S-S波。值得注意的是,图5b中横波穿过溶洞后,幅值明显降低(横波颜色变淡)。

图5 溶洞地层弹性波的波场快照a 0.90ms; b 1.44ms

2.3 地层溶洞对井孔纵波首波的影响

2.3.1 溶洞的位置和半径对纵波首波的影响

图6展示了位于溶洞上方0.4m(源距为2.4m)的接收器所记录的纵波首波波形。其中,R表示溶洞的半径,L表示溶洞距井壁的距离。由图6可知,当溶洞半径较小时(R<20mm),纵波首波幅值变化很小。当溶洞半径较大时(R>20mm),存在一个距离Lp,当溶洞的井旁距离小于Lp时,纵波首波的衰减与Lp呈正相关并在Lp处最大;当溶洞的井旁距离大于Lp时,纵波首波的衰减与Lp呈负相关。经过观察,文中模型的Lp为16cm。当溶洞的井旁距离大于24cm时,井旁距离的变化对纵波首波影响很小,这时溶洞地层与均匀弹性地层的纵波首波比较接近,说明单极子声波测井存在一个最大探测深度Lmax。在本文的模型中,Lmax约为24cm。当溶洞的半径较小(R<30mm)时,溶洞位置的变化对纵波首波的影响很小。

图6 溶洞位置对纵波首波的影响

图7给出了不同半径的溶洞对纵波首波的影响。由图7可知,纵波首波的衰减与溶洞半径呈正相关。当溶洞距井壁的距离<24cm时,纵波的初至时间与溶洞半径呈正相关。当溶洞距井壁的距离为16cm时,纵波首波的衰减最大。

图7 溶洞半径对纵波首波的影响

溶洞上方的接收器记录到的纵波是溶洞引起的散射波与井壁滑行纵波的叠加。当溶洞半径较大,且距井壁的距离不超过声源探测的最大距离时,不同位置的溶洞形成的散射波与井壁滑行纵波的相位、幅值可能不同,故纵波首波的衰减也不同。由于溶洞中充满水,声波在流体中的传播速度小于在地层中的传播速度,故在地层含有溶洞的情况下,纵波首波的初至时间后延。溶洞的半径越大,纵波在流体中的传播时间越长,纵波初至时间越晚。

2.3.2 溶洞的数量对纵波首波的影响

溶洞的数量也会影响声波在井中的传播,我们模拟了等间隔径向排列的多个溶洞对井孔声波的影响,最近的溶洞距井壁1cm,溶洞的间隔为1cm。图8给出了溶洞的数量(N)对纵波首波的影响结果。

图8 溶洞的数量对纵波首波的影响

1) 溶洞的半径为20mm时,纵波首波的衰减与溶洞数量呈正相关,纵波首波的初至时间与溶洞数量呈正相关。

2) 溶洞的半径为40mm时,当溶洞数量≤3时,纵波首波的衰减、初至时间均与溶洞数量呈正相关;当溶洞数量>3时,纵波首波的衰减、初至时间几乎保持不变。

由此可知,若溶洞的井旁距离大于单极子测井的最大探测距离,溶洞已经较难对井孔声波产生影响。当溶洞的半径R为20mm时,横向排列、间隔为1cm的5个溶洞的最大井旁距离为22cm,故此时纵波首波的变化与溶洞数量呈正相关。当溶洞的半径为40mm时,由于第4个和第5个溶洞的井旁距离均大于24cm,纵波首波的衰减几乎不再受第4个和第5个溶洞的影响,此时纵波的变化与前3个溶洞的纵波情况非常接近。

2.4 地层溶洞对井孔横波首波的影响

2.4.1 溶洞的位置和半径对横波首波的影响

图9给出了不同位置的溶洞对横波首波的影响。由图9可知,同样存在一个井旁距离Ls,当溶洞的井旁距离小于Ls时,横波首波衰减与其呈正相关并在Ls处衰减达到最大;当井旁距离大于Ls时,横波首波的衰减较弱。当溶洞的井旁距离>24cm时,井旁距离的变化对横波首波影响很小;当溶洞的半径较小(R<20mm)时,溶洞位置的变化对横波首波的影响很小。对比纵波可知,横波对于溶洞井旁距离的变化更敏感。

图9 溶洞的位置对横波首波的影响

图10给出了不同半径的溶洞对横波首波的影响。由图10可知,溶洞使横波首波明显衰减,横波首波的衰减与溶洞的半径呈正相关。与纵波对比可知,横波对地层溶洞半径的变化更敏感。当L=8cm时,横波的衰减最大。

图10 溶洞的半径对横波首波的影响

地层溶洞内充满流体,横波不能在流体中传播。横波遇到含水溶洞时,溶洞中的流体会抑制横波的传播,一部分横波能量会转化为流体的压缩能量[28]。这部分流体压缩能量可能又转换成其它模式波的能量,继续沿地层传播,最终被接收器接收。因此在溶洞上方的接收器记录到的横波能量会有部分的衰减。溶洞的半径越大,横波的衰减幅度越大。由于单极子声源的探测距离限制,当溶洞的井旁距离超过单极子声波测井的最大探测距离时,溶洞半径对横波首波的影响将会变得很小。

2.4.2 溶洞的数量对横波首波的影响

图11给出了不同的溶洞数量对横波首波的影响。溶洞数量对横波首波的影响与纵波类似。

图11 溶洞数量对横波首波的影响

1) 溶洞的半径为20mm时,横波首波衰减非常明显,并且与溶洞数量呈正相关。

2) 溶洞的半径为40mm时,当溶洞数量≤3时,横波首波衰减与溶洞数量呈正相关;当溶洞数量>3时,横波首波衰减的程度几乎保持不变。

2.4.3 溶洞与横波首波衰减系数的关系

本文将横波首波在溶洞作用下的衰减系数A定义为:

(13)

式中:A1和A2分别为均匀地层和溶洞地层在同一源距横波首波的振幅。图12a为溶洞半径与横波首波衰减系数的关系,其中虚线代表拟合的线性方程。由图12a可知,横波首波衰减系数与溶洞半径存在正线性相关。特别的,当溶洞的井旁距离为4cm时,不仅横波首波衰减系数最大,线性方程的斜率也最大。图12b为溶洞数量与横波首波衰减系数的关系。溶洞数量与横波首波衰减系数之间更像是凸函数,随着溶洞数量的增加,衰减系数增加得越来越缓慢。这是因为声波探测距离有限,一般井旁距离大于24cm的溶洞很难再对井孔声波产生影响,故当溶洞在地层中径向等间距排列时,横波首波衰减系数与有限个溶洞数量呈正相关。图12a和图12b中接收器的源距为2.4m。图12c为跨过溶洞的不同源距与横波首波衰减系数的关系。溶洞的源距为2.0m。当接收器位于溶洞下方时,横波几乎不发生衰减,此时衰减系数接近0。值得注意的是,位于溶洞处的接收器接收到的横波衰减仍然很小(源距为2.0m时)。当接收器逐渐跨越溶洞时,横波在地层中传播会遇到含水溶洞而形成转换纵波和转换横波,横波能量的衰减逐渐增强。靠近溶洞上方的接收器所记录到的波形中转换纵波和转换横波叠加在一起,故横波的波幅较大,首波衰减系数较小。随着源距的增加,转换纵波和转换横波逐渐分离,一部分横波能量就转化为纵波的能量,明显削弱了横波首波幅值,首波衰减系数会增加。随着源距的不断增加,溶洞对横波的影响越来越小,横波能量衰减也越来越小,故衰减系数会逐渐减小并趋于稳定。从图12c中可以看出,当接收器与溶洞距离为0.4m时,转换横波已经与其它模式波分离,当接收器与溶洞距离>0.9m时,衰减系数趋于稳定。

图12 溶洞的半径(a)、数量(b)和源距(c)对横波首波衰减系数的影响

2.5 地层溶洞对散射波的影响

图13a为声波在溶洞地层中的波形。溶洞下方出现了明显的散射波。选取源距为1.6m处的井孔声波(溶洞下方0.4m),我们将均匀地层的波形作为直达波。从图13a中减去直达波,分离出散射波,如图13b和图13c所示。由图13b和图13c可知,分离出的波形主要以散射波为主。图13b中溶洞的半径为10mm,幅值较大的波峰以S-S、P-S散射横波为主,主频约为7.5kHz,幅值较小的波峰由其余散射波组成。图13c中溶洞半径为30mm,自身形成多次复杂的散射波场,其中主频为5～7kHz的大波峰以横波散射波(S-S波等)为主,主频为9kHz的小波峰以纵波散射波(P-P波等)为主。利用傅里叶变换得到散射波的频率谱,如图13d和图13e所示。将每个频率对应数值的平方求和,得到散射波的能量。

图13 溶洞地层中的声波波形(a)、散射波波形(b,c)以及散射波频率谱(d,e)

图14a为散射波能量与溶洞位置的关系曲线。其中半径为6～20mm的溶洞均可视为点散射体,此时散射波能量与井旁距离呈负相关;半径为30～60mm的溶洞可形成多次散射,当井旁距离较大时,散射波能量与井旁距离呈负相关。当溶洞的井旁距离超过单极子声波的最大探测距离时,任何尺寸的溶洞形成的散射波都可忽略。图14b为下行散射波能量与溶洞半径的关系曲线。当溶洞半径可视为点散射体时,形成的下行散射波能量与溶洞半径呈正相关;当溶洞半径较大可形成自身散射时,下行散射波能量降低,井旁距离越小,降低越明显。随着溶洞半径增大,下行散射波能量变化不是很剧烈。图14c展示了溶洞数量对下行散射波能量的影响,当溶洞在地层中径向等间距排列时,下行散射波能量与有限个溶洞数量呈正相关。由于单极子声波探测距离的限制,半径为30mm和40mm的溶洞数量达到2时,下行散射波的能量就不再增加。

图14 下行散射波能量与溶洞位置(a)、半径(b)和数量(c)的关系曲线

3 结论

通过对溶洞地层的单极子声波有限差分方法的数值模拟,研究了溶洞对井孔声波的初至和能量的影响规律,得到以下结论:

1) 单极子声源测井的探测距离是有限的。当接收器源距为2.4m时,若溶洞的距离大于单极子声波的探测距离,很难再对井孔声波产生影响。文中模型声源的最大探测距离约24cm。

2) 由于散射波的叠加,纵波的幅度会发生变化。横波的传播机制决定了横波能量总是有所衰减。相对于纵波,横波对溶洞更敏感,利用横波能更有效地识别溶洞。随着源距的增加,横波的衰减系数会先增大后减小。在单极子声波的最大有效探测距离内,溶洞的数量越多,对声波的影响越明显。纵波和横波对溶洞的距离没有很好的响应规律。值得注意的是,当溶洞离井孔的距离为某个值时,纵波和横波的幅度变化会达到最大。

3) 可视为点散射体的小尺寸溶洞和可形成自身散射的大尺寸溶洞形成的下行散射波的规律明显不同。值得注意的是,离井壁很近的大尺寸溶洞的散射波与溶洞半径或溶洞的井旁距离没有明显的相关性。

溶洞对纵波、横波的影响与裂缝有些类似,上述结论能够为利用声波测井评价溶洞地层提供一些理论依据。由于文中的溶洞位于井旁没有与井孔相交,故不能直接对井孔斯通利波产生影响。本文的理论结果虽然能在一定程度上反映出地层溶洞对井孔声波的影响,但也存在不足之处。本文模拟的模型是二维的理想简化模型,但实际地层中的溶洞无论是数量、形状和填充度等方面都十分复杂,今后需要利用三维复杂地质模型更加系统地研究声波测井的响应特征。