全文快速搜索:   高级搜索

  中国石油大学学报(自然科学版)  2017, Vol. 41 Issue (6): 71-79  DOI:10.3969/j.issn.1673-5005.2017.06.008
0

引用本文 [复制中英文]

雍鹏, 黄建平, 李振春, 等. 优化的时空域等效交错网格有限差分正演模拟[J]. 中国石油大学学报(自然科学版), 2017, 41(6): 71-79. DOI: 10.3969/j.issn.1673-5005.2017.06.008.
[复制中文]
YONG Peng, HUANG Jianping, LI Zhenchun, et al. Forward modeling by optimized equivalent staggered-grid finite-differencemethod for time-space domain[J]. Journal of China University of Petroleum (Edition of Natural Science), 2017, 41(6): 71-79. DOI: 10.3969/j.issn.1673-5005.2017.06.008.
[复制英文]

基金项目

山东省重大创新工程(2017CXGC1608,2017CXGC1602);中国科学院战略性先导科技专项(XDA14010303);国际合作重点项目(41720104006);国家油气重大专项课题(2016ZX05006-004,2016ZX05014001,2016ZX05002);中央高校基本科研业务费专项(R1701023A)

作者简介

雍鹏(1992-), 男, 博士研究生, 研究方向为地震波正演及波形反演。E-mail: yongpeng1991@sina.com

通讯作者

黄建平(1982-), 男, 教授, 博士, 研究方向为地震波传播与成像技术。E-mail: jphuang@upc.edu.cn

文章历史

收稿日期:2017-02-24
优化的时空域等效交错网格有限差分正演模拟
雍鹏 , 黄建平 , 李振春 , 李庆洋 , 刘培君 , 杨明伟 , 袁双齐 , 任英俊     
中国石油大学地球科学与技术学院, 山东青岛 266580
摘要: 地震波有限差分数值模拟中, 使用优化的差分系数可以有效地减少数值频散(相速度误差)。为了同时压制空间频散与时间频散, 提出一种时空优化的等效交错网格有效差分正演模拟方法。采用以减少给定波数范围内所有传播方向的相对频散误差为目标函数, 并借助牛顿法快速获取优化的差分系数。频散分析表明等效交错网格正演在低速介质中容易出现空间频散, 在高速介质中高波数分量容易出现时间频散, 数值模拟显示该方法既能有效地压制空间频散, 也可以减少时间频散。Marmousi模型测试表明该方法适用于复杂地质模型的地震波场正演模拟。相比传统差分方法, 该时空优化差分方法在不增加计算量的前提下, 可以有效地提高数值模拟精度。
关键词: 数值频散    时空优化差分系数    等效交错网格    地震波    正演模拟    牛顿法    
Forward modeling by optimized equivalent staggered-grid finite-differencemethod for time-space domain
YONG Peng , HUANG Jianping , LI Zhenchun , LI Qingyang , LIU Peijun , YANG Mingwei , YUAN Shuangqi , REN Yingjun     
School of Geosciences in China University of Petroleum, Qingdao 266580, China
Abstract: Optimized finite-difference (FD) coefficients can effectively suppress numerical dispersion errors (phase velocity errors) in seismic wave modeling. Most optimized FD schemes used to solve seismic wave equation suppress only spatial but not temporal dispersion errors. In this paper, we develop an optimized FD scheme with high spatial and temporal accuracy for numerical scalar-wave modeling based on equivalent staggered grids. The final aim of optimized FD coefficients is to reduce phase velocity errors. We design an objective function for minimizing relative temporal and spatial dispersion errors of waves propagating in all directions. In addition, Newton method is applied to quickly solve this nonlinear problem. Dispersion analysis has demonstrated that space dispersion easily appear in the low-velocity media, while time dispersion has more probability to occur in the high-velocity media. Modeling tests indicate that the proposed method ensures high accuracy not only in the time domain but also in the space domain. Compared with the traditional methods, the proposed method can improve numerical accuracy without sacrificing computation efficiency.
Keywords: numerical dispersion    time-space domain optimized finite-difference operators    equivalent staggered grid finite-difference    seismic wave    forward modeling    Newton method    

