CN104656140B - 一种基于垂直地震衰减规律的中值滤波方法以及系统 - Google Patents

Abstract

本发明提供一种基于垂直地震衰减规律的中值滤波方法及系统，所述方法包括：采集垂直地震剖面VSP数据；根据所述的VSP数据分选出单分量数据；对所述的单分量数据进行零相位化处理，得到零相位化处理后的剖面；根据所述零相位化处理后的剖面得到时距曲线；根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子；将零相位化处理后的剖面沿所述时距曲线拉齐，得到拉齐后的剖面；根据所述VSP数据的TAR因子对所述拉齐后的剖面进行定向加权中值滤波。解决了用于井控的VSP数据处理中没有可靠波场分离方法的问题。

Description

一种基于垂直地震衰减规律的中值滤波方法以及系统

技术领域

本发明关于地球物理勘探技术领域，特别是关于地球物理勘探中垂直地震数据处理领技术，具体的讲是一种基于垂直地震衰减规律的中值滤波方法及系统。

背景技术

波场信息丰富是垂直地震剖面VSP数据的一个重要特点，在剖面中可以看到多种波场相互交织的现象。通常，VSP数据又不具备地面地震较高的覆盖次数，故而很难通过叠加的方法强化或衰减某一种波场。因此，波场分离方法在VSP数据的处理中的应用就显示尤为重要。

现有技术中，VSP数据的处理方案中并没有针对性的波场分离方法，其波场分离手段都是借用地面地震数据处理的方法，如频率—波数滤波(FK滤波)、均值滤波和中值滤波等。其中，FK滤波容易产生假频，在覆盖次数低、质量要求高的VSP数据波场中的应用十分有限；均值滤波和中值滤波方法都属于时间域的多道速度滤波，中值滤波的抗噪能力优于均值滤波，在实际生产中应用较多，并且派生出自适应去异常、边界变道数等中值滤波方法。在VSP数据处理中，虽然这些波场分离方法并没有考虑到VSP数据自身的特点，但在传统的处理流程中，波场分离步骤是在振幅补偿之后进行的，此时剖面振幅能量相对统一，是可以利用中值方法进行波场分离的。

近年来，VSP数据在井控技术方面得到了快速发展，其振幅衰减参数、球面扩散补偿参数(TAR)、地层吸收补偿参数(Q)等在地面地震处理中大量应用，起到了很好的井控约束作用。然而，这些井控参数的提取都需要在VSP数据振幅相对保真前提下完成，要求在振幅随VSP数据记录深度呈指数变化的剖面中进行波场分离。通常，一个同相轴的不同位置存在着非常大的振幅能量差异，这将极大的影响中值滤波方法的可靠性。

因此，如何提出一种新的滤波方法，进而解决用于井控的VSP数据数据处理中没有可靠波场分离是本领域亟待解决的技术难题。

发明内容

为了解决现有技术中用于井控的VSP数据数据处理中没有可靠波场分离方法的问题，本发明在自适应去异常的、边界变道数的中值滤波基础上，提出一种基于垂直地震衰减规律的中值滤波方法及系统，是一种精确的基于垂直地震衰减规律的中值滤波的方案，由VSP数据获得单分量数据，对单分量数据进行零相位化处理，结合时距曲线确定球面扩散补偿(TAR)因子，并沿时距曲线拉齐剖面，进行定向加权中值滤波，极大提高了中值滤波的可靠性。

本发明的目的之一是，提供一种基于垂直地震衰减规律的中值滤波方法，包括：采集垂直地震剖面VSP数据；根据所述的VSP数据分选出单分量数据；对所述的单分量数据进行零相位化处理，得到零相位化处理后的剖面；根据所述零相位化处理后的剖面得到时距曲线；根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子；将零相位化处理后的剖面沿所述时距曲线拉齐，得到拉齐后的剖面；根据所述VSP数据的TAR因子对所述拉齐后的剖面进行定向加权中值滤波。

本发明的目的之一是，提供了一种基于垂直地震衰减规律的中值滤波的系统，包括：剖面采集装置，用于采集垂直地震剖面VSP数据；单分量数据分选装置，用于根据所述的VSP数据分选出单分量数据；零相位化处理装置，用于对所述的单分量数据进行零相位化处理，得到零相位化处理后的剖面；时距曲线确定装置，用于根据所述零相位化处理后的剖面得到时距曲线；TAR因子确定装置，用于根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子；拉齐装置，用于将零相位化处理后的剖面沿所述时距曲线拉齐，得到拉齐后的剖面；中值滤波装置，用于根据所述VSP数据的TAR因子对所述的拉齐后的剖面进行定向加权中值滤波。

本发明的有益效果在于，在自适应去异常的、边界变道数的中值滤波基础上，提供了一种基于垂直地震衰减规律的中值滤波方法及系统，是一种精确的基于垂直地震衰减规律的中值滤波的方案，由VSP数据分选出单分量数据，对单分量数据进行零相位化处理，结合时距曲线确定TAR因子，并沿时距曲线拉齐剖面，进行定向加权中值滤波，极大提高了中值滤波的可靠性，为用于井控的VSP数据数据处理提供了可靠的波场分离方案。

为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波方法的流程图；

图2为图1中的步骤S104的具体流程图；

图3为图1中的步骤S105的具体流程图；

图4为图1中的步骤S107的具体流程图；

图5为图4中的步骤S404的具体流程图；

图6为图4中的步骤S406的具体流程图；

图7为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波系统的结构框图；

图8为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波系统中的时距曲线确定装置的具体结构框图；

图9为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波系统中的TAR因子确定装置的具体结构框图；

图10为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波系统中的中值滤波装置的具体结构框图；

