CN101334483A - 一种在地震数据处理中衰减瑞雷波散射噪声的方法 - Google Patents

Abstract

提供一种在地震数据处理中衰减Rayleigh波散射噪声的方法，本方法只在主要包含散射波场的频段进行散射噪声衰减处理，不影响其它频段的原始波场特征，是一种相对保真的去噪方法。本方法综合利用线性波场变换、小波变换和α－trimmed滤波技术提取Rayleigh波散射入射波场并计算子波，然后通过散射波场反演技术，实现散射波场的最佳估计，再从原始波场中减去估计的散射噪声，从而有效地提高了地震资料的信噪比。

Description

一种在地震数据处理中衰减瑞雷波散射噪声的方法

技术领域

本发明涉及一种地震数据处理方法，具体涉及地球物理勘探中地震数据处理领域中衰减瑞雷(Rayleigh)波散射噪声的方法。

背景技术

地震勘探是一种利用人工地震技术激发地震波，研究地震波在地层中传播的情况，以查明地下的地质构造或岩性特征，为寻找油气藏或其它勘探目的服务的一种物探方法。在现有勘探石油的各种物探方法中，地震勘探是最有效的方法。

地震勘探的生产工作，基本上可以分地震资料采集、地震资料处理、地震资料解释三个环节：

第一环节是地震资料采集。这个阶段的主要任务是在地质工作和其它物探工作初步确定的油气勘探有利区带，按照一定的观测系统布置测线，使用爆炸震源或可控震源激发地震波，并用地震仪把地震波场记录下来。

第二环节是地震资料处理。这个阶段的主要任务是根据地震波的传播理论，利用计算机和相应地震资料处理软件，对野外采集的地震数据进行各种加工处理，以获得能够反映地下地层结构的“地震剖面图”和反映地下地层岩性变化的地震波振幅、频率、传播速度等信息。

第三个环节是地震资料解释。经过计算机处理得到的地震剖面，虽然已经能够反映地下地质构造的一些特点，但是由于地下情况的复杂性，地震剖面上的许多现象，既可能反映地下的真实情况，也有可能有某些假象。地震解释的主要任务就是运用地震波传播理论和石油地质学原理，综合地质、钻井和其它物探资料，对地下的构造或岩性作出正确的判断，并最终提出钻井井位。

随着勘探开发工作的不断深入，为了满足地震解释及地震储层预测要求，在地震资料处理时有高信噪比、高分辨率和高保真度的要求原则。高信噪比是地震资料高分辨率和高保真度的基础，由于野外采集得到的地震波场除了反射波之外，还包括面波、声波、折射波、多次波、各种次生散射、地面振动源干扰及地面微震等，在成像处理前要尽可能去除掉这些干扰。

地震勘探中遇到的Rayleigh面波噪声，它的特点是频率低，一般为几赫兹到30赫兹；速度低，一般在100m/s到1000m/s之间，以200m/s到500m/s最为常见，目前在地震处理中，现有的Rayleigh面波噪声压制方法主要为滤波法。滤波法可以有效地去除面波干扰，但是滤波法最大的缺陷就是在去除噪声的同时也损失了部分有效波。

目前在提高地震资料信噪比的处理中还存在诸多技术难题，散射波干扰问题就是其中之一。尤其在复杂构造区，由近地表非均匀性引起的Rayleigh波散射噪声，严重地影响了地震资料的质量，提高这些地震资料的信噪比已成为迫切需要。

发明内容

为了解决上面背景技术中存在的问题，本发明提出了一种在地震数据处理中衰减Rayleigh波散射噪声的方法。该方法综合利用线性波场变换、小波变换和α-trimmed滤波技术提取Rayleigh波散射的入射波场并计算子波，然后通过散射波场反演得到散射波场的最佳估计，再从原始波场中减去估计的散射波场，从而有效地提高了地震资料的信噪比。

依据本发明，提供一种在地震数据处理中衰减Rayleigh波散射噪声的方法，该方法包括以下步骤：

(1)在地质工作和其它物探工作初步确定的油气勘探有利区带上，按照一定的观测系统布置测线和炮位，使用爆炸震源或可控震源激发地震波，并用地震仪把地震波场记录下来；

