馬彥彥， 周 輝， 張洪禮

( 中國(guó)石油大學(xué)(北京) 油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249)

?

基于時(shí)間和空間方向差分稀疏約束的疊前AVO反演

馬彥彥， 周 輝， 張洪禮

( 中國(guó)石油大學(xué)(北京) 油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249)

在疊前AVO反演中，常規(guī)的吉洪諾夫正則化方法不能突出異常體在縱向上的邊界，并且忽略了地層的橫向聯(lián)系，導(dǎo)致反演結(jié)果產(chǎn)生邊緣模糊現(xiàn)象。為了解決這一問題，這里基于Zeoppritz方程的Fatti近似公式和貝葉斯理論框架，以模型參數(shù)柯西分布作為先驗(yàn)約束，通過引入時(shí)間和空間方向的差分稀疏約束保護(hù)地層邊界，并最終采用迭代重加權(quán)最小二乘算法實(shí)現(xiàn)縱、橫波阻抗的同時(shí)多道反演。模型反演結(jié)果具有明顯的地層邊界和分層特征，表明該方法能夠解決邊緣模糊問題、保護(hù)不連續(xù)體的橫向邊界，大大提高了油氣預(yù)測(cè)的準(zhǔn)確度。

疊前AVO反演； 差分稀疏約束； 邊界； 同時(shí)多道反演

0 引言

地震疊前AVO反演方法能夠利用不同入射角度疊加數(shù)據(jù)和測(cè)井?dāng)?shù)據(jù)，得到直接反映儲(chǔ)層巖性、物性參數(shù)的多種彈性參數(shù)，提高預(yù)測(cè)精度，日漸成為儲(chǔ)層預(yù)測(cè)的主流方法之一。然而由于地震數(shù)據(jù)頻帶的限制，以及噪聲、正演近似等因素的影響，導(dǎo)致該反演方法高度不適定[1]。許多學(xué)者通過在目標(biāo)函數(shù)中引入先驗(yàn)信息作為正則化約束，例如吉洪諾夫正則化方法，以提高反演問題的穩(wěn)定性。但是這種正則化方法常常會(huì)模糊地層巖性發(fā)生變化的位置和地質(zhì)構(gòu)造的邊界等。隨著石油勘探開發(fā)逐漸向隱蔽油氣藏，特別是向砂巖油氣藏方向轉(zhuǎn)變，迫切需要準(zhǔn)確刻畫砂體的空間分布特征。因此研究基于邊緣保護(hù)的反演方法能夠更好地識(shí)別巖性，描述砂體分布。

邊緣保護(hù)的思想最早起源于圖像處理領(lǐng)域，主要采用全變分作為正則化約束項(xiàng)[2-3]。近年來，許多學(xué)者開始將該思想引入地球物理領(lǐng)域，采用邊緣保護(hù)方法準(zhǔn)確描述地下地質(zhì)構(gòu)造。Valenciano等[4]在柯西約束下運(yùn)用邊緣保護(hù)正則化方法進(jìn)行層速度估計(jì)；Youzwishen等[5]在疊后地震資料中運(yùn)用邊緣保護(hù)正則化算法來估計(jì)縱波速度；Zhang等[6]提出一種自適應(yīng)波阻抗模型重建方法，通過直接引入先驗(yàn)波阻抗信息和聲波測(cè)井信息作為絕對(duì)約束，間接得到能夠表征地層、斷層等地質(zhì)構(gòu)造特征的相對(duì)約束的方式，解決疊后波阻抗反演的不適定及帶寬有限問題；張賽民[7]應(yīng)用似柯西分布函數(shù)作為正則化項(xiàng)，開展了基于邊界保持塊約束的疊后波阻抗反演和疊前AVA三參數(shù)同步反演研究；Anagaw等[8]提出一種基于邊緣保護(hù)的地震成像方法，能夠同時(shí)檢測(cè)地層結(jié)構(gòu)的縱向、橫向不連續(xù)問題；Guitton[9]通過在全波形反演中引入塊狀正則化約束，獲得塊狀速度模型，消除照明和反演假象，由于對(duì)速度模型采用水平導(dǎo)數(shù)稀疏約束，該方法還可以衰減高頻噪聲；Yuan等[10]基于反演的思想，以柯西分布作為先驗(yàn)分布，在貝葉斯理論框架下研究了邊緣保護(hù)隨機(jī)噪聲衰減方法，取得了良好的效果；Tian等[11]以馬爾科夫隨機(jī)場(chǎng)的吉布斯分布形式作為先驗(yàn)約束，用Huber函數(shù)作為吉布斯分布的能量函數(shù)，在地層邊緣得到保護(hù)的前提下反演縱波速度、橫波速度和密度三個(gè)參數(shù)；Yuan等[12]提出了保構(gòu)造的三維同時(shí)多道波阻抗反演方法，能夠穩(wěn)定地震反演，減小高波數(shù)噪聲對(duì)反演結(jié)果的影響，探測(cè)地質(zhì)構(gòu)造的空間連續(xù)性，提升了現(xiàn)有地震波阻抗反演方法的精度。

