時空域與τ—p域預測反褶積壓制多次波方法對比研究

席自彬，麻志國，王海昆，徐 強

（中海油田服務股份有限公司物探事業(yè)部，天津 300451）

1 引言

多次波的產(chǎn)生需要良好的反射界面，只有在反射系數(shù)較大的反射界面上產(chǎn)生的多次反射，才能形成較強的多次波，并被檢波器記錄下來［1］。由于海洋地震資料采集特性，海水—空氣、海水—海底都是良好的反射界面，所以海洋地震資料中多次波比較發(fā)育。常見的海上地震資料中，既有在海水與海底之間多次震蕩產(chǎn)生的鳴震、交混回響等多次波，又有由于海水橫向深度變化劇烈（尤其深水）產(chǎn)生的繞射多次波，以及層間多次波等［2，3］。因此，對于海洋地震資料，關于多次波的去除是一項重要的工作。目前，去除多次波的方法有多種?；跒V波類的多次波去除方法［4，5］依據(jù)的是一次波與多次波信號特征上的差異，即周期性差異或者速度性差異；基于波動方程的多次波預測方法［6－8］是從波場的角度預測，并從原始地震記錄中減去多次波。

預測反褶積方法是工業(yè)生產(chǎn)中常用的多次波壓制方法，它是濾波類多次波去除方法的一種。該方法主要根據(jù)多次波在時間上呈現(xiàn)的規(guī)律性，不需要速度等輔助信息，同時預測步長又可以通過自相關或者步長掃描得到，兼具參數(shù)設置簡單和運行速度快的特點。預測反褶積方法自1954年Robinson首次提出以后［9］，經(jīng)歷了幾十年的發(fā)展與演進，已成為地震資料處理中最基本、最常規(guī)的處理手段［10－12］。唐博文等［13］等于2010年在頻率域設計提出一種新的譜模擬反褶積方法，該方法可以在不做多項式擬合時準確提取子波振幅譜。2015年，王元君等［14］在動態(tài)褶積模型基礎上，提出一種基于時頻域的動態(tài)反褶積方法，該方法是基于S變換的變Q值情形估算反射系數(shù)來提高地震資料信噪比。2021年，張聯(lián)海等［15］研究了基于深度卷積神經(jīng)網(wǎng)絡的稀疏反褶積方法，該方法可以不用多次設置參數(shù)來提高計算的速度和精度。上述研究方法主要是通過不同反褶積方法求取地震反射系數(shù)及提高地震分辨率，對反褶積方法壓制多次波沒有進行深入的探討。在預測反褶積壓制多次波研究方面，2015年，徐云霞等［16］研究了f—x域預測反褶積在共偏移距衰減多次波；2018年，張紅梅［17］將線性拉東域預測反褶積方法應用于海上地震資料處理；張向輝［18］將該方法推廣應用到陸上復雜地質(zhì)成像中，均取得了較好的多次波壓制效果；何林幫等［19］提出雙曲Radon域預測反褶積方法，該方法可以壓制速度逆轉的多次波及中長期多次波，缺點是局部存在損失有效信號現(xiàn)象。但是上述研究方法更偏向于實際資料處理，對壓制多次波時影響參數(shù)和保幅性等問題沒有進行深入的研究。時空域與τ—p域預測反褶積是目前工業(yè)生產(chǎn)中普遍使用且較為成熟的反褶積方法，相較其他方法，具有參數(shù)設置簡便、計算速度快及計算精度高的優(yōu)點。本文首先闡述了時空域預測反褶積和τ—p域預測反褶積方法的基本原理，分析了兩種方法的影響參數(shù)、頻譜特性，同時分析了兩種方法的相對保幅性質(zhì)，在實際海洋地震資料處理中分析兩種方法的優(yōu)缺點。大量實際資料證明，預測反褶積方法在消除海上地震資料中虛反射、交混回響、鳴震或者其他形式簡單的多次波方面有較好的效果。

2 預測反褶積方法原理

2.1 時空域預測反褶積基本原理

在實際地震勘探中，一次反射波對應著實際地層，在同一點不同的時間，一次反射波基本上是不同的，沒有規(guī)律性，即它在時間上不可預測。對于鳴震等多次波，它在地震剖面上出現(xiàn)的時間有規(guī)律性，即可以通過前面的波形預測出后面的波形。從地震記錄中減去多次波干擾，得到一次波反射，即為預測誤差。假設地震子波a（t）為最小相位子波，反射系數(shù)ξ（t）為白噪聲，則褶積模型定義為：

t＋l時刻的輸出為：

而預測誤差為：

從式（3）可以看出，式（2）中第二項為預測值x＾（t＋l）（多次波），第一項為預測誤差e（t＋l）（一次波）。因此，

