基于貪婪—快速迭代收縮閾值的Radon變換及其在多次波壓制中的應用

張 全 雷 芩 林柏櫟 彭 博*② 劉書妍②

(①西南石油大學計算機科學學院,四川成都 610500；②油氣藏地質(zhì)及開發(fā)工程國家重點實驗室(西南石油大學)，四川成都 610500；③電子科技大學信息與通信工程學院，四川成都 611731)

0 引言

1917年，奧地利數(shù)學家Radon[1]首次提出Radon變換，經(jīng)過眾多學者的研究，將其引入各領域，并得到了廣泛的使用。20世紀70年代，Claerbout等[2]首次將Radon變換引入地震數(shù)據(jù)處理，為Radon變換在地震勘探中的發(fā)展奠定了基礎。

Radon變換將地震數(shù)據(jù)沿特定路徑進行曲線積分，目的是將數(shù)據(jù)中的信號疊加為Radon域的稀疏信號，便于識別與分離。為了提高Radon變換壓制多次波的精度，不少學者進行了相關的研究。Thorson等[3]提出了Radon變換最小二乘法和隨機反演的理論，將Radon變換建模為一個稀疏逆問題，使Radon變換的分辨率得到一定程度的提高。李元欽等[4]結合平面映射，提出雙曲Radon正、反變換，減少變換噪音和能量彌散，提高了地震資料的處理質(zhì)量。Sacchi等[5-6]在頻率域利用柯西稀疏約束提高了Radon變換的分辨率。牛濱華等[7]提出次數(shù)為“2”的多項式Radon變換。劉喜武等[8]采用最小二乘反演法實現(xiàn)高分辨率拋物線Radon變換和雙曲線Radon變換，提高了精度和分辨率。Ethan等[9]提出加權最小二乘Radon變換，在壓制多次波過程中能較好地保持振幅。劉保童等[10]提出一種頻率域Radon變換方法，有效地提高了地震數(shù)據(jù)在Radon域的分辨率。Schonewille等[11]證明了Radon變換在時間域比在頻率域更加稀疏，故Radon變換在時間域的分辨率比在頻率域高，提高解復用效率。薛亞茹等[12]提出基于Radon變換和正交多項式變換的多方向正交多項式變換的方法，提高了對一次波與多次波的剩余時差的分辨率，在有效壓制多次波的同時保留了一次波的AVO特性，增強了Radon變換的保幅特性。Lu[13]提出基于迭代二維模型收縮的混合時頻域快速稀疏時不變Radon變換方法，利用混合時頻域的優(yōu)勢，該算法既在時間域中提高了Radon模型的稀疏性，又在頻率域得到更高的計算效率。薛亞茹等[14]提出高分辨率稀疏Radon變換和正交多項式變換結合的高階稀疏Radon變換，能有效分離一次波和多次波。Sun等[15]提出基于光滑的 L0范數(shù)的高階、高分辨率頻率-曲率域Radon變換，將正交多項式與Radon變換相結合，該算法具有更好的聚焦特性，并在壓制多次波的同時能更好地保存一次波的AVO特性，且有更高的計算效率。薛亞茹等[16]提出應用加權迭代軟閾值算法的高分辨率Radon變換，將迭代重加權最小二乘法的加權矩陣思想引入迭代軟閾值算法中，提高了Radon變換分辨率。孫寧娜等[17]提出基于多窗口聯(lián)合優(yōu)化的多次波自適應相減方法，與傳統(tǒng)的基于單窗口匹配的多次波自適應相減方法相比，該方法可更有效地壓制殘余多次波并保護一次波，且計算效率與基于小窗口匹配的傳統(tǒng)多次波自適應相減方法相當。馬繼濤等[18]提出三維高精度保幅Radon變換多次波壓制方法，該方法可獲得高分辨率的模型域數(shù)據(jù)，能有效分離一次波與多次波，同時多項式擬合可保護有效波的振幅，高保真地實現(xiàn)多次波的壓制。在Radon變換壓制多次波的速度提升方面，很多學者也進行了相關的研究。Hampson[19]提出拋物線Radon變換壓制多次波，進一步提升了計算效率。熊登等[20]提出一種新的混合域高分辨率拋物線Radon變換，對于多次波的壓制既有很高的計算效率又有較高的分辨率。Brahim等[21]提出一種基于奇異值分解的改進的拋物線Radon變換，克服了頻域Radon變換實現(xiàn)所需的繁重計算，具有更快的計算速度。張全等[22]提出CPU-GPU異構平臺的拋物線Radon變換并行算法，大幅提高了Radon變換壓制多次波的計算效率。隨著深度學習的發(fā)展，不少學者將多次波壓制與深度學習相結合。劉小舟等[23]提出數(shù)據(jù)增廣的編解碼卷積網(wǎng)絡地震層間多次波壓制方法，結合去噪卷積神經(jīng)網(wǎng)絡(DnCNN)和U形全卷積神經(jīng)網(wǎng)絡(U-Net)的優(yōu)勢，搭建了適合層間多次波壓制的深層編、解碼網(wǎng)絡，并且利用數(shù)據(jù)增廣提高網(wǎng)絡的泛化能力，該網(wǎng)絡能夠很好地壓制層間多次波、保護一次波，提高了多次波的壓制效率。

