基于SL0的高分辨率Radon變換及數(shù)據(jù)重建

薛亞茹 王 敏 陳小宏

(中國石油大學(xué)(北京)地球物理與信息工程學(xué)院，北京 102249)

1 引言

由于地震數(shù)據(jù)在Radon變換域的物理特征直觀明確，利于對比分析，因此Radon變換在勘探地震學(xué)領(lǐng)域，如平面波分解[1]、去噪[2]、數(shù)據(jù)恢復(fù)[3]、多次波壓制[4-6]等方面得到廣泛應(yīng)用。Radon變換算子的非正交性導(dǎo)致直接進行Radon正、反變換時能量的不對等性，于是基于最小范數(shù)反演思想的Radon變換被提出。該方法雖然減少了拖尾現(xiàn)象，但會產(chǎn)生平滑效應(yīng)，不能集中足夠能量，在Radon域分辨率較低[7]。因此，要想提高Radon變換的分辨率，則必須對反演稀疏約束的方法進行改進。Sacchi等[8，9]提出在頻率域通過柯西稀疏約束提高Radon變換的分辨率，主要通過數(shù)據(jù)的先驗信息修正反演算子，使變換后的能量在Radon域得到收斂，這樣不僅有效地消除了空間褶積(或有限孔徑)的影響，同時也進一步提高了解的穩(wěn)定性。高分辨率Radon變換法是一種稀疏約束反演算法[10]，目前常用的最小范數(shù)約束方法是在頻率域通過柯西或L1范數(shù)稀疏約束提高Radon變換的分辨率，但作為L0范數(shù)的凸逼近形式，L1范數(shù)不具有真正意義上的稀疏，它是一種凸松弛算法[11，12]。L0范數(shù)是度量數(shù)據(jù)稀疏度的最好方法，所以為了達到真正的稀疏，要尋求一種嚴格的L0范數(shù)來約束Radon變換。

因求解L0范數(shù)計算量巨大，故在實際應(yīng)用中不斷推出新的稀疏優(yōu)化方法[13-19]。Mohimani等[13，14]提出平滑L0范數(shù)(Smoothed L0Norm，SL0)算法，這是一種基于過完備稀疏分解的快速算法，能直接實現(xiàn)L0范數(shù)最小化。該算法具有重構(gòu)前不需知道信號的稀疏度、計算量小、高匹配度以及重建時間短等優(yōu)點[20]，可有效地應(yīng)用于稀疏信號重構(gòu)。基于SL0范數(shù)算法，Zayyani等[16]、Ghalehjegh等[17]和林婉娟等[18]相繼提出TSL0(Thresholded SL0)算法、BSL0(Block SL0)算法和NSL0算法。曹靜杰等[21]還提出在Curvelet域的零范數(shù)稀疏最優(yōu)化方法，很好地實現(xiàn)了地震數(shù)據(jù)的重構(gòu)，明顯提高了計算效率。

本文將SL0范數(shù)稀疏約束引入Radon變換，可提高地震數(shù)據(jù)在Radon域的能量聚焦效果。利用該方法實現(xiàn)缺失地震道重構(gòu)[22]，通過理論模型和實際數(shù)據(jù)實驗與柯西稀疏約束的高分辨率Radon變換方法作比較，證明L0范數(shù)稀疏約束的Radon變換更具有優(yōu)越性。

2 基本原理及實現(xiàn)

2.1 SL0范數(shù)定義

一個向量的L0范數(shù)是這個向量的不連續(xù)函數(shù)，直接使用L0范數(shù)存在很多問題。為此，通過構(gòu)造一個平滑的連續(xù)函數(shù)逼近這個不連續(xù)的函數(shù)(即SL0范數(shù))。

通常選用標準高斯函數(shù)

(1)

式中參數(shù)σ決定該函數(shù)逼近效果。定義向量函數(shù)

(2)

式中K是向量維數(shù)。由式(1)和式(2)可知，當σ取值較小時，則向量s的L0范數(shù)近似為‖s‖0≈K-Fσ(s)。因此，可通過下式來優(yōu)化L0范數(shù)解

maxFσ(s) s.t.As=x

