本發(fā)明屬于地震勘探中數(shù)據(jù)處理技術(shù)領(lǐng)域，特別涉及一種地震采集數(shù)據(jù)中諧波噪聲壓制的方法。

背景技術(shù)：

隨著全球經(jīng)濟的快速發(fā)展以及環(huán)保意識的提高，傳統(tǒng)的炸藥震源逐漸被可控震源替代?；瑒訏呙枋强煽卣鹪词┕ぶ谐Ｓ玫牟杉绞?，提高了采集效率。但滑動掃描存在采集地震記錄受到較強諧波噪聲干擾的問題，有效信號往往掩藏在諧波噪聲中，對數(shù)據(jù)資料的質(zhì)量造成了嚴重影響。目前壓制諧波噪聲有如下方法：一是純相移濾波法，但該方法處理相關(guān)前地震記錄，數(shù)據(jù)量較大、計算效率低；二是分頻異常振幅壓制法，但該方法對有效信號會造成損傷，保真性不夠；三是地面力信號法，但需要采集地面力信號，采集成本較高，且必須從未受影響的最后一炮數(shù)據(jù)處理，無法進行單炮集處理。

隨著信號稀疏性理論的發(fā)展，Starck等人提出了形態(tài)成分分析的混合信號分解方法。形態(tài)成分分析是指，根據(jù)復(fù)雜信號的組成成分波形形態(tài)特征，將兩種具有不同原子特征的變換字典構(gòu)成超完備字典，實現(xiàn)對復(fù)雜信號更稀疏的表示方式以及更有效的信息識別能力，實現(xiàn)兩種成分的分離。字典通常是根據(jù)經(jīng)驗從已知的數(shù)學(xué)變換中選擇或者構(gòu)造，而選出的字典是否充分滿足形態(tài)成分分析的假設(shè)，是形態(tài)成分分析方法能否成功的關(guān)鍵。使用形態(tài)成分分析的方法對信號與諧波噪聲實現(xiàn)分離來達到壓制諧波噪聲的目的，需要根據(jù)有效信號與諧波噪聲的波形形態(tài)特征構(gòu)造合適的稀疏變換。已有文獻研究選擇使用連續(xù)小波變換來稀疏表示有效信號比較合適。但諧波噪聲較為復(fù)雜，需要構(gòu)造較為合適匹配的變換來稀疏表示，并選擇確定變換的具體參數(shù)，目前相關(guān)研究較少。

針對上述問題，本發(fā)明利用形態(tài)成分分析理論，給出了針對諧波噪聲的波形形態(tài)特征，構(gòu)造合適的稀疏表示變換即Chirplet變換，與連續(xù)小波變換構(gòu)成超完備字典，并根據(jù)數(shù)據(jù)特征選擇確定Chirplet變換參數(shù)，以及使用坐標分塊松弛算法方法實現(xiàn)信噪分離達到壓制諧波噪聲的解決方案。

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

本發(fā)明的目的在于提供了一種基于波形形態(tài)特征稀疏化建模的諧波噪聲壓制方法。針對諧波噪聲的波形形態(tài)特征，構(gòu)造合適的稀疏表示變換即Chirplet變換，與連續(xù)小波變換聯(lián)合構(gòu)成超完備字典，并根據(jù)數(shù)據(jù)特征選擇確定Chirplet變換參數(shù)，以及使用坐標分塊松弛算法方法實現(xiàn)信噪分離以達到有效壓制可控震源滑動掃描記錄的地震數(shù)據(jù)中的諧波噪聲的目的。

為了實現(xiàn)上述目的，本發(fā)明采用如下技術(shù)方案：

基于波形形態(tài)特征稀疏化建模的諧波噪聲壓制方法，包括以下步驟：

1)、根據(jù)可控震源滑動掃描采集地震記錄中諧波噪聲波形形態(tài)特征構(gòu)造Chirplet變換，并與連續(xù)小波變換聯(lián)合構(gòu)成超完備字典；

