起伏地表QR径向基函数有限差分及其在弹性波逆时偏移中的应用

段沛然, 谷丙洛, 李振春, 李青阳

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

2 深层油气重点实验室, 青岛 266580

3 中国石油长庆油田分公司勘探开发研究院, 西安 710018

0 引言

目前,油气地震勘探的重点已由老油田、老区块的简单构造逐渐向山地、丘陵等复杂地表区域转移(刘红伟等,2011;岳玉波等,2012).这些区域勘探面临地表起伏大、近地表横向速度变化剧烈等难题,给地震资料偏移成像带来巨大的挑战(裴正林,2004;黄建平等,2014;杨宏伟等,2022).对于复杂起伏地表地震资料的偏移成像方法主要有两类：一是在偏移前采用高程基准面静校正,将地震资料校正到浮动或者固定基准面上,然后再实施偏移成像(Berryhill,1979;Yilmaz and Lucas,1986;李振春等,2007);二是基于起伏地表直接进行偏移成像(Wiggins,1984;Beasley and Lynn,1992;Margrave and Yao,2000;何英等,2002;Shragge and Sava,2005;程玖兵等,2006;Liu et al.,2011).当基岩出露,地表高程变化大,近地表横向速度变化剧烈时,静校正往往会扭曲波场,降低地震成像的质量.从起伏地表直接进行偏移成像的方法,其算法设计隐含静校正过程,具有自然适应复杂地表条件的能力,一直是是业内研究的热点和目标.

由起伏地表直接偏移的成像方法包括射线类偏移、单程波方程偏移和双程波方程偏移.与射线类(Beylkin,1985;Gray and Bleistein,2009)和单程波方程(Claerbout,1971;Mulder and Plessix,2004)偏移方法相比,基于双程波方程的逆时偏移不受倾角和偏移孔径限制,避免了上下行波分离,可对多种波场(一次反射波、多次波、棱柱波等)准确成像,因此在复杂地下构造成像方面具有巨大的理论优势(Whitmore,1983;Baysal et al.,1983;Gu et al.,2015).波场延拓作为逆时偏移算法的核心环节,其关键在于高精度的地震波数值模拟方法.目前,起伏地表地震波数值模拟方法研究主要分为网格划分和自由地表边界条件两个部分.其中,网格划分可分为规则网格、曲变网格、非结构网格以及无网格.基于规则网格的有限差分方法(Finite Difference Method,FDM)(Kelly et al.,1976;Marfurt,1984;Virieux,1986)具有计算效率高、易于实现、适用于任意复杂介质等优点,是目前应用最为广泛的地震波数值模拟方法(Ren et al.,2017;吴国忱等,2018; 吴建鲁和吴国忱,2018; 段沛然等,2020).为了解决FDM中的起伏地表问题,Jih等(1988)利用具有不同角度的线段拟合地表起伏,采用单边近似法处理不规则自由表面条件实现起伏地表地震波数值模拟.该方法需采用精细网格剖分实现线段贴合起伏地表,因此计算效率较低.Hayashi等(2001)在FDM中引入局部变网格剖分模拟起伏地表自由表面的瑞雷波.变网格由于依赖于网格形变程度,仅能减少杂散衍射,无法完全消除矩形网格阶梯近似所引起的衍射(孙成禹等,2008; 李振春等,2014).Tessmer等(1992)将变网格思想引入交错网格FDM中,但其仅有二阶精度(Hestholm and Ruud,2002),难以应用于实际.Zhang和Chen(2006)根据实际地表起伏划分曲线网格,实现起伏地表贴体网格FDM.de la Puente等(2014)采用基于完全交错网格的垂向变换网格和模拟地震波建模的网格变形策略,实现了起伏地表条件下的弹性波数值模拟.Shragge 和Tapley(2017)采用广义结构网格法求解张量声波方程.Chen 等(2022)将浸入边界法与平面波透射边界相结合,提出一种适用于起伏地表的浸入吸收边界,该方法基于常规的矩形网格剖分,简单高效.贴体网格是一类曲变网格,适用于同位网格和交错网格,使用范围广且较易生成.但是贴体网格要求地形起伏的函数处处可导,因此不适用于地表起伏剧烈的区域(如崎岖地表、近垂直断崖等地质条件).非结构网格地震波数值模拟的实现方法主要有限元法(Käser and Iske,2005;董兴鹏和杨顶辉,2017)、边界元法(符力耘和牟永光,1994; 何彦锋等,2013)、谱元法(Komatitsch et al.,2000;Cristini and Komatitsch,2012)、间断伽辽金法(Cockburn et al.,2000;Käser and Dumbser,2006)、格子法(徐义和张剑锋,2008;孙辉等,2019)等.这些方法均基于非结构网格或非规则网格算子的全空间离散,能够很好的贴合起伏地表,但是这些方法大多采用全局积分方程弱形式,计算量大,长时间模拟存在不稳定问题.