(3)

滿足上式的參數(shù)Fσ要取最小值，當σ>0時上述連續(xù)函數(shù)與L0范數(shù)最接近。所以在σ取最小值的情況下通過Fσ(s)最大化得到SL0范數(shù)最小的解。參數(shù)σ的值決定了函數(shù)Fσ的光滑程度，σ值越大，F(xiàn)σ越光滑，且局部極值較少，可很容易地實現(xiàn)它的最大化，但缺點是不能很好地逼近最優(yōu)解； 相反地，σ值越小，逼近最優(yōu)解的效果越好，但是σ取值太小，F(xiàn)σ高度不平滑，導(dǎo)致包含大量局部最大值，很難實現(xiàn)Fσ的最大化。因此提出在算法中動態(tài)調(diào)整參數(shù)σ的策略。選取一個σ的遞減序列{σ1，σ2，…，σJ}(其中J為遞減序列元素個數(shù))，對每一個σ的值都用最速下降法不斷逼近目標函數(shù)，求得Fσ的最大值。由于σ的變化幅度不是很大，每次計算得到Fσ的最優(yōu)解都與上一次得到的最優(yōu)解很接近，這樣就能避免陷入局部最大值問題，從而得到σ值很小時Fσ的最大值，最終求得使SL0范數(shù)最小的解。

2.2 SL0范數(shù)約束的Radon變換

拋物Radon變換在時域定義為

(4)

式中：d(t，x)為時域地震數(shù)據(jù)；x為炮檢距；m(τ，q)為Radon域數(shù)據(jù)；q為曲率參數(shù)；τ為截距時間。通過傅里葉變換將式(4)變到頻率域，其表達式為

(5)

式中：M(ω，q)和D(ω，x)分別是m(τ，q)和d(t，x)的傅里葉變換；ω為頻率；N為q的取值個數(shù)。寫成矩陣形式為

d=Lm

(6)

式中L為Radon算子。

由上式可知，Radon變換其實就是求解方程組d=Lm的過程，可采用L1范數(shù)的優(yōu)化問題求該方程組的稀疏解。L1范數(shù)稀疏方法雖收斂性較好，但只是L0范數(shù)的凸逼近形式，L0范數(shù)才是最能體現(xiàn)數(shù)據(jù)稀疏的方法。故將Radon變換求稀疏解問題轉(zhuǎn)化為基于SL0范數(shù)的稀疏優(yōu)化問題，其數(shù)學(xué)模型為

min‖m‖0 s.t.Lm=d

(7)

根據(jù)SL0范數(shù)的基本定義，選取高斯函數(shù)作為平滑連續(xù)函數(shù)來逼近L0范數(shù)，可得

‖m‖0≈N-Fσ(m)

對式(7)的優(yōu)化問題可采用最速下降法和梯度投影原理求解，其中最速下降法是由迭代形式

構(gòu)成的。隨著σ的減小，μk值也應(yīng)減小。這是因為σ越小，函數(shù)Fσ的“波動”越大，因此最大化Fσ要用小步長[20]。通過改變σ值觀察不同尺度下的相同曲線，其中尺度是與σ成比例的，不難發(fā)現(xiàn)μk與σ2是成正比的。因此令μk=μσ2，其中μ是一個常數(shù)，其值通過經(jīng)驗設(shè)定，由此得到

式中

為梯度方向，它是通過對函數(shù)

求導(dǎo)得到。利用該梯度方向搜索最優(yōu)值，直至得到Fσ(m)的最優(yōu)解，即L0范數(shù)的最小量。

參數(shù)σ的取值非常重要，關(guān)系到最優(yōu)值的搜索方向，當σ取值足夠大時，滿足

的最優(yōu)解就是Lm=d的最小L2范數(shù)解。下面利用拉格朗日乘數(shù)法驗證該觀點。

設(shè)拉格朗日函數(shù)L(m,λ)=Fσ(m)-λT(Lm-d)，分別對m和λ(拉格朗日乘數(shù))求導(dǎo)，并令其導(dǎo)數(shù)都等于零，可得

(8)

定義式中λ1=-σ2λ。另外,Lm=d的最小L2范數(shù)解可用最小化

