基于表層多次波的地震數(shù)據(jù)插值方法研究

郭書娟，方伍寶，徐兆濤

（1.中國石油化工股份有限公司石油物探技術研究院，江蘇南京211103；2.同濟大學海洋與地球科學學院，上海200092）

由于野外采集的條件限制，原始地震數(shù)據(jù)中往往會出現(xiàn)缺道、壞道和數(shù)據(jù)分布不規(guī)則等情況，從而嚴重影響地震資料處理中波動方程偏移等多道處理算法的精度。因此，如何根據(jù)已獲取的有用信息來恢復缺失的地震信息，成為地震資料處理中急需解決的一個問題。

為了解決缺失地震數(shù)據(jù)的重建問題，有關專家學者研究了多種地震道插值技術：①波動方程插值方法［1－2］，包 括 NMO 和 反 DMO 交替進 行 的 數(shù) 據(jù)規(guī)則化、AMO方法等；②空間域插值方法［3］，如t－x域預測誤差濾波插值等；③頻率域插值方法［4－6］，如f－k域道內(nèi)插值，f－x域預測濾波插值等；④基于Radon變換的插值方法［7－8］，如采用拋物線 Radon變換、雙曲線Radon變換等構(gòu)建插值算子，進行不規(guī)則采樣數(shù)據(jù)重建等；⑤基于稀疏反演的插值方法［9－11］，包括稀疏約束最小二乘反演、非二次模型權函數(shù)的稀疏反演傅里葉重建方法等。這些方法雖然算法不同，但主要都是利用缺失道的相鄰道信息，通過變換或插值外推來填補空白道。

地震數(shù)據(jù)中包含有一次反射波和多次反射波，而通常被看作干擾波的多次波數(shù)據(jù)中也蘊含有可利用的地震信息。Shan等［12－13］提出通過一次波和表層多次波互相關的方法生成準一次波；Berkhout等［14－15］和 Verschuur等［16］提出通過聚焦轉(zhuǎn)換的方法將高階多次波轉(zhuǎn)換為低階多次波，其中包括將一階多次波轉(zhuǎn)換為虛擬一次波；Curry等［17－18］利用非平穩(wěn)預測誤差濾波方法，將準一次波作為輸入訓練數(shù)據(jù)，然后將其用于數(shù)據(jù)插值。這些方法都是利用提取出來的多次波和一次波進行插值，但是實際處理中將多次波提取出來并不容易。

我們提出一種利用本來蘊含在表層多次波數(shù)據(jù)中的地震信息實現(xiàn)疊前炮集資料近道缺失數(shù)據(jù)插值的新方法，不用預先提取多次波，而是直接利用包含有表層多次波的原始數(shù)據(jù)，通過互相關構(gòu)建虛擬一次波，并在滑動時間－空間窗內(nèi)結(jié)合最小二乘匹配濾波（LMF）和均方根（RMS）振幅校正方法對擬地震道進行校正，而后用于地震數(shù)據(jù)插值重建。雖然這種方法構(gòu)建的準一次波有一些噪聲，但是經(jīng)過校正后用于彌補空缺道是完全可以勝任的。

1 擬地震道構(gòu)建

采用互相關方法將地震數(shù)據(jù)多次波路徑部分包含的地震信息提取出來，構(gòu)建“擬地震道”。

如圖1a所示，震源S點激發(fā)的波在地下反射點X反射后傳播至接收點A，得到一次波；此一次波在自由表面反射后繼續(xù)向下傳播，在地下反射點Y反射后傳播至接收點B，得到表層多次波。那么，A點處地震波的旅行時可表示為τSX＋τXA，B點處地震波的旅行時可表示為τSX＋τXA＋τAY＋τYB。在頻率域，一次波地震記錄A和表層多次波地震記錄B可以分別用（1）式和（2）式表示

式中：R是地下反射界面的反射系數(shù)。

通過一次波地震道A（圖1b）與多次波地震道B（圖1a）互相關，可以從多次波地震道B中將共同地震相iω（τSX＋τXA）消去，得到

（3）式即相當于將B道地震數(shù)據(jù)中包含的地震多次波信息－R3eiω（τAY＋τYB）提取了出來（圖1c）。

圖1 擬地震道構(gòu)建原理

以原始含有表層多次波的地震道為數(shù)據(jù)基礎，將所有炮的A道與B道地震數(shù)據(jù)互相關所得結(jié)果相疊加：

由此構(gòu)建出的地震道與實際由A點激發(fā)在地下反射點Y反射后從B點接收的地震道具有相似的運動學特征，包含擬一次波和擬多次波信息，通常被稱作“擬地震道”；實際接收點A點和B點則分別稱為“擬炮點”和“擬接收點”。

