一种适用于AVAF反演的频率域地震道集生成方法与流程

本发明属于石油地球物理勘探领域，具体涉及一种适用于AVAF反演的频率域地震道集生成方法，利用高效的时频分析工具来生成适合于AVAF反演的频率域地震道集。该地震道集不同于时频分析中频谱的概念，是完整的地震道集。

背景技术：

AVAF反演技术是近几年发展起来的粘弹性介质中参数反演的方法，它类似于常规的AVA反演。常规的AVA反演求解的是弹性参数，而AVAF反演需要再反演两个粘弹性参数：P波和S波衰减值，这需要用到频率信息，也就是波形的信息。常规的分频技术包括了傅里叶变换、小波变换、加窗傅里叶变换等等，它们都存在着一个共同的缺陷，分频产生的谱不具有波形的信息，同时会产生相位缠绕现象，这不适用于AVAF反演。

除了这么一个共同的缺点，常规的傅里叶变换、小波变换和固定时窗的加窗傅里叶变换也存在着各自的缺陷。常规的傅里叶变换(FT)假设信号或者图像是稳定的，这不能适应于分析非平稳信号；短时傅里叶变换(STFT)是其中一种固定时窗的时频变换方法，不能在时间和频率上做一个合理的折衷；小波变换(WT)引入了多尺度的概念，使得在对信号进行小波分解时高频具有好的时间分辨率，低频具有比较好的频率分辨率；然而小波变换存变换域不再是频率域，而只是尺度的概念；其次，由于小波变换基的特性导致了相位的概念在小波域中是局部的。这两点与AVAF分析的需要是矛盾的。

技术实现要素：

本公开的目的在于解决上述现有技术中存在的难题，提供一种适用于AVAF反演的频率域地震道集生成方法，通过将地震道集分解成频率域的地震道集，更好地利用了地震道集中的波形信息，也就是频率的信息。这使得下一步的AVAF反演在频率域能够实现，同时用本公开方法生成的数据里可以在频率域或者分道集进行显示，分频道集比常规的分频剖面更加直观，具有波形信息，使得储层预测和含油气性预测更加可信。

根据本公开的一个方面，提供一种适用于AVAF反演的频率域地震道集生成方法，其特征在于，包括以下步骤：获取地震数据；对所述地震数据进行S变换，得到地震数据的S谱；针对反演需求设计窗函数，与S谱相乘再积分；对相乘的结果进行傅里叶反变换，得到频率域道集。

进一步地，所述S变换的公式如下：

$S(\mathrm{\τ},v)=\underset{-\∞}{\overset{+\∞}{\∫}}g\left(t\right)\frac{\left|v\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{{(\mathrm{\τ}-t)}^2{v}^2}{2}}{e}^{-i2\mathrm{\πvt}}\mathrm{dt}---\left(1\right)$

其中τ和v分别是S变换之后的时间坐标和频率坐标。

进一步地，由于离散S变换(DST)的计算量非常大，因为它不是一个正交变换，同时S变换得到了很多的冗余信息。所述对所述地震数据进行S变换使用有限频率域快速S变换方法，选择地震频带作为需要处理的频率进行分解，则地震数据的S变换可以从α域经过反傅里叶变换得到：

$S(\mathrm{\τ},v)=\underset{-\∞}{\overset{+\∞}{\∫}}\mathrm{\α}(v^{\′},v)\·e^{i2\mathrm{\π}{v}^{\′}\mathrm{\τ}}{\mathrm{dv}}^{\′}---\left(2\right).$

进一步地，所述针对反演需求设计窗函数，与S谱相乘再积分的步骤中，高斯窗的选择依赖于地震数据，同时与反演子波的高斯窗函数保持一致；高斯窗函数的个数依赖于反演的频率数，当反演的频率数为N时，窗函数的个数也是N。

所述针对反演需求设计窗函数，与S谱相乘再积分的步骤，具体公式如下式：

$g_{v_0}\left(t\right)={\∫}_{-\∞}^{\∞}\left\{{\∫}_{-\∞}^{\∞}S(\mathrm{\τ},v)\mathrm{d\τ}\right\}{e}^{i2\mathrm{\πvt}}W(v_0,{\mathrm{\σ}}^2)\mathrm{dv}---\left(3\right)$

