可变交错网格优化差分系数法地震波正演模拟

李世中,孙成禹,彭鹏鹏

(中国石油大学(华东)地球科学与技术学院,山东青岛266580)

有限差分方法广泛应用于地震数值模拟中。20世纪70年代,有限差分法首次应用于波动方程正演模拟[1]。之后,VIRIEUX[2]提出了稳定的二阶交错网格弹性波有限差分格式。当用固定步长的有限差分格式对地下低速带或者复杂构造进行数值模拟时,必须选择较小的网格步长来满足计算精度要求或者避免漏掉重要的波场信息,因而不可避免地增加了模拟的计算成本。1989年,MOCZO[3]首次提出采用可变网格思想来提高有限差分数值模拟精度。随后,JASTRAM等[4-5]提出了网格步长可取任意倍数的可变网格算法。李胜军等[6]比较了传统网格与变网格有限差分算法的内存需求、计算效率等,总结了变网格算法的优点。李振春等[7]从理论上推导了变网格界面处虚假反射误差的数学表达式,并通过引入Lanczos滤波算子,实现了一种基于交错网格的稳定且高精度的双变网格正演模拟算法。孙成禹等[8]证实了变网格有限差分算法能提高模拟结果的分辨率,同时降低内存需求量,减少计算时间。孟凡顺等[9]将可变网格有限差分方法引入到交错网格高阶差分数值模拟中,并将可变网格步长倍数变化范围推广至任意奇数倍。KANG等[10]基于一阶速度应力方程,进一步改进有限差分的变时间步长非连续空间交错网格方法。LI等[11]实现了一种灵活的非均匀时间步长不连续曲线网格有限差分方法,用于地震波数值模拟。孙林洁等[12]实现了空间与时间步长任意整数倍变化的双相介质弹性波可变网格正演模拟。马光克等[13]将可变网格有限差分算法引入到逆时偏移中,并分析了可变网格对有限差分算法数值频散的影响。

采用有限差分法进行数值模拟时,由于利用离散网格点的差分来代替空间偏导数,不可避免地产生数值频散现象。许多学者对此做了大量研究,旨在提高模拟的精度。ALFORD等[14]在研究声波方程有限差分模拟精度时指出,为消除数值频散,对于时间2阶、空间4阶精度的有限差分法,在源子波的Nyquist频率所对应的一个波长内不能少于5.5个网格节点。董良国等[15]指出影响数值频散的因素为地震波传播方向、差分精度和一个波长内离散点数。为减少数值频散,可以通过增加有限差分算子的长度去计算空间偏导数[16],但会增加计算成本。TAM等[17]给出了显式频散关系保存方案。LIU等[18]通过使用Taylor级数展开法来最小化时空域频散关系的误差,推导出有限差分系数,以此来减少数值频散。杜启振等[19]通过引入强约束条件和弱约束条件,构造了不同的Lagrange函数,然后通过求取条件极值得到优化差分算子。梁文全等[20]采用前人提出的新的有限差分模板(在保持相同精度的情况下增大了时间步长),基于非线性优化确定时间-空间域隐格式有限差分系数。YANG等[21]基于频散关系和最小二乘理论推导出一阶空间导数的差分系数。YAN等[22]提出一种基于最小二乘理论的优化差分算子进行声波方程正演和逆时偏移。IGEL等[23]使用高斯窗截断得到有限差分算子系数。雍鹏等[24]给出一种显式时间递推格式,并采用共轭梯度法得到精确时间递推匹配系数,实现时空差分算子的同时优化。梁文全等[25]提出使用线性方法压制声波方程矩形网格有限差分算子的数值频散。此外,还有一些关于有限差分算子的优化方法,如牛顿法[26]、模拟退火法[27]、样点逼近法[28]等,都能不同程度地降低较大波数范围内的数值频散,提高地震波数值模拟精度。

以上差分系数优化都是针对固定步长的网格,关于变网格差分系数优化的研究相对较少。本文在前人研究的基础上,提出一种基于最小二乘理论计算可变交错网格优化差分系数算法。通过频散分析和正演模拟算例表明,与Taylor系数展开所获得的差分系数相比,采用优化后的可变交错网格差分系数,能更有效地压制频散现象。

1 方法原理

二维各向同性介质中的一阶应力-速度声波方程可表示为:

(1)

式中:vx,vz分别为质点在x和z方向上的速度分量;t为时间;P为应力;vP为纵波速度;ρ为介质密度。

