CN111856562A - 一种广义高阶同步挤压地震信号时频分解与重构方法 - Google Patents

Abstract

本发明公开了一种广义高阶同步挤压地震信号时频分解与重构方法，包括输入一维信地震号；选择阶数和时窗函数，计算在窗函数下的短时傅里叶变换值，设置阈值，计算K阶矩阵L_K(t,ω)并根据阈值得到m值；构建m阶矩阵L_m(t,ω)和M_m(t,ω)；计算在二维时频谱上的K阶瞬时频率估计值

再用K阶瞬时频率估计值对短时傅里叶变换值进行挤压。本发明根据信号的特征灵活地选择合适的时窗函数，得到聚集的短时傅里叶变换时频谱，且能够直接求解任意阶估计值，有效提高了计算效率，提升了瞬时频率估计值的精确度，使得高阶同步挤压变换的时频谱能量更为集中，极大地提高了信号的时间和频率分辨率，同时在挤压过程中也保持了该变换的可逆性，能够高精度的重构原信号。

Description

一种广义高阶同步挤压地震信号时频分解与重构方法

技术领域

本发明涉及一种地震信号处理方法，尤其涉及一种广义高阶同步挤压地震信号时频分解与重构方法。

背景技术

作为非平稳信号处理领域的一个重要分支，时频分析一直是现代信号处理的研究热点之一。常用的时频分析方法有短时傅里叶变换(STFT)、小波变换(CWT)、S变换(ST)及广义S变换(GST)等。虽然这些时频分析方法的分辨率都有所提高，但由于他们都是借助某一窗函数(不断改进的窗函数)并在该窗函数的区间内对信号的频率成分进行估计，使得其时频能量总是分布在一个以信号真实的瞬时频率为中心的范围内。因此，这些时频变换得到的时频谱的能量聚焦性不够集中。同时受Heisenberg-Gabor不定问题影响，其时间分辨率和频率分辨率不能同时达到最优。

同步挤压变换(SST)是Daubechies等在小波变换的基础上提出的一种新的时频变换方法。它通过严格的数学推导，把小波变换结果在一定频率范围内的时频能量“挤压”到信号真实的瞬时频率附近，极大地提高了时频分辨率。而在同步挤压变换处理过程中，对信号真实的瞬时频率估计的准确性则显得尤为重要，在很大程度上会影响最终的“挤压”效果。

发明内容

本发明的目的就在于提供一种解决上述问题，能根据信号的特征灵活地选择恰当的阶数和时窗函数，使得到的时频谱的能量更为集中，极大地提高地震信号时频分辨率的一种广义高阶同步挤压地震信号时频分解与重构方法。

为了实现上述目的，本发明采用的技术方案是这样的：一种广义高阶同步挤压地震信号时频分解与重构方法，包括以下步骤：

(1)获取一地震信号x(t)；

(2)预设一阶数K和时窗函数g(t)，计算在窗函数t^kg(t)下的短时傅里叶变换值

其中K≤6，k＝0,1,…,max{2K-2,1}；

(3)预设一阈值ε，根据下式构建K阶矩阵L_K(t_,ω)，并查找矩阵L_K(t_,ω)中所有顺序主子式det(L_m(t_,ω))绝对值大于ε的最大值m；

m＝max{j:|det(L_j(t,ω))|＞ε,j＝1,2,…,K}

式中det表示求解矩阵行列式，|·|表示取绝对值；

(4)构建m阶矩阵L_m(t,ω)和M_m(t,ω)，其中，

(5)根据下式计算二维时频谱上的K阶瞬时频率估计值

式中，Re表示取复数的实部；

(6)根据下式对短时傅里叶变换值

进行挤压，得到高阶同步挤压变换值

作为优选：还包括步骤(7)，对步骤(6)中的

取模操作得到最终的高阶同步挤压变换时频谱

作为优选：还包括步骤(7)，采用下式对步骤(6)中的

进行逆变换，重构地震信号x(t)；

作为优选：所述步骤(2)中，采用下式计算

式中，t为时间，ω为频率，σ为时窗函数的调节因子。