2)、Chirplet正、反變換的快速實現(xiàn)；

3)、基于相關(guān)后數(shù)據(jù)時頻分布特征確定Chirplet變換參數(shù)；

4)、根據(jù)參考掃描信號的起始頻率確定諧波噪聲的濾波截止頻率，實現(xiàn)有效信號與諧波噪聲保真分離。

進一步的，根據(jù)可控震源滑動掃描采集地震記錄中諧波噪聲波形形態(tài)特征構(gòu)造Chirplet變換，并與連續(xù)小波變換構(gòu)成超完備字典，包括：

形態(tài)成分分析的對象是含有兩種具有不同形態(tài)特征的成分：

式中：表示待分析信號；表示信號中的兩種成分，具有不同的形態(tài)特征；形態(tài)成分分析的目標是分別提取出兩種成分；假設(shè)和可以分別由字典A₁和A₂有效的稀疏表示，但是用A₂稀疏表示和用A₁稀疏表示時稀疏性差；

在地震數(shù)據(jù)處理中，一般選擇連續(xù)小波變換作為稀疏表示有效信號成分的變換字典，其中連續(xù)小波變換為：

式中WT_x(A,τ)為變換系數(shù)，A表示尺度因子，x(t)表示待分析信號，ψ(t)表示Morlet母小波；其中t為時間，τ為平移量，*表示共軛；

連續(xù)小波變換的反變換為：

式中常數(shù)C_Ψ＜∞為其容許條件；

本發(fā)明中，構(gòu)成合適的Chirplet變換作為稀疏表示諧波噪聲的字典，其中Chirplet正變換的定義為：

CT_x(t,f,a,c,d)＝∫x(τ)h^*(τ-t,f,a,c,d)dτ，

式中CT_x(t,f,a,c,d)表示待分析信號的Chirplet變換值，x(t)表示待分析信號，a為伸縮參數(shù)，c為時間上的線性調(diào)頻率，f為頻率，d為頻率上的線性調(diào)頻率；τ為時間平移量，*表示共軛；h(t,f,a,c,d)為Chirplet變換的核函數(shù)，其表達式如下：

式中g(shù)(t)為高斯函數(shù)；v為平移量；

忽略頻率上的線性調(diào)頻操作，核函數(shù)退化為：

根據(jù)形態(tài)成分分析理論，用上面選定的字典A₁即連續(xù)小波變換和A₂即Chirplet變換，構(gòu)成超完備字典，稀疏表示信號計算稀疏表示系數(shù)：

式中：x₁為重構(gòu)系數(shù)中與A₁對應(yīng)的部分；x₂為重構(gòu)系數(shù)中與A₂對應(yīng)的部分。為拉格朗日乘子；該優(yōu)化問題可以通過坐標分塊松弛算法求解。坐標分塊松弛算法的基本思想是交替迭代的計算x₁和x₂。其主要內(nèi)容步驟為：

初始化：初始迭代步數(shù)k＝0，初始解

迭代：每步迭代k增加1，并計算：

式中，T_λ為硬閾值函數(shù)；與A₁構(gòu)成一對正反變換，與A₂構(gòu)成一對正反變換；

終止條件：當小于預(yù)設(shè)的值時，即繼續(xù)迭代對結(jié)果的影響足夠小時，迭代終止。

輸出：

進一步的，步驟2)即Chirplet正、反變換的快速實現(xiàn)中，具體包括：

步驟101：從原始數(shù)據(jù)中讀出單道時域信號，記為x；

步驟102：計算得到信號x的Chirplet正變換結(jié)果，記為y，其實現(xiàn)方法類似于短時傅里葉變換；

步驟103：計算y的Chirplet反變換結(jié)果，記為

由于步驟103中Chirplet反變換得到的重構(gòu)結(jié)果與原始信號相比會在幅度上存在一定的差別，因此接下來使用一個模型信號計算修正系數(shù)來進行幅度修正；

