CN112925023B - 全波场反演地震数据多次波压制方法 - Google Patents

Abstract

本发明涉及一种全波场反演地震数据多次波压制方法，该方法基于Tau‑p变换，利用双平面波域对地震数据的压缩特性，通过多次迭代反演实现地震数据多次波的压制和一次波的估计；本发明可减小计算量，无需地下介质信息，可以应对复杂的地下数据，且公式推导简单，结果有效。

Description

全波场反演地震数据多次波压制方法

技术领域

本发明涉及地球物理勘探技术领域，具体是关于一种全波场反演地震数据多次波压制方法。

背景技术

海洋地震勘探中，受到海水表面强反射系数的影响，地震数据中存在强烈的多次波干扰；多次波的存在使得地震波波场复杂化，对反射波的振幅和能量造成严重干扰，限制地震数据频带宽度，模糊对地下地质构造的认识，是地震数据处理和解释的严重障碍。

因此，研究多次波的压制机理与理论，对于提高地震处理水平和地震成像精度，丰富和发展地震资料处理理论有着重要的意义。

发明内容

针对上述问题，本发明的目的是提供一种全波场反演地震数据多次波压制方法，该方法基于Tau-p变换，利用全波场反演进行多次波压制和一次波估计的优势，可以一定程度上减小计算量；该方法无需地下介质任何信息，可以应对任意复杂的地下数据，且公式推导结果简单，计算过程简单有效。

为实现上述目的，本发明采取以下技术方案：

本发明所述的全波场反演地震数据多次波压制方法，包括如下步骤：

1)对于给定的时空域地震数据D(x_s,x_r,t)，对共炮点道集进行Tau-p变换，得到Tau-p域地震数据D_tp(x_s,p_r,τ)；

2)对变换得到的Tau-p域地震数据D_tp(x_s,p_r,τ)抽取共射线参数道集，并对共射线参数道集进行高分辨率Tau-p变换，得到双Tau-p域地震数据D_tp(p_s,p_r,τ)；

3)对变换得到的双Tau-p域地震数据D_tp(p_s,p_r,τ)做线性映射变换，得到双平面波域地震数据D_pl(p_d,p₀,τ)；

4)对线性映射变换得到的双平面波域地震数据D_pl(p_d,p₀,τ)，进行全波场反演求取一次波D₀(p_d,p₀,τ)，i次迭代后得到双平面波域一次波地震数据D_0，i(p_d,p₀,τ)；

5)对双平面波域一次波地震数据D_0，i(p_d,p₀,τ)进行反线性映射变换，得到双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)；

6)对双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)进行反变换，得到Tau-p域一次波地震数据D_0.i(x_s,p_r,τ)；

7)对Tau-p域一次波地震数据D_0.i(x_s,p_r,τ)进行反Tau-p变换，得到时空域一次波地震数据D_0,i(x_s,x_r,t)；

8)D_0,i(x_s,x_r,t)即经过i次迭代压制多次波后得到的时空域一次波数据。

所述的全波场反演地震数据多次波压制方法，优选地，所述步骤1)中，由时空域向Tau-p域变换所依据的公式为：

式中，x_r为地震数据检波器的偏移距，x_s为地震数据炮点位置坐标，p_r为检波器射线参数；t和τ为均为时间。

所述的全波场反演地震数据多次波压制方法，优选地，所述步骤2)中，对Tau-p域地震数据D_tp(x_s,p_r,τ)依据公式(2)进行高分辨率Tau-p变换，得到双Tau-p域地震数据D_tp(p_s,p_r,τ)：

式中，p_s为震源射线参数，x_s为地震数据炮点位置坐标，p_r为检波器射线参数；t和τ均为时间。

所述的全波场反演地震数据多次波压制方法，优选地，所述步骤3)中，对双Tau-p域地震数据D_tp(p_s,p_r,τ)依据公式(3)进行线性映射变换，得到双平面波域地震数据D_pl(p_d,p₀,τ)：

D_pl(p_d,p_o,τ)＝D_tp(p_r-p_s,p_r,τ) (3)

式中，p_s为震源射线参数，p_r为检波器射线参数；p_d为下行射线参数，p_o

为上行射线参数，τ为时间。