(2)对记载有地震波场的地震资料按炮号和道号从小到大进行排序；对排好序的炮集数据做振幅补偿；对补偿后的数据做振幅增益控制，并保存增益系数；提取主要包含散射波场频段的信号，并把剩余部分保存起来；提取入射波场和计算子波；进行Rayleigh波散射衰减处理；利用α-trimmed滤波技术去除散射衰减处理结果中极低频段的噪声信号；将去除极低频段噪声后的波场和提取主要包含散射波场频段信号时的剩余波场相加；对相加后的波场进行振幅恢复处理，使地震数据在散射噪声衰减处理前后振幅保持一致；对经过散射噪声衰减的地震数据进行其它常规处理。

其中(2)具体采用如下步骤：

(a)将地震资料按炮号和道号从小到大进行排序；

(b)对排序好的炮集数据做常规振幅补偿处理。

(c)对经过振幅补偿的炮集数据做常规振幅增益控制处理，并保存每一采样点的振幅增益系数。振幅增益系数是这样计算的：对于同一采样点，用增益后的振幅值除以增益前的振幅值，就得到了该采样点的增益系数，如果增益前该采样点振幅值为零，则令该采样点振幅增益系数为零。

(d)选取其中一炮记录，利用α-trimmed滤波技术提取主要包含散射波场频段的信号，记为波场S，将剩余部分记为波场R；

(e)在波场S中，综合利用线性波场变换、小波变换和α-trimmed滤波技术提取入射波场并计算子波；

(f)对波场S进行Rayleigh波散射衰减处理；

(g)利用α-trimmed滤波技术去除(f)步结果中的极低频段的噪声信号；

(h)将(g)步输出的波场与(d)步输出的波场R相加；

(i)对(h)步相加结果进行振幅恢复处理。对于(h)步输出数据的每一个采样点，用该采样点的振幅值除以该采样点的振幅增益系数，对于增益系数为零的采样点，直接令该采样点的振幅值为零；至此完成单炮Rayleigh波散射噪声衰减处理过程；

(j)对每一炮均重复(d)到(i)各步操作，完成整个地震资料的近地表Rayleigh波散射噪声衰减处理；

(k)将去噪后的炮记录进行静校正、动校正、反褶积、叠加、偏移等常规处理，获得改进后的叠加、偏移剖面数据，通过构造、岩性解释来提供精确可靠的井位。

其中为了避免在反演求取散射波场过程中，由于振幅不平衡而出现的奇异性，对原始记录做振幅补偿和振幅增益控制，并保存振幅增益系数，在去除散射噪声以后对振幅进行恢复，使散射噪声衰减处理前后振幅保持一致。

其中为了不改变主要有效波频段成分的波场特征，使用α-trimmed滤波技术提取主要含Rayleigh波散射噪声频段信号，在衰减散射噪声时避开了主要为有效波的频段。

其中在上述(e)步中：综合利用线性波场变换、小波变换和α-trimmed滤波技术提取入射波场及计算子波，具体实现步骤如下：

(41)根据提供的面波视速度，利用线性正变换将Rayleigh波面波同相轴变换呈水平状；

如图1a所示，L1，L2，L3，L4为四条曲线，它们与纵轴的交点分别为a_1，0，a_2，0，a_3，0，a_4，0，图1b为经过线性变换后的结果，L1，L2，L3，L4为四条曲线变为互相平行的直线。

(42)采用小波变换对线性变换后的波场进行分解；

(43)在存在Rayleigh面波的各个小波尺度域利用α-trimmed滤波技术，将具有水平特征的面波分离出来；

(44)利用小波反变换对面波数据重构；

(45)对重构数据进行线性反变换得到入射波场；

(46)在提取的入射波场中选取近炮点的某道，根据给定的子波长度，在消除了Rayleigh面波波速引起的时移后，作为计算子波。

其中在(f)步中：进行散射波场衰减时，使用的时间域离散化散射波场正演公式和反演公式如下：

时间域总的散射波场正演计算公式为：

${\stackrel{\‾}{v}}_3(x,x_s)=v_3^{\mathrm{inc}}(x,x_s)+v_3^{\mathrm{sc}}(x,x_s)$

其中：

$v_3^{\mathrm{sc}}(x,x_s)=\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}\left[\mathrm{shape}(v_3^{\mathrm{inc}}\left(x\right|x^{\′})*v_3^{\mathrm{inc}}\left(x^{\′}\right|x_s))\right]\mathrm{\ξ}\left(r^{\′}\right)\mathrm{\ΔA}\left(r^{\′}\right)$

时间域散射波反演公式为：