无网格法是近些年来兴起的一类数值计算新方法,它是有限元等传统方法的重要补充和发展.基于无网格的径向基函数有限差分方法(Radius-Basis-Function Finite Difference Method,RBF-FDM)能够有效避免虚假反射,具有节点剖分随模型速度变化,复杂地质体剖分灵活等优点(刘立彬等,2020).RBF是RBF-FDM中的一类径向对称函数,用于计算差分权重,并求解不同类型的偏微分方程(Tolstykh and Shirobokov,2003; Fornberg et al.,2013; Fornberg and Flyer,2015;Dehghan and Mohammadi,2019).在地震领域,Martin 等(2013, 2015)首次将RBF-FDM应用于声波方程,实现了二维非均匀介质的数值模拟.Takekawa等(2014a,b)研究了RBF-FDM中的自由表面条件,并将哈密顿粒子法引入RBF-FDM中,进行了起伏地表声波数值模拟.Takekawa和Mikada(2018)在RBF-FDM中提出了一种适用于无网格的自由表面边界条件,并将其应用于频域弹性波数值模拟.Li等(2017)进一步将含时-空域优化算子的RBF-FDM应用于声波数值模拟,同时采用混合吸收边界条件压制边界反射.Takekawa等(2016,2018)通过增加相邻节点数目来提高数值模拟精度,将多参数泰勒展开应用于非均匀介质RBF-FDM.Liu等(2019)利用RBF-FDM进行频率域声波方程数值模拟,采用完全匹配层(Perfectly Matched Layers,PML)边界条件来压制边界反射.Duan等(2021)提出了一种基于频散关系和稳定性条件约束节点间隔的自适应节点剖分方法,实现了灵活稳定的RBF-FDM声波数值模拟.Wu等(2021,2022)实现了基于RBF-FDM的声波逆时偏移和全波形反演.

除了起伏地表地震波数值模拟外,国内外已有多位学者对起伏地表标量波场逆时偏移成像方法进行了研究(张晓波等,2022;孙志广等,2023;姚宗惠等,2023).徐义(2008)提出了一种适用于起伏地表的格子法叠前逆时度偏移方法.Zhang 等(2010)提出了一种基于地表起伏非反射递归算法的叠后逆时偏移算法.刘红伟等(2011)给出了一种实现起伏地表叠前逆时偏移的方法,该方法通过简化自由边界条件,结合GPU加速算法实现了高效稳定的起伏地表叠前逆时偏移.Lan等(2014)基于贴体网格地形“平化”策略较好地拟合平滑地表起伏,实现了贴体网格的起伏地表逆时偏移.Liu等(2017)提出了一种适用于起伏地表的非结构网格法叠前逆时偏移.然而,实际地震波是包含纵横波的矢量场,将纵横波的耦合特性综合应用于基于矢量波场的弹性波逆时偏移,弥补了标量波成像在岩性识别和构造解释等方面存在多解性的缺陷.曲英铭等(2015)提出一种基于曲坐标系的分层坐标变换法,实现起伏地表条件下的弹性波叠前逆时偏移.Qu等(2021)将此方法推广到具有起伏地表的棱柱波弹性波逆时偏移方法中.Naveed和Chen(2019)利用曲网格进行波场延拓,利用贴体网格处理起伏地表实现了弹性波逆时偏移.Zhong等(2021)基于Zhang和McMechan(2010)提出的波场分离方法,利用地形起伏相关的滤波器实现了起伏地表弹性波逆时偏移.上述起伏地表弹性波逆时偏移方法均以结构网格剖分为基础,计算量较大.面对当前复杂的勘探目标,进一步研究起伏地表弹性波逆时偏移方法仍具有重要意义.

综上所述,基于RBF-FDM的地震波数值模拟方法研究主要局限于水平地表声学介质,而针对起伏地表弹性介质的RBF-FDM数值模拟的理论和方法研究尚不完善.另外,针对起伏地表弹性波逆时偏移,目前主要采用曲坐标系下的贴体网格等方法进行波场延拓,尽管曲网格能较好的贴合地表起伏,但对于起伏地表的自由边界条件适用性不足.本文在声波RBF-FDM数值模拟基础上,提出一种优化的起伏地表无网格弹性波RTM方法.通过定量分析基于RBF-FDM的弹性波数值模拟的频散关系与稳定性条件,本文发展了高精度的QR径向基函数有限差分方法(QR-RBF-FDM),同时提出一种优化的起伏地表自适应节点剖分技术,形成了新的基于无网格的起伏地表弹性波数值模拟方法.由于起伏地表无网格边界条件加载困难,本文提出了适用于无网格的精确自由边界条件和弹性波混合吸收边界条件.此外,本文将此无网格RBF-FD应用于精确的纵横波场矢量分解公式(Gu et al.,2019),实现了起伏地表弹性波逆时偏移成像.通过对高斯山丘模型,起伏凹陷模型和起伏地表Marmousi-2模型的数值试算验证了本文提出方法在起伏地表弹性波数值模拟与逆时偏移成像中的可行性与有效性.

1 径向基函数有限差分方法弹性波数值模拟

1.1 径向基函数弹性波有限差分方法

各向同性介质中,弹性波二阶位移方程可表示为

(1)

其中,u(x,z,t)和w(x,z,t)表示在空间x、z方向上的位移波场分量,ρ(x,z)表示介质密度,λ(x,z)和μ(x,z)为拉梅常数.在起伏地表条件下,传统的FDM很难满足要求.基于无网格节点离散模型的RBF-FDM能够有效避免虚假反射,具有剖分节点随模型参数变化,起伏界面处理灵活等优势(刘立彬等,2020).本文将RBF-FDM引入弹性波波动方程数值模拟中,以便获得高精度的弹性波地震波场.

RBF-FDM是基于一组局部无序的离散节点,根据与节点距离相关的径向基函数构造差分算子矩阵实现偏微分方程.在任一时刻,若给定Ω区域内存在波场U(x)及无序离散节点x={x1,x2,…,xn},节点位置及参数由点生成方法可知,波场U(x)可由径向基函数φ(r)近似表示为

(2)

