基于三維拋物線Radon變換的多次波壓制算法對(duì)比分析

馬繼濤

(中國石油大學(xué)(北京)地球物理學(xué)院物探系,北京102249)

Radon變換是地震數(shù)據(jù)處理中常用的數(shù)學(xué)算法之一,廣泛用于地震數(shù)據(jù)速度分析、多次波壓制、波場分離和插值等領(lǐng)域。該變換可分為線性、拋物線和雙曲線變換3類,其中線性Radon變換也稱τ-p變換。Radon變換由CLAERBOUT[1]為首的斯坦福地球物理團(tuán)隊(duì)引入到地球物理領(lǐng)域,并初步使用雙曲線Radon變換進(jìn)行速度分析和反演;當(dāng)時(shí)的算法是在時(shí)間域進(jìn)行的,需要求取大型矩陣的逆,計(jì)算量大;而HAMPSON[2]和BEYLKIN等[3]重新定義了拋物線Radon變換,給出了最小二乘解,并成功應(yīng)用于地震數(shù)據(jù)多次波壓制中。拋物線Radon變換在應(yīng)用時(shí),假定動(dòng)校正后的共中心點(diǎn)道集中,一次波是被校正平,由于校正所用的速度偏大,多次波同相軸校正不足,存在一定剩余時(shí)差,即多次波同相軸是向時(shí)間的正方向彎曲。數(shù)學(xué)推導(dǎo)表明,動(dòng)校正后多次波時(shí)距曲線的形狀可以近似用數(shù)學(xué)上的拋物曲線公式表示。因此,可以通過拋物線Radon變換,將CMP道集中拉平的一次波和彎曲的多次波變換到Radon域。Radon域中的橫坐標(biāo)為剩余時(shí)差大小,一次波和多次波剩余時(shí)差存在的差異,使得二者在Radon域能夠更加有效地分離。為更好地保持有效波的能量,多次波壓制的一般做法是在Radon域切除一次波,并反變換回時(shí)空域得到多次波模型,最后將該多次波模型從原數(shù)據(jù)中減去。

與地震數(shù)據(jù)處理中常用的其它變換,如傅里葉變換和小波變換等不同,Radon變換算子并不是正交的,是不可逆的。對(duì)同一個(gè)數(shù)據(jù)做正、反變換后,反變換后的數(shù)據(jù)會(huì)有能量的損失。最常用的解決手段是反演,即在一定限制條件下,將原始數(shù)據(jù)和變換回來的數(shù)據(jù)之間的差異最小化。最常用的一種最小化策略就是HAMPSON[2]和BEYLKIN等[3]給出的最小二乘范數(shù)解,即讓原始數(shù)據(jù)和變換回來的數(shù)據(jù)在最小二乘意義下差最小。這種算法在線性和拋物線Radon變換中最為常用。