1.1 变网格差分算法

常规网格采用固定的网格步长对模拟区域进行剖分,而变网格一般采用不同空间步长来对模拟区域进行离散,如图1所示。图1中红色线框区域为小步长网格区域,黑色线框区域为大步长网格区域。当对大小步长网格边界点做插值处理时,一方面会降低数值模拟效率,另一方面也会积累计算误差,影响模拟精度。在大小网格边界处引入过渡带,能避免因插值计算引起的误差,使得数值模拟时,网格长度能自然过渡,增加稳定性。

为了保证过渡带的所有半节点导数不用插值就可以求出高阶差分数值,大小步长网格的纵向步长比需为奇整数倍。以纵向上大小网格步长比3∶1为例,图2给出了可变交错网格在过渡带的差分处理方式[13]。

图1 差分网格分布

图2 过渡带网格变化及波场计算

以一阶导数∂P/∂z的可变交错网格差分形式为例,大网格区域的高阶差分计算公式为:

(2)

小网格区域差分格式:

(3)

网格过渡带差分格式:

(4)

1.2 优化差分算子

对于可变交错网格声波波动方程,一般采用时间二阶、空间高阶有限差分算子来提高数值模拟的精度。对于大小网格过渡带的波场f(x)来说,其空间一阶导数交错网格差分格式可近似地表示为:

(5)

假设f(x)为沿x方向传播的单频波,有:

(6)

(7)

其中β=kh/2,且β∈(0,0.5π)。运用Taylor公式将(7)式的三角函数展开为多项式,可得:

(8)

比较β的系数,得:

(9)

首先对(7)式构建基于频散关系的平方误差函数,有:

(10)

式中:b是积分上限,其取值与K有关。在给定的积分上限取值范围[d,u]内,均匀分布着N个采样点,即b1,b2,…,bN。当K确定时,积分上限b取对应的采样点值bK。

(10)式表示可变交错网格有限差分法求取空间导数时在给定间隔[0,b]内引入的平方误差,由最小二乘理论,求取误差的最小值。此外,加上近零波数约束条件以提高可变交错网格有限差分法的数值模拟精度[21]。对于可变交错网格,有:

(11)

令

(12)

由(10)式和(12)式可得极值的必要条件为:

(13)

式中:λ为Lagrange乘数。将(10)式、(12)式代入(13)式中,最终得到:

(14)

其中,

(15)

1.3 频散分析

根据公式(7),定义以下公式描述可变交错网格有限差分的频散关系:

(16)

当δ(β)的值越接近0时,数值频散就越小;当δ(β)的值越远离0时,数值频散就越大。为了对比可变交错网格优化差分系数和Taylor级数展开所得系数的频散关系,以N=5为例,取积分上限取值范围u=0.4π,将可变交错网格优化差分系数和Taylor级数展开所得系数分别代入(16)式中,得到不同大小网格步长比l下的δ(β)值随β变化的曲线,如图3所示。附录A给出了N=5的可变交错网格优化差分系数和Taylor级数展开的差分系数。

从图3可以看出,相同的可变交错网格差分算子长度下,不同的大小网格步长比l对应的积分上限取值范围d值也不同,l越大,d越小。随着过渡带中小网格内的点逐渐远离大小网格交界处,优化的效果越好。因此,可变交错网格优化差分算子相比基于Taylor级数展开的差分算子能够减少差分近似带来的误差,更好地压制数值频散。下面将对大小网格步长比l=5时的可变交错网格优化差分系数和基于Taylor级数展开的差分系数分别进行数值模拟和对比验证。

图3 大小网格边界区域差分精度分析及对比a l=3,d=0.18π; b l=5,d=0.12π; c l=9,d=0.07π; d l=15,d=0.05π

2 数值模拟

2.1 低速夹层模型

如图4所示的一个低速夹层模型宽为4500m,高为4500m。由地表往下各层的速度分别为3800,2500,4200m/s。第2层为低速层,上、下界面深度分别为2355m和2400m。对模型分别利用可变交错网格优化差分系数法和Taylor级数展开法进行声波方程正演模拟,模型虚线框外为大网格区域,步长15m×15m,虚线框内为小网格区域,步长Δx=15m,Δz=3m。再用常规交错网格(15m×15m,3m×3m)的Taylor级数展开法对低速夹层模型进行声波方程正演模拟。介质密度从上往下依次为1700,1500,1800kg/m3,时间采样间隔均为0.1ms,采用时间2阶、空间10阶有限差分算子进行数值模拟。震源点坐标为(2250m,1800m),子波为雷克子波,主频为40Hz。

