張 華 刁 塑 溫建亮 黃光南 朱雯婷 白 敏

(①東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西南昌 330013；②山西省煤炭地質(zhì)物探測(cè)繪院，山西晉中030600；③華北水利水電大學(xué)資源與環(huán)境學(xué)院，河南鄭州 450046)

0 引言

由于表層激發(fā)條件的限制以及其他各種不可抗干擾因素的影響，僅在野外采集階段完全壓制隨機(jī)噪聲不太現(xiàn)實(shí)，需要在室內(nèi)用有效的數(shù)學(xué)去噪方法提高疊前地震資料的信噪比，以滿足不同處理環(huán)節(jié)的需求[1-4]。目前，有許多行之有效的隨機(jī)噪聲壓制方法，其中稀疏變換應(yīng)用較為廣泛，包括傅里葉變換[5-8]、小波變換[9-10]、Radon變換[11-12]、Contourlet變換[13]、曲波變換[14]以及Seislet 變換[15]等。該類方法根據(jù)隨機(jī)噪聲和有效波在稀疏變換域中能量分配的差異，對(duì)稀疏系數(shù)采用適當(dāng)?shù)拈撝等コS機(jī)噪聲的影響。為了取得較好的去噪效果，稀疏基需要盡可能地捕獲地震波的有效信息，并且少數(shù)較大的稀疏系數(shù)能夠代表信號(hào)的主要特征，而大部分較小、被濾除的系數(shù)不影響原始數(shù)據(jù)的主要特征[16]。曲波變換能夠有效地表示地震數(shù)據(jù)的尺度性和方向性，更加稀疏地表示地震波場(chǎng)局部細(xì)節(jié)特征，所以許多學(xué)者選用曲波變換方法去噪，且效果顯著[17-20]。

現(xiàn)有曲波變換去噪方法的前提是地震數(shù)據(jù)為均勻網(wǎng)格采樣，而對(duì)于非均勻網(wǎng)格采樣的含噪數(shù)據(jù)則效果不佳。由于野外地形條件及施工環(huán)境的制約，在很多情況下，地震數(shù)據(jù)常為非均勻采樣[21]。如果在去噪過(guò)程中將非均勻采樣視為均勻采樣進(jìn)行處理，就得不到連續(xù)的地震波場(chǎng)，自然也不能充分壓制噪聲干擾，從而會(huì)影響到后續(xù)其他處理方法的應(yīng)用。Hennenfent等[22]提出基于非均勻采樣的曲波變換重建方法，但未討論該方法在噪聲壓制中的應(yīng)用。

本文在前人研究的基礎(chǔ)上，針對(duì)非均勻采樣地震數(shù)據(jù)，在傳統(tǒng)曲波變換的基礎(chǔ)上，引入非均勻傅里葉變換，建立均勻曲波系數(shù)與空間非均勻采樣地震道之間的規(guī)則化反演算子；然后選擇合適的噪聲估計(jì)值，使用線性Bregman方法反演，在迭代過(guò)程中采用軟閾值對(duì)曲波系數(shù)去噪，由反演得到無(wú)噪聲的均勻曲波系數(shù)；再進(jìn)行常規(guī)曲波反變換，得到去噪后的地震數(shù)據(jù)，由此形成基于非均勻曲波變換和線性Bregman方法的隨機(jī)噪聲壓制方法。模型數(shù)據(jù)和實(shí)際數(shù)據(jù)的實(shí)驗(yàn)結(jié)果表明了本文方法的有效性。

1 非均勻曲波變換

Candès等[23]提出了第二代曲波變換，從而使曲波變換更容易被理解、運(yùn)算效率更高、實(shí)現(xiàn)更簡(jiǎn)單，能夠?yàn)榈卣鹦盘?hào)提供最優(yōu)的稀疏表示方式。實(shí)際上，實(shí)現(xiàn)快速離散曲波變換主要包括兩個(gè)步驟：①對(duì)地震數(shù)據(jù)應(yīng)用二維傅里葉變換，得到頻率—波數(shù)域系數(shù)；②在頻率—波數(shù)域形成角度楔形，將每一個(gè)楔形圍繞到原點(diǎn)重新裝配，并對(duì)每一個(gè)裝配好的楔形應(yīng)用二維傅里葉反變換，得到離散曲波系數(shù)。參考文獻(xiàn)[22]定義曲波正變換算子為

A=TF

(1)

