利用抗假頻凸集投影算法的規(guī)則缺失地震數(shù)據(jù)重建

刁 塑 張 華 張恒琪 張博泓 余 政 龐 洋

（東華理工大學放射性地質與勘探技術國防重點學科實驗室，江西南昌330013）

0 引言

在地震數(shù)據(jù)現(xiàn)場采集中，由于需要考慮客觀存在的障礙物、禁采區(qū)、海洋拖纜羽狀漂移和采集經(jīng)濟成本等因素，地震數(shù)據(jù)沿空間方向通常呈現(xiàn)不規(guī)則分布或規(guī)則缺失分布；在地震數(shù)據(jù)預處理中，必要的剔除廢炮、廢道等操作也會導致地震道缺失。而后續(xù)的常規(guī)數(shù)據(jù)處理，如速度分析、疊加、偏移等，應建立在完整的地震數(shù)據(jù)基礎上。因此，亟待探尋可行的重建完整數(shù)據(jù)的方法。

現(xiàn)行的地震數(shù)據(jù)重建算法可分為四大類：①基于濾波器方法［1］，主要是依據(jù)建立卷積插值濾波器實現(xiàn)數(shù)據(jù)重建；②基于波場延拓算子重建方法［2-3］，如常見的利用Kirchhoff積分算子；③基于快速降秩插值方法［4-6］，即是將插值問題轉化為圖像填充問題；④基于數(shù)學變換方法，如常用的傅里葉變換、拉東變換、小波變換、曲波變換等，先變換到數(shù)字域進行處理，再反變換回時空域［7-14］。本文方法屬第四類，即通過改進算法，在重建的同時消除伴生假頻。

空間均勻網(wǎng)格下缺失的地震數(shù)據(jù)分為規(guī)則缺失和非均勻缺失。一般非均勻缺失數(shù)據(jù)轉換到頻域，它引起的假頻信息轉化為低幅值不相干隨機噪聲，可設置其閾值門限，通過稀疏迭代去噪法進行消除。而規(guī)則缺失數(shù)據(jù)在頻域中產(chǎn)生的假頻與真實頻譜相近，基于傳統(tǒng)的迭代稀疏促進求解方法很難取得高分辨率、高保真、高信噪比的重建結果。為此，人們提出了較多數(shù)據(jù)重建方法以壓制假頻干擾，如傾角掃描法，此類方法需先掃描同相軸傾角［15］，然后沿傾角方向加權、內插地震道；Sinc地震道插值方法［16］是對滿足采樣定理的數(shù)據(jù)做加密采樣點處理，不需獲取同相軸傾角，運算速度快，但無法正確內插含有空間假頻的地震道；小波變換地震道插值方法是利用小波變換的時頻分析特性［17］，多級重構實現(xiàn)地震道插值，但該插值算法較復雜，且精度較低。

對于不規(guī)則空間帶限二維地震數(shù)據(jù)，結合非均勻快速傅里葉變換和最小二乘反演算法，都取得了較好的重建效果［18-23］；此方法也可應用于規(guī)則欠采樣數(shù)據(jù)中，但其抗假頻效果欠佳，且不能有效區(qū)分真實頻譜與假頻信息。凸集投影（Projection onto convex sets，POCS）算法對于隨機的和規(guī)則的欠采樣數(shù)據(jù)都能重建出完整數(shù)據(jù)［24-25］，但對規(guī)則缺失較嚴重數(shù)據(jù)，此方法去假頻效果不甚理想。在f-k域由規(guī)則采樣引起的假頻幅值與真實頻譜相近，為了進行區(qū)分，需對整個頻域范圍內所有頻率進行搜尋。在信號分離和插值方面，與τ-p域相比，f-k域簡捷有效，主要是因為τ-p重建方法計算量大，且重建效果不如f-k域方法，尤其對低信噪比數(shù)據(jù)［26］。同時，f-k域方法可利用所有期望頻率的信息對缺失數(shù)據(jù)施行穩(wěn)定的插值和去噪。

本文采用傅里葉變換，在傳統(tǒng)POCS方法中引入Naghizade［27］提出的選擇函數(shù)，以去除規(guī)則欠采樣引起的假頻信息。將此方法應用于理論數(shù)據(jù)和實際資料處理，均取得了較理想效果。

