基于分治算法的谷值因子斜缆数据鬼波压制方法*

崔明晓 王守东 王豆豆 齐子威

(1.中国石油大学(北京)油气资源与探测国家重点实验室 北京 102249； 2.中国石油大学(北京)海洋石油勘探国家工程实验室 北京 102249 )

海洋地震数据拖缆采集时通常是将采集拖缆放置于海水面以下，由于海水面的反射，检波器除了接收到地下反射波外，还会接收到反射波到达海水面后向下产生的虚反射，该虚反射也被称为鬼波[1]。由于鬼波陷波效应的影响，使得地震数据频带变窄，从而降低了地震剖面的分辨率[2]，增加了地震地质解释的困难，因此必须在地震资料采集和处理时设法消除鬼波的影响[3]。以往常规海洋地震数据采集一般采用等深度缆采集技术[4]。为了压制鬼波，国际上先后发展了OBC[5]、上下缆[6-7]和双检等采集方式[8]，提出了上下缆采集数据处理方法[9]和双检采集数据处理方法[10]，如PGS公司提出了双传感器采集系统中上行波场分离的方法压制虚反射[11],这些方法对于消除鬼波影响，保持地震数据的低频和高频信息都有较好的效果。

近几年，CGG公司推出了变深度缆采集技术[12]。针对变深度缆采集数据，Sablon等[13]提出了镜像偏移压制虚反射的方法，该方法能够有效压制虚反射和拓宽地震频带[14]；Soubaras等[15]提出了联合反褶积压制虚反射的方法，该方法在墨西哥湾取得了良好的应用效果；Wang等[16]提出了一种在F-XY域利用实际炮集和镜像炮集联合反褶积来压制鬼波的算法；许自强 等[17]提出了海上变深度缆数据最优化压制虚反射的方法；叶林 等[18]提出了峰值因子最大化压制鬼波方法。本文推导了利用镜像变深度缆地震数据联合反褶积去鬼波的算法，针对该方法中因鬼波延迟时间不准确而导致不能得到最优压制效果的问题，提出了通过谷值因子来判断最精确鬼波延迟时间的技术思路，并进行了模型试算及实际数据处理，均取得了较好的去鬼波效果。

1 方法原理

1.1 数据最优化联合反褶积去鬼波方法

对于沉放深度为zi处的变深度缆所接收到的实际波场p(x,zi,t)，可以分解为上行波场u(x,zi,t)和下行波场d(x,zi,t)之和的形式，见式(1)，其中上行波场为地震波自激发后向下传播经地下界面反射直接到达检波器的波场，下行波场则是地下反射波先到达海水面后经海水面反射再到达检波器的波场。

p(x,zi,t)=u(x,zi,t)+d(x,zi,t)

(1)

式(1)中：x为偏移距变量;t为时间变量。

把式(1)变换到频率域,可得

P(x,zi,ω)=U(x,zi,ω)+D(x,zi,ω)

(2)

波场传播示意图见图1，其中深度为z1=zr处的实际变深度缆某接收点接收到的总波场为p(x,z1,t)，上行波场为u(x,z1,t)，下行波场为d(x,z1,t)；首先在深度为2zr的位置定义一个假想的海水界面，然后放置一个关于此界面与实际接收点对称的镜像接收点，其深度为z2，且z2=3zr，则假设镜像接收点接收到的总波场为p(x,z2,t)，上行波场为u(x,z2,t)，下行波场为d(x,z2,t)。由式(1)、(2)可得

P(x,z1,ω)=U(x,z1,ω)+D(x,z1,ω)

(3)

P(x,z2,ω)=U(x,z2,ω)+D(x,z2,ω)

(4)

同时使得镜像接收点接收到的下行波场等于实际拖缆接收到的上行波场，即

D(x,z2,ω)=U(x,z1,ω)

(5)

图1 波场传播示意图Fig.1 Schematic of wave propagation

通常情况下，我们可以把海水面的反射系数取值为-1，此时就可以把下行波场写成上行波场和波场延拓算子的形式,即

(6)

将式(5)代入式(6)，可得

U(x,z2,ω)=-U(x,z1,ω)ejk(z2-z1 )

(7)

将式(5)～(7)式代入式(3)、(4)，可得

