表层相关多次波压制算法在CMP道集中的应用

2020-06-10刘仕友马继涛孙万元

物探化探计算技术 2020年2期

刘仕友，马继涛，孙万元

(1. 中海石油(中国)有限公司湛江分公司，湛江 524057；2. 中国石油大学(北京) 地球物理学院，北京 102249)

0 引言

海洋地震勘探中，由于海水表面及海底处的波阻抗差值大，导致海水表面和海底的反射系数较高，从而使得海洋地震数据中存在非常严重的多次波问题。多次波的存在干扰了对有效波地识别，导致速度分析、叠前和叠后偏移的不准确，影响地震成像的真实性和可靠性，对地震解释工作也带来了很大地干扰,有效压制多次波对复杂油气藏区域地震数据的处理和解释有重要的意义。

海洋地震勘探中，海洋自由表面的强反射作用使得地震数据中表层相关的多次波能量占绝大部分，因此，针对海洋地震数据的表层相关多次波压制方法研究是多次波压制的主流,多次波的压制方法可以分为两大类[1]:

1)信号处理方法，是基于有效波和多次波之间不同的某种特性(如周期性和可分离性)。基于多次波和有效波之间的差异，可以设计一个滤波器对其进行滤除，因此这类方法统称为滤波法。但这类方法都是有一定的假设的，如“周期性”假设多次波具有周期性而没有有效波，“可分离性”假设通过一些变换可以将多次波和有效波分离开，进而通过切除的方式压制多次波。比较典型的滤波方法有预测反褶积(周期性)，FK滤波(可分离性)，共中心点叠加(可分离性)，抛物线Radon变换(可分离性)等。滤波法在地下构造比较简单的情况下是有效的，但复杂区域的地震数据有可能无法满足滤波法的假设(随着偏移距的增加，多次波不再具有周期性)；多次波和有效波在变换域的差异越来越小。在这些情况下进行多次波压制有可能出现压制不完全，或者损伤有效波能量的情况。

2)基于波动理论。该类方法以表层相关多次波衰减，即SRME最为知名。该方法由Berkhout[2]提出，他针对表层相关多次波和层间多次波，在考虑反子波的情况下，加入炮点和检波点性质算子的作用，利用多维反演算法初步提出了多次波压制的反馈迭代模型。Verschuur等[3-5]基于Berkhout的理论，采用数据驱动的方式，利用地震数据自身在时空域的二维褶积实现多次波的预测，之后利用自适应相减的手段将预测的多次波和数据中的多次波进行匹配相减，以期达到多次波压制的目的。该类方法不需要地下介质的任何信息，一经提出就引起了广大地球物理学者的兴趣，目前已经成为海洋地震数据表层相关多次波压制最为常用的算法之一。除此之外，Weglein[6]、Faqi Liu[7]等分别基于逆散射级数序列、恒定内插理论推导出了类似的多次波预测相减方法。

常规的SRME是一种基于波动理论的算法，该算法要求地震数据是全波场规则的地震数据，即地震数据由零偏移距至最大偏移距都要有等间距的数据覆盖，且炮点间隔和检波点间隔相等。实际数据无法满足方法要求，因此在利用SRME算法进行实际数据处理之前，或处理数据的同时，需要对数据进行规则化和插值外推处理。如果仅仅是二维地震数据，插值外推和规则化的运算量还可以接受，但若是三维地震数据，由于横测线方向数据极为稀疏，通过三维规则化处理得到地震数据存在着运算量巨大和数据质量不可靠等问题。Verschuur[3]和Kelamis等[8-9]提出了在地下介质横向起伏变化不大时，在CMP道集中应用SRME进行多次波压制的方法，该方法可以在一定程度上避免由于数据不规则和采样稀疏等带来的问题。

笔者对SRME应用在CMP道集中的算法进行讨论，该算法是一种表层相关多次波压制的高效算法，假设地下介质或者局部为一维介质，这样数据中的每个CMP道集可以视为一维介质情况下采集到的炮记录。对于真正的一维介质而言，采集到的所有炮记录都是相同的，因此CMP道集的SRME算法可以通过在频率波数域的简单向量相乘，而不是矩阵相乘的形式实现，这样可以大大减小算法对内存的需求，降低算法的复杂度，提高计算效率。

1 方法原理

首先从表层相关多次波压制的基本算法-反馈迭代理论阐述方法的原理。该理论认为，如果没有自由表层的反馈效应，地表接收到的地震数据P0(z0)可以表示为：

P0(z0)=D(z0)X0(z0，z0)S(z0)

(1)

