石油地球物理勘探  2020, Vol. 55 Issue (5): 1047-1059  DOI: 10.13810/j.cnki.issn.1000-7210.2020.05.013
0
文章快速检索     高级检索

引用本文 

邹强, 黄建平, 雍鹏, 李振春. 平面波优化差分算子弹性波逆时偏移. 石油地球物理勘探, 2020, 55(5): 1047-1059. DOI: 10.13810/j.cnki.issn.1000-7210.2020.05.013.
ZOU Qiang, HUANG Jianping, YONG Peng, LI Zhenchun. An elastic-wave reverse-time migration method based on optimal finite-difference operators using a new plane wave solution. Oil Geophysical Prospecting, 2020, 55(5): 1047-1059. DOI: 10.13810/j.cnki.issn.1000-7210.2020.05.013.

本项研究受国家重点研发计划项目"塔里木盆地盐下超深层大型油气田形成机理与分布预测"(2019YFC0605503)、国家自然科学基金项目"面向深部储层的时空域自适应高斯束成像理论方法及优化"(41874149)和中石油重大科技合作项目"深层复杂高陡构造成像关键技术"(ZD2019-183-003)联合资助

作者简介

邹强  硕士研究生, 1994年生; 2017年获长江大学勘查技术与工程专业学士学位; 现在中国石油大学(华东)攻读地质资源与地质工程专业硕士学位, 主要从事有限差分正演及偏移成像方法方面的学习和研究

黄建平, 山东省青岛市黄岛区长江西路66号中国石油大学(华东)地球科学与技术学院, 266580。Email:jphuang@upc.edu.cn

文章历史

本文于2019年9月12日收到,最终修改稿于2020年6月16日收到
平面波优化差分算子弹性波逆时偏移
邹强①② , 黄建平 , 雍鹏 , 李振春     
① 中国石油大学(华东)地球科学与技术学院, 山东青岛 266580;
② 中国石油塔里木油田公司, 新疆库尔勒 841000
摘要:弹性波逆时偏移基于双程波动方程理论,不受倾角限制,对复杂模型具有较强的适应性,但高额的计算量和内存占用是其用于巨量地震资料处理的主要瓶颈。为此,推导了一种基于平面波优化的等效交错网格弹性波逆时偏移成像方法,通过优化差分系数降低正演波场的频散误差,达到使用低阶优化差分算子替代高阶传统差分算子、降低计算量的目的。在内存需求方面,采用二阶位移弹性波方程,避免了中间变量的出现,且引入了有效边界存储的策略,只利用边界附近区域的部分波场值控制逆时偏移的波场存储方式。模型测试表明:该方法6、10阶空间差分算子逆时偏移结果成像精度分别相当于10、18阶传统方法,有效节约了计算成本;相比于常规弹性波逆时偏移的波场存储方式,有效边界存储在增加部分计算量的前提下,能节约大量内存占用,且使用的差分算子长度越短,内存占用也越少。
关键词弹性波逆时偏移    平面波优化    计算量    内存占用    有效边界存储    
An elastic-wave reverse-time migration method based on optimal finite-difference operators using a new plane wave solution
ZOU Qiang①② , HUANG Jianping , YONG Peng , LI Zhenchun     
① School of Geosciences, China University of Petroleum(East China), Qingdao, Shandong 266580, China;
② PetroChina Tarim Oilfield Company, Korla, Xinjiang 841000, China
Abstract: Elastic-wave revere time migration based on the theory of two-way wave equation is not limited by the dip angle of formation, so it has a strong adaptability to complex models.However, high computing cost and a large memory limit its application for processing big seismic data.This paper deduces an elastic-wave reverse-time migration imaging method using equivalent staggered-grid finite-difference based on time-space domain optimization.To improve imaging accuracy, we reduce the dispersion error by optimizing difference operators, which means that a low-order difference opera-tor can achieve the imaging accuracy of a high-order difference operator, so as to reduce the computing cost.To reduce memory usage, the second-order displacement equation is adopted to avoid intermediate variables.We also adopt an efficient boundary storage strategy, which uses the wavefield value near the boundary to control the reverse-time migration wavefield storage, which effectively reduces the memory use.Model tests show that the new method with a 6th-order spatial difference has a similar imaging accuracy to the conventional staggered-grid scheme with the 10th-order spatial difference, and the new method with a 10th-order corresponds to the conventional method with a 18th-order difference.Compared with conventional elastic wavefield storage strategy, effective boundary storage strategy can save a lot of memory usage on the premise of increasing less computing cost, and this strategy is more beneficial to shorter difference operators.
Keywords: elastic-wave reverse-time migration    plane wave optimization    computational cost    sto-rage usage    efficient boundary storage    
0 引言

逆时偏移的思想起源于数学上的利用有限差分求解波动方程[1],后被广泛应用到地震勘探领域[2-4]。逆时偏移不受倾角限制,能对强横向变速地下结构进行成像。传统地震资料只记录弹性波场的垂直分量,忽略了大部分甚至全部的S波信息[5],相比之下,多分量地震资料的弹性波逆时偏移可以提供更多的地下构造信息。如PS波成像能够补偿P波能量弱的区域的信息,且通过比较PP与PS成像剖面可以验证亮点反射[6]。目前海上多分量地震资料越来越多[7],随着地震采集技术的发展和计算机性能的提高,弹性波逆时偏移的应用也会越来越广泛。

实现逆时偏移成像一般需要三个步骤:震源波场的获取;接收波场的重构;利用成像条件提取成像值[8-9]。因此,地震波场的精确求解是影响成像剖面质量的直接因素[10-11]。由于通常采用有限差分求解波动方程,正演波场会存在数值频散[12-14],对后期成像和反演皆不利。压制数值频散的途径主要有三类:增加差分算子长度[15-17]、通量校正传输[18]、优化差分系数[19-23]。前两类方法都能在一定程度上压制数值频散,但会引入额外的计算量,且不能保证算法的稳定性[24]。优化差分系数的方法在不增加计算量、不破坏算法的稳定性和收敛性的前提下,从数值频散的本质出发,通过合理布设各个网格的权值分配,使数值解更接近于真实解,是压制数值频散的有效方法。传统上,常采用泰勒近似求取差分系数[25],由于直接忽略高阶项,导致其模拟结果数值频散严重。Liu等[26]利用平面波原理推导了一种隐式的有限差分交错网格策略,提高了波场模拟的精度,但隐式格式需要求解更多的方程,难以用于逆时偏移;Yang等[27]在空间—波数域构建了新的频散关系,在最小二乘意义下,求得一阶交错网格差分系数的最优解,有效地压制了空间频散,但该方法比传统方法的时间频散更严重;雍鹏等[28]提出一种基于声波方程的等效交错网格方法,能同时压制空间频散和时间频散。

