一种利用地表一致性统计频谱分析技术处理海量地震数据的方法

摘要

本发明是石油地球物理勘探数据处理中对数据处理，涉及利用地表一致性统计频谱分析技术处理海量地震数据的方法。先将数据分选为炮集，在分析的时窗段上取100ms以内的扰动，形成数据段，在时窗的两端加线性斜坡，并补零至原始数据的道长，对扰动后的数据求平均，得到振幅谱，完成所有炮记录的谱分析和所有检波点集记录的谱分析，将整条测线的谱估计结果与采集的地表结构相对应得到规律性，选择针对性的海量数据处理方法，可提高处理精度与速度、降低噪音，得到高精度的地层剖面及确定含油地层位置。

说明书

技术领域

本发明属于石油地球物理勘探技术，涉及高密度地震勘探数据处理中对数据频率进行快速分析的方法，具体涉及利用地表一致性统计频谱分析技术处理海量地震数据的方法。

背景技术

众所周知，在石油地球物理勘探技术中，提高地震勘探精度、解决陆相沉积盆地复杂的地质问题是一项系统工程，涉及采集方法、静校正技术、噪音压制技术、精细成像技术等许多方面。高密度空间采样技术是地震勘探领域新的发展方向，对于改善资料信噪比和分辨率、提高资料品质、获得更丰富的地下反射信息有独特的作用，因此成为当今地震勘探的热点之一。高密度空间采样技术是指在野外实施点源激发、点接收或小面积组合接收、小道距或小面元观测，对地震波场进行充分采样，对信号和噪声实行“宽进宽出”，避免采集过程中因检波器组合而使反射信息受到伤害，尽量保持采样波场的原始性和丰富性。

但是，对于野外通过高密度空间采样获得的庞大数据体，如何采用快速有效的方法对数据特征进行分析是处理人员首先面临的问题；其次面临的问题的是，在处理过程中，处理人员如何将主要精力放在正确认识数据的特征、客观分析数据存在的问题，在该问题中如何采取适当的技术手段、选择合理的参数对海量数据实施处理，使海量数据的质量不断提高，使处理成果满足地质综合研究的需要；再者面临的问题的是，如何掌握分析手段、质控资料或数据图件监控处理中的质量问题，研究这些问题于近地表特征是否有相关性。拥有的快速地表一致性分析手段是地震数据处理中心综合技术实力的体现。采用地表一致性分析技术掌握海量数据的特征，使用定量分析方法对海量数据处理过程中各个环节的处理质量进行及时、有效地监控，观察数据的细节，可以及时消除隐患，提高工作效率，一直是石油地球物理勘探技术中亟需解决的问题。

众所周知，频率是地震数据的重要特征之一。为了解决上面所述技术问题，在目前常用的处理系统中，多采用交互的单道或多道谱分析方法对地震数据进行频率分析，但是在现有的鲜见的海量数据处理技术中，存在着各种各样的不足或缺陷。下面以OMEGA地震数据处理系统为例加以说明(参见图1)。在图1中所示的图形结构中，交互谱估计的显示分三个区域：(a)输入数据(T-X)区；(b)对应道F-X谱；(c)多道振幅谱。具体方法步骤为：(1)选择输入数据；(2)进行Fourier变换；(3)输出交互谱(F-X谱和振幅谱)。

对于野外高密度采集的海量地震数据而言，常规的谱分析技术(如OMEGA地震数据处理系统)主要在以下几方面存在不足：

(1)当需要处理的数据量大，且数据的频率特征在时间、空间方向存在差异时，常规交互方式谱分析方法需要人工不断修改输入数据，输出的速度降低，不利于快速对比分析；

(2)没有考虑地震勘探中的地表一致性假设前提，不易得到激发、接收条件变化对数据特征造成的影响；

(3)F-X谱受所取数据长度的限制，显示方式无法调整，精度较低；

(4)未能提取频谱上的特征参数，无法实现量化对比分析。

另外，由于地震数据处理的特殊性，即地震数据处理是在全面了解原始数据整体分布特征的基础上开展的针对性处理，数据量和处理精度要求、数据分布特征、分析的难度和分析范围均与其他非海量数据处理方法不同；又由于对高密度地震勘探来说，采集的数据量又增加若干倍，更需要有与之相适应的地表一致性数据分析技术。

因此，对于高密度地震勘探技术来说，目前急需一种合适的海量地震数据的方法，以此来提高海量地震处理精度与速度、降低噪音，籍此通过识别有效区带来得到高精度的地层剖面及确定含油地层位置。

发明内容

