非均勻介質(zhì)地震AVO響應(yīng)模擬及分析

劉 財(cái)，裴思嘉，郭智奇，符 偉，張宇生，劉喜武

1.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026 2.中石油東方地球物理公司，河北 涿州 072751 3.中國石化頁巖油氣勘探開發(fā)重點(diǎn)實(shí)驗(yàn)室，北京 100083 4.中國石化石油勘探開發(fā)研究院，北京 100083

0 引言

地震反射系數(shù)反映地下界面物性特征的變化。對(duì)于各向同性介質(zhì)分界面，平面波反射、透射系數(shù)的精確解可由Zoeppritz[1]方程計(jì)算。許多學(xué)者對(duì)Zoeppritz方程進(jìn)行了不同形式的簡(jiǎn)化，以此開展地震AVO(amplitude versus offset)技術(shù)研究，通過反射波振幅隨偏移距的變化規(guī)律研究反射界面兩側(cè)介質(zhì)的物性特征。傳統(tǒng)的AVO分析都是基于Zoeppritz方程的各自線性近似式提出的[2-4]。1955年，Koefoed[5]不僅簡(jiǎn)化了Zoeppritz方程，而且在實(shí)驗(yàn)中假設(shè)上下界面的泊松比不同，所得的縱波反射系數(shù)隨入射角的改變而變化范圍大大增大，這對(duì)之后的AVO技術(shù)發(fā)展具有重大意義。1980年，Aki等[6]將Zoeppritz方程簡(jiǎn)化為入射角與縱波振幅相關(guān)的近似公式。1985年，Shuey[7]根據(jù)該簡(jiǎn)化公式做了進(jìn)一步研究，將其改寫為入射角小角度項(xiàng)、中角度項(xiàng)和廣角度項(xiàng)三部分之和的近似公式。

然而，層狀介質(zhì)的AVO特征與單一界面是不同的，因?yàn)榉瓷湎禂?shù)與入射波的頻率相關(guān)。實(shí)際地層為非均勻介質(zhì)，地震波在地下傳播時(shí)具有傳播效應(yīng)和層間多次波的調(diào)諧干涉效應(yīng)。對(duì)此，早在1971年，F(xiàn)uchs等[8]就應(yīng)用反射率方法計(jì)算了多層介質(zhì)對(duì)點(diǎn)震源的響應(yīng)。1980年，Aki等[6]提出傳播矩陣算法，在理論上可以更好地描述這種層狀介質(zhì)模型。1983年，Kennett[9]給出了地震波在層狀各項(xiàng)異性介質(zhì)中傳播的描述。1984年，F(xiàn)ryer等[10]又將傳播矩陣用于各向異性介質(zhì)，制作了合成地震記錄。在此之后，有許多學(xué)者對(duì)層狀介質(zhì)進(jìn)行了AVO正演響應(yīng)研究。2001年，Carcione[11]計(jì)算了富有機(jī)質(zhì)頁巖層的反射系數(shù)。2003年，Liu等[12]研究了泊松比和層厚度對(duì)薄層AVO響應(yīng)的影響。2005—2016年，郭智奇等[13-15]基于廣義傳播矩陣?yán)碚撚?jì)算地震波頻變反射系數(shù)，完善了頻變AVO算法，為含油氣儲(chǔ)層頻變AVO響應(yīng)的模擬和分析提供了方法。2014年，蘭慧田[16]應(yīng)用廣義傳播矩陣?yán)碚撚?jì)算了上覆介質(zhì)與裂縫性孔隙介質(zhì)分界面的反射系數(shù)，為裂縫檢測(cè)和裂縫參數(shù)反演提供地震波振幅信息。

本文針對(duì)非均勻地下介質(zhì)模型，選取薄互層結(jié)構(gòu)和頁巖氣儲(chǔ)層兩組實(shí)際測(cè)井?dāng)?shù)據(jù)，分別比較精確的Zoeppritz方程、Shuey二項(xiàng)近似方程和Shuey三項(xiàng)近似方程以及傳播矩陣算法得到的正演模擬響應(yīng)計(jì)算效果，并與實(shí)際的地震數(shù)據(jù)進(jìn)行對(duì)比分析，討論存在調(diào)諧干涉和傳播效應(yīng)的情況下不同方法的精確性和適用性。