波动方程数值模拟是地震波传播与成像研究的重要组成部分, 有限差分因其计算速度快、内存占用小以及易于实现的特点而被广泛应用于地震波正演模拟[1-4]、偏移成像[5-7]及波形反演[8-9]中。有限差分数值模拟中, 空间导数项通常采用高阶差分进行逼近; 考虑到计算内存, 时间导数项通常采用二阶差分。通过在空间域进行泰勒展开获取差分系数的方法得到了广泛应用[10-11], 虽然这样获取的系数在低波数段具有高精度, 但随着波数的增加精度会快速降低[12]。为了降低数值频散(相速度误差), 国内外学者进行了大量的研究, 主要包括空间差分优化, 时空差分同时优化以及网格节点优化3个方面。空间差分优化方面:Holberg[13]通过在一定波数范围内控制空间频散误差获取优化的空间差分系数; Chu等[14]使用二次窗函数来减少空间频散; Zhang等[15-16]采用模拟退火算法给出误差上限求取空间或时间差分系数。使用有限差分法求解波动方程时, 既对空间偏导进行逼近又对时间偏导进行差分近似。时空差分同时优化方面:Etgen[17]以直接减少相速度误差为目标函数, 通过迭代算法获取时空同时优化的中心差分系数; Liu[18]研究中心差分优化时发现以相对时空频散误差为目标函数也可达到同样的精度; Ren等[19]以相对时空频散误差为目标函数并使用最小二乘算法对弹性波差分系数展开了时空优化。一些学者也从网格节点选择的角度开展了研究:Liu等[20]对二阶常密度声波方程进行差分模拟时, 指出在传统网格节点基础上增加少量网格节点可有效地提高时间差分精度; Tan等[21]进一步将这种网格节点应用到变密度交错网格正演模拟中。有限差分声波模拟中, 中心差分通常应用于二阶常密度声波方程, 当介质存在密度变化时, 常使用一阶压强-速度方程并采用交错网格进行模拟。事实上, 变密度二阶声波方程也可采用交错网格进行正演模拟, 并可达到与一阶声波方程交错网格相同的精度, Bartolo等[22]称之为等效交错网格。针对等效交错网格正演模拟, Liu等[23]依据时空频散关系固定特定传播角度, 并在波数域进行泰勒展开得到时空域优化的差分系数, 实现了对多个传播角度的优化。Wang等[24]以减少时空频散差值误差为目标函数, 并使用梯度类算法求取所有传播方向优化的系数。此外, 他们的研究也指出当优化的波数范围较小时可能存在不稳定问题, 故将正则化算子引入系数求取中。由于直接以减少相速度误差为目标函数在求解过程中可能存在不稳定, 笔者首先对频散差值误差、相速度误差以及相对频散误差的函数性质进行对比分析, 选择以减小相对时空频散误差为目标函数, 考虑到优化的波数范围较小时可能存在不稳定问题, 在确定展开优化的波数区间的同时考虑子波的频谱与传统差分系数下的误差函数本身的变化情况后, 使用牛顿法快速求解时空优化的差分系数。

1 等效交错网格正演模拟原理 1.1 等效交错网格差分格式

首先考虑三维常数介质情况下二阶声波方程:

$ \frac{\partial }{{\partial x}}\left( {\frac{{\partial p}}{{\rho \partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {\frac{{\partial zp}}{{\rho \partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {\frac{{\partial p}}{{\rho \partial z}}} \right) = \frac{1}{K}\frac{{{\partial ^2}p}}{{\partial {t^2}}}. $ (1)

其中

$ K = \lambda + 2\mu = \rho {v^2}. $

式中, v为介质速度; ρ为介质密度; p(x, t)为声波压强场; t为时间; x=(x, y, z)为空间位置。对沿着xyz方向的空间一阶偏导数可采用高阶交错网格进行逼近:

$ \begin{array}{l} \frac{{\partial p}}{{\partial x}} \approx \frac{1}{{\Delta x}}\sum\limits_{m = 1}^M {{a_m}\left( {p_{m - 1/2,0,0}^0 - p_{ - m + 1/2,0,0}^0} \right)} ,\\ \frac{{\partial p}}{{\partial y}} \approx \frac{1}{{\Delta y}}\sum\limits_{m = 1}^M {{a_m}\left( {p_{0,m - 1/2,0}^0 - p_{0, - m + 1/2,0}^0} \right)} ,\\ \frac{{\partial p}}{{\partial z}} \approx \frac{1}{{\Delta z}}\sum\limits_{m = 1}^M {{a_m}\left( {p_{0,0,m - 1/2}^0 - p_{0,0, - m + 1/2}^0} \right)} . \end{array} $ (2)

其中

$ p_{i,j,l}^n = p\left( {x + i\Delta x,y + i\Delta y,z + l\Delta z,t + n\Delta t} \right). $

式中, Δx、Δy和Δzxyz方向空间采样间隔; ijl为对应方向网格节点的编号; Δt为时间步长; M为差分算子长度的1/2;am为差分系数。考虑计算内存, 二阶时间导数通常采用二阶差分进行逼近:

$ \frac{{{\partial ^2}p}}{{\partial {t^2}}} = \frac{{p_{0,0,0}^1 - 2p_{0,0,0}^0 + p_{0,0,0}^{ - 1}}}{{\Delta {t^2}}}. $ (3)

本文中将空间采样步长统一为h。将公式(2)与公式(3)代入波动方程(1), 化简可得[24]

$ \begin{array}{l} p_{0,0,0}^1 - 2p_{0,0,0}^0 + p_{0,0,0}^{ - 1} = \\ {\rho _{0,0,0}}r_{0,0,0}^2\sum\limits_{m = 1}^N {\sum\limits_{n = 1}^N {{a_m}{a_n}\left[ {\frac{{\left( {p_{m + n - 1,0,0}^0 - p_{m - n,0,0}^0} \right)}}{{{\rho _{m - 1/2,0,0}}}} - } \right.} } \\ \frac{{\left( {p_{m + n,0,0}^0 - p_{m - n + 1,0,0}^0} \right)}}{{{\rho _{ - m + 1/2,0,0}}}} + \frac{{\left( {p_{0,m + n - 1,0}^0 - p_{0,m - n,0}^0} \right)}}{{{\rho _{0,m - 1/2,0}}}} - \\ \frac{{\left( {p_{0,m + n,0}^0 - p_{0,m - n + 1,0}^0} \right)}}{{{\rho _{0, - m + 1/2,0}}}} + \frac{{\left( {p_{0,0,m + n - 1}^0 - p_{0,0,m - n}^0} \right)}}{{{\rho _{0,0,m - 1/2}}}} - \\ \left. {\frac{{\left( {p_{0,0,m + n}^0 - p_{0,0,m - n + 1}^0} \right)}}{{{\rho _{0,0, - m + 1/2}}}}} \right]. \end{array} $ (4)

$ r = \frac{{v\Delta t}}{h} $为CFL(Courant-Friedrichs-Lewy)计算常数。传统的交错网格正演模拟建立在一阶速度-压强方程基础上, 通常压强定义在整网格点处, 偏振速度定义在半网格点上。然而, 等效交错网格建立在二阶声波方程基础上, 尽管使用交错差分逼近一阶偏导, 由于二阶偏导需要两次使用交错差分, 经过网格节点重新调整, 参与计算的网格节点均位于整网格点上。Bartolo等[22]证明等效交错网格与传统的交错网格有相同的计算精度, 但等效交错网格在二维、三维情况下可分别节省33%、50%的计算内存。

1.2 优化的交错网格算子

将平面波解$ p_{m, n, l}^{k\Delta t} = \exp \left\{ {{\rm{i}}\left[{\left( {x + m\Delta x} \right) + \left( {y + n\Delta y} \right) + \left( {z + l\Delta z} \right)} \right]} \right\}\exp \left\{ { - {\rm{i}}\left[{w\left( {t + k\Delta t} \right)} \right]} \right\} $代入式(4)可得

$ \begin{array}{l} {\left[ {{r^{ - 1}}\sin \left( {0.5\omega \Delta t} \right)} \right]^2} = {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_x}h} \right]} } \right\}^2} + \\ {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_y}h} \right]} } \right\}^2} + \\ {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_z}h} \right]} } \right\}^2}. \end{array} $ (5)

其中

$ \omega = kv,k = \sqrt {k_x^2 + k_y^2 + k_z^2} ,{k_x} = k\cos \theta \cos \varphi , $
$ {k_y} = k\cos \theta \sin \varphi ,{k_z} = k\sin \theta . $

式中, ω为角频率; θ为平面波在X-O-Y平面投影后中与x轴的夹角; φ为平面波的方位角。Liu和Sen[24]通过固定θ=0和φ=π/8, 并在波数域进行泰勒展开得到特定角度优化的交错网格差分系数。另一种思路是, 采用函数拟合优化的方式, 通过调整差分系数使等式(5)右边项在给定的波数范围内和所有传播方向上尽可能逼近左边项:

$ \eta \left( {\theta ,\varphi } \right) = {r^2}q - {\sin ^2}\left( {0.5\omega \Delta t} \right). $ (6)

其中

$ \begin{array}{l} q = {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_x}h} \right]} } \right\}^2} + \\ {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_y}h} \right]} } \right\}^2} + \\ {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_z}h} \right]} } \right\}^2}. \end{array} $

Wang等[24]以时空频散差值误差η$ {\ell_2} $范数建立目标函数, 并引入正则化技术稳定地求解该非线性问题, 得到所有传播角度优化的差分系数。考虑到差分系数优化的最终目的是减少相速度频散, 以减少频散差值误差为目标函数得到的优化解并不能代表相速度误差最小。以相速度误差为目标函数是最直接的选择。首先给出相速度频散误差测量公式[23]:

$ \delta \left( {\theta ,\varphi } \right) = \frac{{{v_{{\rm{FD}}}}}}{v} - 1 = \frac{2}{{kv\Delta t}}\arcsin \left( {r\sqrt q } \right) - 1. $ (7)

式中, vFD为数值计算的相速度; δ(θ, φ)越逼近0, 代表频散误差越小。以减少相速度误差δ$ {\ell_2} $范数为目标函数并采用最优化算法进行求解, 这样可以更直接的达到优化目的。但需要说明的是该问题具有很强的非线性, 同时由于反三角函数的存在, 在求解过程中可能会出现不稳定。为了便于求解, 下面考虑用相对频散误差函数(时间频散与空间频散的比值)逼近式(7):

