一種基于光滑L1范數(shù)的地震數(shù)據(jù)插值方法

李　欣　楊　婷　孫文博　王貝貝

(中海油研究總院，北京 100028)

1　引言

受現(xiàn)場(chǎng)采集因素影響，原始地震數(shù)據(jù)通常不滿足采樣定理，這會(huì)影響去噪、AVO分析、多次波去除以及偏移的效果[1，2]。通過插值處理能得到高密度地震數(shù)據(jù)[3]，從而減少低頻陰影并提高成像精度[4]。目前主流的插值方法都是依托信號(hào)處理[1-3，5]，其中基于稀疏變換類方法是假設(shè)地震數(shù)據(jù)在變換域內(nèi)具有稀疏性，即將地震插值當(dāng)作稀疏反問題?；诘卣鹦盘?hào)在空間上的有限帶寬性，傅里葉變換和線性拉冬變換可用于恢復(fù)不規(guī)則數(shù)據(jù)[1，5-9]。張良等[10]將Shearlet變換應(yīng)用于地震數(shù)據(jù)重建。Herrmann 等[11]提出利用曲波變換壓縮地震數(shù)據(jù)，該變換比其他類型變換能更稀疏地表征地震數(shù)據(jù)。最近，魏小強(qiáng)等[12]將矩陣完備理論用于地震數(shù)據(jù)插值。Kreimer等[4]提出適用于高維數(shù)據(jù)插值的張量完備算法?；谙∈璺囱莸姆椒ㄐ枨蠼獠贿m定性問題，因此可采用稀疏優(yōu)化法尋找理想解法。隨著高維插值方法在地震數(shù)據(jù)中的推廣應(yīng)用，需探尋針對(duì)地震數(shù)據(jù)處理中大規(guī)模問題的快速求解方法。在地震插值方面，Amba等[13]提出求解地震插值的凸集投影(POCS)法。Zwartijes等[14]采用加權(quán)迭代最小二乘法解插值問題。Herrmann等[11]把迭代軟閾值(IST)法引入地震數(shù)據(jù)恢復(fù)中。Van den Berg等[15]提出穩(wěn)健的譜梯度投影(SPGL1)方法。Wang等[16]提出L1范數(shù)信賴域方法實(shí)現(xiàn)地震數(shù)據(jù)的插值。這些方法都不同程度地改善了插值效果。但隨著地震數(shù)據(jù)量的迅猛增大，亟待研發(fā)更快、更穩(wěn)健的方法以降低CPU時(shí)間并獲取可靠的插值數(shù)據(jù)。利用曲波變換[17]，曹靜杰等[18]提出基于光滑L0范數(shù)的梯度投影地震插值方法，大幅度縮短了運(yùn)算時(shí)間；隨后，還提出了基于非凸Lp范數(shù)的重建方法并應(yīng)用于地震反褶積[19]。

在分析、對(duì)比各種稀疏優(yōu)化地震插值方法的時(shí)效性和穩(wěn)健性的基礎(chǔ)上，本文選定Huber函數(shù)作為L1范數(shù)的光滑逼近函數(shù)，并構(gòu)建了快速平滑的L1范數(shù)正則化模型，提出了一種快速、穩(wěn)健的基于Huber函數(shù)正則化梯度投影法。文中選擇曲波變換作為稀疏變換方法，充分利用其正交性可加快計(jì)算速度。數(shù)值計(jì)算結(jié)果證明該方法具有較好的插值結(jié)果和較高的計(jì)算效率；同時(shí)，指出了稀疏優(yōu)化在地震數(shù)據(jù)處理中的適用領(lǐng)域。

2　地震數(shù)據(jù)插值的稀疏優(yōu)化模型

地震插值可被視為一個(gè)反演問題，其對(duì)應(yīng)的正演問題就是地震采集，可寫成如下公式

Φx=b

(1)