P(x,z1,ω)=U(x,z1,ω)[1-e-jk(z2-z1 )]

(8)

P(x,z2,ω)=U(x,z1,ω)[1-ejk(z2-z1)]

(9)

从式(8)、(9)可以看出，实际拖缆和镜像接收点处接收到的总波场都可以用实际拖缆处接收到的上行波场和一个算子计算得到，联立式(8)、(9)，可以得出镜像记录是在实际原始记录的基础上作用了一个算子S(x,ω)，即

P(x,z2,ω)=P(x,z1,ω)S(x,ω)

(10)

其中，镜像算子S(x,ω)为

S(x,ω)=-ejk(z2-z1)

(11)

因为原始记录和镜像记录又可以用实际拖缆z1处的上行波场U(x,z1,ω)和鬼波算子G(x,ω)表示，即

(12)

其中

(13)

式(13)中：Δt为上行波与下行波之间的延迟时间。从而频率空间域的镜像记录可以表示为

P(x,z2,ω)=P(x,z1,ω)(-ejωΔt)

(14)

这样即可推导出用原始记录和带有鬼波延迟时间Δt的算子来求取镜像记录的公式。

利用式(14)可以由原始地震记录计算出镜像地震记录，压制鬼波可以通过原始记录与镜像记录联合反褶积法来实现，具体方法可以通过求解如下最小平方问题来实现，即

U(x,ω)=

(15)

通过求解式(15)可以得到上行波场频率域数据U(x,ω)，再通过反傅里叶变换即可得到时间域上行波场u(x,t)，也就是压制鬼波后的地震数据。

1.2 谷值因子提出

以上推导中，镜像记录求解时需要已知鬼波延迟时间Δt，但实际中Δt的影响因素有很多，包括检波器深度、海水面、地震波传播角度等，因此，往往不能事先得出精确的Δt，使得计算得到的镜像记录不够准确，进而导致通过数据联合反褶积方法得不到最优的鬼波压制效果。为了求得最精确的鬼波延迟时间Δt，本文采用了谷值因子的思路。

当鬼波延迟时间不准确时，原始记录和镜像记录做单道反褶积得到的波形会产生震荡，与精确结果相比波形复杂。为了刻画波形的形态，用得到的单道地震数据最大值和最小值之差和该道记录的均方根作商，假设代入在第i个鬼波延迟时间求得的鬼波数据为di(x,z,t)，由此定义了1个参数V的计算方法，即

(16)

其中，rms[di(x,z,t)]为该下行波地震道数据的均方根，则代入的每一个Δt均可得到一个V，把所有得到的V值中最小的一个称为谷值因子。谷值因子所对应的Δt即为该地震道最精确的鬼波延迟时间，用该延迟时间得到的去鬼波数据是最精确的。因此,搜索最精确的鬼波延迟时间Δt的过程即为搜索谷值因子的过程。谷值因子可通过建立如下目标函数来确定,即

s=minV

(17)

1.3 分治算法搜索谷值因子

根据以上方法，通过搜索谷值因子来确定最精确鬼波延迟时间。为了提高搜索速度，本文采用分治算法搜索，其思想是：当在搜索一组数据中的最小值时，如果参与搜索的数据非常多，根本无法直接找出其中的最小值，使得计算过程非常复杂，可以把参与搜索的数据分成几组，即把一个大的问题分解成几个子问题[19]，如果这些子问题可以直接得出结果，也就是分解产生的数据能够直接判断最小值，则问题就很容易解决了；如果被划分成的每个数据组依然不容易判断最小值，可以把子问题继续划分为更小的子问题，以此类推继续划分数据，直至可以判断最小值，并得出谷值因子为止。

通过分治算法可以有效搜索出谷值因子，得出最精确的鬼波延迟时间，从而得出最优的去鬼波结果。基于分治算法搜索谷值因子来消除鬼波的过程可以分成以下4个步骤：

1) 确定鬼波延迟时间的搜索范围和步长，把该范围内的每一个鬼波延迟时间代入式(14)，计算对应的镜像数据。

2) 由原始数据和镜像数据，应用数据最优化联合反褶积方法求得上行波场和下行波场，利用时间域下行波场分别计算对应的参数V。

