基于迭代閾值收縮的高分辨率Radon變換方法效果對比

馬繼濤，廖震，齊嬌，遲麟

(中國石油大學(北京) 地球物理學院，北京 102249)

0 引言

Radon 變換是數(shù)學領(lǐng)域的一種變換方法，在地震數(shù)據(jù)處理的速度分析、波場分離和多次波壓制等方面均有應用。但受地震數(shù)據(jù)采集空間有限及方法自身原因限制，其常規(guī)方法的變換域存在平滑效應，且分辨率較低，嚴重影響方法的效果，如何提高變換域分辨率一直是地球物理學家們研究的重點[1]。

Radon變換由Claerbout[2]為首的斯坦福大學地球物理專家于20世紀80年代將其引入地球物理領(lǐng)域。Thorson和Claerbout[3]使用雙曲Radon變換進行速度分析和反演，該方法是在時間—空間域進行求解，需要求取大型矩陣的逆，計算量大；Hampson[4]對拋物線Radon變換的反變換進行了重新定義，并應用最小二乘方法求得地震數(shù)據(jù)的拋物線Radon變換，該方法被首次應用到多次波壓制中，但由于地震數(shù)據(jù)采集的孔徑效應[5]，變換域的分辨率較低。Beylkin[6]對離散Radon變換進行了系統(tǒng)闡述，提出在頻域的最小二乘(LS)計算方法，避免了大型矩陣的求逆過程。

為提高Radon變換的分辨率，Scales等[7]提出用迭代重加權(quán)最小二乘算法[8](IRLS)來解決Thorson和Claerbout構(gòu)建的稀疏逆問題RT模型，Sacchi和Ulrych[9]使用概率密度函數(shù)，基于貝葉斯原理提出了一種在頻域的迭代性稀疏Radon算法FSRT[10](frequency sparse radon transform)，該方法通過增強模型的平滑度來避免放大觀測中產(chǎn)生的隨機誤差，具有較強的抗噪能力，被廣泛的用于地震信號的處理；但該方法在處理小時窗中時差十分接近的同相軸，和具有顯著AVO效應的數(shù)據(jù)時會生成假象；Abbad等[11-12]改進了Sacchi和Ulrych提出的FSRT方法，通過引入新的變量消除了變換算子對頻率的依賴性，加快了迭代速度，提高了計算效率；時間域稀疏Radon變換TSRT(time sparse radon transform)是一種比FSRT分辨率更高的變換方法，此法沿著曲率和截距時間方向上施加稀疏約束[13]，能獲得更具稀疏性的Radon模型；混合時頻域稀疏Radon變換TSRTMD[14](time sparse radon transform mixed domain)在頻率和時間混合域中進行稀疏Radon變換，一般的TSRTMD中稀疏Radon變換通過對最小二乘算法得到的時間域Radon模型進行收斂運算，正向和反向Radon變換均在頻域中實現(xiàn)，Lu[15]在此法基礎(chǔ)上引入迭代對計算進行加速，通過對最小二乘算法得到的Radon模型進行迭代收斂處理，進一步使Radon模型更加接近原始值，以達到提高Radon變換分辨率的效果。經(jīng)模擬數(shù)據(jù)和實際數(shù)據(jù)驗證，表明此方法能夠在很大程度上增強Radon變換的分辨率。該類迭代收縮算法涉及到矩陣的向量乘法運算，能夠在每次迭代后進行收縮步長處理。

迭代收斂系列算法通過在時頻域中進行稀疏Radon變換[16]來分離噪聲信號。概括地說，此類迭代方法中，每次迭代過程中，先對地震信號的Radon域模型乘以變換算子，然后用頻域的信號與其結(jié)果作差，之后乘以變換算子的伴隨矩陣，最后對所得結(jié)果的標量進行收縮處理步驟，實現(xiàn)迭代過程中的Radon模型向原始值靠近，收斂Radon域能量團，以達到提高Radon變換分辨率的效果。迭代算法區(qū)別于一般的時域和頻域Radon變換，由于其是在頻域進行正向和反向Radon變換，時域上進行相關(guān)的迭代收縮步驟，所以它具有時域RT和頻域RT兩者的優(yōu)點，比單一域的RT也更具稀疏性[17]，而且不同于單純的LS算法處理，迭代算法在時域上對LS得到的Radon模型進行迭代，在每次迭代后能得到更為收斂的Radon模型，這對提高Radon域的數(shù)據(jù)精度方面非常有效。

