劉　財　鄧馨卉　郭智奇*　劉喜武　劉宇巍

(①吉林大學地球探測科學與技術學院，吉林長春 130026； ②頁巖油氣富集機理與有效開發(fā)國家重點實驗室，北京 100083；③中國石化頁巖油氣勘探開發(fā)重點實驗室，北京 100083； ④中國石化石油勘探開發(fā)研究院，北京 100083)

1　引言

由于頁巖礦物組分和微觀結構復雜，因此頁巖儲層評價對于地球物理方法和技術的要求很高[1-5]。天然發(fā)育的微觀孔隙和裂縫是頁巖儲層中油氣儲存和運移的通道，對儲層評價以及甜點區(qū)識別具有重要意義。

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

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

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

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

2　各向異性AVO響應模擬

2.1　Backus平均理論

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

(1)

2.2　各向異性傳播矩陣理論

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

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

(2)

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

2.3　Ruger理論

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

=P+Bsin2θ+Asin2θtan2θ

(3)

(4)

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

(5)

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

(6)

(7)

(8)

(9)

(10)

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

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

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

(11)

則

(12)

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

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

(13)

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

(14)

由

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

(15)

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

ε1≈2(P-A)

(16)

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

4　應用實例

4.1　頁巖油儲層測井數(shù)據(jù)分析

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

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

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

4.2　正演模擬及井震標定

圖3為羅家地區(qū)A井及等效介質測井曲線。將尺度粗化后的數(shù)據(jù)作為輸入模型，應用各向異性傳播矩陣理論計算合成地震記錄，使用主頻為25Hz的雷克子波，目標層對應的時間約為2408ms。圖4為正演模擬與井震標定結果，可見目標層位(紅線)下方對應正反射系數(shù)。

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

圖3　羅家地區(qū)A井及等效介質測井曲線

圖4　正演模擬與井震標定結果

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

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

圖5　過A井的疊前AVO角道集

圖6　反演時窗范圍選擇

圖7　過A井剖面目標層AVO曲線

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

圖9　不同ε值的反射振幅隨入射角變化曲線

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

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

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

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

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

5　結束語

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

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

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

附錄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的下標P、S分別對應準縱波、準橫波，1、2分別對應上、下層介質； pv(·)意為取復數(shù)的主值； 參數(shù)pij為6×6階彈性剛度系數(shù)矩陣的系數(shù)。對于γ，符號“+”對應準P波，符號“-”對應準S波；sx為水平波慢度

(A-8)

其中

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

(A-9)

sz為垂直波慢度

(A-10)

(+,-)對應向下傳播準P波， (+,+)對應向下傳播準S波， (-,-)對應向上傳播準P波， (-,+)對應向上傳播準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]丁文龍,許長春,久凱等.頁巖裂縫研究進展.地球科學進展,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]王健,石萬忠,舒志國等.富有機質頁巖TOC含量的地球物理定量化預測.石油地球物理勘探,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ī)妥尩?頁巖各向異性參數(shù)間關系及其影響因素分析.石油地球物理勘探,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]蘇建龍,屈大鵬,陳超等.疊前地震反演方法對比分析——焦石壩頁巖氣藏勘探實例.石油地球物理勘探,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]牛聰,劉志斌,王彥春等.應用地球物理技術定量評價遼西凹陷沙河街組烴源巖.石油地球物理勘探,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ū)泥質巖裂縫性儲集層綜合研究.錄井工程,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ū)沙三段泥質巖裂縫特征研究.斷塊油氣田,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.