在海上数据中，多次波有可能是充分发育的.在过去的研究过程中，一般将多次波视为噪音，试图将其去除(Verschuur et al., 1992; Weglein et al., 1997)，但是在成像过程中也可以利用多次波信息(Liu et al., 2011；郭书娟等，2011；叶月明等，2014；刘伊克等，2015).在他们提出的多次波成像方法中，往往需要复杂的预测多次波过程，为此Wang等(2014a)提出了不需要预测多次波的同时利用一次波和多次波的偏移方法.在此基础上，为了避免繁杂的子波预测，Wang等(2014b)和郑忆康等(2015)提出了直接在偏移中利用多次波的数据自相关偏移方法.在这种多次波偏移过程中，上行波场和下行波场延拓的都是含有多次波的地震数据.多次波偏移结果中的假象主要来自于不同地震事件之间的互相关，由于这种互相关是满足成像原理的，导致在成像结果中很难将其直接去除.但是在角度域共成像点道集(angle domain common image gathers, ADCIGs)中, 假象的道集是不平的，可以很好地加以区分.Fei等(2015)和王一博等(2016)指出双程波算子在偏移速度中存在剧烈变化的速度梯度时，会产生偏移假象，因而在数值实现上本文采用的是单程波偏移算法.

Sava和Guitton (2005)提出了采用抛物线Radon变换在角度道集内压制一次波偏移中的残余多次波成像噪音.多次波偏移中的部分成像假象与一次波偏移中的残余多次波成像噪音具有相同的曲率，因此也可以采用抛物线型Radon变换在角道集内压制多次波偏移假象.我们的压制假象方法主要分为以下几步：(1) 在多次波偏移过程中直接获取ADCIGs；(2) 对获得的ADCIGs应用高分辨率抛物线型Radon变换，对于正确的成像事件，其能量应集中在零值附近；(3) 切除Radon域中非零区域的成分; (4) 应用Radon逆变换，获得假象被压制的ADCIGs；(5) 叠加不同角度的ADCIGs，得到最终的成像结果.

为了获得高精度的分离结果，常规Radon变换不能满足高分辨处理的需求，为此本文采用了高精度的抛物线型Radon变换.Zhou和Greenhalgh (1994)将抛物线型Radon变化在频率空间域实现并用最小二乘方法求解，避免了大矩阵求逆并且获得得了高精度解，但是这种解存在低频段的稳定性问题. Sacchi和Ulrych (1995)提出了一种高精度的Radon变换算法，用迭代逼近的方法求解大型线性稀疏算子.这种算法更为稳定，更适合应用到ADCIGs去假象过程中.经过高精度的抛物线型Radon变化，ADCIGs的假象能量会集中在Radon域的曲率非零位置，从而可以有效地被去除.

本文首先简要介绍数据自相关偏移的原理，然后阐述在单程波多次波偏移中移除假象的主要流程，接下来在数值算例中以数据自相关偏移为例验证这种方法的有效性，最后给出相应结论.

数据自相关偏移的成像条件可以表示为

(1)

这里x, z代表笛卡尔坐标，ω代表频率，x_s代表震源位置，Image (x, z)代表偏移结果，下标down表示下行波场，up表示上行波场.P(x, z, ω, x_s)和M(x, z, ω, x_s)分别表示一次波场和地表多次波场，上标*表示共轭运算.地震数据D(x, z, ω, x_s)在深度域上下行延拓，同样的的地震数据D(x, z, ω, x_s)在深度域上上行延拓.

在波场延拓中利用的是傅里叶有限差分法(Fourier finite difference，FFD).这个方法首先是由Ristow和Rühl (1994)提出的，主要思路是把波场延拓分为背景速度和速度校正.背景速度v_ref取为当前深度上的最小速度，用背景速度先进行波场延拓后，再加上波场校正项，以适应真实速度.深度z+Δz处的波场函数p(k_x, ω, z+Δz)由深度z处的波场沿深度方向延拓获得

