多步法逆时偏移是Guan H在2008年提出的方法, 其主要思想是根据偏移模型的地质特点, 将速度模型分成2~3块不等的深度区, 分别进行逆时偏移或多种偏移方法的处理(图1)。首先在区域I应用RTM, 并将震源记录和检波点记录外推至区域I和区域II的交界深度处作为区域II的输入数据, 然后在区域II进行RTM; 同样将震源记录和检波点记录外推至区域II和区域III的交界深度处作为区域III的输入数据, 然后在区域III进行RTM。

图1

	Figure Option ViewDownloadNew Window
	图1 多步法逆时偏移速度分区示意[15]

多步法逆时偏移提高了计算效率, 减少了内存需求。因为偏移的网格尺寸是由偏移速度场的最小速度决定的, 因此随着区域深度的增加, 可以采用更大的偏移网格。较大的偏移网格在提高计算效率的同时, 还可以减少存储震源波场数据的大小, 同时减少内存存储量。

多步法偏移方法还可以在不同的速度区域使用不同的成像方法, 笔者采用的步骤如下。

(1)根据东濮地下构造的特点将偏移速度场分成两个区域:表层构造平缓、速度偏小的区域和中深层断块发育、速度较高的区域; 深度分界线为1 200 m, 1 200 m以上为区域I采用10 m网格剖分, 1 200 m以下为区域II采用20 m网格剖分。

(2)对区域I进行偏移成像, 进行常规克希霍夫积分偏移。偏移后将接收面下延至分界面处, 得到在分界面处接收的野外炮集记录。

(3)对区域II进行常规逆时偏移成像。

(4)叠加两部分偏移结果, 得到完整的偏移结果。

GPGPUs是最新的一种高性能计算平台, 目前最高端的GPU(Tesla K20x)拥有2 688个核, 文中使用的M2095型号显卡也拥有512个核心。GPU中大量的计算单元使其具有海量计算的能力, 可极大地提升计算效率。

逆时偏移成像部分主要利用高阶有限差分方法求解声波波动方程, 三维声波方程的一阶应力— 速度方程形式可表示为

∂P∂t=-ρc2(∂vx∂x+∂vy∂y+∂vz∂z), ∂vx∂t=1ρ∂P∂x, ∂vy∂t=1ρ∂P∂y, ∂vz∂t=1ρ∂P∂z。(1)

其中: v_x, v_y, v_z是质点速度, P为法向应力, c为速度, ρ 为密度。从方程组中可以看出, 主要的工作量也是最耗时的工作量集中在一阶波场导数计算上^{[16, 17, 18]}。因此将计算波场变量的空间导数部分交给GPU进行细粒度并行计算。

使用空间高阶差分法需要大量的内存读写, 以三维8阶差分格式为例, 每计算一个网格点的值需要读取网格点周围25个网格点的数据, 内存读取冗余度很高。针对这一问题, 利用GPU片内访问速度最快的共享内存代替全局内存来减少I/O, 可以极大地节省计算时间。

在开辟共享内存数组时, 一般将数组尺寸定义为16× 16, 但是在有限差分求解中, 每一个网格点的计算都需要周围几个点的参与, 因此需要在数组外进行扩边, 扩边的大小和有限差分阶数有关。假设有限差分求解阶数为2× m(m为偶数), 则开辟共享内存数组尺寸为[16+2× m, 16+2× m]。以8阶差分格式为例, 需开辟共享内存数组尺寸为[16+2× 4, 16+2× 4], 设定计算线程块(Block)大小为16× 16, 并根据计算数组大小得到Grid大小, 即Block的维度。通过这种三维结构体的索引, 将计算数组和计算线程联系起来, 完成向共享内存数组的数据拷贝和扩边。如图2所示, 数组主体部分为当前线程块需要计算的数据, 扩边区域为两个线程块之间共享的数据, 扩边宽度为差分格式长度的1/2。这种共享内存的利用降低了反复从全局内存中读取数据造成的读取冗余。

图2

	Figure Option ViewDownloadNew Window
	图2 16× 16的共享数组和8阶差分格式的扩边[19]

逆时偏移使用的互相关成像条件需要同一时刻的震源正向外推波场和检波点逆向外推波场。这两个波场一个是在时间的正向上进行, 一个是在时间的反向上进行, 为了能同时得到这两个过程的波场, 必须存储其中一个过程的波场。这就需要消耗极大的存储资源, 尤其在实际三维操作中是难以满足的。Eric^[20]归纳总结了5种计算策略, 震源重构计算策略就是其中一种。

