基于Zoeppritz方程的疊前和疊后混合多參數(shù)非線性地震反演

張凌遠 張宏兵* 尚作萍 嚴立志 任 權(quán)

(①河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇南京 211100；②河海大學(xué)力學(xué)與材料學(xué)院，江蘇南京 211100)

0 引言

最近十幾年，基于Zoeppritz方程的疊前反演可以同時獲得縱、橫波阻抗及速度、縱橫波速度比、密度、泊松阻抗、泊松比等儲層參數(shù)，已經(jīng)成為必不可少的特殊處理技術(shù)。盡管如此，疊后地震反演也不可替代，尤其是在一些早期勘探或疊前地震資料品質(zhì)不高的探區(qū)，疊后地震反演可以作為一個必要的補充，需要結(jié)合疊前和疊后地震數(shù)據(jù)開展混合反演。王秀娟等[1]應(yīng)用混合反演技術(shù)分析了中國南海北部水合物賦存情況；李磊等[2]分析了影響混合反演效果的層速度、角度道集等數(shù)據(jù)，提高了疊前反演的實用性，降低了疊后反演的多解性。

地震反演是一項復(fù)雜的系統(tǒng)工程[3]，涉及正演物理數(shù)學(xué)模型、反演目標函數(shù)、適定性問題、搜索算法等。在正演物理數(shù)學(xué)模型方面，全波形反演依據(jù)的波動方程無疑具有巨大的理論優(yōu)勢[4-7]，但是由于計算效率問題其商業(yè)化應(yīng)用效果不高。因此，實際生產(chǎn)中使用的還是基于褶積模型的疊前反演和疊后反演。由于精確Zoeppritz方程的表達式較復(fù)雜，現(xiàn)有商業(yè)化軟件多使用其近似式，如Aki-Richard近似式[8]、Shuey近似式[9]和Fatti近似式[10]等。但是這些近似式是建立在一定的假設(shè)條件之下，其應(yīng)用受到不同程度的限制，因此人們嘗試使用精確Zoeppritz方程開展疊前反演?？v橫波速度比或泊松比對建立氣水、油水識別模型至關(guān)重要，目前通過Shuey近似式或Gray近似式及基追蹤理論反演泊松比[11]，也可以依據(jù)疊前反演獲得的縱、橫波速度計算縱橫波速度比，但不可避免地降低了參數(shù)反演精度。為此，需要嘗試通過精確的Zoeppritz方程直接反演縱橫波速度比或泊松比。

關(guān)于反演目標函數(shù)及適定性問題，疊后反演或疊前反演均使用褶積模型誤差項的極小化或后驗概率的極大化構(gòu)建地震反演目標函數(shù)[12]。地震反演受不適定性所困擾，因此反演是一個病態(tài)問題，體現(xiàn)在正演模型的非線性、反演結(jié)果的非唯一性和數(shù)據(jù)誤差引起的不穩(wěn)定性等方面。目前，解決地震反演不適定性的主要方法是正則化方法，傳統(tǒng)的正則化方法基于Tikhonov正則化思想[13-14]，通過在目標函數(shù)中添加模型參數(shù)的零階或高階先驗項穩(wěn)定反演結(jié)果，但是反演結(jié)果可能過于平滑[15]。為此，人們提出了邊界保護正則化的思想用于地震反演[16-18]。此外，在疊后或疊前地震多參數(shù)同步反演中，各參數(shù)之間數(shù)量級的巨大差異也會引起多參數(shù)同步反演結(jié)果的不穩(wěn)定。為了控制不同數(shù)量級參數(shù)同步反演的穩(wěn)定性，業(yè)界開展了一些有意義的研究，如在改進Fatti方程中引入密度與縱波速度、橫波速度與縱波速度兩個擬合式以及對應(yīng)的密度偏差量、橫波速度偏差量提高反演精度[19-20]。

