陳小春 陳 輝 喻 勤 龍 林

（①成都理工大學，四川成都610066；②成都理工大學數(shù)字地質(zhì)四川省重點實驗室，四川成都610066；③中國石化西南油氣分公司勘探開發(fā)研究院，四川成都610041）

反假頻POCS數(shù)據(jù)規(guī)則化及其在偏移成像中的應(yīng)用

陳小春*①陳 輝①②喻 勤③龍 林①

（①成都理工大學，四川成都610066；②成都理工大學數(shù)字地質(zhì)四川省重點實驗室，四川成都610066；③中國石化西南油氣分公司勘探開發(fā)研究院，四川成都610041）

受地形條件和障礙物等因素的影響，陸上地震數(shù)據(jù)采集容易出現(xiàn)地震道丟失現(xiàn)象，對地震資料的后續(xù)處理和解釋造成不良影響。文中提出一種反假頻凸集投影（POCS）插值方法，借助f-k域低頻成分預測高頻無假頻成分，剔除混疊假頻，達到恢復缺失地震道的目的。Marmousi模型數(shù)據(jù)試驗說明凸集投影反假頻插值能提高偏移成像品質(zhì)。

反假頻 凸集投影（POCS） 偏移 規(guī)則化

1 引言

在野外地震數(shù)據(jù)采集過程中，受復雜環(huán)境和施工條件等影響可能導致地震數(shù)據(jù)道的缺失［1］。不規(guī)則數(shù)據(jù)會對疊加、多次波去除產(chǎn)生影響，產(chǎn)生采集腳印和畫弧等問題，降低偏移成像質(zhì)量。因此，在波動方程偏移等高精度處理之前，需要對不規(guī)則地震數(shù)據(jù)進行規(guī)則化處理，完成缺失數(shù)據(jù)重構(gòu)［2］。

針對數(shù)據(jù)規(guī)則化問題，業(yè)界已提出了多種解決方法。Spitz［3］依據(jù)線性同相軸在f-x域的可預測性，提出f-x域預測濾波反假頻插值方法，該方法能較好地處理假頻問題，但僅限于規(guī)則采樣的地震數(shù)據(jù)；國九英等［4］在Spitz方法的基礎(chǔ)上對f-x域地震道插值方法做了改進，并提出f-x域等道距道內(nèi)插方法，利用f-x域預測算子包含了所有反射波同相軸的傾角信息這一特性實現(xiàn)道內(nèi)插，具有速度快、精度高、內(nèi)插出的地震道波形自然、不受空間假頻影響等優(yōu)點，但前提是要求道內(nèi)插之前道間距必須相等；為了使f-x域道內(nèi)插技術(shù)更趨于合理與完善，周興元［5］討論了在f-x域具有固定道時間差信號的廣義預測算子，并給出了道內(nèi)插實現(xiàn)方法；Gulunay等［6］提出抗假頻f-k域地震道插值方法（UFKI），用f-k域復數(shù)除法替代線性方程組求解，實現(xiàn)道加密插值，提高了計算效率，但只適合奇數(shù)道插值偶數(shù)道，或者偶數(shù)道插值奇數(shù)道，即兩道之間只能插值出1道；國九英等［7］提出了f-k域等道距道內(nèi)插方法，該方法具有抗空間假頻等優(yōu)點；Gulunay等［8，9］提出了廣義的f-k域地震道插值方法（GFKI），可以在兩道之間插值出多道；李學聰?shù)龋?0］提出了f-k偏移和反偏移插值法，該方法計算速度快、沒有經(jīng)過近似處理、精度高，但插值道的振幅相比原始道弱；Zwartjes等［11］提出的稀疏傅里葉重建方法具有較高的插值精度，但計算成本較高；Abma等［12］將反假頻凸集投影（POCS）方法從圖像處理領(lǐng)域引入地震數(shù)據(jù)規(guī)則化，該方法易于實現(xiàn)且重建效果良好，但是不能抗假頻，對隔道缺失的不完整地震數(shù)據(jù)插值無效，Gao等［13］通過引入傾角掃描策略解決了這一難題。

本文首先闡述了POCS數(shù)據(jù)規(guī)則化方法和反假頻技術(shù)的原理及實現(xiàn)，在POCS算法基礎(chǔ)上提出一種新的預測反假頻數(shù)據(jù)規(guī)則化方法。該方法應(yīng)用f-k域無假頻低頻數(shù)據(jù)建立無假頻蒙板，再經(jīng)過線性放大逐步預測高頻數(shù)據(jù)，剔除混疊假頻，達到恢復缺失地震道的目的。方法繼承了POCS算法的高效性，同時具有很強的反假頻能力。數(shù)值實驗驗證了方法的數(shù)據(jù)重建能力和在偏移成像中的作用。