除了非正交問題外,Radon變換還受采集空間和時(shí)間有限、離散采樣等因素的影響,存在假頻和分辨率低等問題。變換域能量無法有效聚焦,使得多次波和一次波的分離困難,導(dǎo)致在壓制結(jié)果中存在多次波的殘余,或者損傷一次波。因此,提高Radon變換的分辨率一直是該方法研究的熱點(diǎn)。提高Radon域分辨率的算法,大致可以分為3類。第1類是在時(shí)間域提高分辨率,該類算法可以很好地衰減變換域的假象,得到清晰的高分辨率結(jié)果,但由于需要引入時(shí)間域的大型稀疏算子,計(jì)算效率低,因此很少采用,在此不作詳細(xì)介紹。第2類是在頻率域迭代提高分辨率的算法,該算法將Radon變換視為一個(gè)稀疏反演問題,將上一次計(jì)算得到的模型結(jié)果,作為下一次迭代計(jì)算的加權(quán)算子,以達(dá)到提高變換域分辨率的目的。但頻率域提高分辨率的算法,對(duì)同一個(gè)地震道的所有時(shí)間施加了相同的約束,因此加權(quán)項(xiàng)對(duì)所有時(shí)間是耦合在一起的,會(huì)使得能量較大的同相軸出現(xiàn)假象[4]。頻率域的算法有:SACCHI等[5]提出的基于柯西范數(shù)概率密度函數(shù)Bayes理論的高分辨率算法、SACCHI等[6]提出的基于迭代共軛梯度逐次更新阻尼參數(shù)的高分辨率算法、HERRMANN等[7]提出的基于低頻約束的去假頻高分辨率Radon變換算法和CHEN等[8]提出的能夠保護(hù)地震數(shù)據(jù)弱信號(hào)、去除假頻的基于地震數(shù)據(jù)特定頻率解約束的非迭代高分辨率算法等。第3類提高Radon變換分辨率的方法是頻率時(shí)間域的混合算法,該算法一般是對(duì)最小二乘算法得到的時(shí)間域模型進(jìn)行收斂處理,而其中的正向和反向Radon變換均在頻域進(jìn)行。該算法綜合了頻率和時(shí)間域兩種算法的優(yōu)勢(shì),既提高了分辨率,又保持了數(shù)據(jù)的波形,同時(shí)還沒有明顯增加計(jì)算時(shí)間。如LU[9]提出的基于加速稀疏時(shí)不變的頻率-時(shí)間混合域迭代收縮高分辨率Radon變換算法,在沒有明顯增加計(jì)算時(shí)間的同時(shí),有效提高了分辨率;馬繼濤等[10]對(duì)二維情況下不同迭代閾值收縮高分辨率Radon變換算法進(jìn)行了對(duì)比分析;薛亞茹等[11]對(duì)加權(quán)迭代軟閾值算法的高分辨率Radon變換算法進(jìn)行了探討;羅騰騰等[12]對(duì)混合域高分辨率Radon變換在繞射波分離和成像中的應(yīng)用進(jìn)行了探討。

三維地震數(shù)據(jù),尤其寬方位數(shù)據(jù),形成的三維小面元道集,受到地層各向異性等因素的影響,同相軸存在抖動(dòng)現(xiàn)象,二維Radon變換無法解決該道集的多次波壓制問題。需要考慮三維地震波場的傳播特征,利用三維Radon理論對(duì)不同方位的數(shù)據(jù)進(jìn)行精確表征,能更好地對(duì)其進(jìn)行處理。三維Radon變換考慮了波場三維傳播的特點(diǎn),利用縱測(cè)線和橫測(cè)線兩個(gè)方向的曲率參數(shù),在Radon域?qū)?dòng)校正后的三維小面元道集進(jìn)行表征,可以實(shí)現(xiàn)三維數(shù)據(jù)多次波的有效壓制。很多學(xué)者針對(duì)三維地震波場傳播的特點(diǎn),對(duì)三維Radon算法展開過較為深入地研究。如DONATI等[13]利用三維τ-p變換對(duì)地震數(shù)據(jù)進(jìn)行插值重建處理;HUGONNET等[14-15]給出了三維拋物線Radon變換算法并應(yīng)用在寬方位地震數(shù)據(jù)處理中;如ZHANG等[16]給出了一種時(shí)間和頻率混合域加速三維稀疏時(shí)不變?chǔ)?p變換算法,即基于迭代閾值收縮的高分辨率τ-p變換算法,并應(yīng)用在三維疊前道集的插值重建中;CAO等[17]給出了一種基于匹配追蹤的高分辨率三維τ-p變換算法,有效解決了空間假頻問題,提高了分辨率,取得了較好的插值效果;SUN等[18]給出了一種用于線性噪聲壓制的三維圓錐線性Radon變換算法,用于壓制三維道集中的面波。針對(duì)三維Radon變換過程中的保幅問題,唐歡歡等[19-20]給出了一種高階高分辨率拋物線保幅Radon變換算法并應(yīng)用于三維地震數(shù)據(jù)重建;薛亞茹等[21]給出了一種高階3D變換地震數(shù)據(jù)重建算法,有效解決了地震數(shù)據(jù)重建過程中的保幅問題;MA等[22]給出了一種高階高分辨率三維Radon變換算法,考慮了三維Radon變換的保幅性,并基于低頻約束的方法提高地震數(shù)據(jù)的分辨率。

