基于改進匹配追蹤傅里葉插值的地震數(shù)據(jù)規(guī)則化重構(gòu)

張春楓 杜啟振，2 張富源 符力耘，2 魏國華

（1. 中國石油大學（華東）山東省深層油氣重點實驗室，山東 青島 266580；2. 青島海洋科學與技術(shù)試點國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東 青島 266237；3. 中國石化勝利油田分公司物探研究院，山東 東營 257022）

0 引 言

隨著中國油氣資源勘探開發(fā)逐步深化，對野外采集的地震資料質(zhì)量提出了更高的要求。理論上常常假定，地震數(shù)據(jù)沿著空間方向上按等距離均勻采樣。在野外地震數(shù)據(jù)采集過程中，地震資料在時間方向上可以做到均勻、密集采樣，受制于實際施工環(huán)境等因素的影響，在空間方向上的采樣往往是不規(guī)則甚至是稀疏的，這會對后續(xù)的地震數(shù)據(jù)處理流程（地震數(shù)據(jù)去噪、CMP 道集疊加、偏移成像等）產(chǎn)生不良影響［1］。為使得基于規(guī)則采樣假設(shè)條件的常規(guī)地震數(shù)據(jù)處理流程得以順利進行，不規(guī)則地震數(shù)據(jù)的插值與規(guī)則化重構(gòu)具有十分重要的意義。

地震數(shù)據(jù)規(guī)則化的方法主要分為以下6 類：第1 類是采用線性預(yù)測濾波的方法進行重構(gòu)［2-5］；第2類是以波動方程為基礎(chǔ)的重構(gòu)方法［6］；第3 類是矩陣降秩的重構(gòu)方法［7-8］；第4 類是相干傾角的重構(gòu)方法［9］；第5 類是基于數(shù)學變換域的重構(gòu)方法；第6 類是基于人工智能的重構(gòu)方法［10-11］。目前，基于傅里葉變換域的重構(gòu)方法受到頗多關(guān)注，其優(yōu)點在于不受空間采樣影響，可以很好地拓展到高維領(lǐng)域。A.J.W.Duijndam 等［12］將傅里葉變換應(yīng)用于二維不規(guī)則地震數(shù)據(jù)的恢復(fù)與重構(gòu)；K.Hindriks 等［13］又將其拓展到三維領(lǐng)域；B.Liu 等［14］采用最小加權(quán)范數(shù)對不規(guī)則地震數(shù)據(jù)進行插值重構(gòu)；S.Xu 等［15-16］結(jié)合了非均勻離散傅里葉變換，提出了基于抗泄露傅里葉變換（Anti-leakage Fourier Transform，ALFT）的不規(guī)則地震數(shù)據(jù)重構(gòu)，并將其應(yīng)用于三維地震數(shù)據(jù)重構(gòu)；在S.Xu 等研究的基礎(chǔ)上，ALFT 算法中傅里葉系數(shù)計算部分也被改進與優(yōu)化［17-19］，包括引入非均勻快速傅里葉變換、GPU 加速算法等。傅里葉變換方法的重構(gòu)效率受傅里葉系數(shù)計算時長的影響，加快傅里葉系數(shù)的計算速度是提升不規(guī)則地震數(shù)據(jù)重構(gòu)效率的關(guān)鍵之一。

為解決傅里葉系數(shù)計算耗時高這一問題，本文采用與ALFT 原理相近的匹配追蹤傅里葉插值（Matching Pursuit Fourier Interpolation，MPFI） 算法，在MPFI 算法的基礎(chǔ)上對其迭代結(jié)構(gòu)進行改進，在不規(guī)則地震數(shù)據(jù)重構(gòu)前計算其傅里葉系數(shù)矩陣。和傳統(tǒng)MPFI 算法相比，改進后的算法優(yōu)化了循環(huán)中大量復(fù)雜的計算，減少迭代所需時間，從而提高不規(guī)則地震數(shù)據(jù)的重構(gòu)效率。

1 地震數(shù)據(jù)重構(gòu)原理

