基于L0范數稀疏約束的地震數據反褶積

梁東輝，陳生昌

(浙江大學地球科學系，浙江杭州310027)

隨著油氣勘探進程的不斷推進，需要尋找地質構造復雜、規模較小及埋藏較深的油氣藏，這就要求地震勘探資料要有更高的分辨率。反褶積是地震數據處理中的一個基本環節，它的作用是壓縮地震記錄中的地震子波，同時可以壓制鳴震和多次波，提高地震記錄的縱向分辨率[1]。實際地震記錄中不可避免的包含有未知噪聲，所以反褶積是一個病態的反問題，需要采取正則化方法來求解。求解此類問題時常用的阻尼最小二乘法(預白化法)在子波自相關矩陣的對角線上加一個小的擾動量，等價于零階二次正則化，雖然這樣能夠很好地提高算法的穩定性，但這是一種較寬泛的二次約束，求得的解是連續的、光滑的，分辨率不夠高。為了提高地震勘探資料的分辨率，可運用稀疏優化手段，使反褶積結果稀疏化，這就是稀疏約束反褶積。傳統的反褶積方法需要假設子波為最小相位和反射系數是白噪聲，而實際地震資料中子波通常不是最小相位的，反射系數也是非白噪的[2-3]，這就導致傳統的反褶積效果不理想。稀疏約束反褶積則避免了這兩種假設，它利用反射系數由稀疏脈沖序列組成這一先驗信息，從帶有噪聲的地震記錄中提取出反射系數的時間位置和振幅值，同時去除噪聲,可得到比常規反褶積具有更高分辨率的反演結果。

稀疏約束反褶積方法有很多，前人已做了很多研究，如L1范數稀疏約束反褶積[4-5]，Lp范數反褶積[6]，L1-L2范數聯合約束稀疏脈沖反演的應用[7]，基于柯西準則的預條件共軛梯度法[8-11]，地震反褶積中的重加權策略[12]等?？挛鳒蕜t約束稀疏反褶積對于較大的反射系數能準確提取，但對于較小的反射系數會有壓制作用[13]。L1范數約束放松了對解的稀疏性的約束，所求的并不是最稀疏的解，并且對于較長的信號其計算復雜度較高[14]。

稀疏約束反褶積由兩部分組成。一是對噪聲的約束:實際地震記錄中的噪聲一般認為是高斯分布的，子波與所求得的反射系數褶積得到正演地震記錄，對噪聲的約束即是使實際地震記錄與正演地震記錄的殘差平方和最小。二是對反射系數的約束:這部分用來限制反射系數的先驗概率分布情況。對于一個向量來說，L0范數就是其中非零元素的個數，因此是向量稀疏性的最佳度量，L0范數稀疏優化已在地震數據壓縮重構中得到很好的應用[15-16]。層狀地層模型的地下地層反射系數是一系列稀疏脈沖，從理論上來講，對反射系數采取L0范數約束能求得最為稀疏的解，可以有效地提高地震資料的分辨率，有利于后期的處理和解釋。

我們將L0范數稀疏約束引入到地震資料反褶積處理中，并通過兩組模型試驗，對比柯西準則，L1范數和L0范數稀疏約束反褶積的效果。

1 方法原理

由震源激發的脈沖信號經過大地濾波器的作用變成一個長度為幾十毫秒的波形，稱為地震子波。根據褶積模型，地震記錄相當于子波和反射系數序列的褶積，又因為野外實際地震數據采集時不可避免地會含有噪聲，故地震記錄可表示為

(1)

式中:d表示地震記錄;ω0表示地震子波;r表示反射系數序列;n表示噪聲?？蓪懗扇缦戮仃囆问?

(2)

這里ω是由地震子波ω0組成的褶積矩陣。

由于地震子波的影響，來自相鄰界面的反射波重疊在一起，難以區分開來，因此需要壓縮子波來提高地震勘探資料的縱向分辨率。傳統的反褶積方法是設計一個與子波性質相反的濾波因子，然后與地震記錄褶積，得到期望的窄脈沖，這個過程稱為反褶積。

1.1 約束優化

稀疏約束反褶積運用正則化方法直接求解方程(2)，假設子波已知或已準確提取出，優化目標函數的一般形式可表示為

(3)

其中，λ>0，是反射系數稀疏性與準確性之間的權衡因子。該目標函數前一部分(J0)對殘差進行約束，保證所求的結果與地震記錄數據符合，也就是對噪聲的約束,這部分一般采用最小平方約束;后一部分(λJr)稱為正則化項或罰函數，是為了對反射系數進行稀疏約束，保證求解結果稀疏，不同的稀疏約束方法選取不同的罰函數。