式中,r节点间距(中心点与其最近邻近节点的距离),rj=‖x-xj‖表示邻近节点xj距离中心节点x的距离,‖·‖表示欧几里得范数;N为中心节点x差分模板所需的邻近节点个数,根据邻近当前任意点的距离选取差分模板的范围,选择N个最近点,记为差分模板,差分模板中的点按距离远近排序;{c1,c2,…,cN}为待定的RBF-FD系数.根据无网格RBF-FD基本原理,对于任意线性偏微分算子L,中心点节点x0处波场U(x)|x=x0的任意空间偏导数可以表示为

(3)

(4)

其中,径向基函数系数矩阵的条件数是影响RBF-FD差分系数求解的重要因素,为保证差分方程的稳定,条件数的取值范围应在106～108(Flyer et al.,2014;Martin et al.,2015;Li et al.,2017).实际中,该条件数由当前差分模板的节点分布,差分格式,节点间隔和径向基函数类型及形变参数等相关.根据条件数经验取值范围(106～108),在类均匀分布(如图1)的模型中得到了一组节点数与节点间隔和形变参数的最佳取值范围(见表2,参见Li et al., 2017).

表1 径向基函数有限差分方法中常用的径向基函数Table 1 Commonly used radial basis function in RBF-FDM

采用RBF-FDM求解弹性波方程,首先通过时间二阶差分来代替式(1)中的时间微分算子如下：

表2 类均匀分布情况下RBF-FD最优差分节点数、节点间隔和形变参数的选取Table 2 Optimal ranges of parameters for RBF-FD with different numbers of neighbors in quasi-uniform distribution

(5)

其中,上标l为当前时刻,τ为时间步长.采用RBF-FDM差分算子对式(1)中u分量的空间导数进行离散：

(6)

(7)

本文采用RBF-FD中常用的K-D算法(KNN-Search)寻点,能够不占用计算量找到计算域中任一中心节点差分模板中的N个邻近差分节点.

1.2 数值频散分析与稳定性条件

在平面波条件下假设下,弹性波场可以表示为

(8)

(9)

(10)

(11)

对式(11)中两式相乘得：

+β2(A+B)]=0,

(12)

根据式(12)可得：

(13)

(14)

将式(10)中A,B代入式(11)化简可得RBF-FD频散关系为

(15)

其中：

(16)

其中,θ表示相角.

稳定性条件也是衡量微分方程数值求解算法有效性的重要指标.本文根据傅里叶谱分析法(Von Neumann,1950)推导RBF-FD情况下的稳定性条件.此处,将式(1)表示为矩阵形式：

(17)

为了方便起见,进一步将式(17)简写为

(18)

式中U(t)=[u(t),w(t)]T为多分量波场矢量,Q为二阶矩阵.

应用RBF-FD的思想求解弹性波方程,对式(18)中的时间导数进行二阶有限差分离散可得：

+O(Δt2M),

(19)

将式(18)代入式(19)得：

(20)

设G(kx,kz)是Q(kx,kz)的傅里叶变换对,对式(20)进行空间傅里叶变换可得：

U(t+Δt,kx,kz)-2U(t,kx,kz)+U(t-Δt,kx,kz)

(21)

将U(t,kx,kz)的系数矩阵记为

(22)

则式(7)改写为

Un+1=2Un-Un-1+AUn,

(23)

式中n表示当前时刻,根据傅里叶分析方法,令：

(24)

联立式(23)、(24)得：

(25)

式中,E为单位矩阵,B为增长矩阵.

根据Von Neumann稳定性条件可知,当增长矩阵特征值的绝对值小于或等于1时,递归是稳定的,即时间稳定性条件是增长矩阵B的谱半径不大于1.由式(25)可知,B的特征值β应满足：

|βE-B|=|β2E-βA-2E|=0,

(26)

由式(26)可知,β的取值与矩阵A有关,因此先求取的A特征值γ：

(27)

根据式(22),若G的特征值为η,求取A特征值γ需要首先求取G的特征值η.根据矩阵Q的表达形式,η应满足：

(28)

(29)

根据行列式规则式(28)可化简为

+Dkzz)2=0,

(30)

利用2次多项式的求解方法计算得到η的2个特征值为

(31)

由式(27)可知,矩阵A特征值γ为

(32)

由于A为2阶方阵,且有2个互异的特征值,因此矩阵A可对角化,有：

A=PΛP-1,

(33)

式中P为可逆矩阵,Λ为A的特征值构成的对角阵,具体表示为

Λ=diag{γ1,γ2},

(34)

由式(27)可得：

|β2E-βΛ-2E|=0,

(35)

由式(35)可得B的特征值β和A的特征值γ的关系式为

β2-βγ-2=0,

(36)

即：

(37)

(38)

若系数矩阵A的特征值γ的绝对值小于或等于1,则满足稳定性条件.对比|γ1|和|γ2|表达式可知,RBF-FD差分格式需满足的稳定性条件是：

(39)

2 QR径向基函数有限差分方法

在常规RBF-FDM中,形变参数ε趋近于0时,径向基函数趋于平坦而拟合精度越高.然而平坦的径向基函数会导致差分矩阵的病态,影响波场延拓的稳定性.Flyer等(2012)研究了含泰勒多项式的RBF-FDM,提高了数值模拟的稳定性,但该方法依赖模型参数的取值,在处理复杂构造模型时,构造泰勒多项式所增加的参数将会导致差分矩阵的病态,使得波场传播不稳定.为了提高模拟精度且不落入病态问题,Fornberg等(2013)提出一种无需形变参数的高斯QR基函数,将其成功用于插值问题.本文将高斯QR基函数的思想引入RBF-FD中,在笛卡尔坐标系下构造差分矩阵,将高斯径向基函数转换至极坐标系,并展开为切比雪夫多项式的形式,利用QR分解法求解由切比雪夫多项式构成的差分矩阵,有效解决了径向基函数的形变参数取值和矩阵病态问题,在保证精度的同时增强波场传播的稳定性.

