高效高精度地震波初至行走计算

高效高精度地震波初至行走计算

0走时计算方法计算各向异性偏移时间是基于全波形反演、色度成像和其他各向异性地震模拟和成像技术的基础。走时计算方法的选择需要综合考虑算法的精度、效率及稳定性等因素。层析成像是一种层状介质速度和各向异性参数建模方法,其结果往往用于为全波形反演提供初始模型,但与波形反演中的旅行时信息不同,其核心之一就是模型初至波走时计算。层析成像是一种迭代反演算法,需要通过正演结果与实际数据的残差不断调整初始模型。走时计算的精度、效率与稳定性和层析成像计算的精度、效率与稳定性息息相关。因此,高效高精度且稳定的走时计算方法是各向异性介质偏移和层析成像等建模技术实现的关键。传统的走时计算方法是射线追踪方法,包括试射法(shootingtrace)和弯曲法(bendingtrace)射线追踪技术结合图形理论和复杂网格技术,在局部进行射线追踪。主要衍生出的方法有插值法和最短路径射线追踪算法等。插值法走时计算的另一种方法是波前面追踪,其包括两种方法。一种方法是基于运动学射线追踪系统的波前构建法;另一种方法是通过求解波动方程WKB(Wentzel,Kramers和Brillouin)展开的零阶近似后来衍生出的方法根据波前扩展方法的不同可分为迎风差分法FMM是由Sethian国内对于各向异性走时计算的方法研究起步较晚,其中以何樵登在各向异性介质中,走时计算不只依赖速度场,随着各向异性参数的增多,构建迎风差分格式求取走时也越来越困难。本文主要利用扰动理论及泰勒公式将VTI介质程函方程展开,完成了对VTI介质程函方程的降维处理。对降维后的VTI介质程函方程构建迎风差分格式并分别进行求解,完成单点走时计算。结合窄带推进理论,构成了各向异性的波前面追踪走时计算方法。首先对均匀弱VTI各向异性介质进行计算并与解析值相减评估其误差;然后为验证其稳定性和有效性,又计算了VTI各向异性层状介质模型的走时;最后以盐丘模型作为测试模型,进行了走时正演试算。1算法的原理和实现1.1vti介质程函方程离散方程在VTI介质中,声波近似的程函方程为以下格式式中:τ(x,y,z)是震源到坐标点(x,y,z)的走时;v此四次方程的非线性程度强于各向同性和椭圆各向异性的程函方程,因此VTI介质程函方程需要更为复杂的有限差分近似。为解此方程,我们可以利用扰动理论式中,τ方程(4)是椭圆各向异性介质的程函方程,独立参数η相同阶数项的系数相等。首先计算一次η项的系数,将一次η项系数相等的公式通过式(4)简化后得到同理,可得到η平方项的系数公式:1.2基于椭圆异性介质模型的背景模型为求取VTI介质走时的解析解,首先从均匀介质背景速度模型出发。为求取式(4)的解析解,设η=0,此时背景模型是椭圆各向异性介质。在这种情况下,震源位于x=0、z=0处的走时解析解如下:该结果满足方程(4)。通过式(7)可求出同理,计算将式(7)、(8)、(9)代入式(3)中,即可求取该点走时的解析值。1.3非均匀介质走时计算在均匀介质中,上述公式的推导是基于VTI介质声波近似程函方程(1)进行的,因此得到的近似结果描述的是介质的群速度。而在非均匀介质正方形网格条件下,由于无法确定指定网格中的震源位置,因此上述公式不能简单地应用于非均匀介质走时计算中。对于此种情况,我们引入FMM。FMM的基本原理是将程函方程的时间梯度项近似为满足波前传播“熵守恒”理论的数值离散形式,通过数值计算来求取黏滞解。这种离散形式为式中:D该方法对于波前传播过程中尖角拐角和拓扑结构变化等问题的有效性在文献然后,将式(12)代入式(4)中,通过计算可以求出点(i,j)处的τ1.4窄带时走时计算FMM的实现策略主要为“熵守恒”理论,该理论可以这样理解:传播的初至波波前可以看做在纸上燃烧的火焰。如果纸上的某个位置被火焰燃烧了,那么这一位置将一直保持着被烧过的状态,后续不会再被燃烧。类比于初至波在空间中的传播,即初至波波前不会两次经过同一个点,每个点的走时是关于位置的单值函数。同时由于波是单向传播的,总是由走时较小点传播至走时较大点,因此,1.3节中所述单点走时计算的循环起点应位于走时最小值点。文中采用了窄带技术作为波前模拟技术,通过排序选择窄带内走时最小的点作为延拓点,再通过对网格节点属性的设定来“储存”已经完成计算的网格节点,最后通过窄带内网格节点的走时更新来模拟波前的演化。整体步骤可分为初始化和循环两部分。12启动式中:l为网格节点的间隔长度;vi2确定窄带以确定深度循环过程分为以下4个步骤。(1)从窄带内选取最小走时点τ(i(2)将点i((3)对点i((4)判断窄带是否为空,如果窄带为空则终止循环,否则跳回(1)继续计算。2算法的精度和稳定性分析2.1走时计算与解析值对比由于复杂界面非均匀模型的走时解析值较难求取,为了便于分析本文方法的正确性及精度,我们首先用一个水平界面均匀弱VTI各向异性模型计算走时,并用走时的解析值与之对比(由于后文复杂模型通过速度及各向异性参数δ和ε定义,因此本文中所有模型均用v水平界面均匀介质模型的大小为4000m×1200m,介质速度为4000m/s,各向异性参数δ=0.1,ε=0.1,网格间距为5m,震源位置为(2000m,600m)。使用本文方法的走时计算结果见图1a。走时解析值通过1.2节中叙述的方法求取,将走时计算值(图1a)与解析值相减后得到误差,误差占解析值的比例(相对误差)见图1b。由图1b可知:在震源附近走时计算结果的相对误差较为明显,主要分布在1.5%～4.0%之间,其原因是差分网格长度太大,对震源处的差分网格进行进一步细分,该现象会有所改善;在远偏处误差趋于稳定,相对误差稳定在0.5%以下。2.2算法的稳定性为了验证本算法对复杂介质模型的稳定性和有效性,分别选取双层、多层和盐丘介质模型进行试算。2.2.1走时模型试算对双层VTI各向异性介质模型进行数据试算。该模型的大小为4000m×1200m,震源点位于坐标(0,0)处。模型参数见表1,模型结构见图2a,使用本文方法试算的走时结果见图2b。地震观测点最先接收到的波称为初至波。在层状介质中,随着偏移距及入射角的增大,下层介质的折射波渐渐取代直达波被率先接收,在某节点产生“毛刺现象”。通过观察图2b可以发现,第一层近偏处与均匀介质传播规律一致,随着偏移距增加,在坐标(500,600)处产生了“毛刺现象”,折射波取代了直达波率先被接收。2.2.2多层模型测试结果对多层VTI各向异性介质模型进行数据试算。该模型的大小为4000m×1200m,震源点位于坐标(0,0)处。模型参数见表2,模型结构见图3a,使用本文方法试算的走时结果见如图3b。由于多层模型的速度变化与双层模型相比较慢,因此其等时线变化较为平缓,近偏处对速度突变层的表现不明显。但在偏移距1000～2500m处,我们可以较清晰地分辨出3个层的折射波对于初至走时的影响,偏移距2500m之后,第二、三层的折射波渐渐被第四层的折射波所取代,并体现于等时线上。2.2.3波场快照的正演模拟对盐丘介质模型进行数据试算。该模型大小为21940m×9144m,模型速度及各向异性参数结构如图4所示。该模型x轴2000～10000m处为一较明显的高速体,15000～19000m处为一断层,两个比较突出的结构体之间是由起伏层状介质连接(图4a)。在层状地质体中,x轴10000～14000m、深度约4000m处有一强各向异性体(图4b、c)。我们在模型x轴方向0、3650、12200、19000m处分别设置炮点,观察其在特殊地形下的试算效果,并与有限差分正演的波场快照进行对比,结果如图5所示。从图5a可以看出,高速体对走时场的影响十分明显。通过与波场快照(图5b)对比可以发现,排除波场快照中界面之间反射波的影响,3000ms等值线与波场快照中的波前面完全吻合。由于图5c中炮点位置位于高速体深度最小的地面投影处,因此选择2000ms的波场快照(图5d)进行对比。可以看出:波场快照0～3000m深度处与等时线吻合;5000m以下受到高速体的影响,波场快照中波前发生分裂,但是最高能量的波前面与等时线依然吻合。排除反射波的影响,从图5e、f可以清晰地看到,不论是在深度0～3000m的层状介质处,还是在深度3000～5000m的高速体和断层处,波场快照都与波前面清晰地吻合。由于图5g中震源点的位置是在断层正上方,为观察到高速体的影响,因此将波场快照的时间设置为5500ms(图5h)。可以清晰地看到深度0～3000m的波前面与等时线吻合;而对于高速体,由于偏移距过大的原因,有限差分正演受到反射波、绕射波等因素的影响,波场快照中最高能量波前面较难分辨。综上所述,经过多个模型的测试,本文提出的方法在不同地下构造及不同偏移距下的稳定性和准确性得到了较好的验证,其计算结果合乎地震波传播规律,能够准确地反映出地下介质的结构。3vti介质声波走时计算方法的适用性1)本文将基于波前面追踪的走时计算方法引入各向异性介质走时计算当中,实现了网格化非均匀VTI介质走时计算。首先对程函方程中的走时参数τ进行展开,成功将四次方程分解为一个二次方程和两个一次方程;然后将各向同性快速匹配走时计算方法引入VTI介质中,对展开的VTI介质程函方程组进行离散化并求解,得到介质的走时。2)大量模型试算检验证明,本文提出的VTI介质声波走时计算方法是有效的,该方法的模型试算结果能够较好地反映出地下复杂介质中地震波的传播规律。虽然本文中用于试算的模型不能涵盖地下介质的所有情况,但本文中的计算方