其中S(τ,ν)是待分频道集的S谱，W(ν₀,σ²)为窗函数，ν₀为窗函数的均值频率，σ²为窗函数的方差。

进一步地，通过公式(3)得到的加窗反变换结果就是最后得到的分频率后的地震道集，该地震道集配合频率域的子波适于进行AVAF反演。

首先使用公式(1)定义的S变换将地震道集变换到时频域S域。在对时频域的数据进行反傅里叶变换时，首先将其在局部时间域上进行全局积分，得到了公式(3)中的然后在频率域进行加窗积分得到$\left\{{\∫}_{-\∞}^{\∞}S(\mathrm{\τ},v)\mathrm{d\τ}\right\}W(v_0,{\mathrm{\σ}}^2).$

相比于已有的所有的时频分析类的方法，该方法拥有独特的优势，那就是使用该发明得到的分频率的地震道集具有波形特征，保留了地震波的相位特征，这极大便利了在频率时间域双域进行AVAF反演。同时该方法配合分频率的地震子波能够完整反演分频率的反射系数，使得后期的油藏定量描述估计成为可能。

本公开的其它特征和优点将在随后的具体实施方式部分予以详细说明。

附图说明

通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。

图1是根据本公开一种实施方式的方法步骤框图。

图2是S谱在频率域的分频段积分示意图。

图3是5Hz的频率地震道集。

图4是35Hz的频率地震道集。

具体实施方式

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

根据本发明的一种实施方式，提供一种适用于AVAF反演的频率域地震道集生成方法，包括：获取地震数据；对所述地震数据进行S变换得到地震数据的S谱；针对反演需求设计窗函数，与S谱相乘再积分；对相乘的结果进行傅里叶反变换，得到频率域道集。

具体地，首先获取地震数据g(t)。然后对地震数据g(t)进行S变换得到地震数据的S谱S(τ,v)。在本实施方式中，S变换正变换的公式如下：

$S(\mathrm{\τ},v)=\underset{-\∞}{\overset{+\∞}{\∫}}g\left(t\right)\frac{\left|v\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{{(\mathrm{\τ}-t)}^2{v}^2}{2}}{e}^{-i2\mathrm{\πvt}}\mathrm{dt}---\left(1\right)$

其中τ和v分别是S变换之后的时间坐标和频率坐标。

由于离散S变换(DST)的计算量非常大，因为它不是一个正交变换，同时S变换得到了很多的冗余信息。优选地，为了提高处理速度，本实施方式可使用有限频率域快速S变换方法。选择需要处理的频率范围内的频率进行分解，该频率范围的选择一般情况下与地震频带重合，则原地震数据的S变换可以从α域经过反傅里叶变换得到：

$S(\mathrm{\τ},v)=\underset{-\∞}{\overset{+\∞}{\∫}}\mathrm{\α}(v^{\′},v)\·e^{i2\mathrm{\π}{v}^{\′}\mathrm{\τ}}{\mathrm{dv}}^{\′}---\left(2\right)$

但是直接用有限频率域快速S变换生成的S谱并不能直接用于AVAF反演，因此需要将其反变换回时间域。如果对单个的频率进行反变换，那么将产生吉布斯效应，导致时间域的地震信号会产生振荡，也不能用于AVAF反演。

本实施方式采用的方法是将S谱与不同均值的高斯窗相乘，再进行S反变换，就能得到适用于AVAF反演的频率域地震道集。高斯窗的选择依赖于地震数据，同时还需要与反演子波的高斯窗函数保持一致，而高斯窗函数的个数依 赖于反演的频率数，当反演的频率数为N时，窗函数的个数也是N。具体公式如下式：

$g_{v_0}\left(t\right)={\∫}_{-\∞}^{\∞}\left\{{\∫}_{-\∞}^{\∞}S(\mathrm{\τ},v)\mathrm{d\τ}\right\}{e}^{i2\mathrm{\πvt}}W(v_0,{\mathrm{\σ}}^2)\mathrm{dv}---\left(3\right)$