2 方法原理

2.1 POCS數(shù)據(jù)規(guī)則化

受地形條件、障礙物、廢道和壞道等因素的影響，實際采集地震數(shù)據(jù)多為非規(guī)則采樣，而非規(guī)則采樣數(shù)據(jù)通常認為是由規(guī)則數(shù)據(jù)采樣不足造成的。實際采集數(shù)據(jù)可以認為是理論觀測系統(tǒng)的數(shù)據(jù)通過稀疏采樣得到的，即可以表示為

式中：dobs為含有缺失道的采集數(shù)據(jù)；S為稀疏采樣算子；d為理論數(shù)據(jù)。地震數(shù)據(jù)重建問題就是由dobs恢復d。對于式（1），通過最小二乘反演可以重建缺失道，以實現(xiàn)觀測系統(tǒng)數(shù)據(jù)的規(guī)則化。根據(jù)L2范數(shù)最小準則建立目標函數(shù)

由于S為稀疏采樣算子，其稀疏性會導致方程欠定、造成多解性。因此，在目標函數(shù)后增加阻尼項、建立罰函數(shù)

式中：C為稀疏變換算子；μ為阻尼項的權(quán)重。最小二乘法計算復雜度高、效率相對低下，對于阻尼最小二乘法來說，阻尼項在約束多解性的同時也降低了算法的收斂性能，進一步影響了最小二乘法的效率。POCS迭代方法具有更強的穩(wěn)健性和收斂性，計算效率也遠高于最小二乘法。在迭代步長滿足收斂條件，即d i接近d的時候，式（1）可改寫為

式中：i為迭代次數(shù)；d i為第i次迭代后得到的重建結(jié)果；E是單位矩陣；F及F-1分別是正、反傅氏變換算子；τi是第i次迭代時的閾值算子。

地震數(shù)據(jù)在f-k域（傅氏變換域）內(nèi)的主要能量集中在其主頻附近。缺失的地震道在觀測系統(tǒng)網(wǎng)格中會形成空道，這些空道在f-k譜中表現(xiàn)為一些幅值相對較小的噪聲，缺失地震道重構(gòu)在f-k譜中則等效于去噪和頻譜能量恢復。將含空道的原始地震數(shù)據(jù)變換到f-k域，若原始數(shù)據(jù)中的空道數(shù)量較少，則f-k譜中有效信號的能量會強于噪聲能量，在這種情況下可以通過設(shè)置一個適當?shù)拈撝祵⒌湍芰康念l譜視為噪聲截去，僅保留強能量部分。然后對所保留的頻譜進行傅里葉反變換，利用頻譜泄露的特性可以在空道上部分恢復出缺失信號。由于閾值的辦法并不能準確區(qū)分有效信號和噪聲，僅通過一次重建無法得到準確的結(jié)果。因此可降低閾值重復上述操作，經(jīng)過數(shù)次迭代后f-k域中的噪聲能量閾值降到遠低于有效信號的水平，這時可認為數(shù)據(jù)被理想重建。

圖1a和圖1b是非規(guī)則數(shù)據(jù)和規(guī)則化后的數(shù)據(jù)，圖1c和圖1d是其對應(yīng)的頻譜。從圖1可以看出，非規(guī)則數(shù)據(jù)在f-k域內(nèi)的噪聲被剔除后，t-x域的空道數(shù)據(jù)得到了很好重建，且與原始數(shù)據(jù)保持同相軸連續(xù)。

在實際數(shù)據(jù)處理過程中，閾值的選擇會影響POCS數(shù)據(jù)規(guī)則化的效果。如果閾值選取過大，會導致有效信號被切除，無法正確恢復有效信號的細節(jié)。如果閾值選取過小，則大量噪聲被當作有效信號予以保留，導致重建信號能量弱。為避免這兩種情況，建議在迭代初期選擇較大的閾值，以凸顯信號同相軸的輪廓，而在后續(xù)的迭代過程中按照一定的步長衰減閾值，逐步恢復低能量部分的細節(jié)，直到閾值遠低于所估計的信噪比。這樣可在不傷及有效信號的前提下盡可能完整地恢復信號。本文采用線性衰減閾值，在模型中給定最大閾值τmax、最小閾值τmin以及迭代次數(shù)N，則第i次迭代時的閾值為

2.2 預測反假頻技術(shù)

