印興耀,劉曉晶,吳國忱,宗兆云

(中國石油大學(華東)地球科學與技術學院,山東青島266580)

?

模型約束基追蹤反演方法

印興耀,劉曉晶,吳國忱,宗兆云

(中國石油大學(華東)地球科學與技術學院,山東青島266580)

摘要:基于反射系數奇偶分解的基追蹤反演方法,補充了地震資料中所缺失的低頻與高頻信息,較好地提高了反演結果對地層的分辨能力。但僅僅使用稀疏約束加入的低頻信息缺乏合理性,可能與工區(qū)的實際地質情況不符,需要進一步改善反演結果的橫向連續(xù)性。因此,提出在基追蹤反演目標函數中加入模型約束,得到模型約束的基追蹤反演目標函數,并使用梯度投影稀疏重構(Gradient Projection for Sparse Reconstruction,GPSR) 算法進行求解。模型約束的加入增強了反演的穩(wěn)定性,使得反演結果中的低頻信息更加符合工區(qū)實際地質背景信息,并且能夠改善反演結果的橫向連續(xù)性。楔形模型和實際數據測試結果表明,模型約束基追蹤反演方法不僅保持了基追蹤反演的稀疏性,地層阻抗呈現塊化,反射界面刻畫清晰,而且反演方法更為穩(wěn)定,反演結果的橫向連續(xù)性得到了改善,從而驗證了該方法的可行性。

關鍵詞：模型約束;基追蹤;地震反演;稀疏性;橫向連續(xù)性

在線性系統(tǒng)假設前提下,ROBINSON[1]提出利用褶積模型描述地震響應,即地震記錄可近似看作是地震子波與地層反射系數的褶積結果。自從TAYLOR等[2]提出確定性反褶積方法以來,基于褶積模型的反演方法成為了獲取儲層參數的重要方法。由于地震反演通常是一個病態(tài)問題,反演不穩(wěn)定,因此需要加入一定的先驗約束將不適定問題轉化為近似適定問題。常規(guī)的反演方法通常是從地層的反射系數出發(fā),假設反射系數是稀疏的,在稀疏約束條件下進行求解。反射系數的稀疏約束條件通常有Lp(p=0,1,2)范數[3],其中,基于L1范數的稀疏脈沖反演方法在儲層預測中得到了廣泛應用;或者假設反射系數服從一定的先驗概率分布(如:HUBER分布[4]、柯西分布[5-6]、改進柯西分布[7-8]),通過貝葉斯理論在后驗概率最大條件下建立目標函數進行反演,DUIJNDAM[9]和TARANTOLA[10]均對基于貝葉斯理論的反演方法進行了詳細闡述。上述方法在增強反演方法穩(wěn)定性與提高反演結果分辨率方面有一定的作用。隨著油氣勘探程度的不斷深入,對反演結果的準確性、合理性以及分辨率提出了更高的要求。常規(guī)反演結果縱向上是平滑的,存在子波旁瓣效應,制約了我們利用反演結果對地層邊界的識別。隨著信號稀疏表示技術的快速發(fā)展,基于稀疏表示理論的反演方法也逐漸被應用到地震反演中,在確?？尚哦鹊幕A上,利用稀疏表示理論反演,將地層的反演結果塊化顯示,反演結果對地層的分辨能力將會得到增強。THEUNE等[11]分析認為塊化反演結果的分辨率更高,對地層的解釋能力更強。ZHANG等[12-13]采納了THEUNE等的觀點，基于反射系數奇偶分解提出了基追蹤反演方法,并應用于疊后與疊前地震資料,得到了塊化的、更有利于地質解釋的反演結果。針對基追蹤反演橫向連續(xù)性差的問題,ZHANG等[14]提出了基于空間正則化的多道基追蹤反演方法,改善了反演結果的橫向連續(xù)性。

