基于小波包变换的频谱成像方法与流程

本发明涉及石油勘探开发技术，尤其涉及一种基于小波包变换的频谱成像方法。

背景技术：

根据薄层调谐原理，不同频率振幅信息反映了不同厚度地层的调谐特征，因此地震记录实质上是不同地层厚度信息的叠合，也就是不同地质体的综合反映。对于薄砂层、河道和砂坝等隐蔽性油气藏，由于受其它地质体信息的干涉，在地震数据体上难以清晰识别，而频谱成像通过时频分析技术对地震数据进行频谱分解，对储层响应较好的频率成分进行成像，从而有利于目标地质体的识别。

频谱成像技术的技术核心是时频分析技术。目前，商业软件采用的时频分析技术基本上是以短时窗傅立叶变换和小波变换为主。众所周知，虽然短时窗傅立叶变换能对信号的局部特征进行分析，但它存在时窗固定，适应性较差的弊病。小波变换在短时窗傅立叶变换基础上进行了发展，它克服了短时窗傅立叶变换对时窗的限制，但它对频率的划分比较粗糙，在成像方面精度不够，从而影响了该项技术的推广应用。

吸收系数往往通过地震子波在介质中传播时的高频衰减进行估算，而已知的只有地震记录的振幅谱，它是地震子波和反射系数的综合结果。在估算地震子波谱时应消除反射系数谱的影响。

单一频率属性通常用于识别特定厚度的储层，而储层厚度在空间上是变化的，因此单一频率属性图在储层预测时存在一定局限。

技术实现要素：

本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种基于小波包变换的频谱成像方法，用于解决针对短时窗傅立叶变换时窗固定和小波变换对频率的划分比较粗糙的问题。

本发明解决其技术问题所采用的技术方案是：一种基于小波包变换的频谱成像方法，包括以下步骤：

1)在偏移后的三维地震数据体上进行加密小波包变换，将地震记录分解为按频带从小到大顺序排列的分频地震记录；

2)根据小波包变换结果计算分频瞬时振幅谱、分频瞬时相位谱、瞬时平均频率、瞬时带宽、瞬时吸收系数五中三维属性数据体；

3)在地震解释层位的控制下，进行属性切片；包括两种可选切片方式：以某一时间间隔δt沿层切片或将两个层位间数据剖分n等份切片；

4)每种属性可以单独分析，也可以对三种不同频率属性进行融合分析。

按上述方案，所述对偏移归位后地震数据进行按频带排序的加密小波包变换；

定义共轭滤波器组hn和gn，满足：

gk＝(-1)^kh1-k(3)

其中，n,k,l∈z，

离散地震信号sk在小波包基下的按频带排序的加密塔式分解方法为：

如果n为偶数，则

如果n为奇数，则

式中，n＝0,1,…2^l-1；l为小波包分解层数；为小波包分解系数，其中：(sk为地震信号)。

通过以上变换，将地震记录分解为按频带从小到大顺序排列的分频地震记录。

按上述方案，步骤2)中分频瞬时振幅谱和分频瞬时相位谱参数确定方法如下：

设s(t,a)是信号h(t)的小波包变换结果，其中a是频率通道号。

则其解析信号为

x(t,a)＝s(t,a)+ih(t,a)(6)

其中，h(t,a)是s(t,a)的希尔伯特变换，i为虚数单位，即

这里，*表示褶积运算；

将x(t,a)写成指数形式为

x(t,a)＝a(t,a)e^iφ(t,a)(8)

其中，a(t,a)为瞬时振幅谱，φ(t,a)为瞬时相位谱，则

φ(t,a)＝arg(x(t,a))+2kπ，k＝0,±1,±2,…(10)

arg(x(t,a))为解析信号x(t,a)的辐角主值。其范围为

-π＜arg(x(t,a))＜π

设k＝0，根据解析信号所在的象限位置按以下变换来计算瞬时相位谱：

在第一象限：s(t,a)＞0，h(t,a)＞0，

在第二象限：s(t,a)＜0，h(t,a)＞0，

在第三象限：s(t,a)＜0，h(t,a)＜0，

在第四象限：s(t,a)＞0，h(t,a)＜0，

按上述方案，步骤2)中从瞬时振幅谱a(t,a)可以得到功率谱密度函数p(t,a)：

p(t,a)＝a²(t,a)(11)

因为a对应着一定频带范围的中心频率，在这里不妨将p(t,a)记为p(t,ω)。

按上述方案，步骤2)中瞬时平均频率和瞬时带宽的确定方法如下:

瞬时平均频率定义为在某一时刻功率谱密度函数的期望值。即：

该值能够反映地层厚度变化和高频能量的衰减情况。

瞬时带宽定义为：

瞬时带宽能够指示总体吸收效果。

按上述方案，步骤2)中瞬时吸收系数的确定方法如下:

瞬时吸收系数的提取包括两个过程：1、地震子波频谱的提取；2、吸收系数的提取。下面分述这两个过程：

1)地震子波频谱的提取

设地震记录为子波与反射系数的褶积，即

s(t)＝w(t)*f(t)(14)

其中，

对上式进行傅里叶变换得:

s(w)＝w(w)f(w)(15)

取对数,得:

ln(s(w))＝lnw(w)+lnf(w)(16)

对上式进行反傅里叶变换得复赛谱:

地震子波相对于反射系数序列而言是比较平滑的、变化较慢的低频部分,所以它的复赛谱主要集中在时间原点附近，而反射系数序列是按时间均匀分布的，通过低通滤波获得地震子波的复赛谱，通过上述过程的逆过程可以由地震子波的复赛谱得到其振幅谱a(w)；

2)吸收系数的提取

在某一时刻，地震子波的振幅谱通常可以用图3表示，

地震子波的振幅谱在谱的高频段(ωd～ωn，ωd-主频，ωn-nyquist频率)谱近似满足指数函数(图3粗黑线条所示)

a(ω)＝cexp(-q(t)ω)(18)

其中：c为常数，q(t)为瞬时吸收系数，ω为角频率。

因此，可以对振幅谱的高频段进行指数拟合来模拟吸收系数q(t)值。

式(18)经过对数变换可以化为线性拟合。

ln(a(ω))＝lnc-q(t)ω(19)

令离散化后的地震子波高频段的(频率，对数振幅谱)对为(f1,a1)，(f2,a2)，(f3,a3)，…，(fn,an)，则

其中：fi为频率，ai为对应的对数振幅谱。

这样就得到了瞬时吸收系数。

本发明产生的有益效果是：本发明通过采用按频带排序的加密小波包变换技术，提高了频带划分精度；研究了分频瞬时吸收系数反演技术，在估算地震子波谱方面方法更为合理；通过将三种不同频率属性叠合显示在一张图上，提升了通过频谱成像技术对储层的识别能力。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例的方法流程图；

图2是本发明的小波包分解算法示意图；

图3是本发明的地震子波的振幅谱；

图4是本发明实施例的沿层切片图；

图5是本发明实施例三频率属性rgb融合显示图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明实施例基于小波包变换的频谱成像方法主要包括以下步骤：

s1：在偏移后三维地震数据体上进行按频带排序的加密小波包变换；

定义共轭滤波器组hn和gn，满足：

gk＝(-1)^kh1-k(3)

其中，n,k,l∈z，

离散地震信号sk在小波包基下的按频带排序的加密塔式分解算法为：

如果n为偶数，则

如果n为奇数，则

式中，n＝0,1,…2^l-1；l为小波包分解层数；为小波包分解系数，其中：(sk为地震信号)。

通过以上变换，将地震记录分解为按频带从小到大顺序排列的分频地震记录。

s2：计算分频瞬时振幅谱、分频瞬时相位谱、瞬时平均频率、瞬时带宽、瞬时吸收系数；

分频瞬时振幅谱和分频瞬时相位谱:

设s(t,a)是信号h(t)的小波包变换结果，其中a是频率通道号。

则其解析信号为

x(t,a)＝s(t,a)+ih(t,a)(6)

其中，h(t,a)是s(t,a)的希尔伯特变换，即

这里，*表示褶积运算。

将x(t,a)写成指数形式为

x(t,a)＝a(t,a)e^iφ(t,a)(8)

其中，a(t,a)为瞬时振幅谱，φ(t,a)为瞬时相位谱，则

φ(t,a)＝arg(x(t,a))+2kπk＝0,±1,±2,…(10)

arg(x(t,a))为解析信号x(t,a)的辐角主值。其范围为

-π＜arg(x(t,a))＜π

通常我们研究φ(t,a)在(-π,π)上的性质，即在(10)式中k＝0时的性质。

arg(x(t,a))可以通过反正切函数表示，而反正切函数是多值函数，其主值定义在上，与前面解析信号x(t,a)的辐角主值范围(-π,π)不一致，也就是说仅仅由反正切函数还不能完全确定解析信号x(t,a)的相位谱。为此，我们根据解析信号所在的象限位置做适当的变换来计算瞬时相位谱。

在第一象限：s(t,a)＞0，h(t,a)＞0，

在第二象限：s(t,a)＜0，h(t,a)＞0，

在第三象限：s(t,a)＜0，h(t,a)＜0，

在第四象限：s(t,a)＞0，h(t,a)＜0，

以上求出的瞬时振幅谱和瞬时相位谱实际上反映了频率通道a的信息，而a对应着以某一频率ω为中心频率的一定频带范围，当频带划分足够精细的情况下，能够近似反映中心频率的信息。

从瞬时振幅谱a(t,a)可以得到功率谱密度函数p(t,a)：