傳統(tǒng)POCS數(shù)據(jù)規(guī)則化方法得以實現(xiàn)的前提是假頻能量明顯小于有效頻譜的能量。當數(shù)據(jù)含有規(guī)則缺失道時，在f-k譜中會產(chǎn)生與有效頻率能量相當?shù)募兕l，而且假頻通常與有效頻譜混疊，此時就無法通過閾值的方法將假頻和有效頻譜分離，因此本文提出一種預測反假頻技術(shù)。

首先，通過頻譜分析在f-k譜中確定無假頻區(qū)域的范圍。將含規(guī)則缺失道的數(shù)據(jù)變換到f-k域，如圖2a所示。從圖中可以看到，在有效頻率兩側(cè)形成嚴重假頻，且在高頻成分假頻與有效頻率混疊嚴重，但在圖2a中灰色的Ⅰ區(qū)，即頻率范圍為［0，f1］，波數(shù)范圍為［-k1，k1］，不存在假頻混疊。

然后對Ⅰ區(qū)的頻譜按照能量二值化，即將譜能量較大的部分置為1，很小的部分置為0，就可得到無假頻蒙板。對該蒙板線性拉伸放大，擴展到Ⅱ區(qū)（圖2b），則蒙板放大前、后非重疊部分即為所預測的高頻有效信號的區(qū)域。

圖1 POCS插值方法對含空道非規(guī)則數(shù)據(jù)的插值重建效果

圖2 預測反假頻技術(shù)實現(xiàn)

由于蒙板本身進行了二值化，再將蒙板與頻譜進行點乘，假頻混疊部分乘零被剔除，只剩下有信號的有效成分。從低頻向高頻部分逐步預測，直到剔除所有混疊的假頻。

2.3 基于POCS的預測反假頻數(shù)據(jù)規(guī)則化

預測反假頻技術(shù)操作簡單而且能很容易地剔除混疊假頻。經(jīng)過拉伸放大后的反假頻蒙板可以用來預測有效頻率的大致范圍，然而采用拉伸放大的方式?jīng)Q定了預測方式是線性的，實際數(shù)據(jù)的同相軸不可能做到準確的線性，這個問題可以通過開窗的辦法來解決。雖然開窗可以使局部區(qū)域的同相軸線性化，但在有噪聲的情況下，這種預測方式獲得蒙板也會因為噪聲而產(chǎn)生誤差。這就難免會在一定程度上造成假頻剔除不干凈，或者傷及有效信號。但是經(jīng)過反假頻處理后，大部分假頻能量已經(jīng)被去除，誤差導致的殘留量已經(jīng)變成隨機噪聲，在反假頻處理之后可結(jié)合傳統(tǒng)的POCS數(shù)據(jù)規(guī)則化方法進一步剔除噪聲從而恢復有效信號。

預測反假頻技術(shù)基于線性插值理論，而實際數(shù)據(jù)的同相軸不可能做到準確的線性，故反假頻蒙板每次線性拉伸的步長會決定預測的精度。拉伸步長越大，預測精度越低；拉伸步長越小，預測精度越高，但是計算量也會相應(yīng)的增加，必須選擇合適的頻率步長從低頻向高頻逐步處理。因此頻率步長就成為權(quán)衡預測精度和計算量的重要參數(shù)。

基于POCS預測的反假頻數(shù)據(jù)規(guī)則化處理步驟（圖3）為：

（1）將窗口數(shù)據(jù)做二維傅氏正變換，由t-x域變換到f-k域；

（2）在f-k域中確定低頻無混疊假頻成分的頻率范圍，記為［0，f1］；

（3）對確定的無假頻成分求模并做線性拉伸，拉伸到頻率范圍為［0，f1+df］；

（4）將拉伸后的能量譜二值化生成反假頻蒙板；

（5）將反假頻蒙板與原始譜相乘剔除混疊假頻后做傳統(tǒng)POCS數(shù)據(jù)規(guī)則化；

（6）判斷是否已經(jīng)處理到最大頻率，若處理到最大頻率計算結(jié)束，否則擴大無假頻區(qū)域進行迭代；

（7）將規(guī)則化后的數(shù)據(jù)反變換到t-x域，將被恢復的道回填至有效信號中。

圖3 預測反假頻數(shù)據(jù)規(guī)則化處理流程

3 數(shù)據(jù)測試