求解。利用拉格朗日乘數(shù)法，最小化結(jié)果可表示為

(9)

比較式(8)和式(9)，可以發(fā)現(xiàn)當σ→∞ (或σ?max{m1,…,mN})時，這兩個方程組是相等的。此時使Fσ(m)最大的最優(yōu)解就是Lm=d的最小L2范數(shù)解。所以用Lm=d的最小L2范數(shù)解作為算法的初始值，是優(yōu)化算法的最好選擇。

綜上所述，基于SL0范數(shù)稀疏約束Radon變換的具體算法分如下四步。

1.對某數(shù)據(jù)d，設(shè)置初值m0為最小二乘解；

2.選取一個合適的σ遞減序列[σ1,σ2,…,σJ]；

3.對于每一個σ的取值都做如下處理：

令σ=σj(j=1,2,…,J)，在解空間M={m|Lm=d}中，利用最速下降法迭代C次使函數(shù)Fσ(m)最大化，然后將得到的新的m值投影到解空間M中，具體步驟如下：

(1)令m=mi-1

(2)對于l=1,…,C

2)更新Radon變換的值m←m-μδ

3)采用梯度投影原理可得m←m-(LHL+μI)-1LH(Lm-d)

(3)mi=m

4.經(jīng)過以上迭代循環(huán)，最終得到Radon域的數(shù)據(jù)m=mJ。

在地震數(shù)據(jù)完整的情況下，數(shù)據(jù)經(jīng)過拋物Radon變換后會在Radon域聚焦得很好。但是，當?shù)卣饠?shù)據(jù)缺失、有空道存在時，經(jīng)過拋物Radon變換得到Radon域的數(shù)據(jù)能量比較發(fā)散，利用拋物Radon變換實現(xiàn)缺失地震數(shù)據(jù)重構(gòu)的方法，實際上就是通過正反拋物Radon變換的多次迭代，使Radon域上的能量重新得到收斂，以恢復(fù)時域上的地震同相軸，達到地震數(shù)據(jù)重構(gòu)的目的。

稀疏重構(gòu)問題實際上就是如何求解欠定方程組的最稀疏解。在缺失地震數(shù)據(jù)重構(gòu)中，利用L0范數(shù)稀疏約束Radon變換來求解Radon域的稀疏解，增強了地震同相軸在Radon域的收斂特性[6]，保持了地震同相軸在時域的連續(xù)性，從而提高了地震數(shù)據(jù)重構(gòu)的精度。

3 理論模型與實際數(shù)據(jù)處理

3.1 理論模型

為了證明本文算法的可行性和有效性，以地震數(shù)據(jù)重建為例對比基于L1范數(shù)的高分辨率Radon變換和SL0Radon變換性能。構(gòu)建了兩個理論模型進行實驗。第一個模型如圖1所示，圖1a是由主頻為30Hz的Ricker子波合成的地震記錄，共64道，道間距為10m，每道采樣點數(shù)為200，采樣間隔為4ms，隨著炮檢距的改變3個同相軸的振幅不發(fā)生變化。圖1b為均勻缺失模型，其中每隔兩道抽取一道作為缺失待重建的數(shù)據(jù)。圖1c和圖1d分別為高分辨率Radon變換和SL0范數(shù)約束的Radon變換方法重構(gòu)的結(jié)果，圖1e和圖1f分別為它們所對應(yīng)的重構(gòu)數(shù)據(jù)誤差剖面圖，從圖中可清楚地看到SL0范數(shù)約束的Radon變換方法將三個同相軸恢復(fù)的很好，而高分辨率Radon變換方法重構(gòu)數(shù)據(jù)殘差中仍然有同相軸部分信息殘留下來。

圖2a為模型一的合成地震記錄的理想Radon譜，圖2b為高分辨率Radon變換譜，圖2c為SL0范數(shù)約束的Radon變換譜，圖2d為高分辨率Radon變換譜殘差，圖2e為SL0范數(shù)約束Radon變換譜殘差，通過對比可看出SL0范數(shù)約束的Radon變換方法比高分辨率Radon變換方法的的分辨率更高，聚焦效果更好。

