一種基于Born-Jordan時(shí)頻分布的地震信號(hào)時(shí)頻峰值濾波方法技術(shù)領(lǐng)域本發(fā)明屬于非平穩(wěn)信號(hào)時(shí)頻分析及地震信號(hào)處理領(lǐng)域，具體涉及一種基于Born-Jordan時(shí)頻分布的地震信號(hào)時(shí)頻峰值濾波方法。

背景技術(shù)：
時(shí)頻峰值濾波是一種地震信號(hào)去噪技術(shù)，它將待處理的含噪聲信號(hào)編碼為解析信號(hào)的瞬時(shí)頻率，通過對(duì)解析信號(hào)進(jìn)行時(shí)頻分析后沿頻率方向?qū)ふ曳逯祦砉烙?jì)解析信號(hào)的瞬時(shí)頻率，從而實(shí)現(xiàn)對(duì)原始信號(hào)的恢復(fù)和去噪。多項(xiàng)研究都證明了該方法在地震勘探資料處理中的實(shí)用性，在提高低信噪比地震勘探資料的信噪比以備后續(xù)分析地質(zhì)結(jié)構(gòu)等工作中起著關(guān)鍵的作用。時(shí)頻峰值濾波(Time-frequencypeakfiltering，TFPF)是M.J.Arnold、M.Roessgen和B.Boashash(1993)在時(shí)頻分析基礎(chǔ)上提出的一種有效壓制噪聲的方法。Boualem.Boashash和Mostefa.Mesbah(2004)采用此算法可恢復(fù)湮沒在加性高斯白噪聲中的人工合成多分量信號(hào)及4進(jìn)制頻移鍵控(4FSK)信號(hào)。同時(shí)可在噪聲程度低至信噪比為-9dB的背景下，清晰的恢復(fù)出新生兒腦電信號(hào)。之后，金雷(2006)將時(shí)頻峰值濾波算法用于濾除地震勘探資料的隨機(jī)噪聲，仿真試驗(yàn)表明信噪比可達(dá)－7dB，證明該方法可以有效地消減地震勘探資料中的隨機(jī)噪聲。TFPF是一維隨機(jī)噪聲壓制方法，無偏估計(jì)信號(hào)的前提條件是信號(hào)近似線性和噪聲為高斯白噪聲。在處理實(shí)際觀測(cè)信號(hào)時(shí)，為增強(qiáng)窗內(nèi)信號(hào)的線性度，Swiercz(2006)等提出一般利用偽魏格納威力分布(pseudo-Wigner-VilleDistribution，PWVD)的窗函數(shù)對(duì)信號(hào)進(jìn)行局部線性化處理，從而實(shí)現(xiàn)對(duì)信號(hào)的無失真恢復(fù)。Liu等(2014)將地震信號(hào)分解為高頻分量和低頻分量以增強(qiáng)其線性并解決了壓制噪聲和保留有效信號(hào)之間的窗長(zhǎng)矛盾問題。Tian等(2014)考慮到各道的相關(guān)性，提出了沿拋物線軌跡(Parabolic-Trace)的時(shí)頻峰值濾波算法，結(jié)果表明可以更好的適應(yīng)實(shí)際地震道變化從而增強(qiáng)信號(hào)線性。林紅波等(2015)提出了一種基于絕對(duì)級(jí)差統(tǒng)計(jì)量(ROAD)的徑向時(shí)頻峰值濾波方法，結(jié)果表明該方法可以壓制空間非平穩(wěn)地震勘探隨機(jī)噪聲且不損害有效信號(hào)，有效抑制了隨機(jī)噪聲空間非平穩(wěn)對(duì)濾波結(jié)果的影響。傳統(tǒng)的和大量改進(jìn)的時(shí)頻峰值濾波都采用加矩形窗的方法來實(shí)現(xiàn)線性化，但是當(dāng)信噪比減少低至某一給定閾值時(shí)，時(shí)頻聚集性和抑制噪聲的能力會(huì)下降，而邵煥(2013)提出的采用定向平滑偽魏格納分布(DirectionallysmoothedpseudoWigner-Villedistribution，DSPWVD)算法則可以更好的去噪。以上方法在如何增強(qiáng)時(shí)窗內(nèi)的地震信號(hào)線性度作了一定的改善，但均使用了基于PWVD對(duì)信號(hào)做時(shí)頻分析，去噪效果不是很理想，而本文方法采用Born-Jordan時(shí)頻分布對(duì)信號(hào)做時(shí)頻分析以進(jìn)一步增強(qiáng)地震信號(hào)的線性性。