式中：Φ為采樣過程；x為反射系數(shù)模型；b為采樣數(shù)據(jù)。由于采樣不足，理論上該方程有無數(shù)多個(gè)解，根據(jù)反演理論，可通過一些先驗(yàn)信息獲得理想的解。地震同相軸可以通過一些變換稀疏地表示，所以變換域解的稀疏性常常被用作先驗(yàn)信息。地震處理中最常用的就是傅里葉變換[1，14]、線性拉冬變換[9]、拋物線拉冬變換[10]和曲波變換[11]等。近年來，高斯束也被用于地震數(shù)據(jù)分解[19]。這種先驗(yàn)信息在信號(hào)處理中也很常見，如信號(hào)壓縮與傳輸、去噪和去卷積[19，20]。如果s=Ψx是稀疏的，其中Ψ是一種變換，那么式(1)可改寫成

Φx=ΦΨ*s=As=b

(2)

式中：s是一個(gè)以向量表示的離散曲波系數(shù)，其稀疏性可用一個(gè)擬范數(shù)‖s‖0度量，該擬范數(shù)表示s的非零元素的個(gè)數(shù)；A=ΦΨ*，是復(fù)合矩陣。于是式(2)即等價(jià)為如下優(yōu)化問題

min‖s‖0s.t.As=b

(3)

上述優(yōu)化問題的含義是在滿足As=b的眾多解中，找到最稀疏的那一個(gè)解。 然而， ‖s‖0不是一個(gè)真

正的范數(shù)，它不連續(xù)且不可求導(dǎo)。為了能采用凸優(yōu)化方法求解，該問題一般變?yōu)榻仆箖?yōu)化問題求解

min‖s‖1s.t.As=b

(4)

這個(gè)問題可通過線性優(yōu)化和內(nèi)點(diǎn)法求解，但是計(jì)算速度低[21，22]。由于L1范數(shù)是不可求導(dǎo)的，所以這類型的問題不能通過共軛梯度法和牛頓型方法直接求解。一種有效的策略是將L1范數(shù)替換成光滑的近似函數(shù)。因此，式(4)可變成

minF(s)s.t.As=b

(5)

(6)

該函數(shù)是光滑的，且當(dāng)a→0時(shí)與|s|非常相似[23]。圖1是a=0.0001時(shí)的Huber函數(shù)形態(tài)及其原點(diǎn)局部放大顯示(在原點(diǎn)為零)。Huber函數(shù)是一個(gè)L1和L2的混合范數(shù)： 當(dāng)變量小時(shí)表現(xiàn)為L2范數(shù)； 當(dāng)變量大時(shí)表現(xiàn)為L1范數(shù)。L2范數(shù)到L1范數(shù)的光滑過度是依靠參數(shù)a。Huber函數(shù)被廣泛地應(yīng)用于地球物理領(lǐng)域，如討論線性擬合的魯棒性問題，但很少應(yīng)用于解的度量方面。此處將其作為正則化項(xiàng)以得到稀疏的解決方案。

圖1　a=0.0001時(shí)的Huber函數(shù)形態(tài)(a)及其局部放大(b)

從上面討論可知： ①Huber函數(shù)中存在一個(gè)參數(shù)控制其擬合度，且是可微的； ②fhuber(s)在原點(diǎn)等于零，在參數(shù)選擇合理的情況下fhuber(s)與|s|非常相似，因此它是L1范數(shù)的最佳近似。下面給出式(5)的以fhuber(s)為目標(biāo)函數(shù)的表達(dá)形式

(7)

根據(jù)凸分析理論， {s|As=b}是一個(gè)凸集，式(7)可用凸集投影方法求解。作為約束優(yōu)化問題，投影梯度法是非常有效的一種方法。梯度投影法是最優(yōu)化算法中解約束優(yōu)化問題的一類重要方法，它基于目標(biāo)函數(shù)的負(fù)梯度方向?qū)Φ膺M(jìn)行初步升級(jí)，然后為了防止升級(jí)后的解不滿足約束條件，通過投影方式將該初步升級(jí)后的解投影到約束條件形成的凸集合中，從而實(shí)現(xiàn)解的一次完整的升級(jí)迭代。這里介紹一種非?？焖俚耐队疤荻确ㄇ蠼馐?7)，該方法要比目前的方法更節(jié)省CPU資源。下面給出求解式(7)的具體步驟。

(1)給出最大迭代次數(shù)L、 Huber函數(shù)參數(shù)a=0.0001和初始迭代次數(shù)k=0； 令初始解s0為As=b的L2范數(shù)解。

(4)輸出最終解s=sk。

3　曲波變換

