唐歡歡，毛偉建*

1 中國(guó)科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院計(jì)算與勘探地球物理研究中心，武漢 430077 2 大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢 430077

0 引言

抗噪性能良好的Radon變換被廣泛應(yīng)用于圖像分析和信號(hào)重構(gòu)等方面，它是一種積分變換，在地震勘探數(shù)據(jù)處理領(lǐng)域，根據(jù)積分路徑的不同，Radon變換可分為線型、雙曲型以及拋物型(Thorson,1978;Hampson,1986;Yilmaz,1989).線性Radon變換(Linear Radon,L_Radon)又稱(chēng)為傾斜疊加，它是將時(shí)空域的同相軸振幅沿著偏移距方向以一定的斜率進(jìn)行疊加運(yùn)算，因此線性同相軸在線性Radon變換域的系數(shù)表現(xiàn)為一個(gè)強(qiáng)能量點(diǎn)，對(duì)應(yīng)的縱橫坐標(biāo)分別表示同相軸的時(shí)間截距及速度信息；而隨機(jī)噪聲經(jīng)過(guò)線性疊加后在變換域分布無(wú)規(guī)律且能量極弱，該特征有利于隨機(jī)噪聲壓制(Turner,1990;鞏向博等，2009；Sabbione et al.,2015)、面波壓制(李衛(wèi)忠等，1998；Wang,1999;張智等，2017)以及波場(chǎng)分離(Moon et al.,1986;王維紅等，2006；李志娜等，2014；Lin and Sacchi,2020).但在油氣勘探地震數(shù)據(jù)中，有效信號(hào)為反射波，其時(shí)距曲線為雙曲型，經(jīng)過(guò)線性Radon變換后，其變換系數(shù)分布并不集中，失去了變換的稀疏性特征且物理意義不明顯.Hampson(1987)將雙曲Radon變換(Hyperbolic Radon,H_Radon)引入到地震數(shù)據(jù)重建中，但雙曲Radon變換在頻率域不解耦，只能在時(shí)間域進(jìn)行計(jì)算導(dǎo)致效率較低(Spitzer et al.,2001;石穎和王維紅，2012).雙曲型同相軸經(jīng)過(guò)部分動(dòng)校正后可近似為拋物型，而拋物Radon變換(Parabolic Radon,P_Radon)在頻率域是解耦的，且在頻率域具有共軛對(duì)稱(chēng)性使其計(jì)算效率較高，因此拋物Radon變換被廣泛地應(yīng)用于多次波壓制(Foster and Mosher,1992;熊登等，2009；鞏向博等，2014；劉菲菲等，2017)及插值重建(Sacchi and Ulrych,1995).Sacchi和Ulrych (1995)提出的利用貝葉斯參數(shù)反演求解拋物Radon變換系數(shù)方法提高了拋物Radon變換的分辨率，有利于提高多次波壓制(Zhang and Wang,2006;Song et al.,2015;Sun et al.,2019)和數(shù)據(jù)重建效果(劉喜武等，2004；Wang et al.,2019).Johansen(1995)利用二階正交多項(xiàng)式對(duì)疊后水平同相軸地震數(shù)據(jù)進(jìn)行擬合，正交多項(xiàng)式系數(shù)可以表征地震數(shù)據(jù)振幅隨偏移距變化情況，該方法缺少同相軸走時(shí)特征的相關(guān)參數(shù)，因此該方法針對(duì)的是水平同相軸數(shù)據(jù)，根據(jù)此特征，薛亞茹等(2014)將正交多項(xiàng)式引入到Radon變換中，得到具有保幅特征的2D高階Radon變換；唐歡歡和毛偉健(2014)、唐歡歡等(2018，2020)將其推廣到3D，提出基于3D高階Radon變換的地震數(shù)據(jù)保幅重建方法.在改善積分路徑方面，牛濱華(2001)提出的多項(xiàng)式Radon變換，將積分路徑從線型或拋物型變?yōu)?階多項(xiàng)式曲線，提高了淺層雙曲型同相軸數(shù)據(jù)在Radon域的稀疏程度，該變換的優(yōu)勢(shì)在于不增加計(jì)算量的情況下可以對(duì)復(fù)雜形態(tài)同相軸數(shù)據(jù)進(jìn)行稀疏表示，但不能同時(shí)對(duì)線性、拋物型及復(fù)雜形態(tài)同相軸數(shù)據(jù)進(jìn)行稀疏表示；Trad等(2001)提出了將雙曲Radon變換算子和線性Radon變換算子相結(jié)合的混合Radon變換，該方法缺少這兩種算子的交叉項(xiàng)無(wú)法對(duì)復(fù)雜同相軸進(jìn)行精確擬合，且在時(shí)間域進(jìn)行計(jì)算效率較低；魯娥和李慶春(2013)將其應(yīng)用到地震數(shù)據(jù)去噪中來(lái)，但該變換與積分路徑相關(guān)的變換參數(shù)只有線性項(xiàng)，導(dǎo)致在變換域出現(xiàn)混疊現(xiàn)象.

