吳庚， 劉財(cái)， 劉殿秘， 劉洋*， 鄭植升

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130026 2 中國(guó)石油吉林油田勘探開(kāi)發(fā)研究院， 松原 138003

0 引言

在地震勘探中，采集地震數(shù)據(jù)的數(shù)量和質(zhì)量是獲得高精度地下參數(shù)反演和解釋的基礎(chǔ).地震數(shù)據(jù)在時(shí)間方向上總是充分采樣的，而由于野外采集的經(jīng)費(fèi)限制以及環(huán)境因素的影響，如陸上勘探地表障礙物的存在以及海上勘探采集時(shí)電纜的羽漂等情況，實(shí)際數(shù)據(jù)在空間方向上是非規(guī)則分布且稀疏的.通過(guò)對(duì)地震數(shù)據(jù)進(jìn)行插值重建，可以使其包含的信息更加充分地反映地下地質(zhì)體的地球物理特征，提高對(duì)復(fù)雜地質(zhì)構(gòu)造和儲(chǔ)層進(jìn)行精細(xì)刻畫(huà)的能力，因此往往需要在地震處理流程的早期階段對(duì)采集到的數(shù)據(jù)進(jìn)行插值或規(guī)則化處理.

缺失地震數(shù)據(jù)插值方法按照地震數(shù)據(jù)的缺失類(lèi)型不同可以分為三類(lèi)：隨機(jī)缺失插值方法、規(guī)則缺失插值方法和連續(xù)缺失插值方法.由于地理環(huán)境因素的影響以及缺道壞道現(xiàn)象的存在，地震數(shù)據(jù)總是存在著隨機(jī)缺失的問(wèn)題.Curry(2003)提出了非平穩(wěn)多尺度預(yù)測(cè)誤差濾波器(PEF)對(duì)隨機(jī)缺失數(shù)據(jù)進(jìn)行插值.Herrmann等(2008)提出了curvelet變換方法對(duì)隨機(jī)缺失地震數(shù)據(jù)進(jìn)行重組.劉洋等(2018)提出了Bregman整形迭代插值方法，該方法能夠利用有效的稀疏變換表征地震波場(chǎng)，迭代求解約束優(yōu)化問(wèn)題，對(duì)隨機(jī)缺失地震道進(jìn)行重建.Liu等(2019)提出了三維t-x-y域多尺度多方向自適應(yīng)PEF對(duì)隨機(jī)缺失數(shù)據(jù)進(jìn)行插值.陳思遠(yuǎn)等(2021)提出了基于頻率-結(jié)構(gòu)雙重加權(quán)閾值的隨機(jī)缺失地震數(shù)據(jù)插值方法，保證了高低頻信息的有效恢復(fù).李昂等(2021)在OVT(Offset Vector Tile)域使用匹配追蹤算法對(duì)五維地震數(shù)據(jù)進(jìn)行插值.賀月等(2021)提出一種與字典學(xué)習(xí)相結(jié)合的凸集投影算法對(duì)地震數(shù)據(jù)插值的同時(shí)去噪.董烈乾等(2022)基于信噪比質(zhì)控因子提出了優(yōu)化迭代反演非規(guī)則數(shù)據(jù)重構(gòu)方法，對(duì)三維隨機(jī)缺失數(shù)據(jù)進(jìn)行恢復(fù).

