具体实施方式

为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

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

本发明提出的一种基于混合范数和小波变换的稀疏时频谱分解方法，首先提出一种完全解析、能够匹配地震子波的母小波，然后根据稀疏表示和混合范数的思想，提出一种具有局域化的时频谱分解方法，最后将该时频谱分解方法用于水合物探测，预测游离气的储层位置。

本发明提出的一种基于数据驱动的稀疏时频谱分解方法包括以下步骤：

步骤1)、获得叠后观测数据

采集原始地震资料，并对原始地震资料进行预处理，得到叠后观测数据，记为其中N是时间采样点数，n∈[1,N]表示当前是第n个采样点。

步骤2)、获得局域化的时频谱分解方法，采用步骤1)获得的叠后观测数据构建标架算子F、叠后观测数据y(n)的小波变换和反变换：

首先，构造小波变换的母小波，母小波在频率域定义如公式(1)所示：

ψ(f)＝U(f)F(f)^α(1-F(f))^β (1)

其中，U(f)为单位阶跃信号；为构造的母小波基函数；f∈[0,f_c]为频率；f_c是单位阶跃信号的截止频率；α,β是调整母小波形态的调节参数，使得母小波更加匹配地震子波，获得高时频局域化的时频谱分解。

其次，在确定母小波的α,β调节参数后，将小波变换构建为紧标架的表示形式，则已知叠后观测数据y(n)小波变换和反变换的表达式如公式(2)和公式(3)所示：

x＝F^*y, (2)

y＝Fx, (3)

其中，x＝[x_j,k],j＝1,2,...,J；k＝1,2,...,K是小波变换的系数，J表示小波变换的尺度采样长度，K表示小波变换时间采样长度；y是由y(n)生成的列向量；F是由母小波ψ(f)生成标架算子，F^*为F的伴随算子，FF*为单位矩阵，已知F的情况下，可以计算F^*，x_j,k表示为小波系数中第(j,k)个系数。

步骤3)、根据步骤2)得到的标价算子F、小波变换与反变换构造基于混合范数的稀疏时频谱分解模型：

首先，根据标架理论，在已知标架算子F和地震信号y时，将求解小波变换的系数x表示为一个带约束的反问题求解，稀疏模型如公式(4)所示：

其中，λ为规则化参数。是惩罚函数。惩罚函数的类型不同，可以获得不同类型的时频谱系数。

为了获得具有更高时频局域化的时频谱系数，在上述稀疏模型中引入混合范数，该混合范数由非凸的稀疏约束和L2范数组成，则新的稀疏时频谱分解模型如公式(5)所示：

式中，y表示由采集到的叠后观测数据生成的列向量；x＝[x_j,k],j＝1,2,...,J；k＝1,2,...,K；x是小波变换的系数，J表示小波变换的尺度采样长度，K表示小波变换时间采样长度；λ_j和λ₂分别表示正则化参数，j表示尺度，k表示时间；||x||₂表示L2范数正则化项，用来防止时频谱系数过于稀疏。φ(x_j,k,a_j)是由已知变量a_j确定的稀疏约束的惩罚函数。与传统的L1范数稀疏约束对比，该惩罚函数是非凸的惩罚函数，能够避免由L1范数不是最稀疏约束而导致的问题，也可以在变量a_j满足一定条件时，确保上述优化问题为凸优化，利用凸优化可以求解上述优化问题。该非凸惩罚函数的种类有很多种，本专利使用类似于反正切形式的非凸惩罚函数，其定义如公式(6)所示：

当公式(6)中的变量a_j满足时，上述稀疏时频谱分解模型(5)是凸的，利用凸优化的理论可以很容易求解该稀疏时频谱分解模型问题。

步骤4)、利用split Bregman迭代算法求解步骤3)中稀疏时频谱分解模型获得具有时频局域化的时频谱系数。

首先，确定稀疏正则化参数λ₁和λ₂以及初始值x₀。引入中间变量u，则上述公式(5)变为公式(7)所示的结果：

