基于曲波變換的面波去噪方法研究

姚恒星，謝凱 （長江大學電子信息學院，湖北 荊州434023）

在地震勘探領域中，噪聲嚴重影響了地震資料的后續處理，因此如何有效地提高地震資料信噪比是地震資料處理的首要任務。面波［1］作為一種常見的規則干擾，在疊前數據上呈掃帚形狀，具有低頻能量強等特點，嚴重影響了地震有效反射波，而且還降低了地震數據的信噪比。常見的去除面波的方法有低通或高通濾波、F－K頻率波數譜濾波［2］等。每種方法有著一定的效果，但這些方法都只側重考慮了面波單一的特性，有著較大的局限性。二維小波變換去噪［3］雖然可以壓制噪聲，但它是以點為元素來描述信號圖像特征，無法有效地表達邊緣信息，并且還會損傷有效信號。為此Candès等提出了曲波變換［4，5］（Curvelet），即基于脊波變換的多尺度多分辨率的幾何分析方法［6］。該變換稱為后小波變換，不僅彌補了小波變換的不足，還可以更加有效的表示二維圖像，具有各向異性、方向性和局部性，可以稀疏表達圖像的平滑區域和邊緣區域［7，8］。把地震數據看成是二維圖像，可以利用面波的方向性［9，10］，用曲波變換去噪的方法去除面波的同時能夠有效地保護有效反射信號。為此，筆者提出基于曲波變換的面波去噪方法研究。

1 曲波變換的基本原理

第1代曲波變換的數字實現很復雜，需要子帶分解、正規化、平滑分塊和脊波分析等許多步驟，并且在進行曲波金字塔分解時帶來了巨大的數據冗余。因此，Strack等在第1代曲波變換的基礎上提出了第2代曲波變換的新框架體系。第2代曲波變換在構造上已經完全不同于第1代曲波變換，實現過程不需要用到脊波變換，直接通過構造曲波窗函數來實現曲波分解，不僅變換意義明確，而且實現起來更加快速和方便。

曲波變換和小波變換屬于稀疏理論的范疇，都是利用基函數與信號作內積來實現信號的稀疏表示。曲波變換可表示為：

式中，φj，l，k表示曲波函數，j，l，k分別表示尺度方向和位置參數。

下面介紹第2代曲波變換的基本原理。曲波變換在頻域內的實現是采用頻域中的窗函數來表示的，首先構造徑向窗和角度窗W（r），r∈和V（t），t∈［－1，1］，那么對所有尺度j，傅里葉頻率窗Uj定義為：

式中，［］表示的整數部分；Uj是受W和V支撐區間的限制而獲得的楔形區域，如圖1所示的陰影區域。

定義在尺度j，方向θl，位置參數k＝（k1，k2）處的連續曲波變換為：

式中，（ω）是二維有效信號的頻域表示表示對φj，l，k（ω）取共軛；Rθl由θl旋轉獲得；＝·2－j，k2·2－j2）。

和小波基礎理論一樣，曲波變換包括粗尺度和精細尺度。粗尺度的曲波變換不具有方向性，所以整個曲波變換是由粗尺度下各向同性的小波和精細尺度下的方向元素組成的。

把笛卡爾坐標系下的二維函數f（t）作為有效信號，離散曲波變換定義：

式中，cD（j，l，k）是曲波變換系數的離散形式（t）是曲波函數的離散形式。

第2代曲波變換的快速算法有2種：USFFT算法和Wrapping算法。筆者采用基于Wrapping的快速離散曲波算法，其核心思想是圍繞原點Wrapping，對任意區域周期化，再一一映射到原點仿射區域。過程如下：

1）對給定的二維函數f［t1，t2］進行二維 FFT，得到頻域表示［n1，n2］，－≤n1，n2≤。

2）在頻域，對每個尺度和角度組（j，l），重采樣［n1，n2］得到［n1，n2－n1tanθl］，（n1，n2）∈Pj。

3）內插后的與窗函數相乘得到［n1，n2］＝［n1，n2－n1tanθl］＾Uj［n1，n2］。

4）圍繞原點 Wrapping局部化［n1，n2］。

5）對，l進行二維IFFT（FFT的逆變換）得到離散曲波系數集合CD（j，l，k）。

Wrapping算法基本思路如下：首先變換到頻域，再在頻域中局部化，最后采用二維IFFT得到曲波系數。此外局部化和二維IFFT可合二為一，即用局部化窗口來乘局部傅里葉基。該算法采用3個參數，使得理解更容易，運算操作更簡便，冗余度更低。

圖1 曲波變換的時域和頻域

2 面波去噪方法

對于地震數據的噪聲去除問題，通常采用的方法是通過變換將信號去噪問題從時域轉換到頻域加以解決。基于曲波變換的面波去噪方法如下：首先對經過預處理的含噪信號進行多尺度分解，然后在各尺度下盡可能提取出有效信號的曲波系數，同時去除噪聲的曲波系數，最后用曲波逆變換重構出地震信號，從而達到去噪目的。

