李　欣　楊　婷　孫文博　王貝貝

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

1　引言

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

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

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

地震插值可被視為一個反演問題，其對應的正演問題就是地震采集，可寫成如下公式

Φx=b

(1)

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

Φx=ΦΨ*s=As=b

(2)

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

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

(3)

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

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

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

(4)

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

minF(s)s.t.As=b

(5)

(6)

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

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

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

(7)

根據(jù)凸分析理論， {s|As=b}是一個凸集，式(7)可用凸集投影方法求解。作為約束優(yōu)化問題，投影梯度法是非常有效的一種方法。梯度投影法是最優(yōu)化算法中解約束優(yōu)化問題的一類重要方法，它基于目標函數(shù)的負梯度方向對迭代解進行初步升級，然后為了防止升級后的解不滿足約束條件，通過投影方式將該初步升級后的解投影到約束條件形成的凸集合中，從而實現(xiàn)解的一次完整的升級迭代。這里介紹一種非常快速的投影梯度法求解式(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　曲波變換

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

(8)

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

4　數(shù)值計算

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

4.1　炮集數(shù)據(jù)實驗

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

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

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

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

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

4.2　疊加剖面實驗

進一步利用疊后數(shù)據(jù)研究光滑L1范數(shù)方法的效率。圖5a是一個疊加剖面，包含130道，道間距是25m， 401個時間采樣點，時間采樣間隔是2ms。不完整采樣數(shù)據(jù)如圖5b所示，其中40%的原始數(shù)據(jù)被隨機地去掉。L1范數(shù)方法的最大迭代次數(shù)是20，圖6a是光滑L1范數(shù)的插值結果； L0范數(shù)的內(nèi)部循環(huán)次數(shù)是2，外部循環(huán)次數(shù)是10，其結果如圖6b； FISTA和SPGL1方法的插值結果見圖7。這些方法的CPU時間、信噪比和相對誤差見表2。與炮集數(shù)據(jù)實驗結果相同，光滑L1范數(shù)方法和光滑L0范數(shù)方法具有相同的速度，比FISTA方法快，所用時間約為SPGL1方法的1/3，因此本文新方法可顯著降低計算成本。由于光滑L0方法是對0范數(shù)的近似，而光滑L1方法是對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)插值結果

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

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

5　結論與建議

本文提出采用Huber函數(shù)作為目標函數(shù)構建地震數(shù)據(jù)重建反演模型。該函數(shù)是一個連續(xù)可微的光滑函數(shù)，可作為L1范數(shù)的近似，因此在求解時能夠采用現(xiàn)有的最優(yōu)化方法求解；針對建立的反演模型提出了一種梯度投影法求解，通過投影方法實現(xiàn)反演模型的快速求解。由于曲波變換能夠直接稀疏地表達曲線形狀的同相軸，不需要對數(shù)據(jù)進行分塊處理，所以本文采用曲波變換作為稀疏變換，該變換的另一個特點是其正交性可用來加速計算。模擬和真實數(shù)據(jù)試驗證明了本文建立模型和求解方法的有效性；數(shù)值結果表明，該新方法比快速迭代軟閾值法更快，約為SPGL1方法計算時間的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)合改進的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]張良，韓立國，徐德鑫等.基于壓縮感知技術的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]魏小強，雷秀麗，馬慶珍.基于多道奇異譜分析的三維地震數(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ù)壓縮重構的正則化與零范數(shù)稀疏最優(yōu)化方法.地球物理學報，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.