科学技术与工程

标量波方程时间高阶广义有限差分法及稳定性条件

PDF下载

袁野 ¹ , 黄健良 ¹ , 陶维祥 ¹ , 柳万春 ¹ , 吴国忱 ²

科学技术与工程 | 论文·天文学、地球科学 2025,25(21): 8796-8804

收起

科学技术与工程 | 论文·天文学、地球科学 2025, 25(21): 8796-8804

标量波方程时间高阶广义有限差分法及稳定性条件

全屏

袁野¹, 黄健良¹, 陶维祥¹, 柳万春¹, 吴国忱²

作者信息

¹ 中国海洋石油国际有限公司海外技术支持中心, 北京 100028

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

袁野(1987—),男,汉族,湖北荆州人,硕士,高级工程师。研究方向:地球物理学、综合地质。E-mail:yuanye3@cnooc.com.cn。

Time High-order Generalized Finite Difference Method and Stability Condition for Scalar Wave Equation

Ye YUAN¹, Jian-liang HUANG¹, Wei-xiang TAO¹, Wan-chun LIU¹, Guo-chen WU²

Affiliations

¹ Global Technology Supporting Center, CNOOC International Ltd., Beijing 100028, China

² School of Geosciences, China University of Petroleum, Qingdao 266580, China

出版时间: 2025-07-28 doi: 10.12404/j.issn.1671-1815.2405936

文章导航

摘要

收起

广义有限差分地震波场数值模拟方法能够适应起伏的地层界面,消除起伏界面造成的阶梯状散射现象,提高正演模拟的准确性。而使用时间二阶广义有限差分法求解波动方程时,由于时间差分精度低,时间采样间隔较大时往往会产生时间频散,影响正演模拟的精度。研究了标量波方程时间四阶广义有限差分正演模拟算法及其稳定性条件和频散特性,通过将时间四阶偏导数转嫁到空间偏导数项上实现时间四阶精度差分,时间频散得到有效的压制。此外相对于时间二阶广义有限差分法,时间四阶广义有限差分可以适应较大的时间采样间隔,一定程度上减少计算量。实验结果表明,所提出的算法能有效压制阶梯状散射和时间频散,具有更高的计算精度,将其应用于逆时偏移,可以获得高质量的成像剖面。

关键词

广义有限差分 / 阶梯状散射 / 高阶差分格式 / 时间频散 / 稳定性条件

Abstract

收起

The generalized finite difference method for seismic wavefields numerical simulation is capable of adapting to undulating stratigraphic interfaces, eliminating the staircase scattering effect caused by such interfaces, and enhancing the accuracy of forward modeling. However, when the second-order generalized finite difference method is used to solve the wave equation, low temporal accuracy can lead to temporal dispersion at larger time intervals, affecting the precision of forward simulation. A fourth-order generalized finite-difference forward modeling algorithm for the scalar wave equation was studied, along with its stability conditions and dispersion characteristics. By transferring the fourth-order time derivative to the spatial derivative term, fourth-order accuracy in time was achieved, effectively suppressing temporal dispersion. Compared to the second-order generalized finite-difference method, the fourth-order approach allows for larger time intervals, thereby reducing computational costs to some extent. Experimental results demonstrate that the proposed algorithm effectively mitigates both staircase scattering and temporal dispersion, yielding higher computational accuracy. When applied to reverse time migration, it produces high-quality imaging profiles.

Key words

generalized finite difference method / stepped scattering / high-order difference scheme / time dispersion / stability conditions

引用本文

袁野, 黄健良, 陶维祥, 柳万春, 吴国忱. 标量波方程时间高阶广义有限差分法及稳定性条件. 科学技术与工程, 2025 , 25 (21) : 8796 -8804 . DOI: 10.12404/j.issn.1671-1815.2405936

Ye YUAN, Jian-liang HUANG, Wei-xiang TAO, Wan-chun LIU, Guo-chen WU. Time High-order Generalized Finite Difference Method and Stability Condition for Scalar Wave Equation[J]. Science Technology and Engineering, 2025 , 25 (21) : 8796 -8804 . DOI: 10.12404/j.issn.1671-1815.2405936

正文

收起

地震波动方程正演模拟经常使用的方法包括伪谱法^[1-2]、有限元法^[3-4]以及有限差分法^[5-8]等多种数值模拟方法,其中有限差分法因简便高效而被广泛使用,在实际资料的应用中也取得了较好的效果^[9-10]。但是常规有限差分方法处理复杂地形存在问题,当地层界面崎岖不平时,使用规则网格对速度模型离散时存在阶梯状网格,在正演模拟过程中会产生虚假的绕射波^[11-12],影响数值模拟结果的准确性。为了消除“阶梯状”虚假反射,研究人员在规则网格有限差分法的基础上提出了变网格有限差分算法^[13]、贴体网格有限差分法^[14-15]等多种差分方法,取得了一定的效果。