基于對(duì)地震數(shù)據(jù)的正交性和優(yōu)良的稀疏表達(dá)能力，本文選擇曲波變換作為稀疏變換。傅里葉變換，拉冬類變換和小波變換都可以用來稀疏地表示地震數(shù)據(jù)，但是對(duì)于基于稀疏反演的插值模型來說，要求地震數(shù)據(jù)在變換域中越稀疏越好，小波變換對(duì)地震數(shù)據(jù)有一定的壓縮性，但是其變換后的稀疏性沒有曲波變換好。曲波變換作為一個(gè)正交、多方向、多尺度、各向異性和局部變換[17]，已被證明是適合地震數(shù)據(jù)的最稀疏變換[11]之一。另一方面，地震數(shù)據(jù)的同相軸大多是曲線形狀的，傅里葉變換只能稀疏地表示直線形狀的同相軸，基于傅里葉變換的地震數(shù)據(jù)插值需要將地震數(shù)據(jù)分塊，使地震數(shù)據(jù)在時(shí)間和空間窗口中近似為線性的。拉冬類變換存在固有的缺點(diǎn)，首先它們不是嚴(yán)格可逆的，另外對(duì)于一些不是規(guī)則線性、雙曲形狀或拋物形狀的同相軸，拉冬類變換對(duì)地震數(shù)據(jù)的稀疏表示結(jié)果并不理想。因此本文選擇存在嚴(yán)格逆變換，并且不需要對(duì)地震數(shù)據(jù)進(jìn)行分塊處理的曲波變換作為稀疏變換，它可以表示成曲波函數(shù)φj,k,l(x)與數(shù)據(jù)f(x)的內(nèi)積

(8)

該變換形式可寫成s=Ψx，Ψ即是曲波變換。曲波正變換和逆變換的計(jì)算成本是O(N2logN)[17]，N表示離散情況下數(shù)據(jù)的長度。離散的曲波變換是一個(gè)緊框架，因此伴隨算子Ψ*等價(jià)于Ψ的偽逆[19]，形式上，曲波變換的逆可寫成x=Ψ*s。

4　數(shù)值計(jì)算

為了測(cè)評(píng)本文方法的計(jì)算效率，選取兩套實(shí)際數(shù)據(jù)進(jìn)行處理，并與常用的先進(jìn)稀疏方法做對(duì)比。首先給出一個(gè)炮集數(shù)據(jù)實(shí)驗(yàn)，第二個(gè)實(shí)驗(yàn)將進(jìn)一步測(cè)試其對(duì)疊后數(shù)據(jù)的效果。本次選用了三種主流方法來做對(duì)比研究。

4.1　炮集數(shù)據(jù)實(shí)驗(yàn)

第一個(gè)炮數(shù)據(jù)的接收點(diǎn)間距是25m，采樣間隔是2ms。數(shù)據(jù)包含了115道，每一道有600個(gè)采樣點(diǎn)，模擬隨機(jī)的采樣 69 道。圖2a是原始數(shù)據(jù)，圖2b是采樣后的數(shù)據(jù)。分別利用光滑L1范數(shù)方法、光滑L0范數(shù)方法[18]、快速迭代閾值法(FISTA)[20]和SPGL1方法[15]進(jìn)行計(jì)算。光滑L1范數(shù)方法的最大迭代次數(shù)是15次； 光滑L0范數(shù)方法有兩個(gè)循環(huán)，每個(gè)外部循環(huán)次數(shù)是2， 內(nèi)部循環(huán)次數(shù)是4；FISTA方法最大迭代次數(shù)是20； SPGL1方法的最大迭代次數(shù)是30?；谶@些參數(shù)，這些方法將獲得相似的插值結(jié)果。這些方法的CPU時(shí)間、信噪比和相對(duì)誤差見表1，基于光滑L1范數(shù)和光滑L0范數(shù)的插值結(jié)果見圖3，F(xiàn)ISTA和SPGL1的插值結(jié)果見圖4。從結(jié)果可以看出，光滑L1范數(shù)方法和光滑L0范數(shù)方法具有相同的計(jì)算速度，都比FIST方法快，約為SPGL1方法的1/3。

表1　光滑L1、光滑L0、FISTA和SPGL1方法炮集數(shù)據(jù)對(duì)比

