利用地应力场模拟确定压裂诱发油气套管变形位置的方法与流程

1.本发明涉及页岩气开发技术领域，具体来讲，涉及一种利用地应力场模拟确定压裂诱发油气套管变形位置的方法。


背景技术：

2.诸如我国四川盆地及其类似结构的高应力地区的页岩气地质构造复杂，储层最大水平主应力接近甚至高于上覆岩石压力，受多种因素影响及多因素叠加，发生套管变形的概率相对较大，压裂过程是最直观的体现。近年来，虽然我国采取并尝试了一些降低套管变形的技术措施，取得了一定成效，但套管变形的概率仍然较高，制约了页岩气的快速开发，同时也增大了页岩气开发成本。
3.现有技术中，国内外虽然开展了相关套管变形的研究，但大多是分析不同区块的套管变形在时间和空间上的分布情况，总结其变化规律，并提出可能的潜在影响因素，随后开展页岩实验和理论分析，明确压裂过程中的套变变形机理，最后提出防治措施。例如，于2019年5月1日公开的名称为页岩气水平井井筒完整性研究、作者为刘奎的博士论文文献中记载了一种根据页岩气井筒附近应力分布规律，建立局部载荷条件下套管应力状态与套管变形的计算模型。该论文文献研究并讨论了局部载荷严重降低了套管的抗挤毁能力，且套管的屈服变形受局部载荷形式、载荷范围、载荷大小、套管外径、套管壁厚等多种因素的影响，径向局部载荷比平行局部载荷更易造成套管屈服损坏，而增加套管壁厚、减小套管外径有利于减小页岩气套管损坏。但该论文文献只研究了套管变形的机理和降低套管变形的理论措施，且作者提出的增加套管壁厚、减小套管外径来减小页岩气套管损坏的防治措施并不适用于改善目前已施工运行的页岩气井段。于2019年5月1日公开的名称为页岩气压裂井套管载荷分析及变形机理研究、作者为郭雪利的博士论文文献中记载了一种基于环空流体物性随温度动态变化规律和页岩实验结果，利用阶段有限元建模方法，建立压裂过程中瞬态温度-压力耦合下井筒组合体数值模型。该论文文献开展了各向异性、水泥环窜槽流体产生压降、压裂施工参数、页岩非均质性及刚度退化等因素下套管载荷动态变化规律研究，得出了页岩气井在大排量多级压裂过程中，井下情况异常复杂，压裂施工参数和页岩性质变化等都会在很大程度上增加套管应力的结论。但该论文文献只结合现场压裂施工数据，开展了页岩气压裂井套管载荷及变形机理研究，没有提出具体防止套管变形的措施，无法为现场施工提供参考。
4.这些研究内容虽然对于套管变形的机理和预防方法提供了理论支撑，能够为预判套管变形风险提供参考，但并不能形成直接确定压裂诱发油气套管变形位置的预测方法，也因此无法形成有效的缓解措施，以使施工现场的套管变形几率大幅降低。


技术实现要素：

