本发明涉及石油天然气资源地球物理勘查技术领域，特别是涉及到一种基于稀疏表示的道集调谐优化增维方法。
背景技术：
我国陆相沉积储层具有厚度相对较薄、空间变化大的特点，会导致地震波在其中的传播发生薄互层的干涉效应。为提高薄互层储层的识别精度，最为重要的解决方法就是拓宽地震资料频带、提高数据主频。针对成像的叠后地震数据，目前常用的处理手段有反q滤波、拓频、分频等拓宽频带和提高主频的方法。反q滤波针对地震波的传播随走时的增大，高频能量吸收衰减和相位畸变的问题，补偿地震波的振幅和相位，从而提高地震资料的主频；拓频方法则是将地震信号的波长压缩，从而达到拓宽地震资料的频带和提高地震数据主频的效果；分频方法通过频谱分解，提取高频信息来实现提高地震资料频率的效果。这些方法对于提高薄互层的识别精度都发挥了一定的作用，但在资料的保真性方面还有待深究。
稀疏性在广泛的信号处理领域已经成为重要的概念之一，无论是在理论上还是实践中，稀疏性始终是一种很有吸引力的信号特征。在地震信号处理领域，稀疏性可以作为先验信息来促进地震数据去噪、反褶积、地震数据重建等线性反演问题的求解。由于自然界中大部分信号通常并不是稀疏的，需要对其进行稀疏表示。稀疏表示就是通过字典中很少量元素的线性组合的形式来表示信号。作为信号的一种特征，无论是地震反射系数序列自身的稀疏性，还是地震数据在过完备字典中的稀疏性，都能够在保证有效信息的前提下降低维度，并有效地促进地震反演问题的求解。从频率域角度来看，数据空间维度要小于模型空间维度，因而地震反演具有多解性。在特定条件下，稀疏约束反问题的解具有唯一、稳定的特征，可以自动选择最优解。
现有的叠前反演道集处理方法主要存在以下问题：①由于叠前道集是高维数据体，数据量大、冗余度高、逐一解释效率低，因而通常采取的方法是利用部分角度叠加的策略，但是叠加会产生干涉效应降低地震分辨率；②传统的叠前反演没有考虑到频率对地震属性的影响，与此同时，传统的时频表征方法也因为时频分辨率不高，不利于对薄互层的精细刻画。为此我们发明了一种新的基于稀疏表示的道集调谐优化增维方法，解决了以上技术问题。
技术实现要素：
本发明的目的是提供一种避免了因为叠加效应而产生的地震波干涉现象，消除了由于频率混叠而造成的失真影响，能够有效地提高后期的反演刻画精度的基于稀疏表示的道集调谐优化增维方法。
本发明的目的可通过如下技术措施来实现：基于稀疏表示的道集调谐优化增维方法，该基于稀疏表示的道集调谐优化增维方法包括：步骤1，构建时频域子波字典；步骤2，进行道集hilbert变换；步骤3，确定匹配子波控制参数；步骤4，进行基于匹配追踪的子波分解；步骤5，进行基于匹配子波的高分辨时频分解。
本发明的目的还可通过如下技术措施来实现：
在步骤1中，在分析地震信号时选择构建ricker子波字典：
其中，ωr(t,fi)表示ricker子波字典，t代表ricker子波的中心时间，fi代表子波字典中第i个ricker子波的主频。
ricker子波具有三个控制参量，旁瓣能量的收敛速度可变,延迟时间可控以及子波库波形丰富，与地震子波匹配较为灵活。
在步骤2中，对于非平稳信号来说，它的瞬时频率为：
其中，
道集内的每一道信号x(t)hilbert变换为：
式中，p为柯西主值，τ表示时移。
因此得到x(t)的解析信号为：
在步骤3中，将匹配子波的三个控制参数，即时移τ，频率因子f，相位因子
在步骤3中，确定了最佳匹配子波ω的中心时间和主频后，根据残差能量最小原则，利用阻尼最小平方法计算，除去这些匹配子波的振幅和相位信息；
假定存在一道带宽有限的地震信号x(t)，通过匹配追踪算法迭代了k-1次，残差信号为rk-1x；当迭代到第k次时，共产生m个匹配子波ωδm(t)，m＝1,2,…,m；对残差信号rk-1x进行hilbert变换构造解析信号：
x＝rk-1x irk-lx(5)
同理构造第k次迭代时的复匹配子波：
wδm(t)＝ωδm(t) iωδm(t)(6)
所以经过k次迭代后产生的残差信号rkx的能量为：
式中x为地震信号的n未列向量，a代表匹配子波的复振幅：
am为实数域匹配子波的振幅，
gδm(t)＝[ωδ(t,f1),ωδ(t,f2),ωδ(t,f3),…,ωδ(t,fm)]为n×m阶匹配子波矩阵，每一列对应一个匹配子波。
为使得地震信号x(t)与匹配子波集ωδm(t)间的相对误差最小，采用阻尼最小平方法计算：
a＝[wtw εi]-1wtx(9)
式中，w为k次迭代后得到的n×m阶复匹配子波矩阵，ε为阻尼因子，i为m×m阶单位矩阵；由式(8)和(9)便可确定匹配子波振幅和相位的控制参数ak和
在步骤4中，匹配子波通过对基本子波ω(t)进行时移、调制和相位变化形成；若基本子波ω(t)为匹配子波库中的满足条件的一个基本子波，令
式中，τ表示时移，f表示为频率因子，
在步骤4中，构建好过完备匹配子波库φ后，尽管φ中的匹配子波并不是全部正交，但仍需满足||ωi||2＝1；从匹配子波库φ中选出与给定信号x(t)最为匹配的ω1，最匹配的要求则是他们的内积<x,ω1>为匹配子波库中所有原子与x(t)内积最大的一个：
<x,ω1>＞<x,ωi>i≠1ωi∈φ(11)
所以，一次迭代后，信号x(t)被分解为：
x(t)＝<x,ω1>ω1 r1x(12)
r1x是第一次迭代后的残差信号；r1x与ω1正交，所以：
<x,ω1>越大，意味着ω1是信号所在hilbert空间中与
rkx＝<x,ωk-1>ωk-1 rk-1x(14)
经过k次迭代分解计算后，原始信号x(t)用k个匹配子波的合成来近似，当残差r1x足够小，即信号的稀疏表示跟信号的近似程度足够好或残差信号满足设定的阈值时，迭代停止。
在步骤5中，根据步骤3、4得到的高精度匹配子波，利用子波集内每一个匹配子波的独立性，选用需要的时频分解方法，对其分别进行时频域展开；再通过线性叠加，得到整个信号的时频谱图；重复对道集内的每一道进行之上的循环，最终得到多分辨去混叠的调谐数据体。
本发明中的基于稀疏表示的道集调谐优化增维方法，主要用于高分辨叠前反演道集处理，进而通过稳定可靠的稀疏表示方法获得高精度多分辨去混叠的调谐数据体。该方法充分利用叠前道集所蕴含的丰富地震信息，考虑道集的地震属性在随偏移距变化而改变的同时也会受到频率的影响，大角度道集信号的叠加会导致波的干涉，影响地震资料的品质，以及在采集的过程中信号会受到频率混叠的影响。对道集内的每一道依次展开，避免了叠加导致的干涉效应，并利用匹配追踪的算法对每一道进行高精度的时频展开，提升了时频分辨率，避免了频率混叠的影响，提高了对薄储层刻画的精度。本发明在反演工作开展之前的道集处理过程中，考虑到了道集的地震属性在随偏移距变化而改变的同时也会受到频率的影响。此外，大角度道集信号的叠加会导致波的干涉，影响地震资料的品质，而在采集的过程中信号会受到频率混叠的影响。为此，本方法对道集中的每一道依次展开，避免了叠加导致的干涉效应，并利用匹配追踪的算法对每一道进行高精度的时频展开，提升了时频分辨率避免了频率混叠的影响，提高了对薄储层刻的画精度。
附图说明
图1为本发明的一具体实施例中rick子波字典的示意图；
图2为本发明的一具体实施例中hilbert变换信号的示意图；
图3为本发明的一具体实施例中信号瞬时包络与瞬时频率的求取的示意图；
图4为本发明的一具体实施例中匹配追踪算法分解的子波集的示意图；
图5为本发明的一具体实施例中匹配追踪算法分解重构对比图；
图6为本发明的一具体实施例中传统时频分析方法与基于匹配追踪算法的时频分析方法时频分辨率对比图；
图7为本发明的一具体实施例中基于匹配追踪算法的时频表征方法区分薄储层的示意图；
图8为本发明的一具体实施例中某道集多分辨去混叠的调谐数据体示意图；
图9为本发明的基于稀疏表示的道集调谐优化增维方法的一具体实施例的流程图。
具体实施方式
为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。
如图9所示，图9为本发明的基于稀疏表示的道集调谐优化增维方法的流程图。
步骤101，时频域子波字典的构建。
为了使稀疏表示的结果能够更好地满足分析信号的要求，过完备原子字典需要满足一定的条件：构造的子波字典应当尽可能与信号的内在结构或性状相匹配；字典中的原子应具有较好的时频聚集性以便更好地适应非平稳信号。由于ricker子波具有波形简单，收敛较快，延迟时间短的特点，在薄互层分析中具有较高的时频分辨率，因而在此分析地震信号时选择构建ricker子波字典。
此外，ricker子波具有三个控制参量，旁瓣能量的收敛速度可变,延迟时间可控以及子波库波形丰富，与地震子波匹配较为灵活。
步骤102，道集hilbert变换。
信号的频域是描述信号的一种有效方法，但它属于全局变换，在时域是完全非局域化的。在物理和工程领域中，对于复杂信号呈周期性变化的现象，瞬时频率是最直观的描述。对于非平稳信号来说，它的瞬时频率为：
其中，
解析信号的实质是将一个窄带实信号转化为复信号，hilbert变换能够对窄带实信号产生唯一的复信号。道集内的每一道信号x(t)hilbert变换为：
式中，p为柯西主值。因此得到x(t)的解析信号为：
步骤103，匹配子波控制参数的确定。
将匹配子波的三个控制参数(时移τ，频率因子f，相位因子
利用得到的解析信号，通过时间序列确定信号采样率和信号长度。利用解析信号求取信号的瞬时包络与瞬时相位，则τk可粗略估计为瞬时包络局部极大值峰值处所对应的时间t。在复数域内对子波库中的原子进行频率参数的扫描，对于某次迭代确定出的若干τm，得到信号在这些位置处的频率先验值fm，假定匹配子波为零相位子波，利用τm和fm确定若干的先验信息点(τm,fm,0)。以这些先验信息点为中心，选取频率扫描范围，确定匹配子波的波峰频率fk。
确定了最佳匹配子波ω的中心时间和主频后，下一步来确定这些匹配子波的振幅和相位信息。根据残差能量最小原则，利用阻尼最小平方法计算，除去这些匹配子波的振幅和相位信息。
假定存在一道带宽有限的地震信号x(t)，通过匹配追踪算法迭代了k-1次，残差信号为rk-1x。当迭代到第k次时，共产生m个匹配子波ωδm(t)，m＝1,2,…,m。对残差信号rk-1x进行hilbert变换构造解析信号：
x＝rk-1x irk-1x(5)
同理构造第k次迭代时的复匹配子波：
wδm(t)＝ωδm(t) iωδm(t)(6)
所以经过k次迭代后产生的残差信号rkx的能量为：
式中a代表匹配子波的复振幅：
am为实数域匹配子波的振幅，
gδm(t)＝[ωδ(t,f1),ωδ(t,f2),ωδ(t,f3),…,ωδ(t,fm)]为n×m阶匹配子波矩阵，每一列对应一个匹配子波。
为使得地震信号x(t)与匹配子波集ωδm(t)间的相对误差最小(即残差信号能量e最小)，采用阻尼最小平方法计算：
a＝[wtw εi]-1wtx(9)
式中，w为k次迭代后得到的n×m阶复匹配子波矩阵，ε为阻尼因子(通常取0.01)，i为m×m阶单位矩阵。由式(8)和(9)便可确定匹配子波振幅和相位的控制参数ak和
步骤104，基于匹配追踪的子波分解。
过完备匹配子波库(冗余字典)是一系列匹配子波构成的子波矩阵，匹配子波通过对基本子波ω(t)进行时移、调制和相位变化形成。若ω(t)为匹配子波库中的满足条件的一个基本子波，令
式中，τ表示时移，f表示为频率因子，
构建好过完备匹配子波库φ后，尽管φ中的匹配子波并不是全部正交，但仍需满足||ωi||2＝1。从匹配子波库φ中选出与给定信号x(t)最为匹配的ω1，最匹配的要求则是他们的内积<x,ω1>为匹配子波库中所有原子与x(t)内积最大的一个：
<x,ω1>＞<x,ωi>i≠1ωi∈φ(11)
所以，一次迭代后，信号x(t)被分解为：
x(t)＝<x,ω1>ω1 r1x(12)
r1x是第一次迭代后的残差信号。r1x与ω1正交，所以：
<x,ω1>越大，意味着ω1是信号所在hilbert空间中与
rkx＝<x,ωk-1>ωk-1 rk-1x(14)
经过k次迭代分解计算后，原始信号x(t)可以用k个匹配子波的合成来近似，当残差r1x足够小，即信号的稀疏表示跟信号的近似程度足够好或残差信号满足设定的阈值时，迭代停止。
步骤105，基于匹配子波的高分辨时频分解。
根据步骤103、104得到的高精度匹配子波，利用子波集内每一个匹配子波的独立性，选用需要的时频分解方法，对其分别进行时频域展开。再通过线性叠加，得到整个信号的时频谱图。重复对道集内的每一道进行之上的循环，最终得到多分辨去混叠的调谐数据体。
在一实施例中，以某叠前共成像点道集为例来说明具体的技术方案：
第一步：时频域子波字典的构建
为了使稀疏表示的结果能够更好地满足分析信号的要求，过完备原子字典需要满足一定的条件：构造的子波字典应当尽可能地与信号的内在结构或性状相匹配；字典中的原子应具有较好的时频聚集性以便更好地适应非平稳信号。利用ricker子波具有波形简单，收敛较快，延迟时间短的特点，在薄互层分析中具有较高的时频分辨率，因而在此分析地震信号时选择构建ricker子波字典。如图1所示，截取展示了ricker子波字典中部分波形。
第二步：道集hilbert变换
频域的信号描述是一种有效方法，但它属于全局变换，在时域是完全非局域化的。在物理和工程领域中，对于复杂信号呈周期性变化的现象，瞬时频率是最直观的描述。而为了求取这些瞬时属性，需要引入解析信号。解析信号的实质是将一个窄带实信号转化为复信号，hilbert变换能够对窄带实信号产生唯一的复信号。我们对道集内的每一道信号x(t)进行hilbert变换，得到了图2所示的解析信号。
第三步：匹配子波控制参数的确定
将匹配子波的三个控制参数所生成的空间离散化，任一离散点
利用得到的解析信号，通过时间序列确定信号采样率和信号长度。利用解析信号求取信号的瞬时包络与瞬时相位，则τk可粗略估计为瞬时包络局部极大值峰值处所对应的时间t。在复数域内对子波库中的原子进行频率参数的扫描，对于某次迭代确定出的若干τm，得到信号在这些位置处的频率先验值fm，假定匹配子波为零相位子波，利用τm和fm确定若干的先验信息点(τm,fm,0)。以这些先验信息点为中心，选取频率扫描范围，确定匹配子波的波峰频率fk。
图3展示了我们通过构建的解析信号求取瞬时包络来固定时间t，在通过对这些点求取瞬时频率，得到信号在这些位置处的频率先验值fm。
确定了最佳匹配子波ω的中心时间和主频后，下一步来确定这些匹配子波的振幅和相位信息。根据残差能量最小原则，可以利用阻尼最小平方法来计算获取这些匹配子波的振幅和相位信息。
第四步：基于匹配追踪的子波分解
过完备匹配子波库(冗余字典)是一系列匹配子波构成的子波矩阵，而匹配子波实质上就是通过对基本子波ω(t)进行时移、调制和相位变化而产生的。构建好过完备匹配子波库φ后，尽管φ中的匹配子波并不是全部正交，但仍需满足||ωi||2＝1。从匹配子波库φ中选出与给定信号x(t)最为匹配的ω1，最匹配的要求则是他们的内积<x,ω1>为匹配子波库中所有原子与x(t)内积最大的一个。经过k次迭代分解计算后，原始信号x(t)可以用k个匹配子波的合成来近似，当残差r1x足够小，即信号的稀疏表示跟信号的近似程度足够好或残差信号满足设定的阈值时，迭代停止。
图4为匹配追踪算法分解的子波集，展示了由信号分解得到的各个匹配ricker子波；图5为匹配追踪算法分解重构对比图，由图可以直观地看出，利用匹配追踪算法分解后的匹配子波线性叠加得到的重构信号与原始信号基本吻合，能够保留应有的信息。
第五步：基于匹配追踪的子波分解
根据步骤(3)、(4)得到的高精度匹配子波，利用子波集内每一个匹配子波的独立性，选用需要的时频分解方法，对其分别进行时频域展开。再通过线性叠加，得到整个信号的时频谱图。重复对道集内的每一道进行之上的循环，最终得到多分辨去混叠的调谐数据体。
图6为传统时频分析方法与基于匹配追踪算法的时频分析方法时频分辨率对比图。图(a)和图(b)分别为利用短时傅里叶变换与维格纳分布所得到的时频图，可以看出stft对于窗口的选择要求很高，如果要提升时间分辨率，窗口必须要尽可能开的小，这必然会损失频率分辨率；反之，若想利用大的窗口提高频率分辨率则也势必需要牺牲时间分辨率。而wvd分析虽然有着较好的时频分辨率，不用选择窗函数，以自身信号的逆序列为分析时窗，将时窗的影响降到最低，时频分布的时长和带宽的乘积达到了最小值，但是存在着严重的交叉项干扰。图(c)和图(d)为利用匹配追踪算法优化后得到的stft和wvd时频图，时频分辨率得到了极大的提升，对频率的刻画更加细致。
图7为基于匹配追踪算法的时频表征方法区分薄储层。传统的时频分析方法无法准确的对薄互层进行识别，而基于匹配追踪算法的高分辨率时频表征手段极大地提升了这种识别能力。
图8为某道集多分辨去混叠的调谐数据体示意图。我们通过对道集内的每一道依次展开，避免了叠加导致的干涉效应，并利用匹配追踪的算法对每一道进行高精度的时频展开，提升了时频分辨率，避免了频率混叠的影响，提高了对薄储层刻的画精度。
本方法在反演工作开展之前的道集处理过程中，考虑到了道集的地震属性在随偏移距变化而改变的同时也会受到频率的影响。此外，大角度道集信号的叠加会导致波的干涉，影响地震资料的品质，而在采集的过程中信号会受到频率混叠的影响。为此，本方法对道集内的每一道依次展开，避免了叠加导致的干涉效应，并利用匹配追踪的算法对每一道进行高精度的时频展开，提升了时频分辨率，避免了频率混叠的影响，提高了对薄储层刻画的精度。
1.基于稀疏表示的道集调谐优化增维方法，其特征在于，该基于稀疏表示的道集调谐优化增维方法包括：
步骤1，构建时频域子波字典；
步骤2，进行道集hilbert变换；
步骤3，确定匹配子波控制参数；
步骤4，进行基于匹配追踪的子波分解；
步骤5，进行基于匹配子波的高分辨时频分解。
2.根据权利要求1所述的基于稀疏表示的道集调谐优化增维方法，其特征在于，在步骤1中，在分析地震信号时选择构建ricker子波字典：
其中，ωr(t,fi)表示ricker子波字典，t代表ricker子波的中心时间，fi代表子波字典中第i个ricker子波的主频；
ricker子波具有三个控制参量，旁瓣能量的收敛速度可变,延迟时间可控以及子波库波形丰富，与地震子波匹配较为灵活。
3.根据权利要求2所述的基于稀疏表示的道集调谐优化增维方法，其特征在于，在步骤2中，对于非平稳信号来说，它的瞬时频率为：
其中，
道集内的每一道信号x(t)hilbert变换为：
式中，p为柯西主值，τ表示时移；
因此得到x(t)的解析信号为：
4.根据权利要求3所述的基于稀疏表示的道集调谐优化增维方法，其特征在于，在步骤3中，将匹配子波的三个控制参数，即时移τ，频率因子f，相位因子
5.根据权利要求4所述的基于稀疏表示的道集调谐优化增维方法，其特征在于，在步骤3中，确定了最佳匹配子波ω的中心时间和主频后，根据残差能量最小原则，利用阻尼最小平方法计算，除去这些匹配子波的振幅和相位信息；
假定存在一道带宽有限的地震信号x(t)，通过匹配追踪算法迭代了k-1次，残差信号为rk-1x；当迭代到第k次时，共产生m个匹配子波ωδm(t)，m＝1,2,…,m；对残差信号rk-1x进行hilbert变换构造解析信号：
同理构造第k次迭代时的复匹配子波：
wδm(t)＝ωδm(t) iωδm(t)(6)
所以经过k次迭代后产生的残差信号rkx的能量为：
式中x为地震信号的n未列向量，a代表匹配子波的复振幅：
am为实数域匹配子波的振幅，
gδm(t)＝[ωδ(t,f1),ωδ(t,f2),ωδ(t,f3),…,ωδ(t,fm)]为n×m阶匹配子波矩阵，每一列对应一个匹配子波；
为使得地震信号x(t)与匹配子波集ωδm(t)间的相对误差最小，采用阻尼最小平方法计算：
a＝[wtw εi]-1wtx(9)
式中，w为k次迭代后得到的n×m阶复匹配子波矩阵，ε为阻尼因子，i为m×m阶单位矩阵；由式(8)和(9)便可确定匹配子波振幅和相位的控制参数ak和
6.根据权利要求1所述的基于稀疏表示的道集调谐优化增维方法，其特征在于，在步骤4中，匹配子波通过对基本子波ω(t)进行时移、调制和相位变化形成；若基本子波ω(t)为匹配子波库中的满足条件的一个基本子波，令
式中，τ表示时移，f表示为频率因子，
7.根据权利要求6所述的基于稀疏表示的道集调谐优化增维方法，其特征在于，在步骤4中，构建好过完备匹配子波库φ后，尽管φ中的匹配子波并不是全部正交，但仍需满足||ωi||2＝1；从匹配子波库φ中选出与给定信号x(t)最为匹配的ω1，最匹配的要求则是他们的内积<x,ω1>为匹配子波库中所有原子与x(t)内积最大的一个：
<x,ω1>><x,ωi>i≠1ωi∈φ(11)
所以，一次迭代后，信号x(t)被分解为：
x(t)＝<x,ω1>ω1 r1x(12)
r1x是第一次迭代后的残差信号；r1x与ω1正交，所以：
<x,ω1>越大，意味着ω1是信号所在hilbert空间中与
rkx＝<x,ωk-1>ωk-1 rk-1x(14)
经过k次迭代分解计算后，原始信号x(t)用k个匹配子波的合成来近似，当残差r1x足够小，即信号的稀疏表示跟信号的近似程度足够好或残差信号满足设定的阈值时，迭代停止。
8.根据权利要求1所述的基于稀疏表示的道集调谐优化增维方法，其特征在于，在步骤5中，根据步骤3、4得到的高精度匹配子波，利用子波集内每一个匹配子波的独立性，选用需要的时频分解方法，对其分别进行时频域展开；再通过线性叠加，得到整个信号的时频谱图；重复对道集内的每一道进行之上的循环，最终得到多分辨去混叠的调谐数据体。
技术总结