Enpemf信號(hào)歸一化stft-wvd時(shí)頻分析方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明設(shè)及一種ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法，屬于地震電磁前兆 研究與地震預(yù)測(cè)領(lǐng)域。
【背景技術(shù)】
[0002] 地震給人類的生活帶來(lái)了巨大的災(zāi)難，據(jù)統(tǒng)計(jì)，全球的自然災(zāi)害之中，地震造成的 死亡人數(shù)占全部自然災(zāi)害死亡人數(shù)的54%，堪稱自然災(zāi)害之最。如何預(yù)測(cè)地震一直W來(lái)都 是一個(gè)熱口敏感的話題。然而，因地震預(yù)測(cè)有著地球內(nèi)部的"不可入性"、大地震的"非頻發(fā) 性"、"地震物理過(guò)程的復(fù)雜性大困難，地震預(yù)測(cè)成為了公認(rèn)的世界性難題，對(duì)于地震前 兆預(yù)測(cè)對(duì)于人類生命安全和社會(huì)財(cái)產(chǎn)安全的具有很大的意義。
[0003] 在現(xiàn)有的孕震信息研究過(guò)程中，STFT與WVD時(shí)頻分析方法是常用的用來(lái)對(duì)采集到 的大量地球天然脈沖電磁場(chǎng)信號(hào)進(jìn)行分析的方法。
[0004] 傳統(tǒng)的STFT-WVD時(shí)頻分析方法會(huì)用到W下公式對(duì)地球天然脈沖電磁場(chǎng)信號(hào)進(jìn)行 處理：
[0007] 傳統(tǒng)的STFT-WVD時(shí)頻分析方法需要設(shè)置閥值CW及幕調(diào)節(jié)系數(shù)a和b。其中C為 STFT譜的交叉項(xiàng)消除閥值。當(dāng)STFT譜的值小于該閥值時(shí)返回0,如果大于該閥值則返回1。WVD中有交叉項(xiàng)對(duì)應(yīng)STFT譜部分的數(shù)值肯定小于該閥值，用數(shù)字0與WVD相乘W消除交叉 項(xiàng)；其中a和b式為幕調(diào)節(jié)系數(shù)，作用是增強(qiáng)兩變換數(shù)值較大部分而消弱有交叉項(xiàng)部分。式 (1)所示方法靈活性差，輸入信號(hào)幅值或能量大小直接影響C值的選擇，目前沒(méi)有人能提出 一種行之有效的自適應(yīng)閥值選擇方法，并且在實(shí)際信號(hào)中有用分量幅值、能量大小往往各 不相同甚至差別較大，因此采用設(shè)置閥值消交叉項(xiàng)容易也將信息項(xiàng)消除；式（2)所示方法 有些許改進(jìn)，但同樣幕調(diào)節(jié)系數(shù)的確定沒(méi)有理論基礎(chǔ)，如何根據(jù)待分析信號(hào)的特征確定幕 調(diào)節(jié)系數(shù)有待進(jìn)一步研究，對(duì)使用造成不便，同時(shí)存在難于避免的交叉項(xiàng)干擾，往往需要進(jìn) 一步利用濾波等方法對(duì)得到的信號(hào)進(jìn)行進(jìn)一步處理，才能得到便于解讀的時(shí)頻圖和譜圖。

【發(fā)明內(nèi)容】