目前有限差分算法主要基于两类网格:规则网格和交错网格。对于弹性波数值模拟,交错网格往往具有更高的模拟精度[29-30],但目前大多数交错网格都基于一阶速度—应力方程,该方程含多个中间变量,在成像和反演过程中未被使用,因此浪费了计算机内存。等效交错网格[31]使用两次一阶交错网格离散二阶空间偏导数,其模拟精度等价于常规的一阶交错网格。因此,使用等效交错网格策略,就可以利用二阶位移方程的交错网格实现弹性波正演、成像和反演,与传统一阶速度—应力方程的交错网格相比,有效节约了计算内存。

另一方面,实现逆时偏移,需要保存全部时刻的震源波场,尤其是三维介质的情况下,需要消耗大量的存储,因此超额的存储量是制约逆时偏移用于巨量地震资料处理的瓶颈[32-33]。近年来许多学者对逆时偏移的存储问题进行了研究,并取得了一定的成效:Symes[34]利用少量的最优检验点控制波场的存储方式,检验点之间的波场通过重新正演获取,实现了计算量和存储量的折中;Clapp[35]提出一种边界存储策略,该方法只需保存边界区域的震源波场,反向传播时,利用边界区域的波场作为边界条件重构之前时刻的波场,该方法虽然增加了计算量,但大幅度降低了逆时偏移的存储需求。

本文基于优化差分系数的思想,提出了一种平面波优化差分系数的等效交错网格算法,主要目的在于利用低阶的优化差分算子替代传统高阶差分算子,在保证精度情况下,减少计算量。研究表明,本文的6、10阶空间差分算子分别能达到传统交错网格法10、18阶的精度,有效降低了计算量。在内存需求方面,将有效边界存储的策略引入弹性波逆时偏移,利用边界附近区域的部分波场值控制逆时偏移的波场存储方式,同时采用二阶位移弹性波方程,避免了中间变量,有效降低了内存占用。模型测试证实了本文方法相比于传统交错网格法在计算量和内存占用方面具有优势。

1 基于平面波优化的等效交错网格有限差分法 1.1 等效交错网格差分格式

在均匀各向同性弹性介质中,用位移表示的二维弹性波动方程为

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}u}}{{\partial {t^2}}} = v_{\rm{P}}^2\frac{{{\partial ^2}u}}{{\partial {x^2}}} + v_{\rm{S}}^2\frac{{{\partial ^2}u}}{{\partial {z^2}}} + (v_{\rm{P}}^2 - v_{\rm{S}}^2)\frac{{{\partial ^2}w}}{{\partial x\partial z}}}\\ {\frac{{{\partial ^2}w}}{{\partial {t^2}}} = v_{\rm{P}}^2\frac{{{\partial ^2}w}}{{\partial {z^2}}} + v_{\rm{S}}^2\frac{{{\partial ^2}w}}{{\partial {x^2}}} + (v_{\rm{P}}^2 - v_{\rm{S}}^2)\frac{{{\partial ^2}u}}{{\partial x\partial z}}} \end{array}} \right. $ (1)

式中:uw分别为在xz方向的位移分量;t表示波传播的时刻;vPvS分别为P波和S波速度。

以式(1)为例,说明等效交错网格有限差分法近似空间偏导的基本原理。二阶空间偏导数的离散格式为

$ \begin{array}{*{20}{l}} {D_{xx}^u = \frac{{{\partial ^2}u_{0,0}^0}}{{\partial {x^2}}} = \frac{1}{{\Delta {x^2}}}\sum\limits_{l = 1}^M {\sum\limits_{m = 1}^M {{a_l}} } {a_m}[(u_{m + l - 1,0}^0 - u_{ - m + l,0}^0) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (u_{m - l,0}^0 - u_{ - m - l + 1,0}^0)]} \end{array} $ (2)
$ \begin{array}{*{20}{l}} {D_{zz}^u = \frac{{{\partial ^2}u_{0,0}^0}}{{\partial {z^2}}} = \frac{1}{{\Delta {z^2}}}\sum\limits_{l = 1}^M {\sum\limits_{m = 1}^M {{b_l}} } {b_m}[(u_{0,m + l - 1}^0 - u_{0, - m + l}^0) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (u_{0,m - l}^0 - u_{0, - m - l + 1}^0)]} \end{array} $ (3)
$ \begin{array}{*{20}{l}} {D_{xz}^w = \frac{{{\partial ^2}w_{0,0}^0}}{{\partial x\partial z}} = \frac{1}{{\Delta x\Delta z}}\sum\limits_{l = 1}^M {\sum\limits_{m = 1}^M {{c_l}} } {c_m}[({w_{m - \frac{1}{2},l - \frac{1}{2}}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} w_{ - m + \frac{1}{2},l - \frac{1}{2}}^0 - (w_{m - \frac{1}{2},l + \frac{1}{2}}^0 - {w_{ - m + \frac{1}{2}, - l + \frac{1}{2}}})]} \end{array} $ (4)

式中:M为空间差分算子长度的一半;ul, m0为波场u在空间网格点(x=lΔxz=mΔz)处0时刻的值,Δx、Δz为空间采样间隔;a=(a1, a2, …, aM)、b=(b1, b2, …, bM)、c=(c1, c2, …, cM)为差分系数。

时间二阶导数的离散格式为

$ D_{tt}^u = \frac{{{\partial ^2}u_{0,0}^0}}{{\partial {t^2}}} = \frac{{u_{0,0}^1 - 2u_{0,0}^0 + u_{0,0}^{ - 1}}}{{\Delta {t^2}}} $ (5)

式中Δt为时间采样间隔。同理,可得其他时间和空间二阶导数DttwDzzwDxxwDxzu。值得注意的是,w分量定义在错开的半个网格点上。式(1)的等效交错网格有限差分离散格式为

$ \left\{ {\begin{array}{*{20}{l}} {D_{tt}^u = v_{\rm{P}}^2D_{xx}^u + v_{\rm{S}}^2D_{zz}^u + (v_{\rm{P}}^2 - v_{\rm{S}}^2)D_{xz}^w}\\ {D_{tt}^w = v_{\rm{P}}^2D_{zz}^w + v_{\rm{S}}^2D_{xx}^w + (v_{\rm{P}}^2 - v_{\rm{S}}^2)D_{xz}^u} \end{array}} \right. $ (6)
1.2 基于平面波优化的差分算子求取

若采用泰勒近似法求解差分系数,则a=b=c,均为传统一阶交错网格的差分系数,可由下式求解

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{1^0}}&{{3^0}}& \cdots &{{{(2M - 1)}^0}}\\ {{1^2}}&{{3^2}}& \cdots &{{{(2M - 1)}^2}}\\ {}&{}& \vdots &{}\\ {{1^{2M - 2}}}&{{3^{2M - 2}}}& \cdots &{{{(2M - 1)}^{2M - 2}}} \end{array}} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {1{a_1}}\\ {3{a_2}}\\ \vdots \\ {(2M - 1){a_M}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1\\ 0\\ \vdots \\ 0 \end{array}} \right] \end{array} $ (7)