采用傅里葉變換將不規(guī)則數(shù)據(jù)規(guī)則化符合信號稀疏分解理論。信號稀疏分解是指具有稀疏性的信號可以在重構(gòu)數(shù)據(jù)最少、重構(gòu)誤差最小的情況下，由有限個最優(yōu)重構(gòu)數(shù)據(jù)的線性加權(quán)組合表示，即使用較少的數(shù)據(jù)重構(gòu)出與原信號最接近或者較為接近的解。對于規(guī)則缺失的地震數(shù)據(jù)，傅里葉變換重構(gòu)十分有效；但是對于不規(guī)則采樣地震數(shù)據(jù)，其傅里葉基函數(shù)不再互相正交，導(dǎo)致各個頻率之間發(fā)生頻譜泄露，即某個傅里葉系數(shù)的能量泄漏到其他的傅里葉系數(shù)上，能量最強的頻率造成的泄漏最大。此時，常規(guī)的傅里葉重構(gòu)方法不再適用，亟待發(fā)展新的方法。S.G.Mallat 等［20］提出了一種基于稀疏計算的匹配追蹤算法（Matching Pursuit，MP）來解決這一問題。任何信號的波形都可以通過一系列正弦波累加得到，可用傅里葉級數(shù)表示為

式中：f(x)——信號；an(n=0，±1，±2，…)——權(quán)重系數(shù)；i——虛數(shù)單位；ω——角頻率；x——空間坐標，m。

地震數(shù)據(jù)規(guī)則化與匹配追蹤算法所解決的問題類似，采用匹配追蹤與傅里葉變換相結(jié)合的方式對地震數(shù)據(jù)作規(guī)則化與重構(gòu)，稱其為匹配追蹤傅里葉插值。MPFI 算法以傅里葉變換為基礎(chǔ)，其基本思想是先計算不規(guī)則數(shù)據(jù)的傅里葉系數(shù)，將其保存到系數(shù)譜中，循環(huán)多次后再把重構(gòu)的傅里葉系數(shù)譜經(jīng)過反傅里葉變換，輸出到期望的空間位置上。

實現(xiàn)MPFI 算法的前提條件是不規(guī)則數(shù)據(jù)具有稀疏性。實際上野外采集得到的t-x域地震數(shù)據(jù)并不具有稀疏性。野外采集地震數(shù)據(jù)在時間方向上是規(guī)則采樣，采樣間隔相同，采用FFT 或者DFT就能將其從時間域變換到頻率域；受各種因素的影響，空間方向上的檢波點排布通常是不規(guī)則的，需要在地震數(shù)據(jù)的空間方向上作非均勻離散傅里葉變換，這樣才能在f-k域中通過迭代計算、更新傅里葉系數(shù)譜。變換后，地震數(shù)據(jù)中大多數(shù)波數(shù)接近于0，滿足不規(guī)則數(shù)據(jù)的稀疏性條件。而且變換后傅里葉系數(shù)譜上的能量更加集中，有利于后續(xù)傅里葉系數(shù)的計算與選取。將最終選取的一系列傅里葉系數(shù)組成新的系數(shù)譜，再將其經(jīng)過反傅里葉變換后重構(gòu)成規(guī)則采樣數(shù)據(jù)。

常規(guī)的MPFI 算法針對不規(guī)則數(shù)據(jù)進行多次迭代，能夠?qū)⑵湫孤┑哪芰炕貧w，從而恢復(fù)原始數(shù)據(jù)。為了壓制能量泄露，計算地震數(shù)據(jù)中其他頻率成分的系數(shù)之前需要對輸入數(shù)據(jù)進行更新。絕對值最大的傅里葉系數(shù)對頻譜泄露的影響最強。故在更新不規(guī)則數(shù)據(jù)時應(yīng)從絕對值最大的頻率成分開始，在原數(shù)據(jù)中減去該頻率在f-x域的貢獻，以此類推，逐漸遞歸到絕對值較小的頻率。

MPFI 與ALFT 算法在計算和尋找傅里葉系數(shù)方面基本相同，T.Nguyen 等［21］開發(fā)了一種快速MPFI 算法，其中NDFT 在每個時間頻率下僅計算一次，從而大大降低了MPFI 算法的計算時間。

對于變換后的f-x域地震數(shù)據(jù)f(xl)，可以用離散和的形式表示為

其逆變換為

式中：f s(xl)——第s次迭代中的f-x域地震數(shù)據(jù)；s——迭代次數(shù)；xl(l= 1，2，…，N)——檢波器坐標，m；N——檢波器數(shù)量；q——空間波數(shù)索引；Nf——空間頻率數(shù)量；(kq)——第s次迭代中的傅里葉系數(shù)；kq——空間波數(shù)。

