基于时空域自适应高阶有限差分的声波叠前逆时偏移

基于时空域自适应高阶有限差分的声波叠前逆时偏移

1叠前逆时偏移思想20世纪80年代的逆偏移技术开始了，该技术可以恢复和返回复杂路径上的信号，以便能够准确地绘制复杂的地下结构。因此，作为目前偏移方法中最精确的成像方法之一，反偏移是一种基于波动理论的深度偏移方法。常用的波动方程包括有限差分法、有限差分法、伪光谱法等。逆时位移思想的第一个起源是利用有限差分法求解波动方程。叠前逆时偏移实现过程中的主要步骤是波场的正向延拓和反向延拓、成像条件以及成像后的低频噪声的消除,其中,波场正向延拓和反向延拓的本质就是求解波动方程.求解波动方程是为了描述地震波场的传播特征,理论上利用全波场的信息对地下构造进行准确成像,所以精确、快速地求解波动方程至关重要.Karazincir和Gerrard然而,常规的有限差分法求解波动方程时,其求空间导数时的差分系数一般都是采用传统的高阶有限差分系数2方法原理2.1数值模拟结果分析逆时偏移的核心问题之一是求解波动方程.在纵波勘探中,一般采用声波方程近似描述地震波在地下介质中的传播,二维声波方程为其中,p为标量波场,v为速度,t为时间,x、z为空间坐标.一般情况下,计算时间二阶导数时,如果采用高阶精度差分形式,则计算量会很大,并且很容易造成不稳定其中,τ为时间采样间隔,p方程(1)中空间二阶导数的2M阶精度差分形式可以写为其中,a将式(2)和式(3)代入声波方程(1)中,可以得到:逆时偏移的正向延拓就是差分方程正演模拟的过程,由式(4)可以推导出正向延拓的表达式:其中,f逆时偏移波场的反向延拓过程就是已知后一时刻的波场值来求前一时刻的波场值,具体计算过程与正向延拓类似,其计算表达式为Liu和Sen其中,令f(θ)=cos为了计算和分析基于时空域有限差分的计算精度,我们分析差分的频散,由式(7)定义频散速度v用频散速度与真实速度的比值来衡量二维声波方程的频散现象,则可以得到二维声波方程的频散分析表达式由于δ(θ)=δ(θ+π/2),所以δ(θ)是一个以π/2为周期的函数.同时考虑到δ(θ)=δ(π/2-θ),所以在计算δ(θ)时,θ的变化是从0到π/4.从式(11)中可以看出,如果δ=1,则模拟时不存在数值频散,如果δ越是偏离1,则存在越大的数值频散.由式(7)变形可得到:从式(12)可以看出,由于由于a其中M设计一个具体模型进行精度分析:介质的速度为3000m/s,模型网格大小为10m×10m,时间步长为0.001s,空间导数精度阶数2M分别为4、8、16、32.图1(a,b)是波在传播方向θ=π/8上,分别采用传统有限差分和时空域有限差分数值模拟时的频散分析曲线图.从图1中可以看出:在相同参数情况下,时空域有限差分的频散程度比传统有限差分的频散程度要小,即时空域有限差分的模拟精度要高于传统有限差分如果空间导数的精度阶数2M固定为16,而波的传播角度θ分别取0、π/16、2π/16、3π/16、4π/16,其按顺序对应的标识为angle1、angle2、angle3、angle4、angle5.图2是有限差分数值模拟在不同方向上频散分析曲线图.从图2(a,b)的频散分析曲线图对比也可以看出:在精度阶数相同时,在波的各个传播方向上,时空域有限差分的模拟精度比传统有限差分高为比较传统有限差分和基于时空域有限差分的稳定性,由式(12)可令稳定因子s为图3是声波方程有限差分在不同的精度阶数(M分别为2、4、8、16)时随数值模拟参数变化的稳定性分析曲线图.从图3中可以看到,当精度阶数越高时,越容易造成不稳定,且时空域有限差分的稳定因子与精度阶数2M和数值模拟参数r有关,而传统有限差分的稳定因子s只与精度阶数2M有关,同时对比两类差分的稳定性分析图,可以看到,在保持求解方程稳定的前提下(s≥r),时空域的有限差分对应的r可取更大更广的范围.所以图3中表明,时空域的有限差分的稳性比传统有限差分的稳定性更好2.2自适应差分算子长度有限差分策略在给定一个地质模型后,采用高阶有限差分进行数值模拟时,通常都是采用固定长度的空间差分算子来计算空间导数.差分算子长度的选择必须同时满足稳定性、精度、计算量等的要求,而这些需要满足的要求又与介质模型的速度有关.对于给定的速度模型,如果整体采用一个固定长度的空间差分算子来进行数值模拟,则这个固定的空间差分算子必须同时满足模型中所有速度各自应该满足的条件.因此,固定的空间差分算子长度一般会要求比较长,导致求解波动方程耗费大量计算时间.Liu和Sen求取自适应差分算子长度时,先定义有限差分数值模拟的误差计算公式给定最大频率f设计一个速度范围为1500~5500m/s,网格大小为10m×10m,时间步长为0.001s,其最大频率为50Hz,而计算误差范围分别控制在102.3成像结果(s)成像条件也是逆时偏移处理过程中的一个重要步骤,成像条件的好坏直接影响着最终成像剖面的质量.目前叠前逆时偏移中常用的成像条件主要有三种类型:激发时间成像条件其中,I(x,z)是点(x,z)的成像结果;s(x,z,t)表示时间域的正向延拓波场;r(x,z,ω)表示时间域的反向延拓波场;ε为无穷小实数,其作用是为了保证计算的稳定性.由于时间域计算相互关的速度相对比较慢,为了提高其计算效率,本文在实现过程中将波场转换到频率域计算互相关,式(16)可以转换为其中,ω表示角频率,S(x,z,ω)表示频率域的正向延拓波场,R(x,z,ω)表示频率域的反向延拓波场,ue494(x,z,ω)表示S(x,z,ω)的共轭.2.4拉普拉斯算子去噪叠前逆时偏移的过程中,不可避免地会产生低频噪声,所以要得到最终的清晰成像结果,还要对偏移的初始剖面进行去噪处理.拉普拉斯算子去噪方法简单且易操作,一般将其写成二阶差分形式,对逆时偏移成像结果进行滤波,对低频噪声去除效果比较明显拉普拉斯算子表达式为其在波数域表示为其中,k应用余弦定理,由波数域矢量计算方法可以得到:式(20)中θ是入射角,k同时,从式(21)也可以看出,在拉普拉斯算子去噪的同时也会损害成像的有用信息,如分子中的频率项应该在成像前的数据上进行补偿.本文采用的是时间域积分方法对成像前数据进行补偿,利用的关系表达式为其中,g(t)是时间域原始信号.3计算值的示例3.1介质模型数值模拟为了进一步验证时空域有限差分相对传统有限差分有更好的数值模拟效果,设计一个大小为201×201(网格数)、网格大小为10m×10m、速度为3000m/s的均匀介质模型.数值模拟时,震源采用一个周期的正弦子波,且其频率为50Hz,时间采样间隔为0.001s,接收的记录长度为2s.震源位于(1000m,1000m)处,分别采用传统有限差分和时空域有限差分,以时间上二阶精度、空间上二十阶精度进行数值模拟,在数值模拟边界处理时,采用完全匹配层(PML)吸收边界条件图5(a,b)分别是采用传统有限差分和时空域有限差分模拟得到的0.3s时刻的波场快照.从快照波前面的形状来看,传统有限差数值分模拟得到的波场快照的波形发生了畸变,有较强频散,而时空域有限差分模拟得到的波场快照的波形保持得比较好,频散相对较弱.图6(a,b)分别是采用传统有限差分和时空域有限差分模拟得到的单炮记录,从(a)、(b)图对比也可以看出,时空域有限差分模拟得到记录的频散更小,波形保持得更好.因此介质模型中的数值模拟结果表明,时空域有限差分数值模拟精度要高于传统有限差分的数值模拟精度3.2数值模拟对比断层模型如图7所示.模型大小为301×201(网格数),网格大小为10m×10m,断层上、下部分的速度分别为1600m/s和2300m/s.接收记录长度为3s,时间步长为0.001s,震源采用一个周期的25Hz的正弦子波,一共58炮,炮点位置在模型90m的深度面从距离最左边100m处开始向右移动,炮间距为50m,每一炮都是全地表等距(10m)接收.在数值模拟和逆时偏移中,都以时间的二阶精度、空间的十二阶精度有限差分求解方程.图8(a,b)是震源在(1500m,90m)时分别采用传统有限差分和时空域有限差分模拟得到的1.5s时刻的波场快照.从图8中对比传统有限差分和时空域有限差分数值模拟的结果表明,时空域有限差分模拟得到的波场的频散更小,精度更高.图9和图10分别是采用传统有限差分和时空域有限差分对断层速度模型进行逆时偏移而得到的剖面.对比去噪前后的偏移剖面,可以看到,高通滤波和拉普拉斯算子去噪都能有效地去除低频噪声,但是高通滤波在去除噪声的同时,也把断层的垂直部分滤掉了,即对陡倾角信号有损害;而拉普拉斯算子去噪能在去除噪声的同时保留好断层的垂直部分,所以拉普拉斯算子去噪在逆时偏移的成像噪声消除中,优于传统高通滤波方法.从图9和图10中也能看到,逆时偏移使水平地层和垂直断层都能较好地成像,且与模型相比,界面成像位置准确,尤其是时空域有限差分逆时偏移得到的最终剖面,断层的垂直部分成像更清楚、更准确.3.3自适应差分算子长度有限差分的计算效率图11是Sigsbee2A偏移速度模型.该模型的特点是拥有较为复杂的盐丘形状以及顶部构造,因此经常被用来验证逆时偏移的效果以及比较噪声去除效果.模型大小为3201×1201(网格数),网格大小为7.62m×7.62m,速度范围为1437~4511m/s.接收记录的时间长度为12s,时间步长为0.0005s,震源采用一个周期的正弦子波,其频率为30Hz.一共200炮,炮点位置是从距离模型最左边1524m处开始向右移动,炮间距为106.68m,每一炮都是全地表等距(7.62m)接收,震源和接收点位于同一深度面15.24m.将求解波动方程的误差控制在10数值模拟和逆时偏移的过程中,基于时空域有限差分方法,分别采用固定差分算子长度和自适应差分算子长度求解波动方程.将自适应差分算子长度应用到偏移的正向延拓和反向延拓计算时,从理论上分析,可以看到自适应差分算子长度平均上要短一些,能有效改进偏移效率.我们在处理平台为ThinkPadR400笔记本电脑(Core2双核,2.0GHz)上做了一个计算效率测试:统计5个炮点正向延拓的耗时,如表1所示,从表中统计可以看出,自适应差分算子长度有限差分的计算效率相对固定差分算子长度有限差分可以节约用时33%左右;统计了5个炮点各自单炮逆时偏移成像整个过程的耗时,如表2所示,从表中也可以看到,在逆时偏移的过程中,自适应差分算子长度有限差分策略能提高效率23%左右.图13(a,b)是分别采用两类算子长度正演模拟得到的波场快照(4s时刻).从图13可以看出,采用两类算子长度正演模拟的结果基本一致.因此采用自适应有限差分法,是一种在能保证精度要求的前提下,提高计算效率的策略图14a是Sigsbee2A速度模型采用时空域自适应有限差分方法逆时偏移的结果.虽然该速度模型对速度场做了平滑,但在盐丘的边界处还是存在着强烈的速度间断,会对逆时偏移带来较强烈的低频噪声,可以看到很强的低频噪声集中在盐丘边界反射处,如图14a所示,因此难以直接分析逆时偏移的效果.图14b是采用高通滤波对逆时偏移结果进行去噪后的结果,从该结果中可以看到,经过高通滤波,低频噪声得到了有效的压制,逆时偏移的效果也显现出来了,但还是没能很干净地去除噪声,损害陡倾角信息.图14c是应用拉普拉斯算子去噪后的逆时偏移剖面,可以看到低频噪声得到了很好的消除,与模型对比,偏移成像结果比较好,成像位置准确,特别是在盐丘边界与盐丘底部成像清晰.由于受观测系统和地下的地质构造的影响,可能会导致在岩丘底部有些区域照明度不够,而要克服此问题,则在偏移过程中需要进行照明度补偿3.4边界条件、自适应差分和自适应差分逆时偏移的应用效果图15是Marmousi速度模型,对其进行重新采样后,模型大小为921×300(网格数),网格大小为10m×10m.该模型的速度范围1500.0~5500.0m/s,总层数超过100层,有大量的反射界面,并含有一些陡倾角、逆冲断层、角度不整合面、地层隆起以及强烈的横向、纵向的速度变化.震源采用一个周期的正弦子波,其频率为30Hz,一共240炮,炮点位置在从距离最左边的2500m处开始向右移动,炮间距为20m,右边放炮,左边接收,每一炮都是120道等距(20m)接收.接收记录长度为3s,时间步长为0.001s.将求解波动方程的误差控制在10分别采用固定差分算子长度和自适应差分算子长度的时空域有限差分法求解波动方程和进行逆时偏移成像.表3是统计两类差分算子长度进行逆时偏移的计算效率,可以看到,自适应差分算子长度有限差分逆时偏移能节时24%左右,有效提高了计算效率.图17a是基于时空域自适应高阶有限差分逆时偏移后的原始结果,由于Marmousi模型存在大量的反射界面,所以在逆时偏移之后,其整个剖面上存在较强的低频噪声.图17b是应用拉普拉斯算子去噪方法对图17a进行处理后的成像剖面,低频噪声得到了比较好的压制.从图17b中也可看出,逆时偏移对模型的主要构造都进行了较好的成像,特别是在中浅部的大倾角构造、断裂及深部的角度不整合地层、目的层(地层隆起)均得到了比较好的成像.4数值模拟结果分析本文采用基于时空域有限差分的差分系数求解声波方程,并分析了其数值频散和稳