基于紧致交错差分格式的二维声波及黏滞声波方程数值模拟

邓光校, 鲍羽, 汪勇*

(1.中国地质大学(武汉)地球物理与空间信息学院, 武汉 430074; 2.中国石化碳酸盐岩缝洞型油藏提高采收率重点实验室,乌鲁木齐 830011; 3.油气资源与勘探技术教育部重点实验室(长江大学), 武汉 430100)

地震波场正演模拟是针对地下介质运用数值计算方法去模拟地震波在地下的传播过程,进而得到位于地面或地下观测点的地震记录,是地震勘探以及地震学的重要基础。基于有限差分数值模拟的发展以及实际生产的需求,针对如何去提高有限差分的计算时的效率[1]、模拟精度[2]、算法稳定性[3-5]、处理复杂介质[6-7]、吸收边界条件[8-9]和压制数值频散[10-11]等研究方向上,研究人员陆续提出了多种方法,也做了大量研究,也取得较为重大的研究成果。

差分计算时增加网格点数是提高数值,模拟精度最直接的方法,但网格点增加必然导致计算时计算量以及存储量的增加。然而紧致差分就可以很好地解决这个问题,紧致差分是一种隐式差分格式,稳定性也较好,也正是具备这些优势也使它成为目前研究较多的有限差分方法之一。Du等[12]为了降低频散将交错网格与紧致有限差分结合,并将其应用于垂直横向各向同性(vertical transversely isotropic, VTI)介质的一阶速度应力方程的数值模拟。Liu等[13-14]证明了在不考虑访问数组的额外成本时,高阶显式差分能够被相同阶的隐式差分所替代。并且将其运用到了声波以及弹性波的数值模拟进行了对比。杨宽德等[15]把紧致差分运用于包含流体的多孔隙各向异性介质中,并且发现通量校正传输(flux-corrected transport,FCT)紧致差分方法能有效压制粗网格条件下模拟弹性波场所引起的数值频散以及源噪声。周成尧等[16]基于五点八阶超紧致差分格式在黏弹性介质中实现了声波方程的正演模拟。紧致差分格式后也逐渐被研究人员依次应用于声波、弹性波和复杂介质等地震波场数值模拟中。交错网格差分格式最早是由Madariaga[17]提出的,交错差分提高了数模模拟时的局部精度,也使得收敛速度更快。因此将交错网格技术引入紧致差分,可以进一步提高数值模拟的精度,更好地压制数值频散。

由于地层的黏弹性所导致的地层对地震波产生的吸收、衰减作用,使得地震波在地下的传播规律变得复杂,所以研究地震波在黏弹性介质中的传播规律对地震勘探有着重要的意义。研究人员也基于不同方法对不同的黏弹介质模型进行了数值模拟研究,如李晓波等[18]基于斑状饱和介质的黏弹特性进行了地震波模拟;罗文山等[19]针对黏滞波动方程推导了分数阶拉普拉斯算子黏滞波动方程的一阶速度-压力形式,在时间域进行了正演模拟。姚振岸等[20-21]在黏弹各向异性介质中进行了微地震波场模拟。

关于声波方程在紧致交错差分格式下的稳定性条件以及该差分格式运用于黏滞声波方程的数值模拟尚未有相关方向的文献发表。现分析和对比它与传统的中心差分和交错差分格式,在模拟精度、数值频散和稳定性条件3个方面的区别,并将其应用于对一阶速度-应力声波方程组和黏滞声波方程组的数值模拟。

1 二维黏滞声波方程

在研究介质的非完全弹性时,地震勘探中常采用开尔文-佛各特(Kelvin-Voigt)体的模型,认为地震波吸收系数与地震信号频率的平方成正比,这也与实际的地质情况更加贴近。该介质中黏滞声波方程的一阶速度-应力波动方程组为

(1)

2 紧致交错有限差分方法

Madariaga[17]最早提出波动方程交错网格数值模拟方法,差分精度为O(Δt2+Δx2),进而Levander[22]将交错网格方法应用于求解P-SV波方程,差分精度也被提高到O(Δt2+Δx4),董良国等[4]进一步提出了一阶弹性波方程的交错网格高阶差分解法,使得差分精度达到O(Δt2M+Δx2N)。现如今,交错有限差分方法已经发展成一种高效且被广泛应用的数值模拟方法。

假设函数f(x)连续,则一阶导数的2N阶精度交错网格差分可以表示为

(2)