針對(duì)常規(guī)三維Radon變換算法分辨率低的問題,很多學(xué)者提出了提高分辨率的算法,包括時(shí)間域、頻率域和混合域的算法,但沒有對(duì)這些算法進(jìn)行系統(tǒng)的對(duì)比分析。我們基于Radon變換多次波壓制,對(duì)常用的頻率域和時(shí)間頻率混合域提高分辨率的算法進(jìn)行了系統(tǒng)的對(duì)比分析。首先,詳細(xì)介紹了基于最小二乘的三維Radon變換算法;然后,分別給出了基于迭代重加權(quán)的高分辨率算法,利用低頻解約束高頻運(yùn)算的高分辨率算法和基于迭代閾值收縮的高分辨率算法;最后,利用模擬數(shù)據(jù)對(duì)每一種算法多次波壓制的有效性進(jìn)行驗(yàn)證和對(duì)比分析,并分析各個(gè)算法的抗噪性能。

1 三維拋物線Radon變換原理

1.1 基于最小二乘的三維拋物線Radon變換

三維拋物線Radon變換沿某個(gè)拋物面對(duì)地震數(shù)據(jù)求和,得到三維Radon域數(shù)據(jù)。設(shè)時(shí)空域三維地震數(shù)據(jù)為d(t,x,y),變換后三維Radon域數(shù)據(jù)為m(τ,qx,qy),其中,x和y分別為縱測(cè)線和聯(lián)絡(luò)測(cè)線兩個(gè)方向的偏移距;qx和qy分別為沿著縱測(cè)線和聯(lián)絡(luò)測(cè)線方向的兩個(gè)曲率參數(shù);t和τ分別為時(shí)空域和Radon域的時(shí)間。根據(jù)Radon變換的定義,沿拋物面路徑對(duì)Radon域數(shù)據(jù)進(jìn)行疊加求和,可得到時(shí)空域地震數(shù)據(jù),則有:

d(t,x,y)=?m(τ=t-qxx2-qyy2,qx,
qy)dqxdqy

(1)

公式(1)可離散為:

(2)

式中:nqx和nqy為曲率qx和qy的個(gè)數(shù)。對(duì)公式(2)兩側(cè)關(guān)于時(shí)間做傅里葉變換,將其轉(zhuǎn)變到頻率域,可得:

(3)

式中:D為頻率域的地震數(shù)據(jù);M為頻率域的Radon數(shù)據(jù)。對(duì)某個(gè)特定的xxj和yyk,可寫出D的矩陣計(jì)算式,如下:

(4)

公式(4)可以擴(kuò)充至所有的數(shù)據(jù)D,

即:

(5)

令M兩側(cè)與qxx2相關(guān)的矩陣為Lx,與qyy2相關(guān)的矩陣為Ly,則式(5)可以寫為如下的算子形式:

D=LxMLy

(6)

式中:Lx為與qx、偏移距x有關(guān)的算子,大小為(nx,nqx);Ly為與qy、偏移距y有關(guān)的算子,大小為(nqy,ny)。在nx≠nqx,ny≠nqy情況下,Lx與Ly均不是方陣,可采用與二維Radon變換類似的方式求取Lx和Ly矩陣的最小二乘逆。在實(shí)際數(shù)據(jù)處理時(shí),為保證處理精度,總是使變換所用曲線的數(shù)目大于地震道的道數(shù),即nqx>nx,nqy>ny,這就使得矩陣Lx為欠定矩陣,Ly為超定矩陣。同時(shí),考慮到公式(6)各矩陣的大小,對(duì)Lx和Ly采用不同的方式求解最小二乘逆,如下式:

(7)

式中:I是單位對(duì)角矩陣,μ是為提高矩陣求逆運(yùn)算穩(wěn)定性施加的阻尼因子,取值一般在0.01～1.00之間;H代表矩陣的共軛轉(zhuǎn)置。

1.2 基于迭代重加權(quán)的三維高分辨率Radon變換

公式(7)中,阻尼因子μ對(duì)所有頻率都是相同的,使得矩陣求逆更加穩(wěn)定,但同時(shí)也降低了Radon變換結(jié)果的分辨率。為提高分辨率,SACCHI等[6]提出了基于迭代重加權(quán)的高分辨率Radon變換算法,該算法利用上一次迭代計(jì)算的模型域結(jié)果,在信號(hào)區(qū)域給μ賦一個(gè)較小的值,在非信號(hào)區(qū)賦一個(gè)較大的值,從而將Radon域能量強(qiáng)的點(diǎn)加強(qiáng),能量弱的點(diǎn)減弱,使Radon域能量聚焦,達(dá)到提高Radon域數(shù)據(jù)分辨率的效果,該變換的表達(dá)式如下。

