馬彥彥， 周 輝， 張洪禮

( 中國石油大學(北京) 油氣資源與探測國家重點實驗室, 北京 102249)

?

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

馬彥彥， 周 輝， 張洪禮

( 中國石油大學(北京) 油氣資源與探測國家重點實驗室, 北京 102249)

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

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

0 引言

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

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

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

1 方法原理

1.1 疊前正演公式

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

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

(1)

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

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

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

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

d=Gm=WADm

(2)

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

因此地震記錄與模型參數的關系可表示為

(3)

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

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

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

d=Gm+n

(4)

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

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

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

(5)

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

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

(6)

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

(7)

對于地下地層而言，縱向波阻抗是稀疏的，具有明顯的分塊特征。當地層橫向連續(xù)性較好時，可認為相鄰地震道的波阻抗差值接近于“0”；當地層存在斷層、尖滅、裂隙等地質構造時，會導致地震同相軸產生錯斷，在構造邊界處，相鄰地震道的波阻抗差值不再為“0”。因此可認為地層橫向波阻抗差值是稀疏的。為此定義兩個方向性一階差分算子，其中水平(空間)方向記為C1，時間方向記為C2，具體可表示為如下稀疏陣形式：

(8)

(9)

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

(10)。

其中：μ為尺度因子。

最終，后驗概率的表達式為式(11)。

(11)

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

(12)

略去常數項，則最大化后驗概率函數等價于使得式(13)最小。

(13)

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

J=J0+λJm

(14)

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

對式(13)目標函數J求導，并令其等于零，得到

GT(Gm-d)+λCTQCm=0

(15)

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

2 例子

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

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

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

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

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

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

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

圖1 Marmousi2真實模型局部顯示

圖2 初始模型

圖3 無邊緣保護時反演結果

圖4 邊緣保護反演結果

圖5 含20%噪聲時邊緣保護反演結果

3 結論

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

[1] 楊培杰, 印興耀. 非線性二次規(guī)劃貝葉斯疊前反演[J]. 地球物理學報,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三參數同步反演研究[D]. 湖南: 中南大學, 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

國家自然科學基金項目(41174117)

馬彥彥(1983-)，女，博士，主要研究方向為地震資料處理及儲層預測，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