地球物理学进展  2015, Vol. 30 Issue (2): 878-888   PDF    
2D声波频率域数值模拟中几种有限差分方法的对比分析
马晓娜1,2, 李志远1, 谷丙洛1,2, 柯沛3, 梁光河1     
1. 中国科学院地质与地球物理研究所, 北京 100029;
2. 中国科学院大学, 北京 100049;
3. 中国石油集团东方地球物理勘探有限责任公司新兴物探开发处, 涿州 072751
摘要:频率域数值模拟是频率域全波形反演的基础, 在地震波场数值模拟中占有重要地位.相对于时间域数值模拟, 频率域数值模拟具有两个明显的优势:没有时间累计误差, 适合于并行计算.然而, 严重的数值频散和巨大的内存损耗是阻碍其应用的两大瓶颈.为解决这两个问题, 基于有限差分方法, 学者提出了多种差分格式, 如优化9点、15点、17点以及25点差分格式.本文从频散关系、计算效率和存储量三个方面, 对比、分析了以上四种差分方法.基于2D声波方程, 通过在均匀模型、层状模型以及Marmousi模型上的应用效果, 对每种方法的优缺点进行了总结, 为高精度数值模拟和声波频率域全波形反演提供方法选择上的参考.
关键词频率域数值模拟     有限差分     优化差分格式    
Comparisons and analysis of several optimization finite-differencing schemes in 2D acoustic frequency-domain numerical modeling
MA Xiao-na1,2, LI Zhi-yuan1, GU Bing-luo1,2, KE pei3, LIANG Guang-he1     
1. Institute of Geology and Geophysics, the Chinese Academy of Sciences, Beijing 100029, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. BGP Inc., China National Petroleum Corporation, Zhuozhou 072751, China
Abstract: Frequency-domain numerical modeling is the basis of full waveform inversion, which is very important in seismic wave modeling. Compared with conventional time domain forward modeling, frequency-space domain forward modeling has two advantages: no time accumulated error; suitable for parallel cmputation. However, there are two problems with this modeling method: severe grid dispersion and huge memory consumption which have restricted its application to a large extent. Finite-differencing method is effective to achieve the simulation in the frequency-domain. To overcome the disadvantages, based on the finite-differencing, scholars proposed optimal 9-point, 15-point, 17-point and 25-point scheme. Based on the frequency domain 2D acoustic equation, these methods are applied in the homogeneous model, the single layer model and Marmousi model, then the results are compared and analyzed. By the comparisons of above methods from the point of dispersion, efficiency and computer memory of numerical simulation, we summarize the advantages and disadvantages of the methods, which gives the selection and reference to the following high precision numerical simulation and frequency domain full waveform inversion.
Key words: frequency-domain numerical modeling     finite-differencing     optimal differencing scheme    
0 引 言

全波形反演是近些年一个热门的研究方向,该方法可以直接使用地震波场的运动学和动力学信息重建地下速度结构(Virieux et al., 2009刘国峰等,2012).全波形反演既可以在时间域实现,也可在频率域实现.Tarantola(1984)最早通过炮点正向传播波场和接收点逆时传播的残差波场互相关生成梯度方向,利用梯度寻优更新模型参数,实现了二维时间域全波形反演.Pratt和Worthington(1990)率先使用频率-空间域有限差分模拟实现全波形反演.频率域全波形反演相对于时间域而言,具有计算效率高、数据选择灵活等优点,特别是在多炮情况下,不同炮可以使用相同的LU分解因子(曹书红和陈景波,2012).

正演模拟是全波形反演的基础,而有效的求解方式可以保证反演结果的计算精度和效率(刘有山等,2014).国内外许多学者对频率域正演有过深入的研究.Lysmer和Drake(1972)最初提出使用有限元方法实现频率域正演模拟.之后,Marfurt(1984)Marfurt和Shin(1989)使用有限差分方法实现频率域数值模拟,并指出频率域数值模拟不存在稳定性问题,而且适合大量炮点的同时模拟.Pratt(1990 ab)、Pratt和Worthington(1990)Shin和Sohn(1998)等学者对此做了进一步的研究.Alford等(1974)对二维声波标量方程有限差分模拟结果和解析解做了分析,认为有限差分方法可以实现精确的数值模拟.目前,有限差分方法被广泛应用于地震波数值模拟(Gu et al., 2012).