所述的全波场反演地震数据多次波压制方法，优选地，所述步骤4)中对双平面波域地震数据D_pl(p_d,p₀,τ)依据如下步骤进行全波场反演压制多次波，求取一次波D₀(p_d,p₀,τ)：

a.双平面波域中含有多次波的地震数据表示为：

D_pl＝(I+R)^-1D₀ (4)

式中，D_pl为含有多次波的双平面波域地震数据，R为反射系数，I为单位矩阵，D₀为不含多次波的一次波数据，由公式(5)得到：

D₀＝D_pl-AD_M (5)

式中，A＝S^-1，S为震源子波，通过多次波自适应相减得到，D_M为常规方法预测得到的多次波数据，且D₀满足以下关系：

D₀＝RA^-1 (6)

因此，D₀还表示为：

D₀＝D_pl+D₀AD_pl (7)

b.建立拟合残差函数V_i，表示为：

V_i＝D_pl+D_0,iAD_pl-D_0,i (8)

c.根据残差函数建立目标函数J：

对目标函数求取D₀的最速下降方向，即负梯度方向-ΔD₀，有：

式中Tr为矩阵的迹；

再按照计算后的更新方向更新D₀，有：

D_0,i+1＝D_0,i+αΔD_0,i (11)

其中，α为更新步长，且应满足目标函数J达到极小值条件，即V_i关于α的导数应为0；根据此条件求得α：

且令：

V＝D+D₀AD-D₀ K＝-ΔD₀(AD-I)

＝D-D₀(I-AD)＝ΔD₀(I-AD)

则式(12)依据化简为：

其中，V^HK＝K^HV，则：

式中Re为取计算结果的实部；

然后用求得的α更新每次迭代后的D_0，i；

d.由更新后得到的D_0，i与D_pl做时空域的二维褶积得到新的多次波数据D_M,i；并自适应相减得到算子A_i；将A_i,D_0,i代入式(10)；

e.重复步骤c，d，直到满足最大迭代次数i时，结束循环，得到完成i次迭代后的D_0，i。

所述的全波场反演地震数据多次波压制方法，优选地，所述步骤5)中，对双平面波域一次波地震数据D_0，i(p_d,p_o,τ)依据公式(15)进行反线性映射，得到双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)：

D_0,i(p_s,p_r,τ)＝D_0,i(p_d-p_o,p_o,τ) (15)。

所述的全波场反演地震数据多次波压制方法，优选地，所述步骤6)中，对D_0,i(p_s,p_r,τ)依据公式(16)进行双平面波域向Tau-p域数据的反变换，得到D_0,i(x_s,p_r,τ)：

所述的全波场反演地震数据多次波压制方法，优选地，所述步骤7)中，对Tau-p域数据D_0,i(x_s,p_r,τ)依据公式(17)进行反Tau-p变换，得到时空域数据D_0,i(x_s,x_r,t)：

本发明由于采取以上技术方案，其具有以下优点：

1、本发明基于Tau-p变换，利用双平面波域对地震数据的压缩特性，通过多次迭代反演实现地震数据多次波的压制和一次波的估计；

2、本发明在地震数据横向存在起伏的情况下，将地震数据变换到双平面波域可以有效对地震数据进行压缩，压缩后地震数据量减小，即可以利用全波场反演进行多次波压制和一次波估计的优势，也可以一定程度上减小计算量；

3.双平面波域全波场反演多次波压制和一次波估计，公式推导结果简单，计算过程简单有效；

4.本发明的多次波压制和一次波估计方法无需地下介质任何信息，可以应对任意复杂的地下数据。

附图说明

图1为本发明的整体流程示意图；

图2为在理论数据测试中双平面波域原始地震数据示意图；

图3为在理论数据测试中本发明方法对图2数据进行常规方法估计得到的双平面波域多次波数据示意图；

图4为在理论数据测试中本发明方法对图2数据进行迭代更新后得到的双平面波域多次波数据示意图；

图5为在理论数据测试中本发明方法最终估计得到的双平面波域一次波数据示意图；

图6为多次波压制前后的时空域地震数据示意图。

具体实施方式