近些年来,无网格广义有限差分法逐渐受到研究者们的关注,相关的研究工作也很多。广义有限差分法(generalized finite difference method,GFDM)的理论基础是多元函数泰勒函数展开和加权最小二乘拟合,使用相邻节点函数值的线性组合表示待求参数的偏导数,克服了传统方法对网格的依赖性^[16-18],节点剖分灵活度高,适用于复杂模型的正演模拟。Ureña等^[19]使用广义有限差分法求解弹性波动方程,讨论了弹性波方程广义有限差分格式的频散及稳定性条件。Benito等^[20]将其用于研究地震波在各向异性介质中的传播特征。贾宗锋等^[21]将广义有限差分法用于标量波方程的正演模拟,消除了起伏界面阶梯状虚假散射。Li等^[22]进行流-固耦合波动方程广义有限差分法正演模拟,并基于此实现了全波形反演速度建模,取得了良好的效果。但是目前关于广义有限差分法在地震勘探领域的研究多以理论为主,在实际资料处理方面的研究还很少见到。

波动方程数值模拟不仅存在空间频散,也存在时间频散,在地震波速度、子波主频不变的情况下,可以通过减小空间网格压制空间频散,减小时间步长压制时间频散。但是小的空间网格间距和小的时间步长就需要更长的模拟时间,意味着计算效率的降低,因此可以采用时间高阶差分格式求解波动方程来压制时间频散,且高阶时间差分格式可以允许更大的时间步长,一定程度上提高计算效率^[23-25]。

现将广义有限差分算法推导至时间四阶,通过将时间四阶项转嫁到空间项,建立标量波方程时间四阶广义有限差分格式,并分析差分格式的稳定性条件和频散特性。同时,结合贴体节点生成方法,消除常规有限差分处理起伏界面存在的虚假反射问题,提高正演模拟精度。通过不同的模型进行试算,验证时间四阶广义有限差分正演模拟方法的稳定性和准确性。同时,将所提出的时间四阶广义有限差分正演模拟算法应用于逆时偏移成像,获得精度较高的成像剖面。以期达到压制时间频散、提高正演模拟稳定性的目的。

1 广义有限差分法原理

收起

广义有限差分法是一种无网格方法,和常规有限差分法不同,其节点离散方式是任意的,不依赖于网格(图1)。广义有限差分法求解偏微分方程时,先将求解区域内的节点进行多元函数泰勒展开,然后构建残差函数,残差函数对偏导数向量求导并令偏导数值为零,可得差分系数矩阵方程,求解矩阵方程可以得到广义有限差分系数及各项偏导数的表达式。

根据多元函数泰勒展开,求解域内的任意一点u_s在点u₀附近的任意阶泰勒展开式和残差函数为

(1)$ u_{s}=u_{0}+\sum_{r=1}^{R} \frac{1}{r!} \sum_{a=0}^{r} C_{r}^{a} \frac{\partial^{r} u_{0}}{\partial x^{r-a} \partial z^{a}} \Delta x_{s}^{r-a} \Delta z_{s}^{a} $

(2)$ \begin{aligned} B(u)= & \sum_{s=1}^{N_{s}}\left(u_{0}-u_{s}+\sum_{r=1}^{R} \frac{1}{r!} \sum_{a=0}^{r} C_{r}^{a} \times\right. \\ & \left.\frac{\partial^{r} u_{0}}{\partial x^{r-a} \partial z^{a}} \Delta x_{s}^{r-a} \Delta z_{s}^{a}\right)^{2} w_{s}^{2}\left(\Delta x_{s}, \Delta z_{s}\right) \end{aligned} $

式中:r为泰勒展开式的阶数;Δx_s=x_s-x₀,Δz_s=z_s-z₀,分别表示差分周围节点u_s到中心节点u₀的水平方向和垂直方向上的距离;w_s(Δx_s,Δz_s)为权重函数,表示周围节点u_s对中心节点u₀的影响程度。权重函数值越大,表明该点距离中心点越近,对中心点的影响越大。

记中心点u₀处的偏导数向量D为

(3)$ \boldsymbol{D}=\left(\frac{\partial u_{0}}{\partial x} \frac{\partial u_{0}}{\partial z} \frac{\partial^{2} u_{0}}{\partial x^{2}} \frac{\partial^{2} u_{0}}{\partial x \partial z} \cdots \frac{\partial^{r} u_{0}}{\partial x^{r-a} \partial z^{a}}\right)^{\mathrm{T}} $

利用残差函数对偏导数向量D求导,并令偏导数值为零,可得

(4)

∂ B (u) ∂ D

差分系数矩阵方程为

(5)AD=b

通过解系数矩阵方程可求得各阶偏导数的值(具体求解过程可见文献[16 -22])。广义有限差分的空间项偏导数可以表示为

(6)$ \boldsymbol{D}=\boldsymbol{A}^{-1} \boldsymbol{b}=-\boldsymbol{M}_{0}^{q} u_{0}+\sum_{s=1}^{N_{s}} \boldsymbol{M}_{s}^{q} u_{s} $

