井洪亮 張少華 方云峰* 馬光凱 張 茜 楊雪霖

(①東方地球物理公司物探技術(shù)研究中心，河北涿州 072751； ②東方地球物理公司，河北涿州 072750； ③東方地球物理公司研究院處理中心，河北涿州 072751)

0 引言

自20世紀(jì)70年代起，Radon變換引入地球物理勘探，經(jīng)過幾十年的發(fā)展，已逐步用于地震資料處理的諸多方面，如地震道重建、VSP上、下波場分離及多次波壓制等[1-6]。尤其在多次波壓制方面，Radon變換因其應(yīng)用簡便、操作靈活，受到業(yè)界的廣泛關(guān)注。Hampson[7]首次利用τ-p變換壓制多次波，但由于變換域能量聚焦不足，存在嚴(yán)重的能量“拖尾”現(xiàn)象，多次波壓制效果并不理想； 隨后，Kostov[8]、Sacchi[9]提出了高分辨率Radon變換，有效地解決了“剪刀狀”能量發(fā)散問題，提高了多次波壓制精度； Wood等[10]推出一種混合域高分辨率Radon變換，充分利用時間域的高分辨率和頻率域可解耦的特點，進(jìn)一步提高了Radon域內(nèi)解的稀疏度和分辨率，但計算成本偏高； Abbad等[11]基于奇異值分解技術(shù)提出了頻率域快速求解算法，即λ-f域Radon變換，通過引入變量λ消除變換算子與頻率的關(guān)系，顯著提高了計算效率。

中國學(xué)者也深入研究了Radon變換多次波壓制方法。劉喜武等[12]利用稀疏約束共軛梯度算法實現(xiàn)高分辨率Radon變換，并成功用于多次波壓制； 王維紅等[13]提出了一種加權(quán)拋物Radon變換，利用變化的權(quán)系數(shù)提高抗假頻能力； 熊登等[14]利用混合域Radon變換在時間域進(jìn)行稀疏約束，提高了Radon域內(nèi)能量團(tuán)的分辨率； 石穎等[15-18]基于二維λ-f域Radon變換優(yōu)化和改進(jìn)地震道重建、多次波壓制和波場分離等方法。近年來，在Radon變換的保幅處理方面也取得一定突破。薛亞茹等[19-22]基于正交多項式變換和常規(guī)Radon變換的不同特征，將二者結(jié)合形成高階Radon變換，通過構(gòu)建多階Radon子區(qū)間，使每個子區(qū)間保留不同的振幅變化信息，從而具有AVO特性； 唐歡歡等[23]將高階二維Radon變換成功拓展到三維，并取得很好的保幅效果。王亮亮等[24]實現(xiàn)了快速Radon變換算法； 張全等[25]提出了CPU-GPU異構(gòu)并行算法，加速比達(dá)到約30倍，提高了Radon變換計算效率。

在實際資料處理中，二維拋物Radon變換壓制多次波的方法已經(jīng)相當(dāng)成熟，但是沒有考慮地震子波的空間傳播特性，忽略了波場隨空間曲率的變化信息，在處理寬方位三維地震數(shù)據(jù)時，往往多次波壓制效果不佳，存在一定的多次波殘余。為此，本文基于二維拋物Radon變換，借鑒Abbad等[11]的方法，深入研究了常規(guī)三維拋物Radon變換，并首次提出了λ-f域三維拋物Radon變換多次波壓制方法。結(jié)果表明：常規(guī)三維拋物Radon變換多次波壓制方法的計算效率過低，不利于大規(guī)模實際應(yīng)用；λ-f域三維拋物Radon變換多次波壓制方法重新定義了變換算子的構(gòu)建方式，消除了算子對地震數(shù)據(jù)頻率的依賴，使整個變換過程僅需計算一次變換算子及其廣義逆算子，大大減少了矩陣求逆運算的次數(shù)，明顯提高了計算效率。同時，根據(jù)一次波和多次波映射到λ-f域的空間分布特征，設(shè)計了三維椎體濾波器，使一次波能量投影到保留區(qū)，多次波能量投影到切除區(qū)，避免了人工定義切除區(qū)域的過程，更便于實際應(yīng)用。因此，針對海量三維地震數(shù)據(jù)處理，本文所提方法可以極大地節(jié)省人工成本。

1 方法原理

1.1 常規(guī)三維拋物Radon變換原理