式中：F實(shí)現(xiàn)了上述離散曲波變換第①步：T則實(shí)現(xiàn)了第②步。

定義曲波反變換算子為

AH=FHTH

(2)

式中：上標(biāo)“H”表示共軛轉(zhuǎn)置；FH表示二維傅里葉反變換，將頻率—波數(shù)域轉(zhuǎn)換到時(shí)間—空間域中；TH表示曲波平鋪算子，即將曲波系數(shù)變換到頻率—波數(shù)域的過(guò)程。

(3)

(4)

該算子描述了非均勻采樣下離散曲波系數(shù)與地震數(shù)據(jù)之間的關(guān)系。

由于快速離散曲波變換具有緊支撐性，正、反曲波變換前、后能量無(wú)損失，即

AAH=I

(5)

式中I為單位矩陣。對(duì)于新的非均勻曲波反變換算子C，由于包含非均勻快速傅里葉變換，那么算子C就不滿足式(5)。為此定義非均勻曲波正變換算子C+為

滿足 ‖d-Cx‖2≤σ

(6)

2 線性Bregman方法

式(6)最優(yōu)化問(wèn)題可轉(zhuǎn)為求解下述基追蹤規(guī)則化問(wèn)題[25]

(7)

式中λ是一個(gè)閾值權(quán)衡因子，平衡L1范數(shù)和L2范數(shù)的比重，可以在第一次迭代過(guò)程計(jì)算。線性Bregman方法求解迭代式為

(8)

式中：zk表示第k次迭代得到的曲波系數(shù)向量；lk為動(dòng)態(tài)步長(zhǎng)，定義為

(9)

軟閾值函數(shù)為

Sλ(x)=sign(x)max(|x|-λ,0)

(10)

(11)

(12)

該式可以保證在每次迭代過(guò)程中對(duì)含噪地震數(shù)據(jù)進(jìn)行噪聲壓制，從中也可以看出線性Bregman方法實(shí)現(xiàn)非常簡(jiǎn)單，沒(méi)有過(guò)多的調(diào)節(jié)參數(shù)。

3 噪聲壓制原理

一個(gè)含噪聲二維地震信號(hào)模型可以表示成為

d(i,j)=s(i,j)+e(i,j)

(13)

式中：s為真實(shí)信號(hào)；e為噪聲。去噪過(guò)程就是從非均勻含噪地震信號(hào)d中提取真實(shí)信號(hào)s，去除噪聲干擾信號(hào)e(圖1)。

圖1 非均勻曲波變換噪聲壓制流程

(14)

式中AH是常規(guī)的曲波反變換算子。通過(guò)對(duì)反演后的無(wú)噪聲均勻曲波系數(shù)應(yīng)用常規(guī)曲波反變換，可得到均勻采樣后的去噪地震數(shù)據(jù)。

4 數(shù)值模擬

信噪比定義為[29]

(15)

式中：s0表示原始均勻不含噪聲數(shù)據(jù)；s表示去噪后的均勻地震數(shù)據(jù)。R越高，表示去噪效果越好。

圖2a為合成的256道地震數(shù)據(jù)，有4個(gè)反射波同相軸，各反射同相軸能量有差異，采樣間隔為1ms，道間距為5m，每道1024個(gè)采樣點(diǎn)。加入標(biāo)準(zhǔn)差為0.05的高斯白噪聲，如圖2b所示，其信噪比為-5.11dB。對(duì)含噪數(shù)據(jù)進(jìn)行空間均勻傅里葉變換，然后再進(jìn)行空間非均勻傅里葉反變換，得到新的空間非均勻采樣下的256道地震數(shù)據(jù)，如圖2c所示，信噪比為-5.25dB，道間距范圍為0～10m。如果直接用均勻曲波變換法對(duì)圖2c數(shù)據(jù)去噪，則扭曲的同相軸得不到校正，導(dǎo)致去噪后的地震記錄誤差較大。因此首先采用非均勻曲波變換方法對(duì)圖2c進(jìn)行規(guī)則化處理，此時(shí)噪聲估計(jì)值σ=0，曲波變換所選擇的尺度數(shù)為5，在第二個(gè)最粗尺度上的角度數(shù)為32。規(guī)則化結(jié)果如圖2d所示，可以看見(jiàn)規(guī)則化后的地震波場(chǎng)非常連續(xù)，幾乎沒(méi)有視覺(jué)上的差異，并且從誤差剖面(圖2e)也可以看出規(guī)則化前后幾乎沒(méi)有誤差，說(shuō)明非均勻曲波變換方法規(guī)則化效果好。