式中，u＝[u_j,k]j＝1,2,...,J；k＝1,2,...,K为引进的中间变量，该变量具有与小波系数x想同的维度，u_j,k表示中间变量u的第(j,k)个元素。

然后，根据凸优化理论，将有约束优化问题变为无约束优化问题：

根据变量分割原理，将上述无约束优化问题(8)可以变为两个子优化问题如公式(9)所示：

式中，k表示迭代次数；u＝[u_j,k]j＝1,2,...,J；k＝1,2,...,K为引进的中间变量；b＝[b_j,k]j＝1,2,...,J；k＝1,2,...,K也表示引进的中间变量；J表示小波变换的尺度采样长度，K表示小波变换时间采样长度；

针对上述公式(9)中的两个子优化问题，分别求解子优化问题，通过两个优化问题之间交替迭代，求解最终最优的解。

一般有两种方式停止迭代：第一是达到最大迭代次数；第二是当当前的迭代结果和上一次迭代结果的误差小于某个门限值的时候，本申请使用的是达到最大迭代次数。

数值仿真结果-合成地震记录数据：

首先，利用一个无噪的合成地震信号来验证本发明的有效性，如图1(a)所示。时间采样间隔为1ms，时间采样次数为512。第一个子波是在0.05s处添加一个主频为60Hz的Ricker子波生成的。第二个子波是在0.15s处给一个负反射系数褶积一个主频为50hz的Ricker子波而生成的。第三个子波由是由0..25s的正40hz Ricker小波和0.275s处的负40hz Ricker子波组成。最后一个子波由3个具有相同主频的Ricker子波组成，主要包括两个30Hz的正Ricker子波(0.35和0.41s)和一个30Hz的负Ricker子波(0.38s)。本实验分别对比了小波变换、短时傅里叶变换、挤压小波变换、基于L1约束的时频谱分解方法和本发明提出的时频谱分解方法。设基于L1约束的时频谱分解方法和本发明提出的时频谱分解方法的最大迭代次数为50。图1(b)-图(f)为这几个时频变换方法的时频谱结果。

如图1(b)和图1(c)所示，与短时傅里叶变换对比，小波变换具有更好的时频局域化。但是，这两个变换都受不确定原理的限制，时频局域化受到限制。挤压小波变换如图1(d)所示虽然具有较好的频率分辨率，但没有考虑时间分辨率。与图1(e)中基于L1约束的时频谱分解结果对比，本发明提出的时频谱分解方法如图1(f)所示具有更稀疏的时频谱分解结果。

其次，图2给出了含噪合成地震时频谱分解结果。图2(a)是给图1(a)添加高斯白噪声生成的，其SNR为10dB。由于传统的时频谱分解方法抗噪性都比较差，所以此实验只验证了小波变换、基于L1约束的稀疏时频谱分解方法和本发明提出的时频谱分解方法，其结果如图2(b)-图2(d)。显然，小波变换的抗噪性能比较差，而基于反演的两种方法具有较好的抗噪性能。与基于L1约束的稀疏时频谱分解方法对比，本发明提出的时频谱分解方法具有更好的抗噪性和稀疏性。

实际地震资料剖面：利用三维的水合物地震数据来进一步验证本发明的有效性。图3显示了三维地震数据，其中Inline编号为95,Xline编号为1306。每个轨迹的采样时间间隔为2ms。如图4所示，在地震剖面Inline 11中，由于数据质量高，可以很容易地追踪到BSR，即红线。由于游离气具有衰减特性，因此利用本发明提出的时频谱分解方法可以检测游离气储层的位置。利用本发明提出的时频谱分解方法计算的沿BSR的衰减剖面结果如图5所示，黑色为衰减严重的地方，即为游离气储层的位置。显然，利用本发明提出的时频谱分解方法可以有效的检测游离气的位置，为水合物储层检测提供一种有效的检测方法。

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