基于Tsvankin参数页岩储层各向异性强度预测

李 翔, 尹 陈, 吴永宏, 熊晓军, 刘 阳, 简世凯

(1.油气藏地质及开发工程国家重点实验室 (成都理工大学)，成都 610059;2.中国石油川庆钻探工程有限公司 地球物理勘探公司，成都 610213 3.中国石油西南油气田分公司 重庆气矿，重庆 400021)

基于Tsvankin参数页岩储层各向异性强度预测

李 翔1, 尹 陈2, 吴永宏3, 熊晓军1, 刘 阳1, 简世凯1

(1.油气藏地质及开发工程国家重点实验室 (成都理工大学)，成都 610059;2.中国石油川庆钻探工程有限公司 地球物理勘探公司，成都 610213 3.中国石油西南油气田分公司 重庆气矿，重庆 400021)

在页岩气勘探开发过程中，各向异性强弱与储层品质、地震波速度、岩石脆性、岩石弹性参数等都有关系，所以对页岩各向异性特征的定量描述具有重要意义。综合考虑页岩页片状层理大量发育，基于VTI介质模型对页岩进行模拟，通过对Thomsen参数简化而利用Tsvankin参数进行页岩各向异性表征，然后以叠前偏移结果为基准，修正均方根速度使得叠前偏移后的CRP道集层拉平效果明显，迭代计算出使得叠前偏移效果最优时的各向异性强度的Tsvankin参数值，并利用Tsvankin参数对叠前道集剩余时差进行校正，最后输出Tsvankin参数数据体，分析Tsvankin参数与页岩层脆性发育及储层性质的关系。将该方法应用于四川长宁“宁201”井区的下志留统龙马溪组页岩，综合地震资料分析发现Tsvankin参数谱值域范围与页岩脆性概率分布趋势正相关，能够预测页岩储层各向异性强度并能反映页岩地层脆性强弱，在“甜点”预测中有一定的指导意义。

Tsvankin参数；页岩；各向异性；脆性

页岩气赋存于具有超低孔隙度的页岩中，是一种在海相、陆相及过渡相地层均有发育的很重要的非常规资源。页岩由于其构成矿物颗粒在成分、形状及大小上均有差异，加之结构上大量发育层理及片理，所以页岩储层各向异性特征十分明显。四川盆地作为中国最重要的页岩气勘探区域，探究各向异性特征对页岩储层“甜点”预测及后期开发意义重大。

为了探究页岩储层各向异性发育情况，邓继新等[1]在实验室超声条件下进行试验，探究了泥、页岩纵横波速度及各向异性参数随压力的变化过程；Carl等[2]探究了Thomsen参数的适用性并结合SEM资料研究了页岩微观物性对其各向异性特征的影响；王倩等[3]根据岩石物理实验探究了页岩弹性模量、泊松比与层理面的关系，设计了一种页岩各向异性指数；谢剑勇等[4]基于岩样测试研究了页岩各向异性参数间的关系并对影响各向异性的因素进行分析。这些研究证明了页岩由于其层理发育可利用横向各向同性(TI)介质进行模拟，在岩石物理领域，通过Thomsen参数评价页岩各向异性特征具有较好的适用性[5]。

以上学者从实验的角度对页岩各向异性特征进行研究；然而在页岩气勘探过程中，岩心资料珍贵而数量极少，实验分析较难以满足生产需求。如何从理论计算的角度进行各向异性参数预测是一个亟待解决的问题。张广智等[6]基于测井资料，采用岩石物理建模方法，反演计算了碳酸盐岩地区横波速度及各向异性参数；孙伟家等[7]从正演的角度对VTI介质强各向异性特征进行数值模拟。而如何充分利用地震资料进行各向异性参数预测的研究成果较少。本文基于页岩工区叠前地震资料，设计了一套针对页岩储层的Tsvankin各向异性参数计算方法。该参数能够对页岩工区各向异性强度进行定量描述，从而对勘探开发有一定的指导意义[8-16]。

1 方法原理

1.1 各向异性强度计算方法

页岩由粒径lt;0.003 9 mm的细粒碎屑、黏土及有机质组成，页岩中广泛发育页状及薄片状的层理，是一种典型的各向异性介质。我们利用表征薄层旋回沉积的VTI介质来模拟页岩储层页片状层理，VTI介质的弹性刚度矩阵为

(1)

VTI介质及HTI介质的弹性刚度矩阵均含5个独立的弹性参数。为便于理论研究，L.Thomsen[17]于1986年提出了一套新的弹性参数表达式，将5个弹性参数简化为3个，定义如下

(2)

其中：ε值与介质的纵波各向异性正相关，当ε=0时，介质的纵波无各向异性;δ是连接vP0和vP90(0和90表示入射角为0°和90°)之间的一种随入射角变化而变化的过渡参数。