图4 低速夹层模型

图5是低速夹层模型正演模拟在0.36s的波场快照,可以看出,当采用固定步长的大网格模拟时,由于低速层的影响,与图5a,图5b和图5d相比,图5c出现非常明显的数值频散。再对比图5a,图5b和图5d 可以看出,即使采用了基于Taylor级数展开的变网格有限差分数值模拟,也出现了较严重的数值频散,干扰了正常波场。而采用可变交错网格优化差分系数法进行数值模拟时,精度更高,压制了频散。

图5 低速夹层模型声波方程波场快照a 可变交错网格Taylor级数展开法(3m/15m); b 可变交错网格优化差分系数法(3m/15m); c 常规大网格Taylor级数展开法(15m×15m); d 常规小网格Taylor级数展开法(3m×3m)

图6是低速夹层模型正演模拟在检波点(2250m,1800m)的部分接收记录,以常规小网格的波形记录为准确值,可以算出可变交错网格优化差分系数法波形记录与准确值的均方根误差为0.002289;而可变交错网格Taylor级数展开法波形记录与准确值的均方根误差为0.008089。均方根误差定义为:

图6 低速夹层模型正演模拟在检波点(2250m,1800m)的部分接收记录

(17)

2.2 起伏界面模型

如图8所示的起伏界面模型宽、高均为4000m,分为4层,第1层速度为1500m/s,第2层速度为2500m/s,第3层速度为2800m/s,第4层速度为3400m/s。第1层为低速层,且第1层和第2层的分界面是起伏的。因此,进行可变交错网格有限差分数值模拟时,令深度600m以上为小网格区域,其步长Δx=10m,Δz=2m,深度600m以下为大网格区域,步长10m×10m。再用常规网格(10m×10m,2m×2m)的Taylor级数展开法对起伏界面模型进行声波方程正演模拟。介质密度设为常数2000kg/m3,时间采样间隔均为0.1ms。震源点坐标为(2000m,10m),子波为雷克子波,主频为24Hz。

图7 低速夹层模型声波方程地震记录a 可变交错网格Taylor级数展开法(3m/15m); b 可变交错网格优化差分系数法(3m/15m); c 常规大网格Taylor级数展开法(15m×15m); d 常规小网格Taylor级数展开法(3m×3m)

图8 起伏界面模型

图9是空间10阶有限差分算子起伏界面模型的正演模拟地震记录,可以明显看出,采用常规大网格进行正演模拟得到的地震记录数值频散最严重,采用可变交错网格Taylor级数展开法得到的地震记录也出现了相当程度的数值频散,不利于实际的波场分析;而在同样条件下采用可变交错网格优化差分系数法得到的地震记录出现的数值频散较少。因此,在采用相同的空间差分算子长度下,可变交错网格优化差分系数法相比Taylor级数展开法能有效降低数值模拟的频散误差。此外,稳定性条件是有限差分数值模拟中十分重要的问题。对于可变交错网格,若是网格步长变化梯度太剧烈,不仅会引起人为的数值干扰,造成计算误差,甚至会使程序计算溢出,无法运行。故空间网格出现太大的步长变化时,往往需要平缓的、经过几次步长变化的处理来保证长时程计算的稳定性。这里长时程的地震记录也间接证明了可变交错网格优化差分系数法计算的稳定性,适用于大区域数值模拟。

图9 空间10阶差分算子起伏界面模型正演模拟地震记录a 可变交错网格Taylor级数展开法(2m/10m); b 可变交错网格优化差分系数法(2m/10m); c 常规大网格Taylor级数展开法(10m×10m); d 常规小网格Taylor级数展开法(2m×2m)

由于过渡带的优化差分算子与空间偏导数算子之间仍存在误差,因此,如果需要进一步降低频散,可以通过提高空间差分阶数来实现,但提高差分阶数将会增大数值模拟计算量,从而增加计算成本。图10是空间18阶差分算子起伏界面模型可变交错网格优化差分系数法地震记录,与图9b相比,数值频散得到了进一步压制。图11是起伏界面模型数值模拟计算成本对比结果,可以看出,当空间差分阶数相同时,可变交错网格的优化差分系数法与可变交错网格Taylor级数展开法计算时间几乎一样,表明用优化差分算子解波动方程不会增加额外的计算成本;而随着空间差分阶数的增大,计算成本也会增加;采用常规小网格做数值模拟,计算成本将急剧增加。