地震反演算法的研究成果眾多，如寬帶約束反演、廣義逆等線性反演算法[21-22]、最速下降法、共軛梯度法、擬牛頓法等局部非線性優(yōu)化算法[23-24]以及粒子群算法、遺傳算法、模擬退火算法等全局非線性優(yōu)化算法[25-27]。現(xiàn)有的研究成果表明，全局非線性優(yōu)化算法具有不依賴初始模擬、可以獲得全局最優(yōu)估計解等優(yōu)點，但是其計算量遠大于線性優(yōu)化算法和局部非線性優(yōu)化算法。目前，商業(yè)應(yīng)用以線性反演算法和局部非線性優(yōu)化算法為主。

本文重點關(guān)注疊前和疊后地震混合反演以及基于精確Zoeppritz方程直接反演縱橫波速度比。首先，建立基于邊界保護正則化的疊前或疊后地震反演目標函數(shù)，結(jié)合快速模擬退火算法進行非線性反演。利用疊后地震反演同步獲得縱波速度和密度，基于精確Zoeppritz方程的疊前地震反演獲得橫波速度或縱橫波速度比。最后，通過實際資料測試、驗證方法的有效性。

1 疊前和疊后地震混合多參數(shù)反演

1.1 疊前和疊后地震響應(yīng)數(shù)學(xué)模型

通常使用褶積模型表征疊后地震響應(yīng)的數(shù)學(xué)模型，即

Y=R*W+N

(1)

式中：Y為疊后地震記錄；R為地震反射系數(shù)序列；W為疊后子波；N為隨機噪聲。一般情況下，假設(shè)N服從高斯分布，其數(shù)學(xué)期望為零，協(xié)方差為σ。垂直入射時R的離散形式為

(2)

式中：下標i為介質(zhì)層號；Z為波阻抗；vP為縱波速度；ρ為密度。疊后地震反演通常反演波阻抗，如果同步反演縱波速度和密度，需要考慮反演的穩(wěn)定性問題。

與疊后地震響應(yīng)類似，疊前地震角度道集的數(shù)學(xué)模型可表示為

Y(θ)=R(θ)*W(θ)+N

(3)

式中：Y(θ)為疊前地震記錄，θ為入射角；R(θ)為地震反射系數(shù)序列；W(θ)為疊前角度子波。R(θ)可以由Zoeppritz方程或其近似式模擬，理論基礎(chǔ)是AVO理論，即描述平面波反射和透射系數(shù)相對入射角變化的Zoeppritz方程

(4)

其中

(5)

式中：T為透射系數(shù)；R為反射系數(shù)；v為速度；下標P和S分別對應(yīng)縱波和轉(zhuǎn)換橫波(為方便表達和閱讀，下文均用“橫波”代替“轉(zhuǎn)換橫波”)，下標兩個字母中第一個字母表示入射波型，第二個字母表示反射或透射波型；θ為P波反射角或透射角，下標1和2分別指示反射界面上覆介質(zhì)和下伏介質(zhì)；φ為S波反射角或透射角。由式(5)可見：RPP不僅與介質(zhì)的彈性參數(shù)有關(guān)，還與入射角θ1有關(guān)，即式(3)中的R(θ)=RPP；當垂直入射(θ1=0)時，RPP與式(2)的Ri相同。

考慮到式(5)較復(fù)雜，物理意義不明確，為此人們提出了一些簡化的近似公式，如Aki近似式[8]、Shuey近似式[9]、Fatti近似式[10]等。這些近似公式物理意義明確，但是受假設(shè)條件(如相鄰地層介質(zhì)彈性參數(shù)變化較小、縱橫波速度比為2.0以及小角度入射等)限制。為此，這里直接使用精確Zoeppritz方程求取反射系數(shù)，以便開展高精度非線性疊前反演。需要說明的是，針對不同的角道集，需要使用不同的角度子波。

1.2 基于邊界保護正則化的反演目標函數(shù)

針對式(1)，將先驗信息引入最小二乘問題，則疊后或疊前地震反演的目標函數(shù)為

