一種地震共反射點道集數(shù)據(jù)的疊前反演方法

張繁昌,印海燕,翁 斌,王保麗

(1.中國石油大學(xué) 地球資源與信息學(xué)院,山東東營 257061;2.中國海洋石油研究總院,北京 100027)

一種地震共反射點道集數(shù)據(jù)的疊前反演方法

張繁昌1,印海燕2,翁 斌2,王保麗1

(1.中國石油大學(xué) 地球資源與信息學(xué)院,山東東營 257061;2.中國海洋石油研究總院,北京 100027)

疊前地震資料含有地層的縱波、橫波速度和密度等信息。利用疊前反演獲得隱藏在地震數(shù)據(jù)中的這些基本參數(shù)后,即可揭示大量巖性及孔隙流體性質(zhì)的信息。這里推導(dǎo)了平面波在層狀彈性介質(zhì)中傳播的正演算子,提出了一種基于B rent正交搜索方向組的疊前三參數(shù)反演方法,該方法不需要求解龐大而復(fù)雜的導(dǎo)數(shù)矩陣。通過自適應(yīng)退火因子和罰函數(shù)來處理約束條件,提高了算法的穩(wěn)定性。將K-L變換引入到方向置換過程,有效防止了搜索方向組的線性相關(guān)。經(jīng)理論模型和油田實際數(shù)據(jù)的反演結(jié)果表明,該反演方法是一種利用疊前地震數(shù)據(jù)進行儲層預(yù)測的有效手段。

疊前反演;地震道集;B rent算法;正演算子;K-L變換

0 前言

縱波的傳播速度隨著巖石中孔隙流體類型的變化而變化,而橫波速度則主要取決于巖石骨架的礦物成份。因此,同時利用縱波、橫波信息,可以有效地進行儲層識別和烴類檢測[1、2]。我們知道,地震波在沿著非法向路徑傳播時,每遇到一個地層分界面,就會產(chǎn)生反射縱波和反射橫波,最后被地面檢波器接收[3]。換句話說,疊前地震數(shù)據(jù)不僅包含有縱波和密度信息,而且攜帶了地層橫波信息。疊前反演的目的,就是利用疊前地震數(shù)據(jù)來得到地下的縱波、橫波等信息[4]。

疊前反演主要有波動方程反演和AVO反演二個方向。疊前波動方程反演在理論上雖然比較成熟,但由于其運算量極大,對噪音敏感,在實際中并沒有得到應(yīng)用。而基于Zoepp ritz方程的疊前道集AVO反演由于運算簡單,抗噪性好而被廣泛應(yīng)用[5]。由于三參數(shù)AVO反演問題是病態(tài)的[6],業(yè)界普遍采用二參數(shù)方程進行反演,以達(dá)到穩(wěn)定反演問題的目的,但同時也引入了較大誤差[7]。另外,一般反演方法都要利用梯度信息,即在反演過程中須求解目標(biāo)函數(shù)的一階導(dǎo)數(shù)矩陣。由于疊前反演數(shù)據(jù)量大,待求解的地層參數(shù)較多,梯度矩陣將占用相當(dāng)大的內(nèi)存空間,從而限制了疊前反演的數(shù)據(jù)規(guī)模。

作者在本文中,將采用一種基于B rent正交搜索方向組[8]的方法,以實現(xiàn)疊前道集的三參數(shù)同時反演。該方法的優(yōu)勢是在反演過程中不需要導(dǎo)數(shù)信息。疊前反演目標(biāo)函數(shù)采用一種退火罰函數(shù)方法構(gòu)建,克服了三參數(shù)反演問題的病態(tài)性質(zhì),提高了該算法的穩(wěn)定性。在反演迭代過程中,將KL變換[9]構(gòu)造的正交方向作為新的搜索方向組,則有效地避免了搜索方向的線性相關(guān),提高了算法的效能。

1 疊前反演問題描述

在非線性疊前地震反演中,實際地震道集dobs與由地質(zhì)模型參數(shù)m合成的道集dmod=A(m),見

式(1)。

其中 矩陣A為正演算子;e為殘差噪音。

可以建立式(2)的疊前反演目標(biāo)函數(shù):

