黏声波方程波场模拟的通用格式自适应系数频域有限差分方法

徐文豪,巴晶*,J.M.Carcione,朱鹤松

1 河海大学地球科学与工程学院,南京 211100 2 National Institute of Oceanography and Applied Geophysics-OGS,Borgo Grotta Gigante 42/c,Trieste,Sgonico 34010,Italy

0 引言

由于地下介质的黏弹性,地震波在地下传播时存在着衰减效应.黏声波方程可以有效表征地震波传播的衰减特征,因此基于黏声波方程的波场数值模拟在地震资料处理、反演研究中具有重要应用,如基于黏声波方程的全波形反演(Malinowski et al.,2011;李海山等,2018;Chen et al.,2021)和逆时偏移(Yang and Zhu,2018;周彤等,2018;Qu et al.,2020;谷丙洛等,2023)等.黏声波模拟的主要方法包括时域有限差分方法(赵建国等,2014;Guo et al.,2016)、频域有限差分(FDFD)方法(张立彬和王华忠,2010;Aghamiry et al.,2022)、有限元法(Zyserman and Gauzellino,2004;Groby and Tsogka,2006)、伪谱法(Zhu et al.,2014;吴玉等,2017)、近似解析离散法(宋国杰等,2019)、傅里叶有限差分法(张金海等,2007,2008)等.在这些方法中,频域有限差分方法具有可灵活处理衰减效应、高效模拟多炮或窄带地震资料、避免时间频散、并行效率高等优点,因而在黏声波模拟中具有较多的应用案例(Arntsen et al.,1998;李添才等,2014;Keating and Innanen,2019).用FDFD方法模拟黏声波方程的主要挑战在于其需要求解大型的稀疏线性方程组,而大型稀疏线性方程组的求解往往需要花费较多的计算量和内存.克服这一挑战的一个常用方法是发展所需一个波长内的网格点数(number of points per wavelength,NPPW)更少的FDFD方法,从而减少波场模拟所需的网格点数,进而减少相应大型稀疏线性方程组的阶数.

作为使用FDFD方法求解黏声波方程的特例以及发展黏声波FDFD方法的基础,关于声波方程的FDFD方法,前人已经进行了较多的研究.针对二维声波方程,Pratt和Worthington(1990)提出了经典的五点FDFD方法,在相速度误差不超过1%的精度要求下,该方法所需的NPPW为13.Liu(2014)提出了一个优化的五点FDFD方法,其精度明显优于传统的五点方法.Jo等(1996)提出一个适用于正方形网格的优化的九点FDFD方法,当x方向和z方向的空间采样间隔相等的时候,该方法能将所需的NPPW减少到4.通过引入平均导数方法(average-derivative method,ADM),Chen(2012)将Jo等(1996)的FDFD方法推广到长方形网格情况下.刘璐等(2013)提出一个优化的十五点FDFD方法,将所需的NPPW减少到2.97,同时相应的稀疏系数矩阵具有可接受的带宽.Shin和Sohn(1998)提出一个优化二十五点FDFD方法,其所需的NPPW为2.5.基于ADM方法,张衡等(2014)将Shin和Sohn(1998)的方法推广到了长方形网格情况.Fan等(2017)提出一个优化通用二十五点FDFD方法,该方法将所需的NPPW减少到2.13,而且很多常用的二维FDFD格式都可以被视作该方法的特殊情况.Xu和Gao(2018)提出一个自适应系数的九点FDFD方法,其在正方形网格和非正方形网格情况下所需的NPPW分别为2.12和2.44,同时所需的计算量和内存需求与优化九点FDFD方法相当.针对三维声波方程,在相速度误差不超过1%的精度要求下,传统的二阶七点FDFD方法所需的NPPW为13(Chen,2014).Operto等(2007)发展了一个适用于立方体网格的优化二十七点方法,将所需的NPPW减少到了4.Chen(2014)发展了一个优化ADM二十七点方法,该方法对一般的长方体网格所需的NPPW为4.Xu等(2021a)发展了自适应系数的三维二十七点FDFD方法,并给出了一种更容易使用和推广的自适应系数FDFD通用格式.自适应系数的二十七点FDFD方法在立方体网格和非立方体网格下所需的NPPW分别为2.2和2.7.

