張 全 雷 芩 林柏櫟 彭 博*② 劉書(shū)妍②

(①西南石油大學(xué)計(jì)算機(jī)科學(xué)學(xué)院,四川成都 610500；②油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室(西南石油大學(xué))，四川成都 610500；③電子科技大學(xué)信息與通信工程學(xué)院，四川成都 611731)

0 引言

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

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

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

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

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

d=Lm

(1)

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

(2)

L通常是病態(tài)的，且LS法不適用于求解病態(tài)方程，常添加正則化項(xiàng)解決此問(wèn)題。基于L1范數(shù)的稀疏正則化是目前常用的提高Radon變換分辨率的方法

(3)

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

d=F-1[LF(m)]

(4)

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

(5)

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

Radon變換壓制多次波在混合域的目標(biāo)方程為

(6)

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

2 Radon變換算法及其改進(jìn)

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

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

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

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

(2)初始化模型m0，當(dāng)前迭代次數(shù)k=0，迭代步長(zhǎng)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。

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

(7)

(8)

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

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

ISTA每次迭代只需考慮當(dāng)前點(diǎn)的信息更新迭代起點(diǎn)，導(dǎo)致算法迭代速度較慢。FISTA在ISTA的基礎(chǔ)上考慮當(dāng)前點(diǎn)和前一點(diǎn)更新迭代起點(diǎn)，比ISTA具有更快的收斂速度。Li等[26]將該算法引入到Radon變換，提出了基于快速迭代收縮的混合域快速稀疏時(shí)不變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í)空域地震數(shù)據(jù)d，最大迭代次數(shù)K；

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

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

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

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

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

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

ISTA計(jì)算簡(jiǎn)單，但在收斂速度上較慢；FISTA的收斂速度雖然快于ISTA，但對(duì)于FISTA，當(dāng)ak→1時(shí)，會(huì)引起振蕩現(xiàn)象，對(duì)收斂速度造成一定的減緩。海量地震數(shù)據(jù)的處理，對(duì)算法的處理效率提出了更高的要求。Liang等[24]提出的Greedy FISTA在保持算法計(jì)算簡(jiǎn)單的同時(shí)，全局收斂速度比ISTA、FISTA更快，有更好的收斂性能，本文將該算法應(yīng)用于稀疏Radon變換算法中，提出一種基于貪婪的快速迭代收縮的混合域快速稀疏時(shí)不變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í)空域地震數(shù)據(jù)d，最大迭代次數(shù)K。

(2)初始化模型m0，當(dāng)前迭代次數(shù)k=0，迭代步長(zhǎng)t>0,y1=m0,步長(zhǎng)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引起的振蕩問(wèn)題，Greedy FISTA對(duì)此引入了重啟動(dòng)方法，將pk重置為1，迫使ak從0開(kāi)始增加，當(dāng)ak→1引起下一次振蕩時(shí)，進(jìn)行重啟，如此循環(huán)。為了縮短重啟動(dòng)的間隔，減小振蕩周期，以此提高收斂速率，將ak取為常數(shù)1，但由于ak控制迭代步長(zhǎng)t的大小，ak為常數(shù)將導(dǎo)致初始迭代步長(zhǎng)過(guò)大，容易造成算法發(fā)散，因此Greedy FISTA算法還有一個(gè)保證收斂的條件。

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

3 實(shí)驗(yàn)數(shù)據(jù)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論和認(rèn)識(shí)

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

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