若設a′（t）＝（a（0），a（1），…，a（l－1）），即取地震子波的前部，則式（4）變?yōu)椋?/p>

式（5）說明：預測誤差為反射系數(shù)和子波前部進行褶積處理，這也是預測反褶積方法能夠壓縮子波長度的本質(zhì)，從而達到提高縱向分辨率的目的。從子波角度看，預測反褶積相當于對子波進行截尾處理，截取子波為a′（t）＝（a（0），a（1），…，a（l－1））。

下面求取預測因子和反褶積因子。眾所周知，預測反褶積是維納濾波，即基于線性預測理論的最小平方濾波。假設b（t）＝（b（0），b（1），…，b（m））為預測因子，則在最小平方情況下，可推出求解預測因子的方程：

式（6）中，rxx為地震記錄的自相關；λ為白噪系數(shù)。計算出預測因子后，通過和地震記錄褶積得到預測值，再由式（3）得到去噪結果。

2.2 τ—p域預測反褶積基本原理

τ—p域預測反褶積也就是線性Rodon域預測反褶積，它的實現(xiàn)需要經(jīng)過三個步驟：首先是將時空域（t－x域）道集經(jīng)過線性Rodon變換轉換到τ—p域（Rodon域也稱慢度—截距時間域）；第二步，在τ—p域對每一個慢度道做預測反褶積，去除多次波；第三步，經(jīng)過線性Rodon反變換轉換到時空域，得到去除了多次波的地震記錄［20］。

由地震反射理論可知，

而射線參數(shù)p為：

所以由式（7）和式（8）可得：

即

式（7）～式（10）中：x為偏移距（m）；t為雙程旅行時（s）；t0為垂直反射時間（s）；v為介質(zhì)速度（m／s）；p為射線參數(shù)；τ為p＝0時的截距時間。通過式（10）進行時空域和τ—p域互相轉換。假設時空域的道集為f（x，t），F(xiàn)（τ，p）是時空域道集的Rodon變換，n為τ—p變換所參與運算的道數(shù)，則將時空域道集經(jīng)式（10）轉換到τ—p域，可得

然后在τ—p域進行預測反褶積，過程和時空域預測反褶積原理一樣，將周期性的多次波去除，得到τ—p域衰減多次波后的一次波記錄Fa（τ，p），然后由式（10）進行線性Rodon反變換，得到時空域衰減多次波后的一次波fa（x，t），即

由圖1可以看出，對于同一個p值，信號雖然可以在很多偏移距接收到，但卻對應著同一個入射角，因此多次波有著嚴格的周期性，多次波的周期不會受到偏移距變化的影響，這也是τ—p域多次波周期性好于時空域多次波周期性的根本原因。

圖2 主要展示道集在時空域和τ—p域之間的相互映射關系。由式（10）可以知道，直達波B（圖中經(jīng)過原點直線）轉換到τ—p域就是一個點B′。對于入射角小于臨界角的折射波A（圖中截距不為0的直線）和反射D，轉換到τ—p域為低P值區(qū)A′和D′。而對于寬角反射C，則映射到高p值區(qū)C′［21］。

圖2 時空域道集和τ—p域道集的相互映射關系Fig.2 Mapping relationship between temporal spatial gather andτ—pgather

3 影響參數(shù)分析

3.1 時空域預測反褶積影響參數(shù)

3.1.1 預測步長的影響

圖3為水平水層多次波模型。由圖3可知，在時空域，低一階次與高一階次反射波時距關系為：

圖3 水平水層多次波模型Fig.3 Multiple wave model of horizontal water layer

式（13）中，n為表面多次波的階數(shù)；h為激發(fā)點處的法向深度（m）。當x＝0時，有Δtn＝。因此在時空域（t－x域），只有在零炮檢距記錄即垂直入射（x＝0）的情況下，多次波才有嚴格的周期性。隨著多次波偏移距的增大和深度的加深，多次波的時差會減小，即多次波在時空域是時空變的。因此，在實際處理中選取預測步長時，如果預測步長取值過小，不能準確預測并有效去除多次波；而且不能以近道的自相關次極值為標準，要以遠道的自相關次極值為參照，才能有效地去除整個道集的多次波。

3.1.2 算子長度的影響

算子長度較易求取，在式（6）中，只要滿足自相關中包括多次波和一次波的主要能量，算子長度就能滿足預測反褶積的要求。但是，算子長度必須大于多次波的一個周期，如果預測算子長度太小，可能會帶來假的能量。在實際處理中，可以適當?shù)亟o大算子長度，對結果沒有影響。

3.2 τ—p域預測反褶積影響參數(shù)

由于τ—p域預測反褶積包括τ—p正、反變換和預測反褶積兩大步，因此影響因素也要充分考慮這兩方面參數(shù)的影響。對于理論模型，τ—p域預測反褶積可以很理想地去除多次波，但是受采集方式和海底復雜情況的影響，實際情況非常復雜。

