基于单频走时敏感度核函数的三维波动方程初至层析

冯波，许荣伟，王华忠，罗飞

1 同济大学海洋与地球科学学院波现象与智能反演成像研究组，上海 200092 2 中国石油化工股份有限公司石油物探技术研究院, 南京 211103

0 引言

在地震勘探中，走时层析技术是宏观速度建模的主要手段.对于陆上地震资料，尤其是地下和地表“双复杂”的山前带探区，原始单炮道集中的反射地震同相轴基本不可见，有效反射信号的振幅远小于近地表散射噪声(如各阶面波及近地表散射体波).此时，地震初至是唯一能够直接识别的有效信号，可以用于反演近地表模型参数(如速度、各向异性参数、甚至品质因子等).

刘玉柱等(2009)指出，菲涅尔体的边界仅能在匀速模型中解析表达.对于变速模型，只有波动方程才能精确描述有限频带地震波的一阶绕射效应(Wu and Toksöz, 1987).Luo和Schuster(1991)基于Born近似和互相关时差测量准则，给出了基于波动方程的走时Fréchet导数(又称为走时敏感度核函数，Marquering et al., 1999；Dahlen et al.，2000；Hung et al.，2000；Zhao et al.，2000).理论分析可知，Born近似成立条件十分苛刻，要求速度异常体的尺度及扰动强度都很小(Wu and Zheng, 2014).相比于Born近似，Rytov近似直接对散射波的相位扰动进行线性化近似，能够更好地描述速度扰动引起的前向散射波的相位扰动(Wu，2003).因此，也可以通过Rytov近似构造单频谐波走时扰动(又称为相位延迟)与模型扰动的线性关系，进而构造有限频走时敏感度核函数(Spetzler and Snieder，2001，2004；Jocker et al.，2006；Xie and Yang，2008，刘玉柱等，2009；Xu et al.，2015；冯波等，2019).

然而，Born近似与Rytov近似均隐含了“弱散射”假设，要求速度扰动远小于背景速度(Beydoun and Tarantola, 1988).为克服弱散射近似的制约，Feng等(2019)提出了适用于强扰动模型的广义Rytov近似(Generalized Rytov Approximation，GRA)，更适用于走时反演.冯波等(2021)将GRA理论应用于透射波走时敏感度核函数的构建，并采用二维模型数据验证了基于GRA的波动方程走时层析方法的有效性.

在数值计算方面，目标函数(Luo and Schuster, 1991; Tromp et al., 2005)的梯度可以采用伴随状态理论推导，通常可以表示为一个入射场与伴随场的零延迟互相关(Lailly, 1983; Liu and Tromp, 2006, 2008; Fichtne et al., 2006)，这与全波形反演(FWI，Tarantola, 1984)的梯度计算方式一致，仅区别于伴随震源的定义不同.考虑到入射场和伴随场的时间外推方向相反(入射场沿时间正向传播，而伴随波场逆时传播)，为同时获得某一时刻的入射场和伴随场从而计算梯度，需要采取一些策略从而兼顾计算和存储：如最优检查点方法(Symes, 2007)，随机边界条件(Clapp, 2009)，波场重构策略(Symes et al., 2008; Feng and Wang, 2012)等.上述方法中，波场重构策略在计算精度、重计算比和数据I/O量三个方面取得了较好的平衡，因此广泛应用于逆时偏移(RTM)及FWI中.

对于三维波动方程走时层析，即使采用了波场重构策略计算梯度，其计算量仍然比射线层析提高了几个数量级，严重限制了该技术的实用化.本文提出用单频梯度代替传统的有限频梯度，降低内存需求并提高计算效率，从而解决三维波动方程走时层析中遇到的计算和存储问题.本文结构如下：在第1节中，首先引入了单频走时敏感度核函数的概念，结合随机边界条件以及离散Fourier数值积分，给出了单频走时核函数的高效计算方法；接着，将单频走时核函数应用于目标函数的梯度计算，在显著降低内存需求的同时提高了计算效率.在第2节中，首先对比了匀速模型中单频走时核函数和有限频核函数的区别；接着，采用三维Overthrust模型数据，分别计算了传统的走时目标函数梯度(即有限频梯度)以及单频梯度.数值测试结果表明，单频梯度具有良好的背景模型更新能力.随后，用三维Overthrust模型数据验证了单频梯度反演的有效性.文章的第3节针对走时层析中的反演策略、数值计算和存储方面展开了讨论，第4节给出了结论.

