地球物理学报  2016, Vol. 59 Issue (5): 1776-1789   PDF    
裂缝诱导双相HTI介质地震波场错格伪谱法模拟与波场特征分析
刘财1,2, 迟唤昭2, 高炜1, 鹿琪1, 兰慧田3    
1. 吉林大学地球探测科学与技术学院, 长春 130026;
2. 吉林大学地球科学学院, 长春 130061;
3. 大庆油田有限责任公司勘探开发研究院, 黑龙江大庆 163712
摘要: 裂缝诱导的双相具有水平对称轴的横向各向同性(HTI)介质模型是由一组平行排列的垂直裂缝嵌入到统计各向同性的流体饱和多孔隙岩石中而组成的,它综合考虑了裂缝型储层岩石的各向异性和孔隙性.高精度的地震波场数值模拟技术是研究该介质中地震波传播规律的主要方法.本文结合错格伪谱法和时间分裂法,求解描述该介质中地震波传播的一阶速度-应力方程.模拟了单层和双层模型中的地震波场,并对其进行了特征分析.研究结果表明:错格伪谱法能有效消除标准网格伪谱法波场模拟结果中出现的数值伪影现象,与时间分裂法结合能够获得稳定的、高精度的模拟结果;裂缝诱导双相HTI介质中的地震波场兼具裂缝各向异性介质和双相介质中传播的地震波的波场特征.
关键词: 裂缝诱导     双相HTI介质     错格伪谱法     时间分裂法     波场特征    
Seismic wavefield simulation and feature analysis for a fracture-induced two-phase HTI medium based on the staggered-grid pseudo-spectral method
LIU Cai1,2, CHI Huan-Zhao2, GAO Wei1, LU Qi1, LAN Hui-Tian3    
1. College of Geo-exploration Science and Technology, Jilin University, Changchun 130026, China;
2. College of Earth Sciences, Jilin University, Changchun 130061, China;
3. Exploration and Development Research Institute of Daqing Oilfield Company Ltd., Heilongjiang Daqing 163712, China
Abstract: A model for a fracture-induced horizontal transverse isotropic (HTI) double-porosity medium has been established by embedding a set of vertical fractures into the isotropic porous background rock in previous work, which considers the fractured reservoir rock anisotropy and porosity. High-precision seismic wavefield simulation is the main technique for the investigation of waves propagating in this media. In this paper, seismic wavefield is simulated using a single-layer model and two-layer model by solving the one-order velocity-stress equations describing seismic wave propagation in this medium by combining staggered-grid pseudo-spectral method and time-splitting method. Then the wavefield features are analyzed. The results of numerical simulation indicate that the artifacts in the wavefield simulation results by the standard-grid pseudo-spectral method can be removed by the staggered-grid pseudo-spectral method, and the stable and high-precision simulations results can be obtained by combining it with the time-splitting method. The seismic wavefield in the fractured-induced two-phase HTI medium has the characteristics of wavefield of seismic waves propagating in a fracture anisotropic medium and two-phase medium.
Key words: Fracture-induced     Two-phase HTI media     Staggered-grid pseudo-spectral method     Time-splitting method     Wavefield characteristics    
1 引言

天然裂缝型油气藏作为一类重要的油气藏类型已经引起了越来越多的关注.基于地震各向异性理论发展的相关技术是裂缝性储层识别和预测的一类重要手段,而地震各向异性与裂缝型储层参数的联系是通过等效介质理论模型来建立的,因而,等效介质模型间接地对裂缝型油气藏的勘探起到重要作用.理论模型与实际储层介质的接近程度,在一定程度上,决定了以此为理论基础所发展起来的储层预测和流体识别技术的有效性.传统的裂缝介质等效介质模型基于单相介质理论建立,主要包括Hudson模型(Hudson,1981)和线性滑动模型(Schoenberg and Sayers,1995),它们都没有考虑储层岩石的孔隙性和渗透性.而实际上,裂缝型储层岩石通常是既包含裂缝又包含孔隙的裂缝型孔隙介质,而且,随着研究的深入,人们逐渐认识到裂缝型油气藏中裂缝定向排列所引起的各向异性以及孔隙流体的存在所引起的双相或多相性对地震波的传播都有很大影响(Chapman,2003;Gurevich,2003;Gurevich et al.,2009;Sil et al.,2011).所以,建立能够综合考虑裂缝型储层岩石的各向异性和孔隙型的等效介质理论模型对于裂缝型储层的精细勘探尤为重要.