J(X)=J1(X)+λ·J2(X)

(6)

式中：J1(X)為數(shù)據(jù)項，X為待反演的參數(shù)；J2(X)為先驗項；D(X)為X的梯度值；Ck為Markov隨機域(MRF)鄰域系統(tǒng)的數(shù)據(jù)點集，k=1, 2 ,3為平滑階數(shù)；Φ為勢函數(shù)；λ為權(quán)系數(shù)；δ為刻度參數(shù)，在被檢測的不連續(xù)處調(diào)節(jié)梯度值。值得注意的是，多參數(shù)(縱、橫波速度和密度)同步反演的J2(X)由三部分組成，即

J2=J2P(vP)+J2S(vS)+J2D(ρ)

(7)

式中J2P(vP)、J2S(vS)和J2D(ρ)分別為vP、vS和ρ的先驗項。于是J2(X)可以重寫為

(8)

式中刻度參數(shù)δP、δS和δρ的取值不同。

需要說明的是，由于本文采用疊后與疊前混合反演，因此在疊后和疊前反演過程中J2(X)是不同的。在疊后反演中

∑CkΦ[DCk(ρ)/δρ]

(9)

在疊前單參數(shù)(如橫波速度)反演中

(10)

如果反演縱橫波速度比，則

(11)

式中：γ=vP/vS為縱橫波速度比；δγ是關(guān)于γ反演的刻度參數(shù)。

有關(guān)勢函數(shù)，主要考慮其邊界保護特性，詳細內(nèi)容見文獻[28]。在實際資料反演中，僅使用具有邊界保護特性的勢函數(shù)ΦGM，勢函數(shù)的具體形式見文獻[14]。此外，為了加快反演收斂速度和改善反演結(jié)果，在反演迭代過程中不斷調(diào)節(jié)正則參數(shù)，即要求λ逐漸減小，δP、δS、δρ和δγ逐漸增大。

1.3 基于模擬退火的反演搜索方法

疊前或疊后地震多參數(shù)同步反演可以看作是使目標函數(shù)(式(6))極小化的優(yōu)化問題。由于式(6)中模型參數(shù)與測量數(shù)據(jù)之間的非線性關(guān)系非常復(fù)雜，為了獲得一個全局最優(yōu)解，使用快速模擬退火算法(FSA)求解目標函數(shù)極小值。FSA需要解決3個問題：①下一個候選點的擾動值；②接受概率；③退火過程中冷卻進度表，包括初始溫度、終止溫度和溫度衰減系數(shù)。本文接受概率使用廣義Gibbs分布函數(shù)，在數(shù)據(jù)測試中確定冷卻進度表，尤其是溫度衰減系數(shù)需要折中反演精度和計算效率。vP、ρ(疊后反演)、vS、γ(疊前單參數(shù)反演)的擾動值分別為

(12)

(13)

γ(m+1)=γ(m)+T(m)sign(ζ4-0.5)×

m=0,1,…

(14)

式中：vP(m)、vS(m)、ρ(m)和γ(m)表示當前值；vP(m+1)、vS(m+1)、ρ(m+1)和γ(m+1)為擾動后值；T(m)為當前溫度值；[vPmin，vPmax]、[vSmin，vSmax]、[ρmin，ρmax]和[γmin，γmax]分別為4個反演參數(shù)的變化范圍；ξ和ζ為0到1的隨機數(shù)(下標1、2、3、4表示取不同的隨機數(shù))；sign(·)為符號函數(shù)。

2 實測疊前和疊后地震數(shù)據(jù)混合反演

圖1為研究區(qū)過井二維地震剖面A。由圖可見，目的層橫向變化很大，地震同相軸連續(xù)性較差。提取了疊前地震數(shù)據(jù)角度道集(圖2)的小、中、大角度子波(圖3)。使用3°～9°、18°～24°和33°～39°等3組角度道集數(shù)據(jù)進行疊前地震反演。圖4為vP、ρ、vS和vP/vS初始模型。本文的疊前和疊后混合反演步驟如下。