基于成本考慮，地震數(shù)據(jù)在道間距或炮間距方向往往采樣不足，可以通過(guò)規(guī)則地震數(shù)據(jù)缺失插值方法進(jìn)行處理.Spitz(1991)提出了頻率域預(yù)測(cè)濾波器(PF)，并用于地震道間插值.Porsani(1999)通過(guò)引入半步PF改進(jìn)了Spitz的方法，實(shí)現(xiàn)更加有效的地震道加密.Wang(2002)進(jìn)一步將預(yù)測(cè)插值從f-x域擴(kuò)展到了f-x-y域，并對(duì)三維地震數(shù)據(jù)進(jìn)行道間插值.Claerbout(1992)首次在地震道間插值中使用t-x域PEF.Crawley等(1999)設(shè)計(jì)了一種平滑變化的PEFs，通過(guò)“階梯濾波器”控制濾波器的平滑來(lái)重建地震道間數(shù)據(jù).Liu和Fomel(2011)提出t-x域正則化非平穩(wěn)自回歸(RNA)方法對(duì)規(guī)則和不規(guī)則缺失的地震數(shù)據(jù)進(jìn)行插值，其方法主要使用整形正則化來(lái)控制自適應(yīng)PEF的平滑性.Jia和Ma(2017)首次把機(jī)器學(xué)習(xí)應(yīng)用到地震數(shù)據(jù)道間插值中，一定程度上克服了傳統(tǒng)方法的物理假設(shè)條件.Wang等(2019)使用8層殘差神經(jīng)網(wǎng)絡(luò)對(duì)規(guī)則缺失地震數(shù)據(jù)進(jìn)行插值.Shao等(2021)使用Radon域地震干涉插值方法對(duì)規(guī)則數(shù)據(jù)進(jìn)行重建.Shao和Wang(2022)使用稀疏Radon變換和動(dòng)態(tài)標(biāo)記函數(shù)對(duì)地震數(shù)據(jù)進(jìn)行加密.

野外觀測(cè)系統(tǒng)布置時(shí)，可能會(huì)遇到無(wú)法布置炮點(diǎn)或者無(wú)法埋置檢波器的地段，如湖泊、工礦場(chǎng)和城鎮(zhèn)等.這種情況下，觀測(cè)系統(tǒng)中的測(cè)線連續(xù)性會(huì)受到影響，造成地震反射剖面間斷現(xiàn)象，形成連續(xù)缺失，盡管實(shí)際施工中地震數(shù)據(jù)采集的變觀方法能夠一定程度解決這種問(wèn)題，但是往往增加施工的復(fù)雜性.另外，由于在野外實(shí)際施工時(shí)，炮點(diǎn)位置附近無(wú)法布置檢波器，總是會(huì)造成近炮檢距缺失現(xiàn)象，這種情況在海上數(shù)據(jù)中尤為嚴(yán)重，影響后續(xù)的數(shù)據(jù)處理流程.因此，提出有效的連續(xù)缺失數(shù)據(jù)插值方法具有重要的實(shí)際意義.Trad等(2002)使用高分辨率時(shí)變Radon變換對(duì)近炮檢距缺失數(shù)據(jù)進(jìn)行重建.Naghizadeh和Sacchi(2007)提出了多步自回歸地震數(shù)據(jù)插值方法，對(duì)規(guī)則、不規(guī)則和連續(xù)缺失數(shù)據(jù)進(jìn)行重建.Van Groenestijn和Verschuur(2009)提出了稀疏反演一次波估計(jì)方法，并應(yīng)用于地震近炮檢距缺失數(shù)據(jù)上，克服了近偏移距數(shù)據(jù)缺失的影響.Wang等(2009)使用干涉插值方法對(duì)近炮檢距缺失數(shù)據(jù)進(jìn)行重建.Curry和Shan(2010)使用多重頻率域PEFs對(duì)近炮檢距缺失數(shù)據(jù)進(jìn)行插值.Baumstein和Neelamani(2010)采用統(tǒng)計(jì)學(xué)和物理意義相結(jié)合的限定條件，將多種類(lèi)型的先驗(yàn)信息和描述不同類(lèi)型數(shù)據(jù)之間映射的方法相結(jié)合對(duì)連續(xù)缺失數(shù)據(jù)進(jìn)行插值.Wang等(2010)提出了以模型為基礎(chǔ)的干涉插值方法，對(duì)近炮檢距缺失數(shù)據(jù)進(jìn)行重建.Naghizadeh和Sacchi(2010)通過(guò)定義峰值的位置估計(jì)數(shù)據(jù)的空間譜，并使用自回歸方法對(duì)規(guī)則和連續(xù)缺失數(shù)據(jù)進(jìn)行插值.Naghizadeh(2012)通過(guò)定義稀疏分布的主傾角在頻率波數(shù)域?qū)B續(xù)缺失數(shù)據(jù)進(jìn)行插值.Zu等(2016)使用傾斜條件限定地震數(shù)據(jù)插值反問(wèn)題的病態(tài)方程，對(duì)連續(xù)缺失數(shù)據(jù)進(jìn)行重建.Wang等(2022)基于深度學(xué)習(xí)方法對(duì)規(guī)則炮缺失數(shù)據(jù)和近偏移距缺失數(shù)據(jù)進(jìn)行高精度重建.

