三维似French模型高精度弹性波正演数值模拟及其波场特征

刘阿男,刘思彤,陈可洋,杨 微

(1.东北石油大学地球科学学院,黑龙江大庆 163318; 2.大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712; 3.大庆油田有限责任公司第六采油厂,黑龙江大庆 163712)

三维似French模型高精度弹性波正演数值模拟及其波场特征

刘阿男1,刘思彤1,陈可洋2,杨 微3

(1.东北石油大学地球科学学院,黑龙江大庆 163318; 2.大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712; 3.大庆油田有限责任公司第六采油厂,黑龙江大庆 163712)

基于三维一阶双曲型各向异性介质速度-应力弹性波动方程,给出各向异性介质向各向同性介质转化的条件,采用高阶交错网格有限差分法和PML吸收边界条件,实现高精度三维弹性波正演数值模拟.以均匀介质为例,研究弹性波在三维空间的波场特征;以三维似French模型为例,合成不同时刻的三分量波场快照和共炮点记录,并对比分析三维各向同性介质和各向异性介质在该模型中的波场响应特征.数值结果验证了方法的正确性和有效性,同时具有精度高和边界吸收效果明显的优点.

三维弹性介质;高精度正演数值模拟;PML吸收边界条件;似French模型;高阶交错网格有限差分法

0 引言

弹性波正演数值模拟在多分量地震资料处理、观测系统采集设计、处理参数试验、解释方案验证等方面发挥越来越重要作用,是地震学和地震勘探学研究领域中的前沿课题之一,也是难度较大的课题之一[1].近些年,随计算机计算能力的增强该技术得到飞跃发展.弹性波正演数值模拟技术大多是基于二维波场观测研究地震勘探问题,原因在于二维弹性波场的计算速度快、算法易实现、剖面易显示,而实际的地球介质是三维空间,二维波场观测忽略三维构造体的侧向效应等很多重要信息,无法准确、真实地刻画三维空间的弹性波场响应特征.只有通过三维弹性波正演数值模拟,才能全面地认识弹性介质与弹性波场之间的耦合关系及弹性波场的空间变化特征,后者对于研究和认识地球介质的各向异性特征尤为重要,然而三维弹性波场正演问题研究困难,其原因在于:三维空间的计算变量较多;空间网格节点配置关系复杂;边界吸收问题处理较难;剖面立体显示困难;计算量和存储量巨大,三维声波正演数值模拟的算法要比弹性波正演数值模拟的算法更为简单、易实现.

目前,对三维声波正演问题已得到丰富的研究成果[2-7],如频散问题[8-15](精度问题)、边界吸收问题[16]、稳定性问题[17]等,而对三维弹性波正演问题的研究相对较少,且大多是基于均匀单层介质研究弹性波场的传播特征(波场响应),如IgelH等[18]讨论三维各向异性介质中弹性波数值模拟的有限差分实现方法;张文生等[19]对三维正交各向异性介质弹性波场进行有限差分模拟;裴正林[20-21]采用高阶交错网格有限差分法对三维各向异性介质中弹性波场进行高精度数值模拟;董良国等[17]提出三维任意各向异性介质弹性波动方程的稳定性条件;夏凡等[22-23]提出三维均匀介质弹性波动方程的吸收边界条件及其高精度数值计算方法.

笔者借助Linux集群计算机的高速运算能力和大容量存储能力,采用高阶交错网格有限差分法求解三维各向异性介质弹性波动方程,通过数值实例分析弹性波在三维空间经典似French速度模型中的三分量波场响应特征,给出三维各向异性介质向各向同性介质转化的条件,同时使用任意方位波场切片技术(包括波场快照切片、模拟记录切片等),实现三分量弹性波场全方位可视化显示,研究三分量弹性波场的传播规律,以此验证方法的准确有效性.

1 基本理论

通常三维二阶各向异性介质(VTI介质)弹性波动方程可表示为

