地震数据空变广义子波的提取方法与流程

本发明涉及一种地震数据子波的提取方法，特别是涉及一种地震数据空变广义子波的提取方法。

背景技术：

地震子波应当反映地震波在地层传播过程中受地层粘弹性介质的吸收衰减作用所导致的波形畸变等效应。地层粘弹性介质对地震信号具有滤波作用，地震信号在地下介质传播过程之中能量逐渐被地层介质吸收而衰减，而且信号的不同频率分量其吸收衰减程度也不同。尤为重要的是地震信号在能量衰减的同时还产生了相位畸变，因为信号的不同频率分量其传播速度也不相同，即存在频散现象或相位延迟。因此，迫切需要从地震数据中提取能够真实反映上述地层滤波效应的空变地震子波。

现有的基于零相位和常相位的子波估算方法所提取的地震子波不能准确表征上述波形畸变等效应。零相位子波提取时对地震数据的功率谱进行平滑处理得到子波的功率谱，并设子波的相位为零，得到时间轴上对称的零相位子波。但是，实际地震信号受频散作用的影响，从来就不是零相位子波。常数相位子波提取的有效方法有峰态扫描法。在统计学中，峰态可以估量地震序列偏离高斯正态分布的程度，因此根据最大峰态值可以确定地震子波的相位。但是，常相位子波仍不能真实反映地震信号在地层传播过程中的吸收和频散作用。

在没有测井记录作为反射系数模型的情况下，混合相位子波的提取需要利用高阶统计分析方法。地震数据的高阶统计量可以精确反映相位信息，因此利用高阶统计量包括三阶累积量和四阶累积量，可以提高地震子波估算的精度（mixed-phasewaveletestimationbyiterativelinearinversionofhigh-orderstatistics,journalofgeophysicsandengineering,vol:4/2007/本发明人;cn102096101a/2010;cn103645500a/2013）。但是，由于需要计算高阶统计量，并需要对高阶统计量进行数值分析运算，该子波估算方法计算量巨大，而且稳定性较差，严重限制了该方法在实际地震数据分析之中的应用。因此，本发明提出一种从实际地震数据中提取空变地震子波的方法，真正反映地震信号在粘弹性空间的变化特征，兼具运算过程的抗噪能力和稳定性。

技术实现要素：

本发明目的在于克服现有技术的上述缺陷，提供一种地震数据空变广义子波的提取方法，所提取的地震子波能够真正反映地震信号的能量衰减和相位畸变及其在空间上的变化特征。

本发明提出地震数据空变广义子波的提取方法，其特别之处在于步骤有：（1）根据地震数据功率谱估算空变子波的初始参数；（2）确定全局化的广义子波参考频率；（3）估算广义子波的空变分数阶值；（4）优化全局化的广义子波参考频率；（5）优化广义子波的空变分数阶值；（6）构建空变地震广义子波。

本发明提出从实际地震数据中提取空变地震子波的方法。本发明人2015年在国际上首次提出“广义地震子波”的理论（generalizedseismicwavelets,geophysicaljournalinternational,vol:203/2015/本发明人），下简称wang’s广义子波理论，其特点是仅需两个参数即可定义非对称子波。因此，本发明在此提出采纳上述广义子波理论表征实际地震数据中的空变子波。本发明的第一有益效果在于其抗噪能力，兼运算过程的稳定性。本发明的第二有益效果则在于地震子波的空变特性，所提取的子波能够真正反映地震信号的能量衰减和相位畸变及其在空间上的变化特征。

作为优化，步骤（1）根据地震数据功率谱估算空变子波的初始参数，分为三步实现；首先，在地震数据体的不同位置，选取时间-空间窗口，计算地震数据的功率谱；其次，设定地震数据功率谱的频段代表子波的功率谱，计算其功率谱的基本统计特性，包括均值频率和标准方差；然后，根据每个窗口功率谱的基本统计特性确定对应的广义子波分数阶值和参考频率。

作为优化，步骤（2）中确定适应全局的广义地震子波的参考频率；依据wang’s广义子波参考频率的物理意义即为地震震源的固有频率，对上述各窗口的广义子波参考频率进行中值滤波，取得能够代表全局的参考频率（ƒ0）。

作为优化，步骤（3）中依据全局化的参考频率以及各窗口的均值频率和标准方差，估算与各窗口对应的分数阶值。

作为优化，广义子波分数阶值（u）的计算方式：

该计算方式同步运用均值频率（ƒm）和标准方差（ƒσ），降低在计算分数阶值时可能的计算误差。

作为优化，步骤（4）中通过地震频谱的匹配实现对全局广义子波参考频率的优化，提出匹配系数如下：