由于基追蹤反演中使用了L1范數的稀疏約束,稀疏約束能夠補充地震資料缺少的低頻與高頻成分。但是,這樣的低頻成分缺乏實際工區(qū)的背景信息,并且基追蹤反演方法的穩(wěn)定性需要進一步加強?；谀Ｐ图s束的地震反演方法能夠較好地解決這兩個問題。YIN等[15]和ZONG等[16-18]分別利用模型點約束與平滑模型約束增強反演的穩(wěn)定性,并使得反演結果中的低頻信息與實際工區(qū)地質背景相一致,從而使反演結果更為合理。為了增強基追蹤反演方法的穩(wěn)定性,改善其橫向連續(xù)性,本文采用在Zhang等提出的基追蹤反演的基礎上,采納Zong等提出的平滑模型約束的思想,在基追蹤反演目標函數中加入模型約束性,在模型約束與稀疏約束的聯合約束下,利用梯度投影稀疏重構(Gradient Projection for Sparse Reconstruction,GPSR)算法[19]進行求解。在確保反演結果稀疏性的基礎上,增強反演的穩(wěn)定性,改善反演結果的橫向連續(xù)性,使得反演結果更為合理。

1反射系數奇偶分解

(1)

式中:re(i,j,Δt)為偶分量;ro(i,j,Δt)為奇分量;Δt為采樣間隔;i為地層頂界面所對應的采樣點位置；j為地層底界面所對應的采樣點位置;a,b分別表示奇偶分解后偶分量與奇分量的系數。其中,奇偶反射系數對可以分別寫作:

(2)

(3)

圖1　稀疏脈沖反演與稀疏層反演研究對象a 稀疏脈沖反演研究對象; b 稀疏層反演研究對象

圖2　反射系數對分量分解(根據Zhang等[13]修改)

ZHANG等[13]在應用反射系數奇偶分解的基礎上考慮地層可能的厚度,構成奇偶反射對的楔形字典,并將單道反射系數序列用奇偶分量楔形字典表示:

(4)

式中:r為反射系數序列;ai,j為偶分量楔形字典的系數序列;bi,j為奇分量楔形字典的系數序列;I為單道采樣點個數;J為地層最大可能厚度對應的采樣點個數;D=[rero]為奇偶分量構成的楔形字典;m=[aTbT]T為楔形字典中對應的系數。

2模型約束基追蹤反演

假設地層為水平層狀介質,每一層為均勻各向同性完全彈性介質,根據Robinson褶積模型,地震資料中每一道地震數據由地震子波與地層反射系數序列褶積得到,考慮有噪聲的情形,寫成矩陣形式:

d=Wr+n

(5)

式中:d表示地震記錄;W為子波矩陣;r為反射系數序列;n為噪聲序列。

將(4)式代入(5)式可以得到:

d=WDm+n=Gm+n

(6)

式中:G為子波矩陣與楔形字典的乘積,G=WD。顯然,矩陣G的列數遠大于其行數,為欠定矩陣。根據稀疏表示理論,m為反射系數在矩陣G中的稀疏表示系數,即待反演參數。

ZHANG等[13]基于上述反射系數的奇偶分解,考慮地層的稀疏性,使用L1范數(即‖m‖1)對反演目標函數進行稀疏性約束,得到基于基追蹤反演的目標函數:

(7)

在L1范數稀疏約束下,能夠補充地震頻帶所缺失的低頻與高頻信息(圖3),提高地震反演的分辨率。然而,PéREZ等[20]認為這些信息是在數學假設的條件下獲取的,缺乏實際工區(qū)的地質背景信息,并且ZHANG等[14]指出(7)式的反演結果缺乏橫向連續(xù)性。利用測井資料建立低頻模型對(7)式進行約束,可以有效補充低頻信息,由于低頻信息來自于實際工區(qū)的測井資料,因此,最終反演結果包含了工區(qū)內更多的真實背景,反演結果更為合理,并且模型約束能夠有效地改善反演結果橫向的連續(xù)性。

圖3　地震反演補充低頻與高頻信息

在t時刻所對應的采樣點上,縱波阻抗與反射系數之間的關系為:

(8)

將(8)式寫作離散化的矩陣-向量形式:

(9)

式中:C為積分矩陣;ξ為相對波阻抗序列。其表達式分別為:

(10)

(11)

因為反射系數在楔形字典分解如(4)式所示,所以(9)式可被寫為:

CDm=ξ

(12)

利用拉格朗日乘子法,將模型約束加入到(7)式表示的基追蹤反演的目標函數中,得到模型約束下反演的目標函數為:

(13)

基追蹤極小化問題是一個凸優(yōu)化問題,基追蹤算法從字典中挑選的向量都是線性無關的,可以被理解為一種最佳基算法[21],(7)式是其標準形式。(13)式是稀疏約束與模型約束的聯合約束下的目標函數,通過替換計算可以將其寫作基追蹤的標準形式。令:

(14)

則:

(15)

通過計算我們可以發(fā)現(13)式與(15)式等價,并且(15)式為基追蹤的標準形式,因此我們將(15)式作為最終的模型約束下基追蹤反演的目標函數。

3GPSR方法求解

隨著稀疏表示技術的興起,CHEN等[22]提出基追蹤算法至今,其求解基追蹤問題算法已成為研究的熱點,大量文獻中給出不同的求解算法[19,23-24]。本文選擇GPSR算法對(15)式進行求解。

首先將待反演參數分解為正數序列和負數序列,對負數序列取相反數,則待反演參數m寫作:m=u-v(u≥0,v>0),因此,待反演參數m的L1范數可以寫作(16)式的形式：

(16)

進而將(16)式代入(15)式,則求解(15)式轉為求解目標函數:

(17)

(18)

可以得到新的目標函數:

(19)

最終通過以下步驟進行迭代,可以求得(19)式的解。

1) 初始化,k=0,給定初始值z(0),選擇參數β∈(0,1)與μ∈(0,0.5)。

2) 利用(20)式與(21)式計算α0,為了防止α0值過大,選取合適的αmin與αmax,并且有0<αmin<αmax,令α0=min{αmin,α0,αmax}。

(20)

(21)

(22)

4) 測試是否收斂,如果滿足精度需求,則迭代終止;如果不滿足精度需求,賦值k=k+1,返回步驟2)。

將反演得到的稀疏表示參數進一步利用(4)式轉化為反射系數信息,最終利用(23)式將地層反射信息轉化為地層的阻抗信息。

(23)

4模型測試

4.1楔形模型測試

模型約束基追蹤反演結果呈現塊化,對地層具有較好的分辨能力。首先測試模型約束基追蹤反演對分辨能力的影響,將模型約束下基追蹤反演方法應用于楔形模型。如圖4所示,楔形模型的時間厚度為1～60ms。圖4a為縱波阻抗模型剖面,上層縱波速度為2.2km/s,密度為2.36g/cm3;中間楔形體部分縱波速度為2.4km/s,密度為2.4g/cm3;楔形體下伏地層的縱波速度、密度值與楔形體上覆地層相同。反射系數是由反射界面兩側的阻抗差異引起的,圖4b為利用圖4a中相關參數計算得到的楔形體界面反射系數。利用35Hz雷克子波正演得到的無噪聲合成記錄如圖4c所示。無噪聲情況下反演的反射系數與阻抗分別如圖4d與圖4e所示,可以看出,反演得到的反射系數稀疏性好、連續(xù)性強,并且縱波阻抗呈塊化顯示,能夠較好地識別楔形體的頂、底界面。圖4中紅色橢圓內,小于調諧厚度時反演得到的反射系數與阻抗值與真實模型存在細微的差別,但不影響我們對地層界面的識別。進一步驗證該方法對數據中噪聲的敏感性,在合成記錄中加入服從高斯分布的隨機噪聲(圖4f),信噪比為1。在含噪聲情況下反演得到的反射系數和縱波阻抗分別如圖4g和圖4h所示。從圖4g和圖4h 可以看出,盡管在薄層處存在一定誤差,但反演結果與真實模型仍具有較高的吻合度,反射系數所對應的界面連續(xù)性較好,縱波阻抗對地層的反射界面仍有很好的識別效果,無論是反射系數還是縱波阻抗所識別的反射界面與真實模型都相一致。因此,在原有的基追蹤反演方法基礎上加入模型約束不會改變反演結果對地層界面的識別能力,且能較好地改善反演結果的連續(xù)性。

