在勘探地震学研究领域，地震偏移成像技术是获取地下构造形态最有效的手段之一.地震偏移成像分为叠加前偏移和叠加后偏移两种类型，根据成像域的不同还可分为时间域偏移和深度域偏移.在叠前深度域偏移方法中，目前石油地震勘探领域一般采用基于射线理论的波动方程积分解法和基于波动理论的微分波动方程单程波解法，这些常规的方法精度有限，而且振幅和相位不准，难以处理横向速度变化剧烈和高陡倾角构造的成像问题^{[1, 2]}.近年来发展起来的基于波动理论的逆时偏移方法是基于双程波算子，精度高，相位准确，不受横向变速和高陡倾角影响，并且可以利用回转波来成像，解决了常规偏移所面临的成像难题^[3~8].逆时偏移方法也分为叠加后和叠加前两类方法，叠前逆时偏移成像原理是基于1971年Claerbout提出的反射地震图成像原理演变而来^[9].假设源波场代表下行波场，接收波场代表上行波场，即可得出叠前逆时偏移零延迟互相关成像条件^{[4, 10~12]}

(1)

其中P_s（X，t）表示源波场，P_g（X，t）表示接收波场.叠前逆时偏移成像效果受到许多因素影响，其中主要的影响因素有：计算速度和存储空间的节省、初始速度模型的建立、震源子波的选择、数值模型边界条件定义、假像的消除等等.对于计算速度和存储量大的问题，随着计算机硬件的快速发展，将会不断得到改善，同时可以采取一些计算技术和存储策略来加以缓解^[13~17].初始速度模型越准确，偏移成像效果越好.震源子波影响偏移成像相位和成像分辨率，子波选择合适，能增加相位准确度和提高成像分辨率.选择一个合适的模型吸收边界，能减弱边界反射对偏移成像的影响，从而可以提高成像质量.对于假像问题，主要由首波、反散射波和回折波等造成，这类波严重干扰偏移成像质量.对于直达波造成的地表假像可以通过切除直达波来消除；对于反散射波和回折波等造成的传播路径上的假像消除，许多学者进行了尝试性的研究和探索^{[4, 18~20]}.本文主要针对以上这些问题进行研究分析.

这里主要分析初始速度模型、震源子波、模型边界、直达波和反射波对叠前逆时偏移成像精度的影响.下面选用一个简单的两层水平层状介质模型进行试验和分析，如图 1.模型大小300×300（网格），网格大小10m×10m，上层介质速度4000m/s，下层介质速度5000m/s，记录长度1s，采样点500，采样间隔2ms，震源采用Ricker子波，主频30Hz，位于模型表面正中央，采用双边接收，检波器遍布整个模型表面.

图 1

Fig. 1

图 1 两层水平层状介质模型 Fig. 1 The two horizontal layered model

叠前逆时偏移成像结果依赖于初始速度模型，初始速度模型越准确，偏移成像结果越好.这里分析了初始速度大小和界面位置存在偏差时的叠前逆时偏移成像结果（白色线条代表界面成像准确位置）.图 2是初始速度大小存在偏差时的叠前逆时偏移成像结果.图 2（a1）和（a2）分别是初始速度比真实速度低2%和5%时的叠前逆时偏移成像结果，从图中可以看出，初始速度偏低时，成像位置上移，同时传播路径上的假像不断减弱.图 2（b1）和（b2）分别是初始速度比真实速度高2%和5%时的叠前逆时偏移成像结果，从图中可以看出，初始速度偏高时，成像位置下移，同时传播路径上的假像不断扩散到边缘.图 3是初始速度模型界面存在偏差及对应的叠前逆时偏移成像结果.图 3（a1）、（a2）和（a3）分别是界面位置上移3.33%、6.67%和16.67%时的初始速度模型，这等效于模型上半部的速度增大.图 3（b1）、（b2）和（b3）分别是对应图 3（a1）、（a2）和（a3）的叠前逆时偏移成像结果，从图中可以看出，随着初始速度模型界面不断上移，成像位置不断下移，同时传播路径上的假像不断扩散.图 3（a4）、（a5）和（a6）分别是界面位置下移3.33%、6.67%和16.67%时的初始速度模型，此时模型上半部的速度没有改变.图 3（b4）、（b5）和（b6）分别是对应图 3（a4）、（a5）和（a6）的叠前逆时偏移成像结果，从图中可以看出，随着初始速度模型界面不断下移，成像位置影响较小，同时传播路径上的假像也逐渐减弱.