(8)

式中:Cmx和Cmy是和上一次迭代運(yùn)算結(jié)果相關(guān)的加權(quán)矩陣。每次迭代都會(huì)生成新的加權(quán)矩陣Cmx和Cmy,以得到高分辨率的三維Radon域結(jié)果。以Cmx為例,加權(quán)矩陣可通過式(9)計(jì)算得到:

(9)

式中:ε為使得求倒數(shù)穩(wěn)定所加的穩(wěn)定因子,diag(·)是將向量變?yōu)閷?duì)角矩陣的算子。如果變換域M的值較小,則計(jì)算出的約束項(xiàng)使得下一次迭代求解出的M的值更小,亦即算法會(huì)使得變換中的假象噪聲變小;反之算法會(huì)使得變換中的信號(hào)增強(qiáng),達(dá)到提高分辨率的目的。

1.3 基于低頻約束的三維高分辨率Radon變換

基于迭代重加權(quán)的高分辨率算法在處理稀疏采樣數(shù)據(jù)時(shí),由于間距大,空間采樣會(huì)產(chǎn)生空間假頻,導(dǎo)致變換域出現(xiàn)假象,降低分辨率。為解決這一問題,HERRMANN[7]提出了一種不需要迭代的高分辨率Radon變換算法,利用數(shù)據(jù)低頻部分的解對(duì)數(shù)據(jù)高頻部分的運(yùn)算進(jìn)行約束,即用多次的重加權(quán)約束Radon域數(shù)據(jù),使其變得稀疏,每一個(gè)加權(quán)矩陣都是由上一個(gè)頻率的計(jì)算結(jié)果生成。地震數(shù)據(jù)的低頻部分波長長,變換時(shí)不會(huì)出現(xiàn)空間假頻,因此低頻約束的方法可以阻止假頻的出現(xiàn);此外,方法無迭代處理過程,可以提高計(jì)算效率;方法計(jì)算過程中,除阻尼參數(shù)外,沒有其他可調(diào)整的參數(shù),運(yùn)算相對(duì)簡單,可一定程度上減輕處理人員對(duì)方法的測(cè)試過程。

基于低頻約束的三維高分辨率Radon變換可表示為:

(10)

式中:Mn和Dn分別為第n個(gè)頻率的三維Radon域和地震數(shù)據(jù);Wn和Vn分別是縱測(cè)線和聯(lián)絡(luò)測(cè)線方向算子求逆所對(duì)應(yīng)的對(duì)角加權(quán)矩陣。該加權(quán)矩陣是由上一次(n-1)頻率的結(jié)果計(jì)算得到。以Lx算子的加權(quán)矩陣Wn為例,可由下式計(jì)算得到:

Wii(ωn)=‖Mi(ωn-1)‖

(11)

式中:Mi(ωn-1)為在上一次計(jì)算頻率運(yùn)算得到的結(jié)果,此處W對(duì)角加權(quán)矩陣i的范圍為1～nqx;若為聯(lián)絡(luò)測(cè)線方向的對(duì)角加權(quán)矩陣V,則i的范圍為1～nqy。在沒有多次波模型先驗(yàn)信息的情況下,該方法利用數(shù)據(jù)的低頻計(jì)算結(jié)果對(duì)高頻計(jì)算進(jìn)行約束,低頻數(shù)據(jù)在運(yùn)算中不產(chǎn)生假頻,在采樣稀疏的情況下仍能取得很好的高分辨率效果。

1.4 基于迭代閾值收縮的三維高分辨率Radon變換

迭代閾值收縮算法是在時(shí)間-頻率混合域利用反變換后數(shù)據(jù)和原數(shù)據(jù)的殘差更新模型,進(jìn)行提高分辨率的處理。每次迭代在時(shí)間域進(jìn)行閾值收縮處理,增強(qiáng)變換域能量團(tuán)中心的能量,減弱邊緣的發(fā)散能量,可以使得模型域能量團(tuán)更加集中,提高分辨率。