图10 空间18阶差分算子起伏界面模型正演模拟地震记录(积分上限取值范围d=0.12π,u=0.42π)

图11 起伏界面模型数值模拟计算成本对比

3 结束语

在常规交错网格差分系数优化的基础上,推导出可变交错网格的优化差分算子系数,既保证了变网格的计算效率,又能进一步压制数值频散,便于波场特征分析。通过对变网格过渡区域优化差分系数的求取和采用优化系数进行数值模拟,并将模拟结果与Taylor级数展开法的模拟结果对比,可以得到以下认识。

1) 可变交错网格差分算子的精度取决于差分系数和差分阶数,提高差分阶数能压制数值频散,但也导致计算成本增加。与基于Taylor级数展开的可变交错网格差分算子相比,基于最小二乘理论的可变交错网格优化差分算子在不增加计算量的情况下能有效地压制数值频散。

2) 可变交错网格优化差分算子的长度不同时,对应的积分上限取值范围u值也不同。相同算子长度下,不同的大小网格步长比l对应的积分上限取值范围d值也不同,一般来说,当差分算子长度相同时,l越大,d越小。

此外,由于时间步长未采用可变算法,为进一步提高可变交错网格优化差分系数法数值模拟的稳定性,可以采用变空间步长和变时间步长相结合的方法。

参 考 文 献

[1] ALTERMAN Z S,LOEWENTHAL D.Seismic wave in a quarter and three quarter plane[J].Geophysical Journal of the Royal Astronomical Society,1970,20(1):101-126

[2] VIRIEUX J.P-SV wave propagation in heterogeneous media:velocity-stress finite difference method[J].Geophysics,1986,51(4):889-901

[3] MOCZO P.Finite-difference technique for SH-waves in 2-D media using irregular grids:application to the seismic response problem[J].Geophysical Journal International,1989,99(2):321-329

[4] JASTRAM C,BEHLE A.Acoustic modeling on a grid of vertically varying spacing[J].Geophysical Prospecting,1992,40(1):157-169

[5] JASTRAM C,TESSMER E.Elastic modeling on a grid with vertically varying spacing[J].Geophysical Prospecting,1994,42(2):357-370

[6] 李胜军,孙成禹,倪长宽,等.声波方程有限差分数值模拟的变网格步长算法[J].工程地球物理学报,2007,4(3):207-212

LI S J,SUN C Y,NI C K,et al.Acoustic equation numerical modeling on a grid of varying spacing[J].Chinese Journal of Engineering Geophysics,2007,4(3):207-212

[7] 李振春,李庆洋,黄建平,等.一种稳定的高精度双变网格正演模拟与逆时偏移方法[J].石油物探,2014,53(2):127-136

LI Z C,LI Q Y,HUANG J P,et al.A stable and high-precision dual-variable grid forward modeling and reverse time migration method[J].Geophysical Prospecting for Petroleum,2014,53(2):127-136

[8] 孙成禹,李胜军,倪长宽,等.波动方程变网格步长有限差分数值模拟[J].石油物探,2008,47(2):123-128

SUN C Y,LI S J,NI C K,et al.Wave equation numerical modeling by finite difference method with varying grid spacing[J].Geophysical Prospecting for Petroleum,2008,47(2):123-128

[9] 孟凡顺,李洋森,王璐,等.基于可变交错网格的地震波数值模拟[J].地球物理学进展,2013,28(6):2936-2943

MENG F S,LI Y S,WANG L,et al.Based on variable and staggered grids seismic numerical simulation[J].Progress in Geophysics,2013,28(6):2936-2943

[10] KANG T S.Finite-difference seismic simulation combining discontinuous grids with locally variable timesteps[J].Bulletin of the Seismological Society of America,2004,94(1):207-219

[11] LI H,ZHANG W,ZHANG Z,et al.Elastic wave finite-difference simulation using discontinuous curvilinear grid with non-uniform time step:two-dimensional case[J].Geophysical Journal International,2015,202(1):102-118

[12] 孙林洁,刘春成,张世鑫.PML边界条件下孔隙介质弹性波可变网格正演模拟方法研究[J].石油物探,2015,54(6):652-664