其中 dobsij表示第i次覆蓋、第j個樣點的實際觀測地震數(shù)據(jù);dmodij表示對應(yīng)位置的合成道集數(shù)據(jù)。

實際上,式(2)定義的是預(yù)測誤差向量e的歐幾里的長度平方,可以用任意范數(shù)來求解地層模型參數(shù)。由于合成數(shù)據(jù)與實際數(shù)據(jù)差異很大的奇異點對低階范數(shù)影響不大,而對高階范數(shù)影響甚大,為了減小少數(shù)奇異點對目標(biāo)函數(shù)的影響,將式(2)用L1模代替,如式(3)所示:

2 模型參數(shù)與正演算子的表示

作者在進行常規(guī)地震疊后反演認(rèn)為,地震波是法向入射到界面上的,但事實上常規(guī)疊加道是不同入射角(或不同偏移距)地震記錄的水平疊加,不能真正代表法線入射的地震記錄。所以,疊加已經(jīng)破壞了真實的振幅關(guān)系。而利用疊前道集反演,要基于Zoepp ritz方程的三項線性簡化方程[10、11],通常用地層彈性參數(shù)的變化率表示(見式4):

其中 Rpp(θj)為縱波反射系數(shù),其大小隨入射角θj而變化,產(chǎn)生AVA效應(yīng);Vp、Vs及ρ分別表示地層分界面二側(cè)的局部平均縱波、橫波速度和密度;γ=4(Vs/Vp)2;縱波、橫波速度和密度的相對變化率分別用ΔVp/Vp、ΔVs/Vs、Δρ/ρ表示,這些參數(shù)的相對變化量稱該參數(shù)在地層界面產(chǎn)生的反射系數(shù)。

由式(4)可以看出,非零偏移距地震道的反射振幅,除了與縱波速度及密度有關(guān)外,還含有橫波速度信息,從而得到泊松比等其它地層彈性參數(shù)。記:

將方程(4)重新整理,并寫成矩陣形式(5)。

其中

方程(5)描述的是單個地層界面,在某一入射角或偏移距下的情形。將其擴展到N個地層界面、K個偏移距的一般情形,可以用分塊矩陣表示為式(6)。

式中 向量Rpk是偏移距為k的所有界面的反射系數(shù);rp、rs、rρ分別是縱波、橫波速度和密度反射系數(shù)向量;分塊矩陣Ek、Fk及Gk均為對角矩陣,如

由于帶限疊前地震數(shù)據(jù)是隨偏移距變化的反射系數(shù)和子波的褶積,設(shè)Wk是偏移距為k時的子波矩陣,dk為相應(yīng)的地震數(shù)據(jù),則有式(7)。

作為解釋人員,更加關(guān)心的是地層參數(shù)信息,而不是地層分界面的反射系數(shù)。由于反射系數(shù)和地層參數(shù)的對數(shù)值之間存在微分關(guān)系,如rp=?lnVp/2?t,以差分形式表示即為:

以矩陣形式表示為rp=D lnVp,其中D為差分矩陣:

類似地:

最后得到疊前道集的正演方程,如式(8)。

其中 系數(shù)矩陣即為式(1)中的正演算子A;地層參數(shù)m=[VpVsρ]T;地震道集數(shù)據(jù)為dmod=[d1,d2,…,dK]T。

在地震數(shù)據(jù)有N個時間采樣點,K次覆蓋的情形下,向量m有N×3個元素,向量d mod有N×K個元素,正演算子A為N×K行、N×3列的大型矩陣。

3 罰函數(shù)的構(gòu)建

由于疊前地震資料普遍信噪比較低,且觀測孔徑不足,導(dǎo)致了問題的不適定性。僅依賴合成道集與實際道集的誤差大小得到的地層參數(shù)是非唯一的,甚至不具有實際物理意義和地質(zhì)意義,這也是造成三參數(shù)反演病態(tài)性的主要原因。為了解決這個問題,得到具有地質(zhì)意義的解,所言反演的地層模型參數(shù)Vp、Vs和ρ,除了要使式(3)最小外,還要滿足約束條件:

這樣,引入先驗信息對解空間進行約束,可以使模型參數(shù)在先驗值附近進行搜索,減小尋優(yōu)范圍,加快收斂速度。