作为优选：所述阈值ε，用于削弱步骤(5)对瞬时频率估计所产生的误差，使步骤(6)获得稳定的高度聚焦的时频表征结果。

本发明的思路为：

首先，获取一待测的地震信号，该信号为随时间变化的一维信号x(t)。

第二，本发明通过选择阶数和时窗函数，构建出一个合适的窗函数t^kg(t)，来对地震信号进行短时傅里叶变换，可以得到分辨率更高的短时傅里叶变换值。

第三，本发明设置了阈值ε，这是由于离散化的处理过程和奇异矩阵会影响对瞬时频率估计的准确性，为此我们设置了一个阈值来控制估计结果的误差；该阈值为一个大于0的实数，主要目的是为了控制瞬时频率估计结果的误差，以此来保证挤压结果的高度聚焦性。通过设置阈值，我们能找到一个满足阈值条件的m阶矩阵，在该矩阵下得到的瞬时频率估计值的误差是相对比较小的，且挤压结果的时频能量分布也比较聚焦。

第四，在m的基础上，构建两个m阶矩阵L_m(t,ω)和M_m(t,ω)；来计算m阶矩阵L_m(t,ω)和M_m(t,ω₎；其目的是为了更加准确地对地震信号瞬时频率进行估计，从而得到K阶瞬时频率，这样得到的K阶瞬时频率的误差相对较小，且能保证挤压后的时频结果的能量分布高度聚焦。

第五，我们用的K阶瞬时频率对短时傅里叶变换值进行挤压，可以获得一个比较稳定的能量高度聚焦的时频表征结果。

与现有技术相比，本发明的优点在于：

(1)本发明以短时傅里叶变换和同步挤压变换为基础，提出了一种新的时频分析方法-广义高阶同步挤压变换。该发明采用短时傅里叶变换，能够根据不同的地震信号特性选择合适的时窗函数，得到时频分辨率更高的短时傅里叶变换变换时频谱。

(2)本发明利用泰勒公式对地震信号瞬时频率进行高阶估计，以数学理论为基础，直接从短时傅里叶变换式出发推导出地震信号瞬时频率高阶估计的统一求解式。与现有的高阶同步挤压变换相比，该求解式摒弃了传统求解式的迭代算法，不再依赖于低阶求解式的求解结果，能够直接求解任意阶估计值，有效提高了计算效率，高阶估计的结果使得对地震信号瞬时频率估计的精确度得以提高。

(3)利用高精度的瞬时频率估计值对短时傅里叶变换时频谱进行挤压操作，可保证挤压后的时频面能量聚集在更窄的时频区间上，使得地震信号的时频能量更加集中，时频分辨率得到更进一步的提升。

(4)基于高时频分辨率的高阶同步挤压变换，我们能实现对原地震信号进行高精度的重构。

综上本发明能使得到的时频谱的能量更为集中，极大地提高地震信号时频分辨率，使地震信号的时间分辨率和频率分辨率能同时达到最优，从而有助于更好地揭示地震信号在时频域的响应特征，在后续对地震信号的分析过程中，有助于更准确的分析，且本发明方法具有可重构性。

附图说明

图1为本发明流程图；

图2为信号的理想时频谱图；

图3为信号的短时傅氏变换谱图；

图4为经本发明方法处理后的一阶同步挤压变换时频谱图；

图5为经本发明方法处理后的二阶同步挤压变换时频谱图；

图6为经本发明方法处理后的三阶同步挤压变换时频谱图；

图7为经本发明方法处理后的四阶同步挤压变换时频谱图。

具体实施方式

下面将结合附图对本发明作进一步说明。

实施例1：参见图1，一种广义高阶同步挤压地震信号时频分解与重构方法，包括以下步骤：

(1)获取一地震信号x(t)；

(2)预设一阶数K和时窗函数g(t)，计算在窗函数t^kg(t)下的短时傅里叶变换值

其中K≤6，k＝0,1,…,max{2K-2,1}；我们采用下式计算

式中，t为时间，ω为频率，σ为时窗函数的调节因子；

(3)预设一阈值ε，根据下式构建K阶矩阵L_K(t,ω)，并查找矩阵L_K(t,ω)中所有顺序主子式det(L_m(t,ω))绝对值大于ε的最大值m；

m＝max{j:|det(L_j(t,ω))|＞ε,j＝1,2,…,K}

式中det表示求解矩阵行列式，|·|表示取绝对值；所述阈值ε，用于削弱步骤(5)对瞬时频率估计所产生的误差，使步骤(6)获得稳定的高度聚焦的时频表征结果；

(4)构建m阶矩阵L_m(t_,ω)和M_m(t_,ω)，其中，

(5)根据下式计算二维时频谱上的K阶瞬时频率估计值

式中，Re表示取复数的实部；

(6)根据下式对短时傅里叶变换值

进行挤压，得到高阶同步挤压变换值

实施例2：参见图1，在实施例1的基础上，还包括步骤(7)，对步骤(6)中的

取模操作得到最终的高阶同步挤压变换时频谱

本实施例是用于输出高阶同步挤压变换时频谱。

实施例3，参见图1，在实施例1的基础上，还包括步骤(7)，采用下式对步骤(6)中的

进行逆变换，重构地震信号x(t)；

本实施例用于重构地震信号。

实施例4：参见图1到图7，在实施例1的基础上，我们具体以一种信号为例，该信号的理想时频谱图如图2所示。图3-图7分别是信号的短时傅氏变换谱、一阶同步挤压变换时频谱、二阶同步挤压变换时频谱、三阶同步挤压变换时频谱和四阶同步挤压变换时频谱。图中，横坐标time(s)表示时间，纵坐标Fre(Hz)表示频率，右侧灰度条表示能量。

通过对比发现，经本发明方法处理后的时频谱图明显优于短时傅氏变换谱，且随着阶数的增加，时频分辨率更高，能量更加集中，更加逼近理想的时频谱图。经本发明方法处理后的时频谱图明显优于图3所示的短时傅氏变换谱，且随着阶数的增加，时频分辨率更高，能量更加集中，更加逼近于图2所示的理想时频谱图。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种广义高阶同步挤压地震信号时频分解与重构方法，其特征在于：包括以下步骤：

(1)获取一地震信号x(t)；

(2)预设一阶数K和时窗函数g(t)，计算在窗函数t^kg(t)下的短时傅里叶变换值

其中K≤6，k＝0,1,…,max{2K-2,1}；

(3)预设一阈值ε，根据下式构建K阶矩阵L_K(t,ω)，并查找矩阵L_K(t,ω)中所有顺序主子式det(L_m(t,ω))绝对值大于ε的最大值m；

m＝max{j:|det(L_j(t,ω))|＞ε,j＝1,2,…,K}

式中det表示求解矩阵行列式，|·|表示取绝对值；

(4)构建m阶矩阵L_m(t,ω)和M_m(t,ω)，其中，

(5)根据下式计算二维时频谱上的K阶瞬时频率估计值

式中，Re表示取复数的实部；

(6)根据下式对短时傅里叶变换值

进行挤压，得到高阶同步挤压变换值

2.根据权利要求1所述的一种广义高阶同步挤压地震信号时频分解与重构方法，其特征在于：还包括步骤(7)，对步骤(6)中的

取模操作得到最终的高阶同步挤压变换时频谱

3.根据权利要求1所述的一种广义高阶同步挤压地震信号时频分解与重构方法，其特征在于：还包括步骤(7)，采用下式对步骤(6)中的

进行逆变换，重构地震信号x(t)；

4.根据权利要求1所述的一种广义高阶同步挤压地震信号时频分解与重构方法，其特征在于：所述步骤(2)中，采用下式计算

式中，t为时间，ω为频率，σ为时窗函数的调节因子。

5.根据权利要求1所述的一种广义高阶同步挤压地震信号时频分解与重构方法，其特征在于：所述阈值ε，用于削弱步骤(5)对瞬时频率估计所产生的误差，使步骤(6)获得稳定的高度聚焦的时频表征结果。

Images