一种基于挤压加窗傅里叶变换的高铁震源地震信号时频分析方法

技术领域

本发明属于勘探地球物理领域，特别涉及一种基于挤压加窗傅里叶变换的高铁震源地震信号时频分析方法。

背景技术

截止2017年底，我国高铁营业里程达2.5万公里，占世界高铁总量的66％。每天有数千趟高铁列车高速运行在分布范围很广的高铁线路上。如此庞大的数量的高铁列车高速运行在高铁线路上，不但会引起高铁列车及路基的振动，而且会将振动以各种类型的地震波传播出去。在高铁线路两侧几十米范围上埋置检波器可接收到高铁列车运行所引起的地震波，即高铁震源地震信号。对检波器接收到的信号进行分析，不但可能分析高铁列车的运行状态，而且有望对高铁线路附近的地下结构进行成像。然而，面对这种全新的高铁震源地震数据，采用何种频谱分析手段来检测频谱成分的变化是极为关键的。目前针对分析这种高铁运行引起地震信号手段极为有限，主要包含：

现有技术1：离散傅里叶变换

该类方法对检波器接收到的数字信号进行快速傅里叶变换或者离散傅里叶变换，可获得信号的振幅谱。

现有技术1的特点：

优点：简单易行、计算量小且不受人为因素干扰。

缺点：1、只能分析频率成分，无法分析各种频率成分的起始和终止时刻；2、无法分析各种频率成分随时间的变化规律。

现有技术2：短时傅里叶变换

该类方法在每个时间点利用窗函数截取一段信号，然后对截取后的信号进行傅里叶变化以获得局部的频率成分。

现有技术2的特点：

优点：实现较为简单、计算量小且不受人为因素干扰。

缺点：频率分辨率较差，受限于不确定性原理的约束。

现有技术3：连续小波变换

该类方法把母小波进行伸缩和平移形成一系列小波族，然后将小波族与信号以此做內积得到一系列连续小波变换系数。

现有技术3的特点：

优点：较为简单、计算量较小。

缺点：1、频率分辨率较差，受限于不确定性原理的约束；2、母小波的选择较为困难；3、所形成的结果为时间-尺度域，而尺度需要转换才能变成频率。

发明内容

本发明的目的在于提供一种基于挤压加窗傅里叶变换的高铁震源地震信号时频分析方法，以解决上述技术问题。本发明能够分析高铁震源地震信号的频率成分随时间的变化，并采用挤压加窗傅里叶变换得到高铁震源地震信号的时频分布，为后续判断列车运行状态提供数据。

为了实现上述目的，本发明采用如下技术方案：

一种基于挤压加窗傅里叶变换的高铁震源地震信号时频分析方法，包括以下步骤：

步骤01：对采集的高铁震源地震信号进行进行去均值处理；

步骤02：对步骤01进行去均值处理后的高铁震源地震信号进行傅里叶变换以获得信号的能量谱分布；

步骤03：利用信号的能量谱分布自适应地确定插值倍数；

步骤04：根据插值倍数对经过步骤1获得的去均值高铁震源地震信号进行插值以获得更高采样率的地震信号；

步骤05：利用窗函数对插值后的地震信号进行加窗傅里叶变换；

步骤06：计算信号加窗傅里叶变换系数相位的关于时间差分并获得挤压位置，然后将信号加窗傅里叶变换系数累积到挤压位置；

步骤07：获得最终的时频分析结果并求模值获得最终的时频分布。

进一步的，步骤01具体包括：

用s[m]表示一维地震信号，共有M个采样点，其时间采样间隔为Δt，m表示信号在时间方向的指标；通过如下方式去均值：

进一步的，步骤02具体包括：

对信号s[m]进行离散傅里叶变换得到S[k]为：

其中，M表示一维信号的采样点数，而Δf为频率域采样间隔，且k表示频率指标，其范围从0到M-1；由S[kΔf]得该一维信号的沿频率域的能量谱EF[k]为：

进一步的，步骤03具体包括：

在频率域计算该一维信号的能量E:

然后计算EF[k]的能量累积函数ACCU_EF[q]

在能量累积函数ACCU_EF[q]中按照如下准则寻找频率上限指标Q：