图 2

Fig. 2

图 2 叠前逆时偏移成像结果(速度有偏差) (1)叠前逆时偏移成像结果(速度低2%); (b1)叠前逆时偏移成像结果(速度高2%); (2)叠前逆时偏移成像结果(速度低5%); (b2)叠前逆时偏移成像结果(速度高5%). Fig. 2 The pre-stack RTM imaging results (velocity with error) (a1) The pre-stack RTM imaging result (lower velocity with error 2%); (b1) The pre-stack RTM imaging result (higher velocity with error 2%); a2) The pre-stack RTM imaging result (lower velocity with error 5%); (b2) The pre-stack RTM imaging result (higher velocity with error 5%).

图 3

Fig. 3

图 3 速度模型和相应的叠前逆时偏移成像结果(界面有偏差) (1)界面上移3.33%的速度模型; (b1)叠前逆时偏移成像结果(对应a1); (2)界面上移6.67%的速度模型; (b2)叠前逆时偏移成像结果(对应a2); (3)界面上移16.67%的速度模型; (b3)叠前逆时偏移成像结果(对应a3); (4)界面下移3.33%的速度模型; (b4)叠前逆时偏移成像结果(对应a4); (5)界面下移6.67%的速度模型; (b5)叠前逆时偏移成像结果(对应a5); (6)界面下移16.67%的速度模型; b6)叠前逆时偏移成像结果(对应a6). Fig. 3 Velocity model and the corresponding pre-stack RTM imaging results (interface with error) (a1) Velocity model (interface upward with error 3.33%); (b1) The pre-stack RTM imaging result (correspond to a1); (a2) Velocity model (interface upward with error 6.67%); (b2) The pre-stack RTM imaging result (correspond to a2); (a3) Velocity model (interface upward with error 16.67%); (b3) The pre-stack RTM imaging result (correspond to a3); (a4) Velocity model (interface downward with error 3.33%); (b4) The pre-stack RTM imaging result (correspond to a4); (a5) Velocity model (interface downward with error 6.67%), (b5) The pre-stack RTM imaging result (correspond to a5); (a6) Velocity model (interface downward with error 16.67%); (b6) The pre-stack RTM imaging result (correspond to a6).

震源子波对叠前逆时偏移成像具有一定的影响，它影响成像相位和成像分辨率.这里选取Ricker子波，并分析了子波相位对偏移成像结果的影响.图 4（a1）、（a2）、（a3）、（a4）和（a5）分别是主频为10、20、30、40Hz和50Hz的Ricker子波，子波主瓣宽度由宽变窄，同时子波相位不断靠近起始点.图 4（b1）、（b2）、（b3）、（b4）和（b5）分别是对应图 5（a1）、（a2）、（a3）、（a4）和（a5）时的叠前逆时偏移成像结果，从图中可以看出，在主频为30Hz的时候，成像位置最准确，因为此时子波与理论地震记录所用子波相同，相位一致.而主频改变时，子波相位偏离真实相位（相对固定起始点），从而引起成像位置偏离真实位置.因此，震源子波越接近地震记录子波，相应的叠前逆时偏移成像位置越准确.

图 4

Fig. 4