为求解一般的黏声波方程,传统的二阶五点FDFD方法所需NPPW较多(Amini and Javaherian,2011).高阶FDFD方法可以有效减少所需的NPPW(Xu et al.,2021b),但其对应的大型稀疏线性方程组含有较大的带宽和较多的非零元素,因此所需的计算量和内存较大(Takekawa and Mikada,2016).Arntsen等(1998)将Holberg(1987)采用的十字形优化有限差分系数应用于黏声波方程的FDFD求解,并针对均匀模型应用解析解验证了其有效性,十字形优化FDFD相对高阶FDFD可以减少所需空间差分的长度,但相对紧凑格式的FDFD方法所需计算量和内存仍然较大,同时十字形优化FDFD在低波数区域的精度相对高阶FDFD会有所减少(Zhang and Yao,2013).通过将声波方程FDFD中的实波数替换为黏声波方程的复波数,Amini和Javaherian(2011)将Jo等(1996)关于声波方程的优化九点FDFD直接应用于黏声波方程,并给出了在一般黏声波模型中的测试,然而Jo等(1996)的FDFD只适用于正方形网格,同时Amini和Javaherian(2011)并没有通过黏声波方程的解析解或参考解验证相应FDFD方法在一般黏声波模型中的有效性.Chen(2012,2014)在提出声波方程ADM优化FDFD方法的同时,给出了在黏声波方程中的相应推广,但并未给出一般黏声波模型的数值算例.在将自适应系数FDFD应用到声波全波形反演的同时,Aghamiry等(2022)给出了自适应系数FDFD方法在黏声波方程的推广,但同样未给出一般黏声波模型的数值算例,同时Aghamiry等(2022)所采用的基于旋转坐标系的自适应系数FDFD方法只适用于空间采样间隔相等的情形(Chen,2012).

为了减少不同空间采样间隔之比下的黏声波波场模拟所需的NPPW,本文针对黏声波波场模拟发展了一种通用格式自适应系数FDFD方法.同时,为验证自适应系数FDFD对一般黏声波模型的有效性,本文针对三个典型的黏声波模型,分别使用解析解和高阶FDFD参考解验证所提方法的有效性.在本方法中,所需的自适应系数FDFD格式可通过在传统的二阶FDFD格式的基础上加入相关的系数随NPPW和空间采样间隔之比变化的校正项得到,其中校正项按网格点与中心点的距离进行分类选取.本文首先给出了基于黏声波方程的通用格式自适应系数FDFD方法,然后分别针对均匀黏声波模型、层状黏声波模型、Marmousi黏声波模型,将所提FDFD与二阶五点FDFD、ADM优化九点FDFD的波场模拟结果进行对比,并分别通过解析解和基于高阶FDFD的参考解验证所提方法的有效性.

1 通用格式自适应系数FDFD方法

1.1 自适应系数FDFD的通用格式

频率域二维无源黏声波方程可表示为

(1)

其中u表示频率域波场,ω表示角频率,vc表示复速度.关于复速度vc,本文采用如下常用的Kolsky-Futterman模型进行描述(Kolsky,1956; Futterman,1962):

(2)

其中vr表示在参考角频率ωr下的参考波传播速度,Q表示品质因子,sgn表示符号函数,i表示虚数单位.据Toverud和Ursin (2005)以及Aghamiry等(2020),本文将参考角频率对应的频率值设为50 Hz.需指出的是,公式(2)对应的一维平面波解的形式是u0e-i(ω/vc)x,如果对应的一维平面波解的形式是u0ei(ω/vc)x,则式(2)虚部的符号将改变.

(3a)

r=Δx/Δz,

(3b)

(4)

其中um,n=u(mΔx,nΔz),(vc)m,n=vc(mΔx,nΔz),ηi(i=1,2,3)表示通用格式中的自适应FDFD系数.

1.2 自适应FDFD系数的确定

据Operto等(2007)、Chen (2012,2014)、Aghamiry等(2022),黏声波FDFD的优化或自适应FDFD系数可通过先假设品质因子趋于无穷大,再极小化相应平面波解的代入误差得到.参考已有的关于黏声波方程的FDFD方法(Chen,2012; Aghamiry et al.,2022),本文假设黏声波FDFD的误差主要来源于数值频散,因此在确定自适应FDFD系数的过程中,只考虑品质因子趋于无穷大的情况(即Q→∞).此时黏声波方程的平面波解与声波方程的平面波解具有相同形式,可表示为