在引入先驗信息后,無約束反演問題就變成了有約束反演問題。為了能夠利用B rent方法進行求解,作者設(shè)計了一種退火罰函數(shù)法,將有約束問題轉(zhuǎn)化為無約束問題。設(shè)模型參數(shù)的罰函數(shù)為:

由式(9)可以看出,當(dāng)反演的模型參數(shù)值偏離先驗值太多,大于門檻值η時,該約束項才起作用;而小于η時,不予懲罰。懲罰因子σ吸取模擬退火[12]的思想:令σ=T(k),T(k)是模擬退火算法中的冷卻進度表,其管理方程為T(k)=αT(k-1),其中0<α<1,k為迭代次數(shù)。這樣,在反演初期階段,先驗約束的作用比較大,能夠很好地防止反演參數(shù)的過度調(diào)整而產(chǎn)生的劇烈抖動。隨著反演過程的不斷進行,T逐漸下降,罰函數(shù)的作用逐漸減小,最終得到具有地質(zhì)意義的解。

4 反演問題的求解

在引入先驗信息約束后,疊前反演所要求解的問題為式(10):

對式(10)的反演有多種方法,業(yè)界普遍使用的是基于梯度信息的反演方法[13],在反演過程中需要計算式(10)的一階導(dǎo)數(shù)。由于疊前反演數(shù)據(jù)量大,待求解的地層參數(shù)較多,由一階導(dǎo)數(shù)構(gòu)成的梯度矩陣將占用相當(dāng)大的內(nèi)存。

在不使用導(dǎo)數(shù)信息的反演方法中,B ren t方法是一種相當(dāng)有效的方法,其基本思想是用任意n個線性無關(guān)的方向ξ1、…、ξn作為初始搜索方向組,經(jīng)過n次迭代,陸續(xù)構(gòu)造出n個互相共扼的搜索方向組u1、…、un。在每次迭代中更新搜索方向組時,都保留已經(jīng)得到的共扼方向,而用新構(gòu)造出的方向替換使目標(biāo)函數(shù)值下降最大的方向。

在目標(biāo)函數(shù)有多個極值點的情況下,B ren t方法會出現(xiàn)這樣的問題:在第k次迭代中,由于搜索方向的改變,使u1、…、uk變得線性相關(guān)或接近線性相關(guān)[14],從而嚴(yán)重限制了該方法在多極值反演問題中的應(yīng)用。為避免搜索方向出現(xiàn)線性相關(guān),作者在方向置換策略過程中引入K-L變換,以保證搜索方向線性無關(guān)。

在疊前反演中,設(shè)地層參數(shù)解空間的初始搜索方向組為ui(i=1,...,n),其中u1、…、uk是n維空間中k個關(guān)于A共軛的方向,令:

則有:

(1)置初始迭代次數(shù)k=0,地層參數(shù)的初始估計mk∈Rn,初始搜索方向組ui=ξi∈Rn(i=1,...,n)。

(2)記mk,0=mk為Rn空間中的初始點,按次序i=1、…、n沿各搜索方向?qū)?yōu),即求 βi,使Q(mk,i-1+βiui)滿足式(10),并置:

(3)找到使目標(biāo)函數(shù)下降最快的p個方向,記錄在集合Φ={l1,l2,…,lp}?{1,2,…,n}中,利用K-L變換對下標(biāo)集中的方向組進行置換。

(4)置mk+1=mk,n,若滿足式(10)規(guī)定的門檻值,則mk+1為最優(yōu)解;否則,置k=k+1,返回步驟(2)繼續(xù)迭代。

5 方法應(yīng)用

5.1 方法測試

圖1 無噪聲時的道集對比Fig.1 Gather com parisonw ithou t no ise

模型試驗所采用的疊前道集(見圖1(a)),由實際井曲線根據(jù)Zoepp ritz方程(如下頁圖2中的虛線所示)產(chǎn)生,地震子波為主頻35 Hz的R icker子波,共188m s,1m s采樣,最小偏移距350m,最大偏移距3 200m,共二十次覆蓋。將采用K-L變換進行方向置換的新算法與原算法進行了對比,搜索方向組的線性相關(guān)指標(biāo)(越大表示越線性相關(guān))如下頁圖3所示。由此可見,原算法中搜索方向組隨著迭代次數(shù)的增加,變得越來越線性相關(guān)(如圖3中實線所示),致使算法不收斂。而新的方向置換算法,有效地防止了線性相關(guān)(如圖3中虛線所示)現(xiàn)象。