图 4 不同主频震源子波和相应的叠前逆时偏移成像结果 (1)震源子波(主频10Hz); (b1)叠前逆时偏移成像结果(对应a1);a2)震源子波(主频20Hz); (b1)叠前逆时偏移成像结果(对应a2); (3)震源子波(主频30 Hz); (b1)叠前逆时偏移成像结果(对应a3); (4)震源子波(主频40Hz); (b1)叠前逆时偏移成像结果(对应a4); (a5)震源子波(王频50 Hz); (b1)叠前逆时偏移成像结果(对应a5). Fig. 4 Source wavelet with different dominant frequency and the corresponding pre-stack RTM imaging results (a1) Source wavelet with 10 Hz dominant frequency; (b1) The pre-stack RTM imaging resutt (correspond to a1); (a2) Source wavelet with 20 Hz dominant frequency; (b2) The pre-stack RTM imaging resutt (correspond to a2).(a3) Source wavelet with 30 Hz dominant frequency, (b3) the pre-stack RTM imaging result (correspond to a3); (a4) Source wavelet with 40 Hz dominant frequency; (b4) The pre-stackRTM imaging result (correspond to a 4); (a5) Source wavelet with 5 0 Hz dominant frequency; (b5) The pre-stack RTM imaging resutt (correspond to a5).

地震勘探中实际产生的地震波是在地下半无限空间介质中传播，但在进行数值模拟时，由于模型大小的限制，人为地引入了边界.地震波在这些边界上会产生反射，对波场成像产生一定的干扰.这里进行了一个简单试验，如图 5（a1）是边界吸收不太好的地震记录，图 5（b1）是相应的叠前逆时偏移成像结果，从图中可以看出，成像位置准确，除了传播路径上反射波造成的假像外，在箭头所指处出现了边界反射造成的假像，但和有效振幅相比，能量比较微弱.图 5（a2）是边界吸收较好时的地震记录，图 5（b2）是相应的叠前逆时偏移成像结果，从图中可以看出，成像位置准确，除了传播路径上反射波造成的假像外，没有边界造成的假像.因此在选择吸收边界时，尽量选择简单易用、计算效率高、并具有一定吸收效果的吸收边界，就可以达到较好的叠前逆时偏移成像效果.

图 5

Fig. 5

图 5 地震记录和相应的叠前逆时偏移成像结果(不同程度边界吸收效果) (1)地震记录(边界吸收不太好); b1)叠前逆时偏移成像结果(对应a1);(2)地震记录(边界吸收较好); b2)叠前逆时偏移成像结果(对应a2). Fig. 5 Seismic records and the corresponding pre-stack RTM imaging results (with absorption boundary condition at different absorption levels) (a1) Seismic record (at bad absorption level); (b1) The pre-stack RTM imaging resutt (correspond to a1); (a2) Seismic record (at good absorption level); (b2) The pre-stack RTM imaging resutt (correspond to a2).

在叠前逆时偏移方法中，直达波形成的假像较容易消除.图 6是直达波切除前后地震记录与对应的叠前逆时偏移成像.图 6（a1）是含有直达波的地震记录，图 6（b1）是相应的叠前逆时偏移成像，从图中可以看出，成像位置准确，除了传播路径上的假像，在地表也出现了严重的假像.图 6（a2）是切除直达波的地震记录，图 6（b2）是相应的叠前逆时偏移成像，从图中可以看出，成像位置准确，除了传播路径上的假像，地表处严重的假像也基本消除.因此可以通过切除直达波来消除地表浅层假像.

图 6

Fig. 6

图 6 地震记录和相应的叠前逆时偏移成像结果 (a1)地震记录(含直达波); (b1)叠前逆时偏移成像结果(对应a1); (2)地震记录(切除直达波); (b2)叠前逆时偏移成像结果(对应a2). Fig. 6 Seismic records and the corresponding pre-stack RTM imaging results (1) Seismic record (with direct wave); (b1) The pre-stack RTM imaging resutt (correspond to a1); (a2) Seismic record (without direct wave); (b2) The pre-stack RTM imaging resutt (correspond to a2).

从图 6（b2）可知，直达波形成的假像消除后，仍然可以看到成像结果中存在较严重的假像，这些假像是在波场延拓过程中波传播的路径上产生的.本研究提出了一种振幅补偿滤波方法，用来消除这种假像.

