一种共反射点道集优化处理方法及系统与流程

本发明属于地震勘探资料处理中的地震反演成像处理领域，具体涉及一种共反射点道集优化处理方法及系统。

背景技术：

叠前偏移得到的共反射点(CRP)道集由反映地下同一点的反射信息组成，理论上基本消除了界面弯曲和倾斜地层的影响，对一系列共反射点道集的全叠加即得到偏移剖面。对于精确构造成像，由于受地震资料信噪比、道集不平、子波动校拉伸等因素的影响，共反射点道集的叠加并不能得到理想的成像结果；对于基于CRP道集的叠前反演，对道集的处理提出了更高的要求，除了要提高地震数据的信噪比和压制噪声的同时必须保持反射波振幅的相对强弱关系，如何获得高信噪比且同相轴水平的叠前资料直接影响着反演结果的稳定性和准确性。因此，无论是对于构造成像还是叠前反演，对共反射点道集进行一定的优化处理，显得十分必要。

针对共反射点道集中存在的信噪比问题，传统的f-k域倾向滤波、f-x域预测滤波、中值滤波等滤波方法由于没有方向特性，易造成对倾斜同相轴振幅的压制，同时也会模糊小的断层和裂缝等不连续结构；针对共反射点道集中存在的道集拉平问题，传统做法是借助道集叠加(全叠加或近-中道部分叠加)来建立标准模型道，然后以模型道作为标准对处理道集进行剩余时差校正，但是当道集中剩余校正量太大时，会在共反射点道集中引入新的误差。

技术实现要素：

本发明所要解决的技术问题是，针对现有技术中存在的难题，提供一种能有效地解决传统方法中存在的不足，在提高道集整体信噪比的同时使得剩余时差被校正，提高道集整体质量的共反射点道集优化处理方法及系统。

本发明解决上述技术问题的技术方案如下：一种共反射点道集优化处理方法，具体包括以下步骤：

a.获取共反射点道集，开设滑动时窗，在滑动窗口内，对共反射点道集中的地震数据采用常规的二次多项式拟合获得一组样点数据，根据样点数据走向进行估计地层倾向；

b.根据所述二次多项式所在方向确定滑动时窗内的同相轴方向；在滑动时窗内以样点数据中的中心样点的同相轴所在方向运用中值滤波，获得去噪共反射点道集；

c.对去噪共反射点道集进行叠加，得到叠加剖面，在叠加剖面上求取模型道；

d.对模型道和处理道求校正量函数，得到共反射点道集的剩余校正量模型，并对剩余校正量模型做平滑处理；

e.基于平滑处理后的剩余校正量模型对去噪共反射点道集进行校正处理，得到优化处理后的结果道集。

本发明的有益效果是：理论模型及实际地震数据试验结果表明：本发明方法对地震数据进行拟合来估计地层倾向并结合结构中值滤波来提高CRP道集信噪比以及对道集进行剩余时差校正的优化处理，有效提高了道集信噪比，道集被拉平，提高了道集整体质量。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述步骤a中所述的二次多项式为：

t＝t₀+c₁x+c₂x² 公式(1)，

其中，t为地震数据S(t)对应的离散时间采样点，x为拟合中心道的相对 道号，x有正负方向，其取值为x＝-M,-M+1,...,-1,0,1...,M-1,M，M为被拟合道集的道数的一半，t₀为拟合中心点的时间，C1和C2为拟合系数(其取值根据实际数据中同相轴倾角范围特征给定，一般在(-6,+6)范围内取值)。

进一步，所述公式(1)中采用最大互相关法则求取最佳拟合系数，即确定公式(1)中C1和C2的取值，在给定系数初始范围内，对不同系数求拟合道互相关，当互相关值最大时的系数即为最佳拟合系数。

进一步，所述互相关目标函数为：

$R\left(i\right)=\frac{1}{N}\frac{\underset{t=-l/2}{\overset{l/2}{\Σ}}{(\underset{k=-M}{\overset{M}{\Σ}}{S}_k(t_0+c_1(i)\×k+t))}^2}{\underset{t=-l/2}{\overset{l/2}{\Σ}}\underset{j=-M}{\overset{M}{\Σ}}{S}_j^2(t_0+c_1\left(i\right)\×j+t)}$ 公式(2)，

