CN113281809B - 一种地震信号的谱分析方法 - Google Patents

Abstract

本发明公开了一种地震信号的谱分析方法，该方法首先设计高斯窗函数，对于输入信号进行滑窗傅里叶变换，得到广义时频变换，对得到的广义时频变换，进行时频浓度测量，得到时频聚集度CM，计算时频能量聚集度的四个参数，遍历限制条件内所有参数获取CM集，再进行参数优化，筛选出CM值最大时的窗口参数。最后，应用优化后的参数获取最优广义时频变换。经实验证明，该方法拥有更高的分辨率和能量聚集性以及更好的抗噪性能，便于分析信号的特征信息。

Description

一种地震信号的谱分析方法

技术领域

本发明涉及信号处理领域，特别涉及一种地震信号的谱分析方法，可应用于地震资料的处理与解释、地震薄层结构的探究以及油气储层的预测等。

背景技术

在地下勘探的探索与资料收集过程中，常用人工激发地震波的方式，根据反馈的信号特征分析地下介质等信息。但是由于地下介质的分布不均匀且种类繁多等特征，信号的传播路径，振动幅度都会受到不同程度的改变，加上采集地震波时的各种噪声影响，所以收集到的信号都是时变非平稳信号。

对于地震信号的分析是地震资料解释与处理的必要步骤。传统的傅里叶变换是一种全局性质的信号分析方法，不能对信号局部细节进行分析，由于地震信号的复杂性，该方法已经不能满足实际的需要。因此使用时频联合分布的方法对信号进行分析是必须的。

短时傅里叶变换将信号从一维时间域分析变换到二维时间频率联合分布进行分析，在时域用是窗函数对输入信号进行截取，对被截取的信号进行快速傅里叶变换，然后移动窗口，重复这个过程直到获取整个时间轴上信号的频谱。短时傅里叶变换能很好地对地震信号进行局部特征的精细刻画。但是，短时傅里叶变换的局限在与其窗函数选定后对于非平稳信号来说是固定的，会导致分辨率固定单一。

Wigner-Ville分布是一种高精度时频分析方法，因为没有使用窗函数，所以也不受海森堡不确定性原理的限制，拥有较高的时频聚集度。但是其缺点是存在严重的交叉干扰项，在实际应用中，由于干扰项的存在，得到的频谱不能准确地分析信号的特征信息。

发明内容

针对传统方法存在的缺陷或不足，本发明的目的在于，提供一种新的地震信号的谱分析方法，该方法通过引入四个参数根据需要灵活的控制时窗宽度，从而获得更高的时频分辨率。

为了实现上述任务，本发明采取如下的技术解决方案：

一种地震信号的谱分析方法，其特征在于，包括以下步骤：

步骤一，对于输入信号长度为x(n),n＝0,1,2,...,N-1,N的表达式，N取偶数，设计高斯窗函数如下式(1)所示：

式中，f为频率，

其中，T为采样后数据时间长度A,B,C,D为可调节的控制窗口形态的参数；为保证控制窗口有一个合适的形态，规定0≤A,B,C,D≤2；

步骤二，基于高斯窗函数，对输入信号进行滑窗Fourier变换，得到广义时频变换如下式(2)所示：

式(2)即为非平稳信号x(n)的新型广义时频变换的离散形式，它给出了信号在n＝m附近的一段时间内的时频信息，V(n,f)的相位谱和幅度谱均由相位函数e^-j2πfm定位；

步骤三，对得到的广义时频变换，进行时频浓度测量，得到时频聚集度CM如下：

其中，V(n,f)的模归一化后为

它的模定义为：

步骤四，计算时频能量聚集度最优的四个参数如下式所示：

{A_opt,B_opt,C_opt,D_opt}＝argmax{CM^ABCD} (5)

步骤五，应用最优的参数A_opt,B_opt,C_opt,D_opt，得到最优窗函数ω_opt(n)如下式所示：

利用最优窗函数ω_opt(n)进行滑窗傅里叶变换，得到非平稳信号x(n)在n＝m处最优的广义时频变换，如下下式所示：

步骤六，应用最优的广义时频变换对地震型号进行分析，获得最终结果V_opt(n,f)并提取瞬时频率等信号特征信息，绘制时频联合分布效果图，描述信号在不同时间和频率的能量强度；基于对图像的分析，以实现储层流体预测，油气检测，薄层辨析等与地震信号相关的数据处理。

