本發(fā)明屬于地球物理技術(shù)領(lǐng)域，尤其是指針對包含直線型同相軸的地震記錄提出了一種基于多層分?jǐn)?shù)階傅里葉霍夫變換的自適應(yīng)時頻峰值濾波消噪方法，該方法通過一種新型的霍夫變換獲取地震數(shù)據(jù)中的同相軸包絡(luò)斜率，進(jìn)而建立具有該斜率的濾波軌線方程并沿著這些軌線對含噪地震數(shù)據(jù)進(jìn)行樣本重采樣，其目的是使得采樣后數(shù)據(jù)中有效信號線性度充分降低，最后對重采樣后數(shù)據(jù)進(jìn)行時頻峰值濾波有效壓制噪聲，為地震勘探信號消噪處理及地震數(shù)據(jù)中有效信號識別提供了新方法。

背景技術(shù)：

地震勘探作為一種重要的物理勘探手段，在油田的開發(fā)中起著重要的作用。野外地震數(shù)據(jù)采集和地震資料的數(shù)字處理技術(shù)是地震勘探的兩個重要環(huán)節(jié)。只有得到高質(zhì)量的地震剖面圖，才能為石油勘探提供可靠的依據(jù)。在利用檢波器進(jìn)行實(shí)際的野外數(shù)據(jù)采集過程中，有效地震信號經(jīng)常會受到各種噪聲的干擾(規(guī)則噪聲和隨機(jī)噪聲)，這些都直接影響地震解釋工作的準(zhǔn)確性。相比之下，隨機(jī)噪聲是不規(guī)則的、沒有規(guī)律的，且相鄰道間彼此是互不相關(guān)的，它們沒有固定的頻率，幾乎分布于整個頻帶，嚴(yán)重影響著地震記錄的信噪比。因此隨機(jī)噪聲壓制是地震數(shù)據(jù)處理中的重點(diǎn)和難點(diǎn)。

近幾年來，對于隨機(jī)噪聲壓制，除了目前發(fā)展較為成熟的幾種消噪算法，如自適應(yīng)濾波、經(jīng)驗(yàn)?zāi)B(tài)分解、多項(xiàng)式擬合、f-x預(yù)測濾波、小波變換、維納濾波等，時頻峰值濾波算法被成功地應(yīng)用于陸地地震勘探隨機(jī)噪聲壓制，并先后發(fā)展出多種改進(jìn)算法，如局部徑向軌線時頻峰值濾波、二次軌線時頻峰值濾波及分段軌線時頻峰值濾波等。時頻峰值濾波方法較上述成熟算法在強(qiáng)隨機(jī)噪聲、非平穩(wěn)信號及先驗(yàn)信息未知等條件下存在明顯的優(yōu)勢，但現(xiàn)有時頻峰值濾波方法中軌線參數(shù)多針對某一地震記錄人為設(shè)置，靈活性較差。

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

本發(fā)明提供一種基于改進(jìn)霍夫變換的自適應(yīng)時頻峰值濾波消噪方法，以解決現(xiàn)有軌線型時頻峰值濾波方法中存在的軌線參數(shù)及形狀固定，靈活性和適應(yīng)性較差問題。目的在于針對包含直線型同相軸的地震記錄，提出一種改進(jìn)的霍夫變換，并將其應(yīng)用與地震記錄中的同相軸形態(tài)識別，從而實(shí)現(xiàn)時頻峰值濾波方法中最優(yōu)軌線的自適應(yīng)選取。

本發(fā)明采取的技術(shù)方案是，包括下列步驟：

