基于Lp擬范數(shù)稀疏約束和交替方向乘子算法的波阻抗反演

張雨強(qiáng),文曉濤,吳 昊,劉 軍,劉 煬

(1.成都理工大學(xué)地球物理學(xué)院,四川成都 610059;2.成都理工大學(xué)油氣藏地質(zhì)及開發(fā)工程重點(diǎn)實(shí)驗(yàn)室,四川成都 610059;3.中國(guó)地質(zhì)大學(xué)構(gòu)造與油氣資源教育部重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430074)

地震數(shù)據(jù)中含有豐富的關(guān)于儲(chǔ)層地層學(xué)、沉積學(xué)和巖石性質(zhì)的信息,地震反演是地質(zhì)解釋中挖掘上述信息的有效途徑[1-2]。根據(jù)反演所利用的基礎(chǔ)數(shù)據(jù),可將反演分為疊后反演和疊前反演。相對(duì)而言,疊前反演信息量更豐富,反演結(jié)果的分辨率更高。但疊后反演也有其自身的優(yōu)勢(shì),如可以相對(duì)快速地將地震波形信息轉(zhuǎn)換為具有地質(zhì)意義的波阻抗信息等,指導(dǎo)儲(chǔ)層的識(shí)別與刻畫[3]。因此,從反演技術(shù)的發(fā)展趨勢(shì)來看,疊后地震反演技術(shù)在油氣勘探中仍然發(fā)揮著重要的作用。

在研究反演的過程中,如何提高反演的精度一直都是國(guó)內(nèi)外學(xué)者重點(diǎn)關(guān)注的問題,比較常用的方法是增加先驗(yàn)信息?；贚2范數(shù)約束的Tikhonov正則化反演方法,一定程度上提高了反演的精度。例如,廣義線性反演[4]和稀疏尖峰反演[5]。但是,使用L2范數(shù)進(jìn)行正則化約束會(huì)導(dǎo)致反演的結(jié)果過于光滑,反射界面不明顯。鑒于反射系數(shù)可看作由大量零值組成的稀疏信號(hào),可利用該稀疏性構(gòu)造稀疏正則項(xiàng),進(jìn)一步對(duì)反演結(jié)果進(jìn)行約束[6]。

稀疏正則化的構(gòu)建,通常是采用L1范數(shù)進(jìn)行正則化約束。在地球物理領(lǐng)域,較早使用L1范數(shù)正則化的是在反褶積[7]和數(shù)據(jù)重構(gòu)[8]等方面。但由于當(dāng)時(shí)沒有合適的算法對(duì)L1范數(shù)約束的優(yōu)化問題進(jìn)行準(zhǔn)確求解,使得L1范數(shù)正則化的反演算法未能得到進(jìn)一步的推廣。近十多年來,壓縮感知理論快速發(fā)展,使得L1范數(shù)正則化再次引起了學(xué)者們的關(guān)注。

ZHANG等[9-10]使用基追蹤算法求解L1范數(shù)約束問題,通過對(duì)疊后地震數(shù)據(jù)反演得到反射系數(shù),為進(jìn)一步增強(qiáng)橫向連續(xù)性,引入“Z型”空間導(dǎo)數(shù)進(jìn)行多道反演。LIU等[11]使用L1范數(shù)正則化反演方法進(jìn)行波阻抗反演,并將該方法與阻尼最小二乘阻抗反演方法對(duì)比,對(duì)比結(jié)果顯示L1范數(shù)正則化反演方法具有更高的反演精度。張繁昌等[12]對(duì)基追蹤反演算法進(jìn)行改進(jìn),將目標(biāo)函數(shù)中保真項(xiàng)的L2范數(shù)約束替換為L(zhǎng)1范數(shù)約束,并使用對(duì)偶對(duì)數(shù)障礙規(guī)劃算法進(jìn)行反演,反演精度相較于傳統(tǒng)基追蹤算法得到進(jìn)一步提高。石戰(zhàn)戰(zhàn)等[13]為了增加在反演中噪聲的魯棒性,在目標(biāo)函數(shù)中引入L1范數(shù)擬合項(xiàng),提出一種L1-L1范數(shù)稀疏約束地震反演方法。ZHOU等[14]提出了一種多道基追蹤疊后反演算法,該算法在L1范數(shù)約束的基礎(chǔ)上引入馬爾可夫過程描述相鄰地震道間的關(guān)系,提高了對(duì)薄層的識(shí)別能力。WU等[15]將L1范數(shù)稀疏約束的拓展形式交疊組稀疏用于波阻抗反演,通過交疊組稀疏的結(jié)構(gòu)特性消除了局部異常值,進(jìn)一步提高反演精度。WANG等[16]提出了一種數(shù)據(jù)驅(qū)動(dòng)的地震波阻抗反演方法,該方法利用L1范數(shù)對(duì)地震空間連續(xù)信息進(jìn)行約束并參與反演,進(jìn)一步提高反演精度。