常規(guī)的在疊前反演目標(biāo)函數(shù)中引入的吉洪諾夫正則化約束存在一定的局限性：①該方法為時(shí)間方向全局光滑，不能突出異常體在縱向上的邊界；②該方法忽略了地層結(jié)構(gòu)的橫向聯(lián)系；③會(huì)受高波數(shù)噪聲的影響。為了解決邊緣模糊問題，同時(shí)保持抗噪作用，作者基于Zeoppritz方程的Fatti近似公式，從貝葉斯理論框架出發(fā)，以柯西分布作為先驗(yàn)約束，基于邊界保持的思想，同時(shí)引入時(shí)間和空間方向差分稀疏約束，反演縱、橫波阻抗，為儲(chǔ)層預(yù)測(cè)和流體識(shí)別提供可靠數(shù)據(jù)。

1 方法原理

1.1 疊前正演公式

Zoeppritz方程的Fatti近似方程式[13]可表示為

R(θ)=(1+tan2θ)Rp-8γ2sin2θRs- (tan2θ-4γ2sin2θ)Rρ

(1)

根據(jù)Walker等的研究[14]，在反射系數(shù)小于0.3的情況下，反射系數(shù)可以用波阻抗對(duì)數(shù)表示。令模型參數(shù)

m=[lnZp,lnZs,lnρ]T，

其中：Zp、Zs、ρ分別為縱波阻抗、橫波阻抗和密度向量。

根據(jù)Buland(2003)的研究[15]，可將正演地震記錄的褶積過程表示為式(2)。

d=Gm=WADm

(2)

其中：W為子波褶積矩陣，A為Fatti近似方程式(1)中與角度有關(guān)的系數(shù)矩陣，D為縱向上的一階差分矩陣。

因此地震記錄與模型參數(shù)的關(guān)系可表示為

(3)

以疊前三參數(shù)同步反演為例，假定待反演的參數(shù)維度為N×M(M為道數(shù)，N為每道采樣點(diǎn)數(shù))，如果參與反演的角度道集個(gè)數(shù)為K，則式(2)中，觀測(cè)數(shù)據(jù)d的維度為KN×M，正演算子矩陣G的維度為KN×3N，模型參數(shù)矩陣m的維度為3N×M。

1.2 基于時(shí)間和空間方向差分稀疏約束的疊前AVO反演

由前文可知，含噪地震記錄可表示為式(4)。

d=Gm+n

(4)

其中：d=[d1,d2,…,dK]T為觀測(cè)地震數(shù)據(jù)；G=WAD為正演算子；m=[lnZp,lnZs,lnρ]T為包含縱波阻抗、橫波阻抗及密度三個(gè)分量的模型參數(shù)矩陣；n為觀測(cè)噪聲。

基于隨機(jī)變量表達(dá)形式的貝葉斯公式可表示為式(5)。

p(m|d)∝p(m)p(d|m)

(5)

其中:p(m)表示模型參數(shù)m的先驗(yàn)概率分布；后驗(yàn)概率分布p(m|d)描述了對(duì)于給定觀測(cè)數(shù)據(jù)d，模型參數(shù)的概率；p(d|m)為似然函數(shù)，描述觀測(cè)數(shù)據(jù)與理論模型之間的差異。