其中：R(i)为时窗内道集互相关值，S_k，S_j为拟合道样点值，l为相关时窗内样点数。相关系数R(i)在0-1间变化，理想情况下，当时窗内为同相轴方向时各道相关性最强，此时R(i)趋于1，R(i)越小则相关性越差。

进一步，求取最佳拟合系数的过程包括：首先假定C2的取值为零，通过扫描C1系数，当相关值R(i)最大时的系数即为最佳系数C1的取值；再将最佳系数C1的取值代入公式(1)，通过扫描C2系数，当相关值R(i)最大时的系数即为最佳系数C2的取值。

进一步，所述步骤b中滑动时窗依次沿时间T轴和X轴方形滑动，获得去噪共反射点道集。

进一步，所述步骤c中在叠加剖面上求取模型道的过程包括：

在叠加剖面上对计算样点所在时窗采用常规的二次多项式拟合估计地层倾向；

在时窗内取倾向方向中值作为该计算样点的值，在叠加剖面上完成对叠加道的处理获得该叠加道对应的共反射点道集的模型道；

进一步，所述步骤d中采用常规最优化方法计算校正量函数；

对剩余校正量模型沿时间T轴方向和X轴方向做平滑处理。

进一步，所述步骤e中校正处理采用以下公式(3)处理：

其中，G(t,x)为原始道集或去噪后道集，优化处理后的结果道集，Δm(t,x)为步骤d获取的校正量函数，x为拟合中心道的相对道号，x有正负方向，其取值为x＝-M,-M+1,...,-1,0,1...,M-1,M。

校正量函数Δm(t,x)的计算方法是将去噪后共反射点道集中的每一道和模型道做相关，当相关值最大时的时移量道集即为当前共反射点道集对应的校正量函数的取值，对一个共反射点道集中所有道重复此处理，得到一个与共反射点道集具有相同道数的道集，此即为校正量函数Δm(t,x)。

进一步，在求取校正量函数过程中，设定的最大允许时移量应随着偏移距的增大而增大；并给定互相关阀值，当互相关值小于阀值则做舍弃处理，以排除由于干扰造成的误差。

技术思路及技术实现：

针对共反射点道集中存在的信噪比和道集拉平问题，本发明提出基于结构中值滤波对CRP道集进行优化处理。首先采用离散扫描法对地层倾向进行估计，为了使得倾向的估计更加准确，本发明采用常规二次多项式对地震数据进行拟合来估计地层倾向，然后在滑动时窗内以中心样点同相轴所在方向运用中值滤波，滤波结果作为最终中心样点的输出即结构中值滤波，处理后提高了道集信噪比，最后在共反射点道集叠加剖面上建立模型道，对去噪后共反射点道集进行剩余时差校正。由于地层倾向方向具有更强的信息相关性，因此本方法能够有效地解决传统方法中存在的不足，在提高道集整体信噪比的同时使得剩余时差被校正，表现为道集被拉平，提高了道集整体质量，并能够较好地保持有效信号和断层、裂缝等地层边缘和细节特征。

本发明解决上述技术问题的技术方案如下：一种共反射点道集优化处理 系统，包括倾向估计模块、滤波模块、叠加模块、平滑处理模块和校正模块；

所述倾向估计模块用于获取共反射点道集，对共反射点道集中的地震数据采用常规的二次多项式拟合获得一组样点数据，根据样点数据走向进行估计地层倾向；

所述滤波模块用于开设滑动时窗，根据所述二次多项式所在方向确定滑动时窗内的同相轴方向；在滑动时窗内以样点数据中的中心样点的同相轴所在方向运用中值滤波，获得去噪共反射点道集；

所述叠加模块用于对去噪共反射点道集进行叠加，得到叠加剖面，在叠加剖面上求取模型道；

所述平滑处理模块用于对模型道和处理道求校正量函数，得到共反射点道集的剩余校正量模型，并对剩余校正量模型做平滑处理；

所述校正模块用于基于平滑处理后的剩余校正量模型对去噪共反射点道集进行校正处理，得到优化处理后的结果道集。