针对现有海量地震数据系统的不足，本发明针对高密度采集中叠前数据量大的特点，本发明提出一种利用地表一致性统计频谱分析技术处理海量地震数据的方法，该方法能够实施快速有效的地震数据谱估计且可根据需要完成地表一致性统计数据的谱估计。

依据本发明的利用地表一致性统计频谱分析技术处理海量地震数据的方法包括以下步骤：

(1)在野外初步确定的有含油气希望的地区，布置测线和炮位，使用爆炸源或地面可控震源激发地震波，并用地震仪把地震波的传播情况记录下来；

(2)首先研究炮集数据上的频率变化规律：将数据分选为炮集，根据研究的目的在一炮数据上选择时窗，地震数据为x_k(i，j)；i为时间序号，j为道序号，k为炮序号，其中每一个单道记录为x(i)；

(3)在分析的时窗段上，取100ms以内的扰动，形成m个数据段，m一般取3-5；根据时窗长度(Wb-Wa+10×n)/SI的样点数(n为0-10之间的随机数，由软件自动产生，保证扰动在100ms以内；SI为采样间隔)，在时窗的两端加线性斜坡，并补零至原始数据的道长；对扰动后的m个时窗段的数据做Fourier变换，采用公式(1)计算单道数据的振幅谱，对扰动后的数据求平均，得到x(i)的振幅谱A(f)；

$A\left(f\right)=\frac{1}{m}\underset{j=\mathrm{1,2},....m}{\Sigma }\left[\mathrm{FFT}\right|\left(x_j\left(t_i\right)\right|)]$(1)

式中：x_j(t_i)为输入地震道，

      t_i为时间，

      j＝1，2，……m为随机扰动后的道序号，

      A(f)为谱估计的结果，

      f为频率；

据此计算出炮集记录在所选范围内各道的Ak(f，j)，计算每炮数据的振幅谱：

描述集总特征的统计平均，即数学期望，如果数据呈正态分布，可采用下面方式(2)

A_k(f)＝E[A_k(f，j)]        (2)

式中：A_k(f，j)为第k炮、第j道的谱估计结果，

      k为炮序号，

      j为道序号，

      A_k(f)为谱估计的结果，

      f为频率。

(4)将一炮记录的谱估计的结果存放在结果道中；

(5)取下一炮记录，重复(3)、(4)，直到完成所有炮记录的谱分析；

(6)其次，将数据分选到检波点集，研究检波点集上数据频率的变化情况：根据分析的目的在同一检波点数据上选择时窗，地震数据为x₁(i，j)；i为时间序号，j为道序号，1为检波点序号，

(7)检波点集上数据谱估计的过程同(3)-(4)；

(8)取下一检波点集数据，重复(3)-(4)，直到完成所有检波点集记录的谱分析；

(9)对于每炮时间域记录(由若干道组成)或每个检波点时间域记录(由若干道组成)，通过上述计算后，该炮点或该检波点的频率信息就可以用1道来表示，压缩了数据量。就可以很方便地将整条测线或整个工区的谱估计结果与野外采集时的地表结构相对应，易于发现其规律性，可快速检查采集中的问题。在处理过程中的各个环节，也采用地表一致性的谱分析方式检查过程中频率的变化、结果的保真度、异常现象出现的区域等，提高数据检测的速度和准确度。

其中在上述公式(2)中，可选择相对稳健的中值，如果数据中有异常值分布，可采用公式(3)

A_k(f)＝median[A_k(f，j)]         (3)

亦可选择灵活的α-TRIM估计，如果数据中有异常值分布，又要取平均，可采用公式(4)

A_k(f)＝α-trim[A_k(f，j)]        (4)

α-TRIM估计中，如果参数选择为20％，即将一炮数据所选范围内的各道的谱在同一频率样点上排序，去掉最小的10％，去掉最大的10％，将中间的80％数据平均作为该频率样点上的值，逐个频率计算，直至截止频率，得到该炮数据的频率值；

在上述利用地表一致性统计频谱分析技术处理海量地震数据的方法中，地表一致性统计谱估计的计算结果以数据的形式存放在地震道中，这样就可以对谱数据进行处理，以便提取特征参数。

对数据进行处理方法：沿频率方向的积分滤波或中值滤波，避免异常值对特征参数的影响。

A′_k(f)＝integralfilter[A_k(f)]        (5)

A′_k(f)＝medianfilter[A_k(f)]          (6)

式中：A′_k(f)为滤波后的数据，

      A_k(f)为输入数据。