式中:U、V、W分别为x、y、z方向质点的振动位移;C为三维弹性张量(矩阵组合形式C也可称为物性矩阵,一般情况下沿主对角线方向是对称正定的);ρ为弹性介质密度;t为时间.参数均为空间坐标的函数.

当各向异性介质对称轴与观测系统不一致时,需要进行观测坐标旋转或者物性矩阵旋转,以达到合理观测目的.此时需要用到Bond矩阵旋转变换,对物性矩阵旋转运算方程为C′=MCMT,其中M为正交正定矩阵,C′为变换后的物性矩阵.

根据Thomsen参数的定义,准纵波相速度Vp0=C33/ρ;准横波相速度Vs0=C44/ρ;度量准纵波各向异性强度的参数ε=(C11-C33)/(2C33);度量准横波各向异性强度的参数γ=(C66-C44)/(2C44);影响垂直于VTI介质对称轴方向附加准纵波速度的参数δ=[(C13+C44)2-(C33-C44)2]/[2C33(C33-C44)].

为实现三维各向异性介质向各向同性介质转化,借用各向同性介质中的泊松比ν概念,建立三维VTI介质中准纵波和准横波相速度之间的关系式,即

同时可推得

为便于应用高阶交错网格有限差分法,将三维二阶各向异性弹性波动方程(式(1))变换为一阶双曲型速度—应力弹性波动方程的微分形式,即

采用高阶交错网格有限差分法(时间2阶、空间10阶差分精度)对式(2)进行高精度差分离散(形成差分行式),对离散方程中的各变量进行交错网格节点配置,可得到三维高精度弹性波正演递推算子.在人工截断边界处,采用统一格式的最佳匹配层吸收边界条件(吸收层厚度取为20个节点),消除或削弱截断边界处的反射波.要保证边界计算过程稳定,还需满足差分迭代所需的稳定性条件[17].

2 波场特征

2.1 均匀介质三维三分量

三维均匀各向同性介质模型尺寸为500m×500m×500m,3个方向空间网格长度为5m,采用纯纵波震源(在3个方向正应力场中加入等量外力扰动)并置震源于模型中央,其最大频率为40Hz(波形为雷克子波),介质的纵波速度为2000m/s;泊松比为0.25;密度为1g/cm3;时间步长取0.5ms.满足稳定性条件,波场快照的记录时间为0.125s.波场任意方向切片显示方式的定义:(x=A,yoz)代表沿x方向深度A处且在yoz平面内的波场切片(其中y为图片的横坐标;z为图片的纵坐标);(t=B,yoz)代表在时间B节点处且在yoz平面内的波场切片,其他情况依此类推.

三维任一分量弹性波场在3个相互垂直方向上的波场快照切片见图1.由图1可知:在纯纵波震源激励下,x方向质点速度场中纵波波形与yoz平面对称;y方向质点速度场中纵波波形与xoz平面对称;z方向质点速度场中纵波波形与yox平面对称,且对称面两侧极性相反.这验证在均匀各向同性介质中,三维纯纵波波场具有球对称的特点.纯横波三分量弹性波场在任意方向的质点振动速度和与其相垂直的平面存在对称关系,且极性相一致(见图2,在x分量质点的振动速度处加入外力扰动).数值计算结果表明计算精度较高、频散现象很弱.该三维三分量弹性波场特征分析方法同样适合于表征三维各向异性介质准纵波和准横波波场的特点,受各向异性参数的影响,准纵波和准横波波形不再是球对称关系,而是椭圆形甚至是具有不规则特征的对称关系.

2.2 三维似French模型

图1 0.125ms时刻纵波震源情况下三维各向同性介质弹性波波场快照

采用似French模型的速度模型,包含1个水平层界面、1个倾斜断层、1个隆起构造、1个凹陷构造(见图3).模型尺寸为500m×500m×500m,3个方向空间网格长度为5m,震源位于模型地表中央位置,其最大频率为50Hz,每层密度均为1g/cm3,时间步长为0.5ms,满足稳定性条件,接收时间为1s.对于三维各向同性介质模型,上层介质纵波速度为2000m/s,下层介质纵波速度为2500m/s,泊松比为0.25;对于三维各向异性介质模型,上层介质准纵波速度为2000m/s,下层介质准纵波速度为2500m/s,泊松比为0.25,ε为0.1,δ为0.2,γ为0.