圖1 研究區(qū)過井二維地震剖面A井1和井2分別位于CDP 2226和CDP 918。疊后地震數(shù)據(jù)包含2542個CDP，時間窗口為2200～2600ms，采樣率為2ms

圖2 疊前地震數(shù)據(jù)角度道集每個CDP有15個角度道集，角度范圍為3°～45°，角度間隔為3°

圖3 小、中、大角度子波

(1)疊后反演vP和ρ。使用vP(圖4a)、ρ(圖4b)初始模型反演vP和ρ，在搜索時利用式(12)進行擾動。

(2)利用疊后反演獲得的vP和ρ開展疊前反演。將vP和ρ作為已知值代入Zoeppritz方程，分別使用vS初始模型(圖4c)、vP/vS初始模型(圖4d)反演vS、vP/vS，在搜索時分別使用式(13)、式(14)進行擾動。

圖4 vP(a)、ρ(b)、vS(c)和vP/vS(d)初始模型

2.1 vP、vS和ρ的混合反演結(jié)果

采用疊后和疊前多參數(shù)混合反演的理由為：①疊后地震資料的品質(zhì)好于疊前地震資料；②由于各參數(shù)之間的數(shù)量級差異，疊前反演同步獲得3參數(shù)的精度很難保證。為此，這里嘗試使用疊后地震數(shù)據(jù)同時反演vP和ρ，在此基礎(chǔ)上使用疊前角度道集地震數(shù)據(jù)進行單參數(shù)(vS、vP/vS、泊松比等)反演。

反演中為了提高穩(wěn)定性，需要引入一些約束條件。通過在改進Fatti方程中引入ρ與vP、vS與vP兩個擬合式以及對應(yīng)的ρ偏差量、vS偏差量提高反演精度。疊后反演使用ρ與vP的擬合式

ρ=1.036+0.0003825vP

(14)

式中ρ偏離擬合線的范圍為(-0.125g/cm3，0.125g/cm3)。

疊前反演使用vS與vP的擬合式

vS=1775.2+0.08818vP

(15)

式中vS偏離擬合線的范圍為(-200m/s，200 m/s)。

圖5為疊后地震反演的vP和ρ，其中疊后反演的反射系數(shù)為垂直入射反射系數(shù)。由圖可見，橫向連續(xù)性和縱向分層特點證明vP(圖5a)和ρ(圖5b)總體反演效果均較好，但在井1周圍連續(xù)性很差，這與該處地層發(fā)育特點及地震資料品質(zhì)有關(guān)。

圖5 疊后地震反演的vP(a)和ρ(b)

vP的值域為(3200m/s，4100m/s)，ρ的值域為(2.25g/cm3，2.60g/cm3)。模擬退火算法的初始溫度為0.5，終止溫度為0.00001，溫度衰減系數(shù)為0.9。正則化勢函數(shù)為ΦGM，正則參數(shù)λ=0.3，δ=(δPδSδρδγ)=(150.0,110.0,0.15,0.11)。利用兩口井的測井參數(shù)約束反演在疊后地震反演獲得的vP和ρ基礎(chǔ)上，進行疊前地震單參數(shù)(vS)反演(圖6)，其中反射系數(shù)直接由Zoeppritz方程獲得?？梢?，中(圖6b)、大(圖6c)角度道集的反演效果好于小角度道集(圖6a)，這是由于小角度道集地震資料品質(zhì)較差所致。圖7為18°～24°角度道集混合反演結(jié)果與井2、井1的測井數(shù)據(jù)對比。由圖可見，疊后反演的ρ與井數(shù)據(jù)吻合很好，vP與井數(shù)據(jù)基本吻合。郭強等[29]認為，中角度道集疊前地震反演的vP、vS的參數(shù)敏感性總體優(yōu)于大角度道集。因此后文重點分析中角度道集反演結(jié)果。

圖6 不同角度疊前地震反演vS(a)3°～9°；(b)18°～24°；(c)33°～39°