圖4　楔形模型疊后合成記錄與反演結果a 楔形模型阻抗; b 真實反射系數; c 合成記錄(無噪聲); d 無噪聲時的反射系數反演結果; e 無噪聲時的阻抗反演結果; f 合成記錄(信噪比為1); g 信噪比為1時反射系數反演結果; h 信噪比為1時阻抗反演結果

4.2二維模型資料測試

為了進一步驗證本方法對實際資料的有效性和實用性,對更接近于實際情況的二維地震資料進行了相應的算法測試。在此,選用某二維地震記錄作為測試數據,其疊加剖面如圖5a所示,圖中黑色的豎線為井所在位置。圖5b為反演中所應用的平滑低頻阻抗模型。圖5c與圖5d分別為反演得到的反射系數與縱波阻抗剖面。從圖5中可以看出,反射系數具有較強的稀疏性,縱波阻抗的結果縱向上表現為塊化,對地層分界面識別能力較強。

模型約束基追蹤反演方法保持了反射系數反演的稀疏性,反演地層縱波阻抗以塊化顯示,地層界面識別明顯,并且由于平滑模型的加入增強了反演的穩(wěn)定性,改善了基追蹤反演結果的橫向連續(xù)性。如圖5中箭頭處所示,圖5a中地震同相軸有較為明顯的間斷現象,圖5c和圖5d中反射系數和縱波阻抗反演結果均較為連續(xù)。

圖6為由圖5中提取的井旁道與常規(guī)反演結果的對比,從左至右依次為原始地震記錄井旁道、井旁道反演結果合成記錄、測井曲線計算的縱波阻抗曲線、模型約束下基追蹤反演的縱波阻抗以及常規(guī)反演方法得到的縱波阻抗。

圖5　Hampson-Russell軟件中數據反演測試a 地震記錄; b 阻抗模型; c 反演的反射系數; d 反演的阻抗

從圖6可以看出,利用模型約束基追蹤反演結果的合成記錄與實際地震記錄有較高的相似系數,驗證了該方法反演結果的正確性。圖6中,0.324～0.336s兩條紅線之間為一氣層,比較模型約束基追蹤反演方法與傳統(tǒng)方法的反演結果可以看出,傳統(tǒng)方法的反演結果是平滑的,對儲層的邊界刻畫不夠清晰,而本文方法的反演結果中,氣層邊界清晰,有很強的儲層邊界識別性能(圖中綠色箭頭處)。

圖6　井旁道反演結果對比

5實際資料應用分析

將本文提出的模型約束基追蹤反演方法應用到某實際工區(qū)的疊后地震資料中,疊后地震記錄如圖7a所示。圖7b為反演過程中所使用的縱波阻抗平滑模型。圖7c和圖7d分別為采用本文提出的模型約束基追蹤反演方法反演得到的反射系數與縱波阻抗。

從圖7可以看出,反射系數稀疏,縱波阻抗呈現塊化,對地層界面識別顯著。圖7a中黑色橢圓內為氣藏發(fā)育區(qū),但同相軸出現了較為明顯的間斷現象,經過模型約束基追蹤反演方法反演得到的反射系數與縱波阻抗在橢圓框內地層連續(xù),驗證模型約束基追蹤反演增強了反演的穩(wěn)定性、改善了反演結果的連續(xù)性。

圖7　實際工區(qū)地震記錄與反演結果a 地震記錄; b 縱波阻抗平滑模型; c 反射系數反演結果; d 縱波阻抗反演結果

6結論

本文在Zhang等提出的基追蹤反演的基礎上加入模型約束,得到模型約束條件下的基追蹤反演目標函數,使得反演結果中的低頻信息更加符合工區(qū)內的實際地質背景信息,從而使得反演結果更加合理,增強了反演的穩(wěn)定性,在一定程度上改善了反演結果的橫向連續(xù)性。