QR径向基函数法利用切比雪夫多项式T(x)对高斯径向基函数φ(r)=e-i(εr)2级数展开(Fornberg et al., 2011),并将其变换至极坐标系(r,θ)得到：

(40)

其中切比雪夫多项式表示为

(41)

其中,j≥0,0≤m≤(j-p)/2;当p=0时,j为偶数,当p=1时,j为奇数;公式(40)中尺度因子dj,m为

(42)

当形变参数ε适当取值时,尺度因子dj,m的取值范围为[0,1];公式(40)中参数cj,m和sj,m为

(43)

其中参数b0=1,bm=2,m>0;参数t0=1/2,tj=1,j>0;F2为超越方程组,其中参数αj,m=(j-2m+p+1)/2,βj,m=[j-2m+1,(j-2m+p+1)/2];由此,公式(40)可改写为矩阵形式：

(44)

其中,矩阵C由参数cj,m和sj,m构成,三角阵D由尺度因子dj,m构成.

利用QR分解对矩阵C进行分解代入公式(44)得：

φ(x)=Q·R·D·T(x),

(45)

将公式(45)代入公式(40)得：

(46)

(47)

求解公式(47)可得到优化的差分系数.由于切比雪夫多项式的定义域为r∈[-1,1],计算差分系数时需将当前点的差分模板预先缩放到[-1,1]的范围内.

图1 类均匀节点剖分示意图

这里采用一个均匀介质模型研究QR径向基函数有限差分方法(QR-RBF-FDM)的拟合精度、稳定性及其频谱分析.介质纵波速度为3000 m·s-1, 横波速度为1732 m·s-1,模型大小为1000 m×1000 m,节点间隔为10 m,时间步长为0.5 ms,最大记录时间为1.0 s.采用20 Hz雷克子波作为纵波震源,位于近(500 m, 20 m)处.图1给出了类均匀节点分布图,红点和蓝点分别表示任意中心差分点及其对应的无网格差分模板包含的N个差分节点.图2给出不同方法得到的点(500 m, 0 m)处的波形图及误差曲线图.其中,解析解(analytical_solution)由Cagniard-De Hoop法(De Hoop, 1960)得到,FDM表示时间2阶空间10阶精度的有限差分.由图2a单道记录与图2b振幅误差曲线图可见：QR-RBF-FDM 和 FDM的波形和振幅与解析解基本吻合,而RBF-FDM的模拟结果取决于形变参数取值.当形变参数为0.01时,RBF-FDM的波形和振幅与解析解基本吻合,而当形变参数为0.1或0.001时,RBF-FDM的模拟结果与解析解偏差很大.图2c为单道记录对应的频谱图,QR-RBF-FDM能够有效抑制形变参数导致的差分矩阵病态问题,获得与FDM和解析解相近的频谱;常规的RBF-FDM采用不适合的形变参数将引起波形震荡,时间空间频散分别出现在波前的浅部和尾部,同时引起了频谱畸变,对频谱影响较大.由此可知,QR-RBF-FDM比常规的RBF-FDM具有更高的模拟精度,且形变参数的取值不影响QR-RBF-FDM的结果.图3为QR-RBF-FDM、常规RBF-FDM与FDM的计算效率对比.与FDM相比,QR-RBF-FDM与RBF-FDM的计算成本主要来源于差分系数计算和差分方程迭代两个方面.其中,差分方程迭代过程与FDM基本一致,而尽管二者仅需在整个计算过程中计算一次差分系数,但由于每个节点进行差分系数矩阵求解仍然需要消耗巨大的计算成本,因此QR-RBF-FDM与RBF-FDM在计算成本和计算时间方面为FDM的4倍左右(如图3).QR-RBF-FDM求解差分系数过程计算相对常规的RBF-FDM更加复杂,因此计算时间较长,但考虑到无需选择形变参数,因此可避免因形变参数选择不当而产生的计算冗余.

图2 均匀模型(500 m, 0 m)处不同方法(a) 单道记录波形; (b) 振幅误差曲线; (c) 对应的频谱.

图3 均匀模型不同方法归一化计算时间对比图

3 起伏地表自适应节点剖分方法

节点剖分是RBF-FDM的第一步,也是连接速度模型与数值模拟算法的关键.Duan等(2021)提出的自适应节点剖分法利用频散关系和稳定性条件约束节点间隔的取值范围,能够在矩形计算域内生成适应模型速度的非均匀节点分布.但考虑起伏地表的不规则边界与地下介质的不均匀性,Duan等(2021)的方法将不再适用.因此,在Duan等(2021)的基础上,本文提出一种自适应起伏地表节点剖分法,边界节点贴体起伏地表,并在内部节点剖分过程中,对起伏地表和地下构造变换剧烈的区域依据速度变化的梯度大小进行自适应调节.起伏地表节点剖分方法(见图4)采用如下的步骤：

(2)利用自适应节点剖分法进行内部节点剖分,剖分方式如图4bi—iii所示.(i) 设定最下层边界节点为潜在点位.(ii) 根据速度模型参数确定各起始点位的节点间距,设置以起始点位为中心,半径为节点间距的圆域.(iii) 删除此圆域内所有其他的潜在点位,在删除节点的圆周范围内放置等距的新点位.由下层潜在点位迭代直至与边界节点距离小于节点间距的一半d/2.记录内部节点的初始坐标(ξ,η)对应的速度模型参数,如图4b(iv)所示.