$ \varepsilon \left( {\theta ,\varphi } \right) = \frac{{{r^2}q}}{{{{\sin }^2}\left( {0.5\omega \Delta t} \right)}} - 1. $ (8)

若式(8)与式(7)的函数形态较为接近, 那么可以用以减少ε$ {\ell_2} $范数为目标函数。图 1为一维情况下使用相同差分系数(由空间域泰勒展开获得)、不同差分阶数下对应数值曲线。图 1(a)为式(6)在一维情况下的差值误差函数曲线, 图 1(b)为式(7)所示的相速度误差函数曲线, 图 1(c)为式(8)所表示比值误差函数曲线。通过对比可以看到, 差值误差函数与真实的频散误差在形态上具有较大的差异, 而通过比值的方式可以非常好地逼近频散误差。本文中采用如下目标函数进行差分系数优化:

$ F\left( a \right) = \frac{1}{2}\sum\limits_{k = 0}^{{k_{\rm{e}}}} {\sum\limits_{\theta = 0}^{{\rm{ \mathsf{ π} }}/4} {\sum\limits_{\varphi = 0}^{{\rm{ \mathsf{ π} }}/4} {\left\| {\frac{{{r^2}q\left( a \right)}}{{{{\sin }^2}\left( {0.5\omega \Delta t} \right)}} - 1} \right\|_2^2} } } . $ (9)
图 1 一维介质中不同差分阶数下3种误差函数对比 Fig.1 Comparison of three error function in 1D case with offerent difference order

式中,ke为有效地震波波数。考虑到传播方向的对称性, θφ的优化范围设定为0到π/4。ke的取值决定了优化波数的范围, 当ke取值较大时, 较广的波数范围得到优化, 但会降低低波数段的精度[21]; 值得注意的是当ke取值较小时, 在求解过程中, 可能会出现不稳定[24]。在地震波传播模拟中, 地震子波通常是一个带宽信号[25], 子波的最大频率与地震波的最大波数有关系$ {k_{\max }} = \frac{{2{\rm{ \mathsf{ π} }}{f_{\max }}}}{v} $。在地震波模拟中, 通常使用雷克子波作为震源函数, 给定子波主频f0, 考虑能量分布, 最高频率通常可设为fmax=2.5f0。由于正演模拟中, 模拟精度受到差分阶数的影响, 当阶数较小时可适当减低最高频率优先确保低波数段的精度。直接令$ {k_{\rm{e}}} = \frac{{2{\rm{ \mathsf{ π} }}{f_{\max }}}}{v} $, 当最高频率较小时求解中容易出现不稳定。本文中给出一种改进的优化区间确定方案:

$ {k_{\rm{e}}} = \max \left\{ {{k_{\max }},{k_0}} \right\}. $ (10)

其中, k0由传统差分系数情况下相速度误差δ(k0) > -0.005决定, 如果想优化更广的波数范围, 可以将δ0设为更小的值, 但这可能降低低波数段的精度。分析图 1(b)可以看到不同阶数下得到的k0通常是不同的。当差分阶数较大时, k0取值较大; 当差分阶数较小时, k0取值相对较小, 这也符合常理。通过式(10)确定优化的波数区间可有效避免区间上限较小的情况, 从而确保了解的稳定性。此外, 在二维、三维介质中, 各传播方向的频散通常不同的, 可以分别求取不同方向对应的k0

1.3 牛顿法优化差分系数

Liu[18]在研究规则网格中心差分优化时提出通过频散比值误差逼近相速度误差, 其面对的是一个线性优化问题, 可以直接使用最小二乘函数拟合法得到一个满意的结果。本文中等效交错网格优化面对的是一个非线性问题, 无法直接使用最小二乘函数拟合法。考虑到牛顿法收敛速度快的特点, 采用牛顿法求解该非线性问题[26](式(9)), 使用牛顿法涉及到目标函数的一阶导数(梯度)以及二阶导数(Hessian矩阵)。首先给出式(9)关于差分系数am的一阶导数

$ \frac{{\partial F\left( a \right)}}{{\partial {a_m}}} = \sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{{\rm{ \mathsf{ π} }}/4} {\sum\limits_{\varphi = 0}^{{\rm{ \mathsf{ π} }}/4} {\left[ {f\left( a \right)\frac{{\partial f\left( a \right)}}{{\partial {a_m}}}} \right]} } } . $ (11)

其中

$ f\left( a \right) = \frac{{{r^2}q\left( a \right)}}{{{{\sin }^2}\left( {0.5\omega \Delta t} \right)}} - 1, $
$ \frac{{\partial f\left( a \right)}}{{\partial {a_m}}} = \frac{{{{\left( {v\Delta t} \right)}^2}}}{{{{\sin }^2}\left( {0.5\omega \Delta t} \right)}}\frac{{\partial q\left( a \right)}}{{\partial {a_m}}}, $
$ \begin{array}{l} \frac{{\partial q\left( a \right)}}{{\partial {a_m}}} = \frac{2}{{{h^2}}}\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_x}h} \right]} } \right\}\sin \left[ {\left( {m - } \right.} \right.\\ \left. {\left. {0.5} \right){k_x}h} \right] + \frac{2}{{{h^2}}}\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_y}h} \right]} } \right\} \times \\ \sin \left[ {\left( {m - 0.5} \right){k_y}h} \right] + \frac{2}{{{h^2}}}\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - } \right.} \right.} } \right.\\ \left. {\left. {\left. {0.5} \right){k_z}h} \right]} \right\}\sin \left[ {\left( {m - 0.5} \right){k_z}h} \right]. \end{array} $

