地球物理学报  2020, Vol. 63 Issue (10): 3838-3848   PDF    
利用矩阵分解理论的双程波方程叠前深度偏移方法
尤加春, 曹俊兴, 王俊     
成都理工大学地球物理学院, 成都 610059
摘要:叠前深度偏移理论及方法一直是地震数据成像中研究的热点问题.业界对单程波叠前深度偏移方法和逆时深度偏移开展了深入的研究,但对双程波方程波场深度延拓理论及成像方法的研究还鲜有报道.本文以地表记录的波场值为基础,利用单程波传播算子估计波场对深度的偏导数,为在深度域求解双程波方程提供充分的边界条件,并提出利用矩阵分解理论实现双程波方程的波场深度外推.通过对强速度变化介质中传播波场的计算,与传统的单程波偏移方法相比,本文提出的偏移方法计算的波场与常规有限差分技术计算的波场相一致,证明了本方法计算的准确性.通过对SEAM模型的成像,在相同的成像参数下,与传统的单程波偏移算法和逆时深度偏移算法方法相比,本文提出的偏移方法能够提供更少的虚假成像和更清晰的成像结果.本文所提偏移算法具有深度偏移和双程波偏移的双重特色,推动和发展了双程波叠前深度偏移的理论和实践.
关键词: 双程波方程      边界条件      波场深度外推      矩阵分解     
Two-way wave equation prestack depth migration using the matrix decomposition theory
YOU JiaChun, CAO JunXing, WANG Jun     
School of Geophysics, Chengdu University of Technology, Chengdu 610059, China
Abstract: Study on the theory and method of prestack depth migration has always been a focused topic in seismic data imaging. The one-way wave equation depth migration and reverse time migration methods have been deeply studied, while less research on wavefield depth extrapolation and imaging based on two-way wave equation is reported. Based on the wavefield recorded at the surface, we use a one-way wave propagator to estimate the derivative wavefield, which provides sufficient boundary conditions for solving two-way wave equation in the depth domain and propose a matrix decomposition scheme to perform wavefield depth extrapolation based on two-way wave equation. Through the calculations of wavefields in the medium with strong velocity changes, it is proved that the down-going waves obtained by the proposed method are consistent with that achieved by the conventional finite difference technique. Compared with the one-way GSP method, this numerical experiment verifies the accuracy and stability of the proposed two-way wave depth migration algorithm. In imaging the SEAM model, we can see that the proposed imaging method is able to provide a clear result with less artifacts compared with the results produced by the conventional GSP and RTM methods. Our proposed method has the characteristics of depth migration and two-way wave equation migration, which is of great significance in theory and practice.
Keywords: Two-way wave equation    Boundary condition    Depth extrapolation    Matrix decomposition    
0 引言

