劉 財(cái),裴思嘉,郭智奇,符 偉,張宇生,劉喜武
1.吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院,長春 130026 2.中石油東方地球物理公司,河北 涿州 072751 3.中國石化頁巖油氣勘探開發(fā)重點(diǎn)實(shí)驗(yàn)室,北京 100083 4.中國石化石油勘探開發(fā)研究院,北京 100083
地震反射系數(shù)反映地下界面物性特征的變化。對于各向同性介質(zhì)分界面,平面波反射、透射系數(shù)的精確解可由Zoeppritz[1]方程計(jì)算。許多學(xué)者對Zoeppritz方程進(jìn)行了不同形式的簡化,以此開展地震AVO(amplitude versus offset)技術(shù)研究,通過反射波振幅隨偏移距的變化規(guī)律研究反射界面兩側(cè)介質(zhì)的物性特征。傳統(tǒng)的AVO分析都是基于Zoeppritz方程的各自線性近似式提出的[2-4]。1955年,Koefoed[5]不僅簡化了Zoeppritz方程,而且在實(shí)驗(yàn)中假設(shè)上下界面的泊松比不同,所得的縱波反射系數(shù)隨入射角的改變而變化范圍大大增大,這對之后的AVO技術(shù)發(fā)展具有重大意義。1980年,Aki等[6]將Zoeppritz方程簡化為入射角與縱波振幅相關(guān)的近似公式。1985年,Shuey[7]根據(jù)該簡化公式做了進(jìn)一步研究,將其改寫為入射角小角度項(xiàng)、中角度項(xiàng)和廣角度項(xiàng)三部分之和的近似公式。
然而,層狀介質(zhì)的AVO特征與單一界面是不同的,因?yàn)榉瓷湎禂?shù)與入射波的頻率相關(guān)。實(shí)際地層為非均勻介質(zhì),地震波在地下傳播時(shí)具有傳播效應(yīng)和層間多次波的調(diào)諧干涉效應(yīng)。對此,早在1971年,F(xiàn)uchs等[8]就應(yīng)用反射率方法計(jì)算了多層介質(zhì)對點(diǎn)震源的響應(yīng)。1980年,Aki等[6]提出傳播矩陣算法,在理論上可以更好地描述這種層狀介質(zhì)模型。1983年,Kennett[9]給出了地震波在層狀各項(xiàng)異性介質(zhì)中傳播的描述。1984年,F(xiàn)ryer等[10]又將傳播矩陣用于各向異性介質(zhì),制作了合成地震記錄。在此之后,有許多學(xué)者對層狀介質(zhì)進(jìn)行了AVO正演響應(yīng)研究。2001年,Carcione[11]計(jì)算了富有機(jī)質(zhì)頁巖層的反射系數(shù)。2003年,Liu等[12]研究了泊松比和層厚度對薄層AVO響應(yīng)的影響。2005—2016年,郭智奇等[13-15]基于廣義傳播矩陣?yán)碚撚?jì)算地震波頻變反射系數(shù),完善了頻變AVO算法,為含油氣儲層頻變AVO響應(yīng)的模擬和分析提供了方法。2014年,蘭慧田[16]應(yīng)用廣義傳播矩陣?yán)碚撚?jì)算了上覆介質(zhì)與裂縫性孔隙介質(zhì)分界面的反射系數(shù),為裂縫檢測和裂縫參數(shù)反演提供地震波振幅信息。
本文針對非均勻地下介質(zhì)模型,選取薄互層結(jié)構(gòu)和頁巖氣儲層兩組實(shí)際測井?dāng)?shù)據(jù),分別比較精確的Zoeppritz方程、Shuey二項(xiàng)近似方程和Shuey三項(xiàng)近似方程以及傳播矩陣算法得到的正演模擬響應(yīng)計(jì)算效果,并與實(shí)際的地震數(shù)據(jù)進(jìn)行對比分析,討論存在調(diào)諧干涉和傳播效應(yīng)的情況下不同方法的精確性和適用性。
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í),將會產(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
Shuey[7]將Zoeppritz方程簡化為關(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)近似方程。
對于非均勻地下介質(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],對于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ù)的主值;對于γ,符號“+”對應(yīng)準(zhǔn)縱波,符號“-”對應(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á)式符號約定為:(+,-)表示向下傳播準(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)
對Uf做反傅里葉變換則可以獲得時(shí)間域的反射波波形ut:
(25)
式中:f為頻率;t為時(shí)間。
圖2a為砂泥巖薄互層結(jié)構(gòu)的實(shí)際數(shù)據(jù)測井曲線,包括縱、橫波速度和密度。基于傳播矩陣?yán)碚撚?jì)算該非均勻地層的地震響應(yīng),圖2b為計(jì)算得到的頻變反射系數(shù)??紤]入射波為40 Hz的Ricker子波,進(jìn)一步計(jì)算反射波頻譜如圖2c所示。
對圖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方程更為接近,對于此數(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 砂泥巖薄互層儲層測井曲線(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
圖5a所示為某頁巖氣儲層的縱、橫波速度和密度測井曲線。從圖5a中可以看到,縱、橫波速度曲線的相關(guān)性很強(qiáng),圖中2 330~2 415 m為頁巖地層,含氣儲層為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í)際觀測數(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í)有明顯振幅減小的趨勢,而其他方法振幅沒有該趨勢。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)碚撛诶碚撋系耐陚湫?。與基于單界面反射的Zoeppritz方程算法相比,傳播矩陣?yán)碚摮浞挚紤]了地震波的傳播效應(yīng),正演結(jié)果與實(shí)際數(shù)據(jù)更接近,為進(jìn)一步研究非均勻儲層地震響應(yīng)提供了精確的模擬方法。為了更加直觀地顯示以上幾種正演算法之間的差異,對4種正演算法進(jìn)行歸一化并計(jì)算出差剖面來進(jìn)行誤差分析。如圖8所示,傳播矩陣算法在T1、T2標(biāo)志層處振幅小于其他幾種算法。振幅趨勢與實(shí)際數(shù)據(jù)振幅趨勢一致。
通過以上分析結(jié)果,進(jìn)一步對實(shí)際數(shù)據(jù)得到的地震剖面與傳播矩陣算法、Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程得到的地震剖面進(jìn)行對比分析。圖9為T1、T2處提取的振幅曲線。由圖9a、b可知,在T1標(biāo)志層處,傳播矩陣算法得到的振幅曲線趨勢與實(shí)際數(shù)據(jù)得到的振幅曲線趨勢基本一致,而Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程沒能合理預(yù)測振幅的變化趨勢;由圖9c、d可知,在T2標(biāo)志層處,這4種算法得到的振幅曲線趨勢均與實(shí)際數(shù)據(jù)得到的振幅曲線基本一致,但也存在一定差異。正演算法得到的振幅提取曲線與實(shí)際數(shù)據(jù)得到的振幅提取曲線趨勢吻合度越高,說明正演結(jié)果越好。Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程在T1處預(yù)測失誤可能由于選取的研究區(qū)上覆層存在較強(qiáng)的反射界面,而Zoeppritz方程是基于單界面的反射,沒有考慮地震波傳播時(shí)層間多次波傳播效應(yīng)和調(diào)諧干涉效應(yīng)。
圖5 頁巖氣儲層測井曲線(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 頁巖氣儲層井震標(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 頁巖氣儲層合成地震記錄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
1)本文應(yīng)用傳播矩陣、Zoeppritz方程、Shuey二項(xiàng)近似方程及Shuey三項(xiàng)近似方程4種正演算法,計(jì)算砂泥巖薄互層和頁巖儲層非均勻地層的高精度地震合成記錄。應(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)用于頁巖儲層進(jìn)行對比研究得出,在近垂直入射時(shí)(0°~20°),Zoeppritz方程算法及Shuey二項(xiàng)近似方程、Shuey三項(xiàng)近似方程與傳播矩陣的結(jié)果基本相同,但是在中、遠(yuǎn)偏移距(入射角>20°)的情況下,傳播矩陣算法與實(shí)際觀測數(shù)據(jù)吻合度最高。因此,基于傳播矩陣?yán)碚撚?jì)算高精度合成地震記錄,可以充分考慮地震波動(dòng)力學(xué)特征,為理解非均勻地層地震反射特征、進(jìn)一步開發(fā)AVO反演和分析技術(shù)提供方法。