综合考虑基于常规声波地震资料的纵波勘探方法，T.Alkhalifah等[18]重新定义了一个有效的各向异性参数“η”(内含Thomsen参数)，将弹性参数进一步简化为

(3)

在常规各向同性介质速度分析过程中，反射波时距曲线满足双曲线特征。TI介质中，时距曲线不再满足简单的双曲线规律，由T.Alkhalifah等[18]重新定义的时距方程

(4)

式中：x为偏移距；t0为零偏移距的双程旅行时；vNM0为入射角为0°时的动校正速度。(4)式充分地描述了长距离传播的纵波正常时差与零入射角时的动校正速度vNM0及各向异性参数η的关系。方程右端第三项不同于常规的双曲线时距方程，其主要控制时差方程中的非双曲线形态，若该式中η=0则该方程与常规双曲线时距方程一致。分母中增加x校正因子，是为了将泰勒展开式中的截断误差进行弥补。由(4)式能够看出，当介质各向异性较强且偏移距较大时，反射波时距关系会偏离各向同性双曲线关系。

由此，计算各向异性η参数可以通过进行速度分析得到均方根速度，并依此速度对地震数据进行叠前时间偏移。由于存在各向异性段的道集，其同相轴形态受(4)式中的η值影响，形态不会是标准的双曲线；因而在常规叠前偏移后其层拉平效果不明显。当偏移迭代计算结果变化不大且层拉平效果不明显时，基于(4)式引入η参数重新计算地震波时差，通过修改η值在CMP道集动校拉平效果最优时输出此时的各向异性η参数作为该道的各向异性强度检测值。

1.2 各向异性参数计算流程

综合考虑上述方法，本文设计了一套针对页岩工区叠前各向异性强度计算流程，其具体步骤如下。

步骤1： 基于页岩气地区三维叠前地震数据，对数据进行叠前道集保边去噪处理。通过引入最小标准差准则，在压制资料噪声的同时能很好地保留边缘等细节信息。在实际资料处理过程中发现该步骤对计算结果影响显著，基于保边去噪后品质较好的地震数据能计算得到更加稳定的η参数。

步骤2： 对叠前道集进行速度分析，计算RMS速度(均方根速度)，基于得到的RMS速度进行叠前时间偏移，并观察偏移后结果的道集拉平效果，可反复迭代计算至偏移后叠前道集位置变化不大时结束计算，作为计算η参数的基础。

步骤3： 给定初始的各向异性参数η场，将η值代入(4)式重新计算地震波走时。利用计算得到的非双曲线形态地震波走时对叠前道集剩余时差进行修正，观察修正后叠前道集层拉平效果，当拉平效果满足要求时输出。迭代计算完成整块三维叠前数据得到三维η数据体，作为工区各向异性强度预测结果。

图1为上述步骤的技术流程。

图1 各向异性参数η场计算流程Fig.1 The flowchart of calculation of anisotropic parameter field

通过以上3个步骤计算得到页岩工区的各向异性参数η数据体。由此对η数据体进行分析，绘制研究区目的层η平面展布图，可结合地震资料分析各向异性强度与地震属性的联系，并基于η参数进一步实现各向异性叠前偏移。

2 应用实例

四川长宁地区页岩气资源丰富，是中国石油天然气股份有限公司页岩气建产核心区域之一[19]。其中宁201-H1井日产气gt;0.15×106m3，其勘探目的层为下志留统龙马溪组页岩。

将本文提出的各向异性强度预测方法应用在长宁的“宁201”井区三维叠前地震数据中，基于常规叠前时间偏移结果为约束，迭代计算使得层拉平效果最优时的各向异性参数谱，并根据得到的η谱对叠前数据剩余时差进行校正。工区实践如图2所示，图2-A中计算所得的η参数谱能量越集中且偏离中线越远表示各向异性发育强度越大，图2-B为未经各向异性参数校正的叠前道集，图2-C为经各向异性校正后的叠前道集。对比发现：经过各向异性参数校正后，浅层远偏移距处的层拉平效果改善明显。与此同时，图2-A中拾取的η参数谱在1.4 s处偏离中线位置，且能量较强，对应各向异性参数校正前后的CRP道集剖面；如图2-B和图2-C中红色箭头所指，叠前道集同相轴的层拉平效果有一定程度的改善，说明各向异性参数对CRP道集剩余时差校正有积极作用，计算得到的各向异性参数数据体对后期的叠前偏移意义同样重大。

图2 各向异性参数对叠前CRP道集校正结果Fig.2 The results of anisotropic parameter for the pre-stack gathers(A)各向异性参数谱剖面; (B)各向异性参数校正前的叠前道集; (C)各向异性参数校正后的叠前道集

