地球物理学报  2015, Vol. 58 Issue (3): 993-1001   PDF    
引用本文
郑忆康, 王一博, 徐嘉亮, 常旭, 姚振兴. 2015. 数据自相关多次波偏移成像. 地球物理学报, 58(3): 993-1001, doi: 10.6038/cjg20150324   
ZHENG Yi-Kang, WANG Yi-Bo, XU Jia-Liang, CHANG Xu, YAO Zhen-Xing. 2015. Imaging multiples by data to data migration. Chinese J. Geophys. (in Chinese), 58(3): 993-1001, doi: 10.6038/cjg20150324   
数据自相关多次波偏移成像
郑忆康1,2, 王一博1, 徐嘉亮1,2, 常旭1, 姚振兴1    
1. 中国科学院地质与地球物理研究所, 北京 100029;
2. 中国科学院大学, 北京 100049
投稿时间:2014-06-09    最后修改时间: 2015-01-30 .
基金项目：国家科技重大专项子课题(2011ZX05030-004-001)资助.
作者简介：郑忆康,男,博士研究生,主要从事地震波成像与速度建模方法研究.E-mail:zhengyk@mail.iggcas.ac.cn
通讯作者：王一博,研究员.E-mail:wangyibo@mail.iggcas.ac.cn
摘要：在常规偏移方法中一般都需要压制地震数据中的多次波,仅利用一次波信息成像,把自由表面反射的多次波视为噪声,但是在多次波中也包含着地下结构信息,应该将其充分利用到成像中来.事实上,已经有不少成像方法试图利用多次波信息,但是大部分方法都需要对多次波进行预测.本文提出了基于傅里叶有限差分偏移算子的数据自相关偏移方法.在这种偏移方法中,对含有一次波和多次波的地震数据,分别进行下行和上行延拓,然后直接利用常规的互相关成像条件成像.由于波场延拓采用了傅里叶有限差分算子,其计算效率高,能够很好地对复杂介质中的地震数据进行延拓.在数值试验中,使用了一个含散射点的三层模型和Marmousi模型.合成数据测试结果表明,这种方法可以对更大范围的地下构造成像,比常规的只利用一次波的傅里叶有限差分法照明度更好,并且在浅层可以提供更高的分辨率.我们提出的数据自相关策略易于实现且避免了繁杂的多次波预测,这对于复杂地下构造成像可能有着重大意义.
关键词： 多次波成像     数据自相关     傅里叶有限差分    
Imaging multiples by data to data migration
ZHENG Yi-Kang1,2, WANG Yi-Bo1, XU Jia-Liang1,2, CHANG Xu1, YAO Zhen-Xing1    
1. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: In conventional migration methods, free surface related multiples are regarded as noise and only the primary wave is utilized in imaging. However, the multiples also contain information of subsurface structures and can offer extra illumination in migration procedure. There are already numerous methods to take advantage of multiples, but most of them need multiples prediction, which is time consuming. To avoid multiples prediction and wavelet estimation, this paper presents a new multiples migration method called data to data migration.
Based on the migration method that uses primaries and multiples simultaneously, this paper presents one-way wave equation migration of data to data migration. The FFD (Fourier finite difference) method is employed in wave field extrapolation. For the data containing primaries and free-surface related multiples, the data to data migration replaces the source wavelet function and the primaries in conventional FFD migration with the recorded data. The same crosscorrelation imaging condition is used. In the FFD migration method, the field data is transformed into the frequency domain. Compared to RTM(reverse time migration), it has higher computation efficiency and less low frequency noise. The FFD is more suitable for the implementation of data to data migration.
Compared to conventional FFD migration that uses primaries only, the images generated by data to data migration based on FFD operator have better illumination and higher resolution in the shallow zone. The images of data to data migration contain artifacts produced by the corsscorrelation of different seismic events. In the deep zone, the migration results are inferior to conventional migration. One reason is the influence of artifacts; another is that multiples need longer receiving time. In the dataset, higher order multiples may not yet be received completely. The results of the conventional migration and data to data migration have different polarities. Assumed the wavelet and the primaries have positive polarity, and then the conventional migration result has a positive polarity. The first order multiple is reflected by the free surface once, so it has negative polarity. It is crosscorrelated with the primary that has positive polarity and generates the result with negative polarity. The second order multiple that has positive polarity is crosscorrelated with the first order multiple that has negative polarity and generates the result with negative polarity. And the m+1 order multiple is crosscorrelated with the m order multiple which have different polarity and generates the result with negative polarity. Therefore the final result of data to data migration has negative polarity and it has different polarity with conventional migration. The computational speed of FFD is much faster than RTM which makes FFD appropriate for the multiples data with long record length.
Recorded data without direct wave are utilized in data to data migration procedure. The numerical examples verify its effectiveness. The FFD operator improves the computational efficiency greatly. The proposed approach has three key advantages: (1) Compared to conventional migration method, data to data migration has wider illumination area and better resolution of scattering points. (2) It can generate subsurface image without multiples prediction, which is time consuming. (3) It can offer fairly good image of shallow reflectors. Another point is wavelet estimation. For the field data, how to extract wavelet function is intricate and prone to error. Data to data migration avoids this complicated process,which may be significant for imaging subsurface structures in the processing of field data. We should also notice that although most energy are imaging correctly, there are still artifacts in the final result.Least squares migration and wide azimuth acquisition technology will partly decrease the artifacts.
Key words: Migration of multiples     Data to data migration     FFD    