以上國(guó)內(nèi)外學(xué)者分別從積分路徑的形態(tài)、計(jì)算效率、分辨率以及保幅性方面對(duì)Radon變換進(jìn)行改進(jìn)，但現(xiàn)有的Radon變換應(yīng)用于地震數(shù)據(jù)處理中仍然存在一個(gè)問(wèn)題，即當(dāng)線性同相軸、雙曲型同相軸以及復(fù)雜形態(tài)同相軸同時(shí)存在時(shí)，不同形態(tài)同相軸在Radon域系數(shù)不能同時(shí)達(dá)到最佳稀疏狀態(tài)，只有同相軸形態(tài)和積分路徑一致的數(shù)據(jù)在Radon域系數(shù)才能集中在一點(diǎn)，而其他的系數(shù)則是發(fā)散狀態(tài)，在空間上有一定的延展性，這將降低隨機(jī)噪聲壓制效果以及波型分離精度.

針對(duì)該問(wèn)題，本文提出了多路徑Radon變換(Multi-path Radon,M_Radon)，這里我們把積分路徑僅為線性、雙曲型或拋物型的Radon變換稱(chēng)為單路徑Radon變換，多路徑Radon變換是指對(duì)每個(gè)同相軸都有線型和拋物型的積分路徑，且具有兩種不同類(lèi)型的變換參數(shù)來(lái)與積分路徑相對(duì)應(yīng)，然后利用最小二乘稀疏反演方法對(duì)不同形態(tài)同相軸自適應(yīng)匹配最佳積分路徑，得到在Radon域同時(shí)聚焦的高分辨率系數(shù).該變換有利于提高隨機(jī)噪聲的壓制、面波壓制以及波型分離效果.

1 方法原理

在地震數(shù)據(jù)處理中，經(jīng)典拋物Radon正反變換如公式(1a)、(1b)所示(Hampson,1986;Beylkin,1987)，是廣義Radon變換的一種特殊形式，公式中d為地震數(shù)據(jù)，x為檢波點(diǎn)坐標(biāo)，q為拋物型積分路徑的曲率，即速度平方的倒數(shù),m為數(shù)據(jù)在Radon域系數(shù)；將積分公式中x2變?yōu)閤，如公式(2a)、(2b)所示，其中p為斜率，即速度的倒數(shù)，公式為線性Radon正反變換，也稱(chēng)為傾斜疊加.公式(1a)、(1b)、(2a)、(2b)分別為

(1a)

(1b)

(2a)

(2b)

通過(guò)線性Radon或拋物Radon變換，理論上可以將t-x域的線性同相軸或拋物型同相軸變換為Radon域的一個(gè)點(diǎn)，如圖1所示，圖1a為一合成的線性同相軸，64個(gè)接收道，道間距為20 m，4 ms時(shí)間采樣率，512個(gè)時(shí)間采樣點(diǎn)，合成該數(shù)據(jù)的速度為v=4000 m·s-1，經(jīng)過(guò)線性Radon變換后其Radon域系數(shù)如圖1b所示；圖1c為一拋物型同相軸，速度為v=1800 m·s-1，經(jīng)過(guò)拋物Radon變換后其Radon域系數(shù)如圖1d所示.從圖中可以看出，t-x域的地震數(shù)據(jù)經(jīng)過(guò)合適類(lèi)型的Radon變換后，其Radon域系數(shù)都各自集中在一個(gè)點(diǎn)附近，從圖1b、d上很容易找到這兩點(diǎn)對(duì)應(yīng)的變換參數(shù)分別為p=0.00025 s·m-1，q=3.08×10-7s2·m-2；再利用公式p=1/v和q=1/v2，可以反向計(jì)算出t-x域地震數(shù)據(jù)線性同相軸和拋物型同相軸的速度分別為v=4000 m·s-1和v=1800 m·s-1，與合成模型數(shù)據(jù)的速度一致，即Radon域系數(shù)分布位置具有明確的物理意義.因此，可以利用Radon變換的稀疏性和明確的物理意義特征進(jìn)行隨機(jī)噪聲壓制、波型分離、多次波壓制以及數(shù)據(jù)重建.