本发明的有益效果是：理论模型及实际地震数据试验结果表明：本发明方法对地震数据进行拟合来估计地层倾向并结合结构中值滤波来提高CRP道集信噪比以及对道集进行剩余时差校正的优化处理，有效提高了道集信噪比，道集被拉平，提高了道集整体质量。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述倾向估计模块中所述的二次多项式为：

t＝t₀+c₁x+c₂x² 公式(1)，

其中，t为地震数据S(t)对应的离散时间采样点，x为拟合中心道的相对道号，x有正负方向，其取值为x＝-M,-M+1,...,-1,0,1...,M-1,M，M为被拟合道集的道数的一半，t₀为拟合中心点的时间，C1和C2为拟合系数(其取值根据实际数据中同相轴倾角范围特征给定，一般在(-6,+6)范围内取值)。

进一步，所述公式(1)中采用最大互相关法则求取最佳拟合系数，即 确定公式(1)中C1和C2的取值，在给定系数初始范围内，对不同系数求拟合道互相关，当互相关值最大时的系数即为最佳拟合系数。

进一步，所述互相关目标函数为：

$R\left(i\right)=\frac{1}{N}\frac{\underset{t=-l/2}{\overset{l/2}{\Σ}}{\left(\underset{k=-M}{\overset{M}{\Σ}}{S}_k\right(t_0+c_1\left(i\right)\×k+t\left)\right)}^2}{\underset{t=-l/2}{\overset{l/2}{\Σ}}\underset{j=-M}{\overset{M}{\Σ}}{S}_j^2(t_0+c_1\left(i\right)\×j+t)}$ 公式(2)，

其中：R(i)为时窗内道集互相关值，S_k，S_j为拟合道样点值，l为相关时窗内样点数。相关系数R(i)在0-1间变化，理想情况下，当时窗内为同相轴方向时各道相关性最强，此时R(i)趋于1，R(i)越小则相关性越差。

进一步，求取最佳拟合系数的过程包括：首先假定C2的取值为零，通过扫描C1系数，当相关值R(i)最大时的系数即为最佳系数C1的取值；再将最佳系数C1的取值代入公式(1)，通过扫描C2系数，当相关值R(i)最大时的系数即为最佳系数C2的取值。

进一步，所述滤波模块中滑动时窗依次沿时间T轴和X轴方形滑动，获得去噪共反射点道集。

进一步，所述叠加模块中在叠加剖面上求取模型道的过程包括：

在叠加剖面上对计算样点所在时窗采用常规的二次多项式拟合估计地层倾向；

在时窗内取倾向方向中值作为该计算样点的值，在叠加剖面上完成对叠加道的处理获得该叠加道对应的共反射点道集的模型道。

进一步，所述平滑处理模块中采用常规最优化方法计算校正量函数；

对剩余校正量模型沿时间T轴方向和X轴方向做平滑处理。

进一步，所述校正模块中校正处理采用以下公式(3)处理：

其中，G(t,x)为原始道集或去噪后道集，优化处理后的结果道集， Δm(t,x)为步骤d获取的校正量函数，x为拟合中心道的相对道号，x有正负方向，其取值为x＝-M,-M+1,...,-1,0,1...,M-1,M。

进一步，在求取校正量函数过程中，设定的最大允许时移量应随着偏移距的增大而增大；并给定互相关阀值，当互相关值小于阀值则做舍弃处理，以排除由于干扰造成的误差。

附图说明

图1为本发明所述的一种共反射点道集优化处理方法流程图；

图2为本发明所述的一种共反射点道集优化处理系统结构框图；

图3为本发明离散扫描法估计地层倾向示意图；

图4a为选取主频为30Hz的含有噪声及剩余时差的模型数据；

图4b为结构中值滤波后的共反射点道集；

图5a为共反射点道集剩余校正量模型；

图5b为对图4b的共反射点道集应用剩余校正量模型校正后的道集；

图6a为实施例2的一个共反射点道集；

图6b为对图6a所示的共反射点道集进行滤波后的共反射点道集；

图7a为对图6b所示的共反射点道集进行优化后得到的剩余校正量模型；

图7b为实施例2的采用本发明所述方法最终处理得到的结果道集；

图8a为实施例2的未经本发明所述方法处理的剖面图；

图8b为实施例2的采用本发明所述方法处理后的剖面图；