步驟一：一張含噪的、包含直線型同相軸的地震記錄，記錄中的反射同相軸就相當(dāng)于圖像中的邊緣信息，為了便于對地震數(shù)據(jù)進(jìn)行形態(tài)識別，首先對地震數(shù)據(jù)進(jìn)行灰度處理，將數(shù)據(jù)值轉(zhuǎn)換到0-255之間，處理后的地震數(shù)據(jù)為一張圖像數(shù)據(jù)，接著，利用改進(jìn)的霍夫變換實(shí)現(xiàn)地震記錄中同相軸形態(tài)識別，具體過程為：先利用公式對灰度處理后的圖像數(shù)據(jù)進(jìn)行多層分?jǐn)?shù)階傅里葉變換，將地震數(shù)據(jù)變換到時域和頻域之間的一個中間域，將其旋轉(zhuǎn)到從時間軸指向頻率軸的一個角度，得到多組不同階數(shù)的傅里葉變換結(jié)果；再將變換結(jié)果轉(zhuǎn)換成由模值ρ和相角θ構(gòu)成的極坐標(biāo)形式，在極坐標(biāo)內(nèi)，對數(shù)據(jù)進(jìn)行共軛對稱，并沿模值方向進(jìn)行一維傅里葉反變換便可得到一張由多條正弦線構(gòu)成的正弦譜圖；實(shí)際上，從多層分?jǐn)?shù)階傅里葉變換到轉(zhuǎn)換成極坐標(biāo)，再到進(jìn)行一維傅里葉反變換，整個過程相當(dāng)于對地震數(shù)據(jù)進(jìn)行了一次霍夫變換，得到的正弦譜圖中多條正弦線的交點(diǎn)，即正弦譜矩陣中的數(shù)據(jù)最大值就對應(yīng)了地震數(shù)據(jù)中的同相軸；對于直線型同相軸的地震記錄，通過搜索正弦譜圖中的峰值便可得到與其對應(yīng)的地震數(shù)據(jù)中的某條同相軸在極坐標(biāo)內(nèi)的模值ρpeak和相角θpeak，再根據(jù)極坐標(biāo)與原圖像坐標(biāo)之間的關(guān)系式ρ＝xcosθ+ysinθ便可以得到地震數(shù)據(jù)內(nèi)同相軸包絡(luò)線的斜率kl＝ctanθpeak；于是沿著多條斜率為kl的直線型濾波軌線對含噪地震記錄進(jìn)行樣本重采樣，有效提高采樣后信號的線性度，降低濾波偏差；

步驟二：根據(jù)檢測到的同相軸包絡(luò)斜率kl建立與同相軸走勢接近的濾波軌線方程yl＝klxl+bl，其中xl為地震數(shù)據(jù)點(diǎn)的橫坐標(biāo)，代表炮檢距方向，yl為地震數(shù)據(jù)點(diǎn)的縱坐標(biāo)，代表時間方向，kl為軌線斜率，bl為截距，決定軌線與地震數(shù)據(jù)縱坐標(biāo)的交點(diǎn)位置，改變截距bl的數(shù)值便會生成多條具有相同斜率的直線軌線，沿著這些軌線對地震數(shù)據(jù)進(jìn)行樣本重采樣，便可實(shí)現(xiàn)對地震數(shù)據(jù)沿同相軸形態(tài)方向的數(shù)據(jù)提取；

步驟三：對采樣后的數(shù)據(jù)樣本，利用時頻峰值濾波方法進(jìn)行噪聲壓制，由于對地震數(shù)據(jù)進(jìn)行樣本重采樣時所采用的軌線斜率是通過檢測同相軸形態(tài)得出的，因此軌線與同相軸趨勢接近，沿著這些軌線重采樣后的地震數(shù)據(jù)有效信號線性度充分提高，主頻充分降低，因此對其進(jìn)行時頻峰值濾波時所采用的濾波窗長比傳統(tǒng)時頻峰值濾波方法中長很多，較長的窗長在完整保持有效成分的同時可以有效地壓制隨機(jī)噪聲。

本發(fā)明所述步驟一中對于直線型同相軸，地震數(shù)據(jù)中的某一點(diǎn)(x,y)與極坐標(biāo)內(nèi)的參數(shù)對(ρ，θ)滿足如下關(guān)系：ρ＝xcosθ+ysinθ,θ∈[0,π]，其中θ為直線法線與x軸的夾角，ρ為圖像空間中原點(diǎn)到直線的距離。

