低渗透油藏为胜利油田的重要油藏类型。胜利油区低渗透油探明储量10.10亿t，浊积砂低渗透油藏为其重要组成部分。然而，浊积砂低渗透储层厚度、岩相、岩性、物性变化大，储层内部应力场分布规律复杂^[1-3]，严重制约着“仿水平井”等先进开发技术的有效实施和开发井网的合理部署。研究高精度浊积砂低渗透储层的构造应力场模拟方法，对提高低渗透油藏的开发效果和油田开发效益意义重大。

由于三维建模复杂繁琐、三维岩石参数求取难度大等原因，目前构造应力场的模拟计算基本上从以下两方面进行：1) 二维模拟分析^[4-7]；2) 三维模拟分析，即将一定区块的地层视为均一的地质体，以计算该区块内参考井的岩石力学参数代替该区块的岩石力学参数来进行三维模拟^[8-11]。然而对于构造形态复杂、储层非均质性强的浊积砂低渗透储层，用这些简化方式计算的结果往往不能体现出由于浊积砂体的厚度变化和储层内部岩石力学性质差异而引起的局部构造应力场异常，而该异常区域是制约油田低渗透油藏开发效果进一步提高的关键。

本次研究采用PETREL与ANSYS联合建模技术，通过构建高精度的几何模型和三维层速度模型，全面计算研究区浊积砂储层的三维岩石力学参数，以提高浊积砂低渗透储层三维构造应力场模拟精度，为强非均质储层的高精度构造应力场模拟提供参考。

1 研究区概况及技术路线 1.1 研究区概况

樊142块是胜利油田低渗透油藏综合利用示范基地重点实验开发区块，自2010年开始采用“仿水平井”技术开发。仿水平井开发技术需要通过大型压裂改造形成200 m以上的水平长裂缝，对长裂缝的定向性要求高，长裂缝的延伸方向直接关系到整个井网的开发效果。而构造应力场的分布控制着压裂裂缝的延伸方向，因此，构造应力场研究的精细程度直接影响着研究区开发效果。

樊142块位于高青县城东约15 km的唐坊镇境内，面积39 km²，四周被博兴、正理庄、金家、大芦湖油田包围。构造位置位于济阳坳陷东营凹陷博兴洼陷南斜坡中部的金家—正理庄—樊家鼻状构造带的北东端。目的层为沙三下亚段，其中上部发育的一套滑塌浊积砂体是本次构造应力场模拟的主体。该砂体厚度变化大，构造形态复杂；主体部位为呈复合韵律特征的一厚层，边侧部位呈沙泥薄互层特征，渗透率为(0.11~6.01)×10^-3 μm²，平均渗透率为1.50×10^-3 μm²，在平面和垂向上都表现出较强的非均质性(表 1)。

樊142块自2003年起已进入滚动勘探阶段，地震资料完全覆盖，测井资料齐全。已完成岩心分析井8口，并有27口压裂井裂缝监测资料，为前期三维地质模型的准确建立及后期模拟计算结果的验证提供了依据。

1.2 技术路线

本次研究采用PETREL与ANSYS联合建模技术和利用三维层速度模型计算研究区三维岩石力学参数的方法。首先，利用PETREL软件进行相控储层地质建模，将浊积砂体的岩相、岩性、物性等非均质特征完整地包含在地质模型之中，为开展浊积砂体内构造应力场模拟奠定基础；然后，通过三维点云逆向工程将PETREL构造模型转换为ANSYS几何模型，将地质信息完整地转换到ANSYS计算模型中，大大提高ANSYS几何模型的精度；同时，以井点岩石力学参数为基础，根据PETREL三维层速度场模型将岩石力学参数拓展到ANSYS三维力学模型之中，使力学模型具有地质属性，为精细准确刻画地质体内各点的应力状态提供保障；最后，参考压裂井点处的应力状态确定边界条件，开展三维构造应力场数值模拟计算，进一步应用井点处应力数据反复检验并修正模拟结果，以获得符合地质实际的储层构造应力场分布结果。

2 PETREL三维地质建模

准确构造模型是精确模拟构造应力场的前提，三维层速度模型是计算三维岩石力学参数的数据基础，PETREL三维地质建模可以得到精确的构造模型和三维层速度模型。本次研究中，首先，综合利用钻井数据、测井资料、地震资料和区域地质资料在PETREL中建立研究区三维构造模型；然后，在构造模型的基础上，结合单井沉积微相划分以及平面沉积微相数据，建立沉积微相模拟；最后，在沉积微相的控制下，利用测井数据计算层速度，建立研究区的三维层速度模型^[12-14]。

2.1 构造模型

三维构造模型是由层面模型和断层模型构成的。层面模型的建立主要有两种方法，一种是直接利用地震构造解释所得的构造图拟合成面，另一种是根据地层划分对比得到的分层数据插值成面，且两种方法得到的层面都可以用井点的分层数据加以调整。断层模型主要根据单井断点数据以及地震构造解释所得的平面及剖面的断层面建立。模型构建过程中要根据研究区的构造特征及时修改和调整，实现构造模型与测井资料、地震资料及地质规律的拟合。