对于二维声波方程和弹性波方程,Pratt和Worthington等(1990)最初使用传统的二阶差分格式实现数值模拟,每个波长需要使用13个网格点以保证相速度误差不超过1%,波阻抗矩阵存储需要巨大的计算机内存.Luo和Schuster(1990)提出省时交错网格方法,该方法后期被证明与Pratt的方法一致.由于频散严重和存储大等问题,最初的这些频率域数值模拟方法没有得到广泛的应用.随后,很多学者对频率域数值模拟做了改进.基于旋转坐标系,Jo等(1996)对声波方程使用优化9点差分格式来近似拉普拉斯算子和质量加速项,在减小计算量和存储量的同时,频散现象得到明显改善.随后,在Jo等(1996)的研究基础上,Shin和Sohn(1989)将坐标系进行三次旋转,利用0°、26.6°、45°和63.4°四个方向的差分项,得出频率域声波方程的25点差分格式,虽然计算精度明显提高,但是系数矩阵带宽的增加导致了计算效率降低.基于Luo等(1990)的省时交错网格方法和Jo等(1996)的旋转坐标系方法,Hustedt 等(2004)提出了9点混合优化算子.

近几年,针对频率域数值模拟计算精度和效率等问题,国内很多学者做了相应的研究.曹书红和陈景波(2012)在传统四阶精度9点格式和旋转坐标系的基础上,构造了四阶精度17点格式,在克服频散方面具有一定的优越性.刘璐等(2013)提出优化15点差分格式,在不明显增加计算量的前提下,较好的压制了频散.Liu(2014)在传统五点方法的基础上,利用泰勒展开和共轭梯度法计算有限差分系数,数值模拟精度提高,保留了传统五点差分格式计算效率高的优点.针对计算效率问题,刘红伟等(2010)使用CPU/GPU协同并行处理对全波形反演加速.廖建平等(2011)采用稀疏矩阵中的嵌套剖分法,减少了计算机内存的需求,加快了计算速度.魏哲枫等(2014)采用基于非规则三角网格的地震波正演方法,应用了最优化剖分技术来离散模型空间,取得较好的效果.传统的频率域正演算法只能适用于相同的纵横向空间采样间隔(Chen,2012张衡等,2014).对此,Chen(2012)提出了平均导数方法(ADM),此方法将声波频率域方程中空间导数项的差分近似表示为正交方向上3个网格点的加权平均,适用于不同的横纵向空间采样间隔.

从国内外学者的研究可以发现,频率-空间域地震波正演具有很多优点(Pratt和Worthington,1990Pratt,1990b卞爱飞等,2010Gu et al., 2012):

(1)使用隐式差分算子,不存在时间累计误差;

(2)多震源进行模拟时,使用相同的LU分解结果,计算效率很高;

(3)各单频分量相互独立,非常适合于并行计算;

(4)模拟衰减因素的影响可以通过直接计算衰减因子关于频率的函数,更灵活;

(5)波动方程反演中,频带选择灵活.

虽然频率-空间域有限差分数值模拟方法有很多优点,但是它本身也有一些缺点.如计算精度与效率之间的矛盾,就相同的计算精度而言,与时间域数值模拟相比,其每个波长需要的网格点数更多(Jo et al., 1996).在对线性方程组求解中,波阻抗矩阵的存储需要很大的内存,这些都是频率-空间域方法发展需要面对的难题,在三维地震勘探中波阻抗矩阵求解问题更加突出(任浩然等,2009).

本文从频散关系、计 算速度和存储量三个方面,对2D声波频率域数值模拟的优化9点(Jo et al., 1996)、25点(Shin和Sohn,1989)、17点(曹书红和陈景波,2012)、15点(刘璐等,2013)差分格式做出了对比分析,并将四种差分格式应用到简单模型和复杂模型的正演中,并总结了优缺点,为后续的频率域数值模拟和全波形反演的研究提供方法基 础.

1 频率域数值模拟方法原理

频率域数值模拟可以表示为

其中,L(w)是波阻抗矩阵,P(w)为频率域波场,s(w)为频率域震源.波阻抗矩阵 L(w)是大型高度稀疏、非对称的矩阵,其元素依赖于频率和传播介质的性质.通常情况下,波阻抗矩阵不适合直接求逆,而是使用直接法和迭代法求解方程.直接法是对系数矩阵直接做LU分解,系数矩阵分解为上、下三角矩阵,分解结果回代求解方程,可以多炮共用分解结果,在二维情况下计算效率较高(Pratt,1990a).本文采用直接法求解方程(1),使用有限差分方法实现数值模 拟.

1.1 声波优化9点差分格式

笛卡尔坐标系下,频率-空间域二维标量声波波动方程为

其中,P(x,z,ω)是频率域波场值,ω是角频率,v(x,z)是传播速度,f(ω)是频率域震源项,δ(x)是狄利克雷(Dirichlet)函数,xs和zs是震源的水平和垂直位置.

优化9点差分格式的构造主要包括三个方面(图 1a):45°旋转坐标系,平均质量加速项和优化系数(Jo et al., 1996).首先对坐标系旋转45°,拉普拉斯算子用0°坐标系和45°坐标系表示为

图 1 四种有限差分格式网格点分布图
(a)9点有限差分格式网格点;(b)15点有限差分格式网格点;(c)17点有限差分格式网格点;(d)25点有限差分格式网格点.
Fig. 1 The figure of four finite difference computational grid
(a)9-point finite difference computational grid;(b)15-point finite difference computational grid;(c)17-point finite difference computational grid;(d)25-point finite difference computational grid.

