CN105334532A - 一种地震子波估计方法 - Google Patents

Abstract

本发明公开了一种地震子波估计方法，所述方法包括以下步骤：获取地震记录并从中选取特定的地震记录作为模型道；基于所述模型道构造地震子波集；基于所述地震子波集求取所述模型道在不同地震子波情况下的稀疏脉冲反射系数；根据所述稀疏脉冲反射系数从所述地震子波集中确定真实地震子波的最佳近似解并作为最终估计的地震子波。与现有技术相比，本发明的方法无需设置太多的初始参数控制，计算量小，方便实用，对窄带地震信号适应性强，具有较强的鲁棒性。

Description

一种地震子波估计方法

技术领域

本发明涉及地质勘探领域，具体说涉及一种地震子波估计方法。

背景技术

在地震勘探领域，地震子波估计的准确性会严重影响子波反褶积或地层参数反演技术的最终应用效果。现有的地震子波估计方法按反射系数是否已知可分为确定性方法和统计性方法两大类。前者需要大量可靠的测井反射系数序列，后者主要针对缺少井资料或井资料无法利用的探区，通过对地层反射系数的分布规律进行假设，然后利用统计类数学方法从偏移后的地震资料中估计出地震子波。受勘探程度和测井反射系数可靠性的制约，统计性地震子波估计方法相对确定性子波估计方法更加实用。

统计性子波估计通常是利用信号分析理论的统计量方法来估计地震子波的振幅谱和相位谱。目前针对缺少井或无井资料探区的地震资料，现有的地震子波估计主要存在以下几个问题：

(1)基于地震记录(绝对)峰度准则这一评价函数的子波相位估计方法不适应于中心频率大于有效频带带宽的窄带地震数据，多解性较强，稳定性较差。

(2)基于高阶统计量的混合相位子波估计法对于含噪数据还难以解决高阶谱相位在2π周期范围内的卷绕问题，难以提取其真实值，很难直接被用于地震子波相位的估算。此外，高阶谱法计算量较大，并且容易受初始值的选取不当而陷入局部解。

(3)目前统计类方法基本没有结合地震反射界面的稀疏分布几何特征和地震稀疏反演的优势来提高地震子波相位估计的可靠性。

因此，针对现有的地震子波估计方法存在的问题，为了在无井约束的情况下从偏移后的地震资料中估计出准确性较高的地震子波，需要一种新的地震子波估计方法。

发明内容

针对现有的地震子波估计方法存在的问题，本发明提供了一种地震子波估计方法，所述方法包括以下步骤：

获取地震记录并从中选取特定的地震记录作为模型道；

基于所述模型道构造地震子波集；

基于所述地震子波集求取所述模型道在不同地震子波情况下的稀疏脉冲反射系数；

根据所述稀疏脉冲反射系数从所述地震子波集中确定真实地震子波的最佳近似解并作为最终估计的地震子波。

在一实施例中，选取多道信噪比以及波形时空一致性满足特定要求的所述地震记录作为所述模型道。

在一实施例中，基于谱模拟和常相位旋转技术构建所述地震子波集。

在一实施例中，在构造所述地震子波集的过程中：

对所述模型道进行谱模拟以求取地震子波的振幅谱；

基于所述振幅谱利用相位扫描法构建所述地震子波集。

在一实施例中，利用所述常相位旋转技术改变所述地震子波的相位谱，形成一系列所述振幅谱相同所述相位谱不同的地震子波序列以构建地震子波集。

在一实施例中，利用基追踪算法求取所述稀疏脉冲反射系数。

在一实施例中，将用于所述谱模拟的所述模型道作为稀疏反射系数反演的模型参考道。

在一实施例中，针对不同的所述地震子波应用非线性稀疏反射系数谱反演技术求取所述模型参考道在L₁范数极小下的稀疏脉冲反射系数集。

在一实施例中，利用稀疏脉冲反射系数的所述L₁范数稀疏准则确定所述最佳近似解。

在一实施例中，在所述稀疏脉冲反射系数集中，最接近所述真实地震子波的所述地震子波反演得到的反射系数序列越稀疏，其对应的所述L₁范数最小。

与现有技术相比，本发明的方法无需设置太多的初始参数控制，计算量小，方便实用，对窄带地震信号适应性强，具有较强的鲁棒性。

本发明的其它特征或优点将在随后的说明书中阐述。并且，本发明的部分特征或优点将通过说明书而变得显而易见，或者通过实施本发明而被了解。本发明的目的和部分优点可通过在说明书、权利要求书以及附图中所特别指出的步骤来实现或获得。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1是根据本发明一实施例的执行流程图；