該地震數(shù)據(jù)f-k譜如圖3a所示。圖3b是其均勻缺失地震數(shù)據(jù)f-k譜，與原始數(shù)據(jù)f-k譜(圖3a)相比，發(fā)現(xiàn)數(shù)據(jù)缺失后出現(xiàn)了明顯假頻現(xiàn)象。圖3c和圖3d分別是高分辨率及SL0兩種方法對應(yīng)的f-k譜；圖3e和圖3f是對應(yīng)兩種方法重建數(shù)據(jù)f-k譜分別與真實頻譜之間的誤差，可看出高分辨率Radon變換方法在重構(gòu)過程中有能量丟失，這與圖1e中殘留的同相軸信息相對應(yīng)，而SL0范數(shù)約束的Radon變換方法在保留原始數(shù)據(jù)能量的同時，還壓制了假頻，重構(gòu)效果明顯優(yōu)于高分辨Radon變換。

圖4a是模型二地震合成記錄，其子波主頻為30Hz，共64道，其中道間距為10m，每道采樣點300個，采樣間隔為4ms，三個同相軸的振幅隨著炮檢距的改變而改變。圖4b是圖4a中數(shù)據(jù)隨機缺失26道所得的缺失數(shù)據(jù)。圖4c和圖4d分別為高分辨率Radon變換及SL0范數(shù)約束的Radon變換重構(gòu)結(jié)果，圖4e和圖4f分別為其所對應(yīng)的誤差剖面，通過誤差可看出高分辨率Radon變換方法的重構(gòu)結(jié)果，不管是在近炮檢距還是在遠炮檢距，都有同相軸信息的丟失，而SL0范數(shù)約束的Radon變換方法重構(gòu)結(jié)果更精確，數(shù)據(jù)的AVO特性恢復(fù)較好。

圖1 振幅不變的均勻缺失數(shù)據(jù)重構(gòu)模型 (a)原始地震數(shù)據(jù)； (b)均勻缺失數(shù)據(jù)； (c)高分辨率方法重構(gòu)結(jié)果； (d)SL0 Radon變換方法重構(gòu)結(jié)果； (e)高分辨率方法誤差剖面(放大5倍)； (f)SL0 Radon變換方法誤差剖面(放大5倍)

圖2 振幅不變均勻缺失時兩種Radon變換方法的Radon譜 (a)合成記錄理想Radon譜； (b)高分辨率Radon變換譜； (c)SL0范數(shù)約束Radon 變換譜； (d)高分辨率Radon變換譜誤差； (e)SL0范數(shù)約束Radon變換譜誤差

圖3 振幅不變的均勻缺失數(shù)據(jù)重建f-k譜比較 (a)原始數(shù)據(jù)f-k譜； (b)缺失數(shù)據(jù)f-k譜； (c)高分辨率方法重建數(shù)據(jù)f-k譜； (d)SL0 Radon變換方法重建數(shù)據(jù)f-k譜； (e)高分辨率方法重建數(shù)據(jù)f-k譜誤差(放大十倍)； (f)SL0 Radon變換方法重建數(shù)據(jù)f-k譜誤差(放大十倍)

圖4 振幅變化的隨機缺失數(shù)據(jù)重構(gòu)模型 (a)原始數(shù)據(jù)； (b)隨機缺失數(shù)據(jù)； (c)高分辨率方法重構(gòu)結(jié)果； (d)SL0 Radon變換方法重構(gòu)結(jié)果； (e)高分辨率方法重構(gòu)誤差剖面； (f)SL0 Radon變換方法重構(gòu)誤差剖面

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