圖2　原始炮集數(shù)據(jù)(a)及其重采樣數(shù)據(jù)(b)

圖3　炮集數(shù)據(jù)的光滑L1方法(a)與光滑L0方法(b)的插值結(jié)果

圖4　炮集數(shù)據(jù)的FISTA方法(a)與SPGL1方法(b)的插值結(jié)果

4.2　疊加剖面實(shí)驗(yàn)

進(jìn)一步利用疊后數(shù)據(jù)研究光滑L1范數(shù)方法的效率。圖5a是一個(gè)疊加剖面，包含130道，道間距是25m， 401個(gè)時(shí)間采樣點(diǎn)，時(shí)間采樣間隔是2ms。不完整采樣數(shù)據(jù)如圖5b所示，其中40%的原始數(shù)據(jù)被隨機(jī)地去掉。L1范數(shù)方法的最大迭代次數(shù)是20，圖6a是光滑L1范數(shù)的插值結(jié)果； L0范數(shù)的內(nèi)部循環(huán)次數(shù)是2，外部循環(huán)次數(shù)是10，其結(jié)果如圖6b； FISTA和SPGL1方法的插值結(jié)果見圖7。這些方法的CPU時(shí)間、信噪比和相對(duì)誤差見表2。與炮集數(shù)據(jù)實(shí)驗(yàn)結(jié)果相同，光滑L1范數(shù)方法和光滑L0范數(shù)方法具有相同的速度，比FISTA方法快，所用時(shí)間約為SPGL1方法的1/3，因此本文新方法可顯著降低計(jì)算成本。由于光滑L0方法是對(duì)0范數(shù)的近似，而光滑L1方法是對(duì)1范數(shù)的近似，光滑L0方法能夠更加逼近0范數(shù)優(yōu)化問題，因此數(shù)值效果要優(yōu)于光滑L1方法。本文的方法的優(yōu)越性在于，迭代只需一層循環(huán)，需調(diào)節(jié)的參數(shù)比光滑L0方法少，且更易調(diào)節(jié)，因此本文算法更加簡便，容易操作。

圖5　原始疊加剖面(a)及其重采樣數(shù)據(jù)(b)

圖6　疊加剖面的光滑L1方法(a)與光滑L0方法(b)插值結(jié)果

圖7　疊加剖面的FISTA方法(a)與SPGL1方法插值結(jié)果(b)

光滑L1光滑L0FISTASPGL1CPU時(shí)間/s565480163信噪比22．180523．787922．709422．9518相對(duì)誤差/(%)7．786．477．327．12

5　結(jié)論與建議

本文提出采用Huber函數(shù)作為目標(biāo)函數(shù)構(gòu)建地震數(shù)據(jù)重建反演模型。該函數(shù)是一個(gè)連續(xù)可微的光滑函數(shù)，可作為L1范數(shù)的近似，因此在求解時(shí)能夠采用現(xiàn)有的最優(yōu)化方法求解；針對(duì)建立的反演模型提出了一種梯度投影法求解，通過投影方法實(shí)現(xiàn)反演模型的快速求解。由于曲波變換能夠直接稀疏地表達(dá)曲線形狀的同相軸，不需要對(duì)數(shù)據(jù)進(jìn)行分塊處理，所以本文采用曲波變換作為稀疏變換，該變換的另一個(gè)特點(diǎn)是其正交性可用來加速計(jì)算。模擬和真實(shí)數(shù)據(jù)試驗(yàn)證明了本文建立模型和求解方法的有效性；數(shù)值結(jié)果表明，該新方法比快速迭代軟閾值法更快，約為SPGL1方法計(jì)算時(shí)間的1/3。