擬炮點A和擬接收點B可以定義在任意接收點位置，因此，只要地震數(shù)據(jù)中蘊含小炮檢距的多次波信息，就可以靈活地構(gòu)建出任意小炮檢距的擬地震道，其在原地震數(shù)據(jù)中通常是缺失的地震信息。若要構(gòu)建出大炮檢距擬地震道，通常需要地下存在特殊地質(zhì)構(gòu)造，如陡傾角反射層［13］。擬炮點和擬接收點位置選定的靈活性是構(gòu)建擬地震道用于地震數(shù)據(jù)插值的基礎。

2 擬地震道校正處理和插值

通過上述互相關方法構(gòu)建出的擬地震道與原始數(shù)據(jù)在動力學特征（如振幅和子波等）方面存在差異，不能直接用于原始地震數(shù)據(jù)的插值重建。為了解決這個問題，采用最小二乘匹配濾波（LMF）和均方根（RMS）振幅校正的方法，在疊前共炮集數(shù)據(jù)域滑動時間－空間窗內(nèi)對擬地震道數(shù)據(jù)進行相位及波形特征方面的匹配和振幅能量的校正。

如圖2，在原始記錄空白道附近開一個數(shù)據(jù)窗，分別選取原始記錄和擬地震記錄中與缺失道鄰近的幾道（如圖2a窗內(nèi)實線和圖2b窗內(nèi)相對應部分），基于LMF和RMS方法構(gòu)建校正算子，用于對圖2b數(shù)據(jù)窗內(nèi)與圖2a缺失道對應的擬地震道進行校正，校正后的擬地震記錄即可補充到原始記錄中的空白道處（如圖2a窗內(nèi)虛線部分）?？v向上根據(jù)淺、中、深層選取不同的時間窗（如圖2a和圖2b所示），時間窗重合的部分進行線性加權疊加。橫向上空間窗也順次向前推進，如圖2c和圖2d中實線窗向虛線窗推進。

為了避免數(shù)據(jù)跨度大影響匹配效果，數(shù)據(jù)窗的選擇要在可以接受的工作效率下盡量窄，但是也不能選取太少，比如臨近空白道的橫向1道或2道，這樣不利于數(shù)據(jù)的穩(wěn)定。

2.1 最小二乘匹配濾波（LMF）

分別在缺道的原始數(shù)據(jù)和構(gòu)建出來的擬地震記錄中某一個數(shù)據(jù)窗內(nèi)臨近缺失道的位置，各選擇出對應位置相同的N道數(shù)據(jù)y（t）和p（t），用于匹配算子的構(gòu)建。將選擇出來的N道原始數(shù)據(jù)表示為yi（t）（i＝1，2，…，N），擬地震數(shù)據(jù)表示為pi（t）（i＝1，2，…，N），i為道號，yi（t）和pi（t）具有相對應的炮檢點位置。

設計一個疊前匹配濾波算子mi（t）作用于地震道pi（t），使pi（t）經(jīng)過匹配濾波后逼近地震道yi（t）。假設疊前濾波器的實際輸出pi（t）＊mi（t）與期望輸出yi（t）的誤差為pi（t）＊mi（t）－yi（t），則總誤差能量為

圖2 滑動時間窗和滑動空間窗

理想狀態(tài)是使得E最小。根據(jù)最小二乘原理，令總誤差能量E對mi（t）的偏導數(shù)等于0，則經(jīng)過整理可以得到如下托布里茲矩陣方程：

式中：Rpp表示pi（t）的自相關函數(shù)矩陣；Ryp表示期望輸出道yi（t）和輸入道pi（t）的互相關函數(shù)向量；M是待求的匹配濾波算子向量。

求解（5）式所示的托布里茲矩陣方程，可以得到第i道的匹配濾波算子mi（t）。對于i＝1，2，…，N，可以得到N個算子。求取這N個算子的平均，即可得到我們所需的匹配算子m（t），即

將匹配濾波算子m（t）作用于數(shù)據(jù)窗內(nèi)缺失道對應的擬地震道，即與擬地震道相褶積，即可完成對該窗內(nèi)擬地震道的匹配濾波處理。

經(jīng)過匹配濾波處理后的擬地震道與原始地震記錄在相位及波形特征方面都很接近，但是振幅能量方面還有一些差異，需要進一步校正后才能用于原始地震數(shù)據(jù)缺失道的插值重建。

2.2 均方根振幅校正

利用均方根（RMS）振幅校正方法對經(jīng)過最小二乘匹配濾波后的擬地震道與原始記錄在振幅能量方面的差異進行校正。

分別計算上述數(shù)據(jù)窗內(nèi)選取出的N道原始數(shù)據(jù)yi（t）和經(jīng)匹配濾波后的擬地震道Pi（t）的平均RMS振幅ay和aP，即

