一种基于谱面积变化的地震衰减估计方法

摘要

本发明公开了一种基于谱面积变化的地震衰减估计方法，包括步骤1：利用常相位子波来估计实际震源子波及其参数；步骤2：计算震源子波的傅里叶谱面积；步骤3：计算衰减子波的傅里叶谱面积；步骤4：计算震源子波和衰减子波的傅里叶谱面积变化量；步骤5：计算相对Q值，进行地震数据衰减估计。本发明利用常相位子波估计震源子波并得到估计参数，利用估计得到的震源子波和衰减子波之间傅里叶谱面积的变化量，计算得到相对Q值，进行地震衰减估计和衰减储层刻画。通过对无噪的和含噪的合成地震记录进行衰减估计，验证了本方法的准确性和稳定性。将该方法应用于实际资料，其衰减估计结果能够很好地刻画含气储层。

说明书

技术领域

本发明属于地球物理勘探中的信号处理、数据处理和解释领域，涉及一种基于谱面积变化的地震衰减估计方法。

背景技术

地下介质是黏弹性介质，介质的黏弹性使得在地下介质中传播的地震子波的波形展宽，主频降低，频带变窄，能量逐渐减弱并最终转化为热能。地震波能量的损耗由多种因素引起，具有不同饱和度、渗透率、孔隙度的介质对地震波能量的损耗程度不同。品质因子Q值是度量介质衰减特性的重要参数，可用于指示地层含油气性。由地震资料估计的地层衰减，可用于岩性及含流体分析、储层识别和烃类检测等。

Tonn(1991)在其综述文章中比较了多种Q值估计方法，包括时间域的子波上升时间法、脉冲振幅衰减法、解析信号法、子波/相位模拟法等，以及频率域的匹配法、频谱模拟法、谱比法等。Tonn指出，没有一种方法能适用于所有的情况，根据噪声水平和接收记录的类型，不同的方法适用于不同的情况，例如当VSP资料的震源可控且资料保幅时，解析信号法估计的结果优于其他方法的估计结果，而当资料不保幅时，谱模拟法更适用；在不含噪的情况下，谱比法可以得到可信的结果，但在含噪情况下该方法不再适用。谱比法是最早提出的方法之一，该方法原理简单，操作方便，在信噪比较高的情况下可以取得较好的结果，得到了广泛的应用。Singleton等(2006)使用Gabor-Morlet谱比法估计Q值并进行Q补偿。Parra等(2006)基于谱比法从直达P波中估计了固有衰减，估计的Q曲线结合其他测井曲线，可以区分异常是与岩性有关还是与含油气饱和度有关。上升时间法和脉冲振幅衰减法对采样率和信噪比依赖性较大，估计结果通常误差较大，而匹配技术法与对数谱比法一样，适用于信噪比较高的资料。

除了Tonn列出的方法外，时间域方法还有脉冲宽度展宽法等，频率域有质心频率偏移法(CFS)、峰值频率移动法(PFS)、子波包络峰值瞬时频率法(WEPIF)等，以及时频域方法，如小波域峰值尺度法、小波域能量衰减法、Lipschitz指数法、时频谱谱比法等。脉冲宽度展宽法由Kjartansson(1979)提出，该方法利用脉冲宽度、Q值、走时之间的经验公式估计Q值，一定程度上减少了噪声干扰；质心频率偏移法与峰值频率移动法分别用质心频率及峰值频率与走时之间的解析关系估计Q值。质心频率的定义是用振幅谱加权的平均频率，质心频率偏移法适用于频带较宽的地震记录，而峰值频率定义为频谱的最大值。子波包络峰值瞬时频率法利用子波的瞬时振幅包络的峰值处对应的瞬时频率与Q值的关系估计Q值，该方法适用于叠后地面反射资料和零偏VSP资料，操作方便，抗噪性强，纵向分辨率高。此外，朱定等(2006)利用主频向低频移动的偏移量与Q值之间的比例关系估计衰减，该方法适用于VSP资料，计算速度快。在频率域中估计Q值是通过提取不同的频率信息如主频、峰值频率、中心频率、瞬时频率、包络峰值处的瞬时频率等，并寻找这些特定的部分频率信息与Q值之间的关系来求取Q值的，不同的频率信息具有不同的稳定性、抗噪性及提取方便性。时频域方法中，李宏兵等用小波变换域的尺度因子(即能量衰减因子)来刻画衰减，从峰值尺度的变化量中直接估计衰减，Innanen(2003)利用局部信号正则化和Lipschitz指数估计衰减，给出了Q值估计的经验公式。

