一种基于变相位雷克子波匹配追踪的理想地震谱分解方法

摘要

本发明公开了一种基于变相位雷克子波匹配追踪的理想地震谱分解方法，首先构建变相位雷克子波库；然后根据匹配追踪方法求取匹配子波及对应的振幅系数，并将振幅系数与子波原子相乘作为分解得到的子信号；接着对单反射界面对应的子信号做傅里叶变换，计算该界面反射地震波的振幅谱和旋转相位谱，将傅里叶变换谱放置在对应的延迟时间处，最终得到理想的时频谱分解结果。该时频分解结果包括振幅谱和旋转相位谱两部分，旋转相位谱能够提供更为丰富的地质和储层信息。振幅谱分解结果精确反映了地震子波到达界面的延迟时间及其能量随频率变化的分布，不受测不准原理极限的限制，具有极高的时间分辨率和频率分辨率，对于地震信号来说是一种理想时频谱。

说明书

技术领域

本发明属于油气地球物理勘探领域，具体涉及一种基于变相位雷克子波匹配追踪的理想地谱分解方法。

背景技术

地震谱分解技术将时间域地震信号转换到时频域进行分析，能够从地震数据中获取更为丰富的地层信息，在确定地层厚度、沉积相解释、油气检测等方面得到了广泛应用。

常见的地震谱分解方法主要包括线性类时频分析、双线性类时频分析和自适应时频分析方法。线性类时频分析方法即核函数分解法，利用窗函数实现了信号的时频局部性，主要包括短时傅里叶变换(STFT)、连续小波变换(CWT)、 S变换(ST)和广义S变换(GST)。线性类时频分析方法受Heisenberg不确定性原理限制，其时间分辨率和频率分辨率相互制约和牵制，无法同时使时间分辨率和频率分辨率达到任意小。双线性类时频分析方法得到的是信号能量在时频域的分布，主要包括Cohen类时频分布和仿射类时频分布。常见的Wigner-Ville分布(WVD)即为一种Cohen类时频分布，相对于线性时频分析方法来说，WVD以自身信号为窗，实现了自适应的匹配性，不受测不准原理限制，具有较高的时频分辨率。WVD分解中存在交叉项，对频谱分解结果产生了严重干扰，而应用平滑函数对WVD进行二维低通滤波能够减小交叉干扰项。在时频域的平滑程度越高，其干扰项的影响越小，时频谱的分辨率就越低，需要在联合时频分辨率和相干项影响之间折衷。

基于匹配追踪(MP)的地震谱分解方法是一种自适应时频分析方法，首先构造一个过完备的子波库，然后基于匹配追踪实现信号的稀疏分解，最后利用 WVD求取分解后的各子成分信号的时频分布，并叠加求和。基于匹配追踪的地震谱分解结果具有很高的时频分辨率，并且避免了交叉项的影响，因此在地震资料解释与反演领域得到了广泛应用。匹配追踪算法中常用的小波原子包括 Gabor小波，Morlet小波和雷克子波。基于雷克子波的匹配追踪算法通常选取零相位的雷克子波建立子波库，将匹配子波的振幅谱和相位谱作为一个整体来进行计算，这种方法得到的相位谱是由于时间延迟引起的线性相位改变。然而，零相位雷克子波库相对于实际地震子波来说冗余度较差，需要构建变相位雷克子波。同时，在实际地震记录中不但存在由时间延迟造成的线性相位改变，还存在由薄层、速度频散和能量衰减等造成的相位旋转，而常规匹配追踪(不论是零相位或是非零相位雷克子波，实数域或复数域匹配追踪)都没有单独考虑旋转相位谱。

发明内容

为了克服上述现有技术的不足，本发明提供了一种基于变相位雷克子波匹配追踪的理想地震谱分解方法，不但具有很高的时频分辨率，还能够得到旋转相位谱，而旋转相位谱中包含了更为丰富的地质信息，预计能够在储层预测和油气检测中发挥重要作用。