二維拋物Radon變換是沿著拋物線路徑對地震數(shù)據(jù)積分、求和的過程。類似地，三維拋物Radon變換的物理意義是沿著拋物面對地震數(shù)據(jù)求和處理，其正、反變換的離散形式分別為

(1)

(2)

式中：d(t,x,y)為三維地震數(shù)據(jù)體，t為時間坐標(biāo)，x、y分別為橫、縱向炮檢距；m(qx,qy,τ)為變換域的模型空間，qx和qy分別為橫、縱向離散曲率參數(shù)，τ為變換域零炮檢距截距時間。

對式(1)、式(2)兩邊進(jìn)行傅里葉變換，分別得到頻率域的正、反變換

(3)

(4)

通常情況下，應(yīng)用Radon變換壓制多次波，多是對正變換的模型空間進(jìn)行處理，再反變換回時空域，而Radon域能量團(tuán)的收斂性直接影響多次波去除效果。因此，本文著重推導(dǎo)了三維拋物Radon正變換的實現(xiàn)過程。根據(jù)式(4)可知，在頻率域內(nèi)，不同頻率切片的三維地震數(shù)據(jù)具有解耦性，對某個頻率，其矩陣簡化形式為

(5)

式中：Lqx=e-ifqxx2為與x有關(guān)的算子；Lqy=e-ifqyy2為與y有關(guān)的算子。采用最小二乘方法，使式(5)的目標(biāo)函數(shù)最小化，即

(6)

分別求取矩陣Lqx和Lqy的逆，可以得到Radon域數(shù)據(jù)M。因此在最小二乘意義下，式(5)的解為

(7)

上標(biāo)“H”表示取共軛。為了防止奇異值對矩陣求逆過程產(chǎn)生異常影響，引入阻尼因子μ(通常取值為0.01～1.00)，主要用于增強矩陣求逆的穩(wěn)定性。最終解的表達(dá)式為

(8)

1.2 λ-f域三維拋物Radon變換原理

針對常規(guī)三維拋物Radon變換計算效率偏低的問題，本文借鑒Abbad等[11]提出的改進(jìn)二維Radon變換方法，重新定義變換算子的構(gòu)建方式，引入兩個數(shù)學(xué)運算變量λx=qxf和λy=qyf(其物理意義為曲率與頻率的乘積)，可將常規(guī)三維拋物Radon變換的f-qx-qy域轉(zhuǎn)換到一個全新的λx-λy-f域。該法繼承了二維λ-f域Radon變換思路，創(chuàng)新實現(xiàn)了三維算法，故簡稱為λ-f域三維拋物Radon變換方法。

針對任意頻率分量，將變量λx和λy引入變換公式，則改進(jìn)三維拋物Radon正、反變換為

(9)

(10)

類似于常規(guī)三維拋物Radon變換方法，針對某個頻率切片，將式(10)的矩陣形式簡化為

(11)

式中Lλx、Lλy為重構(gòu)的變換算子，具體形式為

式中：m和p分別為三維地震道集的縱測線和非縱測線的道數(shù)；n和k為λ-f域的三維模型空間的維度。通過對變換算子Lλx、Lλy重構(gòu)，Lλx由炮檢距xm和λxn(λx的離散值)計算獲得，Lλy由炮檢距yp和λyk(λy的離散值)計算獲得，從而消除了對頻率f的依賴。

采用最小二乘求解方式得到λ-f域三維拋物Radon變換的解

(12)

以模擬的三維地震記錄為例，對比、分析常規(guī)三維拋物Radon正變換及λ-f域三維拋物Radon正變換的譜分布特征及聚焦性(圖1)。利用常規(guī)三維拋物Radon變換將合成數(shù)據(jù)(圖1a)變換到三維Radon域，并抽取四條同相軸對應(yīng)最大能量團(tuán)的垂向(圖1b)和橫向(圖1c)切片。理論上，當(dāng)截距時間和介質(zhì)速度橫向變化不大時，同相軸映射到變換域應(yīng)呈“點”狀分布，但受采集孔徑限制，形成空間截斷效應(yīng)，能量團(tuán)存在“拖尾”現(xiàn)象，能量泄露造成收斂不足(圖1b、圖1c)。對圖1a進(jìn)行λ-f域三維拋物Radon正變換，抽取變換域內(nèi)不同位置的縱向切片(圖1d)，由于地震數(shù)據(jù)中每條同相軸的曲率qx、qy基本固定，因此λx和λy與頻率f成正比，四條拋物同相軸映射到λ-f域呈以不同曲率為斜率的“窄直線”狀分布，且當(dāng)計算頻率極低時，λx和λy趨于0，即所有直線近似過原點。抽取λ-f域內(nèi)不同頻率的橫向切片(圖1e)發(fā)現(xiàn)，四條同相軸在λ-f域能量團(tuán)較收斂，且隨著頻率增加，更容易區(qū)分能量團(tuán)的相對位置關(guān)系，便于后續(xù)的多次波分離。