实际上,泰勒近似的方法所确定的差分系数并非针对方程本身,而只是近似空间偏导数,然后代入波动方程求解。目前有不少的优化方法是针对方程本身构建频散关系,从而优化差分系数,然而它们大多基于声波方程[36],然后直接过渡到弹性波方程[37-38],造成了不必要的精度损失。

由式(1)可以看出,对于弹性波方程,空间偏导数与速度项耦合在一起,与声波方程不同的是,这些速度项不能被单独分离,因此用有限差分去近似空间偏导数时,P、S波速度影响了网格点权值的分配。为解决这一问题,本文定义abc,通过构建三种时空频散关系求解这三组差分系数。此外,声波方程的平面波解并不适用于弹性波方程,要想获取更为准确的差分系数,必须从弹性波方程本身出发,通过矩阵分解的方式求取式(1)的平面波解,然后构建新的时空频散关系,从而获得优化的差分算子。

1.2.1 弹性波的平面波解

二阶位移弹性波方程在波数域可表示为

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}\tilde u}}{{\partial {t^2}}} = - [v_{\rm{P}}^2k_x^2\tilde u + v_{\rm{S}}^2k_z^2\tilde u + (v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}\tilde w]}\\ {\frac{{{\partial ^2}\tilde w}}{{\partial {t^2}}} = - [v_{\rm{P}}^2k_z^2\tilde w + v_{\rm{S}}^2k_x^2\tilde w + (v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}\tilde u]} \end{array}} \right. $ (8)

式中(kx, kz)=(kcosθ, ksinθ)=k为波数矢量,其中θ为波的传播方向与水平面的夹角。将上式写成矩阵形式

$ \left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}\tilde u}}{{\partial {t^2}}}}\\ {\frac{{{\partial ^2}\tilde w}}{{\partial {t^2}}}} \end{array}} \right] = - \left[ {\begin{array}{*{20}{c}} {v_{\rm{P}}^2k_x^2 + v_{\rm{S}}^2k_z^2}&{(v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}}\\ {(v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}}&{v_{\rm{S}}^2k_x^2 + v_{\rm{P}}^2k_z^2} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\tilde u}\\ {\tilde W} \end{array}} \right] $ (9)

上式右端系数为实对称矩阵,根据实对称矩阵的性质可写为

$ \left[ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}\tilde u}}{{\partial {t^2}}}}\\ {\frac{{{\partial ^2}\tilde w}}{{\partial {t^2}}}} \end{array}} \right] = \mathit{\boldsymbol{P}}\left[ {\begin{array}{*{20}{l}} { - v_{\rm{P}}^2k}&{}\\ {}&{ - v_{\rm{S}}^2k} \end{array}} \right]{\mathit{\boldsymbol{P}}^{ - 1}}\left[ {\begin{array}{*{20}{l}} {\tilde u}\\ {\tilde \varpi } \end{array}} \right] $ (10)

式中

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{l}} {\frac{{{k_x}}}{k}}&{ - \frac{{{k_z}}}{k}}\\ {\frac{{{k_z}}}{k}}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{k_x}}}{k}} \end{array}} \right]}\\ {{\mathit{\boldsymbol{P}}^{ - 1}} = {\mathit{\boldsymbol{P}}^{\rm{T}}}} \end{array}} \right. $ (11)

式(10)的通解可以表示为

$ \begin{array}{l} \left[ {\begin{array}{*{20}{l}} {\tilde u(\mathit{\boldsymbol{k}},t)}\\ {\tilde w(\mathit{\boldsymbol{k}},t)} \end{array}} \right] = \mathit{\boldsymbol{P}}\left[ {\begin{array}{*{20}{l}} {{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt)}&{}\\ {}&{{\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)} \end{array}} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{P}}^{ - 1}}\left[ {\begin{array}{*{20}{l}} {\tilde u(\mathit{\boldsymbol{k}},0)}\\ {\tilde w(\mathit{\boldsymbol{k}},0)} \end{array}} \right] \end{array} $ (12)

式中$\tilde u$(k, 0)、$\tilde w$(k, 0)为初始值。将上式展开,可得弹性波动方程的平面波解

$ \left\{ \begin{array}{l} \tilde u(\mathit{\boldsymbol{k}},t)\\ \begin{array}{*{20}{c}} { = \frac{1}{{{k^2}}}[k_x^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) + k_z^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]\tilde u(\mathit{\boldsymbol{k}},0) + }\\ {\frac{1}{{{k^2}}}[{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) - {\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]{k_x}{k_z}\tilde w(\mathit{\boldsymbol{k}},0)} \end{array}\\ \tilde w(\mathit{\boldsymbol{k}},t)\\ \begin{array}{*{20}{c}} { = \frac{1}{{{k^2}}}[k_z^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) + k_x^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]\tilde w(\mathit{\boldsymbol{k}},0) + }\\ {\frac{1}{{{k^2}}}[{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) - {\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]{k_z}{k_x}\tilde u(\mathit{\boldsymbol{k}},0)} \end{array} \end{array} \right. $ (13)
1.2.2 时空频散关系

对式(6)左端项做时间Fourier变换,然后将式(13)代入,可得

$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{tt}^u] = \frac{2}{{{k^2}\Delta {t^2}}}\{ k_x^2[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - 1] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k_z^2[{\rm{cos}}({v_{\rm{S}}}k\Delta t) - 1]\} \tilde u(\mathit{\boldsymbol{k}},t) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2{k_x}{k_z}}}{{{k^2}\Delta {t^2}}}[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - {\rm{cos}}({v_{\rm{S}}}k\Delta t)]\tilde w(\mathit{\boldsymbol{k}},t)} \end{array} $ (14)
$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{tt}^w] = \frac{2}{{{k^2}\Delta {t^2}}}\{ k_z^2[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - 1] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k_x^2[{\rm{cos}}({v_{\rm{S}}}k\Delta t) - 1]\} \tilde w(\mathit{\boldsymbol{k}},t) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2{k_x}{k_z}}}{{{k^2}\Delta {t^2}}}[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - {\rm{cos}}({v_{\rm{S}}}k\Delta t)]\tilde u(\mathit{\boldsymbol{k}},t)} \end{array} $ (15)