SUN L J,LIU C C,ZHANG S X.A variable grid finite-difference scheme with PML boundary condition in elastic wave of porous media[J].Geophysical Prospecting for Petroleum,2015,54(6):652-664

[13] 马光克,李洋森,孙万元,等.可变网格高阶有限差分法逆时偏移研究[J].石油物探,2016,55(5):719-736

MA G K,LI Y S,SUN W Y,et al.Acoustic pre-stack reverse time migration using variable grid finite-difference method[J].Geophysical Prospecting for Petroleum,2016,55(5):719-736

[14] ALFORD R M,KEILY K R,BORE D M.Accuracy of finite difference modeling of the acoustic wave equation[J].Geophysics,1974,39(6):834-842

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

DONG L G,LI P M.Dispersion problem in seismic wave propagation numerical modeling[J].Natural Gas Industry,2004,24(6):53-56

[16] DU Q Z,LI B,HOU B.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered-grid finite difference scheme[J].Applied Geophysics,2009,6(1):42-49

[17] TAM C K W,WEBB J C.Dispersion-relation-preserving finite difference schemes for computational acoustics[J].Journal of Computational Physics,1993,107(2):262-281

[18] LIU Y,SEN M K.Scalar wave equation modeling with time-space domain dispersion-relation-based staggered-grid finite-difference schemes[J].Bulletin of the Seismological Society of America,2011,101(101):141-159

[19] 杜启振,白清云,李宾.横向各向同性介质优化差分系数法地震波场数值模拟[J].石油地球物理勘探,2010,45(2):170-176

DU Q Z,BAI Q Y,LI B.Optimized difference coefficient method seismic wave field numerical simulation in vertical transverse isotropic medium[J].Oil Geophysical Prospecting,2010,45(2):170-176

[20] 梁文全,王彦飞,杨长春.基于优化方法的时间-空间域隐格式有限差分算子确定方法[J].石油物探,2015,54(3):254-259

LIANG W Q,WANG Y F,YANG C C.Determination on the implicit finite-difference operator based on optimization method in time-space domain[J].Geophysical Prospecting for Petroleum,2015,54(3):254-259

[21] YANG L,YAN H,LIU H.Least squares staggered-grid finite-difference for elastic wave modelling[J].Exploration Geophysics,2014,45(4):255-260

[22] YAN H Y,YANG L,LIU H.Acoustic reverse-time migration using optimal staggered-grid finite-difference operator based on least squares[J].Acta Geophysica,2015,63(3):715-734

[23] IGEL H,MORA P,RIOLLET B.Anisotropic wave propagation through finite-difference grids[J].Geophysics,1995,60(4):1203-1216

[24] 雍鹏,黄建平,李振春,等.一种时空域优化的高精度交错网格差分算子与正演模拟[J].地球物理学报,2016,59(11):4223-4233

YONG P,HUANG J P,LI Z C,et al.Optimized staggered-grid finite-difference method in time-space domain based on exact time evolution schemes[J].Chinese Journal of Geophysics,2016,59(11):4223-4233

[25] 梁文全,王彦飞,杨长春.声波方程数值模拟矩形网格有限差分系数确定法[J].石油地球物理勘探,2017,52(1):56-62

LIANG W Q,WANG Y F,YANG C C.Acoustic wave equation modeling with rectangle grid finite difference operator and its linear time space domain solution[J].Oil Geophysical Prospecting,2017,52(1):56-62

[26] KINDELAN M,KAMEL A,SGUAZZERO P.On the construction and efficiency of staggered numerical differentiators for the wave equation[J].Geophysics,1990,55(1):107-110

[27] ZHANG J,YAO Z.Optimized explicit finite-difference schemes for spatial derivatives using maximum norm[J].Journal of Computational Physics,2013,250(10):511-526

[28] YANG L,YAN H,LIU H.Optimal implicit staggered-grid finite-difference schemes based on the sampling approximation method for seismic modeling[J].Geophysical Prospecting,2016,64(3):595-610

附录A 可变交错网格过渡带差分系数

表A-2 大小网格步长比l=5的可变交错网格差分系数(d=0.12π,u=0.4π)

表A-3 大小网格步长比l=9的可变交错网格差分系数(d=0.07π,u=0.4π)

表A-4 大小网格步长比l=15的可变交错网格差分系数(d=0.05π,u=0.4π)