图11为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波系统中的加权数据序列确定模块的具体结构框图；

图12为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波系统中的去异常值处理模块的具体结构框图；

图13为本发明提供的具体实施例中原始VSP数据Z分量示意图；

图14为本发明提供的具体实施例中地震子波零相位化示意图；

图15为本发明提供的具体实施例中VSP数据初至拾取示意图；

图16为本发明提供的具体实施例中利用VSP数据计算的TAR因子的计算结果图；

图17为本发明提供的具体实施例中TAR的计算差异曲线图；

图18为本发明提供的具体实施例中沿下行波时距曲线拉平后的VSP数据示意图；

图19为本发明提供的具体实施例中一个窗内的样点序列示意图；

图20为本发明提供的具体实施例中依VSP数据规律的加权系数序列示意图；

图21为本发明提供的具体实施例中加权之后一个窗内的样点序列示意图；

图22为本发明提供的具体实施例中自适应去异常后的样点序列示意图；

图23为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图1；

图24为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图2。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明涉及地球物理勘探中垂直地震数据处理领域，属于运用地震资料处理手段进行波场分离的范畴，是一种在自适应去异常和边界变道数中值滤波基础上的、应用垂直地震数据衰减规律进行加权的中值滤波方法。

为了解决用于井控的VSP数据数据处理中没有可靠波场分离方法的问题，本发明在自适应去异常的、边界变道数的中值滤波基础上，提出一种利用VSP数据衰减规律进行加权的全新中值滤波方法。图1为本发明实施例提供的一种基于垂直地震衰减规律的中值滤波方法的流程图，由图1可知，所述的方法包括：

S101：采集垂直地震剖面VSP数据。

S102：根据所述的VSP数据分选出单分量数据。

在具体的实施例中，利用现有方法和地震处理软件对VSP原始采集数据进行整理，获得准备进行波场分离的单分量数据，如图13所示，为本发明提供的具体实施例中原始VSP数据Z分量示意图。

所述的数据整理不包括振幅补偿处理、道均衡处理、自动增益(AGC)处理或其它可能改变振幅或能量相对关系的处理步骤。

S103：对所述的单分量数据进行零相位化处理，得到零相位化处理后的剖面。

在具体的实施例中，利用现有软件进行零相位化处理。零相位化处理的目的是将混合相位的复杂子波转化为零相位的子波，使子波的波峰或波谷与地层界面相对应，如图14所示，为本发明提供的具体实施例中地震子波零相位化示意图。

S104：根据所述零相位化处理后的剖面得到时距曲线。

图2为步骤S104的具体流程图，由图2可知，根据所述零相位化处理后的剖面得到时距曲线包括：

S201：将所述零相位化处理后的剖面中下行直达波的波峰位置作为拾取位置；

S202：以所述拾取位置拾取所述下行直达波，得到时距曲线。

在具体的实施方式中，利用现有地震处理软件拾取下行直达波初至，得到时间-深度关系序列，该序列在平面上构成一条时间与距离的关系曲线，以下简称为时距曲线，如图15所示，为本发明提供的具体实施例中VSP数据初至拾取示意图。

所述初至拾取位置为剖面中下行直达波的波峰位置。

所述下行波时距曲线同时作为后续球面扩散补偿(TAR)因子的时距曲线。

由图1可知，该方法还包括：

S105：根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子。图3为步骤S105的具体流程图，由图3可知，该步骤具体包括：

S301：在所述零相位化处理后的剖面中，拾取所述时距曲线的振幅值，得到时间-振幅关系序列；

S302：对所述时间-振幅关系序列中的时间与振幅分别求取自然对数，得到新序列；

S303：对所述的新序列进行线性拟合，得到拟合直线；

S304：确定所述拟合直线的斜率；

S305：对所述斜率取相反数，得到所述VSP数据的球面扩散补偿TAR因子。

在具体的实施方式中，在步骤S103所得到的剖面中，拾取步骤S104时距曲线位置的振幅值，获得时间-振幅关系序列，对该序列中时间与振幅分别求取自然对数，然后对所得到的新序列进行线性拟合，计算得到拟合直线的斜率，并对该斜率值取相反数，即得到垂直地震的TAR因子，如图17所示，为本发明提供的具体实施例中利用VSP数据计算的TAR因子的计算结果图。

由图1可知，该方法还包括：

S106：将零相位化处理后的剖面沿所述时距曲线拉齐，得到拉齐后的剖面。

在具体的实施方式中，利用现有软件将步骤S103所得到的数据沿步骤S104所得到的时距曲线拉齐，获得准备分离下行波场的拉齐后的剖面，如图18所示，为本发明提供的具体实施例中沿下行波时距曲线拉平后的VSP数据示意图。

S107：根据所述VSP数据的球面扩散补偿TAR因子对所述拉齐后的剖面进行定向加权中值滤波。图4为步骤S107的具体流程图，由图4可知，该步骤具体包括：

S401：将所述拉齐后的剖面中地震道从最浅一级检波器记录向最深一级检波器记录排序；

S402：选定所述拉齐后的剖面中每次参与运算的地震道数，称为中值滤波窗。所述中值滤波窗在靠近拉齐后的剖面的边界时逐步减小，以保证中值滤波窗内不出现空样点值，使得后续处理结果相对保真。

S403：将所述中值滤波窗内的振幅值提取出来，组成离散数据序列；

在具体的实施方式中，确定中值滤波窗(跨度)，即每次处理运算的最大样点数，在VSP数据剖面上也就是选定参与每一次统计计算的地震道数。

这些地震道上的同时刻样点值组成离散的数据序列。所述离散数据序列，要求从最浅一级检波器对应地震道向最深一级检波器对应的地震道方向排序。

