一种极值约束的三参数扫描子波分解方法及系统与流程

本发明属于地球物理勘探领域，具体涉及一种极值约束的三参数扫描子波分解方法及系统。

背景技术：

频谱成像技术是近年来发展起来的一项基于地震频谱分解的储层预测特色解释技术，是地震属性分析中重要组成部分。该方法已逐渐成为研究复杂油气区域的一种有价值的后期处理技术，该技术主要利用频谱调谐原理描述储层厚度与分布，也可用于描述沉积相和沉积环境，检测河道和砂体，提取各种时频属性等。

地震频谱分解的核心为信号的时频分析技术，时频分析是分析非平稳信号的常规方法。常用的时频分析方法包括线性时频分析方法、双线性时频分析方法和参数化时频分析方法等。传统的stft、cwt和gst为典型的线性时频分析方法，这类方法受不确定性原理的制约，时频道集不能同时达到较高的时频分辨率。wigner-ville分布为典型的双线性时频分析方法，分解单分量平稳信号具有较高的时频分辨率；但是wvd为双线性的，在分解多分量非平稳信号时会出现严重的交叉项。

匹配追踪(又称子波分解)为典型的参数化时频分析方法，该方法由s.mallat与z.zhang于1993年提出，它可将地震信号分解为一系列子波的集合，对单个子波利用wvd计算时频谱，分解非平稳信号也能获得较高的时频分辨率。该技术自提出以来已应用到地球物理领域的各个方面，如反射系数反演、提高分辨率处理、去噪、去强反射、薄砂体预测、地震反演、地震沉积解释、地质体检测等。但是，该技术是通过不断的搜索局部最佳原子来实现地震信号的分解，计算效率较低。为了提高计算效率，yanghuawang(2007)研究了在迭代过程中使用全局粗粒度预测和局部最优相结合的匹配追踪方法，计算效率得到提高，但局部最优子波搜索过程中没有考虑参数间的关系，分解精度得不到保障；张繁昌(2010)提出了双参数动态扫描技术，提高了分解效率，但没有考虑子波尺度影响，子波尺度为一个定值，不符合实际情况。

因此，目前子波分解中局部最佳原子搜索效率和精度的问题仍亟待解决。

技术实现要素：

本发明的特征和优点在下文的描述中部分地陈述，或者可从该描述显而易见，或者可通过实践本发明而学习。

为克服现有技术的问题，本发明提供一种极值约束的三参数扫描子波分解方法，包括：

s1、对当前地震信号进行复地震道分析，计算瞬时属性；

s2、根据所述瞬时属性计算局部最优原子的初始参数；

s3、根据所述初始参数搜索局部最优原子的相位、频率、尺度并提取局部最优原子的时移；

s4、获取局部最优原子及地震信号的残差xi+1(t)，并将所述地震信号的残差xi+1(t)记为当前地震信号，重复上述步骤s1至s4，直至满足预设条件，得到一系列子波的组合。

可选地，所述瞬时属性包括瞬时振幅、瞬时频率和瞬时相位，所述初始参数包括初始时移初始频率初始相位初始尺度所述步骤s2包括：

搜索所述瞬时振幅的最大值，所述初始时移为所述瞬时振幅的最大值处对应的时间，所述初始频率为该时间对应的瞬时频率，所述初始相位为该时间对应的瞬时相位；

利用初始时移初始相位初始频率和最大匹配投影原理获得初始尺度

可选地，所述步骤s3包括：

根据最大匹配投影原理在给定的各参数变化范围内搜索局部最优原子的相位φi、频率fi和尺度σi；

在给定的时移搜索区间内搜索瞬时相位的单调区间，在单调区间内寻找与局部最优相位φi最接近的相位对应的位置作为局部最佳的时移ui。

可选地，所述步骤s4包括：

根据所述最优原子的相位、频率、尺度、时移获取局部最优原子的振幅ai，所述地震信号的残差xi+1(t)为地震信号xi(t)减去局部最优原子与振幅ai的乘积。可选地，通过下述公式分别计算局部最优原子mi(t)和局部最优原子的振幅ai：

其中，为由第i次迭代的最优原子，r⁽ⁱ⁾x为第i次迭代的地震信号。

可选地，所述步骤s4中的预设条件为：迭代次数达到预设阈值或所述地震信号的残差xi+1(t)的能量小于原始信号x(t)能量的p％。

可选地，所述预设阈值小于地震信号长度的1/3，所述p％小于20％。

本发明提供一种极值约束的三参数扫描子波分解系统，包括：

分析单元，用于对当前地震信号进行复地震道分析，计算瞬时属性；

初始参数计算单元，用于根据所述瞬时属性计算局部最优原子的初始参数；

最优参数计算单元，用于根据所述初始参数搜索局部最优原子的相位、频率、尺度并提取局部最优原子的时移；