三维构造模型共包括4个层面和11条断层。由构造模型可知，樊142块为断层复杂化的北倾单斜构造，沙三下亚段埋深为2 460~2 890 m，平均地层厚度达285 m，砂体的构造形态复杂，厚度变化大，单井钻遇砂体厚度5~16 m；受近东西向挤压作用，发育近东西向的11条正断层，断距10~180 m。

2.2 层速度模型

层速度大小与沉积相分布存在着良好的相关性。在区域变差函数分析的基础上，采用相控建模的方法实现层速度的建模，可以使模型与地质规律更好地吻合。

首先，在区域构造、沉积特征研究的基础上，根据樊142块岩心、测井、录井和分析化验资料，对樊142块沙三下亚段进行沉积微相划分。研究认为，樊142块沙三下亚段总体为深湖相，可进一步划分为物源砂体亚相、滑塌浊积砂体亚相和深湖泥(灰)岩亚相。物源砂体亚相和滑塌浊积砂体亚相可进一步划分为主干水道微相、水道侧缘微相和席状砂微相。

然后，根据沉积微相研究所得的单井相及平面相沉积微相划分成果，采用确定性建模的方法，建立樊142块沙三下亚段沉积微相模型(图 1)。物源位于研究区西南方向，研究区西南部为物源砂体亚相，由物源砂体向东北方向滑塌而成的滑塌浊积砂体亚相主要位于研究区的中部和西北部。

最后，在沉积微相模型的基础上，利用单井的层速度数据，采用相控建模的方法，建立樊142块三维层速度模型。樊142块沙三下亚段地层速度主要为2 700~3 500 m/s，平均速度约为3 170 m/s。其中，滑塌浊积砂体的速度最大，为3 400~3 800 m/s。沙三下亚段上部岩性为泥岩、白云岩、油页岩薄互层，层速度较大，主要为3 300~3 500 m/s；沙三下亚段中部除滑塌浊积砂体外，主要为油页岩、泥岩、灰质白云岩，层速度较小，为2 700~3 200 m/s；沙三下亚段下部岩性以厚层泥岩、油页岩为主，夹薄层白云岩，层速度较大，为3 300~3 500 m/s。

3 樊142块岩石力学参数模型

地质体的非均质性是通过不同的岩石力学参数体现在ANSYS力学模型中的，岩石力学参数与层速度有定量的转换关系^[15-23]。在三维层速度模型的基础上，结合纵横波转换关系和动静态杨氏模量、动静态泊松比的转换关系，可以得到研究区三维岩石力学参数展布^[15-22]。

不同区域不同层段的纵横波转换关系存在一定的差异。根据研究区位置、目的层段及前人的研究成果，选取了济阳坳陷古近系的纵横波速度转换关系。马中高^[23]于2008年利用济阳坳陷古近系不同层系、不同岩性、不同沉积环境的344块样品拟合出了济阳坳陷古近系不同岩性的纵横波速度关系式，相关系数都在0.95以上。

济阳坳陷古近系泥岩纵横波速度转换关系式为

式中：v_P为纵波速度，km/s；v_S为横波速度，km/s。

济阳坳陷古近系砂岩纵横波速度转换关系式为

根据弹性理论，利用由式(1)和式(2)计算得出的砂、泥岩纵横波速度，可计算得到樊142块沙三下亚段的动态杨氏模量和泊松比：

式中：E_d为动态杨氏模量，MPa；μ_d为动态泊松比；ρ为岩石密度，g/cm³。

根据前人在本地区对实测岩石力学参数与常规测井资料计算岩石力学参数的回归计算^[23-25]，得到本地区的动静态岩石力学参数转换公式：

式中：E_s为静态杨氏模量；μ_s为静态泊松比。

根据式(5)和式(6)，可将樊142块沙三下亚段的动态杨氏模量和泊松比转化为静态的杨氏模量和泊松比，从而得到樊142块沙三下亚段的岩石力学参数模型。

4 ANSYS三维构造应力场模拟 4.1 力学模型构建及边界条件确定

如何将PETREL三维构造模型全真转换为ANSYS力学模型，是准确模拟构造应力场的关键。本次层面模型转换的思路是，首先从PETREL构造模型中的层面中提取点云数据，然后通过三维点云逆向工程将点云数据拟合为ANSYS实体，来实现层面由PETREL构造模型到ANSYS力学模型的转换。断层模型则根据PETREL三维构造模型中断层的平面形状、走向和位置以及剖面上的倾向、倾角数据在ANSYS中直接建模，并根据断层断距的大小设置断层的厚度，使ANSYS中的断层模型更加符合地质规律。

考虑到实际地质模型边界的复杂性，难以在实际目的层段模型的外部施加复杂的载荷边界条件；再者，为减小边界条件对研究区的影响，在研究区模型外部建立三维尺寸均大于目的层的立方体。ANSYS几何模型由三个部分组成：3个地层实体(沙三下亚段被浊积砂体一分为三)、11条断层以及外部立方体。

4.2 模型离散化及岩石力学参数加载