前常用的Q值提取方法，如质心频率偏移法(CFS)和峰值频率移动法(PFS)，都是基于单一的、特定的频率变化，来进行Q值估计。在无噪声的理想情况下，可以得到精确的、稳定的衰减估计结果；但是实际地震数据不可避免地含有噪声，这些常规方法很难得到一个稳定的Q值估计结果。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种基于谱面积变化的地震衰减估计方法

为达到上述目的，本发明采用以下技术方案予以实现：

一种基于谱面积变化的地震衰减估计方法，包括以下步骤：

步骤1：利用常相位子波来估计实际震源子波及其参数；

步骤2：计算震源子波的傅里叶谱面积；

步骤3：计算衰减子波的傅里叶谱面积；

步骤4：计算震源子波和衰减子波的傅里叶谱面积变化量；

步骤5：计算相对Q值，进行地震数据衰减估计。

本发明进一步的改进在于：

步骤1中估计震源子波的具体方法如下：

常相位子波的定义式为：

式中，σ₀和δ₀表示常相位子波的主频和标准差；相位φ₀为常数，根据不同的地震数据选取不同的值；i为虚数单位；假设实际震源子波的振幅谱为根据式(2)估计得到估计的震源子波：

式中，为估计的震源子波。

步骤2中计算震源子波的傅里叶谱面积的具体方法如下：

对式(1)两边关于频率积分，得估计的震源子波的傅里叶谱面积如下

其中，f为频率，e为指数函数。

步骤3中计算衰减子波的谱面积的具体方法如下：

设子波在品质因子Q的衰减介质中传播时间t秒，接收子波的振幅谱为

其中，t为子波在衰减介质中传播的双程旅行时；

式(4)两边同样对频率积分，得衰减子波的谱面积为

步骤4中计算震源子波和衰减子波的谱面积差异的具体方法如下：

对式(3)中估计的震源子波傅里叶谱面积和式(5)中衰减子波的傅里叶谱面积作差，得谱面积变化量为

步骤5中计算相对Q值的具体方法如下：

利用式(6)中估计的震源子波的谱面积和衰减子波的谱面积变化，求解得到相对Q值，表达式如下：

选取S变换来计算时变地震数据的傅里叶谱；实地震信号s(t)的S变换的定义为

式中，t和f为子波在衰减介质中传播的双程旅行时和傅里叶频率；τ为时间平移参数，即S变换中高斯窗的中心；根据帕斯瓦尔定理，式(8)中S变换的定义改写为：

式中，S(α)为地震信号s(t)的傅里叶变换；根据式(9)，利用快速傅里叶变换，快速地计算S变换；利用式(9)中的S变换，计算地震数据的傅里叶谱面积，进行地震数据衰减估计。

与现有技术相比，本发明具有以下有益效果：

本发明提出了一种基于傅里叶谱面积变化的地震衰减估计方法。首先该方法利用常相位子波估计震源子波并得到估计参数，利用估计得到的震源子波和衰减子波之间傅里叶谱面积的变化量，计算得到相对Q值，进行地震衰减估计和衰减储层刻画。将该发明应用于合成地震数据和实际地震数据，证明了该方法可以精确地刻画地震数据衰减，并且对含噪数据可以得到一个稳定的、准确的衰减估计结果。本发明充分考虑地震数据的傅里叶谱面积的变化，来进行地震数据的衰减估计，可以得到一个更为稳定的估计结果。本发明通过对无噪的和含噪的合成地震记录进行衰减估计，验证了本方法的准确性和稳定性。将该方法应用于实际资料，其衰减估计结果能够很好地刻画含气储层。