由于反演目標函數中存在L1范數稀疏約束,利用GPSR基追蹤反演算法能夠在反演結果合理的條件下,確保反射系數反演結果的稀疏性。在反射系數稀疏的條件下,阻抗反演結果表現為塊化,能夠有效地消除子波旁瓣的影響,對地層邊界的識別更有優(yōu)勢。

參考文獻

[1]ROBINSON E.Predictive decomposition of seismic traces[J].Geophysics,1957,22(4):767-778

[2]TAYLOR H,BANKS S,MC COY J.Deconvolution with theL1norm[J].Geophysics,1979,44(1):39-52

[3]DEBREYE H W J,RIEL P V.Lp-norm deconvolution[J].Geophysical Prospecting,1990,38(4):381-403

[4]GUITTON A,SYMES W W.Robust inversion of seismic data using the Huber norm[J].Geophysics,2003,68(4):1310-1319

[5]張世鑫,印興耀,張繁昌.基于三變量柯西分布先驗約束的疊前三參數反演方法[J].石油地球物理勘探,2011,46(5):737-743

ZHANG S X,YIN Y X,ZHANG F C.Prestack three term inversion method based on Trivariate Cauchy distribution prior constraint [J].Oil Geophysical Propecting,2011,46(5):737-743

[6]雷文文,印興耀,張繁昌.三元柯西約束疊前地震橫波速度反演[C]∥中國地球物理學會.中國地球物理(2011).合肥:中國科學技術大學出版社，2011：649

LEI W W,YIN X Y,ZHANG F C.Seismic prestack shear wave velocity inversion via a Trivariate Cauchy probability distribution[C]∥Chinese Geophysical Society.The Chinese geophysical (2011).Hefei:University of Science and Technology of China Press,2011:649

[7]宗兆云,印興耀,張繁昌.基于改進柯西約束和點約束的稀疏脈沖彈性阻抗反演[C]∥中國地球物理學會.中國地球物理(2009).合肥:中國科學技術大學出版社,2009:204

ZONG Z Y,YIN X Y,ZHANG F C.Sparse spike elastic inversion based on new Cauchy and points constraint[C]∥Chinese Geophysical Society.The Chinese geophysical(2009).Hefei:University of Science and Technology of China Press,2009:204

[8]楊培杰,印興耀.非線性二次規(guī)劃貝葉斯疊前反演[J].地球物理學報,2008,51(6):1876-1882

YANG P J,YIN X Y.Non-linear quadratic programming Bayesian prestack inversion[J].Chinese Journal of Geophysics,2008,51(6):1876-1882

[9]DUIJNDAM A J W.Detailed Bayesian inversion of seismic data[D].The Netherlands:Technische Universiteit Delft,1987

[10]TARANTOLA A.Inverse problem theory and methods for model parameter estimation[M].Philadelphia:Society for Industrial and Applied Mathematics,2005:1-37

[11]THEUNE U,JENS?S I ?,EIDSVIK J.Analysis of prior models for a blocky inversion of seismic AVA data[J].Geophysics,2010,75(3):C25-C35

[12]ZHANG R,SEN M K,SRINIVASAN S.A prestack basis pursuit seismic inversion[J].Geophysics,2013,78(1):R1-R11

[13]ZHANG R,CASTAGNA J.Seismic sparse-layer reflectivity inversion using basis pursuit decomposition[J].Geophysics,2011,76(6):R147-R158

[14]ZHANG R,SEN M K,SRINIVASAN S.Multi-trace basis pursuit inversion with spatial regularization[J].Journal of Geophysics and Engineering,2013,10(3):35012