步驟201：構(gòu)造一個線性升頻信號模型，記為s，其長度與信號x相同；

步驟202：計算得到信號s的Chirplet正變換結(jié)果，記為s₀，其實現(xiàn)方法類似于短時傅里葉變換；

步驟203：計算得到s₀的Chirplet反變換結(jié)果，記為

步驟204：根據(jù)模型反變換結(jié)果計算修正系數(shù)系數(shù)，即其中||.||₂表示二范數(shù)；

步驟104：將修正系數(shù)r作用到單道時域信號x的Chirplet反變換得到修正的Chirplet反變換結(jié)果，記為則

進一步的，步驟3)即相關(guān)后數(shù)據(jù)的時頻分布特征計算Chirplet變換的參數(shù)包括：

地震資料數(shù)據(jù)是陣列信號，利用形態(tài)成分分析實現(xiàn)有效信號與諧波噪聲分離的方法基于單道數(shù)據(jù)處理，而由于各道數(shù)據(jù)的諧波噪聲強弱不同，因此需要根據(jù)數(shù)據(jù)計算確定Chirplet變換的參數(shù)才能保證對諧波噪聲稀疏，更好的實現(xiàn)信噪分離。其中具體實現(xiàn)步驟如下：

步驟301：將單道時域信號x轉(zhuǎn)換到頻域，記為F；

步驟302：確定頻率域的高、低頻分界，記為κ；

步驟303：計算高頻(頻率分界以上)能量占總頻帶能量的比值，即振幅譜比值，記為b；

步驟304：根據(jù)振幅譜比值計算Chirplet變換的伸縮參數(shù)，記為a，線性調(diào)頻參數(shù)，記為c以及閾值權(quán)系數(shù)參數(shù)，記為λ：

a＝k₁×b，λ＝k₂-b，c＝-k₃×b×a²，

其中，b為振幅譜比值；k₁,k₂,k₃分別為參數(shù)a、c與λ的修正系數(shù)；高低頻分界以及修正系數(shù)需在分析地震資料數(shù)據(jù)的特點之后合理給出。

進一步的，步驟4)即根據(jù)參考掃描信號的起始頻率確諧波噪聲的濾波截止頻率，實現(xiàn)有效信號與諧波噪聲保真分離具體包括：

步驟401：根據(jù)參考掃描信號的起始頻率確定諧波噪聲的起始頻率f_b；

步驟402：諧波噪聲做截止頻率為f_b的高通濾波；

步驟403：原始數(shù)據(jù)減去濾波后的諧波噪聲得到有效信號，實現(xiàn)有效信號與諧波噪聲的保真分離。

相對于現(xiàn)有技術(shù)，本發(fā)明具有以下有益效果：

1)本發(fā)明方法根據(jù)諧波噪聲的時頻分布特征構(gòu)造合適的變換，保證了變換對諧波噪聲的稀疏性；

2)Chirplet正、反變換的快速實現(xiàn)方式保證了變換的計算效率，使其可用于海量實際資料處理；

3)Chirplet反變換中使用幅度修正系數(shù)使得反變換重構(gòu)的準確性得到保證；

4)根據(jù)數(shù)據(jù)驅(qū)動自動確定Chirplet變換的參數(shù)，具有很強的適應(yīng)性，且單道計算，易于并行實現(xiàn)；

5)本發(fā)明對有效信號具有高保真性，能夠有效的保護有效波的高、低頻成分。

基于本發(fā)明的時頻域稀疏優(yōu)化去噪方法解決了可控震源滑動掃描采集地震數(shù)據(jù)中諧波噪聲干擾的問題，達到了提高了地震資料信噪比的目的。

附圖說明

圖1為諧波噪聲的時域波形圖；

圖2為Chirplet變換原子示意圖；

圖3為Chirplet正、反變換的快速實現(xiàn)流程圖；