特征参数的提取为：对滤波处理后的数据在不同的频带范围内进行特征参数的提取。频带范围可分为：干扰波频带范围、反射波的频带范围、优势频带范围、低频段、高频段等，根据需要提取某一频带内的平均能量(Aeng)和峰值频率(fmax)，配以关键字存于道头中，绘制出量化特征曲线或曲面。

频带的平均能量(Aeng)曲线能够清楚地反映该频带内的能量横向变化，分析能量变化与近地表结构变化的对应关系，帮助处理人员选择试验点，确定试验方案。峰值频率(fmax)是所选频带内最大值对应的频率，这一参数使处理人员能快速地了解数据的主频分布以及数据的主频与近地表结构的对应关系，在数据处理的不同环节应用本参数可以研究处理过程中数据的频率变化规律，研究近地表条件的变化对地震数据的影响是否已消除，使得成果数据的频率信息反映的是地下储层的特征或含油气的情况。如果是三维数据，通过本参数的平面显示或立体显示快速分析全区数据的空间分布情况，以便：

(1)对比分析地震数据的频率分布与地表高程、低速带速度、低速带的厚度的关系，研究激发、接收条件对地震数据品质的影响；

(2)根据地表一致性统计谱特征参数分布情况，确定处理中的空变参数和频变参数；

(3)监控处理过程中频率的变化。

附图说明

图1为OMEGA系统中的谱分析的示意图；

图2为使用本发明方法的谱分析结果；

图3(a)为原始数据炮集统计谱分析的示意图；

图3(b)为叠前去噪后炮集统计谱分析的示意图；

图3(c)为振幅处理后炮集统计谱分析的示意图；

图3(d)为反褶积后炮集统计谱分析的示意图；

图4(a)为原始检波点集统计谱分析的示意图；

图4(b)为叠前去噪后检波点集统计谱分析的示意图；

图4(c)为反褶积后检波点集统计谱分析的示意图；

图4(d)为最终检波点集统计谱分析的示意图；

图5(a)为原始检波点集谱分析的示意图，图上曲线为40-60Hz的峰值频率；

图5(b)为处理后的检波点集谱分析的示意图，图上曲线为40-60Hz的峰值频率；

图6(a)为初步叠加剖面的示意图；

图6(b)为最终叠加剖面的示意图。

具体实施方式

下面结合附图，对本发明提出的利用地表一致性统计频谱分析技术处理海量地震数据的方法进行详细描述。

首先，在野外初步确定的有含油气希望的地区，布置测线和炮位，使用爆炸源或地面可控震源激发地震波，并用地震仪把地震波的传播情况记录下来。也就是，使用常规技术手段，在野外采集将要处理的海量地震数据。

首先研究炮集数据上的频率变化规律：将数据分选为炮集，根据研究的目的在一炮数据上选择时窗，地震数据为x_k(i，j)，单道数据为x(i)；

其次，在分析的时窗段上，取100ms以内的扰动，形成m个数据段，m一般取3-5；根据时窗长度(Wb-Wa+10×n)/SI的样点数(n为0-10之间的随机数，Wb-Wa为时窗范围，SI为采样间隔)，在时窗的两端加线性斜坡，并补零至原始数据的道长；对扰动后的m个时窗段的数据做Fourier变换，采用公式(1)计算单道数据的振幅谱，对扰动后的数据求平均，得到x(i)的振幅谱A(f)；

$A\left(f\right)=\frac{1}{m}\underset{j=\mathrm{1,2},....m}{\Sigma }\left[\mathrm{FFT}\right|\left(x_j\left(t_i\right)\right|)]$(1)

式中：x_j(t_i)为输入地震道，

      t_i为时间，

      j＝1，2，……m为随机扰动后的道序号，

      A(f)为谱估计的结果，

      f为频率。

再者，计算出炮集记录在所选范围内各道的A_k(f，j)，可根据数据的特征采用一炮数据谱的估计方式：

描述集总特征的统计平均，即数学期望：(如果数据呈正态分布，可采用本方式)

A_k(f)＝E[A_k(f，j)]        (2)

式中：A_k(f，j)为第k炮、第j道的谱估计结果，

      k为炮序号，

      j为道序号，

      A_k(f)为谱估计的结果，

      f为频率。

可选择相对稳健的中值：(如果数据中有异常值分布，可采用此方式)

A_k(f)＝median[A_k(f，j)]          (3)

亦可选择灵活的α-TRIM估计：(如果数据中有异常值分布，又要取平均，可采用此方式)

A_k(f)＝α-trim[A_k(f，j)]        (4)