迭代收斂算法系列中有各種迭代收縮方法，它們之間各有差異。例如迭代收縮閾值算法[18](iterative shrinkage thresholding,IST)，IST基礎(chǔ)上的快速迭代收縮閾值算法(fast iterative shrinkage thresholding,FIST)，稀疏加速時不變迭代收縮閾值算法(sparse radon transform iterative shrinkage,SRTIS)等。對IST、FIST、SRTIS這3種迭代方法的系統(tǒng)闡述，了解該3種方法的收縮作用原理以及實現(xiàn)方式；利用模擬數(shù)據(jù)驗證3種算法，并對3種算法的收斂速度等進行對比分析；利用實際地震資料驗證3種迭代收縮算法的收斂效率和處理效果。

1 迭代收縮算法原理

設m為Radon域模型，d為時空域地震數(shù)據(jù)，x為偏移距，q為描述拋物線曲線彎曲程度的曲線參數(shù)，t、τ為時間，對地震道集d(x,t)，可將其視為由m沿特定曲率的拋物線疊加求和[19]而得，該變換的離散形式如下：

(1)

式(1)可寫成算子形式D=LM，利用最小二乘方法求解M，可得：

M=(LLH+μI)-1LD，

(2)

式中，M為Radon域模型，L為變換算子，μ為阻尼因子，D為時空域數(shù)據(jù)。

由于地震數(shù)據(jù)采集空間有限，求和過程只能在有限的時空域進行，最小二乘法求解得到的Radon域數(shù)據(jù)存在剪刀狀發(fā)散的假象，一次波和多次波無法聚焦，使得能量團之間存在相互重疊，從而會導致在Radon域進行切除時對有效信號造成損傷，導致地震數(shù)據(jù)失真。因此提高Radon域數(shù)據(jù)的分辨率一直是該方法研究的熱點。常規(guī)提高分辨率的算法是在頻率域進行的，本文在此討論對比幾種基于迭代收縮閾值的時間域提高分辨率的算法。

1.1 迭代收縮閾值算法原理

常規(guī)的迭代收縮閾值算法(IST)原理如下：

mk+1=

Sλ{mk+αF-1[(LTL)-1LT(F[d]-LF[mk])]}(3)

(4)

其中，mk+1是第k+1次迭代后的Radon模型，mk是第k次迭代后的Radon模型，α是步長，Sλ為閾值算子，F(xiàn)和F-1為傅立葉正、反變換算子，λ為閾值運算所對應的參數(shù)。|m|為Radon模型的絕對值。

從式(3)、(4)可以看出，與最小二乘算法相比(式2)，IST在時間域多了一個迭代的步驟，且每一步迭代僅包含1次簡單的閾值收縮操作，因此能夠大大提高計算效率。在每次閾值處理過程中，增強能量團中心，減弱邊緣的發(fā)散能量，可以使得模型本身能量團更為集中，由此消除平滑效應帶來的分辨率降低的影響。

1.2 快速迭代收縮閾值算法原理

迭代收縮閾值算法盡管可以提高分辨率，但其收斂速度較為緩慢。因此，在常規(guī)迭代收縮閾值算法的基礎(chǔ)上，Beck、Teboulle等[20]提出了一種改進的快速迭代收縮閾值算法FIST。此改進方法大大加快了迭代的收斂速度，其公式與IST相比增加了一迭代計算流程，如式(5)：

mk+1=

(6)

式中，初始值m0=0，t0=1，并且tk由下式計算得到：

(7)

從以上迭代步驟可以看出，F(xiàn)IST在原有算法的基礎(chǔ)上增加了一重迭代計算，增加的算式旨在對Radon域數(shù)據(jù)模型迭代收斂前先進行殘差校正，因此能使相同迭代次數(shù)下FIST加快減小與原始數(shù)據(jù)模型的差距，增大收斂速度。

1.3 稀疏加速時不變迭代收縮閾值算法

FIST雖然在收斂速度上提升了不少，但由于此方法只添加了一重迭代步驟，并沒有從本質(zhì)上改變收斂效果?；赥SRTMD，Lu[21]給出了一種改進的SRTIS算法，此法具備更高的迭代收縮效率，不僅能進一步加快收斂速度，而且也能極大地改善收斂效果，使得變換域的分辨率更高。

迭代收縮算法通常應用全局閾值，在Radon變換中，Radon域模型通常沿著時間維度衰減，因此為遵循模型衰減，將模型收縮用二維方式來展示：

mk=Tα{mk-1+2tF-1[(LTL)-1LT(F[d]-LF[mk-1])]}

(8)

其中，t是步長，k為當前迭代次數(shù)，Tα為一收縮算子，計算公式如下：