圖4為根據(jù)相關(guān)后數(shù)據(jù)時頻分布特征確定Chirplet變換參數(shù)流程圖；

圖5為基于參考掃描信號起始頻率的有效信號與諧波噪聲保真分離流程圖；

圖6A為未受諧波干擾模擬數(shù)據(jù)；圖6B為受諧波干擾模擬數(shù)據(jù)；

圖7A本發(fā)明方法從圖6B中分離的有效信號；圖7B本發(fā)明方法從圖6B中分離的諧波噪聲；

圖8A為實際滑動掃描地震記錄圖；

圖8B為本發(fā)明方法壓制圖8A諧波噪聲后剖面；圖8C為本發(fā)明方法獲得的圖8A諧波噪聲剖面；

圖9A為圖8A局部放大顯示結(jié)果；圖9B為圖8B局部放大顯示結(jié)果；圖9C為圖8C局部放大顯示結(jié)果。

具體實施方式

為了使本發(fā)明的目的、技術(shù)方案更加清楚明白，下面結(jié)合附圖和具體實施方式對本發(fā)明做進一步詳細的說明。在此，本發(fā)明的示意性實施方式及其說明用于解釋本發(fā)明，但并不作為本發(fā)明的限定。

在本發(fā)明實施例中，提出了一種基于波形形態(tài)特征稀疏化建模的諧波噪聲壓制方法，包括以下內(nèi)容：

1)、根據(jù)可控震源滑動掃描采集地震記錄中諧波噪聲波形形態(tài)特征構(gòu)造Chirplet變換，并與連續(xù)小波變換構(gòu)成超完備字典；

2)、Chirplet正、反變換的快速實現(xiàn)方法；

3)、基于相關(guān)后諧波噪聲時頻分布特征確定Chirplet變換參數(shù)；

4)、根據(jù)參考掃描信號的起始頻率確定諧波噪聲的濾波截止頻率，實現(xiàn)有效信號與諧波噪聲保真分離。

步驟1)中根據(jù)可控震源滑動掃描采集地震記錄中諧波噪聲波形形態(tài)特征構(gòu)造Chirplet變換，并與連續(xù)小波變換構(gòu)成超完備字典包括：

選擇連續(xù)小波變換作為稀疏表示有效信號成分的變換字典：

式中WT_x(A,τ)為變換系數(shù)，A表示尺度因子，x(t)表示待分析信號，ψ(t)表示Morlet母小波；其中t為時間，τ為平移量，*表示共軛。

連續(xù)小波變換的反變換為：

式中常數(shù)C_Ψ＜∞為其容許條件。

Chirplet正變換的定義為：

CT_x(t,f,a,c,d)＝∫x(τ)h^*(τ-t,f,a,c,d)dτ，

式中CT_x(t,f,a,c,d)表示待分析信號的Chirplet變換值，x(t)表示待分析信號，a為伸縮參數(shù)，c為時間上的線性調(diào)頻率，f為頻率，d為頻率上的線性調(diào)頻率；τ為時間平移量，*表示共軛；h(t,f,a,c,d)為Chirplet變換的核函數(shù)，其表達式如下：

式中g(shù)(t)為高斯函數(shù)；v為平移量。

忽略頻率上的線性調(diào)頻操作，核函數(shù)退化為：

選擇構(gòu)造的Chirplet變換作為稀疏表示諧波噪聲的字典。

如圖1與圖2所示，退化后的Chirplet變換原子與諧波噪聲波形形態(tài)特征比較匹配，能更有效表示諧波噪聲。

根據(jù)形態(tài)成分分析理論，用上面選定的字典A₁即連續(xù)小波變換和A₂即Chirplet變換構(gòu)成超完備字典稀疏表示信號計算稀疏表示系數(shù)：

式中：x₁為重構(gòu)系數(shù)中與A₁對應(yīng)的部分；x₂為重構(gòu)系數(shù)中與A₂對應(yīng)的部分。λ為拉格朗日乘子。該優(yōu)化問題可以通過坐標分塊松弛算法求解。