对三维工区计算其各向异性参数体，针对目的层绘制平面图(图3)。我们发现“宁201”井处于各向异性强度中高值段(图中红黄色段)，与实钻高产气结论吻合，说明各向异性发育强度与页岩储层品质有关，可以一定程度上指导储层“甜点”预测。龙马溪组页岩脆性概率呈西南高东北低的趋势，这与各向异性参数数据强弱趋势范围完全一致，说明各向异性强度与页岩脆性发育存在一定程度的正相关关系，图中红黄色段是脆性概率较大区域，是有利的“甜点”区。

图3 沿目的层向下50 ms各向异性参数平面图Fig.3 The layout of anisotropy parameter downward for 50 ms along the target horizon

工区正东Inline 476线(图3中黑线)叠后地震剖面如图4，目的层段龙马溪组(图中蓝线)地震同相轴连续性较差，可以观察到小断裂较为发育，说明此处地层普遍具有各向异性，与根据地震剩余时差计算所得的高各向异性参数结论吻合，该处页岩层物性较差，为非有效储层段。将Inline 476线投影在各向异性参数平面图上，可以看到其处于强各向异性发育段(图3中白色区域高值区)，证明各向异性参数对评价页岩储层物性有一定意义。

图4 目的层段龙马溪组页岩Inline 476线地震剖面Fig.4 The seismic profile along Inline 476 for Longmaxi Formation

3 结 论

a.页岩中广泛发育各向异性，利用Tsvankin各向异性参数可以独立地表征纵波在页岩中传播时各向异性特征强弱，通过迭代可以得到一个最优各向异性参数。该参数能够进一步校正常规叠前偏移后的剩余时差，为精细刻画页岩断裂及储层性质提供依据。

b.得到的各向异性参数数据体可以为页岩工区“甜点”预测提供一定的依据。将该方法应用到长宁地区，分析认为“宁201”井位于各向异性参数中高值区域，较优的储层发育着中等强度的各向异性特征。龙马溪组页岩脆性概率呈西南高东北低的趋势，与各向异性参数的场数据强弱趋势一致并呈现正相关关系。在页岩气勘探中各向异性较强的区域其岩石脆性较强，是有利的“甜点”区。

c.在工区正东部各向异性场很强，呈片状分布，结合地震资料分析得出：高各向异性可能由断裂广泛发育引起，地震同相轴连续性较差的地区，是勘探中应当回避的非“甜点”区。

[1] 邓继新,史謌,刘瑞珣,等. 泥岩、页岩声速各向异性及其影响因素分析[J].地球物理学报,2004,47(5):862-868.

Deng J X, Shi G, Liu R X,etal. Analysis of the velocity anisotropy and its affection factors in shale and mudstone[J]. Chinese Journal of Geophysics, 2004, 47(5): 862-868. (in Chinese)

[2] Carl H, Chandra S. Elastic anisotropy of shales[J]. The Leading Edge, 2011, 30(3): 324-331.

[3] 王倩,王鹏,项德贵,等.页岩力学参数各向异性研究[J].天然气工业,2012,32(12):62-65.

Wang Q, Wang P, Xiang D G,etal. Anisotropy property of mechanical parameters of shales[J]. Natural Gas Industry, 2012, 32(12): 62-65. (in Chinese)

[4] 谢剑勇,魏建新,狄帮让,等.页岩各向异性参数间关系及其影响因素分析[J].石油地球物理勘探,2015,50(6):1141-1145.

Xie J Y, Wei J X, Di B R,etal. Correlation of shale anisotropy parameters and their influence factors[J]. Oil Geophysical Prospecting, 2015, 50(6):1141-1145. (in Chinese)

[5] Mavko G, Mukerji T, Dvorkin J. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media[M]. Cambridge: Cambridge University Press, 2009.

[6] 张广智,陈怀震,王琪,等.基于碳酸盐岩裂缝岩石物理模型的横波速度和各向异性参数预测[J].地球物理学报,2013,56(5):1707-1715.

Zhang G Z, Chen H Z, Wang Q,etal. Prediction for in-situ formation stress of shale based on rock physics equivalent model[J]. Chinese Journal of Geophysics, 2013, 56(5): 1707-1715. (in Chinese)

[7] 孙伟家,符力耘,管西竹,等.页岩气地震勘探中页岩各向异性的地震模拟研究[J].地球物理学报,2013,56(3):961-970.

Sun W J, Fu L Y, Guan X Z,etal. A study on anisotropy of shale using seismic forward modeling in shale gas exploration[J]. Chinese Journal of Geophysics, 2013, 56(3): 961-970. (in Chinese)