以下将结合附图对本发明的较佳实施例进行详细说明，以便更清楚理解本发明的目的、特点和优点。应理解的是，附图所示的实施例并不是对本发明范围的限制，而只是为了说明本发明技术方案的实质精神。

如图1所示，本发明提供的全波场反演地震数据多次波压制方法，包括如下步骤：

1)对于给定的时空域地震数据D(x_s,x_r,t)，对共炮点道集进行Tau-p变换，得到Tau-p域地震数据D_tp(x_s,p_r,τ)；其中，由时空域向Tau-p域变换所依据的公式为：

式中，x_r为地震数据检波器的偏移距，x_s为地震数据炮点位置坐标，p_r为检波器射线参数；t为和τ为均为时间。

2)对变换得到的Tau-p域地震数据D_tp(x_s,p_r,τ)，抽取共射线参数道集，并依据公式(2)对共射线参数道集进行高分辨率Tau-p变换，得到双Tau-p域地震数据D_tp(p_s,p_r,τ)：

式中，p_s为震源射线参数，x_s为地震数据炮点位置坐标，p_r为检波器射线参数；t为变换前的时间；τ为变换后的时间。

3)对变换得到的双Tau-p域地震数据D_tp(p_s,p_r,τ)，依据公式(3)做线性映射，得到双平面波域地震数据D_pl(p_d,p_o,τ)：

D_pl(p_d,p_o,τ)＝D_tp(p_r-p_s,p_r,τ) (3)

式中，p_s为震源射线参数，p_r为检波器射线参数；p_d为下行射线参数，p_o

为上行射线参数。

4)对线性映射变换得到的双平面波域地震数据D_pl(p_d,p_o,τ)，依据以下步骤进行全波场反演求取一次波D₀(p_d,p_o,τ)，i次迭代后得到双平面波域一次波地震数据D_0，i(p_d,p_o,τ)：

a.双平面波域中含有多次波的地震数据可表示为：

D_pl＝(I+R)^-1D₀ (4)

式中，D_pl为含有多次波的双平面波域地震数据，R为反射系数，I为单位矩阵，D₀为不含多次波的一次波数据，由公式(5)得到：

D₀＝D_pl-AD_M (5)

式中，A为S^-1，S为震源子波，可通过多次波自适应相减得到，D_M为常规方法预测得到的多次波数据，且D₀满足以下关系：

D₀＝RA^-1 (6)

将式(6)代入式(4)，可得到D₀的另一种表达形式：

D₀＝(I+R)D_pl

D₀＝(I+D₀A)D_pl (7)

D₀＝D_pl+D₀AD_pl

b.通过式(7)用最小二乘拟合可建立拟合残差函数V_i，表示为：

V_i＝D_pl+D_0,iAD_pl-D_0,i (8)

c.根据残差函数建立目标函数J：

式中，ω为角频率；

对目标函数求取D_pl的最速下降方向，即负梯度方向-ΔD₀，有：

式中Tr为矩阵的迹；

再按照计算后的更新方向更新D₀，有：

D_0,i+1＝D_0,i+αΔD_0,i (11)

式中，α为更新步长，且应满足目标函数J达到极小值条件，即V_i关于α的导数应为0。根据此条件可求得α：

且令：

V＝D+D₀AD-D₀ K＝-ΔD₀(AD-I)

＝D-D₀(I-AD)＝ΔD₀(I-AD)

则等式依据可化简为：

其中，V^HK＝K^HV，则：

式中Re为取计算结果的实部；

然后用求得的α更新每次迭代后的D_0，i。

d.由更新后得到的D_0，i与D_pl做时空域的二维褶积得到新的多次波数据D_M,i，并自适应相减得到算子A_i；将A_i,D_0,i代入式(10)；

e.重复步骤c，d，直到满足最大迭代次数i时，结束循环，得到完成i次迭代后的D_0，i。

5)依据公式(15)对双平面波域一次波地震数据D_0，i(p_d,p_o,τ)进行反线性映射，得到双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)：

D_0,i(p_s,p_r,τ)＝D_0,i(p_d-p_o,p_o,τ) (15)

式中，p_s为震源射线参数；p_r为检波器射线参数；p_d为下行射线参数；p_o为上行射线参数；τ为时间。