对方程(3)使用二阶有限差分,我们可以得到

其中,Δ=Δx=Δz为网格间距,Pm,n代表 [xm,zn]=[x0+(m-1)Δx,z0+(n-1)Δz]处的波场值.

根据Marfurt(1984)的研究,质量加速项 P可以表示为差分格式中9点的加权平均,并依据各点到中心点的距离划分,可以得到:

其中,b+4c+4d=1.

方程(4)、(5)带入(2),即可得到差分方程.

这种多方向的数值差分形式被证实关于传播角度有相对较小的数值各向异性(Jo et al., 1996).

由于,计算机存储容量和计算速度等运行环境的限制,在地震波场数值模拟时,只能在有限尺寸的模型上求解,因此需要对人工截断边界产生的虚假反射进行特殊处理.本文使用最佳匹配层吸收边界条件(PML).

1.2 声波优化25点差分格式

基于旋转坐标系,将拉普拉斯算子表示为0°、26.6°、45°和63.4°四个方向坐标系的组合,25点有限差分格式可以表示为(如图 1d):

其中,Δ=Δx=Δz为网格间距,=1,∇2P(0,Δ),∇2P(0,2Δ),∇2P(45,2 Δ),∇2P(45,2 2 Δ),∇2P(26.6,5 Δ),∇2(63.4,5 Δ)Pm,n的计算方式与上述优化9点差分格式的相应算子类似,此处不再赘述.

质量加速项 P 可以表示为

其中,

此差分格式使用四个不同的旋转方向来计算加权拉普拉斯算子和质量加速项,各项异性减小,频散情况得到明显改善(Shin和Sohn,1989).

1.3 声波优化四阶精度17点差分格式

为了满足高精度地震成像的需要,曹书红和陈景波(2012)结合经典四阶精度9点格式和优化9点格式,构造了新的17点差分格式,是声波优化9点差分格式的高阶表达,如图 1c所示.

将拉普拉斯算子表示成0°坐标系和45°>坐标系下四阶差分格式的组合,并将质量加速项表示成相应点的结合.则方程(2)的四阶精度17点差分方程可写为

其中,Δ=Δx =Δz,Δ 2P (0,Δ),Δ2P(45,2 Δ)采用四阶精度差分格式,计算方式与上文类似.质量加速项 P的计算类似于优化25点差分格式,只不过相应的计算点数减少.

1.4 声波优化15点差分格式

刘璐等(2013)提出优化15点差分格式,在垂直方向上采用3个层的网格点来近似偏导数算子,而在水平方向本文采用5个层的网格点来近似偏导数算子(图 1b).根据加权平均算子和平均加速度项,构建了由12个参数约束的15点差分格式,水平方向空间差分算子表示为

其中,Δ=Δx=Δz,b1和b2通过优化各项系数使频散关系式

与1最接近求取,在下文中详细介绍.同理,可以得到垂直方向的空间差分算子.

平均质量加速项表示为差分格式中15点的加权平均.因此,将空间差分算子和平均质量加速项带入方程(2)即可得到声波优化15点差分格式.

2 优化系数

优化系数的计算是频率域数值模拟中重要部分,优化系数的选取直接关系到数值模拟的频散情况.频率域数值频散关系与传播角度θ、每个波长内的网格点数G以及各项系数有关.基本思路是优化各项系数使频散关系式与1最接近,也就是二者的残差最小.

Jo等(1996)提出使用最速下降法处理频散关系式,得到优化系数.Shin和Sohn(1998)提出一种新的方法,首先假定归一化相速度对于所有θ和G等于1,推导出以优化系数作为未知数的超定方程,然后使用奇异值分解对超定方程求解,得到优化系数.刘璐等(2013)使用Gauss-Newton法得到优化系数.Liu(2014)使用共轭梯度法得到优化系数.我们给定角度θ取0~ ,间隔选为5° 的取值范围为 ~,间隔选为0.001.

本文以声波优化9点差分方法为例,使用最速下降法求取优化系数,目标函数为

其中,D为各个优化系数,k为迭代次数,J是E 的一阶导数.给定优化系数的初始值,一次迭代后的新参数表示为

D=D+ΔD,反复迭代直至满足精度要求即可.

3 数值频散、计算速度及存储量对比

3.1 频散分析

数值频散压制是频率域波动方程模拟的重要研究内容(卞爱飞等,2010).本文对四种方法分别进行了数值频散分析.求取频散关系的思路基本相同,本文以声波优化9点差分格式为例介绍.

平面波方程的解为

其中,θ为平面波传播角.

把(15)式代入到(6)式中,并且令G= ,用来表示每个波长内的网格数(Jo et al., 1996),求出 .根据相速度的定义,,即得出此差分格式所对应的频散关系为