其中S(τ,ν)是待分频道集(输入的地震数据)的S谱，W(ν₀,σ²)为窗函数，ν₀为窗函数的均值频率，σ²为窗函数的方差。均值频率的选择依赖于需要反演的频率参数，窗函数的方差依赖于经验和实际数据。

通过公式(3)得到的加窗反变换结果就是本实施方式最后得到的分频率后的地震道集，该道集配合频率域的子波可以进行AVAF反演，为后续的粘弹性参数的估计提供输入数据。

下面参照附图详细描述本发明的具体实施例。图1是根据本公开一种实施方式的方法步骤框图。针对每个地震道集可以分为以下几个步骤来实现：

步骤S101：输入一个地震道集数据g(t)。

步骤S102：对地震道集数据g(t)进行S变换得到时间频率域的S谱S(τ,ν)。在本实施例中，S变换正变换的公式如下：

$S(\mathrm{\τ},v)=\underset{-\∞}{\overset{+\∞}{\∫}}g\left(t\right)\frac{\left|v\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{{(\mathrm{\τ}-t)}^2{v}^2}{2}}{e}^{-i2\mathrm{\πvt}}\mathrm{dt}---\left(1\right)$

其中τ和v分别是S变换之后的时间坐标和频率坐标。

作为优选实施例，对地震道集数据g(t)进行有限频率域快速S变换，选择需要处理的频率范围与地震频带重合，则从α域经过反傅里叶变换得到原地震道集数据g(t)的S变换：

$S(\mathrm{\τ},v)=\underset{-\∞}{\overset{+\∞}{\∫}}\mathrm{\α}(v^{\′},v)\·e^{i2\mathrm{\π}{v}^{\′}\mathrm{\τ}}{\mathrm{dv}}^{\′}---\left(2\right)$

步骤S103：针对反演需要设计一系列的高斯窗函数并确定窗函数的个数(例如用于反演的频率个数为5个则窗函数的个数也是5个)，用高斯窗函数分别乘以时间频率域的S谱S(τ,ν)，再进行积分。

$g_{v_0}\left(t\right)={\∫}_{-\∞}^{\∞}\left\{{\∫}_{-\∞}^{\∞}S(\mathrm{\τ},v)\mathrm{d\τ}\right\}{e}^{i2\mathrm{\πvt}}W(v_0,{\mathrm{\σ}}^2)\mathrm{dv}---\left(3\right)$

其中S(τ,ν)是待分频道集的S谱，W(ν₀,σ²)为窗函数，ν₀为窗函数的均值频率，σ²为窗函数的方差，其中均值频率的选择依赖于需要反演的频率参数，窗函数的方差依赖于经验和实际数据。

步骤S104：对步骤3中的结果进行反傅里叶变换，这样就得到了在特定均值频率f₀时候的频率道集。

将本发明应用于一个实际地震数据上，首先使用公式(1)定义的S变换将地震道集变换到时频域S域，如图2所示，这是某一道的一个时间点上的频谱图，在对时频域的数据体进行反变换时，首先将其在局部时间域上进行全局积分，得到了公式(3)中的然后在频率域进行加窗积分就得到了图2中展示了四个高斯窗函数(四条虚线所示)，它们具有相同的展布范围，不同的均值频率。

下一步将加窗后的数据进行傅里叶反变换，得到了分频率的地震道集，如图3和图4所示。图3展示了5Hz分频地震数据，可以看到它的同相轴比图4中35Hz分频地震数据更加粗，同时低频时在深层时衰减较小，而35Hz的高频剖面中，深层能量要明显弱于表层能量，这表明这种分频率地震道集生成方法不仅具有适用于AVAF反演的波形特征，还保留了地震数据衰减的自然属性。

本公开的方法通过将地震道集分解成分频率的地震道集，使得粘弹性介质中的AVAF反演成为可能。本发明首先用有限频率域快速S变换算法将地震道集变换到时频域，由于时频域的数据不能直接用于AVAF反演。然后将时频域的数据乘以一个频率域的高斯窗函数，这不仅能够保留分频率的特征，而且极大程度降低单个频率造成的吉布斯现象。最后通过反变换得到分频率的地震道集，该道集保留了地震信号的相位特征，可以用于后期的反演与含气性的定性估计。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的 情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中技术的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。