逆时偏移剖面假象产生机制及其消除方法

康智清 （中石化河南油田地球物理勘探公司物探研究所，河南 南阳473132）

刘恩 （中石化中原油田地质录井公司，河南 濮阳457000）

司杰戈 （中石化胜利油田物探公司，山东 东营257000）

陈康 （中国石油大学 （华东）地球科学与技术学院，山东 青岛266555）

逆时偏移是地震偏移方法的重要发展。与传统偏移方法不同，逆时偏移是在时间轴上实现外推，可以看作是沿时间反方向的正演模拟过程。传统的沿深度方向的偏移方法基于单程波方程，而逆时偏移则基于全波方程，允许地震波在全方位传播，因而不存在倾角限制。随着计算机技术的发展和对复杂构造成像的更高要求，逆时偏移技术的研究也渐渐深入，由以前的二维发展到三维、声波发展到弹性波、各向同性发展到各向异性。但是偏移假象仍然是影响逆时偏移剖面效果的重要方面［1～5］。

整体看来，逆时偏移需要解决的问题主要包括几方面。①正演方面的问题：包括波动方程的改进；正演模拟的精度；边界吸收处理；数值频散；稳定性问题以及弹性波的波场分离问题。②逆时延拓过程中的问题：逆时延拓中产生的次生干扰；反射和透射损失问题。③成像条件的问题：如何获得正确的保幅的成像条件；如何通过成像条件消除各种偏移假象。④偏移假象问题：假象产生的原因和机制；直达波、边界反射、回转波等的影响造成的偏移假象。⑤算法的效率问题：通过改进算法提高运算效率；通过GPU和并行运算的方法提高运算效率。笔者主要研究第④方面的问题，关于逆时偏移去噪的研究，主要包括修改波动方程、修改成像条件、成像后进行滤波处理等。具体地说，Valenciano等［6］根据反演理论提出了反褶积成像条件。Mulder等［5］运用一个空间域的低通滤波器来消除噪声。Yoon等［3］提出了在零延迟互相关条件中加入Poynting矢量来消除成像噪声。Liu等［4］把全波场分解成单程波分量，并运用成像条件去结合这些波场分量达到消除成像噪声的目的。还有方向性衰减、速度平滑、波场分离成像、无反射波动方程逆时偏移成像等方法［7～12］。

笔者从偏移假象产生的机制出发，重点分析了逆时偏移假象产生的各种原因，针对这些原因提出和总结了解决的方法。特别针对成像条件引入的低频噪声，笔者采用了拉普拉斯滤波算法，通过模型测试结果可知该算法效果较好。

1 逆时偏移的基本原理

逆时偏移主要包括波动方程的正演模拟、逆时外推和成像条件的确定3个步骤［13］，在进行正演模拟和逆时外推过程中分别保存其波场，然后运用成像条件进行求和，得到局部成像数据体。最后将所有炮集的逆时偏移结果进行叠加，得到最终的叠前深度偏移成像结果。三维逆时偏移的定解问题可以描述为：

利用微分和差分关系，对式（1）进行差分离散，得到用于正演模拟和逆时深度偏移的高阶差分方程的初始方程：

式中：u为不同时刻的波长值，m；x、y、z为三维坐标；v为速度，m／s；t为时间，s；Δt为时间采样间隔，s；O为高阶项。式 （2）代表正演过程，式 （3）代表逆时外推过程。将得到的正演波场和逆时外推波场进行互相关就可以得到成像数据体，最后进行叠加就可以得到成像剖面。

2 成像前噪声的产生机制和消除方法

2．1 直达波的影响

在逆时偏移第一个阶段，炮点波场在正向延拓过程中会产生直达波，如果将含有直达波的波场进行互相关就会在成像剖面顶端产生很严重的偏移假象。因此，必须在正演记录中将直达波切除，消除直达波的影响。图1（a）为速度模型，模型大小为200（道）×200（采样点），网格大小为10m×10m，在第100道激发，各道接收，可以得到图1（b），此时为没有去除直达波的正演记录。用图1（b）的正演记录进行逆时外推，最终进行互相关成像就可以得到图1（c），可以看到，在成像结果顶端产生偏移假象。图1（d）为去除了直达波的正演记录。图1（e）为对应的单炮逆时偏移结果，可以看出，去除直达波以后，顶端的偏移假象得到了消除。