(3)删除与边界节点间距小于d/2的内部节点,如图4c所示.

(4)根据自适应节点控制方程(式(51)),将内部节点由初始坐标调整到新的坐标(x,z),通常,需要四次以上的迭代获得满足精度要求的节点分布,得到每个节点的坐标及对应的速度模型参数,如图4d所示.

图4 起伏地表自适应节点剖分过程示意图

为了计算起伏地表上节点的法线方向,步骤1中地表节点的坐标变换到极坐标系：

(48)

式中(x,z)表示笛卡尔坐标系中地表节点的位置向量,(r,θ)表示极坐标系下地表节点的位置向量,取速度模型左上顶点作为极点,则关于极角θ的梯度为

(49)

由此可知,法向量的分量nx和nz可表示为

(50)

其中sgn(·) 表示常用符号函数,g表示梯度.在步骤4中,引入自适应网格理论(Saltzman and Brackbill, 1982)中的控制方程对自适应节点剖分法生成的无网格节点位置进一步调整,控制方程为

(51)

(52)

J=xξzη-xηzξ,

(53)

其中J表示雅可比多项式,其中一阶导数∂x/∂ξ,∂x/∂η,∂z/∂ξ和∂z/∂η简写为xξ,xη,zξ和zη.

无网格与基于坐标变换法的贴体网格和曲网格各有优劣.(1)贴体网格和曲网格是起伏地表条件下的地震波模拟中的重要方法,其网格易于生成且在曲坐标系下可直接使用水平地表假设的自由边界条件.与这两种网格相比,无网格的优势在于无需坐标变换,直接采用弹性波、声波等方程直接进行差分,后期研究中更适用于复杂的黏性介质或各向异性方程.(2)贴体网格和曲网格在生成过程中,仅能通过一些密度函数进行疏密调节,笔者认为此种调节疏密通常是不可控的,因为网格数不变的情况下调节一个部分疏密必将使得另一部分密疏,一定程度上将导致数值频散以及稳定性问题.与这两种方法相比,无网格以模型速度、密度等参数为依据,直接生成满足模型剖分条件的节点分布,不存在调节后一部分点过密一部分点稀疏的现象,因此相对更加具有数值稳定性.(3)贴体网格和曲网格本质上仍然是固定的网格间隔,这就造成在深层、中深层进行数值模拟,网格过密造成冗余.通过文献调研可知,贴体网格往往需要与变网格相互结合以提高其计算效率.与之相比,无网格节点在剖分过程中根据速度、密度等参数可自适应调节节点间距,因此在深层、中深层等复杂构造中可由细网格自然过渡到粗网格.

4 与地形相关的边界条件

4.1 无网格自由边界条件

RBF-FDM的优势在于处理起伏自由边界.在自由表面是固体地球与空气层之间的物性突变面,常规自由边界条件通常假设垂直自由表面的应力为0(Lan and Zhang, 2011),二维水平自由表面(z=0)处的表达式为

(54)

地形起伏会明显影响地震波场的传播,此时水平自由边界条件的假设将不适用.为了在RBF-FDM中精确处理起伏地表,需要将无网格节点自然放置在起伏的曲面地形上,利用起伏地表自适应节点剖分方法获得起伏界面的法线分量,在此基础上重新推导起伏界面的精确自由边界条件.图5显示了规则网格和无网格的剖分示意图.

图5 网格剖分示意图(a) 规则网格; (b) 无网格节点.红色曲线表示起伏地表, 蓝色节点表示内部节点,红色节点表示地表节点,箭头表示地表节点的法向量方向.

因此,本文采用Dirichlet边界条件推导出无网格节点的自由曲面边界条件,将其应用于贴合实际地形的RBF-FDM具有更高的精度.Dirichlet边界条件如下：

(55)

其中σxx和σzz为正应力,σxz为切应力.在此基础上,通过自适应节点生成方法求出了与起伏地表垂直的外法向x分量nx和z分量nz.因此,应力分量可以表示为

(56)

将式(56)代入式(55),经过化简可以得到起伏地表自由边界条件：

(57)

图6 差分模板示意图红点为中心点,蓝点为中心点对应的差分点,椭圆区域为边界节点差分模板,圆形区域为内部节点差分模板.

其中,地表界面上的法向量的分量nx和nz可式(56)计算得到.在规则网格FDM中,通常需在自由表面上方引入虚拟点以满足自由表面处边界条件的差分格式.与之相比,在RBF-FDM中,计算节点可以放在任意位置进行差分,边界节点完全贴体起伏地表分布.因此无需特殊处理,在边界节点处采用如图6椭圆区域所示差分模板,利用QR-RBF-FDM获得边界节点处的空间导数并运用于式(57)实现无网格的自由边界条件.其中一阶空间导数可由RBF-FDM近似得到：

(58)

其中x0为地表节点的坐标,mj和nj表示由QR-RBF-FDM计算得到的一阶差分系数,同理可得垂直分量w的一阶空间导数的表示形式.本文由Dirichlet边界条件出发,推导了基于RBF-FD的自由边界条件,相对于与应力镜像法等常规方法而言,该方法无需设置虚拟层,在笛卡尔坐标系下直接求解边界条件,处理方式简单且计算精度高.本文在数值算例中将对此进行验证.

4.2 无网格弹性波混合吸收边界条件