其中：z0为震源和检波器的垂向位置坐标；S(z0)为震源函数；D(z0)为检波点性质矩阵；X0(z0,z0)为有效波反射矩阵。

若式(1)中加入自由表层的反馈效应，则接收到的地震数据可以表示为：

P(z0)=P0(z0)+[P0(z0)A(ω)]P0(z0)+

[P0(z0)A(ω)]2P0(z0)+…

(2)

(3)

称其为表层算子。

由式(2)可以推导出有效波波场(含层间多次波)，即表层相关多次波压制的表达式：

P0(z0)=P(z0)-A(ω)P2(z0)+

A2(ω)P3(z0)-A3(ω)P4(z0)+…

(4)

如果地下介质是一维的，采集到的炮记录都是相同的，则表层相关多次波压制可以利用式(5)进行。

P0(kx,ω;z0)=P(kx,ω;z0)-

A(ω)P2(kx,ω;z0)+

A2(ω)P3(kx,ω;z0)-

A3(ω)P4(kx,ω;z0)+…

(5)

其中:kx、ω分别为波数和频率；P0(kx,ω;z0)是炮记录对应的单频切片，在一维情况下，因为所有的炮记录都是相同的，所以P0(kx,ω;z0)对应的是P(z0)矩阵的一列。式(4)中的矩阵相乘可以由式(5)中的向量相乘代替。

在实际情况下，地下介质并不完全满足一维的假设，所有的炮记录不是完全相同的，但是如果用CMP道集代替炮记录，会更好地满足一维介质的假设。在产生多次波的地层，即介质上部构造起伏较为平缓的情况下，CMP道集的多次波压制可以取得很好的结果。

需要注意的是，在CMP道集进行表层相关多次波压制的步骤和炮集预测类似，都要进行数据的插值外推的处理，即如果CMP道集最小偏移距不为零，需要利用一定的方法对缺失的数据进行插值外推处理。在动校正后的道集进行插值，之后再反动校正，一般就能满足方法的需求。如果CMP道集的地震道个数太少，可以将相邻的几个CMP道集合并起来，形成CMP超道集，对该超道集进行规则化、插值和多次波预测处理。将CMP道集的多次波预测出来之后，用炮集的SRME方法类似，需要用自适应相减的方法将预测的多次波从原始数据中减去。

因此，利用本文方法进行表层相关多次波压制的步骤如下：

1)读取地震数据D，得到CMP道集的偏移距信息，如地震数据近偏移距不为零，则利用抛物线Radon变换对数据进行近炮检距外推处理，得到数据Dexp。

2)对步骤1)得到的CMP道集Dexp，沿时间和偏移距方向做傅里叶变换，得到频率-波数域的地震数据Dexp_fk。

3)对步骤2)得到的频率-波数域地震数据进行向量相乘运算，得到频率-波数域的多次波记录Mexp_fk。

4)对步骤3)得到的频率波数域多次波记录Mexp_fk进行反傅里叶变换，得到时间-空间域的多次波记录Mexp；若数据进行过外推处理，则将外推的地震道记录进行剔除，得到对应的多次波记录M。

基于最小能量准则，将D和M进行自适应相减，得到多次波压制后的结果P。

2 数据示例

我们分别用一个平层速度模型、一个横向平缓起伏的速度模型生成的模拟地震数据和一个实际海洋数据,对方法的可行性和有效性进行了验证。

2.1 模拟数据算例

第一个模型是平层速度模型，如图1所示。

图1 平层速度模型Fig.1 Flat layer model

在图1中从2 000 m至8 000 m，每隔20 m放置一个震源，并在震源左侧放置101个检波器接收，模拟出201个炮记录。对模拟出的共炮点道集，根据数据震源和检波点坐标计算出每道的CMP数值，并抽取数据共CMP道集。在此只展示其中的3个CMP道集，因为此数据为平层速度模型数据，因此3个CMP道集都是相同的(图2)。

图2 平层数据所对应的3个CMP道集Fig.2 The 3 CMP gathers corresponding to the flat layer data

为了对比常规SRME方法和CMP道集多次波预测方法的方法原理，在此给出了常规SRME方法和本文方法的数据矩阵，分别为频率-空间域和频率-波数域的数据矩阵，如图3所示。

由图3可以看出，常规SRME方法的数据矩阵，每一列对应的为一个共炮点道集，其元素都是相同的；而频率波数域的数据矩阵，其能量几乎都聚焦在了反对角线上，因此数据矩阵的相乘可以用反对角线方向对应元素的点乘代替，可以大大减小计算量。

图3 数据矩阵Fig.3 Data matrices(a)平层数据常规SRME方法数据矩阵；(b)对应的频率-波数域数据矩阵

