基于巖石物理的頁(yè)巖儲(chǔ)層各向異性表征

劉　財(cái)　鄧馨卉　郭智奇*　劉喜武　劉宇巍

(①吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林長(zhǎng)春 130026； ②頁(yè)巖油氣富集機(jī)理與有效開(kāi)發(fā)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100083；③中國(guó)石化頁(yè)巖油氣勘探開(kāi)發(fā)重點(diǎn)實(shí)驗(yàn)室，北京 100083； ④中國(guó)石化石油勘探開(kāi)發(fā)研究院，北京 100083)

1　引言

由于頁(yè)巖礦物組分和微觀(guān)結(jié)構(gòu)復(fù)雜，因此頁(yè)巖儲(chǔ)層評(píng)價(jià)對(duì)于地球物理方法和技術(shù)的要求很高[1-5]。天然發(fā)育的微觀(guān)孔隙和裂縫是頁(yè)巖儲(chǔ)層中油氣儲(chǔ)存和運(yùn)移的通道，對(duì)儲(chǔ)層評(píng)價(jià)以及甜點(diǎn)區(qū)識(shí)別具有重要意義。

頁(yè)巖的微觀(guān)物性特征一般表現(xiàn)為宏觀(guān)VTI性質(zhì)。黏土礦物的定向排列和水平縫的存在是影響儲(chǔ)層VTI特性的主要因素[6-10]。目前，對(duì)黏土礦物的研究較多，而針對(duì)頁(yè)巖水平縫的巖石物理研究尚不充分。

前人主要通過(guò)分析巖石物理模型正演結(jié)果的各向異性特征描述頁(yè)巖中的裂縫特性。Hudson[11，12]和Schoenberg[13]通過(guò)建立等效介質(zhì)巖石物理模型描述裂縫介質(zhì)的地震各向異性特征；Schoenberg等[14，15]在長(zhǎng)波長(zhǎng)條件下，研究了在具有水平縫的各向異性介質(zhì)中彈性波通過(guò)滑動(dòng)裂縫界面的響應(yīng)特征，認(rèn)為裂縫介質(zhì)的地震響應(yīng)可以等效為各向同性背景和水平縫引起的各向異性擾動(dòng)的綜合作用結(jié)果；Liu等[16]應(yīng)用地震各向異性理論描述裂縫特征取得了較好效果。

在儲(chǔ)層各向異性與地震波傳播特征的研究方面，Thomsen[17]提出了弱各向異性理論，并用五個(gè)參數(shù)表征各向異性的影響； Ruger[18，19]研究了VTI介質(zhì)的反射特征，并推導(dǎo)了VTI介質(zhì)反射系數(shù)近似方程； Carcione[20]研究了頁(yè)巖各向異性參數(shù)隨孔隙壓力和干酪根含量的變化規(guī)律； Ran等[21]通過(guò)重構(gòu)裂縫儲(chǔ)層的彈性參數(shù)和各向異性參數(shù)，實(shí)現(xiàn)了彈性參數(shù)和裂縫參數(shù)的疊前地震反演。

本文針對(duì)羅家地區(qū)A井測(cè)井?dāng)?shù)據(jù)，應(yīng)用Backus平均方法[22]將模型從測(cè)井尺度粗化至地震尺度，應(yīng)用各向異性傳播矩陣模擬頁(yè)巖含油儲(chǔ)層的地震AVO響應(yīng)，并結(jié)合實(shí)際地震資料進(jìn)行井震標(biāo)定。之后，基于VTI介質(zhì)反射系數(shù)理論公式[19]進(jìn)行VTI介質(zhì)參數(shù)AVO反演，并應(yīng)用井中巖石物理反演結(jié)果進(jìn)一步計(jì)算頁(yè)巖含油儲(chǔ)層的裂縫空間分布，反演結(jié)果證明了文中方法的有效性。

2　各向異性AVO響應(yīng)模擬

2.1　Backus平均理論

Backus[22]將層狀介質(zhì)在長(zhǎng)波長(zhǎng)條件下等效為均勻TI介質(zhì)。該理論是一種對(duì)測(cè)井曲線(xiàn)進(jìn)行尺度粗化的有效理論，同時(shí)考慮了層狀介質(zhì)在長(zhǎng)波長(zhǎng)條件下的各向異性特征[23，24]。因此，在地震數(shù)據(jù)分析過(guò)程中，為了有效進(jìn)行測(cè)井尺度粗化，目前主要采用Backus平均方法，Schoenberg等[14]進(jìn)一步將其推廣到具有任意各向異性的層狀介質(zhì)情況。Kumar[25]將Backus平均理論用于具有低對(duì)稱(chēng)各向異性的裂縫型頁(yè)巖含油儲(chǔ)層中。對(duì)于VTI介質(zhì)，由Backus平均理論計(jì)算得到地震尺度的等效VTI介質(zhì)的五個(gè)獨(dú)立的模量

