劉仕友， 馬繼濤， 孫萬(wàn)元， 應(yīng)明雄

(1. 中海石油(中國(guó))有限公司 湛江分公司，湛江 524057；2. 中國(guó)石油大學(xué)(北京) 地球物理學(xué)院物探系，北京 102249)

0 引言

多次波是海洋地震數(shù)據(jù)中常見(jiàn)的噪音，有效去除多次波是地震資料處理中的一個(gè)重要研究?jī)?nèi)容。Radon變換是工業(yè)界多次波壓制最為常用的有效手段之一，同時(shí)也被廣泛用于地震數(shù)據(jù)重建，波場(chǎng)分離，速度分析等。

奧地利數(shù)學(xué)家Radon在1917年提出Radon變換，美國(guó)斯坦福大學(xué)地球物理小組在20世紀(jì)70年代對(duì)Radon變換進(jìn)行了重點(diǎn)研究，并取得了顯著進(jìn)展，為其在地球物理領(lǐng)域的應(yīng)用做出了重要貢獻(xiàn)。Hampson[1]將線性Radon變換改進(jìn)為拋物線Radon變換，并將其應(yīng)用于多次波的壓制中；拋物線Radon變換利用一次波和多次波動(dòng)校正速度的差異進(jìn)行多次波壓制。在動(dòng)校正后，地震數(shù)據(jù)的一次波基本被校平，多次波同相軸則由雙曲線被部分動(dòng)校為拋物線，通過(guò)拋物線Radon變換可以把時(shí)空域中的曲線映射為Radon域中的點(diǎn)，多次波和一次波對(duì)應(yīng)曲線曲率具有較大的差異，從而在Radon域可以實(shí)現(xiàn)一次波和多次波的分離，達(dá)到去除多次波的目的。Beylkin[2]對(duì)基于最小二乘理論的Radon變換方法進(jìn)行了討論，由于Radon變換算法自身和地震數(shù)據(jù)采集空間和時(shí)間有限等原因，Radon域的數(shù)據(jù)往往會(huì)出現(xiàn)能量發(fā)散的情況，這嚴(yán)重影響了該方法多次波壓制的效果，因此發(fā)展高分辨率Radon變換算法是十分必要的。Sacchi等[3]提出高精度的Radon變換，奠定了高精度Radon變換的理論基礎(chǔ)。迭代高分辨率算法在處理稀疏采樣數(shù)據(jù)時(shí)分辨率降低，其原因在于有限的空間采樣會(huì)限制Radon變換的分辨率，產(chǎn)生空間假頻，同時(shí)該算法進(jìn)行了多次迭代運(yùn)算，計(jì)算量較大。Herrmann[4]提出一種不需要迭代高分辨率算法，用數(shù)據(jù)的低頻部分來(lái)約束數(shù)據(jù)的拋物線分解，低頻結(jié)果在運(yùn)算中不會(huì)產(chǎn)生假頻，能夠在采樣稀疏的情況下仍取得比較好的高分辨率效果。

國(guó)內(nèi)的專家學(xué)者對(duì)高分辨率Radon變換也展開(kāi)過(guò)較為深入詳細(xì)的研究。吳律[5]系統(tǒng)闡述了Tau-p的原理及應(yīng)用情況；牛濱華[6]等提出了基于多項(xiàng)式理論的Radon變換；許多學(xué)者對(duì)Radon變換及其提高分辨率的方法做了系統(tǒng)深入地研究[7-13]。

筆者對(duì)拋物線Radon變換、基于迭代的高分辨率的Radon變換，及基于低頻約束的高分辨率Radon變換方法進(jìn)行系統(tǒng)闡述，并對(duì)比了這些方法在地震數(shù)據(jù)多次波壓制中的效果。模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)的對(duì)比可以看出，相比于普通最小二乘Radon變換和迭代高分辨率算法，基于低頻約束的高分辨率Radon變換算法分辨率更高，多次波壓制更加徹底。