1 AVO正演模擬基本理論

1.1 Zoeppritz方程

Zoeppritz方程[1]是AVO分析的基礎(chǔ)，該方程通過平面波的入射角和反射界面上下介質(zhì)的6個(gè)彈性參數(shù)求取界面的反射系數(shù)和透射系數(shù)，利用界面兩側(cè)應(yīng)力和位移守恒的原則推導(dǎo)反射波和透射波振幅。反射、透射模式如圖1所示。當(dāng)平面波入射到界面時(shí)，波前面與y軸是平行的，所以與方程、y軸無關(guān)。當(dāng)P波以一定入射角入射到界面時(shí)，將會(huì)產(chǎn)生反射P波、反射S波、透射P波和透射S波，其Zoeppritz方程表達(dá)式為

(1)

式中：RPP、RPS分別為縱波和橫波的反射系數(shù)；TPP、TPS分別為縱波和橫波的透射系數(shù)。

vP1、vP2分別為反射界面上下層介質(zhì)的縱波速度；vS1、vS2分別為反射界面上下層介質(zhì)的橫波速度；ρ1、ρ2分別為反射界面上下層介質(zhì)的密度。圖1 入射波、反射波和透射波的關(guān)系Fig.1 Relationship amongincident wave, reflection wave and transmission wave

1.2 Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程

Shuey[7]將Zoeppritz方程簡(jiǎn)化為關(guān)于入射角θ的函數(shù)，并分為小角度項(xiàng)、中角度項(xiàng)和大角度項(xiàng)三部分，利用泊松比來代替橫波速度，推導(dǎo)出縱波反射系數(shù)(R(θ))公式。

(2)

其中，

(3)

(4)

(5)

(6)

則縱波反射系數(shù)公式為

(7)

式中：ΔvP、Δρ和Δσ分別為反射界面兩側(cè)介質(zhì)縱波速度、密度和泊松比的差值；vP、ρ、σ分別為反射界面兩側(cè)介質(zhì)縱波速度、密度和泊松比的平均值。Shuey三項(xiàng)近似方程(式(7))第一項(xiàng)為垂直入射時(shí)的反射系數(shù)；第二項(xiàng)主要表示中等角度入射時(shí)的反射系數(shù)；第三項(xiàng)為大角度入射時(shí)的反射系數(shù)。當(dāng)入射角較小、并且vP/vS≈2時(shí)，第三項(xiàng)可以忽略(vS為反射界面兩側(cè)介質(zhì)橫波速度的平均值)，此時(shí)式(7)只有前兩項(xiàng)，稱為Shuey二項(xiàng)近似方程。

1.3 傳播矩陣?yán)碚?/h3>

對(duì)于非均勻地下介質(zhì)，地震波的反射特征不僅與入射角度、物性差異有關(guān)，還與入射波頻率、地層厚度、薄互層結(jié)構(gòu)、地層的不均勻性等因素有關(guān)[9,12,17-18]。

由于干涉、調(diào)諧等現(xiàn)象的存在，來自薄互層的反射地震波呈現(xiàn)復(fù)合波型，與基于界面模型的Zeoppritz理論相比，傳播矩陣?yán)碚摳m用于非均勻介質(zhì)情況。根據(jù)傳播矩陣?yán)碚揫19]，對(duì)于P波入射，地層的反射、透射系數(shù)向量r=[RPP,RPS,TPP,TPS]T為

(8)

式中：矩陣A1與A2分別為與上、下層介質(zhì)物性參數(shù)有關(guān)的傳播矩陣；Bα(α=1,…,N)為具有N層結(jié)構(gòu)的中間薄互層的傳播矩陣；iP為P波入射向量，與入射介質(zhì)物性參數(shù)有關(guān)。上述矩陣和向量都是入射波頻率和波慢度的函數(shù)。