假定噪聲服從高斯分布，方差為σ2。若噪聲不相關(guān)且均值為零，并假設(shè)每種噪聲樣本的方差為常數(shù)，則噪聲向量的總概率為式(6)。

(6)

似然函數(shù)可以用噪聲的概率分布表示，又n=d-Gm，則似然函數(shù)為式(7)。

(7)

對(duì)于地下地層而言，縱向波阻抗是稀疏的，具有明顯的分塊特征。當(dāng)?shù)貙訖M向連續(xù)性較好時(shí)，可認(rèn)為相鄰地震道的波阻抗差值接近于“0”；當(dāng)?shù)貙哟嬖跀鄬?、尖滅、裂隙等地質(zhì)構(gòu)造時(shí)，會(huì)導(dǎo)致地震同相軸產(chǎn)生錯(cuò)斷，在構(gòu)造邊界處，相鄰地震道的波阻抗差值不再為“0”。因此可認(rèn)為地層橫向波阻抗差值是稀疏的。為此定義兩個(gè)方向性一階差分算子，其中水平(空間)方向記為C1，時(shí)間方向記為C2，具體可表示為如下稀疏陣形式：

(8)

(9)

假設(shè)模型參數(shù)m服從柯西分布，則其先驗(yàn)概率為式(10)。

(10)。

其中：μ為尺度因子。

最終，后驗(yàn)概率的表達(dá)式為式(11)。

(11)

給定后驗(yàn)概率分布p(m|d)，求取模型參數(shù)最常用的方式為尋找使得p(m|d)最優(yōu)的最大后驗(yàn)概率估計(jì)，方程(5)改寫為式(12)。

(12)

略去常數(shù)項(xiàng)，則最大化后驗(yàn)概率函數(shù)等價(jià)于使得式(13)最小。

(13)

其中：λ=2σ2。式(13)可表示為以下一般形式：

J=J0+λJm

(14)

其中：第一項(xiàng)J0反映了模型參數(shù)與正演算子合成地震記錄與實(shí)際地震記錄的相近程度；第二項(xiàng)Jm為模型參數(shù)約束項(xiàng)，基于柯西分布先驗(yàn)約束，通過引入橫向約束和縱向約束，起到了邊緣保護(hù)的作用。λ為加權(quán)因子，控制模型參數(shù)稀疏特征對(duì)反演結(jié)果的貢獻(xiàn)。

對(duì)式(13)目標(biāo)函數(shù)J求導(dǎo)，并令其等于零，得到

GT(Gm-d)+λCTQCm=0

(15)

顯然，式(15)屬于非線性求解問題，Youzwishen指出,對(duì)于這類問題可以采用迭代重加權(quán)最小二乘算法求解[16]。

2 例子

根據(jù)式(3)所示地震記錄褶積過程及式(13)所示目標(biāo)函數(shù)，這里采用Marmousi2局部模型進(jìn)行方法測(cè)試。如圖1所示，該局部模型地質(zhì)構(gòu)造復(fù)雜，速度變化劇烈，縱向上分塊性較明顯，發(fā)育許多薄層，橫向上具有明顯的地質(zhì)構(gòu)造邊界，發(fā)育一條傾斜生長(zhǎng)斷層(黑色箭頭)以及一個(gè)不整合面(紅色箭頭)，存在地層尖滅(橢圓區(qū)域)等地質(zhì)現(xiàn)象，有利于檢驗(yàn)本文方法的邊緣保護(hù)作用。

地震觀測(cè)數(shù)據(jù)采用合成記錄方式得到。首先利用Fatti近似公式得到反射系數(shù)，然后將其與地震子波褶積，得到角度域合成記錄。采用的地震子波為主頻35Hz的零相位雷克子波，共生成5°、10°、15°、20°、25°五個(gè)角度的共角度域道集。

反演參數(shù)的低頻信息對(duì)反演結(jié)果具有重要影響，有利于解決反演多解性問題，文本將低頻信息作為初始模型代入求解公式。對(duì)圖1所示真實(shí)模型進(jìn)行0Hz～10Hz的帶通濾波，得到待反演參數(shù)的低頻初始模型如圖2所示。