上述算法大多是基于L1范數(shù)對(duì)地震數(shù)據(jù)的稀疏信息進(jìn)行挖掘,但是L1范數(shù)只是L0范數(shù)的凸松弛函數(shù),得到的稀疏先驗(yàn)信息有限,于是直接使用L0范數(shù)進(jìn)行反演[17-19]。L0范數(shù)雖然是最稀疏的范數(shù),但由于目前主流的求解算法是貪婪迭代策略方法[20],這種方法容易在迭代時(shí)陷入局部極值。最近,有學(xué)者提出Lp擬范數(shù)約束,也就是范數(shù)的p次冪(0

選擇比L1范數(shù)更為稀疏的Lp擬范數(shù)(0

1 基本原理

1.1 正演模型

地震信號(hào)可以由反射系數(shù)與地震子波褶積表示:

S=w*R+N

(1)

式中:S表示地震記錄;R表示反射系數(shù);w表示子波;N為隨機(jī)噪聲;“*”代表卷積運(yùn)算符號(hào)。將卷積運(yùn)算轉(zhuǎn)化為矩陣相乘,得到公式:

S=WR+N

(2)

W的定義為:

(3)

當(dāng)反射系數(shù)較小時(shí),反射系數(shù)與波阻抗之間存在以下關(guān)系:

(4)

式中:Z為波阻抗;i為采樣點(diǎn)序號(hào);j為道數(shù)序號(hào)。假設(shè)L=lnZ,反射系數(shù)與波阻抗之間的關(guān)系表示為:

R=DL

(5)

式中:D為差分矩陣。D定義為:

(6)

結(jié)合以上公式,則公式可以表示為:

S=WDL+N

(7)

1.2 目標(biāo)函數(shù)

在上述正演模型基礎(chǔ)上,傳統(tǒng)的目標(biāo)函數(shù)可以表示為

(8)

為了得到更高精度反演結(jié)果,采用Lp擬范數(shù)對(duì)反射系數(shù)進(jìn)行稀疏約束。Lp擬范數(shù)的具體表達(dá)式為:

(9)

式中:0

(10)

(11)

圖1 反射系數(shù)優(yōu)化問題a L1范數(shù); b Lp范數(shù)

結(jié)合(8)式,并在目標(biāo)函數(shù)中加入初始模型約束,得到新的目標(biāo)函數(shù)為:

(12)

式中:L0為初始模型;μ為初始模型約束參數(shù);λ為正則項(xiàng)的約束參數(shù)。

1.3 反演方法

在上述目標(biāo)函數(shù)中,不僅存在Frobenius范數(shù),而且還存在Lp擬范數(shù)。因此,無(wú)法對(duì)該問題直接求解。在此引入交替方向乘子算法,將目標(biāo)函數(shù)中的不同范數(shù)約束轉(zhuǎn)化為可以直接進(jìn)行求解的多個(gè)子函數(shù),交替進(jìn)行求解。具體過程為:首先,在(12)式中引入拉格朗日乘子項(xiàng),用R來代替DL,即

(13)

引入對(duì)偶項(xiàng)將上述目標(biāo)函數(shù)轉(zhuǎn)化為無(wú)約束優(yōu)化問題:

(14)

式中:C表示R的對(duì)偶項(xiàng);η為拉格朗日乘子的正則化參數(shù)項(xiàng)。

基于交替方向乘子算法(ADMM)流程,將目標(biāo)函數(shù)分解為分別與L,R,C有關(guān)的子函數(shù)。其中,與L有關(guān)的目標(biāo)函數(shù)為:

(15)

(15)式為常規(guī)的優(yōu)化問題,可直接采用凸優(yōu)化算法的求解方法直接得到L的更新方式,即

Li+1=(DTWTWD+μE+ηDTD)-1·
(DTWTS+μL0+ηDTRi-ηDTCi)

(16)

式中:E表示單位矩陣。在得到更新之后的L時(shí),則可以對(duì)R的子函數(shù)進(jìn)行求解,即

(17)

由(17)式中可知,該目標(biāo)函數(shù)為一個(gè)基于Lp擬范數(shù)的優(yōu)化問題,在此引入軟閾值收縮算法進(jìn)行求解[25],得到結(jié)果為:

(18)

其中,sign函數(shù)表示為:

(19)

而關(guān)于對(duì)偶項(xiàng)C的子函數(shù)可以表示為:

(20)

該問題也屬于凸優(yōu)化問題,可利用其求梯度進(jìn)行求解:

Ci+1=Ci+DLi+1-Ri+1

(21)

綜合上述內(nèi)容,得到反演技術(shù)流程框架:

算法:基于Lp稀疏約束和交替方向乘子算法的波阻抗反演輸入:地震記錄S,初始模型L0,地震子波w,參數(shù)項(xiàng)μ,λ,η,p,差分矩陣D,收斂誤差ε,道數(shù)trace輸出:地震波阻抗Z初始化:i=0,t=0,L=L0,R=0,C=01. While t

2 模型測(cè)試與實(shí)際應(yīng)用

2.1 模型測(cè)試

采用Marmousi2模型的部分?jǐn)?shù)據(jù)構(gòu)成理論模型(圖2),對(duì)反演方法進(jìn)行測(cè)試。圖2a為理論阻抗模型的二維剖面,該模型共有500道地震數(shù)據(jù),每道共有450個(gè)采樣點(diǎn),在采樣點(diǎn)200至250,道數(shù)250至450之間存在透鏡體形狀的含氣儲(chǔ)層。利用該波阻抗模型,基于公式(4)得到反射系數(shù)(圖2b),然后將得到的反射系數(shù)與30Hz的Ricker子波進(jìn)行褶積得到地震記錄(圖2c)。圖2d為使用高斯濾波器對(duì)原始波阻抗進(jìn)行低通濾波得到的初始模型。