在具体的实施方式中，确定用于中值滤波处理的最大样点数设为n，即选定参与每一次统计计算的地震道数，这些地震道上的同时刻样点值组成离散的数据序列X，如(1)式所示。

X＝{a₁,a₂,a₃,…a_n} (1)

式中：X为离散的振幅序列；a表示单个样点的振幅；n为最大样点数。

由时深关系可知，这n个样点对应的时间序列T如下式：

T＝{t₁,t₂,t₃,…t_n} (2)

式中：T为离散的时间序列；t表示单个样点的时间，单位为秒(s)为；n最大样点数。

在具体的实施方式中，所述用于统计计算的地震道数取值范围为介于7和15之间的奇数，该参数的确定要综合分析波场分离效果和数据保真两个方面，具体数值视信剖面信噪比情况而确定，在信噪比相对比较低的地震数据中选用7作为样点数，在信噪比较高的情况下选用15作为样点数，信噪比介于两者之间时选用奇数9、11或13。该取值在靠近拉齐后的剖面的边界时逐步减小，以保证窗口内不出现空样点值，使得后续处理结果相对保真。

S404：根据所述VSP数据的球面扩散补偿TAR因子以及所述中值滤波窗得到加权数据序列。在具体的实施方式中，利用步骤S105中所计算得到的球面扩散补偿TAR因子，以及步骤S402中所确定的滑动窗长度，设计加权函数，得到递增的加权数据序列。

图5为步骤S404的具体流程图。由图5可知，该步骤具体包括：

S501：根据所述VSP数据的TAR因子确定任意一个中值滤波窗内的TAR因子。利用步骤S105中所计算得到的TAR因子k，可以得到所述地震数据的剖面中任意一个时刻t的球面扩散补偿，如下式所示：

A_t＝A₀t^k (3)

式中：A_t为任意一个时刻的振幅；A₀为零时刻的振幅；t表示样点所在的时刻，单位为秒(s)；k为TAR因子。

S502：根据地震数据中任意一个中值滤波窗内的TAR因子对所述离散数据序列对应的中值滤波窗内的样点值进行相对振幅关系补偿。

在具体的实施方式中，由于本发明致力于保护原始的振幅信息，因此不能基于公式(3)进行整个剖面的补偿，但可以利用该式对步骤S502窗中的样点进行相对振幅关系的补偿。具体方法是选择窗的第个样点为基准点进行补偿处理，如(4)式所示。所选样点即为原始数据上准备替换的样点。该步骤是本发明的关键步骤，可以使原始震幅得到最好的保护。

对(4)式进行整理得到：

S503：对相对振幅关系补偿后的样点值进行基准点归一化处理。

在具体的实施方式中，以第个样点振幅值为进行基准点归一化处理，得到：

在(4)至(6)式中：A_i是在窗内以第个样点为基准点的理论补偿振幅；A′_i为A_i相对于基准点归一化后的值；t表示样点所在的时刻，单位为秒(s)；k为TAR因子；n为窗内最大样点数。

S504：根据基准点归一化处理得到的样点值确定加权数据序列。

递增的加权数据序列Y，可以用下式表示。

式中，Y为振幅加权序列；t表示样点所在的时刻；k为TAR因子；n为最大样点数。

由图4可知，该步骤还包括：

S405：对所述的离散数据序列以及加权数据序列进行定向加权，得到定向加权后的数据序列。

在具体的实施方式中，对序列X和序列Y进行样点对应相乘，即完成定向加权计算，得到新的数据序列Z，如式(8)所示。

以A′表示(8)式中各项，记作：

Z＝{A′₁,A′₂,A′₃,…A′_n} (9)

在(8)式和(9)式中：Z为加权后新的振幅序列；X为原始的振幅序列；Y为振幅加权序列；a表示原始序列中单个样点的振幅；t表示样点所在的时刻；A′表示加权后单个样点的振幅；n为最大样点数。

S406：对定向加权后的数据序列进行去异常值处理，得到去异常值处理后的数据序列。图6为步骤S406的具体流程图，由图6可知，该步骤具体包括：

S601：确定定向加权后的数据序列中各个样点的绝对值的平均值；

S602：根据所述的平均值确定去异常门槛值。在具体的实施方式中，取数列中各个样点平均值的绝对值，并将该值加减20％作为一对去异常门槛值。

S603：从定向加权后的数据序列的各个样点中筛选出绝对值在所述去异常门槛值以外的样点，作为异常样点。在具体的实施方式中，将原数列中各个样点的绝对值与去异常门槛值进行逐一比较，筛选出绝对值在所述去异常门槛值以外的样点，作为异常样点。

S604：以定向加权后的数据序列中各个样点的绝对值的平均值替换所述异常样点，得到异常值处理后的数据序列。

也即本发明的去异常值处理，具体方法是取数列中各个样点平均值的绝对值，并将该值加减20％作为一对去异常的门槛值，再将原数列中各个样点的绝对值与门槛值进行逐一比较，并以平均值替换门槛值范围以外的样点值。如式(10)所示。

此时，自适应去异常后的振幅序列Z′为：

Z′＝{A″₁,A″₂,A″₃,…A″_n} (11)

在(10)式和(11)式中：Z′为自适应去异常后的振幅序列；A′表示加权后单个样点的振幅；A″表示自适应去异常后单个样点的振幅；n为最大样点数。

图19至图22中示意了一个样点窗的加权、自适应去异常的过程。其中，图19为本发明提供的具体实施例中一个窗内的样点序列示意图，图20为本发明提供的具体实施例中依VSP数据规律的加权系数序列示意图，图21为本发明提供的具体实施例中加权之后一个窗内的样点序列示意图，图22为本发明提供的具体实施例中自适应去异常后的样点序列示意图。