图9a为实施例2未经本发明所述方法处理的剖面图(波形变面积显示)；

图9b为实施例2采用本发明所述方法处理后的剖面图(波形变面积显示)。

附图中，各标号所代表的部件列表如下：

1、倾向估计模块，2、滤波模块，3、叠加模块，4、平滑处理模块，5、校正模块。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

如图1所示，为本发明所述的一种共反射点道集优化处理方法，具体包括以下步骤：

a.获取共反射点道集，开设滑动时窗，在滑动窗口内，对共反射点道集中的地震数据采用常规的二次多项式拟合获得一组样点数据，根据样点数据走向进行估计地层倾向；

b.根据所述二次多项式所在方向确定滑动时窗内的同相轴方向；在滑动时窗内以样点数据中的中心样点的同相轴所在方向运用中值滤波，获得去噪共反射点道集；

c.对去噪共反射点道集进行叠加，得到叠加剖面，在叠加剖面上求取模型道；

d.对模型道和处理道求校正量函数，得到共反射点道集的剩余校正量模型，并对剩余校正量模型做平滑处理；

e.基于平滑处理后的剩余校正量模型对去噪共反射点道集进行校正处理，得到优化处理后的结果道集。

如图2所示，为本发明所述的一种共反射点道集优化处理系统，包括倾向估计模块1、滤波模块2、叠加模块3、平滑处理模块4和校正模块5；

所述倾向估计模块1用于获取共反射点道集，对共反射点道集中的地震数据采用常规的二次多项式拟合获得一组样点数据，根据样点数据走向进行估计地层倾向；

所述滤波模块2用于开设滑动时窗，根据所述二次多项式所在方向确定滑动时窗内的同相轴方向；在滑动时窗内以样点数据中的中心样点的同相轴所在方向运用中值滤波，获得去噪共反射点道集；

所述叠加模块3用于对去噪共反射点道集进行叠加，得到叠加剖面，在叠加剖面上求取模型道；

所述平滑处理模块4用于对模型道和处理道求校正量函数，得到共反射点道集的剩余校正量模型，并对剩余校正量模型做平滑处理；

所述校正模块5用于基于平滑处理后的剩余校正量模型对去噪共反射点道集进行校正处理，得到优化处理后的结果道集。

针对共反射点道集中存在的信噪比和道集拉平问题，本发明提出利用二次多项式对地震数据进行拟合来估计地层倾向并结合结构中值滤波来提高CRP道集信噪比以及对道集进行剩余时差校正的优化处理。由于地层倾向方向具有更强的信息相关性，因此，本发明能够有效地解决传统方法中存在的不足，在提高道集整体信噪比的同时使得剩余时差被校正，表现为道集被拉平，提高了道集整体质量，并能够较好地保持有效信号和断层、裂缝等地层边缘和细节特征。

本发明针对共反射点道集中存在的信噪比和道集拉平问题提供了一种有效的优化处理方法,由地层倾向估计、结构中值滤波、道集拉平三部分组成。

地层倾向估计：

采用离散扫描法对地层倾向进行估计是一种简单且计算效率高的方法，该方法基于同相轴在局部满足线性的假设，对离散点在计算时窗内求相关。为了使得倾向的估计更加准确，本发明采用常规二次多项式对地震数据进行拟合来估计地层倾向，此方法是对正交多项式拟合方法的改进，克服了正交多项式拟合方法的不足。为了减少计算量，可以先对一次项系数进行扫描，然后进行二次项系数扫描，对于道集中的每一个样点，在计算时窗内通过给 定不同的系数(倾向方向)计算其相关值，当相关值最大时的系数即为当前样点的最佳一次项拟合系数。如图1所示，拟合道集的中心道为参考道，上、下两线间为一次项拟合系数扫描范围，当扫描方向为计算样点所在同相轴方向时(中间线所示)相关值为最优值，一次项系数确定后，再对二次项系数进行扫描来进一步提高拟合精度。

结构中值滤波：