附图说明

图1无噪合成地震数据；共100道数据，时间采样间隔为1毫秒，该合成地震数据含有两水平层；震源子波为主频为40Hz的零相位Ricker子波；其中，第一层Q值为200，第二层Q值分别为100(第一道到第40道)，30(第41道到第60道)，60(第61道到第100道)；

图2实际震源子波(主频为40Hz的Ricker子波)和零相位子波估计得到的震源子波；其中，(a)实际震源子波和估计的震源子波的时域波形，(b)实际震源子波和估计震源子波的傅里叶谱；

图3图2中地震数据振幅谱的变化和本发明估计的相对Q值；其中，(a)不同层、不同道的地震数据的傅里叶谱，暗实线表示第一层第20道数据的傅里叶谱，深实线、虚线、短划线分别表示第二层第20道、第50道和第80道数据的傅里叶谱；(b)本发明估计得到的相对Q值；

图4含噪合成地震数据；图1中的合成地震数据，加入信噪比为5dB的高斯白噪声后得到该含噪合成地震数据；其他参数与图1中相同；

图5图4中含噪地震数据振幅谱的变化和本发明估计的相对Q值；其中，(a)不同层、不同道地震数据的傅里叶谱，暗线表示第一层第20道数据的傅里叶谱，深实线、虚线、短划线分别表示第二层第20道、第50道和第80道地震数据的傅里叶谱；(b)本发明估计得到的相对Q值；

图6实际地震数据剖面；井1和井2的CMP号为290和700；井1为高产井，井2为干井；

图7时间域和频率域的震源子波和本发明估计的相对Q值；其中，(a)实际震源子波(由浅层数据提取得到)和估计的震源子波的时域波形，(b)实际震源子波和估计震源子波的傅里叶谱，(c)质心频率偏移法(实线)和谱比法(虚线)估计得到的品质因子Q值，(d)本发明估计得到的相对Q值。

具体实施方式

下面结合附图对本发明做进一步详细描述：

参见图，一种基于谱面积变化的地震衰减估计方法

(1)震源子波估计

在实际地震数据处理和衰减估计中，往往震源子波是未知的，因此本发明中，我们首先利用常相位子波来估计实际震源子波，并得到估计的震源子波的参数，进而进行地震数据衰减估计。常相位子波的定义式为

式中，σ₀和δ₀表示常相位子波的主频和标准差。相位φ₀为常数，可以根据不同的地震数据选取不同的值。假设实际震源子波的振幅谱为我们根据式(2)可估计得到估计的震源子波

式中，为本发明估计的震源子波。

(2)估计得到的震源子波的谱面积计算

对式(1)两边关于频率积分，可得估计的震源子波的傅里叶谱面积如下

(3)衰减子波的谱面积计算

设子波在品质因子Q的衰减介质中传播时间t秒，接收子波的振幅谱为

式(4)两边同样对频率积分，可得衰减子波的谱面积为

(4)震源子波和衰减子波的谱面积差异计算

对式(3)中估计的震源子波傅里叶谱面积和式(5)中衰减子波的傅里叶谱面积作差，可得谱面积变化量为

式中，t为子波在衰减介质中传播的双程旅行时。

(5)相对Q值计算

利用式(6)中估计的震源子波的谱面积和衰减子波的谱面积变化，可以求解得到相对Q值，表达式如下

由于估计震源子波和估计参数的误差，该Q值只是一个相对Q值，并不是真正的品质因子，但是本发明认为，该相对Q值仍然可以较好地刻画地震数据的衰减特性。本发明中的合成地震数据实验和实际地震数据实验可证明本结论。

另外，本发明中我们选取S变换来计算时变地震数据的傅里叶谱。实地震信号s(t)的S变换的定义为