所述门槛值设定为均值的20％，该值远低于其它地震数据自适应去异常方法中的门槛值。这是由于VSP数据检波器深入井下，接收环境相对安静，因而出现异常高值或低值的地震野值情况较少；并且经过样点加权处理后，正常值范围内振幅变得更加均衡。

在本发明的其他实施例中，当窗滑动至边界位置时，自动减小窗长，以满足计算需要，回避边界的异常处理效应，该技术为公知技术，已应用于很多商业地震处理软件中。

如上即完成了准备工作，剖面中每一点便对应一个窗。

由图4可知，该步骤还包括：

S407：对去异常值处理后的数据序列进行中值滤波处理，得到定向加权中值滤波处理数据。

在本发明的具体实施例中，该步骤的具体算法是将窗中所有样点按照从大到小的规则重新排序，并取得新序列的中间一个样点的数值，用以代替剖面中的该点，从而完成自适应去异常的、边界变道数的、依照VSP数据规律的定向加权中值滤波处理，如图23为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图1，图24为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图2。

循环步骤S401至步骤S407的处理流程，可完成多组下行波中值滤波处理。对上行波进行中值滤波处理时，拾取上行波时距曲线，同时使用等比递减的加权序列，该变化是由VSP数据方法中地震波传播特点所决定的，其上行波振幅随接收检波器深度的变浅而减小，但衰减规律不变，不需要重新计算TAR因子。

如上即是本发明在自适应去异常的、边界变道数的中值滤波基础上，提出一种利用VSP数据衰减规律进行加权的全新中值滤波方法，解决用于井控的VSP数据数据处理中没有可靠波场分离方法的问题。

图7为本发明实施例提供的一种基于垂直地震衰减的中值滤波基于垂直地震衰减规律的中值滤波的系统的结构框图，由图7可知，所述的系统包括：

剖面采集装置101，用于采集垂直地震剖面VSP数据。

单分量数据分选装置102，用于根据所述的VSP数据分选出单分量数据。

在具体的实施例中，利用现有方法和地震处理软件对VSP原始采集数据进行整理，获得准备进行波场分离的单分量数据，如图13所示，为本发明提供的具体实施例中原始VSP数据Z分量示意图。

所述的数据整理不包括振幅补偿处理、道均衡处理、自动增益(AGC)处理或其它可能改变振幅或能量相对关系的处理步骤。

零相位化处理装置103，用于对所述的单分量数据进行零相位化处理，得到零相位化处理后的剖面。

在具体的实施例中，利用现有软件进行零相位化处理。零相位化处理的目的是将混合相位的复杂子波转化为零相位的子波，使子波的波峰或波谷与地层界面相对应，如图14所示，为本发明提供的具体实施例中地震子波零相位化示意图。

时距曲线确定装置104，用于根据所述零相位化处理后的剖面得到时距曲线。

图8为时距曲线确定装置的具体结构框图，由图8可知，该装置包括：

拾取位置确定模块201，用于将所述零相位化处理后的剖面中下行直达波的波峰位置作为拾取位置；

时距曲线确定模块202，用于以所述拾取位置拾取所述下行直达波，得到时距曲线。

在具体的实施方式中，利用现有地震处理软件拾取下行直达波初至，得到时间-深度关系序列，该序列在平面上构成一条时间与距离的关系曲线，以下简称为时距曲线，如图15所示，为本发明提供的具体实施例中VSP数据初至拾取示意图。

所述初至拾取位置为剖面中下行直达波的波峰位置。

所述下行波时距曲线同时作为后续TAR因子的时距曲线。

由图7可知，该系统还包括：

TAR因子确定装置105，用于根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子。图9为TAR因子确定装置的具体结构框图，由图9可知，该装置具体包括：

时间-振幅关系序列确定模块301，用于在所述零相位化处理后的剖面中，拾取所述时距曲线的振幅值，得到时间-振幅关系序列；

新序列确定模块302，用于对所述时间-振幅关系序列中的时间与振幅分别求取自然对数，得到新序列；

拟合直线确定模块303，用于对所述的新序列进行线性拟合，得到拟合直线；

拟合直线斜率确定模块304，用于确定所述拟合直线的斜率；

TAR因子确定模块305，用于对所述斜率取相反数，得到所述VSP数据的TAR因子。

在具体的实施方式中，在步骤S103所得到的剖面中，拾取步骤S104时距曲线位置的振幅值，获得时间-振幅关系序列，对该序列中时间与振幅分别求取自然对数，然后对所得到的新序列进行线性拟合，计算得到拟合直线的斜率，并对该斜率值取相反数，即得到垂直地震的TAR因子，如图17所示，为本发明提供的具体实施例中利用VSP数据计算的TAR因子的计算结果图。

由图7可知，该系统还包括：

拉齐装置106，用于将零相位化处理后的剖面沿所述时距曲线拉齐，得到拉齐后的剖面。

在具体的实施方式中，利用现有软件将步骤S103所得到的数据沿步骤S104所得到的时距曲线拉齐，获得准备分离下行波场的拉齐后的剖面，如图20所示，为本发明提供的具体实施例中沿下行波时距曲线拉平后的VSP数据示意图。

中值滤波装置107，用于根据所述VSP数据的TAR因子对所述拉齐后的剖面进行定向加权中值滤波。图10为中值滤波装置的具体结构框图，由图10可知，该装置具体包括：

地震道排序模块401，用于将所述拉齐后的剖面中地震道从最浅一级检波器记录向最深一级检波器记录排序；