地震偏移通常是指利用特定的数学模型(通常使用声波或弹性波方程作为数学模型)构建合理的波场外推方程,并应用建立的外推方程计算在已知宏观速度场中从地表观测到的波场到地下任意深度的波场(波场反向传播),然后引入成像条件计算成像振幅(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 基于双程波动方程的波场深度外推

在二维介质中,频率域声波方程可以表示为

(1)

其中,xz分别表示水平坐标轴和垂直坐标轴.v(x, z)是介质的速度,是频率域的波场表示,ω是角频率.

从方程(1)可以看到, 频率域声波方程是关于空间变量的二阶偏微分方程.从理论上讲, 需要提供两个边界条件方能求解二阶偏微分声波方程, 即地表处记录的波场值及其偏导数,频率域声波方程及其边界条件可以写为

(2)

其中φ1, φ2为已知值,为在地表处记录的波场及其偏导数值.

在均匀介质中,方程(2)有其通解形式,通解可表示为

(3)

其中,C1, C2为两个待定系数,是水平波数,c是均匀介质中常速度.

在传统的地震数据采集系统中,由于只能提供一个边界条件,因此理论上只能求解一个待定系数,这样就形成了单程波传播算子.在本文的研究中,利用方程(2)中的两个边界条件可以确定两个系数C1, C2,则根据地表处(z=0 m)记录的波场及其偏导数计算下一个深度(zzm)的波场及其偏导数的深度延拓公式为

(4)

在方程(4)中,可以发现新的波场延拓方程包含有两个单程波传播算子eikzΔz和e-ikzΔz用以根据当前深度的波场值及波场偏导数值计算下一个深度的波场值及波场偏导数值.由于两个单程波传播算子的存在,方程(4)实现了波场的正向传播和反向传播,并将正向传播和反向传播的波场求和得到下一个深度的总波场,具有波场双程传播的特点,因此本文所提的波场深度延拓算法属于双程波波场延拓范畴.

将欧拉公式应用于方程(4),可得

(5)

在方程(5)中,本文构建的双程波波场延拓算法涉及到垂直波数及其三角函数的计算.方程(5)中的波场深度延拓公式仍然是在速度为常数的假设前提下推导建立的.为了实现在非均匀介质中波场的深度延拓计算,现将方程(5)推广到非均匀介质的波场延拓,可表示为:

(6)

其中垂直波数使用方程(2)中提供的两个边界条件,本文实现了对波场及其偏导数从深度z到深度zz的延拓计算,如何精确求解垂直波数和三角函数是解决双程波方程波场延拓问题的核心.

对震源波场和检波点波场进行深度延拓计算之后,本文采用常规的互相关成像条件计算每一层的成像振幅.

本文从成像条件理论上分析深度偏移方法和逆时深度偏移方法的差异.逆时深度偏移方法是在时间域内求解声波方程,在全速度模型空间内开展波场计算.Liu等(2011)提出将检波点波场和震源波场分别分解为上行波和下行波,波场分解表达为:

(7)

其中,sd(x, t)和su(x, t)分别表示震源波场的上行波场和上行波场;rd(x, t)和ru(x, t)分别表示检波点波场的下行波场和上行波场.s(x, t)和r(x, t)分别表示总震源波场和检波点波场

经过波场分解之后的互相关成像条件写为

(8)

其中,计算的是震源下行波场与检波点上行波场的互相关成像结果,计算的是震源上行波与检波点下行波场的互相关成像结果;是震源波场和检波点波场的同向传播波场的互相关成像结果,为逆时深度偏移成像中的低频噪声.

在深度偏移算法中,双程波方程波场深度延拓方法是在深度域求解声波方程,将波场传播的能量主要集中在深度传播方向上.理论上,双程波深度偏移的互相关成像条件可以写为

(9)

对比成像条件(8)和(9)不难发现,双程波深度偏移的互相关成像条件与单程波偏移方法使用的成像条件在结构上一致,从成像条件上分析,深度偏移比逆时偏移产生更少的成像噪声.这种成像差异主要体现在两种偏移算法执行原理上.

2 利用矩阵分解理论实现波场深度外推

在本部分中,本文详细分析如何准确地计算方程(6)中垂直波数及其三角函数实现波场的深度外推.假设速度模型为各向同性或水平层状介质,则可以采用相移法(Phase Shift, PS)进行数值求解(Gazdag, 1978).然而,实际情况很难满足上述假设.一般而言,速度模型在水平和垂直方向上均有变化,在某些情况下(如盐模型),这种速度变化还非常剧烈.

为了解决波场准确延拓的问题,本文以Helmhotz算子为出发点,将在任意速度变化介质中的Helmhotz算子(采用符号L表示)以矩阵形式离散表示为

(10)

其中:

L算子的离散化中,本文采用二阶有限差分算子.根据有限差分数值模拟的经验,采用高阶有限差分算子可以获得较高的精度和更稳健的频散压制.

显然,L算子是一个对称矩阵,根据矩阵分解理论,L算子可以进行特征值-特征向量分解(Strang, 2006张贤达,2015),表示为

(11)

其中对角矩阵Λ包含L算子的特征值,Q表示其特征向量.上标T表示转置.进一步,根据矩阵分解理论,在矩阵函数特征值-特征向量的分解中,矩阵函数的特征值为原矩阵特质值的函数,其特征向量保持不变,则包含L算子函数的特征值及其特征向量可以表示为

(12)

因此,方程(6)中kz及其三角函数可以写成:

(13)

(14)

(15)

(16)

在矩阵分解方程(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深度处的波场及其偏导数波场计算zz深度处的波场及其偏导数波场,只记录了沿深度增加方向的波场及其偏导数值,因此只观测到下行波场.本文以有限差分法计算的正演模拟波场作为参考标准.当波通过倾斜界面时,由于界面两侧的速度差异,入射波和透射波发生分离,如图 2a中的红色箭头所示.本文提出的偏移方法计算的外推下行波与有限差分法计算的下行波结果相同.利用单程波广义屏算法计算的外推下行波场中,入射波和透射波在界面处不发生分离,而是黏合在一起的,如图 2b中的红色箭头所示.显然,利用常规的广义屏偏移方法计算的外推波场与理论分析存在明显差异.从刻画波场传播的精度上对比,本文所提偏移算法计算的准确波场为地震深度成像奠定了基础.

图 1 倾斜界面速度模型(界面倾斜角度30°) Fig. 1 A velocity model with a dipping interface (The angle of the interface is 30 degrees)
图 2 使用不同偏移方法计算在t=2.5 s时刻的波场 (a)常规的有限差分法;(b)单程波广义屏方法;(c)本文提出的基于矩阵分解的双程波波场深度外推方法. Fig. 2 Impulse responses achieved by using different methods at t=2.5 s (a) The conventional finite difference method; (b) one-way GSP method; (c) our proposed two-way wave wavefield depth extrapolation scheme based on matrix decomposition.

除了对比三种算法计算的波场之外,本文统计了三种偏移方法完成整个波场计算所需时间,见表 1.为了尽可能地覆盖有效振幅,频率域偏移方法(包括单程波广义屏方法和本文所提偏移方法)的计算频带为1~20 Hz.在本算例中,CPU并行测试平台搭载的是56核心Intel Xeon 8180 CPU@ 2.5G Hz.从表 1中可见,常规单程波偏移算法所需的计算时间最少,CPU并行版本的有限差分方法次之,本文所提偏移算法用时最多.这是由于本文所涉及的偏移算法涉及到大量的矩阵分解及矩阵乘法运算,常规的频率域偏移算法一般是逐个频率依次计算,这样极大的降低了算法的计算效率.因此,采用频率并行的方式可在一定程度上提高计算效率,本文所提偏移算法的CPU并行程序比无并行的程序计算效率

表 1 三种偏移算法计算时间统计 Table 1 The running time of three migration methods

提升了大约3倍.而且矩阵乘法和矩阵分解运算都具备良好的并行化的操作潜力,也可进一步提升算法的计算效率.

3.2 强速度变化介质中波场计算

对强速度变化介质中波现象的描述一直是地震速度建模和成像研究的重点.准确描述波场的传播是地震成像和反演的基础,因此有必要对波场计算进行深入研究,特别是对复杂介质中波场的刻画显得尤为重要.在数值实验中,本文使用有限差分方法、传统的单程波广义屏方法和本文所提偏移算法计算盐丘模型中波场传播.在波场计算中,采用了两个速度模型; 一个速度模型是原始速度模型,如图 3a所示;另一种速度模型是使用5×5窗口的高斯滤波器生成的平滑速度模型,如图 3b所示,使用平滑速度模型的目的是测试不同方法在速度不准确的情况下计算波场的传播情况.在图 3中,盐丘体速度(4500 m·s-1)大概是围岩速度(平均2200 m·s-1)的两倍,速度模型表现出强速度变化.

图 3 速度模型 (a)原始速度模型; (b)平滑速度模型. Fig. 3 The velocity models (a) Original velocity model; (b) Smoothed velocity model.

在本试验中,在z=0 m和x=6800 m的位置放置了一个10 Hz雷克子波作为激发震源.速度模型在水平方向和垂直方向的网格间距均为20.0 m.本文使用上述三种方法计算波场在盐丘模型中的传播,并在准确速度模型中记录了t=2.0 s时刻的波场,结果如图 4所示.

图 4 基于准确速度模型使用不同偏移方法计算t=2.0 s时刻的波场 (a)常规的有限差分法;(b)单程波广义屏方法;(c)基于矩阵分解的双程波深度外推方法 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.

在精确速度模型的基础上,本文比较了三种偏移算法计算的波场差异.当波通过盐丘时,单程波广义屏方法计算的外推下行波波场与有限差分法计算的下行波波场相比产生较大的误差,如图 4ab中箭头所示,单程波广义屏方法并未准确地计算波场的传播,这将不利于盐丘体下方构造的成像;而本文提出的偏移算法计算的外推下行波与有限差分方法正演模拟计算的下行波较一致,如图 4c中箭头所示.通过与有限差分法和单程波偏移算法计算的波场相比较,本文提出的偏移方法在强速度变化介质中表现出稳健的波场计算性能.

全波形反演中常常使用的是一个相对平滑的速度模型作为反演的初始模型,利用偏移算法开展波场在平滑速度模型中传播规律的研究对偏移成像和全波形反演也至关重要,可为使用本文所提双程波波场深度延拓技术开展全波形反演研究做理论准备.在平滑速度模型中,本文采用上述三种偏移算法计算了t=2.0 s处的波场,如图 5所示.由于采用了平滑的速度模型,平滑速度模型中的断层与围岩的速度差异变小,使得有限差分法正演模拟计算的波场中反射波几乎不可见.同时,在平滑速度模型中,本文提出的偏移方法计算的外推下行波波场与有限差分方法正演模拟计算的下行波波场仍然保持高度一致,而单程波广义屏偏移算法计算的外推下行波波场和理论解释有一定差异,特别是图 5中箭头所指示的波场.此外,与利用准确的速度模型计算的波场相对比,采用平滑的速度模型,本文所提的偏移算法计算的外推波场似乎含有更少的干扰噪声.

图 5 基于平滑速度模型使用不同偏移方法计算t=2.0 s时刻的波场 (a)常规的有限差分法;(b)单程波广义屏方法;(c)基于矩阵分解的双程波深度外推方法. 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.

图 6 速度模型 (a)原始速度模型; (b)平滑速度模型. Fig. 6 The velocity models (a) Original velocity model; (b) Smoothed velocity model.

从单程波广义屏偏移方法计算的成像结果可以看出,单程波广义屏偏移方法对盐丘体的整体成像比较模糊,特别是对盐丘体底界面的成像几乎不可见,如图 7a中箭头所示,并且该偏移结果中存在明显的成像噪声.从整体构造成像效果来看,基于双程波方程的偏移方法(包括逆时深度偏移方法)比单程波偏移方法在成像复杂模型或速度变化剧烈介质时具有明显的优势.对比逆时深度偏移和本文所提偏移算法计算的成像结果,不难发现在相同的参数下,逆时深度偏移算法计算的成像结果出现了较明显的频散噪声,而利用本论文所提的偏移算法计算的成像结果更加稳定、清晰.理论上,可以通过使用更小的网格间隔来压制逆时深度偏移算法中产生的频散,但这将带来更大的计算成本并需要更多的存储空间.为了进行详细的对分析,本文局部放大显示了图 7中虚线框内的成像剖面,如图 8所示.图 8a为速度模型中的层位结构,对比逆时深度偏移算法和本文所提偏移算法计算的局部成像结果(图 8(bc)),可见本文所提偏移算法提供了更加清晰的层位位置和层位接触关系.由于逆时深度偏移算法和本文提出的偏移算法都是基于双程波动方程的,在构造成像中很难看出两种方法的明显区别.本实验的目的是为了证明本文所提的双程波方程波场深度延拓的偏移算法可以和逆时深度偏移算法一样,对复杂介质进行准确的成像,为双程波偏移方法提供另一种选择,而不是代替逆时深度偏移算法.此外,由于本文提出的偏移方法仅在深度方向进行波场外推,因此几乎不对某些特殊结构(例如图 7c中椭圆)进行成像,但这种垂直结构可以通过回转波进行成像(Hale et al.1992; Jia and Wu, 2009).逆时深度偏移方法能在整个速度模型空间范围内对包括回转波在内的各种波场进行正向和反向传播,这使垂直结构能部分被成像,即成像条件计算公式(8)中I1+I2项,如图 7b中的椭圆所示.这也是利用双程波方程开展波场深度延拓及成像需要继续研究的方向之一.

图 7 使用不同的偏移方法计算的成像结果 (a)单程波广义屏方法; (b)逆时深度偏移方法; (c)本文提出的双程波深度偏移方法. Fig. 7 Migration sections calculated by using different methods (a) One-way GSP method; (b) RTM method; (c) Our proposed two-way depth migration method.
图 8 图 7中虚线框中成像剖面局部放大图 (a)速度模型;(b)逆时深度偏移算法;(c)本文所提双程波深度偏移算法. Fig. 8 Partial enlarged imaging sections in dotted box from Fig. 7 (a) Velocity model; (b) RTM method; (c) Our proposed two-way depth migration method.
4 总结

本文以在深度域求解双程波方程所需的边界条件入手,利用单程波偏移算法估计波场对深度的偏导数.基于地表记录的波场及其估计的偏导数波场,建立了含有垂直波数及其三角函数的双程波波场深度外推方程.根据Helmhotz算子与垂直波数之间的数学联系,本文提出使用矩阵分解理论实现双程波方程的波场深度外推.通过利用单程波广义屏方法、有限差分法和本文所提的方法对倾斜界面速度模型和盐丘模型中传播波场的计算,证明了本文所提方法与有限差分方法计算的下行波场精度相同,而常规单程波广义屏方法在强变化介质中计算的外推波场与理论分析存在一定误差.利用常规的逆时深度偏移算法、单程波广义屏方法和本文所提的偏移方法对SEAM模型开展叠前深度偏移成像研究,与单程波广义屏算法计算的结果相比,双程波深度偏移算法在复杂介质成像方面具有很强的优势.此外,在相同的参数下,本文提出的偏移方法比逆时深度偏移方法获得了更清晰的成像结果,频散噪声更小,但逆时深度偏移方法比本文所提深度偏移方法具备一些特殊特征,例如能实现回转波成像.总而言之,本文所提深度偏移方法相较于单程波偏移方法具备更好成像复杂构造的能力,相较于逆时深度偏移算法具有自己的特色和不足,而这些不足还有待进一步的研究.

致谢  特别向各位审稿人对本文提出的建设性意见表示感谢.
References
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.
Gazdag J. 1978. Wave equation migration with the phase-shift method. Geophysics, 43(7): 1342-1351. DOI:10.1190/1.1440899
Gazdag J, Sguazzero P. 1984. Migration of seismic data by phase shift plus interpolation. Geophysics, 49(2): 124-131.
Gray S H, Etgen J, Dellinger J, et al. 2001. Seismic migration problems and solutions. Geophysics, 66(5): 1622-1640.
Grimbergen J L T, Dessing F J, Wapenaar K. 1998. Modal expansion of one-way operators in laterally varying media. Geophysics, 63(3): 995-1005.
Hale D, Hill N R, Stefani J. 1992. Imaging salt with turning seismic waves. Geophysics, 57(11): 1453-1462. DOI:10.1190/1.1443213
Jia X F, Wu R S. 2009. Superwide-angle one-way wave propagator and its application in imaging steep salt flanks. Geophysics, 74(4): S75-S83.
Kosloff D D, Baysal E. 1983. Migration with the full acoustic wave equation. Geophysics, 48(6): 677-687.
Le Rousseau J H, De Hoop M V. 2001. Modeling and imaging with the scalar generalized-screen algorithms in isotropic media. Geophysics, 66(5): 1551-1568.
Li B, Liu H W, Liu G F, et al. 2010. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU. Chinese Journal of Geophysics (in Chinese), 53(12): 2938-2943. DOI:10.3969/j.issn.0001-5733.2010.12.017
Liu F Q, Zhang G Q, Morton S A, et al. 2011. An effective imaging condition for reverse-time migration using wavefield decomposition. Geophysics, 76(1): S29-S39.
Liu H W, Liu H, Zou Z, et al. 2010. The problems of denoise and storage in seismic reverse time migration. Chinese Journal of Geophysics (in Chinese), 53(9): 2171-2180. DOI:10.3969/j.issn.0001-5733.2010.09.017
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.
Ristow D, Rühl T. 1994. Fourier finite-difference migration. Geophysics, 59(12): 1882-1893. DOI:10.1190/1.1443575
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 B Y, Wu R S, Gao J H. 2012. Preliminary investigation of wavefield depth extrapolation by two-way wave equations. International Journal of Geophysics, 2012: 968090.
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.
You J C, Li G C, Liu X W, et al. 2016. Full-wave-equation depth extrapolation for true amplitude migration based on a dual-sensor seismic acquisition system. Geophysical Journal International, 204(3): 1462-1476. DOI:10.1093/gji/ggv535
You J C, Wu R S, Liu X W, et al. 2018a. Two-way wave equation-based depth migration using one-way propagators on a bilayer sensor seismic acquisition system. Geophysics, 83(3): S271-S278.
You J C, Wu R S, Liu X W. 2018b. One-way true-amplitude migration using matrix decomposition. Geophysics, 83(5): S387-S398.
Zhang J F, Liu L N. 2007. Optimum split-step Fourier 3D depth migration:Developments and practical aspects. Geophysics, 72(3): S167-S175.
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.
Zhang Y, Zhang G Q, Bleistein N. 2003. True amplitude wave equation migration arising from true amplitude one-way wave equations. Inverse Problems, 19(5): 1113-1138. DOI:10.1088/0266-5611/19/5/307
Zhang Y, Zhang G, Bleistein N. 2005. Theory of true-amplitude one-way wave equations and true-amplitude common-shot migration. Geophysics, 70(4): E1-E10.
李博, 刘红伟, 刘国峰, 等. 2010. 地震叠前逆时偏移算法的CPU/GPU实施对策. 地球物理学报, 53(12): 2938-2943. DOI:10.3969/j.issn.0001-5733.2010.12.017
刘红伟, 刘洪, 邹振, 等. 2010. 地震叠前逆时偏移中的去噪与存储. 地球物理学报, 53(9): 2171-2180. DOI:10.3969/j.issn.0001-5733.2010.09.017
张贤达. 2015. 矩阵分析与应用. 2版. 北京: 清华大学出版社.