圖1 合成炮記錄Radon正變換結(jié)果

1.3 λ-f域離散采樣條件

為了避免λ-f域三維拋物Radon變換出現(xiàn)假頻而影響算法精度，需要將所有同相軸對應(yīng)的Nλ(λx和λy的個數(shù))值準(zhǔn)確映射到變換域。若Nλ不足使同相軸能量映射不完整，造成能量泄露； 若Nλ過大會引入一定的隨機噪聲，降低資料信噪比。選取合適的Nλ值是必要的，決定該值的兩個重要參數(shù)是采樣間隔和掃描范圍。

根據(jù)常規(guī)三維拋物Radon變換的曲率選取準(zhǔn)則，由λx=qxf、λy=qyf可知，λ-f域重采樣間隔Δλx和Δλy應(yīng)滿足

(13)

(14)

式中：xmax、xmin和ymax、ymin分別為橫向和縱向最大、最小炮檢距； Δqx、Δqy分別為橫、縱向曲率采樣間隔。

根據(jù)式(11)，重構(gòu)的Radon算子Lλx、Lλy的相位誤差必須滿足

(15)

(16)

因此，根據(jù)qx、qy及地震數(shù)據(jù)最大頻率fmax，可以確定參數(shù)λx和λy的掃描范圍

(17)

(18)

由式(13)～式(18)可見，Δλx、Δλy及λx、λy的掃描范圍僅與縱、橫向最大、最小炮檢距以及空間采樣間隔Δx、Δy有關(guān)，而與頻率無關(guān)。

三維拋物Radon變換采樣間隔Δλx、Δλy和λx、λy的掃描范圍類似于常規(guī)二維Radon變換(圖2)[11]，即：在任意頻率范圍內(nèi)，所有同相軸的曲率q的掃描范圍相同(圖2a)，對于低頻同相軸而言，在有限的曲率范圍內(nèi)，頻率過低可能造成采樣不足，無法滿足采樣定理要求，出現(xiàn)低頻同相軸能量泄露、聚焦不穩(wěn)定的問題； 基于λ-f域拋物Radon變換方法，參數(shù)λ的采樣間隔和掃描范圍僅與炮檢距、道間距有關(guān)，在一定的λ取值范圍內(nèi)，曲率和頻率呈反比，即頻率越低，曲率范圍越大，低頻區(qū)的覆蓋區(qū)范圍更寬，增加了低頻段到變換域的完整映射幾率(圖2b)。

圖2 參數(shù)變量掃描范圍示意圖

1.4 λ-f域三維椎體濾波壓制多次波

在實際資料處理中，利用拋物Radon變換壓制多次波首先要對地震數(shù)據(jù)做動校正(NMO)，然后根據(jù)一次波和多次波能量在變換域內(nèi)分布位置的不同，通過手動方式定義切除區(qū)域，從而分離一次波和多次波能量團(tuán)。對于海量三維地震數(shù)據(jù)而言，采用常規(guī)三維拋物Radon變換人為定義三維切除區(qū)，計算成本偏高。本文充分利用NMO后的一次波和多次波能量在λ-f域呈“窄直線”狀分布特點，提出了一種椎體濾波方法，可以有效分離一次波和多次波能量，避免了人工定義切除區(qū)域的過程。實際Radon變換壓制多次波的經(jīng)驗表明，為了滿足淺層非多次波干擾區(qū)地震數(shù)據(jù)的保真性，通常在變換域內(nèi)提取多次波能量，將其反變換回時空域，得到多次波模型數(shù)據(jù)，再采用自適應(yīng)相減法從原始數(shù)據(jù)中去除多次波模型數(shù)據(jù)，以降低對有效波的傷害?；诖怂悸?，在設(shè)計三維椎體濾波器時，以提取多次波能量為目標(biāo)。設(shè)地震數(shù)據(jù)d(x,y,t)的λ-f域三維拋物Radon正變換為M(λx,λy,f)，定義椎體濾波器為G(λx,λy,f)，不同曲率的同相軸會映射到不同的椎體濾波區(qū)(圖3)。假設(shè)一次波同相軸經(jīng)NMO后均被校平，視速度v趨近于∞，由公式q≈1/(2t0v2)可知，一次波的曲率近似為0，則一次波能量從椎體頂點附近垂向分布到椎體內(nèi)部的切除區(qū)域。地震數(shù)據(jù)中多次波同相軸經(jīng)NMO后變?yōu)榫哂幸欢ㄇ?范圍為qmin～qmax)的拋物面，多次波能量從頂點附近呈直線輻射至椎體保留區(qū)域，則