圖2 理論模型a 理論阻抗模型; b 反射系數(shù); c 合成地震記錄; d 初始模型

為了驗(yàn)證新方法的可行性及優(yōu)勢(shì),使用基于L1范數(shù)約束的基追蹤反演算法與Lp擬范數(shù)約束的反演算法對(duì)上述模型進(jìn)行反演,并比較兩者與理論模型的差異。為了進(jìn)一步驗(yàn)證新方法的抗噪性,對(duì)合成地震記錄分別加入20%和50%的高斯隨機(jī)噪聲后,再次反演。

首先,選擇L1范數(shù)約束與Lp擬范數(shù)約束反演的第350道結(jié)果進(jìn)行對(duì)比。不同噪聲下單道反演結(jié)果對(duì)比如圖3所示,黑色線條為理論阻抗模型;綠色線條為初始模型;紅色線條為基于L1范數(shù)約束的單道反演結(jié)果;藍(lán)色線條為新方法的單道反演結(jié)果。由圖3可見,在不同噪聲情況下的單道反演結(jié)果對(duì)比中,采用新方法得到的反演結(jié)果更接近理論模型,該現(xiàn)象在紅色和藍(lán)色橢圓中較為明顯。圖4,圖5和圖6 分別展示了不同噪聲情況下的反演結(jié)果,可以看出,Lp擬范數(shù)約束反演結(jié)果的橫向連續(xù)性優(yōu)于L1范數(shù)約束反演結(jié)果,在縱向展布上的值也更接近理論模型。為了定量比較,計(jì)算了兩種反演結(jié)果與理論模型的信噪比(signal noise ratio,SNR)與均方誤差(root mean square error,RMSE),具體公式為:

圖3 不同噪聲下單道反演結(jié)果對(duì)比a 無(wú)噪反演結(jié)果單道對(duì)比; b 含20%噪聲單道反演結(jié)果對(duì)比; c 含50%噪聲單道反演結(jié)果對(duì)比

圖4 無(wú)噪聲反演結(jié)果a L1范數(shù)約束反演結(jié)果; b Lp擬范數(shù)約束反演結(jié)果

圖5 含20%噪聲反演結(jié)果a L1范數(shù)約束反演結(jié)果; b Lp擬范數(shù)約束反演結(jié)果

圖6 含50%噪聲反演結(jié)果a L1范數(shù)約束反演結(jié)果; b Lp擬范數(shù)約束反演結(jié)果