(10)

(11)

當均值濾波器大小與模型m相同時，二維收縮將達到最好，閾值將降低到全局閾值。

SRTIS是一個混合的時頻域Radon變換，該法與其他迭代算法類似，均在頻域中完成對地震數(shù)據(jù)的正向和反向Radon變換，在時域中對LS求得的Radon模型進行迭代收縮步驟，其不同點在于式(9)的收縮算法中，SRTIS先對初始Radon模型m0的絕對值作二維均值濾波處理，再將得到的結(jié)果歸一化后作為收縮算子Tα的參數(shù)值，這種處理方式不同于上述的IST和FIST中預先設定的固定閾值參數(shù)，它利用每次迭代后Radon模型的更新來進行動態(tài)閾值，能夠保證每次迭代后的收縮處理更加貼合原始數(shù)據(jù)，因此其收斂效率和計算效果大大優(yōu)于IST。

2 方法應用效果

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

本文利用一個簡單的模擬地震數(shù)據(jù)對比3種方法的處理效果，該數(shù)據(jù)有80道，道間距為20 m，每道采樣點為200個，采樣率為4 ms。圖1a給出了時空域的模擬數(shù)據(jù)，圖中可以看出，數(shù)據(jù)中5條曲線振幅較強，振幅值為1，兩條曲線振幅較弱，振幅值為0.1。圖1b為對應的最小二乘Radon變換方法的變換結(jié)果。

圖1 時空域模擬數(shù)據(jù)(a)及對應的最小二乘Radon變換結(jié)果(b)Fig.1 The simulated data in the time domain(a) and the corresponding result of the LS Radon transform(b)

2.1.1 IST方法處理結(jié)果

對最小二乘Radon變換的結(jié)果，基于IST方法，利用式(3)、式(4)對其進行迭代收縮處理，迭代次數(shù)為40、80、160次時的結(jié)果如圖2所示。

從圖2中可以看出，IST對Radon域的能量團起到了一定的收斂作用，但收斂速度緩慢。

之后對IST方法迭代得到的結(jié)果進行多次波壓制，即將動校時差小于0.04 s的Radon域數(shù)據(jù)(一次波)切除掉，得到多次波壓制結(jié)果，如圖3所示。

比較圖3中的多次波壓制結(jié)果，IST反變換得到的Radon域模型仍與原數(shù)據(jù)存在較大差異，導致多次波壓制結(jié)果存在殘余。

從式(4)可以看出，λ對mk的輸出結(jié)果有所影響，λ數(shù)值應比mk小兩個數(shù)量級，且一定范圍內(nèi)，λ越大，Radon域收斂性越強，具有更高的分辨率，但會導致時間域多次波壓制后的殘余更嚴重。因此選取λ為0.001、0.001 1、0.001 3、0.001 5、0.001 7、0.001 9，迭代次數(shù)為160次時進行對比，且此數(shù)據(jù)中，mk的數(shù)量級為10-2。結(jié)果顯示，將λ值取到0.001 9時能達到較好效果。

a—迭代次數(shù)為40次;b—迭代次數(shù)為80次;c—迭代次數(shù)為160次a—the 40 iterations;b—the 80 iterations;c—the 160 iterations圖2 IST不同迭代次數(shù)結(jié)果對比Fig.2 The iterative results with different iterations of the IST transform

a—最小二乘方法;b—迭代次數(shù)為40次;c—迭代次數(shù)為160次a—the LS transform;b—the 40 iterations;c—the 160 iterations圖3 IST方法多次波壓制結(jié)果對比Fig.3 The comparison of multiple suppression results using IST method

2.1.2 FIST方法處理結(jié)果

FIST方法能降低估計Radon域模型的計算成本，同時加快迭代收斂的速度。基于FIST方法，利用式(5)～(7)對最小二乘方法的Radon變換結(jié)果進行迭代收縮處理，迭代次數(shù)為10、20、40次時的結(jié)果如圖4所示。從圖中可見，與IST 方法相比，F(xiàn)IST收斂速度確實加快了。

FIST方法還在一定程度上消除了IST方法中多次波壓制殘余的現(xiàn)象，在圖5中可以看出，隨著迭代次數(shù)的增加，壓制結(jié)果中的多次波殘余逐漸減小，但仍有部分殘余存在。

2.1.3 SRTIS方法處理結(jié)果

基于SRTIS方法，利用式(8)～(11)對最小二乘的Radon變換結(jié)果進行迭代收縮處理，迭代次數(shù)為5、15、30次時的結(jié)果如圖6所示，對應的多次波壓制結(jié)果圖7所示。

