本發(fā)明涉及一種ENPEMF數(shù)據(jù)的格林函數(shù)-SPWVD時(shí)頻分析方法，屬于非平穩(wěn)數(shù)據(jù)的時(shí)頻分析處理領(lǐng)域。

背景技術(shù)：

地震給人類的生活帶來(lái)了巨大的災(zāi)難，據(jù)統(tǒng)計(jì)，全球的自然災(zāi)害之中，地震造成的死亡人數(shù)占全部自然災(zāi)害死亡人數(shù)的54％，堪稱自然災(zāi)害之最。如何預(yù)測(cè)地震一直以來(lái)都是一個(gè)熱門敏感的話題。然而，因地震預(yù)測(cè)有著地球內(nèi)部的“不可入性”、大地震的“非頻發(fā)性”、“地震物理過(guò)程的復(fù)雜性”三大困難，地震預(yù)測(cè)成為了公認(rèn)的世界性難題，對(duì)于地震前兆預(yù)測(cè)對(duì)于人類生命安全和社會(huì)財(cái)產(chǎn)安全的具有很大的意義。

在現(xiàn)有的孕震信息研究過(guò)程中，由現(xiàn)有的信號(hào)采集裝置采集到的地球天然脈沖電磁場(chǎng)(ENPEMF)信號(hào)由于其本身存在時(shí)頻分布不準(zhǔn)確的問(wèn)題，需要經(jīng)過(guò)繁雜的處理才能輸出為穩(wěn)定的二維時(shí)頻譜，再做進(jìn)一步分析。

在現(xiàn)有的孕震信息研究過(guò)程中，STFT與WVD時(shí)頻分析方法是常用的用來(lái)對(duì)采集到的大量地球天然脈沖電磁場(chǎng)信號(hào)進(jìn)行分析的方法。傳統(tǒng)的STFT或WVD時(shí)頻分析方法會(huì)用到以下公式對(duì)地球天然脈沖電磁場(chǎng)信號(hào)進(jìn)行處理：