圖1 單軸數(shù)據(jù)Radon變換效果(a)t-x域線性同相軸；(b)圖(a)數(shù)據(jù)線性Radon變換結(jié)果；(c)t-x域拋物型同相軸；(d)圖(c)數(shù)據(jù)拋物Radon變換結(jié)果.Fig.1 Radon transform result of single-axis data(a)Linear in-phase axis in t-x domain；(b)Linear Radon transform of data in (a)；(c)Parabolic in-phase axis in t-x domain；(d)Radon transform of data in (c).

但是，在實(shí)際地震數(shù)據(jù)中，直達(dá)波、折射波和面波的時(shí)距曲線為線性同相軸，反射波的時(shí)距曲線為雙曲型(經(jīng)部分NMO后可近似為拋物型)，在使用線性Radon或拋物Radon變換時(shí)，不能使所有波型數(shù)據(jù)在變換域同時(shí)達(dá)到最佳稀疏狀態(tài)，如圖2合成數(shù)據(jù)所示，為t-x域地震數(shù)據(jù)，包含線型和拋物型同相軸；其線性Radon變換域系數(shù)如圖3a所示，拋物Radon變換域系數(shù)如圖3d所示；可以看出，在使用線性Radon變換時(shí)，僅有線性同相軸的Radon域系數(shù)是聚焦、稀疏的，其他拋物型同相軸Radon域系數(shù)并不稀疏；同樣地，使用拋物Radon變換時(shí)，稀疏性?xún)H對(duì)拋物型同相軸有效；這是因?yàn)楫?dāng)?shù)卣饠?shù)據(jù)同相軸形態(tài)和積分路徑形態(tài)不一致時(shí)，需要一系列的變換參數(shù)和積分路徑組成的基函數(shù)來(lái)對(duì)地震數(shù)據(jù)擬合，從而在Radon域剖面上一系列連續(xù)的p值或q值位置都有系數(shù)分布，導(dǎo)致系數(shù)不是稀疏的，且不能從Radon域剖面得到所有同相軸的速度信息.圖3b、e分別是從Radon域系數(shù)反變換回時(shí)空域的數(shù)據(jù)；圖3c、f分別為經(jīng)過(guò)不同變換后的數(shù)據(jù)與原始數(shù)據(jù)的差剖面，可以看出拋物Radon變換對(duì)線性同相軸有一定的正反變換誤差，特別是在同相軸的起始位置；用公式(3)來(lái)計(jì)算這兩種變換的相對(duì)均方誤差：

圖2 線性和拋物型同相軸同時(shí)存在的地震數(shù)據(jù)Fig.2 Seismic data with linear and parabolic events

圖3 圖2中地震數(shù)據(jù)Radon變換效果(a)線性Radon正變換系數(shù)；(b)線性Radon反變換數(shù)據(jù)；(c)線性Radon變換誤差剖面；(d)拋物Radon正變換系數(shù)；(e)拋物Radon反變換數(shù)據(jù)；(f)拋物Radon變換誤差剖面.Fig.3 Radon transform of seismic data in Fig.2(a)Coefficients of linear Radon forward transform；(b)Data of inverse linear Radon transform；(c)Error of linear Radon transform；(d)Coefficients of parabolic Radon inverse transform；(e)Data of inverse parabolic Radon transform；(f)Error of parabolic Radon inverse transform.

(3)

本文提出的多路徑Radon變換能解決該問(wèn)題，其正反變換公式為

(4a)

(4b)

式中的積分路徑既包含線性項(xiàng)又包含拋物型項(xiàng)，線性Radon變換參數(shù)p和拋物Radon變換參數(shù)q同時(shí)存在；因此，對(duì)于含有線型的直達(dá)波、折射、面波及雙曲型反射波(可近似為拋物型)的地震數(shù)據(jù)，多路徑Radon變換可以看作線性Radon基函數(shù)和拋物型Radon基函數(shù)的組合.選擇p=0基函數(shù)，公式(4)可以對(duì)拋物型同相軸進(jìn)行稀疏表達(dá)；選擇q=0的基函數(shù)，公式(4)可以對(duì)線性同相軸進(jìn)行稀疏表達(dá)；而p和q都不為0的基函數(shù)，表示公式(4)的積分路徑為多項(xiàng)式，可以對(duì)時(shí)距曲線不是標(biāo)準(zhǔn)的線型或拋物型同相軸進(jìn)行擬合.