(22)

(23)

表1和表2分別為兩種約束條件下反演結(jié)果的信噪比與均方誤差分析結(jié)果。由表1可見,新方法得到的反演結(jié)果的信噪比大于L1范數(shù)得到的反演結(jié)果,而其均方誤差小于L1范數(shù)的反演結(jié)果。由此證明了使用新方法得到的反演結(jié)果優(yōu)于L1范數(shù)的反演結(jié)果。

表1 基于理論模型的信噪比分析結(jié)果

表2 基于理論模型的均方誤差分析結(jié)果

2.2 實(shí)際應(yīng)用

研究區(qū)位于四川盆地南緣,目的層為志留系龍馬溪組,巖性為頁(yè)巖。兩口鉆井(A井和B井)均鉆遇優(yōu)質(zhì)頁(yè)巖氣儲(chǔ)層,但儲(chǔ)層厚度較小(小于10m),常規(guī)的阻抗反演方法由于受到精度和分辨率的限制,很難對(duì)頁(yè)巖氣儲(chǔ)層進(jìn)行有效識(shí)別。為此,選擇該區(qū)域作為試驗(yàn)區(qū),驗(yàn)證新方法的有效性。圖7為過A井和B井的連井剖面。本次反演中,在地震層位信息的指導(dǎo)下,選擇A井波阻抗數(shù)據(jù)進(jìn)行外推插值并做低通濾波處理得到低頻模型,選擇未進(jìn)行反演的B井對(duì)反演結(jié)果進(jìn)行驗(yàn)證。分別使用基于L1范數(shù)約束的反演方法與新方法進(jìn)行反演,得到的結(jié)果如圖8和圖9所示。

圖7 工區(qū)地震剖面

由圖8和圖9可以看到,兩種方法反演的縱波阻抗與A、B兩井測(cè)井曲線吻合較好,表明兩種方法在野外資料應(yīng)用中可行且有效。與基于L1范數(shù)約束反演方法的縱波阻抗剖面相比,新方法能夠反演出更

圖8 基于L1范數(shù)約束反演結(jié)果

多地層信息(圖8和圖9中的黃色圓圈)。同時(shí),新方法獲得的阻抗信息具有更好的橫向連續(xù)性(圖8和圖9 中的紅色圓圈)。為了進(jìn)一步說明新方法的優(yōu)越性,利用兩種方法的反演結(jié)果重新合成地震剖面(圖10 和圖11),再分別與實(shí)際地震剖面進(jìn)行殘差分析(圖12和圖13)。由圖12和圖13可見,新方法的合成地震剖面與實(shí)際地震剖面殘差更小,即新方法反演結(jié)果與實(shí)際數(shù)據(jù)更吻合。

圖9 基于Lp擬范數(shù)約束反演結(jié)果

圖10 基于L1范數(shù)約束反演結(jié)果合成的地震剖面

圖11 基于Lp擬范數(shù)約束反演結(jié)果合成的地震剖面

圖12 基于L1范數(shù)約束反演結(jié)果合成的地震剖面與實(shí)際地震剖面的殘差

圖13 基于Lp擬范數(shù)約束反演結(jié)果合成的地震剖面與實(shí)際地震剖面的殘差

3 結(jié)論與展望

為了提高反演結(jié)果的精度,提出了一種基于Lp擬范數(shù)稀疏約束和交替方向乘子算法的波阻抗反演方法,理論模型測(cè)試以及實(shí)際數(shù)據(jù)應(yīng)用,結(jié)果表明:

1) 針對(duì)傳統(tǒng)L1范數(shù)挖掘地震數(shù)據(jù)的稀疏先驗(yàn)信息有限的問題,引入了更為稀疏的Lp擬范數(shù)(0

2) 針對(duì)基于Lp擬范數(shù)稀疏約束的非凸優(yōu)化問題,采用交替方向乘子算法能夠有效地進(jìn)行求解。

新方法反演結(jié)果的精度相較于常規(guī)方法有較大提升,但由于算法是針對(duì)單道數(shù)據(jù)的反演,反演結(jié)果的橫向連續(xù)性有待提高。下一步考慮引入全變分正則化,對(duì)目標(biāo)函數(shù)進(jìn)行橫向約束,實(shí)現(xiàn)多道同時(shí)反演,以提高反演結(jié)果的橫向連續(xù)性。