图 2 四种差分格式频散曲线
(a)9点差分格式频散曲线;(b)15点差分格式频散曲线;(c)17点差分格式频散曲线;(d)25点差分格式频散曲线.
Fig. 2 Dispersion curves of the four different scheme
(a)The dispersion curve of 9-point difference scheme;(b)The dispersion curve of 15-point difference scheme;(c)The dispersion curve of 17-point difference scheme;(d)The dispersion curve of 25-point difference scheme.

将上面求得的最优化系数带入到频散关系中,计算优化9点(a)、15点(b)、17点(c)以及25点(d)差分格式的频散曲线(图 2).通过对比以上四种差分格式的频散曲线,可以看出,优化9点差分格式的频散现象最严重,以相速度误差1%为标准,最小波长内的网格点数为4.相对于优化9点差分格式,优化25点差分格式频散的频散情况明显改善,最小波长内的网格点数为2.8个;优化15点差分格式也有效的控制了频散情况,最小波长内的网格点数为3个,已非常接近优化25点差分格式的模拟精度;优化17点差分格式的频散情况也明显改善,最小波长内的网格点数仅为2.6个(表 1).

表 1 波阻抗矩阵带宽和频散误差为1%的单位波长内的网格点数(G) Table 1 The b and with of impedance matrix the value of G when the dispersion error is 1%
3.2 计算速度和存储量

频率域正演的效率在很大程度上取决于波阻抗矩阵LU分解的快慢(Husted等,2004),而LU分解的效率则由波阻抗矩阵带宽决定,所以使用较小的带宽来获取较高的精度是频率域数值模拟要解决的关键问题(刘璐等,2013).此外,频率域数值模拟需要消耗巨大的计算机内存,而消耗内存最大的是波阻抗矩阵的存储.在相同计算精度条件下,影响波阻抗矩阵存储量的关键因素有两个,一是最小波长内所需的网格点数,二是矩阵的带宽(Gu et al., 2012).

如果传统二阶有限差分方法需要Nx×Nz个网格点来满足精度要求,那么9点差分需要 个网格点,15点差分格式需要个网格点,17点差分格式需要个网格点,25点差分格式需要个网格点.可以看出,对于相同的精度要求,17点差分格式所需要的网格点数最少,15点差分格式所需要的网格点数接近于25点(表 1).

下面我们对比以上四种差分格式的带宽,以8×8的网格为例显示对应的波阻抗矩阵带宽,如图 3所示.9点差分格式非零元素带宽是2×nx+3,上带宽和下带宽均为nx+2;15点差分格式非零元素带宽是2×nx+5,上带宽和下带宽均为nx+3;17点差分格式带宽为4×nx+3,上带宽和下带宽均为2×nx+2;25点差分格式带宽为4×nx+5,上带宽和下带宽均为2×nx+3.可以看出,15点差分格式相对于9点差分格式,只增加了两条带宽;17点和25点差分格式系数矩阵带宽相对9点和15点差分格式,带宽增加了近两倍(表 1).带宽增加,相应的计算效率随之下降,需要的存储空间变大.

图 3 系数矩阵结构
(a)优化9点差分格式系数矩阵;(b)优化15点差分格式系数矩阵;(c)优化17点差分格式系数矩阵;(d)优化5点差分格式2系数矩阵.
Fig. 3 The structure of the impedance matrix
(a)Impedance matrix of 9-point difference scheme;(b)Impedance matrix of 15-point difference scheme;(c)Impedance matrix of 17-point difference scheme;(d)Impedance matrix of 25-point difference scheme.
4 模型计算 4.1 均匀模型

本文采用均匀模型的单道记录来测试不同差分格式数值频散效果和计算精度.模型网格大小为100×100,空间采样间隔为10 m,时间采样间隔为1 ms,记录时间为0.5 s,震源子波使用雷克子波,主频率为30 Hz,速度为3000 m/s,震源位置为(500 m,20 m),检波器置于z=20 m的界面上.均匀模型测时采用优化9点、15点、17点、25点差分格式进行频率域数值模拟并与解析方法对比,解析解计算使Alford等(1974)提出的公式为

其中,F(ω)为频率域震源,H(2)0为零阶二类Hankel函数,r为观测点到震源点的距离.

取四种差分方法的单道记录并与解析方法进行对比,如图 4.根据Alford等(1974)描述的频散具有延迟的、变宽的、带有振荡尾巴的特征,可以看出9点差分格式的地震记录出现了较严重的频散,与解析解相比出现较大的误差.17点差分格式与15点差分格式的单道地震记录中出现了较轻微的频散现象.25点差分格式可以实现精确模拟,没有出现明显的数值频散误差,与解析解对应良好.

图 4 均匀模型优化差分格式9点、15点、17点、25点与解析法单道记录对比 Fig. 4 Comparison of single traces obtained by optimal difference scheme 9-point,15-point,17-point,25-point and analytical method on homogeneous model
4.2 单层介质模型

