基于L0范數(shù)稀疏約束的地震數(shù)據(jù)反褶積

梁東輝，陳生昌

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

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

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

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

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

1 方法原理

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

(1)

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

(2)

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

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

1.1 約束優(yōu)化

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

(3)

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

也可以從貝葉斯定理出發(fā)，給定噪聲的分布概率和反射系數(shù)的先驗分布概率，然后推導(dǎo)出相應(yīng)的稀疏約束形式。貝葉斯定理可以表示為

(4)

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

假設(shè)噪聲滿足高斯分布，且均值為零，標(biāo)準(zhǔn)差為σn，單獨一個噪聲分量的概率密度函數(shù)為

(5)

噪聲整體的分布概率為

(6)

又因為

(7)

所以似然函數(shù)

(8)

(9)

反射系數(shù)總體的分布概率為

(10)

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

(11)

兩邊同時取對數(shù)，得

(12)

使后驗概率分布最大，等價于求解優(yōu)化函數(shù):

(13)

(13)式就是柯西準(zhǔn)則稀疏約束反褶積的優(yōu)化目標(biāo)函數(shù)。

同理，如果我們假設(shè)噪聲滿足高斯分布，均值為0，反射系數(shù)滿足均值為0的拉普拉斯分布，可根據(jù)貝葉斯定理推導(dǎo)出L1范數(shù)稀疏約束反褶積的優(yōu)化目標(biāo)函數(shù):

(14)

1.2 L0范數(shù)約束反褶積

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

(15)

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

(16)

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

醫(yī)學(xué)和藥學(xué)是實踐性非常強(qiáng)的學(xué)科，課堂學(xué)生內(nèi)需要接受的知識繁多枯燥，以教師、課堂、教材為中心的教學(xué)方式無法滿足學(xué)生可持續(xù)發(fā)展的需要。如果沒有課外科技活動作為學(xué)習(xí)內(nèi)容和形式的補(bǔ)充，學(xué)生的學(xué)習(xí)積極性會受到一定的影響。因此，在醫(yī)藥院校中，學(xué)生課外科技活動在人才培養(yǎng)中占有重要地位，它能在校內(nèi)提供學(xué)生理論聯(lián)系實際的機(jī)會，培養(yǎng)學(xué)生的創(chuàng)新思維，提高人才培養(yǎng)質(zhì)量，是教育行業(yè)“供給側(cè)結(jié)構(gòu)性改革”的途徑之一。

(17)

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

(18)

Ha為硬閾值函數(shù):

(19)

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

(20)

每次迭代后降低閾值，再進(jìn)行下一次迭代。

因為實際地震記錄中的噪聲可認(rèn)為是均值為0，標(biāo)準(zhǔn)差為σn的高斯分布噪聲。

(21)

迭代終止的條件可設(shè)定為

(22)

其中，Nd表示觀測值d的個數(shù)。

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

1) 給定反射系數(shù)初值r0為0，估計噪聲的標(biāo)準(zhǔn)差σn，迭代次數(shù)n=0;

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

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

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

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

2 模型試驗

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

圖1 試驗中所使用的子波

2.1 模型試驗一

構(gòu)建一個簡單的稀疏分布反射系數(shù)序列(圖2a)，與圖1中的子波褶積并加入標(biāo)準(zhǔn)差為0.01的高斯噪聲，得到相應(yīng)的合成地震記錄(圖2b)。

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

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

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

2.2 模型試驗二

構(gòu)建一個復(fù)雜的反射系數(shù)模型(圖4a)，與圖1所示子波褶積并加入標(biāo)準(zhǔn)差為0.01的高斯噪聲，得到相應(yīng)的地震記錄(圖4b)。

分別采用柯西準(zhǔn)則約束、L1范數(shù)約束和L0范數(shù)約束進(jìn)行反褶積計算，結(jié)果如圖5a，圖6a和圖7a 所示。紅色表示正的反射系數(shù)，藍(lán)色表示負(fù)的反射系數(shù)。3張圖中黑色方框?qū)?yīng)的放大圖示分別為圖5b，圖6b和圖7b,圖中箭頭所指處有3條較弱的反射系數(shù)同相軸，在圖5b中基本沒有顯示，在圖6b 和圖7b中有顯示，說明柯西準(zhǔn)則約束對弱反射系數(shù)有壓制作用，L1范數(shù)和L0范數(shù)約束能夠較好地保護(hù)弱反射系數(shù)。

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

表1為模型試驗二中3種反褶積方法的迭代次數(shù)與計算時間?？梢钥吹剑琇0范數(shù)約束具有迭代次數(shù)少、計算時間短的優(yōu)點，有利于處理大規(guī)模地震勘探數(shù)據(jù)。柯西準(zhǔn)則約束只針對一維的地震記錄，因此需要逐道進(jìn)行反褶積運(yùn)算，而利用L1范數(shù)和L0范數(shù)約束進(jìn)行反褶積運(yùn)算時，可以直接處理一個二維的地震記錄，各道之間互不干擾。

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

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

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

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

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

圖8 復(fù)雜的反射系數(shù)模型3種稀疏約束反褶積結(jié)果的波形剖面a 柯西準(zhǔn)則約束; b L1范數(shù)約束; c L0范數(shù)約束

3 結(jié)束語

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

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

參 考 文 獻(xiàn)

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

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

[2] 吳常玉,楊瑞娟,鮑崢,等.基于負(fù)熵的地震盲反褶積方法及其應(yīng)用[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] 鄔世英,孫贊東,朱興卉.動態(tài)反褶積中的反射系數(shù)序列時頻特征研究[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范數(shù)聯(lián)合約束稀疏脈沖反演的應(yīng)用[J].地球科學(xué):中國地質(zhì)大學(xué)學(xué)報,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] 劉喜武,劉洪.實現(xiàn)稀疏反褶積的預(yù)條件雙共軛梯度法[J].物探化探計算技術(shù),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] 朱振宇,劉洪.稀疏反褶積方法及其應(yīng)用[J].石油大學(xué)學(xué)報:自然科學(xué)版,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] 孟小紅,吳何珍,劉國峰.盲源反褶積方法與應(yīng)用研究[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].電子學(xué)報,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] 曹靜杰,王彥飛,楊長春.地震數(shù)據(jù)壓縮重構(gòu)的正則化與零范數(shù)稀疏最優(yōu)化方法[J].地球物理學(xué)報,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 范數(shù)最小化的地球物理數(shù)據(jù)稀疏重構(gòu)[J].應(yīng)用地球物理，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