α-TRIM估计中，如果参数选择为20％，即将一炮数据所选范围内的各道的谱在同一频率样点上排序，去掉最小的10％，去掉最大的10％，将中间的80％数据平均作为该频率样点上的值，逐个频率计算，直至截止频率，得到该炮数据的频率值。

然后，将一炮记录的谱估计的结果存放在结果道中；

对每一炮记录都重复上述运算，直到完成所有炮记录的谱分析；

将数据分选到检波点集，研究检波点集上数据频率的变化情况：根据分析的目的在同一检波点数据上选择时窗，地震数据为x₁(i，j)；

检波点集上数据谱估计的过程与炮集相同；

对每一检波点集记录都重复上述运算，直到完成所有检波点集记录的谱分析。

对于每炮时间域记录(由若干道组成)或每个检波点时间域记录(由若干道组成)，通过上述计算后，该炮点或该检波点的频率信息就可以用1道来表示，压缩了数据量。就可以很方便地将整条测线或整个工区的谱估计结果与野外采集时的地表结构相对应，易于发现其规律性，可快速检查采集中的问题。在处理过程中的各个环节，也采用地表一致性的方式检查过程中频率的变化、结果的保真度、异常现象出现的区域等，提高检测的准确度。

在地表一致性统计谱估计中，计算结果以数据的形式存放在地震道中，这样就可以对数据进行处理，以便提取特征参数。

对数据进行处理方法为：沿频率方向的积分滤波或中值滤波，避免异常值对特征参数的影响。

A′_k(f)＝integralfilter[A_k(f)]    (5)

A′_k(f)＝medianfilter[A_k(f)]      (6)

式中：A′_k(f)为滤波后的数据，

      A_k(f)为输入数据。

特征参数的提取为：对滤波处理后的数据在不同的频带范围内进行特征参数的提取。频带范围可分为：干扰波频带范围、反射波的频带范围、优势频带范围、低频段、高频段等，根据需要提取某一频带内的平均能量(Aeng)和峰值频率(fmax)，配以关键字存于道头中，绘制出量化特征曲线或曲面。

某个频带的平均能量(Aeng)曲线可以清楚地反映该频带内的能量横向变化，分析能量变化与近地表结构变化的对应关系，帮助处理人员选择试验点，确定试验方案。峰值频率(fmax)是所选频带内最大值对应的频率，这一参数使处理人员能快速地了解数据的主频分布，在数据处理的不同环节应用本参数可以研究处理过程中数据的频率变化规律，如果是三维数据，通过本参数的平面显示快速分析全区数据的空间分布情况，确定处理中的空变参数和频变参数，监控处理过程中频率的变化。

另外，本发明的方法还采用以下技术。

1.时窗的选取

在地表一致性统计谱估计中，炮集或检波点集上时窗的选取原则：(1)研究面波的频率特征：时窗取在小炮检距道上，可取不同的时间段研究面波的频散现象，一般取400-800ms，1200-2200ms，2500-3800ms，4000-6000ms；(2)研究反射波的频率特征：时窗避开面波区和线性干扰区，选窗的范围视记录中反射波出现的位置而定，应是时间和炮检距的函数；(3)目的层段的特征：根据探区的目的层所在的位置，确定时变、空变的时窗，研究目的层段上有效信号和干扰波的频率分布状态。

2.频带范围的选取

采用本发明进行特征参数提取时，通常在一定的频带内进行，频带的选取视基础数据的频带宽度和数据中地震波的类型而定。如分析面波的频率分布规律时，频带范围可选在3-20Hz；但在分析有效反射波的频率特征时，可以倍频程的方式选择频带范围，即：15-30Hz，20-40Hz，30-60Hz，40-80Hz，50-100Hz等；也可以选择在优势频带的范围内分析，即：20-60Hz，30-80Hz等；对于特殊的干扰，可有针对性地选择频带范围，如高频干扰可取100-180hz；对于工业电干扰，可取45-55Hz或40-60Hz。

使用本发明方法得到的地震数据处理的谱结果具有以下优点：

(1)谱精度高，因为采用了随机扰动的方式统计计算谱值，并将振幅谱值以数据的方式存于地震道中(注意采样间隔的不同)，便于采用彩色显示或可视化的方式与原始T-X域的数据对比分析；

(2)提取的特征参数稳健，能代表所分析数据集(炮集或检波点集)的综合特征，且以数值的方式存放于道头中，可以曲线、曲面或彩色平面的方式显示，以便与其他特征参数(如地表高程曲线、低降速带速度及厚度、静校正量曲线、能量分布曲线等)对比，进行综合研究。