本發(fā)明以提高時頻峰值濾波方法的適用性和靈活性，有效壓制地震記錄中隨機(jī)噪聲為目的，針對包含直線型同相軸地震記錄，基于多層分?jǐn)?shù)階傅里葉提出了一種改進(jìn)的霍夫變換，并將其應(yīng)用于同相軸形態(tài)檢測，實(shí)現(xiàn)了時頻峰值濾波方法中最優(yōu)濾波軌線的自適應(yīng)選取。改進(jìn)的霍夫變換與現(xiàn)有霍夫變換相比具有較高的分辨率，沒有梯度信息取向的損失，忽略了從笛卡爾到極坐標(biāo)映射的零填充，在同相軸形態(tài)檢測方面計算速度快，準(zhǔn)確性高。

多層分?jǐn)?shù)階傅里葉變換在頻域內(nèi)具有較高的分辨率，沒有梯度信息取向的損失，忽略了從笛卡爾到極坐標(biāo)映射的零填充，計算速度明顯提高。而現(xiàn)有的隨機(jī)霍夫變換有時不能檢測到短線段；線段檢測器對噪聲比較敏感，在嘈雜的圖像中檢測性能明顯下降。

在重采樣過程中由于位移-時間域內(nèi)的樣本格點(diǎn)與新作用域內(nèi)的樣本格點(diǎn)不一致，部分樣本點(diǎn)不能落到具體的地震道上。如果丟棄這些“脫道”的點(diǎn)，采樣得到的序列中就會出現(xiàn)一些斷點(diǎn)和尖峰脈沖。本發(fā)明通過在采樣前在相鄰地震道間插入一些虛擬地震道的方式對地震數(shù)據(jù)進(jìn)行擴(kuò)充，使得原本“脫道”的點(diǎn)能夠落在具體的地震道上，從而有效降低采樣誤差。

本發(fā)明的有益效果是提出了一種基于多層分?jǐn)?shù)階傅里葉變換實(shí)現(xiàn)的改進(jìn)霍夫變換方法，并將其應(yīng)用于直線型地震記錄中同相軸形態(tài)檢測，獲取同相軸斜率。改進(jìn)的霍夫變換與現(xiàn)有霍夫變換相比具有較高的分辨率，沒有梯度信息取向的損失，忽略了從笛卡爾到極坐標(biāo)映射的零填充，在同相軸形態(tài)檢測方面計算速度快，準(zhǔn)確性高。同相軸形態(tài)的準(zhǔn)確識別軌線型時頻峰值濾波方法中最優(yōu)濾波軌線的自適應(yīng)選取為提供了先決條件。此外，沿濾波軌線進(jìn)行的插值重采樣方案通過在相鄰地震道間插入虛擬地震道的方式，有效避免了采樣點(diǎn)的脫道現(xiàn)象，降低了重采樣過程造成的誤差。以上方面均對充分降低濾波偏差，完整保持有效信號具有重要意義。

附圖說明

圖1(a)圖像空間內(nèi)直線各參數(shù)情況，其中點(diǎn)p1、p2和p3為直線上的三個點(diǎn)，直線與坐標(biāo)原點(diǎn)距離為ρ，直線與y軸夾角為θ；

圖1(b)霍夫變換參數(shù)空間內(nèi)正弦曲線譜圖，其中三條正弦曲線與點(diǎn)p1、p2和p3三點(diǎn)對應(yīng)，三條線相交于一點(diǎn)，該點(diǎn)的模值和相角即對應(yīng)了圖1a中直線與坐標(biāo)原點(diǎn)的距離ρ及其與y軸夾角θ；

圖2(a)具有不同噪聲強(qiáng)度的含噪地震記錄，該記錄中包含三條直線型同相軸，其中區(qū)域r1信噪比較低，有效信號被淹沒其中難以辨識；區(qū)域r2包含較強(qiáng)的脈沖噪聲；區(qū)域e1存在交叉軸；

圖2(b)變換后得到的參數(shù)空間的正弦譜圖，其中存在三個明顯的交點(diǎn)(黑色方框標(biāo)出)，分別對應(yīng)原地震數(shù)據(jù)中的三條直線型同相軸；

