缺失声波条件下的页岩储层地应力测井解释方法

时 贤 程远方 常 鑫 蒋 恕 王 欣

1.中国石油大学（华东）石油工程学院 2.美国犹他大学能源与地球科学研究院 3.中国石油勘探开发研究院廊坊分院

北美页岩气商业化开发的巨大突破受益于水平分段压裂等技术的快速进步，衍生出旨在创造复杂缝网的压裂增产理论［1－3］。地应力是分析水力裂缝扩展规律，进行射孔方案及压裂决策与设计的核心参数［4－7］。对页岩储层而言，因层理面发育等地质特征，造成页岩层表现出较强的非均质性和各向异性，传统地应力模型存在较多不适应性。国内外学者Hornby、Schatz、Cicotti、Xu、Barree等通过实验和理论论证了各向异性引 起 的 地 层 应 力 改 变［8－14］。Jaeger、Deenadayalu、Khan、George、Higgins、Olson等则分析了页岩各向异性力学参数在压裂设计中的作用，并将储层力学参数和地应力各向异性看作是复杂裂缝网络形成的原因之一［15－21］。利用声波测井数据是计算地应力的常用方法，较单点深度的室内岩石力学实验而言更具应用价值，但因声波测井费用昂贵，实际工程中进行双极声波测井的井数还不到总施工井数的1%，所以地应力计算数据资料十分有限［22］。

以上客观条件为页岩储层的地应力解释带来巨大不便，所以需要在缺失声波时差条件下，形成一套现场实用的页岩层各向异性地应力解释方法。基于权重分析思想，可参照临井声波数据，并通过中子、密度、伽马等测井资料来反求未进行双极测井油气井的相关声波数据，表现出较好应用价值［22－23］。笔者利用该方法，实现了对缺失声波数据的回归，同时基于该数据分析了某区块页岩储层的各向异性，并分别通过H法（Higgnis）和B法（Blanton）计算了考虑各向异性的地应力大小。笔者旨在系统介绍一种缺失声波测井数据条件下如何进行各向异性地应力计算的实用方法，更好的服务页岩气增产开发工作。

1 基于综合测井资料的声波回归解释方法

通过声波测井解释地应力避免了取心对地层资料原始性质破坏，连续性好、可靠性较高、被广泛用于现场。通过现场应用发现，虽然一个区域内仅在几口井进行了声波测井，但这些声波测井数据和临井其他测井资料，如电阻率、伽马、中子孔隙度和人工神经网络回归的孔隙度等其他数据存在一定关系，所以，可以以此作为基础，实现对缺失声波资料的油气井声波数据的反演解释［22］。以纵波数据为例，如果存在传导率、中子孔隙度和密度孔隙度等测井数据，则需要首先绘制传导率、中子孔隙度和密度孔隙度等倒数与已有临井纵波时差的交汇图，再线性拟合出回归方程，然后将各数据带入回归方程中，即可获得初始声波数据。为了进一步提高精度，需要对初始声波时差与临井声波时差取平均值，再根据平均加权后获得权重因子，从而建立起精度更好的邻井纵波时差的解释模型。同样，横波时差数据也可通过相同计算步骤获取。该方法的优势在于在有限声波测井资料的条件下，可以拟合出高精度的声波数据。并且，提出的权重因子可使回归声波速度结果避免井眼冲蚀和气层存在的影响。下面以纵波时差为例介绍该模型的计算公式，在已知中子孔隙度，电阻率及人工神经网络回归孔隙度等测井资料的条件下，参照邻井实测声波时差（dvp）计算缺失井纵波时差（dvpcom）的表达式为：

2 页岩各向异性成因及本构特征

岩石具有各向异性的原因被认为是地层沉积、压实、胶结等复杂物理和化学过程综合作用的结果，其中层理发育是主因。除了岩石层理结构的因素外，天然裂缝、构造运动、应力加载历史和应力集中效应、孔隙压力改变等也都有可能导致岩石各向异性的产生。

对于页岩层而言，其各向异性通常表现在垂直和水平两个平面上的差异。所以，可以通过横观各向同性模型（Transverse Isotropic）对此进行描述（图1）。

图1 横观各向同性模型中各向异性的单元体图

该模型认为页岩在水平方向是各向同性的，并以垂向和横向上不同的弹性模量Eh和Ev对垂向的各向异性特征进行定量描述［10，24］。根据线性材料应变与应力本构关系，需要四阶张量，共81个独立分量组成的矩阵进行表达。横观各向同性模型从胡克定律出发，针对各向同性特点，应力应变关系仅用5个独立的弹性常数进行描述，模型的本构形式如（3）所示。即