[8] Meadows M, Coen S. Exact inversion of plane layered isotropic and anisotropic elastic media by the state-space approach [J]. Geophysics, 1986, 51(11): 2031-2050.

[9] Meadows M A, Abriel W L. 3-D poststack phase-shift migration in transversely isotropic media[J]. Seg Technical Program Expanded Abstracts, 1994, 13(1): 1205.

[10] Sena A G. Seismic travel-time equations for azimuthally anisotropic and isotropic media: Estimation of interval elastic properties [J]. Geophysics, 1991, 56(12): 2090-2101.

[11] 孙长赞,王建民,冯文霞,等.各向异性叠前时间偏移技术在大庆探区的应用[J].石油地球物理勘探, 2010, 45(增刊1): 71-73.

Sun C Z, Wang J M, Feng W X,etal. Application of anisotropic pre-stack time migration technique in Daqing exploration area[J]. Oil Geophysical Prospecting, 2010, 45(S1): 71-73. (in Chinese)

[12] 万志超,滕吉文,张秉明.各向异性介质中地震波速度分析的研究现状[J].地球物理学进展,1997,12(3):35-44.

Wan Z C, Teng J W, Zhang B M. Velocity analysis in anisotropic media[J]. Chinese Journal of Geophysics, 1997, 12(3): 35-44. (in Chinese)

[13] Johnston J E, Christensen N I. Seismic anisotropy of shales[J]. Journal of Geophysical Research Solid Earth, 1995, 100(B4): 5991-6003.

[14] Sayers C M. The effect of kerogen on the elastic anisotropy of organic-rich shales[J]. Geophysics, 2013, 78: 65-74.

[15] Alkhalifah T. Velocity analysis for transversely isotropic media[J]. Seg Technical Program Expanded Abstracts, 1994, 60(9): 1550-1566.

[16] 王康宁,孙赞东,侯昕晔.基于傅里叶级数展开的纵波方位各向异性裂缝预测[J].石油物探,2015,54(6):755-761.

Wang K N, Sun Z D, Hou X Y. Fracture prediction by P-wave azimuthal anisotropy based on Fourier series decomposition[J]. Geophysical Prospecting For Petroleum, 2015, 54(6): 755-761. (in Chinese)

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

[18] Alkhalifah T, Tsvankin I. Velocity analysis for transversely isotropic media[J]. Geophysics, 1995, 60: 1550-1566.

[19] 兰洁.西南油气田四川长宁区块页岩气产量超6000×104m3[J].天然气与石油,2014,32(5):40.

Lan J. Southwest oil and gas fields Sichuan Chang-ning block shale gas production over 6000×104m3[J]. Natural Gas and Oil, 2014, 32(5): 40. (in Chinese)

PredictionmethodforshalereservoiranisotropybasedontheTsvankinparameter

LI Xiang1, YIN Chen2, WU Yonghong3, XIONG Xiaojun1, LIU Yang1, JIAN Shikai1

1.StateKeyLaboratoryofOilandGasReservoirGeologyandExploitation,ChengduUniversityofTechnology,Chengdu610059,China;2.GeophysicalExplorationCompany,ChuanqingDrillingEngineeringCo.Ltd.,CNPC,Chengdu610213,China;3.ChongqingGasField,SouthwestOilamp;GasFieldCompary,CNPC,Chongqing420021,China

The anisotropy of rock is related to seismic velocity, brittleness, rock elastic parameters and so on. Quantitative description of anisotropy for shale is very important in the early exploration and post-development process of shale gas. VTI media model is used to simulate the shale depending on the constraint with pre-stack migration results and iteratively calculate the anisotropic strength Tsvankin parameter, and the Tsvankin parameter is used to correct the remaining time difference in pre-stack gather. Finally, the output volume of the Tsvankin parameter data is obtained. The anisotropy parameter from Longmaxi Formation shale is used to discuss the significance of shale anisotropy of reservoir exploration and to describe the shale brittleness. It reveals that the Tsvankin parameter of the Longmaxi Formation shale has positive correlation with probability distribution of shale brittleness. It also can be used to predict the intensity of shale anisotropy and reflect the strong and weak of formation brittleness.

Tsvankin parameter; shale; anisotropy; brittleness

P631.4; TE132

A

10.3969/j.issn.1671-9727.2017.06.11

1671-9727(2017)06-0738-06

2016-11-30。

国家科技重大专项(2016ZX05023004-001-003); 国家自然科学基金项目(41274130)。

李翔(1993-)，男，硕士研究生，研究方向：地球物理勘探, E-mail：lixiang_cdut@163.com。

熊晓军(1980-)，男，博士，教授，研究方向：地球物理勘探, E-mail:xiongxiaojun07@cdut.cn。