本发明解决其技术问题所采用的技术方案是：构造一种基于变相位雷克子波匹配追踪的理想地震谱分解方法，包含如下步骤：

(1)输入一道地震数据，当作初始残差信号；

(2)构建当前残差信号的复分析信号，然后计算复分析信号的瞬时包络，根据所述瞬时包络确定匹配子波原子的延迟时间；

(3)根据所述匹配子波原子的延迟时间构建零相位雷克子波，对零相位雷克子波的主频和相位进行扫描，构建变相位雷克子波库；

(4)从变相位雷克子波库提取一个与当前残差信号最匹配的子波原子，求取该最匹配的子波原子对应的振幅系数，将所述振幅系数与所述最匹配的子波原子相乘作为分解得到的子信号，然后从当前残差信号中减去此次分解得到的子信号，将得到的残差当作新的残差信号；

(5)重复步骤(2)到(4)直到迭代次数到达预设值或者剩余残差能量小于某一阈值为止；

(6)迭代结束后，得到若干变相位雷克匹配子波及对应的振幅系数，分别将匹配子波和振幅系数对应相乘，得到单界面反射地震信号，作为子信号集；

(7)对所述子信号集中每一个子信号做傅里叶变换，并将得到的傅里叶变换谱对应到界面对应的延迟时间上，最终得到理想的时频谱分解结果。

进一步地，在本发明的基于变相位雷克子波匹配追踪的理想地震谱分解方法中，步骤(2)具体包括：

构建残差信号s(t)的复分析信号：

S(t)＝s(t)+is_H(t)，

其中s_H(t)表示实信号s(t)的Hilbert变换，复分析信号的瞬时振幅包络为：

取瞬时振幅包络最大值所对应的时间为匹配子波原子的延迟时间。

进一步地，在本发明的基于变相位雷克子波匹配追踪的理想地震谱分解方法中，步骤(3)具体包括：

(31)根据匹配子波原子的延迟时间构建零相位雷克子波：

其中，τ为匹配子波原子的延迟时间，f_d为雷克子波的主频；

(32)通过扫描时频参数γ＝{f_d,φ}构建变相位雷克子波原子库，子波主频的扫描范围为0到奈奎斯特频率，相位的扫描范围为0°到180°，所构建的变相位雷克子波库为：

w_γ(t)＝w(t-τ,f_d)cos(φ)-w_H(t-τ,f_d)sin(φ)，

其中，φ为雷克子波的相位，w_H(t)是零相位雷克子波w(t)的Hilbert变换。

进一步地，在本发明的基于变相位雷克子波匹配追踪的理想地震谱分解方法中，步骤(4)具体包括：

对于第任意次迭代次数k，将变相位雷克子波库中的所有子波原子与当前残差信号s_res^k-1(t)做相关计算，得出相关系数corr：

其中，表示当前残差信号与子波原子w_γ(t)的内积，表示子波原子w_γ(t)的信号归一化算子；

然后，求取相关系数的最大值及其对应的子波原子w_γk(t)，

该子波原子w_γk(t)即为与当前残差信号最匹配的子波原子，然后计算得到该子波原子w_γk(t)对应的振幅系数C_k：

第k次迭代分解出来的子信号即为C_kw_γk(t)，从当前残差中减去第k次分解出来的子信号，得到新的残差信号s_res^k(t)：

s_res^k(t)＝s_res^k-1(t)-C_kw_γk(t)。

进一步地，在本发明的基于变相位雷克子波匹配追踪的理想地震谱分解方法中，步骤(6)具体包括：

将s_k(t)＝C_kw_γk(t)作为第k次迭代生成的子信号，将所有子信号构成了一个子信号集。

进一步地，在本发明的基于变相位雷克子波匹配追踪的理想地震谱分解方法中，步骤(7)具体包括：

求取子信号s_k(t)的傅里叶变换谱S_k(f)：

其中，f表示傅里叶变换的频率；