震源重构计算策略是指, 利用GPU计算正向波场的过程中, 将各个边界波场值存储至CPU上。在利用GPU计算检波点波场的过程中, 同时逆向外推震源波场, 并在对应时刻将存储在CPU上的边界波场值应用到各个边界上, 以保证该过程得到的每个时刻的波场值和原来正向外推该时刻的波场一致。震源重构计算策略可以在外推过程中边外推边成像, 大幅降低存储要求。

本次处理涉及2个三维工区:前梨园北三维工区和文北高精度三维工区。前梨园北工区为6线24炮采集, 道间距50 m, 面元大小为25 m× 25 m, 满覆盖次数30次。文北工区为高精度采集数据, 32线10炮采集, 道间距50 m, 面元大小为5 m× 10 m, 满覆盖次数32次。工区位置见图3所示, 表1为前梨园洼陷目标处理工区野外地震采集的一些主要参数。

图4、图5分别为复杂小断块区域和兰聊大断层区域的叠前深度克希霍夫和逆时偏移成像效果对比。克希霍夫积分深度偏移是建立在射线理论基础上的偏移算法, 存在焦散、多重路径、干涉现象和走时计算不准确的问题, 造成速度分析的多解性, 无法实现真正的共反射点叠加, 从而对复杂构造难以准确成像。从两种成像方法得到的对比剖面中可以看出, 图4b中逆时偏移对复杂小断块成像更为精确, 断块边界刻画清晰, 断点归位准确。另外, 逆时偏移对东濮兰聊大断层的准确归位是一个很好的解决方法, 图5b中逆时偏移对兰聊大断层的成像更精确, 断面归位准确, 断层两侧地层接触关系清晰。

图3

	Figure Option ViewDownloadNew Window
	图3 工区位置

表1

表1

表1 前梨园洼陷目标处理工区野外地震采集部分参数 工区名 前梨园北 文北 年度 2000~2001 2008~2009 观测方式 6线24炮 32线10炮 记录道数 960道 4096道 检波线距 400m 100m 道间距 50m 50m 炮线距 350m 80m 炮点距 50m 80m 观测方式 端点放炮 端点放炮 面元大小 25m× 25m 5m× 10m 覆盖次数 10× 3 8× 4 采样间隔 2ms 1ms 记录长度 8s 6s

表1 前梨园洼陷目标处理工区野外地震采集部分参数

逆时偏移在生产中无法大规模推广应用的一个重要原因是偏移成本太高。逆时偏移依赖于基于全波动方程的波场外推技术, 需要满足一定的稳定性条件和频散条件, 因此计算网格受限于偏移频率和最小偏移速度。如果计算网格选取过小会极大地增加偏移耗时, 如果计算网格选取过大则会降低最终的成像精度, 因此在偏移成像过程中, 需要平衡成像结果精度和偏移计算量之间的关系。笔者选择了使用MPI+CUDA并行的逆时偏移成像方法提高计算效率, 并利用内存优化减小数据I/O。

图4

	Figure Option ViewDownloadNew Window
	图4 复杂小断块区域叠前深度克希霍夫(左)和逆时偏移(右)成像效果对比

图5

	Figure Option ViewDownloadNew Window
	图5 兰聊大断层区域叠前深度克希霍夫(左)和逆时偏移(右)成像效果对比

从东濮前犁园三维地震数据中选取单炮, 测试并行后的逆时偏移算法和并行后进行内存优化的逆时偏移算法效率, 选取的单炮炮检关系如图6所示。其中CPU计算和GPU计算均使用单节点。GPU使用Tesla M2090显卡, 含有512个计算核心; CPU型号为 Intel(R)Core(TM)i7-2600, 8核心。考虑到CPU的计算能力, 采用较低的偏移主频, 可以在一定程度上减小计算量, 但相应的偏移成像精度也会受到影响。在子波主频20 Hz, 20 m空间采样, 1 ms时间采样下, 采用3种策略进行计算, 不同处理器计算耗时如下。

图6

	Figure Option ViewDownloadNew Window
	图6 第50267炮炮检关系

策略1:CPU单炮偏移耗时为:43 h 28 min 37 s, 约2 608.6 min;

策略2:GPU单炮偏移耗时为:2 558.474 s, 约42.65 min;

策略3:利用共享内存优化后, GPU单炮偏移耗时为:2 084.961 s, 约34.75 min。

从数据对比中得到, CPU单炮偏移耗时约为GPU单炮偏移耗时的61.16倍, 是优化内存后GPU单炮耗时的75.07倍, 即文中使用的基于GPU的三维叠前逆时偏移成像方法在东濮前梨园地区计算效率是基于CPU的逆时偏移成像方法的75倍。