与常规中值滤波方法的不同之处在于其滤波具有方向性。基于常规二次多项式对同相轴追踪后，在滑动时窗内以中心样点同相轴所在方向运用中值滤波，滤波结果作为最终中心样点的输出。结构中值滤波相对于常规随机噪声压制方法具有以下优点：(1)在一定倾向方向上进行滤波处理；(2)无论同相轴是否水平，都能对随机噪声和异常振幅进行有效压制；(3)滤波的同时保护断层边界，能够保护断层的不连续性。因此，结构中值滤波不仅能够有效地衰减随机噪声、提高剖面的信噪比，而且能够较好地保持有效信号和断层、裂缝等地层边缘和细节特征。

道集拉平：

静校正法道集拉平中，借助道集叠加(全叠加或近-中道部分叠加)来建立标准模型道，然后以模型道作为标准对处理道集进行剩余时差校正，因此标准模型道的求取是校正效果好坏的关键。借助道集叠加来建立标准模型道时要求道集中剩余校正量不能太大，否则会在共反射点道集中引入新误差。考虑到相关性在叠加剖面中要优于存在剩余时差的道集间，因此本方法在共反射点道集叠加剖面上建立模型道，能充分利用叠加剖面上同相轴相关性优于单道记录特性，增强稳定性，避免了叠加道中引入新的误差。

本发明具体实施例1具体实现步骤为：

步骤一：采用常规二次多项式对地震数据(互相关目标函数中的S_k，S_j即为实际离散的地震数据，拟合关系式t＝t₀+c₁x+c₂x²求得的时间t为地震数 据S(t)对应的离散时间采样点)进行拟合来估计地层倾向，如图3所示，为本发明离散扫描法估计地层倾向示意图，在计算滑动时窗内，局部同相轴满足以下二次拟合多项式：

t＝t₀+c₁x+c₂x² 公式(1)

其中：x为拟合中心道的相对道号，x有正负方向，其取值为x＝-M,-M+1,...,-1,0,1...,M-1,M，M为拟合道数的一半，t₀为拟合中心点的时间。一次项最佳拟合系数的确定采用最大互相关法则，在给定系数初始范围(C1和C2为系数初始扫描范围，根据实际数据中同相轴倾角范围特征给定，一般在(-6,+6)范围内取值)内，对不同系数(若给定C1的范围为(-3，+3)，步长为dc，那么这些系数即为：-3，-3+dc，-3+2dc，…，3-2dc，3-dc，3)求拟合道互相关，当互相关值最大时的系数即为最佳拟合系数(是指C1和C2的值，本发明采用分步求系数C1和C2，第一步求C1时，假定C2始终为零，通过扫描C1系数，当相关值R(i)最大时的系数即为最佳系数C1；第二步求C2时，将第一步中求出的C1作为常量代入t＝t₀+c₁x+c₂x²中，以同样的扫描方法求取最佳系数C2)，互相关目标函数为：

$R\left(i\right)=\frac{1}{N}\frac{\underset{t=-l/2}{\overset{l/2}{\Σ}}{\left(\underset{k=-M}{\overset{M}{\Σ}}{S}_k\right(t_0+c_1\left(i\right)\×k+t\left)\right)}^2}{\underset{t=-l/2}{\overset{l/2}{\Σ}}\underset{j=-M}{\overset{M}{\Σ}}{S}_j^2(t_0+c_1\left(i\right)\×j+t)}$ 公式(2)

其中：R(i)为时窗内道集互相关值，S_k，S_j为拟合道样点值，l为相关时窗内样点数。相关系数R(i)在0～1间变化，理想情况下，当时窗内为同相轴方向时各道相关性最强，此时R(i)趋于1，R(i)越小则相关性越差；

步骤二：获得同相轴方向(步骤一中确定最佳系数C1、C2后，则二次多项式t＝t₀+c₁x+c₂x²所在方向确定下来，此方向即为当前滑动时窗内的同相轴方向)后，在滑动时窗内以中心样点同相轴所在方向运用中值滤波，滤波结果作为最终中心样点的输出。中值a_m的求取表达式为：

$a_m=\underset{a_m\∈S_i}{\mathrm{argmin}}\underset{i=1}{\overset{N}{\Σ}}\left|\right|a_m-a_i|{|}_L$

