CN104950335B - Enpemf信号归一化stft‑wvd时频分析方法 - Google Patents

Abstract

本发明提供了一种ENPEMF信号归一化STFT‑WVD时频分析方法，对地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner‑Ville分布，得到STFT数组和WVD数组，然后选取STFT数组中的最大值归一化得到数组STFT_1，记录1值位置及最小值，将其中0值的数用最小值替换，选取WVD数组中相同位置数归一化得到临时数组A，点除以STFT‑1得到临时数组B，将其中大于设置值x的位置的数与WVD数组中对应位置的数全部置为0，输出临时数组B和WVD数组。本发明较好地消除交叉项干扰，沿袭WVD较高的时频分辨率，克服输入信号改变从而需要重新调节阀值与幂调节系数的缺点，结果理想，使用灵活。

Description

ENPEMF信号归一化STFT-WVD时频分析方法

技术领域

本发明涉及一种ENPEMF信号归一化STFT-WVD时频分析方法，属于地震电磁前兆研究与地震预测领域。

背景技术

地震给人类的生活带来了巨大的灾难，据统计，全球的自然灾害之中，地震造成的死亡人数占全部自然灾害死亡人数的54％，堪称自然灾害之最。如何预测地震一直以来都是一个热门敏感的话题。然而，因地震预测有着地球内部的“不可入性”、大地震的“非频发性”、“地震物理过程的复杂性”三大困难，地震预测成为了公认的世界性难题，对于地震前兆预测对于人类生命安全和社会财产安全的具有很大的意义。

在现有的孕震信息研究过程中，STFT与WVD时频分析方法是常用的用来对采集到的大量地球天然脉冲电磁场信号进行分析的方法。

传统的STFT-WVD时频分析方法会用到以下公式对地球天然脉冲电磁场信号进行处理：

SW_x(t,f)＝W_x(t,f){|S_STFTx(t,f)|>c}………………(1)

SW_x(t,f)＝S_STFTx ^a(t,f)W_x ^b(t,f)…………………(2)

传统的STFT-WVD时频分析方法需要设置阀值c以及幂调节系数a和b。其中c为STFT谱的交叉项消除阀值。当STFT谱的值小于该阀值时返回0，如果大于该阀值则返回1。WVD中有交叉项对应STFT谱部分的数值肯定小于该阀值，用数字0与WVD相乘以消除交叉项；其中a和b式为幂调节系数，作用是增强两变换数值较大部分而消弱有交叉项部分。式(1)所示方法灵活性差，输入信号幅值或能量大小直接影响c值的选择，目前没有人能提出一种行之有效的自适应阀值选择方法，并且在实际信号中有用分量幅值、能量大小往往各不相同甚至差别较大，因此采用设置阀值消交叉项容易也将信息项消除；式(2)所示方法有些许改进，但同样幂调节系数的确定没有理论基础，如何根据待分析信号的特征确定幂调节系数有待进一步研究，对使用造成不便，同时存在难于避免的交叉项干扰，往往需要进一步利用滤波等方法对得到的信号进行进一步处理，才能得到便于解读的时频图和谱图。

发明内容

为了解决现有技术的不足，本发明提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，通过结合短时Fourier变换和Wigner-Ville分布各自特点较好地消除了交叉项的干扰，同时沿袭了WVD较高的时频分辨率。

本发明为解决其技术问题所采用的技术方案是：提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，包括以下步骤：

(1)将地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布化，分别得到STFT数组和WVD数组；

(2)选取STFT数组的最大值max_st，将STFT数组中的各个数除以max_st以对STFT数组进行归一化，得到归一化后的数组STFT_1；

(3)记录数组STFT_1数值为1的数所在的位置(i，j)；记录数组STFT_1中非0值的最小值min_1；

(4)将数组STFT_1中的值为0的数全部用min_1替换；

(5)选取WVD数组中位置为(i，j)的数max_wvd，将WVD数组中的各个数除以max_wvd以对WVD数组进行归一化，得到临时数组A；

(6)临时数组A点除以数组STFT_1得到临时数组B，设置矩阵倍数比值上限值x，x的范围为1和2之间；选取临时数组B中大于x的数并记录它们的位置，将临时数组B中大于x的数全部置为0，将WVD数组中与临时数组B中大于x的数的相同位置的数全部置为0；

(7)输出经步骤(6)置零后的临时数组B以及WVD数组。

步骤(6)中矩阵倍数比值上限值x根据地球天然脉冲电磁场信号的幅值强弱进行设置，x随地球天然脉冲电磁场信号的幅值变强而增大。

步骤(7)输出的WVD数组进行二维时频谱显示。

本发明基于其技术方案所具有的有益效果在于：

本发明的ENPEMF信号归一化STFT-WVD时频分析方法结合短时Fourier变换和Wigner-Ville分布各自特点较好地消除了交叉项的干扰，同时沿袭了WVD较高的时频分辨率，与传统的STFT-WVD算法相比，本发明的ENPEMF信号归一化STFT-WVD时频分析方法不需要设置阀值或者调节幂调节系数，克服了输入信号改变从而需要重新调节阀值与幂调节系数的缺点，结果更加理想，使用更加灵活。

附图说明

图1是本发明的流程示意图。

图2是线性调频信号Wigner-Ville分布二维时频能量分布示意图。

图3是线性调频信号短时Fourier变换二维时频谱示意图。

图4是线性调频信号归一化STFT_WVD二维时频能量分布示意图。