$\left|\right|v_3(x,x_s)-[v_3^{\mathrm{inc}}(x,x_s)+\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}[\mathrm{shape}(v_3^{\mathrm{inc}}\left(x\right|x^{\′})*v_3^{\mathrm{inc}}\left(x^{\′}\right|x_s))\left]\mathrm{\ξ}\left(r^{\′}\right)\mathrm{\ΔA}\left(r^{\′}\right)\right]\left|\right|=\underset{\mathrm{\ξ}\left(r^{\′}\right)}{\min }$

以上公式中：x为检波点位置，x_s为震源的位置，v₃(x，x_s)为时域总散射波场，，v₃ ^inc(x，x_s)为时域入射波场，v₃ ^sc(x，x_s)为由入射场激励的散射场，v₃(x，x_s)为实际炮集资料。shape(～)表示子波反褶积运算，

表示微分运算，^*表示褶积，v₃ ^inc(x|x′)和v₃ ^inc(x′|x_s)为网格点x′上的Green速度张量，ξ(r′)为网格点上的阻抗差函数，||～||表示取2范数，r′为炮点或检波点到网格点的距离，

表示对所有r′求和，ΔA(r′)为离散化的散射点单元面积，为关于ξ(r′)的反演目标函数。

其中在上述(f)步中使用共轭梯度法进行时间域的Rayleigh波散射噪声衰减处理，具体实现步骤如下：

(61)从炮集资料获取测线及其排列等相关信息；

(62)确定浅层散射网格，并给网格上定义待求的波阻抗差函数分布；

(63)根据上面所述的散射波反演公式，建立关于波阻抗差的目标函数，应用共轭梯度法求解目标函数，得到波阻抗差函数最优解；

(64)利用求解的最佳波阻抗差函数，根据上述的总散射场正演计算公式，计算得到总的散射波场；

(65)从波场S中减去总的散射波场。

优选地，在(g)步中：使用α-trimmed滤波技术去除(f)步输出波场中的极低频噪声。

使用所述地震数据处理中衰减Rayleigh波散射噪声的方法，有利于衰减复杂构造区由于地表非均匀性引起的Rayleigh波散射噪声，提高地震数据信噪比。并且该方法还具有以下明显优点：在对Rayleigh波散射噪声进行衰减的过程中避开了主要为有效波的频段，不改变主要有效波频段的原始信号特征；在振幅增益时记录了增益系数，Rayleigh波散射噪声衰减后对振幅作进一步的恢复处理，既保证反演算法的稳定性，又保证处理前后振幅一致性。该方法是一种相对保真的去噪方法。

附图简要说明

图1a为线性正变换前的面波示意图；

图1b为线性正变换后的面波示意图；

图2为原始单炮记录示意图；

图3为振幅补偿后的波场示意图；

图4a为振幅增益控制处理后的波场示意图；

图4b为增益系数示意图；

图5a为从图4a中提取的主要含散射噪声频段波场示意图；

图5b为提取后剩余波场示意图；

图6a为从图5a中提取的入射波场示意图；

图6b为计算的子波示意图；

图7a为反演得到的散射波场示意图；

图7b为图5a对应波场与图7a对应波场相减后的波场示意图；

图8a为去除的极低频噪声示意图；

图8b为去除极低频噪声后剩余波场示意图；

图9为图8b对应波场与图5b对应波场相加后波场示意图；

图10为振幅恢复处理后的波场示意图；

图11a为Rayleigh波散射噪声压制前的初叠加剖面图(3600-5000ms)；

图11b为Rayleigh波散射噪声压制后的初叠加剖面图(3600-5000ms)。

具体实施方式

在本发明中为了提高地震资料的信噪比，提出了一种在地震数据处理中衰减Rayleigh波散射噪声的方法。该方法首先综合利用线性波场变换、小波变换和α-trimmed滤波技术提取Rayleigh波散射的入射波场并计算子波，然后通过散射波反演得到散射波场的最佳估计，再从原始波场中减去估计的散射波场，从而有效地提高了地震资料的信噪比。

在地震数据处理中衰减Rayleigh波散射噪声的方法包括以下步骤：

(1)在地质工作和其它物探工作初步确定的油气勘探有利区带，按照一定的观测系统布置测线和炮位，使用爆炸震源或可控震源激发地震波，并用地震仪把地震波场记录下来；