鉴于Solid186是ANSYS中一个高阶3维20节点固体结构单元，具有二次位移模式，可以更好地模拟不规则的网格模型，因此岩体区域和断层均采用Solid186单元。根据研究精度的要求，兼顾计算效率，将外部立方体离散为尺寸较大的单元(棱长小于300 m)，目的层樊142块沙三下亚段离散为小尺寸单元(棱长小于30 m)。由于断层及其两盘的应力场分布规律较为复杂^[24]，因此对断层离散单元尺寸进行了再次细化(棱长小于10 m)，模型共划分出3 200 674个单元、4 385 999个节点(图 2)。

完成模型的离散化之后，采用按距离扫描的方法(即逐个计算某单元与其周边岩石力学参数的距离，选择与其距离最小的岩石力学参数加载到该单元上)，将三维岩石力学参数加载到ANSYS有限元 分析模型的每一个单元中(图 3,当岩石力学参数多于ANSYS中颜色种类时，其默认用特定的颜色序列循环对不同编号的力学参数着色)。本次模拟计算共加载19 367种岩石力学参数，将浊积砂储层垂向和平面的非均质性全面体现到ANSYS模型中。

4.3 边界条件反演及模拟结果检验

樊142块浊积砂储层内部应力分布十分不均匀，应力大小和方向变化范围都较大。根据樊142块已压裂施工22口井的压裂施工曲线及裂缝监测报告，可以得到该22口井压裂段的水平最大、最小主应力的大小和方向：井点处水平最小主应力为25.8~33.8 MPa，水平最大主应力为28.82~42.71 MPa；水平最大主应力方向为NE66.3°—NE92.9°。樊142块沙三下亚段现今地应力场以近东西向挤压为主。以22口井处实测的应力大小和方向为约束，经过多次试算之后，最终确定边界条件为东部边界施加梯度压力p_x=35~25 MPa，南部边界上施加梯度挤压载荷p_y=45~25 MPa，垂直方向根据深度施加等效的岩体重力载荷^[26-28]。

通过有限元模拟计算后，可以得到樊142块现今应力分布特征，包括最大主应力、中间主应力、最小主应力以及最大剪应力的大小和方向等。提取11口井井点处的应力模拟结果(表 2)与实测的应力大小和方向进行对比，最大水平主应力值标准误差为4.35%，最小水平主应力值标准误差为3.26%，水平最大主应力方向的标准误差为5.6%。

5 结果分析

樊142块构造应力场模拟结果表明，樊142块第一主应力是垂向应力，第二主应力是水平最大主应力，第三主应力是水平最小主应力。最大水平主应力大小整体由西南向东北逐渐增大，最大值为45 MPa，最小值为25 MPa。纵向上表现为随深度增加应力值增大，下亚段顶底面应力差值约为1 MPa。樊142块大部分区域以近东西方向水平挤压为主，最大水平主应力方向由西南部约NE60°向西北部区域NE75°—NE85°过渡。受构造及岩石力学参数影响，构造应力场局部存在异常(图 4)。

断层及周边应力值大小和方向受断层的影响显著。断层的不同部位、规模、组合形式及其走向，都会引起断层周边的应力场异常(图 5)。距离断层端部越近、断层规模越大、断层走向与主应力方向差异越大，其引起的断层周边应力场异常程度与范围越大。随着与断层距离的减小，断层两盘的最大主应力值都表现出先缓慢增大、后急剧减小的特点，且断层上的应力值最小(图 6)；断层的最大主应力方向与其走向一致，且随着与断层距离的减小，断层两盘的最大主应力方向都趋近于断层的最大主应力方向(图 7)。

浊积砂体的岩石力学参数变化较大，水平最大主应力值的压力梯度分界线也因此参差不齐。在区域水平最大主应力低值的背景下，存在局部的应力高值，如F142-X309井、F142-X302井附近(图 5)。

6 结论

1) 通过三维点云逆向工程将PETREL构造模型转换为ANSYS几何模型的方法，能够实现形态复杂的浊积砂储层模型与地质资料的吻合。

2) 根据三维层速度模型计算三维岩石力学参数，用距离扫描法将其加载到ANSYS有限元分析模型的每个单元中，能够较全面地反映储层内部岩石力学性质的差异。

3) PETREL与ANSYS联合建模技术能充分结合PETREL地质建模和ANSYS力学模拟方面的优势，将浊积砂体的岩相、岩性、物性等非均质特征完整地包含在ANSYS力学模型中。

4) 断层及周边应力值大小和方向受断层的影响显著，断层的不同部位、规模、组合形式及其走向都会引起断层周边应力场异常；岩石力学参数对局部应力场的异常也有较显著的影响。

5) 樊142块沙三段构造应力场模拟结果表明，对于构造形态复杂、厚度变化大、储层非均质性强的浊积砂低渗透储层，本文提供的方法能够准确地模拟其现今构造应力场，并且体现出受浊积砂储层构造形态和储层非均质性影响而引起的局部构造应力场的异常，为油田开发中的井网布置、水力压裂和优化开采提供一定的参考。