图2 0.125ms时刻混合源情况下三维各向同性介质弹性波波场快照(x分量)

图3 三维似French模型(单位:m)

在0.20,0.25s时,三维各向同性介质和各向异性介质波场快照中的x、y、z方向的中间位置,且沿平行于yoz、xoz、xoy平面的波场快照切片(只给出三维垂直分量的波场快照)见图4和图5.考虑三维似French速度模型与三维各向同性介质弹性波场快照,经过对比容易识别直达P波、倾斜界面的反射PP波和反射PS波、水平层界面引起的反射PP波和反射PS波、隆起构造的反射PP波和反射PS波、凹陷构造引起的回转反射PP波和反射PS波;弹性波传播至下层介质时形成的透射PP波和透射PS波,以及弹性波在隆起构造和凹陷构造之间相互作用形成的复杂波(包括转换波能量).与三维各向异性介质弹性波场快照对比,弹性波场更为复杂,波场快照中包含由各向异性特性引起的直达S波,模型中构造形成的反射SP波、反射SS波、透射SP波、透射SS波,以及这些波在构造之间相互作用形成的复杂波形,直达S波的波场能量比其他反射波能量要强.由图4可知,在不同的切片中,凹陷构造和隆起构造的侧向传播效应较为明显.三维各向同性介质和各向异性介质数值模拟记录的在不同位置处的波场切片(见图5),两者最明显的特征在于各向异性介质中存在直达S波及其转换波,其波场情况明显比各向同性介质更为复杂.由此可知,三维弹性波在三维空间中进行传播,仅考虑二维空间的弹性波传播不够准确和全面.另外,地震波到达模型边界时无边界反射波形成,三维弹性波场的传播过程中无明显的频散现象,说明数值计算方法具有较高精度.

图4 三维弹性波场快照切片显示

图5 三维弹性波模拟记录切片显示

3 结论

(1)采用高阶交错网格有限差分法求解三维一阶双曲型各向异性介质速度—应力弹性波动方程及统一格式的PML吸收边界条件方程构建三维高精度弹性波正演递推算子,给出三维VTI介质与三维各向同性介质之间的转换关系,实现了三维均匀介质和似French模型的弹性波高精度正演数值模拟.

(2)使用任意方位弹性波场切片技术实现三维三分量弹性波场全方位可视化,研究三分量弹性波场在三维空间中的波场特征及其传播规律,纵波或准纵波具有平面对称、极性相反的特点;横波或准横波具有平面对称、极性相同的特点.

(3)三维弹性波场正演数值模拟技术可以更加全面地表征弹性波在全三维空间中的传播过程,计算方法准确有效.

[1] 裴正林,牟永光.地震波传播的数值模拟[J].地球物理学进展,2004,19(4):933-941.

[2] 刘金平,张向君,董建平.波动方程显式短算子三维正演模拟[J].吉林大学学报:地球科学版,2004,34(1):138-141.

[3] 王润秋,张明,费建博,等.精细积分法三维地震波正演模拟[J].勘探地球物理进展,2006,29(6):394-397.

[4] 何兵寿,魏修成,刘洋.三维波动方程的数值频散关系及其叠前和叠后数值模拟[J].石油大学学报:自然科学版,2001,25(1):67-71.

[5] 熊晓军,贺振华,黄德济.三维波动方程正演及模型应用研究[J].石油物探,2005,44(6):554-556.

[6] 熊高君,张学工,贺振华,等.三维定位原理与三维反射波场模拟[J].矿物岩石,2002,22(3):93-97.

[7] 熊高君,贺振华,张毅祥.三维混合延拓一步法波动方程正演模拟[J].成都理工学院学报,1999,26(1):48-51.

[8] Alford R M,Kelly K R,Boore D M.Accuracy of finite-difference modeling of the acoustic wave equation[J].Geophysics,1974,39 (6):834-842.