本文在Lu[13]的研究基礎上，將Liang等[24]提出的貪婪—快速迭代收縮閾值算法(Greedy Fast Ite-rative Shrinkage Thresholding Algorithm，Greedy FISTA)應用于混合域拋物線Radon變換，構建了一種快速稀疏時不變Radon變換進行地震多次波壓制。實驗結果表明，本文方法有效地提高了Radon變換壓制多次波的精度和效率。

1 Radon變換多次波壓制基本原理

在地震數(shù)據(jù)處理中，常采用拋物線Radon變換實現(xiàn)一次波與多次波分離。動校正后，道集中的一次波被拉平，多次波呈“拋物線”形態(tài)，根據(jù)二者在Radon域的形態(tài)差異實現(xiàn)多次波壓制。時域與Radon域地震數(shù)據(jù)可用數(shù)學模型表示為

d=Lm

(1)

式中：d為時空域地震數(shù)據(jù)；m為地震數(shù)據(jù)在Radon域的表示；L為Radon變換算子矩陣。通常情況下，d為已知數(shù)據(jù)，m為未知數(shù)據(jù)，L通常不可逆，常用共軛轉(zhuǎn)置矩陣LH近似L的逆矩陣，故式(1)在數(shù)學上為一逆問題。求解逆問題的經(jīng)典方法為最小二乘(LS)法

(2)

L通常是病態(tài)的，且LS法不適用于求解病態(tài)方程，常添加正則化項解決此問題?；贚1范數(shù)的稀疏正則化是目前常用的提高Radon變換分辨率的方法

(3)

式中λ>0為正則化系數(shù)。在時域和頻域構成的混合域(簡稱混合域)實現(xiàn)Radon變換

d=F-1[LF(m)]

(4)

式中：F(·)為傅里葉變換算子；F-1(·)為傅里葉逆變換算子。頻率域的Radon變換算子為[26]

(5)

式中：i=0，1，…，Nq，其中Nq為曲率q的個數(shù)；j=0，1，…，Nx，其中Nx為地震數(shù)據(jù)道數(shù)；f為頻率；xj為第j道的炮檢距。

Radon變換壓制多次波在混合域的目標方程為

(6)

該目標方程的求解是影響多次波壓制效率的重要因素。

2 Radon變換算法及其改進

對式(6)的求解，傳統(tǒng)的迭代重加權最小二乘(Iterative Reweighted Least Squares，IRLS)算法涉及大規(guī)模的求逆運算，耗時長。而基于梯度的迭代算法，計算復雜度小，且算法結構簡單。本文主要將當前應用較廣的迭代收縮閾值算法(Iterative Shrinkage-Thresholding Algorithm，ISTA)、快速迭代收縮閾值算法(Fast Iterative Shrinkage-Thresho-lding Algorithm，F(xiàn)ISTA)以及Greety FISTA與拋物線Radon變換結合對多次波進行壓制。

2.1 基于迭代收縮閾值的Radon變換算法

ISTA是經(jīng)典梯度算法的擴展，計算簡單，是目前最常用的方法之一。Lu[13]將二維ISTA引入Radon變換，提出了基于迭代二維模型收縮的混合域快速稀疏時不變Radon變換(an accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on iterative 2D model shrinkage，SRTIS)。SRTIS算法流程如下。

(1)輸入時空域地震數(shù)據(jù)d，最大迭代次數(shù)K；

(2)初始化模型m0，當前迭代次數(shù)k=0，迭代步長t>0；

(3)若k≤K，迭代更新

mk+1=Tα{mk+2tF-1{(LHL)-1×

LH[F(d)-LF(mk)]}}

k+1→k

(4)輸出 Radon域地震數(shù)據(jù)mk。

算法中，一般要求迭代步長t∈(0,1/λmax(LHL))，其中λmax(·)表示求最大特征值，Tα為近似算子，定義為

(7)

(8)

相較于一些傳統(tǒng)方法，SRTIS提高了Radon模型的稀疏性，從而提高了多次波的壓制精度，同時還有更快的收斂速度。

2.2 基于快速迭代收縮閾值的Radon變換算法