(1)

2.2　各向異性傳播矩陣?yán)碚?/h3>

地震尺度的模型可等效為由一系列薄層組成的互層結(jié)構(gòu)，由于干涉、調(diào)諧等現(xiàn)象的存在，其反射波呈復(fù)合波形。影響薄互層介質(zhì)地震反射特征的因素有入射角、物性差異、入射波頻率、地層厚度、薄互層結(jié)構(gòu)等。Rokhlin等[26]給出了層狀介質(zhì)反射系數(shù)表達(dá)式。Carcione[27]進(jìn)一步提出了用于計(jì)算薄互層反射系數(shù)的廣義傳播矩陣?yán)碚摗８鶕?jù)該理論，對(duì)于P波入射，在薄互層介質(zhì)中的反射、透射系數(shù)可以表示為

r=[RPP,RPS,TPP,TPS]T

(2)

式中：RPP、RPS、TPP、TPS分別為PP波、PS波的反射、透射系數(shù)；A1與A2分別為與上、下層介質(zhì)物性參數(shù)有關(guān)的傳播矩陣；Bα(α=1,…,N)為具有N層結(jié)構(gòu)的中間薄互層的傳播矩陣；iP為P波入射向量。詳細(xì)公式與計(jì)算流程見(jiàn)附錄A。

2.3　Ruger理論

基于反射界面兩側(cè)的彈性參數(shù)的弱不連續(xù)性假設(shè)，Ruger[19]提出了VTI介質(zhì)PP波反射系數(shù)近似公式

=P+Bsin2θ+Asin2θtan2θ

(3)

(4)

界面兩側(cè)彈性參數(shù)的弱不連續(xù)性和弱各向異性由Thomsen參數(shù)

(5)

所表示。式(3)中的弱各向異性VTI介質(zhì)的Thomsen參數(shù)表達(dá)式為

(6)

(7)

(8)

(9)

(10)

式中cij為目標(biāo)介質(zhì)的彈性剛度系數(shù)。

3　水平縫密度疊前地震反演

根據(jù)上述VTI介質(zhì)反射系數(shù)理論，本文提出VTI介質(zhì)各向異性參數(shù)反演方法。式(3)的矩陣形式表達(dá)為

(11)

則

(12)

式中：P為與聲阻抗差異相關(guān)的零角度反射系數(shù);B反映了非零角度橫波速度的影響；A為近臨界角的振幅曲率。

根據(jù)Ruger的VTI介質(zhì)反射系數(shù)理論公式，P、B和A分別為

(13)

式中VP01和VP02分別為目標(biāo)反射界面兩側(cè)的P波垂直方向速度。當(dāng)假設(shè)目標(biāo)層反射界面兩側(cè)密度比趨于穩(wěn)定不變時(shí)，P可以近似為

(14)

由

Δε=ε2-ε1≈2(A-P)

(15)

可以得到目標(biāo)層兩側(cè)P波各向異性參數(shù)的變化。當(dāng)下伏巖層為各向同性時(shí)，上覆巖層的各向異性參數(shù)可表達(dá)為

ε1≈2(P-A)

(16)

因此通過(guò)計(jì)算P和A，可求得目標(biāo)層頁(yè)巖的P波各向異性參數(shù)，并根據(jù)井中巖石物理反演得到的各向異性參數(shù)與水平縫密度之間的關(guān)系進(jìn)一步計(jì)算頁(yè)巖水平縫的空間分布。

4　應(yīng)用實(shí)例

4.1　頁(yè)巖油儲(chǔ)層測(cè)井?dāng)?shù)據(jù)分析

巖心觀(guān)測(cè)表明，頁(yè)巖含油儲(chǔ)層裂縫系統(tǒng)發(fā)育。圖1為羅家地區(qū)A井頁(yè)巖含油儲(chǔ)層的測(cè)井曲線(xiàn)。由圖可見(jiàn)： 在約3060m深度處存在彈性性質(zhì)突變界面，測(cè)井解釋表明該界面以上為富有機(jī)質(zhì)頁(yè)巖含油儲(chǔ)層，具有高電阻率、低密度(圖1c)和低速度(圖1b)特征，界面以下為泥質(zhì)灰?guī)r，具有相對(duì)較高的速度和密度以及較低的孔隙度(圖1d)； 在新、老地層過(guò)渡位置，隨深度增加儲(chǔ)層巖石中方解石和白云石含量增加，而黏土和石英含量減少(圖1e)，導(dǎo)致伽馬值逐漸降低。