然后傅里叶变换谱S_k(f)对应到界面对应的延迟时间上，得到理想时频谱>s(t,f)：

得到界面的反射地震波的振幅谱和相位谱，

进一步地，在本发明的基于变相位雷克子波匹配追踪的理想地震谱分解方法中，理想时频谱IDFT_s(t,f)的振幅谱A(t,f)和相位谱Φ(t,f)分别为：

进一步地，在本发明的基于变相位雷克子波匹配追踪的地震谱分解方法中，步骤(7)中，还包括得到理想时频谱IDFT_s(t,f)的旋转相位谱Θ(t,f)：

其中φ_k(f)是传播时间为τ_k的最佳匹配雷克子波的相位。

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

振幅谱是单界面固定传播时间子信号的傅里叶变换组成的，在时间上来说反映了单个界面的反射子波的能量，而且匹配追踪是一种稀疏分解方法，因此该振幅谱具有理想的时间分辨率；在频率上来说，直接对与地震记录长度相同的子信号做傅里叶变换，不会因为窗口长度的影响降低其频率分辨率，因此具有理想的频率分辨率。

相位谱中包含了两种相位信息，第一种是时间延迟造成的线性相位改变，第二种是薄层、速度频散和能量衰减造成的相位旋转。常规的基于零相位雷克子波的匹配追踪得到的相位谱只反映了时间延迟造成的线性相位改变，基于非零相位雷克子波的匹配追踪得到的相位谱反映了线性相位改变和相位旋转两种影响的综合作用。时间延迟造成的线性相位改变与岩性和储层没有直接关系，因此对于常规时频分析，通常只利用振幅谱信息，较少利用相位谱信息。事实上，在地震勘探中，由于旋转相位中包含了薄层、速度频散和能量衰减的信息，因而，旋转相位谱能够在地层识别和储层预测中发挥重要作用。

附图说明

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

图1为本发明基于变相位雷克子波匹配追踪的理想地震谱分解方法流程图；

图2为实施例1模型地震道；

图3为实施例1模型地震道基于变相位雷克子波匹配追踪分解的子信号集；

图4(a)为实施例1模型地震道利用本发明基于变相位雷克子波匹配追踪的理想地震谱分解方法得到的理想地震谱分解结果中的振幅谱；

图4(b)为理想振幅谱的曲线表示；

图4(c)为理想旋转相位谱；

图5(a)、5(b)为实施例1模型地震道S变换时频分析方法得到的振幅谱和相位谱；

图6为实施例1模型地震道匹配追踪-WVD时频分解结果。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

参考图1，图1为本发明基于变相位雷克子波匹配追踪的理想地震谱分解方法的具体实施例的流程图，本方法包括以下步骤：

S1、输入单道地震数据，当作初始残差信号。

S2、构建当前残差信号s(t)的复分析信号：

S(t)＝s(t)+is_H(t)

其中，s_H(t)表示实信号s(t)的Hilbert变换，则复分析信号的瞬时振幅A(t)>

取瞬时振幅包络最大值所对应的时间为匹配子波原子的延迟时间τ。

S3、根据步骤S2中计算得到的匹配子波原子的延迟时间τ，构建零相位雷克子波

w(t)＝{1-2[πf_d(t-τ)]²}exp{-[πf_d(t-τ)]²}

其中，匹配子波原子的延迟时间τ即为雷克子波在地下介质的传播时间， f_d为匹配雷克子波的主频。为了得到任意常相位的雷克子波，求取零相位雷克子波的解析信号表达式

W(t)＝w(t)+iw_H(t)

其中w_H(t)是零相位雷克子波的Hilbert变换。将其相位旋转φ后，新的非零相位复数域雷克子波为

W'(t)＝W(t)exp(iφ)

即

W'(t)＝[w(t)+iw_H(t)][cos(φ)+isin(φ)]

则常相位雷克子波时频原子可以表示为