3) 在得到的所有参数V中通过分治算法快速寻找最小值Vmin，即为谷值因子，其对应的Δt即为最精确的鬼波延迟时间。

4) 根据最精确的鬼波延迟时间计算最准确的镜像记录，进而得到最精确的去鬼波地震数据。

2 模型试算

2.1 验证谷值因子的有效性

为了验证谷值因子对应着最佳鬼波延迟时间的正确性，在单层介质的情况下给定一个主频为50 Hz的单道信号(图2a)，再以50 ms为标准延迟时间给定一个鬼波信号(图2b)，则2个信号相加即可认为是带有鬼波的原始信号(图2c)。设置鬼波延迟时间的搜索范围Δt∈(25 ms，75 ms)，搜索步长Δτ=1 ms，则可知在此范围内共可搜索50次，分别得到Δt从25 ms到75 ms的鬼波记录和每一个鬼波波记录对应的参数V。图2d、e、f、g、h分别为Δt=30、45、50、55、75 ms时得到的鬼波信号，可以看出图2f即Δt=50 ms时的鬼波信号与真实的鬼波信号最接近。

把搜索50次中每次得到的V值计算出来，并绘制变化曲线图(图3)，发现当Δt=50 ms，即Δt为精确的鬼波延迟时间时所对应的V值最小，即为谷值因子，这时提取的鬼波正好最接近真实鬼波。

图2单层介质某单道信号不同鬼波延迟时间得到的鬼波波形

Fig.2Ghostwaveformwithdifferentghostwavedelayofsinglechannelsignalinsinglelayermedia

图3 单层介质参数V值随鬼波延迟时间Δt的变化曲线Fig.3 Curve of the V change with the ghost wave delay Δt in single layer media

以上结果验证了谷值因子在单层介质信号上的有效性。再对多层介质的情况进行验证，给出1个三层介质下的一次波某单道信号(图4a)，然后再以鬼波延迟时间为50 ms得到鬼波信号(图4b)，2个信号相加即可得到带有鬼波的原始信号(图4c)。与单道信号一样，在搜索范围Δt∈(25 ms，75 ms)以步长为1 ms依次代入其中的每一个鬼波延迟时间，分别得到下行波场，并由其计算每一个对应的参数V。图4d、e、f、g、h分别为Δt=30、45、50、55、75 ms时得到的鬼波信号，可以看出图4f即Δt=50 ms时的鬼波信号与真实的鬼波信号最接近。

图4 多层介质某单道信号不同鬼波延迟时间得到的鬼波波形Fig.4 Ghost waveform with different ghost wave delay of single channel signal in multilayer media

把每次搜索对应的参数V值随延迟时间Δt变化曲线画出来(图5)，可以看出多层介质情况下最精确的鬼波延迟时间对应的参数V也是最小的，即多层介质中谷值因子同样可以作为最优鬼波延迟时间Δt的判断标准。

图5 多层介质参数V值随鬼波延迟时间Δt变化曲线Fig.5 Curve of the V change with the ghost wave delay Δt in multilayer media

由此可见，谷值因子所对应的鬼波延迟时间Δt即是最精确的鬼波延迟时间，因此，为了得到最优化的去鬼波结果，计算搜索谷值因子是一种比较有效的方法。

2.2 简单模型试算

建立的简单速度模型如图6a所示。横向采样间隔1.5 m，纵向采样间隔1.5 m，水层速度1 500 m/s，采用70 Hz雷克子波。声波方程有限差分正演模拟，单边接收，每炮645道。震源放于左侧水面上，检波器深度4.5～50.0 m，图6b为缆深变化，图6c为第100道V值大小随鬼波延迟时间的变化曲线。图7a为带有鬼波的原始地震记录，图7b为对应的镜像地震记录，图7c为运用本文方法去除鬼波后的地震记录。图8a、b分别为对应的去鬼波前、后的频谱图。由简单模型的地震记录和频谱图可以看出，运用本文方法对带有鬼波的原始记录进行去鬼波处理，在剖面上可以有效的压制鬼波，提高地震数据的分辨率；在频率域可以消除陷波效应，提高频带宽度。