其中λ为阈值，取值大于或等于0.999；然后根据如下公式确定插值倍数R：

其中表示上取整。

进一步的，步骤04具体包括：

构造一个新的序列SS[k]，共有RM个点，其与S[k]的关系如下：

然后对该新序列SS[k]做反傅里叶变换，可得到重采样后的新时间序列ss[m]：

其中real{}表示取复数的实部。

进一步的，步骤05具体包括：

选取的窗函数为g[m]，对插值后地震信号ss[m]做加窗傅里叶变换得到结果WFT[m,k]为：

其中l为临时时间指标，m为时间指标，k为频率指标。同时假设最终的挤压变换结果为

同时假设最终的挤压变换结果为SWFT[m,k]，并将其全部初始化为0。

进一步的，步骤06具体包括：

WFT_phs[m,k]表示加窗傅里叶变换系数WFT[m,k]的相位，相位WFT_phs[m,k]沿m指标计算差分可得到:

其中

其中img()表示取复数的虚部，real()表示取复数的实部，abs()表示取复数模值。对dif[m,k]通过如下方式得到对应的挤压位置k₁:

其中round()表示对浮点数取整；将相应的加窗傅里叶变换系数WFT[m,k]累积到新的位置：

SWFT[m,k₁]＝SWFT[m,k₁]+WFT[m,k]。

进一步的，步骤07具体包括：

对挤压变换结果SWFT[m,k]计算模值得到高铁震源地震信号的时频分布如下：

SWFT_E[m,k]＝abs(SWFT[m,k])。

进一步的，根据步骤07获得的时频分布，检测高铁列车的运动状态。

相对于现有技术，本发明具有以下有益效果：本发明为针对高铁运行所引起的地震信号的一种时频分析方法，主要用于刻画高铁震源地震信号的频率成分随时间的变化，为一种快速数据处理方法；本发明首先对高铁震源地震信号进行去零均值处理，然后对信号进行频谱分析确定信号插值倍数并插值以提高差分精度，最后对信号的加窗傅里叶变换系数在频率方向上进行挤压操作以获得高精度的时频分布。相比较于常规时频分布，本发明获得的高精度时频分布可以精确刻画高铁震源地震信号中频率成分随时间的变化，可用于检测高铁列车的运行状态(匀速或者加速)。

附图表说明

图1为本发明流程图；

图2为列车1经过时单个检波器接收到的一道高铁震源地震信号；

图3为取均值后的高铁震源信号；

图4为去均值后高铁震源信号的振幅谱；

图5为去均值后高铁震源信号的能量谱；

图6为能量累积函数；

图7为插值后的一道高铁震源信号；

图8为利用加窗傅里叶变换得到的列车1经过时单个检波器接收到的高铁震源地震信号的时频分布。

图9为利用挤压加窗傅里叶变换得到的列车1经过时单个检波器接收到的高铁震源地震信号的时频分布。

图10为利用挤压加窗傅里叶变换得到的列车2经过时单个检波器接收到的高铁震源地震信号的时频分布。

具体实施方式

下面结合附图及对本发明做进一步详细的说明。

本发明为针对高铁经过时所引起地震信号的一种时频分析方法。本发明首先对采集到的地震信号进行傅里叶变换以获得信号的能量谱，然后利用能量谱确定信号的插值倍数并对信号进行插值，接着计算插值后信号的加窗傅里叶变换，利用信号加窗傅里叶变换的相位谱差分确定频率成分的挤压位置，最后将相应的系数累积到挤压位置获得高精度的时频谱图，为分析各时间段频率成分提供一个高精度的时频谱图。

请参阅图1所示，本发明提供一种基于挤压加窗傅里叶变换的高铁震源地震信号时频分析方法，包括以下步骤：

步骤01：对采集的高铁震源地震信号进行进行去均值处理；

步骤02：对步骤01进行去均值处理后的高铁震源地震信号进行傅里叶变换，获得信号的能量谱分布；

步骤03：利用信号的能量谱分布自适应地确定插值倍数；

步骤04：根据插值倍数对经过步骤1获得的去均值高铁震源地震信号进行插值以获得更高采样率的地震信号；

步骤05：利用窗函数对插值后的地震信号进行加窗傅里叶变换；