近年来,研究者们提出了一些综合考虑裂缝各向异性、孔隙性以及裂缝与孔隙之间流体交换作用的介质模型,并对弹性波的传播特征进行了研究.Thomsen(1995)Gurevich(2003)假设裂缝与球形等径孔隙之间流体压力平衡,提出了考虑孔隙和裂缝液压连通的裂缝介质模型;Sil等(2011)唐杰等(2015)将Gurevich模型应用于多孔隙HTI介质的流体替换分析中.Hudson等(19962001)考虑裂缝与孔隙之间的流体压力不完全平衡,推导了与频率相关的动态等效刚度张量.Chapman(2003)Jakobsen等(2003)Gurevich等(2009)提出了含有中观尺度(远大于孔隙尺寸,同时又小于弹性波波长)裂缝的裂缝介质模型,并通过波诱导的流体中观流动机制解释了一些裂缝型储层中地震频段内所观测到的频散和衰减现象.Parra(2000)结合BISQ模型(Dvorkin and Nur,1993)和Thomsen模型,提出了一种考虑微裂缝与等径孔隙之间流体局部喷射流动的等效介质模型,并研究了渗透率各向异性对地震波频散和衰减的影响.结合BISQ模型和Hudson模型,轩义华等(2006)张显文等(2010)分别建立了双相HTI介质和双相裂隙正交各向异性介质模型.刘财等(2013)结合改进的BISQ模型和Gurevich模型提出了又一种考虑微裂缝与等径孔隙之间流体局部喷射流动的等效介质模型,并采用伪谱法进行了波场数值模拟.杜启振等(2009)Du等(2012)孔丽云等(2012)结合孔隙弹性Biot模型和线性滑动模型建立了裂缝诱导HTI双孔隙介质模型,并进行了波场数值模拟与裂缝参数分析.高炜等(2014)结合Biot模型和Gurevich模型建立了一种裂缝诱导双相HTI介质模型,给出了弹性波的传播方程.不同模型在地震频段的适用性主要取决于岩石中流体的流动性(Batzle et al.,2006).

波动方程数值模拟是认识裂缝型储层中地震波传播规律的重要手段.在诸多波动方程数值解法中,伪谱法以其精度高、压制网格数值频散效果好的特点,被广泛应用于复杂介质的波场数值模拟中(刘洋和李承楚,2000;郭智奇等,2007;刘财等,2007;单启铜和乐友喜,2007巴晶等,2008;Ba et al.,2008张军舵等,2008李红星和陶春辉,2009;Carcione and Gurevich,2011),尤其适用于黏弹性介质、孔隙介质等耗散介质的波场模拟问题,因为在耗散介质中传播的弹性波会因能量衰减而伴随产生速度频散,波场数值模拟中采用伪谱法则可最大程度地避免任何频段内这种与介质性质相关的物理频散与数值计算带来的频散的混淆.然而,研究发现,采用伪谱法求解一阶速度-应力方程时,存在Nyquist误差问题,导致波场模拟结果中出现数值伪影现象.Witte和Richards(1987)提出了一种交错网格伪谱法(简称错格伪谱法)来克服这一问题.Özdenvar和McMechan(19961997)详细分析了Nyquist误差问题,从数学上论证了错格伪谱法较标准网格伪谱法的改进,并将其应用到标量方程、声波方程、弹性波方程和孔隙弹性方程的数值求解中.之后,该方法被其他一些学者应用于波场模拟中,例如:Chen(1996)采用错格伪谱法进行了黏声介质中波场模拟;Carcione等(1999)对黏弹各向异性介质和孔隙黏弹性介质中的波场进行了数值模拟;刘炯等(2008;Liu et al.,2010)模拟了缝洞型储层模型、随机介质中的地震波场;巴晶等(2010)模拟了礁、滩相储层的地震波场;杜增利等(2010)吴宝年等(2012)模拟了声学介质和弹性介质中的地震波场.

与前人的模型相比(杜启振等,2009;Du et al.,2012;孔丽云等,2012),裂缝诱导双相HTI介质模型具有物理意义明确、弹性参数计算更为简洁的优点.本文对裂缝诱导双相HTI介质中的地震波场进行研究,为了解决该类介质中弹性波传播一阶速度-应力方程的刚性问题并获得高精度的波场模拟结果,将错格伪谱法与时间分裂法相结合对波传播方程进行数值求解,并对模拟波场进行特征分析,为进一步深入认识实际裂缝型储层介质的地震波传播规律奠定理论基础.

2 裂缝诱导双相HTI介质模型及弹性波传播方程

实际裂缝型储层岩石通常既包含裂缝又包含孔隙,前者孔隙度很小(通常<0.1%)主要影响储层的渗透性,而后者主要存储储层中的流体.在波作用下,裂缝与孔隙之间产生流体压力梯度,流体发生流动以使压力达到新的平衡,这一过程对波的传播有重要影响.对于某些高渗透性裂隙岩石,在地震频段内,流体压力有充足时间在裂缝与孔隙之间达到平衡.基于这一考虑,杜启振等(2009,2012)将一组垂直定向排列的裂缝嵌入到流体饱和各向同性孔隙背景介质中,建立了裂缝诱导HTI双孔隙介质模型.采用类似的方式,高炜等(2014)建立了裂缝诱导双相HTI介质模型(图 1).

图 1 裂缝诱导双相HTI介质模型示意Fig. 1 Diagram of fracture-induced two-phase HTI medium

该模型是基于以下假设建立的(高炜,2014): ①干燥的裂缝型孔隙岩石由一组平行排列的垂直裂缝嵌入到满足统计各向同性的孔隙背景岩石中构成,孔隙之间以及孔隙与裂缝之间相互连通,裂缝和孔隙的最大尺寸远小于弹性波波长;②孔隙和裂缝由一种各向同性、具有黏滞性和可压缩性的流体所饱和;③在波传播方向上,固体骨架和流体之间存在相对位移,流体相对于固体的流动属于Poiseuille型流动,流体渗流遵守达西定律;孔隙和裂缝之间的流体压力在波的半个周期内有足够时间达到完全平衡;④忽略热弹性效应,孔隙流体与岩石基质也不发生化学作用.