同理,对式(6)右端项做空间Fourier变换,可得

$ {\rm{FT}}[D_{xx}^u] = - {\left\{ {\sum\limits_{m = 1}^M 2 {a_m}{\rm{sin}}[(m - 0.5){k_x}h]} \right\}^2}\tilde u(\mathit{\boldsymbol{k}},t) $ (16)
$ {\rm{FT}}[D_{zz}^u] = - {\left\{ {\sum\limits_{m = 1}^M 2 {b_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}^2}\tilde u(\mathit{\boldsymbol{k}},t) $ (17)
$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{xz}^w] = - \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\} \times }\\ {\quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}\tilde w(\mathit{\boldsymbol{k}},t)} \end{array} $ (18)
$ {\rm{FT}}[D_{zz}^w] = - {\left\{ {\sum\limits_{m = 1}^M 2 {a_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}^2}\tilde w(\mathit{\boldsymbol{k}},t) $ (19)
$ FT[D_{xx}^w] = - \left\{ {\sum\limits_{m = 1}^M 2 {b_m}{\rm{sin}}{{[(m - 0.5){k_x}h]}^2}} \right\}\tilde w(\mathit{\boldsymbol{k}},t) $ (20)
$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{xz}^u] = - \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_x}h]} \right\} \times }\\ {\quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}\tilde u(\mathit{\boldsymbol{k}},t)} \end{array} $ (21)

式中hxz为空间采样间隔。

联立式(14)~式(21),化简可得P、S和PS波时空频散关系

$ \begin{array}{l} {\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)\\ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{v_{\rm{P}}^2\Delta {t^2}}}{{{h^2}}}\left\{ {{{\left[ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}((m - 0.5){k_x}h)} \right]}^2} + } \right.}\\ {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left[ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}((m - 0.5){k_z}h)} \right]}^2}} \right\}} \end{array} \end{array} $ (22)
$ \begin{array}{l} {\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{S}}}k\Delta t)\\ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{v_{\rm{S}}^2\Delta {t^2}}}{{{h^2}}}\left\{ {{{\left[ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}((m - 0.5){k_x}h)} \right]}^2} + } \right.}\\ {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left[ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}((m - 0.5){k_z}h)} \right]}^2}} \right\}} \end{array} \end{array} $ (23)
$ \begin{array}{l} \frac{{{k_x}{k_z}}}{{{k^2}}}[{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t) - {\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{S}}}k\Delta t)]\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{(v_{\rm{P}}^2 - v_{\rm{S}}^2)\Delta {t^2}}}{{{h^2}}}\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_x}h]} \right\} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_z}h]} \right\} \end{array} $ (24)

对于式(22)~式(24),本文采用牛顿迭代法使其误差最小化,以求解差分系数abc

1.2.3 牛顿迭代法

式(22)~式(24)是非线性的,无法直接求解,只能依靠迭代逼近方法求解。首先以式(22)为例,利用二阶范数构建目标泛函

$ \varPhi (\mathit{\boldsymbol{a}}) = \frac{1}{2}\sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\| {\frac{{r_1^2q(\mathit{\boldsymbol{a}})}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}} - 1} \right\|_2^2} } $ (25)

式中:ke=π/h为最大优化波数;r1=vPΔt/h

$ \begin{array}{*{20}{l}} {q(\mathit{\boldsymbol{a}}) = {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}[(m - 0.5){k_x}h]} \right\}}^2} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}[(m - 0.5){k_z}h]} \right\}}^2}} \end{array} $ (26)

式(25)的牛顿迭代法解可写为

$ {\mathit{\boldsymbol{a}}^{i + 1}} = {\mathit{\boldsymbol{a}}^i} - {({\mathit{\boldsymbol{H}}^i})^{ - 1}}{\mathit{\boldsymbol{g}}^i} $ (27)

式中:i为迭代次数;H为Hessian矩阵;g为梯度向量。gH可通过对式(25)分别求取一阶和二阶导数获得

$ \begin{array}{l} \frac{{\partial \varPhi (\mathit{\boldsymbol{a}})}}{{\partial {a_m}}} = \sum\limits_{k = 0}^{{k_{\rm{e}}}} {\mathop \sum \limits_{\theta = 0}^{\pi /4} } \left[ {\frac{{r_1^2q(\mathit{\boldsymbol{a}})}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}} - 1} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{r_1^2}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}}\frac{{\partial q(\mathit{\boldsymbol{a}})}}{{\partial {a_m}}} \end{array} $ (28)
$ \begin{array}{l} \frac{{\partial \varPhi (\mathit{\boldsymbol{a}})}}{{\partial {a_m}\partial {a_n}}} = \sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\{ {\left[ {\frac{{r_1^2q(a)}}{{{{\sin }^2}(0.5{v_P}k\Delta t)}} - 1} \right] \times } \right.} } \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {\frac{{r_1^2}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}}\frac{{{\partial ^2}q(\mathit{\boldsymbol{a}})}}{{\partial {a_m}\partial {a_n}}} + }\\ {\left. {{{\left[ {\frac{{r_1^2}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}}} \right]}^2}\frac{{\partial q(\mathit{\boldsymbol{a}})}}{{\partial {a_n}}}\frac{{\partial q(\mathit{\boldsymbol{a}})}}{{\partial {a_m}}}} \right\}} \end{array} \end{array} $ (29)

为了加快收敛速度,将迭代初始值设为基于泰勒展开法的一阶交错网格差分系数,误差阈值设为1%,由于牛顿法的快速收敛性,通常只需要迭代几次即可收敛。因此,就可利用式(27)获得P波项的差分系数a。按照同样的方法,将vP换成vS即可获得差分系数b

当求解频散方程式(24)时,如果直接将左端的sin2(0.5vPkΔt)-sin2(0.5vSkΔt)除到方程的右端,则可能会使分母为零,进而导致迭代算法不收敛。为了解决这一问题,可将sin2(0.5vPkΔt)-sin2(0.5vSkΔt)拆成两项,然后分别与P波速度和S波速度匹配。构建的目标函数为

$ \begin{array}{*{20}{l}} {\varPhi (\mathit{\boldsymbol{c}}) = \frac{1}{2}\sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\| {\frac{{v_{\rm{P}}^2{k^2}\Delta {t^2}q(\mathit{\boldsymbol{c}})}}{{{h^2}{k_x}{k_z}{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}} - 1} \right\|_2^2 + } } }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\| {\frac{{v_{\rm{S}}^2{k^2}\Delta {t^2}q(\mathit{\boldsymbol{c}})}}{{{h^2}{k_x}{k_z}{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{S}}}k\Delta t)}} - 1} \right\|_2^2} } } \end{array} $ (30)