[000引為了解決現(xiàn)有技術(shù)的不足，本發(fā)明提供了一種ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻 分析方法，通過(guò)結(jié)合短時(shí)化urier變換和Wigner-Vi 1 le分布各自特點(diǎn)較好地消除了交叉項(xiàng) 的干擾，同時(shí)沿襲了WVD較高的時(shí)頻分辨率。
[0009] 本發(fā)明為解決其技術(shù)問(wèn)題所采用的技術(shù)方案是；提供了一種ENPEMF信號(hào)歸一化 STFT-WVD時(shí)頻分析方法，包括W下步驟：
[0010] (1)將地球天然脈沖電磁場(chǎng)信號(hào)分別作短時(shí)Fourier變換和Wigner-Ville分布 化，分別得到STFT數(shù)組和WVD數(shù)組；
[001U(2)選取STFT數(shù)組的最大值max_st，將STFT數(shù)組中的各個(gè)數(shù)除wmax_stw對(duì)STFT數(shù)組進(jìn)行歸一化，得到歸一化后的數(shù)組STFT_1;
[0012] (3)記錄數(shù)組STFT_1數(shù)值為1的數(shù)所在的位置（i，j);記錄數(shù)組STFT_1中非0值 的最小值min_l;
[001引 （4)將數(shù)組STFT_1中的值為0的數(shù)全部用min_l替換；
[0014] 妨選取WVD數(shù)組中位置為的數(shù)max_wvd，將WVD數(shù)組中的各個(gè)數(shù)除Wmax_ wvdW對(duì)WVD數(shù)組進(jìn)行歸一化，得到臨時(shí)數(shù)組A;
[00巧]做臨時(shí)數(shù)組A點(diǎn)除W數(shù)組STFT_1得到臨時(shí)數(shù)組B，設(shè)置矩陣倍數(shù)比值上限值X，X的范圍為1和2之間；選取臨時(shí)數(shù)組B中大于X的數(shù)并記錄它們的位置，將臨時(shí)數(shù)組B中 大于X的數(shù)全部置為0,將WVD數(shù)組中與臨時(shí)數(shù)組B中大于X的數(shù)的相同位置的數(shù)全部置為 0 ；
[0016] (7)輸出經(jīng)步驟做置零后的臨時(shí)數(shù)組BW及WVD數(shù)組。
[0017] 步驟做中矩陣倍數(shù)比值上限值X根據(jù)地球天然脈沖電磁場(chǎng)信號(hào)的幅值強(qiáng)弱進(jìn)行 設(shè)置，X隨地球天然脈沖電磁場(chǎng)信號(hào)的幅值變強(qiáng)而增大。
[0018] 步驟（7)輸出的WVD數(shù)組進(jìn)行二維時(shí)頻譜顯示。
[0019] 本發(fā)明基于其技術(shù)方案所具有的有益效果在于：
[0020] 本發(fā)明的ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法結(jié)合短時(shí)化urier變換和 Wigner-Ville分布各自特點(diǎn)較好地消除了交叉項(xiàng)的干擾，同時(shí)沿襲了WVD較高的時(shí)頻分辨 率，與傳統(tǒng)的STFT-WVD算法相比，本發(fā)明的ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法不 需要設(shè)置閥值或者調(diào)節(jié)幕調(diào)節(jié)系數(shù)，克服了輸入信號(hào)改變從而需要重新調(diào)節(jié)閥值與幕調(diào)節(jié) 系數(shù)的缺點(diǎn)，結(jié)果更加理想，使用更加靈活。
【附圖說(shuō)明】
[0021] 圖1是本發(fā)明的流程示意圖。
[0022] 圖2是線性調(diào)頻信號(hào)Wigner-Ville分布二維時(shí)頻能量分布示意圖。
[0023] 圖3是線性調(diào)頻信號(hào)短時(shí)化urier變換二維時(shí)頻譜示意圖。
[0024] 圖4是線性調(diào)頻信號(hào)歸一化STFT_WVD二維時(shí)頻能量分布示意圖。
[0025] 圖5是第18日通道3NH數(shù)據(jù)短時(shí)化urier變換二維時(shí)頻譜示意圖。
[0026] 圖6是第18日通道3NH數(shù)據(jù)Wigner-Ville分布二維時(shí)頻能量分布示意圖。
[0027] 圖7是第18日通道3NH數(shù)據(jù)歸一化STFT_WVD二維時(shí)頻能量分布。
【具體實(shí)施方式】
[002引下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明作進(jìn)一步說(shuō)明。
[0029] 本發(fā)明提供了一種ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法，結(jié)合圖1，包括W 下步驟：
[0030] (1)利用matl油，將地球天然脈沖電磁場(chǎng)信號(hào)分別作短時(shí)化urier變換和 Wigner-Ville分布化，所述兩種變換的函數(shù)分別為：
[0031] [ST,Ts,円=stft(X,Nw,nstep,h,dt)
[0032] [tfr,t,f] =wignerVille(X,fs,fre_bins)
[0033] 分別得到STFT數(shù)組和WVD數(shù)組，其中ST和t化分別表示STFT數(shù)組和WVD數(shù)組；
[0034] 似選取STFT數(shù)組的最大值max_st;
[0035]max_st=max(max(ST))
[0036] 將STFT數(shù)組中的各個(gè)數(shù)除Wmax_stW對(duì)STFT數(shù)組進(jìn)行歸一化，得到歸一化后的 數(shù)組STFT_1;
[0037]STFT_1 =ST/max_st
[0038] (3)記錄數(shù)組STFT_1數(shù)值為1的數(shù)所在的位置（i，j);
[0039] [i,j,v] =find(STFT_l== 1)
[0040] 記錄數(shù)組STFT_1中非0值的最小值min_l; 閨][a,b,C]=find(STFT-l)
[0042]min_l=mink)
[00創(chuàng) （4)將數(shù)組STFT_1中的值為0的數(shù)全部用min_l替換；
[0044]STFT-1 (STFT-1 = =0)=min_l
[0045] 妨選取WVD數(shù)組中位置為（i，_]'）的數(shù)max_wvd;
[0046] max_wvd=tfr(i,j)
[0047] 將WVD數(shù)組中的各個(gè)數(shù)除Wmax_wvd W對(duì)WVD數(shù)組進(jìn)行歸一化，得到臨時(shí)數(shù)組A :
[0048]A=tfr/max_wvd
[0049] 做臨時(shí)數(shù)組A點(diǎn)除W數(shù)組STFT_1得到臨時(shí)數(shù)組B:
[0050] B=A. /STFT-1
[0化1] 根據(jù)地球天然脈沖電磁場(chǎng)信號(hào)的幅值強(qiáng)弱設(shè)置陣倍數(shù)比值上限值X，X隨地球天 然脈沖電磁場(chǎng)信號(hào)的幅值變強(qiáng)而增大，X可W設(shè)置為1. 8,選取臨時(shí)數(shù)組B中大于X的數(shù)并 記錄它們的位置，將臨時(shí)數(shù)組B中大于X的數(shù)全部置為0 :
[0化2] tfHB〉1.8)=0
[0053]將WVD數(shù)組中與臨時(shí)數(shù)組B中大于X的數(shù)的相同位置的數(shù)全部置為0 :
[0化4] B炬〉1.8)=0
[005引 （7)輸出經(jīng)步驟做置零后的臨時(shí)數(shù)組BW及WVD數(shù)組。
[0056] 步驟（7)輸出的WVD數(shù)組進(jìn)行二維時(shí)頻譜顯示。
[0057] 對(duì)于任何一種時(shí)頻分布方法，有一個(gè)公認(rèn)的觀點(diǎn)：如果該時(shí)頻分布方法對(duì)線性調(diào) 頻信號(hào)不能提供好的時(shí)頻聚集性，那么它便不適合用作非平穩(wěn)信號(hào)時(shí)頻分析的工具。圖2 所示為僅經(jīng)過(guò)Wigner-Ville分布化的線性調(diào)頻信號(hào)的二維時(shí)頻能量分布，圖3所示為僅 經(jīng)過(guò)短時(shí)化urier變換后的線性調(diào)頻信號(hào)二維時(shí)頻譜示意圖，圖4是利用本發(fā)明的ENPEMF 信號(hào)歸一化STFT-WVD時(shí)頻分析方法得到的二維時(shí)頻能量分布示意圖。與圖2相比利用本 發(fā)明的ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法得到的二維時(shí)頻能量分布示意圖克服 了短時(shí)化urier變換較差的時(shí)頻聚集特性的缺點(diǎn)，和圖3相比本發(fā)明的ENPEMF信號(hào)歸一化 STFT-WVD時(shí)頻分析方法解決了Wigner-Ville分布固有交叉項(xiàng)的干擾的問(wèn)題，而且在對(duì)信 號(hào)處理前后不需要對(duì)參數(shù)進(jìn)行設(shè)置、調(diào)整，因此本發(fā)明的ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻 分析方法是一種行之有效的時(shí)頻分析工具。
[0化8] 18日蘆山Ms7.0地震前通道3的NH數(shù)據(jù)（每秒接收的地球天然脈沖電磁場(chǎng)信號(hào) 脈沖個(gè)數(shù)）短時(shí)化urier變換的二維時(shí)頻譜如圖5所示，經(jīng)過(guò)Wigner-Ville分布的二維時(shí) 頻能量分布示意圖如圖6所示，采用本發(fā)明的ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法 得到的二維時(shí)頻能量分布示意圖如圖7所示?？蒞看到與使用短時(shí)化urier變換得到的二 維時(shí)頻譜相比采用歸一化STFT_WVD方法得到的二維時(shí)頻能量分布圖更加清晰的描述了信 號(hào)在不同時(shí)間和頻率上的能量、強(qiáng)度和信號(hào)能量集中的時(shí)間和頻段；與使用Wigner-Ville 分布得到的二維時(shí)頻能量分布圖相比采用歸一化STFT_WVD方法得到的二維時(shí)頻能量分布 圖排除了交叉項(xiàng)信號(hào)疊加對(duì)真實(shí)信息項(xiàng)在不同時(shí)間和頻率上的能量的反映的干擾，對(duì)信號(hào) 有用信息分量的提取和判斷起重要作用。
【主權(quán)項(xiàng)】
1. 一種ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法，其特征在于包括以下步驟： (1) 將地球天然脈沖電磁場(chǎng)信號(hào)分別作短時(shí)Fourier變換和Wigner-Ville分布化，分 別得到STFT數(shù)組和WVD數(shù)組； (2) 選取STFT數(shù)組的最大值max_st，將STFT數(shù)組中的各個(gè)數(shù)除以max_st以對(duì)STFT 數(shù)組進(jìn)行歸一化，得到歸一化后的數(shù)組STFT_1 ; (3) 記錄數(shù)組STFT_1數(shù)值為1的數(shù)所在的位置（i，j);記錄數(shù)組STFT_1中非0值的最 小值min_l ; (4) 將數(shù)組STFT_1中的值為0的數(shù)全部用min_l替換； (5) 選取WVD數(shù)組中位置為（i，j)的數(shù)max_wvd，將WVD數(shù)組中的各個(gè)數(shù)除以max_wvd 以對(duì)WVD數(shù)組進(jìn)行歸一化，得到臨時(shí)數(shù)組A ; (6) 臨時(shí)數(shù)組A點(diǎn)除以數(shù)組STFT_1得到臨時(shí)數(shù)組B，設(shè)置矩陣倍數(shù)比值上限值x，X的 范圍為1和2之間；選取臨時(shí)數(shù)組B中大于X的數(shù)并記錄它們的位置，將臨時(shí)數(shù)組B中大于 X的數(shù)全部置為0,將WVD數(shù)組中與臨時(shí)數(shù)組B中大于X的數(shù)的相同位置的數(shù)全部置為0 ; (7) 輸出經(jīng)步驟（6)置零后的臨時(shí)數(shù)組B以及WVD數(shù)組。2. 根據(jù)權(quán)利要求1所述的ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法，其特征在于：步 驟（6)中矩陣倍數(shù)比值上限值X根據(jù)地球天然脈沖電磁場(chǎng)信號(hào)的幅值強(qiáng)弱進(jìn)行設(shè)置，X隨 地球天然脈沖電磁場(chǎng)信號(hào)的幅值變強(qiáng)而增大。3. 根據(jù)權(quán)利要求1所述的ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法，其特征在于： 步驟（7)輸出的WVD數(shù)組進(jìn)行二維時(shí)頻譜顯示。
【專利摘要】本發(fā)明提供了一種ENPEMF信號(hào)歸一化STFT-WVD時(shí)頻分析方法，對(duì)地球天然脈沖電磁場(chǎng)信號(hào)分別作短時(shí)Fourier變換和Wigner-Ville分布，得到STFT數(shù)組和WVD數(shù)組，然后選取STFT數(shù)組中的最大值歸一化得到數(shù)組STFT_1，記錄1值位置及最小值，將其中0值的數(shù)用最小值替換，選取WVD數(shù)組中相同位置數(shù)歸一化得到臨時(shí)數(shù)組A，點(diǎn)除以STFT-1得到臨時(shí)數(shù)組B，將其中大于設(shè)置值x的位置的數(shù)與WVD數(shù)組中對(duì)應(yīng)位置的數(shù)全部置為0，輸出臨時(shí)數(shù)組B和WVD數(shù)組。本發(fā)明較好地消除交叉項(xiàng)干擾，沿襲WVD較高的時(shí)頻分辨率，克服輸入信號(hào)改變從而需要重新調(diào)節(jié)閥值與冪調(diào)節(jié)系數(shù)的缺點(diǎn)，結(jié)果理想，使用靈活。
【IPC分類】G01V3/08, G06F17/14, G01R23/16
【公開號(hào)】CN104950335
【申請(qǐng)?zhí)枴緾N201510209896
【發(fā)明人】郝國(guó)成, 陳忠昌, 趙娟, 康坊, 白雨曉, 宋佳珍
【申請(qǐng)人】中國(guó)地質(zhì)大學(xué)（武漢）
【公開日】2015年9月30日
【申請(qǐng)日】2015年4月28日