以往的方法雖然在缺失地震道插值中有著較好的插值精度，但往往需要迭代計(jì)算，會(huì)導(dǎo)致增加內(nèi)存占用和計(jì)算時(shí)間.為了解決這個(gè)問(wèn)題，在Fomel和Claerbout(2016)提出的非迭代流式預(yù)測(cè)誤差濾波器(SPEF)的基礎(chǔ)上，Liu等(2022)和Zheng等(2022)分別在時(shí)間域和頻率域提出流式預(yù)測(cè)濾波(SPF)方法進(jìn)行地震道的插值重建.但這兩種方法在大空缺地震道位置存在著濾波器系數(shù)無(wú)法更新的情況，難以處理連續(xù)缺失的地震數(shù)據(jù)插值問(wèn)題.基于以上問(wèn)題，本文提出了時(shí)間域非迭代的高階流式預(yù)測(cè)濾波(HOSPF)方法對(duì)連續(xù)缺失地震數(shù)據(jù)進(jìn)行插值，提出的高階限定條件能夠增加局部平滑性，改善以往SPF和SPEF在插值過(guò)程中濾波器難以更新的情況，從而提高連續(xù)缺失地震數(shù)據(jù)的插值精度.同時(shí)設(shè)計(jì)了與時(shí)間域高階SPF相匹配的空間非因果預(yù)測(cè)濾波器和蛇形插值路徑來(lái)減小插值誤差.通過(guò)對(duì)比工業(yè)主流的傅里葉凸集投影(POCS)迭代插值方法，理論模型和實(shí)際數(shù)據(jù)測(cè)試結(jié)果表明，本文提出的方法能夠在節(jié)省計(jì)算資源的同時(shí)保證連續(xù)缺失地震數(shù)據(jù)具有更高的插值精度，更加適合“兩寬一高”地震勘探技術(shù)的需求.

1 基本理論

預(yù)測(cè)濾波是地震數(shù)據(jù)處理的重要理論之一，可以有效解決地震數(shù)據(jù)去噪和插值問(wèn)題，其具體形式包含預(yù)測(cè)誤差濾波器(PEF)和預(yù)測(cè)濾波器(PF)，PEF和相應(yīng)的PF具有不同的濾波器結(jié)構(gòu)和系數(shù)，PEF涉及到預(yù)測(cè)數(shù)據(jù)點(diǎn)時(shí)間方向上的因果預(yù)測(cè)系數(shù)，而PF只有空間方向的濾波器系數(shù)，同時(shí)，PEF的預(yù)測(cè)結(jié)果是殘差，而PF的預(yù)測(cè)結(jié)果是數(shù)據(jù)本身.由于地震數(shù)據(jù)缺失是整道缺失，缺失數(shù)據(jù)在其對(duì)應(yīng)的時(shí)間方向上沒(méi)有預(yù)測(cè)關(guān)系，因此選用僅表征空間預(yù)測(cè)能力的PF在地震數(shù)據(jù)插值和抗噪聲干擾能力上會(huì)優(yōu)于PEF，因此本文選擇PF進(jìn)行地震數(shù)據(jù)插值.Claerbout(1992)提出PF對(duì)缺失數(shù)據(jù)進(jìn)行插值的兩個(gè)步驟.第一步從采集到的已知數(shù)據(jù)計(jì)算濾波器系數(shù)，第二步用計(jì)算的濾波器系數(shù)進(jìn)行缺失數(shù)據(jù)插值.大多數(shù)的自適應(yīng)PF都需要迭代計(jì)算來(lái)處理時(shí)空變化的地震同相軸信息，但是迭代方法增加了計(jì)算成本和存儲(chǔ)量，因此本文提出一種擴(kuò)展的非迭代時(shí)間域自適應(yīng)PF插值方法.