1 方法原理

1.1 建立選擇函數(shù)

首先對所有頻率或在任意頻率范圍內進行掃描搜索，以確定主要能量的斜率或路徑。

設d（t，x）為原始地震數(shù)據(jù)，D（f，k）為其相應的f-k域振幅譜。為了簡化，將頻率和波數(shù)做歸一化處理，令時間軸和距離軸上采樣間隔都為1，且f-k頻譜關于頻率軸對稱，只需考慮正頻率的頻譜，得到歸一化頻率范圍f∈（0，0.5）、歸一化波數(shù)范圍k∈（－0.5，0.5）（圖1）。設定每條角度射線起始位置為頻域原點（f，k）∈（0，0），從圖1看出，可沿著角度射線方向求取射線上振幅能量之和。

圖1 f-k域角度射線掃描示意圖

則每條射線所掃描的能量值為

式中：p為角度射線掃描斜率，其取值范圍依據(jù)不同數(shù)據(jù)的頻譜分布而定；n為歸一化頻率f的索引；算子表示沿數(shù)值較小方向取整數(shù)值。在某一區(qū)間，若p值對應的函數(shù)M值大于鄰近值，則稱M為該范圍對應的峰值。而在實際處理過程中，對于諸多峰值，可通過設置閾值的方法，識別函數(shù)中的幾個較主要峰值，代表主要反射波同相軸的振幅。假設已識別L個峰值，如p1，p2，…，pL。

第二步，依據(jù)函數(shù)M（p）的峰值和掃描斜率范圍，初始化一個零矩陣H，其維度與數(shù)據(jù)頻譜維度相同，目的是把p1，p2，…，pL射線上反射波能量準確分配到矩陣H中（令L個角度射線上的值為1，其余各處為0）。則有

式中：j為角度射線p的索引；n、f、k同式（1）。

最后，考慮到鄰近各角度射線兩側的能量應為同相軸能量，為拾取更完整同相軸信息，將所得H與一維方脈沖函數(shù)B（1，Lb）做卷積，沿波數(shù)軸方向拓寬每個角度射線，最終得出選擇函數(shù)

式中：Lb是函數(shù)B沿波數(shù)軸方向上的長度，代表各個角度射線掃描的寬度，且函數(shù)B沿頻率軸方向無值。從式（2）和式（3）可知，在選擇函數(shù)W中，任何大于1的值都設定為1，小于1的值都設為0。

1.2 POCS算法

傳統(tǒng)POCS重建方法的迭代表達式為

式中：di（t，x）為第i次迭代所得重建數(shù)據(jù)；I為單位矩陣；yobs為原始數(shù)據(jù)，即d0（t，x）；F和FT分別為傅里葉正、反變換算子；N為滿足精度要求時最小迭代次數(shù)；S（t，x）為采樣矩陣；Ti為硬閾值算子，滿足

式中：ci表示第i次迭代得到的傅里葉系數(shù)，滿足ci＝Fdi（t，x）；λi表示第i迭代的閾值，本文選用的閾值［28］為

式中：Cmax為傅里葉系數(shù)最大值；ε為靠近零的正值，不同數(shù)據(jù)ε值有所差別，與數(shù)據(jù)中噪聲能量有關，計算中通過多次測試人工選取。

為體現(xiàn)本文閾值參數(shù)的優(yōu)勢，將其與文獻［28］中的線性閾值與指數(shù)閾值做重建對比，其重建后的信噪比公式采用，其中x0表示原始模型數(shù)據(jù)，x表示重建結果。顯然信噪比越高，表明重建效果越好。計算結果如圖2所示?？梢娫谌〉孟嗤旁氡惹闆r下，本文閾值參數(shù)所用迭代次數(shù)少，計算效率高；在相同迭代次數(shù)下，本文閾值參數(shù)重建信噪比最高。因此，本文閾值明顯優(yōu)于線性閾值和指數(shù)閾值。

圖2 不同閾值函數(shù)重建結果的信噪比曲線

將選擇函數(shù)W代入到傳統(tǒng)的POCS迭代重建過程中，即把式（3）代入式（4）中，最終的抗假頻POCS重建算法為

