王坤喜, 毛偉建

1 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院計(jì)算與勘探地球物理研究中心; 大地測量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 武漢 430077 2 中國科學(xué)院大學(xué)， 北京 100049

0 引言

同時(shí)震源采集技術(shù)，也稱多震源地震數(shù)據(jù)混合采集技術(shù)，可以在常規(guī)采集一個(gè)單炮地震記錄的時(shí)間內(nèi)采集多個(gè)炮的地震記錄.當(dāng)施工的周期不變時(shí)，通過增加激發(fā)的震源數(shù)量，提高對地下構(gòu)造的覆蓋次數(shù)，改善對地下構(gòu)造的照明質(zhì)量；當(dāng)對目標(biāo)區(qū)域的覆蓋次數(shù)不變時(shí)，多個(gè)分布在不同區(qū)域的震源同時(shí)激發(fā)，可以大大縮短施工周期，提高采集效率(Li et al.，2013；Zu et al.，2017；張攀和毛偉建，2018).由于相鄰震源之間的干擾會(huì)產(chǎn)生“混合噪聲”，同時(shí)震源數(shù)據(jù)不能直接用于傳統(tǒng)地震數(shù)據(jù)處理流程，需要先將混疊波場進(jìn)行分離處理，得到常規(guī)的單炮數(shù)據(jù).分離方法通?？梢苑譃楸粍?dòng)分離和主動(dòng)分離兩大類(Akerberg et al.，2008；Liu et al.，2014；Magnussen，2015).

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

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

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

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

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

本文利用Contourlet變換良好的稀疏表達(dá)能力，通過稀疏域閾值收縮來迭代反演分離同時(shí)震源數(shù)據(jù).我們在子帶一致的Wiener閾值基礎(chǔ)上做出改進(jìn)，提出具有尺度與空間自適應(yīng)的Wiener閾值.通過兩個(gè)模擬數(shù)據(jù)和兩個(gè)實(shí)際數(shù)據(jù)的應(yīng)用結(jié)果驗(yàn)證本文提出的自適應(yīng)Wiener閾值方法的有效性.

1 方法

1.1 稀疏迭代反演分離

在時(shí)間域，每個(gè)檢波點(diǎn)接收到的同時(shí)震源混合記錄可表示為

(1)

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

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

(2)

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

同時(shí)震源數(shù)據(jù)的稀疏反演分離問題可以表示如下：

(3)

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

(4)

S=C-1TC,

(5)

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

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

1.2 子帶一致Wiener閾值

在每一次迭代下使用閾值T進(jìn)行收縮去噪前，主炮信號(hào)一般是連續(xù)相干的，干擾炮信號(hào)在延遲時(shí)間作用下會(huì)以串?dāng)_噪聲的形式存在，所以此時(shí)分離問題可以看作噪聲去除問題：

Y=X+N,

(6)

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

(7)

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

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

(8)

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

(9)

(10)

其中Num為Yi,j所處子帶的全部系數(shù)個(gè)數(shù).

1.3 自適應(yīng)Wiener閾值

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

(11)

(12)

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

(13)

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

2 算例

2.1 Marmousi模型數(shù)據(jù)

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