技術(shù)實(shí)現(xiàn)要素：
本發(fā)明提供了一種基于Born-Jordan時(shí)頻分布的時(shí)頻峰值濾波算法，旨在增強(qiáng)時(shí)窗內(nèi)的地震信號(hào)線性度并在更大程度上抑制噪聲，從而對(duì)實(shí)際地震信號(hào)進(jìn)行更合理的去噪濾波，提高信號(hào)恢復(fù)的精確度和進(jìn)一步提高地震資料的信噪比。為解決以上技術(shù)問題，達(dá)到上述目的，本發(fā)明采用如下技術(shù)方案：一種基于Born-Jordan時(shí)頻分布的時(shí)頻峰值濾波算法包括以下步驟：步驟1：輸入待處理的含噪聲原始地震剖面；步驟2：將含噪聲的每一道地震記錄都進(jìn)行振幅值的規(guī)格化處理；步驟3：將規(guī)格化處理后的每一道地震信號(hào)都進(jìn)行編碼，作為單位振幅解析信號(hào)的瞬時(shí)頻率，得到待處理的解析信號(hào)；步驟4：對(duì)解析信號(hào)作Born-Jordan時(shí)頻分析，得到時(shí)頻分布；步驟5：在時(shí)頻分布上沿頻率方向?qū)ふ視r(shí)頻分布的峰值所對(duì)應(yīng)的解析信號(hào)的頻率，對(duì)找到的頻率進(jìn)行幅度的反規(guī)整化處理后作為有效信號(hào)的估計(jì)值；步驟6：對(duì)每一道地震信號(hào)重復(fù)步驟2-步驟5，得到去噪后的整個(gè)地震剖面的估計(jì)值；步驟7：對(duì)去噪后的地震數(shù)據(jù)進(jìn)行信噪比(SNR)分析，若信噪比不滿足需求，重復(fù)步驟2-步驟6進(jìn)行多次迭代濾波，直到滿足預(yù)先設(shè)定的閾值SNR>β,或達(dá)到預(yù)先規(guī)定的迭代次數(shù)N。其中步驟2涉及振幅值的規(guī)格化處理，具體如下：含噪聲信號(hào)為s(t)，則振幅值的規(guī)格化處理計(jì)算公式為：其中a和b兩個(gè)參數(shù)的范圍定義為0.5≥a＝max[sc(t)]>b＝min[sc(t)]≥0，在該范圍內(nèi)，選擇合適的參數(shù)使得信號(hào)不失真，其中max[]和min[]分別為取最大值函數(shù)和取最小值函數(shù)。其中步驟3涉及將地震信號(hào)編碼為單位振幅解析信號(hào)的瞬時(shí)頻率，具體理論如下：對(duì)規(guī)格化處理后的噪聲信號(hào)進(jìn)行編碼，振幅設(shè)定為單位振幅，得到瞬時(shí)頻率為含噪聲信號(hào)的待處理的解析信號(hào)：其中，為避免信號(hào)自變量t和積分變量t混淆，將sc(t)改為用sc(λ)表示，即sc(λ)為規(guī)格化處理后的含噪聲信號(hào)，μ為換算因數(shù)類同于調(diào)制指數(shù)，使用時(shí)可以設(shè)置為1。其中步驟4涉及對(duì)解析信號(hào)作時(shí)頻分析，其計(jì)算如下：對(duì)解析信號(hào)作Born-Jordan時(shí)頻分析，可以得到：其中，z的具體表達(dá)式為(2)式，為避免信號(hào)自變量t和積分變量t混淆，將z(t)改用z(u)表示，τ為時(shí)移變量，a為參數(shù)一般取為*表示共軛函數(shù)，j為虛數(shù)單位，f為頻率。其中步驟5涉及在時(shí)頻分布上尋找峰值并進(jìn)行振幅值的反規(guī)格化處理，恢復(fù)原始振幅值。最后的結(jié)果作為有效信號(hào)的估計(jì)，具體原理如下：瞬時(shí)頻率的估計(jì)是在時(shí)頻分布中沿著頻率方向?qū)ふ視r(shí)頻分布上的峰值所對(duì)應(yīng)的頻率，即：其中，代表沿頻率方向?qū)ふ視r(shí)頻分布峰值所對(duì)應(yīng)的頻率。再對(duì)得到的瞬時(shí)頻率作反規(guī)格化處理，便可得到濾波后理想的數(shù)據(jù)估計(jì)值。反規(guī)格化計(jì)算公式為：式中fz(t)為由時(shí)頻分布得到的峰值所對(duì)應(yīng)的頻率，而是所需要的對(duì)有效信號(hào)的估計(jì)值。與現(xiàn)有技術(shù)相比，本發(fā)明具有如下優(yōu)勢(shì)：本發(fā)明將傳統(tǒng)的時(shí)頻峰值濾波算法中采用偽魏格納時(shí)頻分布(PWVD)改為采用Born-Jordan時(shí)頻分布。PWVD中，對(duì)信號(hào)在時(shí)域加矩形窗，由此在頻域帶來振鈴效應(yīng)而影響了濾波去噪的效果，而本發(fā)明消除了頻域的振鈴現(xiàn)象。同時(shí)由于Born-Jordan時(shí)頻分布的核函數(shù)為sinc函數(shù)，對(duì)旁瓣有抑制作用，且為理想的低通濾波函數(shù)，故可以對(duì)噪聲起到更好的壓制效果，去噪的同時(shí)保留了有效信號(hào)，在實(shí)際應(yīng)用中效果更好。附圖說明圖1為方法流程圖，tracenumber為道數(shù)，n濾波次數(shù)，N為設(shè)定的濾波次數(shù)；圖2為含噪聲的原始地震剖面；圖3為地震數(shù)據(jù)中某一道在幅值規(guī)整化前的波形；圖4為地震數(shù)據(jù)中某一道在幅值規(guī)整化后的波形；圖5為地震數(shù)據(jù)中某一道經(jīng)時(shí)頻峰值濾波后未做幅值規(guī)整化處理的結(jié)果；圖6為地震數(shù)據(jù)中某一道經(jīng)時(shí)頻峰值濾波后做幅值反規(guī)整化處理后的結(jié)果；圖7為經(jīng)時(shí)頻逢值濾波算法去噪后的地震數(shù)據(jù)。具體實(shí)施方式下面結(jié)合附圖對(duì)本發(fā)明進(jìn)行具體說明：步驟1：輸入待處理的含噪聲的地震剖面(如圖2所示)，循環(huán)讀取每一道地震數(shù)據(jù)。步驟2：將讀入的含噪聲的地震剖面的每一道地震道數(shù)據(jù)都進(jìn)行以下幅值規(guī)格化處理，其計(jì)算公式為式中，s(t)即為讀取的含噪聲的每一道地震數(shù)據(jù)，計(jì)算時(shí)設(shè)定a＝0.5，b＝0，選取第20道為例進(jìn)行幅值規(guī)格化處理，其結(jié)果如圖3和圖4所示(幅值規(guī)整化前如圖3所示，幅值規(guī)整化后如圖4所示)。步驟3：將該含噪聲的地震數(shù)據(jù)的每一道都進(jìn)行編碼，作為單位振幅解析信號(hào)的瞬時(shí)頻率，得到待處理的解析信號(hào)，計(jì)算公式為將步驟2處理后的第20道地震數(shù)據(jù)進(jìn)行編碼，其中換算因數(shù)μ設(shè)定為1。步驟4：對(duì)每一道編碼后的解析信號(hào)z(t)都作Born-Jordan時(shí)頻分析，得到地震數(shù)據(jù)每一道的時(shí)頻分布，其計(jì)算公式為式中，Wz(t,f)為每一道地震數(shù)據(jù)所得到的時(shí)頻分布，參數(shù)a設(shè)定為1/2，τ為時(shí)移變量，*表示步驟3中的解析信號(hào)的共軛函數(shù)，j為虛數(shù)單位，f為頻率。對(duì)第20道地震數(shù)據(jù)做時(shí)頻分析，得到其時(shí)頻分布。步驟5：在時(shí)頻分布Wz(t,f)上沿著頻率方向?qū)ふ曳逯邓鶎?duì)應(yīng)的頻率作為每一道數(shù)據(jù)的估計(jì)值，其尋找原理為再對(duì)找到的峰值頻率fz(t)做幅值反規(guī)格化處理，得到去噪后的有效信號(hào)，計(jì)算公式為選取之前的第20道地震數(shù)據(jù)的時(shí)頻分布為例，找到該道的估計(jì)值并對(duì)其做反規(guī)整化處理。步驟6：對(duì)該道濾波結(jié)果進(jìn)行信噪比分析，若信噪比或迭代次數(shù)不滿足需求，重復(fù)步驟2-步驟5。根據(jù)多次實(shí)驗(yàn)情況，設(shè)定對(duì)該道的濾波次數(shù)為2次，經(jīng)過2次時(shí)頻峰值濾波后，得到最終去噪后單道地震記錄，其結(jié)果如圖5和圖6所示(未做幅值反規(guī)整化處理波形如圖5所示，經(jīng)幅值反規(guī)整化處理后波形如圖6所示)。步驟7：循環(huán)讀取所有地震道數(shù)據(jù)，重復(fù)步驟2-6，直到tracenum>＝TotalNum(總道數(shù))。步驟8：輸出最終經(jīng)時(shí)頻峰值濾波后的地震剖面圖(如圖7所示)。