選擇函數(shù)W能消除規(guī)則欠采樣產(chǎn)生的假頻干擾，使每次迭代過程中有效波能量聚焦，因此修改后的POCS算法能成功地重建規(guī)則缺失道。

2 數(shù)值模擬

采用主頻為40Hz雷克子波合成81道地震記錄（圖3a），每道501個采樣點，該記錄總共有3個地震反射同相軸，采樣間隔為2ms，道距為5m。

圖3 原始模型采樣及其f-k譜

為了更清晰地驗證本文抗假頻重建方法的可行性，本文通過對比重建后的振幅譜觀察假頻是否消除，并對比重建后的信噪比。首先，對模型數(shù)據(jù)（圖3a）進行50%規(guī)則欠采樣，即去除偶數(shù)道數(shù)據(jù)，只保留奇數(shù)道數(shù)據(jù)，得到欠采樣數(shù)據(jù)（圖3c）及其頻率譜（圖3d）。可明顯看到，在對應f-k域中（圖3d）因規(guī)則欠采樣導致的假頻與真實頻譜（圖3b）較接近。因此，在之后的重建過程中，需消除假頻的影響，這樣最終才可得到完整而規(guī)則的地震數(shù)據(jù)。

為驗證本文方法的有效性，首先采用傳統(tǒng)傅里葉變換和POCS算法進行重建，選擇迭代次數(shù)N＝20，閾值函數(shù)中的ε值為0.001。觀察重建的結果（圖4a）及其對應的振幅譜（圖4b），可見圖4a中偶數(shù)道上有效波能量并未得到任何恢復、在圖4b中假頻仍然存在，且重建結果的信噪比為2.7dB，說明傳統(tǒng)方法無法重建缺失地震道。

再采用本文抗假頻的POCS算法進行重建，選定與傳統(tǒng)重建相同的參數(shù)和迭代次數(shù)。依據(jù)方法原理，選擇函數(shù)以圖3d中的原點為起點，對頻域中的能量進行掃描，保留3個能量峰值對應的傾角，并在波數(shù)方向拓寬掃描范圍，利用式（7）得到無假頻的重建數(shù)據(jù)。觀察此時所得重建地震數(shù)據(jù)（圖4c），其偶數(shù)道得到很好的重建，同相軸連續(xù)完整；且對應振幅譜（圖4d）上也完全看不到假頻信息。對比圖4d與圖3b，表明本文所提抗假頻POCS算法能有效地重建出規(guī)則缺失的地震道數(shù)據(jù)，且消除了規(guī)則欠采樣所導致的假頻，重建結果的信噪比增至14.5dB。

圖4 兩種重建方法對比結果及其f-k頻譜

圖5 M（p）范圍內能量分布圖

詳細展示本文方法如何得到圖4c重建結果。首先，利用式（1）分別求取原始數(shù)據(jù)和規(guī)則缺失50%數(shù)據(jù)對應的M（p）峰值（圖5），掃描斜率p的范圍設為－8～8，對比圖5a與圖5b，其主要能量峰值對應的p值一致。依據(jù)峰值范圍，設定合適閾值，保留對應同相軸頻譜的3個峰值，并分配到矩陣H中，再利用式（3）建立選擇函數(shù)，結果如圖6。其中大于閾值區(qū)域的值為1，小于閾值區(qū)域的值為0?？芍獔D6b中選擇函數(shù)與圖6a選擇函數(shù)能量分布一致，說明選擇函數(shù)識別出了真實頻譜。然后將該選擇函數(shù)引入到POCS方法中，使得在每次POCS算法迭代過程中，該選擇函數(shù)可很精確地拾取真實頻譜，再參與下一次迭代計算，使得原始信號頻譜能很快聚焦。

為更進一步驗證本文方法的可行性，針對原始模型數(shù)據(jù)進行25%規(guī)則采樣，即得到規(guī)則缺失75%數(shù)據(jù)（圖7a）及其對應的頻譜（圖7b），可見真實頻譜與假頻信息相互交錯纏繞，且能量非常相似，傳統(tǒng)方法很難去除假頻信息。采用本文方法，選定合適的迭代次數(shù)和閾值參數(shù)ε值，得到重建結果（圖7c）及對應f-k譜（圖7d）。該重建結果的信噪比為8.7dB。