图 7a是反射波引起的传播路径上的假像，频率较低，主要均匀分布在波的传播路径上.如何消除假像，许多学者进行了一些研究和探索，其中滤波是一种较好的方法，能适应任何复杂介质.但直接滤波，往往会损失有效振幅.如何有效消除假像且尽可能不损失有效振幅呢？本文给出了振幅补偿滤波方法，即在滤波前后，尽可能保持有效振幅，这实质上是保结构思想的延伸.主要思路是构建一个滤波算子L和新成像条件I°（X），使得滤波算子L不仅能够消除这些假像，而且能够保持有效振幅.用数学公式表述如下：

图 7

Fig. 7

图 7 叠前逆时偏移成像结果 (a)叠前逆时偏移成像结果(未滤波); b)叠前逆时偏移成像结果(振幅补偿滤波). Fig. 7 The pre-stack RTM imaging results (a) The pre-stack RTM imaging resutt (without filter); (b) The pre-stack RTM resutt (with amplitude compensation filtering).

定义新成像条件为：

(2)

使得

(3)

式（3）表明了在滤波前后保持了有效振幅，此即振幅补偿滤波方法.其中P°s是振幅补偿后的源波场，当然也可以对接收波场进行振幅补偿.新的成像条件依赖于滤波算子，对于简单算子，可以求出新成像条件的解析形式.对于复杂算子而言，很难得到新成像条件的解析形式，可以通过其他手段得到近似结果.这里对滤波算子L=c²∇²，通过推导求出新成像条件为

(4)

具体推导过程参见文献[21].

图 7b是应用振幅补偿滤波方法得到的叠前逆时偏移成像结果，从图中可以看出，不仅传播路径上的假像得到消除，而且有效振幅保持得较好.

在上文所述研究工作的基础上，应用本文提出的振幅补偿滤波方法，对勘探地球物理学界给出的SEG/EAGE二维盐丘模型、Marmousi模型和本研究设计的崎岖海底模型进行了叠前逆时偏移成像.

图 8是SEG/EAGE二维盐丘模型和相应的叠前逆时偏移成像结果.图 8a由于盐丘体形状复杂，造成地层介质横向和纵向变速很大.除盐丘体外，主要由一些平缓和过压地层以及高陡断层组成^[22].模型大小649×150（网格），网格大小24m×24m，介质速度1524~4480.56m/s，记录长度5s，采样点626，采样间隔8ms，震源采用Ricker子波，共325炮，炮点位置从模型表面最左边向右移动，炮间距48m，道间距24m，每炮176道，采用右边放炮左边接收观测系统.图 8b是没有经过滤波的叠前逆时偏移成像结果，从图中可以看出，由于假像的严重干扰，造成偏移成像结果非常不好，地层分界面成像模糊不清，基本看不出来.图 8c是应用本文中给出的振幅补偿滤波方法得到的叠前逆时偏移成像结果，从图中可以看出，成像结果较好，成像位置准确，特别是高陡倾角断层也获得了较好的成像.图 8d是图 8c的盐下成像放大部分，从图中可以看出，盐下高陡倾角断层也获得了一定的成像.图 8e是压缩震源子波后的叠前逆时偏移成像结果，图 8f是图 8e的盐下成像放大部分，从图中可以看出，通过压缩震源子波，可以提高成像分辨率，同时盐下可以获得较好的偏移成像，甚至盐下的断层也能显现出来，而常规偏移方法很难得到盐下的清晰成像.当然压缩震源子波后造成了一些频散现象，可以通过细分网格或提高计算精度来压制.

图 8

Fig. 8

图 8 SEG/EAGE二维盐丘模型和叠前逆时偏移成像结果 (a) SEG/EAGE二维盐丘模型; (b)叠前逆时偏移成像(未滤波); (c)叠前逆时偏移成像(振幅补偿滤波); (d)叠前逆时偏移成像(对应图 8c的盐下部分); (e)叠前逆时偏移成像(压缩震源子波); (f)叠前逆时偏移成像(对应图 8e的盐下部分). Fig. 8 SEG/EAGE 2-D satt model and the corresponding Pre-stack RTM imaging results (a) SEG/EAGE 2-D satt model; (b) The Pre-stack RTM imaging resutt (without filtering); (c) The Pre-stack RTM imaging resutt (with amplitude compensation filtering); (d) The Pre-stack RTM imaging resutt (enlarged subsatt part corresponding to c); (e) The Pre-stack RTM imaging resutt (with suppressing source wavelet); (f) The Pre-stack RTM imaging resutt (enlarged subsatt part corresponding to e).