傳播矩陣A1和A2分別為：

(9)

(10)

β=

(11)

γ=

(12)

W=p55(γsx+βsz)，

(13)

Z=βp13sx+γp33sz。

(14)

式中：p11、p33、p55和p13為相應(yīng)層位各項(xiàng)異性彈性參數(shù)；p.v.意為取復(fù)數(shù)的主值；對(duì)于γ，符號(hào)“+”對(duì)應(yīng)準(zhǔn)縱波，符號(hào)“-”對(duì)應(yīng)準(zhǔn)橫波；sx為水平波慢度；sz為垂直波慢度。

p33cos2θ±E)-1/2，

(15)

E={[(p33-p55)cos2θ-(p11-p55)·

sin2θ]2+(p13+p55)2sin22θ}1/2；

(16)

(17)

K1=

(18)

(19)

(20)

其中，sz表達(dá)式符號(hào)約定為：(+,-)表示向下傳播準(zhǔn)縱波；(+,+)表示向下傳播準(zhǔn)橫波；(-,-)表示向上傳播準(zhǔn)縱波；(-,+)表示向上傳播準(zhǔn)橫波。

傳播矩陣Bα為

Bα=T(0)T-1(hα)。

(21)

其中，

(22)

P波入射向量iP為

iP=iω[βqP1,γqP1,-ZqP1,-WqP1]T。

(23)

由薄互層反射系數(shù)的傳播矩陣?yán)碚?，可?jì)算各頻率下的反射、透射系數(shù)向量r相應(yīng)反射波的頻變反射系數(shù)Rf。將頻變反射系數(shù)與頻率域的地震子波Wf相乘可得相應(yīng)反射波的振幅譜Uf，即

Uf=Wf×Rf。

(24)

對(duì)Uf做反傅里葉變換則可以獲得時(shí)間域的反射波波形ut：

(25)

式中：f為頻率；t為時(shí)間。

2 地震AVO模擬及分析

2.1 砂泥巖薄互層

圖2a為砂泥巖薄互層結(jié)構(gòu)的實(shí)際數(shù)據(jù)測(cè)井曲線，包括縱、橫波速度和密度?；趥鞑ゾ仃?yán)碚撚?jì)算該非均勻地層的地震響應(yīng)，圖2b為計(jì)算得到的頻變反射系數(shù)。考慮入射波為40 Hz的Ricker子波，進(jìn)一步計(jì)算反射波頻譜如圖2c所示。

對(duì)圖2c所示的反射波頻譜進(jìn)行反傅里葉變換，得到如圖3所示的時(shí)間域合成地震記錄，從中選取兩個(gè)強(qiáng)反射界面T1、T2進(jìn)行井震標(biāo)定。圖4比較了基于傳播矩陣算法、Zoeppritz方程以及Shuey二項(xiàng)近似方程和Shuey三項(xiàng)近似方程計(jì)算得到的合成地震記錄。從圖4中可以看到，4種方法得到的AVO地震響應(yīng)在不同入射角范圍下波形變化也不同，盡管Zoeppritz算法和傳播矩陣方法計(jì)算的AVO波形在小角度入射時(shí)(入射角<20°)比較接近，但是中、大角度(入射角>20°)入射時(shí)存在明顯差別。如圖4中A、B、C三處可以清晰地看到兩種波形在振幅和走向上都有明顯的不同，A處Zoeppritz方程計(jì)算的振幅明顯強(qiáng)于傳播矩陣，B和C處兩種方法得到的合成記錄強(qiáng)反射界面位置和振幅都有差異。分析可知，相比于傳播矩陣算法，Shuey二項(xiàng)近似方程和Shuey三項(xiàng)近似方程在大角度入射時(shí)與Zoeppritz方程更為接近，對(duì)于此數(shù)據(jù)，入射角為0°～50°時(shí)，相比于Shuey二項(xiàng)近似方程，Shuey三項(xiàng)近似方程與Zoeppritz方程吻合度也更高；傳播矩陣算法與其他幾種方法相比，由于考慮了地震反射的傳播效應(yīng)、調(diào)諧效應(yīng)等動(dòng)力學(xué)特征，可以更加清晰地反映薄互層結(jié)構(gòu)的AVO特征變化，由圖4也可看出其垂向分辨率在這幾種正演算法中最高。