式(2)中:fi为节点函数值;f′i为一阶导数值;Δx为网格大小;cn为差分系数,如表1所示。

表1 一阶导数的差分系数表

对式(1)中的vx(t)分别在t+Δt/2时刻和t-Δt/2时刻通过泰勒公式展开相减,略去高次项,可得vx(t)的时间二阶精度近似式为

(3)

代入黏滞声波方程[式(1)],将vx对时间的导数转换为P对空间的导数,可以得到vx的二阶精度两层显示差分格式,同理也可以得到vz和P的时间二阶精度递推式为

(4)

利用式(4)就可以进行黏滞声波方程的地震波场时间层的推进。

由式(2)可以得到,常规交错差分格式达到2N阶空间差分精度,需要用到2N个节点的函数值。如果想要提高差分格式数值的计算精度,最直接有效的方法就是在差分计算时增加网格点的数目,但这也使得计算量以及所需的存储空间增加。紧致差分格式很好地解决了这个问题,早期的紧致差分格式是建立在Hermite多项式构造的基础之上。

将紧致差分格式与交错网格技术结合,最早由Nagarajan等[23]提出,并且将其用于大涡模拟(large eddy simulation)问题,紧致交错差分仅用4个网格节点得到6阶空间精度的模拟效果。随后,Boersma[24]提出了最高到12阶空间精度的紧致交错差分格式,用于可压缩流体的Navier-Stokes方程的数值模拟。

若函数f(x)是连续的,则一阶导数的2N(N=2～5)阶精度紧致交错差分格式表示为

(5)

式(5)中:a、bn为差分系数,如表1所示。

从式(2)和式(5)可以看出,二者均是利用半网格点上的变量值计算整网格点上的一阶导。区别在于,常规交错差分格式(staggered finite difference,SD)在求取中心点处的导数值时,仅需要知道周围网格点的函数值,而紧致交错差分格式(staggered compact finite difference,SCD)还需使用相邻点的导数值,所以可以认为SD格式是SCD格式的特例,即式(5)中的a=0。由于紧致差分格式是空间隐式差分,若要达到2N阶空间差分精度,仅需2N-2个节点函数值,能够节省计算量和存储空间。

根据式(5)所表示的SCD格式可以求得式(4)中的各偏导数的值。对二维地震波场离散化后假设在x方向有n个节点,z方向有m个节点。紧致交错差分格式在进行一阶导数的求取时,是把变量的值和导数值分别置于两套网格中,因此可以把差分格式定义成向前和向后差分两种类型。以变量P为例,将导数值固定在整网格上,变量P的值固定在半网格上。然后根据式(5)求∂P/∂x和∂P/∂z的值,则差分格式可以写成矩阵形式为

(6)

式(6)中:

向前差分时:

向后差分时:

根据式(6),由应力场P求它的空间一阶导数∂P/∂x和∂P/∂z可表示为

(7)

表1给出了矩阵中的差分系数a和bn(n=1～4),可以求得4～10阶的空间差分精度近似值。

紧致交错差分格式与常规交错差分格式定义变量类似,利用相似的方法,将应力P固定于时间整网格上,速度vx和vz固定于时间半网格,可得到时间为二阶精度的一阶速度-应力黏滞声波方程的紧致交错差分格式表示为

(8)

为了说明SD和SCD格式的区别,这里直接给出一阶速度-应力黏滞声波方程交错网格时间二阶精度差分格式为

(9)

式(9)中:E1和E2分别为x方向的向前和向后差分系数矩阵;F1和F2分别为z方向的向前和向后差分系数矩阵。E1和E2的表达式如下,其中的ci(i=1～5)为表1中的差分系数。

2.1 精度分析

当时间差分精度相同时进行数值模拟,不论是利用常规交错网格差分(SD)格式,还是利用紧致交错网格差分(SCD)格式,它们在时间层的推进方式是相同的。唯一区别在于SD和SCD空间偏导数求取时的差分格式不同,例如SCD格式表示的向前差分为∂P/∂x=PB1A-1/Δx,SD格式则为∂P/∂x=PE1/Δx。在进行二者近似精度的比较时,需要获取它们一阶导数的截断误差去进行比较,如表2所示。不难看出,在同样空间近似精度的条件下,两种格式的截断误差有较大的差别。SD格式计算一阶偏导数4～10阶差分精度的截断误差约是SCD格式的1.46、3.18、9.71和8.7倍,数据表明SCD与SD相比,其截断误差更小和差分精度也更高。

表2 一阶导数SD和SCD方法的截断误差主项系数比较