式中：ayi和aPi分別表示單道yi（t）和Pi（t）的均方根振幅；k表示要用于RMS振幅校正的時間長度，即所在的時間窗的長度。ay／aP＝a即為我們所需對經(jīng)過LMF校正后的擬地震道再進行校正的均方根振幅校正因子。

經(jīng)過LMF和RMS振幅校正后的擬地震道在相位、波形特征以及振幅能量方面都得到了較好的匹配和校正，能夠適用于原始地震數(shù)據(jù)缺失道的插值重建。

3 模型測試與實際資料試處理

3.1 Sigsbee 2B模型數(shù)據(jù)測試

應用Sigsbee 2B模擬數(shù)據(jù)對本文提出的基于表層多次波數(shù)據(jù)的地震數(shù)據(jù)插值方法進行測試。Sigsbee 2B模型數(shù)據(jù)是Subsalt Multiples Attenuation And Reduction Team（SMAART）協(xié)會推出的用于表層多次波數(shù)據(jù)處理的標準模型數(shù)據(jù)。

圖3展示了應用本文方法構(gòu)建擬地震道彌補近炮檢距缺失數(shù)據(jù)的效果。為了更清晰地展示插值效果，只截取了標準模型數(shù)據(jù)記錄2.5～6.5s部分。圖3a為Sigsbee 2B模型數(shù)據(jù)含表層多次波的完整原始單炮記錄；圖3b為圖3a中炮檢距小于2 000×0.304 8m近炮檢距（27道）數(shù)據(jù)缺失的原始記錄；圖3c為用互相關方法生成的整張擬地震單炮記錄，缺失的近炮檢距數(shù)據(jù)已被構(gòu)建，但擬地震記錄中包含有噪聲，這是非成對低階和高階多次波互相關帶來的噪聲；圖3d是將圖3c擬地震記錄的前27道數(shù)據(jù)未經(jīng)校正直接用于圖3b缺道原始記錄插值的結(jié)果，可以看出插入的擬地震道和原始記錄的同相軸特征不相符；圖3e是應用本文方法校正后的擬地震道用于圖3b缺道原始記錄插值的結(jié)果。對比圖3e與圖3a可以看出本文方法對于缺失近炮檢距模擬數(shù)據(jù)的插值效果。

為了更好地展示本文方法的插值效果，將圖3e最終插值結(jié)果和圖3a原始記錄的前27道數(shù)據(jù)放大后進行對比，分別見圖4a和圖4b，兩者的誤差見圖4c?？梢钥闯?，兩者主要的同相軸特征基本一致，缺失的近炮檢距主要信息已得到了較好的補充；插值結(jié)果中除了有很小的能量差異之外，還殘留了一些擬地震道構(gòu)建過程中帶來的噪聲。

圖4 插值結(jié)果和原始數(shù)據(jù)對比

圖5a和5b分別給出了原始炮集和插值結(jié)果的頻譜，從中可以看出兩者的帶寬和主頻基本一致，插值結(jié)果中沒有出現(xiàn)假頻。

模型數(shù)據(jù)測試結(jié)果驗證了本文方法的正確性和有效性。

圖5 原始數(shù)據(jù)（a）和插值結(jié)果（b）頻譜分析

3.2 實際資料試處理

為了進一步測試本文方法的實用效果，對某油田的一塊實際地震資料進行了試處理。該區(qū)由于存在強波阻抗界面，極易產(chǎn)生多次波，具有利用多次波信息進行地震數(shù)據(jù)插值的數(shù)據(jù)基礎。

圖6a給出了缺失炮檢距小于1 125m近炮檢距數(shù)據(jù)的含多次波實際地震道集；圖6b為用互相關方法生成的擬地震數(shù)據(jù)道集，可以看出缺失的近炮檢距數(shù)據(jù)已在擬地震道集中構(gòu)建出來，其中包含了擬一次波和擬多次波信息。對圖6b中炮檢距小于1 125m的擬地震道進行LMF和RMS振幅校正后補充到圖6a原始數(shù)據(jù)中，結(jié)果見圖6c所示。由圖6c可見，淺層1.4s附近的同相軸是一次波信息，實際原始地震道集中缺失的近道地震數(shù)據(jù)得到了較好的補充，插值后的同相軸整體比較連貫。圖6c中2.6s左右的同相軸是多次波，缺失的近道信息也得到了補充，但插值部分的地震道能量稍弱，同相軸不是很連貫。這是因為在擬地震道構(gòu)建過程中，擬一次波是由一階多次波降階得到的，擬多次波則來自于更高階多次波，而高階多次波由于傳播時間長，路徑復雜，能量比較弱，致使構(gòu)建出的擬多次波信息能量較弱。在實際資料處理中，多次波部分通常不是最終偏移處理所需要的成像能量。