圖2 砂泥巖薄互層儲(chǔ)層測(cè)井曲線(a)、頻變反射系數(shù)譜(b)和反射波頻譜(c)Fig.2 Sand-shale thin interbed logging curve(a), reflection coefficients spectrum (b) and amplitude-spectrum of reflection wave (c) in frequency domain of sand-shale thin interbed

圖3 砂泥巖薄互層結(jié)構(gòu)數(shù)據(jù)井震標(biāo)定Fig.3 Well to seismic calibration of sand-shale thin interbed data

a.傳播矩陣算法；b. Zoeppritz方程；c. Shuey二項(xiàng)近似方程；d. Shuey三項(xiàng)近似方程。圖4 砂泥巖薄互層合成地震記錄Fig.4 Synthetic seismogram of sand-shale thin interbed data

2.2 頁巖氣儲(chǔ)層

圖5a所示為某頁巖氣儲(chǔ)層的縱、橫波速度和密度測(cè)井曲線。從圖5a中可以看到，縱、橫波速度曲線的相關(guān)性很強(qiáng)，圖中2 330～2 415 m為頁巖地層，含氣儲(chǔ)層為2 380～2 415 m[20]。基于傳播矩陣算法計(jì)算的頻變反射系數(shù)和反射波頻譜分別如圖5b和圖5c所示。圖6所示為井震標(biāo)定結(jié)果，選取T1和T2兩個(gè)標(biāo)志層進(jìn)行井震標(biāo)定。

圖7為實(shí)際觀測(cè)數(shù)據(jù)和基于傳播矩陣?yán)碚?、Zoeppritz方程、Shuey二項(xiàng)近似方程、Shuey三項(xiàng)近似方程的高精度疊前正演合成記錄。與Zoeppritz方程、Shuey二項(xiàng)近似方程、Shuey三項(xiàng)近似方程相比，傳播矩陣算法在入射角為20°～35°(圖7中T1、T2標(biāo)志層)時(shí)有明顯振幅減小的趨勢(shì)，而其他方法振幅沒有該趨勢(shì)。Zoeppritz方程與Shuey二項(xiàng)近似方程、Shuey三項(xiàng)近似方程在該段入射角范圍內(nèi)區(qū)別不是很明顯。傳播矩陣算法與實(shí)際數(shù)據(jù)吻合度比Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程的計(jì)算結(jié)果更高，如T1標(biāo)志層所示，傳播矩陣方法中模擬出與實(shí)際數(shù)據(jù)一致的反射振幅隨入射角減弱的現(xiàn)象，而Zoeppritz方程計(jì)算結(jié)果中振幅隨入射角增強(qiáng)，驗(yàn)證了傳播矩陣?yán)碚撛诶碚撋系耐陚湫浴Ｅc基于單界面反射的Zoeppritz方程算法相比，傳播矩陣?yán)碚摮浞挚紤]了地震波的傳播效應(yīng)，正演結(jié)果與實(shí)際數(shù)據(jù)更接近，為進(jìn)一步研究非均勻儲(chǔ)層地震響應(yīng)提供了精確的模擬方法。為了更加直觀地顯示以上幾種正演算法之間的差異，對(duì)4種正演算法進(jìn)行歸一化并計(jì)算出差剖面來進(jìn)行誤差分析。如圖8所示，傳播矩陣算法在T1、T2標(biāo)志層處振幅小于其他幾種算法。振幅趨勢(shì)與實(shí)際數(shù)據(jù)振幅趨勢(shì)一致。