1 引言

在常规的地震偏移中，多次波往往被视为噪音，并由此发展出了很多压制多次波的方法.例如f-k域的波场外推法，建立多次波模型再去除多次波(Wiggins，1988)；f-x域的多次迭代预测法，即 SRME(Adaptive surface-related multiple elimination)(Verschuur et al., 1992)； 逆散射级数法压制多次波(Weglein et al., 1997)等.但是在多次波中也包含着地下结构信息，应该将其充分利用到成像中来.实际上，如何在在偏移中利用地表多次波已经被许多研究者探讨过.Verschuur和Berkhout(2005)提出将地表多次波转换为一次波再将其利用到常规偏移方法中.他们的方法的核心问题就是如何从记录的数据中恢复一次波成分.地震干涉法(Schuster，2009; Wapenaar，2006b; Curry and Shan， 2006; Wang et al., 2009)理论上保证了地表多次波能够通过互相关和求和被转换为一次波，其恢复的数据能够在常规地震偏移中加以利用(Yu and Schuster 2002; Schuster et al., 2004; Jiang et al., 2007; He et al., 2007).也有其他的方法，例如在偏移过程中将自由表面多次波直接视为面积震源项(Guitton，2002; Shan，2003; Muijs et al., 2007)，相似的方法有Artman(2006)提出的将记录中的无源成分直接用于地下成像.Liu等(2011)用地震记录和预测出的多次波实现多次波成像.他们的方法体现了在偏移中利用多次波的有效性，但是大部分方法仍然需要做多次波预测，而这是极为耗时而且容易产生误差的.

与传统的Kirchhoff积分法相比，傅里叶有限差分法(Fourier finite difference，FFD)没有采用高频近似，对复杂构造的处理能力更强.与RTM(reverse time migration)相比，在非高陡倾角情形下二者成像效果相当，但是FFD计算效率更高，偏移剖面噪声更少，在实际地震资料处理中有广泛的应用前景.

本文在Wang等(2014a)同时利用一次波和多次波成像方法的基础上，根据FFD成像方法的原理，实现了单程波的数据自相关偏移成像.假设已经有包含了一次波和地表多次波地震记录数据(直达波已经被切除)，数据自相关偏移将常规FFD中的震源项和一次波都替换为完整的记录地震数据，并采用与常规FFD方法中一致的互相关成像条件.数值实验表明提出的方法能够有效地对散射点和浅层反射点成像.这种自相关偏移方法并不需要做多次波预测，这有可能使其在一般的复杂地下构造成像中得到广泛利用. 2 基于FFD的数据自相关偏移方法的基本原理

假设一个震源在地表被激发，记录的地震数据能够在频率域(Wang et al., 2014b)被近似地表达为

这里z0是震源和检波器的深度，R表示真实的地下响应，S(z0)表示真实的震源项，D(z0)表示记录的包含有一次波和地表多次波的地震数据，它能够被表示为：

这里P(z0)表示一次波，M(z0)表示地表多次波.用W(z0)表示我们模拟的震源子波，如果S(z0)能够被W(z0)很好地近似，同时假设我们已经获得了一个偏移速度模型X，那么就能利用傅里叶有限差分或者其他建模方法获得正向传播(下行波场)的合成数据PF(z0)，如下式所示：

这个合成数据PF(z0)能够被用来和记录中一次波P(z0)一起成像，传统的FFD互相关成像条件，即计算上行波场(检波器波场)和下行波场(震源波场)的互相关值作为反射系数.对于震源的估计和地表多次波的去除是传统FFD方法中的关键问题.