目标函数的梯度可表示为$ \mathit{\boldsymbol{g}} = {\left[{\frac{{\partial F\left( a \right)}}{{\partial {a_1}}} \cdots \frac{{\partial F\left( a \right)}}{{\partial {a_m}}} \cdots \frac{{\partial F\left( a \right)}}{{\partial {a_M}}}} \right]^{\rm{T}}} $。使用牛顿法还需计算Hessian矩阵, 目标函数对应的Hessian矩阵的维数为M×M, 元素Hmn可表示为

$ \begin{array}{l} \frac{{{\partial ^2}F\left( a \right)}}{{\partial {a_m}\partial {a_n}}} = \sum\limits_{k = 0}^{{k_{\rm{e}}}} {\sum\limits_{\theta = 0}^{{\rm{ \mathsf{ π} }}/4} {\sum\limits_{\varphi = 0}^{{\rm{ \mathsf{ π} }}/4} {\left[ {f\left( a \right)\frac{{{\partial ^2}f\left( a \right)}}{{\partial {a_m}\partial {a_n}}} + \frac{{\partial f\left( a \right)}}{{\partial {a_m}}} \times } \right.} } } \\ \left. {\frac{{\partial f\left( a \right)}}{{\partial {a_n}}}} \right]. \end{array} $ (12)

其中

$ \frac{{{\partial ^2}f\left( a \right)}}{{\partial {a_m}\partial {a_n}}} = \frac{{{{\left( {v\Delta t} \right)}^2}}}{{{{\sin }^2}\left( {0.5\omega \Delta t} \right)}}\frac{{{\partial ^2}q\left( a \right)}}{{\partial {a_m}\partial {a_n}}}, $
$ \begin{array}{l} \frac{{{\partial ^2}q\left( a \right)}}{{\partial {a_m}\partial {a_n}}} = \frac{2}{{{h^2}}}\sin \left[ {\left( {m - 0.5} \right){k_x}h} \right]\sin \left[ {\left( {n - 0.5} \right){k_x}h} \right] + \\ \frac{2}{{{h^2}}}\sin \left[ {\left( {m - 0.5} \right){k_y}h} \right]\sin \left[ {\left( {n - 0.5} \right){k_y}h} \right] + \\ \frac{2}{{{h^2}}}\sin \left[ {\left( {m - 0.5} \right){k_z}h} \right]\sin \left[ {\left( {n - 0.5} \right){k_z}h} \right]. \end{array} $

在得到目标函数的梯度g与Hessian矩阵H后, 可以使用如下的迭代公式进行求解差分系数:

$ {\mathit{\boldsymbol{a}}^{k + 1}} = {\mathit{\boldsymbol{a}}^k} - {\left( {{\mathit{\boldsymbol{H}}^k}} \right)^{ - 1}}{\mathit{\boldsymbol{g}}^k}. $ (13)

式中, k为迭代次数。设定当误差达到初始误差1%作为迭代终止条件, 由于牛顿法具有很快的收敛速度, 通常只需几次迭代就可得到满意的结果。需要说明的是, 若直接以子波的频带作为优化的波数区间, 当子波的最高频率较小, 求H的逆时容易出现不稳定。通过改进的优化波数选择方案, 在系数计算中可以有效地避开不稳定问题。

2 数值频散分析 2.1 一维频散

一维声波方程等效交错网格数值模拟频散可定义为

$ \delta = \frac{{{v_{{\rm{FD}}}}}}{v}\frac{2}{{kv\Delta t}}\arcsin \left( {r\sqrt q } \right). $ (14)

其中

$ q = {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right)kh} \right]} } \right\}^2}. $

式中, v为声介质传播速度; vFD为地震波数值模拟传播速度。δ大于1代表数值模拟的相速度大于真实速度, 波场快照中出现波前面超前现象(时间频散); δ小于1代表数值模拟的相速度小于真实速度, 波场出现“拖尾”现象(空间频散)。

图 2(a)为低速介质中传统方法(红色线)与本文方法(蓝色线)获取系数的不同差分阶数情况下频散对比。计算参数:空间采样10 m, 时间步长1 ms, 介质速度为1 500 m/s。从图中可以看到, 低速介质中高波数段容易出现空间频散; 随着空间差分阶数的增加, 两种系数在高波段的精度都有所提高; 相比传统系数, 通过本文优化方法得到的系数可在更大的波数范围保持高精度。