模型网格大小为100×100,空间采样间隔为6 m,时间采样间隔为2 ms,记录时间为1.4 s,震源子波主频率为25 Hz,第一层纵波速度为1000 m/s,第二层纵波速度为2000 m/s,震源位置在(300 m,18 m),检波器置于z=18 m的界面上(如图 5).对此模型分别进行优化9点、15点、17点、25点差分格式的数值模拟,并分别对正演中LU分解耗时进行记录.四种差分格式的地震记录如图 6所示,从所标示的区域可以看出,9点差分格式的正演结果频散明显,15点差分格式的频散现象有所消除,17点和25差分格式正演结果都没有观测到频散现象.

图 5 单层速度模型 Fig. 5 The velocity model of single layer

图 6 单层模型正演记录
(a)优化9点差分格式求得的地震记录;(b)优化15点差分格式求得的地震记录;(c)优化17点差分格式求得的地震记录;

(d)优化25点差分格式求得的地震记录.
Fig. 6 Forward modeling result of single layer model
(a)Wavefield acquired by 9-point difference scheme;(b)Wavefield acquired by 15-point difference scheme;(c)Wavefield acquired by 17-point difference scheme;(d)Wavefield acquired by 25-point difference scheme.

表 2看出,9点差分格式LU分解耗时为449.21 s;15点差分格式LU分解耗时为4733.73 s;17点差分格式LU分解耗时为15434.09 s;25点差分格式LU分解耗时为15996.11 s.9点差分格式计算耗时最小,25点耗时最大.15点差分格式系数矩阵带宽比9点差分格式多两条,LU分解耗时相差十倍;17点与25点差分格式LU分解耗时相差500 s左右,相差不大.17点和25点差分格式与9点差分格式相比,矩阵带宽相差一倍,LU分解耗时相差近34倍.

4.3 复杂模型

复杂模型为Marmous模型的一部分,网格大小为150×300,空间采样间隔为7 m,时间采样间隔为1 ms,记录时间为0.7 s,震源子波主频率为25 Hz,最大速度为4450 m/s,最小速度为2282 m/s,平均速度为2904 m/s,如图 7所示.震源位置在(900 m,30 m),检波器置于z=30 m的界面上.从地震记录可以看出,四种优化差分方法均适用于复杂模型的数值模拟计算,优化9点与15点差分格式的出现轻微的频散情况,优化17点和25点差分格式的模拟结果没有明显的频散情况出现(图 8).

图 7 Marmousi速度模型 Fig. 7 The velocity model of Marmousi model

图 8 Marmousi模型正演记录
(a)优化9点差分格式求得的地震记录;(b)优化15点差分格式求得的地震记录;(c)优化17点差分格式求得的地震记录;(d)优化25点差分格式求得的地震记录.
Fig. 8 Forward modeling result of Marmousi model
(a)Wavefield acquired by 9-point difference scheme;(b)Wavefield acquired by 15-point difference scheme;(c)Wavefield acquired by 17-point difference scheme;(d)Wavefield acquired by 25-point difference scheme.

对于复杂模型,优化15点差分格式的计算时间约为9点差分格式的46倍,优化17点差分格式的计算时间约为9点差分格式的179倍,约为15点差分格式的3.8倍,约为25点差分格式的18倍(表 2).

表 2 LU分解耗时对比 Table 2 The comparison of LU decomposition time

表 2记录的简单模型和复杂模型的LU分解耗时可以发现,系数矩阵带宽的数量直接影响LU分解的耗时,进而影响整个频率域数值模拟的计算效率.从简单模型到复杂模型的数值模拟可以看出,15点相对于9点差分格式的LU分解耗时比翻了5倍;同时,17点相对于15点差分格的LU分解耗时比翻了5倍.此外,随着模型复杂程度的增加,17点差分格式的耗时超过25点差分格式. 4.4 全波形反演算例

本小结将四种正演差分格式分别应用于全波形反演中,并对反演效果进行对比分析.这里选用图 9a的凹陷模型作为初始模型,网格大小为60×60,纵横网格间距为10 m,共59炮,炮点位于z=20 m的界面上,炮间间隔为10 m,检波器位于地表,每炮60道,震源采用雷克子波,主频为20 Hz.图 9b是用10×10的网格光滑之后的初始速度模型.反演频率以3 Hz为步长从1 Hz选取到50 Hz,共17个反演频率,每个频率的最大迭代次数为30次.本文使用BFGS算法进行全波形反演,目标函数使用残差的L2范数(董良国等,2013).除正演所采用的差分格式不同之外所有参数完全相同.

图 10中可以看出,使用四种正演差分格式进行反演基本都能得到模型的形态,浅部反演的模型值与理论值吻合较好.正演中使用优化15点、17点、25点差分格式的模型反演结果略好于优化9点差分格式,特别是在模型下部两层交界处;就前三种差分格式使用的反演耗时来看,优化15点差分格式计算速度最快(表 3).此外,优化25点差分格式在模型凹陷处的反演结果有些波动.