1 理论

1.1 单频走时敏感度核函数

基于广义Rytov近似(Feng et al., 2019)，冯波等(2021)给出了如下有限频带地震信号的走时扰动与慢度模型扰动的线性关系：

(1)

其中，xr和xs分别为检波器和震源坐标，Δs(x)=s(x)-s0(x)为慢度扰动，s0(x)=v0(x)-1，s(x)=v(x)-1分别为背景慢度模型和扰动后的慢度模型，K(x;xr,xs)称为基于广义Rytov近似的带限走时敏感度核函数(或称为有限频走时敏感度核函数)，表示为：

(2)

其中，u0(x,t;xs)为入射波场，走时伴随场λ0(x,T-t;xr)由走时伴随震源(冯波等，2021)逆时传播生成.

(2)式中的核函数可以通过波场重构策略(表1)高效计算，需要一次波场正演模拟(计算入射场)以及两次波场逆时传播(震源波场逆时重构及伴随场逆时传播).对于大规模三维问题，即使仅存储6个面上的边界波场，依然会占用数十甚至上百GB的内存空间，严重影响梯度的计算效率.

表1 基于波场重构策略的有限频敏感度核函数计算流程

利用Parseval定理，可以将(2)式中对时间的积分转化为对频率的积分：

×(λ0(x,f;xr))*df}.

(3)

只要对频率的采样满足采样定理，(3)式与(2)式是等价的.例如，假定地震波正演模拟的时间步长dt=2 ms，记录时长T=4 s(道长ns=2000).根据采样定理，频率采样间隔df满足：df=1/T=0.25 Hz.若地震数据的高截频fmax=50 Hz，此时需要nf=fmax/df=200个单频波场即可实现核函数计算.Wang和Dong(2021)通过计算频率域波场(多个频率)，实现了基于高斯-牛顿反演方法的二维波动方程初至层析.Wang等(2021，2022)进一步将该方法推广至各向异性介质，提出了单频反演策略，并用二维数值实验验证了方法的有效性.然而，对于三维问题，存储全部或部分的频率域三维波场仍然需要占用巨大的存储空间，难以全部载入内存.

考虑到影响地震波一阶散射效应的主要区域集中在连接炮检的中心射线邻域的第一菲涅尔带内，且在光滑缓变模型中地震波主频对应的单频菲涅尔带内，单频核函数与带限核函数基本一致(刘玉柱等，2009)，本文提出用单频核函数代替带限核函数：

Kf(x,fm;xr,xs)=2s0(x)Re{((i2πfm)2u0(x,fm;xs))

×(λ0(x,fm;xr))*},

(4)

其中，fm代表地震波的主频.

对比单频及带限核函数的空间形态可知，用单频核函数近似带限核函数的误差主要来自第一菲涅尔带外的空间振荡.对于波动方程初至层析，地震观测系统的多次覆盖特性有助于第一菲涅尔带外的能量相互干涉(详见图6)，因此用单频核函数实现背景速度反演在理论上具有可行性.

为了更好地实现层析核函数在第一菲涅尔带外的相互干涉，本文引入随机边界条件用于入射场的正演模拟和伴随场的逆时传播.通过在随机边界内生成随机模型(互相关函数为高斯或指数函数)，并令其特征尺度(Baig et al., 2003)与入射波长相仿，使得入射场和伴随场在随机边界内的空间相关性显著降低.

在数值计算方面，随机速度模型的引入仅仅修改了模型边界附近内的速度模型，无需引入额外的数值吸收边界(如PML(perfectly matched layer，完美匹配层)边界条件等，对于三维问题，采用随机边界后波场模拟无需单独处理6个面和8个角点的反射).表2中给出了将随机边界条件用于单频核函数(4)的具体计算流程.由于无需单独处理吸收边界，新的单频核函数计算方法的计算量进一步降低，且算法更加简洁(更适用于GPU加速).

表2 基于随机边界条件的单频敏感度核函数计算流程

1.2 走时层析反问题

基于走时残差L2范数的误差泛函可以表示为：

(5)

其中，m∈M为模型空间M中的向量(本文中m为慢度模型)，Δt∈D为数据空间D中的向量.

