三维平面波叠前时间偏移实用化研究

段心标，方伍宝

（中国石油化工股份有限公司 石油物探技术研究院，江苏南京 210014）

三维平面波叠前时间偏移实用化研究

段心标，方伍宝

（中国石油化工股份有限公司 石油物探技术研究院，江苏南京 210014）

三维平面波叠前时间偏移技术实用性不强，浅层成像信噪比较低，在计算过程中内存开支大，并行节点有效利用率低，且不具有断点保护功能，无法适用于大规模叠前偏移成像处理。为了提高三维平面波叠前时间偏移的实用化水平，从平面波叠前时间偏移原理和实现方案出发，针对存在的问题，提出了相应的解决方案，包括通过倾角滤波方法压制浅层的偏移成像噪声，提高浅层成像信噪比；采用分层成像方案，减小计算过程中的内存开支，提高并行节点的有效利用率；开发断点保护功能，考虑了长时间偏移计算中的断点风险。在应用于某工区实例资料处理后表明，经优化后的平面波叠前时间偏移技术较原平面波偏移具有更好的成像精度，并且具备了大规模数据处理能力。

三维平面波叠前时间偏移；倾角滤波；分层成像；断点保护

0 前言

目前在生产中，常用的Kirchhoff积分法叠前时间偏移［1～4］，对速度横向缓慢变化介质的适应性不强，并且该方法通常是不保幅的。为了发展适应弱横向变速介质，并具有相对保幅特性的叠前时间偏移技术，王华忠［5］在Mosher［6］平面波偏移思想的基础上，提出了二维偏移距域平面波有限差分叠前时间偏移方法。Feng［7］把该方法推广到三维情况，并提出了一种与方位角无关的有限孔径倾斜叠加平面波分解［8、9］方法。这种三维平面波叠前时间偏移方法，基于时间域层速度采用有限差分算法对平面波波场进行外推，可以很好地实现复杂构造的偏移成像。但是，该偏移方法的实用性并不强，主要表现以下为以下几点：

（1）浅层成像信噪比较低。

（2）计算过程中内存开支大，并行节点有效利用率低，不能适用于大规模数据处理。

（3）不具备有断点保护功能，没有考虑长时间偏移计算中的断点风险。

为了提高三维平面波叠前时间偏移的实用化水平，作者在本文中，从平面波叠前时间偏移原理和实现方案出发，针对上述几个问题采用有效的解决手段，有效地压制了浅层偏移噪声，并且使得偏移算法具备了大规模数据处理的能力。

1 平面波偏移

1.1 偏移原理

由时间域双平方根算子的频散关系，可以推导出如下方程：

其中 kτ为双程旅行时对应的波数；γ为炮检点射线的半张角为常速度。

式（1）在频率～空间域可以表示为：

其中 v＝v（x，y，τ）；U珦＝U珦（x，y，τ；ω）。

γ与平面波矢量p珝h有如下关系式：

进而得三维偏移距平面波方程：

对每一个平面波数据求解方程（4），即可实现偏移成像。叠加所有平面波的成像结果，可以得到平面波叠前时间偏移最终成像剖面。

1.2 实现方案

一个单ph平面波数据的偏移成像，类似于一个叠后波动方程偏移，在逐次对所有平面波数据进行偏移处理后，即可得到成像数据体。平面波叠前时间偏移流程如图1所示。

（1）在单ph平面波偏移中，主进程从磁盘中读取三维平面波数据，并通过傅里叶变换，把数据变换到频率空间域中，然后按频率片把数据发送给从进程。

（2）从进程接收到单频数据后，逐时间片对波场进行向下延拓和成像，并累加本进程其它频率成像结果。

（3）当所有频率成像结束后，主进程归约累加所有进程的成像值，可以得到当前单ph平面波的偏移结果，并写入数据磁盘。

作者在平面波偏移算法实现中用到了MPI主从模式并行技术，主进程负责分发作业，从进程负责对每个单频平面波数据进行偏移计算。一般而言，偏移中频率数为300～500，因此该实现方案可以有效地利用计算机集群的运算能力。

