基于时空域交错网格有限差分法的应力速度声波方程数值模拟

2022-01-28彭更新郭念民胡自多徐凯驰裴广平

石油物探 2022年1期

彭更新,刘威,郭念民,胡自多,徐凯驰,裴广平

(1.中国石油天然气集团公司塔里木油田分公司,新疆库尔勒841000;2.中国石油勘探开发研究院西北分院,甘肃兰州730020)

有限差分法具有占用内存小和编程实现简单等优点,是应用最广泛的一种波动方程数值模拟方法[1-3]。有限差分法利用差分算子近似表示波动方程中的时间和空间偏微分算子,进而通过迭代求解差分离散波动方程实现波动方程数值模拟。差分算子近似微分算子,会导致地震波的传播速度与真实速度不相等,并且不同频率成分的传播速度不相同,即出现数值频散现象[4-5]。固有的数值频散严重影响有限差分法的数值模拟精度;同时,相比基于射线理论[6]和高斯束理论[7]的数值模拟方法,有限差分法的计算效率相对较低。因此,提高模拟精度,同时保持甚至提高计算效率是有限差分数值模拟方法面临的重要挑战。

波动方程数值模拟广泛采用的有限差分法可以分成两类,即常规网格有限差分法和交错网格有限差分法。常规网格有限差分法主要应用于二阶标量波动方程数值模拟[8-10];交错网格有限差分法主要应用于一阶应力-速度声波[11-13]或弹性波动方程数值模拟[14-16]。相比常规网格差分法,交错网格差分法通常具有更高的模拟精度,并且对速度对比度(介质最大速度与最小速度之比)较大的介质数值模拟稳定性更高[17],同时交错网格差分法更便于施加应力和位移震源。

KANE[18]首次将交错网格有限差分应用于电磁波方程数值模拟,VIRIEUX[14-15]将时间二阶和空间二阶交错网格有限差分法应用于地震波场的数值模拟。为了提高交错网格差分法的模拟精度,LEVANDER[19]在VIRIEUX[14-15]研究成果的基础上,进一步发展了时间二阶和空间四阶交错网格差分法。通过增大空间差分算子长度构建出空间任意偶数阶交错网格差分法,相应的差分系数通常基于空间域频散关系和泰勒展开求解[20]。本文将这种时间二阶、空间2M阶并基于空间域频散关系和泰勒展开计算差分系数的交错网格差分法称为空间域交错网格有限差分法(SD-SFD)。差分离散波动方程在时空域同时迭代求解,而空间域交错网格有限差分法仅基于空间域频散关系计算差分系数,导致空间域交错网格有限差分法容易出现数值频散,模拟精度低。

本文旨在改进空间域交错网格有限差分法的差分系数算法,进而提高模拟精度。在保持时间二阶、空间2M阶差分格式不变的情况下,提出利用时空域频散关系和泰勒展开计算差分系数,称这种改进差分系数算法的交错网格差分法为时空域交错网格有限差分法(TSD-SFD)。首先,给出时间二阶、空间2M阶交错网格差分法对二维应力-速度声波方程的差分离散方程;其次,推导基于时空域频散关系和泰勒展开的差分系数算法;再次,分析和对比时空域交错网格有限差分法和空间域交错网格有限差分法的数值频散特性和稳定性特性;最后,利用层状介质模型和塔里木盆地典型复杂构造模型进行数值模拟,对比两种交错网格差分法的模拟精度和计算效率。

1 交错网格有限差分法

二维应力速度声波方程可以表示为:

(1)

式中:P=P(x,z,t)为压力场,u=u(x,z,t)和w=w(x,z,t)分别为质点震动速度场的x和z分量,κ为体积模量,ρ为介质密度。

应力速度声波方程数值模拟普遍采用交错网格有限差分法。交错网格有限差分法的基本思想是将压力场P、质点震动速度场u和w定义在交错的网格系统中,如图1所示。

图1 二维交错网格差分法示意

交错网格有限差分法通常采用时间二阶、空间2M阶差分法。对波动方程中的一阶时间偏导数∂P/∂t、∂u/∂t和∂w/∂t采用二阶差分近似,可以得到:

(2)

(3)

(4)

对波动方程中的一阶空间偏导数采用2M阶差分近似,可以得到:

(5)

(6)

(7)

(8)

式中:am(m=1,2,…,M)为差分系数。

将方程(2)～(8)代入方程(1),可以得到:

(9)

方程(9)为时间二阶、空间2M阶交错网格差分法对应力速度声波方程的差分离散波动方程。空间域交错网格有限差分法和时空域交错网格有限差分法均通过迭代求解方程(9)实现波动方程数值模拟,只是二者采用的差分系数算法不同。

2 差分系数计算

差分系数直接影响交错网格有限差分法的数值模拟精度和稳定性,因此,差分系数算法是交错网格有限差分法的重要研究内容。

空间域交错网格有限差分法仅在空间域计算差分系数,相应的差分系数算法可概括为:首先,根据平面波理论和空间差分算子表达式推导空间域频散关系;其次,对空间域频散关系中的三角函数进行泰勒级数展开,建立差分系数求解方程组,再通过解方程组计算差分系数。

时空域交错网格有限差分法在时空域计算差分系数,相应的差分系数算法可概括为:首先,根据平面波理论和差分离散波动方程推导时空域频散关系;其次,对时空域频散关系中的三角函数进行泰勒级数展开,建立差分系数求解方程组,通过解方程组计算差分系数。

空间域交错网格有限差分法和时空域交错网格有限差分法的差分系数计算过程均采用了平面波理论。在均匀介质中,应力速度声波方程具有如下形式的离散平面波解:

(10)

(11)

(12)

kx=kcosθkz=ksinθ

(13)

式中:AP、Au和Aw为平面波的振幅,k为波数,θ为平面波传播方向与x轴的夹角,ω为角频率。

2.1 空间域差分系数算法

空间域交错网格有限差分法仅在空间域计算差分系数,将离散平面波解方程(10)代入空间差分算子方程(5)得到:

(14)

方程(14)为空间2M阶交错网格差分法的空间域频散关系,泰勒展开其中的正弦函数得到:

(15)

(16)

解方程(16)得到:

(17)

m=1,2,…,M

方程(17)为空间域交错网格有限差分法的差分系数通解,可以看出其差分系数仅与M的取值有关,而与地震波的传播速度无关。

2.2 时空域差分系数算法

时空域交错网格有限差分法在时空域计算差分系数,将平面波解方程(10)～(13)代入差分离散波动方程(9)得到:

(18a)

(18b)

(18c)

消除方程18a～18c中的AP,Au和Aw,并且考虑到ω=υk,κ=ρυ2得到:

(19)

式中:r=υτ/h为Courant条件数,表示单位时间步长内地震波的传播距离与空间采样步长之比。

方程(19)为时间二阶和空间2M阶交错网格差分法的时空域频散关系。利用时空域频散关系和泰勒展开求解差分系数,需要选择一个特定θ值。为了尽可能简化差分系数计算过程,取θ=0或θ=π/2,可以得到:

(20)

泰勒展开式(20)中的正弦函数得到:

(21)

使方程(21)左右两边k2j-1h2j-2的系数对应相等得到:

(22)

j=1,2,…,M

解方程(22)得到:

(23)

m=1,2,…,M

方程(23)为时空域交错网格有限差分法的差分系数通解。取r=0时,方程(23)与方程(17)完全等价,因此,空间域交错网格有限差分法的差分系数是时空域交错网格有限差分法的差分系数的一个特解。同时,方程(23)表明时空域交错网格有限差分法的差分系数不仅与M的取值相关,还与r=υτ/h取值相关,数值模拟过程中时间步长τ和空间步长h取值固定,差分系数随速度υ自适应变化,这是时空域交错网格有限差分法比空间域交错网格有限差分法具有更高模拟精度的根本原因。

3 数值频散分析

数值频散的严重程度直接影响交错网格有限差分法的数值模拟精度。本文采用归一化相速度误差εph(θ)度量数值频散的大小。根据时间二阶、空间2M阶交错网格差分法的时空域频散关系方程(19),可导出εph(θ)的表达式为:

(24)

(25)

当εph(θ)=0时,相速度等于真实速度,无数值频散;εph(θ)>0时,相速度大于真实速度,表现出时间频散;εph(θ)<0时,相速度小于真实速度,表现出空间频散。

空间域交错网格有限差分法和时空域交错网格有限差分法的归一化相速度误差εph(θ)的表达式完全相同,均由方程(24)和(25)表示,但是二者的差分系数am(m=1,2,…,M)不同,从而表现出不同的数值频散特性。

根据εph(θ)的方程(24)和方程(25),结合两种差分法各自的差分系数,可绘制两种差分法的相速度频散曲线,进而分析数值频散特性。

图2给出了r=0.2和r=0.3时,空间域交错网格有限差分法(M=10)和时空域交错网格有限差分法(M=10)的相速度数值频散曲线,分析对比可以看出:

图2 相速度数值频散曲线

1)r取值增大,空间域交错网格有限差分法(M=10)和时空域交错网格有限差分法(M=10)的相速度数值频散均明显增大;

2)在波动方程有限差分法数值模拟中,通常保证单位波长内有4个采样点,对应kh/π的取值为0.5,在0≤kh/π≤0.5范围内,r取值相同时,时空域交错网格有限差分法(M=10)的相速度数值频散明显小于空间域交错网格有限差分法(M=10),因此,时空域交错网格有限差分法(M=10)压制数值频散的效果更好,模拟精度更高;

3)在0≤kh/π≤0.5范围内,空间域交错网格有限差分法(M=10,r=0.2)和时空域交错网格有限差分法(M=10,r=0.3)的相速度数值频散大小基本相当,但时空域交错网格有限差分法(M=10)采用的r取值更大,意味着能采用更大的时间步长,从而提高计算效率,同时不损失模拟精度;

4)当kh/π≥0.6时,时空域交错网格有限差分法(M=10)存在一定的空间数值频散,这时相比空间域交错网格有限差分法在压制数值频散,提高模拟精度方面的优势会降低。

4 稳定性分析

根据时间二阶、空间2M阶交错网格差分法的时空域频散方程得到:

(26)

取空间波数kx=kz=π/h,并且根据三角函数有界性sin2(rkh/2)≤1,可得到:

(27)

式中:S为稳定性因子。

空间域交错网格有限差分法和时空域交错网格有限差分法的稳定性条件均由方程(27)所描述,但两种方法的差分系数不同。

图3a和图3b分别给出了两种差分法的稳定性因子S随r取值变化的曲线。空间域的差分系数与r取值无关,其稳定性因子S不随r变化,呈现为一条平行于r坐标轴的直线;时空域的差分系数为r的函数,其稳定性因子S表现为随r取值变化的曲线。

根据稳定性条件不等式r≤S,稳定性因子S曲线(直线)与直线S=r的交点为稳定性条件约束下的最大r取值。图3c给出了两种差分法稳定性条件不等式约束下的最大r取值随M的变化曲线,称为稳定性曲线。稳定性条件不等式约束下的最大r取值越大,稳定性越强,反之,则稳定性越弱。图3c表明,空间域交错网格有限差分法和时空域交错网格有限差分法的稳定性均随M增大而降低;当M≥1时,时空域交错网格有限差分法的稳定性明显强于空间域交错网格有限差分法,因此利用时空域交错网格有限差分法进行数值模拟时,能采用更大的r取值,即采用更大的时间步长,从而提高计算效率。

图3 稳定性曲线

5 数值模拟实例

利用层状介质模型和塔里木盆地典型复杂构造模型开展数值模拟实验,进一步分析、对比两种方法的数值模拟精度和计算效率。

5.1 层状介质模型

图4给出了一个网格尺寸为Δx=Δz=h=15 m,网格数为nz×nx=601×801的层状介质模型。主频为25 Hz的Ricker子波作为震源,位于网格点(401,3)。空间域交错网格有限差分法(M=10)和时空域交错网格有限差分法(M=10)分别采用时间步长τ=1.5 ms和τ=1.0 ms进行数值模拟。图5和图6分别给出了利用两种差分法(M=10)进行数值模拟得到的2.1 s时刻压力场P和质点震动速度场x分量u的波场快照。

图4 层状介质速度模型

图5 层状介质模型数值模拟2.1 s时刻波场快照(时间步长τ=1.5 ms)

图6 层状介质模型数值模拟2.1 s时刻波场快照(时间步长τ=1.0 ms)

对比两组波场快照,可以看出:采用时间步长τ=1.5 ms时,空间域交错网格有限差分法(M=10)存在严重的数值频散,而时空域交错网格有限差分法(M=10)仅表现轻微的数值频散;减小时间步长至τ=1.0 ms,空间域交错网格有限差分法(M=10)的数值频散明显减弱,但仍存在较明显的数值频散,时空域交错网格有限差分法(M=10)的数值频散进一步减弱,无明显的数值频散。因此,采用相同时间步长时,即计算效率基本相同时,时空域交错网格有限差分法比空间域交错网格有限差分法能更有效地压制数值频散,获得更高的模拟精度。

进一步分析可以看出,时空域交错网格有限差分法(M=10)采用时间步长τ=1.5 ms比空间域交错网格有限差分法(M=10)采用时间步长τ=1.0 ms数值频散更弱,模拟精度更高。因此,相比空间域交错网格有限差分法,时空域交错网格有限差分法能采用更大的时间步长以提高计算效率,同时达到更高的模拟精度。

5.2 塔里木盆地典型复杂构造模型

图7为塔里木盆地典型复杂构造模型,模型网格规模为nz×nx=1 051×1 201,网格尺寸为Δx=Δz=h=15 m。主频25 Hz的Ricker子波作为震源,位于网格点(501,3)。两种差分法均采用时间步长τ=1.0 ms进行数值模拟。图8为两种差分法数值模拟3.0 s时刻压力场P的波场快照,图9为两种差分法数值模拟得到的压力场P的炮集记录。

图7 塔里木盆地典型复杂构造速度模型

图8 塔里木盆地复杂构造模型数值模拟3.0 s时刻压力场P的波场快照(时间步长τ=1.0 ms)

图9 塔里木盆地复杂构造模型数值模拟压力场P的炮集记录(时间步长τ=1.0 ms)

对比图8和图9可以看出,空间域交错网格有限差分法(M=10)数值模拟得到的波场快照和炮集记录均表现出明显的数值频散,模拟精度低;时空域交错网格有限差分法(M=10)数值模拟的波场快照和炮集记录均无明显的数值频散,模拟精度高。塔里木盆地典型复杂构造模型模拟结果进一步证明了时空域交错网格有限差分法能更有效地压制数值频散,从而提高模拟精度。