ISTA每次迭代只需考慮當前點的信息更新迭代起點，導致算法迭代速度較慢。FISTA在ISTA的基礎上考慮當前點和前一點更新迭代起點，比ISTA具有更快的收斂速度。Li等[26]將該算法引入到Radon變換，提出了基于快速迭代收縮的混合域快速稀疏時不變Radon變換(an accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on fast iterative shrinkage-thresholding algorithm，SRTFIS)。SRTFIS算法流程如下。

(1)輸入時空域地震數(shù)據(jù)d，最大迭代次數(shù)K；

(2)初始化模型m0，當前迭代次數(shù)k=0，迭代步長t>0，p0=1,y=m0；

(3)若k≤K，迭代更新

(4)輸出Radon域地震數(shù)據(jù)mk。

算法中一般要求t=1/λmax(LHL)。

與SRTIS相比，SRTFIS使用yk代替mk，即在迭代步驟中，SRTFIS算法的迭代點mk不僅僅依賴于前一個迭代點mk-1，還在yk上使用了前兩點{mk-1,mk-2}的線性組合，提高了收斂速度。慣性參數(shù)ak用于控制mk-mk-1的增長速度。

2.3 基于貪婪—快速迭代收縮閾值的Radon變換算法

ISTA計算簡單，但在收斂速度上較慢；FISTA的收斂速度雖然快于ISTA，但對于FISTA，當ak→1時，會引起振蕩現(xiàn)象，對收斂速度造成一定的減緩。海量地震數(shù)據(jù)的處理，對算法的處理效率提出了更高的要求。Liang等[24]提出的Greedy FISTA在保持算法計算簡單的同時，全局收斂速度比ISTA、FISTA更快，有更好的收斂性能，本文將該算法應用于稀疏Radon變換算法中，提出一種基于貪婪的快速迭代收縮的混合域快速稀疏時不變Radon變換(an accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on greedy fast iterative shrinkage-thresholding algorithm，SRTGFIS)。SRTGFIS算法流程如下。

(1)輸入時空域地震數(shù)據(jù)d，最大迭代次數(shù)K。

(2)初始化模型m0，當前迭代次數(shù)k=0，迭代步長t>0,y1=m0,步長t的收縮因子ξ<1,算法收斂的控制因子S>1。

(3)若k≤K，迭代更新

yk=(mk-mk-1)

mk+1=Tα{yk+2tF-1{(LHL)-1×

LH[F(d)-LF(yk)]}}

k+1→k

重啟條件：若(yk-mk+1)T(mk+1-mk)≥0，則

yk=mk；

收斂條件：若‖mk+1-mk‖≥S‖m1-m0‖，則

t=max{ξt,t}

k+1→k

(4)輸出Radon域地震數(shù)據(jù)mk。

該算法中，一般要求t∈[1/λmax(LHL),2/λmax(LHL)]。

為了解決FISTA中ak→1引起的振蕩問題，Greedy FISTA對此引入了重啟動方法，將pk重置為1，迫使ak從0開始增加，當ak→1引起下一次振蕩時，進行重啟，如此循環(huán)。為了縮短重啟動的間隔，減小振蕩周期，以此提高收斂速率，將ak取為常數(shù)1，但由于ak控制迭代步長t的大小，ak為常數(shù)將導致初始迭代步長過大，容易造成算法發(fā)散，因此Greedy FISTA算法還有一個保證收斂的條件。

整個重啟動框架在保證收斂的情況下，縮短了重啟間隔和振蕩周期，相較于ISTA、FISTA算法，該算法提高了收斂速度。

3 實驗數(shù)據(jù)

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

本文實驗的模擬數(shù)據(jù)來自SeismicLab工具包合成的多次波模擬數(shù)據(jù)(圖1)。該模擬數(shù)據(jù)只有一個地震道集，共49道，每道有1001個樣點，采樣間隔為4ms。分別將頻率域最小二乘法Radon變換(LSRT)、SRTIS、SRTFIS、SRTGFIS四種算法用于模擬數(shù)據(jù)進行比較(圖2)。

圖1 合成的多次波全波場模擬數(shù)據(jù)

比較四種算法對模擬數(shù)據(jù)多次波的壓制效果(圖2)，其中SRTIS、SRTFIS和SRTGFIS算法的K均為20。由圖可見，LSRT的壓制效果一般，處理之后還有多次波的殘余(圖2a紅色箭頭所指)，而SRTIS(圖2b)、SRTFIS(圖2c)和SRTGFIS(圖2d)均能完全壓制多次波。

圖2 不同算法對模擬數(shù)據(jù)的多次波壓制效果對比

圖3 三種算法的收斂速度比較