式中

$ \begin{array}{*{20}{l}} {q(\mathit{\boldsymbol{c}}) = \left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_x}h]} \right\} \times }\\ {\quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_z}h]} \right\}} \end{array} $ (31)

按照求解差分系数a的方法,以式(30)为目标泛函,可求得差分系数c

此外,通过特征值分析方法[39],可求得本文方法的稳定性条件为

$ \begin{array}{l} {\rm{max}}\left[ {{v_{{\rm{max}}}}\Delta t\sqrt {\frac{1}{{\Delta {x^2}}}\sum\limits_{m = 1}^M | {a_m}{|^2} + \frac{1}{{\Delta {z^2}}}\sum\limits_{m = 1}^M | {b_m}{|^2}} ,} \right.\\ {\kern 1pt} \left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {v_{\max }}\Delta t\sqrt {\left( {\frac{1}{{\Delta {x^2}}} + \frac{1}{{\Delta {z^2}}}} \right)\sum\limits_{m = 1}^M | {c_m}{|^2}} } \right] \le 1 \end{array} $ (32)
1.3 频散分析

定义P波、S波和PS波相速度误差分别为

$ {E_{\rm{P}}} = \frac{{v_{\rm{P}}^2\Delta {t^2}s(\mathit{\boldsymbol{a}})}}{{{h^2}{\rm{si}}{{\rm{n}}^2}(\pi f\Delta t)}} $ (33)
$ {{E_{\rm{S}}} = \frac{{v_{\rm{S}}^2\Delta {t^2}s(\mathit{\boldsymbol{b}})}}{{{h^2}{\rm{si}}{{\rm{n}}^2}(\pi f\Delta t)}}} $ (34)
$ {{E_{{\rm{PS}}}} = \frac{{v_{\rm{P}}^2\Delta {t^2}{s_1}(c) + v_{\rm{S}}^2\Delta {t^2}{s_2}(\mathit{\boldsymbol{c}})}}{{{h^2}{\rm{sin}}2\theta {\kern 1pt} {\rm{si}}{{\rm{n}}^2}(\pi f\Delta t)}}} $ (35)

式中

$ \begin{array}{*{20}{l}} {s(\mathit{\boldsymbol{a}}) = {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $ (36)
$ \begin{array}{*{20}{l}} {s(\mathit{\boldsymbol{b}}) = {{\left\{ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $ (37)
$ \begin{array}{*{20}{l}} {{s_1}(\mathit{\boldsymbol{c}}) = {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $ (38)
$ \begin{array}{*{20}{l}} {{s_2}(\mathit{\boldsymbol{c}}) = {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $ (39)

式(33)~式(35)中利用关系k=2πf/v将自变量设为频率f,根据采样定理可知f∈[0,vS/(2h)]。

图 1为均匀低速介质中传统方法和本文方法的三种波的频散误差曲线,其中模拟参数为:vP=2500m/s、vS=1443m/s、h=15m、M=5、Δt=0.5ms。可以发现,三种波的频散误差曲线形态不同,间接地证明了对于弹性波方程,有必要采取三组不同的差分系数单独地分析这三种波的频散。由图 1a可以看出:随着频率的增大,传统方法的频散误差值逐渐大于1,说明出现了时间频散,而本文的方法能够在一定程度上压制时间频散,但由于平均效应,会损失一定的空间精度;随着传播角度的增大,这种空间精度的损失会越来越小。由图 1b可以看出,S波产生了严重的空间频散,本文方法对S波的频散有很好的压制效果。图 1c为PS波频散误差曲线,其频散误差介于P波与S波之间,与传统方法相比,本文方法也能有效压制PS波频散。

图 1 低速介质中本文方法与传统方法的频散曲线对比 (a)P波;(b)S波;(c)PS波

为了验证本文方法在高速介质中的优势,将速度增大为:vP=4000m/s、vS=2309m/s,其余参数不变,计算三种波的频散曲线(图 2)。可以看出,在高速介质中,三种波都出现了很严重的时间频散,P波时间频散最为严重,随着频率的增大,频散误差越来越大。与传统方法相比,本文方法能够有效压制时间频散,尤其是在高频段。另外,对比图 2图 1可以看出,图 2的有效频带范围要宽于图 1(红圈所示),说明随着速度的增大,有效频带范围变宽。

图 2 高速介质中本文方法与传统方法的频散曲线对比 (a)P波;(b)S波;(c)PS波
2 弹性波逆时偏移成像策略 2.1 矢量波分离和成像条件

由于弹性波中的P波和S波是耦合在一起的,在成像的过程会存在串扰假象,严重影响成像剖面的质量,因此在成像之前需要对弹性波进行分离。为了便于推导,将式(1)写成矢量形式

$ \frac{{{\partial ^2}\mathit{\boldsymbol{u}}}}{{\partial {t^2}}} = v_{\rm{P}}^2\nabla (\nabla \cdot \mathit{\boldsymbol{u}}) + v_{\rm{S}}^2\nabla \times \nabla \times \mathit{\boldsymbol{u}} $ (40)

式中u=(u, w)表示位移矢量。P波与S波的耦合关系可以表示为

$ \mathit{\boldsymbol{u}} = \mathit{\boldsymbol{p}} + \mathit{\boldsymbol{s}} $ (41)

式中:p表示P波分量;s为S波分量。将式(41)代入式(40),有

$ \frac{{{\partial ^2}\mathit{\boldsymbol{p}}}}{{\partial {t^2}}} + \frac{{{\partial ^2}\mathit{\boldsymbol{s}}}}{{\partial {t^2}}} = v_{\rm{P}}^2\nabla (\nabla \cdot \mathit{\boldsymbol{u}}) + v_{\rm{S}}^2\nabla \times \nabla \times \mathit{\boldsymbol{u}} $ (42)

由于P波是无旋场,S波是无散场,因此对式(42)两边同时求散度和旋度,可得

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}\mathit{\boldsymbol{p}}}}{{\partial {t^2}}} = v_{\rm{P}}^2\nabla (\nabla \cdot \mathit{\boldsymbol{u}})}\\ {\frac{{{\partial ^2}\mathit{\boldsymbol{s}}}}{{\partial {t^2}}} = v_{\rm{S}}^2\nabla \times \nabla \times \mathit{\boldsymbol{u}}} \end{array}} \right. $ (43)

u=px+sxw=pz+sz,将式(43)展开为笛卡尔坐标(xz)的形式,有

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}{p_x}}}{{\partial {t^2}}} = v_{\rm{P}}^2\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right)}\\ {\frac{{{\partial ^2}{p_z}}}{{\partial {t^2}}} = v_{\rm{P}}^2\left( {\frac{{{\partial ^2}w}}{{\partial {z^2}}} + \frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right)}\\ {\frac{{{\partial ^2}{s_x}}}{{\partial {t^2}}} = v_{\rm{S}}^2\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}} - \frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right)}\\ {\frac{{{\partial ^2}{s_z}}}{{\partial {t^2}}} = v_{\rm{S}}^2\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}} - \frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right)} \end{array}} \right. $ (44)

