1 引言

逆时偏移的历史可以追溯到1978年.1978年Hemon^[1]阐述了用有限差分法求解波动方程实现逆时偏移的思想，但没有强调这种方法的潜力.而真正将逆时偏移引入地震勘探这一领域的是1983年几位科学家独立的科研工作：Baysal等^[2]采用拟谱法实现逆时偏移，他的方法规定了波场传播算子的符号，成像角度限制在90°以内；McMechan等^[3]采用时间和空间域二阶有限差分法求解全声波方程，突破了倾角的限制；Loewenthal等^[4]在频率波数域利用指数解析解求解波动方程，计算速度快，但是只适合常速介质；Whitmore等^[5]利用逆时偏移结果拾取层位，更新速度模型，迭代求取速度.所有这些方法最初都是用做叠后偏移.由于逆时偏移计算量及存储量巨大，无法适应工业生产的需要，逆时偏移一直没有在工业界得到广泛应用.相反，Kirchhoff积分法以及单程波方法由于其在计算量以及存储量方面的优势一直占据着地震偏移领域的主导地位.近年来，由于计算机技术的飞速发展以及对复杂构造油气藏成像要求的提高，逆时偏移重新被认识开发.目前，逆时偏移技术已经从叠后走向叠前，从二维走向三维^[6]，从声波方程走向弹性波方程^[7]，从各向同性介质走向TTI（Tilted Transversely Isotropic，倾斜横向各向同性）介质^[8~11].

国内对逆时偏移的研究可追溯到20世纪90年代，邓玉琼等^[12]将弹性波有限元反时偏转法用于同时对合成地表反射数据和合成VSP反射数据进行偏移，不但可以改正因使用范围有限的地表数据而导致的弯曲畸变，反射点也可被恰当地成像于他们的真实深度处；王山山等^[13]采用混合法求解双程无反射波动方程，这个方法沿x方向使用伪谱法，沿z方向采用有限元法，沿时间方向有限差分，是算法灵活性，计算效率以及内存需求的一种折中的选择；底青云等^[14]开展了弹性波有限元逆时偏移技术研究；张美根等^[15]将弹性波有限元叠前逆时偏移技术应用到各向异性介质；王童奎等^[16]用谱元法求解波动方程实现叠前逆时偏移；薛东川等^[17]采用最大振幅成像条件替代最小走时成像条件；杜启振等^[18]将叠前逆时偏移技术应用到横向各向同性介质弹性波多分量数据，陈可洋^[19]用高阶有限差分法求解波动方程实现叠前逆时偏移.

近年来，国内外的逆时偏移研究成功的案例并不多见，重要的原因是逆时偏移研究耗费巨大，30余年来，缺乏一种快速的、适合工业生产规模的逆时偏移实现方法.正如前文提到的，逆时偏移的一大瓶颈是计算量太大.与单程波方法相比，单程波方法沿着深度方向递推，只需要计算地震数据有限频带内的信息；逆时偏移沿时间方向递推，由于算法稳定性以及频散关系等要求，时间递推步数要远大于单程波深度递推时所需要的频率数目.如何提高逆时偏移的计算效率成为逆时偏移能否大规模应用的重要环节.目前，逆时偏移提速主要有两种方法：从算法角度减少计算量和高性能并行计算.减少计算量的方法主要有：Yue Wang^[20]提出一种变时间步长和空间网格大小的方法；Guan等^[21]提出了一种分段多步偏移方法，将偏移区域在深度方向按照构造特征划分为几个连续的空间，对每个空间依次偏移，并在空间拼接区域保存波场信息用于下一个空间的偏移，这种做法的优点是对每个偏移空间所选择的空间网格以及时间步长可以不同，而且不同的空间可以选择不同的偏移方法，可以显著降低计算量；Zhang等^[22]通过引入复波场以及平方根算子将二阶声波方程降为类似于单程波的时间一阶导数的偏微分方程，采用可分近似求解，算法突破了稳定性及频散的影响，时间步长可以取到Nyquist频率，是一种高效率高精度的算法；Robert Soubaras等^[23]指出单步法计算过程中拟微分算子难以合成的缺点并采用双步法方程克服了这一缺点，同时保持了单步法在时间步长方面的高效性.逆时偏移提速的另一种方法依靠计算机硬件性能的提高：Matthew H Karazincir等^[24]比较了伪谱法和高阶有限差分两种方法，高阶有限差分计算量小，不存在Fourier变换产生的噪音，从稳定性的角度比较，时间间隔可以更大，并且便于并行实现，Villarreal等^[25]利用这一思想，采用大规模集群并行实现波动方程有限差分法正演模拟；多核处理器（如Cell，FPGA以及GPU）的出现，为RTM提供了新的机遇，Araya-Polo M等^[26]首先用Cell处理器实现了三维逆时偏移，Micikevicius P^[27]给出了利用GPU实现高阶有限差分的算法；Darren Foltinek等^[28]在利用GPU实现逆时偏移方面做了有益的尝试.本文的一项主要工作就是围绕利用GPU加速逆时偏移展开.