(2)对记载有地震波场的地震资料按炮号和道号从小到大进行排序；对排序好的炮集数据做振幅补偿；对补偿后的数据做振幅增益控制，并保存增益系数；提取主要包含散射波场频段的信号，并把剩余部分保存起来；提取入射波场和计算子波；进行Rayleigh波散射衰减处理；利用α-trimmed滤波技术去除散射衰减处理结果中极低频段的噪声信号；将去除极低频段噪声后的波场和提取主要包含散射波场频段信号时的剩余波场相加；对相加后的波场进行振幅恢复处理，使地震数据在散射噪声衰减处理前后振幅保持一致；对经过散射噪声衰减的地震数据进行其它常规处理。

其中(2)具体采用如下步骤：

(a)将地震资料按炮号和道号从小到大进行排序；

(b)对排序好的炮集数据做常规振幅补偿处理。

(c)对经过振幅补偿的炮集数据做常规振幅增益控制处理，并保存每一采样点的振幅增益系数。振幅增益系数是这样计算的：对于同一采样点，用增益后的振幅值除以增益前的振幅值，就得到了该采样点的增益系数，如果增益前该采样点振幅值为零，则令该采样点振幅增益系数为零。

(d)选取其中一炮记录，利用α-trimmed滤波技术提取主要包含散射波场频段的信号，记为波场S，将剩余部分记为波场R；

(e)在波场S中，综合利用线性波场变换、小波变换和α-trimmed滤波技术提取入射波场并计算子波；

(f)对波场S进行Rayleigh波散射衰减处理；

(g)利用α-trimmed滤波技术去除(f)步结果中的极低频段的噪声信号；

(h)将(g)步输出的波场与(d)步输出的波场R相加；

(i)对(h)步相加结果进行振幅恢复处理。对于(h)步输出数据的每一个采杆点，用该采样点的振幅值除以该采样点的振幅增益系数，对于增益系数为零的采样点，直接令该采样点的振幅值为零。至此完成单炮Rayleigh波散射噪声衰减处理过程；

(j)对每一炮均重复(d)到(i)各步操作，完成整个地震资料的近地表Rayleigh波散射噪声衰减处理；

(k)将去噪后的炮记录进行静校正、动校正、反褶积、叠加、偏移等常规处理，获得改进后的叠加、偏移剖面，通过构造、岩性解释提供精确可靠的井位。

其中为了避免在反演求取散射波场过程中，由于振幅不平衡而出现的奇异性，对原始记录做振幅补偿和振幅增益控制，并保存振幅增益系数，在去除散射噪声以后对振幅进行恢复，使散射噪声衰减处理前后振幅保持一致。

其中为了不改变主要有效波频段成分的波场特征，使用α-trimmed滤波技术提取主要含Rayleigh波散射噪声频段信号，在衰减散射噪声时避开了主要为有效波的频段。

其中在上述(e)步中：综合利用线性波场变换、小波变换和α-trimmed滤波技术提取入射波场及计算子波，具体实现步骤如下：

(41)根据提供的面波视速度，利用线性正变换将Rayleigh波面波同相轴变换呈水平状；

如图1a所示，L1，L2，L3，L4为四条曲线，它们与纵轴的交点分别为a_1，0，a_2，0，a_3，0，a_4，0，图1b为经过线性变换后的结果，L1，L2，L3，L4为四条曲线变为互相平行的直线。

(42)采用小波变换对线性变换后的波场进行分解；

(43)在存在Rayleigh面波的各个小波尺度域利用α-trimmed滤波技术，将具有水平特征的面波分离出来；

(44)利用小波反变换对面波数据重构；

(45)对重构数据进行线性反变换得到入射波场；

(46)在提取的入射波场中选取近炮点的某道，根据给定的子波长度，在消除了Rayleigh面波波速引起的时移后，作为计算子波。

其中在(f)步中：进行散射波场衰减时，使用的时间域离散化散射波场正演公式和反演公式如下：

时间域总的散射波场正演计算公式为：

${\stackrel{\‾}{v}}_3(x,x_s)=v_3^{\mathrm{inc}}(x,x_s)+v_3^{\mathrm{sc}}(x,x_s)$

其中：

$v_3^{\mathrm{sc}}(x,x_s)=\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}\left[\mathrm{shape}(v_3^{\mathrm{inc}}\left(x\right|x^{\′})*v_3^{\mathrm{inc}}\left(x^{\′}\right|x_s))\right]\mathrm{\ξ}\left(r^{\′}\right)\mathrm{\ΔA}\left(r^{\′}\right)$

时间域散射波反演公式为：