1 方法原理

1.1 Radon變換原理

二維拋物線Radon變換是沿著拋物線對(duì)地震數(shù)據(jù)進(jìn)行疊加得到Radon域中的數(shù)據(jù)，二維Radon變換可表示為：

(1)

其中:m為Radon域數(shù)據(jù)；d為時(shí)空域地震數(shù)據(jù)；t、τ為時(shí)間(一個(gè)是Radon域的時(shí)間，一個(gè)是時(shí)空域的時(shí)間)；q為描述曲線彎曲程度的曲線參數(shù)；x為偏移距。對(duì)式(1)進(jìn)行離散化，可得式(2)。

(2)

其中：Nq為Radon域數(shù)據(jù)的曲率個(gè)數(shù)。對(duì)兩邊做傅立葉變換，轉(zhuǎn)換到頻率域中進(jìn)行計(jì)算，

(3)

將式(3)寫為算子形式，則拋物線Radon變換可以表示為：

D=LM

(4)

其中，算子L可以寫為：

(5)

因此，為得到Radon域的數(shù)據(jù)M，可以建立如式(6)目標(biāo)函數(shù)進(jìn)行求解。

(6)

將式(6)最小化，可以得到M的最小二乘解，即：

(7)

其中,μ是阻尼因子，取值一般在0.1～1之間。由式(7)可以看出，在進(jìn)行Radon域M的計(jì)算時(shí)，對(duì)所有頻率μ均為相同的定值，其目的是為了提高矩陣運(yùn)算的穩(wěn)定性，但該因子同時(shí)也對(duì)Radon域的數(shù)據(jù)進(jìn)行了平滑，降低了該域數(shù)據(jù)的分辨率。

1.2 基于迭代的高分辨率Radon變換

因?yàn)槌Ｒ?guī)Radon變換在Radon域數(shù)據(jù)能量發(fā)散，影響多次波壓制的效果，因此Sacchi等發(fā)展了基于迭代的高分辨率Radon變換，該方法建立了如下的目標(biāo)函數(shù)：

(8)

其中，R(M)是施加在Radon域數(shù)據(jù)M上的規(guī)則化因子。對(duì)式(8)最小化，可得到Radon域數(shù)據(jù)的計(jì)算公式如下：

M(ω)=(L(ω)HL(ω)+μQ(M))-1L(ω)HD(ω)

(9)

其中，施加在矩陣LHL上的規(guī)則化因子與Radon域數(shù)據(jù)M存在如下的關(guān)系：

(10)

即某個(gè)頻率的加權(quán)矩陣是通過(guò)某規(guī)則，由Radon域的計(jì)算結(jié)果M給出的；如果Radon域計(jì)算結(jié)果M中的某部分能量較強(qiáng)，則加權(quán)矩陣Q較小，如果Radon域計(jì)算結(jié)果M中的某部分能量較弱，則加權(quán)矩陣較大，這樣通過(guò)式(9)中的矩陣求逆過(guò)程，會(huì)使得Radon域能量強(qiáng)的點(diǎn)更強(qiáng)，能量弱的點(diǎn)減弱，從而使得Radon域能量聚焦，達(dá)到高分辨率的效果。該算法是一種基于貝葉斯原理的空間稀疏約束算法，通過(guò)迭代運(yùn)算對(duì)Radon域內(nèi)的能量進(jìn)行聚焦，實(shí)際計(jì)算時(shí)，可采取如下的迭代公式進(jìn)行運(yùn)算：

Mk(ω)=(L(ω)HL(ω)+

λQ(Mk-1(ω)))-1L(ω)HD(ω)

(11)

其中，k為迭代次數(shù)?？梢钥闯觯炒蔚募訖?quán)矩陣Q，是由上一次迭代的運(yùn)算結(jié)果M給出的。加權(quán)矩陣可以根據(jù)上一次迭代結(jié)果的Radon域數(shù)據(jù)M的大小直接賦值，例如，