基于以上假设,二维裂缝诱导双相HTI介质一阶速度-应力弹性波传播方程为(高炜,2014):

式中,

表示时间微分算子,则表示空间微分算子;τij(i.j=x,z)是总应力张量的分量,p为流体压力;vi是固体质点速度矢量v=(vx,vz)的分量,qi(i=x,z)是流体质点相对于固体质点的速度矢量q=(qx,qz)在δTi*方向上的分量,q=Φət(U-u),u=(ux,uz)和U=(Ux,Uz)分别是固体质点和流体质点的位移矢量;Φ是孔隙度;ρ=(1-Φ)ρs+Φρf是含流体孔隙介质的密度,ρs和ρf分别是固体和流体的密度;mi=Tiρf/Φ(i=1,3),Ti表示主轴方向的弯曲度;η表示流体的黏滞系数;k11k33分别表示x方向和z方向的渗透率;cIJSat(I,J=1,…,6)是流体饱和孔隙介质刚度张量的分量,可以写为裂缝参数、背景孔隙介质参数和流体性质参数的函数,具体计算表达式见相关文献(高炜等,2014).

3 波传播方程的数值解法3.1 刚性问题

由于双相介质中存在传播速度与单相介质相似的快纵波和传播速度非常慢的慢纵波,传播矩阵具有两种大小相差很大的不同特征值,表明一阶速度-应力微分方程组是刚性的.在数值模拟时需要采用非常小的时间步长才能满足稳定性条件,这将大大增加计算量.为了解决刚性问题,以便可以采用大步长的显式时间积分法对其进行数值求解,Carcione和Quiroga-Goode(1995)引入了一种时间分裂解法.在时间分裂法中,刚性一阶偏微分方程组被分成两部分,一部分为刚性,另一部分为非刚性,交替求解两个方程组,其中一个方程组的解作为另一个方程组解的初始值.

根据时间分裂法,可以将方程(1)分裂为刚性方程组和非刚性方程组,其中刚性方程组为

该刚性方程组可以解析求解,给定(n-1/2)Δt时刻的速度场值qxn-1/2qzn-1/2vxn-1/2vzn-1/2,通过方程组(2)式,可以得到(n+1/2)Δt时刻速度场值的中间解为

其中,λs(i)=-(η/ki)22(i),当η=0时,vx*=vxn-1/2vz*=vzn-1/2qx*=qxn-1/2qz*=qzn-1/2中间波场矢量解W=[vx*vz*qx*qz*τxxnτzznτxznpn]T作为非刚性方程组(4)式数值求解过程的初始值,采用错格伪谱法数值求解非刚性方程组(4)式所得到的波场矢量:

即为方程(1)式的数值解.

3.2 错格伪谱法离散递推格式

对波传播平面的空间区域和传播时间进行离散化,取x=iΔxx=(i±1/2)Δxz=jΔzz=(j±1/2)Δzt=nΔtt=(n±1/2)Δt,其中Δx、Δz和Δt分别表示空间离散采样步长和时间采样步长,i,j为整数,分别表示xz方向的空间离散采样点号,n表示离散时间采样点号.场分量(速度或应力分量)和介质物性参数以图 2的方式分配在整网格点或半网格点上.半网格点处的介质物性参数通过整网格点处的平均值来计算,则(i+1/2,j)和(i,j+1/2)处的ρρfT,Φ和η/k值分别为

图 2 二维裂缝诱导双相HTI介质速度-应力方程交错网格示意Fig. 2 Diagram of velocity-stress equation staggered grid in 2D fracture-induced two-phase HTI medium

这里a代表上述介质物性参数.而在(i+1/2,j+1/2)处的c55的值为

场量的一阶空间微分由错格傅里叶伪谱微分算子计算,若以φ表示场量,则其沿x1方向的一阶空间微分为

其中,φ是场量φ的傅里叶变换,kx(N)是Nyquist波数.对于时间导数采用二阶精度的中心差分近似.

因此,非刚性方程组(4)式的错格伪谱法离散格式为

4 数值算例与分析4.1 单层模型

设计一个单层裂缝诱导双相HTI介质模型,模型网格数为255×255,网格大小为10 m×10 m.震源采用流体体积源加垂直分量固相源,即将震源能量加载在流体压力和固相垂直分量上.震源中心点位置为(1270 m,1270 m),震源时间函数采用Ricker子波,震源函数形式如下:

式中,f0为震源中心频率(震源主频),t0=1/f0为震源子波延迟时间,(x0z0)表示震源中心位置,λ表示震源力空间作用的集中系数.震源子波主频为40 Hz,延迟时间为0.025 s.时间采样间隔为1 ms.

为消除人工边界反射,采用简便实用的Cerjan衰减边界(Cerjan et al.,1985),即沿人工边界向外扩充N个网格点,进入扩充区域内的波场通过乘以下面的因子

逐渐衰减为0,其中a为衰减系数,可通过试验确定最佳值.

