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

（東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西南昌330013）

0 引言

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

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

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

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

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

1 方法原理

1.1 建立選擇函數(shù)

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

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

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

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

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

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

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

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

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

1.2 POCS算法

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

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

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

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

為體現(xiàn)本文閾值參數(shù)的優(yōu)勢(shì)，將其與文獻(xiàn)［28］中的線性閾值與指數(shù)閾值做重建對(duì)比，其重建后的信噪比公式采用，其中x0表示原始模型數(shù)據(jù)，x表示重建結(jié)果。顯然信噪比越高，表明重建效果越好。計(jì)算結(jié)果如圖2所示。可見(jiàn)在取得相同信噪比情況下，本文閾值參數(shù)所用迭代次數(shù)少，計(jì)算效率高；在相同迭代次數(shù)下，本文閾值參數(shù)重建信噪比最高。因此，本文閾值明顯優(yōu)于線性閾值和指數(shù)閾值。

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

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

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

2 數(shù)值模擬

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

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

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

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

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

圖4 兩種重建方法對(duì)比結(jié)果及其f-k頻譜

圖5 M（p）范圍內(nèi)能量分布圖

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 重建結(jié)果及其f-k譜

4 結(jié)論與討論

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

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

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