[15]YIN X Y,YANG P J,ZHANG G Z.A novel prestack AVO inversion and its application[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008,2041-2045

[16]ZONG Z Y,YIN X Y,WU G C.AVO inversion and poroelasticity with P- and S-wave moduli[J].Geophysics,2012,77(6):N17-N24

[17]ZONG Z,YIN X,WU G.Direct inversion for a fluid factor and its application in heterogeneous reservoirs[J].Geophysical Prospecting,2013,61(5):998-1005

[18]ZONG Z,YIN X,WU G.Elastic impedance variation with angle inversion for elastic parameters[J].Journal of Geophysics and Engineering,2012,9(3):247-260

[19]FIGUEIREDO M A T,NOWAK R D,WRIGHT S J.Gradient projection for sparse reconstruction:application to compressed sensing and other inverse problems[J].IEEE Journal of Selected Topics in Signal Processing,2007,1(4):586-597

[20]PéREZ D,VELIS D,SACCHI M.High-resolution prestack seismic inversion using a hybrid FISTA least-squares strategy[J].Geophysics,2013,78(5):R185-R195

[21]MALLAT S.A wavelet tour of signal processing:the sparse way[M].3rded.Salt Lake City:Academid Press,2008:611-696

[22]CHEN S S,DONOHO D L,SAUNDERS M A.Atomic decomposition by basis pursuit[J].Siam Review,2001,43(1):129-159

[23]BERG E V D,FRIEDLANDER M P.Probing the Pareto frontier for basis pursuit solutions[J].Siam Journal on Scientific Computing,2008,31(2):890-912

[24]YANG J,ZHANG Y.Alternating direction algorithms forL1-problems in compressive sensing[J].Eprint Arxiv,2009,33(1):250-278

(編輯：顧石慶)

Basis pursuit inversion method under model constraint

YIN Xingyao,LIU Xiaojing,WU Guochen,ZONG Zhaoyun

(SchoolofGeosciences,ChinaUniversityofPetroleum,Qingdao266580,China)

Abstract:The basis pursuit inversion (BPI) method based on the dipole decomposition of reflection coefficients compensates the low-frequency and high-frequency information lacked in the seismic data,which can well improve the identification capability of inversion results on formations.However,it might not be reasonable to add the low frequency information just by the sparse constraint,and the information might not coincide with the true geological background.Therefore,we derived the model constrained BPI objective function by putting the model constraint to the origin objective function of BPI,and utilized the gradient projection for sparse reconstruction (GPSR) algorithm to solve the problem.It would not only enhance the inversion stability and make the low frequency information contained in inversion results more consistent with the actual work area,but also improve the lateral continuity of the inversion results.The results of the wedge model and real field data testing demonstrate that the model constrained BPI method not only keep the sparsity of the inversion results and shows the impedance in blocky,which characterizes the interface more clearly and improves the reliability of the inversion and the lateral continuity of the inversion results.

Keywords:model constraint,basis pursuit,seismic inversion,sparsity,lateral continuity

文章編號：1000-1441(2016)01-0115-08

DOI:10.3969/j.issn.1000-1441.2016.01.015

中圖分類號：P631

文獻標識碼：A

基金項目：國家自然科學基金石油化工聯合基金重點項目(U1562215)、國家重點基礎研究發(fā)展計劃(973計劃)項目(2013CB228604)與國家科技重大專項(2011ZX05030-004-002,2011ZX05006-002,2011ZX05009-003)聯合資助。

作者簡介：印興耀(1962—),男,教授,博士生導師,從事勘探地球物理理論與方法的教學與科研工作。

收稿日期：2014-02-25;改回日期：2015-06-28。

印興耀,劉曉晶,吳國忱,等.模型約束基追蹤反演方法[J].石油物探,2016,55(1):-122

YIN Xingyao,LIU Xiaojing,WU Guochen,et al.Basis pursuit inversion method under model constraint[J].Geophysical Prospecting for Petroleum,2016,55(1):-122

This research is financially supported by the Petrochemical Joint Foundation Major Project of National Natural Science Foundation of China (Grant No.U1562215),the National Key Basic Research Program of China (973 Program) (Grant No.2013CB228604) and the National Science and Technology Major Project of China (Grant Nos.2011ZX05030-004-002,2011ZX05006-002,2011ZX05009-003).