$\left|\right|v_3(x,x_s)-[v_3^{\mathrm{inc}}(x,x_s)+\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}[\mathrm{shape}(v_3^{\mathrm{inc}}\left(x\right|x^{\′})*v_3^{\mathrm{inc}}\left(x^{\′}\right|x_s))]\mathrm{\ξ}\left(r^{\′}\right)\mathrm{\ΔA}\left(r^{\′}\right)\left]\right||=\underset{\mathrm{\ξ}\left(r^{\′}\right)}{\min }$

以上公式中：x为检波点位置，x_s为震源的位置，v₃(x，x_s)为时域总散射波场，，v₃ ^inc(x，x_s)为时域入射波场，v₃ ^sc(x，x_s)为由入射场激励的散射场，v₃(x，x_s)为实际炮集资料。shape(～)表示子波反褶积运算，

表示微分运算，^*表示褶积，v₃ ^inc(x|x′)和v₃ ^inc(x′|x_s)为网格点x′上的Green速度张量，ξ(r′)为网格点上的阻抗差函数，||～||表示取2范数，r′为炮点或检波点到网格点的距离，

表示对所有r′求和，ΔA(r′)为离散化的散射点单元面积，

为关于ξ(r′)的反演目标函数。

其中在上述(f)步中使用共轭梯度法进行时间域的Rayleigh波散射噪声衰减处理，具体实现步骤如下：

(61)从炮集资料获取测线及其排列等相关信息；

(62)确定浅层散射网格，并给网格上定义待求的波阻抗差函数分布；

(63)根据上面所述的散射波反演公式，建立关于波阻抗差的目标函数，应用共轭梯度法求解目标函数，得到波阻抗差函数最优解；

(64)利用求解的最佳波阻抗差函数，根据上述的总散射场正演计算公式，计算得到总的散射波场；

(65)从波场S中减去总的散射波场。

另外，在(g)步中使用α-trimmed滤波技术去除(f)步输出波场中的极低频噪声。

为了进一步说明本发明，下面结合实例和附图以实际单炮数据为例说明本发明的步骤和过程：

第一步，对地震资料做预处理，将炮记录按炮号和道号从小到大进行排序。图2为原始单炮记录。这一步是针对非炮集资料进行的预处理，如果已经是按照炮号和道号排好顺序的炮集资料，就可以直接进行第二步处理。

第二步，对炮第一步生成的数据做振幅补偿。

由于地震波在传播的过程中，随着传播距离不断增大，波前面也不断增大，单位波前面上的振幅能量就越来越小，因此需要对振幅做几何扩散补偿，使深层和浅层单位波前面上的振幅能量基本保持一致。在补偿时一般用距离(均方根速度乘以时间)作为补偿因子。在此均方根速度设为常速1600m/s，图3为经过振幅补偿的单炮记录示意图。

第三步，对第二步生成的数据进行振幅增益控制，并保存增益系数。

振幅增益控制的作用就是要保证整个记录上最终的平均振幅值是均匀的。由于振幅增益对振幅具有破坏作用，所以要记录下增益系数，以便以后对振幅作恢复处理。振幅增益控制时在给定时窗内进行，在此时窗长度选取0.5ms。图4a为增益后的单炮记录示意图，图4b为记录的增益系数示意图。

第四步，利用α-trimmed滤波技术在增益后的波场中提取主要含散射噪声频段的信号。

为了不改变有效波频段的波场特征，使用α-trimmed滤波技术提取主要包含散射波场频段的信号。在此α-trimmed滤波裁剪比例取0.3，滤波统计窗口长度为21个采样点。图5a为提取的主要含散射噪声波场示意图，图5b为提取后剩余波场示意图。

第五步，综合利用线性波场变换，小波变换和α-trimmed滤波提取入射波场及计算子波。提取入射波场及计算子波的步骤如下：

(1)线性波场变换。

在炮集记录剖面上估计面波速度范围，在此估计给出三组面波速度，三组面波速度下限分别为1830m/s，405m/s，275m/s，速度上限分别为1850m/s，415m/s，285m/s，子波长度均为220ms。根据估计的面波速度计算出面波到达时间，利用计算的面波到达时间进行线性正变换，把面波同相轴校平。

(2)利用小波分解对线性变换后的数据进行分解。

利用小波变换可以实现对地震数据的精细描述，小波变换有很多种算法，在此选用Mallat算法将输入数据变换到不同尺度域，小波分解级别为7级。

(3)在存在面波的尺度域利用α-trimmed滤波技术提取Rayleigh波面波。