记录的地震数据D(z0)能够被用来当作重建的面积震源，如果我们有偏移的速度模型X，那么得到的地表相关多次波数据MF(z0)能够被表示为

一次波数据经历地表反射，重新经过地下响应被接收，数据中包含了地下介质信息.这里的多次波数据MF(z0)能够和记录数据中的M(z0)一起成像，这就是利用多次波偏移的基本思路.但是这种方法需要提前做地表多次波预测，耗时量大且难度高.

事实上，公式(3)中的震源项W(z0)和(4)中的面积震源D(z0)能够被结合起来组合成一个新的震源函数，这样得到的合成正向传播数据能够被表示为

这就是Wang等(2014a)同时利用多次波和一次波成像的公式，即正传一个震源子波和含有多次波的地震数据组合成的混合震源，反传地震数据，然后应用成像条件.

公式(5)可以进一步拆分为两部分.第一部分如下：

这一部分包含了常规FFD方法中一次波和震源子波的准确成像结果以及由于多次波和子波成像带来的噪音.第二部分可以表示为

这就是我们所称的数据自相关偏移，正传与反传的都是含有多次波的地震数据.实际上第一部分的有效成分就是传统的FFD偏移方法，多次波贡献体现在第二部分，也就是数据自相关偏移中.

这种方法的基本思路是下行的一次波场正向外推和上行的一阶地表多次波反向外推场互相关，下行的一阶地表多次波正向外推和上行的二阶地表多次波互相关，下行的二阶地表多次波正向外推和上行的三阶地表多次波互相关…，叠加之后得到反射剖面.利用这些规则，一次波和多次波能够被同时用来偏移而不需要多次波预测.我们的方法中有三个关键点：(1)记录的地震数据(包含着一次波和地表多次波，不包含直达波)被视为一个面积震源用来下行延拓；(2)同样的的地震数据用来上行延拓得到检 波器波场；(3)在选定的频率和波场位置上对震源波场和检波器波场应用成像条件.

Wang等(2014a，b)文中采用的是RTM方法，这个方法的优点是成像精度高，但是对内存和计算资源要求高.多次波数据往往是长时间记录，使用RTM方法需要大量的计算时间.FFD将地震数据转换到频率域内，计算效率更高，低频噪声相比于RTM更少，更加适用于我们的数据自相关偏移方法.

常规FFD中广泛使用的互相关成像条件(Claerbout，1971)也能在我们的方法中加以利用.这种成像条件可以表示为

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

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

这里M(x，z，ω，xs)代表地表多次波，上标表明了不同阶的多次波.利用公式(9)，公式(8)可以表示为

进一步展开得到

为了清楚地表明公式(9)右侧表达式的不同成分，我们将其写为

这里，第一部分是准确的成像结果，第二和第三部分产生了人为偏移误差.第一部分中有两种成像类型：(1)上行延拓的一次波场和下行延拓的一阶地表多次波场成像；(2)上行延拓的n阶地表多次波场和下行延拓的(n+1)阶多次波场成像.第二部分包含了两种成像误差：(1)上行延拓的一次波场和下行延拓的m阶多次波场互相关产生的误差(m≥2);(2)上行延拓的n阶和下行延拓的m阶多次波场互相关产生的误差(m≥n+2).第三部分的成像误差是下行延拓一次波场和n阶多次波场与上行延拓的m(m≤n)阶多次波场互相关产生的.由于高阶的多次波经历多次反射，其能量低于低阶多次波能量，在误差中低阶多次波贡献也更大.

本文中利用的是FFD来实现波场延拓(Ristow and Rühl，1994)，即将延拓分成背景常速度vref加上速度校正来完成，首先在[ω，kx]域利用参考速度vref高效地完成简单的相移，然后在[ω，x]域采用合适的校正来补偿垂向传播的误差，最后在[ω，x]域内考虑横向速度变化，进行类似于剩余波场传播的校正. 3 数值算例

我们首先采用图 1所示的含有多个散射点的三层速度模型进行测试.此三层模型在水平方向上有900个网格点，在垂直方向上有401个网格点，网格间距都是10 m.炮点间距为40 m，每一炮有900道记录，地表全覆盖接收来进行数值模拟.检波器间距为10 m，同时震源和检波器深度都是10 m.时间记录的长度是8 s，采样间距为1 ms.第一个和最后一个震源位置分别为2010 m和7010 m.