與二維算法類似,三維算法是提前計(jì)算各頻率Lx和Ly矩陣的最小二乘逆,并存儲(chǔ)起來在計(jì)算中直接調(diào)用,以犧牲小的存儲(chǔ)內(nèi)存為代價(jià)提高計(jì)算效率;利用Radon域模型反變換得到數(shù)據(jù)和原始數(shù)據(jù)的殘差,通過逐次迭代收縮閾值更新Radon域數(shù)據(jù),達(dá)到提高分辨率的效果。我們將LU[9]給出的二維算法擴(kuò)展到三維,設(shè)算子Lx和Ly的最小二乘逆分別為Lx_inv和Ly_inv,即令:

(12)

則可通過公式(13)計(jì)算更新變換域模型:

mk=Tα{mk-1+βF-1[Lx_inv·(F[d]-
LxF[mk-1]Ly)·Ly_inv]}

(13)

式中:Lx_inv和Ly_inv在每次迭代中直接調(diào)用;mk是k次迭代后時(shí)間域的Radon模型數(shù)據(jù);Tα代表對(duì)括號(hào)內(nèi)數(shù)據(jù)的閾值處理;β為步長;d為時(shí)空域的地震數(shù)據(jù);F和F-1是傅里葉正變換和反變換算子。由公式(13)可以看出,迭代閾值收縮算法是將地震數(shù)據(jù)和Radon域模型反變換后得到的地震數(shù)據(jù)之差,以一定比例疊加到Radon域,更新時(shí)空域的模型。之后基于公式(14)的收縮函數(shù),對(duì)Radon域模型數(shù)據(jù)進(jìn)行收縮和提高分辨率的運(yùn)算。

(14)

2 模型數(shù)據(jù)試算對(duì)比分析

利用模擬數(shù)據(jù)對(duì)方法的有效性和精確性進(jìn)行驗(yàn)證,分別展示三維常規(guī)最小二乘Radon變換、迭代重加權(quán)的高分辨率Radon變換、基于低頻約束的高分辨率Radon變換以及基于迭代閾值收縮的高分辨率Radon變換結(jié)果。為方便陳述,以LS(least-squares)代表最小二乘變換算法;以IRHR(iterated reweight high-resolution)代表迭代重加權(quán)高分辨率算法;以LOWHR(lower frequencies constrained high resolution)代表低頻約束高分辨率變換算法;以SRTIS(sparse radon transform iterative shrinkage)代表迭代閾值收縮高分辨率算法,對(duì)變換域和多次波壓制結(jié)果進(jìn)行對(duì)比分析。

首先,模擬生成三維地震數(shù)據(jù)CMP道集。假定縱測(cè)線偏移距范圍為-4800～4800m,每個(gè)道集數(shù)據(jù)有50道;聯(lián)絡(luò)測(cè)線方向偏移距范圍為-2400～2400m,測(cè)線間距為200m。數(shù)據(jù)中有4個(gè)一次波和4個(gè)具有不同剩余時(shí)差的多次波,一次波的剩余時(shí)差為0,多次波的剩余時(shí)差分別為600ms,400ms,200ms和40ms,最下方的多次波與一次波之間的時(shí)差非常小,與一次波所在點(diǎn)橫向只有0.04s,即2.3個(gè)樣點(diǎn)的距離。多次波能量有強(qiáng)有弱,其振幅由上到下依次為-1.0,0.5,0.1,-0.9。三維小面元CMP道集如圖1所示。利用本文算法對(duì)此三維CMP道集進(jìn)行多次波壓制測(cè)試。

圖1 三維小面元CMP道集

為更好地分析對(duì)比多次波壓制的效果,還抽取了原始數(shù)據(jù)和一次波的第1,7,13,19和25個(gè)小面元道集(圖2)?？梢钥闯?不同位置的小道集中,多次波和一次波具有不同的時(shí)差差異特點(diǎn),因而利用相同參數(shù)對(duì)所有小道集進(jìn)行處理不可靠。圖3展示的是與圖1對(duì)應(yīng)的按標(biāo)量偏移距大小排列的三維CMP道集,在該道集中可以看出,一次波和多次波受方位各向異性的影響,同相軸有抖動(dòng)?；诙SRadon變換無法直接對(duì)其進(jìn)行處理,需要考慮基于三維拋物線Radon理論的多次波壓制方法。

圖2 5個(gè)小面元道集a 含多次波數(shù)據(jù)道集;b 一次波道集