為了檢驗預測反假頻數(shù)據(jù)規(guī)則化方法的效果，利用簡單線性模型數(shù)據(jù)進行測試。該模型是由主頻為60Hz的Ricker子波合成，它包含8條斜率不同的直線同相軸，空間方向有110道，道間距為15m，時間方向共620個采樣點，采樣率為2ms，在該模型中數(shù)據(jù)被嚴重抽稀，只保留了三分之一的數(shù)據(jù)（圖4a）。在實驗中進行傳統(tǒng)POCS數(shù)據(jù)規(guī)則化處理時給定的最大閾值為0.4，而最小的閾值為0.001，反假頻數(shù)據(jù)規(guī)則化的迭代頻率步長為6 Hz。圖4c為其頻譜，從中可以看出因數(shù)據(jù)道缺失而產(chǎn)生的幅值與主頻相當?shù)募兕l，且在60Hz以上出現(xiàn)假頻混疊。圖4d為重建后的頻譜，經(jīng)過預測反假頻數(shù)據(jù)規(guī)則化處理后假頻得到有效的壓制，特別是在混疊部分既壓制了假頻，又很好地保護了有效信號。重建后的頻譜反變換回t-x域，如圖4b所示?？盏罃?shù)據(jù)被重建，從同相軸的連續(xù)性上看線性模型數(shù)據(jù)得到了正確的重建。

為驗證該方法在噪聲干擾下的效果，對上述數(shù)據(jù)加上信噪比為30%的隨機噪聲。經(jīng)過同樣的處理（圖5）可以看出，在隨機噪聲和假頻的同時干擾下，本文提出的方法依然可以有效地重建缺失道信號。

為了更加準確地分析重建效果，提取原始數(shù)據(jù)、不含噪聲的重建數(shù)據(jù)和含有噪聲的重建數(shù)據(jù)的第76道、0.40～0.50s的波形進行疊加顯示（圖6），可以看出三組數(shù)據(jù)的走時和振幅都基本吻合，其中原始數(shù)據(jù)和不含噪聲的重建數(shù)據(jù)吻合度更高，當受到噪聲影響的時候重建結(jié)果的波形和理論值存在一定的差異。

在實際數(shù)據(jù)處理過程中可能會遇到一些極端情況。比如低頻部分依然存在假頻混疊，就意味著數(shù)據(jù)非常稀疏，或者隨機噪聲過強導致無法通過閾值辨認。數(shù)據(jù)質(zhì)量太差時，算法選擇往往不起作用，只能通過重新采集解決。

在數(shù)據(jù)缺失非常嚴重的情況下，難以完全精確重建缺失的地震道，僅根據(jù)同相軸的連續(xù)性無法判斷重建算法對最終處理結(jié)果的影響。本文通過Marmousi模型數(shù)據(jù)規(guī)則化前、后道集及偏移成像結(jié)果驗證數(shù)據(jù)規(guī)則化對成像效果的改善能力。為了更接近實際情況，本文沒有使用理論速度模型，而是用Seis＠Home軟件的偏移和速度分析模塊構(gòu)建時間域的偏移速度模型，如圖7所示。

圖4 簡單線性模型理論數(shù)據(jù)測試結(jié)果

圖5 簡單線性模型含隨機噪聲數(shù)據(jù)測試結(jié)果

圖6 三組數(shù)據(jù)第76道0.40～0.50s的波形對比

Marmousi模型數(shù)據(jù)的時域樣點數(shù)為726個，時間采樣率為4ms。圖8a為理論的共炮檢距道集數(shù)據(jù)，圖8b為稀疏采樣二分之一后的數(shù)據(jù)，圖8c為數(shù)據(jù)規(guī)則化重建后的數(shù)據(jù)。由于Marmousi模型數(shù)據(jù)本身就比較稀疏，在大傾角構(gòu)造區(qū)域肉眼可見明顯的鋸齒狀，減少二分之一的道數(shù)據(jù)后稀疏程度更加嚴重，所以重建效果也不可能完全與原始數(shù)據(jù)相同。但是從Kirchhoff疊前時間偏移結(jié)果（圖9）看，抽稀的數(shù)據(jù)（所有共炮檢距數(shù)據(jù)都進行相同的稀疏采樣）做偏移處理（圖9b）后，淺層出現(xiàn)嚴重的偏移畫弧現(xiàn)象，影響成像結(jié)果解釋。然而規(guī)則化后，偏移成像（圖9c）整體得到改善，特別是在中淺層效果尤為明顯，偏移畫弧問題得到解決。雖然在一些高陡構(gòu)造區(qū)域成像能量較弱，但是有效消除了偏移噪聲，突出了有效信息。

圖7 時間域偏移速度建模

圖8 Marmousi模型數(shù)據(jù)規(guī)則化效果