在程序执行中，主进程要为一个复数类型的三维频率空间域平面波数据和二个实数类型的三维时间空间域成像数据开辟内存空间，从进程要为一个实数类型的三维时间空间域成像数据开辟内存空间。一般情况下，在一个满覆盖面积为100km2的地震数据平面波叠前时间偏移处理中，主进程需要5G的内存空间，从进程需要2G的内存空间。通常而言，计算机集群无法满足如此大的内存需求，这就制约了平面波叠前时间偏移的大规模数据处理应用，使得该偏移技术的实用性不高。

此外，由于偏移成像需要大量的计算时间，在长时间计算中不可避免的存在着计算机节点死机、节点通信不稳定等风险。而在平面波偏移实现流程中，并没考虑到这些风险带来的危害，程序也不具备断点保护功能。如果出现断点情况，则需要重新偏移，会浪费大量的计算时间，这也是该偏移技术实用化不高的原因之一。

图1 平面波偏移主体流程Fig.1 Flow chart of plane－wave migration

2 实用化策略

平面波叠前时间偏移实用性不强，首先表现为浅层偏移噪声严重，成像信噪比较低。图2（a）是某资料的平面波偏移成像结果，图2（b）为其积分法偏移结果。通过对比可见，平面波偏移信噪比低，浅层同相轴不清楚。此外，还有二个阻碍平面波叠前时间偏移实用化的影响因素，如前面分析它们分别是程序内存开支过大和不具有断点保护功能。

下面，作者将从压制浅层偏移噪声、减小内存开支和增加断点保护等三个方面，对平面波叠前时间偏移进行优化。

2.1 浅层偏移噪声压制

平面波偏移浅层噪声问题是由浅层的偏移孔径过大引起的，所以要想控制三维波动方程偏移的偏移孔径，需要对单道偏移响应进行倾角滤波［10、11］。

在偏移计算时，浅层的滤波角度应该较小，随着深度的变化，滤波角度也要逐渐变大。

深度Z方向的视倾角为：

其中 kz为Z方向的视波数，且又有：

由公式（1）可得：

把公式（6）和公式（7）代入到公式（5）中，则有：

令视倾角截断值为cosα，则当Kz小于cosα时，把波场值赋为0，即可实现偏移响应倾角滤波。

图3是ph＝0时一道信号的偏移响应，图3（a）为无倾角滤波，图3（b）为有倾角滤波，且滤波角度从浅至深逐渐变大。通过倾角滤波，平面波叠前时间偏移的偏移孔径得到控制。采用此方法对图2资料做平面波偏移处理，滤波角度α设定为40°，偏移成像结果如图4所示。

从图4可以看出，浅层噪声得到了较好的压制，信噪比较图2（a）有明显提高。

图2 某资料浅层偏移结果Fig.2 The migration results on shallow profile of a seismic data

2.2 分层成像

由于需要为成像数据体开辟很大内存空间，原平面波偏移不具备大规模数据处理能力。为了减小内存开支，我们提出了分层成像偏移技术。

（1）首先由节点最大可用内存，自动计算分层成像时每个层的厚度和成像层数。

（2）主进程在读取平面波数据后，并将其转换到频率空间域，按频率逐个发送给从进程。

（3）从进程每次对一定厚度的层进行成像，在本层计算结束后主进程归约计算结果写入数据盘并且从进程把延拓后的波场传回到主进程。

（4）然后主进程把频率空间域的波场值重新发送出去，进行下一个层的成像。

分层成像实现方案见下页图5。

假设偏移延拓总步数为Nall，每层时间点厚度为Np，那么实数类型的三维时间空间域成像数据内存大小，仅为原来的Np/Nall倍，这样就有效地减小了内存开支，提高了节点利用率。

图5 分层成像平面波偏移主体流程Fig.5 Flow chart of layer－by－layer plane－wave migration

2.3 断点保护

断点保护策略是当某ph偏移完成后，记录相关信息到日志文件，如果断点后重新提交作业，则从已计算结束的下一个ph开始偏移。该断点保护方案不用写或读临时数据，只需主进程写或读日志文件，记录或获取偏移运行信息。断点后重新提交作业时，主进程将从日志文件中获取的信息广播告知从进程，由从进程按此信息进行偏移计算。该断点保护方案使得并行偏移程序在节点选择时具有较强的灵活性，断点保护前后可以使用完全不同的节点。