采用均勻曲波變換方法對(duì)非均勻含噪數(shù)據(jù)(圖2c)直接去噪，對(duì)分解后的曲波系數(shù)進(jìn)行軟閾值處理，將保留下來(lái)的曲波系數(shù)進(jìn)行曲波反變換，從而得到最終的去噪結(jié)果，如圖3a所示，其信噪比為4.45dB。圖3b為去除的噪聲干擾，可以看出，其中含有部分有效波信號(hào)，主要原因是常規(guī)曲波變換不能直接處理非均勻采樣數(shù)據(jù)，而是把非均勻采樣數(shù)據(jù)當(dāng)成均勻采樣數(shù)據(jù)進(jìn)行處理。根據(jù)噪聲能量計(jì)算出噪聲估計(jì)值σ=2.5，采用本文方法在非均勻采樣數(shù)據(jù)規(guī)則化過(guò)程中同時(shí)去噪，結(jié)果如圖3c所示，其信噪比為5.11dB。去除的噪聲如圖3d所示?？梢钥闯?，本文方法不僅可以將非均勻采樣地震數(shù)據(jù)內(nèi)插為均勻采樣數(shù)據(jù)，還可以同時(shí)有效地去除隨機(jī)噪聲干擾，去噪后的地震同相軸更加連續(xù)、清晰，大幅度地提高了信噪比，且去除的噪聲中幾乎不含有效波能量，表明對(duì)于非均勻采樣地震數(shù)據(jù)本文方法具有較好的去噪效果。

(a)合成理論地震數(shù)據(jù)；

(b)加噪地震數(shù)據(jù)；

(c)非均勻采樣含噪數(shù)據(jù)；

(d)規(guī)則化數(shù)據(jù)；

(e)圖d與圖b數(shù)據(jù)的差

(a)常規(guī)曲波變換方法的去噪結(jié)果；

(b)常規(guī)曲波變換方法去除的噪聲；

(c)本文方法的去噪結(jié)果；

(d)本文方法去除的噪聲

為了詳細(xì)對(duì)比去噪效果，將原始理論地震數(shù)據(jù)、非均勻地震數(shù)據(jù)、含噪非均勻地震數(shù)據(jù)、常規(guī)曲波變換去噪結(jié)果和本文方法去噪結(jié)果分別進(jìn)行局部放大，如圖4所示?？梢钥吹椒蔷鶆虻卣饠?shù)據(jù)顯示在均勻網(wǎng)格上使有效波同相軸扭曲錯(cuò)動(dòng)明顯，再加上噪聲的影響，這種扭曲現(xiàn)象更為嚴(yán)重，使有效波能量和噪聲能量相互纏繞在一起，降低了整個(gè)地震數(shù)據(jù)的信噪比。從去噪結(jié)果的局部放大顯示來(lái)看，由于常規(guī)曲波變換的去噪前提條件是均勻采樣，所以在壓制非均勻地震數(shù)據(jù)噪聲時(shí)不能有效調(diào)整非均勻采樣點(diǎn)位置，使去噪后的有效波同相軸仍然扭曲錯(cuò)動(dòng)，并且噪聲壓制也不徹底。而本文方法去噪后的有效波與原始地震數(shù)據(jù)有效波非常接近，去噪后的同相軸光滑連續(xù)，這進(jìn)一步說(shuō)明本文方法在去除隨機(jī)噪聲的同時(shí)，也可以將非均勻地震數(shù)據(jù)調(diào)整為均勻采樣數(shù)據(jù)。

為了對(duì)比不同噪聲水平下本文方法的去噪效果，對(duì)原始理論數(shù)據(jù)加入標(biāo)準(zhǔn)差分別為0.075和0.100的高斯白噪聲，然后采用非均勻傅里葉反變換獲得非均勻地震數(shù)據(jù)，如圖5a和圖5b所示，信噪比分別為-6.99dB和-8.23dB?？梢钥闯鲇行Рㄐ盘?hào)被噪聲淹沒(méi)，難以識(shí)別，并且有效波同相軸局部扭曲，模糊不清，尤其是圖5b，幾乎看不出有效信號(hào)。采用本文方法同時(shí)進(jìn)行數(shù)據(jù)規(guī)則化和噪聲壓制，噪聲估計(jì)值分別為σ=3.8和σ=5.1，結(jié)果如圖5c和圖5d所示。去噪后的信噪比分別為3.34dB和1.99dB，可以看出本文方法能夠有效地去除非均勻采樣數(shù)據(jù)中的隨機(jī)噪聲，并且將其內(nèi)插為均勻采樣，極大地提高了原始含噪數(shù)據(jù)的信噪比，但隨著噪聲能量大幅度增強(qiáng)，有效波損傷也會(huì)增加。