基于上式,就可以利用有限差分法外推P波和S波两个方向的波场,再利用内积成像公式提取PP波、PS波成像值

$ \left\{ \begin{array}{l} {I_{{\rm{PP}}}}({\rm{ }}\mathit{\boldsymbol{x}}{\rm{ }}){\rm{ }} = \int {[{p_x}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){p_x}({\mathit{\boldsymbol{x}}_{\rm{r}}},t) + {p_z}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){p_z}({\mathit{\boldsymbol{x}}_{\rm{r}}},t)]{\rm{d}}t} \\ {I_{{\rm{PS}}}}({\rm{ }}\mathit{\boldsymbol{x}}{\rm{ }}){\rm{ }} = \int {[{p_x}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){s_x}({\mathit{\boldsymbol{x}}_{\rm{r}}},t) + {p_z}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){s_z}({\mathit{\boldsymbol{x}}_{\rm{r}}},t)]{\rm{d}}t} {\rm{ }} \end{array} \right. $ (45)

式中:xs为炮点坐标;xr为检波点坐标。由于地震勘探主要使用P波震源,本文没有给出SP波和SS波成像过程。

2.2 有效边界存储策略

逆时偏移技术在求解波动方程的时候需要存储全部时刻的震源波场和接收波场,需大量的内存而难以处理巨量数据。将平衡计算量和存储量的有效边界存储策略引入到弹性波逆时偏移中,可有效节约存储空间,如图 3所示。值得注意的是,为了保证空间差分算子长度的边界区域波场值的正确性,需要进行两次波场反传,因此会增加部分计算量。

图 3 有效边界存储策略示意图 红框内:计算区域;红框外:边界区域;黄色区域:有效边界区域

边界存储的基本做法是:先进行一次震源波场的正向延拓,但只存储边界区域波场,由于接收波场的重构是震源波场的反过程,因此在反传时,可以将边界区域的波场作为边界条件重构之前时刻的波场,这样会增加一些计算量,但大大地降低了存储量,而且边界区域的波场无需重新计算,直接取出替换即可,有效缓解了计算负担。

研究发现,要重构反传波场,无需存储全部边界区域的波场,只需要保证部分边界区域的波场值是正确的,就可以进行波场反向延拓,即有效边界存储策略,该区域的厚度与有限差分空间阶数有关。如采用10阶空间有限差分,根据传统有限差分法的延拓公式,只需要保证边界区域厚度为10以及计算区域边界厚度为10的区域波场即可,如图 3黄色区域所示。

如果不采用边界存储策略,弹性波逆时偏移需要存储量为

$ {S_1} = ({N_x} + 2 \times {N_{\rm{b}}}) \times ({N_z} + 2 \times {N_{\rm{b}}}) \times {N_t} \times 4{\rm{B}} $ (46)

如果采用边界存储策略,需要的存储量为

$ {S_2} = 2 \times {N_{\rm{b}}} \times ({N_x} + {N_z} + 2 \times {N_{\rm{b}}}) \times {N_t} \times 4{\rm{B}} $ (47)

如采用有效边界存储策略,需要的存储量为

$ {S_3} = 4 \times N \times ({N_x} + {N_z}) \times {N_t} \times 4{\rm{B}} $ (48)

式中:NxNz分别为横纵向网格点数;Nt为时间采样点数;Nb为边界层厚度;N=2M为有限差分空间阶数。假设模型尺寸为201×201,时间样点数为1000,边界层厚度为30,有限差分阶次为10,则三种方法所需的存储量分别为0.254、0.103和0.060GB。当模型更大时,有效边界存储的优势会更加明显。

3 模型试算 3.1 洼陷模型

首先选择比较简单的洼陷模型进行测试,详细讨论本文方法在计算量、内存需求和模拟精度三个方面的优劣。模型尺寸为301×301,时间采样间隔为0.5ms,空间采样间隔为15m,总接收时长为3.0s,采用全接收方式,布置14炮,炮间距为300m,道间距为15m,纵横波速度比为1.73,P波速度模型如图 4所示。选用PML边界条件,厚度为30个空间采样间隔。采用主频为30Hz的Ricker子波进行波场模拟。

图 4 洼陷P波速度模型

图 5展示了传统一阶交错网格法(简称传统方法)与本文平面波优化的等效交错网格法(简称本文方法)水平分量单炮地震记录,有限差分空间阶数均为10,可见本文方法能有效压制数值频散(图 5红框所示)。

图 5 洼陷模型两种方法单炮记录对比 (a)传统方法;(b)本文方法

图 6a图 6b分别为传统方法和本文方法的PP波成像剖面,可以看出,本文方法PP波成像剖面的偏移噪声能量更弱(黑色箭头所示),且对同相轴的刻画更精细(红框所示)。图 6c图 6d为两种方法的PS波成像剖面,实施了与图 6a图 6b相同的增益和滤波处理,与PP波成像结果相比,PS波成像剖面存在更为严重的噪声和假象,这与频散分析结论一致,但PS波成像分辨率要高于PP波成像。对比图 6c图 6d可见,本文方法的PS成像剖面质量要优于传统方法(黑色箭头和红圈所示)。

图 6 洼陷模型两种方法成像剖面对比 (a)传统方法,PP波成像;(b)本文方法,PP波成像;(c)传统方法,PS波成像;(d)本文方法,PS波成像

为了验证本文方法利用低阶的空间差分能实现传统方法高阶差分的成像效果,抽取两种方法、不同差分阶次的PP波(图 7)、PS波(图 8)成像单道进行对比,并计算与参考道的残差

图 7 洼陷模型x=3.0km处两种方法、不同阶次有限差分逆时偏移PP波成像单道对比 (a)本文方法,2M=6;(b)本文方法,2M=10;(c)传统方法,2M=10;(d)传统方法,2M=18。图中右下部分是矩形框中局部放大

图 8 洼陷模型x=3.0km处两种方法、不同阶次有限差分逆时偏移PS波成像单道对比 (a)本文方法,2M=6;(b)本文方法,2M=10;(c)传统方法,2M=10;(d)传统方法,2M=18。图中右下部分是矩形框中局部放大
$ E = \sum\limits_{n = 1}^{{N_z}} {(|\mathit{\boldsymbol{I}}_n^{{\rm{cal}}} - \mathit{\boldsymbol{I}}_n^{{\rm{ref}}}|)} $ (49)