圖2(c)變換后得到的參數(shù)空間的正弦譜圖的三維圖，其中的三個尖峰就是圖2b正弦譜圖中的三個交點(diǎn)，通過搜索峰值便可找到正弦譜圖中的交點(diǎn)位置；

圖2(d)含噪記錄及同相軸檢測結(jié)果，圖中黑色線給出了三條同相軸的檢測結(jié)果，可以看到檢測結(jié)果與同相軸幾乎完全吻合；

圖3重采樣過程中虛擬地震道插值示意圖，其中x方向?yàn)榈卣鸬婪较?，xi-1，xi，xi+1，xi+2為實(shí)際的地震道，其間虛線為插入的虛擬地震道；t方向?yàn)闀r間方向，點(diǎn)p1、p2和p3為與時刻ti-1，ti，ti+1對應(yīng)的某一濾波軌線上的三個點(diǎn)，可見點(diǎn)p2沒有落到具體的地震道上，而點(diǎn)p4為對點(diǎn)p2橫坐標(biāo)四舍五入對應(yīng)的點(diǎn)，插入虛擬地震道后，點(diǎn)p2可以落到某一具體的虛擬地震道上，有效提高了數(shù)據(jù)采樣的準(zhǔn)確性；

圖4(a)實(shí)際的含噪地震記錄，中包含的同相軸呈現(xiàn)直線形態(tài)；

圖4(b)傳統(tǒng)時頻峰值濾波消噪結(jié)果，其濾波方向沿著時間方向，選取的窗長為9，濾波后同相軸可見；

圖4(c)基于形態(tài)檢測的自適應(yīng)時頻峰值濾波消噪結(jié)果，其濾波方向沿著檢測到的同相軸方向，濾波窗長為32，濾波后隨機(jī)噪聲得到了有效壓制，同相軸清晰連貫。

具體實(shí)施方式

包括下列步驟：

步驟一：一張含噪的地震記錄可以看成是一張圖像數(shù)據(jù)，記錄中的反射同相軸就相當(dāng)于圖像中的邊緣信息。為了便于對地震數(shù)據(jù)進(jìn)行形態(tài)識別，首先對地震數(shù)據(jù)進(jìn)行灰度處理，將數(shù)據(jù)值轉(zhuǎn)換到0-255之間，處理后的地震數(shù)據(jù)可看成是一張圖像數(shù)據(jù)；接著利用改進(jìn)霍夫變換對灰度化的圖像數(shù)據(jù)中的同相軸形態(tài)進(jìn)行識別，首先，利用式(1)中的分?jǐn)?shù)階傅里葉變換公式將地震數(shù)據(jù)變換到時域和頻域之間的一個中間域，將其旋轉(zhuǎn)到從時間軸指向頻率軸的一個角度；

其中f^α(k1,k2)的n×n個頻率分散在[-απ,απ]×[-απ,απ]區(qū)間內(nèi)，給定一個分?jǐn)?shù)系數(shù)α，就會產(chǎn)生一系列的分?jǐn)?shù)階傅里葉頻率，于是有

pi{(αik1,αik2)|-n/2≤k1,k2≤n/2}，(2)

其中i＝1,2,…,l，l是層的個數(shù)，整個頻域內(nèi)的數(shù)據(jù)點(diǎn)為l層數(shù)據(jù)點(diǎn)的合并，

l越大，得到的頻率分辨率越高，對含噪地震數(shù)據(jù)圖像進(jìn)行多層分?jǐn)?shù)階傅里葉變換后，得到多個不同階數(shù)的傅里葉變換頻域散點(diǎn)圖；

接著，將變換得到的散點(diǎn)圖轉(zhuǎn)換成由模值ρ和相角θ構(gòu)成的極坐標(biāo)形式，在極坐標(biāo)內(nèi)，對數(shù)據(jù)進(jìn)行共軛對稱，并沿模值方向進(jìn)行一維傅里葉反變換便可得到一張由多條正弦線構(gòu)成的正弦譜圖；實(shí)際上，從多層分?jǐn)?shù)階傅里葉變換到轉(zhuǎn)換成極坐標(biāo)，再到進(jìn)行一維傅里葉反變換，整個過程相當(dāng)于對地震數(shù)據(jù)進(jìn)行了一次霍夫變換過程；

