一种基于匹配追踪的Wigner高阶谱地震信号谱分解方法
【技术领域】
[0001] 本发明属于非平稳信号时频分析及地震信号处理领域，具体涉及一种基于匹配追 踪的Wigner高阶谱地震信号谱分解方法。
【背景技术】
[0002] 谱分解是一种地震信号解释技术，它将地震数据分解到时频域，由此揭示出时频 域中包含有用的含油气信息。多项研究都证明了该方法在储层解释和预测方面的实用性， 例如进行储层厚度估计，地层解释以及流体识别等。传统的谱分解方法，如短时傅里叶变 换，它通过引入窗函数来计算信号的时频分布，故而其时频分辨率被相应的窗函数所限制， 不能满足高精度地震勘探对精细储层预测的要求。Ville (1948)将Wigner分布（1932)引 入了信号处理领域，提出著名的 Wigner-Ville 分布（Wigner-Ville Distribution，WVD)， 它在时频平面上直接定义能量密度，不需要受分辨率的限制。业已普遍承认，没有任何一种 时频联合分布的时频分辨率能出其右，然而在实际应用中，该变换的二次性产生了干扰，即 WVD的交叉项问题。
[0003] 在众多谱分解方法中，由Mallat和Zhang(1993)提出的匹配追踪（Matching pursuit，MP)算法具有很高的时频分辨率，得到了广泛的应用。MP是一种灵活的自适应信 号分解方法，它结合信号的先验信息构造合适的冗余字典，将信号在该字典上分解，自适应 地得到最能够匹配信号实际结构的分解表达式。同时，由于字典中的原子具有很好的时频 聚集性，将其与信号的WVD相结合，一方面可以利用WVD的最佳时频分辨率提高MP时频谱 的分辨率，另一方面可以解决WVD的交叉项问题。
[0004] Chakraboty和0kaya(1995)首先将MP算法引入地震信号谱分解中，取得了良好的 效果，自此，MP被广泛地应用于地震信号低频阴影检测，时频属性提取等。但MP是一个密集 计算过程，计算效率低，并且其结果依赖于字典的构建。Liu等（2004, 2005)引入Morlet小 波和Ricker子波来构建字典，并将复数道分析加入地震信号的MP分解中，提出利用Morlet 小波的动态MP，提高了 MP的效率。Wang (2007)给出了完整的基于复数域动态最优搜索的MP 分解方法，并导出一种快速计算表达式。此后，该方法被广泛应用于地震信号时频分析中， 并得到了一定程度的发展和改进。张繁昌等（2010,2013)提出了利用动态子波库的双参数 和单参数扫描算法，使算法效率得到大幅度提尚。黄桿东等（2012)在Morlet小波中引入 能量衰减因子，改进了算法灵活性和重构精度。Zhao等（2012)提出将字典扩展为Ricker 子波、Morlet小波和多相位地震子波的集合，表明多成分字典更能够有效地反映地震信号 中包含的信息。
[0005] 随着信号处理技术的发展，信号高阶统计量具有一些重要的性质，使得人们对信 号高阶统计量的研究逐渐重视起来，许多学者从高阶统计量的角度来讨论WVD，即高阶统计 量与WVD相结合的一种产物--Wigner高阶谱时频分布。Gerr (1988)首先提出三阶Wigner 分布。Swami (1992)对时变高阶谱的定义、性质及应用作了进一步研究，取得了较大的进展。 Wigner高阶谱是WVD在高阶谱域的延伸，既继承了时频分布同时反映信号频谱成分和时间 之间关系的优点，还引入了高阶谱对高斯噪声的良好抑制能力，同时可以获得时频聚集性 更好的时频谱。然而，Wigner高阶谱与WVD-样同样存在交叉干扰项的问题。并且随着阶 数的增加，计算量大大提升，因此Wigner高阶谱的计算主要是Wigner双谱和Wigner三谱。
【发明内容】
[0006] 本发明提供了一种基于匹配追踪的Wigner高阶谱地震信号谱分解方法，旨在去 除Wigner高阶矩谱的交叉项，并获得高时频聚集性的谱分解结果，提高储层预测的精度。 具体的技术方案如下所述。
[0007] 本发明为了解决上述技术问题，采用以下技术方案：
[0008] 本发明中所述方法是为了克服上述现有技术的缺点，主要针对去除Wigner高阶 矩谱的交叉项，并获得高时频聚集性的谱分解结果，提高储层预测的精度的问题，提出了一 种基于匹配追踪的Wigner高阶谱地震信号谱分解方法，包括以下步骤：
[0009] 步骤1 :读入地震剖面的q道数据，选定原子类型；
[0010] 步骤2 :设置初始取变量i = 1 ;
[0011] 步骤3 :读取第i道地震数据Xl (t)，设置匹配追踪分解次数N，将第η次分解的信 号记为Rn⑴，用Xi⑴对札⑴赋值，即札(t) = Xi (t)，n e [1，2,…，Ν];
[0012] 步骤4 :，设置初始取变量η = 1 ;
[0013] 步骤5 :对分解信号Rn(t)进行复数道分析，确定Rn(t)对应原子的初始时间延迟 u。、初始频率ω。和初始相位1%三个参数；
[0014] 步骤6 :对尺度因子进行全局搜索，确定分解信号Rn(t)对应原子的初始尺度因子 σ。，由此得到初始参数集％: = 叫，％丨；
[0015] 步骤7 :对初始参数集Ad%，丨进行局部搜索，找到与信号Rn(t)最匹配 的原子f);
[0016] 步骤8 :计算信号Rn(t)的最匹配原子S';.,.⑴的k阶Wigner高阶谱的对角切片谱 SWkgJt,f)；
[0017] 步骤9 :计算信号比(〇在最匹配原子^⑴方向上投影后的残差，并将残差视为新 的信号Rn+1 (t)，令n = n+1，重复步骤5-8,直到达到最大分解次数N ;
[0018] 步骤10 :对η所有取值下的信号Rn⑴的最匹配原子:仏,⑴的k阶Wigner高阶谱 的对角切片谱求和，作为该道地震数据的Wigner高阶谱时频谱截取 单频切片SFliV;
[0019] 步骤11 :令变量i = i+1，重复步骤3-10,直到取完整个地震剖面的数据，即直到i =q，由所有道地震数据的不同单频切片可以构成整个二维剖面的一系列单频属性，即谱分 解的结果SFV= {SFliV，1彡i彡q}。
[0020] 上述技术方案中，所述步骤5包括以下几个步骤：
[0021] 步骤 5. 1 :对 Rn(t)做 Hilbert 变换：
[0023] 其中，t是时间变量，P表示柯西主值，τ表示时间积分变量。
[0024] 步骤5. 2 :令X (t) = Rn⑴，确定Rn⑴的解析信号z⑴：
[0026] 其中，j是虚数单位，e是自然常数，a表示瞬时振幅；
[0027] 步骤5. 3 :确定瞬时振幅a (t):
[0029] 步骤5. 4 :选取瞬时振幅a (t)达到最大的时刻u。作为原子的时间延迟初始值：
[0030] iir, =argmaxi7(/)
[0031] 步骤5. 5 :确定原子的初始相位即该时刻的瞬时相位炉:
[0034] 步骤5. 6 :确定原子的初始频率ω。，即该时刻的瞬时频率ω (u。）：
[0037] 上述技术方案中，所述步骤6包括以下几个步骤：
[0038] 步骤6. 1 :设置尺度因子〇的全局搜索范围[σ σ 2]
[0039] 步骤6. 2 :根据尺度因子〇的数值计算原子：
L0041」 具中，gY⑴是通迓对基本汲形g⑴做Ζ 树雙化而得到的原子波形表 达式，t是时间变量，j是虚数单位，e是自然常数；u为时间延迟，控制基本波形的时间中心； ω为频率调制因子，控制基本波形的频率中心；识为相位调制，控制波形的形状；〇为尺度 因子，控制其在时间域的跨度；
[0042] 步骤6. 3 :计算信号Rn(t)与原子⑴的内积<Rn(t), (t)> :
[0044] 步骤6.4:找到使〈&(〇4"〇>达到最大值的尺度因子〇作为原子的初始尺度 因子〇。：
[0046] 步骤6. 5 :得到初始参数集γ。：
[0048] 其中，γ表示四个参数的集合用以描述一个原子。
[0049] 上述技术方案中，所述步骤7包括以下几个步骤：
[0050] 步骤7. 1 :设置四参数集γ ξ的局部搜索范围为[γ (J-Δ γ, yQ+A γ];
[0051] 步骤7. 2 :找到与信号Rn(t)最匹配的原子4,.(0:
[0053] 其中，Α ω是对应γ ξ参数集的原子，γ n表示第η次分解信号下的最佳参数集。
[0054] 上述技术方案中，所述步骤8包括以下步骤：
[0055] 步骤8. 1 :计算信号Rn(t)的最匹配原子及...⑴的k阶Wigner高阶谱
[0057] 其中，t是时间变量，f是频率变量，j是虚数单位，e是自然常数，τ表示时间积 分变量，r、m和s表示整数变量。
[0058] 步骤8. 2 :计算信号Rn(t)的最匹配原子⑴的k阶Wigner高阶谱的对角切片谱
[0060] 上述技术方案中，所述步骤9中计算信号Rn(t)在最匹配原子方向上投影后的残 差，并将残差视为新的信号Rn+1 (t)，包括以下步骤：
[0061] 步骤9. 1 :计算信号Rn(t)在最匹配原子方向上投影后的残差，并将残差视为新的 十目号Rn+i⑴：
[0063] 上述技术方案中，所述步骤10包括以下步骤：
[0064] 步骤10. 1 :对第i道地震数据的所有原子的k阶Wigner高阶谱的对角切片谱求 和^ (〖,/):
[0066] 步骤10. 2 :将I (?,/)作为第i道地震数据的Wigner高阶谱时频谱；
[0067] 步骤10. 3 :截取单频切片SFliV:
[0069] 其中，v表示单频切片的频率值。
[0070] 本发明是一种基于匹配追踪的Wigner高阶谱地震信号谱分解方法，并将该模型 应用于非平稳信号时频分析及地震信号处理领域。
[0071] 综上所述，由于采用了上述技术方案，本发明的有益效果是：
[0072] 本发明将基于匹配追踪的Wigner-Ville时频分布扩展到基于匹配追踪的Wigner 高阶谱，利用匹配追踪方法去除Wigner高阶谱的交叉项，可以获得时频聚集性更高的地震 谱分解结果，为后续地震储层预测和流体识别提供更精确的信息，具有较高的实用性。
【附图说明】
[0073] 图1为方法流程图。
[0074] 图2为一个原始地震剖面；
[0075] 图3为一个Morlet原子；
[0076] 图4为第15道地震数据；
[0077] 图5为第一次分解的最佳匹配原子的波形；
[0078] 图6为第一次分解的最佳匹配原子的时频谱，6a为Wigner-Ville时频谱，6b为 Wigner双谱（k = 2)对角切片谱，6c为Wigner三谱（k = 3)对角切片谱；
[0079] 图7为第15道地震数据的Wigner高阶谱时频谱，7a为Wigner-Ville时频谱，7b 为Wigner双谱（k = 2)对角切片谱，7c为Wigner三谱（k = 3)对角切片谱；
[0080] 图8为地震数据的谱分解结果，8a为Wigner-Ville时频谱得到的45Hz单频剖面， 8b为Wigner双谱（k = 2)对角切片谱得到的45Hz单频剖面，8c为Wigner三谱（k = 3) 对角切片谱得到的45Hz单频剖面。
【具体实施方式】
[0081] 本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥 的特征和/或步骤以外，均可以以任何方式组合。
[0082] 下