傳統(tǒng)的STFT或WVD時(shí)頻分析方法需要設(shè)置閥值c以及冪調(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相乘以消除交叉項(xiàng)；其中a和b式為冪調(diào)節(jié)系數(shù)，作用是增強(qiáng)兩變換數(shù)值較大部分而消弱有交叉項(xiàng)部分。式(a)所示方法靈活性差，輸入信號(hào)幅值或能量大小直接影響c值的選擇，目前沒(méi)有人能提出一種行之有效的自適應(yīng)閥值選擇方法，并且在實(shí)際信號(hào)中有用分量幅值、能量大小往往各不相同甚至差別較大，因此采用設(shè)置閥值消交叉項(xiàng)容易也將信息項(xiàng)消除；式(b)所示方法有些許改進(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í)頻圖和譜圖。

技術(shù)實(shí)現(xiàn)要素：

為了解決現(xiàn)有技術(shù)的不足，本發(fā)明提供了一種ENPEMF數(shù)據(jù)的格林函數(shù)-SPWVD時(shí)頻分析方法，可以處理地球天然脈沖電磁場(chǎng)信號(hào)數(shù)據(jù)，獲取較為明顯的震前時(shí)頻能量譜的分布異常，得出基于地球天然脈沖電磁場(chǎng)分析孕震信息的有效性，該方法的結(jié)構(gòu)新穎，效果良好，滿足科學(xué)分析研究的需要，對(duì)地震前兆研究以及非平穩(wěn)信號(hào)的研究有積極的意義。

本發(fā)明為解決其技術(shù)問(wèn)題所采用的技術(shù)方案是：提供了一種ENPEMF數(shù)據(jù)的格林函數(shù)-SPWVD時(shí)頻分析方法，包括以下步驟：

(1)對(duì)原始ENPEMF數(shù)據(jù)進(jìn)行平滑壓縮預(yù)處理；

(2)針對(duì)步驟(1)中平滑壓縮后的ENPEMF數(shù)據(jù)，采用基于泊松方程的格林函數(shù)進(jìn)行處理，得到格林函數(shù)化的ENPEMF數(shù)據(jù)；

(3)對(duì)經(jīng)步驟(2)處理得到的格林函數(shù)化的ENPEMF數(shù)據(jù)，采用平滑偽WVD算法進(jìn)行掃頻處理；

(4)對(duì)經(jīng)步驟(3)處理得到的數(shù)據(jù)，進(jìn)行時(shí)頻分析匯總，得出原數(shù)據(jù)的時(shí)間-頻率聯(lián)合分布圖。

步驟(1)所述平滑壓縮預(yù)處理采用平方平均算法。

步驟(2)所述針對(duì)步驟(1)中平滑壓縮后的ENPEMF數(shù)據(jù)，采用基于泊松方程的格林函數(shù)進(jìn)行處理，得到格林函數(shù)化的ENPEMF數(shù)據(jù)，具體包括將基于泊松方程的格林函數(shù)正規(guī)化，再將平滑壓縮后的ENPEMF數(shù)據(jù)代入正規(guī)化后的格林函數(shù)得到格林函數(shù)化的ENPEMF數(shù)據(jù)。

步驟(3)所述對(duì)經(jīng)步驟(2)處理得到的格林函數(shù)化的ENPEMF數(shù)據(jù)，采用平滑偽WVD算法進(jìn)行掃頻處理，具體過(guò)程為：設(shè)格林函數(shù)化的ENPEMF數(shù)據(jù)用以下疊加信號(hào)形式表示：

其中，φ₁(t)和φ₂(t)為相位；則利用以下公式對(duì)格林函數(shù)化的ENPEMF數(shù)據(jù)進(jìn)行時(shí)頻分析變換：

其中，t和f分別為格林函數(shù)化的ENPEMF數(shù)據(jù)的時(shí)域變量和頻域變量，是格林函數(shù)化的ENPEMF數(shù)據(jù)在時(shí)頻和頻域的各自成分特點(diǎn)描述；u和τ分別表示格林函數(shù)化的ENPEMF數(shù)據(jù)的時(shí)域積分變量和頻域積分變量，g(u)和h(τ)為設(shè)置的實(shí)偶窗函數(shù)。

實(shí)偶窗函數(shù)g(u)、h(τ)均采用高斯窗。

本發(fā)明基于其技術(shù)方案所具有的有益效果在于：

(1)本發(fā)明充分利用了格林函數(shù)的疊加思想：格林函數(shù)是表示一種特定的場(chǎng)和產(chǎn)生這種場(chǎng)的源之間的數(shù)學(xué)物理方程，即當(dāng)一個(gè)源被分解成很多點(diǎn)源的疊加時(shí)，如果能知道分解后每個(gè)點(diǎn)源產(chǎn)生的場(chǎng)，利用疊加原理就可以求出同樣邊界條件下任意源的場(chǎng)；本發(fā)明基于格林函數(shù)的疊加思想，將原信號(hào)首先經(jīng)過(guò)格林函數(shù)分解，能夠加強(qiáng)頻率的分布效果；再對(duì)經(jīng)格林函數(shù)分解后的信號(hào)進(jìn)行SPWVD算法的掃頻分解，最后得到原始信號(hào)的時(shí)間-頻率聯(lián)合分布；

(2)本發(fā)明可以較好的降噪和降低WVD掃頻時(shí)產(chǎn)生的交叉項(xiàng)，并且結(jié)構(gòu)新穎，步驟簡(jiǎn)單，在針對(duì)大量繁雜地球天然脈沖電磁場(chǎng)(ENPEMF)數(shù)據(jù)分析時(shí)能夠提高處理效率，具有良好的應(yīng)用前景；

(3)本發(fā)明針對(duì)特殊的非平穩(wěn)ENPEMF數(shù)據(jù)采用平方平均算法進(jìn)行平滑壓縮處理，與算術(shù)平均、幾何平均、調(diào)和平均等平滑壓縮處理方法相比，平方平均算法壓縮后的數(shù)據(jù)包絡(luò)與原數(shù)據(jù)包絡(luò)一致性最高，效果最優(yōu)；

(4)經(jīng)大量仿真實(shí)驗(yàn)論證，本發(fā)明與傳統(tǒng)的STFT方法、小波變換(Wavelet)方法相比，能夠得到效果最優(yōu)的時(shí)頻分布圖。

附圖說(shuō)明

圖1是本發(fā)明的一種ENPEMF數(shù)據(jù)的格林函數(shù)-SPWVD時(shí)頻分析方法的流程示意圖。

圖2是本發(fā)明實(shí)施例信號(hào)G1的FFT頻率分布圖，其中圖2(1)是為時(shí)域圖，圖2(2)為頻域圖。

圖3是本發(fā)明實(shí)施例信號(hào)G1的STFT時(shí)頻圖。

圖4是本發(fā)明實(shí)施例信號(hào)G1的Wavelet時(shí)頻圖。

圖5是本發(fā)明實(shí)施例信號(hào)G1的WVD時(shí)頻圖。

圖6是本發(fā)明實(shí)施例信號(hào)x1的直接SPWVD時(shí)頻圖。

圖7是本發(fā)明實(shí)施例信號(hào)G1的SPWVD時(shí)頻圖。

圖8是本發(fā)明實(shí)施例信號(hào)G2的FFT頻域圖。

圖9是本發(fā)明實(shí)施例信號(hào)G2的STFT時(shí)頻圖。

圖10是本發(fā)明實(shí)施例信號(hào)G2的Wavelet時(shí)頻圖。

圖11是本發(fā)明實(shí)施例信號(hào)G2的WVD時(shí)頻圖。

圖12是本發(fā)明實(shí)施例信號(hào)x2直接WVD時(shí)頻圖。

圖13是本發(fā)明實(shí)施例信號(hào)x2直接SPWVD時(shí)頻圖。

圖14是本發(fā)明實(shí)施例信號(hào)G2的SPWVD時(shí)頻圖。

圖15是本發(fā)明實(shí)施例420CN2AH的FFT時(shí)頻圖。

圖16是本發(fā)明實(shí)施例420CN2AH的格林函數(shù)化后的SPWVD時(shí)頻圖。

圖17是原始ENPEMF數(shù)據(jù)1的時(shí)間-脈沖數(shù)示意圖。

圖18是原始ENPEMF數(shù)據(jù)1的不同壓縮方法處理效果比較示意圖。

圖19是原始ENPEMF數(shù)據(jù)2的日期-脈沖數(shù)示意圖。

圖20是原始ENPEMF數(shù)據(jù)2的不同壓縮方法處理效果比較示意圖。

具體實(shí)施方式

下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明作進(jìn)一步說(shuō)明。

本發(fā)明提供了一種ENPEMF數(shù)據(jù)的格林函數(shù)-SPWVD時(shí)頻分析方法，參照?qǐng)D1，包括以下步驟：

