基于低信噪比條件下新型Seislet變換的閾值去噪方法

劉 財，崔芳姿，劉 洋，王 典，劉殿秘，張 鵬

1.吉林大學地球探測科學與技術學院，長春 1300262.中國石油吉林油田公司地球物理勘探研究院，吉林 松原 138000

?

基于低信噪比條件下新型Seislet變換的閾值去噪方法

劉 財1，崔芳姿1，劉 洋1，王 典1，劉殿秘2，張 鵬1

1.吉林大學地球探測科學與技術學院，長春 130026
2.中國石油吉林油田公司地球物理勘探研究院，吉林 松原 138000

Seislet變換是一種類小波變換方法，主要根據小波基沿地震同相軸的局部傾角方向來分析數據，其中，局部地震傾角的表征是該方法的核心。局部傾角的求取方法有很多種，但是往往在低信噪比的條件下存在著一些局限。根據共中心點道集中基于時距關系的地震傾角定義，提出一種適應于低信噪比條件下的傾角求取方法。對比基于時距關系與平面波分解濾波器計算出的局部地震傾角，結果證明，該方法能更加準確地表征低信噪比條件下同相軸的傾角信息。將基于時距關系的局部地震傾角用于Seislet框架，建立表征低信噪比數據的新型Seislet變換方法。在地震數據處理中，引入語音信號中改進的閾值方法。結合新型Seislet變換，提出閾值去噪方法，此方法不但適于地震數據，而且在提高信噪比方面也優于傳統的閾值去噪方法。實際數據處理的結果驗證了新型Seislet變換與改進閾值去噪方法的組合能夠有效地解決低信噪比條件下的信號提取任務。

Seislet變換；低信噪比；局部地震傾角；閾值；去噪

0 引言

小波去噪是一類有效分離信號與隨機噪聲的方法。在地震隨機噪聲消除中:張波等[1]利用小波分析和反演方法聯合去除隨機噪聲；王典等[2]基于提升算法和百分位數軟閾值消除地震隨機噪聲，并取得明顯的效果；樊計昌等[3]提出了小波包節點域和空間域傾角掃描高階相關去噪技術；Zhang等[4]基于非線性小波閾值去噪和遞歸最優閾值除去地震隨機噪聲，保持了處理信號的平滑度和完整性；王志強等[5]利用檢波器組合響應來提高地震資料的信噪比，壓制隨機噪聲。然而，小波變換對非平穩信號的表征存在很大的局限性。

類小波變換在地震圖像處理和數據分析中具有很廣泛的應用[6-8]，其主要利用數據的方向特性[9-11]。Fomel等[12-13]開發了一種針對地震數據方向性的數學變換方法，命名為“Seislet變換”，Seislet變換的定義結合了小波提升算法和局部平面波分解技術。一般來說，沿測線方向的離散小波變換是一種特殊的零傾角Seislet變換，Seislet變換能夠更有效地壓縮地震數據[14]。其中，局部地震傾角的求取是該方法的核心。局部傾角的計算有多種方法。相干分析中的歸一化互相干、任意多道相干、三維地震數據體相干算法都可以用來計算地震同相軸的局部傾角[15]。平面波分解濾波器(PWD)是一種特殊的預測誤差濾波器，通過局部平面波方程的表征，使分解殘差達到最小來確定局部傾角[16]。也可以通過希爾伯特變換用瞬時空間頻率與瞬時時間頻率的比值來獲取局部傾角[17]。利用曲波變換在曲波域可以實現同時對多個傾角的估計[18]。然而，在低信噪比的條件下，上述方法存在一些局限。

Donoho等[19]于1994年提出了信號去噪的軟閾值方法和硬閾值方法，這種方法廣泛應用于各個領域。然而，硬閾值函數不連續，軟閾值函數不可導。為了克服硬閾值函數和軟閾值函數的缺點，很多學者對這一問題進行了研究。Gao等[20]提出了半軟閾值方法進行改進；李沖泥等[21]于1999年提出一種新的閾值函數，該函數連續可導，彌補了原函數的分段性和不連續性；許麗群[22]于2010年提出了軟硬閾值折中法，既克服了硬閾值函數不連續的缺點，同時又克服了軟閾值函數中估計小波系數與分解小波系數之間存在著恒定偏差的缺陷。

在本次研究中，筆者從共中心點道集的反射時距曲線出發，建立了一種基于地震波運動方程且適應于低信噪比條件下傾角求取的方法，并將求出的局部地震傾角用于構建新型Seislet變換，將語音信號中改進的閾值方法引入地震數據處理中，通過理論模型和實際數據的處理結果驗證本次研究方法的有效性。

1 理論基礎

1.1 基于時距關系地震傾角的新型Seislet變換

Seislet變換是把小波提升算法同局部平面波分解技術相結合(見附錄A)，改進了小波提升算法中的預測算子P(公式(A-4))和更新算子U(公式(A-7))，使之能夠沿著同相軸的方向進行預測和更新。

以往Seislet變換在低信噪比的條件下存在著一些局限。在本文中，筆者從共中心點道集中地震傾角的定義出發，提出一種適應于低信噪比條件下的傾角求取方法，并用于定義新型Seislet變換。經典反射波雙曲時距關系[24]為

式中：t(x)是炮檢距為x時的雙程旅行時；t0是零炮檢距雙程旅行時；vn(t0,x)是第n層的疊加或均方根速度。根據地震同相軸傾角的定義p=dt/dx，由方程(3)得到

由傾角的公式(公式(4))可以看出，如果能夠求出對應于每一點(t，x)的速度vn(t0,x)，就能得到隨炮檢距x和時間t而逐漸變化的傾角。在共中心點道集中，速度掃描是一種能夠在低信噪比條件下穩定求取均方根速度的方法，利用得到的掃描速度，根據公式(4)能夠準確地計算出地震局部傾角。

由于經典反射波雙曲時距關系是在水平界面和近炮檢距假設基礎上提出的，因此通過公式(4)計算的同相軸傾角也需要滿足這2個條件。在不同地層產狀的地震模型下，只要滿足該假設條件，本方法能夠得到準確的時空變化傾角，在低信噪比條件下也具有較好的適應性。對于大炮檢距和各向異性所引起的非雙曲反射同相軸，可以選取非雙曲時距曲線方程[25]進行推導，此時將會增加額外參數，多參數掃描方法能夠提供穩定的求解方案。

1.2 閾值去噪

閾值去噪是一種應用稀疏變換進行信噪分離的方法。信號在變換域內能量主要集中在有限的幾個系數中，而噪聲的能量分布于整個變換域內，因此，經稀疏變換后，信號的變換系數幅值大于噪聲的變換系數幅值。如果找到一個合適的閾值，當變換系數小于該閾值時，認為這時的系數主要是由噪聲引起的，當變換系數大于該閾值時，則認為其主要是由信號引起的；當對變換系數進行閾值處理后，反變換就可以達到去除噪聲而保留有用信號的目的。

1.2.1 硬閾值與軟閾值

硬閾值處理是把信號稀疏變換系數的絕對值與閾值比較，小于閾值的變換系數變為0，大于閾值的變換系數不變，再根據反稀疏變換進行信號重建，如圖1a所示。

其中：w是閾值處理后的變換域系數；d是含噪信號的變換域系數；λ是閾值。

軟閾值處理是把稀疏變換系數中大于閾值的變為該值與閾值的差值，如圖1b所示。

硬閾值與軟閾值函數的優缺點比較如表1所示。

1.2.2 改進的閾值方法

本文引入語音信號中改進的新閾值函數[26]：

w=

表1 硬閾值與軟閾值特點

a. 硬閾值方法；b. 軟閾值方法；c, d. 改進的閾值方法。圖1 不同閾值方法比較Fig.1 Comparison of different threshold methods

其中：α為形狀因子并且0≤α≤1.0；μ為壓縮系數。比較不同形狀因子和壓縮系數的閾值曲線(圖1c，d)可以看出：當α=0.5，μ=10 000時過渡比較平滑，故本文采用α=0.5，μ=10 000進行處理。另外，由圖1c和圖1d還可以看出：改進閾值去噪方法的閾值函數是連續的，噪聲分量隨著變換域系數增大而逐漸減小；并且在|d|<λ時，對小波系數進行壓縮處理，而不是簡單地置為零，避免了將有用信號當作噪聲濾除掉。

2 理論模型分析

2.1 基于時距關系地震傾角的新型Seislet變換

為了驗證本文提出的局部地震傾角的有效性，將基于時距關系和平面波分解濾波器計算出的局部地震傾角進行比較。首先建立一個經典的反射波雙曲時距關系模型(圖2a)，道間距為20 m，共128道。圖2b是無噪聲條件下用平面波分解濾波器計算出的傾角。可以看出，平面波分解濾波器在高信噪比條件下能夠獲得準確的傾角。

接下來，在圖2a的模型中加入隨機噪聲，使得信噪比(SNR)為-8.20 dB(圖2c)。一般認為，信號能量與噪聲能量之比大于3∶1，即SNR>4.77 dB時為高信噪比；相反，信號能量與噪聲能量之比小于3∶1，即SNR<4.77 dB時為低信噪比。圖2d是低信噪比條件下用平面波分解濾波器計算出的傾角。此時，噪聲對于平面波分解濾波器產生較大的影響，傾角值有很大的畸變。

對含噪模型(圖2c)進行速度分析，結果如圖3a所示。可以看出：在共中心點道集中，速度掃描能夠在低信噪比條件下穩定求取均方根速度。圖3b為基于時距關系求得的地震局部傾角。

將求得的傾角(圖2d和圖3b)分別用于Seislet變換，得到2種情況的變換域系數(圖3c和圖3d)。對比兩圖可以看出：雙曲型反射同相軸都得到了較好的壓縮；但是PWD對應的變換域系數相比基于時距關系的變換域系數更加集中在小尺度范圍內。這主要是因為，圖3c是由圖2c中的信號與隨機噪聲同時沿圖2d中畸變的傾角進行了壓縮而得到的，而圖3d是由圖2c中的信號沿圖3b中準確的傾角進行壓縮、隨機噪聲無法沿圖3b中準確的傾角進行壓縮而得到的。

a. 理論模型；b. 圖2a的傾角；c. 低信噪比模型；d. 圖2c的傾角。p為慢度。圖2 不同信噪比條件下平面波分解濾波器求得的傾角比較Fig.2 Comparison of slopes of PWD method with different SNR

a. 速度譜；b. 基于時距關系的局部地震傾角；c. 基于PWD的Seislet變換系數；d. 基于時距關系的Seislet變換系數。圖3 低信噪比模型不同方法變換系數比較Fig.3 Comparison of transform coefficients of different methods with low SNR

2.2 閾值去噪

首先選取在尺度方向的簡單閾值方法，即將基于PWD和時距曲線關系的Seislet變換域中尺度大于5的Seislet變換系數置零(圖4a，b)，然后做Seislet反變換得到去噪結果，如圖4c和圖4d所示。可以看到：基于PWD的Seislet反變換去噪結果同相軸產生較大的扭曲，基于時距曲線關系的Seislet反變換去噪結果比較合理；但是2幅圖都含有一定程度的噪聲，考慮用改進的閾值去噪方法對圖像進行改進。

a. 置零后基于PWD的Seislet變換系數；b. 置零后基于時距關系的Seislet變換系數；c. 基于PWD的Seislet去噪結果；d. 基于時距關系的Seislet去噪結果。圖4 低信噪比模型不同傾角的Seislet反變換效果比較Fig.4 Comparison of inverse Seislet transform of different slopes with low SNR

a. 硬閾值方法；b. 軟閾值方法；a. 改進的閾值方法。圖5 不同閾值方法的處理效果比較Fig.5 Comparison of denoising results by using different threshold methods

3 實際數據處理

在實際數據處理中，使用墨西哥灣某海洋數據[28]來評估新方法的去噪效果。如圖6a所示，在時間-共中心點-炮檢距域中，左下角剖面代表炮檢距為0.134 km的共中心點道集，右下角剖面表示共中心點為11.892 km的共炮檢距剖面，左上角剖面代表雙程旅行時為2 s的時間切片。由速度分析(圖6b)的結果可以看出：該數據體層間多次波影響較小，速度譜中能量團比較集中。圖6c為提取的均方根速度，可以看出，該區域速度場比較平穩。圖6d為基于時距關系求得的局部地震傾角。將此傾角用于Seislet變換，如圖7a所示，反射信息主要集中在很小的尺度范圍內，噪聲卻分布在整個變換域內。通過閾值估計模型得到閾值λ=1 251 348.25。然后用改進的閾值方法對基于時距關系的Seislet變換域系數做處理，Seislet反變換得到的去噪結果如圖7b所示。可以看出：改進的閾值方法去除了大部分隨機噪聲，提高了信噪比。

a. 實際地震數據；b. 速度譜；c. 均方根速度；d. 新方法求得的局部地震傾角。圖6 實際地震數據的局部地震傾角Fig.6 Local event slopes of field seismic data

a. 基于時距關系的Seislet變換；b. 改進的閾值方法處理效果。圖7 改進的閾值方法的處理效果Fig.7 Denoising results by using improved threshold method

Table 2SNRof the results by using different threshold methods

方法SNR/dB低信噪比模型-8.20硬閾值方法3.70軟閾值方法3.11改進的閾值方法4.84

4 結語

本文介紹了一種基于共中心點道集反射波時距關系的傾角計算方法，通過與平面波分解濾波器計算出的局部地震傾角對比，該方法能夠求取出更適合低信噪比條件下的傾角，將這種傾角與Seislet變換框架結合，能夠獲得適用于低信噪比條件下的改進型Seislet變換方法。引入語音信號中改進閾值方法并與傳統的硬閾值以及軟閾值方法進行對比，理論模型的處理結果和去噪效果衡量指標證明這種改進的閾值方法不但適于地震信號，而且在提高信噪比方面優于傳統閾值方法，是一種更加適合低信噪比條件下Seislet的變換閾值去噪方法。通過對實際數據的處理，驗證了本文方法的正確性和優越性。

附錄A 小波變換提升算法

提升算法的基本思想在于通過一個基本小波，逐步構建出一個具有更加良好性質的新的小波。一個規范的提升算法有3個步驟：分解、預測、更新。

1)將數據e0分解成為2個小的子集e-1和r-1。假定相鄰的數據間有最大的相關性(在實際中也往往是這種情況)，按照數據的奇偶序號對數據列進行間隔采樣：

2)設一個與數據無關的預測函數P，可以用r-1和e-1的預測值之間的差值來代替r-1。如果使用合理的預測，那么差值將包含比原來r-1少得多的信息，產生的小波系數r-1中保留的是原信號的高頻分量，也即信號中的隨機信息分量。記

r-1?

例如最簡單的預測算子使用偶數序列中2個相鄰數的均值作為它們之間奇數序列的預測值：

3)預測出原信號中的高頻分量后，還需要求出信號中的低頻分量，也即信號中的相關信息分量。通過設計更新算子U對e-1更新，就可以求出表示低頻分量的小波變換尺度系數：

e-1?

U(r-1,k)=(r-1,k+r-1,k-1)/4=

此時，已經可以用較小的數據序列e-1和小波系數r-1來代替原始的數據。如果有一個好的預測，那么2個子集{e-1,r-1}將產生比原來的序列e0更為緊湊的表示。接下來對這個算法進行周期重復。將e-1抽取成兩個序列e-2和r-2，然后用r-2和r-2預測值之間的差值來代替r-2，用r-2對e-2進行更新，可以用小波表示{e-2,r-2,r-1}取代e0。這樣經過n步，就可以用小波表示{e-n,r-n, …,r-1}來取代原來的數據序列。假使小波系數是由基于某種相關性的預測模型得到的，那么，{e-n,r-n, …,r-1}將給出一個比原序列更為緊湊的表示。然而，地震數據空間方向的數據分布特點并不滿足小波變換的假設。

[1] 張波，印興耀，吳國忱，等.用小波分析和反演方法聯合去除隨機噪聲的研究[J].石油地球物理勘探，1999，34(6)：635-641. Zhang Bo, Yin Xingyao, Wu Guochen, et al. Random Noise Elimination Using Wavelet Analysis and Inversion[J]. Oil Geophysical Prospecting, 1999, 34(6): 635-641.

[2] 王典，劉財，劉洋，等.基于提升算法和百分位數軟閾值的小波去噪技術[J]. 地球物理學進展，2008， 23(4)：1124-1130. Wang Dian, Liu Cai, Liu Yang, et al. Application of Wavelet Transform Based on Lifting Scheme and Percentiles Soft Threshold to Elimination of Seismic Random Noise[J]. Progress in Geophysics, 2008, 23(4): 1124-1130.

[3] 樊計昌，劉明軍，王夫運，等.小波包節點域和空間域傾角掃描高階相關去噪技術[J]. 石油地球物理勘探，2009，44(6)：695-699. Fan Jichang, Liu Mingjun, Wang Fuyun, et al. Dip Scanning High Order Correlation Denoise Technique in Wavelet Packet Node Field and Space Domain[J]. Oil Geophysical Prospecting, 2009, 44(6): 695-699.

[4] Zhang Yi, Cheng Lizhi. A New Wavelet Denoising Method of Seismic Signals Based on a Recursive Optimal Thresholding[J]. Acta Scientiarum Naturalium Universitatis Nankaiensis, 2011, 4: 14.

[5] 王志強，韓立國，鞏向博，等. 起伏地表檢波器組合響應[J]. 吉林大學學報:地球科學版，2014,44(2)：694-703. Wang Zhiqiang，Han Liguo，Gong Xiangbo，et al．Geophone Array Response on Undulating Surface[J]．Journal of Jilin University：Earth Science Edition，2014，44(2)：694-703.

[6] Chauris H, Nguyen T. Seismic Demigration/Migration in the Curvelet Domain[J]. Geophysics, 2008, 73(2): 35-46.

[7] Herrmann F J, Wang D, Hennenlent G, et al. Curvelet-Based Seismic Data Processing: A Multiscale and Nonlinear Approach[J]. Geophysics, 2008, 73(1): 1-5.

[8] Pennec E L, Mallat S. Sparse Geometric Image Representation with Bandelets[J]. IEEE Trans Image Process, 2005, 14(4): 423-438.

[9] Do M N, Vetterli M. The Contourlet Transform: An Efficient Directional Multiresolution Image Representation[J]. IEEE Trans Process, 2005, 14(12): 2091-2106.

[10] Starck J L, Candes E J, Donoho D L. The Curvelet Transform for Image Denoising[J]. IEEE Trans Image Process, 2002, 11(6): 670-684.

[11] Velisavljevic V. Directionlets: Anisotropic Multi-Directional Representation with Separable Filtering[D]. Lausanne: Ecole Polytechnique Federale de Lausanne, 2005.

[12] Fomel S. Towards the Seislet Transform[J]. SEG Extended Abstracts, 2006, 25(1):2847-2851.

[13] Fomel S, Liu Y. Seislet Transform and Seislet Frame[J]. Geophysics, 2010, 75(3): 25-38.

[14] 劉洋，Fomel S，劉財，等．高階Seislet變換及其在隨機噪聲消除中的應用[J]．地球物理學報, 2009, 52(8): 2142-2151. Liu Yang, Fomel S, Liu Cai, et al. High-Order Seislet Transform and Its Application of Random Noise Attenuation[J]. Chinese Journal of Geophysics, 2009, 52(8): 2142-2151.

[15] 張鵬，陸文凱. 利用局部傾角的地震成像研究[C]//中國地球物理學會第二十六屆年會、中國地震學會第十三次學術大會論文集. 寧波：中國地球物理學會，中國地震學會，2010. Zhang Peng, Lu Wenkai. Study of Seismic Imaging by Using Local Event Slopes[C]//Papers of the Twenty-Sixth Annual Meeting of Chinese Geophysical Society and the Thirteenth Academic Conference of Seismological Society of China. Ningbo: Chinese Geophysical Society, Seismological Society of China, 2010.

[16] Fomel S. Applications of Plane-Wave Destruction Filter[J]. Geophysics, 2002, 67(6):1946-1960.

[17] Bona A. Velocity-Less Migration Using Horizontal Slownesses[J]. EGU General Assembly, 2009, 11:EGU2009-3883.

[18] Douma H. Leading-Order Seismic Imaging Using Curvelets[J]. Geophysics, 2007, 72(6): 231-248.

[19] Donoho D L, Johnstone J M. Ideal Spatial Adaptation by Wavelet Shrinkage[J]. Biometrika, 1994, 81(3): 425-455.

[20] Gao H Y, Bruce A G. WaveShrink with Firm Shrinkage[J]. Statistica Sinica, 1997, 7(4): 855-874.

[21] 李沖泥，胡光銳.一種改進的子波域語音增強方法[J].通信學報，1999，20(4)：88-91. Li Chongni, Hu Guangrui. A Modified Wavelet Domain Speech Enhancement Method[J]. Journal of China Institute of Communications, 1999, 20(4) : 88-91.

[22] 許麗群.小波閾值去噪改進算法研究[J].電子測量技術， 2010，33(8)：43-45. Xu Liqun.Research of Improved Algorithm of De-Noising in Wavelet Threshold[J].Electronic Measurement Technology, 2010, 33(8): 43-45.

[23] 張新超. Seislet變換去除疊前線性干擾[D]. 成都：成都理工大學，2009. Zhang Xinchao. Removal of Prestack Linear Noise by Seislet Transform[D]. Chengdu: Chengdu University of Technology, 2009.

[24] Fomel S. Velocity-Independent Time-Domain Seismic Imaging Using Local Event Slopes[J]. Geophysics, 2007, 72(3): 139-147.

[25] Fomel S, Grechka V. Nonhyperbolic Reflection Moveout of P Waves: An Overview and Comparison of Reasons[M]. Golden:Center for Wave Phenomena Colorado School of Mines, 2001.

[26] 鄧玉娟.基于小波變換的語音閾值去噪算法研究[D].重慶：重慶大學，2009. Deng Yujuan. Study of Threshold De-Noising Algorithm of Speech Signal Based on Wavelet Transformation[D]. Chongqing: Chongqing University, 2009.

[27] 余晃晶.小波降噪閾值選取的研究[J].紹興文理學院學報：自然科學版，2004，24(9)：34-38. Yu Huangjing. Slection of a Wavelet Noise-Reduction Threshold[J]. Journal of Shaoxing University: Natural Science Edition, 2004,24(9): 34-38.

[28] Claerbout, J F. Basic Earth Imaging: Stanford Exploration Project[EB/OL].[2012-10-10]. http:// www. stanford. edu/sep/prof/.

Threshold Denoising Method Based on New Seislet Transform in the Condition of LowSNR

Liu Cai1, Cui Fangzi1, Liu Yang1, Wang Dian1, Liu Dianmi2, Zhang Peng1

1.CollegeofGeoExplorationScienceandTechnology,JilinUniversity,Changchun130026,China
2.ResearchInstituteofGeophysicalProspecting,JilinOilfieldCompany,CNPC,Songyuan138000,Jilin,China

Seislet transform is a wavelet-like transform based on wavelet base that analyzes seismic data along variable local event slopes. The calculation of the local event slopes is the key step of this method.There are many kinds of methods to calculate the local slopes. But in the condition of lowSNR, these methods have some limits. We propose a method that is suitable for lowSNRbased on the definition of the slopes oft-xrelationship in CMP traces. Comparing the slopes oft-xrelationship in CMP traces with PWD method, we can see that the new method is more accurate in calculating the local slopes. Apply these slopes to Seislet frame to establish a new Seislet transform that characterize lowSNRdata. In seismic data processing, we introduce one method which works well for acoustic signals. We propose the new threshold method with Seislet transform. It is not only suitable for seismic signals, but also better than the traditional threshold method to improveSNR. The results show that the combination of the new Seislet transform and improved threshold method can extract signal in the condition of lowSNReffectively.

Seislet transform; lowSNR; local event slope; threshold; denoising

10.13278/j.cnki.jjuese.201501303.

2014-01-16

國家自然科學基金項目(41274119, 41174080)；國家“863”重大計劃項目(2012AA09A20103)

劉財(1963——)，男，教授，博士生導師，主要從事地震波場正反演理論、綜合地球物理等研究，E-mail:liucai@jlu.edu.cn

劉洋(1979——)，男，教授，主要從事地震數據處理，海洋電磁數據處理和地質-地球物理綜合研究等工作，E-mail:yangliu1979@jlu.edu.cn。

10.13278/j.cnki.jjuese.201501303

P631.4

A

劉財，崔芳姿，劉洋，等. 基于低信噪比條件下新型Seislet變換的閾值去噪方法.吉林大學學報:地球科學版,2015,45(1):293-301.

Liu Cai, Cui Fangzi, Liu Yang, et al. Threshold Denoising Method Based on New Seislet Transform in the Condition of LowSNR.Journal of Jilin University:Earth Science Edition,2015,45(1):293-301.doi:10.13278/j.cnki.jjuese.201501303.