利用本文中的反演算法,對圖1(a)的疊前道集進行縱橫波速度及密度三參數(shù),同時反演的結(jié)果如下頁圖2中的實線所示。通過與實際井曲線的對比可見,二者基本重合。需說明的是,本反演方法是按照地震數(shù)據(jù)的采樣間隔進行反演的,事先并不知道地層界面的位置。反演結(jié)果合成的疊前地震道集如圖1(b)所示,圖1(c)為前二者的殘差,可見反演結(jié)果的合成道集與真實模型參數(shù)的道集極為相似,殘差基本為零。

為測試地震數(shù)據(jù)中隨機噪聲對反演結(jié)果的影響,在圖1(a)的地震道集中,加入不同程度的隨機噪聲。

后面圖4(a)是加入30%隨機噪聲時的地震道集,用本文中的方法對其進行疊前反演,反演結(jié)果如后面圖5所示。

由此可見,縱波、橫波曲線與實際井曲線能夠很好地吻合,雖然密度曲線的質(zhì)量有所下降,但仍具有較好的可比性。反演結(jié)果的合成道集如圖4(b)所示(見下頁),圖4(c)所示(見下頁)的殘差數(shù)據(jù)基本上是隨機噪聲。

當(dāng)?shù)卣鸬兰械碾S機噪聲達(dá)到50%時(見后面圖6(a)),有效信號已淹沒在噪聲中,后面圖6(b)所示反演結(jié)果的合成道集與真實地層模型的地震道集(如圖1(a)所示)已差別很大,后面圖6(c)的殘差記錄中還隱約殘留有效信號同相軸的影子。圖6(a)地震道集的反演結(jié)果如后面圖7所示。

與前面的反演結(jié)果相比,各參數(shù)與實際井曲線的偏差明顯增大,但從總體上看,縱波速度與原測井曲線仍具有較好的相似性,橫波和密度參數(shù)的反演質(zhì)量下降較大,可見縱波速度較其它二個參數(shù)的反演具有更強的抗噪能力,原因可能是由于縱波對地震AVO響應(yīng)的貢獻(xiàn)最大,這從方程(4)中各參數(shù)前面的系數(shù)不難看出。

5.2 實際數(shù)據(jù)的應(yīng)用

后面的圖8(a)是麗水凹陷過L2井的疊前地震道集,其中在1.87 s附近(實心箭頭所指)為氣層,1.95 s(空心箭頭所指)處為假亮點。L2井已有實測的縱波速度、橫波速度和密度曲線,對該位置進行疊前反演,是為了檢驗方法的有效性,然后再對工區(qū)內(nèi)沒有橫波速度的井位,利用本方法得到相應(yīng)的橫波速度。

圖2 無噪聲時反演結(jié)果與實際井曲線的對比Fig.2 Comparison between inverted and realwell logw ithoutnoise(Dashed line is realwell log,so lid line is inverted resu lt)

圖3 線性相關(guān)指標(biāo)值隨迭代次數(shù)的變化Fig.3 Changesof linear-dependen t ind icato rw ith iteration

圖9(見后面)從左到右是L2井位置處反演的縱波速度、橫波速度及密度曲線,與該井實際測井曲線的對比,其中實線為反演曲線,虛線為L2的實測井曲線。

圖8(b)(見后面)是由反演結(jié)果生成的道集,與圖8(a)的同相軸位置和振幅變化關(guān)系都很一致。

圖4 含較小噪聲時的道集對比Fig.4 Gather comparisonw ith less noise

圖5 圖4(a)所示地震道集的反演結(jié)果與實際井曲線的對比Fig.5 Comparison between inverted resultof Fig 4(a)and realwell log

圖6 含較大噪聲時的道集對比Fig.6 Gather comparisonw ith large noise

圖7 圖6(a)所示地震道集的反演結(jié)果與實際井曲線的對比Fig.7 Comparison between inverted resultof Fig 6(a)and realwell log