圖1　羅家地區(qū)A井頁(yè)巖含油儲(chǔ)層測(cè)井?dāng)?shù)據(jù)(a)自然伽馬； (b)縱波速度； (c)密度； (d)孔隙度； (e)礦物體積含量

圖2為羅家地區(qū)A井頁(yè)巖含油儲(chǔ)層的巖心掃描電鏡(SEM)與地層微電阻掃描成像(FMI)結(jié)果，可觀(guān)察到頁(yè)巖在水平背景下水平縫較發(fā)育。

4.2　正演模擬及井震標(biāo)定

圖3為羅家地區(qū)A井及等效介質(zhì)測(cè)井曲線(xiàn)。將尺度粗化后的數(shù)據(jù)作為輸入模型，應(yīng)用各向異性傳播矩陣?yán)碚撚?jì)算合成地震記錄，使用主頻為25Hz的雷克子波，目標(biāo)層對(duì)應(yīng)的時(shí)間約為2408ms。圖4為正演模擬與井震標(biāo)定結(jié)果，可見(jiàn)目標(biāo)層位(紅線(xiàn))下方對(duì)應(yīng)正反射系數(shù)。

圖2　羅家地區(qū)A井頁(yè)巖含油儲(chǔ)層的SEM(a)與FMI(b)結(jié)果

圖3　羅家地區(qū)A井及等效介質(zhì)測(cè)井曲線(xiàn)

圖4　正演模擬與井震標(biāo)定結(jié)果

4.3　基于Ruger理論的水平縫密度反演

圖5為過(guò)A井的疊前AVO角道集，圖6為目標(biāo)層時(shí)窗范圍，圖7為過(guò)A井剖面目標(biāo)層AVO曲線(xiàn)。由圖可見(jiàn)，目標(biāo)層具有I類(lèi)AVO特征，即零炮檢距振幅強(qiáng)且為正極性，AVO呈遞減趨勢(shì)，當(dāng)入射角足夠大時(shí)可觀(guān)察到極性反轉(zhuǎn)。圖8為根據(jù)A井的數(shù)據(jù)設(shè)計(jì)的目標(biāo)層反射模型。由圖可見(jiàn)，反射模型為以目標(biāo)層為反射界面的單界面層狀模型，上覆層為低速VTI頁(yè)巖，下伏層為高速各向同性灰?guī)r。

圖5　過(guò)A井的疊前AVO角道集

圖6　反演時(shí)窗范圍選擇

圖7　過(guò)A井剖面目標(biāo)層AVO曲線(xiàn)

圖8　根據(jù)A井的數(shù)據(jù)設(shè)計(jì)的目標(biāo)層反射模型

圖9　不同ε值的反射振幅隨入射角變化曲線(xiàn)

應(yīng)用Ruger理論，改變P波各向異性參數(shù)ε值計(jì)算該反射模型的AVO反射振幅曲線(xiàn)，得到不同ε值的反射振幅隨入射角的變化曲線(xiàn)(圖9)，可見(jiàn)隨著ε的增加，反射振幅在遠(yuǎn)炮檢距逐漸減小，為頁(yè)巖各向異性參數(shù)反演提供了依據(jù)。地震反演方法首先要求解式(12)中的超定矩陣。圖10為反演的目標(biāo)層參數(shù)。由圖可見(jiàn)，在A(yíng)、B和 C三口井位處，各向異性參數(shù)ε分別為0.0923、0.3038和0.2526(圖10c)，水平縫密度分別為0.0428、0.0072和0.0547(圖10d)。

圖11為各向異性參數(shù)反演結(jié)果，圖12為各向異性參數(shù)ε與水平縫密度εH交會(huì)圖，由ε、εH的關(guān)系可獲得目標(biāo)層水平縫密度的空間分布(圖10d)。 油氣開(kāi)發(fā)資料表明，羅家地區(qū)C井的產(chǎn)量大于A(yíng)井和B井，與C井的水平縫密度最大相吻合(圖10d)，證明利用反演的水平縫密度進(jìn)行儲(chǔ)層預(yù)測(cè)的有效性。

圖10　反演的目標(biāo)層參數(shù)

圖11　各向異性參數(shù)反演結(jié)果[28]

圖12　各向異性參數(shù)ε與水平縫密度εH交會(huì)圖