基于稀疏柯西分布先驗(yàn)信息約束，通過在反演目標(biāo)函數(shù)引入橫向空間約束C1和縱向時(shí)間約束C2，保護(hù)地層邊界信息。在計(jì)算過程中，涉及到加權(quán)因子的選取問題，該值的大小主要通過試驗(yàn)的方法求取。

圖3為不加約束項(xiàng)，即約束項(xiàng)權(quán)重為λ=0時(shí)的縱波阻抗、橫波阻抗反演結(jié)果。與圖1對(duì)比可以看出，反演結(jié)果的邊界模糊，觀察不到地層尖滅特征，橫波阻抗反演結(jié)果剖面上幾乎無法識(shí)別斷層的斷點(diǎn)位置(黑色箭頭所指)。此外，對(duì)模型下方500ms以下的傾斜地層而言，未加邊緣保護(hù)約束項(xiàng)時(shí)，反演方法幾乎完全不能恢復(fù)其地層形態(tài)。

圖4是選取最優(yōu)邊緣保護(hù)約束項(xiàng)權(quán)重時(shí)的縱波阻抗、橫波阻抗反演結(jié)果。與圖1、圖3進(jìn)行對(duì)比可以看出，反演結(jié)果得到很大改善，橫向上的斷層邊界(黑色箭頭所指)及地層尖滅(黑色橢圓所示)、縱向上的地層界面刻畫地更加清楚，薄層(枚紅色箭頭所指)也得到很好恢復(fù)。

為了驗(yàn)證方法的抗噪性，對(duì)含噪合成記錄進(jìn)行反演測(cè)試。在合成的五個(gè)角道集地震記錄中加入全頻帶隨機(jī)噪聲，研究約束項(xiàng)對(duì)于噪聲壓制及邊緣保護(hù)的作用。加入的噪聲占地震記錄能量的20%。圖5為最終反演結(jié)果，盡管較圖4反演結(jié)果分辨率有所降低，但依然較好地保持了構(gòu)造邊界，地層尖滅、薄層等細(xì)節(jié)得到了較好恢復(fù)。

圖1 Marmousi2真實(shí)模型局部顯示

圖2 初始模型

圖3 無邊緣保護(hù)時(shí)反演結(jié)果

圖4 邊緣保護(hù)反演結(jié)果

圖5 含20%噪聲時(shí)邊緣保護(hù)反演結(jié)果

3 結(jié)論

作者提出了一種基于時(shí)間和空間方向的差分稀疏約束疊前AVO反演方法。該方法以模型參數(shù)柯西分布作為先驗(yàn)信息，通過引入縱向上的塊約束和橫向上相鄰道之間的約束，實(shí)現(xiàn)同時(shí)多道疊前阻抗反演，解決常規(guī)正則化約束疊前方法帶來的邊緣模糊問題，能夠保護(hù)地層的縱向邊界，描述地層的橫向連續(xù)性，刻畫地質(zhì)構(gòu)造、地層巖性在橫向上的邊界。對(duì)模型資料進(jìn)行反演處理，驗(yàn)證了該方法的反演效果和穩(wěn)健性。

[1] 楊培杰, 印興耀. 非線性二次規(guī)劃貝葉斯疊前反演[J]. 地球物理學(xué)報(bào),2008, 51(6):1876-1882.YANGPJ,YINXY.Non-linearquadraticprogrammingBayesianprestackinversion[J].ChineseJournalofGeophysics, 2008, 51(6): 1876-1882.(InChinese)