利用RBF-FDM进行起伏地表弹性波逆时偏移时,必须考虑起伏地表情况下的边界吸收问题.基于Li 等(2017)提出的分段光滑曲线边界节点生成方法,本文在现有边界节点基础上,基于弹性波Clayton和Engquist(CE)单程波方程推导了无网格起伏地表条件下的混合吸收边界条件.起伏地表条件下的弹性波CE单程波方程表示为(任志明,2016)：

(59)

(60)

将式(59)写为位移分量方程的形式可得：

(61)

通过引入边界区域每个节点对应的最近内点,为单程波方程提供稳定的差分：

(62)

其中,位移分量u和w的上标-1,0,1表示前一时刻,当前时刻与后一时刻,边界区域的节点(表示为下标B)与其对应的最近内点(表示为下标I)构成的差分项Al和Bl(l=1,2,3,4)为

(63)

(64)

(65)

在此基础上,在边界区域的同时采用双程波方程与上述单程波方程,即可构成起伏地表条件下的无网格混合吸收边界.混合吸收边界条件实现步骤(Liu and Sen, 2010)如下：

(3)对于边界区域,将双程波波场和单程波波场加权平均得到最终的波场：

(66)

其中,加权系数为

aBi=(i-1)/N.

(67)

5 基于径向基函数有限差分方法的弹性波逆时偏移

在无网格模型参数化的基础上,我们开展了基于RBF-FDM的弹性波逆时偏移成像方法应用研究.逆时偏移的理论基础是Claerbout(1971)提出的时间一致性原理,即反射面存在于地层内下行波波至时间与某一上行波波至时间相一致之处.因此QR-RBF-FDM弹性波逆时偏移成像可分为三个部分：一是基于QR-RBF-FDM震源波场的正向延拓,二是基于QR-RBF-FDM检波点波场的逆向延拓,三是应用成像条件获得偏移成像结果.波场分离是弹性波RTM至关重要的环节.Gu等(2019)解决了亥姆霍兹分解导致的波场振幅和相位畸变,建立了精确的波场分解和重新组合方程,并实现了3D标量波弹性波逆时偏移,该方法能够直接获得具有正确振幅、相位和物理单位的偏移结果.本文基于Gu等(2019)的方法,我们推导得到在无网格模型参数化条件下通过RBF-FDM实现的精确纵横波波场分离方程,在正向和逆向波场延拓过程中,利用无网格震源波场和检波点波场进行精确的纵横波波场分离,并应用互相关成像条件得到成像结果.

根据Gu等(2019)的方法推导,精确的波场分离方程如下：

(68)

(69)

其中,vp和vs为P-和S-速度.由于无网格法的本质是局部无序的离散节点,无法直接求解式(68)、(69)中的积分.考虑到矢量场V仅有子波项与时间有关,因此在正向和反向波场计算过程中,对震源子波做二次积分提取式(68)、(69)中的积分,矢量场中P和S势则已包含积分,即可去掉式(68)、(69)中的积分.将式(68)代入式(69)可得：

(70)

根据散度、旋度和梯度的关系求解式(68)得：

(71)

采用QR-RBF-FDM对式(71)进行离散,代入式(6)二阶空间算子可得：

(72)

(73)

其中x=(x,y,z)为空间向量;其中sm和rn分别表示源波场的m波模式和接收器波场的n波模式.

6 数值算例

6.1 高斯起伏模型

本节通过两类高斯起伏地表模型(凸型模型和凹型模型)验证在自由边界条件情况下本文方法的有效性,如图7a、b所示.模型大小为2 km×1.5 km,其中凸型高斯函数为

(74)

凹型高斯函数为

(75)

我们分别对两类模型采用自适应起伏地表节点剖分法和规则网格剖分法进行离散化.图8展示了自适应起伏地表节点剖分法离散的结果,由于本测试选用均匀速递模型以测试起伏自由地表下QR-RBF-FDM的可靠性和稳定性,在速度层上为真空,地表采用自由地表条件,矩形边界采用无网格弹性波混合吸收条件.模型参数设置如下：纵波速度为3000 m·s-1,横波速度2000 m·s-1,而对应的自适应起伏地表节点剖分结果呈现出节点间隔相同的特征,其间隔大小为10 m(事实上当速度非均匀条件下,自适应起伏地表节点剖分结果的节点间隔会呈现出不同速递区域的差异化现象,我们将在后续模型测试中展示).图9展示了规则网格剖分结果,其正交方向的网格间隔统一为10 m.图9中红线部分为高斯起伏界面位置,由于规则网格剖分无法精确捕捉真实的界面位置,因此实际数值模拟过程中,将选取纵向最接近真实界面位置的点作为该模型的自由边界顶点.