式(6)中:

M 0 q

∑ s = 1 N s M s q

,表示差分模板内各个节点的差分系数。

2 标量波方程时间四阶广义有限差分格式

收起

波动方程数值模拟在空间和时间上均会产生数值频散,因此要在空间和时间两方面同时压制频散才会取得理想的效果^[26],为了满足高频粗网格大时间步长的正演模拟精度需求,除了提高空间差分阶数外,还要采用时间高阶差分。因此建立标量波方程时间四阶广义有限差分格式,提高正演模拟的精度。

针对二维标量波方程有

(7)

1 v 2 ∂ 2 u ∂ t 2

∂ 2 u ∂ x 2

∂ 2 u ∂ z 2

+f(x,z,t)

式(7)中:v为介质速度;f(x,z,t)为震源项。

将泰勒展开式保留到时间四阶项,并将时间四阶项转嫁到空间项上进行求取,可得时间项二阶偏导数的四阶表达形式(Dablain^[23])为

(8) $\begin{aligned} \frac{\partial^{2} u}{\partial t^{2}}= & \frac{1}{(\Delta t)^{2}}[u(t+\Delta t)+u(t-\Delta t)- \\ & \left.2 u(t)-\frac{(v \Delta t)^{4}}{12}\left(\frac{\partial^{4} u}{\partial x^{4}}+\frac{\partial^{4} u}{\partial z^{4}}+2 \frac{\partial^{4} u}{\partial x^{2} \partial z^{2}}\right)\right] \end{aligned}$

空间项偏导数的广义有限差分近似表示如下。

(9)

∂ 2 u 0 ∂ x 2 = - ξ 0 u 0 + ∑ s = 1 N s ξ s u s ∂ 2 u 0 ∂ z 2 = - η 0 u 0 + ∑ s = 1 N s η s u s ∂ 4 u 0 ∂ x 4 = - ξ ξ 0 u 0 + ∑ s = 1 N s ξ ξ s u s ∂ 4 u 0 ∂ z 4 = - η η 0 u 0 + ∑ s = 1 N s η η s u s ∂ 4 u 0 ∂ x 2 ∂ z 2 = - ξ η 0 u 0 + ∑ s = 1 N s ξ η s u s

将式(8)、式(9)代入式(7)中,不考虑震源项,可得标量波方程时间四阶广义有限差分格式为

(10)$ \begin{aligned} u_{0}^{n+1}= & 2 u_{0}^{n}-u_{0}^{n-1}+(v \Delta t)^{2}\left(-\xi_{0} u_{0}^{n}+\sum_{s=1}^{N_{s}} \xi_{s} u_{s}^{n}-\right. \\ & \left.\eta_{0} u_{0}^{n}+\sum_{s=1}^{N_{s}} \eta_{s} u_{s}^{n}\right)+\frac{(v \Delta t)^{4}}{12}\left[-\xi \xi_{0} u_{0}^{n}+\right. \\ & \sum_{s=1}^{N_{s}} \xi \xi_{s} u_{s}^{n}-\eta \eta_{0} u_{0}^{n}+\sum_{s=1}^{N_{s}} \eta \eta_{s} u_{s}^{n}+ \\ & \left.2\left(-\xi \eta_{0} u_{0}^{n}+\sum_{s=1}^{N_{s}} \xi \eta_{s} u_{s}^{n}\right)\right] \end{aligned} $

正演模拟过程中通常会涉及边界反射的问题,此处采用完全匹配层(perfectly matched layer,PML)边界条件对边界反射进行吸收衰减^[27]。时间域PML标量波方程为

(11)

∂ 2 u x ∂ t 2 + 2 d x (x) ∂ u x ∂ t + d x 2 (x) u x = v 2 ∂ 2 u ∂ x 2 ∂ 2 u z ∂ t 2 + 2 d z (z) ∂ u z ∂ t + d z 2 (z) u z = v 2 ∂ 2 u ∂ z 2

式(11)中:d_x(x)、d_z(z)分别为x方向和z方向的衰减系数。

将时间偏导数和空间偏导数代入时间域PML标量波方程中,最终可以得到PML标量波方程时间四阶广义有限差分格式如式(12)、式(13)所示。