也可以從貝葉斯定理出發，給定噪聲的分布概率和反射系數的先驗分布概率，然后推導出相應的稀疏約束形式。貝葉斯定理可以表示為

(4)

其中，P(r|d)表示反射系數的后驗概率分布;P(d|r)表示給定反射系數時觀測數據的符合程度，稱為似然函數，等于噪聲的分布概率;P(r)表示反射系數的先驗概率分布。

假設噪聲滿足高斯分布，且均值為零，標準差為σn，單獨一個噪聲分量的概率密度函數為

(5)

噪聲整體的分布概率為

(6)

又因為

(7)

所以似然函數

(8)

(9)

反射系數總體的分布概率為

(10)

根據(4)式，反射系數的后驗概率分布為

(11)

兩邊同時取對數，得

(12)

使后驗概率分布最大，等價于求解優化函數:

(13)

(13)式就是柯西準則稀疏約束反褶積的優化目標函數。

同理，如果我們假設噪聲滿足高斯分布，均值為0，反射系數滿足均值為0的拉普拉斯分布，可根據貝葉斯定理推導出L1范數稀疏約束反褶積的優化目標函數:

(14)

1.2 L0范數約束反褶積

稀疏約束反褶積希望求得稀疏的反射系數，即較少的非零元素值。反射系數序列的L0范數就是其中非零元素的個數，因此采用L0范數稀疏約束來求解反褶積問題是一種理想的策略，理論上能求得最為稀疏的解。L0范數要解決的優化問題是

(15)

其中，δ是與噪聲有關的參數，表示對噪聲水平的估計量。(15)式表示在殘差滿足一定條件時，要求反射系數序列中非零元素數量最少。通過拉格朗日乘子法可將(15)式轉化為無約束問題，構建L0范數稀疏約束反褶積優化目標函數:

(16)

如果按照貝葉斯定理，假設噪聲滿足高斯分布，反射系數為稀疏分布，同樣可推導出(16)式。解此優化問題的方法很多，如迭代重加權最小二乘(IRLS)算法[17]、基追蹤(BP)算法[18]、貪婪算法[19-20]、迭代硬閾值法等[21]。IRLS算法構建加權的二范數來近似零范數，需要進行矩陣求逆運算，不適合大規模數據的處理;BP算法將L0范數最小優化轉化為L1范數最小優化問題，放松了對稀疏性的限制并且計算量較大;貪婪算法，如匹配追蹤、正交匹配追蹤等，求解精度較低，且不適合處理欠稀疏數據;迭代硬閾值計算形式簡單，計算量小，有利于處理大數據量的地震勘探資料。本文選取迭代硬閾值求解L0范數稀疏約束優化問題，其迭代格式為

醫學和藥學是實踐性非常強的學科，課堂學生內需要接受的知識繁多枯燥，以教師、課堂、教材為中心的教學方式無法滿足學生可持續發展的需要。如果沒有課外科技活動作為學習內容和形式的補充，學生的學習積極性會受到一定的影響。因此，在醫藥院校中，學生課外科技活動在人才培養中占有重要地位，它能在校內提供學生理論聯系實際的機會，培養學生的創新思維，提高人才培養質量，是教育行業“供給側結構性改革”的途徑之一。

(17)

其中，μ為步長。為了保證迭代收斂[22]，步長的選取范圍要滿足

(18)

Ha為硬閾值函數:

(19)

a>0表示閾值，計算時先給定一個較大的閾值:

(20)

每次迭代后降低閾值，再進行下一次迭代。

因為實際地震記錄中的噪聲可認為是均值為0，標準差為σn的高斯分布噪聲。

(21)

迭代終止的條件可設定為

(22)

其中，Nd表示觀測值d的個數。

L0范數稀疏約束反褶積的計算步驟如下:

1) 給定反射系數初值r0為0，估計噪聲的標準差σn，迭代次數n=0;

2) 選取合適的步長μ及閾值初值;

3) 按照(17)式由rn計算rn+1;

4) 閾值下降，迭代次數n增加1;

5) 重復步驟3)和步驟4),直至滿足迭代終止條件，輸出反褶積結果。

2 模型試驗

選取延遲30ms，主頻為25Hz的雷克子波(圖1)作為下面兩個模型試驗中的地震子波，采樣間隔為2ms。

圖1 試驗中所使用的子波

2.1 模型試驗一

構建一個簡單的稀疏分布反射系數序列(圖2a)，與圖1中的子波褶積并加入標準差為0.01的高斯噪聲，得到相應的合成地震記錄(圖2b)。