w_γ(t)＝w(t-τ,f_d)cos(φ)-w_H(t-τ,f_d)sin(φ)

如果信号的采样频率为f_s，根据奈奎斯特定律，雷克子波原子的最大频率为f_s/2，实际上，地震波的主频通常在100Hz以内；雷克子波的相位区间为>d,φ}即可得到雷克子波库。扫描的间隔越小，生成的小波原子越多，所得到的分解结果越稀疏，但是同时也会增加时频分解所需要的计算量。

本发明中，假设雷克子波的相位为φ，通过扫描时频参数γ＝{f_d,φ}对零相位雷克子波的主频和相位进行扫描，构建变相位雷克子波原子库，其中子波主频的扫描范围为0到奈奎斯特频率，相位的扫描范围为0°到180°，所构建的变相位雷克子波库w_γ(t)为：

w_γ(t)＝w(t-τ,f_d)cos(φ)-w_H(t-τ,f_d)sin(φ)

其中w_H(t)是零相位雷克子波w(t)的Hilbert变换。通常，若信号的采样频率为f_s，频率扫描间隔可以取傅里叶变换中的最小频率间隔，即△f＝f_s/N，其中N是数据长度。相位角的扫描间隔可以取π/36弧度。扫描的间隔越小，生成的小波原子越多，所得到的分解结果越稀疏，但是计算需要的时间就越多，需要在计算精度和计算时间两方面进行折衷。

S4、将所有雷克子波库中的子波与当前残差信号s_res^k-1(t)做相关计算，其相关系数为

求取相关系数的最大值及其对应的子波原子

其中，D为匹配雷克子波原子库，该子波原子w_γk(t)即为与当前残差信号最匹配的子波原子，其对应的振幅系数为：

第k次迭代分解出来的信号即为从当前残差中减去第k次分解出来的子信号，得到新的残差信号：

S5、对于新的残差信号，重复步骤S2到S4，直到残差信号能量小于某一阈值或者迭代次数超过设定的最大迭代次数时，迭代结束。

S6、将分解出来的匹配雷克子波与对应的振幅系数相乘，作为子信号：

则原地震信号可以利用匹配追踪分解为若干子信号以及分解后的余项：

其中，K为所进行的总的迭代次数，R_K(t)是匹配追踪分解余项，可忽略不计，将原地震信号近似写成：

所有子信号都可以看作是单界面的反射信号，实际地震记录则可以看作是若干反射界面反射波的叠加。

S7、对于所有分解得到的子信号，利用傅里叶变换，求取其振幅谱和相位谱。对于任一子信号，其传播时间τ_k是根据解析信号瞬时振幅极值计算得到的，因此所有子信号在地下介质的传播时间是已知的。常规的基于匹配追踪的时频分解方法通常是对分解出来的子信号做WVD，得到高分辨率时频谱。事实上，根据叠加原理，地震记录可以看作是若干单界面反射波的叠加，因此，将分解得到的子波做傅里叶变换：

并且将傅里叶变换结果S_k(f)放置在界面对应的传播时间，可以得到理想的时频谱IDFT_s(t,f)：

理想时频谱是一个复数矩阵，将其元素写成指数形式，

其振幅谱A(t,f)、相位谱Φ(t,f)和旋转相位谱Θ(t,f)分别是：

在本发明的一个实施例中，对一个模型地震道(图2)进行匹配追踪(图 3)，然后计算其理想振幅谱和旋转相位谱(见图4(a)-(c))。将基于本发明的频谱分解结果(图4)与其它时频分析方法的频谱分解结果(图5(a)、图 5(b)、图6)相对比，可见：S变换的时频分辨率相对较差，其相位谱中包括了线性相位和旋转相位，实用性较差；基于匹配追踪的WVD分解结果中振幅谱表示的是信号的能量而非振幅，不存在相位谱；本发明的谱分解结果具有更高的时间和频谱分辨率，并且能够得到旋转相位信息，具有较高的应用价值。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。