3.2.1 算子本身的影響

從式（7）和式（11）可以看出，由于公式近似和數(shù)值離散造成空間采樣率不足，會在τ—p正、反變換過程中帶來假頻。為有效防止空間假頻，已知道間距Δx時，可知最大波數(shù)k為，進一步推導最大慢度pmax為：

為有效防止慢度p上的假頻，在已知偏移距范圍xr和最大頻率fmax的情況下，必須滿足相鄰p道間的時差小于一個周期（T），所以有

為有效防止截距時間τ上的假頻，可以與時間采樣定理一致，須滿足：

3.2.2 參數(shù)pmax的影響

一個入射角對應一個射線參數(shù)p，同時對應著τ—p域一道；pmax既對應資料中最大傾角，也對應τ—p域的最大道數(shù)。若pmax過大，即在時空域道集上并不存在與之對應的同相軸，會在τ—p域相應道上產(chǎn)生噪音。產(chǎn)生的噪音也會參與預測反褶積運算，然后通過τ—p反變換轉換到時空域，會在τ—p反變換后的時空域道集上帶來噪音，造成道集信噪比下降。因此，pmax應選擇地震資料中處理感興趣的傾角范圍，不應過大。

3.2.3 其他因素的影響

1）直達波以及折射波對τ—p域預測反褶積有影響，由式（10）知，直達波以及折射波轉換到τ—p域是一個點，不能準確地再轉換回時空域，所以在做τ—p域預測反褶積時，需要切除直達波和折射波，或做相應特殊處理；

2）關于τ—p域采樣psample，當采樣不足時，會產(chǎn)生一些隨機噪音，但是過采樣不會產(chǎn)生任何影響，也不會提高道集的質(zhì)量；

3）關于最大頻率Fmax，主要用于計算p的道數(shù)，必須給足，但是大小不會影響道集的質(zhì)量，過大會影響計算的效率。

4 保幅性分析

假設時空域道集上兩道地震記錄為x1（t）和x2（t），首先經(jīng)過地表一致性處理，消除道集的子波橫向性差異，在此基礎上，假設兩個地震記錄振幅在時窗t1～t2時刻內(nèi)任意采樣點滿足x2（t）＝kx1（t），其中k為任意常數(shù)，預測反褶積后，相應的地震記錄為y1（t）和y2（t），b1（t）和b2（t）為相應的預測反褶積算子，則最小平方反褶積濾波方程為式（6）。假設rxx為地震記錄的自相關，由于rxx（τ）＝，可得rx2x2（τ）＝k2rx1x1（τ），其中，τ＝0，1，…，m。因此，式（6）左邊自相關矩陣滿足R2＝k2R1。當l為預測步長時，預測反褶積可以將rxx（l＋τ）當作rdw（τ）。因此可得，rx2x2（l＋τ）＝k2rx1x1（l＋τ），所以可得b1（t）＝b2（t），則預測值滿足y′2（t）＝ky′1（t），最終可求得反褶積后的振幅值y2（t）＝ky1（t），即預測反褶積前后的振幅比例沒有變化，所以單道的預測反褶積是相對保幅的。

對于多道預測反褶積來說相對簡單，由于相鄰多道采用相同的反褶積因子，即b1（t）＝b2（t），又因為x2（t）＝kx1（t），由反褶積運算可得y2（t）＝ky1（t），所以多道預測反褶積也是相對保幅的。在預測反褶積運算中，使用較多的是多道預測反褶積，因為多道預測反褶積統(tǒng)計誤差小，子波提取比單道預測反褶積穩(wěn)定，所以相對保幅性要好于單道預測反褶積［22，23］。本文采用的就是多道預測反褶積方法。

5 實際海上資料處理應用

本文選擇渤海某工區(qū)實際資料，首先對τ—p域預測反褶積影響參數(shù)進行數(shù)據(jù)測試，然后分別在時空域和τ—p域做預測反褶積，對比反褶積效果。

5.1 τ—p域預測反褶積影響參數(shù)數(shù)據(jù)測試

5.1.1 算子本身的影響

由于計算公式采取近似和數(shù)值離散造成采樣率不足，會在τ—p正、反變換過程中帶來假頻。如圖4所示，圖4（a）為原始炮集，4（b）為經(jīng)過τ—p正、反變換后的炮集，圖4（c）為τ—p正、反變換帶來的假頻。

圖4 τ—p正、反變換前后效果Fig.4 Effect before and afterτ—pForward and reverse transform

5.1.2 參數(shù)pmax的影響