圖9 數(shù)據(jù)規(guī)則化前后偏移成像結(jié)果對比

4 結(jié)論

（1）傳統(tǒng)的POCS數(shù)據(jù)規(guī)則化方法只能重建無混疊假頻的數(shù)據(jù)，而對含有混疊假頻的數(shù)據(jù)不起作用。

（2）基于POCS算法的預測反假頻技術(shù)能對含有混疊假頻的數(shù)據(jù)進行良好的重建，而且只需較小的計算成本就能得到高質(zhì)量的結(jié)果，恢復后的數(shù)據(jù)能量與其相鄰道的有效數(shù)據(jù)相當，并且保持了同相軸的連續(xù)性。

（3）數(shù)據(jù)經(jīng)過數(shù)據(jù)規(guī)則化處理后，偏移結(jié)果中畫弧問題得到明顯改善，更利于資料的解釋。

［1］ 陸基孟.地震勘探原理（第三版）.山東東營：中國石油大學出版社，2011.

［2］ 黃德濟，賀振華，包吉山.地震勘探資料數(shù)字處理.北京：地質(zhì)出版社，1990.

［3］ Spitz S.Seismic trace interpolation in theF-Xdomain.Geophysics，1991，56（6）：785-794.

［4］ 國九英，周興元，俞壽朋.F-X域等道距道內(nèi)插.石油地球物理勘探，1996，31（1）：28-34.Guo Jiuying，Zhou Xingyuan，Yu Shoupeng.Iso-interval trace interpolation inF-Xdomain.OGP，1996，31（1）：28-34.

［5］ 周興元.F-X域地震道內(nèi)插理論及實現(xiàn).石油地球物理勘探，1997，32（2）：154-162.Zhou Xingyuan.Seismic trace interpolation inF-Xdomain.OGP，1997，32（2）：154-162.

［6］ Gulunay N，Chambers E.UnaliasedF-Kdomain trace interpolation.SEG Technical Program Expanded Abstracts，1996，15：1461-1464.

［7］ 國九英，周興元.F-K域等道距道內(nèi)插.石油地球物理勘探，1996，31（2）：211-218.Guo Jiuying，Zhou Xingyuan.Iso-interval trace interpolation inF-Kdomain.OGP，1996，31（2）：211-218.

［8］ Gulunay N，Chambers E.Generalizedf-kdomain trace interpolation.SEG Technical Program Expanded Abstracts，1997，16：1100-1103.

［9］ Gulunay N.Seismic trace interpolation in the Fourier transform domain.Geophysics，2003，68（1）：355-369.

［10］ 李學聰，劉伊克，常旭等.基于F-K偏移和反偏移的地震道插值方法研究.地球物理學進展，2009，24（6）：2020-2029.Li Xuecong，Liu Yike，Chang Xu et al.Interpolation for seismic traces based onF-Kmigration and demigration.Progress in Geophysics，2009，24（6）：2020-2029.

［11］ Zwartjes P，Sacchi M.Fourier reconstruction of nonuniformly sampled，aliased seismic data.Geophysics，2004，72（1）：21-32.

［12］ Abma R，Kabir N.3D interpolation of irregular data with a POCS algorithm.SEG Technical Program Expanded Abstracts，2005，24：2150-2153.

［13］ Gao J，Stanton A，Naghizadeh M et al.Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction.Geophysical Prospecting，2013，61（S1）：138-151.

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.003

陳小春，陳輝，喻勤，龍林.反假頻POCS數(shù)據(jù)規(guī)則化及其在偏移成像中的應(yīng)用.石油地球物理勘探，2017，52（1）：13-19.

1000-7210（2017）01-0013-07

*四川省成都市成都理工大學地球物理學院地球物理系，610066。Email：974609310＠qq.com

本文于2015年12月4日收到，最終修改稿于2016年12月6日收到。

本項研究受國家自然科學基金項目（41304111）、四川省科技廳重點項目（2016JY0200）、四川省教育廳自然科學重點項目（14ZA0061，16ZB0101）、“油氣地球物理勘探”四川省青年科技創(chuàng)新研究團隊專項計劃（2016TD0023）、成都理工大學優(yōu)秀創(chuàng)新團隊培育計劃（KYTD201410）及國家科技重大專項（2016ZX05002-004-05）等聯(lián)合資助。

（本文編輯：宜明理）

陳小春 1994年生；現(xiàn)為成都理工大學地球物理學院地球物理專業(yè)本科生，研究方向為地震信號高分辨率處理方法和重建方法及模塊開發(fā)。