为了比较SD、SCD格式和常规中心差分格式(finite difference,FD)的数值模拟精度,以一维平面简谐波初值问题进行试算。一维平面谐波初值问题数学表达式为

(10)

上述偏微分方程的达朗贝尔解,即精确的解析解为

(11)

设置一维介质模型长度200 m,v=50 m/s,时间步长为0.1 ms,空间网格大小设置为1 m。3种差分格式的空间差分固定为四阶精度进行数值模拟。可得图1所示的1 s时刻的∂u/∂x分量波场快照,图1中黑色的实线部分是套用一维平面简谐波的达朗贝尔解计算得到的精确解析解。不难看出,紧致交错差分的数值模拟结果更加接近精确的解析解,说明其具有更好的模拟精度,常规交错网格次之,有限差分格式的精度最差。由定量计算可以得知,该时刻FD、SD和SCD的数值解与解析解之间的相对误差依次为24%,8.5%和5.2%,导致这一误差的结果正是由于它们在计算空间导数时存在不同截断误差大小,这也与表2的理论分析结果相符。需要说明的是,这里为了比较精度差异,选取的网格较大,导致波形方波化,若减小网格大小,能够使波形足够光滑,模拟精度也更高。

图1 3种差分格式的一维声波方程数值模拟

2.2 频散分析

如果在进行数值模拟的时候采用的空间网格过大,就会得到较大的求解误差,产生数值频散,因此在进行判别数值模拟方法的优劣时,频散分析是一种重要的分析手段,也是进行数值模拟确定网格大小的重要依据。利用一阶导数的紧致交错差分格式[式(5)]进行频散分析,依据数值波数与真波数的比值,分析该方法的适用条件。令

(12)

将式(12)代入式(5)中可得

(13)

利用欧拉公式可得

(14)

将B=(Ik′)A代入式(14)可得

(15)

假设φ=kh,φ′=k′h,那么修正波数与真波数之比就可以定义为

(16)

在理想情况下,如果不存在数值频散则波数比R恒等于1。R偏离1越大,则说明该方法的数值频散越严重,反之则说明该方法能更好地压制数值频散。同样的方法可以得到一阶导数的常规交错差分(SD)和中心差分(FD)格式的波数比,图2为3种格式在不同差分精度情况下的波数比曲线。

图2 3种差分格式的波数比曲线

将φ∈[0,π] 作为坐标系的横轴,可由波数与空间步长相乘得到,单位波长中采样点个数N=2π/φ,因此横坐标亦可以看作N由∞逐渐减小至2。根据图2中的速度比曲线不难得出,交错差分以及紧致交错差分的数值波数相较于常规的中心差分更接近于真波数,直接表明这两种方法在压制数值频散上的效果更好;紧致交错差分格式与常规交错差分格式相比,压制数值频散的特征更加明显,例如图2中6阶SCD曲线接近10阶SD,所以SCD方法能使用更少的网格点,可以节省计算资源;3种方法均在10阶精度的前提下,为保证无数值频散,要求波数比大于99.99%,则FD、SD和SCD方法每个波长内需要采样11.6、5.2和4.2个点。所以紧致交错差分更能适用于粗网格时的数值模拟,能够更好地压制数值频散。

前人研究表明FD格式不利于压制数值频散,这里只比较SD和SCD两种差分方法的差异。模型大小设置为4 000 m×4 000 m,纵波速度给定为4 000 m/s,在模型的中心加载雷克子波震源,峰值频率25 Hz。时间步长2 ms,空间网格20 m,再针对一阶速度-应力声波方程采用SD和SCD格式在不同空间差分精度进行波场模拟,进而得到同一时刻的应力P分量的波场快照(若未加说明,波场快照和地震记录均为应力P分量),如图3所示。

图3 SD和SCD不同差分精度计算的400 ms时刻波场快照

从图3可以看出,同前面的理论分析结果一致,SCD格式在压制数值频散方面比SD格式更具有优势,图3(d)为6阶SCD格式,其波场快照与图3(c)所示的10阶SD格式接近,要远远好于图3(a)所示的6阶SD格式。为了更清楚地观察差异,提取模型x=2 000 m,z=2 000 m处质点的振动曲线,见图4。可以看出,6阶SD的振动曲线存在明显的尾波,即数值频散,而10阶SD、6阶和10阶SCD的振动曲线则基本重合,波形平滑,无尾波现象,这表明SCD格式能用较低的差分阶数,能够使用粗网格计算,提高了运算效率。