两类高斯模型测试过程中,时间间隔统一为0.5 ms,采样点数为2001,总时长为1 s.震源选取主频为20 Hz的雷克子波,震源位置坐标为 (0.4 km, 0 km).与常规FDM相比(如图10所示),图10展示了利用QR-RBF-FDM所得的不同时刻凸型高斯起伏地表模型波场快照,图8d和图11中对地震相做了详细的分析,并在记录上对波的类型做了相应标注.对比可知,常规FDM和本文方法均能够较好地处理起伏地表,波场快照中能够清晰地呈现出直达P波、反射PP波、反射PS波和面波R.但与常规FDM不同,QR-RBF-FDM在起伏构造中没有产生台阶散射,模拟结果更准确,反射波同相轴更清晰.对比二者数值结果可得以下结论：(1)在0.2 s之前,波在水平地表传播,QR-RBF-FDM和常规FDM的模拟结果基本一致;(2)0.2 s至0.6 s之间,波经过了起伏地表,常规FDM模拟的波场中产生大量“波纹状”的次级散射波,这些是由规则网格离散所造成的台阶散射效应.与之不同,QR-RBF-FDM的地表节点分布完全贴合了起伏地形,因此QR-RBF-FDM模拟的波场不存在上述次级散射,展现出了该方法对于复杂区域模拟的优势;(3)0.6 s之后,波穿过起伏地表并在水平地表继续传播,该阶段QR-RBF-FDM和常规FDM的模拟结果基本相同,该结论与0.2 s之前的情况一致.图12和图13展示了凹型高斯起伏地表模型波场模拟结果,对比该结果可以得出与凸型高斯起伏地表模型类似的结论.进一步地,对比图11和图13中波场快照可知,凸起和凹陷模型在传播过程中均会产生直达波P、反射纵波PP、反射横波PS和面波R.与凸起模型相比,0.6时地震波围绕凹陷模型成回转形式向前扩散,无论直达波或是面波在正方向产生了很强的散射.图14和图15分别为两类模型采用QR-RBF-FDM和FDM模拟得到的地震记录,从中可以很清晰地看出QR-RBF-FDM相较于FDM在消除阶梯散射上的优势,且面波同相轴清晰.综上可得,常规FDM和QR-RBF-FDM方法有很大的区别,QR-RBF-FDM得到的波场在起伏界面处非常稳定,并且不会产生大量次级散射,能够准确地对自由边界条件下的起伏地表模型进行数值模拟.

图7 均匀介质起伏地表模型(a) 凸型高斯模型; (b) 凹型高斯模型.

图8 均匀介质起伏地表模型自适应起伏地表节点剖分结果(a) 凸型高斯模型; (b) 凹型高斯模型.

图10 凸型高斯起伏地表模型常规FDM波场快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

图11 凸型高斯起伏地表模型QR-RBF-FDM波场快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

6.2 起伏凹陷模型

本节利用起伏凹陷模型验证非均匀模型下自适应起伏地表节点剖分法和QR-RBF-FDM的稳定性和适应性,此时地表及矩形边界采用无网格弹性波混合吸收条件.该模型见图16,模型大小为6 km×4 km,地表起伏函数为6.1节两类高斯函数的组合.图17展示了自适应起伏地表节点剖分法离散的结果.与图17b中常规节点剖分方法采用统一尺度步长相比,图17a基于自适应起伏地表节点剖分方法中节点间距随着速度增加而变大,范围在10 m至15 m之间.基于该模型进行数值模拟,时间步长为0.5 ms,记录时间长度为3 s,选取主频为20 Hz的雷克子波作为震源时间函数,检波点均匀布设于地表,接收范围为4 km,相邻检波点间距为10 m.

图18给出了不同时刻QR-RBF-FDM模拟的波场快照,震源位于地表横向3 km处.对比不同时刻波场可知,在近地表区域模拟的波场不存台阶散射,波传播至深层区域时(即节点间隔增大区域)仍然稳定且无频散现象.由图19的不同炮点位置(地表x方向坐标1 km, 3 km 和5 km处)激发的地震记录可以观察到,不论在地表任何位置激发的震源,QR-RBF-FDM都能很好地处理起伏地表,采用无网格弹性波混合吸收条件后地震记录中可以清晰地看出直达波,反射波,转换波的发育而无面波与散射波,验证了方法的正确性.综上可知,本文提出的QR-RBF-FDM不仅通过灵活的自适应起伏地表节点剖分法提供一种稳定的离散节点模型,而且与常规RBF-FDM相比无需选取形变参数,大幅提升了无网格算法的计算效率和计算精度.与常规曲坐标系算法相比,QR-RBF-FDM差分形式更为简单,更易于推广.

图12 凹型高斯起伏地表模型常规FDM波场快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

图13 凹型高斯起伏地表模型QR-RBF-FDM波场快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

图14 凸型高斯起伏地表模型QR-RBF-FDM和FDM地震记录(a) QR-RBF-FDM水平分量; (b) QR-RBF-FDM垂直分量; (c) FDM水平分量; (d) FDM垂直分量.

图16 起伏凹陷模型(a) 纵波速度模型; (b) 横波速度模型.

图17 起伏凹陷模型节点剖分结果(a) 自适应起伏地表节点剖分结果; (b) 常规节点剖分结果.

图18 起伏凹陷模型QR-RBF-FDM波场快照(a) 0.5 s水平分量; (b) 1.0 s水平分量; (c) 1.5 s水平分量; (d) 0.5 s垂直分量; (e) 1.0 s垂直分量; (f) 1.5 s垂直分量.

图19 起伏凹陷模型QR-RBF-FDM地震记录炮点位于地表1 km：(a)水平分量;(b)垂直分量; 炮点位于地表3 km：(c)水平分量;(d)垂直分量; 炮点位于地表5 km：(e)水平分量;(f)垂直分量.

图20 起伏凹陷模型多炮成像结果(a) PP成像; (b) PS成像.

图21 起伏地表Marmousi-2模型(a) 纵波速度模型; (b) 横波速度模型.

此处,我们进一步地利用上述方法实现了ERTM,采用60炮主频20 Hz的雷克子波作为震源时间函数在地表附近放炮,中间激发两端接收,炮间距100 m.最大记录时长3 s,每炮400道接收,道间距10 m.图20给出了多炮叠加后的PP和PS剖面,从中可以看出,QR-RBF-RTM可灵活地处理起伏地表,且地下界面位置准确,同相轴清晰连续,该数值结果验证了本文方法的正确性.

6.3 起伏地表的Marmousi模型