图 9 凹陷模型和初始模型
(a)凹陷模型;(b)初始模型.
Fig. 9 The hollow model and initial model
(a)Hollow model;(b)Initial model.
图 10 四种正演差分格式的反演效果对比
(a)优化9点差分格式;(b)优化15点差分格式;(c)优化17点差分格式;(d)优化25点差分格式.
Fig. 10 The comparison of inverse result from four forward modeling difference scheme
(a)9-point difference scheme;(b)15-point difference scheme;(c)17-point difference scheme;(d)25-point difference scheme.

表 3 全波形反演耗时对比 Table 3 The comparison of full waveform inversion times
5 结 论

综上分析,以上四种差分格式中,优化9点和15点差分格式的计算耗时小,而且后者的计算精度高;优化17点和25点差分格式的模拟精度高,但是前者在复杂模型的计算中耗时较其他三种差分格式都大.

通过频散曲线对比分析,可以发现,优化9点差分格式频散情况明显,其他三种差分格式的频散误差得到改善.同时,相速度误差1%的情况下,四种差分格式最小波长内的网格点数都减小到4以下.对于计算机耗时和存储情况,通过对比波阻抗矩阵带宽可以发现,25点和17点差格式的系数矩阵带宽约为9点和15点差分格式的两倍,相应的计算效率低,需要的存储空间大.通过从简单模型到复杂模型的数值模拟实验,可以发现25点和17点差分格式LU分解耗时远大于9点和15点差分格式.而且,随着模型复杂程度的增加,LU分解耗时比呈数倍的增加.其中,对于部分Marmousi模型正演,17点差分格式的耗时最大.

优化15点差分格式兼顾了系数矩阵带宽和精度,改善精度的同时尽可能地减小了矩阵带宽,降低了存储空间,提高了计算效率. 优化17点差分格式实现了四阶精度差分,改变了高精度正演采用“加密网格,低阶实现”(增加存储量、降低计算效率)的现状,但是随着模型复杂程度的增加,计算效率降低.

最后,本文将四种正演差分格式分别应用于全波形反演中,发现使用优化15点差分格式进行正演,得到的反演结果好,而且计算效率高,具有一定的优越性.未来的工作中将重点研究如何改进频率域正演,包括精度和计算效率等,为高精度、高效率的全波形反演打下基础.