图 1 含散射点的三层速度模型 Fig.1 The three layer velocity model with scattering points

图 2是常规的只利用了一次波的FFD结果.此算例采用合成的Ricker子波作为震源，记录的一次波作为检波器数据反向传播.图 3是数据自相关偏移得到的结果.这个成像结果是把含有多次波的记录数据正向延拓作为震源波场，同样的记录数据反传作为检波器波场.比较图 2和图 3，图 2是不含假象的，图 3在浅层有更好的成像结果和覆盖范围，但是其中包含有一定的假象.如图中方框1和3所示，图 3相比于图 2，照明范围更大，能提供更大区域内的地下信息；如方框2所示，对于散射点，图 3的分辨率更高，数据自相关偏移在对于多次波充分发育的数据，成像精度是更高的.在深部区域，图 3成像结果不及图 2的，一个原因是假象的干扰，其次是多次波需要更长的接收时间，这里我们采用的数据记录时间与常规一次波相同，更高阶的多次波尚未完全被接收，影响了成像结果.为了消除成像误差，有三种进一步的改进策略：一是Wang等(2014b)文中提出的在角度域内的共成像点道集上滤波，去除假象；二是采用最小二乘偏移方法，假象产生的反射数据在迭代中消去；三是宽方位角接收，不同方位角的叠加偏移剖面可以压制假象.

图 2 常规一次波FFD偏移结果 Fig. 2 The migration result of conventional FFD

图 3 数据自相关多次波FFD偏移结果 Fig. 3 The result of data to data Fourier finite difference migration using multiples

图 4是采用RTM方法得到的结果，与图 3相比较，其成像分辨率更高，存在着略为明显的低频噪声，需要进一步的去噪处理.这与Guan等(2009)、Higginbotham等(2010)的结论一致.对于这个模型，同样炮数使用同样核数的CPU进行计算，Wang等(2014a，b)使用的RTM方法需要40分27秒，而FFD方法只需要5分47秒，显然FFD的计算效率 是远远优于RTM的，更适合于含长时间记录的多次波数据.

图 4 数据自相关多次波RTM偏移结果(Laplace滤波后) Fig. 4 The result of data to data reverse time migration using multiples(filtered by Laplacian)

值得注意的是，数据自相关偏移与常规偏移所得到的偏移结果相位是相反的(Wang et al., 2014b).考虑子波与一次波相位为正，二者互相关成像结果为正；经地表反射的一阶多次波相位为负，其与相位为正的一次波成像，所得成像结果相位为负; 二阶地表多次波相位为正，其与相位为负的一阶多次波成像，所得成像结果相位为负……叠加之后所得数据互相关偏移结果相位自然为负，与常规偏移结果相位相反.

图 5是一个添加了水层的Marmousi速度模型.这个模型有水平方向的737个网格点和垂直方向的950个网格点，水平向与深度向网格间距分别是12.5 m、4 m.炮点间距是37.5 m，每一炮对应737道记录，检波器间距为12.5 m，炮点和检波点深度均为4 m.记录的时间长度为6s，采样间隔为0.4 ms.

常规的利用一次波的FFD结果如图 6所示，图 7是数据自相关偏移成像结果.相比于图 6，图 7中的浅层反射层由于更高的覆盖次数和更好的照明度成像结果更好.但是也可以看到在偏移结果中存在着由成像误差带来的假象.同样地，图 6与图 7偏移结果相位是相反的.这也是在实际数据处理过程中，使用数据自相关偏移结果和常规偏移结果时所需要注意的一点.

图 5 Mamousi速度模型 Fig. 5 The Marmousi velocity model

图 6 常规FFD偏移结果 Fig. 6 The migration result of conventional FFD

图 7 数据自相关多次波偏移结果 Fig. 7 The result of data to data migration using multiples
4 结论

数据自相关多次波偏移不需要多次波预测而且全部的地震记录(包含一次波和地表多次波并且直达波被切除)都用于偏移过程.数值试验证实了这种数据自相关偏移的可靠性.我们提出的方法有3个主要优点：(1)相比于常规偏移方法，它对于地下有更大的成像范围且能提高散射点的照明度；(2)它能得到相当好的散射点成像结果并且不需要多次波预测，而这种预测是耗时量大且容易产生误差的；(3)浅层反射点的成像效果与常规偏移方法相当.另一点需要指出的是震源子波的估计.对于实际数据，震源项的估计是一个难题，而这种数据自相关偏移方法则避免了繁琐的震源子波估计.这种方法对于实际地震资料处理中的地下复杂构造成像可能有着重要意义.