图6 简单速度模型(a)、缆深变化(b)和第100道V值随鬼波延迟时间Δt变化曲线(c)Fig.6 Simple speed model(a),cable depth change(b)and the curve of V change with the ghost wave delay Δt in the 100th trace(c)

图7 简单模型原始地震记录局部(a)、镜像记录局部(b)和去鬼波后地震记录局部(c)Fig.7 Part of original seismograms(a),mirror seismograms(b) and seismograms after deghost(c)of the simple model

图8 简单模型处理前(a)与处理后(b)频谱图对比Fig.8 Spectrum before(a)and after(b)processing of the simple model

2.3 复杂模型试算

建立的复杂速度模型如图9a所示。横向采样间隔1.5 m，纵向采样间隔1.5 m，水层速度1 500 m/s，采用70 Hz雷克子波。声波方程有限差分正演模拟，单边接收，每炮1 345道。震源放于左侧水面上，检波器深度4.5～50.0 m，图9b为缆深变化，图9c为第150道V值随鬼波延迟时间的变化曲线。图10a为带有鬼波的原始地震记录，图10b为最优鬼波延迟时间计算得到的镜像地震记录，图10c为运用本文方法去除鬼波后的地震记录。图11a、b分别为对应的去鬼波前、后的频谱图。由复杂模型的地震记录和频谱图可以看出，本文方法对去除复杂模型地震记录的鬼波和消除陷波效应也具有良好的效果。

图9 复杂速度模型(a)、缆深变化(b)和第150道V值随鬼波延迟时间Δt变化曲线(c)Fig.9 Complex speed model(a),cable depth change(b)and the curve of V change with the ghost wave delay Δt in the 150th trace(c)

图10 复杂模型原始地震记录镜局部(a)、镜像记录局部(b)和去鬼波后地震记录局部(c)Fig.10 Part of original seismograms(a),mirror seismograms(b)and seismograms after deghost(c)of the complex model

图11 复杂模型处理前(a)与处理后(b)频谱图对比Fig.11 Spectrum before(a)and after(b)processing of the complex model

3 实际数据处理

用本文方法对某海上地震数据进行处理。图12a为原始地震记录，图12b为镜像地震记录，图12c为用本文方法进行处理后得到的地震剖面，图13a、b、c分别为实际资料相应的局部显示。由实际地震数据的处理结果可以看出，本文压制鬼波方法取得了良好的效果。

图12 实际资料原始地震记录(a)、镜像地震记录(b)和压制鬼波后地震记录(c)Fig.12 Original seismograms(a),mirror seismograms(b) and seismograms after deghost(c)of the actual data

图13 实际资料原始局部显示(a)、镜像局部显示(b)和去鬼波后局部显示(c)Fig.13 Part of original seismograms(a),mirror seismograms(b) and seismograms after deghost(c)of the actual data

4 结论

1) 在应用数据最优化联合反褶积方法压制斜缆数据鬼波时，通过搜索谷值因子可以得到最精确的鬼波延迟时间，从而解决鬼波延迟时间不准确的问题，进而得到更优化的斜缆数据鬼波压制效果。

2) 在搜索谷值因子时，分治算法可以作为一种有效的搜索方法，通过该方法可以减少搜索次数，并准确地找到谷值因子，进而确定鬼波延迟时间，得到去鬼波的地震记录。

3) 应用本文方法对模拟数据和实际数据进行处理均得到了较好的去鬼波效果，在剖面上提高了地震资料的分辨率，拓宽了地震频带。

[1] 陆敬安，伍忠良，曾宪军.海洋地震勘探中地震波、鬼波综合效应分析与应用[J].海洋技术，2006,25(4):76-78，98.

LU Jing’an,WU Zhongliang,ZENG Xianjun.The synthesized effect of seismic wave and ghost reflection and its application in marine seismic survey[J].Ocean Technology,2006,25(4):76-78，98.

[2] 李套山，高洪强，刘振夏，等.虚反射信息的采集及其形成机制、频率响应的理论讨论[J].石油物探，1997,36(4)：38-44.

LI Taoshan,GAO Hongqiang,LIU Zhenxia,et al.The acquisition of ghost information and the study of its generation and frequency response[J].GPP,1997,36(4):38-44.