通過以上分析結(jié)果，進(jìn)一步對(duì)實(shí)際數(shù)據(jù)得到的地震剖面與傳播矩陣算法、Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程得到的地震剖面進(jìn)行對(duì)比分析。圖9為T1、T2處提取的振幅曲線。由圖9a、b可知，在T1標(biāo)志層處，傳播矩陣算法得到的振幅曲線趨勢(shì)與實(shí)際數(shù)據(jù)得到的振幅曲線趨勢(shì)基本一致，而Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程沒能合理預(yù)測(cè)振幅的變化趨勢(shì)；由圖9c、d可知，在T2標(biāo)志層處，這4種算法得到的振幅曲線趨勢(shì)均與實(shí)際數(shù)據(jù)得到的振幅曲線基本一致，但也存在一定差異。正演算法得到的振幅提取曲線與實(shí)際數(shù)據(jù)得到的振幅提取曲線趨勢(shì)吻合度越高，說明正演結(jié)果越好。Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程在T1處預(yù)測(cè)失誤可能由于選取的研究區(qū)上覆層存在較強(qiáng)的反射界面，而Zoeppritz方程是基于單界面的反射，沒有考慮地震波傳播時(shí)層間多次波傳播效應(yīng)和調(diào)諧干涉效應(yīng)。

圖5 頁巖氣儲(chǔ)層測(cè)井曲線(a)、頻變反射系數(shù)譜(b)和反射波頻譜(c)Fig.5 Shale reservoirs loging curve(a), reflection coefficients spectrum (b) and amplitude-spectrum of reflection wave (c) in frequency domain of shale reservoirs

圖6 頁巖氣儲(chǔ)層井震標(biāo)定Fig.6 Well to seismic calibration of shale reservoirs data

a.實(shí)際地震數(shù)據(jù)；b. 傳播矩陣算法；c. Zoeppritz方程；d. Shuey二項(xiàng)近似方程；e. Shuey三項(xiàng)近似方程。圖7 頁巖氣儲(chǔ)層合成地震記錄Fig.7 Synthetic seismogram shale reservoirs data

a.Zoeppritz方程與傳播矩陣算法差剖面；b. Shuey二項(xiàng)近似方程與傳播矩陣算法差剖面；c. Shuey三項(xiàng)近似方程與傳播矩陣算法差剖面。圖8 誤差分析Fig.8 Error analysis

a.實(shí)際數(shù)據(jù)T1處振幅曲線和平滑后振幅曲線；b. 4種算法T1處振幅曲線；c.實(shí)際數(shù)據(jù)T2處振幅曲線和平滑后振幅曲線；d. 4種算法T2處振幅曲線。圖9 振幅曲線Fig.9 Amplitude-angle curve

3 結(jié)論

1)本文應(yīng)用傳播矩陣、Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程4種正演算法，計(jì)算砂泥巖薄互層和頁巖儲(chǔ)層非均勻地層的高精度地震合成記錄。應(yīng)用于薄互層實(shí)際數(shù)據(jù)的計(jì)算結(jié)果表明，與Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程相比，傳播矩陣方法考慮了地震反射的傳播效應(yīng)、調(diào)諧效應(yīng)等動(dòng)力學(xué)特征，可以更加清晰地反映薄互層結(jié)構(gòu)的AVO特征變化，并為薄層、薄互層條件下的地震反射振幅信息提供了更為準(zhǔn)確的正演方法。

2)將上述4種正演方法應(yīng)用于頁巖儲(chǔ)層進(jìn)行對(duì)比研究得出，在近垂直入射時(shí)(0°～20°)，Zoeppritz方程算法及Shuey二項(xiàng)近似方程、Shuey三項(xiàng)近似方程與傳播矩陣的結(jié)果基本相同，但是在中、遠(yuǎn)偏移距(入射角>20°)的情況下，傳播矩陣算法與實(shí)際觀測(cè)數(shù)據(jù)吻合度最高。因此，基于傳播矩陣?yán)碚撚?jì)算高精度合成地震記錄，可以充分考慮地震波動(dòng)力學(xué)特征，為理解非均勻地層地震反射特征、進(jìn)一步開發(fā)AVO反演和分析技術(shù)提供方法。