和傳統(tǒng)Radon變換一樣，多路徑Radon變換也將在頻率域通過(guò)反演方法來(lái)求解Radon域系數(shù)m，對(duì)公式(4a)、(4b)作Fourier正變換，即得到頻率域多路徑Radon變換公式為

(5a)

(5b)

根據(jù)指數(shù)函數(shù)的性質(zhì)，公式(5)可以進(jìn)一步分解為

(6a)

(6b)

(7)

(8)

(9)

(10)

可以通過(guò)加權(quán)迭代方法(Sacchi,1997)進(jìn)行近似求解.

(11)

(12)

(13)

2 模型數(shù)據(jù)測(cè)試

2.1 多路徑Radon變換的稀疏表示能力分析

用多路徑Radon變換對(duì)圖2中的數(shù)據(jù)作正變換，線性積分路徑參數(shù)p的最小值p0=0 s·m-1，采樣間隔dp=2×10-5s·m-1，采樣個(gè)數(shù)np=12；拋物型積分路徑參數(shù)q的最小值q0=0 s2·m-2，采樣間隔dq=2×10-8s2·m-2，采樣個(gè)數(shù)nq=21，3次迭代求解多路徑Radon域系數(shù)；系數(shù)分布如圖4所示，圖中有12列，分別表示積分路徑參數(shù)取值為(p0,[q0,q1,…,q20])，(p1,[q0,q1,…,q20])，…,(p11,[q0,q1,…,q20])時(shí)的多路徑Radon域系數(shù)，縱坐標(biāo)表示時(shí)間，橫坐標(biāo)表示拋物型積分路徑參數(shù).從圖4中可以看出，與傳統(tǒng)Radon變換不同(圖3a、d)，多路徑Radon域系數(shù)可以同時(shí)聚焦，且獨(dú)立分布在兩個(gè)不同區(qū)域上，即圖4中紅色方框標(biāo)識(shí)出的第1列及第11列區(qū)域.第一個(gè)區(qū)域是線性積分路徑參數(shù)p=0×dp區(qū)，此時(shí)p=0，表示積分路徑為拋物型，為q剖面，該區(qū)域的Radon系數(shù)有4個(gè)，分別對(duì)應(yīng)了4個(gè)拋物型積分參數(shù)，能夠與圖2中時(shí)空域四個(gè)拋物型同相軸的起始時(shí)間、速度一一對(duì)應(yīng)；第二個(gè)區(qū)域是線性積分路徑參數(shù)p=10×dp區(qū)，該區(qū)域的Radon系數(shù)有1個(gè)，對(duì)應(yīng)的拋物型積分路徑參數(shù)q=0×dq，表示積分路徑為線性，為p剖面，該位置的Radon系數(shù)與圖2中時(shí)空域線性同相軸起始時(shí)間、速度一一對(duì)應(yīng)，且速度可由線性積分路徑參數(shù)p=10×dp計(jì)算得到；因此，通過(guò)多路徑Radon不僅可以使時(shí)空域不同形態(tài)的同相軸在Radon域系數(shù)同時(shí)達(dá)到稀疏狀態(tài)，而且能夠很容易地識(shí)別這些系數(shù)對(duì)應(yīng)的時(shí)空域同相軸形態(tài)、時(shí)間截距以及速度.這些優(yōu)勢(shì)有助于提高利用Radon變換方法進(jìn)行隨機(jī)噪聲壓制、面波壓制以及波型分離等技術(shù)的處理效果.將多路徑Radon域系數(shù)通過(guò)反變換，我們即可得到時(shí)空域地震數(shù)據(jù)，如圖5a所示，變換誤差剖面如圖5b所示，利用公式(3)計(jì)算得到相對(duì)均方誤差為1.36×10-4，誤差極小，這是該變換能夠用于地震數(shù)據(jù)處理的基礎(chǔ).

圖4 多路徑Radon正變換系數(shù)Fig.4 Coefficients in multi-path Radon forward transform