图4 SD和SCD不同差分精度计算的波场记录

2.3 稳定性分析

有限差分数值模拟中另外一个非常重要的问题就是稳定性条件,也是影响差分方法计算效率的重要因素之一。这里参照董良国等[4]和杜启振等[5]使用的Fourier方法对差分格式[式(8)]进行稳定性分析,省略掉推导过程,直接利用式(17)表示的二维声波方程(2M,2N)阶差分精度的紧致交错网格差分格式的稳定性条件。

(17)

定义α=vΔt/h为差分格式的Courant数,并与常规交错网格差分格式的稳定性条件进行对比,对比结果如表3所示。从表3来看,差分格式的时间差分精度越高,差分格式也越稳定,但是更高阶的时间精度在计算时需要转变为更高阶的空间导数,会相应地增加计算量。当差分精度相同时,SCD格式的稳定性条件比SD格式要略微严格,也就是说空间步长相同时,允许采纳的时间步长要略小。

表3 二维声波方程不同差分格式的Courant数

2.4 PML边界条件

数值模拟过程模型大小的限制在计算时的设置都会导致人工边界的存在,倘若不对其进行处理,就会产生边界反射,从而对正常的地震波场造成干扰,因此处理人工边界的好坏也直接影响到了数值模拟的精度和效率。对于边界的处理,采用完全匹配层(perfectly matched layer,PML),按照王守东[25]和刘有山等[26]推导声波PML控制方程的思路和方法,略去推导过程,这里直接给出一阶速度-应力黏滞声波方程的PML边界条件的控制方程组,当η=0时,即可得到声波方程的PML控制方程,即

(18)

式(18)中:P1和P2为引入的中间变量;d(x)和d(z)分别为x和z方向的衰减系数,在这两个方向上起到衰减的作用。高刚等[27]的研究指出,d(x)和d(z)的衰减函数类型和层数会影响边界反射的衰减效果,认为余弦型吸收衰减函数效果较好,所以在模型试算中采用该衰减函数,表示为

(19)

式(19)中:αi为边界层的吸收衰减因子;B为衰减幅度因子,计算时取500;PML为整个吸收边界层的厚度,在进行实际的数值模拟时,可令其为20。为了验证所使用差分格式在PML条件下对边界反射的吸收效果,在2.2节中的模型之上采用SCD(2,10)(时间2阶,空间10阶精度)格式进行数值模拟,空间网格10 m,时间步长1 ms,对比未采用PML边界条件和采用PML边界条件时在600 ms时刻的波场快照,如图5(a)和图5(b)所示,并将二者做差,得到图5(c)只包含边界反射的波场快照。不难看出,SCD格式在PML控制方程作用下,人工边界反射得到了有效吸收,并且有效波没有发生明显改变,处理结果十分可靠。

图5 600 ms时刻波场快照

3 模型试算

第2节通过理论分析和数值模拟讨论了二维声波方程的紧致交错差分格式的计算精度、频散关系、稳定性和边界条件,然后基于SCD(2,10)差分格式对式(1)表示的二维黏滞声波方程进行试算,以验证该方法是否适用于黏滞声波方程的数值模拟。

3.1 均匀介质模型

图6 350 ms时刻波场快照

提取数值模拟时得到的波场快照中x=2 000 m的波形记录,如图7所示。不难看出,波形曲线平滑,并且不存在数值频散,且随着品质因子的减小,能看到地震波振幅发生显著衰减,频率降低,波形展宽,相应的波长也增加了。

图7 不同品质因子时的波形曲线

3.2 水平层状介质模型

设置4层层状介质模型,利用SCD(2,10)格式进行数值模拟,模型参数如表4所示。

表4 多层层状介质模型参数

模型深度和长度均设置成2 000 m,时间步长Δt=0.5 ms,空间网格大小h=10 m,在地表(x=1 000 m,z=0 m)位置加载主频30 Hz雷克子波震源,由地表接收的单炮记录(记录长度1.5 s)如图8所示,为了方便比较显示,对地震记录进行了增益处理,提高深层反射波振幅强度。

图8 地面地震记录

利用四层水平层状介质模型对地面接收到的地震波场进行模拟,发现地震记录清晰可见,不存在数值频散和边界反射情况,反射波显示清晰,充分说明本文算法可以有效地模拟黏滞声波在多层介质模型中的传播。