式中对实际地震频谱wobs(ƒ)和地震子波频谱w(ƒ)分别做了归一化处理，点乘（·）的结果则为互相关系数；式中采用n个互相关系数的平均值表征匹配系数。

作为优化，步骤（4）中地震频谱匹配的优化式如下：

式中k表示的是窗口号，ƒ0^ref是参考频率初始估算值，作为参数优化问题的模型约束，其约束系数设为μ；采用参考频率初始估算值作为模型约束，增强优化过程的抗噪声能力。

作为优化，步骤（5）中通过地震频谱的匹配实现空变广义子波分数阶值的优化，提出优化式如下：

式中u为分数阶值空变函数u(t,x)，纵向上随时间t（或者深度z）而变化，横向上随距离x而变化，u^ref则是初始估算值。上述分数阶值的优化过程之中，在运用初始估算对模型约束的同时，对模型在时间和空间上的变化形态亦予以约束，三项约束的系数分别设为（μ1，μ2，μ3）。

作为优化，步骤（6）分为两步实现；首先，计算每个窗口的广义子波频谱；然后，通过数值求解的计算方式对频谱做傅立叶逆变换得到时间域的广义子波。

本发明提出的地震数据空变广义子波的提取方法，采用上述技术方案之后，具有如下优点：提取方法运算稳定，具有非常强的抗噪能力，提取的地震广义子波形态简洁，所提取的子波能够真实反映地震信号在地层介质传播过程之中的吸收衰减等空变特征。

附图说明

图1是本发明地震数据空变广义子波提取方法的流程框图；图2是本发明地震数据空变广义子波提取方法在不同深度提取的地震广义子波（虚线）与实测地震子波（实线）之间的波形对比图；图3是本发明地震数据空变广义子波提取方法实施例的地震数据剖面，左侧选取上下a、b两个空间-时间窗口，右侧选取上下c、d两个空间-时间窗口；图4是本发明地震数据空变广义子波提取方法的实际地震数据频谱以及提取的广义地震子波波形图。

具体实施方式

本发明地震数据空变广义子波的提取方法的步骤有：（1）根据地震数据功率谱估算空变子波的初始参数；（2）确定全局化的广义子波参考频率；（3）估算广义子波的空变分数阶值；（4）优化全局化的广义子波参考频率；（5）优化广义子波的空变分数阶值；（6）构建空变地震广义子波。本发明在此提出采纳wang’s广义子波理论表征从实际地震数据中提取的空变子波，因为该理论的特点是仅需两个参数即可定义非对称的广义地震子波（generalizedseismicwavelets.geophysicaljournalinternational,vol:203/2015/本发明人）。本发明第一有益效果在于其空变特性，所提取的地震子波能够真正反映地震信号能量衰减和相位畸变及其在空间上的变化特征。本发明第二有益效果则在于其抗噪能力，兼运算过程的稳定性。

步骤（1）根据地震数据功率谱估算空变子波的初始参数，分为三步实现；首先，在地震数据体的不同位置，选取时间-空间窗口，计算地震数据的功率谱；其次，设定地震数据功率谱的频段代表子波的功率谱，计算其功率谱的基本统计特性，包括均值频率和标准方差；然后，根据每个窗口功率谱的基本统计特性确定对应的广义子波分数阶值和参考频率。

步骤（2）中确定适应全局的广义地震子波的参考频率；依据wang’s广义子波参考频率的物理意义即为地震震源的固有频率，对上述各窗口的广义子波参考频率进行中值滤波，取得能够代表全局的参考频率（ƒ0）。

步骤（3）中依据全局化的参考频率以及各窗口的均值频率和标准方差，估算与各窗口对应的分数阶值。广义子波分数阶值（u）的计算方式：

该计算方式同步运用均值频率（ƒm）和标准方差（ƒσ），降低在计算分数阶值时可能的计算误差。

步骤（4）中通过地震频谱的匹配实现对全局广义子波参考频率的优化，提出匹配系数如下：

式中对实际地震频谱w0bs(ƒ)和地震子波频谱w(ƒ)分别做了归一化处理，点乘（·）的结果则为互相关系数；式中采用n个互相关系数的平均值表征匹配系数。步骤（4）中地震频谱匹配的优化式如下：

式中k表示的是窗口号，ƒ0^ref是参考频率初始估算值，作为参数优化问题的模型约束，其约束系数设为μ；采用参考频率初始估算值作为模型约束，增强优化过程的抗噪声能力。

步骤（5）中通过地震频谱的匹配实现空变广义子波分数阶值的优化，提出优化式如下：