中值滤波窗确定模块402，用于选定所述拉齐后的剖面中每次参与运算的地震道数，称为中值滤波窗。所述中值滤波窗在靠近拉齐后的剖面的边界时逐步减小，以保证中值滤波窗内不出现空样点值，使得后续处理结果相对保真。

离散数据序列确定模块403，用于将所述中值滤波窗内的振幅值提取出来，组成离散数据序列；

在具体的实施方式中，确定中值滤波窗(跨度)，即每次处理运算的最大样点数，在VSP数据剖面上也就是选定参与每一次统计计算的地震道数。

这些地震道上的同时刻样点值组成离散的数据序列。所述离散数据序列，要求从最浅一级检波器对应地震道向最深一级检波器对应的地震道方向排序。

在具体的实施方式中，确定用于中值滤波处理的最大样点数设为n，即选定参与每一次统计计算的地震道数，这些地震道上的同时刻样点值组成离散的数据序列X，如(1)式所示。

由时深关系可知，这n个样点对应的时间序列T如(2)式。

在具体的实施方式中，所述用于统计计算的地震道数取值范围为介于7和15之间的奇数，该参数的确定要综合分析波场分离效果和数据保真两个方面，具体数值视信剖面信噪比情况而确定，在信噪比相对比较低的地震数据中选用7作为样点数，在信噪比较高的情况下选用15作为样点数，信噪比介于两者之间时选用奇数9、11或13。该取值在靠近拉齐后的剖面的边界时逐步减小，以保证窗口内不出现空样点值，使得后续处理结果相对保真。

加权数据序列确定模块404，用于根据所述VSP数据的TAR因子以及所述中值滤波窗得到加权数据序列。在具体的实施方式中，利用步骤S105中所计算得到的TAR因子，以及步骤S402中所确定的滑动窗长度，设计加权函数，得到递增的加权数据序列。图11为加权数据序列确定模块的具体结构框图。由图11可知，该模块具体包括：

TAR因子确定单元501，用于根据所述VSP数据的TAR因子确定任意一个中值滤波窗内的TAR因子。利用步骤S105中所计算得到的TAR因子k，可以得到所述地震数据的剖面中任意一个时刻t的球面扩散补偿，如(3)式所示：

式中：A_t为任意一个时刻的振幅；A₀为零时刻的振幅；t表示样点所在的时刻，单位为秒(s)；k为TAR因子。

相对振幅关系补偿单元502，用于根据地震数据中任意一个中值滤波窗内的TAR因子对所述离散数据序列对应的中值滤波窗内的样点值进行相对振幅关系补偿。

在具体的实施方式中，由于本发明致力于保护原始的振幅信息，因此不能基于公式(3)进行整个剖面的补偿，但可以利用该式对步骤S502窗中的样点进行相对振幅关系的补偿。具体方法是选择窗的第个样点为基准点进行补偿处理，如(4)式所示。所选样点即为原始数据上准备替换的样点。该步骤是本发明的关键步骤，可以使原始震幅得到最好的保护。

对(4)式进行整理得到(5)式。

归一化处理单元503，用于对相对振幅关系补偿后的样点值进行基准点归一化处理。

在具体的实施方式中，以第个样点振幅值为进行基准点归一化处理，得到(5)式。

在(4)至(6)式中：A_i是在窗内以第个样点为基准点的理论补偿振幅；A′_i为A_i相对于基准点归一化后的值；t表示样点所在的时刻，单位为秒(s)；k为TAR因子；n为窗内最大样点数。

加权数据序列确定单元504，用于根据基准点归一化处理得到的样点值确定加权数据序列。

递增的加权数据序列Y，可以用(7)式表示。

由图10可知，该模块还包括：

定向加权模块405，用于对所述的离散数据序列以及加权数据序列进行定向加权，得到定向加权后的数据序列。

在具体的实施方式中，对序列X和序列Y进行样点对应相乘，即完成定向加权计算，得到新的数据序列Z，如式(8)所示。

以A′表示(8)式中各项，记作(9)。

在(8)式和(9)式中：Z为加权后新的振幅序列；X为原始的振幅序列；Y为振幅加权序列；a表示原始序列中单个样点的振幅；t表示样点所在的时刻；A′表示加权后单个样点的振幅；n为最大样点数。

去异常值处理模块406，用于对定向加权后的数据序列进行去异常值处理，得到去异常值处理后的数据序列。图12为去异常值处理模块的具体结构框图，由图12可知，该模块具体包括：

平均值确定单元601，用于确定定向加权后的数据序列中各个样点的绝对值的平均值；

去异常门槛值确定单元602，用于根据所述的平均值确定去异常门槛值。在具体的实施方式中，取数列中各个样点平均值的绝对值，并将该值加减20％作为一对去异常门槛值。

异常样点确定单元603，用于从定向加权后的数据序列的各个样点中筛选出绝对值在所述去异常门槛值以外的样点，作为异常样点。在具体的实施方式中，将原数列中各个样点的绝对值与去异常门槛值进行逐一比较，筛选出绝对值在所述去异常门槛值以外的样点，作为异常样点。

异常样点替换单元604，用于以定向加权后的数据序列中各个样点的绝对值的平均值替换所述异常样点，得到异常值处理后的数据序列。

也即本发明的去异常值处理，具体方法是取数列中各个样点平均值的绝对值，并将该值加减20％作为一对去异常的门槛值，再将原数列中各个样点的绝对值与门槛值进行逐一比较，并以平均值替换门槛值范围以外的样点值。如式(10)所示。

此时，自适应去异常后的振幅序列Z′为(11)所示。

在(10)式和(11)式中：Z′为自适应去异常后的振幅序列；A′表示加权后单个样点的振幅；A″表示自适应去异常后单个样点的振幅；n为最大样点数。

