随着我国经济的高速增长，石油消费量和需求也在急剧地增加，陆上老油区产量增长缓慢，所以，进行海洋油气资源勘探与开发尤为重要和迫切.海洋地震勘探是海洋地质与油气资源调查的主要手段，近年来海洋油气勘探的难度也越来越大，诸如渤海古潜山气藏、黄海中生代海相碳酸盐岩、南海崎岖海底等复杂地区的油气勘探对海洋地震资料成像处理精度提出了更高的要求.

传统偏移成像算法都是基于波场一次反射的假设，只利用了地面接收到的一次反射波，其他波场都被认为是噪声而尽可能地衰减掉.一直以来，大家都把注意力放在如何衰减多次波(Verschuur et al., 1992; 刘伊克等，2008; 石颖等，2013).事实上，任何接收到的反射波都是信号，尤其是在海洋地震资料中发育的强多次反射波蕴含着地下特殊构造的反射信息，花费高额成本采集到的多次波反射信号却因无法有效利用而被舍弃掉，无疑是非常可惜的.近年来，伴随着计算机技术的进步和发展，一系列关于如何利用多次波成像的方法成为研究热点(Verschuur, 2014)，受到专家学者的重视.Reiter等(1991)最初尝试了基于射线理论的方法对深水海底检波器接收的多次波数据进行叠前深度偏移，后期的发展主要是基于波动理论的多次波成像.

目前多次波成像方法主要分为三类：第一类是多次波降阶转化为准一次波偏移成像.基于共聚焦点(CFP)技术可以将多次波转变为一次反射波，聚焦变换后的一次波在零时间位置聚焦为一个点，而二阶多次波将转变为一次反射波(Berkhout and Verschuur, 2003).Berkhout和Verschuur(2006)结合自由表面相关多次波衰减(SRME)方法与多道加权互相关(MWCC)技术把多次波转变为一次反射波，将非线性的地震成像过程通过一系列线性步骤完成.Shan(2003)将包含多次波的记录与多次波互相关后得到准一次波，再按照传统一次波偏移方法基于双平方根偏移算子进行成像.第二类是低阶多次波作为高阶多次波震源实现多次波成像，该类偏移方法可以认为是基于面炮偏移的方法(Berkhout, 1992; Rietveld et al., 1992)，将所有可作为多次反射波的一次激发源作为面炮，按照传统偏移算子进行成像.Berkhout和Verschuur(1994)提出以包含多次波的记录为面炮、多次波为面炮记录的方式实现了基于波动方程的多次波成像，并基于“WRW”模型，理论上证明其正确性.Guitton(2002)指出，当以同时包含多次波与一次波数据为面炮记录时会出现很强的干扰噪声，影响对资料的构造解释.随着计算机技术的发展，基于全波方程的逆时偏移技术逐步得到了工业界的推广应用.Liu等(2011)提出了修改了传统逆时偏移初始上、下行波场实现逆时多次波偏移成像，并在Sigsbee2B模型测试上取得了明显的效果，类似的方法在海底地震(OBS)数据中也得到了同样应用(Zhang and Schuster, 2014).干涉噪声是多次波成像方法固有的难题，Wang等(2014a, 2014b)在局部角度域成像道集上，通过曲率差异来衰减干涉噪声，改善多次波成像质量.为了能够同时利用一次波多次波成像，郭书娟等(2011)提出了类似的基于广义的炮偏移方法实现地表多次波和一次波联合成像，Wang等(2014a, 2014b)提出一次波与多次波同时叠前逆时偏移成像.从计算效率和实用性方面考虑，基于单程波偏移算子多次波成像方法也得到了一定的发展，叶月明等(2014)尝试分别单独成像一次波与多次波，再将成像结果匹配叠加而避开子波难以估算的问题.Ye等(2013, 2014)将基于单程波算子的多次波成像技术应用到了海洋地震资料处理中，多次波成像在均衡照明与高分辨率方面展现了它的优势.第三类是基于反演的多次波成像方法，初期的最小二乘一次波与多次波联合偏移成像方法也是基于射线追踪的方法，关键点在于计算多次波的走时信息(He and Schuster, 2003).Brown和Guitton(2005)在地下局部变化缓慢的假设下，提出了共中心点(CMP)域进行的一次波与多次波最小二乘Kirchhoff偏移，这种方法需要预先从叠加剖面上识别出可能产生多次波的反射层位，再基于模型正演模拟出多次反射波，并将惩罚系数作用于干涉噪声来降低其对成像的影响, 在深水地震资料处理中得到了成功应用，Vyas和Brown(2006)修改了正则化算子，使这种方法能够在浅水地震资料中也适用.近年来，基于全波方程的多次波最小二乘偏移技术也得到了发展(Berkhout and Verschuur, 2011).Zhang和Schuster(2014)以预测出的自由表面多次波为上行波场，包含多次波的记录为下行波场实现了多次波的逆时最小二乘偏移，Wong等(2015)提出了基于最小二乘逆时偏移方法的一次波与多次波联合成像，通过最小二乘偏移的多次波迭代来压制噪声, 基于反演类的多次波成像方法具备较高的精度，但是非常耗时，而且非常依赖速度场的准确性.

本文属于第二类多次波成像方法，借鉴一次波与多次波同时叠前逆时偏移成像方法的思路，基于单程波偏移算子，提出了一种可衰减干涉假象的一次波与多次波同时成像方法，压制了一次波与多次波同时叠前偏移成像方法中产生的干涉假象，给出了主要干涉假象的表达式，并采取了“预测与匹配相减”的方式将干涉假象在成像域从一次波与多次波联合成像结果中衰减掉.该方法计算效率高、实用性强，同时利用了一次波与多次波的信息，并且衰减掉了强能量的干涉假象.Sigsbee2B模型测试证明了本方法的正确性，此外，在我国某深海探区的实际资料中得到应用，由于多次波提供了更多小角度入射信息与照明均衡的特点，所以资料成像的垂向分辨率得到了提高，基底照明更佳均衡，成像连续性得到了加强.

基于单程波算子的偏移方法相对于全波逆时偏移具有较高的计算效率，相对于射线类偏移方法具有较强的速度横向变化适应性，对于实际资料应用处理方面有一定的效果和效率的均衡优势，本研究所采用的都是基于单程波算子的偏移方法.图 1a所示的是一次波与多次波同时成像示意图，其中实线部分是自由表面相关多次波的传播路径.其中W是震源子波，下行传播至地下界面X位置处发生反射，反射波上行传播至地表后被记录为一次反射波P；一次反射波下行传播至地下界面X₁位置处发生反射，反射波上行传播至地表后被记录为一阶自由表面相关多次波M₁；一阶自由表面相关多次波下行传播至地下界面X₂位置处发生反射，反射波上行传播至地表后被记录为二阶自由表面相关多次波M₂，同样得到第N阶自由表面相关多次波M_N.利用一次反射波只能够得到X位置处的成像，而多次反射波携带着X₁至X_N位置的反射信息，但是，传统的一次波成像方法并不能够使得这些多次波正确归位.Liu等(2011)通过以一次反射波和多次波为面源，多次波为记录的方式实现多次波逆时偏移成像，得到X₁至X_N位置的正确成像(图 1a).为了能够同时进行一次波与多次波成像，Wang等(2014a, 2014b)以震源子波、一次波和多次波同时为下行波场，一次波和多次波为上行波场的方式进行逆时偏移成像，同时得到X，X₁至X_N位置的成像，但是，还产生了A₁，A₂等位置所示的干涉假象，该成像条件可以表示为：

图 1

Fig. 1

图 1 一次波与多次波成像方法原理图 (a)一次波与多次波同时成像；(b)正向延拓震源子波与逆向延拓多次波互相产生构造假象；(c)一次波与多次波对真实构造成像. Fig. 1 Schematic diagram showing principle of primary and multiple imaging (a) Simultaneous imaging of primaries and multiples; (b) Artifacts caused by cross correlation between forward propagated wavelets and backward propagated multiples; (c) Real structure imaging from primary and multiple.

(1)

其中x, z代表了空间位置，w表示频率，下标F代表了正向延拓，下标B代表了反向延拓，M代表了多次波，包含了1至N阶自由表面相关多次波：

(2)

将(2)式代入(1)式后得到

(3)

式(3)右边包含了两部分，第Ⅰ部分为正向延拓的面炮震源，第Ⅱ部分为反向延拓的面炮记录，展开后得到一次波与多次波同时成像更为具体的表达式：

(4)

(4) 式包含了四部分，第Ⅰ部分是真实构造成像，包含了正向延拓震源子波与逆向延拓一次反射波互相关成像、正向延拓一次反射波与逆向延拓一阶自由表面相关多次波互相关成像及N-1阶正向延拓自由表面相关多次波与N阶逆向延拓自由表面相关多次波互相关成像，其中同时含有一次波与多次波对真实构造成像的贡献，如图 1a中的X，X₁, …, X_N位置；第Ⅱ部分是由于正向延拓震源子波与逆向延拓不同阶次自由表面相关多次波互相关成像后产生的构造假象，如图 1b中的A₁，A₂等，其中正向延拓的震源子波与逆向延拓的一阶自由表面相关多次波产生的构造假象能量较强，是一次波与多次波同时成像中主要的影响因素；第Ⅲ部分是由正向延拓的一次反射波与逆向延拓的二阶及以上阶次多次波、正向延拓多次波与逆向延拓比正向延拓高两阶及以上多次波相关成像后产生的构造假象，该部分波场能量很弱，对整体构造成像并无明显影响；第Ⅳ部分的构成是一些背景噪声，也可以忽略不计.

由式(4)可知，影响一次波与多次波同时成像主要原因是正向延拓震源子波与逆向延拓自由表面相关多次波互相关成像后产生的构造假象，为了衰减这种假象，本文借鉴多次波压制的方法，采用预测与匹配相减的策略.在式(4)中的第Ⅱ部分干涉假象可以表示为

(5)

式(5)第Ⅰ部分是已知的子波，第Ⅱ部分是预测出的自由表面相关多次波，所以，以震源子波为震源波场正向延拓、以预测出的自由表面相关多次波记录波场逆向延拓，互相关成像后得到一次波与多次波同时成像方法中的干涉假象部分如图 1b所示(式4中的第Ⅱ部分).预测出干涉假象后，采用匹配相减的方法将式(5)得到的干涉假象从式(4)中减去，得到干涉假象衰减后的一次波与多次波同时成像Image(x, z)：

(6)

其中，f是基于传统维纳滤波的方法求取匹配因子，在此不做过多赘述，式(6)中的第Ⅲ项与第Ⅳ项可忽略不计，所以衰减干涉假象后的一次波与多次波成像条件可以近似表示为

(7)

式(7)只包含了一次波成像的贡献(第Ⅰ部分)和多次波成像的贡献(第Ⅱ部分)，得到如图 1c所示X，X₁，…，X_N位置处成像的同时，图 1b中的A₁，A₂，…等干涉假象也得到了较好的压制.

为了验证本文方法的有效性，首先对简单的三层模型进行了衰减干涉假象的一次波与多次波同时成像方法测试，最后对更为复杂的Sigsbee2B模型进行验证，取得了明显效果.

该三层模型的速度场如图 2a所示，第一层、第二层和第三层速度分别是1500 m·s^-1、2500 m·s^-1和3500 m·s^-1.采用高阶有限差分正演模拟得到单炮记录，中间放炮、151道接收、道间距20 m、3 s记录时间和4 ms采样，其中图 2b是包含多次波的单炮记录，P1和P2分别是第一层和第二层的一次反射波，M11和M21分别是第一层界面和第二层界面对应的一阶自由表面相关多次波，M12和M22分别是第一层界面和第二层界面对应的二阶自由表面相关多次波，如果接收时间足够长同样会接收到第一层界面和第二层界面对应的二阶以上自由表面相关多次波，分离多次波后得到包含直达波的一次反射波和多次波记录分别如图 2c和2d所示.

图 2

Fig. 2

图 2 三层模型与速度场 (a)三层模型速度场; (b)正演模拟单炮记录; (c)一次反射波; (d)分离出的多次波. Fig. 2 Schematic diagram of primary and multiple imaging on a three-layer model (a) Velocity model; (b) Forward modeling of single shot gathers; (c) Primaries; (d) Separated multiples.

图 3a是无干涉假象衰减的一次波与多次波同时成像结果，真实构造能够得到较好的归位，但是与此同时，在深层位置处出现两个平层(图 3a中箭头所指层位)，由原始模型可知，深部两个平层为构造假象.对于实际资料来说，真实地下构造形态是未知的，当出现类似构造时无法确定其真实性，给解释人员带来困难.图 3b是本文方法预测出的构造假象，除能量差异外(增益影响)，成像位置及构造形态与图 3a中的完全一致，匹配衰减后构造假象得到了较好的压制(图 3c), 与传统一次波成像(图 3d)相比，箭头所指构造假象得到较好压制的同时，椭圆形框内构造边界照明得到了明显加强，第二层成像更为均衡，成像边界下拉的不收敛现象有了一定的提高，这正体现出了多次波成像均衡照明度和高垂向分辨率的特点.

图 3

Fig. 3

图 3 简单模型偏移成像结果 (a)无干涉假象衰减的一次波与多次波同时成像; (b)预测出的干涉假象; (c)干涉假象衰减后的一次波与多次波同时成像; (d)传统一次波成像. Fig. 3 Migration imaging results of a simple model (a) Simultaneous imaging of primaries and multiples without artifacts suppressing; (b) Predicted artifacts; (c) Simultaneous imaging of primaries and multiples with artifacts suppressing; (d) Conventional imaging of primaries.

为了测试本方法对复杂模型的适用性，下面对国际上通用的Sigsbee2B模型进行试算.该模型的海底较深，最深处近3 km，海底下存在高速盐丘构造(图 4a)，海底与盐丘之间发育多个散射体，盐下有多个高陡断层构造.盐丘陡倾角边界由于初至反射波照明不足与盐丘的遮挡屏蔽作用使得成像照明不足.图 4b是含有多次波的记录，图 4c是一次反射波记录，图 4d是采用经典的SRME方法预测出的自由表面相关多次波，记录总时间长度为12 s，在7 s和11 s附近的强同相轴为海底反射的一阶与二阶自由表面相关多次波.图 5a是干涉假象衰减前的一次波与多次波同时成像，整体构造虽然得到了较好的刻画，但是其中箭头所指层位与整体构造趋势明显不一致，其形态与海底类似，这些是由于子波震源的正向延拓与逆向延拓的一阶自由表面相关多次波互相关成像而产生的干涉假象，严重影响了构造成像的真实性.图 5b是本文方法预测出的干涉假象，箭头所指的构造假象在形态上与图 5a中的完全一致，匹配衰减后该假象得到了明显压制(图 5c)，与传统一次波成像(图 5d)相比，一次波与多次波同时成像的高频信息更加丰富、垂向分辨率也更高.为了便于与一次波单独成像对比，选取局部范围(图 5c与图 5d中矩形框)对比，该区域包含了盐丘高陡边界、多个散射体构造和弱照明地层.图 6a是干涉假象衰减后的一次波与多次波同时成像局部放大，与传统一次波成像相比(图 6b)，其优势体现在三个方面，首先是盐丘侧翼高陡倾角边界得到了更好的成像(图 6中的长箭头)；其次是局部反射层成像照明更加均衡(图 6中椭圆形框内)；再次是盐丘上部散射体成像更加收敛(图 6中的短箭头).

图 4

Fig. 4

图 4 Sigsbee2B模型速度场和记录 (a) Sigsbee2B模型速度场；(b)包含多次波单炮记录；(c)自由表面相关多次波衰减后的单炮记录；(d)自由表面相关多次波. Fig. 4 Sigsbee2B velocity model and shot gather (a) Sigsbee2B velocity model; (b) Shot gather with surface related multiples; (c) Shot gather with surface-related multiple suppressing; (d) Surface related multiples.

图 5

Fig. 5

图 5 Sigsbee2B模型偏移 (a)无干涉假象衰减的一次波与多次波同时成像；(b)预测出的干涉假象；(c)干涉假象衰减后的一次波与多次波同时成像；(d)传统一次波成像. Fig. 5 Migration of Sigsbee2B model (a) Simultaneous imaging of primaries and multiples without artifact suppressing; (b) Predicted artifacts; (c) Simultaneous imaging of primaries and multiples with artifact suppressing; (d) Conventional imaging of primaries.

图 6

Fig. 6

图 6 Sigsbee2B模型偏移局部放大 (a)干涉假象衰减后的一次波与多次波同时成像; (b)传统一次波成像. Fig. 6 Local enlargement of Sigsbee2B model migration (a) Simultaneous imaging of primaries and multiples with artifacts suppressing; (b) Conventional imaging of primaries.

模型测试证明了本方法的正确性，为进一步验证对实际资料的处理效果，对某深海探区实际资料进行了试处理.该区位于我国深水水域，盆地内含油气构造发育，地层和构造条件具有良好的油气运聚和成藏体系，是我国深海油气勘探开发的重要目标区域.对于具备稳定海底反射的深海地震资料，海底鸣震不但表现得十分明显，而且近距离的强反射层之间会产生较强的层间多次波，该地区发育强深水海底反射的长周期自由表面相关多次波和部分非规律层间多次波，适于本方法的应用.

图 7a是该地区实际资料的深度域层速度场，海底深度在1 km以上，海水速度1500 m·s^-1，图 7b是含多次波单炮记录，其中，道距是12.5 m，每炮301道接收.图 7c是预测出的自由表面相关多次波，箭头所指为海底一阶与二阶自由表面相关多次波.以包含多次波记录与提取的稳定海底反射子波为起始正向延拓波场、包含多次波记录为起始反向延拓波场，干涉假象衰减前一次波与多次波同时成像如图 8a所示，箭头所示位置出现了连续性较强的同相轴，其构造形态与整体地层区域走势并不一致，图 8b是正向延拓震源子波与逆向延拓多次波相关成像产生的构造假象，匹配相减后构造假象得到了较好压制(图 8c)，与传统一次反射波成像(图 8d)相比，箭头所指浅层与深层成像都有了明显改善和提高，从频谱对比(图 8e)上可以看出，一次波与多次波同时成像后整个频带得以展宽，尤其是高频段信息更加丰富.图 9a与图 9b分别是干涉假象压制后的一次波、多次波同时成像与一次波单独成像的浅层局部放大，多次波参与成像后层位照明更加均衡.图 9c与图 9d分别是二者的深层成像对比，矩形框内的地层垂向分辨率得到提高，箭头所指基底构造成像连续性得到了显著改善.

图 7

Fig. 7

图 7 某深海探区速度场和单炮记录 (a)速度场; (b)含多次波单炮记录; (c)预测出的自由表面相关多次波. Fig. 7 Velocity model and shot gather in a real deep marine area (a) Velocity model; (b) Shot gather with surface related multiples; (c) Predicted surface related multiples.

图 8

Fig. 8

图 8 某深海探区实际资料处理 (a)无干涉假象衰减的一次波与多次波同时成像；(b)预测出的干涉假象；(c)干涉假象衰减后的一次波与多次波同时成像；(d)传统一次波成像; (e)一次波成像与一次波、多次波同时成像频谱对比. Fig. 8 Processing of deep marine area (a) Simultaneous imaging of primaries and multiples without artifacts attenuation; (b) Predicted artifacts; (c) Simultaneous imaging of primaries and multiples with artifacts attenuation; (d) Conventional imaging of primaries; (e) Frequency spectrum comparison between imaging of primaries and simultaneous imaging of primaries and multiples.

图 9

Fig. 9

图 9 某深海探区实际资料成像局部放大 (a)图 8c浅层成像局部放大; (b)图 8d浅层成像局部放大; (c)图 8c深层成像局部放大; (d)图 8d深层成像局部放大. Fig. 9 Local enlargement of deep marine area imaging (a) Shallow area imaging in Fig. 8c; (b) Shallow area imaging in Fig. 8d; (c) Deep area imaging in Fig. 8c; (b) Deep area imaging in Fig. 8d.

多次波普遍存在于海洋地震资料中，利用好多次反射波信息是提高海洋地震资料处理品质的重要手段，本文基于波动方程面炮偏移的思想，提出了可衰减干涉假象的一次波与多次波同时成像方法，发挥自由表面相关多次波对构造成像贡献的同时，压制了由于多次波的引入而产生的构造假象，并且避开了现有一次波与多次波同时成像方法中子波信息无法确定缺陷，Sigsbee2B模型测试证明了其正确性，我国某深海探区资料处理验证了其对海洋实际资料处理的实用性，取得了明显效果并得到以下结论：

(1) 多次反射波对地下的多次照明和较小入射角的特点使得其对成像有均衡照明度和高垂向分辨率的优势，具备稳定海底反射的深海探区资料更适合于自由表面相关多次波预测，有利于准确地预测出干涉假象，更适合于本方法的应用；

(2) 速度场的准确性和实际采集资料的信噪比是影响本方法的重要因素，当速度场不准确或海底起伏较大的情况下会产生较多的干涉噪声或假象难以衰减，在这种情况下本方法的使用需要慎重；

(3) 多次波成像的研究与应用是对一次波成像的有利补充，但并不是代替一次波成像，基于反演类技术的多次波成像方法将通过多次迭代来消除干涉假象的影响，也是多次波成像技术的下一步研究方向.