式中:Ical为PP波或PS波成像结果;Iref为参考道,是传统方法差分算子长度为20时的逆时偏移结果(接近无频散)。

图 7图 8可见,本文方法6阶空间差分能达到传统方法10阶空间差分的成像效果,而10阶空间差分则能达到传统方法18阶的成像效果,因此在相同的精度要求下,本文方法使用的差分阶次低,能有效节约计算成本。

以有效边界存储策略为例,表 1统计了本文方法和传统方法采用不同阶次有限差分进行洼陷模型逆时偏移的用时,可见:差分算子阶次越高,计算效率越低;如利用6阶本文方法代替10阶传统方法,计算效率可提升约37%;如利用10阶本文方法代替18阶传统方法,计算效率可提升约31%。对比洼陷模型不同阶次、不同方法、采用不同存储策略的内存占用(表 1),可见有效边界存储策略能节省大量内存,且空间差分阶数越低节省得越多。

表 1 洼陷模型逆时偏移不同方法内存占用和耗时对比
3.2 大庆模型

选取大庆油田A区块速度模型(图 9a)验证本文方法对复杂模型的适用性。模型尺寸为1250×248,Δx=12m,Δz=8m,纵横波速度比为1.73。震源选取主频为30Hz的Ricker子波,Δt=0.5ms,有限差分算子长度为10。观测系统采用全接收的方式布置51炮,炮间距为240m,道间距为12m,接收时长为4.0s。

图 9 大庆P波速度模型(a)及传统方法(b)与本文方法(c)模拟单炮记录对比

图 9b图 9c分别为传统方法和本文方法模拟的单炮记录,图 10a图 10b分别为图 9b图 9c红框中的局部放大。对比图 10a图 10b红色箭头处可以看出,传统方法产生了严重的数值频散,本文方法利用更优的差分算子,能有效压制数值频散。

图 10 大庆模型不同方法模拟单炮记录局部(图 9b图 9c方框部分)放大对比 (a)传统方法;(b)本文方法

图 11为大庆模型本文方法弹性波逆时偏移结果。可以看出,本文方法能较好地展示成像剖面上的构造层位,对孔洞的几何边界刻画也较为精确,说明了本文方法对复杂模型的有效性。

图 11 大庆模型本文方法弹性波逆时偏移剖面 (a)PP波;(b)PS波

抽取6、10阶本文方法成像结果与10、18阶传统方法进行单道对比,并计算其残差(图 12图 13)。可见,6和10阶本文方法成像精度分别与10、18阶传统方法相当,意味着在相同的精度要求下,利用本文方法进行波场外推所需的差分算子更短,因此能有效降低节约计算量。

图 12 大庆模型x=9.6km处两种方法、不同阶次有限差分逆时偏移PP波成像单道对比 (a)本文方法,2M=6;(b)本文方法,2M=10;(c)传统方法,2M=10;(d)传统方法,2M=18图中左下部分是矩形框中局部放大;参考道是20阶传统方法成像结果

图 13 大庆模型x=9.6km处两种方法、不同阶次有限差分逆时偏移PS波成像单道对比 (a)本文方法,2M=6;(b)本文方法,2M=10;(c)传统方法,2M=10;(d)传统方法,2M =18图中下部是矩形框中局部放大;参考道是20阶传统方法成像结果
4 结论

(1) 本文通过推导二阶弹性波方程的平面波解,构建了新的时空频散关系,利用牛顿迭代法,获得了优化的差分系数。平面波优化的差分算子在低速介质中能压制空间频散,在高速介质中能有效压制时间频散。与一阶速度—应力方程相比,避免了中间变量的出现,因此能减少正演和偏移的内存占用。

(2) 将本文提出的平面波优化的差分算子应用于有效边界存储弹性波逆时偏移,有效地降低了存储需求,减小了计算量。

(3) 洼陷模型和复杂的大庆模型弹性波逆时偏移结果表明,6、10阶本文方法能达到10、18阶传统方法的成像精度,因此能有效减少计算量。