1.1 高階流式預(yù)測(cè)濾波系數(shù)解析估計(jì)

在二維情況下，自適應(yīng)PF能夠表示為：

(1)

濾波器系數(shù)可以通過(guò)求解如下最小二乘問(wèn)題得到：

(2)

(3)

(4)

其中，εt和εx分別是沿著時(shí)間和空間方向的標(biāo)量正則化參數(shù).

公式(4)對(duì)應(yīng)的最小二乘解形式為：

a(t,x)=[d(t,x)d(t,x)T+4ε2I]-1[d(t,x)d(t,x)

(5)

其中：

(6)

(7)

I表示單位矩陣.

公式(5)中逆矩陣可以通過(guò)迭代方法進(jìn)行求解，比如共軛梯度法等，但迭代方法會(huì)增加計(jì)算成本和內(nèi)存占用.由于該矩陣具有特殊的結(jié)構(gòu)，本文使用Sherman-Morrison方法(Hager，1989)直接求解形如(L-MN)-1矩陣的逆，其中矩陣M和N分別代表了列向量和行向量.該矩陣的逆可以表示為：

(L-MN)-1=L-1+L-1M(I-NL-1M)-1NL-1,

(8)

其中，L和I-NL-1M是可逆的.在本文中，上述公式中L、M和N分別代表了4ε2I，-d(t,x)和d(t,x)T.因此，公式(5)可以寫(xiě)成解析解的形式：

(9)

由于公式(9)中當(dāng)前系數(shù)a(t,x)依賴于處理路徑中前面相鄰的濾波器系數(shù)和數(shù)據(jù)值，因此在缺失數(shù)據(jù)中僅需要保證初始位置具有數(shù)據(jù)即可持續(xù)更新濾波器的時(shí)空變系數(shù)，初始的濾波器一般賦值為零，對(duì)于單邊接收的近炮檢距缺失數(shù)據(jù)，可以使用數(shù)據(jù)空間方向鏡像的方式進(jìn)行處理，從而保證HOSPF系數(shù)初始化的準(zhǔn)確性.由于使用了更多時(shí)間和空間方向的已知濾波器進(jìn)行約束，高階流式預(yù)測(cè)濾波器能夠更好地估計(jì)連續(xù)缺失地震數(shù)據(jù)位置的濾波器系數(shù).

預(yù)測(cè)濾波器具有空間因果和空間非因果兩種結(jié)構(gòu)，空間非因果的PF使用了更多目標(biāo)點(diǎn)附近的數(shù)據(jù)，這會(huì)比空間因果的濾波器有著更精準(zhǔn)的插值結(jié)果.因此，本文選用空間非因果流式濾波器，以20個(gè)濾波器系數(shù)的空間非因果PF為例，其中濾波器大小為5(時(shí)間)×4(空間)，表達(dá)式為：

(10)

其中，·表示零值.

本文方法可以直接擴(kuò)展為三維結(jié)構(gòu)，相比于上述二維情況，僅需要額外增加一個(gè)空間方向上的約束.在三維情況下，增加沿著第二個(gè)空間軸y方向上的約束條件，因此三維情況下公式(4)可以寫(xiě)成：

(11)

公式(11)的最小二乘解為：

(12)

其中：

(13)

(14)

對(duì)比二維公式(6)，可以發(fā)現(xiàn)僅僅增加了一個(gè)y軸方向上的濾波器系數(shù)和標(biāo)量正則化項(xiàng).

正則化參數(shù)εt、εx和εy一般選取在待插值數(shù)據(jù)的最大值和最小值之間，并根據(jù)不同方向平滑程度的不同對(duì)三個(gè)正則化參數(shù)進(jìn)行不同的調(diào)整.因此為了方便參數(shù)的設(shè)置，在進(jìn)行插值之前，可以先對(duì)數(shù)據(jù)進(jìn)行歸一化處理.

1.2 非迭代缺失地震數(shù)據(jù)快速重建

在流式預(yù)測(cè)濾波器估計(jì)的同時(shí)，可以利用計(jì)算好的濾波器系數(shù)直接對(duì)缺失的地震道進(jìn)行插值重建.以二維情況為例，對(duì)于HOSPF的殘差r(t,x)，可以表示為：

r(t,x)=d(t,x)-d(t,x)Ta(t,x)

(15)

通過(guò)對(duì)公式(15)進(jìn)行簡(jiǎn)單的數(shù)學(xué)變換，可以得到：

(16)

(17)

當(dāng)使用估計(jì)出的HOSPF系數(shù)對(duì)缺失數(shù)據(jù)進(jìn)行插值時(shí)，r(t,x)取0，則公式(17)進(jìn)一步簡(jiǎn)化為：

(18)

由于r(t,x)取0，公式(16)可簡(jiǎn)寫(xiě)為：

(19)

本文方法需要人為輸入的參數(shù)有濾波器大小和標(biāo)量正則化參數(shù)εn.通常根據(jù)地震同相軸的傾斜程度選擇濾波器大小.εn的選取范圍一般在輸入數(shù)據(jù)最大值和最小值之間，為了方便參數(shù)的選取，可以對(duì)數(shù)據(jù)振幅進(jìn)行歸一化處理.

為了進(jìn)一步提升插值精度，參考Zheng等(2022)的思路，針對(duì)HOSPF設(shè)計(jì)相匹配的時(shí)間-空間域插值路徑.以二維情況為例，插值路徑通常為圖1a所示的沿著空間方向逐行進(jìn)行插值，但這會(huì)導(dǎo)致每一行開(kāi)始時(shí)，濾波器系數(shù)需要重新進(jìn)行初始化計(jì)算，使得濾波器無(wú)法在整個(gè)數(shù)據(jù)中連續(xù)更新.針對(duì)這種情況，在二維時(shí)間-空間域內(nèi)設(shè)計(jì)如圖1b所示的蛇形插值路徑.在逐行插值的基礎(chǔ)上，相鄰行采用空間反方向?qū)?shù)據(jù)進(jìn)行濾波器系數(shù)計(jì)算和缺失數(shù)據(jù)重建，從而使得濾波器系數(shù)不需要對(duì)每行數(shù)據(jù)進(jìn)行初始化操作，進(jìn)一步提高HOSPF在數(shù)據(jù)邊緣位置的插值精度.與Zheng等(2022)方法不同，本文的濾波器是在時(shí)間-空間平面沿著蛇形插值路徑進(jìn)行移動(dòng)，當(dāng)一個(gè)時(shí)空平面計(jì)算完后，移動(dòng)到下一個(gè)時(shí)空平面上，而非在空間x-y平面沿著蛇形路徑移動(dòng).