a—迭代次數(shù)為10次;b—迭代次數(shù)為20次;c—迭代次數(shù)為40次a—the 10 iterations;b—the 20 iterations;c—the 40 iterations圖4 FIST迭代結(jié)果Fig.4 The iterative result of the FIST transform

a—迭代次數(shù)為10次;b—迭代次數(shù)為20次;c—迭代次數(shù)為40次a—the 10 iterations;b—the 20 iterations;c—the 40 iterations圖5 FIST多次波壓制結(jié)果對比Fig.5 The comparison of multiple suppression results using FIST method

a—迭代次數(shù)為5次;b—迭代次數(shù)為15次;c—迭代次數(shù)為30次a—the 5 iterations;b—the 15 iterations;c—the 30 iterations圖6 SRTIS迭代結(jié)果Fig.6 The iterative result of the SRTIS transform

a—迭代次數(shù)為5次;b—迭代次數(shù)為15次;c—迭代次數(shù)為30次a—the 5 iterations;b—the 15 iterations;c—the 30 iterations圖7 SRTIS多次波壓制結(jié)果對比Fig.7 The comparison of multiple suppression results using SRTIS method

圖4和圖6對比可見，SRTIS在收斂速度上略快于FIST，在迭代次數(shù)15次時，SRTIS已經(jīng)能看到清晰的數(shù)據(jù)點，能量充分聚集，平滑效應被消除，但兩個弱振幅點并不能顯示，在迭代次數(shù)達到30次時能看到7個數(shù)據(jù)點。從圖7可見在迭代次數(shù)為30次時多次波殘余基本消失，說明SRTIS的壓制效果最佳。在上述SRTIS結(jié)果中，α和t均選擇同一數(shù)值，α=0.04，t=0.8，此參數(shù)設置為經(jīng)過驗證后的最佳試驗數(shù)值。

2.1.4 3種方法的結(jié)果對比

為對比3種方法的數(shù)據(jù)處理能力，本文對3種迭代算法的收斂效率進行了比較，3種算法1次到20次迭代后的模型mk與原始模型m的誤差能量結(jié)果如圖9所示。圖中可見SRTIS下降速率最快，F(xiàn)IST次之，IST最慢。說明SRTIS方法的收斂性優(yōu)于其他兩種，改進后的FIST優(yōu)于IST。

圖9 3種迭代算法的收斂效率對比Fig.9 The compasion of three iterative algorithms' contraction efficiency

為比較三者之間的收斂速度，用模擬數(shù)據(jù)以及一總長度為7 s、采樣間隔為0.004 s的CDP道集計算達到10-1級別的精度標準時所耗費的時間，結(jié)果見表1所示。

從表1中可見IST方法的收斂速度非常緩慢，而FIST與SRTIS收斂速度沒有很大差別，但是在保幅性能上，F(xiàn)IST并不能很好的達到要求，在實際資料的處理中，較易破壞同相軸，因此綜合考慮，SRTIS的收斂效率較高。

此外，用一相近拋物線試驗這些方法的壓制干擾能力，其動較時差為0.028 s，壓制動校時差大于0.02 s的拋物線，兩者相差僅為0.008 s。4種方法中，最小二乘算法殘余嚴重；IST能在一定程度上切除掉干擾拋物線，但是在尾端仍有明顯的殘余；FIST尾端的殘余現(xiàn)象有得到改善；SRTIS壓制效果十分明顯，壓制后已經(jīng)很難觀測到干擾的拋物線痕跡。

表1 3種迭代算法的收斂結(jié)果對比

以上結(jié)果均說明SRTIS優(yōu)于其他兩種迭代收縮算法，能較大程度上提高Radon變換的分辨率，且能較為徹底地對數(shù)據(jù)中剩余時差較小的多次波進行壓制，因此方法對壓制剩余時差較小的層間多次波較為有利。

2.2 實際地震資料處理結(jié)果

為進一步驗證迭代收縮算法的實用性，本文引入實際地震資料，基于3種迭代收縮方法與傳統(tǒng)最小二乘算法進行對比。實際資料為墨西哥灣的某多次波干擾嚴重的CDP道集。

利用最小二乘方法進行了Radon變換和多次波的壓制，并基于該變換結(jié)果和上文所述的各種迭代收縮方法進行提高分辨率運算，迭代次數(shù)均為40次，得到高分辨率的Radon域結(jié)果；對所得到的Radon域數(shù)據(jù)進行一次波和多次波的分離，并將分理出的多次波變換回時空域，將其由原數(shù)據(jù)中減去，得到多次波壓制后的道集數(shù)據(jù)。原始數(shù)據(jù)、LS、IST、FIST和SRTIS各方法估計的多次波模型如圖10所示，多次波壓制后的結(jié)果如圖11所示。