由于所建立的裂缝诱导双相HTI介质模型考虑了流体的全局流动,根据Biot理论,其波场中可能会包含有慢纵波,但双相介质中流-固界面的黏滞摩擦阻碍流-固相对运动,将导致慢纵波快速衰减.为了研究慢纵波的传播特征,本文除进行含标准黏滞流体的黏滞相界情况下的弹性波传播数值模拟外,也开展流体为黏滞系数非常小的近似理想相界情况下的波场模拟.两种情况下的数值模拟所采用的介质物性参数除黏滞系数η不同以外,其余物性参数相同,见表 1,其中KsμsKf分别表示固体基质的体积模量、剪切模量和流体体积模量,ΔN和ΔT表示裂缝法向和切向弱度参数,ΦpΦc表示背景岩石孔隙度和裂缝孔隙度,其余参数物理意义同前文.在模拟近似理想相界时,取η=10-10Pa·s,而在模拟黏滞相界时,取水的黏滞系数η=0.001 Pa·s.图 3给出了近似理想相界情况下,固体和流体(相对于固体)质点速度垂直分量和水平分量在350 ms时的波场传播快照.相应地,图 4给出了黏滞相界情况下的波场传播快照.图 5是近似理想相界情况下,位于(108,108)网格点处的质点速度垂直分量时间记录图,图 6是黏滞相界情况下,该质点速度垂直分量时间记录图.

表 1 裂缝诱导双相HTI介质物性参数 Table 1 Parameters of fracture-induced two-phase HTI medium

图 3 裂缝诱导双相HTI介质中固相(A)和流相(相对于固相)(B)质点速度水平分量(a1,b1)和垂直分量(a2,b2)在350 ms时的波场传播快照(η=10-10 Pa·s)Fig. 3 Wavefield snapshots of solid phase (A) and fluid phase (relative to solid phase) (B) particle velocity horizontal component (a1,b1) and vertical component (a2,b2) in fracture-induced two-phase HTI medium at 350 ms (η=10-10 Pa·s)

图 4 裂缝诱导双相HTI介质中固相(A)和流相(相对于固相)(B)质点速度水平分量(a1,b1)和垂直分量(a2,b2)在350 ms时的波场传播快照(η=0.001 Pa·s)Fig. 4 Wavefield snapshots of solid phase (A) and fluid phase (relative to solid phase) (B) particle velocity horizontal component (a1,b1) and vertical component (a2,b2) in fracture-induced two-phase HTI medium at 350 ms (η=0.001 Pa·s)

图 5 近似理想相界情况固相(实线)和流相(虚线)质点(108,108)速度垂直分量时间记录(附图为固相质点速度分量时间记录中慢准P波的局部放大图)Fig. 5 Vertical component time record of solid phase(solid line) and fluid phase(dashed line) particle(108,108) in the ideal phase boundary case (top right corner is the enlarged figure of slow quasi-P wave in solid phase particle velocity component time record)

图 6 黏滞相界情况固相(实线)和流相(虚线)质点(108,108)速度垂直分量时间记录Fig. 6 Vertical component time record of solid phase (solid line) and fluid phase(dashed line) particle (108,108) in the viscous phase boundary case

图 3图 4中各分量的波场快照可以看出,波场中各波前面都非常清晰,没有可见网格数值频散的产生,表明了错格伪谱法在裂缝诱导双相HTI介质波场数值模拟中的有效性和精确性.

分析以上波场数值模拟结果可知:(1)震源在裂缝诱导双相HTI介质中激发了三种类型的波,从外向内分别为快准P波、准SV波和慢P波;(2)快准P波和准SV波类似于单相各向异性介质中的准纵波和准SV波,波前面呈椭圆形,准SV波的波前面出现波面尖角现象,这体现了裂缝诱导双相HTI介质中裂缝各向异性所引起的波场特征;(3)慢准P波在不同相界情况下具有不同的特征,当流体为近似理想相界时(图 3),慢纵波呈传播模式,波前面也呈椭圆向外传播,其传播速度明显低于快准P波和准SV波,而当流体为黏滞性时(图 4),慢纵波表现为扩散性,以静态模式出现在震源位置,如图 4中箭头所标示的慢P波;(4)结合图 3图 4的波场快照以及图 5图 6的质点速度垂直分量的时间记录可以看出:快准P波和准SV波引起固相质点和流相质点的同相位振动,慢准P波引起固相质点和流相质点反相位振动;近似理想相界情况下,快准P波和准SV波在固相中的振幅明显大于在流相中的振幅,慢准P波在固相中的振幅明显小于在流相中的振幅,即在流相中更容易观测到慢纵波;在黏滞相界情况下,快准P波和准SV波在固相和流相中振幅差异非常小,慢准P波在固相和流相中均呈扩散性,在震源附近快速衰减.

三种波的波前面传播特征、固相和流相质点的振动特性以及慢准P波的存在共同表明:裂缝诱导双相HTI介质具有裂缝各向异性和双相性,模型中背景孔隙介质主导了介质的双相性,而裂缝系统则诱导了介质的各向异性.

为了分析不同流体类型对裂缝诱导双相HTI介质中波传播特征的影响,将上述模型中的流体由水替换为气,介质模型的物性参数分别为:流体体积模量Kf=0.025 GPa、流体密度ρf=65 kg·m-3、流体黏滞系数η=2×10-5Pa·s,其余参数同表 1.图 7给出了裂缝诱导双相HTI介质含水和含气两种情况下的固相质点速度垂直分量的波场快照,图 8给出了固相质点(108,108)速度垂直分量的时间记录.从图中可以看出,快准P波在含气饱和裂缝诱导双相HTI介质中的传播速度比含水时的传播速度小,准SV波则与之相反.快准P波和准SV波在含不同流体的裂缝诱导双相HTI介质中传播时的速度大小关系,符合Gassmann方程,对于准P波,速度的大小由含流体多孔隙岩石的弹性模量和密度共同决定,因而,流体体积模量和密度共同影响地震波的传播速度,对于低频范围,含不同流体的岩石,其体积模量的差异比密度的差异对速度的影响更大,所以含水时速度大于含气时;而对于准SV波,不同流体密度的差异对速度的影响更大,所以,含气时速度大于含水时.另外,由于两种流体均是黏滞性流体,慢准P波表现出弥散性,快速衰减,波场中无法观测到.