图19至图22中示意了一个样点窗的加权、自适应去异常的过程。其中，图19为本发明提供的具体实施例中一个窗内的样点序列示意图，图20为本发明提供的具体实施例中依VSP数据规律的加权系数序列示意图，图21为本发明提供的具体实施例中加权之后一个窗内的样点序列示意图，图22为本发明提供的具体实施例中自适应去异常后的样点序列示意图。

所述门槛值设定为均值的20％，该值远低于其它地震数据自适应去异常方法中的门槛值。这是由于VSP数据检波器深入井下，接收环境相对安静，因而出现异常高值或低值的地震野值情况较少；并且经过样点加权处理后，正常值范围内振幅变得更加均衡。

在本发明的其他实施例中，当窗滑动至边界位置时，自动减小窗长，以满足计算需要，回避边界的异常处理效应，该技术为公知技术，已应用于很多商业地震处理软件中。

如上即完成了准备工作，剖面中每一点便对应一个窗。

由图10可知，该模块还包括：

中值滤波模块407，用于对去异常值处理后的数据序列进行中值滤波处理，得到定向加权中值滤波处理数据。

在本发明的具体实施例中，该步骤的具体算法是将窗中所有样点按照从大到小的规则重新排序，并取得新序列的中间一个样点的数值，用以代替剖面中的该点，从而完成自适应去异常的、边界变道数的、依照VSP数据规律的定向加权中值滤波处理，如图23为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图1，图24为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图2。

循环步骤S401至步骤S407的处理流程，可完成多组下行波中值滤波处理。对上行波进行中值滤波处理时，拾取上行波时距曲线；同时使用等比递减的加权序列，该变化是由VSP数据方法中地震波传播特点所决定的，其上行波振幅随接收检波器深度的变浅而减小，但衰减规律不变，不需要重新计算TAR因子。

如上即是本发明在自适应去异常的、边界变道数的中值滤波基础上，提出一种利用VSP数据衰减规律进行加权的全新中值滤波系统，解决用于井控的VSP数据数据处理中没有可靠波场分离方法的问题。

下面结合具体的实施方式，详细介绍本发明的技术方案。本发明的目的是提出了一种针对VSP数据衰减规律进行加权的中值滤波方法，解决用于井控的VSP数据处理中没有可靠波场分离方法的问题。是在自适应去异常的、边界变道数的中值滤波基础上，利用VSP数据的衰减规律进行加权处理的全新中值滤波方法。

在具体的实施方式中，该方案具体包括：

1)数据准备：

利用现有方法和地震处理软件对VSP数据原始采集数据进行整理，获得准备进行波场分离的单分量数据，如图13所示，为本发明提供的具体实施例中原始VSP数据Z分量示意图。为了表现VSP数据振幅关系，该数据未加增益显示(下同)。

所述的数据整理不包括振幅补偿处理、道均衡处理、自动增益(AGC)处理或其它可能改变振幅或能量相对关系的处理步骤。

2)零相位化处理：

利用现有软件进行零相位化处理。零相位化处理的目的是将混合相位的复杂子波转化为零相位的子波，使子波的波峰或波谷与地层界面相对应，如图14所示，为本发明提供的具体实施例中地震子波零相位化示意图。

3)时间-深度关系拾取：

利用现有地震处理软件拾取下行直达波初至，得到时间-深度关系序列，该序列在平面上构成一条时间与距离的关系曲线，以下简称为时距曲线，如图15所示，为本发明提供的具体实施例中VSP数据初至拾取示意图。

所述初至拾取位置为剖面中下行直达波的波峰位置。

4)球面扩散补偿因子(TAR因子)计算：

在步骤2)所得到的剖面中，拾取步骤3)时距曲线位置的振幅值，获得时间-振幅关系序列，对该序列中时间与振幅分别求取自然对数，然后对所得到的新序列进行线性拟合，计算得到拟合直线的斜率，并对该斜率值取相反数，即得到垂直地震的TAR因子，如图18为本发明提供的具体实施例中利用VSP数据计算的TAR因子的计算结果图，图17为本发明提供的具体实施例中TAR的计算差异曲线图。用来监测TAR计算的精度。

5)时距曲线拉齐：

利用现有软件将2)步骤所得到的数据沿3)步骤所得到的时距曲线拉齐，获得准备分离下行波场的地震数据，如图18所示，为本发明提供的具体实施例中沿下行波时距曲线拉平后的VSP数据示意图。

6)以下述方式进行中值滤波处理：

(1)确定用于中值滤波处理的最大样点数n，即选定参与每一次统计计算的地震道数，这些地震道上的同时刻样点值组成离散的数据序列X，如(1)式所示。

由步骤3)的时深关系可知，这n个样点对应的时间序列T如(2)式所示。

所述用于统计计算的地震道数取值范围为介于7和15之间的奇数，该参数的确定要综合分析波场分离效果和数据保真两个方面，具体数值视信剖面信噪比情况而确定，在信噪比相对比较低的地震数据中选用7作为样点数，在信噪比较高的情况下选用15作为样点数，信噪比介于两者之间时选用奇数9、11或13。

所述离散数据序列，要求从最浅一级检波器对应地震道向最深一级检波器对应的地震道方向排序。

(2)利用步骤4)中所计算得到的TAR因子k，可以得到剖面中任意一个时刻t的球面扩散补偿，如(3)式所示。

由于本发明致力于保护原始的振幅信息，因此不能基于公式(3)进行整个剖面的补偿，但可以利用该式对(1)步骤窗中的样点进行相对振幅关系的补偿。具体方法是选择窗的第个样点为基准点进行补偿处理，如(4)式所示。所选样点即为原始数据上准备替换的样点。该步骤是本发明的关键步骤，可以使原始震幅得到最好的保护。