圖3 按標(biāo)量偏移距排列的三維CMP道集

2.1 多次波壓制的有效性

分別利用4種算法對(duì)該模擬數(shù)據(jù)進(jìn)行了處理。變換域應(yīng)為三維模型(τ,qx,qy),但由于三維模型并不能直觀顯示各方法的聚焦效果,將變換域中每個(gè)qy面板的數(shù)據(jù)進(jìn)行了累加求和,并基于該結(jié)果對(duì)各方法分辨率進(jìn)行了對(duì)比。

圖4為變換域的結(jié)果,每個(gè)箭頭處都有一次波或多次波聚焦的能量?？梢钥闯?圖4a的LS算法可以將數(shù)據(jù)變換到Radon域,振幅為0.1的弱多次波在變換域也清晰可見,但該算法假象嚴(yán)重;一般而言,能量越強(qiáng),同相軸假象越嚴(yán)重。因此在0.1s和1.0s處的多次波,都有較嚴(yán)重的剪刀狀發(fā)散假象。切除時(shí)會(huì)損傷一次波能量,或?qū)е露啻尾▔褐撇粡氐?。而其?種高分辨率算法都較好地減弱或壓制了剪刀狀發(fā)散假象,變換后的分辨率都得到了很大提升。在弱多次波的識(shí)別方面,LOWHR和SRTIS兩種算法都能看到較為清晰的聚焦點(diǎn),但I(xiàn)RHR算法在多次迭代后,將弱多次波誤認(rèn)為假象噪聲,導(dǎo)致該多次波的能量減弱(圖4b),這會(huì)給變換域多次波的識(shí)別帶來一定干擾。

在圖4中的三維Radon變換結(jié)果中,設(shè)置切除參數(shù)為0.02,即將變換域qx小于0.02的能量切除掉。數(shù)據(jù)下方1.0s處,一次波位于橫向第2個(gè)樣點(diǎn),多次波位于橫向第5個(gè)樣點(diǎn),切除點(diǎn)位于第3個(gè)樣點(diǎn),切除位置位于一次波和多次波中間位置,且與二者相距非常近,因此切除對(duì)方法分辨率的要求非常高。切除一次波能量后,保留多次波并將其變換回時(shí)空域,得到圖5所示的多次波模型,之后將得到的多次波模型從原數(shù)據(jù)中減掉,以最大程度保留一次波的能量,得到圖6和圖7所示的多次波壓制后的小道集和標(biāo)量偏移距道集,其中圖6為與圖2所對(duì)應(yīng)的小道集壓制后的結(jié)果,圖7為與圖3相對(duì)應(yīng)的按標(biāo)量偏移距大小排列的多次波壓制后的結(jié)果。在圖5所示的多次波模型中,LS算法有很嚴(yán)重的一次波能量泄漏問題,而IRHR結(jié)果由于分辨率得到了提高,一次波能量有所減弱;圖5a和圖5b箭頭所指處均為所泄漏的一次波能量。在LOWHR和SRTIS算法結(jié)果中未見一次波能量同相軸,說明這兩種高分辨率算法將一次波和多次波分離得更為徹底。在圖6和圖7的多次波壓制結(jié)果中也可以看出,LS和IRHR算法壓制結(jié)果中有多次波能量的殘余(箭頭所指),同時(shí)兩種算法由于估計(jì)多次波模型中有一次波能量的泄漏,得到的一次波能量也有損傷,即LS和IRHR算法既存在多次波壓制不徹底的問題,也存在一次波能量的損傷問題。SRTIS算法多次波壓制的結(jié)果中(圖6d),箭頭所指處存在多次波能量的殘余,說明該方法雖然在變換域能夠很好的分離開一次波和多次波,但變換回的多次波與數(shù)據(jù)中的多次波振幅不一致,導(dǎo)致壓制結(jié)果有殘余,圖7d箭頭所指也可以看到多次波的殘余;與其它方法相比,LOWHR算法無論在多次波和一次波的分離,還是多次波壓制中均取得了較好的多次波壓制效果。

圖4 幾種三維Radon變換算法的變換結(jié)果a LS;b IRHR;c LOWHR;d SRTIS

圖5 不同變換方法估計(jì)的多次波模型a LS;b IRHR;c LOWHR;d SRTIS