图 7 裂缝诱导双相HTI介质含水(a)和含气(b)时固相质点速度垂直分量350 ms波场快照Fig. 7 Wavefield snapshots of solid phase particle velocity vertical component in fracture-induced two-phase HTI medium at 350 ms when the pore fluid is water (a) and gas (b)

图 8 裂缝诱导双相HTI介质含水(虚线)和含气(实线)时固相质点(108,108)速度垂直分量时间记录Fig. 8 Vertical component time record of solid phase particle(108,108) in fracture-induced two-phase HTI medium when the pore fluid is water(dotted line) and gas(solid line)
4.2 双层模型

通过一个两层模型来分析弹性波在双层裂缝诱导双相HTI介质分界面处的反射和透射特征.模型网格数为255×255,网格大小为10 m×10 m.界面位于1680 m处.上、下层均为裂缝诱导双相HTI介质,物性参数见表 2.震源采用流体体积注入源加垂直分量固相源,震源中心点位置为(1270 m,1480 m),震源函数和震源参数同上.时间采样间隔为0.5 ms.接收排列在深度1180 m处水平分布.为了研究慢纵波在介质分界面处的反射和透射特征,假设上层介质中含有非黏滞性流体,即上层介质为理想相界情况.对于下层介质,考虑含非黏滞性和黏滞性流体两种情况,即除考虑表 2中下层介质黏滞系数η=0外,还考虑黏滞系数η=0.001的情况.

表 2 裂缝诱导双相HTI介质双层模型的物性参数 Table 2 Parameters of fracture-induced HTI two-phase medium two-layer model

为了验证错格伪谱法在强波阻抗界面存在情况下波场模拟的优势,所设计的模型界面两侧波阻抗差较大.分别采用标准网格伪谱法和错格伪谱法进行波场模拟,并将两种方法的模拟结果进行对比分析.

图 9图 10分别是采用标准网格伪谱法和错格伪谱法进行波场模拟所得到的裂缝诱导双相HTI介质固相和流相质点速度分量的波场快照.从图 9可以看出,在标准网格伪谱法模拟得到的波场快照中存在明显的振铃拖尾现象,主要集中在垂直方向上,而在图 10所给出的错格伪谱法波场模拟结果中则不存在着这样的现象.这表明,采用错格伪谱法能够适应地下存在强波阻抗反射界面的波场数值模拟问题,与标准网格伪谱法相比,能够获得更高的精度.

图 9 标准网格伪谱法波场模拟所得的固相(A)和流相(B)质点速度水平分量(a1,b1)和垂直分量(a2,b2)在500 ms时的波场传播快照Fig. 9 Wavefield snapshots of solid phase (A) and fluid phase (B) particle velocity horizontal component (a1,b1) and vertical component (a2,b2) in standard grid pseudo-spectral method at 500 ms

图 10 错格伪谱法波场模拟所得的固相(A)和流相(B)质点速度水平分量(a1,b1)和垂直分量(a2,b2)在500 ms时的波场传播快照Fig. 10 Wavefield snapshots of solid phase (A) and fluid phase (B) particle velocity horizontal component (a1,b1) and vertical component (a2,b2) in staggered-grid pseudo-spectral method at 500 ms

图 11图 12分别给出了当下层介质为理想相界时,双层裂缝诱导双相HTI介质波场模拟的波场快照和合成地震记录,而图 13则给出了下层介质为黏滞相界时的波场快照.

图 11 固相(a1,b1)和流相(相对于固相)(a2,b2)质点速度水平分量(B)和垂直分量(A)在500 ms时的波场传播快照(下层介质为理想相界)
1 直达快准P波; 2 直达准SV波; 3 直达慢准P波; 4 反射快准P波; 5 反射准SV波; 6 反射慢准P波; 7 透射快准P波; 8 透射准SV波; 9 透射慢准P波; 10 快准P波产生的反射转换准SV波; 11 快准P波产生的反射转换慢准P波12 准SV波产生的折射转换快准P波; 13 准SV波产生的反射转换慢准P波; 14 慢准P波产生的反射转换快准P波; 15 慢准P波产生的反射转换准SV波; 16 快准P波产生的透射转换准SV; 17 快准P波产生的透射转换慢准P波; 18 准SV波产生的透射转换快准P波; 19 准SV波产生的透射转换慢准P波; 20 慢准P波产生的透射转换快准P波; 21 慢准P波产生的透射转换准SV波.
Fig. 11 Wavefield snapshots of solid phase (left) and fluid phase (right) particle velocity horizontal component (down) and vertical component (up) at 500 ms (the lower layer medium is the ideal phase boundary)
1 Direct fast quasi-P wave (qP1 wave); 2 Direct quasi-SV wave (qSV wave); 3 Direct slow quasi-P wave (qP2 wave); 4 Reflected qP1 wave; 5 Reflected qSV wave; 6 Reflected qP2 wave; 7 Transmitted qP1 wave; 8 Transmitted qSV wave; 9 Transmitted qP2 wave; 10 Reflected converted quasi-SV wave by fast quasi-P wave (qP1-qSV); 11 Refracted qSV-qP1; 12 Reflected qSV-qP1; 13 Reflected qSV-qP2; 14 Reflected qP2-qP1; 15 Reflected qP2-qSV; 16 Transmitted qP1-qSV; 17 Transmitted qP1-qP2; 18 Transmitted qSV-qP1; 19 Transmitted qSV-qP2; 20 Transmitted qP2-qP1, 21 Transmitted qP2-qSV.

