Enpemf信号归一化stft-wvd时频分析方法
【技术领域】
[0001] 本发明设及一种ENPEMF信号归一化STFT-WVD时频分析方法，属于地震电磁前兆 研究与地震预测领域。
【背景技术】
[0002] 地震给人类的生活带来了巨大的灾难，据统计，全球的自然灾害之中，地震造成的 死亡人数占全部自然灾害死亡人数的54%，堪称自然灾害之最。如何预测地震一直W来都 是一个热口敏感的话题。然而，因地震预测有着地球内部的"不可入性"、大地震的"非频发 性"、"地震物理过程的复杂性大困难，地震预测成为了公认的世界性难题，对于地震前 兆预测对于人类生命安全和社会财产安全的具有很大的意义。
[0003] 在现有的孕震信息研究过程中，STFT与WVD时频分析方法是常用的用来对采集到 的大量地球天然脉冲电磁场信号进行分析的方法。
[0004] 传统的STFT-WVD时频分析方法会用到W下公式对地球天然脉冲电磁场信号进行 处理：
[0007] 传统的STFT-WVD时频分析方法需要设置阀值CW及幕调节系数a和b。其中C为 STFT谱的交叉项消除阀值。当STFT谱的值小于该阀值时返回0,如果大于该阀值则返回1。WVD中有交叉项对应STFT谱部分的数值肯定小于该阀值，用数字0与WVD相乘W消除交叉 项；其中a和b式为幕调节系数，作用是增强两变换数值较大部分而消弱有交叉项部分。式 (1)所示方法灵活性差，输入信号幅值或能量大小直接影响C值的选择，目前没有人能提出 一种行之有效的自适应阀值选择方法，并且在实际信号中有用分量幅值、能量大小往往各 不相同甚至差别较大，因此采用设置阀值消交叉项容易也将信息项消除；式（2)所示方法 有些许改进，但同样幕调节系数的确定没有理论基础，如何根据待分析信号的特征确定幕 调节系数有待进一步研究，对使用造成不便，同时存在难于避免的交叉项干扰，往往需要进 一步利用滤波等方法对得到的信号进行进一步处理，才能得到便于解读的时频图和谱图。
【发明内容】
[000引为了解决现有技术的不足，本发明提供了一种ENPEMF信号归一化STFT-WVD时频 分析方法，通过结合短时化urier变换和Wigner-Vi 1 le分布各自特点较好地消除了交叉项 的干扰，同时沿袭了WVD较高的时频分辨率。
[0009] 本发明为解决其技术问题所采用的技术方案是；提供了一种ENPEMF信号归一化 STFT-WVD时频分析方法，包括W下步骤：
[0010] (1)将地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布 化，分别得到STFT数组和WVD数组；
[001U(2)选取STFT数组的最大值max_st，将STFT数组中的各个数除wmax_stw对STFT数组进行归一化，得到归一化后的数组STFT_1;
[0012] (3)记录数组STFT_1数值为1的数所在的位置（i，j);记录数组STFT_1中非0值 的最小值min_l;
[001引 （4)将数组STFT_1中的值为0的数全部用min_l替换；
[0014] 妨选取WVD数组中位置为的数max_wvd，将WVD数组中的各个数除Wmax_ wvdW对WVD数组进行归一化，得到临时数组A;
[00巧]做临时数组A点除W数组STFT_1得到临时数组B，设置矩阵倍数比值上限值X，X的范围为1和2之间；选取临时数组B中大于X的数并记录它们的位置，将临时数组B中 大于X的数全部置为0,将WVD数组中与临时数组B中大于X的数的相同位置的数全部置为 0 ；
[0016] (7)输出经步骤做置零后的临时数组BW及WVD数组。
[0017] 步骤做中矩阵倍数比值上限值X根据地球天然脉冲电磁场信号的幅值强弱进行 设置，X随地球天然脉冲电磁场信号的幅值变强而增大。
[0018] 步骤（7)输出的WVD数组进行二维时频谱显示。
[0019] 本发明基于其技术方案所具有的有益效果在于：
[0020] 本发明的ENPEMF信号归一化STFT-WVD时频分析方法结合短时化urier变换和 Wigner-Ville分布各自特点较好地消除了交叉项的干扰，同时沿袭了WVD较高的时频分辨 率，与传统的STFT-WVD算法相比，本发明的ENPEMF信号归一化STFT-WVD时频分析方法不 需要设置阀值或者调节幕调节系数，克服了输入信号改变从而需要重新调节阀值与幕调节 系数的缺点，结果更加理想，使用更加灵活。
【附图说明】
[0021] 图1是本发明的流程示意图。
[0022] 图2是线性调频信号Wigner-Ville分布二维时频能量分布示意图。
[0023] 图3是线性调频信号短时化urier变换二维时频谱示意图。
[0024] 图4是线性调频信号归一化STFT_WVD二维时频能量分布示意图。
[0025] 图5是第18日通道3NH数据短时化urier变换二维时频谱示意图。
[0026] 图6是第18日通道3NH数据Wigner-Ville分布二维时频能量分布示意图。
[0027] 图7是第18日通道3NH数据归一化STFT_WVD二维时频能量分布。
【具体实施方式】
[002引下面结合附图和实施例对本发明作进一步说明。
[0029] 本发明提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，结合图1，包括W 下步骤：
[0030] (1)利用matl油，将地球天然脉冲电磁场信号分别作短时化urier变换和 Wigner-Ville分布化，所述两种变换的函数分别为：
[0031] [ST,Ts,円=stft(X,Nw,nstep,h,dt)
[0032] [tfr,t,f] =wignerVille(X,fs,fre_bins)
[0033] 分别得到STFT数组和WVD数组，其中ST和t化分别表示STFT数组和WVD数组；
[0034] 似选取STFT数组的最大值max_st;
[0035]max_st=max(max(ST))
[0036] 将STFT数组中的各个数除Wmax_stW对STFT数组进行归一化，得到归一化后的 数组STFT_1;
[0037]STFT_1 =ST/max_st
[0038] (3)记录数组STFT_1数值为1的数所在的位置（i，j);
[0039] [i,j,v] =find(STFT_l== 1)
[0040] 记录数组STFT_1中非0值的最小值min_l; 闺][a,b,C]=find(STFT-l)
[0042]min_l=mink)
[00创 （4)将数组STFT_1中的值为0的数全部用min_l替换；
[0044]STFT-1 (STFT-1 = =0)=min_l
[0045] 妨选取WVD数组中位置为（i，_]'）的数max_wvd;
[0046] max_wvd=tfr(i,j)
[0047] 将WVD数组中的各个数除Wmax_wvd W对WVD数组进行归一化，得到临时数组A :
[0048]A=tfr/max_wvd
[0049] 做临时数组A点除W数组STFT_1得到临时数组B:
[0050] B=A. /STFT-1
[0化1] 根据地球天然脉冲电磁场信号的幅值强弱设置阵倍数比值上限值X，X随地球天 然脉冲电磁场信号的幅值变强而增大，X可W设置为1. 8,选取临时数组B中大于X的数并 记录它们的位置，将临时数组B中大于X的数全部置为0 :
[0化2] tfHB〉1.8)=0
[0053]将WVD数组中与临时数组B中大于X的数的相同位置的数全部置为0 :
[0化4] B炬〉1.8)=0
[005引 （7)输出经步骤做置零后的临时数组BW及WVD数组。
[0056] 步骤（7)输出的WVD数组进行二维时频谱显示。
[0057] 对于任何一种时频分布方法，有一个公认的观点：如果该时频分布方法对线性调 频信号不能提供好的时频聚集性，那么它便不适合用作非平稳信号时频分析的工具。图2 所示为仅经过Wigner-Ville分布化的线性调频信号的二维时频能量分布，图3所示为仅 经过短时化urier变换后的线性调频信号二维时频谱示意图，图4是利用本发明的ENPEMF 信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图。与图2相比利用本 发明的ENPEMF信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图克服 了短时化urier变换较差的时频聚集特性的缺点，和图3相比本发明的ENPEMF信号归一化 STFT-WVD时频分析方法解决了Wigner-Ville分布固有交叉项的干扰的问题，而且在对信 号处理前后不需要对参数进行设置、调整，因此本发明的ENPEMF信号归一化STFT-WVD时频 分析方法是一种行之有效的时频分析工具。
[0化8] 18日芦山Ms7.0地震前通道3的NH数据（每秒接收的地球天然脉冲电磁场信号 脉冲个数）短时化urier变换的二维时频谱如图5所示，经过Wigner-Ville分布的二维时 频能量分布示意图如图6所示，采用本发明的ENPEMF信号归一化STFT-WVD时频分析方法 得到的二维时频能量分布示意图如图7所示。可W看到与使用短时化urier变换得到的二 维时频谱相比采用归一化STFT_WVD方法得到的二维时频能量分布图更加清晰的描述了信 号在不同时间和频率上的能量、强度和信号能量集中的时间和频段；与使用Wigner-Ville 分布得到的二维时频能量分布图相比采用归一化STFT_WVD方法得到的二维时频能量分布 图排除了交叉项信号叠加对真实信息项在不同时间和频率上的能量的反映的干扰，对信号 有用信息分量的提取和判断起重要作用。
【主权项】
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_l ; (4) 将数组STFT_1中的值为0的数全部用min_l替换； (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数组进行二维时频谱显示。
【专利摘要】本发明提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，对地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布，得到STFT数组和WVD数组，然后选取STFT数组中的最大值归一化得到数组STFT_1，记录1值位置及最小值，将其中0值的数用最小值替换，选取WVD数组中相同位置数归一化得到临时数组A，点除以STFT-1得到临时数组B，将其中大于设置值x的位置的数与WVD数组中对应位置的数全部置为0，输出临时数组B和WVD数组。本发明较好地消除交叉项干扰，沿袭WVD较高的时频分辨率，克服输入信号改变从而需要重新调节阀值与幂调节系数的缺点，结果理想，使用灵活。
【IPC分类】G01V3/08, G06F17/14, G01R23/16
【公开号】CN104950335
【申请号】CN201510209896
【发明人】郝国成, 陈忠昌, 赵娟, 康坊, 白雨晓, 宋佳珍
【申请人】中国地质大学（武汉）
【公开日】2015年9月30日
【申请日】2015年4月28日