式中u为分数阶值空变函数u(t,x)，纵向上随时间t（或者深度z）而变化，横向上随距离x而变化，u^ref则是初始估算值。上述分数阶值的优化过程之中，在运用初始估算对模型约束的同时，对模型在时间和空间上的变化形态亦予以约束，三项约束的系数分别设为（μ1，μ2，μ3）。

步骤（6）分为两步实现；首先，计算每个窗口的广义子波频谱；然后，通过数值求解的计算方式对频谱做傅立叶逆变换得到时间域的广义子波。

如图所示，为使本发明地震数据空变广义子波的提取方法和优点更加清楚明白，下面结合附图对本发明实施方式作详细具体说明。本发明提供了基于实际地震数据功率谱提取空变广义地震子波的方法，图1是该方法的基本流程图。

步骤s101，根据地震数据功率谱估算空变子波的初始参数，分为三步实现。首先，在地震数据体的不同位置选取时间-空间窗口，每个窗口中的地震道构建成一个超长序列，计算该超长序列的功率谱；其次，选定起始频率和截止频率，保证其间的地震数据功率谱能够代表子波的功率谱，计算所截频段功率谱的基本统计特性，包括均值频率和标准方差；然后，根据地震功率谱的基本统计特性确定对应每个窗口的广义子波分数阶值和参考频率。

步骤s102，确定全局化的广义子波参考频率。本发明指出wang’s广义子波参考频率的物理意义是地震震源的固有频率，同一个地震数据体中地震震源一致，其固有频率也一致。基于上述物理意义对各窗口的参考频率进行中值滤波，取得代表全局的广义子波参考频率（ƒ0）。

步骤s103，估算广义子波的空变分数阶值（u），本发明提出如下计算方式：

计算过程中同步运用均值频率（ƒm）和标准方差（ƒσ）计算分数阶值，而非仅仅运用两者之一，用于降低这些基本统计特性的估算误差可能导致的分数阶值计算潜在误差。

步骤s104，通过地震频谱的匹配实现对全局化的广义子波参考频率的优化。本发明提出匹配系数如下：

式中对实际地震数据频谱wobs(ƒ)和地震子波频谱w(ƒ)进行了归一化处理，向量点乘（·）则为互相关系数，式中采用n个互相关系数的平均值表征匹配系数。本发明提出地震频谱的匹配优化方式如下：

式中k是窗口编号，参考频率的初始估算值ƒ0^ref作为参数优化问题的模型约束，其约束系数设为μ，目的是增强参数优化过程的抗噪能力及其稳定性。

步骤s105，优化广义子波的空变分数阶值。通过地震频谱的匹配实现空变广义子波分数阶值的优化如下：

式中u为分数阶值空变函数u(t,x)，纵向上随时间t（或者深度z）而变化，横向上随距离x而变化，u^ref则是初始估算值。上述分数阶值的优化过程之中，在运用初始估算对模型约束的同时，对模型在时间和空间上的变化形态亦予以约束，三项约束的系数分别设为（μ1，μ2，μ3）。

步骤s106，构建地震空变广义子波，分为两步实现。首先，计算每个窗口的广义子波频谱；然后，通过采用数值求解的计算方式对子波频谱做傅立叶逆变换得到时间域的广义子波。

本发明还提供两个实施例，说明本发明的有益效果。实施例图2对比了依据本发明方法提取的广义地震子波（虚线）与实测地震子波（实线）的波形。地震子波随着地层深度z而变化，实测地震的深度位置从600米至1700米。本实施例提取的全局化的广义地震子波参考频率f0为32.4赫兹，分数阶值随着深度增加而逐渐变小。所以，如图2所示，本发明的有益效果在于，所提取的广义地震子波能够真实反映地震信号在地层介质传播过程中的吸收衰减等空变特征。

实施例图3-图4演示本发明方法从实际地震数据中提取空变广义子波的效果。图3展示一个实际地震数据剖面。为了演示本发明方法的效果仅从地震数据剖面的左侧选取上下a、b两个窗口，从其右侧选取上下c、d两个窗口。地震数据剖面的时间段在3.3-5.0秒，选取如此深部的地震可以数据检验本发明方法的抗噪能力及其稳定性。

图4是本发明地震数据空变广义子波提取方法中实际地震数据的频谱以及相对应的广义地震子波的波形图。如图4所示，本发明地震数据空变广义子波的提取方法具有地震子波形态简洁的优点，广义地震子波的参考频率为19.5hz。就实施例图4中选取的两列对比而言，广义子波的分数阶值横向上变化不很明显，但是在纵向上有着非常显著的变化。所以，如图4所示，本发明有益效果还在于，空变广义子波提取方法具有稳定性好、抗噪能力强、以及广义地震子波空变特征提取精确等优点。