(1)對(duì)原始ENPEMF數(shù)據(jù)進(jìn)行平滑壓縮預(yù)處理；所述平滑壓縮預(yù)處理采用平方平均算法，效果最優(yōu)。

如圖17至圖20所示，給出了兩組原始ENPEMF數(shù)據(jù)用不同的壓縮方法處理的效果比較，可以很明顯對(duì)比出平方平均方法處理得到的數(shù)據(jù)最接近原始ENPEMF數(shù)據(jù)的波動(dòng)韻律，再經(jīng)過(guò)后續(xù)步驟進(jìn)行分析，可以最大程度保證數(shù)據(jù)有效性。

(2)針對(duì)步驟(1)中平滑壓縮后的ENPEMF數(shù)據(jù)(即多頻信號(hào)s)，采用基于泊松方程的格林函數(shù)進(jìn)行處理，得到格林函數(shù)化的ENPEMF數(shù)據(jù)；

基于泊松方程的格林函數(shù)為：

為得到正確結(jié)果，將公式1進(jìn)行正規(guī)化得到以下公式：

公式2中，m為引入的一個(gè)無(wú)窮小的質(zhì)量參數(shù)，r＝|x|，此時(shí)再令m→0，得：

公式3即為正規(guī)化后的格林函數(shù)。將多頻信號(hào)s作為變量代入公式3，即得到格林函數(shù)化的ENPEMF數(shù)據(jù)；

(3)對(duì)經(jīng)步驟(2)處理得到的格林函數(shù)化的ENPEMF數(shù)據(jù)，采用平滑偽WVD算法進(jìn)行掃頻處理；

通常，一個(gè)疊加信號(hào)的WVD分布包括以下四項(xiàng)：

前兩項(xiàng)為信號(hào)自身項(xiàng)(Signalterms)，后兩項(xiàng)就是交叉干擾項(xiàng)(Crossterms)，即為需要抑制的由于算法變換產(chǎn)生的噪聲分量。

假設(shè)格林函數(shù)化的ENPEMF數(shù)據(jù)用以下疊加信號(hào)形式表示：

其中，φ₁(t)和φ₂(t)為相位，s(t)確定后φ₁(t)、φ₂(t)以及j即可確定。

將信號(hào)(公式5)代入公式4，得：