图 9a是Marmousi模型，它主要由一些倾斜地层、高角度逆冲断层、角度不整合地层以及地层隆起构成.模型大小737×750（网格），网格大小12.5m×4m，介质速度1500~5500m/s，记录长度3s，采样点750，采样间隔4ms，震源采用Ricker子波，共240炮，炮点位置从模型左边（240，0）处右移，炮间距25m，道间距25m，每炮96道，采用右边放炮左边接收观测系统.图 9b是通过振幅补偿滤波后的叠前逆时偏移成像结果，从图中可以看出，成像结果较好，成像位置准确，构造清晰，振幅对应很好.

图 9

Fig. 9

图 9 Marmousi模型和叠前逆时偏移成像结果 OMarmous模型; (b)叠前逆时偏移成像(振幅补偿滤波). Fig. 9 Marmousi model and the corresponding Pre-stack RTM imaging results (a) Marmouti model; (b) The Pre-RTM imaging resutt (with amplitude compensation filtering).

图 10a是本研究设计的崎岖海底模型，主要由水平海面、崎岖海底、水平地层以及渐变地层组成.模型大小900×400（网格），网格大小12m×12m，介质速度1500~3700m/s，记录长度6s，采样点750，采样间隔8ms，震源采用Ricker子波，共450炮，炮点位置从模型最左边右移，炮间距24m，道间距24m，每炮120道，采用右边放炮左边接收观测系统.图 10b是通过振幅补偿滤波后的叠前逆时偏移成像结果，从图中可以看出，成像结果较好，成像位置准确，构造清晰，振幅对应很好.

图 10

Fig. 10

图 10 本研究设计的崎岖海底模型和叠前逆时偏移成像结果 (a)本研究设计的崎岖海底模型; (b)叠前逆时偏移成像(振幅补偿滤波). Fig. 10 Rugged water-bottom model designed by us and the corresponding Pre-stack RTM image results (a) rugged water-bottom model designed by us, (b) the Pre-stack RTM imaging resutt (with amplitude compensation filtering).

影响叠前逆时偏移成像的因素很多，本文只分析了其中几个影响因素.对于计算和存储问题，我们采取在有限孔径上进行计算和存储，在一定程度上减少了计算和存储量，但对于太大的模型来说，计算和存储还是一个问题.Clapp于2009年提出在叠前逆时偏移中采用随机边界的方法^[17].先让源波场传播到最大时刻，再进行逆传播，使之与接收波场同步，不用存储全部波场，缓解了波场存储问题，但这种方法需要全模型计算，大大增加了计算量，并且只适用于无耗散介质中的声波传播.复杂介质叠前逆时偏移计算和存储问题还需要进一步研究.还有频散问题，它主要和地震波速度、时间步长、空间步长、震源子波、数值计算精度等有关，其中地震波速度、时间步长和空间步长组成一个Courant数，它不仅影响计算稳定性，而且影响频散，合理搭配这些因素，可以有效压制频散，从而提高成像质量.同时在互相关成像的同时，在一定程度上也压制了频散效应.关于数值计算方法，本文采用交错网格有限差分方法，这种方法比传统的中心差分格式计算精度要高.相对有限差分方法而言，伪谱法数值计算频散很小，但计算量太大.如何权衡时间步长、空间步长、震源子波、数值计算精度以及计算效率、频散问题等，需要进一步研究和分析.

本文主要分析和探讨了叠前逆时偏移成像的几个主要的影响成像精度的因素，针对反射波造成的假像，给出了振幅补偿滤波方法.并对SEG/EAGE二维盐丘模型、Marmousi模型和崎岖海底模型进行了叠前逆时偏移成像，均取得了较好的成像效果.