GPU的全称是图形处理器（Graphic Processing Unit），是NVIDIA公司推出的一款用于计算机显示的设备，俗称“显卡”.2006年NVIDIA公司发布了统一计算设备架构平台CUDA（Compute Unified Device Architecture），可让开发商编写指令程序，以便让GPU执行通常由CPU执行的计算任务.由于GPU在科学计算方面的巨大优势以及在计算机中的重要作用，GPU的含义可能扩展为通用处理器（General PurposeUnit）.和传统的集群相比，GPU具有计算核心多，硬件成本低，省电，省空间等优点，目前已经广泛应用于生命科学、医学设备、工业生产、电子设计自动化、制造、金融和通讯等行业.在石油天然气行业，国际上Hess、CGGVeritas、Chevron、Headwave、Acceleware、Seismiccity等公司和组织都在致力于这方面的研究.中国科学院地质与地球物理研究所刘洪研究员与吉星吉达公司共同合作开发的非对称走时叠前时间偏移GPU算法^{[29, 30]}已经在胜利、大庆等油田投入生产，效果良好.GPU已经逐渐成为下一代高性能计算的首选工具.

本文首先给出了波动方程的高阶有限差分解法，并分析了算法稳定性条件以及频散关系；然后介绍了逆时偏移的实现过程以及GPU实现流程；最后对脉冲响应、垂直断层以及BP模型做了测试，通过与单程波方法比较，阐明了逆时偏移方法对陡倾构造成像方面的优势.

2 波动方程高阶有限差分解

三维均匀横向各项同性介质声波方程如下：

(1)

式中速度变量v是空间变量的函数，我们采用规则网格高阶有限差分解法求解声波方程.

2.1 差分格式

规则网格中心差分系数的求解通常借助泰勒展开，生成系数矩阵，利用Cramer求取系数.以二阶导数为例，要求2n阶差分格式的系数，需要求解以下系数矩阵：

(2)

应用Cramer法则可以解此方程，主要涉及两个Vandermonde行列式的求解.差分系数有如下形式：

(3)

以x方向的空间导数为例，差分格式如下：

(4)

其中，

(5)

差分格式的求取还有另外一种更为简洁的方法，Bengt Fornberg ^[31]利用递推法推导了高阶规则网格与交错网格、中心差分以及单边差分格式，通过与求解系数矩阵求得的系数比较，二者完全相同，文中我们采用递推法推导了2~20阶差分格式的系数，从算法稳定性以及频散关系的角度考虑，时间方向采用二阶差分格式，空间采用高阶差分格式.

2.2 稳定性分析

采用状态传递矩阵特征值^[32]方法，可以求得时间域二阶差分，空间域高阶差分格式的稳定性条件为：

一维情形：

(6)

二维情形：

(7)

三维情形：

(8)

将稳定性条件系数作为方程维数及空间差分阶数的函数，如图 1所示：阶数越高，维数越高，则稳定性条件系数越小，相应的时间步长越小，导致计算量与存储量越大.从而可以得出结论：进行数值模拟时，从二维到三维，不仅仅是计算的空间网格点数增加，时间递推的步数也增加.

2.3 频散关系

差分算子对微分算子的近似会造成数值频散，时间方向和空间方向的差分都是如此.下面，我们以二维声波方程空间差分为例推导频散关系.

二维空间高阶差分格式如下：

(9)

将p（t，x，z）=exp[i（wt-kcosθx-ksinθz）]以及ω²=k²v²带入上式，可得频散关系：

(10)