(5)式可以用梯度导引算法迭代求解：

mk+1=mk-αkP(g(mk)),

(6)

其中，k为迭代次数，αk为步长，可以通过线性搜索方法计算，g(mk)为泛函梯度，P为光滑算子，用于提取梯度中的低波数成分.

根据(5)式，泛函梯度为：

(7)

其中，mk为第k次迭代时的慢度模型，Δtk为第k次迭代时的走时残差，K为带限走时敏感度核函数(冯波等，2021).

将泛函梯度(7)中的带限核函数用单频走时敏感度核函数替换，(7)式可以近似为：

(8)

其中，u0(x,fm;xs)和λ0(x,fm;xs)分别为频率为fm的单炮入射波场和单炮走时残差Δtk(xr,xs)逆时传播得到的单频伴随波场.与(7)式中定义的有限频梯度(走时残差向带限核函数反投影)相类比，本文将走时残差向单频核函数的反投影结果(即公式(8))定义为单频梯度.

利用表2中的单频核函数计算流程，可以高效地计算目标泛函的单频梯度.相比于基于波场重构算法的梯度计算策略，单频梯度计算方法对内存需求极低(相对于波动方程正演，仅需要增加两个单频波场以及一个成像体)，且计算量至少减少1/3(仅需要两次波传播，且无需增加数值吸收边界内的额外计算量)，因此更适用于大规模层析反问题求解.表3定量对比三种不同方法计算梯度所需的波场传播次数和存储空间.以存储入射波场策略为例(其他两种方案不再赘述)：时空域的入射波场需要占用Nt·N3个浮点数，逆时传播的伴随波场占用2N3个浮点数，三维速度模型和梯度占用4N3个浮点数.此外，还需要额外处理数值吸收边界内的辅助波场计算和存储，因此内存需求(浮点数目)≥Nt·N3+4N3.

表3 不同梯度计算策略所需的波场传播次数和存储空间对比(假定三维速度模型在x, y, z方向网格数目均为N，时间方向外推Nt步，正演算子为二阶常密度声波方程)

2 数值实验

首先计算三维均匀速度模型中的单频走时敏感度核函数，并与带限走时敏感度核函数的形态相对比；接着，用三维复杂模型分别计算单频梯度和传统的有限频梯度，通过分析单频梯度的特点，验证单频反演的可行性和有效性.

数值实验1：三维均匀速度模型中的单频核函数形态

首先，如图1所示，设计了一个三维均匀速度模型(速度2500 m·s-1，上表面为自由边界，其余5面为随机边界).速度模型在x,y,z方向网格数目分别为801，801，301，三个方向上网格间距均为5 m.震源子波为主频40 Hz的Ricker子波，正演步长0.5 ms，记录时间1 s.震源和检波器放置在地下一个网格深度处，偏移距为2.5 km.图2为基于广义Rytov近似的有限频走时敏感度核函数(冯波等，2021)，呈现中空的形态(Tromp et al., 2005).图3为单频(40 Hz)走时敏感度核函数.对比图2和图3可以看出，在第一菲涅尔体内，有限频核函数与单频核函数形态相近；而在第一菲涅尔体外，单频核函数呈明显的空间振荡特性，符合理论预测，也与前人研究结果(Marquering et al., 1999)相符合.此外，随机边界产生的散射波场在第一菲涅尔体外出现了一定程度上的相互干涉，减弱了单频核函数的规律振荡.

图1 三维匀速模型(地表自由边界+四周随机边界)

图2 传统的有限频率走时敏感度核函数

图3 单频走时敏感度核函数(40 Hz)

数值实验2：三维Overthrust模型

为验证利用单频梯度进行速度反演的可行性，我们采用三维Overthrust速度模型(图4，Aminzadeh et al., 1994)测试走时目标泛函的梯度形态.原始速度模型在x,y,z方向网格数目分别为801，801，187，三个方向上网格间距均为25 m.宽方位观测系统设计为：一共39条炮线，炮线间隔500 m.每条炮线共39炮，炮间隔500 m.总计1521炮，震源深度为25 m.第一炮坐标(sx,sy)=(0.5, 0.5)km.每炮共242条接收线，每束线242道接收.检波器线间距和道间距均为50 m，深度为25 m.单炮观测系统如图5所示，x,y方向的偏移距范围均为-6025～-25 m以及25～6025 m.最小、最大偏移距分别为25 m和8521 m.所有的检波器坐标(gx,gy)满足：0 m≤gx≤20000 m, 0 m≤gy≤20000 m.