6)依据公式(16)对双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)进行反变换，得到Tau-p域一次波地震数据D_0.i(x_s,p_r,τ)：

式中，x_s为地震数据炮点位置坐标；p_s为震源射线参数；p_r为检波器射线参数；t和τ均为时间。

7)依据公式(17)对Tau-p域一次波地震数据D_0,i(x_s,p_r,τ)进行反Tau-p变换，得到时空域一次波地震数据D_0,i(x_s,x_r,t)：

式中，x_s为地震数据炮点位置坐标；x_r为地震数据检波器的偏移距；p_r为检波器射线参数；t和τ均为时间。

8)D_0,i(x_s,x_r,t)为该方法经过i次迭代压制多次波后得到的时空域一次波数据。

如图2所示，图中地震数据采样点为750，采样率为4ms，经过变换到双平面波域，图中的数据伴随有较大程度的多次波干扰，导致数据显示较为混杂，且多次波与一次波形态一致，在时间域中较难去除。

如图3～4所示，在理论数据测试中，首先利用全波场反演的理论通过全波场数据预测得到其中的多次波信息，之后利用本发明算法对该多次波数据进行迭代更新。图3所示的为常规方法预估得到的多次波数据在双平面波域的显示，可以看到，图中多次波数据团的能量较为发散，不够聚焦，因此在往后的压制过程中会造成较大影响。图4为多次波数据经过4次自适应相减迭代更新后的双平面波域显示，与图3相比较，图中的能量团更为聚集，有收敛的趋势，利于多次波压制。

如图5所示，利用本发明算法对数据进行双平面波域多次波压制。由图5的压制结果可以看出，压制后的地震数据中多次波的影响被减弱。且构建的一次波数据中，具有很好的保幅性。

图6为多次波压制前后的时空域地震数据，可以更明显的看出方法在压制多次波的同时，具有较好的保幅性能。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (8)

1.一种全波场反演地震数据多次波压制方法，其特征在于，包括如下步骤：

1)对于给定的时空域地震数据D(x_s,x_r,t)，对共炮点道集进行Tau-p变换，得到Tau-p域地震数据D_tp(x_s,p_r,τ)；

2)对变换得到的Tau-p域地震数据D_tp(x_s,p_r,τ)抽取共射线参数道集，并对共射线参数道集进行高分辨率Tau-p变换，得到双Tau-p域地震数据D_tp(p_s,p_r,τ)；

3)对变换得到的双Tau-p域地震数据D_tp(p_s,p_r,τ)做线性映射变换，得到双平面波域地震数据D_pl(p_d,p₀,τ)；

4)对线性映射变换得到的双平面波域地震数据D_pl(p_d,p₀,τ)，进行全波场反演求取一次波D₀(p_d,p₀,τ)，i次迭代后得到双平面波域一次波地震数据D_0，i(p_d,p₀,τ)；

5)对双平面波域一次波地震数据D_0，i(p_d,p₀,τ)进行反线性映射变换，得到双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)；

6)对双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)进行反变换，得到Tau-p域一次波地震数据D_0.i(x_s,p_r,τ)；

7)对Tau-p域一次波地震数据D_0.i(x_s,p_r,τ)进行反Tau-p变换，得到时空域一次波地震数据D_0,i(x_s,x_r,t)；

8)D_0,i(x_s,x_r,t)即经过i次迭代压制多次波后得到的时空域一次波数据；

其中，x_s为地震数据炮点位置坐标；x_r为地震数据检波器的偏移距；t和τ均为时间；p_s为震源射线参数；p_r为检波器射线参数；p_d为下行射线参数；p_o为上行射线参数。

2.根据权利要求1所述的全波场反演地震数据多次波压制方法，其特征在于，所述步骤1)中，由时空域向Tau-p域变换所依据的公式为：

式中，x_r为地震数据检波器的偏移距，x_s为地震数据炮点位置坐标，p_r为检波器射线参数；t和τ为均为时间。

3.根据权利要求1所述的全波场反演地震数据多次波压制方法，其特征在于，所述步骤2)中，对Tau-p域地震数据D_tp(x_s,p_r,τ)依据公式(2)进行高分辨率Tau-p变换，得到双Tau-p域地震数据D_tp(p_s,p_r,τ)：

