張 華 刁 塑 溫建亮 黃光南 朱雯婷 白 敏
(①東華理工大學放射性地質與勘探技術國防重點學科實驗室,江西南昌 330013;②山西省煤炭地質物探測繪院,山西晉中030600;③華北水利水電大學資源與環境學院,河南鄭州 450046)
由于表層激發條件的限制以及其他各種不可抗干擾因素的影響,僅在野外采集階段完全壓制隨機噪聲不太現實,需要在室內用有效的數學去噪方法提高疊前地震資料的信噪比,以滿足不同處理環節的需求[1-4]。目前,有許多行之有效的隨機噪聲壓制方法,其中稀疏變換應用較為廣泛,包括傅里葉變換[5-8]、小波變換[9-10]、Radon變換[11-12]、Contourlet變換[13]、曲波變換[14]以及Seislet 變換[15]等。該類方法根據隨機噪聲和有效波在稀疏變換域中能量分配的差異,對稀疏系數采用適當的閾值去除隨機噪聲的影響。為了取得較好的去噪效果,稀疏基需要盡可能地捕獲地震波的有效信息,并且少數較大的稀疏系數能夠代表信號的主要特征,而大部分較小、被濾除的系數不影響原始數據的主要特征[16]。曲波變換能夠有效地表示地震數據的尺度性和方向性,更加稀疏地表示地震波場局部細節特征,所以許多學者選用曲波變換方法去噪,且效果顯著[17-20]。
現有曲波變換去噪方法的前提是地震數據為均勻網格采樣,而對于非均勻網格采樣的含噪數據則效果不佳。由于野外地形條件及施工環境的制約,在很多情況下,地震數據常為非均勻采樣[21]。如果在去噪過程中將非均勻采樣視為均勻采樣進行處理,就得不到連續的地震波場,自然也不能充分壓制噪聲干擾,從而會影響到后續其他處理方法的應用。Hennenfent等[22]提出基于非均勻采樣的曲波變換重建方法,但未討論該方法在噪聲壓制中的應用。
本文在前人研究的基礎上,針對非均勻采樣地震數據,在傳統曲波變換的基礎上,引入非均勻傅里葉變換,建立均勻曲波系數與空間非均勻采樣地震道之間的規則化反演算子;然后選擇合適的噪聲估計值,使用線性Bregman方法反演,在迭代過程中采用軟閾值對曲波系數去噪,由反演得到無噪聲的均勻曲波系數;再進行常規曲波反變換,得到去噪后的地震數據,由此形成基于非均勻曲波變換和線性Bregman方法的隨機噪聲壓制方法。模型數據和實際數據的實驗結果表明了本文方法的有效性。
Candès等[23]提出了第二代曲波變換,從而使曲波變換更容易被理解、運算效率更高、實現更簡單,能夠為地震信號提供最優的稀疏表示方式。實際上,實現快速離散曲波變換主要包括兩個步驟:①對地震數據應用二維傅里葉變換,得到頻率—波數域系數;②在頻率—波數域形成角度楔形,將每一個楔形圍繞到原點重新裝配,并對每一個裝配好的楔形應用二維傅里葉反變換,得到離散曲波系數。參考文獻[22]定義曲波正變換算子為
A=TF
(1)
式中:F實現了上述離散曲波變換第①步:T則實現了第②步。
定義曲波反變換算子為
AH=FHTH
(2)
式中:上標“H”表示共軛轉置;FH表示二維傅里葉反變換,將頻率—波數域轉換到時間—空間域中;TH表示曲波平鋪算子,即將曲波系數變換到頻率—波數域的過程。


(3)


(4)
該算子描述了非均勻采樣下離散曲波系數與地震數據之間的關系。
由于快速離散曲波變換具有緊支撐性,正、反曲波變換前、后能量無損失,即
AAH=I
(5)
式中I為單位矩陣。對于新的非均勻曲波反變換算子C,由于包含非均勻快速傅里葉變換,那么算子C就不滿足式(5)。為此定義非均勻曲波正變換算子C+為
滿足 ‖d-Cx‖2≤σ
(6)

式(6)最優化問題可轉為求解下述基追蹤規則化問題[25]
(7)
式中λ是一個閾值權衡因子,平衡L1范數和L2范數的比重,可以在第一次迭代過程計算。線性Bregman方法求解迭代式為
(8)
式中:zk表示第k次迭代得到的曲波系數向量;lk為動態步長,定義為
(9)
軟閾值函數為
Sλ(x)=sign(x)max(|x|-λ,0)
(10)

(11)

(12)
該式可以保證在每次迭代過程中對含噪地震數據進行噪聲壓制,從中也可以看出線性Bregman方法實現非常簡單,沒有過多的調節參數。
一個含噪聲二維地震信號模型可以表示成為
d(i,j)=s(i,j)+e(i,j)
(13)
式中:s為真實信號;e為噪聲。去噪過程就是從非均勻含噪地震信號d中提取真實信號s,去除噪聲干擾信號e(圖1)。

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