p(t,a)＝a²(t,a)(11)

因为a对应着一定频带范围的中心频率，在这里不妨将p(t,a)记为p(t,ω)。

瞬时平均频率:

瞬时平均频率定义为在某一时刻功率谱密度函数的期望值。即：

该值能够反映地层厚度变化和高频能量的衰减情况。

瞬时带宽:

瞬时带宽定义为：

瞬时带宽能够指示总体吸收效果。

瞬时吸收系数:

由robinson褶积模型可知，地震记录的频谱是地震子波的频谱和地层反射系数序列的频谱的乘积。在此，我们主要研究地震子波在传播过程中由于地层的的吸收作用引起的能量衰减情况。因此首先要消除地层反射系数序列频谱的影响。

归结起来，瞬时吸收系数的提取包括两个过程：1、地震子波频谱的提取；2、吸收系数的提取。下面分述这两个过程：

1)地震子波频谱的提取

复赛谱分析是一种非线性滤波,用以分解两个相褶积的过程。要分解的两个过程应具有不相重叠的谱函数或者说两个过程变化率不同,一个变化快,一个变化慢。地震反射波记录可视为地震子波与反射系数的褶积,脉冲地震记录具有白噪性质,变化快,而地震子波一般较为稳定,变化慢。。

设地震记录为子波与反射系数的褶积，即

s(t)＝w(t)*f(t)(14)

对上式进行傅里叶变换得:

s(w)＝w(w)f(w)(15)

取对数,得:

ln(s(w))＝lnw(w)+lnf(w)(16)

对上式进行反傅里叶变换得复赛谱:

地震子波相对于反射系数序列而言是比较平滑的、变化较慢的低频部分,所以它的复赛谱主要集中在时间原点附近，而反射系数序列是按时间均匀分布的。因此可通过低通滤波来获得地震子波的复赛谱。通过上述过程的逆过程可以由地震子波的复赛谱得到其振幅谱a(w)。

2)吸收系数的提取

在某一时刻，地震子波的振幅谱通常可以用图3表示，

在谱的高频段(ωd～ωn，ωd-主频，ωn-nyquist频率)谱近似满足指数函数(图3粗黑线条所示)

a(ω)＝cexp(-q(t)ω)(18)

因此，可以对振幅谱的高频段进行指数拟合来模拟吸收系数q值。

上式经过对数变换可以化为线性拟合。

ln(a(ω))＝lnc-q(t)ω(19)

令离散化后的地震子波高频段的(频率，对数振幅谱)对为(f1,a1)，(f2,a2)，(f3,a3)，…，(fn,an)，则

这样就得到了瞬时吸收系数。

s3：在地震解释层位的控制下，进行属性切片；可以两种方式切片：1、以某一时间间隔δt沿层切片；2、将两个层位间数据剖分n等份切片。

s4：每种属性可以单独分析，也可以对三种不同频率属性进行融合分析。

图4是某工区河道砂体的小波包频谱成像处理结果。从上至下依次为40hz-100hz瞬时振幅谱平面图，右下图是由钻井资料解释得到的砂岩等厚图，反映了物源方向和河道延展形态。从60hz瞬时振幅谱平面图可以看出，河道砂体刻画的非常清晰，与由钻井资料解释得到的砂岩等厚图比较吻合。在高频部分(100hz瞬时振幅谱图上)，河道仍能清晰地显现，说明本方法在高频部分具有较高的分辨能力。由此可以看出，分频属性很好地利用了地震记录的有效高频信息，使对薄层的平面展布研究成为可能。

单一频率属性通常用于识别特定厚度的储层，而储层厚度在空间上是变化的，需要多张分频属性图来反映储层的分布特征。将分频成果中频段互不重叠的低频段、中频段、高频段能量属性以rgb模式混合起来显示，利用颜色混合特有的效果来突出目标，能够有效改善单一频率属性对目标的分辨效果。然后在该数据体上进行储层的刻画与分析。将rgb混合模式引入到多分频能量属性的融合显示中。

图5-b是用33hz、48hz和60hz瞬时振幅谱合成的三原色混频显示图。与原始瞬时振幅谱图(图5-a)相比，提供了更加丰富的信息，除了主河道显示更加清楚外，一些细小构造也能反映出来，为后续做进一步精细解释提供了更多的依据。

综上，本发明基于小波包变换的频谱成像技术主要改进包括三点：(1)按频带排序的加密小波包变换技术，克服了短时窗傅立叶变换时窗固定和小波变换对频率的划分比较粗糙的弊病，分频数据精度更高；(2)在计算瞬时吸收系数时，通过复赛谱技术消除反射系数谱的影响，估算的瞬时吸收系数更为合理。(3)三种不同频率属性的融合显示利用颜色混合特有的效果来突出目标，能够有效改善单一频率属性对目标的分辨效果。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。