應用加權迭代軟閾值算法的高分辨率Radon變換

薛亞茹 郭蒙軍* 馮璐瑜 馬繼濤 陳小宏

(①中國石油大學(北京)信息科學與工程學院，北京 102249； ②中國石油大學(北京)地球物理學院，北京 102249)

0 引言

奧地利數學家Radon于上世紀初提出了Radon變換理論，經過不斷的發展完善，逐漸從數學領域沿用到其他領域。Claerbout等[1]將Radon變換引入地震勘探領域，極大地促進了Radon變換在地震資料處理方面的應用，如多次波壓制、缺失地震道重建、平面波分解以及去噪，且該變換可在時域、頻域及混合域實現。

由于線性Radon變換和拋物Radon變換具有時不變特性，故可在頻域快速求解。Hampson[2]提出拋物Radon變換并將其應用于多次波壓制，且在頻域用最小二乘(Least square，LS)算法求解，但該方法的分辨率并不高； 之后，Sacchi等[3-4]結合貝葉斯原理，通過引入模型的先驗信息提高Radon變換的分辨率。由于高分辨率Radon變換迭代過程涉及矩陣求逆運算，導致其計算量較大，為此采用低頻約束方法，即用上一頻率計算結果約束當前頻率的計算[5-6]。劉喜武等[7]利用稀疏約束共軛梯度算法求解Radon變換，與阻尼LS算法相比，提高了Radon變換的分辨率和計算效率； 王維紅等[8]提出基于Levinson遞推法的加權拋物Radon變換疊前地震數據重建方法，與傳統的高分辨率Radon變換方法相比，計算效率有所提升； 薛亞茹等[9]將傳統拋物Radon變換與正交多項式相結合，克服了空間假頻，并保留地震波的AVO特性； 之后，王亮亮等[10]在此基礎上引入新變量，消除了變換算子對頻率的依賴，提高了計算效率。由于L1范數不具有真正意義上的稀疏性，薛亞茹等[11]引入SL0范數約束，并通過最速下降法和梯度投影原理實施了高分辨率Radon變換。

由于Radon模型在時域比在頻域更具有稀疏性，Cary[12]和Schonewille等[13]在時域實現了Radon變換，提高了其分辨率，但收斂速度較慢； 鞏向博等[14]提出各向異性Radon變換，通過最優相似系數加權Gauss-Seidel迭代算法，提高了Radon變換計算效率。結合Radon變換在時域的稀疏性及在頻域的計算高效性，Trad等[15]在頻域實現Radon變換的正向和反向過程，并利用迭代加權最小二乘法(Iterative reweighted least square，IRLS)在時域求解稀疏Radon變換，充分兼顧時域Radon變換和頻域Radon變換的優點； 此后，Lu[16]用迭代收縮算法(Iterative-shrinkage algorithms)代替IRLS，進一步提高混合域Radon變換的計算效率； 在此基礎上，Wang等[17]針對噪聲偏離高斯分布的Radon變換，引入Bregman方法，提高了混合域Radon變換的魯棒性； 鞏向博等[18]在混合域實現雙曲Radon變換，并引入快速迭代軟閾值算法(Fast iterative shrinkage-thresholding algorithm，FISTA)加快反演收斂速度。隨著Radon變換在地震數據領域的更廣泛應用，近年來還研發出基于字典學習的用于去噪和插值的Radon變換，通過稀疏表征地震數據并結合K-SVD方法實現求解[19-21]。由于高分辨率Radon變換計算量較大且計算結果易受正則化影響，也有學者提出在壓縮感知框架下處理地震數據，利用匹配追蹤等算法，以提高地震數據處理的計算精度和效率[22-26]。

本文分析了迭代軟閾值法(Iterative soft threshold algorithm，ISTA)和IRLS實現Radon變換的基本原理，提出將IRLS的加權矩陣思想引入ISTA方法中，形成加權迭代軟閾值法(Reweighted-Iterative soft threshold algorithm，R-ISTA)，以進一步提高Radon變換分辨率。

1 方法原理

1.1 Radon變換反演基本原理

在地震勘探領域，常用的Radon變換有線性Radon變換、拋物Radon變換和雙曲Radon變換。拋物Radon變換將動校正后的同相軸沿時距曲線進行疊加，得到Radon參數。由于采集到的地震數據為離散數據，一般采用離散Radon變換，Radon正變換離散形式為

(1)

式中：m為Radon域數據；d為時空域地震數據；q為曲率參數；h為炮檢距；t為時間；τ為時空域雙重旅行時截距；Nq為Radon域曲率參數個數。

拋物Radon變換具有時不變性，可將其變換到頻域進行計算，以提高計算效率。對式(1)做傅里葉變換得到對應的頻域Radon正變換公式