图 12 固相(a1,b1)和流相(相对于固相)(a2,b2)速度水平分量(B)和垂直分量(A)合成地震记录Fig. 12 Synthetic seismogram of solid phase (relative to solid phase) (a1,b1) and fluid phase (a2,b2) velocity horizontal component (B) and vertical component (A)

图 13 固相(a1,b1)和流相(相对于固相)(b2,b2)质点速度水平分量(B)和 垂直分量(A)在500 ms时的波场传播快照(下层介质为黏滞相界)Fig. 13 Wavefield snapshots of solid phase (a1,b1) and fluid phase (a2,b2) particle velocity horizontal component (B) and vertical component (A) at 500 ms (the lower layer medium is the viscous phase boundary)

图 11可知:当两层裂缝诱导双相HTI介质中的流体均为非粘滞性时,快准P波、准SV波和慢准P波在介质分界面处均发生了反射和透射现象,并发生了波型转换,波场快照中可同时观测到多种波,这些波大致可以分为以下4类:(1)直达波类:直达快准P波、直达准SV波和直达慢准P波;(2)反射波类:反射快准P波、反射准SV波、反射慢准P波;(3)透射波类:透射快准P波、透射准SV波、透射慢准P波;(4)转换波类:由快准P波在界面反射之后形成的反射转换准SV波和慢准P波,由准SV波在界面折射之后形成的折射转换快准P波,由准SV波在界面反射之后形成的反射转换慢准P波,由慢准P波在界面反射之后形成的反射转换快准P波和准SV波,由快准P波在界面透射之后形成的透射转换准SV波和慢准P波,由准SV波在界面透射之后形成的透射转换快准P波和慢准P波,由慢准P波在界面透射之后形成的透射转换快准P波和准SV波.在合成地震记录中可以观测到反射域的各种弹性波,如图 12所示.从图 13中可以看出,当下层介质中的流体为粘滞性时,透射的慢准P波以及由快准P波和准SV波形成的透射转换慢准P波由于强衰减已观测不到.

5 结论与认识

本文采用错格伪谱法和时间分裂法求解裂缝诱导双相HTI介质的一阶速度-应力方程进行地震波场数值模拟,以此研究了裂缝诱导双相HTI介质单层模型和双层介质模型中的地震波传播规律,取得了以下结论和认识:

(1)错格伪谱法能够有效消除标准网格伪谱法波场模拟结果中存在的数值伪影现象,与时间分裂法相结合,能够获得稳定的、高精度的波场模拟结果,是一种进行裂缝诱导双相HTI介质中波场数值模拟的高精度算法.

(2)快准P波、准SV波和慢准P波在近似理想相界和黏滞相界情况下的波前面传播特征和质点振动特性表明,裂缝诱导双相HTI介质中的地震波场兼具裂缝各向异性介质和双相介质中传播的地震波的波场特征.

(3)各种波在介质分界面处均发生反射和透射现象,并发生波型转换,波场快照中可同时观测到多种波,这些波可以分为直达波、反射波、透射波和转换波4类,快准P波和准SV波与慢P波之间可发生相互转换.对裂缝诱导双相HTI介质分界面处地震波反射和透射现象的研究,为进一步加深对裂缝介质储层地震响应的认识奠定基础.

