2013, Vol. 56
2. 中国石油大学(北京) CNPC物探重点实验室, 北京 102249
2. CNPC Key Lab of Geophysical Exploration, China University of Petroleum, Beijing 102249, China
海底地震仪OBS观测是一种将检波器直接放置于海底的地震观测系统[1].近年研究表明，在海洋地球物理调查中，可以利用OBS探测海底地质构造、天然气水合物和游离气、海底油气资源，调查地下岩层岩性，确定海底的弹性和各向异性等物性参数[2].在油气勘探中，由于OBS造价昂贵，成本较高，OBS站点一般分隔较远，不同于海底电缆地震勘探观测系统[3].在海洋深水区域，海底电缆的应用受到限制，因此，OBS技术在海洋深水勘探中具有重要的作用.OBS数据和其他海洋勘探数据（拖缆数据、海底电缆数据等）一样面临着多次波干扰问题，由于海水表面反射系数大，与表层相关的多次波压制就成为OBS数据处理中最需解决的问题之一[4-5].
多次波的压制方法主要分为两类[6-7]：一类是利用一次波和多次波的特征差异来压制多次波，统称为滤波法，如预测反褶积、抛物Radon变换等；另一类是从地震数据中预测多次波，然后将其从地震数据中减去，该类方法基于波动方程理论，也称为波动方程预测相减法，主要有：波场外推法[8]、反馈迭代法[9]、逆散射级数法[10]等.本文研究涉及的方法属于反馈迭代法的范畴.在反馈迭代方法中，SRME是经典的自由表面多次波压制方法[9].SRME一般分为两步：多次波预测和自适应相减，在多次波与一次波同相轴相交的情况下，基于一次波能量最小假设的自适应相减方法会损伤一次波，黄新武等[7]、薛亚茹等[11]、李鹏等[12]都对该方法做过深入的研究. Berkhout（2006）将SRME理论推广到反数据域压制多次波，该方法与正数据域的SRME相比，理论表达简单，避免了自适应相减过程，因此地震记录可以完全保留有效波信息[13].马继涛等在反数据域压制多次波理论基础上，提出了平面波域反数据压制多次波的方法，提高了矩阵求逆的稳定性和计算有效性[14-15].Van Borselen等[16]，Biersteker[17]在反演理论基础上，通过最小能量约束研究了一次波的反演估计问题，从反演的角度避免了自适应相减的过程.van Groenestijn和Verschuur（2009）[18]提出了EPSI方法，在稀疏约束的基础上，通过多维反演迭代求解估计一次波，该方法可以精确地估计一次波，不需要自适应相减，且计算稳定.
对于不同的地震观测方式，采用SRME的方法略有不同.马继涛等[19]在SRME理论基础上，研究了海底电缆地震数据中自由表面多次波的压制问题，得到了较好的效果.反馈迭代法需要全波场的信息，如果地震数据有缺失，需要对地震数据做插值重建处理.然而对于OBS数据，检波点非常稀疏，插值结果难以保真，因此仅仅采用OBS数据利用反馈迭代法预测多次波具有很大的局限性[3, 5].Verschuur和Neumann（1999）[4]在SRME基础上，利用拖缆数据和OBS数据联合研究了OBS多次波压制问题.本文将SRME方法拓展到EPSI方法，联合海洋常规拖缆数据和OBS数据，研究OBS数据自由表面多次波压制方法.首先，详细介绍了反馈迭代自由表面多次波压制理论；然后在此基础上推导了联合拖缆数据压制OBS多次波的EPSI原理，给出了OBS多次波压制的具体步骤；最后，利用模拟数据验证了本文方法的有效性.
2 反馈迭代模型与多次波预测理论地震数据是炮点坐标xs=[sx，sy，sz]、检波点坐标xr=[rx，ry，rz]和时间t的函数，将地震数据表示为d（xs，xr，t），在频率域地震数据变为d（xs，xr，ω）.省略炮点和检波点坐标，某一频率切片的地震数据可以用矩阵D表示.地震数据可以表示为一次波（指除了自由表面多次波之外的波场信息，即包含层间多次波）和自由表面多次波的和
|
(1) |
其中，P表示一次波，M表示多次波，X为一次波脉冲响应，S为震源矩阵，R为表面反射系数矩阵.通过公式（1），可以推导出压制多次波的三种反馈迭代方法.
（1）SRME方法.由公式M=XRD，可以看出，表面相关多次波可以通过地震数据D作用算子XR预测得到.因此有M=PS-1RD，S-1R是表层相关算子，通过迭代法就可以预测多次波模型，然后利用自适应相减方法消除掉S-1R的影响，这就是经典的SRME方法.
（2）反数据域方法.公式（1）可以写成（I-XR）D=XS，那么
|
(2) |
由于表层相关算子RS-1不含有任何与时间相关的信息，在反数据域是一个位于零时间附近聚焦点，因此在反数据域地震数据变成一次波加上表面相关算子，在时间域可以通过切除零时间附近的表面相关算子，从而得到一次波的估计.该方法需要做矩阵的逆运算，因此需要采用稳定的矩阵求逆方法.
（3）EPSI方法.SRME方法需要自适应相减，反数据域方法在频率域矩阵求逆存在不稳定性问题，而EPSI方法很好地解决了这些问题，该方法直接从方程（1）入手，对时间域一次波脉冲响应施加稀疏约束，采用迭代反演的方法估计一次波与多次波.假设自由表面反射系数为-1，那么方程（1）中的R=-I.采用Claerbout对优化问题的描述方式[20]，EPSI方法可以描述成下面的优化问题：
|
(3a) |
|
(3b) |
|
(3c) |
其中‖A‖2
2表示对矩阵A所有元素求平方和，|A|0表示A的L-0范数.公式（3a）是由公式（1）得到的，采用L-2范数最小描述公式（1）两端的能量差最小；公式（3b）表示时间域一次波脉冲响应是稀疏的，采用零范数最小描述其稀疏性，
上述EPSI方法是针对常规观测系统地震数据的方法，要求检波点与炮点位于同一观测平面上，且只用一个数据体.而OBS数据观测方式不同于常规观测方式，检波点稀疏且放置于海底，炮点和检波点并不位于同一观测平面上，因此本文研究联合拖缆数据压制OBS多次波.图 1显示了OBS波场传播示意图，可以看出，OBS接收到的表面相关多次波（包括上行波场和下行波场）都可以看作是拖缆采集数据Dsur与OBS一次波脉冲响应Xobs的多维褶积，因此我们将公式（1）改写为
|
图 1 OBS波场传播示意图 （a）OBS下行波场中的多次波示意图；（b）OBS上行波场中的多次波示意图.左：OBS多次波；中：拖缆采集数据Dsur；右：OBS一次波（包括直达波）脉冲响应Xobs. Fig. 1 The sketch map of OBS wavefields (a) OBS downgoing wavefield; (b) OBS upgoing wavefield. Left:OBS multiples; middle:streamer data Dsur; right:the primary (including direct wave) impulse response of OBS data Xobs. |
|
(4) |
其中，Dsur是常规拖缆数据，Dobs是OBS数据，Xobs是OBS一次波脉冲响应，S是震源矩阵.在计算过程中，由于OBS数据检波点非常稀疏，我们仅仅需要在共检波点道集上压制多次波，对某一共检波点道集来说，Dobs和Xobs是一维向量，Dsur是拖缆数据的二维矩阵.在采集过程中，拖缆和OBS数据一般采用同样类型的震源，因此这里假设震源特性S是不变的.公式（4）与公式（1）的差别是公式（4）含有拖缆数据和OBS数据两个数据体，但求解过程与公式（1）近似.类似于公式（3），我们将公式（4）写成下面的优化问题
|
(5a) |
|
(5b) |
|
(5c) |
类似于van Groenestijn和Verschuur（2009）[18]给出的公式（3）的求解方法，这里我们给出公式（5）的详细求解方法.
（1）首先，确定Xobs及其梯度ΔXobs.将公式（5a）写成目标函数的形式
|
(6) |
对上述目标函数求导得到Xobs的梯度
|
(7) |
其中（•）H表示（•）的复共轭转置，Dobs -XobsS+ XobsDsur是残差项.（Dobs-XobsS+XobsDsur）SH的作用是将残差中的OBS一次波信息转化为OBS一次波脉冲响应，而（Dobs-XobsS+XobsDsur）Dsur H的作用是将残差中的OBS多次波信息转化为OBS一次波脉冲响应.如果迭代初始值设为：Xobs （0）=0，S（0）=0，那么经过一次迭代的ΔXobs结果是OBS数据与拖缆数据的多维相关，即
|
(8) |
（2）其次，满足时间域Xobs的稀疏性.为了采用稀疏约束，van Groenestijn和Verschuur（2008）[22]研究了直接对X加稀疏约束的方法，结果表明反演收敛速度非常慢.为了提高计算效率，van Groenestijn和Verschuur（2009）[18]将稀疏约束加到X的更新项ΔX上.类似地，本文将约束加到OBS一次波脉冲响应更新项ΔXobs上.具体做法是将ΔXobs变换到时间域，选择一个时间窗口，在时间窗口内保留一个或者几个时间域ΔXobs的最大值，其他设为零，该最大值代表了时间域的同相轴，这样保证了ΔXobs在时间域同相轴的个数比较少，从而满足了稀疏性.随着迭代次数的增加，选择的时间窗口依次增大以保证迭代能够收敛，然后将稀疏的ΔXobs加到一次波脉冲响应上，得到一次波脉冲响应下一次迭代的结果
|
(9) |
其中
|
(10) |
其中i，j是矩阵坐标，“• *”表示矩阵元素对应相乘.第一次迭代的时间窗口选取要适当，不能将多次波同相轴包括进来参与计算.随着迭代次数的增加，残差里面浅层一次波及其产生的多次波能量已经被减掉了，这有助于估计深层的弱一次反射.
（3）最后求解震源特性S.震源特性可以在时间域通过匹配Xobs （i）和Dobs+X obs（i） Dsur得到，也可以通过频率域的除法得到.
至此，已经得到迭代一次后的一次波脉冲响应Xobs和震源特征S，通过多次迭代即可得到最终结果.估计的一次波可以直接通过Xobs和S的乘积得到
|
(11) |
也可以采用保守的方法，即从原始数据中减去估计的多次波
|
(12) |
通过少量的迭代计算浅层对应的多次波就能够估计出来，而深层产生的多次波能量非常弱，因此在精度要求不高的情况下，可以用较少的迭代运算来估计多次波，然后采用保守方法计算一次波，从而提高计算效率.对于海底地震仪不同分量（压力和X，Y，Z分量）的数据，可以分别利用EPSI估计各分量数据中的一次波.本文算例以压力分量为例说明EPSI压制OBS多次波的方法.
4 模型算例用一个模拟的OBS数据压力分量和拖缆数据说明本文提出算法的有效性.图 2a是速度模型，图 2b是模拟的OBS共检波点道集，由于炮点和检波点不在同一平面上，因此直达波具有双曲线形态.为便于比较，这里也给出了吸收边界条件下模拟的不含表面相关多次波的记录（图 2c）.图 3是模拟的拖缆数据.比较图 2b和图 2c可以看出OBS数据多次波干扰非常严重，尤其在深层部分，一次波能量几乎被多次波淹没，例如箭头所指的第四个反射界面，一次波和多次波几乎完全重合，对于这样的同相轴SRME自适应相减会有很大的局限性.
|
图 2 （a）速度模型；（b）自由表面模拟的OBS数据；（c）吸收边界条件模拟的OBS数据 Fig. 2 (a) Velocity model; (b) Synthetic OBS data with free-surface boundary conditions; (c) Synthetic OBS data with absorbing boundary conditions |
|
图 3 模拟的拖缆数据（前5炮数据） Fig. 3 Synthetic streamer data (first five shot gathers) |
图 4-图 7显示了EPSI方法迭代的结果.随着迭代次数的增加，浅层的梯度逐渐衰减，说明在迭代过程中选择时间窗口应逐渐变大，这样才能把新的一次波能量包含进来（图 4）.从图 5可以看出，一次波脉冲响应随着迭代次数的增加，逐渐地被反演出来，当迭代10次时，浅层的一次波脉冲响应已经被反演出来，由于时间窗口还没有达到深层对应的时间，因此深层的一次波脉冲响应没有被反演出来.正是由于迭代过程是从浅层到深层逐步地估计一次波，所以保证了一次波和多次波能够很好的预测和分离.从EPSI估计得到的一次波数据（图 6c）可以看出EPSI得到的一次波非常清晰，很好地保持了有效信号.图 7给出了估计的多次波随着迭代次数的变化情况，随着迭代次数的增加，与深层有关的多次波逐渐被反演出来，由于深层有关的多次波能量非常弱，第一次迭代的结果（图 7a）就能够反映主要的多次波能量.注意箭头所指的深层有关的多次波能量，在第一次迭代结果中并没有被估计出来，随着迭代次数的增加它被很好地估计出来了.由于深层有关的多次反射波能量很弱，利用EPSI方法时，可以选用较少的迭代次数得到多次波，然后根据公式（12）得到一次波，这样可以提高计算效率，当然也会降低计算精度.最后比较EPSI、SRME自适应相减和吸收边界条件模拟得到的一次波记录（图 8）.这里SRME自适应相减采用的是非稳态自回归自适应相减方法[23]，该方法充分考虑了子波的时变性.从图 8可以看出SRME方法由于深层多次波和一次波相交比较严重，效果较差，第四个反射轴与多次波几乎完全重合，SRME自适应相减将其当作多次波能量直接减掉了，第五个反射轴处多次波同相轴非常复杂，SRME自适应相减结果也不理想.而EPSI方法可以很好地估计一次波，与吸收边界条件模拟得到的一次波记录非常相近，充分说明EPSI方法可以很好地压制OBS数据中的多次波.
|
图 4 不同迭代次数的梯度ΔXobs （a）迭代1次；（b）迭代10次；（c）迭代20次. Fig. 4 The gradient ΔXobs during different iterations (a) Iterations i=1;(b) Iterations i=10;(c) Iterations i=20. |
|
图 5 不同迭代次数的一次波脉冲响应Xobs （a）迭代1次；（b）迭代10次；（c）迭代20次. Fig. 5 The primary impulse response Xobs during different iterations (a) Iterations i=1;(b) Iterations i=10;(c) Iterations i=20. |
|
图 6
不同迭代次数估计的一次波数据 |
|
图 7 不同迭代次数估计的自由表面多次波数据XobsDsur （a）迭代1次；（b）迭代10次；（c）迭代20次. Fig. 7 The estimated multiples XobsDsur during different iterations (a) Iterations i=1;(b) Iterations i=10;(c) Iterations i=20. |
|
图 8 不同方法得到的一次波结果比较 （a）SRME自适应相减；（b）EPSI方法由公式（11）得到的一次波；（c）EPSI方法由公式（12）得到的一次波. Fig. 8 The comparison of obtained primaries between different methods (a) SRME; (b) EPSI by Eq.(11);(c) EPSI by Eq.(12). |
在OBS偏移成像过程，由于检波点鬼波数据有较大的照明范围，常常用来做镜像偏移[5]，这时需要用到检波点鬼波信息，在图 8得到的一次波记录中不含有检波点鬼波信息.为了得到含有检波点鬼波信息的一次波记录，首先需要将OBS数据中的直达波（图 9a中虚线以上的部分）切除作为EPSI的输入数据.图 9c显示的多次波是不含有检波点鬼波的记录，通过上下行波的分离可以用于OBS数据镜像偏移.
|
图 9 包含检波点鬼波的一次波估计结果 （a）原始共检波点道集；（b）检波点鬼波+一次波；（c）不含有检波点鬼波的多次波. Fig. 9 The result of estimated primaries with receiver ghost (a) Common receiver gather; (b) Receiver ghost + primaries; (c) Multiples without receiver ghost. |
文中介绍了压制多次波的反馈迭代模型及其衍生的三种方法：经典SRME方法、反数据域方法和EPSI方法.EPSI方法与其他反馈迭代法一样是一种完全数据驱动的方法，不需要地下介质的任何信息.但EPSI与SRME不同，EPSI方法不需要自适应相减步骤，在稀疏约束基础上直接反演估计一次波信息，从而保护了一次有效波反射信息；EPSI与反数据域方法不同，EPSI方法不需要考虑矩阵求逆的稳定性问题，采用矩阵乘法的迭代计算估计一次波.OBS数据采集具有自己的特点，与海底电缆数据不同，海底电缆数据一般检波点较多，可以只用海底电缆数据本身预测和估计一次波，而OBS数据检波点非常稀疏，本文充分考虑了OBS采集的特点，以EPSI方法为理论基础，提出了联合拖缆数据压制OBS多次波方法，利用拖缆数据的信息克服了稀疏检波点OBS多次波预测的困难.模拟算例验证了本文方法的有效性.
EPSI方法需要大量的矩阵乘法运算，计算量较大，但由于OBS数据量相对较少，计算量可以接受.另外，为了提高计算效率可以采用较少的迭代次数，利用公式（12）保守方法估计一次波，这样，在满足精度的要求下可以提高计算效率.本文采用时间域L-0范数稀疏约束的EPSI方法估计OBS一次波，还可以将其扩展到利用其他稀疏变换（如curvelet变换，seislet变换等）稀疏约束的EPSI方法，发展利用其他稀疏变换的EPSI方法压制OBS多次波是下一步的研究方向.
致谢感谢马继涛博士对SRME和反数据域压制多次波方法的有关讨论，感谢荷兰Delft理工大学的Verschuur博士的指导和国家留学基金委的资助.
| [1] | Moghaddam P, Libak A, Keers H, et al. Efficient and accurate modeling of ocean bottom seismometer data using reciprocity. Geophysics , 2012, 77(6): T211-T220. DOI:10.1190/geo2011-0498.1 |
| [2] | 阮爱国, 初凤友, 孟补. 海底天然气水合物地震研究方法及海底地震仪的应用. 天然气工业(地质与勘探) , 2007, 27(4): 46–48. Ruan A G, Chu F Y, Meng B. Seismic methods and OBS application in the study of oceanic gas hydrates. Natural Gas Industry (in Chinese) , 2007, 27(4): 46-48. |
| [3] | Dash R, Spence G, Hyndman R, et al. Wide-area imaging from OBS multiples. Geophysics , 2009, 74(6): Q41-Q47. DOI:10.1190/1.3223623 |
| [4] | Verschuur D J, Neumann E I. Integration of OBS data and surface data for OBS multiple removal. SEG Technical Program Expanded Abstracts, 1999: 1350-1353. |
| [5] | Pica A, Manin M, Granger P, et al. 3D SRME on OBS data using waveform multiple modeling. SEG Technical Program Expanded Abstracts, 2006: 2659-2663. |
| [6] | Weglein A B. Multiple attenuation: an overview of recent advances and the road ahead. The Leading Edge , 1999, 18(1): 40-44. DOI:10.1190/1.1438150 |
| [7] | 黄新武, 孙春岩, 牛滨华, 等. 基于数据一致性预测与压制自由表面多次波--理论研究与试处理. 地球物理学报 , 2005, 48(1): 173–180. Huang X W, Sun C Y, Niu B H, et al. Surface-related multiple prediction and suppression based on data-consistence: A theoretical study and test. Chinese J. Geophys. (in Chinese) , 2005, 48(1): 173-180. |
| [8] | Sava P, Guitton A. Multiple attenuation in the image space. Geophysics , 2005, 70(1): V10-V20. DOI:10.1190/1.1852789 |
| [9] | Verschuur D J, Berkhout A, Wapenaar C. Adaptive surface related multiple elimination. Geophysics , 1992, 57(9): 1166-1177. DOI:10.1190/1.1443330 |
| [10] | 李翔, 胡天跃. 逆散射级数法去除自由表面多次波. 地球物理学报 , 2009, 52(6): 1633–1640. Li X, Hu T Y. Surface-related multiple removal with inverse scattering series method. Chinese J. Geophys. (in Chinese) , 2009, 52(6): 1633-1640. |
| [11] | 薛亚茹, 陈小宏, 陆文凯. 压制多次波的正交多项式谱减法. 地球物理学报 , 2009, 52(3): 817–823. Xue Y R, Chen X H, Lu W K. Orthogonal polynomial spectrum subtraction for multiple attenuation. Chinese J. Geophys. (in Chinese) , 2009, 52(3): 817-823. |
| [12] | 李鹏, 刘伊克, 常旭, 等. 均衡拟多道匹配滤波法在波动方程法压制多次波中的应用. 地球物理学报 , 2007, 50(6): 1844–1853. Li P, Liu Y K, Chang X, et al. Application of the equipoise pseudomultichannel matching filter in multiple elimination using wave-equation method. Chinese J. Geophys. (in Chinese) , 2007, 50(6): 1844-1853. |
| [13] | Berkhout A J. Seismic processing in the inverse data space. Geophysics , 2006, 71(4): A29-A33. DOI:10.1190/1.2217727 |
| [14] | Ma J T, Sen M K, Chen X H. Free-surface multiple attenuation using inverse data processing in the coupled plane-wave domain. Geophysics , 2009, 74(4): V75-V81. DOI:10.1190/1.3137059 |
| [15] | 马继涛, SenM K, 陈小宏. 平面波域反数据处理压制多次波方法研究. 地球物理学报 , 2009, 52(3): 808–816. Ma J T, Sen M K, Chen X H. Multiple attenuation using inverse data processing in the plane-wave domain. Chinese J. Geophys. (in Chinese) , 2009, 52(3): 808-816. |
| [16] | van Borselen R G, Fokkema J T, van den Berg P M. Removal of surface-related wave phenomena-the marine case. Geophysics , 1996, 61(1): 202-210. DOI:10.1190/1.1443940 |
| [17] | Biersteker J. Magic: Shell's surface multiple attenuation technique. SEG Technical Program Expanded Abstracts, 2001: 1301-1304. |
| [18] | van Groenestijn G, Verschuur D. Estimating primaries by sparse inversion and application to near-offset data reconstruction. Geophysics , 2009, 74(3): A23-A28. DOI:10.1190/1.3111115 |
| [19] | 马继涛, SenM, 陈小宏, 等. 海底电缆多次波压制方法研究. 地球物理学报 , 2011, 54(11): 2960–2966. Ma J T, Sen M, Chen X H, et al. OBC multiple attenuation technique using SRME theory. Chinese J. Geophys. (in Chinese) , 2011, 54(11): 2960-2966. |
| [20] | Claerbout J F. Image estimation by example: Geophysical soundings image construction. Stanford Exploration Project, 2012, available on http://sep.stanford.edu/sep/prof/. |
| [21] | Lin T, Herrmann F. Estimating primaries by sparse inversion in a curvelet-like representation domain. The 73rd Conference and Exhibition, EAGE, Extended Abstracts, 2011. |
| [22] | van Groenestijn G, Verschuur D. Towards a new approach for primary estimation. SEG Technical Program Expanded Abstracts, 2008: 2487-2491. |
| [23] | Fomel S. Adaptive multiple subtraction using regularized nonstationary regression. Geophysics , 2009, 74: V25-V33. DOI:10.1190/1.3043447 |