从（10）式可以发现，v < v₀（v₀是源点的速度），也就是说，由于空间离散造成的数值频散在波形上作为尾巴而出现.空间差分引起的数值频散由三个因素决定：一是地震波传播方向，在dx=dz的情况下，角度为45°时数值频散最小，如图 2a和图 3a所示；二是空间差分精度，随着空间差分精度的提高，差分解法产生的数值频散会逐渐减小，如图 2所示.因此，可通过提高差分精度的办法来减小数值频散；三是一个波长内离散点数目，对任意阶空间差分精度，一个波长内离散点数越少，数值频散越严重，如图 2b所示.也就是说，在相同的网格间距情况下，子波频率越高，介质速度越低，频散越严重.因此，海水以及陆地表层低降速带以及高频子波对地震波模拟结果有很大影响，这一结论可以解释文献[20]中变时间步长和空间网格大小的方法.在实际应用中，应根据不同的地质模型以及数据频率信息来选择不同的差分精度.

采用类似的推导方法，可以推导出时间差分引起的数值频散关系如下：

(11)

可以发现v>v₀，由于时间离散造成的数值频散在波形上在正常到达时之前出现.时间差分引起的数值频散由两个因素决定：差分精度越高，一个时间周期内离散点数目越多，由时间离散造成的数值频散越低.也就是说，时间递推步长必须满足稳定性条件与频散关系两方面的要求.在逆时偏移中，由于算法稳定性的要求，dt一般很小，时间方向二阶差分精度足以满足频散关系的要求，地震波传播中的数值频散主要是由空间离散造成的，因此，空间方向要采用高阶差分格式.

图 3是一个简单的常速模型，震源位于地表中央处，模型速度为3000m/s，dx=dz=12.5m，地震波主频为20 Hz，时间导数采用二阶差分.图 3a是空间导数二阶差分的结果，频散严重；图 3b是空间导数十二阶差分的结果，频散得到了有效的压制.

3 叠前逆时偏移实现过程及GPU加速

炮域叠前逆时偏移的实现过程包括三个步骤：首先炮点波场沿时间正方向传播到最大时刻并将中间结果保存下来；然后检波点波场数据沿时间反方向传播，并读取炮点相应时刻的波场值；最后应用合适的成像条件，将所有炮的偏移结果累加得到最终偏移结果.波场传播采用上文中介绍的时间域二阶差分、空间域高阶差分格式.Sandip Chattopadhyay等^[33]对逆时偏移成像条件做了比较系统和全面的介绍，本文中我们选用互相关成像条件，原因是互相关成像条件容易实现，便于并行，不存在稳定性问题，并且可以处理多波至问题，不会丢失波场信息.

(12)

式（9）和式（12）表明，波场传播和应用成像条件时空间中每个网格点都是解耦的，即独立的.所有网格点可以并行计算.并行粒度很小，每个线程只计算一个网格点或者几个网格点的值，与CPU集群适合粗粒度并行不同，GPU计算核心多，对这种细粒度并行的问题更有优势.实验过程中选用的GPU型号是Tesla S1070^[34]，它由四块Tesla1060显卡组成，每块显卡有240个计算核心和4GB的显存，选用的CUDA版本是2.3.1^[35].作为对比，CPU的型号是Pentium（R）Dual-CoreCPU E5200 @ 2.50GHz，配有2GBDDR2内存.

图 4是CPU/GPU协同计算实现逆时偏移的流程图，计算密集的波场传播与成像部分借助GPU并行计算，中间结果需要借助通讯保存在CPU内存或硬盘上.需要指出的是，运用空间高阶差分法需要大量的内存读写，以三维8阶差分格式为例，每计算一个网格点的值需要读取网格点周围25个网格点的数据，内存读取冗余度（read redundancy）很高.针对这一问题，我们采用文献[25]中的方法，运用矩阵分片的思想，借助GPU共享存储器（share memory）机制重用内存数据，可以大幅度降低内存读取冗余度.

4 模型试算

我们针对垂向变速介质的脉冲响应，以及垂直断层和BP模型分别使用单程波偏移和逆时偏移方法进行了测试，测试结果如下.

4.1 脉冲响应

NoletG^[36]指出，对如下线性变速介质：

(13)

t时刻的波前面解析解方程如下：

(14)