分析重建結果的f-k譜（圖7d），本文引入的選擇函數(shù)在重建過程中有效地壓制了假頻干擾，保留有效波能量。同時從建立的選擇函數(shù)（圖6c）中可見，同相軸能量對應的峰值與圖6a一致。但相比之下，其75%規(guī)則缺失數(shù)據(jù)重建效果較差，說明當缺失較為嚴重時，重建結果則難以達到理想效果，周邊已知的地震道不足以恢復大范圍連續(xù)缺失的地震道，因此需發(fā)展高維的重建方法，從另一維度進行重建。同時需說明的是，因基于傅立葉變換的二維重建方法只能處理近似線性同相軸，而對于非線性同相軸則需做分窗口重建，使其滿足近似線性同相軸的要求。

圖6 f-k域選擇函數(shù)值

圖7 25%規(guī)則采樣重建結果及其f-k譜

3 實際數(shù)據(jù)算例

選取實際采集地震數(shù)據(jù)驗證本文方法的適用性。圖8為M區(qū)海上二維原始地震數(shù)據(jù)及其振幅譜，道間距為25m，采樣率為4ms，180道接收。本文方法基于傅里葉變換，只適應線性或近似線性同相軸地震數(shù)據(jù)重建，而對于非線性同相軸，通常采用分窗口重建。

為滿足處理要求，從原始數(shù)據(jù)中截取第100～第140道數(shù)據(jù)，時間軸方向截取0.8～1.6s（圖9a、圖9b）。對該窗口數(shù)據(jù)進行50%的規(guī)則欠采樣（圖9c），可見規(guī)則采樣引起的假頻信息（圖9d）與真實頻譜（圖9b）非常相似。

圖8 現(xiàn)場實測數(shù)據(jù)（a）及其f-k譜（b）

圖9 截取數(shù)據(jù)和規(guī)則缺失數(shù)據(jù)及其f-k譜

針對該實際數(shù)據(jù)體，本次選取選擇函數(shù)中計算傾角p值的范圍是－12～12。選擇函數(shù)以圖9d中原點為起點，對頻域中的能量進行掃描，掃描的同相軸能量分布如圖10所示，此掃描范圍可最大程度地保留反射波同相軸能量、去除規(guī)則采樣產(chǎn)生的假頻信息。同時，迭代次數(shù)和閾值函數(shù)ε值的選取也影響重建效果。

由最終重建結果（圖11a）及對應f-k譜（圖11b）得到重建結果誤差（圖11c）及其頻譜（圖11d），可見規(guī)則采樣導致的假頻信息得到較徹底壓制，并保留了主要反射波同相軸能量，重建結果的信噪比為12.78d B。

圖10 實測數(shù)據(jù)的規(guī)則缺失50%數(shù)據(jù)f-k域選擇函數(shù)

圖11 重建結果及其f-k譜

4 結論與討論

本文結合傅里葉變換與POCS重建算法，研究了傾角掃描選擇函數(shù)，該函數(shù)對全頻段的能量軸進行掃描，確定出主要傾角；并把該函數(shù)引入傳統(tǒng)的重建算法，重建線性或近似線性的同相軸數(shù)據(jù)；在此基礎上，將該方法應用于理論和實際數(shù)據(jù)，重建出規(guī)則缺失的地震數(shù)據(jù)并壓制了假頻干擾。

根據(jù)對傾角掃描選擇函數(shù)的分析，在傳統(tǒng)傅里葉變換和POCS重建算法的應用中，只可重建線性或近似線性的地震數(shù)據(jù)。然而，實際數(shù)據(jù)中缺失的反射波同相軸往往是非線性的，此時需對非線性地震數(shù)據(jù)做分窗口重建，使得每個窗口內地震數(shù)據(jù)滿足近似線性同相軸。

還可嘗試選取另一種數(shù)學變換替代傅里葉變換，并結合傾角掃描選擇函數(shù)，重建非線性同相軸缺失數(shù)據(jù)。這將是后續(xù)研究課題。