分別采用柯西準則約束、L1范數約束和L0范數約束進行反褶積計算，結果如圖3所示。對比圖3a,圖3b,圖3c和圖3d可見，上面兩個橢圓中有幾個較弱的反射系數，柯西準則約束反褶積結果未能顯示，L1范數約束能顯示但振幅值偏弱，L0范數約束全部顯示且較清晰;下面兩個橢圓中各有1個較弱的反射系數，柯西準則和L1范數約束反褶積結果未能顯示，L0范數約束結果中有顯示。說明L0范數對弱反射系數保護能力最強。

圖2 簡單的稀疏反射系數模型(a)及其合成地震記錄(b)

圖3 稀疏反射系數序列(a)和柯西準則(b)，L1范數(c)，L0范數(d)稀疏約束反褶積結果

2.2 模型試驗二

構建一個復雜的反射系數模型(圖4a)，與圖1所示子波褶積并加入標準差為0.01的高斯噪聲，得到相應的地震記錄(圖4b)。

分別采用柯西準則約束、L1范數約束和L0范數約束進行反褶積計算，結果如圖5a，圖6a和圖7a 所示。紅色表示正的反射系數，藍色表示負的反射系數。3張圖中黑色方框對應的放大圖示分別為圖5b，圖6b和圖7b,圖中箭頭所指處有3條較弱的反射系數同相軸，在圖5b中基本沒有顯示，在圖6b 和圖7b中有顯示，說明柯西準則約束對弱反射系數有壓制作用，L1范數和L0范數約束能夠較好地保護弱反射系數。

圖5b，圖6b和圖7b中黑色橢圓內有兩條藍色的反射系數同相軸，在圖5b中顯示不全，圖6b中顯示的振幅較弱，在圖7b中清晰顯示。圖8為3種稀疏約束反褶積結果的波形剖面。圖8a,圖8b 和圖8c中的紅色橢圓內有兩條負的(向左)反射系數同相軸，可以明顯看出，柯西準則稀疏約束反褶積(圖8a)只顯示了左邊一部分;L1范數稀疏約束反褶積(圖8b)顯示較全，但幅值較弱;L0范數稀疏約束反褶積(圖8c) 清晰顯示了這兩條負的反射系數同相軸。

表1為模型試驗二中3種反褶積方法的迭代次數與計算時間?？梢钥吹?，L0范數約束具有迭代次數少、計算時間短的優點，有利于處理大規模地震勘探數據?？挛鳒蕜t約束只針對一維的地震記錄，因此需要逐道進行反褶積運算，而利用L1范數和L0范數約束進行反褶積運算時，可以直接處理一個二維的地震記錄，各道之間互不干擾。

表1 模型試驗二3種方法的迭代次數與計算時間

圖4 復雜的稀疏反射系數模型(a)及其合成地震記錄(b)

圖5 柯西準則稀疏約束反褶積a 反褶積結果; b 圖5a黑色方框部分的放大顯示

圖6 L1范數稀疏約束反褶積a 反褶積結果; b 圖6a黑色方框部分的放大顯示

圖7 L0 范數稀疏約束反褶積a 反褶積結果; b 圖7a黑色方框部分的放大顯示

圖8 復雜的反射系數模型3種稀疏約束反褶積結果的波形剖面a 柯西準則約束; b L1范數約束; c L0范數約束

3 結束語

我們將L0范數稀疏約束引入地震勘探數據的反褶積處理，并與柯西準則約束和L1范數約束反褶積結果進行對比?？挛鳒蕜t約束反褶積對于較弱的反射系數有壓制作用；L1范數約束反褶積對弱反射系數的保護比柯西準則強，但是求得的反射系數振幅可能會偏弱，且迭代次數較多。L0范數約束反褶積對弱反射系數的保護最好，有利于提高地震勘探資料的精度，并且具有迭代次數少、計算時間短的優點。地震勘探需要處理的數據量很大，因此在進行反褶積運算時以采用L0范數約束可以提高處理效果和計算效率。

本文的模型試驗是對合成地震記錄進行的，子波是已知的。在實際數據的反褶積處理過程中，一個至關重要的問題是如何從含噪地震道中提取準確的高精度地震子波，這將是下一步的研究方向。

參 考 文 獻

[1] 牟永光.地震數據處理方法[M].北京:石油工業出版社,2007:57-79

Mou Y G.Processing methods of seismic data[M].Beijing:Petroleum Industry Press,2007:57-79

[2] 吳常玉,楊瑞娟,鮑崢,等.基于負熵的地震盲反褶積方法及其應用[J].石油物探,2009,48(3):232-238

Wu C Y,Yang R J,Bao Z,et al.Seismic blind deconvolution method based on negative entropy and its application[J].Geophysical Prospecting for Petroleum,2009,48(3):232-238

[3] 鄔世英,孫贊東,朱興卉.動態反褶積中的反射系數序列時頻特征研究[J].石油物探,2011,50(4):324-330