圖 5a為實際地震數(shù)據(jù)，共有92道，每道取401個采樣點，采樣間隔為4ms(來自SeismicLab軟件包，http：//www-geo.phys.ualberta.ca/saig/SeismicLab)。將該地震數(shù)據(jù)隨機缺失32道作為缺失數(shù)據(jù)，如圖5b所示。高分辨率 Radon變換及SL0范數(shù)約束的Radon變換方法重構(gòu)結(jié)果分別如圖5c和圖5d所示。圖5e和5f分別為高分辨率Radon變換和SL0范數(shù)約束的Radon變換方法重構(gòu)結(jié)果與實際數(shù)據(jù)的誤差剖面，通過比較可看出,SL0范數(shù)約束的Radon變換在重構(gòu)時誤差較小，較好地保持了地震同相軸的連續(xù)性，恢復(fù)了數(shù)據(jù)的AVO特性。

圖5 實際數(shù)據(jù)隨機缺失重構(gòu) (a)原始數(shù)據(jù)； (b)隨機缺失數(shù)據(jù)； (c)高分辨率方法重構(gòu)結(jié)果 ；(d)SL0 Radon 變換方法重構(gòu)結(jié)果； (e)高分辨率方法重構(gòu)誤差； (f)SL0 Radon變換重構(gòu)誤差

4 結(jié)論

Cauchy準則和L1范數(shù)對數(shù)據(jù)的稀疏都存在局限性，它們不能充分反映數(shù)據(jù)稀疏特征，L0范數(shù)是表征數(shù)據(jù)稀疏度的最好方法，本文將L0范數(shù)稀疏約束應(yīng)用于Radon變換，提出SL0范數(shù)約束的Radon變換方法。該方法比高分辨率Radon變換法的分辨率更高，壓制假頻效果更好。在該方法中參數(shù)σ取值決定了SL0范數(shù)的稀疏度，σ越大，稀疏度越小，Radon變換分辨率越低；σ越小，稀疏度越大，Radon變換分辨率越高，聚焦效果越好。所以在實際應(yīng)用中要適當選取σ值，來提高Radon變換的效果。

地震道缺失重構(gòu)問題對地震數(shù)據(jù)處理意義重大。將該方法應(yīng)用于地震數(shù)據(jù)重構(gòu)中，通過理論模型實驗和實際地震資料處理，表明相比高分辨率Radon變換方法，該方法通過SL0范數(shù)稀疏約束Radon變換更能提高Radon變換的分辨率，更好地保持地震同相軸的連續(xù)性和恢復(fù)地震數(shù)據(jù)的AVO特性，對于不同的缺失數(shù)據(jù)都有更好的重構(gòu)效果。

[1] 王雄文，王華忠.基于壓縮感知的高分辨率平面波分解方法研究.地球物理學(xué)報，2014，57(9)：2946-2960. Wang Xiongwen，Wang Huazhong.A research of high-resolution plane-wave decomposition based on compressed sensing.Chinese Journal of Geophysics，2014，57(9)：2946-2960.

[2] 徐彥凱，曹思遠，何元.雙曲Radon-ASVD 方法壓制疊前地震數(shù)據(jù)隨機噪聲.石油地球物理勘探，2017，52(3)：451-457. Xu Yankai，Cao Siyuan，He Yuan.Hyperbolic Radon-ASVD method for suppressing random noise of pre-stack seismic data.OGP，2017，52(3)：451-457.

[3] 薛亞茹，唐歡歡，陳小宏.高階高分辨率Radon變換地震數(shù)據(jù)重建方法.石油地球物理勘探，2014，49(1)：95-100，131. Xue Yaru，Tang Huanhuan，Chen Xiaohong.Reconstruction method of seismic data using high-order high resolution Radon transform.OGP，2014，49(1)：95-100，131.

[4] 范景文，李振春，宋翔宇等.各向異性高分辨率Radon變換壓制多次波.石油地球物理勘探，2016，51(4)：665-669. Fan Jingwen，Li Zhenchun，Song Xiangyu et al.Suppression of multiple waves by anisotropic high resolution Radon transform.OGP，2016，51(4)：665-669.

[5] 謝俊法，孫成禹，韓文功.迭代拋物Radon變換法分離一次波與多次波.石油地球物理勘探，2014，49(1)：76-81. Xie Junfa，Sun Chengyu，Han Wengong.Separation of primary and multiple waves by iterative parabolic Radon transform.OGP，2014，49(1)：76-81.

[6] 劉喜武，劉洪，李幼銘.高分辨率Radon變換方法及其在地震信號處理中的應(yīng)用.地球物理學(xué)進展，2004，19(1)：8-15. Liu Xiwu，Liu Hong，Li Youming.High resolution radon transform and its application in seismic signal processing.Progress in Geophysics，2004，19(1)：8-15.

[7] 安鵬.基于高分辨率Radon變換的波場分離方法研究[學(xué)位論文].山東東營：中國石油大學(xué)(華東)，2009. An Peng.Wave Field Separation Technology Research Based on High-resolution Radon Transform [D].China University of Petroleum(East China)，Dongying，Shandong，2009.

[8] Sacchi M，Ulrych T.High resolution velocity gathers and offset space reconstruction.Geophysics，1995，60(4)：1169-1177.

[9] Sacchi M，Porsani M.Fast high resolution parabolic Radon transform.SEG Technical Program Expanded Abstracts,1999，18：1477-1480.

[10] Trad D，Ulrych T，Sacchi M.Latest view of the sparse Radon transform.Geophysics，2003，68(1)：386-399.

[11] Figueiredo M A T，Nowak R D，Wright S J.Gradient projection for sparse reconstruction：application to compressed sensing and other inverse problems.IEEE Journal of Selected Topics in Signal Processing，2007，1(4)：586-597.

[12] Chen S S，Donoho D L，Saunders M A.Atomic decomposition by basis pursuit.SIAM Journal of Scientific Computing，1998，20(1)：33-61.

[13] Mohimani H，Babaie-Zadeh M and Jutten C.Fast sparse representation using smoothed L0norm.IEEE Tran-sactions on Signal Processing，2007，46(6)：1-12.

[14] Mohimani H，Babaie-Zadeh M，Jutten C.A fast ap-proach for overcomplete sparse decomposition based on smoothed L0norm.IEEE Transactions on Signal Processing，2009，57(1)：289-301.

[15] Hyder M M and Mahata K.An improved smoothed L0approximation algorithm for sparse representation.IEEE Transactions on Signal Processing，2010，58(4)：2194-2205.

[16] Zayyani H，Babaie-Zadeh M.Thresholded smoothed-L0dictionary learning for sparse representations.IEEE International Conference on Acoustics，Speech and Signal Processing，2009，1825-1828.

[17] Ghalehjegh S H，Babaie-Zadeh M，Jutten C.Fast block-sparse decomposition based on SL0.9th International Conference on Latent Variable Analysis and Signal Separation，2010，426-433.

[18] 林婉娟，趙瑞珍，李浩.用于壓縮感知信號重建的NSL0算法.新型工業(yè)化，2011，1(7)：78-84. Lin Wanjuan，Zhao Ruizhen，Li Hao.The NSL0algorithm for compressive sensing signal reconstruction.Industrialization，2011，1(7)：78-84.

[19] 祁銳，李宏偉，張玉潔.基于改進光滑L0范數(shù)的塊稀疏信號重構(gòu)算法.計算機工程，2015，41(11)：294-298. Qi Rui，Li Hongwei，Zhang Yujie.Block sparse signal reconstruction algorithm based on improved smoothed L0norm.Computer Engineering，2015，41(11)：294-298.

[20] 楊良龍，趙生妹，鄭寶玉等.基于SL0壓縮感知信號重建的改進算法.信號處理，2012，28(6)：834-841. Yang Lianglong，Zhao Shengmei，Zheng Baoyu et al.The improved reconstruction algorithm for compressive sensing on SL0.Signal Processing，2012，28(6)：834-841.

[21] 曹靜杰，王彥飛，楊長春.地震數(shù)據(jù)壓縮重構(gòu)的正則化與零范數(shù)稀疏最優(yōu)化方法.地球物理學(xué)報，2012，55(2)：596-607. Cao Jingjie，Wang Yanfei，Yang Changchun.Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization.Chinese Journal of Geophysics，2012，55(2)：596-607.

[22] 張良，韓立國，劉爭光等.基于壓縮感知和Contourlet變換的地震數(shù)據(jù)重建方法.石油物探，2017，56(6)：804-811. Zhang Liang,Han Liguo,Liu Zhengguang et al.Seismic data reconstruction based on compressed sending and Contourlet transform.GPP,2017，56(6)：804-811.