由于面波具有低频特性，所以经过小波变换后，面波一般只存在于小波时窗较宽对应的尺度域内。所以在实施α-trimmed滤波时，只在这几个尺度域进行即可。在此α-trimmed滤波裁剪比例取0.1，滤波统计窗口长度为31个采样点。

(4)利用小波反变换重构信号。

在重构的过程中只利用上面提取面波时对应的尺度域。在此选取小波重构开始级别为5，小波重构终止级别为7。

(5)对重构数据进行线性反变换得到入射波场。

图6a是经过线性反变换后得到的入射波场示意图。

(6)计算子波。

在第(5)步提取的入射波场中，选取距离炮点最近的道，根据炮检距和估计的面波速度，计算出时移量，把时移量作为子波起始位置，再根据估计的子波长度220ms，在选取道中截取这一部分作为子波。图6b是计算的子波示意图。

第六步，进行Reilaygh散射衰减处理。

散射衰减处理时首先要根据散射波反演公式建立关于波阻抗差的目标函数，再利用共轭梯度法求解目标函数，得到最佳的波阻抗差函数，最后利用散射波正演技术得到散射波场。具体实现步骤如下：

(1)从炮集资料获取测线及其排列等相关信息；

(2)确定浅层散射网格，并在网格上定义待求的阻抗差函数分布ξ；

(3)求解阻抗差函数ξ；

将第五步中提取的入射波场及计算的子波带入下面近地表散射场反演公式建立关于阻抗差ξ的目标函数，并利用共轭梯度法对目标函数进行求解。

散射波场反演公式如下：

$\left|\right|v_3(x,x_s)-[v_3^{\mathrm{inc}}(x,x_s)+\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}[\mathrm{shape}(v_3^{\mathrm{inc}}\left(x\right|x^{\′})*v_3^{\mathrm{inc}}(x^{\′}|x_s))\left]\mathrm{\ξ}(r^{\′})\mathrm{\ΔA}(r^{\′})\right]\left|\right|=\underset{\mathrm{\ξ}\left(r^{\′}\right)}{\min }$

(4)计算浅层散射波场。

将第五步中提取的入射波场及计算的子波以及(3)中求取的阻抗差函数都带入近地表散射波正演公式中，得到浅层的散射波场。图7a为得到的浅层散射波场示意图。

近地表散射波正演公式如下：

${\stackrel{\‾}{v}}_3(x,x_s)=v_3^{\mathrm{inc}}(x,x_s)+v_3^{\mathrm{sc}}(x,x_s)$

其中：

$v_3^{\mathrm{sc}}(x,x_s)=\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}\left[\mathrm{shape}(v_3^{\mathrm{inc}}\left(x\right|x^{\′})*v_3^{\mathrm{inc}}\left(x^{\′}\right|x_s))\right]\mathrm{\ξ}\left(r^{\′}\right)\mathrm{\ΔA}\left(r^{\′}\right)$

(5)用图5a所对应得数据减去图7a所对应的数据，图7b为散射噪声衰减处理后的示意图。

其中(3)和(4)公式中所涉及的参数意义如下：

||～||表示取2范数，x为检波点位置，x_s为震源的位置，v₃(x，x_s)为图5a对应的实际资料，v₃ ^inc(x，x_s)为第五步中提取的时域入射波场，r′为炮点、检波点到网格点的x′距离，

表示对所有r′求和，shape(～)表示子波反褶积运算，

表示微分运算，^*表示褶积，v₃ ^inc(x|x′)和v₃ ^inc(x′|x_s)为网格点x′上的Green速度张量，ξ(r′)为网格点上的阻抗差函数，ΔA(r′)为离散化的散射点单元面积，

为关于ξ(r′)反演目标函数。v₃(x，x_s)为模拟的时域总散射波场，v₃ ^sc(x，x_s)为由入射场激励的散射场。

第七步，去除第六步输出结果中的极低频噪声。

在地震处理中一般认为极低频段信号中很少含有效波成分，因此为了得到更好的去噪效果，利用α-trimmed滤波技术去除极低频部分噪声。在此α-trimmed滤波裁剪比例取为0.1，滤波统计窗口长度为41个采样点。图8a为去除的极低频噪声示意图，图8b为去除极低频噪声后剩余波场示意图，

第八步，将图8b对应波场与图5b对应的波场相加。图9为相加后的波场示意图。

第九步，对相加后的波场进行振幅增益恢复。

将图9对应的数据除以第三步保存的增益系数。图10为经过增益恢复的示意图。经过以上九步操作，就完成了单炮近地表Rayleigh波散射噪声衰减过程。