前兩項(xiàng)自身項(xiàng)不是信號(hào)瞬時(shí)頻率支撐區(qū)內(nèi)的任何分量，來(lái)自自身項(xiàng)的疊加，表現(xiàn)出交叉項(xiàng)的特點(diǎn)。

以式(公式6)的第一項(xiàng)來(lái)分析，將φ₁(t)按照Taylor級(jí)欺展開得到：

由公式7可知，由屬于同一個(gè)信號(hào)分量的不同部分之間相互作用產(chǎn)生的內(nèi)部交叉項(xiàng)(自交叉項(xiàng))也是干擾噪聲源之一，有別于由不同的信號(hào)分量之間相互作用造成的外部交叉項(xiàng)。這兩項(xiàng)都是需要采用合適的方法來(lái)去除或抑制。

在存在多個(gè)信號(hào)分量的情況下，原始信號(hào)的WVD時(shí)頻圖之間產(chǎn)生的交叉項(xiàng)將十分顯著。并且，當(dāng)信號(hào)分量越多時(shí)，其交叉項(xiàng)分布也就越復(fù)雜。為了抑制交叉項(xiàng)，出現(xiàn)一種WVD的改良算法，即平滑偽WVD，簡(jiǎn)稱SPWVD。

在WVD(公式4)的基礎(chǔ)上增加兩個(gè)窗函數(shù)，改良為SPWVD，其定義如下：

公式8中，其中，t和f分別為格林函數(shù)化的ENPEMF數(shù)據(jù)的時(shí)域變量和頻域變量，是格林函數(shù)化的ENPEMF數(shù)據(jù)在時(shí)頻和頻域的各自成分特點(diǎn)描述；u和τ分別表示格林函數(shù)化的ENPEMF數(shù)據(jù)的時(shí)域積分變量和頻域積分變量，格林函數(shù)化的ENPEMF數(shù)據(jù)確定后u和τ即可確定。

g(u)和h(τ)為設(shè)置的實(shí)偶窗函數(shù)，且h(0)＝g(0)＝1，SPWVD的作用實(shí)質(zhì)是相當(dāng)于WVD在時(shí)域和頻域方向同時(shí)進(jìn)行平滑濾波處理。通過(guò)調(diào)整上述兩個(gè)窗g(u)、h(τ)的寬度，可以有效的抑制交叉項(xiàng)。需要注意的是，時(shí)頻圖的分辨率與窗的寬度大小成反比。因此，根據(jù)實(shí)際需求來(lái)選取合適的窗將變得尤為重要。g(u)、h(τ)的形狀可選擇矩形窗、高斯窗、布萊克曼、漢明窗、漢寧窗等，這些形狀的窗各有優(yōu)缺點(diǎn)，針對(duì)ENPEMF數(shù)據(jù)的非周性、非平穩(wěn)性，由于高斯窗具有最好的時(shí)寬-帶寬乘積(等于1/2)，效果相對(duì)較好，這里的兩個(gè)窗函數(shù)都選擇高斯窗。

由上述步驟(2)中得到格林函數(shù)化的ENPEMF數(shù)據(jù)G(x)，然后再對(duì)G(x)信號(hào)按SPWVD方法(公式8)進(jìn)行時(shí)頻分析變換，即完成平滑偽WVD算法進(jìn)行掃頻處理。

(4)對(duì)經(jīng)步驟(3)處理得到的數(shù)據(jù)，進(jìn)行時(shí)頻分析匯總，得出原數(shù)據(jù)的時(shí)間-頻率聯(lián)合分布圖。具體根據(jù)數(shù)據(jù)分析的需要，采用matlab中的畫圖輸出函數(shù)plot()，可輸出二維時(shí)頻圖，也可輸出三維時(shí)頻圖；可繪制單天的數(shù)據(jù)時(shí)頻分析，也可繪制多天的數(shù)據(jù)時(shí)頻分析進(jìn)行對(duì)比研究；可以單幅圖像輸出，也可多幅圖像對(duì)比輸出等等。即根據(jù)分析要求和目的采用有針對(duì)性的時(shí)頻分析匯總。

對(duì)本發(fā)明的效果進(jìn)行仿真驗(yàn)證。

仿真實(shí)驗(yàn)1：假定有一個(gè)4頻分量的信號(hào)x1：

x1＝sin(2*pi*300*t1)+sin(2*pi*150*t1)+sin(2*pi*250*t1)+sin(2*pi*90*t1)。