實際地震資料試處理的結(jié)果表明，應用本文方法可以將多次波作為有效信息來彌補實際原始地震資料中缺失的一次波信息，這進一步證明了本文方法的有效性和實用性。

圖6 某油田實際地震道集資料基于表層多次波的插值方法試處理效果

4 結(jié)束語

我們提出了一種基于表層多次波數(shù)據(jù)實現(xiàn)疊前地震資料近道數(shù)據(jù)插值重建的新方法。此方法利用原始數(shù)據(jù)多次波中蘊含的信息來彌補在采集中缺失的近炮檢距記錄，較適用于表層多次波比較發(fā)育地區(qū)的地震數(shù)據(jù)，為解決含有多次波并需要近炮檢距數(shù)據(jù)補充的地震數(shù)據(jù)處理問題提供了一個很好的思路。

基于表層多次波的插值新方法無需進行表層多次波和一次波提取，直接利用含多次波的原始資料構(gòu)建擬地震數(shù)據(jù)，這樣在構(gòu)建擬地震道的過程中會引入一些噪聲，但是并不會影響后續(xù)經(jīng)校正處理后用于地震數(shù)據(jù)的插值重建。當然，如果可以去除擬地震道中的噪聲，將會帶來更好的數(shù)據(jù)插值效果，這也是本文方法下一步改進和發(fā)展的方向。

［1］Ronen J.Wave－equation trace interpolation［J］.Geophysics，1987，52（7）：973－984

［2］Chemingui N，Biondi B.Handling the irregular geometry in wide azimuth surveys［J］.Expanded Abstracts of 69thAnnual Internat SEG Mtg，1999，32－35

［3］Claerbout J，Nichols D.Interpolation beyond aliasing by（t－x）domain in PEFs［J］.Expanded Abstracts of 53rdAnnual Internat SEG Mtg，1991，1－10

［4］國九英，周興元.F－K域等道距道內(nèi)插［J］.石油地球物理勘探，1996，31（2）：211－218 Guo J Y，Zhou X Y.Iso－interval trace interpolation in F－K domain［J］.Oil Geophysical Prospecting，1996，31（2）：211－218

［5］Spitz S.Seismic trace interpolation in the F－X domain［J］.Geophysics，1991，56（6）：785－794

［6］Ji J.Interpolation using prediction－error filter simulation（PEFs）［J］.Expanded Abstracts of 63rdAn－nual Internat SEG Mtg，1993，1170－1173

［7］Thorson J R，Claerbout J F.Velocity－stack and slant stack stochastic inversion［J］.Geophysics，1985，50（12）：2727－2741

［8］Kabir M M N，Verschuur D J.Restoration of missing offsets by parabolic Radon transform［J］.Geophysical Prospecting，1995，43（3）：347－368

［9］Liu B，Sacchi M.Minimum weighted norm interpolation of seismic data with adaptive weights［J］.Expanded Abstracts of 73rdAnnual Internat SEG Mtg，2001，1921－1924

［10］Wang Y H.Sparseness－constrained least－squares inversion：application to seismic wave reconstruction［J］.Geophysics，2003，68（5）：1633－1638

［11］Zwartjes P M，Gisolf A.Fourier reconstruction with sparse inversion［J］.Geophysical Prospecting，2007，55（2）：199－221

［12］Shan G J.Source－receiver migration of multiple reflections［J］.Expanded Abstracts of 73rdAnnual Internat SEG Mtg，2003，1008－1011

［13］Shan G J，Guitton A.Migration of surface－related multiples：tests on the Sigsbee 2Bdataset［J］.Ex－panded Abstracts of 74thAnnual Internat SEG Mtg，2004，1497－1500

［14］Berkhout A J，Verschuur D J.Transformation of multiples into primary reflections［J］.Expanded Abstracts of 73rdAnnual Internat SEG Mtg，2003，1925－1928

［15］Berkhout A J，Verschuur D J.Imaging of multiple reflections［J］.Geophysics，2006，71（4）：S1209－S1220

［16］Verschuur D J，Berkhout A J.Transforming multiples into primaries：experience with field data［J］.Expanded Abstracts of 75thAnnual Internat SEG Mtg，2005，2103－2106

［17］Curry W，Shan G J.Interpolation of near offsets using multiples and prediction－error filters［J］.Expanded Abstracts of 78thAnnual Internat SEG Mtg，2008，2421－2424

［18］Curry W，Shan G J.Interpolation of near offsets using multiples and prediction－error filters［J］.Geophysics，2010，75（6）：WB153－WB164