5　結(jié)束語(yǔ)

本文針對(duì)頁(yè)巖油儲(chǔ)層裂縫進(jìn)行地震響應(yīng)模擬及分析，基于VTI介質(zhì)反射理論，提出了頁(yè)巖含油儲(chǔ)層各向異性參數(shù)AVO反演方法，并結(jié)合巖石物理分析結(jié)果，進(jìn)一步計(jì)算儲(chǔ)層裂縫的空間分布，得到以下認(rèn)識(shí)。

(1)應(yīng)用Backus平均方法將模型從測(cè)井尺度粗化至地震尺度，應(yīng)用各向異性傳播矩陣方法模擬及分析頁(yè)巖含油儲(chǔ)層的地震AVO響應(yīng)，并結(jié)合實(shí)際地震資料進(jìn)行井震標(biāo)定。同時(shí)，在Ruger理論計(jì)算中發(fā)現(xiàn)，隨著模型中各向異性參數(shù)ε數(shù)值的增大，遠(yuǎn)炮檢距反射振幅逐漸減小。

(2)由地震反演得到的各向異性參數(shù)和水平縫密度與巖石物理的計(jì)算結(jié)果一致，驗(yàn)證了反演方法的正確性。同時(shí)，由反演得到的水平縫的空間分布與油氣產(chǎn)量具有相關(guān)性，說(shuō)明該方法在頁(yè)巖油有利儲(chǔ)層預(yù)測(cè)中的有效性。

附錄A

式(2)中的傳播矩陣A1和A2分別為

(A-1)

A2=iω×

(A-2)

(A-3)

(A-4)

W=p55(γsx+βsz)

(A-5)

Z=βp13sx+γp33sz

(A-6)

而

(A-7)

式中：變量ρ、γ、W和Z的下標(biāo)P、S分別對(duì)應(yīng)準(zhǔn)縱波、準(zhǔn)橫波，1、2分別對(duì)應(yīng)上、下層介質(zhì)； pv(·)意為取復(fù)數(shù)的主值； 參數(shù)pij為6×6階彈性剛度系數(shù)矩陣的系數(shù)。對(duì)于γ，符號(hào)“+”對(duì)應(yīng)準(zhǔn)P波，符號(hào)“-”對(duì)應(yīng)準(zhǔn)S波；sx為水平波慢度

(A-8)

其中

E={[(p33-p55)cos2θ-(p11-p55)sin2θ]2+

(A-9)

sz為垂直波慢度

(A-10)

(+,-)對(duì)應(yīng)向下傳播準(zhǔn)P波， (+,+)對(duì)應(yīng)向下傳播準(zhǔn)S波， (-,-)對(duì)應(yīng)向上傳播準(zhǔn)P波， (-,+)對(duì)應(yīng)向上傳播準(zhǔn)S波。其中

(A-11)

(A-12)

(A-13)

傳播矩陣Bα=T(0)T-1(hα)，其中

(A-14)

P波入射向量為

iP=iω[βP1,γP1, -ZP1,-WP1]T

(A-15)

[1]丁文龍,許長(zhǎng)春,久凱等.頁(yè)巖裂縫研究進(jìn)展.地球科學(xué)進(jìn)展,2011,26(2):135-144.

Ding Wenlong,Xu Changchun,Jiu Kai et al.The research progress of shale fractures.Advances in Earth Science,2011,26(2):135-144.

[2]王健,石萬(wàn)忠,舒志國(guó)等.富有機(jī)質(zhì)頁(yè)巖TOC含量的地球物理定量化預(yù)測(cè).石油地球物理勘探,2016,51(3):596-604.

Wang Jian,Shi Wanzhong,Shu Zhiguo et al.TOC content quantitative prediction in organic-rich shale.OGP,2016,51(3):596-604.

[3]謝劍勇,魏建新,狄?guī)妥尩?頁(yè)巖各向異性參數(shù)間關(guān)系及其影響因素分析.石油地球物理勘探,2015,50(6):1141-1145.

Xie Jianyong,Wei Jianxin,Di Bangrang et al.Correlation of anisotropic parameters and their influence factors.OGP,2015,50(6):1141-1145.

[4]蘇建龍,屈大鵬,陳超等.疊前地震反演方法對(duì)比分析——焦石壩頁(yè)巖氣藏勘探實(shí)例.石油地球物理勘探,2016,51(3):581-588.

Su Jianlong,Qu Dapeng,Chen Chao et al.Application of pre-stack inversion technique for shale gas in Jiaoshiba gas field.OGP,2016,51(3):581-588.