图 2 两种方法不同差分阶数下的一维数值模拟频散分析 Fig.2 1D dispersion analysis by two methods with different difference order

图 2(b)为两种系数在高速介质中频散曲线对比。计算参数:介质速度为4 500 m/s, 空间采样为10 m, 时间步长为1 ms。从图中可以看到, 传统方法(红色线)在高速介质中容易出现时间频散; 由于传统方法求取空间差分系数时未考虑时间频散, 增加空间差分阶数并不能有效解决时间频散。本文方法(蓝色线)将时间频散与空间频散统一起来, 并通过优化算法获取差分系数, 可有效地将时间频散考虑到空间差分求取中, 因而可以有效减少时间频散。

2.2 二维频散

二维声波方程数值等效交错网格数值频散可定义为

$ \delta \left( \theta \right) = \frac{{{v_{{\rm{FD}}}}}}{v} = \frac{2}{{kv\Delta t}}\arcsin \left( {r\sqrt q } \right). $ (15)

其中

$ \begin{array}{l} q = {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_x}h} \right]} } \right\}^2} + \\ {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_y}h} \right]} } \right\}^2}. \end{array} $

二维数值频散不仅与波数有关, 还与传播角度有关, 不同传播方向其频散是不同的。图 3(a)为差分阶数M=5在低速介质v=1 500 m/s中两种方法的频散图, 可以看到:传统方法在0°角方向(x轴或者z轴)有严重的空间频散, 在45°角方向上空间频散相对较小; 本文方法在获取系数时考虑了不同传播方向的影响, 因而可以实现对不同方向的优化。图 3(b)为高速介质v=4 500 m/s中两种方法的频散图, 其中差分阶数M=5。在高速介质中, 传统方法在中波数段出现了时间频散, 而在高波数段不同方向出现的频散类型是不同的。通过本文方法可对多个传播的时间频散与空间频散进行有效压制。

图 3 两种方法不同差分阶数下的二维数值模拟频散分析 Fig.3 2D dispersion analysis by two methods with different difference order
2.3 三维频散

三维声波方程数值等效交错网格数值频散可定义为

$ \delta \left( {\theta ,\varphi } \right) = \frac{{{v_{{\rm{FD}}}}}}{v} = \frac{2}{{kv\Delta t}}\arcsin \left( {r\sqrt q } \right). $ (16)

其中

$ \begin{array}{l} q = {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_x}h} \right]} } \right\}^2} + \\ {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_y}h} \right]} } \right\}^2} + \\ {\left\{ {\sum\limits_{m = 1}^M {{a_m}\sin \left[ {\left( {m - 0.5} \right){k_z}h} \right]} } \right\}^2}. \end{array} $

图 4为传统方法与本文方法在三维声波数值模拟的不同传播方向的频散图。其中, 介质速度为2 500 m/s, 时间步长为1 ms, 空间步长为10 m, 差分阶数M=8。对比图 4(a)图 4(b), 可以看到通过本文方法可以有效地将数值频散误差限定在一个较小的范围中。由于本文中以时空频散关系为基础, 因此本文方法既有效地减少了时间频散, 又压制了空间频散。采用本文优化的系数, 三维正演时各个传播方向的精度可得到有效地提高。

图 4 两种方法不同差分阶数下的三维数值模拟频散分析 Fig.4 3D dispersion analysis by two methods with different difference order
3 模型试算结果

为了验证本文方法的有效性, 对3个典型模型进行了正演模型测试。首先, 使用两个均匀模型分别测试本文方法在压制空间频散与时间频散的情况; 最后, 使用Marmousi模型验证本文方法对复杂速度场的适应性。

3.1 均匀模型

频散分析表明等效交错网格正演方法在低速介质中容易出现空间频散。为了验证本文方法对空间频散的改善, 首先使用等效交错网格模拟地震波在速度为1 500 m/s的均匀介质中的传播。横向与纵向空间采样均为10 m, 时间步长为1 ms, 纵横向网格点数为201。考虑空间采样定理, 采用主频为24 Hz的雷克子波作为震源, 优化系数求取时考虑的最高频率为55 Hz。图 5为3种参数对应的1.2 s波场快照。其中, 图 5(a)采用传统系数M=4的模拟结果, 图 5(b)为传统系数M为8的模拟结果, 图 5(c)使用本文优化后系数M=4的模拟结果。在图 5(a)出现了高波数成分延迟的现象(空间频散), 通过增加差分阶数(M=8)空间频散得到有效压制(图 5(b)), 对比图 5(b)图 5(c)可以看到采用本文方法在不增加差分阶数下可达到传统系数更高阶的正演效果。尽管增加差分阶数可有效压制频散, 但会带来计算量的增加, 采用本文优化方法可在不增计算量的情况下达到更高阶差分的精度。

图 5 1.2 s时刻波场快照 Fig.5 Snapshots at 1.2 s