u(x,z)=u0e-ikr[sin(θ)x+cos(θ)z],

(5)

(6a)

(6b)

(6c)

(6d)

-1],

(6e)

其中η=[η1,η2,η3].极小化平面波解的代入误差,所需的自适应FDFD系数可通过求解如下优化问题得到:

(7)

式(7)中的优化问题可通过Polak-Ribière-Polyak共轭梯度法(Polyak,1969)快速求解,

(8a)

ηk+1=ηk+αkdk,

(8c)

1.3 FDFD系数查找表

表1 基于声波数值频散关系的黏声波方程自适应系数九点FDFD方法的部分系数查找表(r=3)

需指出的是,对于通用格式自适应系数九点FDFD方法,相应的查找表只是一个3×500的实数双精度矩阵,因此几乎不会增加波场模拟所需的内存需求.此外,对网格尺寸为nx×nz的波场模拟,构造自适应系数FDFD系数矩阵相对构造常系数FDFD系数矩阵只需额外增加3nxnz次线性插值,其计算量相对计算FDFD系数矩阵的各分量以及求解FDFD线性方程组而言很小.

1.4 数值频散分析

使用FDFD方法数值求解波动方程时容易导致数值频散,即FDFD方法对应的相速度与模型真实速度有偏差.这里给出通用格式自适应系数九点FDFD方法的数值频散分析结果.令vph表示自适应FDFD方法对应的数值相速度.则自适应系数FDFD方法的平面波解(和波场方程的平面波解不同)可表示为

(9)

为得到通用格式自适应系数九点FDFD的数值相速度vph和模型真实速度vr的偏差,将自适应系数FDFD的平面波解代入自适应九点FDFD格式,即式(4),得到

(10a)

(10b)

(10c)

-1],

(10e)

图1和图2分别给出空间采样间隔比率r=1和r=3时自适应系数九点FDFD和二阶五点FDFD、优化ADM九点FDFD的数值频散曲线对比.图1和图2结果表明,相对于二阶五点FDFD和ADM优化九点FDFD,通用格式自适应系数九点FDFD可以更为有效的压制正方形网格和非正方形网格下的数值频散.以数值频散误差不超过1%为标准,相对于二阶五点FDFD和ADM优化九点FDFD,新方法可将黏声波方程一个波长内所需的网格点数从4减少到2.5.

图1 r=1时不同FDFD方法的数值频散曲线对比

图2 r=3时不同FDFD方法的数值频散曲线对比

1.5 吸收边界条件

本文采用完美匹配层(perfectly matched layer,PML)吸收边界条件(Liu and Tao,1997; Ren and Liu,2013;高静怀等,2018)压制人工边界反射.带有PML吸收边界和震源项的二维黏声波方程可表示为

-f(ω)δ(x-xs)δ(z-zs),

(11)

(12)

其中(xm,zn)=(mΔx,nΔz),δm,ms和δn,ns表示克罗内克函数,(ms,ns)表示震源的网格点下标.

2 数值算例

为验证所提方法的有效性,针对三个典型的黏声波模型,将新方法的模拟结果分别与二阶五点FDFD以及ADM优化九点FDFD的模拟结果进行对比.同时,对均匀模型采用解析解作为参考解,对非均匀模型采用阶数足够高的高阶FDFD作为参考解,高阶FDFD采用的具体阶数参考Zhang和Yao(2013)关于使用不同阶数的有限差分近似空间偏导的波数域误差的研究.对于所有算例,FDFD生成的大型稀疏线性方程组均通过PARDISO直接求解器(Petra et al.,2014)求解,计算平台均为超算平台的30个节点(AMD EPYC 7452 @2.35GHz),每个节点有64核和256 GB的内存,不同频率的波场通过消息传递接口(message passing interface,MPI)多节点并行计算.

2.1 均匀模型算例

