You J C, Cao J X, Wang J. 2020. Two-way wave equation prestack depth migration using the matrix decomposition theory. Chinese J. Geophys . (in Chinese), 63(10): 3838-3848, doi: 10.6038/cjg2020O0027 . 2020-08-14 收修定稿 第一作者简介: 尤加春, 男, 1988年生, 现成都理工大学从事博士后研究工作, 主要从事深度偏移成像研究 摘要 ：叠前深度偏移理论及方法一直是地震数据成像中研究的热点问题.业界对单程波叠前深度偏移方法和逆时深度偏移开展了深入的研究，但对双程波方程波场深度延拓理论及成像方法的研究还鲜有报道.本文以地表记录的波场值为基础，利用单程波传播算子估计波场对深度的偏导数，为在深度域求解双程波方程提供充分的边界条件，并提出利用矩阵分解理论实现双程波方程的波场深度外推.通过对强速度变化介质中传播波场的计算，与传统的单程波偏移方法相比，本文提出的偏移方法计算的波场与常规有限差分技术计算的波场相一致，证明了本方法计算的准确性.通过对SEAM模型的成像，在相同的成像参数下，与传统的单程波偏移算法和逆时深度偏移算法方法相比，本文提出的偏移方法能够提供更少的虚假成像和更清晰的成像结果.本文所提偏移算法具有深度偏移和双程波偏移的双重特色，推动和发展了双程波叠前深度偏移的理论和实践. 地震偏移通常是指利用特定的数学模型(通常使用声波或弹性波方程作为数学模型)构建合理的波场外推方程，并应用建立的外推方程计算在已知宏观速度场中从地表观测到的波场到地下任意深度的波场(波场反向传播)，然后引入成像条件计算成像振幅( Berkhout, 1982 ; Gray et al., 2001 ; Stolt and Weglein, 2012 ).地震偏移的本质是使用数学手段沿着波传播路径反向传播地表观测数据，并消除能量损失或相位畸变等的传播效应，以获得地下结构的准确图像.从地震偏移的基本描述可以概括出偏移成像涉及到的两个关键问题：波场外推及成像条件.本文重点讨论波场的外推方式和方法. 对于波场深度延拓方法，学者更多讨论的是基于单程波方程的深度偏移方法，如已开展深入研究的傅里叶有限差分(Fourier finite difference, FFD)偏移方法和广义屏(Generalized Screen Propagator, GSP)偏移方法( Gazdag and Sguazzero, 1984 ; Ristow and Rühl, 1994 ; Wu, 1994 ; Le Rousseau and De Hoop, 2001 ).单程波方程是通过对双程波方程进行上行波和下行波分解得到的.传统的单程波方程在处理上行波和下行波分解时，舍弃了振幅修正项，因此不具备保幅特性. Zhang等(2003 , 2005 )从理论上推导了具有振幅校正项的单程波真振幅方程.此外，传统的单程波偏移方法往往利用泰勒级数或Pade级数等近似理论计算垂直波数，导致对强横向速度变化介质和复杂介质的成像角度有限，特别是对陡倾角介质的成像.为了提高单程波偏移方法在复杂介质中的成像能力， Grimbergen等(1998) 采用波场模态分解的方法计算单程波传播算子. Wu等(2000) 提出了局部分解的Beamlet传播算子. Zhang等(2010) 采用优化算法优化傅里叶有限差分展开系数提高傅里叶有限差分偏移方法的成像角度. Zhang和Liu(2006 , 2007 )提出采用优化的裂步傅里叶偏移算法实现更加准确的三维深度偏移. You等(2018a) 发展了基于矩阵分解的单程波真振幅偏移算法，提高了单程波真振幅偏移算法在复杂介质中的成像精度. 逆时深度偏移(Reverse time Migration, RTM)方法作为一种基于双程波方程的偏移方法，由于其能描述各种波的传播规律，受到了工业界和学术界的广泛研究( Baysal et al., 1983 ； 李博等，2010 ； 刘红伟等，2010 ).为了准确地对复杂介质进行成像，利用双程波方程或全声波方程开展偏移成像是实现这一目的理论基础，特别是逆时深度偏移方法对复杂介质的成像效果更说明了其优势.显然，单程波偏移方法在成像精度方面并不能满足这一要求.目前，关于双程波方程波场深度延拓及成像方法的理论和实践报道还较少.因此，对利用双程波方程开展波场深度外推的理论和实践进行研究具有理论意义.众所周知，标量声波方程是关于时间和空间变量的二阶偏微分方程.理论上，需要两个边界条件才能在深度域和时间域中准确求解.传统的地震数据采集系统只记录地表的波场，通常只能提供一个边界条件.因此，直接利用常规地震数据采集系统求解双程波方程缺乏充分的边界条件.而对于双程波波场深度延拓及成像方面，一些学者开展了相关研究工作. Kosloff和Baysal(1983) 建立了一个包含Helmhotz算子的一阶方程组，采用经典的四阶Runge-Kutta方法求解该方程组.但对于波场深度延拓过程中的倏逝波问题，Kosloff和Baysal在频率-波数域由每层最大速度构建了一个低通滤波器，达到消除倏逝波的目的； Sandberg和Beylkin(2009) 根据倏逝波和有效波所对应特征值的差异，提出采用谱投影算法解决双程波方程波场深度延拓过程中的倏逝波问题. Wu等(2012) 提出采用单程波Beamlet传播算子实现双程波波场深度延拓及成像，但该偏移算法在强速度变化介质中存在不稳定. Pan(2015) 提出采用类似于单程波偏移算法中相移加插值的数学原理在深度域实现双程波波场延拓.为了解决双程波方程深度延拓的边界条件问题， You等(2016) 提出在不同深度处设置上下两层排列的检波器记录两个深度的波场值，研究了双程波深度偏移算法在保幅计算上的优势.基于双检波器地震数据采集系统， You等(2018b) 提出采用单程波传播算子实现双程波深度偏移成像的思想，数值实验对比了该偏移算法与常规偏移算法的成像性能差异.在执行原理上，时间域逆时深度偏移算法是在时间域上求解声波方程，即在时间轴上实现震源波场的正向计算和检波点波场的反向逆时外推；而利用双程波方程开展波场深度延拓的偏移算法是在深度域求解声波方程，即在深度方向上实现震源波场和检波点波场的延拓计算. 本文以双程波方程波场深度延拓及成像为研究内容，首先为了解决利用声波方程开展波场深度延拓边界条件不足的问题，本文根据地表处记录的波场值使用单程波偏移方法估计波场对深度的偏导数，提出利用矩阵分解理论开展双程波方程的波场深度外推.在数值实验中，本文利用常规的单程波广义屏传播算子、本文提出的偏移方法和有限差分技术计算波场在一个倾斜界面速度模型和盐丘模型中的传播规律，对比三种偏移算法在波场计算上的差异；此外，利用常规的广义屏偏移算法、本文提出的偏移算法和逆时深度偏移算法对SEAM模型进行成像处理，以评估不同算法之间的成像性能差异. 1 基于双程波动方程的波场深度外推 在二维介质中，频率域声波方程可以表示为 在方程(4)中，可以发现新的波场延拓方程包含有两个单程波传播算子e ^{i k _z Δ z} 和e ^{- i k _z Δ z} 用以根据当前深度的波场值及波场偏导数值计算下一个深度的波场值及波场偏导数值.由于两个单程波传播算子的存在，方程(4)实现了波场的正向传播和反向传播，并将正向传播和反向传播的波场求和得到下一个深度的总波场，具有波场双程传播的特点，因此本文所提的波场深度延拓算法属于双程波波场延拓范畴. 将欧拉公式应用于方程(4)，可得 其中垂直波数 使用方程(2)中提供的两个边界条件，本文实现了对波场及其偏导数从深度 z 到深度 z +Δ z 的延拓计算，如何精确求解垂直波数和三角函数是解决双程波方程波场延拓问题的核心. 对震源波场和检波点波场进行深度延拓计算之后，本文采用常规的互相关成像条件计算每一层的成像振幅. 本文从成像条件理论上分析深度偏移方法和逆时深度偏移方法的差异.逆时深度偏移方法是在时间域内求解声波方程，在全速度模型空间内开展波场计算. Liu等(2011) 提出将检波点波场和震源波场分别分解为上行波和下行波，波场分解表达为：