見下頁,綜合圖8和圖9可以看出,在1.87 s處,縱波速度降低,密度明顯減小,而橫波速度變化微弱,因此在圖8中實心箭頭所指的同相軸是由儲層含氣所致。而在1.95 s處的縱波、橫波速度均減小,密度值卻較大,故在圖8中空心箭頭所指的同相軸是由巖性變化引起的,并不是由于儲層含氣引起的。

6 結(jié)論

圖8 L2井反演合成道集與實際地震道集的對比Fig.8 Comparison between synthetic and real gatherofwellL2

圖9 L2井三參數(shù)反演結(jié)果與實際井曲線的對比Fig.9 Comparison between 3-term inverted resultand wellL2 real log

疊前地震共反射點道集中含有地層的縱波速度、橫波速度和密度等信息,在利用疊前反演獲得這些參數(shù)后,可揭示巖性及儲層所含流體性質(zhì)的信息。因此,疊前反演已成為巖性油氣藏勘探和烴類檢測的有力手段。

作者在本文中詳細(xì)推導(dǎo)了平面波在層狀彈性介質(zhì)中傳播的正演算子,進而提出了基于B rent正交搜索方向組的疊前三參數(shù)反演方法,該方法不需要求解和存儲復(fù)雜的導(dǎo)數(shù)矩陣。通過自適應(yīng)退火因子和罰函數(shù)來構(gòu)建約束條件,克服了三參數(shù)反演的病態(tài)性質(zhì),提高了算法的穩(wěn)定性。方向置換則利用K-L變換實現(xiàn),有效地防止了搜索方向組的線性相關(guān)現(xiàn)象。經(jīng)理論模型檢驗和實際數(shù)據(jù)的應(yīng)用表明,本方法反演過程穩(wěn)定,效果良好,是一種利用疊前道集數(shù)據(jù)反演地層縱波、橫波及密度信息的有效手段。

[1]ANDERSON P F,GRAY FD.U sing LMR for dual attribute litho logy identification[J].71 th Ann.Internat M tg.,SEG Expanded Abstracts,2001:201.

[2]周中彪.基于巖石物理模型的測井約束橫波速度計算方法研究[J].物探化探計算技術(shù),2010,32(5):536.

[3]印興耀,韓文功,李振春,等.地震技術(shù)新進展(下)[M].山東東營:中國石油大學(xué)出版社,2006.

[4]肖思和,李曙光,許多,等.疊前彈性波阻抗反演在儲層預(yù)測中的應(yīng)用[J].物探化探計算技術(shù),2010,32(5):476.

[5]王永剛.地震資料綜合解釋方法[M].山東東營:中國石油大學(xué)出版社,2007.

[6]DEBSK IW,TARANTOLA A.Inform ation on elastic param eters obtained from the amp litudes of reflected waves[J].Geophysics,1995,60(5):1426.

[7]陳建江,印興耀.基于貝葉斯理論的AVO三參數(shù)波形反演[J].地球物理學(xué)報,2007,50(4):1251.

[8]BRENTR P.A lgorithm form inim izationw ithoutderivatives[M].Englewood C liffs:Prentice-Hall,1973.

[9]張賢達(dá).矩陣分析與應(yīng)用[M].北京:清華大學(xué)出版社,2004.

[10]AK IK,R ICHARDSPG.Quantitative seismo logy:theory andm ethods[M].Cam bridge:W.H.Freem an and Co.,1980.

[11]CONNOLLY P.Elastic impedance[J].The Leading Edge,1999,18(4):438.

[12]PEID,LOU IE JN,SATISH K P.App lication of sim ulated annealing inversion on high-frequency fundam ental-mode Rayleigh wave dispersion curves[J].Geophysics,2007,72(5):R77.

[13]施光燕,董加禮.最優(yōu)化方法[M].北京:高等教育出版社,2005.

[14]高旅端.正交程度及其在B rent方法中的應(yīng)用[J].北京工業(yè)大學(xué)學(xué)報,1997,23(2):42.

P 631.4

A

1001—1749(2011)01—0011—09

國家“863”項目資助(2006AA 09A 102)

2010-06-17 改回日期:2010-11-02

張繁昌(1972-),男,副教授,博士,主要從事地震反演、儲層解釋方面的研究。