式中，t和f为时间和傅里叶频率。τ为时间平移参数，即S变换中高斯窗的中心。根据帕斯瓦尔定理，式(8)中S变换的定义可改写为

式中，S(α)为地震信号s(t)的傅里叶变换。根据式(9)，利用快速傅里叶变换(FFT)，可以快速地计算S变换。利用式(9)中的S变换，我们可以计算地震数据的傅里叶谱面积，进行地震数据衰减估计。

实施例：

(1)合成地震数据算例

首先，我们用无噪和含噪的合成地震数据验证本发明的有效性。图1是无噪的合成地震数据，共有100道地震数据，其中含有两个水平层。震源子波为主频为40Hz的零相位Ricker子波。时间采样间隔为1毫秒。第一层的Q值为200，第二层的Q值分别为100(第1道到第40道)，30(第41道到第60道)，60(第61道到第100道)。首先我们用式(1)中的常相位子波根据式(2)估计震源子波，其结果如图2所以。图2(a)给出了实际震源子波和估计得到的常相位子波(即估计得到的震源子波)的时域波形，图2(b)给出了实际震源子波和估计得到的常相位子波(即估计得到的震源子波)的傅里叶谱，二者之间的差异可忽略，可见当震源子波未知时，我们可以用常相位子波估计得到震源子波，并得到估计参数。图3(a)给出了不同层、不同道的地震数据的傅里叶谱，可见Q值不同导致衰减子波的傅里叶谱也是不尽相同的，因此本发明利用不同衰减子波的傅里叶谱面积变化来进行Q值估计是有理论根据的。图3(b)给出了本发明估计的相对Q值，虽然该Q值不是理论的、真实的品质因子，但是该相对Q值仍可以准确地刻画地震数据衰减的相对关系，能够很好地刻画不同地层的衰减。

因为实际处理中地震数据不可避免地含有噪声，因此我们尝试将本发明应用于含噪的合成地震数据来进行衰减估计，以证明该方法的有效性和稳定性。图4是在图1的合成地震数据加入信噪比为5dB的高斯白噪声得到的含噪合成地震数据，可见在噪声影响下，地震数据的质量很差，这必然给地震衰减估计带来很大的影响。图5(a)给出了不同层、不同道的含噪地震数据的傅里叶谱，可见噪声对傅里叶谱的严重干扰，如采用质心频率偏移法(CFS)等利用单一的、特定的频率变化来进行衰减估计，必然带来很大的估计误差，很难得到一个准确的估计结果。本发明方法估计得到的相对Q值如图5(b)所示，其结果能够很好地揭示不同层和不同衰减的相对关系，可见本发明在含噪情况下，仍能够得到一个稳定的估计结果，该相对Q值的估计结果仍然可以准确地反映含噪合成地震数据的衰减差异。

(2)实际地震资料算例

我们将本发明的估计方法用于实际地震数据，尝试进行衰减估计。该2D实际地震数据如图6所示，其时间采样间隔为1毫秒，共有600道数据。两口井的CDP号分别为290和700(图6中竖线标注处)，其中井1为高产井，含9.2米厚的含气砂岩，井2为低产井，含有8米厚的砂岩，但测井显示结果并不含气。根据式(2)估计得到震源子波如图7(a)和(b)所示，根据实际震源子波和估计的震源子波之间可以忽略的差异，可见本发明中估计震源子波的正确性。注意，由于真实震源子波未知，该实际震源子波是由浅层地震数据提取而得到。我们用常用的质心频率偏移法(CFS)和谱比法(LRS)作为对比方法，其结果如图7(c)所示，在井1附近并未出现低Q值，即强衰减的现象，在井2附近Q值有所下降，即衰减仍较强，这跟测井结果并不吻合。本发明估计得到的相对Q值结果如图7(c)所示，注意290道附近的低Q值，即强衰减，700道附近的高Q值，即弱衰减。本发明估计得到的相对Q值与测井结果吻合，证明了该方法的正确性和有效性，该方法可用于衰减储层的刻画和预测。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。