圖7 18°～24°角度道集混合反演結(jié)果與井2(左)、井1(右)的測井數(shù)據(jù)對比

vS的值域為(1600m/s，2300m/s)，模擬退火算法的初始溫度為0.5，終止溫度為0.00001，溫度衰減系數(shù)為0.9。正則化勢函數(shù)為ΦGM，正則參數(shù)λ=0.3，δ=(δPδSδρδγ)=(150.0,110.0,0.15,0.11)。利用兩口井的測井參數(shù)約束反演。

需要特別說明的是，由于研究區(qū)小角度道集地震數(shù)據(jù)品質(zhì)不好，而中角度道集反演結(jié)果又略優(yōu)于大角度道集。因此“中角度道集反演效果最好”的認識不具普適性。

2.2 vP/vS直接反演結(jié)果

利用vP/vS能夠有效識別氣層(或氣水同層)與水層。取疊后反演vP(圖5a)與疊前反演vS(圖6)的比值，得到不同角度vP/vS(圖8)——間接反演結(jié)果?？梢?，在較深時間段vP/vS值明顯降低，可能與含氣有關(guān)。

此外，通過疊前地震數(shù)據(jù)可直接反演vP/vS，再結(jié)合vP可以獲得vS，然后由Zoeppritz方程計算反射系數(shù)。vP/vS初始模型(圖4d)由vP初始模型(圖4a)和vS初始模型(圖4c)獲得。圖9為由不同角度疊前地震數(shù)據(jù)直接反演的vP/vS。對比圖9與圖8可見：前者(圖9)時窗上部的vP/vS值大于1.8，時窗下部的vP/vS值為1.7～1.8，都明顯高于中部的值(1.6～1.7)；后者(圖8)時窗上部的vP/vS值偏低，分布于1.65～1.8。

圖8 不同角度vP/vS間接反演結(jié)果(a)18°～24°； (b)33°～39°

圖9 由不同角度疊前地震數(shù)據(jù)直接反演的vP/vS(a)18°～24°； (b)33°～39°

圖10為不同角度反演獲得的vP/vS與井2、井1的測井數(shù)據(jù)對比。由圖可見：vP/vS間接反演結(jié)果與井數(shù)據(jù)存在誤差(與井2吻合較好，與井1吻合略差)，這主要是因為單獨反演vP、vS很難保證vP/vS的值穩(wěn)定；vP/vS直接反演結(jié)果與井數(shù)據(jù)吻合很好，尤其是18°～24°角度道集的結(jié)果(圖10a)。

圖10 不同角度反演獲得的vP/vS與井2(上)、井1(下)的測井數(shù)據(jù)對比(a)18°～24°； (b)33°～39°

3 結(jié)論

(1)直接使用Zoeppritz方程進行反演，可以避免Zoeppritz方程近似式反演引起的誤差，對于大角度道集尤其重要。使用Zoeppritz方程直接反演幾乎并不影響多參數(shù)反演的速度和效果。此外，由于多參數(shù)之間的量級差異以及多參數(shù)同步隨機搜索增大了不確定性，在反演中引入兩個附加約束條件，提高了多參數(shù)反演結(jié)果穩(wěn)定性。

(2)疊后和疊前混合多參數(shù)反演可以克服疊前地震多參數(shù)同步反演中不同參數(shù)數(shù)量級、參數(shù)敏感性、疊前地震資料品質(zhì)差異引起的問題。由于研究區(qū)的實際疊前地震道集品質(zhì)較差，而疊后地震數(shù)據(jù)品質(zhì)較好，故可由疊后反演獲得縱波速度和密度，由疊前反演獲得橫波速度、縱橫波速度比、泊松比等。結(jié)果顯示，疊后反演獲得的縱波速度和密度精度較高，疊前反演獲得的縱橫波速度比的精度高于疊后反演，18°～24°角度道集的反演效果好于3°～9°和33°～39°角度道集。