圖2為原信號(hào)經(jīng)過(guò)格林函數(shù)化后的時(shí)域圖和頻率分布圖，原信號(hào)x1經(jīng)過(guò)格林函數(shù)化后表示為G1。

信號(hào)G1為格林函數(shù)化后的4頻疊加信號(hào)，從圖2至圖5中可以看出，快速傅里葉(FFT)可清晰反映G1信號(hào)的頻率分布，但無(wú)法顯示時(shí)間與頻率的聯(lián)合分布；短時(shí)傅里葉(STFT)時(shí)頻圖模糊，分辨率低；小波(Wavelet)時(shí)頻圖相對(duì)清晰正確，但是頻率的聚集度不好；Wigner-Ville分布(WVD)時(shí)頻圖在四個(gè)頻率之間存在5條交叉項(xiàng)干擾，干擾頻率大致分別為120、170、200、230、270。圖6為信號(hào)x1直接采用平滑偽Wigner-Ville分布(SPWVD)分解，沒(méi)有結(jié)合格林函數(shù)的時(shí)頻圖，在時(shí)間軸兩端會(huì)出現(xiàn)一些噪聲干擾。而圖7為采用格林函數(shù)后的時(shí)頻分析，效果較好。

仿真實(shí)驗(yàn)2：對(duì)于信號(hào)x2＝s1+s2：

s1＝(1+0.3*cos(2*t)).*exp(-t./15).*cos(2*pi*(2.4*t+0.5*t.^1.2+0.3*sin(t)))；

s2＝cos(2*pi*(5.3*t+0.2*t.^1.3))；

t＝linspace(0,10,2048)；x1＝(1+0.2*cos(t)).*cos(2*pi.*(2*t+0.3*cos(t)))。

圖8為信號(hào)x2經(jīng)過(guò)格林函數(shù)化后得到的G2的時(shí)域圖和快速傅里葉變換后的頻域圖。圖9為信號(hào)G2的短時(shí)傅里葉變換(STFT)后的時(shí)頻分布圖，圖10為信號(hào)G2的小波變換(Wavelet)時(shí)頻圖，圖11為信號(hào)G2的Wigner-Ville分布(WVD)時(shí)頻圖，圖12為原始信號(hào)x2的直接的Wigner-Ville分布(WVD)時(shí)頻圖，圖13原始信號(hào)x2的直接平滑偽Wigner-Ville分布(SPWVD)時(shí)頻圖，圖14為本發(fā)明的原始信號(hào)x2經(jīng)過(guò)格林函數(shù)化后G2的SPWVD時(shí)頻圖。

從圖8至圖14可以看出，F(xiàn)FT只能單獨(dú)觀察頻率分布，無(wú)法了解時(shí)間與頻率的聯(lián)合分布，經(jīng)過(guò)格林函數(shù)化后的信號(hào)G2在多種時(shí)頻分析方法的表現(xiàn)都不盡如意：STFT時(shí)頻圖十分模糊，難以分辨，Wavelet時(shí)頻圖較好，f＝30Hz處稍顯模糊，頻率分布的聚焦度較差，WVD時(shí)頻圖較好，但是交叉項(xiàng)干擾嚴(yán)重，出現(xiàn)了原信號(hào)沒(méi)有的新增頻率，影響了對(duì)原始信號(hào)的頻率分析。如果原始信號(hào)不經(jīng)過(guò)格林函數(shù)化，則從圖12和圖13可以看出，頻率分布效果很差。只有采用本發(fā)明的圖14，原始信號(hào)先經(jīng)過(guò)格林函數(shù)化后，再結(jié)合SPWVD時(shí)頻分析方法，最后得到的時(shí)頻聯(lián)合分布的效果最優(yōu)。

將本方法應(yīng)用于地球天然脈沖電磁場(chǎng)(ENPEMF)信號(hào)的時(shí)頻聯(lián)合分析中，可以較好的得到該信號(hào)的時(shí)間-頻率的聯(lián)合分布情況，有助于進(jìn)一步了解地震發(fā)生前ENPEMF信號(hào)的時(shí)頻分布變化，為地震的前兆信號(hào)特性研究提供依據(jù)和幫助。圖15為2013年4月20日第二通道AH數(shù)據(jù)類型(420CN2AH)的信號(hào)時(shí)頻分布圖和快速傅里葉的頻率分布圖。圖16為該ENPEMF數(shù)據(jù)經(jīng)過(guò)格林函數(shù)化后的SPWVD時(shí)頻圖。