图2是根据本发明一实施例地震子波与地震合成记录的二次谱；

图3是真实地震子波与根据本发明一实施例估计的地震子波的振幅谱；

图4是根据本发明一实施例相同振幅谱不同相位谱的子波集示意图；

图5a-5d为根据本发明一实施例稀疏反射系数模型的合成地震记录在不同子波下的反演结果及其L₁范数；

图6a为一层状稀疏地层二维模型的合成地震记录；

图6b为根据本发明一实施的估计子波的振幅谱和真实子波的对比图；

图6c为根据本发明一实施的最终估计的子波波形。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此本发明的实施人员可以充分理解本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程并依据上述实现过程具体实施本发明。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

综上，目前针对缺少井或无井资料探区的地震资料，现有的地震子波估计主要存在以下几个问题：

(1)基于地震记录(绝对)峰度准则这一评价函数的子波相位估计方法不适应于中心频率大于有效频带带宽的窄带地震数据，多解性较强，稳定性较差。

(2)基于高阶统计量的混合相位子波估计法对于含噪数据还难以解决高阶谱相位在2π周期范围内的卷绕问题，难以提取其真实值，很难直接被用于地震子波相位的估算。此外，高阶谱法计算量较大，并且容易受初始值的选取不当而陷入局部解。

(3)目前统计类方法基本没有结合地震反射界面的稀疏分布几何特征和地震稀疏反演的优势来提高地震子波相位估计的可靠性。

因此，针对现有的地震子波估计方法存在的问题，为了在无井约束的情况下从偏移后的地震资料中估计出准确性较高的地震子波，服务于子波反褶积或地层参数反演，识别地下地层或岩性结构，本发明提出了一种地震子波估计方法。

近年来，一些研究结果表明，相同振幅谱不同相位的子波的非线性稀疏反射系数谱反演结果的稀疏性存在明显差异，准确子波的非线性稀疏反射系数谱反演结果的稀疏性最好。

地下岩性地层分界面(地震反射界面)在地震地层学意义下是稀疏层状分布的。这意味着可以利用非线性稀疏反射系数谱反演结果的稀疏性作为评价准则来评价子波是否准确。在数学上，通常利用L₁范数作为向量序列稀疏性的评价准则。此外，研究成果表明非线性稀疏反射系数谱反演对窄带地震信号的适应性较强，反演结果的稳定性和抗噪性较好。因此，利用L₁范数稀疏准则作为子波相位估计是否准确的判断准则，可以增强常相位扫描法对窄带地震信号的适用性，稳定性和抗噪性。

联合常相位扫描法和L₁范数稀疏准则的优点，在基于低通滤波的谱模拟法和非线性稀疏反射系数谱反演法的基础上，发明了一种基于L₁范数稀疏准则的地震子波估计方法。本发明无需井资料和井震标定等作约束，基于地震资料的频谱信息，利用谱模拟和相位旋转技术建立地震子波集，在此基础上，利用稀疏反射系数谱反演技术求取不同子波的稀疏反射系数集，利用稀疏层状地层的数学表征函数(L₁范数)作为子波相位估计的评价准则函数，选取L₁范数最小的稀疏反射系数对应的子波为最终估计的子波。

本发明与传统的基于峰度准则的常相位扫描地震子波估计方法相比，其创新之处在于，抛弃了传统峰度准则直接应用于地震信号的策略，在地震地层学意义稀疏层状地层假设下，通过引入L₁范数稀疏准则来评价相同振幅谱不同相位的地震子波的稀疏反射系数谱反演结果的稀疏性来确定地震子波的相位，将非线性稀疏反射系数谱反演和子波相位估计有机地结合起来，可以增强常相位扫描法对窄带地震信号的适用性和稳定性。

接下来基于附图来详细描述本发明一实施例的执行过程，附图的流程图中出示的步骤可以在包含诸如一组计算机可执行指令的计算机系统中执行。虽然在流程图中示出了各步骤的逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

如图1所示，执行本发明的方法首先执行步骤S110，选取模型道步骤，获取地震记录并从中选取特定的地震记录作为模型道。在步骤S110中，从地震记录中选取多道信噪比以及波形时空一致性满足特定要求的地震记录作为模型道。具体的，对比地震记录的信噪比以及时空一致性，从中选出信噪比以及时空一致性相对较好的地震记录。