第十步，重复第四步到第九步操作，完成所有炮集的近地表Rayleigh波散射噪声衰减。

第十一步，处理人员可根据实际需要对去噪后的炮集数据进行静校正、动校正、反褶积、叠加、偏移等常规处理，获得改进后的叠加、偏移剖面，通过构造、岩性解释最终提供精确可靠的井位。图11a为Rayleigh散射噪声压制前的初叠加剖面图(3600-5000ms)，图11bRayleigh散射噪声压制后的初叠加剖面图(3600-5000ms)。从图11b中可以看出，Rayleigh波散射噪声得到较好的压制。

通过上面表述，可以看出使用所提供的地震数据处理中衰减Rayleigh波散射噪声的方法，有利于衰减复杂构造区由于地表非均匀性引起的Rayleigh波散射噪声，提高地震数据信噪比。并且该方法还具有以下明显优点：对Rayleigh波散射噪声进行衰减的过程中避开了主要为有效波的频段，不改变有效波频段的原始信号特征；在振幅增益时记录了增益系数，Rayleigh波散射噪声衰减处理后对振幅作进一步的恢复处理，既保证反演算法的稳定性，又保证处理前后有效波振幅一致性。该方法是一种相对保真的去噪方法。

进一步地，尽管已经清楚详细地描述了本发明提出的在地震数据处理中衰减Rayleigh波散射噪声的方法，并且本发明的优选实施例详细描述并解释了本发明，但是本领域普通的技术人员可以理解，在不背离所附权利要求定义的本发明的精神和范围的情况下，可以在形式和细节中做出多种修改。

Claims (7)

1、一种在地震数据处理中衰减Rayleigh波散射噪声的方法，该方法包括以下步骤：

(1)在地质工作和其它物探工作初步确定的油气勘探有利区带上，按照一定的观测系统布置测线和炮位，使用爆炸震源或可控震源激发地震波，并用地震仪把地震波场记录下来；

(2)对记载有地震波场的地震资料按炮号和道号从小到大进行排序；对排序好的炮集数据做振幅补偿；对补偿后的数据做振幅增益控制，并保存增益系数；提取主要包含散射波场频段的信号，并把剩余部分保存起来；提取入射波场并计算子波；进行Rayleigh波散射衰减处理；利用α-trimmed滤波技术去除散射衰减处理结果中极低频段的噪声信号；将去除极低频段噪声后的波场和提取主要包含散射波场频段信号时的剩余波场相加；对相加后的波场进行振幅恢复处理，使地震数据在散射噪声衰减处理前后振幅保持一致；对经过散射噪声衰减的地震数据进行其它常规处理；

其中(2)具体采用如下步骤：

(a)将地震资料按炮号和道号从小到大排序；

(b)对排序好的炮集数据做常规振幅补偿处理；

(c)对经过振幅补偿的炮集数据做常规振幅增益控制处理，并保存每一采样点的振幅增益系数。振幅增益系数是这样计算的：对于同一采样点，用增益后的振幅值除以增益前的振幅值，就得到了该采样点的增益系数，如果增益前该采样点振幅值为零，则令该采样点振幅增益系数为零。

(d)选取其中一炮记录，利用α-trimmed滤波技术提取主要包含散射波场频段的信号，记为波场S，将剩余部分记为波场R；

(e)在波场S中，综合利用线性波场变换、小波变换和α-trimmed滤波技术提取入射波场和计算子波；

(f)对波场S进行Rayleigh波散射衰减处理；

(g)利用α-trimmed滤波技术去除(f)步结果中的极低频段的噪声信号；

(h)将(g)步输出的波场与(d)步输出的波场R相加；

(i)对(h)步输出波场进行振幅恢复处理。对于(h)步输出波场的每一个采样点，用该采样点的振幅值除以该采样点的振幅增益系数，对于增益系数为零的采样点，直接令该采样点的振幅值为零；至此完成单炮Rayleigh波散射噪声衰减处理过程；

(j)对每一炮均重复(d)到(i)各步操作，完成整个地震资料的近地表Rayleigh波散射噪声衰减处理；

(k)将去噪后的炮记录进行静校正、动校正、反褶积、叠加、偏移的常规处理，获得改进后的叠加、偏移剖面数据，通过构造、岩性解释来提供精确可靠的井位。