圖2 Marmousi模型中共檢波點(diǎn)道集數(shù)據(jù) (a) 原始主炮數(shù)據(jù); (b) 原始干擾炮數(shù)據(jù); (c) 混疊數(shù)據(jù).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域不同尺度不同方向分解的子帶系數(shù)分布圖 (a) 原始主炮數(shù)據(jù); (b) 混疊數(shù)據(jù).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域的不同尺度不同方向子帶系數(shù)統(tǒng)計(jì)分布直方圖 (a) 原始主炮數(shù)據(jù); (b) 串?dāng)_噪聲數(shù)據(jù).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 不同方法分離主炮的結(jié)果和放大10倍的殘差 (a) 子帶一致閾值法分離結(jié)果; (b) Curvelet域指數(shù)衰減閾值法分離結(jié)果; (c) 本文提出的自適應(yīng)閾值法分離結(jié)果; (d)、(e)和(f)分別是分離結(jié)果(a)、(b)和(c)所對應(yīng)的殘差(已放大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 不同方法分離干擾炮的結(jié)果和放大10倍的殘差 (a) 子帶一致閾值法分離結(jié)果; (b) Curvelet域指數(shù)衰減閾值法分離結(jié)果; (c) 本文提出的自適應(yīng)閾值法分離結(jié)果; (d)、(e)和(f)分別是分離結(jié)果(a)、(b)和(c)所對應(yīng)的殘差(已放大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 原始主炮數(shù)據(jù)、混疊數(shù)據(jù)和用三種分離方法得到的分離后主炮數(shù)據(jù)的單道波形對比圖 右上角子圖為紅色虛線矩形框的放大圖，‘A’表示子帶一致閾值法，‘B’表示Curvelet域指數(shù)衰減閾值法，‘C’表示本文提出的自適應(yīng)閾值法. (a) 圖2a、圖2c和圖5a—c中第110道數(shù)據(jù); (b) 圖2a、圖2c和圖5a—c中第300道數(shù)據(jù).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 三種分離方法得到單道分離誤差對比圖 藍(lán)色虛線表示子帶一致閾值法分離得到的誤差，綠色實(shí)線表示用Curvelet域指數(shù)衰減閾值法分離得到的誤差，紅色實(shí)線表示 本文自適應(yīng)閾值法分離得到的誤差. (a) 圖5d—f中第110道數(shù)據(jù); (b) 圖5d—f中第300道數(shù)據(jù).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模型數(shù)據(jù)的分離信噪比曲線 (a) 子帶一致閾值法; (b) Curvelet域指數(shù)衰減閾值法; (c) 本文提出的自適應(yīng)閾值法.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模型數(shù)據(jù) (a) 原始主炮數(shù)據(jù); (b) 原始干擾炮數(shù)據(jù); (c) 混疊數(shù)據(jù).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 本文提出的自適應(yīng)閾值法分離結(jié)果和殘差 (a) 主炮分離結(jié)果; (b) 干擾炮分離結(jié)果; (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域指數(shù)衰減閾值法(祖紹環(huán)等，2016)和本文提出的自適應(yīng)閾值法分離圖2c中的混疊數(shù)據(jù)，證明自適應(yīng)閾值法的優(yōu)越性.圖5a—c為主炮的分離結(jié)果，圖6a—c為干擾炮的分離結(jié)果.從分離結(jié)果可以看出，這三種方法都取得了較好的分離效果，但是自適應(yīng)閾值法能夠保護(hù)有效信號(hào)特別是深部弱信號(hào)，殘留噪聲更少.為了更好地比較這三種方法的分離效果，分別將圖2a與圖5a—c數(shù)據(jù)相減、將圖2b與圖6a—c數(shù)據(jù)相減得到殘差，然后將殘差放大10倍后得到如圖5d—f、圖6d—f的放大后殘差圖.從殘差圖可以看出，子帶一致閾值法沒能很好地壓制串?dāng)_噪聲，同時(shí)也造成了有效信號(hào)的泄漏，殘差大，分離結(jié)果較差；Curvelet域指數(shù)衰減閾值法對深部弱信號(hào)保留較好，但是串?dāng)_噪聲強(qiáng)的部位不能夠被很好地壓制；本文提出的自適應(yīng)閾值法能夠有效的壓制串?dāng)_噪聲，更好地恢復(fù)和保護(hù)深層弱信號(hào)，分離結(jié)果好，殘差小.為了更加細(xì)致的對比分離效果，我們分別提取原始主炮數(shù)據(jù)(圖2a)、混疊數(shù)據(jù)(圖2c)和分離結(jié)果(圖5a—c)中第110和300道的數(shù)據(jù)，繪制如圖7所示的單道波形對比圖.從單道波形對比圖中可以發(fā)現(xiàn)，這三種方法都能去除混疊數(shù)據(jù)中的噪聲，分離的結(jié)果都與原始未混疊數(shù)據(jù)接近，但從圖中放大部位可以看出本文自適應(yīng)閾值法分離得到的數(shù)據(jù)與原始主炮數(shù)據(jù)更相近.圖8為從殘差圖5d—f中抽取的第110和300道數(shù)據(jù)，反映了分離誤差的大小，能夠進(jìn)一步展示圖7中的三種方法分離得到的結(jié)果與原始主炮數(shù)據(jù)的接近程度.圖8中的數(shù)值越小，證明分離結(jié)果與原始主炮數(shù)據(jù)越接近，分離效果越好.可以看出，自適應(yīng)閾值法分離得到的誤差比子帶一致閾值法和Curvelet域指數(shù)衰減閾值方法分離得到的誤差小很多，在4s以后的弱信號(hào)部位，自適應(yīng)閾值法得到的分離誤差幾乎為0，由此證明了本文自適應(yīng)閾值法具有良好地保護(hù)弱信號(hào)的能力.

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

(14)

2.2 SEAM模型數(shù)據(jù)

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

2.3 第一個(gè)實(shí)際數(shù)據(jù)

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

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

2.4 第二個(gè)實(shí)際數(shù)據(jù)

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

3 討論

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

圖12 本文提出的自適應(yīng)閾值法分離 SEAM模型數(shù)據(jù)的信噪比曲線 其中實(shí)線代表主炮的分離信噪比，虛線代表干擾炮的分離信噪比.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 結(jié)論

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

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

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

圖14 本文提出的自適應(yīng)閾值法分離第一個(gè)實(shí)際數(shù)據(jù)的分離結(jié)果和殘差 (a) 主炮分離結(jié)果; (b) 干擾炮分離結(jié)果; (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 本文提出的自適應(yīng)閾值法分離 第一個(gè)實(shí)際數(shù)據(jù)的信噪比曲線 其中實(shí)線代表主炮的分離信噪比，虛線代表干擾炮的分離信噪比.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對本文自適應(yīng)閾值法 分離效果的影響 其中實(shí)線代表主炮的分離信噪比，虛線代表干擾炮的分離信噪比.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 第二個(gè)實(shí)際數(shù)據(jù)的分離測試 (a) 混疊數(shù)據(jù); (b) 主炮分離結(jié)果; (c) 干擾炮分離結(jié)果.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 衰減系數(shù)α取不同值時(shí)，本文自適應(yīng)閾值法 得到的主炮分離信噪比曲線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 α