本发明的地震信号的谱分析方法，简谐波在时域内的伸缩变换与该窗函数在时域内的平移伸缩变换合作，形成一种新的适用于地震信号的谱分解法。

与现有技术相比具有如下优点：

1、从本质上来说是一种新的时频联合分布的方法，具有更高的时频分辨率，便于人们理解和使用。

2、传统的时频分析方法存在窗函数固定从而导致分辨率单一的限制，而本发明通过改造高斯窗函数，在分子分母共引入四个参数的情况下，可以更加灵活的控制窗口形状且更好的消除因窗口形状引起的误差。

3、控制窗口形状的标准差进行范围约束，不会存在因窗口过宽或过窄而引起的时频分辨率下降等问题。

4、应用的是基于能量浓度的参数优化方法，可以筛选出时频浓度最大值时的参数，然后应用最优参数进行时频分析，从而获得更高的时频分辨率；

5、相比于传统的时频分析方法拥有更好的抗噪性能，在低信噪比的情况下，信号的时频特征仍然清晰可辨。

附图说明

图1为本发明的地震信号的谱分析方法流程图；

图2为参数a，b对时窗形状的影响对比；

图3为参数c,d对时窗形状的影响对比；

图4为合成信号的振幅谱；

图5为对合成信号进行短时傅里叶变换得到的时频图；

图6为对合成信号进行连续小波变换得到的时频图；

图7为对合成信号进行平滑伪魏格纳-维尔分布的时频图；

图8为对合成信号用本发明所提方法得到的时频图；

图9合成信号在1db下应用本发明的时频分析图；

图10为合成楔状模型地震剖面；

图11为楔状模型应用现有Stockwell变换的时频分析图；

图12为楔状模型应用本发明的时频分析图。

以下结合附图和实施例对本发明作进一步的详细说明。

具体实施方式

本实施例给出一种地震信号的谱分析方法，具体包括以下步骤；

步骤一，对于输入信号长度为x(n),n＝0,1,2,...,N-1,N的表达式，N取偶数，设计高斯窗函数如下式(1)所示：

式中，f为频率，

其中，T为采样后数据时间长度A,B,C,D为可调节的控制窗口形态的参数；为保证控制窗口有一个合适的形态，规定0≤A,B,C,D≤2；

步骤二，基于高斯窗函数，对输入信号进行滑窗Fourier变换，得到广义时频变换如下式(2)所示：

式(2)即为非平稳信号x(n)的新型广义时频变换的离散形式，它给出了信号在n＝m附近的一段时间内的时频信息，V(n,f)的相位谱和幅度谱均由相位函数e^-j2πfm定位；

步骤三，对得到的广义时频变换，进行时频浓度测量，得到时频聚集度CM如下：

其中，V(n,f)的模归一化后为

它的模定义为：

步骤四，计算时频能量聚集度的四个参数如下式所示：

{A_opt,B_opt,C_opt,D_opt}＝argmax{CM^ABCD} (5)

步骤五，应用参数A_opt,B_opt,C_opt,D_opt，得到窗函数ω_opt(n)如下式所示：

利用窗函数ω_opt(n)进行滑窗傅里叶变换，得到非平稳信号x(n)在n＝m处最优的广义时频变换，如下下式所示：

步骤六，应用广义时频变换对地震型号进行分析，获得最终结果V_opt(n,f)并提取瞬时频率信号特征信息，绘制时频联合分布效果图，描述信号在不同时间和频率的能量强度；基于对图像的分析，以实现储层流体预测，油气检测，薄层辨析等与地震信号相关的数据处理。

本实施例中，为保证时频分析的能量聚集性，高斯窗不可太宽也不可太窄，因此对影响窗函数的标准差作出以下限制：

L<σ(f)<U

其中，L取nTs，U取lTs，n和l分别取10，1000；f∈[f_min,f_max],f_min,f_max取决于待分析信号，那么，对于所述计算时频能量聚集度的四个参数有如下限制：

0≤A,B,C,D≤2。

所述的高斯窗函数按如下方式进行处理：

归一化参数