图1 直达波对逆时偏移成像剖面的影响

2．2 边界反射的影响

在正演的过程中，边界的处理是很重要的，如果边界问题处理不好，将会引入边界反射，边界反射会将一些干扰波引入到正演记录中，最终直接导致偏移剖面中出现偏移假象问题。图2（a）为层状速度模型，图2（b）代表含有边界反射的正演记录。图2（c）代表含边界反射的单炮偏移剖面，可以明显地看到，在偏移剖面中，存在另外一些同相轴，即是边界反射带来的偏移假象。笔者采用了PML（perfect matched layer） 吸 收 边 界 条件消除边界反射，通过边界反射的消除，图2（d）中的正演记录不含边界反射，图2（e）为消除边界反射后对应的单炮偏移剖面。对比图2（e）和（c）可知，通过边界条件的使用，偏移剖面中由于边界反射引起的偏移假象得到了很好的压制。

图2 直达波对成像剖面的影响

3 成像过程中噪声的产生机制和消除方法

3．1 成像条件引入噪声的机制

目前，常用的成像条件有零时刻成像条件、互相关成像条件、波阻抗成像条件［14］。在逆时偏移成像中，通常使用的是互相关成像条件［15］。由于一些不正确的互相关，导致了成像剖面中存在着偏移噪声，这些噪声主要分布在成像剖面浅部，遮盖地下真实形态，从而降低了成像质量。下面具体以层状模型为例分析噪声的形成机制。

图3（a）为t＝0．2s时刻的炮点波场，波是向下传播的。图3（b）为对应于图3（a）的逆时波场，在对应的点波传播方向是向上的，图3（a）和 （b）的互相关就形成了图3（c）中的一道弧线，即为成像噪声。图3（d）为t＝0．4s时刻的炮点波场，波场存在向下传播的透射波和向上传播的反射波；图3（e）为同时刻的逆时外推波场，波场中也存在逆时入射波和逆时反射波，图3（d）和 （e）的互相关就构成了图3（f）的A、B两点和一条弧线，其中A、B两点是正确的像点，而弧形是互相关引入的噪声。

3．2 拉普拉斯算子去噪原理

对于成像过程中的噪声，目前已经由很多的学者提出很多去除的方法，主要有方向性衰减、速度平滑、波场分离成像、无反射波动方程逆时偏移成像等。但是方向性衰减、波场分离成像、波因廷矢量成像条件等在实现上有很大的难度。无反射波动方程效果不是太好，特别是对于大角度入射时反射还是比较明显。速度平滑是一种比较简单的方法，效果也不错，可是对精确的速度模型来说不是一个好的方法，因为速度模型的平滑引入了速度误差。因此，笔者采用拉普拉斯滤波算法消除成像过程中的噪声，它在去噪的过程中既起到了高通滤波的作用，也有速度平滑的效果。

对于常规的二维拉普拉斯滤波算子：

式中：为拉普拉斯算子；I为波场值。从式（4）中可以看出，拉普拉斯算子可以表示为一个二阶微分的形式。

图3 成像条件引入噪声的形成机制

图4 为不同阶数拉普拉斯算子对低频的压制效果，可以看到，微分有提升高频的作用，同时可以压制低频，且随着阶数越高，微分作用对低频的削弱越明显，对高频成分也是非线性的提升，2阶微分对高频的提升和对低频的压制作用明显强于1阶微分。通过图4可以得出下面的结论：拉普拉斯算子是基于2阶微分的组合算子，具有很好的压制低频和提升高频的作用，把一个信号或者图像通过拉普拉斯算子，其过程相当于一个高通滤波器，相比高通滤波器，拉普拉斯算子滤波还能起到平滑的效果。