(12)

其中，ε為使得求倒數(shù)穩(wěn)定所施加的一個(gè)因子，可以通過(guò)式(13)進(jìn)行計(jì)算。

ε=max(0.01*max|M|,1×10-6)

(13)

該算法對(duì)每個(gè)頻率都需要迭代進(jìn)行，運(yùn)算量較大；如果想要取得較為滿意的結(jié)果，一般每個(gè)頻率都需要迭代數(shù)十次。如果處理的數(shù)據(jù)體很大，需要耗費(fèi)大量的運(yùn)算時(shí)間。

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

由前面可知，式(9)中的加權(quán)因子Q，可以由上一次迭代的計(jì)算結(jié)果給出。Herrmann對(duì)公式(9)進(jìn)行了改進(jìn)，由上一個(gè)頻率的計(jì)算結(jié)果計(jì)算加權(quán)因子Q[4]。該算法的一個(gè)優(yōu)勢(shì)是避免了每個(gè)頻率為得到加權(quán)矩陣進(jìn)行的迭代運(yùn)算，對(duì)于每個(gè)頻率，直接代入由上個(gè)頻率計(jì)算結(jié)果得到的加權(quán)矩陣，就可以得到該頻率的計(jì)算結(jié)果。另外，由于使用了低頻的計(jì)算結(jié)果對(duì)高頻的計(jì)算進(jìn)行了約束，該算法可以抗假頻。

非迭代的，基于低頻約束的高分辨率Radon變換可表示為：

M(ωn)=(L(ωn)HL(ωn)+

μQ(M(ωn-1)))-1L(ωn)HD(ωn)

(14)

其中：Q為對(duì)角加權(quán)矩陣，由上一頻率的計(jì)算結(jié)果得到；角標(biāo)n表示第n個(gè)頻率；式中的阻尼因子μ的取值范圍變化較大，可以通過(guò)測(cè)試得到，一般的取值范圍為0.01～1；Q可表示為：

Qii(ωn)=‖Mi(ωn-1)‖,i=1,2…Nq

(15)

其中，Mi為在上一頻率運(yùn)算得到計(jì)算結(jié)果。從式(14)、式(15)中可以看出，某個(gè)頻率的加權(quán)矩陣，是由上一個(gè)頻率的計(jì)算結(jié)果給出的。在沒(méi)有多次波曲率先驗(yàn)信息的情況下，該方法更看重低頻結(jié)果對(duì)整體結(jié)果的約束，低頻結(jié)果在運(yùn)算中不會(huì)產(chǎn)生假頻，能夠在采樣稀疏的情況下仍取得比較好的高分辨率效果。

基于低頻約束的高分辨率Radon變換進(jìn)行多次波壓制的處理步驟如下：

1)通過(guò)傅立葉變換將地震數(shù)據(jù)轉(zhuǎn)換到f-x域。

2)給初始對(duì)角加權(quán)矩陣賦值為單位矩陣。

3)對(duì)低頻數(shù)據(jù)進(jìn)行處理，將處理結(jié)果在矩陣中保存。

4)利用上一頻率的數(shù)據(jù)計(jì)算結(jié)果計(jì)算新的加權(quán)矩陣，用來(lái)進(jìn)行下一個(gè)頻率數(shù)據(jù)的計(jì)算，直至所有頻率都計(jì)算完成。

圖1 模擬數(shù)據(jù)Fig.1 Simulated data

圖2 常規(guī)最小二乘算法結(jié)果Fig.2 The results of regular least-squares result(a)Radon變換結(jié)果;(b)多次波壓制后的結(jié)果

5)反傅立葉變換，并對(duì)計(jì)算結(jié)果進(jìn)行一次波切除。

6)進(jìn)行Radon反變換，完成相減運(yùn)算，得到多次波壓制結(jié)果。

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

2.1 模擬數(shù)據(jù)測(cè)試