(2)

可寫成如下矩陣形式

D=LM

(3)

其中Radon變換算子矩陣L的具體表達式為

(4)

式中Nh為道數。

由于矩陣L通常不是方陣，無法求逆，其LS解分辨率低，基于L1范數的稀疏正則化是現今常用的高分辨率反演方法。構造高分辨率反演目標函數為

(5)

式中λ為拉格朗日因子，且有λ>0，用于平衡模型的準確度和稀疏性。ISTA是上述優化問題常用求解方法之一[27-30]，該算法求解過程不涉及矩陣求逆，計算效率高，每次迭代過程中經矩陣乘法后通過軟閾值函數逐漸逼近所求變量，其迭代式為

Mn+1=Sσ[Mn+ηLH(D-LMn)]

(6)

式中：n為迭代次數；η為正參數，通常取為LHL最大特征值的倒數；Sσ為軟閾值函數，其表達式為

Sσ(M,σ)=sign(M)max{0,|M|-σ}

(7)

式中： sign為符號函數；σ為閾值，該值與上一次迭代結果有關，即達到自適應閾值的目的，其經驗表達式為

σ=0.01 max(|Mn|)

(8)

由于矩陣L具有非正交性，式(6)中共軛算子LH降低了Radon變換的分辨率，使得上述算法應用于Radon變換反演時收斂速度較慢。

1.2 基于R-ISTA算法的高分辨Radon變換

IRLS是高分辨Radon變換的另一種常用方法[31]。Sacchi等[3-4]基于貝葉斯原理，將模型先驗信息與Radon變換相結合，將前一次迭代結果作為下一次反演的加權矩陣； 經過若干次迭代之后，得到高分辨率的Radon變換。

據式(3)構造如下目標函數

(9)

式中：μ為阻尼因子，取值范圍是0.01～1.00；WM為數據M的協方差矩陣。將式(9)最小化，可得

M=(LHL+μW)-1LHD

(10)

為克服ISTA中Radon共軛算子分辨率低的問題，引入IRLS的Radon變換中權重矩陣思想，將IRLS的加權反演矩陣引入ISTA中，得到R-ISTA公式，即改進的式(6)表達式

Mn+1=Sσ{Mn+ηB-1[LH(D-LMn)]}

(11)

與式(6)所示傳統的ISTA公式相比，本文引入與加權矩陣相關的反演矩陣

B=LHL+μW

(12)

式中對角矩陣W的對角元素與前一次迭代Radon域數據M有關，用來聚焦Radon域能量，即

(13)

式中b是穩定因子。

與傳統ISTA方法相比，本文所提R-ISTA方法通過引入權重矩陣使Radon域能量聚焦，當上一次迭代得到的Radon域數據M部分能量較低時，則相應的權重矩陣較大。類似地，當前一次迭代得到的M數據部分能量較高時，其對應的權重矩陣較小； 隨后經矩陣求逆運算，使得數據M中能量強的部分繼續加強，能量弱的部分進一步減弱，以此達到高分辨率Radon變換的目的。

1.3 基于主頻約束的R-ISTA高分辨Radon變換

迭代加權方法對每個頻率都需計算反演權重矩陣(式(12))，還需做求逆運算，故計算量較大。為了降低運算耗時，可用單個頻率的運行結果約束其他頻率，避免逐個分別對每一頻率計算權重矩陣，這樣就可提高Radon變換計算效率[5-6，32]。

為提高R-ISTA的計算效率，本文設計了適用于多次波壓制的主頻約束R-ISTA。首先選取地震數據的主頻進行迭代訓練，在分辨率達到要求后停止迭代，并保存其迭代得到的最終權重矩陣，然后用該頻率訓練得到的權重矩陣去約束其他所有頻率的反演。基于主頻約束的R-ISTA迭代公式為

(14)

式中Bmain為定值，不隨迭代次數變化，其表達式為

Bmain=LHL+μWmain

(15)

式中Wmain為主頻訓練迭代得到的加權矩陣，此后其值不再隨頻率變化。由于只需進行一次求逆運算，所以計算量大幅度減少。

基于主頻約束的R-ISTA壓制多次波的步驟如下：

(1)設置初始參數，通過傅里葉變換把時空域地震數據d變換到頻域D=F(d)；

(2)選擇主頻對應的頻域數據Dmain，利用式(11)獲得該主頻的Radon參數，保存其加權參數矩陣Wii=(|Mi|2+b2)-1，以此約束其他頻率的反演；