在多次波壓制精度相當?shù)那闆r下，比較三種算法不同迭代次數(shù)的多次波壓制結果(圖4)，對比三種算法對模擬數(shù)據(jù)的計算速度(表1)。圖4a為SRTIS在不同迭代次數(shù)下壓制多次波的結果，當?shù)降?次時，多次波殘余明顯(紅色箭頭所示，下同)，迭代到第16次時，多次波基本被壓制，但依然有殘余；圖4b為SRTFIS在不同迭代次數(shù)下壓制多次波的結果，迭代到第7次時，多次波有少量殘余，迭代到第8次，多次波基本被壓制；圖4c為SRTGFIS在不同迭代次數(shù)下壓制多次波的結果，迭代到第6次時，壓制效果好于SRTFIS，迭代到第7次時，多次波已基本被壓制。

圖4 不同算法不同迭代次數(shù)多次波壓制結果對比

同時，在相同的計算環(huán)境和實驗數(shù)據(jù)下，對程序運行時間進行了測試。當達到相同的壓制精度時，SRTIS、SRTFIS和SRTGFIS三種算法耗時(由各算法分別運行10次求得的平均值)見表1，可見，SRTGFIS較其他兩種算法計算效率更高。

表1 模擬數(shù)據(jù)三種算法達到相同壓制精度所需迭代次數(shù)和時間對比

3.2 實際數(shù)據(jù)

實際數(shù)據(jù)來自SeismicLab工具包，只有1個地震道集，包含92道，每道401個樣點，采樣間隔為4ms。

圖5為實際包含多次波與一次波的全波場數(shù)據(jù)。對于SRTIS、SRTFIS和SRTGFIS算法，當t=1/λmax(LHL)時，步長太小，多次波的壓制效率均太低?？紤]模擬數(shù)據(jù)的測試情況，對于實際數(shù)據(jù)，將步長設置為0.10時，SRTIS、SRTFIS和SRTGFIS三種算法分別需要迭代20、14、10次才能達到較好的多次波壓制效果，故初始化Radon模型相關參數(shù)為：m0=0，t=0.10，K=20。LSRT、SRTIS、SRTFIS、SRTGFIS四種算法多次波的壓制結果如圖6所示。LSRT算法的壓制結果中還有殘留的多次波(圖6a紅色箭頭所指)；SRTIS(圖6b)、SRTFIS(圖6c)和SRTGFIS(圖4d)算法的多次波壓制效果均優(yōu)于LSRT算法。

圖5 實際數(shù)據(jù)

圖6 四種方法對實際數(shù)據(jù)的多次波壓制結果對比

實際數(shù)據(jù)三種算法的收斂曲線(圖7)顯示，SRTIS的收斂速度最慢，SRTGFIS算法的收斂速度最快，由于ak=1導致初始迭代步長過大，引起收斂曲線的波動。

圖7 三種算法的收斂速度比較

在多次波壓制精度相當?shù)那闆r下，分別比較三種算法不同迭代次數(shù)下的多次波壓制效果(圖8～圖10)。SRTIS算法在迭代到第18次時，多次波有少量殘余(圖8中紅色箭頭所指，下同)，迭代到第20次時，多次波得到基本壓制；SRTFIS算法在迭代到第12次時，多次波有少量殘余，迭代到第14次時，多次波基本壓制(圖9)；SRTGFIS算法在迭代到第8次時，多次波有少量殘余，迭代到第10次時，多次波得到很好壓制(圖10)。

圖8 SRTIS 法不同迭代次數(shù)下的多次波壓制結果對比

圖9 SRTFIS 法不同迭代次數(shù)下的多次波壓制結果對比

圖10 SRTGFIS 法不同迭代次數(shù)下的多次波壓制結果對比

在相同的計算環(huán)境和實驗數(shù)據(jù)下，對程序運行所需時間進行了測試。當達到相同的壓制精度時，SRTIS、SRTFIS及SRTGFIS算法耗時如表2所示(各算法分別運行10次求得的平均值)，SRTGFIS算法優(yōu)于前兩種，計算效率更高。

表2 實際數(shù)據(jù)三種算法達到相同壓制精度所需迭代次數(shù)和時間對比

4 結論和認識

本文將Greedy FISTA算法引入到Radon變換壓制多次波的逆問題求解中，構建了SRTGFIS算法。模擬數(shù)據(jù)和實際數(shù)據(jù)的處理結果均表明：該算法對多次波的壓制效果優(yōu)于LSRT算法；與SRTIS和SRTFIS算法相比，當壓制精度相同時，效率更高。

隨著深度學習方法的不斷發(fā)展，在其他學科領域均已展現(xiàn)出了巨大的優(yōu)勢，將其引入到多次波壓制，避免傳統(tǒng)多次波壓制算法中逆問題求解的時間消耗，是下一步提高多次波壓制效率的研究方向。