馬繼濤 劉仕友 廖 震
(①中國石油大學(北京)地球物理學院,北京 102249; ②中海石油(中國)有限公司海南分公司,海南海口 570100; ③中國石油大學(華東)地球科學與技術(shù)學院,山東青島 266580)
地震波從地下巖層反射回地表后,會在地表面發(fā)生下行反射,形成多次波。受海水表面強反射系數(shù)的影響,多次波是海洋地震數(shù)據(jù)中常見的干擾,對其壓制是地震數(shù)據(jù)處理中的一項重要內(nèi)容。Radon變換是業(yè)界多次波壓制最常用的有效手段之一,同時也被廣泛地應(yīng)用于地震數(shù)據(jù)重建、波場分離、速度分析等領(lǐng)域。
拋物線Radon變換由Hampson[1]引入二維地震數(shù)據(jù)處理。Radon變換假設(shè)正常時差校正后的共中心點(CMP)道集中一次波是能被校正平的,剩余時差為零。由于對多次波做校正所用的速度是深層一次波速度,動校速度偏大,導致同相軸校正不足,產(chǎn)生一定的剩余時差。此時多次波時距曲線的形狀由雙曲線轉(zhuǎn)變?yōu)閽佄锞€,拋物線Radon變換可將CMP道集中動校平的一次波和彎曲的多次波變換到Radon域,該域中橫坐標為剩余時差。由于一次波與多次波剩余時差存在差異,二者在變換域成為分離開的聚焦點。為保護有效波能量,壓制多次波的常用方法是在Radon域切除一次波,并將多次波變換回時空域,得到多次波模型,之后將其從原始數(shù)據(jù)中消除。
但常規(guī)Radon變換算法受到采集空間和時間有限、離散采樣等因素的影響,存在假頻、分辨率低等問題,致使在變換域無法徹底分離多次波與一次波; 此外,該變換只是將數(shù)據(jù)沿拋物線路徑進行疊加求和,并未考慮同相軸振幅的橫向變化,使分離得到的多次波與原數(shù)據(jù)中多次波的振幅存在一定差異,這些都會影響最終壓制效果。
在提高變換域內(nèi)一次波與多次波的可分離性方面,很多學者做了研究。如Sacchi等先后提出基于柯西范數(shù)概率密度函數(shù)的Bayes理論的高分辨率算法[2]及基于迭代共軛梯度逐次更新阻尼參數(shù)的高分辨率算法[3],Herrmann等[4]提出基于低頻約束的去假頻高分辨率Radon變換算法,Trad等[5]綜合分析了稀疏Radon變換的算法,Chen等[6]提出了能保護地震數(shù)據(jù)弱信號、去除假頻的基于地震數(shù)據(jù)特定頻率解約束的非迭代高分辨率算法。近年來,也有學者提出頻率—時間混合域提高分辨率的迭代收縮算法。Lu[7]提出的基于加速稀疏時不變的頻率—時間混合域迭代收縮高分辨率Radon變換算法,在未顯著增加計算時間的情況下,有效提高了分辨率; 薛亞茹等[8]研究了基于加權(quán)迭代軟閾值的頻率—時間混合域高分辨率Radon變換算法; 羅騰騰等[9]基于頻率—時間混合域的高分辨率Radon變換算法對繞射波進行分離與成像。
在提高Radon變換的保幅性方面,也有很多學者做了相關(guān)研究。Nowak等[10]、薛昭等[11]分別分析了Radon變換及其去噪方法的保幅性; 薛亞茹等[12-13]、Vyas等[14]分別給出了基于不同正交多項式的高階Radon變換算法,增加了描述振幅變化的梯度和曲率信息,拓展了傳統(tǒng)Radon變換算法,使其在具有振幅變化特征的數(shù)據(jù)處理中具有更高準確度; Wang等[15]給出一種能模擬道集中振幅變化的Radon變換算法,也取得了較好的保幅效果。
對于二維數(shù)據(jù)和部分三維窄方位地震數(shù)據(jù),應(yīng)用二維Radon變換方法可取得較好的效果。隨著地震采集技術(shù)的進步,由于寬方位采集具有方位寬、照明充分等優(yōu)勢,越來越受到人們的青睞。與窄方位數(shù)據(jù)相比,寬方位地震數(shù)據(jù)存在隨方位變化的特點,相同或相近炮檢距、但不同方位的地震道,受不同方向地層傾角或地層各向異性的影響,可能存在較大的時差差異; 當把動校正后的地震數(shù)據(jù),按照標量炮檢距從小到大的順序排列,形成一個與二維情形類似的三維共中心點道集(Common-Cell Gathers)時,不同方位的數(shù)據(jù)由于存在一定的時差差異,相鄰道時差變化劇烈,呈現(xiàn)同相軸樣點跳動的特點,常規(guī)二維拋物線Radon算法無法精確表征該數(shù)據(jù),更難以有效分離有效波與干擾波,即壓制效果差。三維共中心點道集可分別按照x和y兩個方向的炮檢距大小,即矢量炮檢距進行排列,形成由不同方位的小道集組成的共中心點道集; 該三維共中心點道集可更直觀地描述三維波場的傳播; 但不同方位道集形態(tài)可能存在較大差異,加之拖纜漂移、不同方位地層特性差異等的影響,不同方位的多次波時差會存在差異。顯然,利用相同參數(shù)的二維變換對小道集進行多次波壓制也無法取得較好效果。
為此,需要考慮三維地震波場的傳播特征,利用三維Radon理論精確表征不同方位的數(shù)據(jù),以更好地對其進行處理。
三維Radon算法考慮了波場三維傳播的特點,利用x和y兩個方向的曲率參數(shù),在Radon域?qū)有U蟮娜S共中心點道集進行表征,可對時空域的三維地震數(shù)據(jù)做較好描述。基于三維地震波場傳播的特點,很多學者對三維Radon算法開展了較深入的研究。Donati等[16]利用三維τ-p變換對地震數(shù)據(jù)進行插值重建; Hugonnet 等[17-18]給出三維拋物線Radon變換方法并應(yīng)用于寬方位地震數(shù)據(jù)處理。但常規(guī)三維變換方法在Radon域?qū)挿轿坏卣饠?shù)據(jù)進行表征仍存在一些與二維方法共性的問題,這緣于寬方位地震數(shù)據(jù)本身特點。如寬方位數(shù)據(jù)道集測線/炮線間距大,普通變換得到的結(jié)果存在嚴重的假頻問題; 可利用將低頻運算的解作為約束的方法,或迭代閾值收縮的算法去假頻,以提高分辨率。如Zhang等[19]提出一種加速三維稀疏時不變τ-p變換方法,即基于迭代收縮的高分辨率τ-p變換方法,并應(yīng)用于三維疊前道集的插值重建;Cao等[20]給出基于匹配追蹤的高分辨率三維τ-p變換方法,有效應(yīng)對了空間假頻,提高了分辨率,取得了較好插值效果。在保幅處理方面,同樣可考慮利用正交多項式擬合表征地震數(shù)據(jù)的振幅。唐歡歡等[21]提出高階高分辨率拋物線保幅Radon變換方法,并應(yīng)用于三維地震數(shù)據(jù)重建; Wang等[22]也給出一種類似的三維保幅Radon變換數(shù)據(jù)重建方法。
本文系統(tǒng)闡述了三維Radon變換算法,并給出該變換的最小二乘解;為提高變換的分辨率,引入基于迭代閾值收縮的三維高分辨率算法;為提高算法的保幅性,引入正交多項式變換,并形成一種三維高精度保幅拋物線Radon變換算法。文中給出一套橫向振幅變化的模擬數(shù)據(jù),并基于最小二乘和高分辨率算法對其進行處理,分別對比了變換域、多次波壓制后的結(jié)果,并給出了壓制結(jié)果的一次波與真實一次波之間的差,表明該方法具有保幅性。最后給出海洋寬方位地震數(shù)據(jù)實例,進一步驗證該方法兼具分辨率高且保幅的優(yōu)勢。
二維拋物線Radon變換沿拋物線路徑對地震數(shù)據(jù)進行求和,得到Radon域數(shù)據(jù)。類似地,三維拋物線Radon變換沿某個拋物面對地震數(shù)據(jù)進行求和,得到三維Radon域數(shù)據(jù)。設(shè)x、y分別為縱測線和非縱測線兩個方向的炮檢距,qx、qy分別為(變換所用的)曲率參數(shù),時空域三維地震數(shù)據(jù)為d(t,x,y),三維Radon域數(shù)據(jù)為m(τ,qx,qy),則三維Radon反變換可表示為
d(t,x,y)=?m(τ=t-qxx2-qyy2,qx,qy)dqxdqy
(1)
其離散形式寫成
(2)
式中Nqx、Nqy分別為曲率qx和qy的個數(shù)。對該式等號兩邊做傅里葉變換,再將其轉(zhuǎn)變到頻率域,可得
(3)
將其寫成算子形式
D=LxMLy
(4)
式中Lx、Ly分別為與炮檢距x、y有關(guān)的算子。
從式(4)可知,可采用與二維Radon變換類似的最小二乘方法,分別求取矩陣Lx、Ly的逆,以得到Radon域數(shù)據(jù)M。因此,在最小二乘意義下有
(5)
式中:I為單位矩陣;λ是提高矩陣求逆運算穩(wěn)定性的阻尼因子,取值范圍一般是0.1~1.0。
三維迭代收縮閾值方法,是提前計算矩陣Lx和Ly的最小二乘逆,并存儲到大型矩陣中,在后續(xù)運算中直接調(diào)用; 利用Radon域模型變換回的時空域數(shù)據(jù)與原始數(shù)據(jù)的殘差,通過逐次迭代收縮閾值更新Radon域數(shù)據(jù),達到提高分辨率的效果。
設(shè)存儲Lx和Ly最小二乘逆的大型矩陣分別為BLxinv、BLyinv,則有
(6)
mk=Tα{mk-1+βF-1
[BLxinv(F[d]-LxF[mk-1]Ly)BLyinv]}
(7)
式中:Tα代表對括號內(nèi)數(shù)據(jù)的一個閾值處理;mk是k次迭代后時間域的Radon模型數(shù)據(jù);F和F-1分別是傅里葉正變換和反變換算子。由上式可看出,迭代閾值收縮算法,將地震數(shù)據(jù)與Radon域模型反變換后得到的地震數(shù)據(jù)之差,以一定比例投影到Radon域模型中,并在每次迭代處理之后,基于如下函數(shù)
(8)
在本算法的迭代過程中,直接調(diào)用了已有的最小二乘逆計算結(jié)果。即只需犧牲很小的存儲內(nèi)存,就可大大提高計算效率。
Radon變換算法只是沿特定的曲線對地震數(shù)據(jù)進行求和疊加,且假設(shè)地震數(shù)據(jù)具有標準的特定形態(tài)的同相軸特征,變換過程中并未考慮地震數(shù)據(jù)振幅的橫向變化; 而實際上由于受地下介質(zhì)復雜性的影響,地震數(shù)據(jù)一般不具有標準同相軸的特性,這樣會導致變換分辨率的降低和變換的不保幅,使變換回的數(shù)據(jù)與原始數(shù)據(jù)存在一定的振幅差異,導致處理結(jié)果中有多次波殘余。地震數(shù)據(jù)的保幅性是與變換的分辨率、變換算子的可逆性緊密聯(lián)系的。分辨率高,會導致保幅性變差; 變換算子的相對可逆,會導致分辨率降低。因此,研究三維Radon變換的保幅算法對于提高變換精度、提升最終處理效果是非常必要的。
由于CDP道集內(nèi)地震數(shù)據(jù)的振幅變化平緩,Johansen等[23]采用了正交多項式對數(shù)據(jù)的AVO效應(yīng)進行擬合; Xue等[13]將Radon變換與正交多項式擬合相結(jié)合,提出了保幅高階Radon變換算法。本文采用類似思路,利用兩個正交多項式系數(shù)p(x)和p(y),分別對x和y方向的地震數(shù)據(jù)振幅變化進行擬合,并將其融合到三維Radon變換算子中。
以x方向的正交多項式p(x)為例。設(shè)正交多項式系數(shù)的階數(shù)為3,則保幅算子可由下式求取
(9)
(10)
(11)
將保幅算子與Radon變換算子Lx、Ly相結(jié)合,即為新的保幅Radon變換算子
(12)
將式(12)代入式(5)或式(7),即可進行保幅的三維最小二乘和高分辨率Radon變換。
本文分別利用模擬數(shù)據(jù)和實際數(shù)據(jù)驗證方法的精確度和有效性。分別展示三維Radon變換的最小二乘算法、基于迭代閾值收縮的三維高分辨率算法及保幅的三維高分辨率算法等的結(jié)果,并對變換域多次波壓制結(jié)果進行了對比分析。
首先生成模擬數(shù)據(jù)CDP道集。假設(shè)橫測線方向共有13個CDP道集,組成三維CDP道集。利用本文三維Radon算法對該三維CDP道集展開多次波壓制測試。生成模擬CDP道集時,假定縱測線炮檢距范圍是0~4800m,道間距為200m,縱測線方向CDP道集共有25道; 橫測線炮檢距范圍是-2400~2400m,道間距為400m,由此模擬得到三維情形下的CDP道集(圖1)。
圖1 三維CDP道集
從所生成的數(shù)據(jù)可見,一次波和多次波在y方向都有變化,同相軸平的是一次波,其剩余時差接近于0; 其余同相軸都有一定的彎曲,是多次波,四個多次波的剩余時差分別為0.043、0.080、0.080、0.020s。其中,最上部多次波振幅很弱,剩余時差都為0.080s的多次波沿y方向的變化趨勢存在一定差異; 最下部多次波的振幅與一次波相反,且緊鄰一次波。該數(shù)據(jù)中,有與一次波同相軸十分接近、剩余時差小的多次波; 也有時差差異大、能量很弱的多次波。數(shù)據(jù)x方向同相軸振幅存在變化。
為更好地對比分析多次波的壓制效果,還抽取了原始數(shù)據(jù)和一次波的第1、7和13個小面元道集(圖2)。為更直觀地顯示三維情形下的CDP道集,在圖3中展示了與圖1對應(yīng)的沿標量炮檢距方向分布的三維CDP道集。在該道集中可看到,一次波和多次波受方位各向異性的影響,同相軸有很大的抖動。基于二維Radon變換直接對其進行處理,無法得到較理想效果; 同時,從圖2a中的小面元道集也可看出,不同的小道集中多次波與一次波存在不同的時差差異特點,利用相同參數(shù)對所有小道集進行處理是不可靠的。因此,需要利用上文提及的三維拋物線Radon變換對模擬數(shù)據(jù)的多次波做壓制處理。
首先給出基于最小二乘理論的三維拋物線Radon變換的結(jié)果(圖4)。該Radon域結(jié)果分別展示了變換后的三維模型和單個qy切片的變換結(jié)果??梢娨淮尾ㄅc剩余時差為0.020s的多次波存在嚴重的相互干涉,剩余時差為0.080s的兩個多次波能量發(fā)散,但在結(jié)果中可看到能量較強。而剩余時差為0.043s的多次波,由于其能量較弱,在變換域剖面中并未呈現(xiàn)明顯的能量團。定義的切除數(shù)值為0.010s,即在Radon域?qū)⑹S鄷r差小于0.010s的數(shù)據(jù)切除,并將其變換回時空域,得到圖5所示的處理結(jié)果。從中可看到能量較弱的多次波(圖5a),在一次波下方也可看到有一些多次波殘余(圖5b)。
圖2 原始數(shù)據(jù)(a)與一次波(b)的小面元道集
圖3 沿標量炮檢距方向排列的三維CDP道集(a)含有多次波的數(shù)據(jù); (b)真正的一次波
為更好地對壓制結(jié)果進行對比分析,同樣選取其中三個小道集進行放大展示(圖6)。從圖6a可見,由于三維最小二乘算法拋物線Radon變換的變換精度不夠、分辨率低,同相軸能量無法在Radon域聚焦,切除時會產(chǎn)生能量泄露現(xiàn)象,導致變換回的多次波中有一次波的信息(箭頭所指),從而使最終壓制結(jié)果中一次波能量有損失; 同時分辨率較低也導致多次波殘余(圖6b箭頭所指)。以上結(jié)果均說明最小二乘方法無法滿足高精度地震數(shù)據(jù)處理的需求。
為改進多次波壓制效果,利用基于迭代閾值收縮的三維高分辨率拋物線Radon變換方法對此數(shù)據(jù)進行提高分辨率處理。圖7為該方法Radon域的結(jié)果,可見該方法能大幅度提高Radon域的分辨率,在一次波與多次波交叉處,原本離散的兩個點變得更聚焦; 振幅較弱、剩余時差為0.043s的多次波在Radon域也可看到; 但由于數(shù)據(jù)橫向振幅的變化,變換域的聚焦點無法進一步集中,導致剩余時差為0的一次波與剩余時差為0.020s的多次波無法徹底分離,在變換域數(shù)據(jù)圖7上部中間0.010s處做切除,會對多次波壓制帶來一定影響。在Radon域?qū)⒁淮尾▽?yīng)的點切除后,變換回x-t域,并將變換回的多次波(圖8a)與原數(shù)據(jù)做相減,得到多次波壓制后的結(jié)果(圖8b)。從中可見多次波壓制的結(jié)果與最小二乘方法相比有明顯改善。但從圖9放大的小道集剖面中可看出,多次波的壓制損傷了部分有效波的能量(圖9a箭頭所指)。
圖4 最小二乘三維拋物線Radon變換Radon域結(jié)果(a)三維模型; (b)單個qy的變換結(jié)果
圖5 最小二乘方法處理結(jié)果的小道集(a)壓制的多次波; (b)壓制多次波后的一次波
圖6 最小二乘方法處理結(jié)果的三個小道集(a)壓制的多次波; (b)壓制多次波后的一次波
圖7 三維高分辨率拋物線Radon變換后Radon域結(jié)果(a)三維模型; (b)單個qy的變換結(jié)果
圖8 高分辨率方法處理結(jié)果小道集(a)壓制的多次波; (b)壓制多次波后的一次波
地震數(shù)據(jù)橫向振幅的變化降低了高分辨率算法的聚焦性,導致有效波能量損傷,且存在多次波的殘余。必須考慮橫向振幅變化的影響,因此采用保幅的高分辨率算法,引入x和y方向的保幅算子對該數(shù)據(jù)進行處理。圖10為該方法變換域第二階對應(yīng)結(jié)果,從中可見一次波與多次波得到了較好的分離,且振幅較弱的多次波也有很好的聚焦效果,這有助于時空域振幅較弱多次波的有效辨別。與圖8和圖9相比,圖11、圖12中壓制的多次波剖面中已經(jīng)沒有一次波能量; 壓制多次波后的一次波剖面(圖11b、圖12b)上基本沒有多次波殘余,壓制更徹底。由于在分辨率得到提升的同時,對同相軸振幅的橫向變化進行了擬合,因此多次波壓制更徹底,一次波能量得到更充分保護。
圖9 高分辨率方法處理結(jié)果(放大顯示)的三個小道集(a)壓制的多次波; (b)壓制多次波后的一次波
圖10 三維保幅高分辨率拋物線Radon變換Radon域結(jié)果(a)三維模型; (b)單個qy的變換結(jié)果
圖11 三維保幅高分辨率方法處理結(jié)果小道集(a)壓制的多次波; (b)壓制多次波后的一次波
為進一步對比三種方法壓制多次波的效果,將多次波壓制后結(jié)果與真實一次波(圖2b)進行相減(圖13)??擅黠@看出,無論是最小二乘方法還是高分辨率方法,對橫向振幅存在變化的數(shù)據(jù)進行處理是不保幅的,高分辨率的算法盡管會提高變換域的分辨率,但這些方法都不可避免地對橫向振幅變化的一次波同相軸造成損傷。由于時空域近炮檢距數(shù)據(jù)變換到三維Radon域會表現(xiàn)出橫向拖尾假象,遠炮檢距數(shù)據(jù)在三維Radon域表現(xiàn)為縱向拖尾假象,而多次波壓制是在橫向某剩余時差位置對一次波與多次波進行的分離,因此多次波壓制主要對近炮檢距一次波數(shù)據(jù)造成了損傷。這與圖13結(jié)果相對應(yīng),壓制結(jié)果與一次波之間的差值主要表現(xiàn)在近炮檢距中。圖14為標量炮檢距排列的道集壓制結(jié)果對比,也可明顯看出三維保幅高分辨率算法的優(yōu)勢。
表1給出了對模擬數(shù)據(jù)進行處理時,最小二乘、高分辨率以及保幅高分辨率三種算法的計算時間??煽闯霰7呔人惴ㄓ嬎懔枯^大,約為高分辨率算法的5倍。在實際地震數(shù)據(jù)處理應(yīng)用中,應(yīng)權(quán)衡計算效率與計算精度,以選取最適宜算法。
為更好地對比各算法的保幅性能,本文還給出了真實一次波和壓制后結(jié)果一次波的振幅變化曲線(圖15)??梢姳7叻直媛仕惴?HOHR)與原始數(shù)據(jù)(True)曲線最接近,也表明本文算法對一次波能量具有良好保護性能。
利用某海洋實際數(shù)據(jù)CDP道集對方法進行了測試,該數(shù)據(jù)為M淺海OBC數(shù)據(jù),經(jīng)過一系列處理后,中遠炮檢距處有一些動校時差較小的多次波,影響了中深層的有效波信息。分別利用最小二乘、高分辨率和保幅高分辨率三種算法對該道集進行多次波壓制處理。由于高分辨率Radon變換方法結(jié)果與最小二乘方法相比變化不大,在此就未給出最小二乘算法的結(jié)果。壓制之前的道集、高分辨率及保幅高分辨率算法壓制的小道集結(jié)果如圖16所示。
圖13 不同方法多次波壓制結(jié)果與真實一次波的差值(a)最小二乘方法; (b)高分辨率方法; (c)保幅高分辨率方法
圖14 不同方法多次波壓制結(jié)果對比(標量炮檢距排列)(a)最小二乘方法; (b)高分辨率方法; (c)保幅高分辨率方法
圖15 各算法保幅性對比
表1 各算法計算時間對比
從壓制處理前、后的小道集來看,高分辨率Radon變換方法對于遠纜的多次波壓制效果較差。圖16b中第7個小道集為近纜數(shù)據(jù),對應(yīng)的多次波基本已被壓制干凈; 第1和第13個小道集的遠纜數(shù)據(jù)中箭頭所指的多次波卻因振幅與原數(shù)據(jù)存在差異,并未被徹底壓制。保幅高分辨率算法(圖16c)對遠纜小道集的多次波壓制則較徹底。
從圖17中的CDP道集對比也可看出,保幅高分辨率算法較高分辨率算法而言,可以更多地壓制數(shù)據(jù)中由方位各向異性引起的抖動的多次波(向上箭頭所指),即小道集中、遠纜的多次波。
圖16 三維Radon變換壓制前、后小面元道集對比(a)原始數(shù)據(jù); (b)高分辨率方法壓制結(jié)果;(c)保幅高分辨率方法結(jié)果
圖17 多次波壓制前、后CDP道集對比(a)原始數(shù)據(jù); (b)高分辨率方法壓制結(jié)果;(c)保幅高分辨率方法結(jié)果
本文對三維Radon變換多次波壓制理論進行了系統(tǒng)論述; 利用基于迭代閾值收縮的方法在時間和頻率域運算,應(yīng)對方法分辨率低導致的一次波與多次波分離不徹底問題; 利用多項式擬合地震數(shù)據(jù)橫向振幅的變化,應(yīng)對方法的不保幅缺陷。針對模擬和實際數(shù)據(jù)進行了方法驗證,結(jié)果表明保幅高分辨率Radon變換在多次波壓制方面能取得較好效果。通過分析和研究,得出以下認識和結(jié)論:
(1)受采集空間和時間、離散采樣等因素影響,三維最小二乘Radon變換方法的分辨率較低;
(2)迭代閾值收縮算法可提高三維Radon變換域的分辨率,提高一次波與多次波的可分離性,一定程度上提升多次波壓制效果;
(3)地震數(shù)據(jù)橫向振幅的變化也會降低地震數(shù)據(jù)的分辨率,而多項式擬合可保存地震數(shù)據(jù)橫向的振幅變化,高分辨率算法與多項式擬合的結(jié)合可兼顧分辨率和保幅的問題。
(4)在本文所進行的數(shù)據(jù)處理中,為取得更接近實際的處理結(jié)果,需要對保幅高精度算法中涉及提高分辨率的多個參數(shù)進行測試。此外,與常規(guī)最小二乘算法相比,本文算法的計算量偏大,因此需在追求準確處理結(jié)果與采用適量計算耗時之間做一個權(quán)衡。