MPFI 算法中需要不斷地對地震數(shù)據(jù)進行更新，更新量由所選取傅里葉系數(shù)對應(yīng)的貢獻值決定，系數(shù)越大，泄露在頻譜上的能量就越多。每次迭代中最強傅里葉系數(shù)對應(yīng)空間波數(shù)ks可由公式計算得到

式中ks——第s次迭代中所選取的空間波數(shù)。

初始化迭代次數(shù)s= 1，s+ 1 次迭代后，不規(guī)則地震數(shù)據(jù)的剩余量表示為

式中：f s+1(xl)——第s+ 1 次迭代中的f-x域地震數(shù)據(jù)；?(ks)——第s次迭代中的最大傅里葉系數(shù)。

常規(guī)MPFI 算法的重構(gòu)步驟［22］：

（1）在地震數(shù)據(jù)的時間方向上作FFT 變換，得到f-x域地震數(shù)據(jù)；

（2）在空間方向上，對輸入的數(shù)據(jù)作非均勻離散傅里葉變換，計算其傅里葉系數(shù)譜；

（3）選擇能量最強的傅里葉系數(shù)，保存在新系數(shù)譜中；

（4）將選中的系數(shù)按照采樣點的實際空間位置作非均勻反傅里葉變換；

（5）原始不規(guī)則f-x域數(shù)據(jù)減去該次反傅里葉變換的結(jié)果，更新不規(guī)則數(shù)據(jù)，再次循環(huán)；

（6）不斷重復(fù)循環(huán)步驟（2）到（5），直到滿足迭代誤差或者預(yù)定義的迭代次數(shù)時為止；

（7）將保存的傅里葉系數(shù)譜按照期望空間位置進行反傅里葉變換，得到重構(gòu)后規(guī)則地震數(shù)據(jù)。

2 重構(gòu)算法的改進

MPFI 算法的主要不足是數(shù)據(jù)重構(gòu)計算時間長、效率低，這是因為循環(huán)中涉及到大量非均勻離散正反傅里葉變換的計算，特別是循環(huán)中的第一步，傅里葉系數(shù)的計算量巨大。循環(huán)過程中不規(guī)則數(shù)據(jù)每更新一次，其剩余量中的傅里葉系數(shù)都要重新計算一遍。如果能夠在迭代前提前估算出這些傅里葉系數(shù)，那么在循環(huán)中就可以省略大量復(fù)雜運算，不規(guī)則地震數(shù)據(jù)的重構(gòu)效率將會大大提高。計算傅里葉系數(shù)的方法多種多樣。本文采用直接計算傅里葉系數(shù)的方式，對公式（5）兩邊同時作非均勻離散傅里葉逆變換，并結(jié)合公式（3）得到

整理后，可得到不規(guī)則采樣地震數(shù)據(jù)的剩余量在f-k域中的更新迭代公式

式中(kq)——第s+ 1 次迭代中的傅里葉系數(shù)。

公式（7）右邊指數(shù)項的計算量仍然很大，但是其中的空間波數(shù)kq、偏移距xl只與地震數(shù)據(jù)有關(guān)，而每次選中的最強傅里葉系數(shù)所對應(yīng)波數(shù)ks是空間波數(shù)集合kq中的元素，因此改進后迭代公式中的復(fù)雜指數(shù)項可以提前計算得到，并且只計算一次。與迭代公式（5）相比，迭代公式（7）將不規(guī)則地震數(shù)據(jù)剩余量的計算從f-x域轉(zhuǎn)換到f-k域，地震數(shù)據(jù)只在f-k域中更新迭代。此優(yōu)勢是能夠避免迭代中不規(guī)則地震數(shù)據(jù)在f-x域與f-k域之間頻繁轉(zhuǎn)換，減少了循環(huán)中大量非均勻離散傅里葉變換，從而達到提高計算效率的目的。

改進后使用MPFI 算法重構(gòu)地震數(shù)據(jù)的過程變?yōu)椋?/p>

（1）計算地震數(shù)據(jù)的傅里葉分量矩陣；

（2）在地震數(shù)據(jù)的時間方向和空間方向分別作傅里葉變換；

（3）選擇能量最強的傅里葉系數(shù)，并保存在新傅里葉系數(shù)譜中；

（4）使用改進后的地震數(shù)據(jù)剩余量迭代公式在f-k域中更新不規(guī)則地震數(shù)據(jù)；

（5）重復(fù)進行步驟（3）到（4），直到剩余量滿足誤差或者達到預(yù)定義的迭代次數(shù)為止；