(a)合成地震數(shù)據(jù)；

(b)無(wú)噪聲非均勻地震數(shù)據(jù)；

(c)含噪聲非均勻地震數(shù)據(jù)；

(d)常規(guī)曲波變換去噪結(jié)果；

(f)本文方法去噪結(jié)果

(a)標(biāo)準(zhǔn)差為0.075的非均勻含噪數(shù)據(jù)；

(b)標(biāo)準(zhǔn)差為0.100的非均勻含噪數(shù)據(jù)；

(c)圖5a的去噪結(jié)果；

(d)圖5b的去噪結(jié)果

5 應(yīng)用實(shí)例

圖6a為海上單炮地震數(shù)據(jù)，數(shù)據(jù)采集時(shí)設(shè)計(jì)道距為12.5m，然而由于海上拖纜的羽狀漂移導(dǎo)致采集道距不均勻，道距范圍為6～20m。圖6b為其局部放大，可以看出隨機(jī)噪聲較為發(fā)育，并且非均勻采樣扭曲了部分有效波同相軸。采用本文方法進(jìn)行去噪，分解的尺度為5，在第二個(gè)最粗尺度上的角度數(shù)為16，噪聲估計(jì)值σ=8.5。圖6c為本文方法去噪結(jié)果，圖6d為局部放大顯示，可以看出絕大部分噪聲能量得到了有效壓制，去噪后的地震數(shù)據(jù)同相軸較為連續(xù)，幾乎沒(méi)有損失有效波。圖7a和圖7b分別為原始單炮數(shù)據(jù)和本文方法去噪結(jié)果所對(duì)應(yīng)的二維頻譜，也可以看出本文方法去噪相對(duì)徹底，并且將非均勻采樣地震數(shù)據(jù)內(nèi)插為均勻采樣。從去噪前、后振幅譜差(圖7c)可以看出，去除了大部分隨機(jī)噪聲，有效波的損傷較小，大幅提高了信噪比。

(a)海上非均勻采樣道集；

(b)圖a局部放大顯示；

(c)本文方法去噪結(jié)果；

(d)圖c局部放大顯示

(a)原始數(shù)據(jù)；

(b)本文方法去噪結(jié)果；

(c)圖a與圖b數(shù)據(jù)的差

圖7 海上實(shí)際數(shù)據(jù)去噪前、后頻譜對(duì)比

6 結(jié)論與認(rèn)識(shí)

本文在多尺度多方向二維曲波變換基礎(chǔ)上，提出了基于非均勻二維曲波變換和線性Bregman方法的地震數(shù)據(jù)隨機(jī)噪聲壓制方法。由于非均勻曲波變換在正、反變換過(guò)程中具有能量無(wú)損性質(zhì)，而這種性質(zhì)特別適合在曲波域處理非均勻采樣下的含噪地震數(shù)據(jù)，并且將其內(nèi)插為均勻采樣，在此過(guò)程中能有效地壓制隨機(jī)噪聲。理論和實(shí)際資料的處理結(jié)果表明，該方法可以有效地壓制地震數(shù)據(jù)中的隨機(jī)噪聲，且盡可能地保護(hù)了微弱的有效波信號(hào)，使反射波同相軸更加連續(xù)、清晰。

本文方法采用非均勻采樣快速傅里葉變換，需要進(jìn)行褶積和反褶積運(yùn)算，計(jì)算速度遠(yuǎn)遠(yuǎn)低于均勻采樣快速傅里葉變換，而且曲波變換實(shí)現(xiàn)過(guò)程復(fù)雜、冗余度高，因此本文所提出的非均勻曲波變換去噪方法計(jì)算效率相對(duì)較低。本文只進(jìn)行了二維非均勻曲波變換，如推廣到三維數(shù)據(jù)，則計(jì)算時(shí)間會(huì)進(jìn)一步增加，需要發(fā)展其他快速算法提高本文方法的計(jì)算效率。