考虑断点保护的偏移流程如图6所示。

3 实例

作者利用实用化优化以后的平面波叠前时间偏移技术，处理某工区实际资料。该资料满覆盖面积为54.17km2，面元大小为25m×25m，采样时间为6s，采样间隔为2ms，采样点数为3 000，CMP道集数据量为70G。若采用原平面波偏移技术，主进程内存需求为2 457M，每个从进程的内存需求为1 006M。采用实用化优化后的平面波叠前时间偏移时，假设最大可用内存空间为350M，则程序自动对偏移成像深度进行分层，共分为三层，每层厚度为1 000个采样点。

该资料经实用化的平面波偏移技术处理后，成像剖面见下页图7（a）。图7（b）（见下页）是原平面波偏移成像剖面，图7（c）（见下页）是商业软件叠前时间偏移剖面。对比图7（a）和图7（b）可以看出，平面波偏移实用化优化以后，浅层噪音得到很好的改善。对比图7（a）和图7（c）则可知，优化后的平面波偏移技术具有很好的偏移成像精度，在成像效果上与商业软件叠前时间偏移相当。在某些位置平面波偏移成像的层间反射内幕更清晰，同相轴连续性更好，如蓝色圈标识部份。

4 结论

针对平面波叠前时间偏移实用化程度不高，无法适应大规模数据处理的问题，作者给出了具体的解决方案，形成了实用化的平面波叠前时间偏移技术。通过倾角滤波方法压制了浅层的偏移成像噪声，采用分层成像方案，减小了计算过程中内存开支，设计开发断点保护功能，考虑了长时间偏移计算中的断点风险。经实例表明，优化后的平面波叠前时间偏移技术具有很好的成像精度，并且具备了大规模数据处理能力。

图6 断点保护平面波偏移流程Fig.6 Flow chart of plane－wave migration with breakpoint protection

图7 某工区资料偏移成像剖面Fig.7 The imaging profiles of a seismic data

致谢：

感谢同济大学海洋与地球科学学院“波现象与反演成像”研究组提供的支持和帮助。

［1］ SCHNEIDER W A.Integral formulation for migration in two and three dimensions［J］.Geophysics，1978，43（1）：49.

［2］ LUMLEY D，CLAERBOU J，BEVC D.Anti－aliased Kirchhoff 3Dmigration［J］.Expanded Abstracts of 64thAnnual International Meeting，SEG，1994：1282.

［3］ RESHEF M.Depth migration from irregular surfaces with depth extrapolation methods［J］.Geophysics，1991，56（1）：119.

［4］ CAUSSE E，HAUGEN G U，ROMMEL B E.Large－offset approximation to seismic reflection traveltimes［J］.Geophys，2000，48（8）：763.

［5］ MOSHER C C，KEHO T H，WEGLEIN A B，et al.The impact of migration on AVO［J］.Geophysics，1996，61（6）：1603.

［6］ 王华忠，冯波，任浩然.二维offset平面波有限差分法叠前时间偏移［J］.石油物探，2009，48（1）：11.

［7］ BO FENG，HUAZHONG WANG，SHAOYONG LIU.3Doffset plane－wave finite－difference prestack time migration［J］.Expanded Abstracts of 79thAnnual International Meeting，SEG，2009：2934.

［8］ YILMAZ O.Seismic data analysis：Processing，inversion，and interpretation of seismic data［M］.Tulsa Society of Exploration Geophysics，2001.

［9］ MITCHELL A R，KELAMIS P G.Efficient tau－p hyperbolic velocity filtering［J］.Geophysics，1990，55（5）：619.

［10］ZHOU H，马建波.并行应用于3D地震偏移的Butterworth和Chebyshev倾角滤波器［J］.石油物探译丛，2000（5）：29.

［11］苑益军，牛滨华，王焕弟，等.去噪技术在地震资料处理中的应用［J］.东华理工学院学报，2005，48（1）：12

book=1,ebook=1

P 631.4＋43

A

10.3969/j.issn.1001－1749.2012.03.05

段心标（1982－），男，河南夏邑人，硕士主要从事地震成像方面的研究。

国家重大科技专项（2008ZX05014－001－002）

2011－10－19 改回日期：2011－11－25