其中， s _d ( x , t )和 s _u ( x , t )分别表示震源波场的上行波场和上行波场； r _d ( x , t )和 r _u ( x , t )分别表示检波点波场的下行波场和上行波场. s ( x , t )和 r ( x , t )分别表示总震源波场和检波点波场 经过波场分解之后的互相关成像条件写为 其中， 计算的是震源下行波场与检波点上行波场的互相关成像结果， 计算的是震源上行波与检波点下行波场的互相关成像结果； 是震源波场和检波点波场的同向传播波场的互相关成像结果，为逆时深度偏移成像中的低频噪声. 在深度偏移算法中，双程波方程波场深度延拓方法是在深度域求解声波方程，将波场传播的能量主要集中在深度传播方向上.理论上，双程波深度偏移的互相关成像条件可以写为 对比成像条件(8)和(9)不难发现，双程波深度偏移的互相关成像条件与单程波偏移方法使用的成像条件在结构上一致，从成像条件上分析，深度偏移比逆时偏移产生更少的成像噪声.这种成像差异主要体现在两种偏移算法执行原理上. 2 利用矩阵分解理论实现波场深度外推 在本部分中，本文详细分析如何准确地计算方程(6)中垂直波数及其三角函数实现波场的深度外推.假设速度模型为各向同性或水平层状介质，则可以采用相移法(Phase Shift, PS)进行数值求解( Gazdag, 1978 ).然而，实际情况很难满足上述假设.一般而言，速度模型在水平和垂直方向上均有变化，在某些情况下(如盐模型)，这种速度变化还非常剧烈. 为了解决波场准确延拓的问题，本文以Helmhotz算子为出发点，将在任意速度变化介质中的Helmhotz算子(采用符号 L 表示)以矩阵形式离散表示为 在矩阵分解方程(11)中，对角化矩阵 Λ 包含有负特征值和正特征值.由于负特征值的平方根产生虚数，对应的是波场传播中的倏逝波，倏逝波在波场深度延拓过程中会使振幅呈现指数增长，导致波场深度延拓的不稳定，因此在波场深度延拓过程中必须滤掉负特征值. 3 数值实例 在数值实验中，本文设计了一系列的理论模型并对比分析了本文提出的双程波深偏移算法与传统的单程广义屏算法和逆时深度偏移算法在成像性能上的差异.在本数值实验中使用的有限差分方法和逆时深度偏移算法均采用的是10阶交错网格有限差分格式. 3.1 倾斜界面中的脉冲响应 本实验设计了倾角为30°的倾斜界面速度模型.界面两侧速度分别为2000 m·s ^-1 和4000 m·s ^-1 ，倾斜界面两侧表现出较强的横向速度变化，速度模型如 图 1 所示，其垂直方向和水平方向的网格间隔分别为20 m.计算脉冲响应的震源位于 x =4000 m, z =0 m处，子波采用主频率为10 Hz的雷克子波，时间采样率为0.001 s，总计算时长为4.0 s.本文使用有限差分算法、单程波广义屏偏移算法和本文所提偏移算法计算波场传播，并记录 t =2.5 s时刻的波场，如 图 2 所示.虽然本文中所提偏移方法是以双程波波动方程为基础，但是在波场外推时，波场深度延拓算法是根据 z 深度处的波场及其偏导数波场计算 z +Δ z 深度处的波场及其偏导数波场，只记录了沿深度增加方向的波场及其偏导数值，因此只观测到下行波场.本文以有限差分法计算的正演模拟波场作为参考标准.当波通过倾斜界面时，由于界面两侧的速度差异，入射波和透射波发生分离，如 图 2a 中的红色箭头所示.本文提出的偏移方法计算的外推下行波与有限差分法计算的下行波结果相同.利用单程波广义屏算法计算的外推下行波场中，入射波和透射波在界面处不发生分离，而是黏合在一起的，如 图 2b 中的红色箭头所示.显然，利用常规的广义屏偏移方法计算的外推波场与理论分析存在明显差异.从刻画波场传播的精度上对比，本文所提偏移算法计算的准确波场为地震深度成像奠定了基础. 除了对比三种算法计算的波场之外，本文统计了三种偏移方法完成整个波场计算所需时间，见 表 1 .为了尽可能地覆盖有效振幅，频率域偏移方法(包括单程波广义屏方法和本文所提偏移方法)的计算频带为1~20 Hz.在本算例中，CPU并行测试平台搭载的是56核心Intel Xeon 8180 CPU@ 2.5G Hz.从 表 1 中可见，常规单程波偏移算法所需的计算时间最少，CPU并行版本的有限差分方法次之，本文所提偏移算法用时最多.这是由于本文所涉及的偏移算法涉及到大量的矩阵分解及矩阵乘法运算，常规的频率域偏移算法一般是逐个频率依次计算，这样极大的降低了算法的计算效率.因此，采用频率并行的方式可在一定程度上提高计算效率，本文所提偏移算法的CPU并行程序比无并行的程序计算效率 ( Table 1 ) 3.2 强速度变化介质中波场计算 对强速度变化介质中波现象的描述一直是地震速度建模和成像研究的重点.准确描述波场的传播是地震成像和反演的基础，因此有必要对波场计算进行深入研究，特别是对复杂介质中波场的刻画显得尤为重要.在数值实验中，本文使用有限差分方法、传统的单程波广义屏方法和本文所提偏移算法计算盐丘模型中波场传播.在波场计算中，采用了两个速度模型; 一个速度模型是原始速度模型，如 图 3a 所示；另一种速度模型是使用5×5窗口的高斯滤波器生成的平滑速度模型，如 图 3b 所示，使用平滑速度模型的目的是测试不同方法在速度不准确的情况下计算波场的传播情况.在 图 3 中，盐丘体速度(4500 m·s ^-1 )大概是围岩速度(平均2200 m·s ^-1 )的两倍，速度模型表现出强速度变化. Fig. 4 Wavefield snapshots calculated by using different methods based on the accurate velocity model at time t =2.0 s. (a) Conventional finite difference method; (b) One-way GSP method; (c) Our proposed two-way wave wavefield depth extrapolation scheme based on matrix decomposition. 在精确速度模型的基础上，本文比较了三种偏移算法计算的波场差异.当波通过盐丘时，单程波广义屏方法计算的外推下行波波场与有限差分法计算的下行波波场相比产生较大的误差，如 图 4a 和 b 中箭头所示，单程波广义屏方法并未准确地计算波场的传播，这将不利于盐丘体下方构造的成像；而本文提出的偏移算法计算的外推下行波与有限差分方法正演模拟计算的下行波较一致，如 图 4c 中箭头所示.通过与有限差分法和单程波偏移算法计算的波场相比较，本文提出的偏移方法在强速度变化介质中表现出稳健的波场计算性能. 全波形反演中常常使用的是一个相对平滑的速度模型作为反演的初始模型，利用偏移算法开展波场在平滑速度模型中传播规律的研究对偏移成像和全波形反演也至关重要，可为使用本文所提双程波波场深度延拓技术开展全波形反演研究做理论准备.在平滑速度模型中，本文采用上述三种偏移算法计算了 t =2.0 s处的波场，如 图 5 所示.由于采用了平滑的速度模型，平滑速度模型中的断层与围岩的速度差异变小，使得有限差分法正演模拟计算的波场中反射波几乎不可见.同时，在平滑速度模型中，本文提出的偏移方法计算的外推下行波波场与有限差分方法正演模拟计算的下行波波场仍然保持高度一致，而单程波广义屏偏移算法计算的外推下行波波场和理论解释有一定差异，特别是 图 5 中箭头所指示的波场.此外，与利用准确的速度模型计算的波场相对比，采用平滑的速度模型，本文所提的偏移算法计算的外推波场似乎含有更少的干扰噪声. Fig. 5 Wavefield snapshots calculated by using different methods based on the smoothed velocity model at time t =2.0 s (a) Conventional finite difference method; (b) One-way GSP method; (c) Our proposed two-way wave wavefield depth extrapolation scheme based on matrix decomposition. 3.3 SEAM模型的叠前深度偏移 SEAM模型包含一个典型的盐丘体，该模型具有精细的层位信息和充满流体的储层，因此对该模型开展偏移成像研究对于油气资源勘探具有重要意义( Fehler and Keliher, 2011 ).由于盐丘体与周围沉积层的速度差异较大，盐丘体对地震波的散射作用较强，因此对盐丘体及其构造的准确成像一直是地震偏移研究的难点问题.本文采用单程波广义屏方法、常规逆时深度偏移方法和本文提出的深度偏移方法对SEAM模型正演的炮集进行偏移成像处理.在数值实验中，本文利用了两种速度模型，一种是用于生成炮集记录的原始速度模型；另一种是利用5×5的高斯窗口平滑的速度模型，用于偏移成像，速度模型如 图 6 所示.速度模型在水平方向有24 km，深度方向有15 km.SEAM模型在水平方向和垂直方向的网格间距分别为20 m和40 m.采用有限差分技术模拟计算了150炮的地震记录，炮间距为160 m.在每一炮的计算中，采集器均匀分布于 x =0 m到 x =24000 m的范围内，检波器间距20 m. 从单程波广义屏偏移方法计算的成像结果可以看出，单程波广义屏偏移方法对盐丘体的整体成像比较模糊，特别是对盐丘体底界面的成像几乎不可见，如 图 7a 中箭头所示，并且该偏移结果中存在明显的成像噪声.从整体构造成像效果来看，基于双程波方程的偏移方法(包括逆时深度偏移方法)比单程波偏移方法在成像复杂模型或速度变化剧烈介质时具有明显的优势.对比逆时深度偏移和本文所提偏移算法计算的成像结果，不难发现在相同的参数下，逆时深度偏移算法计算的成像结果出现了较明显的频散噪声，而利用本论文所提的偏移算法计算的成像结果更加稳定、清晰.理论上，可以通过使用更小的网格间隔来压制逆时深度偏移算法中产生的频散，但这将带来更大的计算成本并需要更多的存储空间.为了进行详细的对分析，本文局部放大显示了 图 7 中虚线框内的成像剖面，如 图 8 所示. 图 8a 为速度模型中的层位结构，对比逆时深度偏移算法和本文所提偏移算法计算的局部成像结果( 图 8 ( b ， c ))，可见本文所提偏移算法提供了更加清晰的层位位置和层位接触关系.由于逆时深度偏移算法和本文提出的偏移算法都是基于双程波动方程的，在构造成像中很难看出两种方法的明显区别.本实验的目的是为了证明本文所提的双程波方程波场深度延拓的偏移算法可以和逆时深度偏移算法一样，对复杂介质进行准确的成像，为双程波偏移方法提供另一种选择，而不是代替逆时深度偏移算法.此外，由于本文提出的偏移方法仅在深度方向进行波场外推，因此几乎不对某些特殊结构(例如 图 7c 中椭圆)进行成像，但这种垂直结构可以通过回转波进行成像( Hale et al.1992 ; Jia and Wu, 2009 ).逆时深度偏移方法能在整个速度模型空间范围内对包括回转波在内的各种波场进行正向和反向传播，这使垂直结构能部分被成像，即成像条件计算公式(8)中 I ₁ + I ₂ 项，如 图 7b 中的椭圆所示.这也是利用双程波方程开展波场深度延拓及成像需要继续研究的方向之一. 本文以在深度域求解双程波方程所需的边界条件入手，利用单程波偏移算法估计波场对深度的偏导数.基于地表记录的波场及其估计的偏导数波场，建立了含有垂直波数及其三角函数的双程波波场深度外推方程.根据Helmhotz算子与垂直波数之间的数学联系，本文提出使用矩阵分解理论实现双程波方程的波场深度外推.通过利用单程波广义屏方法、有限差分法和本文所提的方法对倾斜界面速度模型和盐丘模型中传播波场的计算，证明了本文所提方法与有限差分方法计算的下行波场精度相同，而常规单程波广义屏方法在强变化介质中计算的外推波场与理论分析存在一定误差.利用常规的逆时深度偏移算法、单程波广义屏方法和本文所提的偏移方法对SEAM模型开展叠前深度偏移成像研究，与单程波广义屏算法计算的结果相比，双程波深度偏移算法在复杂介质成像方面具有很强的优势.此外，在相同的参数下，本文提出的偏移方法比逆时深度偏移方法获得了更清晰的成像结果，频散噪声更小，但逆时深度偏移方法比本文所提深度偏移方法具备一些特殊特征，例如能实现回转波成像.总而言之，本文所提深度偏移方法相较于单程波偏移方法具备更好成像复杂构造的能力，相较于逆时深度偏移算法具有自己的特色和不足，而这些不足还有待进一步的研究. Baysal E, Dan Kosloff D, Sherwood J W C. 1983. Reverse time migration. Geophysics , 48(11): 1514-1524. DOI:10.1190/1.1441434 Berkhout A J. 1982. Seismic Migration:Imaging of Acoustic Energy By Wave Field Extrapolation. 2nd ed. Amsterdam: Elsevier. Fehler M, Keliher P J. 2011. SEAM Phase I:Challenges of Subsalt Imaging in Tertiary Basins, with Emphasis on Deepwater Gulf of Mexico. Tulsa: Society of Exploration Geophysicists. Liu L N, Zhang J F. 2006. 3D wavefield extrapolation with optimum split-step Fourier method. Geophysics , 71(3): T95-T108. Pan J H. 2015. A new upward and downward wavefield continuation for two-way wave equation migration.//85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 4180-4184. Sandberg K, Beylkin G. 2009. Full wave equation depth extrapolation for migration. Geophysics , 74(6): WCA121-WCA128. Stolt H R, Weglein A B. 2012. Seismic Imaging and Inversion:Volume 1:Application of Linear Inverse Theory. New York: Cambridge University Press. Strang G. 2006. Linear Algebra and Its Application. 4th ed. Boston,MA: Cengage Learning. Wu R S. 1994. Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method. Journal of Geophysical Research:Solid Earth , 99(B1): 751-766. DOI:10.1029/93JB02518 Wu R S, Wang Y Z, Gao J H. 2000. Beamlet migration based on local perturbation theory.//70th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1008-1011. Zhang J H, Wang W M, Wang S Q, et al. 2010. Optimized Chebyshev Fourier migration:a wide-angle dual-domain method for media with strong velocity contrasts. Geophysics , 75(2): S23-S34. Zhang X D. 2015. Matrix Analysis and Applications (in Chinese). 2nd ed. Beijing: Tsinghua University Press. 刘红伟, 刘洪, 邹振, 等. 2010. 地震叠前逆时偏移中的去噪与存储. 地球物理学报, 53(9): 2171-2180. DOI:10.3969/j.issn.0001-5733.2010.09.017 张贤达. 2015. 矩阵分析与应用. 2版. 北京: 清华大学出版社.