基于广义S变换的谱相关系数的Q值估计方法与流程

本发明涉及地震信号处理领域，特别涉及一种基于广义s变换的谱相关系数的q值估计方法。

背景技术：

目前，在地震信号处理领域，国内外已有多种q(地震波衰减品质因子)值估计方法，主要包括时间域估计方法、频率域估计方法和反演类q值估计方法。其中，频率域q值估计方法是应用比较多的一种方法，其在q值估计的过程中需要提取振幅谱，而振幅谱的提取可以采用不同的时频分析方法。在近几年，发展出多种时频分析方法，可用于提取振幅谱，主要包括gabor变换、小波变换，s变换和广义s变换等。

现有技术把各种时频分析方法与谱比法等频率域q值估计方法相结合来进行q值估计，并且不断地进行了改进和创新。目前常用的q值估算方法如对数谱比法、峰值频率法等在进行q值估计时都受到加时窗问题的影响，峰值频率法和解析信号法对噪声很敏感。多数改进都是在峰值频率法、谱比法和质心频移法的基础上，结合多种时频分析方法，但仍然不能突破每种方法的局限性，并且会受到所选择的时频分析方法的影响。

在应用时频分析方法时，选用可变时窗的时频变换能够得到更为精确的结果。s变换的窗函数以固定趋势随频率变化，不能根据实际需要进行调整，因此其应用受到了一定限制。因此研究一种能避免频带范围选择的问题而且具有一定的抗噪性，同时能减小时频分析方法的影响的q值估计方法，具有重要意义。

技术实现要素：

本发明的目的是提供一种基于广义s变换的谱相关系数的q值估计方法，其能够避免频带范围选择的问题而且具有一定的抗噪性，同时能减小时频分析方法的影响的q值估计方法。

本发明采用以下解决方案：

一种基于广义s变换的谱相关系数的q值估计方法，包括以下步骤：

步骤1：输入地震数据；

步骤2：对所述地震数据进行广义s变换，得到变换后的地震数据；

步骤3：获得变换后的地震数据的估计振幅谱和实际振幅谱；

步骤4：基于所述估计振幅谱和实际振幅谱，构建谱相关系数；

步骤5：求解所述谱相关系数最小时的q值。

优选地，所述地震数据为u(ω,z)，其中ω是地震波的角频率，z是地震波的传播距离。

优选地，根据以下公式(8)得到变换后的地震数据sn(τ,f)：

其中，f表示频率，t表示地震波从地面开始传播到待分析地层的旅行时间，τ表示所述窗函数w(t-τ,f)的时间位置，λ表示用于调节所述窗函数的时间延续长度的参数，p表示用于调节所述窗函数的衰减趋势的参数，其中所述窗函数w(t-τ,f)表示为以下公式(1)：

优选地，根据以下公式(11)计算所述估计振幅谱

其中，h(f,q,τ)表示衰减传播算子，由以下公式(10)表示：

其中t0和t1分别表示地震波传播到待分析地层的上界面和下界面所用的双程旅行时，as(t0,f)和as(t1,f)分别表示t0和t1时的振幅谱。

优选地，根据以下公式(12)构建所述谱相关系数：

其中，f(q)表示所述谱相关系数，〈·，·〉表示内积算子，║·║表示l²范数。

优选地，针对所述谱相关系数f(q)，在给定范围内改变q值，获得所述谱相关系数f(q)最小时的q值。

与现有技术相比，本发明的有益效果是：

(1)在应用谱比法估计q值时，需要选择一个频带范围，这个频带范围会影响q值估计的精度。应用本发明的方法估计q值时，可以选择全频带进行q值估计，因此不受频带范围选择的影响；

(2)在实际应用中，谱相关系数法本身具有抗噪性，其抗噪性好于谱比法和峰值频率法，精度也较高。

附图说明

通过结合附图对本公开示例性实施例进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显。

图1示出了根据示例性实施例的基于广义s变换的谱相关系数的q值估计方法的流程图；

图2a和图2b分别示出了根据示例性实施例的合成地震记录及其广义s变换；