对于横向各向同性材料而言，Cij可以表达为：

在该模型当中，仅用5个独立变量C11，C12，C13，C55，C66就可以表达弹性模量（E）、泊松比（v）及剪切模量（G）。Ev与Eh为水平和垂直杨氏模量，其表达公式为式（5）和（6）［25］：

vh与vv为水平和垂直方向的泊松比，可由下式进行计算：

其中，C11、C12和C66存在线性关系，但是C33、C44、C66需要通过纵波，横波及管流声波速度分别进行计算，通常需要高级的声波解释工具来进行测量［21］。

通过上述测井声波解释方法，就实现了对各向异性岩石力学参数的测量。另外，实验室也提供了计算各向异性弹性模量等静力学参数的测量方法，但通常式样大小有限，费用昂贵，所以常作为核心目的层测井解释数据的校正。实验测试方法如图2所示。它需要在垂直，水平和角度为45°这3个方向进行分别测试，才可取得各向异性的力学参数结果［19］。

3 页岩各向异性地应力解释模型

横观各向同性地应力模型和常规地应力模型的差异主要在水平地应力的计算上，目前有两种计算模式：①通过提出各向异性参数（K0）的Higgins法来对经典地应力模型的进行修正，简称H法；②考虑岩石岩性差异引起的应变不同，并耦合温度等其他参数的Blanton法，简称 B法［25］。

图2 考虑各向异性的岩石力学实验室取心方法图

3.1 Higgins地应力解释方法

各向异性参数（K0）代表横向和纵向的弹性模量及泊松比对岩石在平行和垂直层理面的两个方向上的对岩石力学性质的不同影响程度，即

所以，根据各向异性参数（K0）建立的最小地应力解释方程为：

根据式（10）可以看出，对于各向异性岩石，水平弹性模量增加会导致Eh／Ev项的增加，从而使最终计算水平最小主应力结果也增加。对于各向同性岩石，Eh／Ev值为1，则与常规地应力模型相同。

Higgins等基于应变理论提出一套新的地应力解释方法，它主要将各向异性对构造应力的影响更加准确地考虑到方程之中，其计算公式为［9］：

目前，没有理想的方法确定εH与εh的标定关系，一般在计算中多假设应变εH为水平最小主应力方向上的应变εh的两倍。

3.2 Blanton地应力解释方法

在进行最小水平主应力计算时，很多时候并没有考虑构造应力，岩石线弹性和热变量的影响。因此，Blanton等提出一种综合系数校正方法来计算的最小水平主应力［8］：

其中构造应变的大小需要通过对目的层进行相关测试来获取，热扩展系数（εT）在岩石拉伸时取负值，而压缩时取正值。通常，αT的取值和岩性密切相关，对于砂岩来说，αT一般为5.56×10－6／F，页岩为5.00×10－6／F，而碳酸盐大约为4.44×10－6／F。对于ΔT而言，则可以通过地区地热梯度和井底温度进行反推。

图3 MB井横、纵声波时差与中子孔隙度和密度孔隙度的交汇图

图4 MC井横、纵声波时差与中子孔隙度和密度孔隙度的交汇图

4 实例分析及模型应用效果分析

以美国某页岩区块为例介绍在缺失声波资料的情况下如何回归岩石力学参数和计算地应力的基本过程，并分析该地区的各向异性程度。该页岩层位于3 200～3 900m的白垩系 Mancos地层，属于海相沉积页岩，油气藏上部由粉砂岩，泥岩及页岩夹层组成，中间有大段较纯的黑色页岩层，其中页岩有机质含量介于3%～10%，孔隙度介于4%～8%，渗透率介于（5～20）×10－6mD，属于极低渗透的页岩气藏。所选区域目前共有3口页岩气井 MA、MB、MC，其中仅有1口气井MA具有裸眼的双极声波测井数据，井MB和MC同时具有中子孔隙度和密度孔隙度测井资料，MC另外具有自然伽马资料，MB具有部分电阻率测井资料。由于只有中子孔隙度和密度孔隙度测井资料。所以，首先将MB和MC的中子孔隙度，密度孔隙度测井资料与井MA的纵波和横波数据分别进行回归解释。图3（a～d）分别表示MB井横波时差和纵波时差分别和密度孔隙度、中子孔隙度的交汇图，图4（a～d）分别表示MC井横波时差和纵波时差分别和密度孔隙度、中子孔隙度的交汇图。