[3] 刘建磊，王修田.海上地震虚反射相位剔除法反褶积[J].海洋地质与第四纪地质，2002,22(4):117-121.

LIU Jianlei,WANG Xiutian.Deconvolution by phase spectrum deghosting in marine seismic survey[J].Marine Geology & Quaternary Geology,2002,22(4):117-121.

[4] 王振峰，李绪深，孙志鹏，等．琼东南盆地深水区油气成藏条件和勘探潜力[J].中国海上油气，2011,23(1):7-13.

WANG Zhenfeng,LI Xushen,SUN Zhipeng,et al.Hydrocarbon accumulation conditions and exploration potential in the deep-water region，Qiongdongnan basin[J].China Offshore Oil and Gas,2011,23(1):7-13.

[5] 王守君.海底电缆地震技术优势及在中国近海的应用效果[J].中国海上油气，2012,24(2):9-12.

WANG Shoujun.Technical advantages of OBC seismic survey and its application effects offshore China[J].China Offshore Oil and Gas,2012,24(2):9-12.

[6] HILL D,COMBED L,BACON J.Over/under acquisition and data processing:the next quantum leap in seismic technology[J].First Break,2006,24(6):81-95.

[7] MOLDOVEANU N,COMBEE J.Over/under towed-streamer acquisition:a method to extend seismic bandwidth to both higher and lower frequencies[J].The Leading Edge,2007,26(1):41-58.

[8] TENGHAMN R,VAAGE S,BORRESEN C.A dual-sensor,towed marine streamer:its viable implementation and initial results[J].SEG Expanded Abstracts,2007(1):989.

[9] 刘春成，刘志斌，顾汉明，等.利用上下缆合并算子确定海上上下缆采集的最优沉放深度组合[J].石油物探，2013,52(6)：623-629.

LIU Chuncheng,LIU Zhibin,GU Hanming,et al.The determination of optimal sinking depths of over/under streamers in offshore survey by merge operator[J].GPP,2013,52(6):623-629.

[10] SOLLNER W,BROX E,WIDMAIER M,et al.Surface related multiple suppression in dual senor towed steamer data[J].SEG Technical Program Expanded Abstracts,2007,26:2540-2544.

[11] CARLSON C,LONG A,SOLLNER W,et al.Increased resolution and penetration from a towed dual-sensor streamer[J].First Break,2007,25(12):71-77.

[12] SOUBARAS R,DOWLE R.Variable-depth streamer,a broadband marine solution[J].First Break,2010,28:89-96.

[13] SABLON R,RUSSIER D,ZURITA O,et al.Multiple attenuation for variable-depth streamer data:from deep to shallow water[C].81st Annual Meeting,SEG,Expanded Abstracts,2011:3505-3509.

[14] SOUBARAS R.Pre-stack deghosting for variable-depth streamer data[C].82nd Annual Meeting,SEG,Expanded Abstracts,2012:780-789.

[15] SOUBARAS R.Deghosting by joint deconvolution of a migration and a mirror migration[C].80th Annual Meeting,SEG,Expanded Abstracts,2010:3406-3410.

[16] WANG P,PENG C.Premigration deghosting for marine streamer data using a bootstrap approach[C].82nd Annual Meeting,SEG,Expanded Abstracts,2012:1-5.

[17] 许自强，方中于，顾汉明，等.海上变深度缆数据最优化压制鬼波方法及其应用[J].石油物探，2015,54(4):404-413.

XU Ziqiang,FANG Zhongyu,GU Hanming,et al,The application of optimal deghosting algorithm on marine variable-depth streamer data[J].GPP,2015,54(4):404-413.

[18] 叶林,韩立国.基于峰值因子最大化的鬼波压制方法[J].世界地质,2016,35(4):1102-1107.

YE Lin,HAN Liguo.Deghosting method based on maximization of peak factor[J].Global Geology,2016,35(4):1102-1107.

[19] 李绍华，王建新，马振宇，等.基于加权分治技术的setpacking精确算法[J].小型微型计算机系统，2010,31(6)：1180-1184.

LI Shaohua,WANG Jianxin,MA Zhenyu,et al.Exact algorithm for set packing based on measure and conquer[J].Chinese Computer Systems,2010,31(6):1180-1184.