使高斯窗成为一个单位面积，从而使得新的广义时频变换得幅度具有与傅里叶变换相同的意义。由于高斯窗是关于频率的函数，因此窗口形态会随着信号频率的变化而变化，那么高斯窗在时间域不仅做平移变换，同时也在做伸缩变换，同样也意味着广义时频变换是一种多分辨率的信号分析方法；当[A B C D]＝[0 1 1 0]时，所述的高斯窗即等于标准的高斯窗。

关于参数对于窗口的影响，参数B和D用于粗调窗口形态，参数A和C用于微调窗口形态。

在时间域，窗口长度随着参数A，B，D的增加而变窄，随着参数C的增加而变长；窗口幅度随着参数A和B的增加而增加，随着参数C和D的增加而减小；其中B和D的影响较为剧烈，用于粗调窗口形态，而参数A和C的影响较为微小，用于微调窗口形态。

本实施例中，聚集度CM计算方法类似于统计学中“峰值”的定义。不同的窗口参数条件下，对应有不同的时频分布，不同的时频分辨率，相应地拥有不同的时频聚集度CM。

当聚集度CM值越大时，意味着时频分布越为集中。因此在本实施例中，将CM作为参数优化的目标函数，选择聚集度最高条件下的窗口参数，以此获得分辨率更高的时频分布，便于分析信号在不同时间对应的频率特征。

以下是发明人给出的具体实施例。

实施例：

针对传统的时频分析方法在地震信号中存在窗口形态固定，分辨率单一，能量聚集性差等问题，本实施例给出一种地震信号的谱分析方法。首先应用新的窗函数对各地震道信号进行滑窗傅里叶变换，然后计算不同参数条件下的聚集度CM值，选择最优参数进行广义时频变换。实验证明，该方法比传统的时频分析方法拥有更高的分辨率和能量密度。更高的时频分辨率意味着对地震信号时频特征更清晰地刻画，便于后续薄层识别和勘探等领域的研究。

具体的实现步骤包括：

(1)输入待分析的原始信号x(t)，合成信号x(t)是由三个子波信号叠加而成，合成信号的振幅谱参见图4所示。

(2)应用高斯窗，对信号x(t)进行滑窗傅里叶变换，新型滑动窗的窗口形态根据参数设置不同而变化；

(3)为获取时频分辨的最佳效果，将时频聚集度CM作为目标函数，并在限定的约束条件内进行参数优化。

首先，计算边界条件[0-2]内所有参数进行滑窗傅里叶变换后的时频聚集度CM，并选择当CM最大时的参数{A_opt,B_opt,C_opt,D_opt}。

应用优化后的参数进行广义时频分析。获得的时频效果图参见图8，优化后参数选择为[A B C D]＝[0.8128 0.5075 1.0211 0.1221]。

对于合成信号，其中x₁(t)是雷克子波主频为10Hz分布在0.1s，0.5s，0.8s的子波信号，x₂(t)是主频为30Hz分布在0.2s，0.3s，0.6s和0.7s的雷克子波信号，x₃(t)是雷克子波主频为80Hz分布在0.1s，0.4s，0.65s的子波信号。

其中信号数学模型如下：

x(t)＝x₁(t)+x₂(t)+x₃(t)

对于高斯窗，作用类似于对被截取的信号进行局部加权，从而突出信号的局部特征。参数对时窗形状的影响参见图2和图3所示，在主频确定的情况下，由图2可知，时窗随着参数A，B的增加而变得窄而尖，从而使时间分辨率提高，频率分辨率降低；反之，由图3可知，时窗随着参数C，D的增加而变得更宽更扁，从而使频率分辨率提高，时间分辨率降低。虽然四参数对时窗控制大体趋势有重复地方，但参数之间造成的影响大小，程度效果各不相同。四参数共同参与的控制使窗口形状更加灵活。

根据图5、图6、图7、图8进行对比分析发现，短时傅里叶变换在低频时能量聚集度较高，时间分辨率较高，但整体时频能量模态混叠，很难准确分辨信号特征。连续小波变换相比较于短时傅里叶变换，整体时频聚集度有所提升，但低中频信号时频特征仍难以分辨。对于平滑伪魏格纳-维尔分布，虽然时频聚集度大幅提升，但由于其方法本身存在严重交叉干扰项的缺陷，难以分辨信号本身特征。本实施例采用的地震信号谱分析方法，对于合成信号效果显著，时频能量集中，时间分辨率高，信号分布清晰可辨。