(12)$ \begin{aligned} u_{x}^{n+1}= & \frac{1}{1+\frac{1}{6} d_{x}^{3}(x) \Delta t^{3}+d_{x}(x) \Delta t}\left\{\left[2-\frac{1}{4} d_{x}^{4}(x) \Delta t^{4}-\right.\right. \\ & \left.d_{x}^{2}(x) \Delta t^{2}\right] u_{x}^{n}+\left[\frac{1}{6} d_{x}^{3}(x) \Delta t^{3}+d_{x}(x) \Delta t-\right. \\ & 1] u_{x}^{n-1}+\left[\frac{1}{4} v^{2} \Delta t^{4} d_{x}^{2}(x)+v^{2} \Delta t^{2}\right] \times \\ & \left(-\xi_{0} u_{0}+\sum_{s=1}^{N_{s}} \xi_{s} u_{s}\right)+\frac{1}{12} v^{4} \Delta t^{4}\left(-\xi \xi_{0} u_{0}+\right. \\ & \left.\sum_{s=1}^{N_{s}} \xi \xi_{s} u_{s}-\xi \eta_{0} u_{0}+\sum_{s=1}^{N_{s}} \xi \eta_{s} u_{s}\right)-\frac{1}{12} v^{4} \Delta t^{5} d_{x}(x) \times \\ & {\left[2 \sum_{t=1}^{n-1}\left(-\xi \xi_{0} u_{0}^{t}+\sum_{s=1}^{N_{s}} \xi \xi_{s} u_{s}^{t}-\xi \eta_{0} u_{0}^{t}+\sum_{s=1}^{N_{s}} \xi \eta_{s} u_{s}^{t}\right)+\right.} \\ & \left.\left.\left(-\xi \xi_{0} u_{0}+\sum_{s=1}^{N_{s}} \xi \xi_{s} u_{s}-\xi \eta_{0} u_{0}+\sum_{s=1}^{N_{s}} \xi \eta_{s} u_{s}\right)\right]\right\} \end{aligned} $

(13)$ \begin{aligned} u_{z}^{n+1}= & \frac{1}{1+\frac{1}{6} d_{z}^{3}(z) \Delta t^{3}+d_{z}(z) \Delta t}\left\{\left[2-\frac{1}{4} d_{z}^{4}(z) \Delta t^{4}-\right.\right. \\ & \left.d_{z}^{2}(z) \Delta t^{2}\right] u_{z}^{n}+\left[\frac{1}{6} d_{z}^{3}(z) \Delta t^{3}+d_{z}(z) \Delta t-\right. \\ & 1] u_{z}^{n-1}+\left(\frac{1}{4} v^{2} \Delta t^{4} d_{z}^{2}(z)+v^{2} \Delta t^{2}\right) \times \\ & \left(-\eta_{0} u_{0}+\sum_{s=1}^{N_{s}} \eta_{s} u_{s}\right)+\frac{1}{12} v^{4} \Delta t^{4}\left(-\xi \eta_{0} u_{0}+\right. \\ & \left.\sum_{s=1}^{N_{s}} \xi \eta_{s} u_{s}-\eta \eta_{0} u_{0}+\sum_{s=1}^{N_{s}} \eta \eta_{s} u_{s}\right)-\frac{1}{12} v^{4} \Delta t^{5} d_{z}(z) \times \\ & {\left[2 \sum_{t=1}^{n-1}\left(-\xi \eta_{0} u_{0}^{t}+\sum_{s=1}^{N_{s}} \xi \eta_{s} u_{s}^{t}-\eta \eta_{0} u_{0}^{t}+\sum_{s=1}^{N_{s}} \eta \eta_{s} u_{s}^{t}\right)+\right.} \\ & \left.\left.\left(-\xi \eta_{0} u_{0}+\sum_{s=1}^{N_{s}} \xi \eta_{s} u_{s}-\eta \eta_{0} u_{0}+\sum_{s=1}^{N_{s}} \eta \eta_{s} u_{s}\right)\right]\right\} \end{aligned} $

3 稳定性条件及频散特性

收起

3.1 稳定性条件

差分格式的稳定性是数值求解波动方程至关重要的一方面,在稳定性条件容许范围选择合适的时间步长既保证了地震波场长时间传播的稳定性,又能减少迭代的次数,故基于平面波解求取标量波方程时间四阶广义有限差分格式的稳定性条件。

二维标量波方程的平面波解为

(14)

u 0 n = A ζ n e - j (k x x 0 + k z z 0) u s n = A ζ n e - j (k x x s + k z z s) ζ = e j w Δ t

式(14)中:k_x=ksinθ,k_z=kcosθ,k=ω/v_ph为波数,v_ph为相速度,θ为波的传播方向与z轴的夹角。

将平面波解代入差分格式[式(10)]中并进行整理可得

(15)$ \begin{aligned} 2 \cos (\omega \Delta t)-2= & (v \Delta t)^{2}\left\{\sum_{s=1}^{N_{s}}\left(\xi_{s}+\eta_{s}\right) \times\right. \\ & {\left[\cos \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)-1\right]+} \\ & \left.\mathrm{j} \sum_{s=1}^{N_{s}}\left(\xi_{s}+\eta_{s}\right) \sin \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)\right\}+ \\ & \frac{(v \Delta t)^{4}}{12}\left\{\sum_{s=1}^{N_{s}}\left(\xi \xi_{s}+\eta \eta_{s}+2 \xi \eta_{s}\right) \times\right. \\ & {\left[\cos \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)-1\right]+} \\ & \mathrm{j} \sum_{s=1}^{N_{s}}\left(\xi \xi_{s}+\eta \eta_{s}+2 \xi \eta_{s}\right) \\ & \left.\sin \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)\right\} \end{aligned} $