对(4)式进行整理得到(5)式，再以第个样点振幅值为进行基准点归一化处理，得到(6)式。在(4)至(6)式中：A_i是在窗内以第个样点为基准点的理论补偿振幅；A′_i为A_i相对于基准点归一化后的值；t表示样点所在的时刻，单位为秒(s)；k为TAR因子；n为窗内最大样点数。

于是递增的加权数据序列Y，如(7)式所示。

(3)对(1)步骤中序列X和(2)步骤中序列Y进行样点对应相乘，即完成定向加权计算，得到新的数据序列Z，如式(8)所示。以A′表示(8)式中各项，记作(9)式。

在(8)式和(9)式中：Z为加权后新的振幅序列；X为原始的振幅序列；Y为振幅加权序列；a表示原始序列中单个样点的振幅；t表示样点所在的时刻；A′表示加权后单个样点的振幅；n为最大样点数。

(4)对上一步所得到的序列进行去异常值处理，具体方法是取数列中各个样点平均值的绝对值，并将该值加减20％作为一对去异常的门槛值，再将原数列中各个样点的绝对值与门槛值进行逐一比较，并以平均值替换门槛值范围以外的样点值。如式(10)所示。

此时，自适应去异常后的振幅序列Z′如(11)式所示。

在(10)式和(11)式中：Z′为自适应去异常后的振幅序列；A′表示加权后单个样点的振幅；A″表示自适应去异常后单个样点的振幅；n为最大样点数。

在图19至图22中示意了一个样点窗的加权、自适应去异常的过程。图19为本发明提供的具体实施例中一个窗内的样点序列示意图，图20为本发明提供的具体实施例中依VSP数据规律的加权系数序列示意图，图21为本发明提供的具体实施例中加权之后一个窗内的样点序列示意图，图22为本发明提供的具体实施例中自适应去异常后的样点序列示意图。

所述门槛值设定为均值的20％，该值远低于其它地震数据自适应去异常方法中的门槛值。这是由于VSP数据检波器深入井下，接收环境相对安静，因而出现异常高值或低值的地震野值情况较少；并且经过样点加权处理后，正常值范围内振幅变得更加均衡。

(5)当窗滑动至边界位置时，自动减小窗长，以满足计算需要，回避边界的异常处理效应，该技术为公知技术，已应用于很多商业地震处理软件中。

(6)完成步骤(1)至步骤(5)的准备工作，剖面中每一点便对应一个窗，最后一步是进行中值滤波计算，具体算法是将窗中所有样点按照从大到小的规则重新排序，并取得新序列的中间一个样点的数值，用以代替剖面中的该点，从而完成自适应去异常的、边界变道数的、依照VSP数据规律的定向加权中值滤波处理，图23为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图1，图24为本发明提供的具体实施例中定向加权的、自适应去异常的、边界变道数的中值滤波后VSP数据示意图2。

7)完成多组波场的分离：

循环步骤3)至步骤6)的处理流程，可完成多组下行波中值滤波处理。如果对上行波进行中值滤波处理，则需要在步骤3)中，拾取上行波时距曲线；同时在步骤6)中使用等比递减的加权序列，该变化是由VSP数据方法中地震波传播特点所决定的，其上行波振幅随接收检波器深度的变浅而减小，但衰减规律不变，不需要重新计算TAR因子。

综上所述，在自适应去异常的、边界变道数的中值滤波基础上，本发明提供了一种基于垂直地震衰减规律的中值滤波方法及系统，是一种精确的基于垂直地震衰减规律的中值滤波的方案，由VSP数据获得单分量数据，对单分量数据进行零相位化处理，结合时距曲线确定TAR因子，并沿时距曲线拉齐剖面，进行定向加权中值滤波，极大提高了中值滤波的可靠性，为用于井控的VSP数据数据处理提供了可靠的波场分离方案。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一般计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random AccessMemory，RAM)等。

本领域技术人员还可以了解到本发明实施例列出的各种功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (12)

1.一种基于垂直地震衰减规律的中值滤波方法，其特征是，所述的方法包括：

采集垂直地震剖面VSP数据；

根据所述的VSP数据分选出单分量数据；

对所述的单分量数据进行零相位化处理，得到零相位化处理后的剖面；

根据所述零相位化处理后的剖面得到时距曲线；

根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子；

将零相位化处理后的剖面沿所述时距曲线拉齐，得到拉齐后的剖面；

根据所述VSP数据的球面扩散补偿TAR因子对所述拉齐后的剖面进行定向加权中值滤波。

2.根据权利要求1所述的方法，其特征是，根据所述零相位化处理后的剖面得到时距曲线包括：

将所述零相位化处理后的剖面中下行直达波的波峰位置作为拾取位置；

以所述拾取位置拾取所述下行直达波，得到时距曲线。

3.根据权利要求2所述的方法，其特征是，根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子包括：

在所述零相位化处理后的剖面中，拾取所述时距曲线的振幅值，得到时间-振幅关系序列；

对所述时间-振幅关系序列中的时间与振幅分别求取自然对数，得到新序列；

对所述的新序列进行线性拟合，得到拟合直线；

确定所述拟合直线的斜率；

对所述斜率取相反数，得到所述VSP数据的球面扩散补偿TAR因子。

4.根据权利要求3所述的方法，其特征是，根据所述VSP数据的球面扩散补偿TAR因子对所述的拉齐后的剖面进行定向加权中值滤波处理包括：

将所述拉齐后的剖面中地震道从最浅一级检波器记录向最深一级检波器记录排序；

选定所述拉齐后的剖面中每次参与运算的地震道数，称为中值滤波窗；