圖1 不同模式的插值路徑(a) 逐行插值路徑； (b) 蛇形插值路徑.Fig.1 Different strategies for interpolation path(a) The row-by-row interpolation path；(b) The snaky interpolation path.

2 理論模型測(cè)試

2.1 二維理論模型測(cè)試

本文將提出的HOSPF插值方法與工業(yè)常用的傅里葉凸集投影(POCS)(Abma and Kabir，2006)方法進(jìn)行對(duì)比.為了定量地評(píng)價(jià)理論模型插值結(jié)果的精度，本文選取全局信噪比作為衡量標(biāo)準(zhǔn)(劉洋等，2017)：

(20)

首先建立如圖2a所示的理論模型，原始完整采樣的數(shù)據(jù)包含了多條雙曲同相軸，并且同相軸的振幅具有時(shí)空變化特征.首先隨機(jī)移除40%地震道，并在此基礎(chǔ)上，為了重點(diǎn)測(cè)試本文方法在連續(xù)缺失數(shù)據(jù)上的處理效果，又在第100道位置的附近連續(xù)移除了25道數(shù)據(jù)，其信噪比為2.72 dB.圖3分別展示了二維傅里葉POCS、低階流式預(yù)測(cè)濾波和本文提出的HOSPF方法的重建結(jié)果和插值誤差，為了保證處理效果，傅里葉POCS方法中使用了6×3個(gè)重疊窗，每個(gè)重疊窗大小為150(時(shí)間方向)×150(空間方向)，經(jīng)過(guò)測(cè)試迭代次數(shù)為575次可以獲得收斂的插值結(jié)果，如圖3a所示.傅里葉POCS方法能夠比較有效地恢復(fù)隨機(jī)缺失位置的地震道，同相軸的連續(xù)性得到提升，但是在連續(xù)缺失位置上效果并不理想，產(chǎn)生了較大的插值誤差(圖3b).低階流式預(yù)測(cè)濾波器方法同樣無(wú)法處理大面積連續(xù)缺失數(shù)據(jù)，同時(shí)在隨機(jī)缺失數(shù)據(jù)位置也產(chǎn)生了一些空間假頻.本文提出的HOSPF的參數(shù)中εt和εx均為0.5，濾波器大小為7(時(shí)間方向)×45(空間方向)，從插值結(jié)果(圖3c)和插值誤差(圖3d)可以看到，本文提出的方法在連續(xù)缺失數(shù)據(jù)上有著更合理的插值精度，只是在同相軸傾角比較大的位置準(zhǔn)確性有所降低，可以通過(guò)進(jìn)一步提高濾波器的階數(shù)來(lái)改善處理效果，但是同時(shí)也會(huì)使計(jì)算效率有所降低.傅里葉POCS、低階流式預(yù)測(cè)濾波器和本文方法插值結(jié)果的信噪比分別為11.25 dB、10.38 dB和18.61 dB，同時(shí)這三種方法在使用相同單核CPU的計(jì)算時(shí)間分別為23.65 s、1.93 s和2.30 s，可見(jiàn)在保持合理插值結(jié)果的同時(shí)，本文提出的方法比傅里葉POCS具有更快的計(jì)算速度，同時(shí)，相比于低階算法并沒(méi)有增加太多的計(jì)算時(shí)間.