設地震信號為s，有效信號為d，噪聲為n，則含噪聲地震數據可表示為：

有效信號可用下述方法估算，即：

式中，C表示曲波變換；C－1表示曲波逆變換；Cs表示對地震信號s作曲波變換后的系數；F表示閾值函數，定義為：

式中，m為大小與尺度有關的常數；σc為曲波域中噪聲標準差；σ為噪聲標準差。

該方法所述的地震信號去除面波方法的流程包括4個步驟：

1）通過對大量含面波的實際疊前地震資料進行分析，確定地震信號中包含的噪聲模型及其相應的各項參數，從許多具有代表性的面波里取其平均作為面波噪聲模型。

2）因為地震信號是二維的，在空間域中分析存在很多局限性。由于曲波變換具有多尺度和多方向性的特點，能夠稀疏地表示二維信號，使得在曲波域中能更精細地分離出面波噪聲和有用地震地震反射信號。該方法對二維地震數據進行曲波分解，選擇合適的尺度，將空域地震信號變換到曲波域中，得到地震信號的各方向各尺度的曲波系數。同時對面波模型也進行曲波分解，得到面波曲波系數，將面波數據變換到曲波域。

3）根據面波噪聲模型的曲波系數在不同尺度和方向上的分布特點來設置去除面波的閾值。接著采用閾值去噪法，低于閾值的系數可以認為是面波的曲波系數，從而將其置零去掉；大于閾值的系數認為是有效反射信號的曲波系數，將其保留。

4）對濾掉面波后的曲波系數進行曲波反變換，重構得到去噪后的地震信號。

3 試驗結果與分析

筆者分別對合成地震資料和疊后實際資料進行了相應處理。其測試結果如圖2所示。

圖2（a）是理論合成數據，該模型有2條主頻為35Hz和45Hz的反射同相軸和1條主頻為15Hz的面波同相軸，2條斜線代表面波，2條曲線代表有效反射波。該模型數據由100道組成，每道包括600個采樣點。圖2（b）是采用筆者所描述的方法去噪后結果，圖2（c）是去除掉的面波。由圖2（c）可見去除的大部分是面波。

圖2 模型地震數據的去噪結果對比

圖3（a）是含面波實際地震數據，共33道，每道1501個采樣點。從圖3（a）中可以看出，原始實際地震數據受面波干擾嚴重，許多有效波同相軸無法識別。筆者設計了包含面波的噪聲模型，將地震數據和面波噪聲模型都進行曲波分解，得到其在曲波域中的表示。根據面波噪聲模型在曲波域中分布的尺度和方向，確定了濾波閾值，然后根據此閾值對地震信號的曲波系數進行處理。圖3（b）是濾波后的重構圖像顯示效果，可以看出有效波的水平同相軸變得更加清晰，部分在處理前無法識別的反射層在處理后顯現出來。圖3（c）是通過該方法處理所去除掉的面波噪聲。

表1是實際地震數據去除面波前后的信噪比對比，抽取的是地震數據的前10道數據。由表1中可以看出在經過曲波變換去面波處理后，地震數據的信噪比有了明顯的提升，由此可見曲波去面波方法具有較強的實用性。

4 結語

通過以上試驗可知，將該方法運用于地震信號的去噪處理，能夠充分利用面波在曲波域中的分布特點，有效地壓制面波干擾并保護反射地震信號，重構后的地震有效反射信號同相軸變得更加清晰，不僅提高了信噪比，還提高了成像質量。

表1 實際地震數據去噪信噪比對比

圖3 實際地震數據的去噪結果對比

［1］李晶 .面波在地震波場中的特性研究及其應用 ［D］.成都：成都理工大學，2006.

［2］李彩芹，張華 .小波變換與F－K聯合濾波在面波分離中的應用 ［J］.中國煤田地質，2007，19（4）：60～61，84.

［3］林椹尠，宋國鄉，薛文 .圖像的幾種小波去噪方法的比較與改進 ［J］.西安電子科技大學學報（自然科學版），2004，31（4）：626～629.

［4］董烈乾，李振春，王德營，等 .第2代Curvelet變換壓制面波方法 ［J］.石油地球物理勘探，2011，46（6）：897～904.

［5］蔡炳煌 .基于曲波分析的圖像處理與應用 ［D］.汕頭：汕頭大學，2007.

［6］才溪 .多尺度圖像融合理論與方法 ［M］.北京：電子工業出版社，2014.

［7］Starck J L，Candes E，Donoho D L.The Curvelet Transform for Image Denoising ［J］.IEEE Transactions on Image Processing，2002，11（6）：670～684.

［8］Kristof De Meersman.Ground Roll polarization filtering with spatial smoothness constraints ［J］.SEG Technical Program Expanded Abstracts，2008（27）：413～416.

［9］代虎 .地震面波資料處理的基本方法 ［J］.黑龍江水利科技，2012，41（9）：28～30.

［10］張恒磊，劉天佑.Curvelet域蒙特卡羅估計的噪聲衰減 ［J］.西南石油大學學報（自然科學版），2011，33（4）：64～68.