圖3 三維數(shù)據(jù)體λ-f域椎體濾波區(qū)域示意圖

(19)

(20)

2 數(shù)據(jù)實驗

2.1 三維模型數(shù)據(jù)測試

為了檢驗λ-f域三維拋物Radon變換方法的效果，以三維正演速度模型(圖4)為例進(jìn)行多次波壓制試算。圖5展示了圖4模型數(shù)據(jù)多次波壓制效果。由圖可見：①NMO三維模型數(shù)據(jù)的一次波同相軸均被校平，多次波同相軸仍為拋物形態(tài)(圖5a)。②常規(guī)三維拋物Radon正變換抽取一階自由表面多次波能量最強位置，當(dāng)一次波和多次波的截距時間較小時，由于存在能量團(tuán)“拖尾”現(xiàn)象，部分能量交織在一起，很難分離多次波能量，在定義切除線過程中會損傷一次波能量或殘留多次波能量(圖5b)。③λ-f域三維拋物Radon正變換結(jié)果的一次波能量垂直分布，而多次波能量具有一定斜率，因此根據(jù)空間位置可分離一次波和多次波(圖5c)。④常規(guī)三維拋物Radon變換和λ-f域三維拋物Radon變換均能有效提取多次波模型數(shù)據(jù)，但是在近炮檢距位置都?xì)埩粢淮尾?圖5d、圖5e紅色箭頭處)。這是由于近炮檢距位置下層的一次波與上層的一階自由表面多次波的NMO時差較小，很難有效分離所致，這也是Radon類方法在近炮檢距位置多次波壓制效果不佳的原因。⑤常規(guī)三維拋物Radon變換多次波壓制結(jié)果明顯殘留多次波，這是由于常規(guī)三維拋物Radon變換的能量團(tuán)收斂性差，在提取多次波能量過程中，為了減少一次波能量損傷，出現(xiàn)部分多次波能量泄露(圖5f白色箭頭處)。⑥由于在λ-f域易于分離一次波和多次波能量，因此λ-f域三維拋物Radon變換較干凈、徹底地壓制了多次波(圖5g)。

圖4 三維正演速度模型

為了進(jìn)一步分析常規(guī)三維拋物Radon變換和λ-f域三維拋物Radon變換的多次波壓制精度，隨機抽取圖5a、圖5f和圖5g同一位置的地震道(圖6)，可見，后者的多次波壓制精度高于前者。

圖5 圖4數(shù)據(jù)多次波壓制效果

圖6 由圖5a、圖5f和圖5g隨機抽取的同一位置的地震道對比

表1展示了利用不同方法壓制多次波的計算效率，可見在采用相同的計算平臺和計算參數(shù)情況下，對同一個節(jié)點的單線程運算，常規(guī)三維拋物Radon變換壓制多次波方法的計算耗時(不包括人為定義切除過程)是λ-f域三維拋物Radon變換的8倍以上。因此，針對海量地震數(shù)據(jù)去噪，后者的計算高效性更明顯，具有實用價值。

表1 利用不同方法壓制多次波的計算效率

理論上，二維拋物Radon變換壓制多次波方法存儲變換算子的矩陣維度不大，對于每個頻率，重新計算變換算子的耗時較短。因此，λ-f域二維拋物Radon變換壓制多次波的效率提升并不明顯。常規(guī)三維拋物Radon變換壓制多次波方法存儲變換算子的維度大，對于每個頻率切片，需要重復(fù)進(jìn)行大量的矩陣運算，尤其是求取廣義逆運算，會嚴(yán)重制約計算效率。因此，λ-f域三維拋物Radon變換壓制多次波充分體現(xiàn)了計算效率優(yōu)勢。

2.2 實際資料測試