步骤06：计算信号加窗傅里叶变换系数相位的关于时间差分并获得挤压位置，然后将信号加窗傅里叶变换系数累积到挤压位置；

步骤07：获得最终的时频分析结果并求模值获得最终的时频分布。

步骤01：对采集的高铁震源地震信号进行去均值处理，具体包括：

用s[m]表示一维地震信号，共有M个采样点，其时间采样间隔为Δt，m表示信号在时间方向的指标。通过如下方式去均值：

步骤02中去均值处理后的高铁震源地震信号进行傅里叶变换以获得信号的能量谱分布，具体包括：

对信号s[m]进行离散傅里叶变换得到S[k]为：

其中，M表示一维信号的采样点数，而Δf为频率域采样间隔，且k表示频率指标，其范围从0到M-1。由S[kΔf]得该一维信号的沿频率域的能量谱EF[k]为：

步骤03中利用信号的能量谱分布自适应地确定插值倍数，具体包括：

在频率域计算该一维信号的能量E:

然后计算EF[k]的能量累积函数ACCU_EF[q]

在能量累积函数ACCU_EF[q]中按照如下准则寻找频率上限指标Q：

其中λ为阈值，取值一般选取为大于0.999的值。然后根据如下公式确定插值倍数R：

其中表示上取整。

步骤04中根据插值倍数对经过步骤1获得的去均值高铁震源地震信号进行插值以获得更高采样率的地震信号，具体包括：

构造一个新的序列SS[k]，共有RM个点，其与S[k]的关系如下：

然后对该新序列SS[k]做反傅里叶变换，得到重采样后的新时间序列ss[m]：

其中real{}表示取复数的实部。

步骤05中窗函数对插值后的地震信号进行加窗傅里叶变换，具体包括：

若选取的窗函数为g[m]，则对插值后地震信号ss[m]做加窗傅里叶变换得到结果WFT[m,k]为：

其中l为临时时间指标，m为时间指标，k为频率指标。同时假设最终的挤压变换结果为SWFT[m,k]，并将SWFT[m,k]全部初始化为0。

步骤06中计算信号加窗傅里叶变换系数相位的关于时间差分并获得挤压位置，然后将信号加窗傅里叶变换系数累积到挤压位置，具体包括：

WFT_phs[m,k]表示加窗傅里叶变换系数相位，WFT_phs[m,k]沿m指标计算差分得到:

其中phs_img_dif[m,k]及phs_real_dif[m,k]如下：

其中img()表示取复数的虚部，real()表示取复数的实部，abs()表示取复数模值。对dif[m,k]通过如下方式得到对应的挤压频率位置k₁:

其中round()表示对浮点数取整。将相应的加窗傅里叶变换系数WFT[m,k]累积到新的位置[m,k₁]：

SWFT[m,k₁]＝SWFT[m,k₁]+WFT[m,k]

步骤07中获得最终的时频分析结果并求模值获得最终的时频分布，具体如下：

对挤压变换结果SWFT[m,k]计算模值得到高铁震源地震信号的时频分布如下：

SWFT_E[m,k]＝abs(SWFT[m,k])

以高铁经过时距高铁线路30m的单个低频检波器所接收到的信号为例。图2为列车1经过时单个检波器所接收到的高铁震源所引起的振动信号，采样间隔为5ms，共有3001个采样点。图3为去均值后的结果，图4为对应的振幅谱，图5为对应的能量谱，图6为能量累积函数。取为λ为0.999，由图6所示的能量谱得Q为1147，由此可以得到插值倍数R为14倍。图7为14倍插值后的信号，而图8为列车1经过时利用常规加窗傅里叶变换得到的高铁震源地震信号的时频分布，而图9为列车1经过时利用挤压加窗傅里叶变换得到的高铁震源地震信号的时频分布，各个分立谱的频率随时间保持不变，说明列车经过检波器时处于匀速行驶状态。图10为列车经过时基于挤压加窗傅里叶变换的高铁震源地震信号时频分布，各个分立谱的频率随时间增加而变大，说明列车经过检波器时处于加速行驶状态。

最后需要说明的是，以上算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等(例如更改窗函数、替换插值方案、替换差分方法等)，均应在本发明的保护范围内。