参考文献
[1] Ba J, Lu M H, Hu B, et al. 2008. The skeleton-relaxed model for poroviscoelastic media. Chinese J. Geophys. (in Chinese), 51(5): 1527-1537, doi: 10.3321/j.issn:0001-5733.2008.05.027.
[2] Ba J, Nie J X, Cao H, et al. 2008. Mesoscopic fluid flow simulation in double-porosity rocks. Geophysical Research Letters, 35(4): L04303.
[3] Ba J, Cao H, Yao F C, et al. 2010. Application of staggerred pseudospectrum method in simulating seismic wave field of reef beach facies reservoirs. OGP (in Chinese), 45(2): 177-184.
[4] Batzle M L, Han D H, Hofmann R. 2006. Fluid mobility and frequency-dependent seismic velocity-Direct measurements. Geophysics, 71(1): N1-N9.
[5] Carcione J M, Quiroga-Goode G. 1995. Some aspects of the physics and numerical modeling of Biot compressional waves. Journal of Computational Acoustics, 3(34): 261-280.
[6] Carcione J M, Helle H B. 1999. Numerical solution of the poroviscoelastic wave equation on a staggered mesh. J. Comput. Phys., 154(2): 520-527.
[7] Carcione J M. 1999. Staggered mesh for the anisotropic and viscoelastic wave equation. Geophysics, 64(6): 1863-1866.
[8] Carcione, J M., Gurevich, B., 2011. Differential form and numerical implementation of Biot's poroelasticity equations with squirt dissipation. Geophysics, 76(6),: N55-N64.
[9] Cerjan C, Kosloff D, Kosloff R, et al. 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equation. Geophysics, 50(4): 705-708.
[10] Chapman M. 2003. Frequency-dependent anisotropy due to meso-scale fractures in the presence of equant porosity. Geophysical Prospecting, 51(5): 369-379.
[11] Chen H W. 1996. Staggered-grid pseudospectral viscoacoustic wave field simulation in two-dimensional media. J. Acoust. Soc. Am., 100(1): 120-131.
[12] Du Q Z, Kong L Y, Han S C. 2009. Wavefield propagation characteristics in the fracture-induced anisotropic double-porosity medium. Chinese J. Geophys. (in Chinese), 52(4): 1049-1058, doi: 10.3969/j.issn.0001-5733.2009.04.022.
[13] Du Q Z, Wang X M, Ba J, et al. 2012. An equivalent medium model for wave simulation in fractured porous rocks. Geophysical Prospecting, 60(5): 940-956.
[14] Du Z L, Xu F, Gao H L. 2010. Pseudo spectral seismic wavefield simulation with staggered grid. GPP (in Chinese), 49(5): 430-437.
[15] Dvorkin J, Nur A. 1993. Dynamic poroelasticity: A unified model with the squirt and the Biot mechanisms. Geophysics, 58(4): 524-533.
[16] Gao W, Liu C, Guo Z Q, et al. 2014. Fracture-induced two-phase HTI medium model and its elastic wave propagation equations. Global Geology (in Chinese), 33(4): 904-909.
[17] Guo Z Q, Liu C, Yang B J, et al. 2007. Seismic wave fields modeling and feature in viscoelastic anisotropic media. Progress in Geophysics (in Chinese), 22(3): 804-810.
[18] Gurevich B. 2003. Elastic properties of saturated porous rocks with aligned fractures. Journal of Applied Geophysics, 54(3-4): 203-218.
[19] Gurevich B, Brajanovski R, Galvin R J, et al. 2009. P-wave dispersion and attenuation in fractured and porous reservoirs—poroelasticity approach. Geophysical Prospecting, 57(2): 225-237.
[20] Hudson J A. 1981. Wave speeds and attenuation of elastic waves in material containing cracks. Geophysical Journal of the Royal Astronomical SocietyInternational, 64(1): 133-150.
[21] Hudson J A, Liu E, Crampin S. 1996. The mechanical properties of materials with interconnected cracks and pores. Geophysical Journal International, 124(1): 105-112.
[22] Hudson J A, Pointer T, Liu E. 2001. Effective-medium theories for fluid-saturated materials with aligned cracks. Geophysical Prospecting, 49(5): 509-522.
[23] Jakobsen M, Johansen T A, McCann C. 2003. The acoustic signature of fluid flow in complex porous media. Journal of Applied Geophysics, 54(3-4): 219-246.
[24] Kong L Y, Wang Y B, Yang H Z. 2012. Fracture parameters analyses in fracture-induced HTI double-porosity medium. Chinese J. Geophys. (in Chinese), 55(1): 189-196, doi: 10.6038/j.issn.0001-5733.2012.01.018.
[25] Li H X, Tao C H. 2009. Features analysis of seismic wave field in two-phase anisotropic random medium with the pseudo-spectral method. Acta Physica Sinica (in Chinese), 58(4): 2836-2842.
[26] Liu C, Guo Z Q, Yang B J, et al. 2007. Analysis of reflection and transmission problems of waves in viscoelastic anisotropic media. Chinese J. Geophys.(in Chinese), 50(4): 1216-1224.
[27] Liu C, Lan H T, Guo Z Q, et al. 2013. Pseudo-spectral modeling and feature analysis of wave propagation in two-phase HTI medium based on reformulated BISQ mechanism. Chinese J. Geophys. (in Chinese), 56(10): 3461-3473, doi: 10.6038/cjg20131021.
[28] Liu J, Ma J W, Yang H Z, et al. 2008. Staggered-grid pseudo-spectrum simulation of fracture and cave reservoir. OGP (in Chinese), 43(6): 723-727.
[29] Liu J, Ba J, Ma J W, et al. 2010. An analysis of seismic attenuation in random porous media.Science. China Ser G-Phys. Mech. Astron., 53(4): 628-637.
[30] Liu Y, Li C C. 2000. The study on elastic wave propagation with the pseudo-spectral method in two-phase anisotropic medium. Acta Seismologica Sinica (in Chinese), 22(2): 132-138.
[31] Özdenvar T, McMechan G A. 1996. Causes and reduction of numerical arteifacts in pseudo-spectral wavefield extrapolation. Geophy. J. Int., 126(3): 819-828.
[32] Özdenvar T, McMechan G A. 1997. Algorithms for staggered-grid computations for poroelastic, elastic, acoustic, and scalar wave equations. Geophysical Prospecting, 45(3): 403-420.
[33] Parra J O. 2000. Poroelastic model to relate seismic wave attenuation and dispersion to permeability anisotropy. Geophysics, 65(1): 202-210.
[34] Schoenberg M, Sayers C M. 1995. Seismic anisotropy of fractured rock. Geophysics, 60(1): 204-211.
[35] Shan Q T, Yue Y X. 2007. Wavefield simulation of 2-D viscoelastic medium in Perfectly Matched Layer boundary. GPP (in Chinese), 46(2): 126-130.
[36] Sil S, Sen M K, Gurevich B. 2011. Analysis of fluid substitution in a porous and fractured medium. Geophysics, 76(3): WA157-WA166.
[37] Tang J, Fang B, Sun C Y, et al. 2015. Study of seismic wave propagation characteristics based on anisotropic fluid substitution in fractured medium. Geophysical Prospecting for Petroleum (in Chinese), 54(1): 1-8.
[38] Thomsen L. 1995. Elastic anisotropy due to aligned cracks in porous rock. Geophysical Prospecting, 43(6): 805-829.
[39] Witte D, Richards P G. 1987. Contribution to the pseudo-spectral method for computing synthetic seismograms.//Expanded Abstracts of 57th SEG Annual Int. Mtg., 517-519.
[40] Wu B N, Wu S Q, Huang Y, et al. 2012. Application of staggered grids in elastic wavefield modeling by pseudo-spectral method. GPP (in Chinese), 51(5): 440-445.
[41] Xuan Y H, He J D, Meng Q S, et al. 2006. Wavefield analysis of biphase EDA medium based on BISQ mechanism. OGP (in Chinese), 41(5): 550-556.
[42] Zhang J D, Yue Y X, Wang Y X. 2008. Numerical simulation wave field by pseudo-spectrum method in isotropic two-phase media. GPP (in Chinese), 47(4): 338-345.
[43] Zhang X W, Wang D L, Wang Z J, et al. 2010. The study on azimuth echaraecterlistiecs of attenuation and dispersion in 3D two-phase orthotropic crack medium based on BISQ mechanism. Chinese J. Geophys. (in Chinese), 53(10): 2452-2459, doi: 10.3969/j.issn.0001-5733.2010.10.019.
[44] 巴晶, 卢明辉, 胡彬等. 2008. 黏弹双相介质中的松弛骨架模型. 地球物理学报, 51(5): 1527-1537, doi: 10.3321/j.issn:0001-5733.2008.05.027.
[45] 巴晶, 曹宏, 姚逢昌等. 2010. 利用错格虚谱法模拟礁、滩相储层地震波场. 石油地球物理勘探, 45(2): 177-184.
[46] 杜启振, 孔丽云, 韩世春. 2009. 裂缝诱导各向异性双孔隙介质波场传播特征. 地球物理学报, 52(4): 1049-1058, doi: 10.3969/j.issn.0001-5733.2009.04.022.
[47] 杜增利, 徐峰, 高宏亮等. 2010. 虚谱法交错网格地震波场数值模拟. 石油物探, 49(5): 430-437.
[48] 高炜, 刘财, 郭智奇等. 2014. 裂缝诱导双相HTI介质模型及其弹性波传播方程. 世界地质, 33(4): 904-909.
[49] 郭智奇, 刘财, 杨宝俊等. 2007. 粘黏弹各向异性介质中地震波场模拟与特征. 地球物理学进展, 22(3): 804-810.
[50] 孔丽云, 王一博, 杨慧珠. 2012. 裂缝诱导HTI双孔隙介质中的裂缝参数分析. 地球物理学报, 55(1): 189-196, doi: 10.6038/j.issn.0001-5733.2012.01.018.
[51] 李红星, 陶春辉. 2009. 双相各向异性随机介质伪谱法地震波场特征分析. 物理学报, 58(4): 2836-2842.
[52] 刘财, 郭智奇, 杨宝俊等. 2007. 黏弹各向异性介质中波的反射与透射问题分析. 地球物理学报, 50(4): 1216-1224
[53] 刘财, 兰慧田, 郭智奇等. 2013. 基于改进BISQ机制的双相HTI介质波传播伪谱法模拟与特征分析. 地球物理学报, 56(10): 3461-3473, doi: 10.6038/cjg20131021.
[54] 刘炯, 马坚伟, 杨慧珠等. 2008. 缝洞型储层错格伪谱法模拟. 石油地球物理勘探, 43(6): 723-727.
[55] 刘洋, 李承楚. 2000. 双相各向异性介质中弹性波传播伪谱法数值模拟研究. 地震学报, 22(2): 132-138.
[56] 单启铜, 乐友喜. 2007. PML边界条件下二维粘弹性介质波场模拟. 石油物探, 46(2): 126-130.
[57] 唐杰, 方兵, 孙成禹等. 2015. 基于各向异性流体替换的裂隙介质波传播特征研究. 石油物探, 54(1): 1-8.
[58] 吴宝年, 吴肃琴, 黄义等. 2012. 交错网格在伪谱法弹性波场数值模拟中的应用. 石油物探, 51(5): 430440-437445.
[59] 轩义华, 何樵登, 孟庆生等. 2006. 基于BISQ机制的双相EDA介质的波场分析. 石油地球物理勘探, 41(5): 550-556.
[60] 张军舵, 乐友喜, 王艳香. 2008. 双相各向同性介质伪谱法地震波场数值模拟. 石油物探, 47(4): 338-345.
[61] 张显文, 王德利, 王者江等. 2010. 基于BISQ机制三维双相正交裂隙各向异性介质衰减及频散方位特性研究. 地球物理学报, 53(10): 2452-2459, doi: 10.3969/j.issn.0001-5733.2010.10.019.