分析表明在高速介质中高波数分量容易出现相速度大于真实速度的情况, 为了检测本文方法对时间频散的改善, 模拟地震波在速度为4 000 m/s的均匀介质中的传播, 横向与纵向空间采样均为10 m, 纵横向网格点数为501。为了便于从波场快照观察时间频散, 采用主频为80 Hz的雷克子波作为震源, 优化系数求取时考虑的最高频率为165 Hz。数值模拟中, 差分阶数M=5。图 6(a)为传统方法时间步长1 ms对应的1.0 s时刻的波场快照, 图 6(b)为传统方法时间步长dt=0.5 ms对应的相同时刻的波场快照, 图 6(c)为本文优化方法时间步长1 ms在1.0 s时刻的波场快照。在图 6(a)中出现了较为严重的时间频散(尤其是对角线方向), 而在横向与纵向上存在轻微的空间频散; 通过减少时间步长, 时间频散得到很好的压制, 但空间频散并未得到很好解决(图 6(b)); 由于本文方法以减少相对时空频散为目标函数, 因而采用本文优化方法可有效地减缓对角线上的时间频散, 同时对横向与纵向上的频散也有一定的改善。

图 6 1.0 s时刻波场快照 Fig.6 Snapshots at 1.0 s
3.2 Marmousi模型

进一步测试本文方法对复杂模型的适应性, 为此取SEG提供的国际Marmousi模型的一部分(图 7(a))进行正演测试处理。考虑到人工边界反射, 使用吸收边界法[27]压制边界反射。模型网格大小为1 001×751, 空间采样间隔为Δxz=15 m, 模型速度为1 500~5 500 m/s, 密度如图 7(b)所示。本文方法在正演模拟前须求取不同速度下的差分系数, 由于速度相差不大时, 计算出的差分系数相差很小, 对于复杂模型本文采取间隔5 m/s求取系数, 总共计算了801组差分系数。计算参数为:差分阶数M=4, 采用主频为13 Hz的雷克子波作为震源, 计算时间步长为Δt=1 ms。正演模拟总共进行了5 001次迭代, 检波器位于地表, 地震记录时长3 s(图 8)。为了方便与传统方法对比, 抽取出图 7中的一部分进行放大(图 9(a))。图 9(b)为传统方法M=4的正演地震记录局部放大图, 此外, 图 9(c)为传统方法M=6的正演地震记录局部放大图。在图 9(b)中白色箭头所指处存在明显的空间频散, 而在图 9(a)图 9(b)对应位置处频散得到了有效的压制。正演测试表明本文方法适用于复杂模型并且以较低差分阶数达到传统高阶差分的模拟精度, 然而传统方法以计算量换取计算精度, 本文方法在几乎不增加计算量的情况可有效地提高计算精度。因此本文方法可应用于逆时偏移、全波形反演以及最小二乘逆时偏移中提高其计算效率。

图 7 Marmousi速度与密度模型 Fig.7 Marmousi P-wave velocity model and density model
图 8 本文方法模拟的地震记录 Fig.8 Seismogram for Marmousi model by presented method
图 9 地震记录局部放大 Fig.9 Zoom of seismograms
4 结束语

本文中提出了一种时空域等效交错网格差分系数确定方法, 使用牛顿法减少所有传播方向上相对时空频散关系误差获取优化的差分系数, 为了稳定求取有效的差分系数, 在决定优化展开的波数区间同时参考子波的频带范围与传统系数下误差函数的形态。通过频散分析发现本文方法可有效地压制时间频散与空间频散。此外, 数值模拟也表明在空间频散方面, 本文方法能以较低的差分阶数获得传统高阶差分正演的精度; 而针对时间频散, 本文方法可使用较大的时间步长达到传统方法小时间步长的精度, 可有效地提高计算效率。此外, 本文方法可有效地应用于复杂模型正演模拟, 并易于推广到非对称网格逆时偏移成像、全波形反演及最小二乘偏移中。