图5是第18日通道3NH数据短时Fourier变换二维时频谱示意图。

图6是第18日通道3NH数据Wigner-Ville分布二维时频能量分布示意图。

图7是第18日通道3NH数据归一化STFT_WVD二维时频能量分布。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

本发明提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，结合图1，包括以下步骤：

(1)利用matlab，将地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布化，所述两种变换的函数分别为：

[ST,Ts,F]＝stft(x,Nw,nstep,h,dt)

[tfr,t,f]＝wignerVille(x,fs,fre_bins)

分别得到STFT数组和WVD数组，其中ST和tfr分别表示STFT数组和WVD数组；

(2)选取STFT数组的最大值max_st：

max_st＝max(max(ST))

将STFT数组中的各个数除以max_st以对STFT数组进行归一化，得到归一化后的数组STFT_1：

STFT_1＝ST/max_st

(3)记录数组STFT_1数值为1的数所在的位置(i，j)：

[i,j,v]＝find(STFT_1＝＝1)

记录数组STFT_1中非0值的最小值min_1：

[a,b,c]＝find(STFT_1)

min_1＝min(c)

(4)将数组STFT_1中的值为0的数全部用min_1替换：

STFT_1(STFT_1＝＝0)＝min_1

(5)选取WVD数组中位置为(i，j)的数max_wvd：

max_wvd＝tfr(i,j)

将WVD数组中的各个数除以max_wvd以对WVD数组进行归一化，得到临时数组A：

A＝tfr/max_wvd

(6)临时数组A点除以数组STFT_1得到临时数组B：

B＝A./STFT_1

根据地球天然脉冲电磁场信号的幅值强弱设置阵倍数比值上限值x，x随地球天然脉冲电磁场信号的幅值变强而增大，x可以设置为1.8，选取临时数组B中大于x的数并记录它们的位置，将临时数组B中大于x的数全部置为0：

tfr(B>1.8)＝0

将WVD数组中与临时数组B中大于x的数的相同位置的数全部置为0：

B(B>1.8)＝0

(7)输出经步骤(6)置零后的临时数组B以及WVD数组。

步骤(7)输出的WVD数组进行二维时频谱显示。

对于任何一种时频分布方法，有一个公认的观点：如果该时频分布方法对线性调频信号不能提供好的时频聚集性，那么它便不适合用作非平稳信号时频分析的工具。图2所示为仅经过Wigner-Ville分布化的线性调频信号的二维时频能量分布，图3所示为仅经过短时Fourier变换后的线性调频信号二维时频谱示意图，图4是利用本发明的ENPEMF信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图。与图3相比利用本发明的ENPEMF信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图克服了短时Fourier变换较差的时频聚集特性的缺点，和图2相比本发明的ENPEMF信号归一化STFT-WVD时频分析方法解决了Wigner-Ville分布固有交叉项的干扰的问题，而且在对信号处理前后不需要对参数进行设置、调整，因此本发明的ENPEMF信号归一化STFT-WVD时频分析方法是一种行之有效的时频分析工具。

18日芦山Ms7.0地震前通道3的NH数据(每秒接收的地球天然脉冲电磁场信号脉冲个数)短时Fourier变换的二维时频谱如图5所示，经过Wigner-Ville分布的二维时频能量分布示意图如图6所示，采用本发明的ENPEMF信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图如图7所示。可以看到与使用短时Fourier变换得到的二维时频谱相比采用归一化STFT_WVD方法得到的二维时频能量分布图更加清晰的描述了信号在不同时间和频率上的能量、强度和信号能量集中的时间和频段；与使用Wigner-Ville分布得到的二维时频能量分布图相比采用归一化STFT_WVD方法得到的二维时频能量分布图排除了交叉项信号叠加对真实信息项在不同时间和频率上的能量的反映的干扰，对信号有用信息分量的提取和判断起重要作用。

Claims (3)

1.一种ENPEMF信号归一化STFT-WVD时频分析方法，其特征在于包括以下步骤：

(1)将地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布化，分别得到STFT数组和WVD数组；

(2)选取STFT数组的最大值max_st，将STFT数组中的各个数除以max_st以对STFT数组进行归一化，得到归一化后的数组STFT_1；

(3)记录数组STFT_1数值为1的数所在的位置(i，j)；记录数组STFT_1中非0值的最小值min_1；

(4)将数组STFT_1中的值为0的数全部用min_1替换；

(5)选取WVD数组中位置为(i，j)的数max_wvd，将WVD数组中的各个数除以max_wvd以对WVD数组进行归一化，得到临时数组A；

(6)临时数组A点除以数组STFT_1得到临时数组B，设置矩阵倍数比值上限值x，x的范围为1和2之间；选取临时数组B中大于x的数并记录它们的位置，将临时数组B中大于x的数全部置为0，将WVD数组中与临时数组B中大于x的数的相同位置的数全部置为0；

(7)输出经步骤(6)置零后的临时数组B以及WVD数组。

2.根据权利要求1所述的ENPEMF信号归一化STFT-WVD时频分析方法，其特征在于：步骤(6)中矩阵倍数比值上限值x根据地球天然脉冲电磁场信号的幅值强弱进行设置，x随地球天然脉冲电磁场信号的幅值变强而增大。

3.根据权利要求1所述的ENPEMF信号归一化STFT-WVD时频分析方法，其特征在于：步骤(7)输出的WVD数组进行二维时频谱显示。