图4 平层数据CMP道集所对应的多次波结果Fig.4 The multiple result corresponding to the CMP gathers of flat layer data

图5 平层数据多次波压制后的CMP道集Fig.5 The CMP gathers after multiple attenuation of the flat layer data

图6 共偏移距道集Fig.6 Common offset gather(a)多次波压制前;(b)多次波压制后

遵循CMP道集多次波预测的方法步骤，对每个CMP道集做二维傅里叶变换，并对变换后的频率-波数域数据进行向量点乘运算，得到该CMP道集对应的多次波。图4展示的是图2的三个CMP道集所对应的三个多次波记录，由图4可以看出，预测的多次波结果与原数据中的多次波在旅行时和同相轴形态上都较为吻合。将预测出的多次波从原数据中自适应减去，得到图5即多次波压制后的结果。

多次波压制前后的共偏移距剖面如图6所示，可以看出，对于平层数据，由于所有位置处的CMP道集都是相同的，利用本文方法进行多次波预测和压制，和常规SRME方法进行多次波压制的效果是一样的。

利用横向平缓起伏变化的速度模型(图7)正演模拟出的共炮点道集，抽取其CMP道集(图8)进行多次波预测的工作。同样给出了常规SRME方法频率-空间域数据矩阵和本文方法的频率-波数域数据矩阵(图9)。在图9中可以看出，由于此数据横向有了起伏变化，每一列元素所对应的共炮点道集也有了些微的变化，但频率-波数域的数据矩阵能量仍然集中在反对角线元素上，因此可以同样利用向量的点乘运算预测出数据的多次波。

图7 一个横向平缓起伏的速度模型Fig.7 A horizontal gentle varying velocity model

图8 横向平缓起伏数据的三个CMP道集Fig.8 3 CMP gathers of the horizontal gentle varying velocity model

图9 横向平缓起伏数据的数据矩阵Fig.9 The data matrix of the horizontal gently varying model(a)频率-空间数据矩阵；(b)频率-波数数据矩阵

图10 起伏数据预测的多次波Fig.10 Predicted multiple of the horizontal gently varying model

图11 自适应相减之后的结果Fig.11 The result after adaptive subtraction

图12 横向平缓起伏数据共偏移距剖面Fig.12 Common offset gather of the horizontal gently varying model(a)多次波压制前；(b)多次波压制后

对图8中的CMP道集沿时间和空间方向分别进行傅里叶变换，得到频率-波数域的CMP道集，并将其进行点乘运算，得到频率-波数域的多次波，之后对其进行反傅里叶变换，得到时间-空间域的多次波(图10)；将预测得到的多次波从原数据中自适应减去，得到多次波压制后的结果(图11)。从图11中可以看出，单个CMP道集的多次波压制可以取得较好的效果，但在远偏移距处，由于地形起伏使得预测的多次波带有一定的误差，从而导致多次波压制有一定的残余，此数据所对应的多次波压制前后的共偏移距剖面如图12所示。从图12中也可以看出，除了在地形起伏变化较为剧烈的个别区域之外，对该数据进行多次波压制也可以取得较好的效果。

图13 实际数据的CMP道集Fig.13 CMP gathers of real data(a)原数据；(b)预测的多次波；(c)自适应相减之后的结果

图14 实际数据数据矩阵Fig.14 Data matrix of real data(a)频率-空间数据矩阵；(b) 频率-波数域数据矩阵

2.2 实际数据算例

利用墨西哥湾的一个实际数据对方法进行了进一步验证。图13(a)给出的是墨西哥湾数据的一个CMP道集，该CMP道集的近偏移距为﹣100 m，首先利用抛物线Radon变换对该数据的近偏移距进行了插值外推处理，之后对插值外推的数据沿空间和时间进行傅里叶变换，之后对该频率-波数域数据进行点乘并做反傅里叶变换，得到该CMP道集对应的多次波(图13(b))，将该多次波从原CMP道集中自适应减去，得到多次波压制后的结果(图13(c))，可以看出对于实际数据，本文方法仍然可以取得较好的效果。

从图14中可以看出，即使本实际数据比较复杂，其频率-波数域数据矩阵的能量仍然大部分聚焦在了反对角线附近。因此利用两次傅里叶变换之后数据的点乘运算进行多次波预测对于实际数据仍然是有一定的意义的。

图15为实际数据对应的共偏移距剖面，图16为地形起伏区域及中部平缓区域多次波压制后剖面的局部放大图。从图15和图16中可以明显看出，除了在地形起伏变化较大，不满足本方法假设的区域(箭头所指)，本文方法均能够较好地预测并压制多次波。