参考文献
[1] Alford R M, Kelly K R, Boore D M. 1974. Accuracy of finite-difference modeling of the acoustic wave equation [J]. Geophysics, 39(6): 834-842, doi:10.1190/1.1440470.
[2] Berenger J-P. 1994. A perfectly matched layer for the absorption of electromagnetic waves [J]. Journal of Computational Physics, 114(2): 185-200, doi:10.1006/jcph.1994.1159.
[3] Bian A F, Yu W H, Zhou H W. Progress in the frequency-domain full waveform inversion method[J]. Progress in Geophys. (in Chinese), 25(3): 982-993, doi:10.3969/j.issn.1004-2903.2010.03.037.
[4] Cao S H, Chen J B. 2012. A 17-point scheme and its numerical implementation for high-accuracy modeling of frequency-domain acoustic equation [J]. Chinese J. Geophys. (in Chinese), 55(10): 3440-3449, doi:10.6038/j.issn.0001-5733.2012.10.27.
[5] Chen J B. 2012. An average-derivative optimal scheme for frequency-domain scalar wave equation [J]. Geophysics, 77(6): T201-T210, doi: 10.1190/Geo2011-0389.1.
[6] Chen J B. 2013. A generalized optimal 9-point scheme for frequency-domain scalar wave equation [J]. Journal of Applied Geophysics, 92: 1-7, doi: 10.1016/j.jappgeo.2013.02.008.
[7] Chen J B. 2014. A 27-point scheme for a 3D frequency-domain scalar wave equation based on an average-derivative method [J]. Geophysical Prospecting, 62(2): 258-277, doi: 10.1111/1365-2478.12090.
[8] Clapp R G. 2009. Reverse time migration with random boundaries [C].// 79th Annual International Meeting, SEG, Expanded Abstracts, 2809-2813.
[9] Ding J C, Chang X, Liu Y K, et al. 2007. Layer by layer waveform inversion of seismic reflection data[J]. Chinese J. Geophys. (in Chinese), 50(2): 574-580, doi:10.3321/j.issn:0001-5733.2007.02.031.
[10] Dong L G, Chi B X, Tao J X, et al. 2013. Objective function behavior in acoustic full-waveform inversion[J]. Chinese J. Geophys. (in Chinese), 56(10): 3445-3460, doi: 10.6038/cjg20131020.
[11] GAO Feng-Xia, LIU Cai, FENG Xuan, et al. 2013. Comparisons and analyses of several optimization methods in the application of frequency-domain full waveform inversion [J]. Progress in Geophys. (in Chinese), 28(4): 2060-2068, doi: 10.6038/pg20130450.
[12] Gu B L, Liang G H, Li Z Y. 2012. A 21-point finite difference scheme for 2D frequency-domain elastic wave modelling [J]. Exploration Geophysics, 44(3): 156-166.
[13] Husted B, Operto S, Virieux J. 2004. Mixed-grid and staggered grid finite-difference methods for frequency-domain acoustic wave modelling[J]. Geophysical Journal International, 157(3): 1269-1296, doi: 10.1111/j.1365-264X.2004.02289.x.
[14] Jo C-H, Shin C S, Suh J H. 1996. An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator [J]. Geophysics, 61(2): 529-537, doi: 10.1190/1.1443979.
[15] Liao J P, Liu H X, Wang H Z, et al. 2011. Study on rapid highly accurate acoustic wave numerical simulation in frequency space domain [J]. Progress in Geophysics (in Chinese), 26(4): 1359-1363, doi: 10.3969/j.issn.1004-2903.2011.04.029.
[16] Liu G F, Liu H, Meng X H, et al. 2012. Frequency-related factors analysis in frequency domain waveform inversion [J]. Chinese J. Geophys. (in Chinese), 55(4): 1345-1353, doi: 10.6038/j.issn.0001-5733.2012.04.030.
[17] Liu H W, Li B, Liu H, et al. 2010. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation [J]. Chinese J. Geophys. (in Chinese), 53(7): 1725-1733, doi: 10.3969/j.issn.0001-5733.2010.07.024.
[18] Liu L, Liu H, Liu H W. 2013. Optimal 15-point finite difference forward modeling in frequency-space domain [J]. Chinese J. Geophys. (in Chinese), 56(2): 644-652, doi: 10.6038/ cjg20130228.
[19] Liu Y. 2014. An optimal 5-point scheme for frequency-domain scalar wave equation [J]. Journal of Applied Geophysics, 108(3): 19-24, doi: 10.1016/j.jappgeo.2014.06.006.
[20] LIU You-Shan, TENG Ji-Wen, Xu Tao, et al. 2014. Numerical modeling of seismic wavefield with the SEM based on Triangles [J]. Progress in Geophysics (in Chinese), 29(4): 1715-1726, doi: 10.6038/pg20140430.
[21] Lysmer J, Drake L A. 1972. A finite-element method for seismology[A].// Bolt B A ed. Methods in Computational Physics, Volume 11: Seismology: Surface Waves and Earth Oscillations[M]. New York: Academic Press Inc., 181-216.
[22] Mallick S, Frazer L N. 1987. Practical aspects of reflectivity modeling [J]. Geophysics, 52(10): 1355-1364, doi: 10.1190/1.1442248.
[23] Marfurt K, Shin C. 1989. The future of iterative modelling in geophysical exploration[A]. //Eisner E ed. Handbook of Geophysical Exploration: I-Seismic Exploration, 21-Supercomputers in Seismic Exploration [M]. Pergamon Press, 203-228.
[24] Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations [J]. Geophysics, 49(5): 533-549, doi: 10.1190/1.1441689.
[25] Operto S, Virieux J, Amestoy P, et al. 2007. 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study [J]. Geophysics, 72(5): SM195-SM211, doi: 10.1190/1.2759835.
[26] Pratt R G. 1990a. Frequency-domain elastic wave modeling by finite differences: A tool for cross-hole seismic imaging [J]. Geophysics, 55(5): 626-632.
[27] Pratt R G. 1990b. Inverse theory applied to multi-source cross-hole tomography, Part 2: Elastic wave-equation method [J]. Geophysical Prospecting, 38(3): 311-329.
[28] Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography, Part 1: Acoustic wave equation method [J]. Geophysical Prospecting, 38(3): 287-310.
[29] Pratt R G, Shin C, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency- space seismic waveform inversion [J]. Geophysics, 133(2): 341-362.
[30] Ren H R, Wang H Z, Gong T. 2009. Seismic modeling of scalar seismic wave propagation with finite-difference scheme in frequency-space domain[J]. Geophysical Prospecting for Petroleum (in Chinese), 48(1): 20-26.
[31] Ren H R, Huang G H, Wang H Z, et al. 2013. A research on the Hessian operator in seismic inversion imaging[J]. Chinese J. Geophys. (in Chinese), 56(7): 2429-2436, doi: 10.6038/cjg20130728.
[32] Shi Y M, Zhang Y, Yao F C, et al. 2014. Methodology of seismic imaging for hydrocarbon reservoirs based on acoustic full waveform inversion [J]. Chinese J. Geophys., 57(5): 1599-1611, doi: 10.6038/cjg20140523.
[33] Shin C, Sohn H. 1998. A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator[J]. Geophysics, 63(1): 289-296. doi: 10.1190/1.1444323.
[34] tekl I, Pratt R. 1998. Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators [J]. Geophysics, 63(5): 1779-1794, doi: 10.1190/1.1444472.
[35] Symes W M. 2007. Reverse time migration with optimal chekpointing[J]. Geophysics, 72(5): SM213-SM-221, doi: 10.1190/1.2742686.
[36] Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation [J]. Geophysics, 49(8): 1259-1266.
[37] Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics [J]. Geophysics, 74(6): WCC1-WCC26, doi: 10.1190/1.3238367.
[38] Wei Z F, Gao H W, Zhang J F. 2014. Time-domain full waveform inversion based on an irregular-grid acoustic modeling method[J]. Chinese J. Geophys. (in Chinese), 57(2): 586-594, doi: 10.6038/cjg20140222.
[39] Wu G C, Luo C M, Liang K. 2007. Frequency-space domain finite difference numerical simulation of elastic wave in TTI media [J]. Journal of Jilin University (Earth Science Edition) (in Chinese), 37(5): 1023-1033, doi: 10.3969/j.issn.1671-5888.2007.05.030.
[40] Zhang H, Liu H, Liu L, et al. 2014. Frequency domain acoustic equation high-order modeling based on an average-derivative method [J]. Chinese J. Geophys. (in Chinese), 57(5): 1599-1611, doi: 10.6038/cjg20140523.
[41] 卞爱飞, 於文辉, 周华伟. 2010. 频率域全波形反演方法研究进展. 地球物理学进展[J]. 25(3): 982-993, doi: 10.3969/j.issn.1004-2903.2010.03.037.
[42] 曹书红, 陈景波. 2012. 声波方程频率域高精度正演的17点格式及数值实现[J]. 地球物理学报, 55(10): 3440-3449, doi: 10.6038/j.issn.0001-5733.2012.10.027.
[43] 丁继才, 常旭, 刘伊克,等. 2007. 反射地震数据的逐层波形反演[J]. 地球物理学报, 50(2): 574-580, doi: 10.3321/j.issn:0001-5733.2007.02.031.
[44] 董良国, 迟本鑫, 陶纪霞,等. 2013. 声波全波形反演目标函数性态[J]. 地球物理学报, 56(10): 3445-3460, doi: 10.6038/cjg20131020.
[45] 高凤霞, 刘财, 冯晅等. 2013. 几种优化方法在频率域全波形反演中的应用效果及对比分析研究[J]. 地球物理学进展, 28(4): 2060-2068, doi: 10.6038/pg20130450.
[46] 廖建平, 刘和秀, 王华忠,等. 2011. 快速高精度的频率空间域声波数值模拟方法研究[J]. 地球物理学进展, 26(4): 1359-1363, doi: 10.3969/j.issn.1004-2903.2011.04.029.
[47] 刘国峰, 刘洪, 孟小红,等. 2012. 频率域波形反演中与频率相关的影响因素分析[J]. 地球物理学报, 2012, 55(4): 1345-1353, doi: 10.6038/j.issn.0001-5733.2012.04.030.
[48] 刘红伟, 李博, 刘洪,等. 2010. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 53(7): 1725-1733, doi: 10.3969/j.issn.0001-5733.2010.07.024.
[49] 刘璐, 刘洪, 刘红伟. 2013. 优化15点频率-空间域有限差分正演模拟[J]. 地球物理学报, 56(2): 644-652, doi: 10.6038/cjg20130228.
[50] 刘有山, 滕吉文, 徐涛,等. 2014. 三角网格谱元法地震波场数值模拟[J]. 地球物理学进展, 29(4): 1715-1726, doi: 10.6038/pg20140430.
[51] 任浩然, 王华忠, 龚婷. 2009. 标量地震波频率-空间域有限差分法数值模拟[J]. 石油物探, 48(1): 20-26.
[52] 任浩然, 黄光辉, 王华忠,等. 2013. 地震反演成像中的Hessian 算子研究[J]. 地球物理学报, 56(7): 2429-2436, doi: 10.6038/cjg20130728.
[53] 石玉梅, 张研, 姚逢昌,等. 2014. 基于声学全波形反演的油气藏地震成像方法 [J]. 地球物理学报, 57(2): 607-617, doi: 10.6038/cjg20140224.
[54] 魏哲枫, 高红伟, 张剑锋. 2014. 基于非规则网格声波正演的时间域全波形反演[J]. 地球物理学报, 57(2): 586-594, doi: 10.6038/cjg20140222.
[55] 吴国忱, 罗彩明, 梁楷. 2007. TTI 介质弹性波频率-空间域有限差分数值模拟[J]. 吉林大学学报(地球科学版), 37(5): 1023-1033.
[56] 张衡, 刘洪, 刘璐,等. 2014. 基于平均导数方法的声波方程频率域高阶正演[J]. 地球物理学报, 57(5): 1599-1611, doi: 10.6038/cjg20140523.