其中N为序列个数，L为范数，a_m为需要求取的中值，a_i为参与求取中值的样点序列即同相轴方向上的离散地震数据。滑动时窗依次沿着时间T轴和X轴方向滑动即完成一个共反射点道集的处理，处理后得到的是经过结构中值滤波即随机噪声得到有效压制，信噪比得到提高后的共反射点道集；

步骤三：对结构中值滤波后随机噪声得到有效压制，信噪比得到提高后，得到去噪后CRP道集，同时对道集进行叠加，得到叠加剖面；

步骤四：在叠加剖面上对计算样点所在时窗内进行地层倾向估计(方法和步骤一相同，只是前面是在道集上处理，这里是在叠加剖面上处理)，在时窗内取倾向方向中值作为该计算样点的值，在叠加剖面上，完成一个叠加道的处理(一个共反射点道集经过叠加后在叠加剖面上变成了一道，在叠加剖面上对这一道进行步骤1和步骤2的结构中值滤波后的输出同样也是一道，此道即为对应这个共反射点道集的模型道)即可得到该叠加道对应的共反射点道集的模型道G_m(t,x)，此步骤在叠加剖面上求取的模型道；另一方面去噪后道集也作为需要进行剩余时差校正的处理道输入即G(t,x))；

步骤五：对模型道和处理道以最优化方法求校正量函数(以常规最优化方法求校正量函数，即将处理道集中的每一道和模型道做相关，当相关值最大时，得到一个时移量道即为当前处理道对应的校正量函数道，对一个共反射点道集中所有道重复此处理，得到一个与共反射点道集具有相同道数的道集，即得到此共反射点道集的剩余校正量模型)，得到共反射点道集的剩余校正量模型并沿T方向与X方向做平滑处理(平滑后的剩余校正量模型即为Δm(t,x)，只是公式中是以关于空间和时间的函数给出，其值即对应平滑后的剩余校正量模型中的数值)，以达到准确求取时移量的目的；

步骤六：对CRP道集进行校正处理：

其中，G(t,x)为原始道集或步骤3中得到的去噪后道集，为校正后道集，Δm(t,x)为步骤5中得到的校正量函数。在求取校正量函数过程中，设定的最大允许时移量应随着偏移距的增大而增大；并给定互相关阀值，当互相关值小于阀值则做舍弃处理，以排除由于干扰造成的误差。

通过理论模型及实际地震数据试验，本方法对共反射点道集的优化处理效果明显，具有较强的针对性。

图4a为选取主频为30Hz的雷克子波合成含有随机噪声的共反射点道集，采样点数为401，采样率为4ms，共25道。合成模型中共四个同相轴，作为对比第三个同相轴为完全水平，其余三个同相轴含有不同剩余时差。时窗取5x5,初始系数扫描范围为(-1,1)，图4b为结构中值滤波后道集，对比图4a可知道集中随机噪声得到压制，信噪比明显提高。对于理论模型，以同相轴完全水平时道集的叠加作为模型道，对去噪后道集进行剩余量分析和校正。图5a为道集剩余校正量模型(平滑后)，由于第三个同相轴为完全水平，其剩余校正量均为零值，图5b是对图4b应用剩余校正量模型校正后的道集。

图6-图9为本发明的实施例2，图6a为以实际资料为一个共反射点道集(局部)，道集中存在噪声干扰和剩余时差，道集叠加不能完全聚焦。选取时窗为7x5，初始系数扫描范围为(-3,3)，图6b为结构中值滤波后结果，在叠加剖面中建立对应CRP道集的参考模型道并对滤波后道集(图6b)进行优化处理。图7a为得到的剩余校正量模型，图7b为校正后道集，原始道集(图6a)与最终处理道集(图7b)的对比可知，道集中噪声得到压制并且剩余时差被校正，表现为道集被拉平，提高了道集整体质量。图8为应用本方法处理前后叠加对比剖面，图8a为未经本发明所述方法处理的剖面图；图8b为采用本发明所述方法处理后的剖面图；处理后叠加剖面(图8b)信噪比得到改善并且同相轴连续性增强，整体剖面质量得到明显提高。叠加剖 面的波形变面积显示更清楚地反映了处理前后同相轴的变化，如图9所示，图9a为未经本发明所述方法处理的波形变面积显示剖面图；图9b为采用本发明所述方法处理后的波形变面积显示剖面图。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。