如圖5所示，最大道數(shù)pmax過大時，會在τ—p正變換后產(chǎn)生噪音。從圖6可以看出，產(chǎn)生的噪音也會參與預測反褶積運算，然后通過τ—p反變換轉換到時空域，會在時空域道集上形成噪音，造成道集信噪比下降。因此，對于pmax的大小應選擇地震資料中處理感興趣的傾角范圍，不應過大。由式（14）可知，pmax取值應小于。

圖5 τ—p域不同pmax大小效果Fig.5 Different pmaxeffect inτ—pdomain

圖6 不同pmax值時τ—p域預測反褶積效果Fig.6 Predictive deconvolution effect with different pmaxinτ—pdomain

5.2 兩種預測反褶積方法實際數(shù)據(jù)測試對比

圖7 為原始炮集，圖8為時空域原始炮集的自相關和預測反褶積后的自相關，選取的時窗為1600～3200ms。從自相關上可以看到，時空域近偏移距多次波去除效果較好，但是在遠偏移距上還有少許殘留的多次波。圖9為原始炮集τ—p域自相關和τ—p域預測反褶積后的自相關，選取時窗和時空域相同。從圖9可以看出，τ—p域預測反褶積多次波去除效果比時空域好，但是帶來了一些低頻噪音。圖10為預測反褶積前（紅色）、時空域預測反褶積（綠色）、τ—p域預測反褶積（藍色）三者頻譜對比，選取的時窗為500～2000ms。從圖10可以看出，τ—p域預測反褶積優(yōu)勢頻帶更寬，高頻成分、低頻成分更豐富。圖11為預測反褶積前（紅色）、時空域預測反褶積（綠色）、τ—p域預測反褶積（藍色）三者振幅譜對比。從圖11可以看出，預測反褶積方法相對保幅，時空域預測反褶積振幅值約是τ—p域預測反褶積振幅值的1.2倍，因此，時空域預測反褶積相對保幅效果更好。由于τ—p域預測反褶積的實現(xiàn)要經(jīng)過τ—p域正、反變換和預測反褶積兩步，因此它的計算效率要低于時空域預測反褶積。實驗證明，τ—p域預測反褶積計算時間是時空域預測反褶積的1.5倍。

圖7 原始炮集Fig.7 Original shot gather

圖8 時空域預測反褶積前后炮集自相關Fig.8 Shot gather auto－correlation before and after temporal spatial predictive deconvolution

圖9 τ—p域預測反褶積前后炮集自相關Fig.9 Shot gather auto－correlation before and afterτ—pdomain predictive deconvolution

圖10 不同反褶積方法頻率譜Fig.10 Frequency spectra of different deconvolution methods

圖11 不同反褶積方法振幅譜Fig.11 Amplitude spectra of different deconvolution methods

在兩種預測反褶積均采用最優(yōu)參數(shù)的情況下，對比兩種預測反褶積效果。從圖12疊加剖面可以看出，時空域預測反褶積多次波有少許殘留（圖中黑色箭頭位置），τ—p域預測反褶積對多次波壓制更干凈。但在相對保幅性方面，時空域預測反褶積疊加剖面同向軸更實，連續(xù)性更好（圖中綠色箭頭位置），保幅效果要優(yōu)于τ—p域預測反褶積。

圖12 兩種預測反褶積疊加效果對比Fig.12 Comparison of two predicted deconvolution stack effects

6 結論

時空域和τ—p域預測反褶積在實際資料處理時存在明顯差別：時空域預測反褶積主要受預測步長和算子長度兩項參數(shù)影響較大，在實際處理中選取預測步長過小時，不能準確預測并有效去除多次波。算子長度也必須大于多次波的一個周期，預測算子長度太小，會帶來假的能量。相比而言，多次波在τ—p域有更明顯的規(guī)律性，但是由于該方法本身在計算時存在公式近似和數(shù)值離散，會在τ—p正、反變換過程中因空間采用率不足帶來假頻。并且pmax參數(shù)的選擇至關重要，需要經(jīng)過反復實驗論證，數(shù)值太小會造成有效波缺失，數(shù)值太大會帶來假的噪音。同時，τ—p域預測反褶積還需要關注直達波、折射波、τ—p域采樣間隔、最大頻率等參數(shù)的影響，因此在實際處理時更加需要精細對待。通過本文的對比研究，可以得到以下幾個結論：

1）在時空域，多次波只有在垂直入射即零炮檢距情況下才能保持其周期性，但是在τ—p域，多次波具有嚴格的周期性，不隨炮檢距變化而變化，所以在τ—p域去除多次波效果更理想；

2）τ—p域預測反褶積去除多次波效果優(yōu)于時空域預測反褶積，同時可以拓寬優(yōu)勢頻帶，豐富低頻、高頻信息；

3）預測反褶積方法是相對保幅的，時空域預測反褶積相對保幅性效果要優(yōu)于τ—p域預測反褶積，同時時空域預測反褶積的計算效率是τ—p域預測反褶積的1.5倍。