以海洋拖攬數(shù)據(jù)為例測試多次波壓制效果。圖7為λ-f域三維拋物Radon變換壓制海洋單炮記錄多次波效果。由圖可見：①NMO單炮記錄中廣泛存在一階和二階自由表面多次波，而且能量級別較大，部分有效波同相軸被掩蓋(圖7a)。②λ-f域三維拋物Radon正變換結(jié)果的一次波能量在λx=0、λy=0區(qū)域附近幾乎垂直分布，根據(jù)多次波的最大和最小視速度范圍定義椎體濾波區(qū)域，可以分離有效波和多次波(圖7b)。③λ-f域三維拋物Radon變換有效恢復(fù)了能量強、曲率大的多次波(圖7c)。④經(jīng)多次波壓制，圖7a的強能量多次波均得到衰減，同時顯現(xiàn)了部分能量較弱的一次波(圖7d紅色箭頭處)。

圖7 λ-f域三維拋物Radon變換壓制海洋單炮記錄多次波效果

圖8為三維共炮檢距剖面多次波壓制效果。由圖可見：原始數(shù)據(jù)存在強能量多次波，且地層具有一定傾角(圖8a箭頭標(biāo)注處)； 常規(guī)三維拋物Radon變換有效壓制了大部分較平緩的一階、二階自由表面多次波，但在一些構(gòu)造傾角較大位置仍殘留多次波模型包含兩層水平層狀介質(zhì)，利用有限差分正演方法記錄了上層的一階和二階自由表面多次波。模擬觀測系統(tǒng)為規(guī)則網(wǎng)格，其中，縱、橫向各500個網(wǎng)格點，網(wǎng)格間距為10m(圖8b)； 與圖8b相比，λ-f域三維拋物Radon變換干凈、徹底地壓制了多次波，尤其在一次波和多次波同相軸較接近的位置(t=3.3s)，且未損傷有效信號(圖8c)。在相同的計算環(huán)境和計算參數(shù)情況下，常規(guī)三維拋物Radon變換處理681炮地震數(shù)據(jù)的耗時為684min，而λ-f域三維拋物Radon變換的耗時為71min，即后者的加速比提高了約9倍。

然而，盡管λ-f域三維拋物Radon變換方法的計算精度和效率高于常規(guī)三維拋物Radon變換方法，但在圖8的紅色方框位置仍殘留多次波。為了進(jìn)一步驗證，抽取該位置附近的疊前CMP道集速度譜(圖9)?？梢姡瑑煞N方法均存在明顯的低速能量團(tuán)，尤其在4.0s和4.7s位置附近。分析表明，在淺層(2s位置)存在繞射點，進(jìn)而誘發(fā)強能量繞射多次波，其同相軸呈近似雙曲線形態(tài)，不滿足拋物Radon變換的假設(shè)條件，因此較難徹底壓制。綜上所述，利用λ-f域三維拋物Radon變換多次波壓制方法去除海洋地震資料表面多次波的效果較好。

圖8 三維共炮檢距剖面多次波壓制效果

圖9 疊前CMP道集速度譜

3 結(jié)束語

λ-f域三維拋物Radon變換是對常規(guī)三維 Radon變換的改進(jìn)，本文將其成功應(yīng)用于三維地震資料多次波壓制，理論模型和實際數(shù)據(jù)試算表明方法具有實用價值，適合于海量地震數(shù)據(jù)去噪，獲得以下認(rèn)識。

(1)與常規(guī)三維 Radon變換方法相比，λ-f域三維拋物Radon變換方法壓制表面多次波更徹底，同相軸在λ-f域呈“窄直線”狀分布，緩解了能量團(tuán)“拖尾”現(xiàn)象，在多次波切除過程中降低了空間截斷效應(yīng)的影響。

(2)λ-f域三維拋物Radon變換通過對變換算子進(jìn)行重構(gòu)，消除了變換算子對頻率的依賴，有效減少了矩陣求逆運算次數(shù)，較常規(guī)三維 Radon變換方法的計算效率提高了約8倍以上。

(3)利用λ-f域三維拋物Radon變換方法壓制多次波時，需要根據(jù)多次波同相軸的視速度確定椎體濾波器的母線長度，即視速度越大，濾波器的切除線越向有效信號能量靠近。為了避免損傷有效信號，建議實際應(yīng)用時不易選取過大的最大視速度。

尚需指出，λ-f域三維拋物Radon變換方法壓制多次波時仍然存在Radon類方法的弊端，即當(dāng)一次波和多次波速度差異較小時，在近炮檢距位置無法得到較理想的多次波壓制效果。