由于差分模板中的节点不一定是均匀分布的,虚数部分不能忽略,因此将式(15)拆分成实数和虚数两部分。

实数部分为

(16)$ \begin{aligned} 2 \cos (\omega \Delta t)-2= & (v \Delta t)^{2} \sum_{s=1}\left(\xi_{s}+\eta_{s}\right) \times \\ & {\left[\cos \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)-1\right]+} \\ & \frac{(v \Delta t)^{4}}{12} \sum_{s=1}^{N_{s}}\left(\xi \xi_{s}+\eta \eta_{s}+2 \xi \eta_{s}\right) \times \\ & {\left[\cos \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)-1\right] } \end{aligned} $

虚数部分为

(17)$ \begin{aligned} 0= & (v \Delta t)^{2} \sum_{s=1}^{N_{s}}\left(\xi_{s}+\eta_{s}\right) \sin \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)+ \\ & \frac{(v \Delta t)^{4}}{12} \sum_{s=1}^{N_{s}}\left(\xi \xi_{s}+\eta \eta_{s}+2 \xi \eta_{s}\right) \times \\ & \sin \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right) \end{aligned} $

根据余弦函数的性质,式(16)左边的取值范围为[-4,0],且式(16)右端项小于零,因此对式(16)整理可得

(18)$ \begin{array}{l} \frac{(v \Delta t)^{4}}{12} \sum_{s=1}^{N_{s}}\left(\xi \xi_{s}+\eta \eta_{s}+2 \xi \eta_{s}\right)\left[1-\cos \left(k_{x} \Delta x_{s}+\right.\right. \\ \left.\left.k_{z} \Delta z_{s}\right)\right]+(v \Delta t)^{2} \sum_{s=1}^{N_{s}}\left(\xi_{s}+\eta_{s}\right)[1- \\ \left.\cos \left(k_{x} \Delta x_{s}+k_{z} \Delta z_{s}\right)\right] \leqslant 4 \end{array} $

当cos(k_xΔx_s+k_zΔz_s)=-1时,式(17)右端项sin(k_xΔx_s+k_zΔz_s)=0,不等式(18)左端取最大值,可得标量波方程时间四阶广义有限差分格式的稳定性条件为

(19)$ \Delta t \leqslant \frac{1}{v} \sqrt{\frac{6\left(\sqrt{b^{2}+\frac{2}{3} a}-b\right)}{a}} $

式(19)中:a=ξξ₀+ηη₀+2ξη₀,b=ξ₀+η₀。

3.2 频散特性

采用相速度的相对误差对时间四阶广义有限差分格式的时间频散进行分析,讨论时间网格频散与时间步长及频率成分之间的关系。

标量波方程时间四阶精度的差分格式为

(20)$ \begin{aligned} u_{0}^{n+1}= & 2 u_{0}^{n}-u_{0}^{n-1}+(v \Delta t)^{2}\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial z^{2}}\right)+ \\ & \frac{(v \Delta t)^{4}}{12}\left(\frac{\partial^{4} u}{\partial x^{4}}+2 \frac{\partial^{4} u}{\partial x^{2} \partial z^{2}}+\frac{\partial^{4} u}{\partial z^{4}}\right) \end{aligned} $

将二维标量波方程的平面波解代入式(20)中,并对方程进行整理可得

(21)$ \begin{array}{c} \frac{(\omega \Delta t)^{4}}{12}\left(\frac{v_{\mathrm{ph}}}{v}\right)^{4}-(\omega \Delta t)^{2}\left(\frac{v_{\mathrm{ph}}}{v}\right)^{2}- \\ {[2 \cos (\omega \Delta t)-2]=0} \end{array} $

求解方程(21)可得时间四阶精度差分格式的频散关系为

(22)$ \left(\frac{v_{\mathrm{ph}}}{v}\right)^{2}=\frac{6\{1-\sqrt{1+2 / 3[\cos (\omega \Delta t)-1]}\}}{(\omega \Delta t)^{2}} $

设G=2π/ωΔt,表示一次震动包含的时间网格点数,则式(22)化为

(23)$ \left(\frac{v_{\mathrm{ph}}}{v}\right)^{2}=\frac{6\left\{1-\sqrt{1+2 / 3\left[\cos \left(\frac{2 \pi}{G}\right)-1\right]}\right\}}{\left(\frac{2 \pi}{G}\right)^{2}} $

对比二阶精度差分格式的频散关系和四阶精度差分格式的频散关系可以看出(图2),时间频散的相速度超前于实际地震波速度,在频率成分较低,时间采样间隔较小时,1/G较小,相速度与实际地震波速度的差异较小,随着频率增加,时间采样间隔增大,时间二阶精度的相速度与实际地震波速度的差异增大的较快,远大于时间四阶精度。