图5（a）是一个合成地震记录的深度域表示，由一个倾斜轴和一个水平轴组成，Kz表示Z方向的波数域，Kx表示X方向的波数域。图5（b）是将图5（a）通过拉普拉斯算子后的结果，可以发现，通过拉普拉斯算子，低频成分得到了很好的压制，同时提升和保留了高频成分。

图4 不同阶数拉普拉斯算子对低频的压制效果

图5 滤波前、后波数域分析

3．3 模型测试

图6 、7是对拉普拉斯滤波的效果测试。图6（a）和图7（a）分别是速度模型和sigbee模型，图6（b）是30炮逆时偏移剖面，可以看到成像噪声分布在剖面各个地方，严重影响了剖面的质量；同理可以看到图7（b）剖面中的噪声也是非常严重。图6（c）和图7（c）是将含噪剖面进行拉普拉斯算子滤波的结果，可以看出，通过滤波，成像剖面的质量都得到了很好的改善，成像噪声得到了很好的压制。

图6 速度模型拉普拉斯算子滤波效果

图7 sigbee模型拉普拉斯算子滤波效果

4 结 语

笔者分成像前和成像过程两个阶段分析了逆时偏移剖面中的偏移假象的原因，成像前主要有直达波和边界反射的影响，成像过程主要是成像条件的影响。针对这些问题，笔者提出了有效的解决方法，在成像前分别用切除直达波和采用PML边界条件的方法达到目的；对于成像条件引入的噪声，笔者采用了拉普拉斯算法进行压制，通过模型测试结果可以看出，拉普拉斯算法能够有效地对成像条件引入的噪声进行压制。

［1］Baysal E，Kosloff D D，Sherwood J W C．Reverse－time migration ［J］．Geophysics，1984，48 （4）：1514～1524．

［2］Chang W F，McMechan G A．Elastic reverse－time migration ［J］．Geophysics，1987，52 （4）：1365～1378．

［3］Yoon K，Marfurt K J．Reverse－time migration using the Poynting vector．Exploration Geophysics，2006，37 （1）：102～107．

［4］Liu Faqi，Zhang，Guanquan，Morton S A，et al．An effective imaging condition for reverse－time migration using wavefield decomposition ［J］．Geophysics，2011，76 （1）：S29～S39．

［5］Mulder W A，Plessix R E．A comparison between one way and two way wave equation migtation ［J］．Geophysics，2004，69 （6）：1491～1504．

［6］Valenciano A A，Biondi B．Deconvolution imaging condition for reverse－time migration ［J］．Stanford Exploration Project，2002，Report 112：83～96．

［7］Liu Faqi，Zhang G，Morton S A，et al．Reverse－time migration using one－way wavefield imaging condition ［A］．2007SEG Annual Meeting ［C］．San Antonio，2007－09－23～28．

［8］Guitton A，Kaelin B．Least－squares attenuation of reverse－time migration artifact ［J］．Geophysics，2004，69 （6）：1491～1504．

［9］Rickett J，Sava P．Offset and angle－domain common image－point gathersfor shot－profile migration ［J］．Geophysics，2002，67 （7）：883～889．

［10］Fletcher R P，Fowler P J．Suppressing unwanted internal reflections in prestack reverse－time migration ［J］．Geophysics，2006，71（6）：79～82．

［11］红伟，刘洪，邹振 ．地震叠前逆时偏移中的去噪与存储 ［J］．地球物理学报，2010，53（9）：2171～2180．

［12］张军华，吕宁，田连玉，等 ．地震资料去噪方法技术综述评述 ［J］．地球物理学进展，2006，21（2）：546～553．

［13］杨仁虎，常旭，刘伊克 ．叠前逆时偏移影响因素分析 ［J］．地球物理学报，2010，53（8）：1902～1913．

［14］陈可洋 ．基于高阶有限差分的波动方程叠前逆时偏移方法 ［J］．石油物探，2010，49（2）：153～157．

［15］丁亮，刘洋 ．逆时偏移成像技术研究进展 ［J］．地球物理学进展，2011，26（3）：1085～1100．