(2)

在具体求解上，第一步需要在频率波数域，利用v_ref完成简单高效的相移，第二步是在频率空间域用薄透镜项(Claerbout，1985)补偿垂向误差，第三步也是在频率空间域进行，将非垂向的传播效应考虑进来，用有限差分方程来进行求解.

地表多次波能够被分解为不同的阶数成分：

(3)

这里M(x, z, ω, x_s)代表地表多次波，不同的上标值表明了不同阶数的多次波.将公式(3)代入公式(1)，展开后可以得到

(4)

公式(4)中的第一部分是正确的成像结果，其中有两种成像类型：(1) 下行延拓的一次波场和上行延拓的一阶地表多次波场互相关；(2) 下行延拓的n阶地表多次波场和上行延拓的(n+1)阶多次波场互相关.第二部分和第三部分会产生假象，包含了三种成像误差：(1) 下行延拓的一次波场和上行延拓的m阶多次波场互相关产生的误差(m≥2); (2) 下行延拓的n阶和上行延拓的m阶多次波场互相关产生的误差(m=n+2)；(3) 下行延拓一次波场和n阶多次波场与上行延拓的m(m≤n)阶多次波场互相关产生的以及下行延拓的n阶和上行延拓的m阶多次波场互相关产生的误差(m>n+2).

其中第二部分假象能量在ADCIGs中会表现为抛物线型，这也是我们试图使用高精度抛物线型Radon变换去除的假象能量.使用FFD方法实现波场延拓和互相关成像条件之后，会得到含有假象的成像结果.为了进一步提升成像精度，需要考虑采用合适的方法压制假象.

在偏移过程中可以用多种方法来提取ADCIGs, 例如可以在在频率波数域求取传播角和方位角信息，在给定的成像条件下输出ADCIGs (Xu et al., 2011)；或者采用矢量波动方程进行波场延拓，并用能流密度矢量计算反射角，输出ADCIGs (Yoon and Marfurt，2004；王保利等，2013).Biondi和Symes (2004)提出可以先提取水平偏移距域共成像点道集(HOCIGs)，然后再通过转换公式将HOCIGs由偏移距域变换到角度域，提取ADCIGs.由于数据自相关偏移是多震源波场，为了提高计算效率，可以采用与Biondi和Symes (2004)类似的方法提取ADCIGs，在单程波波场延拓中，分别计算下行波场和上行波场，应用互相关成像条件，得到水平偏移距域共成像点道集(HOCIGs)，然后根据Sava和Fomel (2003)的转换公式将HOCIGs由偏移距域变换到角度域，提取ADCIGs.

Rickett和Sava (2002)提出在互相关成像条件中，将震源(下行)和接收(上行)波场分别添加一个对应的水平偏移距，就可以得到水平偏移距域共成像点道集(HOCIGs)，也就是

(5)

这里S(x, z, ω, x_s)代表下行延拓的震源波场，R(x, z, ω, x_s)是上行延拓的检波器波场，x_h代表的是地下水平半偏移距.由于FFD算子是在深度方向上延拓，根据公式(5)，在计算中可以快速高效地计算出HOCIGs.

Sava和Fomel (2003)提出角度域内的反射角γ和偏移距域内的参数x_h存在着如下的对应关系

(6)

利用这个公式，就可以将HOCIGs转化为ADCIGs，这个变换实际上就是一种倾斜叠加操作，可以表示为

(7)

在角度域内，部分多次波偏移假象与正确的偏移结果的曲率差别较大，可以采用抛物线型Radon变换将噪音分离出来.抛物线型Radon变换和反变换的公式如下(Sava and Guitton, 2005; Alvarez et al., 2007)

(8)

(9)

这里z是角度域内的深度，z₀是抛物线变换后Radon域内的深度, γ是反射角，q是变换曲率. g(γ)是由Biondi和Symes (2004)给出的对原来的Radon变换的优化，即

(10)

为了提高Radon变换的精度，进一步引入了高分辨率的抛物线型Radon变换.频率域内的Radon变换可以通过Fourier变换得到

(11)

其中j是虚数单位，x代表空间域内的坐标位置，q代表Radon域内的变换曲率，ω代表频率，对应的逆变换为

(12)

将公式(12)离散化后得到

(13)

其中l=0, 1, …, L; k=0, 1, …, K，u(x_l, ω)表示每道数据，v(q_k, ω)是Radon变换结果，q_k是曲率参数，一般是均匀离散化的，写成矩阵形式有

(14)

其中R是(L+1)(K+1)的矩阵，元素R_lk=e^-jωq_kx_l2，标准的最小二乘法求解抛物线Radon变换，解的形式为V=(R^HR)^-1R^HU.但是这种求解是不稳定的，Sacchi和Ulrych (1995)为了压制噪声，采用最大熵分布作为约束条件，构造了如下的求解方程

(15)

其中矩阵W为模型空间的加权矩阵，是一个对角阵.其中的元素由上次迭代的V获得，数据分辨矩阵R^HR+W^TW变为Hermit矩阵，可以通过Cholesky分解方法，将对称正定矩阵表示为一个下三角矩阵及其转置的乘积，进而求逆，获得方程组的解.高分辨率Radon变换利用稀疏约束反演，对q_k选择不同的加权系数，来提高变换的分辨率.如果偏移速度是准确的，经过高分辨率Radon变换处理后，在Radon域内，ADCIGs中的真实成像的能量会集中于曲率为零区域，而公式(4)中第二部分假象的能量集中在非零区域内.这样就可以很容易地切除掉假象，再做Radon逆变换就得到假象被压制的ADCIGs.本文在波场延拓中采用的是FFD方法，沿着深度域方向进行，相比于Wang等(2014b)采用的逆时偏移方法(RTM，reverse time migration)方法，可以快速地提取HOCIGs，计算速度可以大幅度地提升，FFD方法不足之处是对陡倾角成像不好，但是陡倾角是很难利用多次波，经过数据自相关成像的.对于长时间记录的多次波数据，使用FFD进行数据自相关偏移处理更为合适.

数值算例采用的是一个添加了水层的Marmousi模型，如图 1所示.这个速度模型横向有737个网格点，网格间距为12.5 m，纵向有950个网格点，网格间距为4 m.炮点间距为37.5 m，每炮有737道.检波器间距为12.5 m，检波点的间距为12.5 m.记录数据的时间长度为6 s，采样间隔为1 ms.图 2是正演获得的含多次波炮集记录，红色箭头标注的是部分多次波，由于Marmousi海底层速度较低，多次波能量不强.

图 1

Fig. 1

图 1 添加了水层的Marmousi速度模型 Fig. 1 The Marmousi velocity model with water layer

图 2

Fig. 2

图 2 含多次波的炮集记录, 红色箭头所示的是部分多次波 Fig. 2 The typical shotgathers generated by forward modeling, part of the multiples are pointed with red arrows

只利用一次波数据的常规FFD偏移中提取的的ADCIGs如图 3所示.这里的正传的是一个Ricker子波，而反传的是记录的一次波数据.对于一个准确的速度模型，RTM与FFD方法提取的角道集基本是一致的.但是RTM的成像精度比FFD更高，对于陡倾角的成像效果更好，表现在角道集上就是道集深度分辨率更高，含有更多的大角度成分.但是RTM提取ADCIGs需要的计算量远大于FFD方法.对于长时间记录的多次波实际数据，在速度模型精度不高的条件下，为了高效率地提取ADCIGs，FFD方法更为实用.

图 3

Fig. 3

图 3 一次波FFD偏移提取的角道集 Fig. 3 ADCIGs of the conventional FFD migration using primaries only

为了压制偏移中的低频噪声，在角度域内切除掉大角度成分，然后将其叠加，得到的偏移结果如图 4所示.可以看到图 4的成像结果中低频噪音被很好地压制，值得注意的是图 4中是没有假象的.

图 4

Fig. 4

图 4 一次波FFD偏移结果 Fig. 4 The migration result of the conventional FFD migration using primaries only

数据自相关偏移的ADCIGs如图 5所示，包含有多次波成分的地震数据正传作为下行波场，而同样的数据反传作为上行波场.采用和图 4一致的窗函数滤除大角度成分，叠加结果如图 6所示，可以看到成像结果中是包含假象的，这就是公式(3)中误差部分互相关所带来的.由于多次波能量不强，假象能量较弱.

图 5

Fig. 5

图 5 数据自相关FFD偏移角道集 Fig. 5 ADCIGs of data to data migration

图 6

Fig. 6

图 6 数据自相关偏移结果 Fig. 6 The preliminary migration result of data to data migration

图 7是将图 5中获得的ADCIGs经过高分辨率抛物线型Radon变换后的结果，正确的偏移图像其能量在经过抛物线型Radon变换后应当聚焦到零值附近，也就是q=0这条直线附近，而假象能量会聚焦到远离零值的区域.加以一个简单的窗函数加以去除，然后应用Radon逆变换回角度域，得到的ADCIGs如图 8所示.叠加各个角度结果，最后得到的自相关偏移结果如图 9所示.显然，相比于图 6，图 9偏移结果的信噪比大幅度的提升了.图 10所示的是图 6和图 9的残差剖面.本文使用的高分辨率Radon变换并不完全保幅，图 6与图 9量级并不一致，图 10中不仅仅只有假象能量，还有其他噪声和振幅不一致带来的成像能量损失.寻求一种高分辨率的保幅抛物线型Radon变换是下一步的工作目标.

图 7

Fig. 7

图 7 数据自相关偏移角道集变化到Radon域内道集 Fig. 7 The ADCIGs in Radon domain after the parabolic Radon transform

图 8

Fig. 8

图 8 Radon域内去噪后变换回角度域内的道集 Fig. 8 The ADCIGs after muting procedure and the adjoint Radon transform

图 9

Fig. 9

图 9 角度域内去噪后数据自相关偏移结果 Fig. 9 The final image after eliminating the migration artifacts in ADCIGs

图 10

Fig. 10

图 10 图 6与图 9的残差剖面 Fig. 10 The residual image of Fig. 6 and Fig. 9

多次波偏移通过角度域内压制假象的办法可以得到与常规FFD方法相近的偏移结果，并且不需要预测子波.对于实际采集数据，这有着重要意义.因为实际数据中震源函数的估计是相当难以处理的.FFD算子沿着深度方向进行延拓，因而在计算中它可以快速高效的计算出HOCIGs.在提取长时间记录的多次波数据的ADCIGs时，内存消耗少，计算速度快的FFD更为适用.

而数据自相关偏移结果中的假象能够被近似视为两部分的结合，第一部分也就是公式(4)中的第二项产生了深部的假象，第二部分也就是公式(4)中的第三项产生了浅层的假象.第一部分的假象，比如由一次波和高阶的地表多次波产生的假象，经过本文提出的的假象压制策略处理之后，可以得到很好的压制. Marmousi模型的数值算例验证了这一点.然而如何压制第二部分互相关产生的假象还需要进一步的探讨.本文中由于这一部分的假象是产生在设置的水层中，所以可以很容易地将其去除.但是对于水层较浅的实际数据这种处理方式就不再适用了.可能需要采用宽方位角采集的方式来部分减少假象.另一个需要考虑的是尽管FFD方法可以快速高效地计算HOCIGs，但是如果需要提取垂直偏移距域共成像点道集(VOCIGs)，与HOCIGs结合生成更稳定的ADCIGs时，FFD方法可能就不再适用，需要做进一步的修正处理.

感谢审稿专家和编辑部的鼎力支持和帮助.