（6）將保存的傅里葉系數(shù)譜按照規(guī)則的空間采樣點位置進行反傅里葉變換，得到規(guī)則地震數(shù)據(jù)。

3 模型測試

3.1 凹陷模型

為了檢驗改進后MPFI 算法的可行性，采用二維凹陷速度模型作時間二階、空間十階聲波正演模擬測試。

二維凹陷速度模型的橫向采樣點為800，縱向采樣點為440，橫向網(wǎng)格寬度為5 m，縱向網(wǎng)格寬度為5 m，上層速度為2 500 m/s，下層速度為2 800 m/s（圖1）。

圖1 凹陷速度模型Fig. 1 Graben velocity model

采用完全匹配層（Perfectly Matched Layer，PML）吸收邊界條件，炮點位于橫向2 000 m、縱向25 m 處，空間方向共160 道，道間距為25 m，時間方向共1 700 個采樣點，采樣間隔時間為1 ms，雷克子波的主頻為25 Hz。正演模擬后得到規(guī)則采樣單炮地震記錄（圖2（a））。

圖2 凹陷模型正演與重構(gòu)結(jié)果Fig. 2 Graben model forward modeling and reconstruction result

一般來說，對于單炮地震記錄，直達波的能量強于地下分界面所反射的有效波的能量。重構(gòu)地震數(shù)據(jù)時，如果直達波與反射波位置較為接近或者疊合在一起，應(yīng)先切除強直達波，否則在選擇傅里葉系數(shù)時往往選取的是直達波的信息而不是有效波的信息，會嚴重影響到后續(xù)的地震數(shù)據(jù)重構(gòu)。對規(guī)則正演地震記錄中的檢波器作隨機缺失（缺失數(shù)量為32），得到欠采樣單炮地震記錄（圖2（b））。相對于規(guī)則采樣數(shù)據(jù)，隨機欠采樣地震數(shù)據(jù)的同相軸出現(xiàn)了中斷、不連續(xù)現(xiàn)象。使用MPFI 算法與改進MPFI 算法重構(gòu)后，隨機缺失的檢波點被重構(gòu)到規(guī)則的期望空間位置上，且重構(gòu)后正演結(jié)果中同相軸變得光滑、連續(xù)，隨機欠采樣地震數(shù)據(jù)被重構(gòu)為規(guī)則采樣地震數(shù)據(jù)（圖2（c）、（d））。

相比于規(guī)則采樣數(shù)據(jù)的頻波譜（圖3（a）），不規(guī)則采樣數(shù)據(jù)的頻波譜上還出現(xiàn)了能量泄露，其能量峰值降低，頻率泄漏到整個頻波譜內(nèi)（圖3（b）），使用MPFI算法與改進MPFI算法重構(gòu)后，頻波譜中泄露的能量回歸到有效頻率中（圖3（c）、（d））。單道地震數(shù)據(jù)對比表明（圖4），使用2種算法重構(gòu)后的單道地震數(shù)據(jù)與原始單道地震數(shù)據(jù)相比，誤差基本相同，2 種方法重構(gòu)后的單道數(shù)據(jù)十分接近，可認為2 種方法重構(gòu)效果一致。上述結(jié)果表明在二維凹陷速度模型中，改進后的重構(gòu)方法具有可行性。

圖5 Marmousi聲波速度模型Fig. 5 Marmousi acoustic velocity model

3.2 復(fù)雜模型

為了檢驗改進后MPFI 算法在復(fù)雜模型中的可行性，采用抽稀Marmousi 速度模型作時間二階、空間十階有限差分聲波正演模擬測試。為了減弱直達波對反射波的干擾，將模型上方水層加厚。如圖5 所示，該速度模型的橫向采樣點數(shù)為530，縱向采樣點數(shù)為227，采用PML 吸收邊界條件，其色標表示介質(zhì)中聲波傳播速度。

正演模擬中設(shè)置的參數(shù)為：地震數(shù)據(jù)的空間方向共有106 道，橫向網(wǎng)格寬度為5 m，縱向網(wǎng)格寬度為5 m，道間距為25 m，時間方向共有1 601 個采樣點，時間采樣間隔為1 ms，炮點坐標位于橫向1 325 m、縱向25 m 處，雷克子波的主頻設(shè)為25 Hz。正演模擬后得到規(guī)則采樣單炮地震記錄（圖6（a））。將正演模擬中的檢波器按照隨機缺失的方式進行數(shù)據(jù)采樣（缺失數(shù)量為24），以模擬實際采集得到的不規(guī)則地震數(shù)據(jù)，可以看出，不規(guī)則采樣地震數(shù)據(jù)中的反射波同相軸出現(xiàn)了隨機中斷，有明顯的不連續(xù)現(xiàn)像（圖6（b）），經(jīng)過重構(gòu)，同相軸缺失信息得到補償，地下波場得到恢復(fù)（圖6（c）、（d））。