(3)對頻域所有地震數據D，由LS計算其迭代初始值M0=(LHL+μI)-1LHD； 將M0代入主頻約束R-ISTA迭代式(式(14))，其中偽逆矩陣由式(15)計算，迭代至收斂，得到Radon域數據Mn；

(4)在Radon域分離一次波與多次波，實現多次波壓制。

2 數據實驗

2.1 模擬數據

為比較ISTA、IRLS及R-ISTA三種方法實現Radon變換的分辨率及壓制多次波的效果，構建如圖1a所示的模擬地震數據：采用主頻為30Hz的Ricker子波，采樣點數為200，采樣間隔為4ms，共計64道； 共含4條同相軸，其中一次波(圖1b)和多次波(圖1c)各兩條。針對圖1a模擬數據對應的原始Radon數據(圖2)，分別用ISTA、主頻約束IRLS及主頻約束R-ISTA三種方法進行處理，反演得到對應的Radon域數據(圖3a～圖3c)及其與原始Radon域數據的誤差(圖3d～圖3f)。其中：ISTA和主頻約束R-ISTA的迭代次數為10； 主頻約束IRLS的主頻迭代也為10次，其他頻率反演無需迭代。

圖1 模擬地震數據

圖2 原始Radon域數據

觀察、對比可見：ISTA反演得到的Radon域數據能量較分散，分辨率較低，無法分離剩余時差差異較小的兩同相軸； 主頻約束IRLS雖可將其分離，但還是存在能量擴散現象； 而主頻約束R-ISTA反演得到的Radon參數能量集中，分辨率高。

進一步分析反演所得結果(圖3)，可知一次波與多次波已分離，去除其中一次波，經Radon反變換可得到時域多次波數據，從原始地震數據減去多次波數據，即可達到去除多次波的效果。分別用ISTA、主頻約束IRLS和主頻約束R-ISTA處理模擬地震數據后，恢復出圖4a～圖4c所示的多次波數據，用原始多次波數據(圖1c)分別減去上述三種方法恢復的多次波數據，得到圖4d～圖4f所示的(放大了5倍的)多次波殘差數據。從原始地震數據(圖1a)減去三種方法恢復的多次波數據(圖4a～圖4c)，可得到三種方法估計的一次波數據(圖5a～圖5c)及其與原始一次波(圖1b)的殘差(圖5e～圖5f，放大了5倍)。

圖3 不同方法反演得到的Radon域數據(上)及其與原始數據的差值(下)

從圖4d～圖4f可觀察到，主頻約束IRLS和主頻約束R-ISTA恢復的多次波更接近于真實值，ISTA恢復的多次波與真實值相差較大。觀察圖5e～圖5f，可發現在一次波同相軸與多次波同相軸相鄰位置，ISTA方法的多次波殘留明顯，且一次波數據有丟失； 主頻約束IRLS和主頻約束R-ISTA兩種方法對多次波壓制效果都優于ISTA，但主頻約束IRLS方法(圖5e)仍殘存明顯一次波數據； 而圖5f中則基本看不到一次波殘留，因此主頻約束R-ISTA恢復的一次波數據更接近實際一次波數據。

圖4 不同方法恢復的多次波數據(上)及其與原始多次波數據的差值(下，放大5倍)

圖5 不同方法恢復的一次波數據(上)及其與原始一次波數據的差值(下，放大5倍)

用信噪比衡量不同方法對多次波的壓制效果，信噪比定義為

(16)

式中：d為原始一次波數據；dm為不同方法各自得到的一次波數據。本次采用的三種方法所用數據的信噪比如表1所示，可見主頻約束R-ISTA得到的信噪比為31.0404dB，遠大于ISTA的7.7926dB，同樣高于主頻約束IRLS的21.6724dB。因此，主頻約束R-ISTA對多次波的壓制效果優于另兩種方法。

為評估R-ISTA的抗噪性能，針對圖1a模擬數據加入高斯白噪聲，得到圖6所示含噪地震數據。根據信噪比定義式(式(16))算出該數據的信噪比為 0dB。分別用ISTA、主頻約束IRLS及主頻約束R-ISTA壓制多次波，得到一次波(圖7a～圖7c)及其誤差數據(圖7d～圖7f)，三種方法去噪后信噪比見表1。從圖7及表1可見，主頻約束R-ISTA的去噪效果略優于主頻約束IRLS，且都遠優于ISTA。

圖7 針對含噪地震數據用不同方法恢復的一次波數據(上)及其與原始一次波數據的差值(下，放大5倍)

表1 加噪前后三種方法得到一次波的信噪比

圖6 含高斯白噪模擬地震數據

2.2 實際數據