得到的正弦譜圖中多條正弦線的交點(diǎn)(即正弦譜矩陣中的數(shù)據(jù)最大值)就對應(yīng)了地震數(shù)據(jù)中的同相軸；于是，對于直線型同相軸，通過搜索正弦譜圖中峰值便可得到與其對應(yīng)的地震數(shù)據(jù)中的某一線條在極坐標(biāo)內(nèi)的模值ρpeak和相角θpeak，再根據(jù)極坐標(biāo)與原圖像坐標(biāo)之間的關(guān)系式ρ＝xcosθ+ysinθ便可以得到地震數(shù)據(jù)內(nèi)同相軸包絡(luò)線的斜率kl＝ctgθpeak，于是就可以沿著多條斜率為kl的直線對含噪地震記錄進(jìn)行樣本重采樣，使得采樣軌跡充分逼近同相軸形態(tài)，從而有效降低采樣后得到信號的線性度，降低濾波偏差；

對于直線型同相軸，地震數(shù)據(jù)中的某一點(diǎn)(x,y)與極坐標(biāo)內(nèi)的參數(shù)對(ρ，θ)滿足如下關(guān)系：

ρ＝xcosθ+ysinθ,θ∈[0,π](4)

其中θ為直線法線與x軸的夾角，ρ為圖像空間中原點(diǎn)到直線的距離，圖1給出了圖像空間與參數(shù)空間的對應(yīng)關(guān)系示意圖。圖1(a)為圖像數(shù)據(jù)空間中某一直線上的點(diǎn)p1、p2和p3，圖1(b)為這三個點(diǎn)對應(yīng)的極坐標(biāo)空間的三條正弦曲線，這三條正弦線交于坐標(biāo)為(ρ1,θ1)的點(diǎn)p，通過點(diǎn)p的模值和相角信息便可得到原圖像數(shù)據(jù)中的直線(即同相軸)的斜率；

為了驗(yàn)證改進(jìn)的多層分?jǐn)?shù)階傅里葉霍夫變換方法在地震記錄同相軸形態(tài)檢測中的有效性，圖2(a)給出了一張40道，每道1024個采樣點(diǎn)的人工合成地震記錄，各道加入不同強(qiáng)度的隨機(jī)噪聲得到信噪比不同的含噪記錄，信噪比范圍由-8db到0db。該記錄中包含三條直線型同相軸，其中兩條相交于e1區(qū)域，在r2區(qū)域存在一段突變噪聲，在r1區(qū)域有效信號淹沒于強(qiáng)隨機(jī)噪聲當(dāng)中，信噪比較低。圖2(b)中的蝶型圖中存在三個明顯的蝶型中心點(diǎn)，分別對應(yīng)地震記錄中的三條同相軸。圖2(c)是蝶型圖的一個三維圖，圖中的三個明顯尖峰就對應(yīng)蝶型圖中的三個蝶型中心點(diǎn)。從參數(shù)空間的三維立體圖可以看出，累加器矩陣存在三個明顯峰值，選取其中峰值最大的前三個，即可恢復(fù)出圖像域內(nèi)相應(yīng)的三條同相軸包絡(luò)線，如圖2(d)中黑色線所示。

步驟二：根據(jù)檢測到的同相軸包絡(luò)斜率kl建立與同相軸形態(tài)接近的濾波軌線方程

yl＝klxl+bl,(5)

其中xl為地震數(shù)據(jù)點(diǎn)的橫坐標(biāo)，代表炮檢距方向，yl為地震數(shù)據(jù)點(diǎn)的縱坐標(biāo)，代表時間方向，kl為軌線斜率，bl為截距，決定軌線與地震數(shù)據(jù)縱坐標(biāo)的交點(diǎn)位置；不斷改變截距bl的數(shù)值便會生成多條具有相同斜率的直線軌線，沿著這些軌線對地震數(shù)據(jù)進(jìn)行樣本重采樣，便可實(shí)現(xiàn)對地震數(shù)據(jù)沿同相軸形態(tài)方向數(shù)據(jù)提?。?/p>

