余昌輝

（中石化海洋石油工程有限公司上海物探分公司，上海 201208）

基于曲波變換的地震隨機噪聲壓制新方法

余昌輝

（中石化海洋石油工程有限公司上海物探分公司，上海 201208）

此文提出了一種通過對曲波系數(shù)進行處理來實現(xiàn)地震隨機噪聲壓制的新方法。地震信號經(jīng)曲波變換得到的曲波系數(shù)在各尺度各方向的分布仍具有同相軸特征，但噪聲部分的曲波系數(shù)的分布是隨機的，故可將曲波系數(shù)作為新的目標(biāo)信號，對其進行“去噪”，最后用經(jīng)過處理的曲波系數(shù)還原地震數(shù)據(jù)實現(xiàn)對隨機噪聲的壓制。數(shù)值算例驗證了該方法的正確性和有效性，并與用離散余弦變換和小波變換兩種方法處理曲波系數(shù)的結(jié)果進行對比。實驗表明用曲波變換處理曲波系數(shù)的去噪效果要優(yōu)于后兩者，且比單純進行一次曲波變換去噪的結(jié)果在信噪比上有明顯提高。

曲波變換；曲波系數(shù)；多尺度多方向；隨機噪聲

噪聲壓制是地震勘探數(shù)據(jù)處理的主要任務(wù)之一。根據(jù)噪聲與有效信號的相干性可將其分為隨機噪聲和相干噪聲，本文主要針對隨機噪聲進行壓制。曲波變換是目前備受關(guān)注的去噪方法之一，該方法由Candes和Donoho于1999年以一個數(shù)學(xué)概念提出，2004年首次被地球物理學(xué)家引入到地震數(shù)據(jù)去噪[1-3]。曲波變換對二階可導(dǎo)函數(shù)具有最優(yōu)逼近階，且具有多方向、多尺度和各向異性的特點，這使得該變換很適合結(jié)構(gòu)有很強幾何性、正則性并表現(xiàn)為二階可導(dǎo)曲線特征的地震信號；另外，曲波變換對地震數(shù)據(jù)的表達具有稀疏性，即用極少的曲波系數(shù)就可以表示原信號中的有效信號。本文在前人工作的基礎(chǔ)上，通過對曲波域內(nèi)的曲波系數(shù)結(jié)構(gòu)及分布進行深入研究，發(fā)現(xiàn)有效信號的曲波系數(shù)仍具有地震同相軸特征，據(jù)此提出了一種新的噪聲壓制方法。

1 方法原理

1.1 曲波變換基本原理

曲波變換的定義是將均方可積函數(shù)f映射為系數(shù)序列c ( j, l, k )的變換，其中j表示尺度參數(shù)，l表示方向參數(shù)，k表示位置參數(shù)。在快速離散曲波變換中，對于尺度為j，空間位置為b = ( k1·2-j, k2·2-[j/2])，直角坐標(biāo)下的曲波函數(shù)為：

曲波變換最核心的關(guān)系是曲波基的支撐區(qū)間滿足：width∝length2，該關(guān)系稱為各向異性尺度關(guān)系，它表明了曲波是一種具有方向性的基原子，這也是曲波變換具有很好的非線性逼近能力的根本原因。通過信號與某一曲波函數(shù)內(nèi)積的形式能實現(xiàn)對信號的稀疏表達[4-8]：

1.2 基于曲波變換的地震資料去噪

利用曲波的稀疏性，我們可將地震數(shù)據(jù)中隨機噪聲的衰減表述為如下的約束優(yōu)化問題[9]：式中：y為帶有噪聲的地震數(shù)據(jù)；A=CT表示反曲波變換；x表示曲波系數(shù)；ε為一任意小量；為最終輸出結(jié)果。

由于該問題的求解較困難，用如下優(yōu)化問題代替了上述約束優(yōu)化問題：

這里加入了一個起到正則化作用的閾值算子λ。問題Pλ是通過迭代閾值的方法求解的。其迭代序列為：

1.3 新方法的原理及實現(xiàn)

將含噪的地震數(shù)據(jù)做曲波變換，得到各尺度各方向上的曲波系數(shù)，發(fā)現(xiàn)曲波系數(shù)在曲波域各區(qū)塊的分布仍具有與地震信號相似的同相軸特征，即地震信號部分的曲波系數(shù)仍呈地震信號特征，而隨機噪聲部分的曲波系數(shù)的分布仍是隨機的，基于此特征我們提出了將各個區(qū)塊內(nèi)的曲波系數(shù)作為處理目標(biāo)的去噪新方法，即對各區(qū)塊內(nèi)曲波系數(shù)分別進行曲波變換并施以閾值，以此實現(xiàn)對曲波系數(shù)的“去噪”。

根據(jù)曲波變換原理，若f為含噪地震信號，進行一次曲波去噪可得到去噪結(jié)果f'：

式中：c' ( j, l, k)為在曲波域經(jīng)過處理的曲波系數(shù)；φ'j,l,k為反曲波函數(shù)。

由（6）式有，新方法去噪的結(jié)果f"可表示為：

式中：c?（j,l,k）為經(jīng)過曲波去噪處理后的c″（j,l,k）；c′（j,l,k）為經(jīng)一次曲波變換去噪后的數(shù)據(jù)

經(jīng)過再次曲波變換后分解到各尺度各方向上的曲波系數(shù)；根據(jù)曲波變換的原理有：

所以有

將式（6）、式（8）、式（9）代入式（7）得新方法去噪結(jié)果 ：

該方法的核心步驟是將一次曲波變換的系數(shù)分解到不同尺度方向的小塊區(qū)域，然后分別對這些區(qū)域系數(shù)再進行曲波變換。分解時的尺度數(shù)和方向數(shù)由被分解數(shù)據(jù)大小決定，圖1a中地震記錄水平采樣點為256個，故最合適的尺度數(shù)為5，基礎(chǔ)方向數(shù)為8。經(jīng)計算第一尺度和最高尺度的方向數(shù)均為1，第二、三尺度方向數(shù)為8，第四尺度方向數(shù)為16。曲波域內(nèi)系數(shù)分布情況如圖1b所示。

圖1 地震數(shù)據(jù)所對應(yīng)曲波系數(shù)在曲波域內(nèi)分布

新方法的實現(xiàn)如下：

對于一地震信號s，設(shè)有效信號為d，n為隨機噪聲，則去噪模型可表示為：

（1）根據(jù)數(shù)據(jù)大小確定尺度數(shù)和方向數(shù)，對s進行曲波變換，對得到的曲波系數(shù)c ( j, l, k )加閾值，再做反曲波變換，得去噪結(jié)果s'；

（2）對s'進行曲波變換，將得到的曲波系數(shù)分解到各尺度各方向，把每一塊曲波系數(shù)作為新的去噪對象；

（3）根據(jù)每個小塊曲波系數(shù)的分布情況確定第二次曲波變換的尺度數(shù)和方向數(shù)，分別對每一塊曲波系數(shù)進行曲波變換，得到二次變換后的曲波系數(shù)c" ( j, l, k )；

（4）對c" ( j, l, k )進行閾值處理，并進行反曲波變換，得到每一塊的去噪后的曲波系數(shù)，將這些曲波系數(shù)合成為新的c'" ( j, l, k )；

（5）對c'" ( j, l, k )做反曲波變換，得到最終的去噪結(jié)果s"。

圖2 合成數(shù)據(jù)實驗

2 數(shù)值實驗

2.1 陡傾構(gòu)造模型

為了驗證新方法的正確性和有效性，首先在合成的共炮點道集數(shù)據(jù)上做數(shù)值實驗。速度模型為陡傾構(gòu)造模型，如圖2a所示。水平網(wǎng)格點數(shù)為3 375，網(wǎng)格間隔為10，深度網(wǎng)格點數(shù)為570，網(wǎng)格間隔為10；記錄共256道，道間距為50 m，時間采樣點數(shù)為1 025，采樣間隔為0.001 ms。共炮點道集如圖2b所示。加入幅值與有效信號最大值在一個數(shù)量級上的范圍在-0.1~0.1之間的隨機噪聲，得到含噪剖面（圖2c），PSNR=50.32 dB。這里用PSNR（峰值信噪比）來衡量去噪水平。其計算公式為：

PSNR=10 log10（max2/MSE） （12）式中：max為數(shù)據(jù)的峰值；MSE為原數(shù)據(jù)與處理數(shù)據(jù)之間的均方誤差。

曲波變換去噪的結(jié)果如圖3a所示，新方法得到的去噪結(jié)果如圖3b所示。為方便對比去噪的效果，分別將原始數(shù)據(jù)、含噪數(shù)據(jù)、曲波變換去噪結(jié)果和新方法去噪結(jié)果的剖面取出相同位置的一部分進行放大，如圖3c、3d、3e、3f所示。計算得常規(guī)曲波變換去噪結(jié)果，PSNR=55.992 dB，新方法去噪結(jié)果，PSNR=61.164 dB。從圖3中可看出，新方法的去噪結(jié)果在較好地保存了地震波同向軸的同時，進一步去掉了噪聲，其結(jié)果優(yōu)于曲波變換的去噪結(jié)果。

圖4可見，經(jīng)處理的曲波系數(shù)得到了明顯改善，對應(yīng)有效信號的曲波系數(shù)留了下來，曲波系數(shù)的稀疏性增強。用這些新的曲波系數(shù)組合起來進行反曲波變換得到的結(jié)果，噪聲必然會受到抑制，因此達到了改進去噪效果的目的。

2.2 對比實驗

為了驗證本文提出的新方法（即用曲波變換處理曲波系數(shù)）的優(yōu)越性，進行了與用小波變換和離散余弦變換處理曲波系數(shù)的對比實驗。其結(jié)果如圖5和表1所示。

由圖5和表1可知，三種方法對隨機噪聲都起到一定的壓制，且均方誤差都較小，其中曲波變換得到的結(jié)果均方誤差最小，信噪比最高；從去掉噪聲多少的角度看，小波變換和曲波變換所得結(jié)果優(yōu)于離散余弦變換，但從有效信號的形態(tài)上講，用離散余弦變換和曲波變換的處理結(jié)果要優(yōu)于小波變換，前兩者更好地保留了有效信號的信息，三者相比，后者在有效信號形態(tài)上失真程度最高。綜合以上因素曲波變換處理結(jié)果最優(yōu)。

圖3 新方法去噪結(jié)果與常規(guī)曲波變換去噪結(jié)果對比

圖4 第三尺度、第三方向的曲波系數(shù)去噪前后對比情況

3 結(jié)論

本文利用了曲波變換的多尺度多方向性，及其對地震信號稀疏表達的特性，通過對曲波系數(shù)分布規(guī)律的研究，發(fā)現(xiàn)有效信號在曲波域內(nèi)的系數(shù)也具有同相軸特征。基于現(xiàn)有的曲波變換去噪方法提出了以曲波系數(shù)為去噪目標(biāo)的新方法，該方法是對傳統(tǒng)方法的改進。經(jīng)實驗證明新方法正確且有效。與傳統(tǒng)曲波變換去噪結(jié)果對比顯示，該方法在提高信噪比上有顯著作用；另外，通過與小波變換和離散余弦變換的對比實驗表明曲波變換無論是對信噪比的提高還是對有效信號的保留都優(yōu)于二者，印證了曲波變換對地震信號的最佳適應(yīng)性，從而證明本文的方法在地震數(shù)據(jù)處理中具有一定價值。

圖5 不同方法處理曲波系數(shù)去噪結(jié)果

表1 不同方法處理曲波系數(shù)去噪結(jié)果量化對比

[1]胡天躍. 地震資料疊前去噪技術(shù)的現(xiàn)狀與未來[J]. 地球物理學(xué)進展， 2002， 17（2）： 218-223.

[2]CANDèS E J， DONOHO D L. Curvelets-a Surprisingly Effective Nonadaptive Representation for Objects with Edges[R]. Technical Report. Stanford： Department of Statistics， Stanford University，1999.

[3]MA J W， PLONKA G. A Review of Curvelets and Recent Applications[J]. IEEE Signal Processing Magazine， 2010， 27（2）： 118-133.

[4]包乾宗， 陳文超， 高靜懷. 基于第二代Curvelet變換的地震資料隨機噪聲衰減[J]. 煤田地質(zhì)與勘探， 2010， 38（1）： 66-70.

[5]仝中飛. Curvelet閾值迭代法在地震數(shù)據(jù)去噪和插值中的應(yīng)用研究[D]. 長春： 吉林大學(xué)， 2009.

[6]韓佳君. Curvelet變換組合法壓制地震隨機噪聲研究[D]. 長春： 吉林大學(xué)， 2010.

[7]薛念. 基于Curvelet變換的地震數(shù)據(jù)插值和去噪[D]. 成都：西南交通大學(xué)， 2010.

[8]吳愛弟， 趙秀玲. 基于曲波變換的地震信號去噪方法[J]. 中國石油大學(xué)學(xué)報（自然科學(xué)版）， 2010， 34（3）： 30-33.

[9]WANG D L， TONG Z F， TANG C， et al. An Iterative Curvelet Thresholding Algorithm for Seismic Random Noise Attenuation[J]. Applied Geophysics， 2010， 7（4）： 315-324.

A New Method to Suppress Random Seismic Noise Based on Curvelet Transform

YU Changhui
（Shanghai Geophysical Division of SINOPEC Offshore Oilfield Engineering Company, Shanghai 201208, China）

This paper proposes a new method to suppress the random seismic noise by disposing the curvelet transform (CVT) coeffcient. The curvelet coeffcients of seismic signal after CVT still remain the characteristics of seismic event on different scales and along different directions, but the curvelet coeffcients of noise distribute randomly. The curvelet coeffcient, therefore, can be taken as the new target signal and de-noised. Then the processed CVT coeffcients are used to restore the seismic data and to suppress the random seismic noises. Numerical experiment verifed the correctness and validity of this method. Compared with the result obtained by discrete cosine transform and wavelet transform, the de-noising effect of CVT is better, and also better than that of only once CVT in signal to noise ratio.

CVT; curvelet coeffcient; multiple scale&multiple direction; random noise

P631.4

A

10.3969/j.issn.1008-2336.2017.01.011

1008-2336（2017）01-0011-05

2016-03-25；改回日期：2016-07-05

余昌輝，男，1988年生，助理工程師，理學(xué)學(xué)士，地球物理專業(yè)。E-mail: yuchanghui@sopgc.com。