图4 三维Overthrust速度模型

图5 观测系统示意图

采用有限差分方法求解三维标量声波方程正演地震记录.震源子波为主频8 Hz的Ricker子波，正演步长2 ms，记录时间4 s.初始速度模型为线性递增的速度场(v(z=0)=2.5 km·s-1,v(z=4.0 km)=5.0 km·s-1).分别用真实模型和初始模型进行波动方程正演，并用初至自动拾取方法(罗飞和王华忠，2021)拾取初至并计算初至时差.由于对观测数据(真实模型正演记录)和合成数据(初始模型正演记录)采用了相同的拾取算法及参数，消除了走时拾取方法对初至时差测量的误差.

分别计算走时目标函数的有限频梯度和单频梯度.图6a展示了以初始模型为背景速度，第16条炮线中第15炮的单频(8 Hz)梯度(内存占用约为1.73 GB).与单个炮检对的单频核函数(图3)相比，由于多个检波器(每炮共58564道)的核函数相互叠加干涉，第一菲涅尔体外的单频核函数的空间振荡特征明显减弱，而第一菲涅尔体内的能量(图6a中蓝色部分)得到了增强.另一方面，传统的有限频单炮梯度(冯波等，2021)如图6b所示(有限频梯度的频带展布范围同震源子波，内存占用约为7.33 GB).虽然传统的有限频梯度形态更加光滑，但单频梯度仍然较好的呈现了有限频梯度的低波数背景部分.图7a为多炮单频梯度(931炮单频梯度求和).与单炮单频梯度相比，经过多炮叠加后的梯度，低波数成分得到明显增强，而高波数振荡部分进一步相互干涉相消.高波数噪声(梯度假象)主要体现在近震源处的采集脚印.图7b为传统的有限频多炮梯度.与图7a中的单频梯度相比，主要的区别在于有限频多炮梯度的高波数噪声得到了更好的压制.

图6 单炮梯度对比

为了有效消除多炮梯度中的高波数干扰(主要是稀疏震源空间采样引起的采集脚印以及未干涉掉的空间高频振荡噪声)，我们用高斯光滑滤波器P(x,σ)=exp(-x/σ)2与图7中的多炮梯度进行褶积，光滑后的多炮单频梯度和多炮有限频梯度如图8所示.可以看出，高斯光滑滤波器完全消除了多炮梯度中的高波数噪声，单频梯度和有限频梯度的光滑背景部分(即梯度的低波数成分)基本一致，均能用于背景速度反演.因此，用光滑后的多炮单频梯度进行背景速度反演是可行的.

图7 多炮梯度对比

图8 高斯平滑(σx=σy=500 m, σz=100 m)后的多炮梯度对比

为验证单频梯度的反演效果，本文设计了一种多偏移距反演策略(偏移距递减策略)：即最大偏移距分别为8000 m，6000 m，4000 m，2000 m(本文称之为4个尺度).每个尺度内采用梯度导引的优化算法(Métivier and Brossier, 2016)，且每个尺度的初始模型为上一个尺度的反演结果.目标函数梯度采用高斯平滑滤波器提取梯度中的低波数成分，步长用线性搜索算法计算.每个尺度内最多迭代5次，每个尺度内反演停止准则为走时相对误差的L2范数小于预设的阈值(10-1).图9中展示了每个尺度内第一次迭代时目标函数的梯度形态，其中，前两个尺度(最大偏移距为8000 m和6000 m)采用了400炮计算梯度(x和y方向上炮间距1000 m，共计400炮，高斯平滑参数取：σx=σy=500 m,σz=100 m)，后两个尺度(最大偏移距为4000 m和2000 m)采用了所有炮计算梯度(x和y方向上炮间距500 m，共计1521炮，高斯平滑参数取：σx=σy=250 m,σz=50 m).每个尺度内的累计速度扰动(当前尺度的反演结果减去上一尺度的反演结果)分别如图10所示.从泛函梯度以及每个尺度内的累计速度更新量可以看出，随着最大偏移距递减，模型的有效反演深度逐渐变浅，但近地表的分辨率逐渐提升.