首先考虑网格点数为201×201的均匀模型.模型在50 Hz参考频率下的速度和品质因子分别为2000 m·s-1和50.模型的采样间隔为Δx=Δz=16.7 m,各侧PML吸收边界层数均为20,震源选为主频20 Hz的雷克子波,震源位置在模型中心.检波器布置在地表所有网格点处.波场模拟的频率范围是0.5～60 Hz且频率采样间隔为0.5 Hz.从频域地震记录转换为时域地震记录的时间范围是0到2 s且时间采样间隔为1 ms.图3给出本算例使用的部分自适应FDFD系数(r=1),表明所使用的自适应FDFD系数是较光滑的.图4给出了0.9 s时二阶五点FDFD、ADM优化九点FDFD、自适应系数九点FDFD的波场快照对比.图5给出了各检波器接收的地震记录对比.图6给出了在(x,z)=(1670 m,0 m)的地表中心处上述三种FDFD接收的地震记录与解析解(Carcione et al.,1988)的对比.图4、图5、图6表明对这个均匀黏声波模型,相对于二阶五点FDFD和ADM优化九点FDFD,本文所提的自适应系数九点FDFD可以更有效地减少黏声波模拟误差.同时,图4、图5、图6也给出了相应速度模型的声波模拟结果,图4和图5体现了黏声波传播相对声波传播主频降低的特征,图6体现了黏声波传播相对声波传播振幅衰减的特征.此外上述三种FDFD在超算30个节点上的计算时间分别为1 s、1 s、1 s,所需内存分别为0.22 GB、0.25 GB、0.25 GB,所需的计算时间和内存相当.

图3 均匀黏声波模型使用的r=1时的自适应FDFD系数

图4 均匀黏声波模型(Q=50)在0.9 s时不同方法得到的波场快照对比

图5 均匀黏声波模型(Q=50)不同方法得到的接收地震记录对比

图6 均匀黏声波模型(Q=50)中声波方程解析解、黏声波方程二阶五点FDFD、黏声波方程ADM优化九点FDFD、黏声波方程通用格式自适应系数九点FDFD波场模拟结果与解析解单道对比

为验证通用格式自适应系数FDFD的误差主要来源于数值频散,我们在上述均匀模型的基础上,进一步考虑品质因子Q值的变化对模拟效果的影响.图7给出了Q值分别取250、50、10的三种情形下,在(x,z)=(1670 m,0 m)地表中心处通用格式自适应系数FDFD接收的地震记录与解析解的对比.图7表明在不同的Q值情况下,通用格式自适应系数FDFD均能与解析解较好符合.

图7 对取不同Q值的均匀黏声波模型,通用格式自适应系数九点FDFD波场模拟结果与解析解单道对比

2.2 层状模型算例

其次考虑一个网格点数为401×401的两层模型.在50 Hz参考频率下,模型在第一层和第二层的速度分别为2000 m·s-1和4000 m·s-1,品质因子分别为50和200.模型的采样间隔为Δx=16.7 m、Δz=8.35 m,沿x和z方向的单侧PML吸收边界层数分别为20和40,震源选为主频20 Hz的雷克子波,震源位置为(xs,zs)=(3340 m,83.5 m).检波器布置在地表的所有网格点处.波场模拟的频率范围是0.25～60 Hz且频率采样间隔为0.25 Hz.从频域地震记录转换为时域地震记录的时间范围是0到4 s且时间采样间隔为1 ms.本算例采用沿x方向72阶、沿z方向12阶的高阶FDFD作为参考解,该高阶FDFD为传统二阶五点FDFD的高阶版本(Xu et al.,2021b).图8给出了本算例使用的部分自适应FDFD系数(r=2).图9给出了1.5 s时上述四种FDFD的波场快照对比.图10给出了各检波器接收的地震记录的对比.图11给出了在(x,z)=(1670 m,0 m)处上述四种FDFD接收的地震记录的对比.图9、图10、图11表明,对这个层状黏声波模型,相对于二阶五点FDFD和ADM优化九点FDFD,自适应系数九点FDFD可更为有效的减少黏声波模拟的误差.同时,图9、图10、图11也给出了相应速度模型的声波模拟结果,图9和图10体现了黏声波传播相对声波传播主频降低的特征,图11体现了黏声波传播相对声波传播振幅衰减的特征.此外上述四种FDFD在超算30个节点上的计算时间分别为3 s、4 s、4 s、216 s,所需内存分别为0.78 GB、0.91 GB、0.91 GB、16.5 GB.通用格式自适应系数九点FDFD的计算时间和所需内存与二阶五点FDFD、ADM优化九点FDFD相当,均远小于参考解所需的计算时间和内存.