图3a和图3b分别示出了示例性实施例中第二层上下界面的广义s变换振幅谱及其拟合后的广义s变换振幅谱；

图4示出了示例性实施例中目标函数f随q值变化的关系曲线；

图5示出了示例性实施例中的实际叠后地震数据；

图6示出了示例性实施例中的实际叠后地震数据的广义s变换结果；以及

图7示出了示例性实施例中估计的实际数据的单道q值。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施例。虽然附图中显示了本公开的优选实施例，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

近几年，随着时频分析方法的发展，将时频分析方法与谱比法等频率域q值估计方法相结合，可以估计地层q值。应用谱相关系数法求振幅谱时可以借助于gabor变换、小波变换、s变换和广义s变换等时频分析方法。reine等人在2009年比较了stft和可变时窗的时频变换，如s变换、连续小波变换等，指出选用可变时窗的时频变换能够得到更为精确的结果。s变换的窗函数以固定趋势随频率变化，不能根据实际需要进行调整，因此其应用受到了一定限制，对s变换进行了改进提出几种不同窗函数的广义s变换。本发明的示例性实施例将广义s变换运用到地震数据中估计q值，先利用广义s变换求得振幅谱，再求得振幅谱相关系数，最后运用迭代法估算q值。

图1示出了根据示例性实施例的基于广义s变换的谱相关系数的q值估计方法的流程图，以下参考图1描述根据示例性实施例的基于广义s变换的谱相关系数的q值估计方法，其包括以下步骤：

步骤1：输入地震数据

地震数据表示为u(ω,z)，其中ω是地震波的角频率，z是地震波的传播距离。

步骤2：对地震数据进行广义s变换，得到变换后的地震数据sn(τ,f)

pinnegar等人提出了既可以调节窗函数标准差、窗口又不对称的广义s变换。广义s变换的窗函数为双曲窗函数，由以下公式(1)表述：

其中，f表示频率，t表示地震波从地面开始传播到待分析地层的旅行时间，τ表示所述窗函数w(t-τ,f)的时间位置，λ表示用于调节所述窗函数的时间延续长度的参数，p表示用于调节所述窗函数的衰减趋势的参数。

基于以上窗函数，信号h(t)的广义s变换表达式如以下公式(2)所示：

其中，窗函数w(t-τ,f)随参数λ和p的改变而变化，因此可以根据实际应用需要合理选择参数λ和p，调节广义s变换的时频分辨率。

在常q介质模型中，认为地震波为平面波，地震波u(ω,z)在一维粘弹性介质中沿z轴的增加方向以角频率ω传播时，传播方程可表示为以下公式(3)：

其中，ω是地震波的角频率，z是地震波的传播距离，c(ω)是地震波的相速度。

假定震源为理想脉冲源，即|u(ω,0)|＝1，并且忽略由于散射引起的地震波衰减，对公式(3)进行简化，可以得到以下公式(4)

公式(4)中的相速度c(ω)可以表示为公式(5)：

其中，c0为角频率ω0处的相速度，地震波传播时间t＝z/c(ω)，则公式(4)可改写为以下公式(6)：

对公式(6)进行广义s变换，可以得到以下公式(7)：

已知则公式(7)可改写成以下公式(8)，通过以下公式(8)即可得到变换后的地震数据sn(τ,f)：

步骤3：获得变换后的地震数据sn(τ,f)的估计振幅谱和实际振幅谱as(t1,f)

根据步骤2得到的变换后的地震数据sn(τ,f)，可以获得其估计振幅谱详细过程如下所述，而实际振幅谱as(t1,f)是通过实际测量获得的。

对于步骤2得到的变换后的地震数据sn(τ,f)，其振幅谱可以表示为以下公式(9)：

其中，t表示地震波从地面开始传播到待分析地层的旅行时间。

对于待分析的地层，假设地震波传播到该层上下界面所用的双程旅行时分别为t0和t1，振幅谱分别为as(t0,f)和as(t1,f)。根据谱相关系数法，由于只考虑振幅谱，衰减传播算子h(f,q,τ)可以表示为以下公式(10)：