参考文献
[1]
Hemon C. Equations d'onde et modeles[J]. Geophy-sical Prospecting, 1978, 26(4): 790-821. DOI:10.1111/j.1365-2478.1978.tb01634.x
[2]
王娟, 李振春, 孙小东, 等. TTI介质逆时偏移成像[J]. 石油地球物理勘探, 2012, 47(4): 573-577.
WANG Juan, LI Zhenchun, SUN Xiaodong, et al. Reverse time migration in tilt transversely isotropic (TTI) media[J]. Oil Geophysical Prospecting, 2012, 47(4): 573-577.
[3]
Baysal E, Kosloff D D, Sherwood J W C. Reverse-time migration[J]. Geophysics, 1983, 48(11): 1514-1524. DOI:10.1190/1.1441434
[4]
Larner K, Beasley C. Cascaded migrations:Improving the accuracy of finite-difference migration[J]. Geophysics, 1987, 52(5): 618-643. DOI:10.1190/1.1442331
[5]
Du Q Z, Gong X F, Zhang M Q, et al. 3D PS-wave i-maging with elastic reverse-time migration[J]. Geophysics, 2014, 79(5): S173-S184. DOI:10.1190/geo2013-0253.1
[6]
Yong P, Huang J P, Liao W Y, et al. Elastic-wave reverse-time migration based on decoupled elastic-wave equations and inner-product imaging condition[J]. Journal of Geophysics and Engineering, 2016, 13(6): 953-963. DOI:10.1088/1742-2132/13/6/953
[7]
Li Y Y, Du Y, Yang J D, et al. Elastic reverse time migration using acoustic propagators[J]. Geophysics, 2018, 83(5): S399-S408. DOI:10.1190/geo2017-0687.1
[8]
Whitmore N D.Iterative depth migration by backward time propagation[C].SEG Technical Program Expanded Abstracts, 1983, 2: 382-385.
[9]
Mcmechan G A. Migration by extrapolation of time-dependent boundary values[J]. Geophysical Prospecting, 2006, 31(3): 413-420.
[10]
柯璇, 石颖. 基于一步法波场延拓的正演模拟和逆时偏移成像[J]. 地球物理学报, 2017, 60(11): 4468-4479.
KE Xuan, SHI Ying. Forward simulation and reverse time migration imaging based on the step wavefield extrapolation[J]. Chinese Journal of Geophysics, 2017, 60(11): 4468-4479.
[11]
严红勇, 刘洋. 基于时空域自适应高阶有限差分的声波叠前逆时偏移[J]. 地球物理学报, 2013, 56(3): 971-984.
YAN Hongyong, LIU Yang. Acoustic prestack reverse time migration using the adaptive high-order finite-difference method in time-space domain[J]. Chinese Journal of Geophysics, 2013, 56(3): 971-984.
[12]
Etgen J T.A tutorial on optimizing time domain finite-difference schemes: "Beyond Holberg"[R].Stanford Exploration Project, 2007, 129: 33-43.
[13]
Dablain M A. The application of high-order differencing to the scalar wave equation[J]. Geophysics, 1986, 51(1): 54-66.
[14]
张志禹, 谭显波, 黄璐瑶, 等. 抗频散有限差分波动方程数值模拟及逆时偏移[J]. 石油地球物理勘探, 2014, 49(6): 1115-1121.
ZHANG Zhiyu, TAN Xianbo, HUANG Luyao, et al. Anti-dispersion finite difference simulation and reverse-time migration for wave equations[J]. Oil Geophysical Prospecting, 2014, 49(6): 1115-1121.
[15]
朱峰, 黄建平, 李振春, 等. 基于构造滤波的VTI介质平面波最小二乘FFD叠前深度偏移[J]. 石油地球物理勘探, 2018, 53(5): 932-944.
ZHU Feng, HUANG Jianping, LI Zhenchun, et al. Structure-constrained plane-wave least-squares FFD prestack depth migration for VTI media[J]. Oil Geophysical Prospecting, 2018, 53(5): 932-944.
[16]
刘红伟, 李博, 刘洪, 等. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 2010, 53(7): 1725-1733.
LIU Hongwei, LI Bo, LIU Hong, et al. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J]. Chinese Journal of Geophysics, 2010, 53(7): 1725-1733.
[17]
Tan S, Huang L. A staggered-grid finite-difference scheme optimized in the time-space domain for mode-ling scalar-wave propagation in geophysical problems[J]. Journal of Computational Physics, 2014, 276(1): 613-634.
[18]
Liu H F, Dai N X, Niu F L, et al. An explicit time evolution method for acoustic wave propagation[J]. Geophysics, 2014, 79(3): T117-T124.
[19]
Fei T, Larner K. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport[J]. Geophysics, 1995, 60(6): 1830-1842. DOI:10.1190/1.1443915
[20]
Chen S, Yang D H, Deng X Y. A weighted Runge-Kutta method with weak numerical dispersion for solving wave equations[J]. Communications in Computational Physics, 2010, 7(5): 1027-1048. DOI:10.4208/cicp.2009.09.088
[21]
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(1): 141-159. DOI:10.1785/0120100041
[22]
Wang Y F, Liang W Q, Nashed Z, et al. Seismic mo-deling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion relationship preserving method[J]. Geophy-sics, 2014, 79(5): T277-T285.
[23]
Ren Z M, Liu Y. Acoustic and elastic modeling by optimal time-space-domain staggered-grid finite-diffe-rence schemes[J]. Geophysics, 2015, 80(1): T17-T40.
[24]
Liu Y. Globally optimal finite-difference schemes based on least squares[J]. Geophysics, 2013, 78(4): T113-T132.
[25]
Ren Y J, Huang J P, Yong P, et al. Window functions and optimized staggered-grid finite-difference operators[J]. Applied Geophysics, 2018, 15(2): 253-260. DOI:10.1007/s11770-018-0668-7
[26]
Liu Y, Sen M K. An implicit staggered-grid finite-difference method for seismic modelling[J]. Geophysical Journal International, 2009, 179(1): 459-474. DOI:10.1111/j.1365-246X.2009.04305.x
[27]
Yang L, Yan H Y, Liu H. Least squares staggered-grid finite-difference for elastic wave modelling[J]. Exploration Geophysics, 2014, 45(4): 255-260. DOI:10.1071/EG13087
[28]
雍鹏, 黄建平, 李振春, 等. 优化的时空域等效交错网格有限差分正演模拟[J]. 中国石油大学学报(自然科学版), 2017, 41(6): 72-79.
YONG Peng, HUANG Jianping, LI Zhenchun, et al. Forward modeling by optimized equivalent staggered-grid finite-difference method for time-space domain[J]. Journal of China University of Petroleum(Edition of Natural Science), 2017, 41(6): 72-79.
[29]
Yong P, Huang J P, Li Z C, et al. Optimized equivalent staggered-grid FD method for elastic wave mo-deling based on plane wave solutions[J]. Geophysical Journal International, 2016, 208(2): 1157-1172.
[30]
Virieux J. SH wave propagation in heterogeneous media:velocity-stress finite-difference method[J]. Geophysics, 1984, 49(11): 1933-1957. DOI:10.1190/1.1441605
[31]
BartoloL D, Dors C, Mansur W J. A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation[J]. Geophysics, 2012, 77(5): T187-T199. DOI:10.1190/geo2011-0345.1
[32]
杨仁虎, 凌云, 瞿立建, 等. 基于边界存储的地震波场重构及逆时偏移成像应用研究[J]. 地球物理学进展, 2017, 32(3): 354-357.
YANG Renhu, LING Yun, QU Lijian, et al. Seismic wavefield reconstruction based on boundary storage and application of reverse time migration[J]. Progress in Geophysics, 2017, 32(3): 354-357.
[33]
李青阳, 吴国忱, 段沛然. 准规则网格高阶有限差分法非均质弹性波波场模拟[J]. 石油地球物理勘探, 2019, 54(3): 539-550.
LI Qingyang, WU Guochen, DUAN Peiran. Elastic wavefield forward modeling in heterogeneous media based on the quasi-regular grid high-order finite difference[J]. Oil Geophysical Prospecting, 2019, 54(3): 539-550.
[34]
Symes W W. Reverse time migration with optimal checkpointing[J]. Geophysics, 2007, 72(5): SM213-SM221. DOI:10.1190/1.2742686
[35]
Clapp R G.Reverse time migration: Saving the boun-daries[R].Stanford Exploration Project, 2008, 136: 136-144.
[36]
Finkelstein B, Kastner R. Finite difference time do-main dispersion reduction schemes[J]. Journal of Computational Physics, 2007, 221(1): 422-438.
[37]
Graves R W. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite diffe-rences[J]. Bulletin of the Seismological Society of America, 1996, 86(4): 1091-1106.
[38]
Xia F, Dong L G, Ma Z T. The numerical modeling of 3-D elastic wave equation using a high-order staggered-grid finite difference scheme[J]. Applied Geophysics, 2004, 1(1): 38-41.
[39]
牟永光, 裴正林. 三维复杂介质地震数值模拟[M]. 北京: 石油工业出版社, 2005.