局部最优原子获取单元，用于获取局部最优原子；

迭代计算单元，用于计算地震信号的残差xi+1(t)，并将所述地震信号的残差xi+1(t)记为当前地震信号；

迭代判断单元，用于判断是否满足预设条件，若是，则终止迭代计算，获取一系列子波的组合。

可选地，所述最优参数计算单元具体用于：根据最大匹配投影原理在给定的各参数变化范围内搜索局部最优原子的相位φi、频率fi和尺度参数σi；在给定的时移搜索区间内搜索瞬时相位的单调区间，在单调区间内寻找与局部最优相位φi最接近的相位对应的位置作为局部最佳的时移ui。

本发明提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机可执行的至少一个程序，所述至少一个程序被所述计算机执行时使所述计算机执行本发明任一实施例提供的方法中的步骤。

本发明提供的极值约束的三参数扫描子波分解方法及系统，充分考虑子波控制参数间的关系和局部范围内相位极值对最佳匹配参数的影响，提高计算效率的同时能够获得更准确的局部最优时频原子，从而实现高精度、高效率的地震子波分解。

附图说明

图1为本发明实施例的极值约束的三参数扫描子波分解方法的步骤流程图。

图2为本发明实施例的极值约束的三参数扫描子波分解方法的步骤流程图。

图3为本发明实施例的极值约束的三参数扫描子波分解系统的结构示意图。

图4为一道实际地震信号。

图5为原子波分解技术分解后的各个子波。

图6为原子波分解技术获得的地震信号的时频道集。

图7为采用本发明分解的各个子波。

图8为采用本发明计算的地震信号的时频道集。

具体实施方式

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

如图1所示，本发明提供一种极值约束的三参数扫描子波分解方法，包括：

s1、对当前地震信号进行复地震道分析，计算瞬时属性；

瞬时属性包括瞬时振幅、瞬时频率和瞬时相位。输入地震信号x(t)，对当前地震信号xi(t)进行复地震道分析，即可计算瞬时振幅、瞬时相位和瞬时频率三条曲线。xi(t)为第i次迭代的信号，第一次迭代即为原始输入地震信号x(t)。对数字信号进行复地震道分析为公知的常用技术，在三瞬属性分析中经常使用，计算过程和三瞬属性提取在公知的教科书中都有详细说明，在此不再赘述。

s2、根据所述瞬时属性计算局部最优原子的初始参数；

更具体地，搜索瞬时振幅的最大值，初始时移即为瞬时振幅最大值处对应的时间，初始频率即为该时间对应的瞬时频率，初始相位即为该时间对应的瞬时相位。

利用初始时移初始相位初始频率和最大匹配投影原理获得初始尺度

s3、根据所述初始参数搜索局部最优原子的相位、频率、尺度并提取局部最优原子的时移；

根据最大匹配投影原理在给定的各参数变化范围内搜索局部最优的相位φi、频率fi和尺度参数σi。

在给定的时移搜索区间内搜索瞬时相位的单调区间，在单调区间内寻找与局部最优相位φi最接近的相位对应的位置作为局部最佳的时移ui。

在具体实施时，初始尺度是在一组固定的下，通过计算最优化公式(1)得到的；

其中，d＝{mr(t)}r∈γ为时频原子字典，是函数r⁽ⁱ⁾x和的内积，且r⁽ⁱ⁾x表示第i次迭代的剩余信号，第1次迭代即为x(t)，mr(t)为由参数r决定的时频原子，时频原子由公式(2)获得，ri＝{ui,σi,fi,φi}，i表示迭代次数。

在初始时移固定的情况下，根据最大匹配投影原理在给定的各参数变化范围内搜索局部最优的相位φi、频率fi和尺度参数σi。

本实施例中，对相位、频率和尺度三个参数进行局部优化是在局域内寻找ri＝{ui0,σi,fi,φi}的最优值，其中初始时移ui0固定不变，搜索范围为[ri-δr,ri+δr]，其中δr＝(δu,δσ,δf,δφ)，即δu＝0为时间偏移量、δσ＝0.1为尺度偏移量、δf＝5hz为频率偏移量、δφ＝50为相位偏移量。在具体实施时，各参数给定的范围内任意取不同的值，由公式(2)获得一组时频原子，使公式(1)达到最大时对应的参数即为局部最优的φi、fi和σi。

根据局部最佳的相位和瞬时相位曲线的关系获得局部最佳的时移。更具体地，在给定的时移搜索区间内搜索瞬时相位的有效单调区间在该区间内瞬时相位为单调递增或单调递减，在该区间内搜索与局部最佳相位φi差值最小的相位位置，该位置即为搜索的局部最佳时移ui。

s4、获取局部最优原子及地震信号的残差，将所述地震信号的残差xi+1(t)记为当前地震信号，重复上述步骤s1至s4进行迭代分解，直至满足预设条件，得到一系列子波的组合。