圖2 二維理論模型測(cè)試(a) 原始完整采樣數(shù)據(jù)； (b) 包含40%隨機(jī)缺失和連續(xù)缺失的數(shù)據(jù).Fig.2 2D synthetic data test(a) Original well-sampled data； (b) Data with 40% randomly removed and continuously missing traces.

圖3 二維模型不同方法的插值結(jié)果和誤差(a)和(b) 傅里葉POCS插值結(jié)果和誤差； (c)和(d) 低階流式預(yù)測(cè)濾波插值結(jié)果和誤差；(e)和(f) 本文提出HOSPF方法插值結(jié)果和誤差.Fig.3 Interpolation results and errors of different methods for 2D synthetic model(a) and (b) Interpolated result and error by using Fourier POCS method； (c) and (d) Interpolation result and error by using low-order streaming prediction filter； (e) and (f) Interpolated result and error by using the proposed HOSPF.

2.2 三維理論模型測(cè)試

為了進(jìn)一步測(cè)試提出方法在三維數(shù)據(jù)上處理的效果，選擇有多個(gè)彎曲反射界面和一個(gè)斷層的三維Qdome模型，如圖4a所示.圖4b展示了在隨機(jī)去除70%地震道的同時(shí)移除空間范圍35道×35道的方形區(qū)域，模擬連續(xù)缺失模型數(shù)據(jù)，相應(yīng)的數(shù)據(jù)信噪比為1.40 dB.選用三維傅里葉POCS和低階流式預(yù)測(cè)濾波器與本文提出的三維HOSPF插值方法進(jìn)行對(duì)比.由于該模型數(shù)據(jù)中同相軸具有較小的斜率，因此未使用分窗處理，當(dāng)?shù)螖?shù)為84次，算法達(dá)到收斂.與二維模型處理結(jié)果類(lèi)似，三維傅里葉POCS方法能夠較好地重建非連續(xù)缺失的隨機(jī)地震道，在近似線性同相軸的位置處具有較高的精度(圖5a)，但是對(duì)于連續(xù)缺失的位置和斷層及附近的同相軸則產(chǎn)生較大的插值誤差，如圖5b所示.低階流式預(yù)測(cè)濾波器同樣無(wú)法對(duì)連續(xù)缺失的地震道進(jìn)行有效地插值，并在隨機(jī)缺失的部分也產(chǎn)生了較大的誤差，如圖5c、d所示.選用本文提出的三維HOSPF進(jìn)行數(shù)據(jù)的插值重建，濾波器參數(shù)包括：εt為0.00005，εx為0.0008，εy為0.00001，濾波器大小為7(時(shí)間方向)×5(空間X方向)×5(空間Y方向)，正則化參數(shù)的數(shù)值較小的原因在于該數(shù)據(jù)振幅值較小.從圖5e所示的插值結(jié)果可以看到，本文提出的方法在三維情況下仍然能對(duì)大面積連續(xù)缺失地震道進(jìn)行有效地恢復(fù)，除了斷層位置和連續(xù)缺失的地震道相比較隨機(jī)地震道位置的插值精度有一定程度的降低外，總體具有更合理的插值效果.傅里葉POCS、低階算法和本文方法的處理結(jié)果信噪比分別提升至13.77 dB、12.61 dB和16.12 dB，三維傅里葉POCS方法在該模型上收斂的比較快，在使用相同單核CPU的計(jì)算時(shí)間為48.01 s，而低階和高階流式預(yù)測(cè)濾波器方法的計(jì)算時(shí)間分別為21.68 s和23.48 s.