4 模型试算及对比分析

收起

4.1 层状起伏模型

首先采用层状起伏模型进行正演模拟来验证本文算法的有效性,速度模型如图3(a)所示,模型采样点数为401×201,纵、横向采样间隔均为15 m,时间采样间隔0.002 s,采样时长2.0 s,震源采用雷克子波,主频22 Hz,震源位于(201,1)处。分别采用时间二阶精度规则网格有限差分法、时间二阶精度广义有限差分法、时间四阶精度广义有限差分法求解标量波方程,得到它们的波场快照和正演记录。

图3(b)~图3(d)、图4(a)~图4(c)分别展示了使用时间二阶规则网格有限差分法、时间二阶广义有限差分法、时间四阶广义有限差分法计算得到的波场快照和正演记录,从中可以看出时间二阶规则网格有限差分法正演模拟结果存在“阶梯状”散射,并且在波形前有明显的时间频散,而无论是采用时间二阶精度还是时间四阶精度,广义有限差分法得到的正演模拟结果均不存在“阶梯状”散射,但是时间二阶精度广义有限差分法的正演模拟结果存在时间频散,而时间四阶精度广义有限差分法则不存在时间频散现象,明显提高了数值模拟的精度。

图5(a)~图5(c)分别展示了使用时间二阶规则网格有限差分法、时间二阶广义有限差分法、时间四阶广义有限差分法得到的逆时偏移成像剖面,可以看出,由于“阶梯状”散射的干扰,时间二阶规则网格有限差分法的成像剖面存在不连续现象,成像效果较差,而时间二阶广义有限差分法的成像剖面则由于时间频散的缘故,分辨率降低,时间四阶广义有限差分法得到的成像效果最好,验证了本文算法的正确性和有效性。

4.2 复杂模型

为进一步验证时间四阶精度广义有限差分算法的正确性和稳定性,采用更为复杂的多层起伏模型进行正演模拟,速度模型如图6所示,模型采样点数为601×301,纵、横向采样间隔均为8 m,时间采样间隔0.001 s,采样时长2.4 s,震源采用雷克子波,主频27 Hz,震源位于(301,1)处。分别采用时间二阶广义有限差分法、时间四阶广义有限差分法求解标量波方程,得到它们的正演模拟记录。

图7(a)、图7(b)分别展示了使用时间二阶精度广义有限差分法、时间四阶精度广义有限差分法得到的正演模拟记录。相较于时间二阶精度广义有限差分法,采用时间四阶精度广义有限差分法进行正演模拟能够有效地压制时间频散现象,得到高精度的正演模拟结果,验证了所提出的时间四阶精度广义有限差分算法对复杂模型的适用性。

图8(a)、图8(b)分别展示了使用时间二阶广义有限差分法、时间四阶广义有限差分法得到的逆时偏移成像剖面,可以看出,时间四阶广义有限差分法的成像剖面好于时间二阶精度的成像剖面,进一步验证了本文算法的正确性及其对复杂模型的适用性。

4.3 实际资料

为了说明方法的准确性与实用性,使用某工区实际速度模型进行正演模拟,图9展示了该工区的某二维测线速度模型,图10是实际资料炮集,实际资料炮记录采用拖缆采集方式,拖缆长度8 km,因此实际单炮使用的正演模型大小为8 km × 7.5 km。

分别采用时间二阶广义有限差分法、时间四阶广义有限差分法进行正演模拟,时间采样间隔0.003 s,采样时长7.5 s,震源采用雷克子波,主频9 Hz,右侧放炮,左侧接收。图11(a)和图11(b)分别是使用时间二阶广义有限差分法、时间四阶广义有限差分法得到的正演模拟记录。从正演模拟记录中可以出,时间四阶广义有限差分法得到的正演模拟记录的时间频散得到了很好的压制,与实际资料炮集更加一致,验证了本文算法的正确性与实用性。

5 结论

收起

推导了标量波方程的时间四阶广义有限差分格式,将时间高阶偏导数项转嫁到空间项上求取,避免了多时间层的运算。并对时间四阶精度广义有限差分格式的稳定性和时间频散特性进行了分析,通过数值模拟验证了差分格式的精度与稳定性,取得了以下结论。

(1) 广义有限差分法节点剖分灵活,不受网格的束缚,可以根据速度模型建立合适的场节点分布形式,避免了常规有限差分数值模拟中网格剖分与速度界面不一致的情况,有效地消除了模拟过程中由于界面起伏产生的“阶梯状”绕射。

(2) 在参数相同的条件下,时间四阶广义有限差分格式的正演模拟精度明显高于时间二阶,有效地压制了由时间方向上差分引起的数值频散。

(3) 时间四阶精度广义有限差分算法能够有效压制“阶梯状”散射和时间频散,基于本文算法的逆时偏移成像可以获得高质量的成像剖面。

基金

收起

中国海洋石油集团有限公司“十四五”重大科技项目(KJGG2022-0903)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Kosloff