根据所述最优原子的相位、频率、尺度、时移获取局部最优原子的振幅ai；再根据所述最优原子的相位、频率、尺度、时移、振幅获取地震信号的残差xi+1(t)。地震信号的残差xi+1(t)为地震信号xi(t)减去局部最优原子与振幅ai的乘积。

更具体地，根据局部最优原子的相位φi、频率fi、尺度σi和时移ui，结合公式(3)获得局部最优原子的振幅ai。公式(3)如下所示：

其中，为由第i次迭代的最优参数ri＝{ui,σi,fi,φi}和公式(2)决定的局部最优原子，r⁽ⁱ⁾x为第i次迭代的地震信号，第1次迭代的r⁽ⁱ⁾x为原始地震信号x(t)。

将局部最优原子的相位φi、频率fi、尺度σi和时移ui参数带入公式(2)获得局部最优的原子mi(t)。将地震信号xi(t)减去局部最优原子与振幅ai的乘积，获得地震信号的残差xi+1(t)。用于计算残差xi+1(t)的公式如下所示：

xi+1(t)＝xi(t)-aimi(t)(4)

预设条件包括迭代次数达到预设阈值或所述地震信号的残差xi+1(t)的能量小于原始信号x(t)能量的p％。其中，预设阈值小于地震信号长度的1/3，p不大于20。本实施例中，退出条件为事先约定的迭代次数167次或所述残差xi+1(t)的能量小于原始信号x(t)能量的10％。

最终原始地震信号被分解为一系列子波的组合，如公式(5)所示：

请参照图2，本发明提供一种极值约束的三参数扫描子波分解方法，包括步骤：

101、开始。

102、输入当前地震信号；

103、复地震道分析；

对当前地震信号xi(t)进行复地震道分析，计算瞬时振幅、瞬时相位和瞬时频率。根据瞬时振幅获得初始时移、初始频率和初始相位；更具体地，搜索瞬时振幅的最大值，初始时移即为瞬时振幅最大值处对应的时间，初始频率即为该时间对应的瞬时频率，初始相位即为该时间对应的瞬时相位。

104、匹配初始尺度；

在初始时移、初始频率和初始相位的基础上获得初始尺度；本实施例中，利用初始时移初始相位初始频率和上述公式(1)和(2)获得初始尺度

105、局部动态扫描；

在初始时移固定的情况下，根据最大匹配投影原理在给定的各参数变化范围内搜索局部最优的相位φi、频率fi和尺度参数σi

106、提取局部最佳原子的时移；

利用“复地震道分析”中计算的瞬时相位属性和“局部动态扫描”得到的局部最佳的相位提取局部最佳原子的时移。更具体地，在给定的时移搜索区间内搜索瞬时相位的有效单调区间在该区间内瞬时相位为单调递增或单调递减，在该区间内搜索与局部最佳相位φi差值最小的相位位置，该位置即为搜索的局部最佳时移ui。

107、计算局部最佳原子；

将局部最优原子的相位φi、频率fi、尺度σi和时移ui参数带入公式(2)获得局部最优的原子mi(t)。

108、判断是否满足预设条件，若是，则进入步骤11，即结束进程；若否，则进入步骤109；

109、去除局部最佳原子获得残差；

根据所述最优原子的相位、频率、尺度、时移结合公式(3)获取局部最优原子的振幅ai；再根据所述最优原子的相位、频率、尺度、时移、振幅获取地震信号的残差xi+1(t)。地震信号的残差xi+1(t)为地震信号xi(t)减去局部最优原子与振幅ai的乘积。

预设条件包括迭代次数达到预设阈值或所述地震信号的残差xi+1(t)的能量小于原始信号x(t)能量的p％。其中，预设阈值小于地震信号长度的1/3，p不大于20。本实施例中，退出条件为事先约定的迭代次数167次或所述残差xi+1(t)的能量小于原始信号x(t)能量的10％。

110、将残差记为当前地震信号，并返回步骤102。

111、结束进程。

最终原始地震信号被分解为如公式(5)所示的一系列子波的组合。

本发明提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机可执行的至少一个程序，所述至少一个程序被所述计算机执行时使所述计算机执行本发明任一实施例提供的方法中的步骤。

如图3所示，本发明提供一种极值约束的三参数扫描子波分解系统，包括：分析单元10、初始参数计算单元20、最优参数计算单元30、局部最优原子获取单元40、迭代计算单元50、迭代判断单元60。其中：

分析单元10用于对当前地震信号进行复地震道分析，计算瞬时属性；更具体地，分析单元用于实现步骤s1，在此不再赘述。

初始参数计算单元20与分析单元10相连，用于根据所述瞬时属性计算局部最优原子的初始参数；更具体地，初始参数计算单元用于实现步骤s2，在此不再赘述。