圖6 不同變換方法多次波壓制后的小道集a LS;b IRHR;c LOWHR;d SRTIS

圖7 不同變換方法多次波壓制后的標(biāo)量偏移距道集a LS;b IRHR;c LOWHR;d SRTIS

表1列出了采用幾種方法處理模擬數(shù)據(jù)所用的時(shí)間。從表1中可以看出,各方法運(yùn)行時(shí)間由大到小順序依次為:LOWHR>IRHR>SRTIS>LS。LS算法耗時(shí)最短,但其效果差;LOWHR算法耗時(shí)最長,但其效果最好。相對(duì)而言,SRITS和LOWHR算法均可以取得較好的壓制效果,但SRTIS算法運(yùn)行時(shí)間更短;SRTIS算法為取得較好的結(jié)果,需要調(diào)試多個(gè)參數(shù),而LOWHR算法與LS算法類似,只需調(diào)整阻尼參數(shù),計(jì)算較為簡單直接。

表1 各算法運(yùn)行時(shí)間對(duì)比

2.2 多次波壓制方法的抗噪性

對(duì)模擬數(shù)據(jù)加入一定量的高斯噪聲,測(cè)試了各算法的抗噪性。結(jié)果表明,噪聲對(duì)各算法的變換分辨率及多次波壓制結(jié)果均存在一定的影響,但調(diào)整阻尼參數(shù)后均可取得較好的結(jié)果。本文給出了信噪比為20的情況下,變換域及多次波壓制的效果。圖8為LS、IRHR、LOWHR和SRTIS算法抗噪性的測(cè)試結(jié)果。其中,LS、LOWHR和SRTIS算法均調(diào)大了阻尼參數(shù),IRHR算法保持原阻尼參數(shù)不變?？梢钥闯?各種高分辨率算法在沒有噪聲的情況下,均有很好的多次波壓制效果,但弱能量多次波在SRTIS算法變換域更加清晰。說明SRTIS算法充分利用了頻率域和時(shí)間域兩種算法的優(yōu)勢(shì),既可以得到高分辨率的變換結(jié)果,又保持了數(shù)據(jù)中弱同相軸的能量,有利于數(shù)據(jù)中多次波的識(shí)別和壓制。

圖8 加噪數(shù)據(jù)不同變換算法的結(jié)果(左:變換域;右:壓制后的小道集)a LS;b IRHR;c LOWHR;d SRTIS

3 結(jié)論與建議

系統(tǒng)介紹了三維Radon變換的理論,給出了求解該問題的最小二乘算法,基于迭代重加權(quán)、低頻約束和迭代閾值收縮的高分辨率算法,并利用模擬數(shù)據(jù)驗(yàn)證了各方法的有效性、抗噪性和對(duì)小時(shí)差多次波的壓制能力。通過對(duì)比和分析研究,可以得出如下結(jié)論。

1)受采集空間和時(shí)間有限、離散采樣的影響,三維Radon變換的最小二乘解存在假象,分辨率低。

2)頻率域高分辨率算法在變換域時(shí)間方向存在噪聲干擾,但可提高橫向分辨率;時(shí)間頻率混合域高分辨率算法可以在雙方向提高變換分辨率。

3)在有噪聲干擾的情況下,通過調(diào)整變換參數(shù),各算法均可以取得較好的多次波壓制效果。

4)在小時(shí)差多次波的壓制方面,迭代收縮閾值高分辨率算法在變換域具有較好的收斂性,可以有效識(shí)別多次波,但分辨率的提升會(huì)影響所估計(jì)多次波模型的保幅性;相對(duì)而言,低頻約束高分辨率算法在小時(shí)差多次波壓制方面效果最佳。

隨著反演精度的提高,對(duì)地震數(shù)據(jù)疊前道集的保幅性要求越來越高,如何在多次波壓制過程中保持有效波的振幅,需要進(jìn)一步的研究。另外,與二維變換算法相比,三維變換算法計(jì)算效率較低,盡管迭代閾值收縮算法在計(jì)算過程中以保存矩陣逆的方式減少了運(yùn)算時(shí)間,然而三維算法的計(jì)算量仍然較大,不利于大規(guī)模生產(chǎn)應(yīng)用。因此,需要進(jìn)一步研究如何提高三維算法的計(jì)算效率。