D D

, Baysal

. Forward modeling by a Fourier method[J]. Geophysics, 1982, 47(10): 1402-1412.

[2]

唐杰, 刘英昌, 韩盛元, 等. 黏弹性VTI介质微地震正演模拟与震源机制全波形反演[J]. 地球物理学报, 2022, 65(1): 301-319.

Tang

Jie

, Liu

Yingchang

, Han

Shengyuan

, et al. Microseimic forward modeling and full-waveform inversion of focal mechanisms for viscoelastic vertically transverse isotropic media[J]. Chinese Journal of Geophysics, 2022, 65(1): 301-319.

[3]

张剑锋. 一种高精度有限元弹性波场正演方法[J]. 石油地球物理勘探, 1994, 29(1): 29-36.

Zhang

Jianfeng

. An accurate method for finite-element elastic wave field modeling[J]. Oil Geophysical Prospecting, 1994, 29(1): 29-36.

[4]

Santos

J E

, Carcione

J M

, Savioli

G B

, et al. Wave propagation in thermo-poroelasticity: a finite-element approach[J]. Geophysics, 2023, 88(1): 161-175.

[5]

Alterman

Z S

, Karal

F C J R

. Propagation of elastic waves in layered media by finite difference methods[J]. Bulletin of the Seismological Society of America, 1968, 58(1): 367-398.

[6]

董良国, 马在田, 曹景忠, 等. 一阶弹性波方程交错网格高阶差分解法[J]. 地球物理学报, 2000, 43(3): 411-419.

Dong

Liangguo

, Ma

Zaitian

, Cao

Jingzhong

, et al. A staggered-grid high-order difference method of one-order elastic wave equation[J]. Chinese Journal of Geophysics, 2000, 43(3): 411-419.

[7]

周诚尧, 汪勇, 桂志先, 等. 二维黏弹介质五点八阶超紧致有限差分声波方程数值模拟[J]. 科学技术与工程, 2020, 20(1): 54-63.

Zhou

Chengyao

, Wang

Yong

, Gui

Zhixian

, et al. Numerical simulation of 2D viscoelastic acoustic wave equation used five point eight order combined compact difference scheme[J]. Science Technology and Engineering, 2020, 20(1): 54-63.

[8]

邓光校, 鲍羽, 汪勇. 基于紧致交错差分格式的二维声波及黏滞声波方程数值模拟[J]. 科学技术与工程, 2023, 23(21): 8966-8979.

Deng

Guangxiao

, Bao

, Wang

Yong

. Numerical simulation of two dimensional acoustic and viscous acoustic wave equations based on compact staggered difference scheme[J]. Science Technology and Engineering, 2023, 23(21): 8966-8979.

[9]

郑雪静, 李俊伦, 熊强青, 等. 主动源全波形反演在内蒙古居延海坳陷多资源勘查中的应用[J]. 地球物理学报, 2024, 67(8): 3120-3135.

Zheng

Xuejing

, Li

Junlun

, Xiong

Qiangqing

, et al. Application of active-source full waveform inversion in multi-resource exploration in Juyanhai depression of Inner Mongolia[J]. Chinese Journal of Geophysics, 2024, 67(8): 3120-3135.

[10]

王璇, 王腾飞, 程玖兵, 等. 超深水区基准面延拓驱动的反射波形反演成像[J]. 石油地球物理勘探, 2024, 59(4): 790-799.

Wang

Xuan

, Wang

Tengfei

, Cheng

Jiubing

, et al. Reflection waveform inversion imaging driven by redatuming in ultra-deepwater area[J]. Oil Geophysical Prospecting, 2024, 59(4): 790-799.

[11]

李振春, 肖建恩, 曲英铭, 等. 时间域起伏自由地表正演模拟综述[J]. 地球物理学进展, 2016, 31(1): 300-309.

Zhenchun

, Xiao

Jianen

, Qu

Yingming

, et al. The review of irregular free-surface forward modeling in time domain[J]. Progress in Geophysics, 2016, 31(1): 300-309.

[12]

Yao

, da

Silva N V

, Debens

H A

, et al. Accurate seabed modeling using finite difference methods[J]. Computational Geosciences, 2018, 22: 469-484.

[13]

曲英铭, 李振春, 黄建平, 等. 基于多尺度双变网格的时间域全波形反演[J]. 石油物探, 2016, 55(2): 241-250.

Yingming

, Li

Zhenchun

, Huang

Jianping

, et al. Full waveform inversion based on multi-scale dual-variable grid in time domain[J]. Geophysical Prospecting for Petroleum, 2016, 55(2): 241-250.

[14]

李庆洋, 黄建平, 李振春, 等. 起伏地表贴体全交错网格仿真型有限差分正演模拟[J]. 石油地球物理勘探, 2015, 50(4): 633-642.

Qingyang

, Huang

Jianping

, Li

Zhenchun