但值得注意的是，尽管大部分自由表面多次波都正确成像了，结果中依然有不同的地震事件互相关产生的偏移假象.进一步的研究工作是如何压制假象，提高成像精度.一种策略是采用宽方位角接收，叠加之后，不同方位角中真实反射点的成像被加强，另一种策略是最小二乘偏移，通过正演数据和实际地震数据的比较，在逐次迭代中压制假象.

参考文献
[1] Artman B. 2006. Imaging passive seismic data. Geophysics, 71(4): SI177-SI187.
[2] Claerbout J F. 1971. Toward a unified theory of reflector mapping. Geophyscis, 36(3): 467-481.
[3] Curry W J, Shan G. 2006. Interpolation of near offsets with multiples and prediction-error filters. 68th Ann. Internat. Mtg, Session: A041, Eur. Assn. Geosci. Eng.
[4] Guan H, Kim Y C, Ji J, et al. 2009. Multistep reverse time migration. The Leading Edge, 28(4): 442-447.
[5] Guitton A. 2002. Shot-profile migration of multiple reflections. 72nd Annual International Meeting, SEG, Expanded Abstracts, 1296-1299.
[6] He R Q, Hornby B, Schuster G. 2007. 3D wave-equation interferometric migration of VSP free-surface multiples. Geophysics, 72(5): 195-203.
[7] Higginbotham J, Brown M, Macesanu C, et al. 2010. Onshore wave-equation depth imaging and velocity model building. The Leading Edge, 29(11): 1386-1392.
[8] Jiang Z Y, Sheng J M, Yu J H, et al. 2007. Migration methods for imaging different-order multiples. Geophysical Prospecting, 55(1): 1-19.
[9] Liu Y K, Chang X, Jin D G, et al. 2011. Reverse time migration of multiples for subsalt imaging. Geophysics, 76(5): WB209-WB216.
[10] Muijs R, Robertsson A J O, Holliger K. 2007. Prestack depth migration of primary and surface-related multiple reflections: Part I—Imaging. Geophysics, 72(2): S59-S69.
[11] Ristow D, Rühl T. 1994. Fourier finite-difference migration. Geophysics, 59(12): 1882-1893.
[12] Shan G J. 2003. Source-receiver migration of multiple reflections. 73rd Annual International Meeting, SEG, Expanded Abstracts, 1008-1011.
[13] Schuster G T, Yu J, Sheng J, et al. 2004. Interferometric/daylight seismic imaging. Geophysical Journal International, 157(2): 838-852
[14] Schuster G T. 2009. Seismic Interferometry. Cambridge: Cambridge Press.
[15] Verschuur D J, Berkhout A J, Wapenaar C P A. 1992. Adaptive surface-related multiple elimination. Geophysics, 57(9): 1166-1177.
[16] Verschuur D J, Berkhout A J. 2005. Transforming multiples into primaries: experience with field data. 75th Annual International Meeting, SEG, Expanded Abstracts, 2103-2106.
[17] Wang Y B, Luo Y, Schuster G T. 2009. Interferometric interpolation of missing seismic data. Geophysics, 74(3): SI37-SI45.
[18] Wang Y B, Chang X, Hu H. 2014a. Simultaneous reverse time migration of primaries and free-surface related multiples without multiple prediction. Geophysics, 79(1): S1-S9.
[19] Wang Y B, Zheng Y K, Zhang L L, et al. 2014b. Reverse time migration of multiples: eliminating migration artifacts in angle domain common image gathers. Geophysics, 79(6): S263-S270.
[20] Wapenaar K. 2006a. Green's function retrieval by cross-correlation in case of one-sided illumination. Geophysical Research Letters, 33(19): L19304, doi: 10.1029/2006GL027747.
[21] Wapenaar K. 2006b. Green's function representations for seismic interferometry. Geophysics, 71(4): SI33-SI46, doi: 10.1190/1.2213955.
[22] Weglein A B, Gasparotto F A, Carvalho P M, et al. 1997. An inverse-scattering series method for attenuating multiples in seismic reflection data. Geophysics, 62(6): 1975-1989.
[23] Wiggins J W. 1988. Attenuation of complex water-bottom multiples by wave-equation-based prediction and subtraction. Geophysics, 53(12): 1527-1539.
[24] Yu J H, Schuster G T. 2002. Joint migration of primary and multiple reflections in RVSP data. 72nd SEG, Expanded Abstracts, 2373-2376.