在获取上述交汇图后，可以回归出相应的声波时差。根据加权平均的思想，通过计算求得MB和MC井纵横波中子孔隙度和密度孔隙度权重因子（表1）。

表1 MB井及MC井权重因子统计表

图5 MB井与MC井垂直与水平弹性模量与泊松比对比图

表2 定井深条件下各向异性系数与岩性解释结果对比表

根据已知A井声波时差资料，可以回归出MB井及MC井的声波数据，然后结合测井解释的其他各向异性参数，可以对动态杨氏模量和泊松比进行分别计算（图5）。表2是通过QEMSCAN的岩性解释结果与岩石力学参数的对比。

运用美国犹他大学的QEMSCAN技术对取心处岩样的岩石矿物组分进行了定量分析，显示了矿物的具体分布结构，实验结果发现所取岩样的胶结物成分主要为伊利石、高岭石、石英以及伊／蒙混层等，图6－a中显示的黄色区域为黄铁矿层，并清晰显示出页岩的层理，图6－b中的层理较不发育，与定井深点的各向异性参数比较如表2所示。通过岩石力学参数的解释结果比较发现，石英含量较高，脆性指数较高的页岩层具有较低的泊松比，而黏土含量高的页岩泊松比较高，主要由其轴向和横向的变形差异较小导致。通过回归的垂直与水平杨氏模量和泊松比的曲线发现，整体上看，水平弹性模量要略高于垂直弹性模量，垂直弹性模量更符合垂向的岩性变化，并且页岩层的各向异性相对较为突出。对于泊松比而言，垂向泊松比值比横向泊松比略大，但整体差异较小。

根据解释的力学参数结果，分别应用传统方法和两种考虑页岩各向异性的方法对地应力进行计算分析，H法主要根据公式（11），而B法则根据公式（14）对最小水平主应力进行计算，计算结果如图7所示。构造应力的计算需要通过小型压裂获取的闭合应力结果进行校正，同时构造应力的大小也和岩性相关，需要根据层位进行分别标定。所以，对于B法中的应变值，一般随岩性不同而变化，对于砂岩来说，一般取0.000 27，而页岩一般取0.000 24。

图6 高精度（微米级精度）QEMSCAN扫描结果图

图7 考虑各向异性的MB井和MC井地应力剖面与常规方法地应力剖面关系图

表3 小型压裂数据与C井回归地应力数据对比表

将所求地层最小水平主应力和小型压裂闭合应力数据进行比较，并做误差分析，结果如表3所示。通过应力结果比较发现，H法和B法预测的结果都比常规方法的地应力值要高，其中B法的值最大，主要在于模型中泊松比各向系数最小。H法和B法计算地应力值较高的原因认为是模型中考虑了不同岩性导致的构造应力差异。通常在构造板块活跃的区域，构造应力的作用忽略会导致计算的最小水平主应力的增加。对于常规方法计算的应力剖面而言，可以看出，如果考虑构造应力项，只能在最终应力剖面上增加一条平滑的曲线，忽略了地层的差异性。考虑各向异性的方法由于将岩石力学性质改变导致的应力改变施加到不同的井深当中，所以结果将更加精确。从压裂施工的角度来看，考虑各向异性计算出最小水平主应力值偏高。因此将更加有利射孔层位的选择及对水力裂缝缝高延伸的约束。

5 结论

1）通过综合测井资料和临井声波资料，可以实现缺失声波数据的反演，大大扩展了声波数据的利用广度。从技术操作上看，该种方法应用简单且精确度较高，可为力学参数与地应力的计算提供技术参数。

2）进行考虑各向异性页岩的力学参数解释时，需要分别计算垂向和横向杨氏模量值和泊松比值。一般来说，水平杨氏模量值比垂直弹性模量值高，并且垂直杨氏模量值更符合岩性规律，而垂向与横向泊松比值相差不多，但垂向泊松比值要比横向泊松比值略大。

3）由于传统地应力模型在构造应力的处理上仅仅做了简单的处理，所以附加在整个深度总应力剖面上的构造应力为平滑直线，无法体现出岩石力学性质对应力的影响作用，同时也不利于区分各层的应力差异。横观各向同性地应力模型考虑了岩石力学性质改变对每个深度上地应力的影响，其中泊松比值越小则最小主应力越大，其中B法预测的最小地应力最大正是因为泊松比系数较低所导致，H法主要受横向泊松比影响，所以计算的最小水平主地应力较大。使用横观各向同性模型将更有助于准确计算设计的裂缝几何形态，但对于各向同性地层，横观各向同性模型和常规模型的结果则基本相同。另外，通过成像测井等技术可对所求地应力结果做进一步验证。