, et al. Undulating surface body-fitted grid seismic modeling based on fully staggered-grid mimetic finite difference[J]. Oil Geophysical Prospecting, 2015, 50(4): 633-642.

[15]

李庆洋, 李振春, 黄建平, 等. 基于贴体全交错网格的起伏地表正演模拟影响因素[J]. 吉林大学学报(地球科学版), 2016, 46(3): 920-929.

Qingyang

, Li

Zhenchun

, Huang

Jianping

, et al. Factor analysis of seismic modeling with topography based on a fully staggered body-fitted grids[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(3): 920-929.

[16]

Benito

J J

, Urena

, Gavete

. Influence of several factors in the generalized finite difference method[J]. Applied Mathematical Modelling, 2001, 25(12): 1039-1053.

[17]

张挺, 任聿飞, 杨志强, 等. 基于广义有限差分法求解非线性自由液面的液体晃动问题[J]. 四川大学学报: 工程科学版, 2016, 48(1): 8-14.

Zhang

Ting

, Ren

Yufei

, Yang

Zhiqiang

, et al. Generalized finite difference method for non-linear free-surface problems of liquid sloshing[J]. Journal of Sichuan University (Engineering Science Edition), 2016, 48(1): 8-14.

[18]

陈剑, 刘春明, 王茂海, 等. 广义有限差分法在静态电磁场计算中的应用[J]. 电工技术学报, 2018, 33(7): 1579-1587.

Chen

Jian

, Liu

Chunming

, Wang

Maohai

, et al. Application of the generalized finite difference method to static electromagnetic problems[J]. Transactions of China Electrotechnical Society, 2018, 33(7): 1579-1587.

[19]

Ureña

, Benito

J J

, Salete

, et al. A note on the application of the generalized finite difference method to seismic wave propagation in 2D[J]. Journal of Computational and Applied Mathematics, 2012, 236(12): 3016-3025.

[20]

Benito

J J

, Ureña

, Gavete

, et al. A GFDM with PML for seismic wave equations in heterogeneous media[J]. Journal of Computational and Applied Mathematics, 2013, 252: 40-51.

[21]

贾宗锋, 吴国忱, 李青阳, 等. 标量波方程广义有限差分正演模拟[J]. 石油地球物理勘探, 2022, 57(1): 101-110.

Jia

Zongfeng

, Wu

Guochen

, Li

Qingyang

, et al. Forward modeling of scalar wave equation with generalized finite difference method[J]. Oil Geophysical Prospecting, 2022, 57(1): 101-110.

[22]

, Wu

, Jia

, et al. Full-waveform inversion in acoustic-elastic coupled media with irregular seafloor based on the generalized finite-difference method[J]. Geophysics, 2023, 88(2): 83-100.

[23]

Dablain

M A

. The application of high-order differencing to the scalar wave equation[J]. Geophysics, 1986, 51(1): 54-66.

[24]

唐怀谷, 何兵寿. 一阶声波方程时间四阶精度差分格式的伪谱法求解[J]. 石油地球物理勘探, 2017, 52(1): 71-80.

Tang

Huaigu

, He

Bingshou

. Pseudo spectrum method of first-order acoustic wave equation finite difference schemes with fourth-order time difference accuracy[J]. Oil Geophysical Prospecting, 2017, 52(1): 71-80.

[25]

李青阳, 吴国忱, 梁展源. 基于PML边界的时间高阶伪谱法弹性波场模拟[J]. 地球物理学进展, 2018, 33(1): 228-235.

Qingyang

, Wu

Guochen

, Liang

Zhanyuan

. Time domain high-order pseudo spectral method based on PML boundary for elastic wavefield simulation[J]. Progress in Geophysics, 2018, 33(1): 228-235.

[26]

吴国忱, 王华忠. 波场模拟中的数值频散分析与校正策略[J]. 地球物理学进展, 2005(1): 58-65.

Guochen

, Wang

Huazhong

. Analysis of numerical dispersion in wave-field simulation[J]. Progress in Geophysics, 2005(1): 58-65.

[27]

姚铭, 汪勇, 杨晓柳, 等. 不同边界条件下的二维声波方程数值模拟[J]. 科学技术与工程, 2017, 17(32): 11-16.

Yao

Ming

, Wang

Yong

, Yang

Xiaoliu

, et al. Numerical simulation of two-dimensional acoustic equation under different boundary conditions[J]. Science Technology and Engineering, 2017, 17(32): 11-16.

2025年第25卷第21期

PDF下载

152

引用本文

BibTeX

文章信息

doi: 10.12404/j.issn.1671-1815.2405936

接收时间：2024-08-07
首发时间：2026-01-13
出版时间：2025-07-28

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-08-07
修回日期：2025-04-09

基金

中国海洋石油集团有限公司“十四五”重大科技项目(KJGG2022-0903)

作者信息

¹ 中国海洋石油国际有限公司海外技术支持中心, 北京 100028

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

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/kxjsygc/CN/10.12404/j.issn.1671-1815.2405936

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT