基于Contourlet域自適應Wiener閾值的同時震源波場分離

王坤喜, 毛偉建

1 中國科學院精密測量科學與技術創新研究院計算與勘探地球物理研究中心; 大地測量與地球動力學國家重點實驗室， 武漢 430077 2 中國科學院大學， 北京 100049

0 引言

同時震源采集技術，也稱多震源地震數據混合采集技術，可以在常規采集一個單炮地震記錄的時間內采集多個炮的地震記錄.當施工的周期不變時，通過增加激發的震源數量，提高對地下構造的覆蓋次數，改善對地下構造的照明質量；當對目標區域的覆蓋次數不變時，多個分布在不同區域的震源同時激發，可以大大縮短施工周期，提高采集效率(Li et al.，2013；Zu et al.，2017；張攀和毛偉建，2018).由于相鄰震源之間的干擾會產生“混合噪聲”，同時震源數據不能直接用于傳統地震數據處理流程，需要先將混疊波場進行分離處理，得到常規的單炮數據.分離方法通常可以分為被動分離和主動分離兩大類(Akerberg et al.，2008；Liu et al.，2014；Magnussen，2015).

被動分離是將分離問題轉換為去噪問題，將經過延遲時間編碼的干擾炮信號視為噪聲.當不同震源之間存在隨機時間差時，混合噪聲在非共炮域(如共檢波點域、共中心點域、共偏移距域等)不再具有連續性，表現為偽隨機脈沖，而有效信號仍然保持連續相干性(Doulgeris et al.，2012).基于上述性質，Moore等(2008)在共偏移距道集進行濾波分離.Huo等(2009)利用多方向矢量中值濾波法對地層的傾角進行了掃描，求得最佳傾角方向后利用中值濾波去除共中心點域的混合噪聲.Kim等(2009)基于數據本身建立了一個噪聲估計模型，然后從采集數據中自適應的提取噪聲，該算法適用于共炮檢距道集.Zhang和Olofsson(2012)采用加權τ-p變換壓制鄰炮干擾.Yang等(2017)提出先用PWD(plane-wave destruction)求取地層的傾角信息，再結合矢量中值濾波對串擾噪聲進行壓制.

另一種分離類方法是主動分離，也被稱為反演類方法，該方法是將分離問題提煉成數學中的反問題，利用反演類方法進行求解.Akerberg等(2008)、Moore(2010)和劉強等(2014)在Radon域進行了稀疏約束下的迭代反演分離.Lin和Herrmann(2009)在Curvelet域進行稀疏表示，利用壓縮感知對混合的波場進行分離.Abma等(2010)在頻率域通過稀疏約束的POCS(projection onto convex sets)方法使不相干信號能量最小化，進而獲取分離信號.Tan等(2013)基于互易定理，利用稀疏反演的方法進行波場分離.Chen等(2014)提出整形正則化反演的方法并在Seislet域利用閾值進行約束反演.祖紹環等(2016)在Curvelet域進行稀疏約束，同樣利用閾值進行約束反演.宋家文等(2019)將地震數據轉換到頻率—波數—波數域(FKK)，提出了一種高效混采數據分離的方法.在稀疏迭代分離框架下，除了使用固定的變換域外，一些學者也考慮采用字典學習來完成分離(Chen，2017；Zhou et al.，2017)，但是字典學習的方法需要很大的計算開銷.

在同時震源稀疏域迭代反演中，選擇一個合適的稀疏基十分重要,它決定了分離效果和迭代收斂的速度.Do和Vetterli(2002)首次提出了Contourlet變換理論，將拉普拉斯金字塔(Burt and Adelson，1983；Do and Vetterli，2003)和二維方向濾波器組串聯在一起組成金字塔多方向濾波器組(Do and Vetterli，2001)，可以將信號數據分解成不同尺度下的方向子帶，很好地抓住數據的幾何結構特征，滿足曲線的各向異性尺度關系，能夠快速并結構化分解采樣信號.Contourlet變換將多尺度分析和方向分析分開進行，首先由拉普拉斯金字塔變換對圖像進行多尺度分解以“捕獲”點奇異，接著由方向濾波器組將分布在同方向上的奇異點合成一個系數.Contourlet變換的最終結果是用類似于線段的基結構來逼近原圖像，相對于Wavelet變換和Curvelet變換等多尺度分析工具，Contourlet變換具有更少的冗余度和更好的逼近性能，因此Contourlet變換已在信號處理、圖像處理等領域得到廣泛應用，如信號去噪(Eslami and Radha，2003)，目標檢測(Wilbur et al.，2009)，圖像壓縮(Chang et al.，2000a)等等.

與數字圖像去噪的目標相似，同時震源的稀疏域迭代反演分離需要在去除串擾噪聲的同時盡可能地保留地震數據中的弱信號和細節信息(Zu et al.，2017)，合理地選取閾值非常重要，直接影響到分離效果(Yang et al.，2013).Donoho和Johnstone(1994)在小波域提出了閾值收縮方法，給出了Universal閾值，并從漸進意義上證明了Universal閾值的最優性(Donoho，1995)，然而Universal閾值是閾值的上限，但不是最佳收縮閾值，它過度“扼殺”了小波系數，丟失信號的細節部分.Donoho和Johnstone(1995)又提出了SURE閾值，但是該閾值偏于保守，在信噪比很小的情況下，其效果不如Universal閾值.Zhang和Desai(1998)改進SURE閾值提出了SureShrink閾值，可以求出理想閾值的估計值，但沒有顯式表達式，閾值的計算往往需要預先知道信號本身，但實際求取中這是不可能的.Chang等(2000a)假設原始不帶噪聲信號在小波域的系數服從廣義高斯分布，在貝葉斯框架下通過最小化貝葉斯風險，得到了著名的BayesShrink閾值，該閾值考慮了信號在小波域的先驗知識，具有很好的去噪效果.在借鑒圖像壓縮編碼中的上下文模型方法后，Chang等(2000b)又提出了具有局部自適應能力的BayesShrink閾值，能夠靈活實現去噪與保留信號之間的平衡.在貝葉斯最大后驗概率估計理論下，Moulin和Liu(1999)假設圖像子帶的小波系數服從拉普拉斯分布或者高斯分布，可以分別推導得到MapShrink閾值與Wiener閾值，這兩個閾值算子都具有較好的去噪效果.

以上閾值方法的提出主要用于數字圖像去噪中，一些學者也在研究適用于同時震源稀疏域迭代反演分離的閾值算子.Chen等(2014)提出在Seislet域采用百分位法(Yang et al.，2013)分段設置硬閾值，取得了一定的分離效果，后來又在PNMO-MF-FK聯合濾波器下使用線性衰減硬閾值也得到了較好的分離結果.該閾值算子在初始迭代時取系數最大值，并隨著迭代次數線性衰減為零(Chen et al.，2015).Li等(2018)將混疊數據變換到Curvelet域后，發現相干信號主要集中在低頻子帶中，串擾噪聲主要集中在高頻子帶中，通過考慮在不同子帶設置不同的閾值，取得了比全局一致方法更好的效果.Zu等(2018)將該閾值方法推廣并運用于三維同時震源數據分離.在同時震源稀疏域迭代反演分離方法中，另一種常見的閾值是指數衰減閾值(Gao et al.，2010；Cheng and Sacchi，2015；Song et al.，2019).與常數和線性衰減閾值相比，指數衰減閾值具有更快的收斂速度，但是該閾值沒有充分利用稀疏域系數的分布情況等先驗信息，不具備自適應去噪的能力.

本文利用Contourlet變換良好的稀疏表達能力，通過稀疏域閾值收縮來迭代反演分離同時震源數據.我們在子帶一致的Wiener閾值基礎上做出改進，提出具有尺度與空間自適應的Wiener閾值.通過兩個模擬數據和兩個實際數據的應用結果驗證本文提出的自適應Wiener閾值方法的有效性.

1 方法

1.1 稀疏迭代反演分離

在時間域，每個檢波點接收到的同時震源混合記錄可表示為

(1)

其中dobs為混疊道集記錄，db為第b個震源的原始道集記錄，m為混疊度.對角矩陣Γb表示第b個震源的延遲時間矩陣，其表達式為

Γb=F-1diag(e-iω tb,1,e-iω tb,2,…,e-iω tb,u,…)F,

(2)

其中，tb,u為第b個震源第u道激發的延遲時間，i為虛數單位，ω為角頻率，F和F-1表示一對正逆傅里葉變換.通常，在同時震源數據分離中，第一個震源(本文稱為主炮)的延遲時間算子(Γ1)為單位矩陣I，得到連續相干的有效信號；其他震源(本文統稱為干擾炮)在延遲時間算子的作用下會形成串擾噪聲.

同時震源數據的稀疏反演分離問題可以表示如下：

(3)

求解式(3)的一個常用迭代框架是非線性整形正則化框架(Chen et al.，2014)：

(4)

S=C-1TC,

(5)

其中T為閾值參數，C和C-1分別代表稀疏變換及其逆變換，本文中C為Contourlet變換.

共檢波點道集或共偏移距道集中的主炮信號在Contourlet域是聚焦的能量點，而干擾炮信號是分散的串擾噪聲，T可以對變換域系數進行收縮從而有效地壓制串擾噪聲進而分離混疊波場.因此，同時震源數據的稀疏迭代反演分離的關鍵是如何確定合適的T.

1.2 子帶一致Wiener閾值

在每一次迭代下使用閾值T進行收縮去噪前，主炮信號一般是連續相干的，干擾炮信號在延遲時間作用下會以串擾噪聲的形式存在，所以此時分離問題可以看作噪聲去除問題：

Y=X+N,

(6)

其中Y、X和N分別表示變換域中包含噪聲的混疊數據、不含噪聲的原始主炮數據和串擾噪聲數據.在每一次迭代中，混疊數據在Contourlet域被分解為K層，每層L個方向，用k=1,2…,K表示分解的第k個尺度，頻率越高k越大；用l=1,2…,L表示第k個尺度的方向數.在Contourlet域，假設原始主炮數據和串擾噪聲數據的各方向子帶系數都服從高斯分布，根據貝葉斯最大后驗概率估計理論，可以求得Wiener濾波公式(Moulin and Liu，1999)

(7)

圖1 Wiener閾值收縮函數形式示意圖Fig.1 Schematic diagram of the Wiener threshold contraction function

由Wiener閾值收縮原理(式(7))，可得Wiener閾值

(8)

式(8)中σN可以采用魯棒中位數估計法(Donoho and Johnstone，1994，1995)求取

(9)

(10)

其中Num為Yi,j所處子帶的全部系數個數.

1.3 自適應Wiener閾值

子帶一致閾值法能夠取得一定的分離效果，但是由于沒有考慮到串擾噪聲在Contourlet域中不同尺度的分布規律以及有效信號的空間分布規律，該方法在串擾噪聲的壓制效果、有效信號尤其是弱信號的保留和迭代收斂速度方面存在不足.基于以上考慮，我們分別從噪聲方差計算和有效信號方差計算兩個方面進行改進，獲得自適應Wiener閾值.

(11)

(12)

根據式(11)與式(12)，改進后的尺度與空間自適應Wiener閾值算子可以表示為

(13)

在Contourlet域使用式(13)的閾值算子迭代反演分離同時震源數據的方法稱為Contourlet域自適應Wiener閾值法，下文簡稱為自適應閾值法.

2 算例

2.1 Marmousi模型數據

Marmousi模型是石油工業界用于方法測試的一個復雜地質模型.圖2a和圖2b分別為主炮和干擾炮激發采集的未混疊共檢波點道集數據，共767道；圖2b數據在0～1 s隨機延遲時間作用下以偽隨機噪聲的形式加到主炮波場中，形成如圖2c所示的混疊波場，混疊度m=2.模擬計算在Intel Xeon E7-4807、內存為128G的IBM服務器上展開，沒有并行計算.Contourlet變換選擇“9-7”塔式分解濾波器和“pkva”方向濾波器進行三層分解，方向數分別為4、8和8.在本次測試的自適應Wiener閾值算子中，α的取值為0.6.圖3a—b分別為圖2a數據與圖2c數據在Contourlet域中不同尺度不同方向的系數分布情況，可以看出在迭代初期，有效信號的系數集中分布在中低頻部分，而串擾噪聲的系數主要分布在中高頻部分.圖4a為圖2a數據在Contourlet域中不同尺度不同方向的系數統計分布直方圖，可以看出系數多集中分布在零點附近，形態上表現為尖峰長尾，可以用高斯分布來逼近，當分解的尺度越低，得到數值較大的系數越多，進一步證明了連續相干的有效信號集中分布于低頻和中頻子帶中；圖4b為串擾噪聲數據在Contourlet域中不同尺度不同方向的系數統計分布直方圖，可以看出串擾噪聲分布滿足高斯模型，且數值較大的系數多集中存在于中高頻子帶中.

圖2 Marmousi模型中共檢波點道集數據 (a) 原始主炮數據; (b) 原始干擾炮數據; (c) 混疊數據.Fig.2 Data of Marmousi model in common receiver domain (a) Original data of the main source; (b) Original data of the interference source; (c) Blended data.

圖3 Contourlet域不同尺度不同方向分解的子帶系數分布圖 (a) 原始主炮數據; (b) 混疊數據.Fig.3 Distribution of subband coefficients in different scales and different directions in Contourlet domain (a) Original data of the main source; (b) Blended data.

圖4 Contourlet域的不同尺度不同方向子帶系數統計分布直方圖 (a) 原始主炮數據; (b) 串擾噪聲數據.Fig.4 Statistical distribution histogram of subband coefficients in different scales and different directions in Contourlet domain (a) Original data of the main source; (b) Crosstalk noise data.

圖5 不同方法分離主炮的結果和放大10倍的殘差 (a) 子帶一致閾值法分離結果; (b) Curvelet域指數衰減閾值法分離結果; (c) 本文提出的自適應閾值法分離結果; (d)、(e)和(f)分別是分離結果(a)、(b)和(c)所對應的殘差(已放大10倍).Fig.5 Deblended results of the main source from different methods and corresponding residuals with 10 times scaled up (a) Deblended result from the subband consistent threshold method; (b) from the exponential attenuation threshold method in Curvelet domain; (c) from the proposed adaptive threshold method; (d), (e) and (f) are the residuals with 10 times scaled up corresponding to the results (a), (b) and (c), respectively.

圖6 不同方法分離干擾炮的結果和放大10倍的殘差 (a) 子帶一致閾值法分離結果; (b) Curvelet域指數衰減閾值法分離結果; (c) 本文提出的自適應閾值法分離結果; (d)、(e)和(f)分別是分離結果(a)、(b)和(c)所對應的殘差(已放大10倍).Fig.6 Deblended results of the interference source from different methods and corresponding residuals with 10 times scaled up (a) Deblended result from the subband consistent threshold method; (b) from the exponential attenuation threshold method in Curvelet domain; (c) from the proposed adaptive threshold method; (d), (e) and (f) are the residuals with 10 times scaled up corresponding to the results (a), (b) and (c), respectively.

圖7 原始主炮數據、混疊數據和用三種分離方法得到的分離后主炮數據的單道波形對比圖 右上角子圖為紅色虛線矩形框的放大圖，‘A’表示子帶一致閾值法，‘B’表示Curvelet域指數衰減閾值法，‘C’表示本文提出的自適應閾值法. (a) 圖2a、圖2c和圖5a—c中第110道數據; (b) 圖2a、圖2c和圖5a—c中第300道數據.Fig.7 The single trace waveform comparison of the original data, blended data, and deblended data from the three deblending methods for the main source The upper right subgraph is the zoomed-in section of the red rectangle box. ‘A’ is for the subband consistent threshold method, ‘B’ for the exponential attenuation threshold method in Curvelet domain, and ‘C’ for the proposed adaptive threshold method. (a) The 110th trace in Fig.2a,Fig.2c and Figs.5a—c; (b) The 300th trace in Fig.2a,Fig.2c and Figs.5a—c.

圖8 三種分離方法得到單道分離誤差對比圖 藍色虛線表示子帶一致閾值法分離得到的誤差，綠色實線表示用Curvelet域指數衰減閾值法分離得到的誤差，紅色實線表示 本文自適應閾值法分離得到的誤差. (a) 圖5d—f中第110道數據; (b) 圖5d—f中第300道數據.Fig.8 A comparison of the errors for deblended single trace obtained from the three deblending methods The blue dotted line is the error for the subband consistent threshold method, the green solid line for the exponential attenuation threshold method in Curvelet domain, and the red solid line for the proposed adaptive threshold method. (a) The 110th trace in Figs.5d—f; (b) The 300th trace in Figs.5d—f.

圖9 不同分離方法處理Marmousi模型數據的分離信噪比曲線 (a) 子帶一致閾值法; (b) Curvelet域指數衰減閾值法; (c) 本文提出的自適應閾值法.Fig.9 The SNR of the Marmousi model data from the different deblending methods (a) The subband consistent threshold method; (b) The exponential attenuation threshold method in Curvelet domain; (c) The proposed adaptive threshold method.

圖10 SEAM模型數據 (a) 原始主炮數據; (b) 原始干擾炮數據; (c) 混疊數據.Fig.10 Data of the SEAM model (a) Original data of the main source; (b) Original data of the interference source; (c) Blended data.

圖11 本文提出的自適應閾值法分離結果和殘差 (a) 主炮分離結果; (b) 干擾炮分離結果; (c) 圖10a與圖11a之間殘差; (d) 圖10b與圖11b之間殘差.Fig.11 Deblended results and the corresponding residuals from the proposed adaptive threshold method (a) Deblended result of the main source; (b) Deblended result of the interference source; (c) Residual betweenFig.10a andFig.11a; (d) Residual betweenFig.10b andFig.11b.

在相同的迭代反演框架下，分別用子帶一致閾值法、Curvelet域指數衰減閾值法(祖紹環等，2016)和本文提出的自適應閾值法分離圖2c中的混疊數據，證明自適應閾值法的優越性.圖5a—c為主炮的分離結果，圖6a—c為干擾炮的分離結果.從分離結果可以看出，這三種方法都取得了較好的分離效果，但是自適應閾值法能夠保護有效信號特別是深部弱信號，殘留噪聲更少.為了更好地比較這三種方法的分離效果，分別將圖2a與圖5a—c數據相減、將圖2b與圖6a—c數據相減得到殘差，然后將殘差放大10倍后得到如圖5d—f、圖6d—f的放大后殘差圖.從殘差圖可以看出，子帶一致閾值法沒能很好地壓制串擾噪聲，同時也造成了有效信號的泄漏，殘差大，分離結果較差；Curvelet域指數衰減閾值法對深部弱信號保留較好，但是串擾噪聲強的部位不能夠被很好地壓制；本文提出的自適應閾值法能夠有效的壓制串擾噪聲，更好地恢復和保護深層弱信號，分離結果好，殘差小.為了更加細致的對比分離效果，我們分別提取原始主炮數據(圖2a)、混疊數據(圖2c)和分離結果(圖5a—c)中第110和300道的數據，繪制如圖7所示的單道波形對比圖.從單道波形對比圖中可以發現，這三種方法都能去除混疊數據中的噪聲，分離的結果都與原始未混疊數據接近，但從圖中放大部位可以看出本文自適應閾值法分離得到的數據與原始主炮數據更相近.圖8為從殘差圖5d—f中抽取的第110和300道數據，反映了分離誤差的大小，能夠進一步展示圖7中的三種方法分離得到的結果與原始主炮數據的接近程度.圖8中的數值越小，證明分離結果與原始主炮數據越接近，分離效果越好.可以看出，自適應閾值法分離得到的誤差比子帶一致閾值法和Curvelet域指數衰減閾值方法分離得到的誤差小很多，在4s以后的弱信號部位，自適應閾值法得到的分離誤差幾乎為0，由此證明了本文自適應閾值法具有良好地保護弱信號的能力.

為了定量分析分離效果，我們定義如下信噪比公式：

(14)

2.2 SEAM模型數據

SEAM模型包含傾斜地層、鹽丘等復雜地質構造，是工業界常用的另一個測試模型.圖10a和圖10b分別為主炮和干擾炮激發采集的未混疊的共檢波點道集數據，共1024道，可以看到鹽丘產生的反射層信號能量很強.在0～1 s隨機延遲時間作用下，干擾炮波場會以偽隨機噪聲的形式疊加到主炮波場中，混疊度m=2，形成如圖10c的混疊波場.本次測試的分離方法是本文提出的自適應閾值法，加速因子λ2中α為0.7，其他測試條件與前文相同.圖11a為主炮分離結果，收斂后的信噪比為27.88 dB，圖11b為干擾炮分離的結果，收斂后的信噪比為31.31 dB，可以看出，使用自適應閾值法分離SEAM模型的同時震源數據同樣具有較高信噪比，有效信號幾乎被完全恢復，位于鹽丘內部的弱反射信號得到了很好的保留，同樣證明了自適應閾值法能夠很好地保護弱信號，具有良好的分離效果；圖11c為圖10a與圖11a之間的殘差，圖11d為圖10b與圖11b之間的殘差，可以看出，除了反射波能量過大的部位殘留少量噪聲外，大部分串擾噪聲得到了很好地壓制，殘差較小.圖12為分離得到的信噪比曲線，通過40次迭代實現收斂，用時93.2 s，效率高.

2.3 第一個實際數據

使用Petroleum Geo-Services(PGS)公司公開的海上實際數據測試本文提出的自適應閾值法的分離效果.為了獲得更加復雜的同時震源數據，我們截取原數據進行了人為混疊.實際數據的炮數和檢波點數都是256，采樣間隔為4 ms，混疊度m=2，加速因子λ2中α為0.7，其他條件與前文相同.圖13a為原始主炮數據；圖13b為原始干擾炮數據；原始干擾炮數據在2 s內的延遲時間編碼下與原始主炮數據混疊得到混疊數據，如圖13c所示.可以看出，在原始主炮數據與原始干擾炮數據中，淺部信號的能量都很強，遠大于深部弱信號.

使用自適應閾值法分離圖13c中的混疊數據.圖14a為主炮分離結果，收斂后的信噪比為23.37 dB，圖14b為干擾炮分離的結果，收斂后的信噪比為21.27 dB，圖14c為圖13a與圖14a之間的殘差，圖14d為圖13b與圖14b之間的殘差.圖15為分離的信噪比曲線.自適應閾值法取得了較高的分離信噪比，能很好地壓制串擾噪聲，保護深層弱信號，由殘差數據(圖14c—d)可見深部幾乎沒有有效信號泄漏，本次測試在第37次迭代時實現了收斂，用時89.2 s，證明了本文自適應閾值法對實際的同時震源數據同樣具有優越的分離效果和高效的分離速度.

2.4 第二個實際數據

使用某工區實際采集的同時震源地震數據進一步驗證自適應閾值法的有效性.檢波點間距和炮點間距均為25 m，每道共1020個采樣點，采樣間隔為4 ms，混疊度m=2.在本次測試中，加速因子λ2中α為0.7，其他測試條件與前文相同.圖16a為混疊的共偏移距道集，串擾噪聲散布在整個道集中，深部弱有效信號被湮沒.使用本文提出的自適應閾值法得到的主炮和干擾炮的分離結果分別如圖16b和圖16c所示.從分離結果可見淺部強反射同相軸連續性更好，深部弱反射同相軸清晰可見，強噪聲被完全壓制，剖面的信噪比顯著提高.本次測試進一步證明了自適應閾值法的有效性.

3 討論

在自適應閾值法中，衰減系數α能夠促進和加快迭代的收斂，我們同樣使用Marmousi模型數據討論α對分離效果的影響.當α的取值分別為0.2，0.3，0.4，0.5，0.6，0.7，0.8和0.9時，使用自適應閾值法分離圖2c中的混疊數據得到如圖18所示的主炮分離信噪比曲線.在圖18中，隨著α的減小，收斂需要的迭代次數減少，收斂速度明顯加快，但是當α小于0.6時，收斂后的信噪比也出現了明顯的下降；當α取值在0.7～0.9之間時，收斂后的信噪比一致且較高；當α為0.9時，收斂后信噪比較高，但是收斂需要的迭代次數明顯多于α為0.8時收斂需要的迭代次數.通過討論可以發現，α越小，收斂速度越快，但是過小的α會影響分離效果，所以我們需要在分離效果和收斂速度之間做出一個平衡，而α取值為0.6～0.8是一個較為合理的選擇.

圖12 本文提出的自適應閾值法分離 SEAM模型數據的信噪比曲線 其中實線代表主炮的分離信噪比，虛線代表干擾炮的分離信噪比.Fig.12 The SNR with respective to the iteration number for SEAM model data using the proposed adaptive threshold method The solid line is for the separated SNR of the main source, and the dashed line for the separated SNR of the interference source.

4 結論

在同時震源數據的迭代反演分離過程中，待處理地震數據經過Contourlet變換后，有效信號和串擾噪聲在不同尺度和不同方向的子帶中具有不同的分布規律，因此使用子帶一致Wiener閾值方法往往得不到好的分離結果.在稀疏域迭代反演框架下，本文提出了一種在Contourlet域分離同時震源數據的自適應Wiener閾值方法.在迭代初期，自適應Wiener閾值方法中的噪聲方差計算具有尺度自適應能力，能夠適當增大高頻子帶中的噪聲方差估計量，相比于使用全局一致的噪聲方差計算方法，更利于達到快速壓制串擾噪聲的目的.在迭代后期，地震數據中的有效信號在稀疏域中多表現為數值較大的能量聚集點，自適應Wiener閾值方法使用局部加窗法計算有效信號方差，比子帶一致的有效信號方差能更加準確地反映實際方差的大小，能夠避免因為有效信號方差計算偏小而導致弱信號泄露.實驗結果證明，本文Contourlet域自適應Wiener閾值方法能夠快速有效地壓制串擾噪聲和保護有效信號，取得了比Contourlet域子帶一致Wiener閾值方法和Curvelet域指數衰減閾值方法更好的分離效果.

在估計噪聲方差時，本文自適應Wiener閾值方法只考慮了噪聲數據在Contourlet域中不同尺度下的分布規律，如何將同一尺度下的不同方向特征也考慮進去需要進一步研究.

圖13 第一個實際數據 (a) 原始主炮數據; (b) 原始干擾炮數據; (c) 混疊數據.Fig.13 The first field data (a) Original data of the main source; (b) Original data of the interference source; (c) Blended data.

圖14 本文提出的自適應閾值法分離第一個實際數據的分離結果和殘差 (a) 主炮分離結果; (b) 干擾炮分離結果; (c) 圖13a與圖14a之間殘差; (d) 圖13b與圖14b之間殘差.Fig.14 Deblended results of the first field data from the proposed adaptive threshold method, and the corresponding residuals (a) Deblended result of the main source; (b) Deblended result of the interference source; (c) Residuals betweenFig.13a andFig.14a; (d) Residuals betweenFig.13b andFig.14b.

圖15 本文提出的自適應閾值法分離 第一個實際數據的信噪比曲線 其中實線代表主炮的分離信噪比，虛線代表干擾炮的分離信噪比.Fig.15 The SNR with respective to the iteration number for the first field data using the proposed adaptive threshold method The solid line is for the SNR of the main source, and the dashed line for the SNR of the interference source.

圖17 不同窗口長度w對本文自適應閾值法 分離效果的影響 其中實線代表主炮的分離信噪比，虛線代表干擾炮的分離信噪比.Fig.17 The effect of different window sizes w on the separation effect from the proposed adaptive threshold method The solid line is for the SNR of the main source, and the dashed line for the SNR of the interference source.

圖16 第二個實際數據的分離測試 (a) 混疊數據; (b) 主炮分離結果; (c) 干擾炮分離結果.Fig.16 The deblending test of the second field data (a) Blended data; (b) Deblended result of the main source; (c) Deblended result of the interference source.

圖18 衰減系數α取不同值時，本文自適應閾值法 得到的主炮分離信噪比曲線Fig.18 The SNR curves during the separation of the main source from the proposed adaptive threshold method by taking different values of attenuation parameter α