式中，v'是速度梯度，α是速度梯度方向与z方向的夹角.波场t时刻波前面是半径为v₀/v'sh（v't），圆心为{x₀-v₀/v'sinα[1-ch（v't）]，z₀ -v₀/v'cosα[1-ch（v't）]}的圆.

测试过程中，我们选用的是仅在垂直方向存在变速的介质，其中v₀=1800 m/s，v'=2 m/s/m，dx=dz=8m，震源放在地表模型的中心位置1036 m处，选用主频为20 Hz的Ricker子波为震源，波场在t=0.45s时刻的波前面是一个半径为924m，圆心为（1036m，390m）的圆，如图 5中的白线所示.分别用高阶有限差分逆时偏移和分裂步相移^[37]单程波偏移方法计算脉冲响应，如图 5所示.对于横向不变速介质，分裂步相移法得到的脉冲响应走时信息是准确的，和理论解吻合很好，但对于角度大于90°的波场无法成像，如图 5b所示，这是传统单程波方法固有的局限；逆时偏移求解的是全声波方程，没有做任何的近似，突破了90°倾角的限制，可以对任意角度成像，如图 5a所示.

4.2 垂直断层

第二个实验是如图 6a所示的垂直断层模型，成像难点是断层的垂直部分.我们使用基于GPU/CPU协同计算的高阶有限差分法正演模拟得到正演数据，图 6b给出了某一时刻的波场快照，箭头方向指示了棱柱波（Prismaticwave）的传播路径，在剖面上是一条与一次反射波相切的同相轴，如图 6c中箭头所示.

对于棱柱波，传统的单程波方法是无法成像的，因为波场传播方向发生了改变，而逆时偏移可以对棱柱波准确成像.图 7a是逆时偏移的结果，断层的垂直部分成像清楚；图 7b是Fourier积分单程波^[38]方法的成像结果，这种方法直接运用准确的单程波算子积分，不做任何的近似，是一种可以精确到90°的单程波方法，但依然无法对断层的垂直部分成像.

4.3 BP盐丘模型

BP^[39]二维模型是2004年6月由Frederic Billette and Sverre Brandsberg-Dahl在巴黎第66届EAGE年会上发布的.这个模型由三部分组成：模型的左侧部分模拟的是穿过墨西哥湾西部（Western Gulf of Mexico）的横切面；中间部分模拟的是墨西哥湾中东部（Eastern/Central Gulf of Mexico）与安哥拉（off-shore Angola）近海岸的地质构造；右侧部分模拟的是里海（Caspian Sea）、北海（North Sea）以及特立尼达岛（Trinidad）的地质构造.图 8是中部的速度模型，模型的难点是盐丘内幕构造成像，需要借助上面提到的棱柱波成像.

图 9（a，b）分别是逆时偏移过程和Fourier积分法单程波偏移过程中炮点波场传播的快照，白色箭头指示了逆时偏移过程中棱柱波的传播过程，单程波方法无法模拟棱柱波，也就无法利用棱柱波成像，因为单程波偏移过程中波场沿z方向的传播方向不能改变.

图 10b是Fourier积分单程波偏移方法的结果，盐丘顶部以及低速区成像清楚，但盐丘内壁以及盐丘右侧垂直边界无法成像；图 10a是逆时偏移的成像结果，盐丘内壁以及盐丘右侧垂直边界成像清楚，逆时偏移优势明显.

逆时偏移测试过程中，我们选用了时间2阶中心差分，空间12阶中心差分算法，分别在CPU和GPU上测试，测试结果表明，Tesla1060的计算速度比CPU快30倍左右.

5 结语

本文介绍了高阶有限差分法求解声波方程的稳定性条件和频散关系，据此给出了叠前逆时偏移的实现过程.通过和单程波方法相比较，逆时偏移突破了90°倾角的限制，对垂直断层，盐丘内幕等陡倾构造成像效果显著，是单程波方法无法比拟的.但是逆时偏移计算量大，计算周期长，为此本文选择GPU实现算法加速，取得了30倍的加速比.

本文关于GPU加速的研究仅限于单个GPU细粒度并行，对于大规模的三维叠前逆时偏移，单个GPU是无法满足生产规模的需要的，为此，需要多GPU间粗粒度并行、单GPU内部细粒度并行的算法，关于这方面的工作正在开展.

计算量大是叠前逆时偏移的瓶颈之一，成像噪音消除以及存储量大等问题是叠前逆时偏移的另外两个瓶颈，将另文阐述.

致谢

对BP公司Frederic Billette提供的二维BP模型数据，谨致谢意，感谢李幼铭研究员对本文的指导和帮助.