致谢：本研究得到了美国犹他大学能源与地学研究院（EGI）提供的软件及资料支持，另外对中国留学基金委（CSC）的海外留学资助也表示衷心感谢。

符 号 说 明

Ed为动态杨氏模量，MPa；vd为泊松比；vp为纵波速度，m／μs；vs为横波速度，m／μs；ρ 为密度，g／cm3；Z 为地层深度，m；ρ（z）为上覆岩层密度函数，kg／m3；g 为重力加速度，m／s2；σv为上覆岩层压力，MPa；pp为孔隙压力，MPa；α为Biot系数，无量纲；αT为构造应力，MPa；vs为岩石静态泊松比，无量纲；dvphin为中子孔隙度回归声波时差，μs／m；dvp为实测声波时差，μs／m；W1～W4为纵波声波时差权重因子，无量纲；dvpcom为中子孔隙度回归声波时差，μs／m；dvs为实测声波时差，μs／m；M1～M4为横波声波时差权重因子，无量纲；σij、εkl为二阶应力与应变张量；Cijkl为四阶刚度张量；ξ为线弹性常数；εH与εh为最大水平主应力和最小水平主应力方向的应变；C1为平面应变模量；C2为岩石杨氏模量。

［1］贾承造，郑民，张永峰.中国非常规油气资源与勘探开发前景［J］.石油勘探与开发，2012，39（2）：129－136.JIA Chengzao，ZHENG Min，ZHANG Yongfeng.Unconventional hydrocarbon resources in China and the prospect of exploration and development［J］.Petroleum Exploration and Development，2012，39（2）：129－136.

［2］林腊梅，张金川，唐玄，等.中国陆相页岩气的形成条件［J］.天然气工业，2013，33（1）：35－39.LIN Lamei，ZHANG Jinchuan，TANG Xuan，et al.Conditions of continental shale gas accumulation in China［J］.Natural Gas Industry，2013，33（1）：35－39.

［3］程远方，李友志，时贤，等.页岩气体积压裂缝网模型分析及应用［J］.天然气工业，2014，34（9）：53－58.CHENG Yuanfang，LI Youzhi，SHI Xian，et al.Analysis and application of fracture network models of volume fracturing in shale gas reservoirs［J］.Natural Gas Industry，2014，34（9）：53－58.

［4］胡永全，贾锁刚，赵金洲，等.缝网压裂控制条件研究［J］.西南石油大学学报：自然科学版，2013，35（4）：126－132.HU Yongquan，JIA Suogang，ZHAO Jinzhou，et al.Study on controlling conditions in network hydraulic fracturing［J］.Journal of Southwest Petroleum University：Science ＆Technology Edition，2013，35（4）：126－132.

［5］李勇明，彭瑶，王中泽.页岩气压裂增产机理与施工技术分析［J］.西南石油大学学报：自然科学版，2013，35（2）：90－96.LI Yongming，PENG Yao，WANG Zhongze.Analysis of shale gas fracture stimulation mechanism and operating techniques［J］.Journal of Southwest Petroleum University：Science ＆ Technology Edition，2013，35（2）：90－96.

［6］赵金洲，任岚，胡永全.页岩储层压裂缝成网延伸的受控因素分析［J］.西南石油大学学报：自然科学版，2013，35（1）：1－9.ZHAO Jinzhou，REN Lan，HU Yongquan.Controlling factors of hydraulic fractures extending into network in shale formations［J］.Journal of Southwest Petroleum University：Science ＆ Technology Edition，2013，35（1）：1－9.

［7］王香增，吴金桥，张军涛.陆相页岩气层的CO2压裂技术应用探讨［J］.天然气工业，2014，34（1）：64－67.WANG Xiangzeng，WU Jinqiao，ZHANG Juntao.Application of CO2fracturing technology for terrestrial shale gas reservoirs［J］.Natural Gas Industry，2014，34（1）：64－67.

［8］BLANTON T L，OLSON J E.Stress magnitudes from logs：Effects of tectonic strains and temperature［J］.SPE Reservoir Evaluation ＆ Engineering，1999，2（1）：62－68.

［9］HIGGINS S M，GOODWIN S A，BRATTON T R，et al.Anisotropic stress models improve completion design in the Baxter Shale［C］∥paper 115736－MS presented at the SPE Annual Technical Conference and Exhibition，21－24September 2008，Denver，Colorado，USA.New York：SPE，2008.