圖4 三維理論模型測(cè)試(a) 三維理論數(shù)據(jù)； (b) 隨機(jī)缺失70%地震道和空間連續(xù)缺失數(shù)據(jù).Fig.4 3D synthetic modeling test(a) 3D synthetic data； (b) The data with 70% randomly removed traces and a large continuous gap.

圖5 三維模型不同方法的插值結(jié)果和誤差(a)和(b) 傅里葉POCS插值結(jié)果和誤差； (c)和(d) 低階流式預(yù)測(cè)濾波插值結(jié)果和誤差； (e)和(f) 本文提出HOSPF方法插值結(jié)果和誤差.Fig.5 Interpolation results and errors of different methods for 3D synthetic model(a) and (b) Interpolated result and error by using Fourier POCS method； (c) and (d) Interpolation result and error by using low-order streaming prediction filter； (e) and (f) Interpolated result and error by using the proposed HOSPF.

3 實(shí)際數(shù)據(jù)測(cè)試

為了驗(yàn)證本文方法在實(shí)際數(shù)據(jù)處理中的適應(yīng)性，選用墨西哥灣海上三維地震數(shù)據(jù)(Claerbout，2008)進(jìn)行測(cè)試.原始數(shù)據(jù)如圖6a所示，拖纜采集的實(shí)施方式造成數(shù)據(jù)近炮檢距缺失.首先對(duì)數(shù)據(jù)進(jìn)行了40%隨機(jī)地震道的去除，并且在共中心點(diǎn)方向模擬30個(gè)共中心點(diǎn)道集的連續(xù)缺失(圖6b).同樣選用三維傅里葉POCS、低階流式預(yù)測(cè)濾波器與本文提出的三維HOSPF插值方法作為對(duì)比，圖7展示了三種方法的處理效果.三維傅里葉POCS在74次迭代時(shí)達(dá)到收斂標(biāo)準(zhǔn)，能夠有效重建隨機(jī)缺失和小范圍連續(xù)缺失的近炮檢距缺失道，但是對(duì)于大范圍連續(xù)缺失的地震道則無(wú)法有效進(jìn)行插值，說(shuō)明該方法的處理能力受限于連續(xù)缺失的地震道數(shù)量(圖7a).低階流式預(yù)測(cè)濾波器在處理連續(xù)缺失和近炮檢距數(shù)據(jù)效果同樣不理想，無(wú)法有效恢復(fù)同相軸的能量(圖7b).接下來(lái)，選取三維HOSPF的濾波器參數(shù)，根據(jù)數(shù)據(jù)的振幅和時(shí)空方向的數(shù)據(jù)平滑性，設(shè)置εt和εy均為0.5，εx為0.01，濾波器大小為15(時(shí)間方向)×40(共中心點(diǎn)方向)×5(炮檢距方向)，本文提出的三維HOSPF相較于傅里葉POCS和傳統(tǒng)流式預(yù)測(cè)濾波器插值方法有更好的插值精度，不僅在大范圍連續(xù)缺失的位置明顯改善插值效果，對(duì)于近炮檢距和隨機(jī)缺失道的重建效果也更加合理，插值后數(shù)據(jù)與相鄰的已知數(shù)據(jù)之間能量更加接近(圖7c).圖8和圖9分別為實(shí)際數(shù)據(jù)炮檢距為0.1675 km位置的共炮檢距剖面插值結(jié)果和誤差的局部放大圖，二者使用了相同的顯示參數(shù)，可以看出在連續(xù)缺失和隨機(jī)缺失位置上，本文提出的方法相較于其他兩種方法插值效果更好，同相軸更加連續(xù)，同時(shí)計(jì)算速度更快，三維傅里葉POCS方法、低階流式預(yù)測(cè)濾波器和本文提出HOSPF方法在使用相同單核CPU的計(jì)算時(shí)間分別為：98.75 s、48.46 s和53.70 s.