(14)
式中AH是常規的曲波反變換算子。通過對反演后的無噪聲均勻曲波系數應用常規曲波反變換,可得到均勻采樣后的去噪地震數據。
信噪比定義為[29]
(15)
式中:s0表示原始均勻不含噪聲數據;s表示去噪后的均勻地震數據。R越高,表示去噪效果越好。
圖2a為合成的256道地震數據,有4個反射波同相軸,各反射同相軸能量有差異,采樣間隔為1ms,道間距為5m,每道1024個采樣點。加入標準差為0.05的高斯白噪聲,如圖2b所示,其信噪比為-5.11dB。對含噪數據進行空間均勻傅里葉變換,然后再進行空間非均勻傅里葉反變換,得到新的空間非均勻采樣下的256道地震數據,如圖2c所示,信噪比為-5.25dB,道間距范圍為0~10m。如果直接用均勻曲波變換法對圖2c數據去噪,則扭曲的同相軸得不到校正,導致去噪后的地震記錄誤差較大。因此首先采用非均勻曲波變換方法對圖2c進行規則化處理,此時噪聲估計值σ=0,曲波變換所選擇的尺度數為5,在第二個最粗尺度上的角度數為32。規則化結果如圖2d所示,可以看見規則化后的地震波場非常連續,幾乎沒有視覺上的差異,并且從誤差剖面(圖2e)也可以看出規則化前后幾乎沒有誤差,說明非均勻曲波變換方法規則化效果好。
采用均勻曲波變換方法對非均勻含噪數據(圖2c)直接去噪,對分解后的曲波系數進行軟閾值處理,將保留下來的曲波系數進行曲波反變換,從而得到最終的去噪結果,如圖3a所示,其信噪比為4.45dB。圖3b為去除的噪聲干擾,可以看出,其中含有部分有效波信號,主要原因是常規曲波變換不能直接處理非均勻采樣數據,而是把非均勻采樣數據當成均勻采樣數據進行處理。根據噪聲能量計算出噪聲估計值σ=2.5,采用本文方法在非均勻采樣數據規則化過程中同時去噪,結果如圖3c所示,其信噪比為5.11dB。去除的噪聲如圖3d所示。可以看出,本文方法不僅可以將非均勻采樣地震數據內插為均勻采樣數據,還可以同時有效地去除隨機噪聲干擾,去噪后的地震同相軸更加連續、清晰,大幅度地提高了信噪比,且去除的噪聲中幾乎不含有效波能量,表明對于非均勻采樣地震數據本文方法具有較好的去噪效果。

(a)合成理論地震數據;

(b)加噪地震數據;

(c)非均勻采樣含噪數據;

(d)規則化數據;

(e)圖d與圖b數據的差

(a)常規曲波變換方法的去噪結果;

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

(c)本文方法的去噪結果;

(d)本文方法去除的噪聲
為了詳細對比去噪效果,將原始理論地震數據、非均勻地震數據、含噪非均勻地震數據、常規曲波變換去噪結果和本文方法去噪結果分別進行局部放大,如圖4所示。可以看到非均勻地震數據顯示在均勻網格上使有效波同相軸扭曲錯動明顯,再加上噪聲的影響,這種扭曲現象更為嚴重,使有效波能量和噪聲能量相互纏繞在一起,降低了整個地震數據的信噪比。從去噪結果的局部放大顯示來看,由于常規曲波變換的去噪前提條件是均勻采樣,所以在壓制非均勻地震數據噪聲時不能有效調整非均勻采樣點位置,使去噪后的有效波同相軸仍然扭曲錯動,并且噪聲壓制也不徹底。而本文方法去噪后的有效波與原始地震數據有效波非常接近,去噪后的同相軸光滑連續,這進一步說明本文方法在去除隨機噪聲的同時,也可以將非均勻地震數據調整為均勻采樣數據。
為了對比不同噪聲水平下本文方法的去噪效果,對原始理論數據加入標準差分別為0.075和0.100的高斯白噪聲,然后采用非均勻傅里葉反變換獲得非均勻地震數據,如圖5a和圖5b所示,信噪比分別為-6.99dB和-8.23dB。可以看出有效波信號被噪聲淹沒,難以識別,并且有效波同相軸局部扭曲,模糊不清,尤其是圖5b,幾乎看不出有效信號。采用本文方法同時進行數據規則化和噪聲壓制,噪聲估計值分別為σ=3.8和σ=5.1,結果如圖5c和圖5d所示。去噪后的信噪比分別為3.34dB和1.99dB,可以看出本文方法能夠有效地去除非均勻采樣數據中的隨機噪聲,并且將其內插為均勻采樣,極大地提高了原始含噪數據的信噪比,但隨著噪聲能量大幅度增強,有效波損傷也會增加。

(a)合成地震數據;

(b)無噪聲非均勻地震數據;

(c)含噪聲非均勻地震數據;

(d)常規曲波變換去噪結果;

(f)本文方法去噪結果

(a)標準差為0.075的非均勻含噪數據;

(b)標準差為0.100的非均勻含噪數據;

(c)圖5a的去噪結果;

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

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

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

(c)本文方法去噪結果;

(d)圖c局部放大顯示

(a)原始數據;

(b)本文方法去噪結果;

(c)圖a與圖b數據的差
圖7 海上實際數據去噪前、后頻譜對比
本文在多尺度多方向二維曲波變換基礎上,提出了基于非均勻二維曲波變換和線性Bregman方法的地震數據隨機噪聲壓制方法。由于非均勻曲波變換在正、反變換過程中具有能量無損性質,而這種性質特別適合在曲波域處理非均勻采樣下的含噪地震數據,并且將其內插為均勻采樣,在此過程中能有效地壓制隨機噪聲。理論和實際資料的處理結果表明,該方法可以有效地壓制地震數據中的隨機噪聲,且盡可能地保護了微弱的有效波信號,使反射波同相軸更加連續、清晰。
本文方法采用非均勻采樣快速傅里葉變換,需要進行褶積和反褶積運算,計算速度遠遠低于均勻采樣快速傅里葉變換,而且曲波變換實現過程復雜、冗余度高,因此本文所提出的非均勻曲波變換去噪方法計算效率相對較低。本文只進行了二維非均勻曲波變換,如推廣到三維數據,則計算時間會進一步增加,需要發展其他快速算法提高本文方法的計算效率。