[2]GEMAND,REYNOLDSG.Constrainedrestorationandtherecoveryofdiscontinuities[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 1992, 14(3): 367-383.

[3]BOUMANC,SAUERK.AgeneralizedGaussianimagemodelforedge-preservingMAPestimation[J].IEEETransactionsonImageProcessing, 1993, 2(7): 296-310.

[4]VALENCIANOAA,BROWNM,GUITTONA,etal.Intervalvelocityestimationusingedge-preservingregularization[C].The77thAnnualInternationalMeeting,SEG,ExpandedAbstracts, 2004, 2431-2434.

[5]YOUZWISHENCF,SACCHIMD.Edgepreservingimaging[J].JournalofSeismicExploration, 2006, 15: 45-58.

[6]ZHANGHB,SHANGZP,YANGCC.Adaptivereconstructionmethodofimpedancemodelwithabsoluteandrelativeconstraints[J].JournalofAppliedGeophysics, 2009, 67: 114-124

[7] 張賽民. 基于邊界保持塊約束疊后波阻抗及疊前AVA三參數(shù)同步反演研究[D]. 湖南: 中南大學(xué), 2011.ZHANGSM.Researchofpost-stackimpedanceinversionandAVAThree-termSimultaneousInversionbasedonedgepreservingblockyconstrain[D].Hunan:CentralSouthUniversity, 2011.(InChinese)

[8]ANAGAWAY,SACCHIMD.Edge-preservingseismicimagingusingthetotalvariationmethod[J].JournalofGeophysicsandEngineering, 2012(9): 138-146.

[9]GUITTONA.Blockyregularizationschemesforfull-waveforminversion[J].GeophysicalProspecting, 2012, 60: 870-884.

[10]YUANSY,WANGSX.Edge-preservingnoisereductionbasedonBayesianinversionwithdirectionaldifferenceconstraints[J].Journalofgeophysicsandengineering, 2013, 10(2): 25001-25010.

[11]TIANYK,ZHOUH,CHENHM,etal.Bayesianprestackseismicinversionwithaself-adaptiveHuber-Markovrandom-fieldedgeprotectionscheme[J].AppliedGeophysics, 2013, 10(4): 453-460.

[12]YUANSY,WANGSX,LUOCM,etal.Simultaneousmultitraceimpedanceinversionwithtransform-domainsparsitypromotion[J].Geophysics, 2015, 80(2): 71-80.

[13]FATTIJL,SMITHGC,VAILPJ,etal.DetectionofgasinsandstonereservoirsusingAVOanalysis:a3DSeismicCaseHistoryUsingtheGeostackTechnique[J].Geophysics, 1994, 59(9): 1362-1376.

[14]WALKERC,ULRYCHTJ.Autoregressiverecoveryoftheacousticimpedance[J].Geophysics, 1983, 48(10): 1338-1350.

[15]BULANDA,OMREH.BayesianlinearizedAVOinversion[J].Geophysics, 2003, 68(1): 185-198.

[16]YOUZWISHENCF.Non-linearsparseandblockconstraintsforseismicinversionproblems[D].Alberta:UniversityofAlberta, 2001.

Prestack AVO inversion with space and time directional difference sparse constraints

MA Yan-yan, ZHOU Hui, ZHANG Hong-li

(China University of Petroleum, State Key Laboratory of Petroleum Resource and Prospecting, Beijing 102249, China)

In seismic prestack AVO inversion, the conventional Tikhonov regularization method can not highlight the vertical edge of geologic anomalies and usually ignores the lateral connection of strata, which will cause fuzzy edges in inversion results. To address this issue, based on the Fatti approximate formula of Zeoppritz equation and Bayesian framework, this paper introduces the space and time directional difference sparse constraints expressed by the Cauchy norm to protect edges, and realizes the simultaneous multitrace inversion of P-impedance and S-impedance with iteratively reweighted least squares algorithm. The inversion results for model have obvious stratigraphic edge and layered characteristics, and it indicates that the new method can solve the edge fuzzy problem, protect the lateral edge of discontinuities, and largely improve the accuracy of predicting reservoirs.

prestack AVO inversion; difference sparse constraint; edge; simultaneous multitrace inversion

2015-07-14改回日期：2015-08-03

國(guó)家自然科學(xué)基金項(xiàng)目(41174117)

馬彥彥(1983-)，女，博士，主要研究方向?yàn)榈卣鹳Y料處理及儲(chǔ)層預(yù)測(cè)，E-mail：cup_mayanyan@163.com。

1001-1749(2015)05-0616-06

P 631.4

A

10.3969/j.issn.1001-1749.2015.05.12