图8 两层黏声波模型使用的r=2时的自适应FDFD系数

图9 两层黏声波模型在1.5 s时不同方法得到的波场快照对比

图10 两层黏声波模型不同方法得到的接收地震记录对比

图11 两层黏声波模型中声波方程高阶FDFD(沿x方向72阶、沿z方向12阶)、黏声波方程二阶五点FDFD、黏声波方程ADM优化九点FDFD、黏声波方程通用格式自适应系数九点FDFD波场模拟结果与参考解(黏声波方程沿x方向72阶、沿z方向12阶的高阶FDFD)的单道对比

2.3 Marmousi模型算例

本算例考虑一个网格点数为737×750的Marmousi模型.在50 Hz参考频率下,模型速度如图12a所示,模型的品质因子(图12b)是基于李氏经验公式通过速度模型换算得到,即Q≈3.516×v2.2×10-6(李庆忠,1993).模型的采样间隔为Δx=12.5 m、Δz=4 m,沿x和z方向的单侧PML吸收边界层数分别为20和40,震源选为主频20 Hz的雷克子波,震源位置为(xs,zs)=(4600 m,40 m).检波器布置在地表所有网格点处.波场模拟频率范围是0.125～60 Hz且频率采样间隔为0.125 Hz.从频域地震记录转换为时域地震记录的时间范围是0～8 s且时间采样间隔为1 ms.本算例采用沿x方向72阶、沿z方向8阶的高阶FDFD作为参考解,该高阶FDFD为传统二阶五点FDFD的高阶版本(Xu et al.,2021b).图13给出了本算例使用的部分自适应FDFD系数(r=3.125).图14给出了各检波器接收的地震记录的对比.图15给出了在(x,z)=(3750 m,0 m)处上述四种FDFD接受的地震记录的对比.图14和图15的结果表明,对该Marmousi黏声波模型,相对于二阶五点FDFD和ADM优化九点FDFD,自适应系数九点FDFD可以更为有效的减少模拟误差.同时,图14和图15也给出了相应速度模型的声波模拟结果,图14体现了黏声波传播相对声波传播主频降低的特征,图15体现了黏声波传播相对声波传播振幅衰减的特征.此外上述四种FDFD在超算30个节点的计算时间分别为20 s、23 s、24 s、974 s,所需内存分别为2.4 GB、3.0 GB、3.0 GB、47.5 GB,通用格式自适应系数九点FDFD的计算时间和所需内存与另两种方法相当,远小于参考解所需时间和内存.

图12 Marmousi黏声波模型

图13 Marmousi黏声波模型使用的r=3.125时的自适应FDFD系数

图14 Marmousi黏声波模型不同方法得到的接收地震记录对比

图15 Marmousi黏声波模型中声波方程高阶FDFD(沿x方向72阶、沿z方向8阶)、黏声波方程二阶五点FDFD、黏声波方程ADM优化九点FDFD、黏声波方程通用格式自适应系数九点FDFD波场模拟结果与参考解(黏声波方程沿x方向72阶、沿z方向8阶的高阶FDFD)的单道对比

3 结论

本文针对黏声波方程给出了一种适用于不同空间采样间隔之比的通用格式自适应系数FDFD方法,通过假设黏声波FDFD的误差主要来源于数值频散,所需的自适应FDFD系数根据声波方程的数值频散关系和查找表快速给出,并用于黏声波方程FDFD的数值频散分析和数值模拟中.针对三个典型的黏声波模型,分别采用解析解和基于高阶FDFD的参考解验证了所提出方法的有效性.数值频散分析表明,在空间采样间隔相等或不等的情况下,以相速度误差不超过1%为标准,通用格式自适应系数FDFD方法所需的一个波长内的采样点数均小于2.5.数值实验表明,在不同的空间采样间隔之比的较粗网格情况下,相对于二阶五点FDFD和ADM优化九点FDFD,新方法均可在相似的计算量和内存需求下,有效提高黏声波数值模拟的精度.