图9给出检波点(x=1 500 m,z=0 m)接收到的振动曲线,可以看出,由于介质的吸收衰减作用,黏滞声波中的反射波振幅比声波模拟要小,且随着传播时间和距离的增加,反射波信号变宽,频率变低,这也更加符合实际地震信号特征。

图9 检波器(x=800 m,z=0 m)接收到的波形曲线

3.3 Marmousi模型

为了检验紧致交错差分格式是否适用于复杂的介质的数值模拟,采用经典的二维Marmousi纵波速度模型进行模型试算。模型如图10所示,当速度v的范围是1 729～5 500 m/s,品质因子Q通过经验公式Q=14v2.2计算得到,其中v的单位是km/s,可以算出品质因子范围为46.7～595.6。模型大小设置为501×501个网格点,空间网格大小固定为8 m,时间步长0.5 ms,纵波震源位于(x=2 000 m,z=0 m),加载30 Hz的Ricker子波,采样时长3 s。

图10 Marmousi模型

用Marmousi模型对黏滞声波方程和声波方程进行数值模拟得到了如图11所示的波场快照,可以看到,波场特征清晰,无明显数值频散以及边界反射,边界的吸收效果良好。阐明了本文方法对复杂介质的适用性。图12分别是由声波方程以及黏滞声波方程得到的地震记录,由于波前扩散导致深部的反射振幅明显减弱,因此对地震记录进行了自动增益控制处理,加强了深层的能量使显示更为清晰,时窗设置为500 ms,通过对比可以发现:黏滞声波的数值模拟结果除了波前扩散造成的振幅衰减,地震波传播过程中还受到了介质的吸收衰减作用,从而反射波的振幅衰减幅度要大于声波方程中的模拟结果,而且地震波的频率也会随之变化,频率会随着地震波的传播距离的增加而减小。

图11 Marmousi模型波场快照及黏滞声波模拟

图12 Marmousi 模型的地面地震记录(z=0 m)

对于进行AGC处理后的地面地震记录,提取单道的地震记录波形图(x=1 000 m)如图13所示。可以看到,黏滞声波记录的地震波振幅明显小于声波,子波的延续时间大于声波,反射波个数也更少,并且不难发现,传播的时间越长这种差异越明显,进而导致了黏滞声波模拟结果从垂向分辨率上小于声波。

图13 检波器(x=1 000 m)接收的地震记录

对模拟记录进行广义S变化,再进行时频分析,得到的结果如图14所示,可以看到两道模拟记录在频谱上的差别。从结果来看,声波的模拟记录[图14(a)]的地震波振幅随着时间的增加而减小,导致这种变化的原因是波前扩散。但是振幅虽然衰减了,反射波的主频却未发生明显改变,几乎与激发主频平直,也是30 Hz左右。而黏滞声波记录[图14(b)]除了振幅衰减,反射波主频也会降低到10 Hz左右。通过对比发现,两个记录的时频谱,声波的数值模拟结果其能量的峰值相对集中也更加清晰,与时间域的主要反射波对应;而黏滞声波记录中,1.4 s和2.5 s位置的能量峰值相对模糊一些,这也从频域说明了地震垂向分辨率的降低。

图14 x=1 000 m处单道地震记录时频分析结果

4 结论与建议

利用声波方程结合紧致交错差分格式,在模拟精度、频散曲线和稳定性分析3个方面,研究认为SCD格式仅需使用2N-2个节点,就能使一阶导数的差分精度达到2N阶精度,少于FD格式(2N+1个)和SD格式(2N个);与SD格式相比,SCD格式具有更高的模拟精度和略微严格的稳定性条件;此外,6阶SCD格式的频散关系几乎可以等同于10阶SD格式,这也说明相同精度下的SCD格式其数值频散更小,更加适用于粗网格下的大尺度的地震波场数值模拟。

其次,基于黏滞声波的PML控制方程,在均匀介质模型、水平层状介质模型和Marmousi模型上利用SCD格式,进行了数值模拟。最终可见地震波场特征清晰,表明本文方法对于复杂介质的适用性,也进一步说明该差分格式不仅可以用于黏滞声波的波场模拟,并且还能够推广到二维或三维的各向异性介质以及双相介质等复杂介质的声波或弹性波数值模拟中,同时为研究地震波传播规律、逆时偏移以及全波形反演等工作提供了一种行之有效的方法。

为减少计算内存和时间,所用的差分格式时间为二阶,在数值模拟时可以通过提高时间差分的阶数,从而达到增加数值模拟计算精度的目的。此外,也可以通过优化差分系数的方法进一步提升该方法的优势。