图9 目标泛函的梯度形态

图10 每个尺度内的累计速度扰动

最终的反演结果如图11所示，可以看出浅层的反演效果较好，有效反演深度约为1.0 km.为进一步对比反演效果，图12展示了深度分别为150 m，200 m和250 m处的等深度切面，无论是低速异常还是局部高速异常，反演结果上都有较好的体现.在150 m和200 m深度切面上，小尺度的速度异常体(红色箭头)也能很好的刻画出来.随着深度的增加，反演结果的分辨率逐渐降低，但仍然可以较好的刻画速度模型的背景变化趋势.图13是等深度面沿x方向在不同y坐标处的抽线对比，在浅层的地方趋势上已经与真实的速度模型接近，局部甚至优于平滑后的真实模型.

图11 最终反演结果与真实结果对比

图12 真实模型和反演结果的等深度切面对比

图13 等深度(200 m)切面沿x方向抽线对比

3 讨论

本文通过两个数值实验，初步验证了基于随机边界条件的单频梯度在计算效率、存储需求方面的优势以及用于三维走时反演的可行性.

首先，随机边界条件的引入降低了入射场和伴随场在随机边界内的空间相关性，且地震观测系统的多次覆盖特性有助于敏感度核函数在第一菲涅尔带外的空间振荡部分相互干涉.单频梯度中残余的高波数假象可以用光滑算子(如三维高斯函数)有效消除，保证了单频梯度更新方向的稳健性，进而有助于实现背景速度反演.

其次，在数值计算方面，本文提出的单频走时敏感度核函数计算方法仅需要一次正传播+一次反传播，而基于波场重构策略的传统梯度计算方法需要一次正传播+两次反传播.此外，波场重构策略中的三次波场传播均需要施加吸收边界条件，引入了额外的计算量及内存消耗.相比之下，本文方法无需采用数值吸收边界且仅需要两次波场模拟，因而计算量严格小于传统方法的2/3.同时，由于随机边界的引入使得正演算法进一步简化(即算法各向同性)，因而更适用于GPU加速.在存储方面，基于波场重构策略的传统梯度计算方法需要存储6个面上的边界波场，存储量仍然非常可观，而随机边界方法仅需存储(复值)单频入射场和伴随场，内存需求大幅降低.因此，本文方法的计算效率高且内存需求低，能够满足三维反演的计算需求.

再次，在反演分辨率方面，均匀模型中单频波的波长决定了能够反演的异常体尺度.对于大尺度背景速度模型反演，可以采用更稀疏的炮间距(如本文第二个数值实验中前两个尺度的反演，仅用了400炮(炮间距为1.0 km)就实现了大尺度的速度更新)从而减低总体计算量.另一方面，为提高近地表反演精度，可以在单频反演结果的基础之上，仅用小偏移距数据作为传统的有限频走时层析的输入，进一步更新近地表速度模型.

最后，在三维波动方程反演实用化方面，计算效率是关键问题.本文仅采用了最基础的最速下降方法用于展示三维单频走时反演的可行性.由于忽略了Hessian逆算子对于振幅均衡的影响，导致有效反演深度不超过1 km.若采用二阶反演算法，如作者提出的基于隐式矩阵-向量乘的高斯-牛顿方法(冯波等，2019)，有望进一步提高深层的反演质量并加速反演收敛.

4 结论

针对波动方程走时反演方法在求解大规模三维问题时遇到的计算和存储瓶颈，提出了基于随机边界的单频梯度高效反演方法.通过将随机边界条件引入三维单频敏感度核函数的构造，有效降低了波传播的算法复杂度(无需单独处理边界)，并大幅降低了单频梯度计算过程对内存的需求.三维复杂模型反演结果验证了本文方法的有效性和可行性.

后续的研究中可以引入随机梯度的思想，将炮点随机采样与多偏移距反演策略结合，进一步降低反演的总体计算量.

致谢感谢中国石油天然气股份有限公司勘探开发研究院和西北分院、中海油研究总院有限责任公司和中海石油(中国)有限公司湛江分公司、中国石油化工股份有限公司石油物探技术研究院和胜利油田分公司对波现象与智能反演成像研究组(WPI)研究工作的资助与支持.