衰减传播算子h(f,q,τ)的频率取值范围与所考虑的地震信号频率范围相同。为了建立振幅谱互相关系数与地震波衰减品质因子q的关系，将衰减传播 算子h(f,q,τ)作用到t0时信号的振幅谱上就得到了对t1时信号的振幅谱的一个估计如以下公式(11)所示：

步骤4：基于估计振幅谱和实际振幅谱，构建谱相关系数f(q)

基于t1时信号的估计振幅谱和实际振幅谱as(t1,f)，构建估计振幅谱和实际振幅谱之间的谱相关系数f(q)，如以下公式(12)所示：

其中,〈·，·〉表示内积算子，║·║表示l²范数，该公式用来定量测量估计振幅谱和实际振幅谱as(t1,f)之间的相似程度。对于任意一个实函数而言，它与它本身最相似且f(q)为0，这样便只要求解f(q)的最小值即可。

步骤5：求解谱相关系数f(q)最小时的q值

利用线性粘弹介质的深度z0+δz处的理论频谱去逼近深度z0+δz处实际记录的频谱，不断的修正衰减传播算子h(f,q,τ)，即不断的修正q值，使理论振幅谱与实际振幅谱达到最佳一致，即可得到介质的地震波衰减品质因子q。

实际处理时，可以给出一个q值搜索范围，在这个搜索范围内改变q值，可以搜索到f(q)最小时的q值。理论上，当f(q)为极值0时的q值即为估计的地震波衰减品质因子值，此时，理论振幅谱与实际振幅谱最为拟合。

应用示例

以下参考附图描述根据示例性实施例的基于广义s变换的谱相关系数的q值估计方法的应用效果。在本示例中，选择一个合成地震记录(见图2a)和一个实际地震记录(见图5)来测试本发明的q值估计的效果。

图2a和图2b分别示出了一个衰减地震道的合成地震记录及其广义s变换，反射界面在0.4s、0.8s、1.2s和1.6s，震源为30hz的ricker子波。加入少量随机噪声，通过广义s变换得到振幅谱。求取第二层的理论q值为50。

图3a和图3b分别示出了第二层上下界面的广义s变换振幅谱及其拟合后的广义s变换振幅谱。首先给定一个q值的范围，在这个范围内对q值进行扫描，得到每个q值对应的目标函数f的值。图4所显示目标函数f随q值变化的关系曲线：当q为49.8时，f达到最大值。此时，估计的广义s变换振幅谱和实际的广义s变换振幅谱最为相似。图3b显示估计q值与理论q值的误差为0.4％，该方法能有效估计q值，精度较高。

为考察该方法的抗噪性，在合成地震记录中加入噪声，分别用对数谱比法、峰值频移法以及本发明的方法估计q值。表1给出了未加入噪声和加入随机噪声情况下三种方法估计的q值，在三种方法中，理论q值均为50。结果表明，当信号中没有加入随机噪声时，三种方法估计的q值误差都较小；当信号中叠加随机噪声时，本发明的方法估计的q值误差最小。与谱比法和峰值频移法相比，本发明的方法具有较强的抗噪性。

表1

在叠后地震数据中应用根据示例性实施例的基于广义s变换的谱相关系数的q值估计方法。图5为某油田的保幅纯波叠后地震数据，其为剩余静校正后，地表一致性反褶积前的叠加数据，图中显示的是其中的第861道。对其进行广义s变换，结果如图6所示。应用本发明的方法估计q值，目的层在1350ms和1450ms之间。在实际地层中，含油气的地层的q值都较小。在设定q值范围时，认为q值大于500，即认为没有衰减，或者衰减可以忽略，因此q值估计范围为1-500。图7为所估计的q值，目的层处的q值较小，能比较好的指示油气储层。

上述技术方案只是本发明的一种实施例，对于本领域内的技术人员而言，在本发明公开的原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施例的描述，因此前面的描述只是优选的，而并不具有限制性的意义。