在重采樣過程中由于位移-時間域內(nèi)的樣本格點(diǎn)與新作用域內(nèi)的樣本格點(diǎn)不一致，部分樣本點(diǎn)不能落到具體的地震道上，如果丟棄這些“脫道”的點(diǎn)，采樣得到的序列中就會出現(xiàn)一些斷點(diǎn)和尖峰脈沖，本發(fā)明通過在采樣前在相鄰地震道間插入一些虛擬地震道的方式對地震數(shù)據(jù)進(jìn)行擴(kuò)充，使得原本“脫道”的點(diǎn)能夠落在具體的地震道上，從而有效降低采樣誤差，如圖3所示，圖中虛線是濾波軌線，橫坐標(biāo)為地震道方向，縱坐標(biāo)為時間方向，點(diǎn)p1、p3為時刻點(diǎn)ti-1和ti+1對應(yīng)的軌線上的樣本點(diǎn)，它們分別落在地震道xi-1和xi+2上，從圖中可以看出，菱形點(diǎn)p2在插值后可以落到了某一虛擬道上，在采樣過程中它就不會被丟失，也不再利用其它點(diǎn)進(jìn)行代替，使得采樣得到的序列更準(zhǔn)確；

采樣過程首先在地震數(shù)據(jù)中選取一個樣本點(diǎn)，設(shè)其縱坐標(biāo)為將其代入式(5)中的軌線方程計算出相應(yīng)的橫坐標(biāo)即可提取一個坐標(biāo)為的數(shù)據(jù)點(diǎn)；再選取縱坐標(biāo)代入式(5)中的軌線方程計算出相應(yīng)的橫坐標(biāo)xli+1，提取第二個樣本點(diǎn)依次類推。直到計算的橫坐標(biāo)達(dá)到最大的地震道數(shù)為止，這樣就可提取得到一個樣本序列。接著依次提取第二個、第三個…樣本序列。重采樣后得到的是一個包含多條序列的數(shù)據(jù)矩陣；

步驟三：選取合適窗長，利用時頻峰值濾波方法對沿多條軌線采樣后的重采樣數(shù)據(jù)進(jìn)行噪聲壓制，由于對地震數(shù)據(jù)進(jìn)行樣本重采樣時所采用的軌線斜率是通過檢測同相軸形態(tài)得出的，因此軌線與同相軸趨勢接近，沿著這些軌線重采樣的地震數(shù)據(jù)中有效信號線性度充分提高，主頻充分降低，因此對其進(jìn)行時頻峰值濾波時所采用的濾波窗長比傳統(tǒng)時頻峰值濾波方法中要長很多，較長的窗長在完整保持有效成分的同時可以有效地壓制隨機(jī)噪聲。選取較大窗長對重采樣后數(shù)據(jù)序列進(jìn)行時頻峰值濾波，實(shí)現(xiàn)對含噪地震記錄的消噪處理。

下面給出具體的實(shí)驗(yàn)例子來說明本發(fā)明的實(shí)驗(yàn)效果：

圖4(a)為某一地震剖面部分記錄，其中包含大量隨機(jī)噪聲，有效同相軸淹沒其中，該記錄中同相軸大都呈現(xiàn)直線型形態(tài)分布。分別利用傳統(tǒng)時頻峰值濾波算法和基于改進(jìn)霍夫變換的自適應(yīng)時頻峰值濾波算法對該含噪地震記錄進(jìn)行噪聲壓制處理，結(jié)果分別如圖4(b)和圖4(c)所示。通過觀察兩種方法的消噪結(jié)果可以看出，兩種方法去噪后，隨機(jī)噪聲都得到了壓制，同相軸變得清晰。但相比之下，本發(fā)明中方法消噪后，同相軸變得更加連貫，噪聲壓制效果更為理想。該實(shí)驗(yàn)結(jié)果驗(yàn)證了自適應(yīng)軌線時頻峰值濾波算法的有效性。