参考文献
[1]
BOORE D M. Finite difference methods for seismic wave propagation in heterogeneous materials[J]. Methods in Computational Physics, 1972, 11: 1-37.
[2]
VIRIEUX J. P-SV wave propagation in heterogeneous me-dia: velocity stress finite difference method[J]. Geophysics, 1986, 51(4): 889-901. DOI:10.1190/1.1442147
[3]
LEVANDER A R. Fourth-order finite-difference P-SV seismograms[J]. Geophysics, 1984, 53(11): 1425-1436.
[4]
方刚, FOMELSergey, 杜启振. 交错网格Lowrank有限差分及其在逆时偏移中的应用[J]. 中国石油大学学报(自然科学版), 2014, 38(2): 44-51.
FANG Gang, FOMEL Sergey, DU Qizhen. Lowrank finite difference on a staggered grid and its application on reverse time migration[J]. Journal of China University of Petroleum (Edition of Natural Science), 2014, 38(2): 44-51.
[5]
李振春, 雍鹏, 黄建平, 等. 基于矢量波场分离弹性波逆时偏移成像[J]. 中国石油大学学报(自然科学版), 2016, 40(1): 42-48.
LI Zhenchun, YONG Peng, HUANG Jianping, et al. Elastic wave reverse time migration based on vector wavefield separation[J]. Journal of China University of Petroleum (Edition of Natural Science), 2016, 40(1): 42-48.
[6]
黄建平, 曹晓莉, 李振春, 等. 最小二乘逆时偏移在近地表高精度成像中的应用[J]. 石油地球物理勘探, 2014, 49(1): 107-112.
HUANG Jianping, CAO Xiaoli, LI Zhenchun, et al. Least square reverse time migration in high resolution imaging of near surface[J]. Oil Geophysical Prospecting (in Chinese), 2014, 49(1): 107-112.
[7]
黄建平, 李振春, 孔雪, 等. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究[J]. 地球物理学报, 2013, 56(5): 1716-1725.
HUANG Jianping, LI Zhenchun, KONG Xue, et al. A study of least-squares migration imaging method for fractured-type carbonate reservoir[J]. Chinese Journal of Geophysics (in Chinese), 2013, 56(5): 1716-1725. DOI:10.6038/cjg20130529
[8]
PRATT R G, SHIN C, HICKS G J. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophysical Journal International, 1998, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x
[9]
VIRIEUX J, OPERTO S. An overview of full-waveform inversion in exploration geophysics[J]. Geophysics, 2009, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367
[10]
董良国, 马在田, 曹景忠, 等. 一阶弹性波方程交错网格高阶差分解法[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 (in Chinese), 2000, 43(3): 411-419.
[11]
LIU Y, SEN M K. A practical implicit finite-difference method: examples from seismic modelling[J]. Journal of Geophysics and Engineering, 2009, 6(3): 231-249. DOI:10.1088/1742-2132/6/3/003
[12]
吴国忱, 王华忠. 波场模拟中的数值频散分析与校正策略[J]. 地球物理学进展, 2005, 20(1): 58-65.
WU Guochen, WANG Huazhong. Analysis of numerical dispersion in wave-field simulation[J]. Progress in Geophysics (in Chinese), 2005, 20(1): 58-65.
[13]
HOLBERG O. Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena[J]. Geophysical Prospecting, 1987, 35(6): 629-655. DOI:10.1111/gpr.1987.35.issue-6
[14]
CHU C, STOFFA P L. Determination of finite-difference weights using scaled binomial windows[J]. Geophysics, 2012, 77(3): W17-W26. DOI:10.1190/geo2011-0336.1
[15]
ZHANG J, YAO Z. Optimized finite-difference operator for broadband seismic wave modeling[J]. Geophysics, 2013, 78(1): A13-A18. DOI:10.1190/geo2012-0277.1
[16]
ZHANG J, YAO Z X. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm[J]. Journal of Computational Physics, 2013, 250: 511-526. DOI:10.1016/j.jcp.2013.04.029
[17]
ETGEN J T. A tutorial on optimizing time domain finite-difference schemes: "Beyond Holberg"[J]. Stanford Exploration Project Report, 2007, 129: 33-43.
[18]
LIU Y. Globally optimal finite-difference schemes based on least squares[J]. Geophysics, 2013, 78(4): T113-T132. DOI:10.1190/geo2012-0480.1
[19]
REN Z, LIU Y. Acoustic and elastic modeling by optimal time-space-domain staggered-grid finite-difference schemes[J]. Geophysics, 2015, 80(1): T17-T40. DOI:10.1190/geo2014-0269.1
[20]
LIU H, DAI N, NIU F, et al. An explicit time evolution method for acoustic wave propagation[J]. Geophysics, 2014, 79(3): T117-T124. DOI:10.1190/geo2013-0073.1
[21]
TAN S, HUANG L. A staggered-grid finite-difference scheme optimized in the time-space domain for modeling scalar-wave propagation in geophysical problems[J]. Journal of Computational Physics, 2014, 276: 613-634. DOI:10.1016/j.jcp.2014.07.044
[22]
BARTOLO L 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
[23]
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
[24]
WANG Y, LIANG W, NASHED Z, et al. Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method[J]. Geophysics, 2014, 79(5): T277-T285. DOI:10.1190/geo2014-0078.1
[25]
WU Z, ALKHALIFAH T. The optimized expansion based low-rank method for wavefield extrapolation[J]. Geophysics, 2014, 79(2): T51-T60. DOI:10.1190/geo2013-0174.1
[26]
YONG P, HUANG J, LI Z, et al. Optimized equivalent staggered-grid FD method for elastic wave modelling based on plane wave solutions[J]. Geophysical Journal International, 2017, 208(2): 1157-1172. DOI:10.1093/gji/ggw447
[27]
CERJAN C, KOSLOFF D, KOSLOFF R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J]. Geophysics, 1985, 50(4): 705-708. DOI:10.1190/1.1441945