圖5 多路徑Radon反變換結(jié)果(a)反變換結(jié)果；(b)誤差剖面.Fig.5 Result of multi-path Radon inverse transform(a)Result of inverse transform；(b)Error.

圖6 三種變換的稀疏表示能力Fig.6 Sparse representation capability of three Radon transforms

2.2 多路徑Radon變換在模型數(shù)據(jù)隨機(jī)噪聲壓制及波型分離中的應(yīng)用

多路徑Radon變換的強(qiáng)稀疏表示能力使其在隨機(jī)噪聲壓制和波型分離問(wèn)題中具有天然優(yōu)勢(shì)，先以模型數(shù)據(jù)為例進(jìn)行展示.圖7a為干凈的原始數(shù)據(jù)，包含待分離的線性和拋物型兩個(gè)同相軸，同相軸在遠(yuǎn)偏移距形態(tài)接近且有重疊.地震數(shù)據(jù)共有64道，道間距為20 m，每道有300個(gè)時(shí)間采樣點(diǎn)，采樣間隔為4 ms.將采取常用的傅里葉變換(FK)、傳統(tǒng)的線性Radon、拋物Radon以及本文提出的多路徑Radon變換方法對(duì)其進(jìn)行分離.圖7b、c分別為該數(shù)據(jù)在線性Radon變換域和拋物Radon變換域的系數(shù)分布圖，可以看出在線性Radon變換中僅有線性同相軸的系數(shù)分布是稀疏、聚焦的，如圖7b中白色箭頭所示，而拋物性同相軸的系數(shù)分布在空間上具有延展性，使二者的系數(shù)有部分重疊，很難將其完全分離開(kāi)來(lái)；類(lèi)似的情況也出現(xiàn)在拋物Radon變換系數(shù)中，如圖7c所示；圖7d、e、f分別為線性同相軸的FK譜、拋物同相軸的FK譜以及二者混合一起后的FK譜，很明顯這兩個(gè)軸的FK譜能量差異極大，重疊在一起后的FK譜中這兩個(gè)軸的能量也有部分重疊，同樣也不能將其完全分離開(kāi)；圖7g為多路徑Radon域系數(shù)，可以看出線性同相軸和拋物同相軸的系數(shù)分布在不同區(qū)域且同時(shí)稀疏、聚焦，很容易將其分離.

圖7 時(shí)空域數(shù)據(jù)在不同變換域系數(shù)(a)原始數(shù)據(jù)；(b)線性Radon變換域系數(shù)；(c)拋物Radon變換域系數(shù)；(d)線性軸FK譜；(e)拋物軸FK譜；(f)原始數(shù)據(jù)FK譜；(g)多路徑Radon變換域系數(shù).Fig.7 Coefficients of temporal-spatial data in different transform domains(a)Original data；(b)Coefficients in linear Radon transform domain；(c)Coefficients in parabolic Radon transform domain；(d)FK spectrum of linear event；(e)FK spectrum of parabolic event；(f)FK spectrum of original data；(g)Coefficients of multi-path Radon transform domain.

表1 四種變換波型分離效果比較Table 1 Comparison of data separation by 4 transforms

為了測(cè)試本文方法抗噪能力，在圖7a原始干凈數(shù)據(jù)上添加了一定的隨機(jī)噪聲，數(shù)據(jù)的信噪比為-16.09 dB，如圖9a所示.根據(jù)上文三種變換的稀疏表示能力曲線來(lái)確定噪聲壓制的閾值，首先將圖7變換域中前30%的強(qiáng)能量系數(shù)保留，其他系數(shù)進(jìn)行壓制，再進(jìn)行下一步的分離.我們比較了本文提出的多路徑Radon變換與傳統(tǒng)線性Radon、拋物Radon變換方法的去噪效果.圖9b、c分別為數(shù)據(jù)在線性Radon、拋物Radon變換域系數(shù)，可以看出，由于隨機(jī)噪聲的存在，使得線性同相軸和拋物同相軸的變換域系數(shù)更加難以區(qū)分，進(jìn)而影響下一步的分離；圖9d為數(shù)據(jù)在多路徑Radon變換域系數(shù)，可以看出，即使存在能量較強(qiáng)的隨機(jī)噪聲，由于不同類(lèi)型同相軸在該變換域系數(shù)都同時(shí)稀疏、聚焦，從而仍然能夠較容易將其區(qū)分并分離；線性Radon變換方法的噪聲壓制及波型分離結(jié)果如圖8e、f所示，拋物Radon變換方法的噪聲壓制及波型分離結(jié)果如圖9g、h所示,本文提出的多路徑Radon變換方法的噪聲壓制及波型分離結(jié)果如圖9i、j所示，從圖中可以很明顯看出，線性Radon變換僅對(duì)線性同相軸的噪聲壓制效果較好，拋物Radon變換僅對(duì)拋物同相軸的噪聲壓制有效，而多路徑Radon變換方法對(duì)線性同相軸和拋物同相軸的分離及噪聲壓制同時(shí)有效，同相軸分離完全、無(wú)振幅損失且噪聲壓制干凈.三種方法噪聲壓制及分離后數(shù)據(jù)的信噪比如表2所示，本文提出的多路徑Radon變換方法在去噪及分離后，信噪比達(dá)到28 dB左右，與原始-16.09 dB含噪數(shù)據(jù)相比，信噪比提高了44 dB，其效果明顯優(yōu)于線性Radon和拋物Radon變換.