隨著壓縮感知理論的發(fā)展，許多稀疏優(yōu)化模型，稀疏變換和求解方法都可用于地震數(shù)據(jù)重建問題。L1范數(shù)是目前常用的一種的稀疏約束，此外，Lp(0

[1]Liu B.Multi-dimensional Reconstruction of Seismic Data [D].University of Alberta，Edmonton，Canada，2004.

[2]Naghizadeh M，Sacchi M.Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data.Geophysics，2010，75(6)：WB189-WB202.

[3]Spitz S.Seismic trace interpolation in the F-X domain.Geophysics，1991，56(6)：785-794.

[4]Kreimer N，Edmonton A，Sacchi M.Reconstruction of seismic data via tensor completion.2012 IEEE Statistical Signal Processing Workshop (SSP)，2012，29-32.

[5]Duijndam A，Schonewille M，Hindriks C.Reconstruction of band-limited signals，irregularly sampled along one spatial direction.Geophysics，1999，64(2)：524-538.

[6]王本鋒，陳小宏，李景葉等.POCS聯(lián)合改進(jìn)的Jitter采樣理論曲波域地震數(shù)據(jù)重建.石油地球物理勘探， 2015，50(1)：20-28.

Wang Benfeng，Chen Xiaohong，Li Jingye et al.Seismic data reconstruction based on POCS and improved Jittered sampling in the curvelet domain.OGP，2015，50(1)：20-28.

[7]Sacchi M，Ulrych T，Walker C.Interpolation and ex-trapolation using a high resolution discrete Fourier transform.IEEE Transactions on Signal Processing，1998，46()：31-38.

[8]Xu S，Zhang Y，Pham D et al.Antileakage Fourier transform for seismic data regularization.Geophysics，2005，70(4):V87-V95.

[9]薛亞茹，唐歡歡，陳小宏.高階高分辨率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.

[10]張良，韓立國，徐德鑫等.基于壓縮感知技術(shù)的Shearlet變換重建地震數(shù)據(jù).石油地球物理勘探，2017，52(2)：220-225.

Zhang Liang，Han Liguo，Xu Dexin et al.Reconstruction of seismic data of Shearlet transform based on compressed sensing technology.OGP，2017，52(2)：220-225.

[11]Herrmann F and Hennenfent G.Non-parametric seismic data recovery with curvelet frames.Geophysical Journal International，2008，173(1)：233-248.

[12]魏小強(qiáng)，雷秀麗，馬慶珍.基于多道奇異譜分析的三維地震數(shù)據(jù)規(guī)則化方法.石油地球物理勘探，2014，49(5)：846-851.

Wei Xiaoqiang，Lei Xiuli，Ma Qingzhen.3D seismic data regularization based on multi-channel singularspectrum analysis.OGP，2014，49(5)：846-851.

[13]Amba R，Kabir N.3D interpolation of irregular data with a POCS algorithm.Geophysics，2006，71(6)：E91-E97.

[14]Zwartijes P，Sacchi M.Fourier reconstruction of nonuniformly sampled，aliased seismic data.Geophysics，2007，72(1)：V21-V32.

[15]Van den Berg E，Michael F.Probing the pareto frontier for basis pursuit solutions.SIAM Journal on Science Computing，2008，31(2)：890-912.

[16]Wang Y F，Cao J J，Yang C C.Recovery of seismic wavefields based on compressive sensing by a l1-norm constrained trust region method and the piecewise random sub-sampling.Geophysical Journal International，2011，187(1)：199-213.

[17]Candes E.Compressive sampling.Proceedings of In-ternational Congress of Mathematicians，European Mathematical Society Publishing House，Madrid，Spain，2006，33-52.

[18]曹靜杰，王彥飛，楊長春.地震數(shù)據(jù)壓縮重構(gòu)的正則化與零范數(shù)稀疏最優(yōu)化方法.地球物理學(xué)報(bào)，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.

[19]曹靜杰.基于廣義高斯分布和非凸Lp范數(shù)正則化的地震稀疏盲反褶積.石油地球物理勘探，2016，51(3)：428-433.

Cao Jingjie.Seismic sparse blind deconvolution based on generalized Gaussian distribution and non-convex Lpnorm regularization.OGP，2016，51(3)：428-433.

[20]Lustig M，Donoho D，Pauly J et al.The application of compressed sensing for rapid MR imaging.Magnetic Resonance on Medicine，2007，58(6)：1182-1195.

[21]Chen S，Donoho D，Saunders M.Atomic decomposi-tion by basis pursuit.Siam Review，2001，43(1)：129-159.

[22]Candes E，Tao T.Decoding by linear programming.IEEE Transactions on Information Theory，2005，51，4203-4215.

[23]Bube K P，NemethT.Fast line searches for the robust solution of linear systems in the hybrid L1/L2 and Huber norms.Geophysics，2007，72(2)：A13-A17.