［10］SCHOENBERG M，MUIR F，SAYERS C.Introducing ANNIE：A simple and three－parameter anisotropic velocity model for shales［J］.Journal of Seismic Exploration，1996，5（1）：35－50.

［11］VERNIK L，LIU X.Velocity anisotropy in shales：A petrophysical study［J］.Exploration Geophysicists，1997，62（2）：521－532.

［12］BULLER D，HUGHES S，MARKET J.Petrophysical evaluation for enhancing hydraulic stimulation in horizontal shale gas wells［C］∥paper 132990－MS presented at the SPE Annual Technical Conference and Exhibition，19－22 September 2010，Florence，Italy.New York：SPE，2010.

［13］FRANQUET J A，RODRIGUEZ E F.Orthotropic horizontal stress characterization from logging and core derived acoustic anisotropies［C］∥paper ARMA－2012－644 presented at the 46thU.S.Rock Mechanics／Geomechanics Symposium，24－27June 2012，Chicago，Illinois，USA.New York：SPE，2012.

［14］RAJAPPA B，CONLET J T.In－situ sress profiling in the Agency Draw Field，Utah［C］∥paper 107818－MS presented at the Rocky Mountain Oil ＆ Gas Technology Symposium，16－18April 2007，Colorado，USA.New York：SPE 2007.

［15］SONG L，HARELAND G.Minimum horizontal stress profile from logging data for Montney Formation of North East British Columbia［C］∥paper 162233－MS presented at the SPE Canadian Unconventional Resources Conference，30October －1November 2012，Calgary，Alberta，Canada.New York：SPE，2012.

［16］CIPPLLA C L，LIU Diabin，KYTE D G.Practical application of in－situ stress profiles［C］∥paper 28607－MS presented at the SPE Annual Technical Conference and Exhibition，25－28September 1994，New Orleans，Louisiana，USA.New York：SPE，1994.

［17］KHAN S，ANSARI S A，HAN Hongxue，et al.Importance of shale anisotropy in estimating in－situ stresses and wellbore stability analysis in Horn River Basin［C］∥paper 149433－MS presented at the Canadian Unconventional Resources Conference，15－17November 2001，Alberta，Canada.New York：SPE，2011.

［18］WATERS G A，LEWIS R E，BENTLEY D.The effect of mechanical properties anisotropy in the generation of hydraulic Fractures in organic shales［C］∥paper 146776－MS presented at the SPE Annual Technical Conference and Exhibition，30October －2November 2011，Denver，Colorado，USA.New York：SPE，2011.

［19］BARREE R D，GILBERT J V，CONWAYA M W.Stress and rock property profiling for unconventional reservoir stimulation［C］∥paper 118703－MS presented at the SPE Hydraulic Fracturing Technology Conference，19－21January 2009，the Woodlands，Texas，USA.New York：SPE，2009.

［20］LARKIN S D，BROWN E K，BAZAN L W.Stimulation design and post fracture production analysis：A tight gas sand case history［C］∥paper 74361－MS presented at the SPE International Petroleum Conference and Exhibition in Mexico，10－12February 2002，Villahermosa，Mexico.New York：SPE，2002.

［21］WALSH J，SINHA B，DONALD A.Formation anisotropy parameters using borehole sonic data［C］∥paper 2006－TT presented at the SPWLA 47thAnnual Logging Symposium，4－7June 2006，Veracruz，Mexico.

［22］MULLEN M J，ROUNDTREE R，TURK G A.A composite determination of mechanical rock properties for stimulation design（what to do when you don＇t have a sonic log）［C］∥paper 108139－MS presented at the Rocky Mountain Oil ＆ Gas Technology Symposium，16－18April 2007，Denver，Colorado，USA.New York：SPE，2007.

［23］RICKMA R，MULLEN M J.A practical use of shale petrophysics for stimulation design optimization：All shale plays are not clones of the Barnett Shale［C］∥paper 115258－MS presented at the SPE Annual Technical Conference and Exhibition，21－24September 2008，Denver，Colorado，USA.New York：SPE，2008.

［24］AMADEI B.Importance of anisotropy when estimating and measuring in situ stresses in rock.［J］.International Journal of Rock Mechanics and Mining Sciences ＆ Geomechanics，1996，33（3）：293－325.

［25］AMADEI B，SAVAGE W Z，SWOLFS H S.Gravitational stresses in anisotropic rock masses［J］.International Journal of Rock Mechanics and Mining Sciences ＆ Geomechanics，1987，24（1）：5－14.