2.如权利要求1所述的在地震数据处理中衰减Rayleigh波散射噪声的方法，其中为了避免在(f)步反演求取散射波场过程中，由于振幅不平衡而出现的奇异性，对原始记录做振幅补偿和振幅增益控制，并保存振幅增益系数，在去除散射噪声以后对处理后的振幅进行恢复，保证处理前后振幅一致。

3.如权利要求1所述的在地震数据处理中衰减Rayleigh波散射噪声的方法，其中为了不改变有效波频段的波场特征，使用α-trimmed滤波技术提取主要包含散射波场频段的信号，在衰减散射噪声时避开了主要为有效波的频段。

4.如权利要求1所述的在地震数据处理中衰减Rayleigh波散射噪声的方法，其中在(e)步中：综合利用线性波场变换、小波变换和α-trimmed滤波技术提取入射波场并计算子波，实现步骤如下：

(41)根据提供的面波视速度，利用线性正变换将Rayleigh面波同相轴变换呈水平状；

(42)采用小波变换对线性变换后的波场进行分解；

(43)在存在Rayleigh面波的各个小波尺度域利用α-trimmed滤波技术，将具有水平特征的面波分离出来；

(44)利用小波反变换对面波数据重构；

(45)对重构数据进行线性反变换得到入射波场；

(46)在提取的入射波场中选取近炮点的某道，根据给定的子波长度，在消除了Rayleigh面波波速引起的时移后，作为计算子波。

5.如权利要求1所述的在地震数据处理中衰减Rayleigh波散射噪声的方法，其中在(f)步中，进行散射波场衰减时，使用的时间域离散化散射波场正演公式和反演公式如下：

时间域总的散射波场正演计算公式为：

${\stackrel{\‾}{v}}_3(x,x_s)=v_3^{\mathrm{inc}}(x,x_s)+v_3^{\mathrm{sc}}(x,x_s)$

其中：

$v_3^{\mathrm{sc}}(x,x_s)=\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}\[\mathrm{shape}\left(v_3^{\mathrm{inc}}{\left(x\right|x^{\′})*v}_3^{\mathrm{inc}}\left(x^{\′}\right|x_s)\right)\]\mathrm{\ξ}\left(r^{\′}\right)\mathrm{\ΔA}\left(r^{\′}\right)$

时间域散射波反演公式为：

$\left|\right|v_3(x,x_s)-\[v_3^{\mathrm{inc}}(x,x_s)+\underset{r^{\′}}{\mathrm{\Σ}}\frac{d}{\mathrm{dt}}\[\mathrm{shape}\left(v_3^{\mathrm{inc}}{\left(x\right|x^{\′})*v}_3^{\mathrm{inc}}\left(x^{\′}\right|x_s)\right)\]\mathrm{\ξ}\left(r^{\′}\right)\mathrm{\ΔA}\left(r^{\′}\right)\]\left|\right|=\underset{\mathrm{\ξ}\left(r^{\′}\right)}{\min }$

以上公式中：x为检波点位置，x_s为震源的位置，v₃(x，x_s)为时域总散射波场，v₃ ^inc(x，x_s)为时域入射波场，v₃ ^sc(x，x_s)为由入射场激励的散射场，v₃(x，x_s)为实际炮集资料。shape(～)表示子波反褶积运算，

表示微分运算，*表示褶积，v₃ ^inc(x|x′)和v₃ ^inc(x′|x_s)为网格点x′上的Green速度张量，ξ(r′)为网格点上的阻抗差函数，||～||表示取2范数，r′为炮点或检波点到网格点的距离，

表示对所有r′求和，ΔA(r′)为离散化的散射点单元面积，

为关于ξ(r′)的反演目标函数。

6.如权利要求1所述的在地震数据处理中衰减散射噪声的方法，其中在(f)步中使用共轭梯度法进行时间域的Rayleigh波散射噪声衰减处理，具体实现步骤如下：

(61)从炮集资料获取测线及其排列等相关信息；

(62)确定浅层散射网格，并给网格上定义待求的波阻抗差函数分布；

(63)根据权利要求(5)中所述的散射波反演公式，建立关于波阻抗差的目标函数，应用共轭梯度法求解目标函数，得到波阻抗差函数最优解；

(64)利用求解的最佳波阻抗差函数，根据权利要求(5)中的总散射场正演计算公式，得到总的散射波场；

(65)从波场S中减去总的散射波场。

7.如权利要求1所述的在地震数据处理中衰减散射噪声的方法，其中在(g)步中：使用α-trimmed滤波技术去除(f)步输出波场中的极低频噪声。

Images