圖6 實(shí)際模型測(cè)試(a) 有近炮檢距缺失的三維實(shí)際數(shù)據(jù)； (b) 包含40%隨機(jī)缺失和30個(gè)CMP連續(xù)缺失的數(shù)據(jù).Fig.6 The field modeling test(a) 3D field data with near-offset missing data； (b) Data with 40% randomly removed and 30 continuously missing CMPs.

圖8 實(shí)際數(shù)據(jù)插值結(jié)果局部放大比較(a) 三維傅里葉POCS方法； (b) 低階流式預(yù)測(cè)濾波方法； (c) 本文提出的三維HOSPF方法.Fig.8 Close-up comparison of field data interpolation results(a) 3D Fourier POCS； (b) Low-order streaming prediction filter; (c) The proposed 3D HOSPF.

圖9 實(shí)際數(shù)據(jù)誤差局部放大比較(a) 三維傅里葉POCS方法； (b) 低階流式預(yù)測(cè)濾波方法； (c) 本文提出的三維HOSPF方法.Fig.9 Close-up comparison of field data errors(a) 3D Fourier POCS； (b) Low-order streaming prediction filter; (c) The proposed 3D HOSPF.

4 結(jié)論與展望

本文提出了一種非迭代的時(shí)間域高階流式預(yù)測(cè)濾波(HOSPF)插值方法.通過(guò)高階約束條件進(jìn)一步強(qiáng)化局部平滑性在缺失數(shù)據(jù)插值中的適用性，改善傳統(tǒng)SPF方法在處理連續(xù)缺失數(shù)據(jù)插值問(wèn)題中濾波器難以連續(xù)更新的問(wèn)題，本文提出的方法能夠有效處理大范圍連續(xù)缺失數(shù)據(jù)的插值重建問(wèn)題，提供更合理的插值精度.由于本方法不需要進(jìn)行迭代，所以有高效的計(jì)算速度和較低的計(jì)算成本.同時(shí)，針對(duì)HOSPF的處理特點(diǎn)，選用空間非因果濾波器并設(shè)計(jì)了新的時(shí)間-空間域蛇形插值路徑，進(jìn)一步提高插值精度并改善邊界效應(yīng).對(duì)理論模型和實(shí)際數(shù)據(jù)進(jìn)行處理和分析，并對(duì)比工業(yè)主流傅里葉POCS方法，本文提出的方法在提供更加合理的連續(xù)缺失數(shù)據(jù)插值精度的前提下，占用更少的計(jì)算資源，尤其適用于高維地震數(shù)據(jù)的插值重建問(wèn)題.

致謝感謝同濟(jì)大學(xué)王本鋒研究員的建議和幫助，感謝兩位匿名審稿專家提出的寶貴意見(jiàn).