接下来需要执行步骤S120，构造地震子波集步骤，基于模型道构造地震子波集。本实施例基于谱模拟和常相位旋转技术构建地震子波集。具体的，在步骤S120中，首先对模型道进行谱模拟求取地震子波的振幅谱。

本实施例采用基于低通滤波的谱模拟法来估计地震子波的振幅谱。该方法计算过程稳定、简单、效率高，估计的振幅谱准确性较高。其基本原理是：在地震子波振幅谱光滑的假设条件下，利用地震信号二次振幅谱能够分离高低频信息的特点，设计低通滤波器将分布在地震信号二次振幅谱低频位置的子波二次谱估计出来，进而通过一次傅立叶逆变换求得子波的振幅谱。

基于上述原理，在本实施例中，首先执行步骤S121，求取地震记录的振幅谱。然后执行步骤S122，求取地震记录的二次振幅谱。具体的，对地震记录连续做两次傅立叶变换可以得到地震记录的二次谱，然后令地震记录的二次谱

|S(t)|⁽²⁾＝FT(|S(f)|)＝|W(t)|⁽²⁾×|R(t)|⁽²⁾(1)。

执行步骤S123，设计低通滤波器b(t)。接着就可以执行步骤S124，求取地震子波的二次振幅谱|W(t)|⁽²⁾。并进一步执行步骤S125，求取地震子波的振幅谱|W(f)|。步骤S124以及步骤S125的计算公式如下：