将所述中值滤波窗内的振幅值提取出来，组成离散数据序列；

根据所述VSP数据的球面扩散补偿TAR因子以及所述中值滤波窗得到加权数据序列；

对所述的离散数据序列以及加权数据序列进行定向加权，得到定向加权后的数据序列；

对定向加权后的数据序列进行去异常值处理，得到去异常值处理后的数据序列；

对去异常值处理后的数据序列进行中值滤波处理，得到定向加权中值滤波处理数据。

5.根据权利要求4所述的方法，其特征是，根据所述VSP数据的球面扩散补偿TAR因子以及所述中值滤波窗得到加权数据序列包括：

根据所述VSP数据的球面扩散补偿TAR因子确定任意一个中值滤波窗内的球面扩散补偿TAR因子；

根据任意一个中值滤波窗内的球面扩散补偿TAR因子对所述离散数据序列对应的中值滤波窗内的样点值进行相对振幅关系补偿；

对相对振幅关系补偿后的样点值进行基准点归一化处理；

根据基准点归一化处理得到的样点值确定加权数据序列。

6.根据权利要求4所述的方法，其特征是，对定向加权后的数据序列进行去异常值处理，得到去异常值处理后的数据序列包括：

确定定向加权后的数据序列中各个样点的绝对值的平均值；

根据所述的平均值确定去异常门槛值；

从定向加权后的数据序列的各个样点中筛选出绝对值在所述去异常门槛值以外的样点，作为异常样点；

以定向加权后的数据序列中各个样点的绝对值的平均值替换所述异常样点，得到异常值处理后的数据序列。

7.一种基于垂直地震衰减规律的中值滤波系统，其特征是，所述的系统包括：

剖面采集装置，用于采集垂直地震剖面VSP数据；

单分量数据分选装置，用于根据所述的VSP数据分选出单分量数据；

零相位化处理装置，用于对所述的单分量数据进行零相位化处理，得到零相位化处理后的剖面；

时距曲线确定装置，用于根据所述零相位化处理后的剖面得到时距曲线；

球面扩散补偿TAR因子确定装置，用于根据所述零相位化处理后的剖面以及所述的时距曲线确定所述VSP数据的球面扩散补偿TAR因子；

拉齐装置，用于将零相位化处理后的剖面沿所述时距曲线拉齐，得到拉齐后的剖面；

中值滤波装置，用于根据所述VSP数据的球面扩散补偿TAR因子对所述的拉齐后的剖面进行定向加权中值滤波。

8.根据权利要求7所述的系统，其特征是，所述的时距曲线确定装置包括：

拾取位置确定模块，用于将所述零相位化处理后的剖面中下行直达波的波峰位置作为拾取位置；

时距曲线确定模块，用于以所述拾取位置拾取所述下行直达波，得到时距曲线。

9.根据权利要求8所述的系统，其特征是，所述的球面扩散补偿TAR因子确定装置包括：

时间-振幅关系序列确定模块，用于在所述零相位化处理后的剖面中，拾取所述时距曲线的振幅值，得到时间-振幅关系序列；

新序列确定模块，用于对所述时间-振幅关系序列中的时间与振幅分别求取自然对数，得到新序列；

拟合直线确定模块，用于对所述的新序列进行线性拟合，得到拟合直线；

拟合直线斜率确定模块，用于计算所述拟合直线的斜率；

球面扩散补偿TAR因子确定模块，用于对所述斜率取相反数，得到所述VSP数据的球面扩散补偿TAR因子。

10.根据权利要求9所述的系统，其特征是，所述的中值滤波装置包括：

地震道排序模块，用于将所述拉齐后的剖面中地震道从最浅一级检波器记录向最深一级检波器记录排序；

中值滤波窗确定模块，用于选定所述拉齐后的剖面中每次参与运算的地震道数，称为中值滤波窗；

离散数据序列确定模块，用于将所述中值滤波窗内的振幅值提取出来，组成离散数据序列；

加权数据序列确定模块，用于根据所述VSP数据的球面扩散补偿TAR因子以及所述中值滤波窗得到加权数据序列；

定向加权模块，用于对所述的离散数据序列以及加权数据序列进行定向加权，得到定向加权后的数据序列；

去异常值处理模块，用于对定向加权后的数据序列进行去异常值处理，得到去异常值处理后的数据序列；

中值滤波模块，用于对去异常值处理后的数据序列进行中值滤波处理，得到定向加权中值滤波处理数据。

11.根据权利要求10所述的系统，其特征是，所述的加权数据序列确定模块包括：

球面扩散补偿TAR因子确定单元，用于根据所述VSP数据的球面扩散补偿TAR因子确定任意一个中值滤波窗内的球面扩散补偿TAR因子；

相对振幅关系补偿单元，用于根据任意一个中值滤波窗内的球面扩散补偿TAR因子对所述离散数据序列对应的中值滤波窗内的样点值进行相对振幅关系补偿；

归一化处理单元，用于对相对振幅关系补偿后的样点值进行基准点归一化处理；

加权数据序列确定单元，用于根据基准点归一化处理得到的样点值确定加权数据序列。

12.根据权利要求10所述的系统，其特征是，所述的去异常值处理模块包括：

平均值确定单元，用于确定定向加权后的数据序列中各个样点的绝对值的平均值；

去异常门槛值确定单元，用于根据所述的平均值确定去异常门槛值；

异常样点确定单元，用于从定向加权后的数据序列的各个样点中筛选出绝对值在所述去异常门槛值以外的样点，作为异常样点；

异常样点替换单元，用于以定向加权后的数据序列中各个样点的绝对值的平均值替换所述异常样点，得到异常值处理后的数据序列。