[9] Marfurt KJ.Accuracy of finite difference and finite element modeling of the scalar and elastic equations[J].Geophysics,1984,49(5): 533-549.

[10] Dablain M A.The application of high-differencing to the scalar wave equation[J].Geophysics,1986,51(1):54-66.

[11] Fei T,Larner K.Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport[J].Geophysics,1995,60(6):1830-1842.

[12] 董良国,李培明.地震波传播数值模拟中的频散问题[J].天然气工业,2004,24(6):53-56.

[13] 吴国忱,王华忠.波场模拟中的数值频散分析与校正策略[J].地球物理学进展,2005,20(1):58-65.

[14] 宁刚,熊章强,陈持.波动方程有限差分正演模拟误差来源分析[J].物探与化探,2008,32(2):203-206.

[15] 董良国,马在田,曹景忠,等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报,2000,43(3):411-419.

[16] 陈可洋,杨微.优化的三维地震波旁轴近似吸收边界条件[J].勘探地球物理进展,2009,32(3):179-181,206.

[17] 董良国,马在田,曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究[J].地球物理学报,2000,43(6):856-864.

[18] Igel H,Morap,Riolletb.Anisotropic wave propagation through finite-difference grids[J].Geophysics,1995,60(4):1203-1216.

[19] 张文生,宋海斌.三维正交各向异性介质三分量高精度有限差分正演模拟[J].石油地球物理勘探,2001,36(4):422-432.

[20] 裴正林.三维各向异性介质中弹性波方程交错网格高阶有限差分法数值模拟[J].石油大学学报:自然科学版,2004,28(5):23-29.

[21] 裴正林.三维各向同性介质弹性波方程交错网格高阶有限差分法模拟[J].石油物探,2005,44(4):308-315.

[22] Xia F,Dong L G,Ma Z T.Absorbing boundary conditions in 3D elastic-wave numerical modeling[J].Chinese J.Geophys,2004,47 (1):132-136.

[23] Xia F,Dong L G,Ma Z T.The numerical modeling of 3-D elastic wave equation using a high-order,staggered-grid,finite difference scheme[J].Applied Geophysics,2004,1(1):38-41.

High-precision elastic w ave forw ard numerical simulation and its w ave field characteristic analysis on 3D like French model/2010,34(4):45-50

LIU A-nan1,LIU Si-tong1,CHEN Ke-yang2,YANG Wei3
(1.Geoscience College,Northeast Petroleum University,Daqing,Heilongjiang163318,China;2. Ex ploration and Development Institute ofDaqing Oilf ield Company L td.,Daqing,Heilongjiang 163712,China;3.Oil Recovery Plant N o.6,Daqing Oilf ield Corp.L td.,Daqing,Heilongjiang 163712,China)

3D elastic wave forward numerical simulation is the key means for studying seismic wave’s wave field characteristic and its propagating rules in complex 3D spatial elastic wave medium.Based on 3D one-order hyperbolic anisotropic medium velocity-stress elastic wave equation,the paper gives out the transformation condition from anisotropic medium to isotropic medium,and uses high-order staggered-grid finite-difference method and PML absorbing boundary condition to achieve high-precision 3D elastic wave forward numerical simulation.Take an isotropic medium as an example to study the wave field characteristic in 3D space,and take 3D like French model as another model to obtain three-component snapshots with different time and common-shot record,compare and analyze the wave field response in this model of 3D isotropic medium and anisotropic medium.The numerical result verifies the method’s correctness and effectiveness,and shows the advantages of high-precision and obvious perfect boundary absorption effect.

3D elastic medium ; high-precision forward numerical simulation ; PML absorbing boundarycondition ; French like model ; high-order staggered-grid finite-difference met hod

book=4,ebook=364

TE19

A

1000-1891(2010)04-0045-06

2010-06-22;审稿人:王云专;编辑:任志平,张兆虹

国家“973”重大基础研究项目(2001CB209104)

刘阿男(1973-),女,硕士,工程师,主要从事油气成藏方面的研究.