(3)谱分析速度快，分析人员确定分析范围、参数后，地表一致性谱估计和统计由程序自动完成，适用于高密度采集的海量数据的频率分析或大面积的三维数据分析。

下面对附图进行详细解释。

图1为OMEGA系统中的谱分析图件，图1(a)中横坐标为道号(TRACE NUM)，纵坐标为时间(Time(ms))，图1(b)中横坐标为频率(Frequency(Hz))，纵坐标为道号(TRACENUM)；图1(c)中横坐标为频率(Frequency(Hz))纵坐标为振幅(Amplitude)。可以统计出输入炮记录的振幅谱，这种方式适合于点分析，不能方便地进行炮与炮之间的对比，振幅谱仅以曲线的方式绘出，炮之间或道之间的频率与能量的量化关系难以体现出来。对于高密度采集的海量数据，叠前分析的工作量很大，需要进行横向和空间方向的对比，快速稳健的谱分析和显示是处理参数分析和选择的依据。

图2为本方法处理后的一条线共检波点集谱分析结果。横坐标是检波点号(SRF_SLOC)，纵坐标是频率(Frequency(Hz))。可以看出在测线上560个检波点的频率变换情况：在检波点2951-2970过水泡子，冬季施工产生严重的冰上干扰，干扰波的频带宽可达200Hz，能量强；检波点3000-3260段接收的数据频带较宽，但也受地表高程变化的影响；检波点3300至测线尾部，激发点和接收点均在高部位，导致地震数据频带变窄、能量减弱。

通过对共接收点集的统计谱分析，全测线上数据频率和能量的变化情况及其与地表条件的对应关系清晰，一目了然，可以达到快速分析的目的。

应用实例：GA地区保真处理中的谱分析

按照该区资料处理任务的要求，地震资料要进行精细的叠前时间偏移处理，要求剖面具有较高的信噪比和分辨率，能满足后续储层预测及烃类检测处理的需要。

由于该区地质任务对数据处理的保真度要求较高，质量监控和数据分析极为重要。常规的分析方法往往找不准问题的症结所在，且分析范围不全面。采用本方法对数据处理过程进行了全程跟踪监控后，通过诊断，找到了处理中各个环节的问题，发现了其他软件不易发现的问题，即某些频带上的特殊干扰影响了正常处理的效果，降低了处理的精度。为此，我们采取了针对性的迭代分频去噪技术，重新调整了处理流程，取得了良好的效果。

图3(a)、(b)、(c)、(d)为处理过程中炮集数据的频谱跟踪分析的结果，横坐标是炮点号(SOU_SLOC)，纵坐标是频率(Frequency(Hz))。从谱上可以清楚地看出2310-2420炮之间严重的50Hz、150Hz干扰波，横向炮间频率、能量的差异可以快速地分析出来。图5为处理过程中检波点集数据的频谱跟踪分析的结果，检波点异常引起的高频干扰比在炮集上表现得更显著。处理中针对检波点集的处理应更为细致，监控分析也应该更加细化。图4(a)、(b)、(c)、(d)为处理过程中共检波点集数据的频谱跟踪分析结果，横坐标是检波点号(SRF_SLOC)，纵坐标是频率(Hz)看出处理过程中数据频率成分的变化。

图5(a)为原始检波点集的频率数据，横坐标是检波点号(SRF_SLOC)，纵坐标是频率(Hz)。图上的曲线是40-60Hz频带内的峰值频率，纵坐标是频率值(Hz)，在峰值频率曲线上可以清楚地看出50Hz工业干扰出现的检波点位置，在2196-2235之间，对这些检波点位置上的数据采用时间域单频波衰减方法压制工业干扰，并通过反褶积展宽数据频带，其结果如5(b)所示。采用本发明中的频带量化分析技术，可以快速反映出不同频带数据的特征，其中存在的问题能够清晰地展现出来，具有较强的针对性。

处理人员结合这些频率图件和量化曲线，采用针对性的处理措施，在保证信噪比的情况下，尽可能提高分辨率，使得地震剖面资料取得了较好的效果，有效频带从10-55Hz提高到8-80Hz(见图6(a)、图6(b)横坐标是CDP号，纵坐标是时间(ms))。

如上述，已经清楚详细地描述了本发明提出的地表一致性海量数据谱估计方法。尽管本发明的优选实施例详细描述并解释了本发明，但是本领域普通的技术人员可以理解，在不背离所附权利要求定义的本发明的精神和范围的情况下，可以在形式和细节中做出多种修改。