圖10 各方法所估計的多次波模型Fig.10 The estimated multiple model by various methods

圖11 各方法多次波壓制后的結(jié)果Fig.11 The multiple suppression result by various methods

上述方法的處理效果中，LS方法利用最小二乘原理取得Radon域的變換結(jié)果，但最小二乘算法僅為局部收斂，因此估計的多次波與原始數(shù)據(jù)仍有較大差異，壓制的時候會影響到一次波的反變換結(jié)果；IST雖然提高了分辨率，但其公式中并未將每次計算迭代后的殘差量帶回，因此導致所估計多次波與原數(shù)據(jù)有一定差異，且IST方法遠偏移距箭頭所指一次波也受到影響；FIST方法雖能進一步提高分辨率，也考慮了殘差量，但因其迭代程度相對過深，導致Radon域數(shù)據(jù)失真；SRTIS因為將每次迭代后的數(shù)據(jù)模型如式(10)歸一化后代入到收縮閾值步驟中，在一定程度上能夠約束模型的收斂方向，使得其更靠近真解，因此估計的多次波較為真實，壓制效果也與實際情況較為吻合，既提高了Radon域分辨率，也得到了相對較好的壓制結(jié)果，近偏移距和遠偏移距效果均不錯。且從圖11壓制效果中可見，圖11e的結(jié)果最好的保持了原數(shù)據(jù)的振幅特性，而圖11c、d中一次波的橫向振幅特性遭到破壞，導致一次波水平同相軸較多，不符實際。

為突出振幅能量，本文將處理后的CDP道集疊加用于對比，疊加道集效果如圖12所示。

對比圖11的單道壓制效果，圖12疊加道集在4.3 s附近的壓制效果較為明顯。圖12中，LS方法的疊加大致能夠?qū)?.5～3.7 s出現(xiàn)的干擾同相軸去除，但仍殘留著多次波的干擾能量，這是該方法的理論公式所帶來的誤差影響，是無法消除的；IST的疊加雖然提高了分辨率，但由于在壓制過程中受到多次波殘余誤差的影響，并不能很好地辨識出一次波的同相軸，信噪比較低，僅能在4～4.3 s模糊地看到3條同相軸，且4.7 s和5.1 s附近的同相軸聚焦并不是很好；FIST方法的疊加進一步提高了分辨率，對多次波的壓制程度也更大，但因其迭代程度相對過深，造成的信噪比損失也更為嚴重，因此4～4.3 s處的同相軸周圍干擾影響更為嚴重；SRTIS的疊加因為其方法在每次迭代過程中都能與原始數(shù)據(jù)模型保持一致性，所以既提高了分辨率，也并沒有對信噪比造成較大破壞，平衡了分辨率和信噪比，在圖中能夠清楚地看到4～4.3 s的3條同相軸，且4.5 s和5.1 s附近的同相軸相對集中，于多次波壓制為較好的結(jié)果。綜上所述，SRTIS在實際地震資料的多次波壓制應用中依然強于其他兩種迭代方法。

3 結(jié)論

本文基于IST、FIST和SRTIS這3種迭代收縮Radon變換算法壓制多次波，并給出了算法實現(xiàn)的計算步驟，模型數(shù)據(jù)和實際資料的處理驗證了這3種方法的正確性和可行性，得出如下結(jié)論：

1)傳統(tǒng)LS方法在用于壓制多次波方面上存在局限性，而迭代收縮算法通過迭代中的閾值收縮操作來逼近原始值，在每次閾值處理過程中，增強能量團中心，削弱邊緣的發(fā)散能量，使模型能量團集中，由此消除平滑效應帶來的降低分辨率的影響。

2)SRTIS的壓制效果均優(yōu)于IST和FIST，其先對初始Radon模型m0的絕對值作二維均值濾波處理，再將得到的結(jié)果歸一化后作為收縮算子Tα的參數(shù)值。不同于IST和FIST的固定閾值參數(shù)，SRTIS利用了每次迭代后Radon模型的更新來進行動態(tài)閾值，能夠保證每次迭代后的收縮處理更加貼合原始數(shù)據(jù)。

3)參數(shù)的選擇對于算法的實現(xiàn)效果有很大影響，因此在使用算法處理實際數(shù)據(jù)時需先對測試數(shù)據(jù)大小進行評估，才能達到較好的收縮閾值效益。