Wu S Y,Sun Z D,Zhu X H.Research on the time-frequency spectrum of reflectivity sequence based on dynamic deconvolution[J].Geophysical Prospecting for Petroleum,2011,50(4):324-330

[4] Taylor H L,Banks S C,Mccoy J F.Deconvolution with the1norm[J].Geophysics,1979,44(1):39-52

[5] Dossal C,Mallat S.Sparse spike deconvolution with minimum scale[C]∥Signal Processing with Adaptive Sparse Structured Representations.Rennes,France:SPARS workshop,2005:1-4

[6] Debeye H,Riel V P.Lp‐norm deconvolution1[J].Geophysical Prospecting,1990,38(4):381-403

[7] 王宇,韓立國,周家雄,等.L1-L2范數聯合約束稀疏脈沖反演的應用[J].地球科學:中國地質大學學報,2009,34(5):835-840

Wang Y,Han L G,Zhou J X,et al.Application of combined norm constrained sparseness spike inverse[J].Earth Science—Journal of China University of Geosciences,2009，34(5):835-840

[8] 劉喜武,劉洪.實現稀疏反褶積的預條件雙共軛梯度法[J].物探化探計算技術,2003,25(3):215-219

Liu X W,Liu H.The preconditional dual conjugate gradient algorithm for sparse deconvolution[J].Computing Techniques for Geophysical and Geochemical Exploration,2003,25(3):215-219

[9] 朱振宇,劉洪.稀疏反褶積方法及其應用[J].石油大學學報:自然科學版,2005,29(6):20-22

Zhu Z Y,Liu H.Sparse deconvolution method and its application[J].Journal of the University of Petroleum,2005,29(6):20-22

[10] 孟小紅,吳何珍,劉國峰.盲源反褶積方法與應用研究[J].石油地球物理勘探,2005,40(6):642-645

Meng X H,Wu H Z,Liu G F.Study of blind deconvolution and application of method[J].Oil Geophysical Prospecting,2005,40(6):642-645

[11] 劉喜武,寧俊瑞,張改蘭.Cauchy稀疏約束Bayesian估計地震盲反褶積框架與算法研究[J].石油物探,2009,48(5):459-464

Liu X W,Ning J R,Zang G L.Cauchy sparse constrained bayesian estimation based seismic blind deconvolution frame and algorithm[J].Geophysical Prospecting for Petroleum,2009,48(5):459-464

[12] Sacchi M D.Reweighting strategies in seismic deconvolution[J].Geophysical Journal International,1997,129(3):651-656

[13] 張繁昌,劉杰,印興耀,等.修正柯西約束地震盲反褶積方法[J].石油地球物理勘探,2008,43(4):391-396

Zang F C,Liu J,Yin X Y,et al.Modified cauchy-constrained seismic blind deconvolution[J].Oil Geophysical Prospecting,2008,43(4):391-396

[14] 焦李成,楊淑媛,劉芳,等.壓縮感知回顧與展望[J].電子學報,2011,39(7):1651-1662

Jiao L C,Yang S Y,Liu F,et al.Development and prospect of compressive sensing[J].Chinese Journal of Electronics,2011,39(7):1651-1662

[15] 曹靜杰,王彥飛,楊長春.地震數據壓縮重構的正則化與零范數稀疏最優化方法[J].地球物理學報,2012,55(2):596-607

Cao J J,Wang Y F,Yang C C.Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization[J].Chinese Journal of Geophysics,2012,55(2):596-607

[16] 陳國新,陳生昌,王漢闖,等.基于 L0 范數最小化的地球物理數據稀疏重構[J].應用地球物理，2013,10(2):181-190

Chen G X,Chen S C,Wang H C,et al.Geophysical data sparse reconstruction based on L0-norm minimization[J].Applied Geophysics,2013,10(2):181-190

[17] Wohlberg B,Rodr guez P.An iteratively reweighted norm algorithm for minimization of total variation functionals[J].Signal Processing Letters,IEEE,2007,14(12):948-951

[18] Chen S S,Donoho D L,Saunders M A.Atomic decomposition by basis pursuit[J].SIAM Journal on Scientific Computing,1998,20(1):33-61

[19] Mallat S G,Zhang Z.Matching pursuits with time-frequency dictionaries[J].IEEE Transactions on Signal Processing,1993,41(12):3397-3415

[20] Tropp J A,Gilbert A C.Signal recovery from random measurements via orthogonal mathing pursuit[J].IEEE Transactions on Information Theory,2007,53(12):4655-4666

[21] Blumensath T,Davies M E.Iterative thresholding for sparse approximations[J].Journal of Fourier Analysis and Applications,2008,14(5):629-654

[22] Beck A,Teboulle M.A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J].SIAM Journal on Imaging Sciences,2009,2(1):183-202