如圖3所示，步驟2)中Chirplet正、反變換的快速實現(xiàn)流程包括如下步驟：

步驟101：從原始數(shù)據(jù)中讀出單道時域信號，記為x；

步驟102：計算得到信號x的Chirplet正變換結(jié)果，記為y，其實現(xiàn)方法類似于短時傅里葉變換；

步驟103：計算y的Chirplet反變換結(jié)果，記為

由于步驟103中Chirplet反變換得到的重構(gòu)結(jié)果與原始信號相比會在幅度上存在一定的差別，因此接下來使用一個模型信號計算修正系數(shù)來進行幅度修正。

步驟201：構(gòu)造一個線性升頻信號模型，記為s，其長度與信號x相同；

步驟202：計算得到信號s的Chirplet正變換結(jié)果，記為s₀，其實現(xiàn)方法類似于短時傅里葉變換；

步驟203：計算得到s₀的Chirplet反變換結(jié)果，記為

步驟204：根據(jù)模型反變換結(jié)果計算修正系數(shù)系數(shù)，即其中||.||₂表示二范數(shù)；

步驟104：將修正系數(shù)r作用到單道時域信號x的Chirplet反變換x，得到修正的Chirplet反變換結(jié)果，記為則

通過本發(fā)明的Chirplet快速正反變換，計算效率較高，且修正系數(shù)使得重構(gòu)精確度更高。

如圖4所示，步驟3)中基于相關(guān)后數(shù)據(jù)時頻分布特征確定Chirplet變換參數(shù)的方法包括如下步驟：

步驟301：將單道時域信號x轉(zhuǎn)換到頻域，記為F；

步驟302：確定頻域的高低頻分界，記為κ；

步驟303：計算高頻(頻率分界以上)能量占總頻帶能量的比值，即振幅譜比值，記為b；

步驟304：根據(jù)振幅譜比值計算Chirplet變換的伸縮參數(shù)，記為a，線性調(diào)頻參數(shù)，記為c以及閾值權(quán)系數(shù)參數(shù)，記為λ：

a＝k₁×b，λ＝k₂-b，c＝-k₃×b×a²，

其中，b為振幅譜比值；k₁,k₂,k₃分別為參數(shù)a、c與λ的修正系數(shù)。

本發(fā)明實施例中，高、頻分界選擇為40Hz，k₁,k₂,k₃分別?。簁₁＝2.0,k₂＝1.5,k₃＝5。

通過本發(fā)明的確定Chirplet變換參數(shù)的方法，能夠各道數(shù)據(jù)得到不同的參數(shù)，從而具有較強的適應(yīng)性，可以處理含不同強弱程度諧波噪聲的地震記錄數(shù)據(jù)。

如圖5所示，步驟4)中基于參考掃描信號起始頻率的諧波噪聲保真分離方法包括如下步驟：

步驟401：根據(jù)參考掃描信號的起始頻率確定諧波噪聲的起始頻率，記為f_b；

步驟402：諧波噪聲做截止頻率為f_b的高通濾波；

步驟403：原始數(shù)據(jù)減去濾波后的諧波噪聲得到有效信號，實現(xiàn)有效信號與諧波噪聲的保真分離。

使用基于本發(fā)明的保真分離方法，能夠更加有效地保護有效信號，降低對有效信號的損傷。

在本例中，提出了一種基于波形形態(tài)特征稀疏化建模的諧波噪聲壓制方法，從而達到了有效壓制可控震源滑動掃描記錄地震數(shù)據(jù)中諧波噪聲的目的。本發(fā)明具有如下有益效果：

1)本發(fā)明方法根據(jù)諧波噪聲的時頻分布特征構(gòu)造合適的變換，保證了稀疏性；

2)Chirplet正、反變換的快速實現(xiàn)方式保證了變換的效率，使其可用于海量實際資料處理；