本节利用经典Marmousi-2模型验证复杂模型下自适应起伏地表节点剖分法和QR-RBF-FDM的稳定性和适应性,此时地表采用无网格弹性波混合吸收条件.模型见图21,模型大小为13.6 km×2.8 km,采用Foothill模型的起伏地形作为Marmousi的地形起伏.图22展示了自适应起伏地表节点剖分法离散的结果.图22基于自适应起伏地表节点剖分方法中节点间距随着速度增加而变大,范围在5～18 m之间.基于该模型进行数值模拟,时间步长为0.5 ms,记录时间长度为4.5 s,选取主频为20 Hz的雷克子波作为震源时间函数,检波点均匀布设于地表,接收范围为10 km,相邻检波点间距为10 m.

图23给出了不同时刻水平和垂直分量的波场快照,震源位于地表横向3 km处.图24为纵横波分离的波场快照,其中左列为纵波波场分量,右列为横波波场分量.图23表明本文方法对于地表高程变化剧烈的复杂地表模型能够实现高精度数值模拟,图24表明QR-RBF-FD应用于分离方程能够有效分离纵横波.由图25地震记录可见,在复杂模型下弹性波场仍然稳定,且无矩形网格剖分引起的虚假反射现象.由图25的不同炮点位置(地表x方向坐标3.2 km, 6.5 km和9.8 km处)激发的地震记录可以观察到,不论在地表任何位置激发的震源,在记录中均可以清晰地看出,直达波、反射波和转换波的发育.由于起伏地表的存在,反射波和转换波的双曲线特性发生规则扭曲,验证了在复杂构造和起伏地表情况下QR-RBF-FDM的适应性.

此处,我们进一步地利用上述方法实现了ERTM,采用135炮主频20 Hz的雷克子波作为震源时间函数在地表附近放炮,中间激发两端接收,炮间距100 m,最大记录时长4.5 s,每炮1000道接收,道间距10 m.图26给出了多炮叠加后的PP和PS成像结果剖面,从中可以看出,成像结果与偏移速度模型的构造界面信息吻合,且地下界面位置准确,无论是PP成像还是PS额成像,同相轴均清晰连续,具有较高的信噪比,表明基于QR-RBF- ERTM方法对起伏地表下得到复杂模型具有良好的适应性.

7 结论

本文提出一种自适应节点剖分法实现起伏地表无网格剖分,在此基础上提出QRRBF-FDM,推导出弹性波模拟中RBF-FDM的频散关系和稳定性条件.同时,提出了一种适用于起伏地表的无网格起伏自由地表条件和弹性波混合边界条件,实现了基于无网格的起伏地表弹性波数值模拟.此外,推导出一种适用于无网格的精确矢量波分离方法,实现了纵横波分离的无网格弹性波逆时偏移方法.通过对高斯起伏模型,起伏凹陷模型和起伏地表的Marmousi-2模型的数值模拟和偏移试算,得出到以下几点认识：

(1)自适应节点剖分法不仅能与地表完全匹配,而且节点分布密度随地下介质的速度变化而变化,在速度大的区域节点分布稀疏,速度小的区域节点分布密集.同时采用自适应节点控制方程对节点位置进行调整,保证节点分布的均匀性.基于该方法,利用无形变参数扰动QR-RBF-FDM进行数值模拟,既保持了传统FDM的高精度,又弥补了传统FDM在自由起伏地表处理中缺乏的灵活性的缺点.

(2)在起伏地表条件下,贴体网格与曲网格是模拟地震波的重要方法,其网格易于生成,而且可利用水平地表假设的自由边界条件,在曲坐标系下直接使用.无网格具有利用弹性波、声波等方程进行直接差分的优点,无需坐标变换,后期研究中更适用于复杂的黏性介质或各向异性方程.需要指出的是,无网格数值模拟的差分模板仍基于FDM的数值差分思想,其传播的稳定性依赖于差分模板的均匀性和连续性,因此在现阶段仍然无法突破梯度非连续模型,这将是后续研究工作的重点与难点.

(3)基于Gu等(2019)的基础上,本文推导并实现了无网格模型参数化条件下的起伏地表弹性波逆时偏移.模型试算结果表明,该方法对地表附近和复杂深层成像效果较好,而且震源的干扰效应能准确反演起伏地表形态.

(4)应用QR-RBF-FDM,以无网格模型参数化为基础,实现逆时偏移,对未来起伏地表复杂构造、复杂介质处理具有广阔的应用前景;同时,QR-RBF-FDM也可直接应用于最小二乘偏移、全波形反演.

图22 起伏地表Marmousi-2模型自适应起伏地表节点剖分结果

图23 起伏地表Marmousi-2模型RBF-FDM波场快照(a) 1.0 s水平分量; (b) 1.0 s垂直分量; (c) 2.0 s水平分量; (d) 2.0 s垂直分量; (e) 3.0 s水平分量; (f) 3.0 s垂直分量.

图24 起伏地表Marmousi-2模型RBF-FDM波场快照(a) 1.0 s纵波分量; (b) 1.0 s横波分量; (c) 2.0 s纵波分量; (d) 2.0 s横波分量; (e) 3.0 s纵波分量; (f) 3.0 s横波分量.

图25 起伏地表Marmousi-2模型RFB-FDM地震记录炮点位于地表3.2 km： (a) 水平分量; (b) 垂直分量; 炮点位于地表6.5 km： (c) 水平分量; (d) 垂直分量; 炮点位于地表9.8 km： (e) 水平分量; (f) 垂直分量.

图26 起伏地表Marmousi-2模型弹性波逆时偏移结果 (a) PP成像; (b) PS成像.