圖6 Marmousi模型正演與重構(gòu)結(jié)果Fig. 6 Marmousi model forward modeling and reconstruction results

本文采用的2 種重構(gòu)方法均是依靠計算傅里葉系數(shù)譜重構(gòu)不規(guī)則地震數(shù)據(jù)，改進后的MPFI 算法是在常規(guī)算法的基礎(chǔ)上對傅里葉系數(shù)的計算方式作出優(yōu)化，對迭代公式兩端同時應(yīng)用傅里葉變換，因此2 種方法在重構(gòu)效果上基本相同，但改進后的重構(gòu)方法計算效率更加高效，這是因為使用MPFI 方法重構(gòu)不規(guī)則地震數(shù)據(jù)的效率還與重建過程中傅里葉變換次數(shù)有關(guān)。

對于相同的地震數(shù)據(jù)，其數(shù)據(jù)量越大、頻率范圍越寬，循環(huán)中所涉及的非均勻離散傅里葉變換就越多。這些非均勻離散傅里葉變換在計算機中需要使用歐拉公式計算，因此重構(gòu)時間就越長。采用MPFI 算法重構(gòu)凹陷不規(guī)則采樣與抽稀Marmousi 不規(guī)則采樣數(shù)據(jù)，在計算傅里葉系數(shù)譜時分別需要進行大量空間方向上的非均勻傅里葉變換，在所占內(nèi)存相同的情況下，使用改進后的MPFI 算法在計算傅里葉系數(shù)譜的過程中不需進行空間方向上的非均勻傅里葉變換，計算量相比于優(yōu)化前大大減少，而且迭代公式（7）中涉及的復(fù)雜運算在循環(huán)前就已經(jīng)一次完成，不需要在每次循環(huán)中去重復(fù)計算傅里葉系數(shù)。此外2 種方法在時間方向上的快速傅里葉變換次數(shù)相同，因此地震數(shù)據(jù)重構(gòu)效率得到提高。

表1 說明了程序運行的環(huán)境配置，以及改進前后所需的總運行時間，可以看出，改進后的MPFI算法的重構(gòu)時間明顯縮短，運行效率顯著提升。

表1 改進前、后重構(gòu)計算時間對比Table 1 Comparison of reconstruction calculation before and after improvement

4 應(yīng)用實例

為了驗證改進后MPFI 算法在實際地震資料中的應(yīng)用，采用某陸上工區(qū)疊前CMP 道集記錄進行測試。由于數(shù)據(jù)量過大，現(xiàn)截取原地震數(shù)據(jù)中某一部分。原始CMP 道集道間距為25 m，共70 道，時間方向共500 個采樣點，時間采樣間隔為2 ms，其波形如圖7（a）所示。不規(guī)則CMP 道集中空間方向為56 道，其采樣方式為隨機欠采樣，占原始道集采樣的80%（圖7（b））。使用改進后MPFI 算法對隨機欠采樣CMP 道集數(shù)據(jù)進行重構(gòu)，重構(gòu)后其同相軸變密，缺失信息得到補償，有效信號明顯得到改善（圖7（c））。

圖7 實際地震資料應(yīng)用Fig. 7 Application of real seismic data

5 結(jié) 論

（1）改進后的重構(gòu)算法在計算、更新傅里葉系數(shù)時，不依賴大量的正反非均勻離散傅里葉變換，大大減少了迭代中的計算量，重構(gòu)效率相比常規(guī)MPFI 算法得到提高。

（2）使用二維凹陷速度模型與抽稀Marmousi速度模型的正演單炮地震記錄對改進后的重構(gòu)算法進行測試，測試結(jié)果表明改進的MPFI 算法可以將不規(guī)則采樣地震數(shù)據(jù)規(guī)則化，并且在相同條件下改進后算法的重構(gòu)效率更高；將改進后的算法應(yīng)用到實際地震資料中，取得了良好的重構(gòu)效果。