[5]牛聰,劉志斌,王彥春等.應(yīng)用地球物理技術(shù)定量評(píng)價(jià)遼西凹陷沙河街組烴源巖.石油地球物理勘探,2017,52(1):131-137.

Niu Cong,Liu Zhibin,Wang Yanchun et al.Quantitative evaluation of Shahejie formation source rocks in Liaoxi Sag with geophysical approaches.OGP,2017,52(1):131-137.

[6]慈興華,劉宗林,王志戰(zhàn).羅家地區(qū)泥質(zhì)巖裂縫性?xún)?chǔ)集層綜合研究.錄井工程,2006,17(1):71-74.

Ci Xinghua,Liu Zonglin,Wang Zhizhan.Composite study on fractured reservoir of argillaceous rock in Luojia Area.Mud Logging Engineering,2006,17(1):71-74.

[7]尹克敏,李勇,慈興華等.羅家地區(qū)沙三段泥質(zhì)巖裂縫特征研究.斷塊油氣田,2002,9(5):24-27.

Yin Kemin,Li Yong,Ci Xinghua et al.Study on characters of muddy fractural reservoirs in Luojia Area.Fault-Block Oil & Gas Field,2002,9(5):24-27.

[8]Cheng C H.Seismic velocities in porous rocks: Direct and inverse problems.Massachusetts Institute of Technology,1978,42(1):143-144.

[9]Cheng C H.Crack models for a transversely anisotropic medium.Journal of Geophysical Research Atmospheres,1993,98(B1):675-684.

[10]Guo Z,Li X Y,Liu C.Anisotropy parameters estimate and rock physics analysis for the Barnett Shale.Journal of Geophysics & Engineering,2014,11(6):65006-65016.

[11]Hudson J A.Overall properties of a crack solid. Mathematical Proceedings of the Cambridge Philosophical Society,1980,88(2):371-384.

[12]Hudson J A.Wave speeds and attenuation of elastic waves in material containing cracks.Geophysical Journal International,1981,64(1):133-150.

[13]Schoenberg M.Reflection of elastic waves from periodically stratified media with interfacial slip.Geophysical Prospecting,1983,31(2):265-292.

[14]Schoenberg M,Muir F.A calculus for finely layered anisotropic media.Geophysics,1989,54(5):581-589.

[15]Schoenberg M,Sayers C M.Seismic anisotropy of fractured rock.Geophysics,1995,60(1):204-211.

[16]Liu E,Martinez A.Seismic Fracture Characterization.EAGE Publication,Netherlands,2012.

[17]Thomsen L.Weak elastic anisotropy.Geophysics,1986,51(10):1954-1966.

[18]Rüger A.Reflection Coefficients and Azimuthal AVO Analysis in Anisotropic Media[D].Colorado School of Mines,Colorado,1996.

[19]Rüger A.P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry.Geophysics,1997,62(3):713-722.

[20]Carcione J M.A model for seismic velocity and attenuation in petroleum source rocks.Geophysics,2000,65(4):1080-1092.

[21]Ran B,Sengupta M,Salama A et al.Reconstruction of the layer anisotropic elastic parameters and high-resolution fracture characterization from P-wave data:a case study using seismic inversion and Bayesian rock physics parameter estimation.Geophysical Prospecting,2009,57(2):253-262.

[22]Backus G E.Long-wave elastic anisotropy produced by horizontal layering.Journal of Geophysics Research,1962,67(11):4427-4440.

[23]Bayuk I O,Ammerman M,Chesnokov E M.Upscaling of elastic properties of anisotropic sedimentary rocks.Geophysical Journal of the Royal Astronomical Society,2008,172(2):842-860.

[24]Lindsay R,Koughnet R V.Sequential Backus averaging:Upscaling well logs to seismic wavelengths.The Leading Edge,2001,20(2):188-191.

[25]Kumar D.Applying Backus averaging for deriving seismic anisotropy of a long-wavelength equivalent medium from well-log data.Journal of Geophysics and Engineering,2013,10(9):1-15.

[26]Rokhlin S I,Wang Y J.Equivalent boundary condi-tions for thin orthotropic layer between two solids:reflection,refraction,and interface waves.Journal of the Acoustical Society of America,1992,91(1):1875-1887.

[27]Carcione J M.AVO effects of a hydrocarbon source-rock layer.Geophysics,2001,66(2):419-427.

[28]Guo Z Q,Liu C,Liu X W et al.Research on anisotropy of shale oil reservoir based on rock physics model.Applied Geophysics,2016,13(2):382-392.