基于低頻約束的高分辨率Radon變換多次波壓制方法研究

劉仕友， 馬繼濤， 孫萬元， 應明雄

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

0 引言

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

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

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

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

1 方法原理

1.1 Radon變換原理

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

(1)

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

(2)

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

(3)

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

D=LM

(4)

其中，算子L可以寫為：

(5)

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

(6)

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

(7)

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

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

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

(8)

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

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

(9)

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

(10)

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

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

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

(11)

其中，k為迭代次數(shù)。可以看出，某次迭代的加權(quán)矩陣Q，是由上一次迭代的運算結(jié)果M給出的。加權(quán)矩陣可以根據(jù)上一次迭代結(jié)果的Radon域數(shù)據(jù)M的大小直接賦值，例如，

(12)

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

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

(13)

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

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

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

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

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

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

(14)

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

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

(15)

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

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

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

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

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

4)利用上一頻率的數(shù)據(jù)計算結(jié)果計算新的加權(quán)矩陣，用來進行下一個頻率數(shù)據(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)反傅立葉變換，并對計算結(jié)果進行一次波切除。

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

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

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

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

利用基于最小二乘的拋物線Radon變換對該數(shù)據(jù)進行多次波壓制，在進行多次波切除時，選擇的Radon域切除的q值為50 ms，以減小對一次波的損害。該方法在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ā)散，導致其壓制結(jié)果具有明顯的多次波殘余；另外在近偏移距處還產(chǎn)生了較為明顯的假象。

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

從圖3中可以看出，基于迭代的高分辨率Radon變換在Radon域，沒有了剪刀狀發(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，該假象仍然存在。因此可以認為該假象有可能是由剪刀狀發(fā)散所引起的，因為該迭代提高分辨率的方法，初始迭代值是最小二乘Radon變換的結(jié)果，該方法在進行迭代時，誤將剪刀狀發(fā)散的能量認為是有效能量進行進一步的聚焦，從而產(chǎn)生了本不存在的能量假象。該方法多次波壓制過程中采用了和最小二乘Radon相同的切除參數(shù)，結(jié)果如圖3(b)所示，可以看出，壓制效果有了明顯改善，但遠偏移距的多次波仍有一些殘余。

采用基于低頻約束的拋物線Radon變換對該數(shù)據(jù)進行了多次波壓制處理(圖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ù)據(jù)炮集Fig.5 The real shot gather(a)帶有多次波的實際數(shù)據(jù)；(b)局部放大圖

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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