5.本发明的目的在于解决现有技术存在的上述不足中的至少一项。例如，本发明的目的之一在于提供一种能够结合套管变形的数值模拟和测井数据直接确定压裂诱发油气
套管变形位置的方法。
6.为了实现上述目的，本发明提供了一种利用地应力场模拟确定压裂诱发油气套管变形位置的方法，所述方法包括以下步骤：按照“井-震结合”的技术方法，在地震波数据的基础上，结合单井层位信息，建立目标区块的三维地质模型；根据目标区块内单井的测井数据计算岩石力学参数，针对其他测井数据和岩石力学参数进行单井地质力学分析，获取单井地质分析结果；综合分析世界应力地图的信息、已有井的单井测量信息、以及目标区块内水平井压裂微地震监测信息，获取目标区块内的最大水平主应力方向角；引入三维地质模型，建立目标区块的地应力场有限元模型，并基于单井地质分析结果和已有井在压裂施工期间的微地震实测结果，验证目标区块地应力场的模拟结果，将验证合格的地应力场数值解构建为区块精细地应力场；基于区块精细地应力场，建立地质-压裂工程-水泥环-套管一体化的有限元模型，对水平井全水平段的套管柱进行三维有限元变形及应力分析，得到套管柱在压裂之前的初始应力分布；对通过单井的目标层段压裂施工后的套管柱进行三维有限元变形及应力分析，计算压裂引起的地应力场变化，找出对套管施加剪切载荷的位置分布；将模拟结果中对套管施加最大剪切载荷的位置预测为压裂诱发油气套管变形的位置。
7.在本发明的一个示例性实施例中，所述在地震波数据的基础上，结合单井层位信息，建立目标区块的三维地质模型，可包括：根据三维地震波数据建立地质模型的网格几何尺寸，划分地层并定义各地层的网格单元大小，其中，储层所在地层的网格单元小于其他地层的网格单元。
8.在本发明的一个示例性实施例中，所述储层所在地层的网格单元厚度与所述其他地层的网格单元厚度之比可为1：18～1：30。
9.在本发明的一个示例性实施例中，所述其他测井数据可包括地层的伽马射线、压缩声波时长和密度，所述岩石力学参数可包括杨氏模量、泊松比、内聚力和内摩擦角。
10.在本发明的一个示例性实施例中，所述单井地质分析结果可包括利用其他测井数据间接分析获得的地应力结果和利用岩石力学参数直接计算获得的地应力结果。
11.在本发明的一个示例性实施例中，所述综合分析世界应力地图的信息、已有井的单井测量信息、以及目标区块内水平井压裂微地震监测信息，获取目标区块内的最大水平主应力方向角，可包括：利用世界应力地图判断目标区块所属区域的最大水平主应力方向角，明确目标区块内最大水平主应力方向角所属的区间范围；根据单井测量的信息，判断目标区块内不同构造部位的最大水平主应力方向角所属的区间范围，明确目标区块内最大水平主应力方向角随地势的变化规律；分析微地震监测信息，并结合微地震监测信息的分析结果验证和修正利用世界应力地图和单井测量信息总结的目标区块内的最大水平主应力方向角的取值，和/或补充修正利用世界应力地图和单井测量信息所遗漏区域的目标区块内的最大水平主应力方向角的取值。
12.在本发明的一个示例性实施例中，所述分析微地震监测信息可包括：将目标区块内水平井的微地震事件点中出现呈单一色彩条状分布的位置判定为同一个压裂段产生的缝网条带，并将缝网条带的方向角判断为该位置的最大水平主应力的方向角；将目标区块内水平井的微地震事件点中出现由多种色彩组成、且事件点分布超出储层范围的彩色条带的位置判定为天然裂缝的位置；将目标区块内水平井的微地震事件点中出现微地震响应点成片状或团状分布的情况判断为该位置的最大水平主应力方向性不明显，最小水平主应力
方向和垂向主应力方向接近。
13.在本发明的一个示例性实施例中，所述建立目标区块的地应力场有限元模型，可包括：将地质模型网格采用的单元设置为三维8节点线性单元，地应力场有限元模型的载荷设为重力载荷，地应力场有限元模型的四边和底部的边界条件设置为法向位移约束，顶部的边界条件设置为自由边界，并将初始地应力参数和初始孔隙压力参数设置为初始条件。
14.在本发明的一个示例性实施例中，所述地应力场有限元模型中储层所在地层采用的单元可设置为c3d8rp位移-孔隙压力耦合单元，储层以外的其他地层采用的单元可设置为c3d8r位移单元。
15.在本发明的一个示例性实施例中，所述初始地应力参数可包括单井地质力学分析获得的三轴地应力主分量和最大水平主应力方向角。
16.在本发明的一个示例性实施例中，所述基于单井地质分析结果和已有井在压裂施工期间的微地震实测结果，验证目标区块地应力场的模拟结果，将验证合格的地应力场数值解构建为区块精细地应力场，可包括：针对目标区块的每个单井，都分别利用其他测井数据间接分析获得的地应力结果和利用岩石力学参数直接计算获得的地应力结果作为初始地应力参数，进行目标区块地应力场的数值模拟，得到对应的地应力场数值解；对比上述两种地应力场数值解与已有井压裂施工的地应力场实际测量值的误差，并将与地应力场实际测量值误差较小的数值解作为区块精细地应力场。
17.在本发明的一个示例性实施例中，在所述对比上述两种地应力场数值解与已有井压裂施工的地应力场实际测量值的误差时，可选择地应力场数值解的最小主应力结果与已有井压裂施工的地应力场实际测量值最接近的数值解作为后续压裂及套管变形模拟的输入数据。
18.与现有技术相比，本发明的有益效果和优点包括以下内容中的至少一项：
19.(1)本发明通过将单井的测井数据和岩石力学参数分别输入地应力场模型进行模拟计算，并将模拟结果与真正实施压裂时的微地震数据进行反复对比验证，从而获得能够较为真实的再现目标区块的地应力场分布和裂缝走向的精细地应力场；
20.(2)本发明将区块精细地应力场作为后续压裂及套管变形模拟的输入数据，消除了现有技术中由于测井数据误差大或实测参数不足造成的模型输入参数的不确定性，保证了输入地应力场的精度，提高了利用地质-压裂工程-水泥环-套管一体化的有限元模型预测套管变形位置的准确性。
21.(3)本发明的预测方法以测井数据为依据和分析基础，具有良好的可操作性。
附图说明
22.图1示出了本发明的一个示例性实施例中的利用地应力场模拟确定压裂诱发油气套管变形位置的计算流程图。
23.图2示出了本发明的一个示例性实施例中的威202区块的三维地质模型的结构示意图。
24.图3a示出了本发明的一个示例性实施例中的威202直井的测井数据图；图3b示出了本发明的一个示例性实施例中的威202直井的岩石力学参数图。
25.图4示出了本发明的一个示例性实施例中的威202直井的单井地质分析结果图。
26.图5示出了本发明的一个示例性实施例中的威远地区13口单井上的最大水平主应力方向角分布图。
27.图6a示出了本发明的一个示例性实施例中的威202井区的微震信息图；图6b示出了本发明的一个示例性实施例中的威204井区的微震信息图；图6c示出了本发明的一个示例性实施例中的威202h10-3水平井的压裂微震信息图；图6d示出了本发明的一个示例性实施例中的威204h1-5的微地震监测信息。
28.图7示出了本发明的一个示例性实施例中的威202区块的地应力场有限元模型网格图。
29.图8示出了本发明的一个示例性实施例中的威202区块五峰组-龙马溪组地层中的最大水平主压应力方向分布图。
30.图9a示出了本发明的一个示例性实施例中的威202区块五峰组-龙马溪地层中的最小水平主应力方向分布图；图9b示出了威202区块五峰组-龙马溪地层中的中间主应力方向分布图；图9c示出了威202区块五峰组-龙马溪地层中的中间主应力的三维矢量方向分布图。
31.图10示出了本发明的一个示例性实施例中的威202直井储层单元地应力场三维数值解分布云图。
具体实施方式
32.在下文中，将结合示例性实施例和附图来详细说明本发明的利用地应力场模拟确定压裂诱发油气套管变形位置的方法。
33.在本发明的一个示例性实施例中，本发明提供了一种利用地应力场模拟确定压裂诱发油气套管变形位置的方法。
34.在本实施例中，图1为利用地应力场模拟确定压裂诱发油气套管变形位置的计算流程图。如图1所示，利用地应力场模拟确定压裂诱发油气套管变形位置的方法可包括：
35.步骤s1、按照“井-震结合”(即测井与地震结合)的技术方法，在地震波数据的基础上，结合单井层位信息，建立目标区块的三维地质模型。
36.具体地，结合目标区块的构造起伏特征，在三维地震波数据的基础上，可构建目标区块整体的几何形状和网格大小，然后结合目标区块已有井的单井层位信息，可划定目标区块的层位信息，确定储层所在地层和储层之外的其他地层的几何形状和网格大小，获得目标区块的三维地质模型。
37.这里，建立目标区块的三维地质模型的具体实施方式，可包括：根据三维地震波数据建立地质模型的网格几何尺寸，划分地层并定义各地层的网格单元大小，其中，储层所在地层的网格单元小于其他地层的网格单元。所述储层所在地层的网格单元厚度与所述其他地层的网格单元厚度之比可为1：18～1：30。例如，其他地层的网格单元厚度可为储层所在地层的网格单元厚度的20倍。
38.以威远地区页岩气三维区块威202区块的实地应用为例进行说明。
39.威远地区页岩气三维区块内地面地腹构造格局一致，构造简单，西北高东南低，轴线近东西向。三维区块中北部主要存在威ⅰ号高点，南部为威远构造的南翼，断层稀少且断层落差较小。从侏罗系中统沙溪庙组叶肢介页岩构造可知：三维区块内构造简单，断层不发
育，中北部为威远构造的威ⅰ号高点，威ⅰ号高点上存在三个次一级高点，轴向近东西向；三维区南部为威远构造南翼。地腹构造格局与地表大体一致，但局部构造细节有一定的变化，褶皱相对增强，断层相对发育。由于断层落差较小，其对初始地应力场分布的影响不大。另外，由于断层落差较小，通过地震波数据识别断层的难度较大，有些落差很小的断层难以识别出来。
40.由现场数据可知，威202区块的长为19km，宽为13.5km。根据三维地震波数据，建立的威202区块的三维地质模型，并结合单井所在地区的地层层序，将地质模型划分为7个层位，从上至下分别为：地面至嘉陵江组底界、嘉陵江组底界至飞仙关组底界、飞仙关组底界至上二叠纪底界、上二叠纪底界至下二叠纪底界、下二叠纪底界至五峰组底界(包括了龙马溪组地层)、五峰组底界至寒武系底界、寒武系底界至海拔-5500米。其中，储层位于五峰组-龙马溪组地层。
41.考虑到有限元计算量的优化，与其他地层的网格相比，在五峰组-龙马溪组地层(也就是储层所在地层)划分了较密的网格。即五峰组-龙门溪组地层的网格为20层，最小单元厚度1米，而其它地层采用较大的单元尺寸。威202区块的地质模型总共采用的单元数为121068。由于地层剥蚀现象的影响，五峰组-龙马溪组之间没有进行区分，两者一起视作一个地层。
42.图2示出了威202区块的三维地质模型示意图(深度单位，米)，图中给出了储层的几何形状和网格，呈西北高-东南低的构造。图中从上到下不同深度颜色的网格分别代表的是地面至嘉陵江组底界、嘉陵江组底界至飞仙关组底界、飞仙关组底界至上二叠纪底界、上二叠纪底界至下二叠纪底界、下二叠纪底界至五峰组底界(包括了龙马溪组地层)、五峰组底界至寒武系底界、寒武系底界至海拔-5500米。
43.步骤s2、根据目标区块内单井的测井数据计算岩石力学参数，针对其他测井数据和岩石力学参数进行单井地质力学分析，获取单井地质分析结果。
44.所述其他测井数据可包括地层的伽马射线、压缩声波时长和密度。所述岩石力学参数包括杨氏模量、泊松比、内聚力(粘聚力)和内摩擦角。采用经验公式和相关计算原理，根据声波测井数据可计算获得杨氏模量等岩石力学参数。例如，采用下面相关公式(式(1)至式(4))可计算杨氏模量等岩石力学参数。
45.e＝103ρb·
[3(vs/v
p
)
2-4]/vs2[(vs/v
p
)
2-1]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(1)
[0046]
v＝0.5[(vs/v
p
)
2-2]/[(vs/v
p
)
2-1]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(2)
[0047]
c＝4.69433
×
10
7vp4
ρb[(1 v)/(1-v)](1-2v)(1 0.78v
sh
)
ꢀꢀꢀꢀꢀꢀ
式(3)
[0048][0049]
式中，e为杨氏模量，mpa；ρb为岩石密度，g/cm3；vs为纵波，us/m；vp为横波，us/m；v为泊松比，无量纲；c为粘聚力，mpa；v
sh
为泥质含量，％；为内摩擦角，度。
[0050]
所述单井地质分析结果可包括利用其他测井数据间接分析获得的地应力结果和利用岩石力学参数直接计算获得的地应力结果。所述利用其他测井数据间接分析获得的地应力结果可包括利用其他测井数据间接分析获得的最大水平主应力、利用其他测井数据间接分析获得的最小水平主应力、利用其他测井数据间接分析获得的垂向主应力、以及利用其他测井数据间接分析获得的地层孔隙压力。所述利用岩石力学参数直接计算获得的地应力分析结果包括利用岩石力学参数直接计算获得的最大水平主应力、利用岩石力学参数直
接计算获得的最小水平主应力、利用岩石力学参数直接计算获得的垂向主应力、以及利用岩石力学参数直接计算获得的地层孔隙压力。
[0051]
仍以威远地区页岩气三维区块威202区块的实地应用为例进行说明。
[0052]
图3a、图3b和图4为采用步骤s2的方法对威202区块进行单井地质力学分析后获得的单井地质分析结果。其中，图3a和图3b分别为威202区块的某一直井的测井数据和根据声波测井数据计算得到的岩石力学参数。这个井的下部井轨迹最大倾角8.6度。图3a的左一栏为地层的伽马射线(单位：api)，左二栏为压缩声波时长dc(单位：微秒/英尺)，左三栏为井径(单位：英时)，右一栏为密度(单位：g/cc)，右二栏为声波时长(单位：毫秒/英尺)。图3b的左一栏为威202下部井段的弹性模量的取值曲线(单位：mpsi)，左二栏为泊松比的取值曲线，右一栏为地层材料的粘结强度(单位：mpa)，右二栏为内摩擦角(单位：度)，其中粘结强度和内摩擦角是用于根据莫尔-库伦条件计算坍塌压力的主要参数。
[0053]
图4为威202直井的单井地质分析结果，其中，图4的左一栏为威202直井分析的地应力主分量结果的工程梯度形式的显示，图4的右一栏为单井的地应力分析结果的全量显示。图中的曲线a对应的图标为线段shg elastic(w202)，曲线b对应的图标为线段shg-shg-obg(w202)，曲线c对应的图标为线段sfg mohr-coulomb(w202)，曲线d对应的图标为线段pp dt e3(w202)，曲线e对应的图标为线段obg rhob(w202)，点a对应的图标为方块shgmin(w202)，点b对应的图标为方块shmax(w202)，点d对应的图标为方块pp(w202)，点e对应的图标为方块obg-m(w202)。
[0054]
图4中各曲线的图标所代表的含义如下：曲线obg rhob(即曲线e)代表计算得到的(也就是利用岩石力学参数直接计算获得的)垂向应力(即上覆岩层压力)，方块obg-m(即点e)代表测量分析得到的(也就是利用其他测井数据间接分析获得的)垂向应力；曲线shg elastic(即曲线a)代表计算得到的(也就是利用岩石力学参数直接计算获得的)最小水平主应力，方块shgmin(即点a)代表水力压裂方法测得的(也就是利用其他测井数据间接分析获得的)最小水平主应力；曲线shg-shg-obg(即曲线b)代表计算得到的(也就是利用岩石力学参数直接计算获得的)最大水平主应力，方块shmax(即点b)代表测量分析得到的(也就是利用其他测井数据间接分析获得的)最大水平主应力；曲线pp dt e3(即曲线d)代表计算得到的(也就是利用岩石力学参数直接计算获得的)地层孔隙压力，方块pp(即点d)代表测量分析得到的(也就是利用其他测井数据间接分析获得的)地层孔隙压力。曲线sfg mohr-coulomb(即曲线c)代表由摩尔库伦塑性屈服条件计算得到的(也就是利用岩石力学参数直接计算获得的)坍塌压力，即泥浆密度窗口下限。
[0055]
图中可以看出，方块obg-m(也就是利用其他测井数据间接分析获得的垂向应力)的数值明显小于曲线obg(也就是利用岩石力学参数直接计算获得的垂向应力)的数值，相对误差大约为10％。曲线shg在垂深大约2560米左右的深度上(也就是利用岩石力学参数直接计算获得的最大水平主应力)的数值与方块shg-max(也就是利用其他测井数据间接分析获得的最大水平主应力)的数值相重合，说明这个shg曲线结果是合理的。曲线shg在垂深大约2560米左右的深度上(也就是利用岩石力学参数直接计算获得的最小水平主应力)的值与方块shgmin(也就是利用其他测井数据间接分析获得的最小水平主应力)的数值相重合，说明这个shg曲线结果是合理的。曲线pp在垂深大约2560米左右的深度上(也就是利用岩石力学参数直接计算获得的地层孔隙压力)的数值与方块pp的(也就是利用其他测井数据间
接分析获得的地层孔隙压力)数值相重合，说明这个pp曲线结果是合理的。
[0056]
根据图4的单井地质力学结果，可以认为：威202直井附近地层中垂向应力幅值在三个主应力分量中处于中间，因此属于“走滑断层应力格式”。
[0057]
步骤s3、综合分析世界应力地图(world stress map，wsm)的信息、已有井的单井测量信息、以及目标区块内水平井压裂微地震监测信息，获取目标区块内的最大水平主应力方向角。
[0058]
这里，所述综合分析世界应力地图的信息、已有井的单井测量信息、以及目标区块内水平井压裂微地震监测信息的具体实施方式，是指：先利用世界应力地图判断目标区块所属区域的最大水平主应力方向角，明确目标区块内最大水平主应力方向角所属的区间范围；然后根据单井测量的信息，判断目标区块内不同构造部位的最大水平主应力方向角所属的区间范围，明确目标区块内最大水平主应力方向角随地势的变化规律；最后再分析微地震监测信息，并结合微地震监测信息的分析结果验证和修正上述两种信息总结的目标区块内的最大水平主应力方向角的取值，或补充上述两种信息遗漏的某些区域的目标区块内的最大水平主应力方向角的取值。
[0059]
在本实施例中，所述分析微地震监测信息可包括：将目标区块内水平井的微地震事件点中出现呈单一色彩条状分布的位置判定为同一个压裂段产生的缝网条带，并将缝网条带的方向角判断为该位置的最大水平主应力的方向角；将目标区块内水平井的微地震事件点中出现由多种色彩组成、且事件点分布超出储层范围的彩色条带的位置判定为天然裂缝的位置；将目标区块内水平井的微地震事件点中出现微地震响应点成片状或团状分布的情况判断为该位置的最大水平主应力方向性不明显，最小水平主应力方向和垂向主应力方向接近。
[0060]
仍以威远地区页岩气三维区块威202区块的实地应用为例进行说明。
[0061]
在世界应力地图wsm中的信息显示：四川盆地的区域地应力格式特点为以逆冲断层及走滑断层应力模式为主。区域内最大水平主应力方向以东-西向为主，其它方向如东北-西南及西北-东南也有存在。
[0062]
图5为威远地区内13口井的单井井位上储层五峰组-龙马溪组的最大主应力方向分布。这些信息综合了微地震监测、交叉偶极阵列声波测井分析等应力方位分析结果。图5中的左边第一条虚线a是指海拔-1400m的位置，第二条虚线b是指海拔-2400m的位置，第三条虚线c是指海拔-3100m的位置。如图5所示，单井测量的信息显示：在区块内不同构造部位，最大水平主应力呈现出一定规律的变化趋势。整个威远页岩气区块可以分为：位于左上边界至虚线a之间的顶部平缓带、位于虚线a和虚线b之间的中部陡斜带、位于虚线b和虚线c之间的右下部缓斜带、以及虚线c与右下边界的近洼平缓带。其中顶部平缓带的最大水平主应力方向角为130
°
，中部陡斜带的最大水平主应力方向角为95
°
～105
°
(变化范围10
°
)，右下部缓斜带的最大水平主应力方向角为85
°
～95
°
(变化范围10
°
)，近洼平缓带的最大水平主应力方向角为65
°
～95
°
(变化范围30
°
)。从以上信息可以看出：威远区块内储层五峰组-龙马溪组中的最大水平地应力变化比较复杂，不仅随着水平位置在变化(从130
°
到65
°
～90
°
)、在储层内随深度也有30
°
的变化。
[0063]
图6a、图6b、图6c和图6d为威远区块部分已有水平井的微震监测信息。其中，图6a为威202井区的微震信息，图6b为威204井区的微震信息图，图6c为威202h10-3水平井的压
裂微震信息，图6d为威204h1-5的微地震监测信息。
[0064]
由于微震事件点的位置取决于最大主应力方向和天然裂缝分布两个因素，在图6a、图6b、图6c和图6d中的不同色彩代表不同时间/压裂段产生的微震事件。每一个色彩呈单一色彩条状分布的时候，就意味着这是同一个压裂段产生的缝网条带，而这个缝网条带的方向角就是这个位置的最大水平主应力的方向角。图6a、图6b、图6c和图6d中的实线段就是按照这个原则画上去的，它代表此处最大水平主应力的方向。当代表微震事件的彩色条带是由多种色彩组成、且事件点分布超出储层范围，这个时候，这个彩色条带代表的就是天然裂缝的位置，在图6a、图6b、图6c和图6d中用虚线段表示。
[0065]
也就是说，在图6a、图6b、图6c和图6d中，实线段表示同一个压裂段产生的缝网条带处的最大水平主应力的方向，虚线段表示天然裂缝的位置。在标注图6a中的虚线段所指代的天然裂缝位置时，参考了图6c和图6d中的微震事件分布。从图6d中可以看出，微震事件点超出了目的层箱体，属于天然裂缝相关事件。天然裂缝的方向角取决于地质构造运动的方向，经常与现在的地层主应力方向角不一致，而且一般没有明确的解析关系。
[0066]
另外，在图6a、图6b、图6c和图6d中，微地震响应点成团状分布，表明最大水平主应力方向性不明显，两个主应力(最小水平主应力和竖向主应力)方向接近。如图所示，目标区块内水平井压裂微地震监测信息显示：区块内中部的最大主应力方向主要沿东西走向；区块内两侧位置上的最大主应力方向分别向各自的两侧上方偏离；在威202区块左侧，方向角最大可到110
°
。
[0067]
图6a中的虚线圆圈中的微震事件点呈片状分布，表明此处的两个水平主应力大小接近，没有明显占优的主应力方向。这与图5中的单井最大水平主应力方向角信息相符合，在图5的主应力方向角变化范围较大的位置，微震事件呈片状/团状分布，没有明显的水平主方向。
[0068]
步骤s4、引入三维地质模型，建立目标区块的地应力场有限元模型，并基于单井地质分析结果和已有井在压裂施工期间的微地震实测结果，验证目标区块地应力场的模拟结果，将验证合格的地应力场数值解构建为区块精细地应力场。
[0069]
具体地，先将步骤s2中单井地质分析结果获得的三轴地应力主分量、以及步骤s3中分析获得的最大水平主应力方向角作为模型设置参数输入至步骤s1中获得的目标区块的三维地质模型，可建立目标区块的地应力场有限元模型；然后进行目标区块的地应力场数值模拟，将地应力场数值模拟结果与已有井在压裂施工期间的微地震实测结果和单井地质分析结果进行对比验证，并调整模型设置参数，将验证合格后的地应力场有限元模型作为能够真实模拟目标区块地应力场的数值模型，其获得目标区块的三维精细地应力场数值解可作为后续压裂及套管变形模拟的输入数据。
[0070]
例如，所述建立目标区块的地应力场有限元模型的具体实施方式，可包括：地质模型网格采用的单元为三维8节点线性单元。其中，储层所在的地层设为c3d8rp-孔隙压力耦合单元，其他地层设为c3d8r位移单元。其中，c3d8rp和c3d8r是指软件abaqus中的单元类型，c表示为实体单元，3d表示为三维，8是这个单元所具有的节点数目，r是指这个单元是缩减积分单元，p是指三线孔隙压力。模型载荷设为重力载荷。模型四边和底部的边界条件设为法向位移约束，顶部的边界条件设为自由边界。初始条件输入初始地应力参数和初始孔隙压力参数，其中，初始地应力参数包括三轴地应力主分量(即最大水平主应力、最小水平
主应力和垂向主应力)和最大水平主应力方向角。
[0071]
针对目标区块的每个单井，都可分别输入利用其他测井数据间接分析获得的地应力结果和利用岩石力学参数直接计算获得的地应力结果作为初始地应力参数，进行目标区块地应力场的数值模拟，得到对应的地应力场数值解。然后对比上述两种地应力场数值解与已有井压裂施工的地应力场实际测量值的误差，并将与地应力场实际测量值误差较小的数值解构建为区块精细地应力场，并作为后续压裂及套管变形模拟的输入数据。
[0072]
例如，可以对比分析两种模拟所得的储层的最大水平主应力方向分布是否符合已有井在压裂施工期间的微地震实测结果、储层所属的应力格式是否符合单井地质力学分析结果、以及储层的地应力场分布规律是否符合已有井在压裂施工期间的微地震实测结果等，将两种模拟中误差更小的地应力分析结果判断为更适合设置为重现实际压裂施工过程的区块地应力场的模拟参数，其得到的地应力场数值解也更适合作为后续压裂及套管变形模拟的输入数据。考虑到模型数值结果主要目的是为了后续的“压裂引起的套管变形预测”，模型参数准备时主要强调了数值解最小主应力结果与实测值接近，即使用了“优先保证数值解的最小主应力结果与实测值最接近、兼顾最大水平主应力方向角数值解与实测值中间值趋势一致”的原则，来保证模型的合理性。
[0073]
需要说明的是，单井地质力学分析的目的是分析裂缝的延伸走向，由于裂缝始终沿着垂直最小水平主应力的方向延伸，确定地应力的方位即可预测裂缝延伸的方位。不同的井段所处的地质局部构造细节有一定的变化，故适合模拟该直井的地应力分析结果也有所不同。通过将利用其他测井数据间接分析获得的地应力结果和利用岩石力学参数直接计算获得的地应力结果作为初始地应力参数，分别进行目标区块地应力场的数值模拟，能够判断出两组地应力数据的模拟准确性，从而挑选出最适合模拟该直井的地应力场模拟初始输入参数。如此针对寻找目标区块的所有直井，反复寻找最适合设置为重现实际压裂施工过程(即误差最小)的区块地应力场的模拟参数的地应力分析结果，能够提高地应力场模型的整体模拟精度，从而建立精细地应力场，保证后续压裂及套管变形模拟时，模拟的裂缝延伸走向接近真实的裂缝状态。
[0074]
仍以威远地区页岩气三维区块威202区块的实地应用为例进行说明。
[0075]
将步骤s2中的单井地质分析结果的和步骤s3中的最大水平主应力方向角作为初始条件输入至本步骤的地应力场模型。
[0076]
表1为综合上述单井地质分析结果和经验之后的模型弹性力学参数取值列表。其中的弹性模量是随深度变化的。计算模型中采用了abaqus用户子程序来实现其tvd深度(tvd是指垂直深度)依赖特征。
[0077]
表1综合单井分析结果和经验之后的模型弹性力学参数取值列表
[0078][0079]
根据单井的测井数据来计算地层材料参数(也就是岩石力学参数)如杨氏模量和初始地应力相关参数，能够减少由于实测参数不足造成的模型输入参数的不确定性，保证精细地应力场模型初始输入参数的准确性。
[0080]
另外，步骤s2中的根据摩尔库伦塑性屈服条件计算得到的坍塌压力(例如图4中的曲线c)可与本步骤的地应力场模型的模拟结果进行对比修正，并不一定要输入地应力场模型。
[0081]
如图7所示为威202区块的地应力场有限元模型网格，来自于步骤s1所建立的地质模型网格。采用的单元为三维8节点线性单元，在储层为c3d8rp位移-孔隙压力耦合单元，在储层以外为c3d8r位移单元。长边方向取为x-轴方向。模型的载荷为重力载荷。边界条件为四边的法向位移约束和底部的法向位移约束；顶部为地面，自由边界。
[0082]
使用图7中的三维有限元网格模型，设置初始条件，建立目标区块的地应力场有限元模型。初始条件包括初始地应力场和初始孔隙压力场，其中孔隙压力场仅在目的层五峰组-龙马溪组地层中存在，威202区块的孔隙压力系数为1.4g/cc。初始地应力场的设置按照图4给出的单井地应力分量进行各个地层的地应力分量参数设置。
[0083]
经过数值计算得到的威202区块五峰组-龙马溪组地层中最大水平主压应力方向分布如图8所示。从图8看出，数值结果中的最大水平主应力在威202区块左上部的平缓带的主应力方向约为130
°
；在区块右下部的最大水平主应力方向逐渐过渡到接近东西方向的90
°
。这个结果符合图5给出的方向角测量分析结果。
[0084]
图9a和图9b分别给出了威202区块五峰组-龙马溪组地层中最小水平主压应力和中间主应力的方向，图9c为中间主应力的三维矢量方向分布图。从图9a、图9b和图9c可以看出，在威202区块的右下/东南位置，中间主应力呈接近竖直的方向。这表明这些位置的竖向应力是中间主应力分量，应力格式属于“走滑断层应力格式”。这个结果符合图4的单井地质力学结果。
[0085]
图10为威202直井储层单元地应力场三维数值解分布云图(tvd＝2550m)。单元的边长为275米。图中的符号约定遵从固体力学的符号约定，即拉伸应力为正、压应力为负。从图中看出，应力的幅值大小从西北往东南逐渐增大。同时，表2给出了威202直井地应力数值解和测量值的比较。从表2可以看出，数值结果中的两个水平主应力的幅值与实测值很接近，竖向应力与实测值误差相对较大，为6.37％。
[0086]
表2威202直井地应力数值解和测量值的比较
[0087][0088]
分析上述模拟结果，发现主应力大小和方向的数值解与实测值有一定的误差。原因之一是由于计算采用的简化模型在局部构造细节方面与实际情况可能有偏离。此外，实测值给出的主方向角度是一个范围，而主应力值是一个平均意义上的数值。这说明实测值本身也是一个平均值，有一个误差范围。
[0089]
考虑到模型数值结果主要目的是为了后面的“压裂引起的套管变形预测”，模型参数准备时主要强调了数值解最小主应力结果与实测值接近，即使用了“优先保证数值解的最小主应力结果与实测值最接近、兼顾最大水平主应力方向角数值解与实测值中间值趋势一致”的原则，来保证模型的合理性。
[0090]
因此，综合分析考虑后，认为上述地应力场数值结果能够构建为威202区块的精细地应力场，从而被用作后续压裂及套管变形模拟的输入数据。
[0091]
步骤s5、基于区块精细地应力场，建立地质-压裂工程-水泥环-套管一体化的有限元模型，对水平井全水平段的套管柱进行三维有限元变形及应力分析，得到套管柱在压裂之前的初始应力分布。
[0092]
具体地：在步骤s4中的区块的地应力场有限元模型的基础上建立井、水泥环、套管的几何模型，并设置压裂施工的模型数据，获得地质-压裂工程-水泥环-套管一体化的有限元模型，并将区块精细地应力场的数值模拟结果作为压裂及套管变形模拟的输入数据。
[0093]
例如，建立地质-压裂工程-水泥环-套管一体化的有限元模型可包括以下内容：先建立井的几何模型，然后设置模型数据，模型数据包括储层厚度分布、水平段测井数据、地应力场分布、设计的压裂施工参数、射孔参数、地层压力系数、地应力参数、岩石力学参数，最后考虑多种材料(包括地层材料、水泥环材料和套管材料)、多种几何形状(包括地层几何形状、水泥环几何形状和套管几何形状)、结构变形和物理渗流两种物理场、流体压力载荷、地应力载荷、重力载荷、弹性力学模型、以及塑性力学模型，建立井的三维有限元模型。
[0094]
对于预测套管变形的地质-压裂工程-水泥环-套管一体化的数学模型，它的基础理论模型为金属的塑性加载屈服准则，包括式(5)的特雷斯卡(tresca)屈服准则和式(6)的米泽斯(von mises)屈服准则。二者都是基于剪切强度的屈服准则。换句话说，金属材料的塑性变形是剪切塑性变形。
[0095]
τ
max
＝k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(5)
[0096]
(σ
1-σ2)2 (σ
2-σ3)2 (σ
3-σ1)2＝2σ
s2
＝6k2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(6)
[0097]
其中，σ1，σ2，σ3为三个主应力，mpa，σs为屈服应力，mpa，k为材料的剪切屈服强度，mpa。
[0098]
套管的拉伸塑性，宏观看上去是拉伸塑性变形，但微观机理是金属晶体的剪切塑性滑移。
[0099]
最终建立的“地质-压裂工程-水泥环-套管一体化”的模型包括了：
[0100]
1)多种材料、多种几何形状；
[0101]
2)结构变形和渗流两种物理场；
[0102]
3)流体压力载荷、地应力载荷、重力载荷；
[0103]
4)弹性力学本构模型、塑性力学本构模型。
[0104]
这个复杂模型的力学行为可以用下述式子来概括性表示：
[0105]
k(u)
·
u＝f
[0106]
这里的变量定义为：
[0107]
k是系统刚度矩阵，代表了模型的材料属性包括弹性本构和塑性本构、以及几何特征；u是运动矢量，包括了模型中各点的位移和孔隙压力；f是载荷矢量，代表了模型中涉及的各种载荷。
[0108]
步骤s6、对通过单井的目标层段压裂施工后的套管柱进行三维有限元变形及应力分析，计算压裂引起的地应力场变化，找出对套管施加剪切载荷的位置分布。
[0109]
具体地，可按照设计的压裂施工参数对所述地质-压裂工程-水泥环-套管一体化的有限元模型中的套管柱的目标层段设置模型参数，并进行储层压裂的数值模拟，计算对该目标层段压裂施工后引起的地应力场变化，并找到对套管施加剪切载荷的位置分布。
[0110]
步骤s7、将模拟结果中对套管施加最大剪切载荷的位置预测为压裂诱发油气套管变形的位置。
[0111]
在压裂导致的地层破碎区的外边缘，存在地层刚度不对称、局部刚度明显变化的情况，位于该位置的套管发生变形的几率较大。所述压裂导致的地层破碎区的外边缘可通过储层压裂的数值模拟确定。在储层压裂的数值模拟结果中可获得压裂引起的剪切局部化带，所述剪切局部化带可以视为压裂导致的地层破碎区。这里，剪切局部化带是指地层注入孔隙压力和地应力场的共同作用下剪切应变集中的区域，位于该区域的套管受到的地层载荷明显高于其他位置的地层载荷，故位于该区域的套管发生变形的几率也高于其他位置。
[0112]
综上所述，本发明的有益效果和优点包括以下内容中的至少一项：
[0113]
(1)本发明通过将单井的测井数据和岩石力学参数分别输入地应力场模型进行模拟计算，并将模拟结果与真正实施压裂时的微地震数据进行反复对比验证，从而获得能够较为真实的再现目标区块的地应力场分布和裂缝走向的精细地应力场。
[0114]
(2)本发明将区块精细地应力场作为后续压裂及套管变形模拟的输入数据，消除了现有技术中由于测井数据误差大或实测参数不足造成的模型输入参数的不确定性，保证了输入地应力场的精度，提高了利用地质-压裂工程-水泥环-套管一体化的有限元模型预测套管变形位置的准确性。
[0115]
(3)本发明的预测方法以测井数据为依据和分析基础，具有良好的可操作性。
[0116]
尽管上面已经结合示例性实施例及附图描述了本发明，但是本领域普通技术人员应该清楚，在不脱离权利要求的精神和范围的情况下，可以对上述实施例进行各种修改。