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

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

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

0 引言

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

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

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

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

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

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

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

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

Y=R*W+N

(1)

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

(2)

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

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

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

(3)

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

(4)

其中

(5)

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

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

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

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

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

(6)

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

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

(7)

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

(8)

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

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

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

(9)

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

(10)

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

(11)

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

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

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

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

(12)

(13)

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

m=0,1,…

(14)

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

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

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

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

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

圖3 小、中、大角度子波

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

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

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

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

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

反演中為了提高穩(wěn)定性，需要引入一些約束條件。通過在改進(jìn)Fatti方程中引入ρ與vP、vS與vP兩個(gè)擬合式以及對(duì)應(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ù)性和縱向分層特點(diǎn)證明vP(圖5a)和ρ(圖5b)總體反演效果均較好，但在井1周圍連續(xù)性很差，這與該處地層發(fā)育特點(diǎn)及地震資料品質(zhì)有關(guān)。

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

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

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

圖7 18°～24°角度道集混合反演結(jié)果與井2(左)、井1(右)的測(cè)井?dāng)?shù)據(jù)對(duì)比

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

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

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

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

此外，通過疊前地震數(shù)據(jù)可直接反演vP/vS，再結(jié)合vP可以獲得vS，然后由Zoeppritz方程計(jì)算反射系數(shù)。vP/vS初始模型(圖4d)由vP初始模型(圖4a)和vS初始模型(圖4c)獲得。圖9為由不同角度疊前地震數(shù)據(jù)直接反演的vP/vS。對(duì)比圖9與圖8可見：前者(圖9)時(shí)窗上部的vP/vS值大于1.8，時(shí)窗下部的vP/vS值為1.7～1.8，都明顯高于中部的值(1.6～1.7)；后者(圖8)時(shí)窗上部的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的測(cè)井?dāng)?shù)據(jù)對(duì)比。由圖可見：vP/vS間接反演結(jié)果與井?dāng)?shù)據(jù)存在誤差(與井2吻合較好，與井1吻合略差)，這主要是因?yàn)閱为?dú)反演vP、vS很難保證vP/vS的值穩(wěn)定；vP/vS直接反演結(jié)果與井?dāng)?shù)據(jù)吻合很好，尤其是18°～24°角度道集的結(jié)果(圖10a)。

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

3 結(jié)論

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

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