利用一個(gè)模擬數(shù)據(jù)對(duì)三種Radon變換方法的分辨率和多次波壓制效果進(jìn)行了對(duì)比。首先生成一個(gè)模擬數(shù)據(jù)，該模擬數(shù)據(jù)由兩個(gè)一次波和兩個(gè)對(duì)應(yīng)的多次波組成，一次波分別位于0.4 s、0.9 s，多次波分別位于0.8 s、1.2 s，兩個(gè)多次波的動(dòng)校時(shí)差量分別為330 ms、480 ms,該模擬數(shù)據(jù)如圖1所示。

利用基于最小二乘的拋物線Radon變換對(duì)該數(shù)據(jù)進(jìn)行多次波壓制，在進(jìn)行多次波切除時(shí)，選擇的Radon域切除的q值為50 ms，以減小對(duì)一次波的損害。該方法在Radon域的變換結(jié)果，以及多次波壓制后的結(jié)果如圖2所示。

圖3 迭代方法提高分辨率結(jié)果Fig.3 The high resolution results based iterative method(a)Radon域結(jié)果；(b)多次波壓制后的結(jié)果

從圖2中可以明顯看出，基于最小二乘理論的拋物線Radon變換在變換域具有明顯的剪刀狀發(fā)散現(xiàn)象，由于剪刀狀能量發(fā)散，導(dǎo)致其壓制結(jié)果具有明顯的多次波殘余；另外在近偏移距處還產(chǎn)生了較為明顯的假象。

我們利用基于迭代的高分辨率Radon變換對(duì)該模擬數(shù)據(jù)進(jìn)行了多次波壓制，其結(jié)果如圖3所示。

從圖3中可以看出，基于迭代的高分辨率Radon變換在Radon域，沒(méi)有了剪刀狀發(fā)散現(xiàn)象，分辨率有了明顯提高，但該結(jié)果在一次波下方的1.4 s處，以及多次波附近的0.9 s和1.2 s處，出現(xiàn)了假象。我們將迭代次數(shù)分別選取為10、20、40次，阻尼因子分別選取為0.01、0.05、0.1，該假象仍然存在。因此可以認(rèn)為該假象有可能是由剪刀狀發(fā)散所引起的，因?yàn)樵摰岣叻直媛实姆椒ǎ跏嫉凳亲钚《薘adon變換的結(jié)果，該方法在進(jìn)行迭代時(shí)，誤將剪刀狀發(fā)散的能量認(rèn)為是有效能量進(jìn)行進(jìn)一步的聚焦，從而產(chǎn)生了本不存在的能量假象。該方法多次波壓制過(guò)程中采用了和最小二乘Radon相同的切除參數(shù)，結(jié)果如圖3(b)所示，可以看出，壓制效果有了明顯改善，但遠(yuǎn)偏移距的多次波仍有一些殘余。

采用基于低頻約束的拋物線Radon變換對(duì)該數(shù)據(jù)進(jìn)行了多次波壓制處理(圖4)。從圖4可以看出，基于低頻約束的拋物線Radon變換變換域的分辨率很高，而且不存在假象；多次波壓制后的結(jié)果優(yōu)于其他兩種方法的結(jié)果。

圖4 低頻約束提高分辨率結(jié)果Fig.4 The high results based on lower frequency constraint(a) Radon域結(jié)果；(b)多次波壓制結(jié)果

圖5 實(shí)際數(shù)據(jù)炮集Fig.5 The real shot gather(a)帶有多次波的實(shí)際數(shù)據(jù)；(b)局部放大圖