最优参数计算单元30与所述分析单元10、初始参数计算单元20相连，最优参数计算单元30用于根据所述初始参数搜索局部最优原子的相位、频率、尺度并提取局部最优原子的时移；更具体地，最优参数计算单元用于实现步骤s30，在此不再赘述。

局部最优原子获取单元40用于获取局部最优原子；将局部最优原子的相位φi、频率fi、尺度σi和时移ui参数带入公式(2)获得局部最优的原子mi(t)。

迭代计算单元50与局部最优原子获取单元40相连，用于计算地震信号的残差xi+1(t)，并将所述地震信号的残差xi+1(t)记为当前地震信号。地震信号的残差xi+1(t)为地震信号xi(t)减去局部最优原子与振幅ai的乘积；根据所述最优原子的相位、频率、尺度、时移结合公式(3)获取局部最优原子的振幅ai。

迭代判断单元60与迭代计算单元50相连，用于判断是否满足预设条件，若是，则终止迭代计算，获取一系列子波的组合。预设条件包括迭代次数达到预设阈值或所述地震信号的残差xi+1(t)的能量小于原始信号x(t)能量的p％。其中，预设阈值小于地震信号长度的1/3，p不大于20。本实施例中，退出条件为事先约定的迭代次数167次或所述残差xi+1(t)的能量小于原始信号x(t)能量的10％。

下面通过一个具体的实例说明本发明实施例的有益效果。

图4示出了一道实际地震信号，从图中可以看出地震波形的最大值和最小值在800ms-1000ms之间，计算的时频谱的能量最大值也应该在800ms-1000ms之间。

图5示出了原有的四参数动态扫描子波分解方法分解得到的各个子波，从左到右依次为每一步迭代输出的子波，迭代次数为80次，迭代时间为381秒。从图中的实线圆圈可以看出在同一位置出现了多个子波，极性相反、且下一个子波的能量大于了上一个子波的能量，出现这种现象是由于在搜索时移时偏离了瞬时振幅最大值的位置，其根本原因在于搜索局部最优原子的时移参数时，在时移扰动范围内瞬时相位出现了跳跃，即从π突然变化到-π，搜索的局部最优原子的相位大小没有问题，但是时移出现了差异，偏离了瞬时振幅的最大位置，导致本次迭代的局部最优的原子不正确，那么计算的残差信号也不正确，直接影响到后续的迭代过程。同时，从图中还能看出，这样情况一般都是成对出现，原因在于第i次迭代时该位置瞬时振幅能量最大，但是拾取局部最优原子出现了偏差，残差信号在该位置附近能量得到加强，导致下一次迭代搜索的子波能量高于上一次迭代，不符合迭代逻辑。

图6示出了原有的四参数动态扫描子波分解方法分解的子波计算的时频谱。由原始地震信号图5所示，地震信号的最大能量应该在800ms-1000ms之间，按照时频分析原理可知，时频分析的最大能量团也应该在800ms-1000ms之间，但是图6中的最大能量(最黑的部位)在1400ms以下，这与图5所示的第一个黑色圆圈中的最大能量子波一致；1000ms以下的强能量团与图5中最右边的圆圈里的子波有关，同一个位置出现了三个子波(极性相反)，直接导致时频谱的能量团能量出现问题。由于分解过程没有考虑相位的影响，导致分解子波不正确，所以计算的时频谱也不正确，直接影响到时频谱后续的应用。

图7为本发明提取的极值约束的三参数扫描子波分解方法分解的子波，迭代次数为60次，迭代时间为19秒。从图中可以看出，子波的能量不断减小，这与实际迭代分解思路一致，且没有两个子波在同一位置的情况，得益于本发明提出的相位极值的约束。第一个子波的能量最大，出现在800ms-1000ms之间，图8给出了本发明计算的时频谱，从时频谱的能量关系也可以看出，时频谱最大能量也出现在800ms-1000ms之间与，与地震信号最大能量位置保持一致。迭代次数相差不大，但是计算时间大大减小。

总之，与现有子波分解技术相比，本发明充分考虑了各时频原子决定参数之间的关系，将现有的四参数动态搜索减少为三参数动态搜索，提高了计算效率；同时，考虑了相位突变对局部最优时移的影响，在瞬时相位的单调区间内搜索最佳相位，由局部最佳相位确定局部最佳时移，提高了子波分解的精度。

本发明提供一种极值约束的三参数扫描子波分解方法及系统，能够实现更高精度的子波分解、提高计算效率。本发明充分考虑了各时频原子决定参数之间的关系，将现有的四参数动态搜索减少为三参数动态搜索，提高了计算效率；同时，考虑了相位突变对局部最优时移的影响，在瞬时相位的单调区间内搜索最佳相位，由局部最佳相位确定局部最佳时移，提高了子波分解的精度。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。