为了验证本实施例的地震信号的谱分析方法的优异抗噪性能，给合成信号加上高斯白噪声，在1db下对信号进行时频变换，时频分析结果参照图9，可以看到，尽管在低信噪比的干扰下，原信号的时频能量仍然较为集中，可以清晰分辨合成信号中主要成分分布，直到1db下才稍有模糊，但时间分辨率仍较为准确。

实验例：

参见图10，图11和图12，设计一个楔状地震信号模型，采集时长共500ms，共60道地震数据，采用主频为30Hz的雷克子波作为震源子波，该模型顶部和底部的反射系数分别为0.4，-0.4。图10为楔状模型合成地震记录。图11为对合成的地震记录应用现有Stockwell变换得到的时频结果图，可以看到在楔状模型的第10道，第30道和第50道地震数据及附近，合成信号之间的干涉较为严重，对薄层的识别效果模糊。图12为对合成地震记录应用对每一道数据进行时频分析，并提取主频，做振幅切片，通过该振幅谱对薄层界面进行分辨，可以看到分辨率显著提升，薄层识别更为精确。

Claims (3)

1.一种地震信号的谱分析方法，其特征在于，包括以下步骤：

步骤一，对于输入信号长度为x(n),n＝0,1,2,...,N-1,N的表达式，N取偶数，设计高斯窗函数如下式(1)所示：

式中，f为频率，

其中，T为采样后数据时间长度A,B,C,D为可调节的控制窗口形态的参数；为保证控制窗口有一个合适的形态，规定0≤A,B,C,D≤2；

步骤二，基于高斯窗函数，对输入信号进行滑窗傅里叶变换，得到广义时频变换如下式(2)所示：

式(2)即为非平稳信号x(n)的广义时频变换的离散形式，它给出了信号在n＝m附近的一段时间内的时频信息，V(n,f)的相位谱和幅度谱均由相位函数e^-j2πfm定位；

步骤三，对得到的广义时频变换，进行时频浓度测量，得到时频聚集度CM如下：

其中，V(n,f)的模归一化后为

它的模定义为：

步骤四，计算时频能量聚集度的四个参数如下式所示：

{A_opt,B_opt,C_opt,D_opt}＝argmax{CM^ABCD} (5)

步骤五，应用参数A_opt,B_opt,C_opt,D_opt，得到窗函数ω_opt(n)如下式所示：

利用窗函数ω_opt(n)进行滑窗傅里叶变换，得到非平稳信号x(n)在n＝m处最优的广义时频变换，如下下式所示：

步骤六，应用广义时频变换对地震型号进行分析，获得最终结果V_opt(n,f)并提取瞬时频率信号特征信息，绘制时频联合分布效果图，描述信号在不同时间和频率的能量强度；基于对图像的分析，以实现储层流体预测，油气检测，薄层辨析，与地震信号相关的数据处理。

2.如权利要求1所述的方法，其特征在于：

为保证时频分析的能量聚集性，高斯窗不可太宽也不可太窄，因此对影响窗函数的标准差作出以下限制：

L<σ(f)<U

其中，L取nTs，U取lTs，n和l分别取10，1000；f∈[f_min,f_max],f_min,f_max取决于待分析信号，那么，对于所述计算时频能量聚集度的四个参数有如下限制：

0≤A,B,C,D≤2。

3.如权利要求1所述的方法，其特征在于，所述的高斯窗函数按如下方式进行处理：

归一化参数

使高斯窗成为一个单位面积，从而使得新的广义时频变换得幅度具有与傅里叶变换相同的意义；

由于高斯窗是关于频率的函数，因此窗口形态会随着信号频率的变化而变化，那么高斯窗在时间域不仅做平移变换，同时也在做伸缩变换，同样也意味着广义时频变换是一种多分辨率的信号分析方法；当[A B C D]＝[0 1 1 0]时，所述的高斯窗即等于标准的高斯窗；

关于参数对于窗口的影响，参数B和D用于粗调窗口形态，参数A和C用于微调窗口形态；

在时间域，窗口长度随着参数A，B，D的增加而变窄，随着参数C的增加而变长；窗口幅度随着参数A和B的增加而增加，随着参数C和D的增加而减小；其中B和D的影响较为剧烈，用于粗调窗口形态，而参数A和C的影响较为微小，用于微调窗口形态。

Images