3)Chirplet反變換使用修正系數(shù)保證了結(jié)果的準確性；

4)根據(jù)數(shù)據(jù)驅(qū)動自動確定Chirplet變換的參數(shù)，具有很強的適應(yīng)性，且單道計算，便于并行處理；

5)本發(fā)明對有效信號具有高保真性，能夠有效的保護有效波的高、低頻成分。

下面結(jié)合一個具體實施例對上述變換構(gòu)造及參數(shù)確定方法進行具體說明，且該實施例僅是為了更好說明本發(fā)明，并不構(gòu)成對本發(fā)明的不當限定。

利用本發(fā)明提供的分析處理方法應(yīng)用到合成模擬與實際滑動掃描地震記錄數(shù)據(jù)中驗證本發(fā)明構(gòu)造Chirplet變換及參數(shù)的方法壓制諧波噪聲的有效性。本發(fā)明的方法不僅能夠有效的壓制諧波噪聲，而且有效信號具有較高保真性。

圖6A-圖6B為針對滑動掃描地震記錄特性，給出的未受諧波干擾的合成記錄以及受諧波干擾的合成記錄。該模型包含三個反射層，所用掃描信號為線性升頻掃描信號，其最低頻率為3Hz，最高頻率為90Hz，掃描信號時長16s，有效信號包含直達波以及反射波，諧波噪聲包含二次諧波和三次諧波。有效信號受到很強的諧波噪聲干擾，尤其是深層反射波受到諧波噪聲干擾更為突出。

如圖7A-圖7B所示，本發(fā)明方法對該含諧波噪聲的炮記錄(圖6B)進行處理，得到對應(yīng)的有效信號剖面和對應(yīng)的諧波噪聲剖面。本發(fā)明方法可以有效的對諧波噪聲進行壓制，諧波噪聲中幾乎沒有殘留有效信號能量，表明本發(fā)明的方法能夠有效壓制諧波噪聲，而且對有效信號具有高保真性。

圖8A為我國西部某油田基于滑動掃描觀測到的相關(guān)后炮集數(shù)據(jù)，接下來驗證本發(fā)明方法處理實際地震資料的有效性。該炮集數(shù)據(jù)總共401道，采樣時間間隔為2ms，記錄長度為6s。從剖面上可以看出，該記錄中的信號受到很強的諧波噪聲干擾，導(dǎo)致有效信號被諧波噪聲覆蓋，對地震資料的分析以及解釋造成嚴重影響。

如圖8B-圖8C所示，應(yīng)用本發(fā)明的方法對該炮集數(shù)據(jù)進行處理獲得其有效信號和諧波噪聲，本發(fā)明方法有效地壓制了地震記錄中的諧波噪聲。并且如圖8B中方框區(qū)域指示處被諧波噪聲覆蓋的有效信號在有效信號剖面中清晰地顯示。

對矩形框區(qū)域進行放大顯示，如圖9A-9C所示，在諧波噪聲剖面中幾乎不含有效信號能量，表明本發(fā)明的方法對有效信號具有較高保真性。

以上的模型和實際資料算例中，利用本發(fā)明的方法構(gòu)造諧波噪聲的稀疏表示變換及確定變換的參數(shù)，對地震資料數(shù)據(jù)進行諧波噪聲壓制，不僅能夠有效的壓制諧波噪聲，而且有效信號具有較高保真性，為后續(xù)資料的分析奠定基礎(chǔ)。

最后需要說明的是，以上模型和實際資料算例對本發(fā)明的目的，技術(shù)方案以及有益效果提供了進一步的驗證，這僅屬于本發(fā)明的具體實施算例，并不用于限定本發(fā)明的保護范圍。凡在本發(fā)明的精神和原則之內(nèi)，所做的任何修改，改進或等同替換等，均應(yīng)在本發(fā)明的保護范圍內(nèi)。