式中，p_s为震源射线参数，x_s为地震数据炮点位置坐标，p_r为检波器射线参数；t 和τ均为时间。

4.根据权利要求1所述的全波场反演地震数据多次波压制方法，其特征在于，所述步骤3)中，对双Tau-p域地震数据D_tp(p_s,p_r,τ)依据公式(3)进行线性映射变换，得到双平面波域地震数据D_pl(p_d,p₀,τ)：

D_pl(p_d,p_o,τ)＝D_tp(p_r-p_s,p_r,τ) (3)

式中，p_s为震源射线参数，p_r为检波器射线参数；p_d为下行射线参数，p_o为上行射线参数，τ为时间。

5.根据权利要求1所述的全波场反演地震数据多次波压制方法，其特征在于，所述步骤4)中对双平面波域地震数据D_pl(p_d,p₀,τ)依据如下步骤进行全波场反演压制多次波，求取一次波D₀(p_d,p₀,τ)，其中p_d为下行射线参数，p_o为上行射线参数，τ为时间：

a.双平面波域中含有多次波的地震数据表示为：

D_pl＝(I+R)^-1D₀ (4)

式中，D_pl为含有多次波的双平面波域地震数据，R为反射系数，I为单位矩阵，D₀为不含多次波的一次波数据，由公式(5)得到：

D₀＝D_pl-AD_M (5)

式中，A＝S^-1，S为震源子波，通过多次波自适应相减得到，D_M为常规方法预测得到的多次波数据，且D₀满足以下关系：

D₀＝RA^-1 (6)

因此，D₀还表示为：

D₀＝D_pl+D₀AD_pl (7)

b.建立拟合残差函数V_i，表示为：

V_i＝D_pl+D_0,iAD_pl-D_0,i (8)

c.根据残差函数建立目标函数J：

式中，ω为角频率；

对目标函数求取D_pl的最速下降方向，即负梯度方向-ΔD₀，有：

式中Tr为矩阵的迹；

再按照计算后的更新方向更新D₀，有：

D_0,i+1＝D_0,i+αΔD_0,i (11)

其中，α为更新步长，且应满足目标函数J达到极小值条件，即V_i关于α的导数应为0；根据此条件求得α：

且令：

则式(12)依据化简为：

其中，V^HK＝K^HV，则：

式中Re为取计算结果的实部；

然后用求得的α更新每次迭代后的D_0，i；

d.由更新后得到的D_0，i与D_pl做时空域的二维褶积得到新的多次波数据D_M,i；并自适应相减得到算子A_i；将A_i,D_0,i代入式(10)；

e.重复步骤c，d，直到满足最大迭代次数i时，结束循环，得到完成i次迭代后的D_0，i。

6.根据权利要求1所述的全波场反演地震数据多次波压制方法，其特征在于，所述步骤5)中，对双平面波域一次波地震数据D_0，i(p_d,p_o,τ)依据公式(15)进行反线性映射，得到双Tau-p域一次波地震数据D_0,i(p_s,p_r,τ)：

D_0,i(p_s,p_r,τ)＝D_0,i(p_d-p_o,p_o,τ) (15)

式中，p_s为震源射线参数；p_r为检波器射线参数；p_d为下行射线参数；p_o为上行射线参数；τ为时间。

7.根据权利要求1所述的全波场反演地震数据多次波压制方法，其特征在于，所述步骤6)中，对D_0,i(p_s,p_r,τ)依据公式(16)进行双平面波域向Tau-p域数据的反变换，得到D_0,i(x_s,p_r,τ)：

式中，x_s为地震数据炮点位置坐标；p_s为震源射线参数；p_r为检波器射线参数；t和τ均为时间。

8.根据权利要求1所述的全波场反演地震数据多次波压制方法，其特征在于，所述步骤7)中，对Tau-p域数据D_0,i(x_s,p_r,τ)依据公式(17)进行反Tau-p变换，得到时空域数据D_0,i(x_s,x_r,t)：

式中，x_s为地震数据炮点位置坐标；x_r为地震数据检波器的偏移距；p_r为检波器射线参数；t和τ均为时间。

Images