在計(jì)算時(shí)間上，對(duì)于此模擬數(shù)據(jù)，基于最小二乘的拋物線Radon變換和基于低頻約束的高分辨率Radon變換分別耗時(shí)0.40 s和0.51 s，基于迭代的高分辨率Radon變換運(yùn)算時(shí)間和所選取的迭代次數(shù)有關(guān)，本文運(yùn)算中選取迭代次數(shù)與運(yùn)算時(shí)間的關(guān)系如表1所示，從表1可以看出，基于迭代的高分辨率Radon變換運(yùn)算時(shí)間在迭代次數(shù)是1時(shí)，運(yùn)算時(shí)間和其他兩種方法相比略大，但隨著迭代次數(shù)增加，運(yùn)算所需時(shí)間也隨著大大增加。對(duì)比可以看出，基于低頻約束的高分辨率Radon變換計(jì)算速度快，計(jì)算結(jié)果最優(yōu)。

表1 基于迭代的高分辨率Radon變換 運(yùn)算時(shí)間和迭代次數(shù)表Tab.1 The iteration number and computation time of iterative based high resolution method

2.2 實(shí)際數(shù)據(jù)測(cè)試

利用墨西哥灣的一個(gè)CDP道集對(duì)三種拋物線Radon變換算法進(jìn)行了測(cè)試(圖5)，從圖5中可以看出，一次波已經(jīng)被校平，未校平的多次波在4 s以下。為更好地對(duì)結(jié)果進(jìn)行對(duì)比，選取了該道集的3.2 s～4.8 s處進(jìn)行了放大。

圖6 實(shí)際數(shù)據(jù)最小二乘結(jié)果Fig.6 The least-squares results of real data(a)Radon域結(jié)果；(b)多次波壓制結(jié)果

圖7 實(shí)際數(shù)據(jù)迭代方法提高分辨率結(jié)果Fig.7 The real data high resolution results based iterative method(a)Radon域結(jié)果；(b)多次波壓制結(jié)果

圖8 實(shí)際數(shù)據(jù)低頻約束提高分辨率結(jié)果Fig.8 The real data high results based on lower frequency constraint(a)Radon域結(jié)果；(b)多次波壓制結(jié)果

我們分別利用三種拋物線Radon變換方法對(duì)CDP道集進(jìn)行了處理(圖6、圖7、圖8)。從圖6～圖8中可以看出，基于低頻約束的拋物線Radon變換在變換域具有更高的分辨率，基于最小二乘理論的拋物線Radon變換，和其他兩種方法相比，其壓制結(jié)果在中遠(yuǎn)偏移距處有少量的多次波殘余(箭頭所指)。基于迭代的高分辨率Radon變換(迭代20次)和基于低頻約束的高分辨率Radon變換具有相近的多次波壓制結(jié)果，但后者計(jì)算速度要快很多，對(duì)于該CDP道集，前者用時(shí)23.68 s，而后者僅用時(shí)3.49 s。在處理大批量實(shí)際數(shù)據(jù)時(shí)，基于低頻約束的高分辨率Radon變換算法具有明顯的優(yōu)勢(shì)。

3 結(jié)論

筆者對(duì)拋物線Radon變換及其提高變換域分辨率的兩種方法進(jìn)行了探討，提出一種基于低頻約束的高分辨率Radon變換算法，并用模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)對(duì)算法進(jìn)行了測(cè)試，通過(guò)分析和研究得出以下的結(jié)論：

1)基于迭代算法的高分辨率Radon變換利用前一次迭代變換域的結(jié)果對(duì)下一次迭代的計(jì)算進(jìn)行約束，而基于低頻約束的高分辨率Radon變換是一種非迭代的算法，它利用前一個(gè)頻率的計(jì)算結(jié)果對(duì)下一個(gè)頻率的計(jì)算進(jìn)行約束。

2)基于低頻約束的高分辨率Radon變換可以取得更高分辨率的結(jié)果和更好的多次波壓制結(jié)果。

3)與以迭代方式提高分辨率的算法相比，低頻約束的高分辨率Radon變換算法計(jì)算時(shí)間大大縮短，在實(shí)際生產(chǎn)中可以極大地提高計(jì)算效率和生產(chǎn)工作效率。