$\left\{\begin{array}{c}{\left|W\left(t\right)\right|}^{\left(2\right)}\≈{\left|S\left(t\right)\right|}^{\left(2\right)}\×b\left(t\right)\\ \left|W\left(f\right)\right|={\mathrm{FT}}^{-1}\[{\left|W\left(t\right)\right|}^{\left(2\right)}\]\end{array}\right.---\left(2\right).$

其中，FT和FT^-1为正反傅立叶变换，地震记录s(t)的振幅谱|S(f)|＝|FT[s(t)]|，低通滤波器b(t)的表达式如下：

$b\left(t\right)=\left\{\begin{array}{cc}1.0& 0\≤\left|t\right|\≤\frac{T}{2}\\ \frac{1}{2}(1+cos(\mathrm{\π}\frac{t-\mathrm{\λ}\frac{T}{2}}{2(1-\mathrm{\λ})\frac{T}{2}}\left)\right)& \mathrm{\λ}\frac{T}{2}\≤\left|t\right|\≤\frac{T}{2}\end{array}\right.---\left(3\right)$

其中，λ∈(0,1]，通常取1；T为低通滤波的高截位置。低通滤波器高截位置同地震子波有效延展长度和地震记录二次谱高低频信息的分离情况有关。高截位置可以通过地震记录的谱估计出地震子波的主频，进而估计出有效延展长度，并结合地震记录二次谱高低频信息的分离情况来确定。

在本发明的一具体应用实例中，如图2所示，横坐标为时间(单位毫秒(ms))，纵坐标为二次振幅谱，虚线代表地震记录，实线代表地震子波。图2中地震子波的主频为25Hz，有效延展长度为79ms，二次谱高低频信息在50ms附近分离开。因此，根据公式3可设计低通滤波器的高截位置在T/2＝40ms附近是合理的。考虑二次谱的对称性，另外半只低通滤波器的高截位置为-T/2＝-40ms，将该低通滤波器应用于图2中所示的例子，即可估计出地震子波的二次谱。

应用傅立叶逆变换处理地震子波的二次谱，可估计出地震子波的振幅谱。如图3所示，横坐标为频率(单位赫兹(Hz))，纵坐标为振幅谱，虚线代表子波振幅谱的估计值，实线代表子波振幅谱的真实值。从图3中可以看出，估计的子波振幅谱同真实子波的振幅谱吻合程度较高，这说明上述基于低通滤波的谱模拟法是有效的。

接下来就可以基于地震子波振幅谱利用相位扫描法构建地震子波集。若地震子波的振幅谱已知，假设地震子波的相位谱与频率无关，则可利用常相位旋转技术改变地震子波的相位谱，形成一系列振幅谱相同相位谱不同的地震子波序列以构建地震子波集。

具体的，首先要执行步骤S126，确定地震子波相位扫描间隔。设子波集的相位范围在若忽略子波的极性，可只取一、四象限作为相位的范围来构建子波集，即令子波相位扫描间隔为则子波集的相位扫描序列为运算符[·]_n代表向下取整。

接下来执行步骤S127，形成地震子波集步骤。在频率域，常相位旋转技术可通过式(4)实现：

其中，i表示虚数单位。

在时间域中，常相位旋转技术可通过式(5)实现：

其中，H[·]表示希尔伯特(Hilbert)变换。w₀(t)为振幅谱|W(f)|对应的零相位子波，w₀(t)＝FT^-1[|W(f)|]。

由此，即可通过式(4)或式(5)获得相同振幅谱不同相位谱的子波集W＝{w_j(t)}。

以图3中所示的振幅谱为例，以相位扫描间隔构建的地震子波集如图4所示(横坐标为相位(单位度(°))，纵坐标为时间(单位毫秒(ms)))。从图4中可以看出这些子波由于具有相同的振幅谱，因此有效延展长度基本一致。

由于相位的不同，在波形形态上存在明显的差异。为了使构建的子波集尽可能的逼近真实子波，就需要加密相位扫描间隔，涵盖更多不同相位的地震子波，通常子波相位扫描间隔的取值范围为(1°～5°)π/180°。常相位旋转技术在频率域和时间域均可实现。

根据步骤S120可建立起包含真实地震子波近似解的地震子波集。接下来，就可以执行步骤S130，求取系数脉冲反射系数步骤。基于地震子波集求取模型道在不同地震子波情况下的稀疏脉冲反射系数。在本实施例中，利用基追踪算法求取稀疏脉冲反射系数。具体的，首先将用于谱模拟的模型道作为稀疏反射系数反演的模型参考道。然后针对不同的地震子波应用非线性稀疏反射系数谱反演技术求取模型参考道在L₁范数极小下的稀疏脉冲反射系数集。

由地震褶积模型和反演理论可知，L₁范数极小下的非线性稀疏反射系数谱反演的目标函数如下：

min‖S(f)-W(f)G(f,t)r(t)‖₂+μ‖r(t)‖₁(6)

其中，f∈[f_low,f_high]，f_low和f_high分别是模型参考地震道频谱S(f)有效频带的低、高截止频率；W(f)为子波的频谱，G(f,t)＝exp(-i2πft)为子波延迟传递函数，r(t)为稀疏地层反射系数，μ为正则化因子。在本实施例中，式(6)采用基追踪算法稳健地求出全局最优解，得到不同子波w_j(t)的稀疏脉冲反射系数r_j(t)。

接下来执行步骤S140，确定最佳近似解，根据稀疏脉冲反射系数从地震子波集中确定真实地震子波的最佳近似解。在本实施例中，利用稀疏脉冲反射系数L₁范数稀疏准则确定最佳近似解。具体的，准确(最接近真实地震子波)的地震子波反演得到的反射系数序列越稀疏，其对应的L₁范数最小。

根据准确子波对应的稀疏脉冲反射系数L₁范数最小的原则，可以利用公式

$w_e\left(t\right)=\arg \underset{w_j\∈W}{min}\left|\right|r_j\left(t\right)|{|}_1---\left(7\right)$

估计出地震子波。

以一具体应用实例为例，图5a-5d为稀疏反射系数模型的合成地震记录在不同子波下的反演结果及其L₁范数。从图5a-5d中可以看出，准确子波反演得到的反射系数序列越稀疏，其L₁范数最小。显然，采用上述稀疏准则，可以从子波集中找到真实的子波或其最佳近似解。

对于如图5a所示的某一稀疏地层模型的合成地震记录(纵坐标为时间(单位毫秒(ms)，横向上分别为稀疏反射系数以及合成地震记录波形图)，可以构建出如图5b所示的一系列相同振幅谱不同相位谱的子波集(横坐标为相位(单位度(°))，纵坐标为时间(单位毫秒(ms)))。若振幅谱估计相对准确，相位离散间隔足够，那么图5b中必然存在一个子波和真实的子波非常逼近。

图5b中灰色子波是用于合成图5a所示的地震记录的真实子波，相位为30°。对于图5b所示的每个相同振幅谱不同相位的子波，应用公式(6)，对图5a所示的合成地震记录进行反演，可以得到如图5c所示的每个子波对应的稀疏反射系数(横坐标为相位(单位度(°))，纵坐标为时间(单位毫秒(ms))。从图5c中可以看出，不同相位的子波的反射系数反演结果的稀疏性存在明显的差别，这就是应用稀疏准则作为相位估计评判函数的基础。

利用公式‖r_j(t)‖₁即可求出如图5d所示的L₁范数稀疏准则的函数曲线(横坐标为相位(单位度(°))，纵坐标为L₁范数)。从图5d中不难看出，L₁范数最小值位于相位30°的位置，和图5b中所示的真实子波位置一致，这说明准确子波(最接近真实子波)的稀疏反射系数的反演结果的稀疏性最佳。应用公式(7)，即可从子波集中筛选出最终估计的子波。此外，L₁范数稀疏准则函数曲线的最小值位置的波谷呈“深洼状”，从数学最优化理论上讲，本发明估计子波相位的稳定性较好，多解性较弱，具有较佳的鲁棒性。

最后执行步骤S150，将上述最佳近似解作为最终估计的地震子波并输出。

为测试基于L₁范数稀疏准则的地震子波估计方法的应用效果。现以一层状稀疏地层二维模型为例进行试验，其合成地震记录如图6a所示(横坐标为道号，纵坐标为时间(单位毫秒(ms))。图6b为估计子波的振幅谱和真实子波的对比图(横坐标为频率(单位赫兹(Hz))，纵坐标为振幅谱)，可以看出，和真实子波的振幅谱非常接近。图6c为最终估计的子波波形(横坐标为时间(单位毫秒(ms)，纵坐标为振幅)，从图6c中可以看出，估计子波的波形和真实子波的波形吻合较好。试验结果表明，研制的基于L₁范数稀疏准则的地震子波估计方法能够较准确地从地震记录中估计出子波反褶积或反演所需的子波。

本发明本质上是一种基于准则函数判别的相位扫描法地震子波估计方法，同传统的峰度准则不同的是，采用了非线性稀疏反射系数谱反演结果的稀疏性这一地震地层学假设作为地震子波相位估计的准则，通过稀疏反射系数的L₁范数极小来判断子波相位是否估计准确。由于非线性稀疏反射系数谱反演方法对窄带信号的适应性和抗噪性较强，使得L₁范数稀疏准则相比传统的峰度准则具有更强的鲁棒性。

与现有技术相比，本发明的方法计算简单，运用方便，估计的子波准确性较好，对窄带地震信号的适应性较强，能够为子波反褶积或模型参数反演等地震资料解释性处理技术的应用提供较为准确的子波。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。本发明所述的方法还可有其他多种实施例。在不背离本发明实质的情况下，熟悉本领域的技术人员当可根据本发明作出各种相应的改变或变形，但这些相应的改变或变形都应属于本发明的权利要求的保护范围。

Claims (10)

1.一种地震子波估计方法，其特征在于，所述方法包括以下步骤：

获取地震记录并从中选取特定的地震记录作为模型道；

基于所述模型道构造地震子波集；

基于所述地震子波集求取所述模型道在不同地震子波情况下的稀疏脉冲反射系数；

根据所述稀疏脉冲反射系数从所述地震子波集中确定真实地震子波的最佳近似解并作为最终估计的地震子波。

2.根据权利要求1所述的方法，其特征在于，选取多道信噪比以及波形时空一致性满足特定要求的所述地震记录作为所述模型道。

3.根据权利要求1所述的方法，其特征在于，基于谱模拟和常相位旋转技术构建所述地震子波集。

4.根据权利要求2所述的方法，其特征在于，在构造所述地震子波集的过程中：

对所述模型道进行谱模拟以求取地震子波的振幅谱；

基于所述振幅谱利用相位扫描法构建所述地震子波集。

5.根据权利要求4所述的方法，其特征在于，利用所述常相位旋转技术改变所述地震子波的相位谱，形成一系列所述振幅谱相同所述相位谱不同的地震子波序列以构建地震子波集。

6.根据权利要求5所述的方法，其特征在于，利用基追踪算法求取所述稀疏脉冲反射系数。

7.根据权利要求6所述的方法，其特征在于，将用于所述谱模拟的所述模型道作为稀疏反射系数反演的模型参考道。

8.根据权利要求7所述的方法，其特征在于，针对不同的所述地震子波应用非线性稀疏反射系数谱反演技术求取所述模型参考道在L₁范数极小下的稀疏脉冲反射系数集。

9.根据权利要求8所述的方法，其特征在于，利用稀疏脉冲反射系数的所述L₁范数稀疏准则确定所述最佳近似解。

10.根据权利要求9所述的方法，其特征在于，在所述稀疏脉冲反射系数集中，最接近所述真实地震子波的所述地震子波反演得到的反射系数序列越稀疏，其对应的所述L₁范数最小。