圖8 波型分離效果(a)(b)傅里葉變換波型分離后的數(shù)據(jù)；(c)(d)線性Radon變換波型分離后的數(shù)據(jù)；(e)(f)拋物Radon變換波型分離后的數(shù)據(jù)；(g)(h)多路徑Radon變換波型分離后的數(shù)據(jù).Fig.8 Waveform separations by different transforms(a)and (b)Fourier transform；(c)and (d)Linear Radon transform；(e)and (f)Parabolic Radon transform；(g)and (h)Multi-path Radon transform.

圖9 去噪及分離效果(a)含噪數(shù)據(jù);(b)線性Radon變換系數(shù);(c)拋物Radon變換系數(shù);(d)多路徑Radon變換系數(shù);(e)(f)線性Radon變換去噪及分離結(jié)果;(g)(h)拋物Radon變換去噪及分離結(jié)果;(i)(j)多路徑Radon變換去噪及分離結(jié)果.Fig.9 Denoised and separated data by multi-path Radon transform(a)Noised data；(b)Coefficients of linear Radon transform；(c)Coefficients of parabolic Radon transform；(d)Coefficients of multi-path Radon transform；(e)and (f)Denoised and separated data by linear Radon transform；(g)and (h)Denoised and separated data by parabolic Radon transform；(i)and (j)Denoised and separated data by multi-path Radon transform.

表2 含噪數(shù)據(jù)噪聲壓制及波型分離效果比較Table 2 Comparison of denoised and separated data by three Radon transforms

2.3 多路徑Radon變換在實(shí)際數(shù)據(jù)噪聲壓制及波型分離中的應(yīng)用

現(xiàn)將該方法應(yīng)用于塔里木某區(qū)域采集的實(shí)際數(shù)據(jù)面波壓制測(cè)試中，如圖10a所示，共71道，道間距為30 m,每道625個(gè)時(shí)間采樣,采樣間隔為4 ms，強(qiáng)烈的面波表現(xiàn)為線性，將近偏移距的有效信號(hào)湮沒(méi).根據(jù)上文模型數(shù)據(jù)的測(cè)試結(jié)果，可以發(fā)現(xiàn)在線性和拋物型同相軸同時(shí)存在的地震數(shù)據(jù)中，線性Radon變換的處理效果優(yōu)于拋物Radon變換，因此在實(shí)際數(shù)據(jù)面波壓制測(cè)試中，將僅提供線性Radon和本文提出的多路徑Radon變換的處理效果對(duì)比.二者面波壓制效果分別如圖10b、c所示，從圖中可以看出這兩種方法都能夠?qū)⒕€性面波壓制，但多路徑Radon變換方法能夠得到分辨率更高的近偏移距反射波同相軸，如圖中A白色橢圓區(qū)域所示；同時(shí)，對(duì)于深層較弱信號(hào)，如2.4 s深度B白色橢圓區(qū)域所示，本文方法恢復(fù)的同相軸更為連續(xù)；從去除的面波數(shù)據(jù)圖10e、f中紅色箭頭處，可以進(jìn)一步發(fā)現(xiàn)線性Radon變換方法去除了更多有效信號(hào).