選取M地區三維實際地震數據(圖8a)，共計13炮，每炮25道，采樣點數為2501，采樣間隔為2ms。圖8b為對原始數據按炮點距離排序后的實際數據。 經動校正后，一次波同相軸基本被校平，曲率較小， 而多次波仍有剩余時差，曲率較大，故可根據曲率的不同分離一次波與多次波。從圖8b可見多次波主要存在于圖像中部(紅色橢圓)區域，為便于對多次波壓制過程做數據分析，選擇僅展示圖8c所示的中間5炮數據。該數據Radon變換曲率參數q取值范圍是-0.4～0.6s。圖8d～圖8f分別為ISTA、主頻約束IRLS和主頻約束R-ISTA三種方法反演的Radon域數據，切除0.2s以下部分(黑線左側區域)，即可獲取多次波對應的Radon域數據； 再經Radon反變換得到多次波數據，用原始數據減去多次波數據即可得三種方法壓制多次波后結果(圖8g～圖8i)。

由于從圖8g～圖8i難以判明三種方法優劣性，故用上述三種方法處理13炮數據后，再按炮點距離排序，得到圖9a～圖9c所示多次波數據； 進而獲得三種方法去除多次波后的數據(圖9d～圖9f)。其中：ISTA迭代20次； 主頻約束R-ISTA迭代5次； 主頻約束IRLS的主頻迭代10次，其他頻率反演無需迭代。觀察圖8d～圖8f中三種方法反演的Radon域數據，可見主頻約束R-ISTA求取的Radon域數據比另兩種方法更稀疏，去除多次波后可保留更多一次波信息； 從圖9d～圖9f中箭頭所指部分可見，與原始數據相比，經三種方法處理后，ISTA對多次波的壓制效果欠佳，另兩種方法的壓制多次波效果較理想，且主頻約束R-ISTA比主頻約束IRLS的一次波同相軸更清晰。

圖8 實際地震數據

圖9 按炮點距離排序地震數據不同方法去除的多次波(上)及剩余數據(下)

2.3 收斂速度

已知ISTA的收斂速度為O1(1/k， 其中k為迭代次數，下同)，即為次線性收斂； FISTA的收斂速度為O2(1/k2)[28]； IRLS為指數收斂[33]。本文以均方誤差為收斂性能的評判指標，其表達式為

(17)

將R-ISTA與收斂速度已知的ISTA、FISTA、IRLS及基礎的LS對比，得到圖10所示的收斂速度圖。可見FISTA收斂速度快于ISTA，且此兩者的精度略高于LS。迭代初始，IRLS和R-ISTA的收斂速度相差不大，且都快于ISTA； 迭代12次之后，IRLS趨于收斂，R-ISTA的均方誤差則繼續下降，迭代20次趨于收斂； 最終，R-ISTA的均方誤差小于另四種算法。

圖10 ISTA、FISTA、LS、IRLS、R-ISTA收斂速度

從表2可見，對于相同迭代次數(1次迭代除外)，ISTA的運行速度明顯快于R-ISTA和IRLS，這是由于ISTA不需進行矩陣求逆運算，而IRLS和R-ISTA的每次迭代過程，都需計算偽逆矩陣B=LHL+μW，當待處理數據較復雜時，其計算量更大。本次計算環境為Intel Core i5，主頻為2.3 GHz ，采用Windows平臺下的Matlab實現。

表2 三種算法的運算時間對比

為減小R-ISTA計算量，本文在R-ISTA基礎上引入主頻約束，此時偽逆矩陣Bmain=LHL+μWmain只需求解一次，計算量顯著減少，主頻約束對運行速度的影響見表3。可見當只進行一次迭代時，主頻約束對運行時間影響不大； 隨著迭代次數的增加，對無主頻約束的R-ISTA而言，偽逆矩陣要隨之計算更新，導致其計算量顯著增大； 而有主頻約束的R-ISTA由于只需求解一次偽逆矩陣，之后不再隨迭代次數增加而更新，且軟閾值函數的相關計算量較小，當處理大量數據時，其計算量增加緩慢。

表3 有/無主頻約束R-ISTA方法的迭代時間

3 結束語

本文將高分辨Radon變換中加權矩陣的思想引入ISTA方法，利用Radon參數的先驗信息約束反演誤差函數，提出R-ISTA方法，解決了傳統軟閾值算法用于Radon變換表現的分辨率低及收斂速度慢問題。模擬數據和實際資料處理結果表明，與ISTA和IRLS相比，所提R-ISTA的Radon變換的分辨率有所提升，且具有更強壓制多次波能力。

由于R-ISTA方法對于每一個頻率仍需計算加權矩陣并進行矩陣求逆運算，計算量較大，因此本文在壓制多次波時，引入主頻約束思想，提高了R-ISTA方法對大型數據的處理能力。