圖10 面波壓制效果(a)原始數(shù)據(jù)；(b)線性Radon變換面波壓制結(jié)果；(c)多路徑Radon變換面波壓制結(jié)果；(d)留白區(qū)域；(e)(f)分別為(b)(c)去除的面波.Fig.10 Suppression of surface waves by different transforms(a)Original data；(b)Linear Radon transform；(c)Multi-path Radon transform；(d)Whitespace；(e)and (f)Surface waves removed by (b)and (c),respectively.

接著將多路徑Radon變換方法應(yīng)用于VSP數(shù)據(jù)的上下行波分離中，原始VSP數(shù)據(jù)如圖11所示，該數(shù)據(jù)共201道，道間距為20 m，每道1905個(gè)時(shí)間采樣，采樣間隔為4 ms，該數(shù)據(jù)較為復(fù)雜，每個(gè)同相軸都不是標(biāo)準(zhǔn)的線性或雙曲型.圖12為三種Radon變換方法的上下行波分離效果：圖12a、b分別為線性Radon變換方法分離的下行波和上行波；圖12c、d分別為拋物Radon變換方法分離的下行波和上行波；圖12e、f分別為多路徑Radon變換方法分離的下行波和上行波.可以看出，多路徑Radon變換方法的分離效果最佳，能夠?qū)崿F(xiàn)上下行波的有效分離；線性Radon變換方法次之，其振幅存在一定的損失，且在上下行波的起始交叉點(diǎn)位置有部分假頻，如圖12a、b中紅色箭頭位置所示；拋物Radon變換方法的效果最差，在近偏移距處傾角較大的同相軸沒(méi)有得到有效恢復(fù)，且存在拋物型腳印，如圖12c、d中紅色虛框處所示；圖13為近偏移距0～1000 m數(shù)據(jù)分離效果的放大圖，可以更明顯地看出多路徑Radon變換方法上下行波分離的更為徹底且無(wú)信號(hào)損失，該實(shí)例進(jìn)一步顯示了本文方法的有效性.

圖11 原始VSP數(shù)據(jù)Fig.11 Original VSP data

圖13 VSP數(shù)據(jù)上下行波分離近偏移距放大效果(a)(b)線性Radon變換近偏移距分離效果；(c)(d)拋物Radon變換近偏移距分離效果；(e)(f)多路徑 Radon變換近偏移距分離效果.Fig.13 Enlargement of separation of near-offset upward and downward waves in VSP data (a)and (b)Linear Radon transform；(c)and (d)Parabolic Radon transform；(e)and (f)Multi-path Radon transform.

3 結(jié)論

本文提出的基于最小二乘稀疏反演多路徑Radon變換方法，能夠使不同形態(tài)同相軸的Radon域系數(shù)自動(dòng)分離，且保持系數(shù)的稀疏性.該方法解決了在地震數(shù)據(jù)中同時(shí)存在線性和拋物型同相軸時(shí)，傳統(tǒng)Radon變換域系數(shù)不能同時(shí)準(zhǔn)確地進(jìn)行稀疏表達(dá)的問(wèn)題.本文從變換的稀疏程度、計(jì)算效率以及處理效果方面做了詳細(xì)的分析，與傳統(tǒng)線性Radon、拋物Radon相比，在不增加額外計(jì)算量的情況下，多路徑Radon變換能夠更好地識(shí)別不同系數(shù)所對(duì)應(yīng)的同相軸形態(tài)、時(shí)間截距以及速度.這些優(yōu)勢(shì)使Radon變換在隨機(jī)噪聲壓制、面波壓制以及波型分離中獲得更佳的處理效果.

需要注意的是，在地質(zhì)結(jié)構(gòu)特別復(fù)雜區(qū)域，地震剖面中的同相軸不滿(mǎn)足線性及雙曲型走時(shí)特征時(shí)，需要采樣間隔較小的p參數(shù)和q參數(shù)對(duì)其進(jìn)行擬合，即p和q的采樣個(gè)數(shù)增加，使多路徑Radon變換算子規(guī)模成倍加大，從而導(dǎo)致該變換的計(jì)算效率降低，且變換域系數(shù)不再集中于某些點(diǎn)而是在一定范圍內(nèi)進(jìn)行展布，失去稀疏特征.因此，綜合考慮處理效果和計(jì)算效率，多路徑Radon變換在同相軸有明顯線性和雙曲型形態(tài)的地震數(shù)據(jù)中更具優(yōu)勢(shì).

致謝感謝評(píng)審專(zhuān)家提出的中肯意見(jiàn)和建議，使論文得以完善和豐富.