引用本文
位云生, 王军磊, 于伟, 齐亚东, 苗继军, 袁贺, 刘楚溪. 基于三维分形裂缝模型的页岩气井智能化产能评价方法[J]. 石油勘探与开发, 2021,48(4): 787-796.
WEI Yunsheng, WANG Junlei, YU Wei, QI Yadong, MIAO Jijun, YUAN He, LIU Chuxi. A smart productivity evaluation method for shale gas wells based on 3D fractal fracture network model[J]. Petroleum Exploration & Development, 2021,48(4): 787-796.  
Doi: 10.11698/PED.2021.04.11
Permissions
Copyright©2021, 《石油勘探与开发》编辑部
《石油勘探与开发》编辑部 所有
基于三维分形裂缝模型的页岩气井智能化产能评价方法
位云生1, 王军磊1, 于伟2,3, 齐亚东1, 苗继军3, 袁贺1, 刘楚溪2
1. 中国石油勘探开发研究院,北京 100083
2. 德克萨斯大学奥斯汀分校,德克萨斯 78712,美国
3. SimTech有限责任公司,休斯顿 77494,美国
联系作者简介:王军磊（1986-）,男,山东文登人,博士,中国石油勘探开发研究院工程师,主要从事复杂气藏渗流仿真模拟、解析/数值解分析和机器学习方面的研究工作。地址：北京市海淀区学院路20号,中国石油勘探开发研究院气田开发研究所,邮政编码：100083。E-mail: wangjunlei@petrochina.com.cn

第一作者简介：位云生（1979-）,男,河南项城人,博士,中国石油勘探开发研究院高级工程师,主要从事气田开发综合研究工作。地址：北京市海淀区学院路20号,中国石油勘探开发研究院气田开发研究所,邮政编码：100083。E-mail: weiys@petrochina.com.cn

收稿日期: 2020-11-23 修回日期: 2021-05-07
基金项目: 国家科技重大专项“昭通页岩气勘探开发示范工程”（2017ZX05063-005）; 中国石油勘探开发研究院科技开发项目“复杂缝网条件下页岩气井生产动态数值模拟程序开发”（YGJ2019-12-04）
摘要

基于分形理论多次信息叠加算法设计了三维天然裂缝生成方法,人工裂缝模型与天然裂缝模型耦合可表征页岩压后复杂多尺度裂缝系统。采用马尔科夫链-蒙特卡洛智能化历史拟合算法,结合嵌入式离散裂缝技术耦合数值模拟器预测页岩气井产能,形成了一体化页岩气井产能评价模型。研究表明：三维天然裂缝生成方法可利用分形参数控制裂缝网络的整体分布,与人工裂缝耦合可表征页岩压后复杂的跨尺度裂缝系统;嵌入式离散裂缝模型具有裂缝网格数量少、运算耗时短的优点,能够灵活表征天然裂缝及人工裂缝属性,在有效降低计算量的同时能够精确地模拟流体在基质-裂缝中的交换过程;嵌入式离散裂缝模型与智能化历史拟合算法相结合能够降低裂缝、储集层等未知参数计算的不确定性,实现储集层、裂缝关键参数的有效反演,并实现气井产能量化预测。经实际应用验证一体化井产能评价模型预测结果具有较高的可信度。 图13 表7 参26

关键词: 分形裂缝网络; 多次信息叠加算法; 嵌入式裂缝模型; 智能化历史拟合; 储集层参数反演; 页岩气; 智能化产能评价
中图分类号:TE329 文献标志码:A 文章编号:1000-0747(2021)04-0787-10
A smart productivity evaluation method for shale gas wells based on 3D fractal fracture network model
WEI Yunsheng1, WANG Junlei1, YU Wei2,3, QI Yadong1, MIAO Jijun3, YUAN He1, LIU Chuxi2
1. Research Institute of Petroleum Exploration & Development, PetroChina, Beijing 100083, China
2. University of Texas at Austin, Austin TX 78712, USA
3. SimTech LLC, Houston TX 77494, USA
Abstract

The generation method of three-dimensional fractal discrete fracture network (FDFN) based on multiplicative cascade process was developed. The complex multi-scale fracture system in shale after fracturing was characterized by coupling the artificial fracture model and the natural fracture model. Based on an assisted history matching (AHM) using multiple-proxy-based Markov chain Monte Carlo algorithm (MCMC), an embedded discrete fracture modeling (EDFM) incorporated with reservoir simulator was used to predict productivity of shale gas well. When using the natural fracture generation method, the distribution of natural fracture network can be controlled by fractal parameters, and the natural fracture network generated coupling with artificial fractures can characterize the complex system of different- scale fractures in shale after fracturing. The EDFM, with fewer grids and less computation time consumption, can characterize the attributes of natural fractures and artificial fractures flexibly, and simulate the details of mass transfer between matrix cells and fractures while reducing computation significantly. The combination of AMH and EDFM can lower the uncertainty of reservoir and fracture parameters, and realize effective inversion of key reservoir and fracture parameters and the productivity forecast of shale gas wells. Application demonstrates the results from the proposed productivity prediction model integrating FDFN, EDFM and AHM have high credibility.

Keyword: fractal discrete fracture network; multiplicative cascade process; embedded discrete fracture model; intelligent history matching; reservoir parameter inversion; shale gas; smart productivity evaluation
0 引言

页岩气储集层体积压裂后形成的裂缝网络十分复杂, 具有跨尺度特征, 裂缝描述及高效模拟难度大, 页岩气井压后产能预测难度大。要获得准确的产能评价结果, 必须具备3个条件：①合理表征裂缝网络的空间展布特征; ②精确模拟流体在裂缝网络中的渗流过程; ③高效准确地拟合气井生产历史数据。

页岩中天然裂缝具有离散裂缝介质的特征, 传统等效连续模型难以适用[1, 2]。离散裂缝网络（DFN）模型将裂缝描述为具有不同尺度和形态的裂缝片, 采用随机模拟方法, 利用蚂蚁体追踪算法结合井震资料获得实际裂缝网络[3, 4, 5], 但该模型受制于随机模拟方法的局限性, 关键表征参数存在着很大的不确定性。由于天然裂缝的几何构型具有分形特征[6, 7, 8]（即裂缝特征参数是关于分形维数的统计概率场）, 相较于DFN模型所依据的均匀随机理论, 分形几何理论使用分数维来描述随机的、无规则的但具有统计意义上自相似性的裂缝系统[7], 所生成的分形离散裂缝模型（FDFN）可以更准确地描述天然裂缝的非均质性及几何形态的复杂性。对于相同的分形参数, 裂缝系统可能有不同的分布规律及几何形态, 但却具有几乎相同的储集层渗流能力, 即相同的缝网等效渗透率[8, 9], 其对应的压力动态曲线也一致[10, 11]。天然裂缝所具有的自相似性和尺度不变性, 使得用分形参数可以很好地控制天然裂缝的整体分布, 在提高刻画精度的同时又保证了内含的关键数学信息不丢失。

受制于局部网格加密[12]、等效连续模型[13, 14]和非结构网格[15, 16]技术在特征尺度、网格剖分量、模拟精度和计算效率等方面的局限性, 近年来兴起的嵌入式离散裂缝模型（EDFM）技术避免了非结构化网格剖分的复杂性以及伴随的算法不收敛性[17, 18], 同时保证了裂缝-裂缝、裂缝-基质间的流动模拟精度[19, 20], 在高效模拟任意复杂裂缝网络产能方面具有明显的技术优势。此外, 基于产能模型的历史拟合是进行参数反演及产能预测的必要环节, 基于普通优化算法（如牛顿优化算法、粒子群优化算法、遗传算法等）的自动历史拟合流程繁琐、耗时、耗力, 难以求取最有代表性的拟合结果[21]。而基于智能化算法的辅助性历史拟合技术能够充分考虑到产能主控因素间的强关联性, 提高历史拟合效率、减少人为主观因素干扰, 但由于模拟本身的复杂性, 相适应的拟合算法研发极具挑战性。

本文提出了从裂缝建模、动态模拟、历史拟合到产能预测的产能评价方法, 包括：①使用分形几何理论建立离散天然裂缝表征模型; ②以生成的三维天然/人工裂缝模型为基础, 结合新一代嵌入式离散裂缝模型, 精确模拟页岩气井全生命周期生产动态; ③通过智能化算法自动取样进行机器学习, 实现关键储集层和裂缝参数的有效反演, 并对裂缝参数及产能进行风险量化。

1 三维分形裂缝模型
1.1 分形离散裂缝网络生成

分形理论分为二维分形与三维分形, 每个维度模型有其特定的参数。裂缝分形维数可以通过露头、成像测井和岩心分析获得。裂缝中心点、开度、长度、走向、分布等参数能够控制裂缝系统的成簇分布。目前国际上均为二维分形裂缝系统, 三维裂缝采用圆盘模型, 但这与实际裂缝形态并不一致[9]。

三维空间对应的欧几里得维数为3。裂缝长度在l~l+dl之间的裂缝条数记为$n\left( l, L \right)\text{d}l$, 满足[6]：

$n\left( l, L \right)\text{d}l={{\alpha }_{3\text{D}}}{{L}^{{{D}_{\text{c3D}}}}}{{l}^{-\left( {{D}_{\text{l3D}}}+1 \right)}}\text{d}l$ （1）

α 3D决定三维空间中的裂缝数量, 与模型尺度无关, 需要计算获得。根据Piggott公式[22], 可获得二维密度常数和三维密度常数间的关系式：

${{\alpha }_{\text{3D}}}=\alpha \sqrt{\text{ }\!\!\pi\!\!\text{ }}\frac{\Gamma \left( 1.5+0.5{{D}_{\text{l}}} \right)}{\Gamma \left( 1.0+0.5{{D}_{\text{l}}} \right)}$ （2）

通常, 引入参数α 3D用以表征长裂缝与短裂缝的比例, Darcel等[22]命名为幂率长度分布指数, 当α 3D=3时长裂缝数量占主导地位, 当α 3D=∞ 时裂缝长度一致且趋近于最小值lmin。其与Dl3D之间满足如下关系式[3]：

${{\alpha }_{\text{3D}}}={{D}_{\text{l3D}}}+1$ （3）

同时二维分形维数与三维分形维数的关系满足[6]：

$\left\{ \begin{align} & {{D}_{\text{c3D}}}={{D}_{\text{c}}}+1 \\ & {{D}_{\text{l3D}}}={{D}_{\text{l}}}+1 \\ \end{align} \right.$ （4）

在（4）式中, Dc取值范围为1~2; Dl取值范围为1~∞ 。对（1）式进行积分, 得到三维模型中所有长度大于lmin的裂缝总数量[9]：

$N\left( L \right)=\frac{{{\alpha }_{\text{3D}}}}{{{D}_{\text{l3D}}}}{{L}^{{{D}_{\text{c3D}}}}}l_{\min }^{-{{D}_{\text{l3D}}}}$ （5）

（5）式同样适用于二维分形裂缝数量计算。裂缝几何中心点在二维平面空间中的分布可以使用分形维数Dc和成对校正函数${{C}_{2}}\left( r \right)$来计算[23]：

${{C}_{2}}\left( r \right)=\frac{2{{N}_{\text{p}}}\left( r \right)}{N\left( N-1 \right)}=c{{r}^{{{D}_{\text{c}}}}}$ （6）

需要强调的是, 适用于二维分形裂缝几何参数生成的多次信息叠加算法也适用于三维裂缝模型, 该算法主要应用于天然裂缝区域性分布概率计算。二维模型经过多次迭代分解而形成多个子区域（7~9次迭代）, 每个区域出现天然裂缝的概率由特定公式约束[6, 22], 所有子区域的概率总和不一定等于1（否则是完全随机分布）。多次信息叠加迭代次数满足[9]：

${{N}_{\text{iter}}}=\frac{\lg \frac{L}{{{l}_{\min }}}}{\lg {{l}_{\text{ratio}}}}$ （7）

本文基于Python语言研发了三维分形裂缝生成模块, 该模块能够生成裂缝属性参数。图1展示了经过多次信息叠加过程生成的天然裂缝中心点的三维分布, 图中每个点代表天然裂缝中心点的空间位置, 具体参数给定为：Dc3D=2.8, Dl3D=3.5, L=1 000 m, α =0.5, lmin=50 m, Niter=6。

图1 三维天然裂缝中心点的位置分布

Kim等[9]给出的三维分形模型, 假设裂缝为圆盘形状, 模拟区为长宽高相等的箱体, 这与实际裂缝形态不一致, 故本文采用等效面积方法将圆盘模型校正为矩形模型。由于天然裂缝的长度遵从幂律分布, 而非正态分布, 因此使用幂律分布模拟天然裂缝的长度[11]：

$r={{\left[ l_{\min }^{-{{D}_{\text{c}}}}-F\left( l_{\min }^{-{{D}_{\text{c}}}}-l_{\max }^{-{{D}_{\text{c}}}} \right) \right]}^{-\frac{1}{{{D}_{\text{c}}}}}}$ （8）

α 3D和Dc3D参数控制着三维空间中裂缝系统的分布特征：①α 3D参数决定长裂缝与短裂缝数量比例（取值3~∞ ）。当α 3D=3.0时模型中长裂缝数量占主导地位, 随α 3D值的增加模型中长裂缝数量占比降低, 当α 3D趋近于无穷大时模型中所有裂缝长度一致且等于最小裂缝长度lmin; ②Dc3D参数决定了裂缝分布的非均质性（取值2.0~3.0）, 当Dc3D=2时裂缝分布不均匀且成簇出现, 空白区域多, 随着Dc3D值增加, 裂缝分布趋于均匀; ③极限条件下, 当α 3D=∞ , Dc3D=3.0时, 裂缝系统长度一致且均匀分布, 此时等效于连续性模型。

图2给出了两组裂缝条件下天然裂缝长度、走向及倾角的统计结果。可以看到裂缝长度基本呈偏正态分布, 而天然裂缝方位角、倾角由高斯正态分布决定。

图2 天然裂缝长度、走向及倾角统计分布图

将三维裂缝描述为矩形, 则对应的裂缝高度由0~1中生成的随机数与天然裂缝的主长度的乘积决定。采用裂缝中心点、长度、高度、走向及倾角等几何参数的计算结果, 以裂缝中心点为基点, 根据走向及角度确定裂缝空间形态, 根据长度及高度确定裂缝几何尺寸, 最终可获得三维分形裂缝几何形态的空间分布图（见图3）。

图3 分形理论产生的三维分形天然裂缝模型

1.2 裂缝生成算法验证

三维分形离散裂缝模型是基于二维模型算法改进的, 为了便于直观分析, 这里以二维模型为例进行验证。控制二维分形离散裂缝模型的参数包括二维裂缝分布的分形维数Dc、二维裂缝长度分布的分形维数Dl、二维裂缝密度常数α 以及裂缝方位角中值。算例中设定生成区域大小为10 m× 10 m, 最小裂缝长度lmin=0.05 m, 其余所需分形参数值见表1中的模型输入值。输入参数, 通过分形裂缝生成模型, 获得分形离散裂缝网络（见图4）, 模型共生成9 326条裂缝。

图4 分形离散裂缝网络

表1 模型输入参数与拟合结果参数对比表

根据图4中生成的裂缝分布情况, 统计裂缝方位角、长度和位置的分布（见图5）, 采用二维条件下的（5）式（形式如$\frac{N\left( L \right)}{{{L}^{{{D}_{\text{c}}}}}}=\frac{\alpha }{{{D}_{\text{l}}}}l_{\min }^{-{{D}_{\text{l}}}}$）、（6）式, 通过数据拟合重新获得裂缝分形参数（见表1）。模拟结果表明, Dc、Dl及裂缝方位角中值的模型输入值和拟合结果基本一致, α 计算结果虽然与输入值有一定差异, 但仍在误差范围之内, 说明本文分形裂缝生成算法是可靠的, 即根据可靠的二维分形离散裂缝生成算法能够获得较精确的三维分形裂缝网络模型[9, 11]。

图5 算例生成分形离散裂缝网络及裂缝属性统计

2 裂缝网络流动模拟
2.1 数值模型及假设条件

在生成的三维分形天然裂缝系统上, 设定分段压裂水平井为基本井型, 表征天然裂缝与人工裂缝耦合的压后复杂裂缝系统（见图6）, 并将裂缝系统与数值模拟模型进行结合。

图6 页岩压后天然-人工裂缝物理模型（二维投影）

采用黑油模型进行流动模拟, 考虑页岩气储集层地质、开发特征, 对模型做如下假设。

①储集层均质、水平、等厚, 平面基质渗透率各向同性, 考虑到层理影响气体纵向流动, 纵向渗透率为平面渗透率的0.1倍。

②气藏采用三维单孔数值模型, 模型包括3个不同尺度空间：基质、天然裂缝和人工裂缝。基质系统采用规则正交网格, 裂缝系统（天然、人工裂缝）使用EDFM产生的虚拟网格描述, 该方法最大优点在于考虑了各类孔隙间的流体传递特征。

③为简化气体解吸过程, 采用瞬时解吸模型描述吸附气解吸过程, 吸附模型选用BET（Brunauer- Emmett-Teller）多层等温吸附模型：

${{v}_{\text{ag}}}=\frac{1-\left( n+1 \right){{\left( \frac{p}{{{p}_{\text{o}}}} \right)}^{n}}+n{{\left( \frac{p}{{{p}_{\text{o}}}} \right)}^{n+1}}}{1+\left( C-1 \right)\frac{p}{{{p}_{\text{o}}}}-C{{\left( \frac{p}{{{p}_{\text{o}}}} \right)}^{n+1}}}\frac{{{v}_{\text{m}}}\frac{Cp}{{{p}_{\text{o}}}}}{1-\frac{p}{{{p}_{\text{o}}}}}$ （9）

上式中, C和n为无因次实验常数。当n=1时, BET多层模型简化为Langmuir等温吸附模型：

${{v}_{\text{ag}}}={{v}_{\text{L}}}\frac{p}{p+{{p}_{\text{L}}}}$ （10）

其中 vL=vm, pL=po/C。

④考虑生产过程中的压裂液返排、裂缝内流体流动满足气水两相流动, 故选用Corey型相渗模型：

${{K}_{\text{rw}}}={{K}_{\text{rwo}}}{{\left( \frac{{{S}_{\text{w}}}-{{S}_{\text{wirr}}}}{1-{{S}_{\text{wirr}}}-{{S}_{\text{gr}}}} \right)}^{{{N}_{\text{w}}}}}$ （11）

${{K}_{\text{rg}}}={{K}_{\text{rgo}}}{{\left( \frac{1-{{S}_{\text{w}}}-{{S}_{\text{gr}}}}{1-{{S}_{\text{wirr}}}-{{S}_{\text{gr}}}} \right)}^{{{N}_{\text{g}}}}}$ （12）

在（11）— （12）式中, Krwo、Krgo、Nw、Ng为实验拟合值。模型采用非平衡态初始化方法, 裂缝内含水饱和度较高, 基质内含水饱和度低于束缚水饱和度。

⑤考虑人工裂缝系统内渗透率应力敏感效应, 渗透率随压力变化, 满足衰减公式：

$K={{K}_{\text{i}}}\exp [-\gamma ({{p}_{\text{i}}}-p)]$ （13）

⑥考虑流体重力及相态间毛管压力。本文中视气、水两相为无传质现象的两个组分。根据达西定律, 组分jgw（jgw=g代表气相, jgw=w代表水相）的体积流速满足：

${{u}_{{{j}_{_{\text{gw}}}}}}={{10}^{-6}}\frac{K{{K}_{\text{r, }{{j}_{\text{gw}}}}}}{{{u}_{{{j}_{\text{gw}}}}}}\left( \nabla {{p}_{{{j}_{\text{gw}}}}}-{{10}^{-6}}{{\gamma }_{{{j}_{\text{gw}}}}}\nabla D \right)$ （14）

裂缝和基质的流动特性通过K和${{K}_{\text{r, }{{j}_{\text{gw}}}}}$来区别表征。

⑦将（9）— （14）式与质量守恒方程结合, 获得气相和水相的压力控制方程：

$\left\{ \begin{align} & \frac{\partial }{\partial t}\left[ \phi {{S}_{\text{g}}}\left( 1-{{c}_{\text{ag}}} \right){{\rho }_{\text{g}}}+\left( 1-\phi \right){{\rho }_{\text{s}}}{{v}_{\text{ag}}}{{\rho }_{\text{ag}}} \right]+\nabla \cdot \left( {{\rho }_{\text{g}}}{{u}_{\text{g}}} \right)=\frac{{{q}_{\text{g}}}}{{{V}_{\text{b}}}} \\ & \frac{\partial }{\partial t}\left( \phi {{S}_{\text{w}}}{{\rho }_{\text{w}}} \right)+\nabla \cdot \left( {{\rho }_{\text{w}}}{{u}_{\text{w}}} \right)=\frac{{{q}_{\text{w}}}}{{{V}_{\text{b}}}} \\ \end{align} \right.$ （15）

使用正交网格对数学模型（（15）式）进行有限差分处理。本文使用自主研发的离散裂缝模型（EDFM）技术对裂缝连接进行处理[24], 在使用正交网格基础上, EDFM技术能够在保持规则网格计算高效性的同时, 通过等效处理方法快速模拟复杂裂缝系统。

为了保证处理前后裂缝体积的一致性, 虚拟的裂缝网格孔隙度应满足：

${{\phi }_{\text{f}}}=\frac{{{V}_{\text{f}}}}{{{V}_{\text{b}}}}=\frac{{{S}_{\text{seg}}}{{w}_{\text{f}}}}{{{V}_{\text{b}}}}$ （16）

EDFM技术核心在于非邻近连接对（NNC）的传导率计算, 主要用于处理物理模型上相邻但在计算域上不相邻网格之间的流量交换。非邻近网格连接对的传导率TNNC计算通式为：

${{T}_{\text{NNC}}}\text{=}\frac{{{K}_{\text{NNC}}}{{A}_{\text{NNC}}}}{{{d}_{\text{NNC}}}}$ （17）

上式中, 当裂缝与基质连接时, dNNC为基质块到裂缝面的平均距离, KNNC为基质渗透率; 当裂缝与裂缝连接时, dNNC为裂缝元之间的法向距离, KNNC为裂缝平均渗透率。

2.2 对比验证

通过对比局部加密网格（LGR）、非结构化网格（PEBI）和嵌入式离散裂缝模型（EDFM）的计算结果验证模型可靠性及计算效率。考虑LGR方法的局限性, 使用规则的分段压裂水平井作为算例, 气藏尺寸1 700 m× 400 m× 20 m, 纵向为单层网格且假设裂缝贯穿, 裂缝开度为0.01 m, 孔隙度为35%, 气井保持1.5 MPa恒定井底压力, 单相气体生产20年, 其他基础参数见表2。

表2 模型输入基础参数

LGR和EDFM模型中基础网格数为170× 40× 1。在LGR模型中, 垂直裂缝面方向采用对数加密原则剖分网格, 沿裂缝面方向采用均匀加密网格; PEBI网格在裂缝附近采用径向网格局部加密, 并根据裂缝周围流场构建网格。EDFM模型中裂缝以虚拟网格的形式与原始基础网格连接, 裂缝真实属性通过显式公式描述, 不需要做等效处理。

图7对比了LGR、PEBI和EDFM方法模拟的气井生产动态数据, 其中EDFM和PEBI结果吻合度极高。因受制于网格剖分程度LGR方法存在一定误差, 但仍在误差范围之内。需要强调的是, LGR模拟结果与网格剖分程度有关, 裂缝附近网格越密, 计算精度越高, 但模型计算量也大幅度增加[19]。

图7 EDFM方法与PEBI、LGR方法计算结果对比

表3对比了EDFM、PEBI和LGR方法的计算效率。LGR方法需要大量裂缝网格描述裂缝周边的流动特征, 由于正交网格本身的局限性, 该方法难以刻画复杂裂缝网络, 而且当储集层渗透率较低时需要更为精细的网格, 大大增加了计算量; PEBI方法虽然整体网格数最少, 但刻画裂缝所需网格量仍较大, 而且当描述复杂裂缝网络时, 需要更多的细化网格来计算裂缝交点间的流体交换, 严重影响计算效率。对比来看, EDFM方法在保证计算精度的同时能够有效降低裂缝网格数, 显著节约模拟成本、提高模拟效率, 尤其对构型复杂的裂缝网络, EDFM方法计算优势更为明显。

表3 EDFM和LGR、PEBI方法计算效率对比
3 智能化历史拟合

气藏自动历史拟合本质是求取目标函数最小值的方法。利用EDFM方法处理复杂裂缝网络的强大能力, 结合马尔科夫链-蒙特卡洛（MCMC）机器学习算法, 开发出基于神经网络（NN）的智能化历史拟合模块（NN-MCMC）[25], 实现了复杂裂缝系统（包括有效缝高、缝长、导流能力、簇效率、裂缝的应力敏感程度等关键参数）的高效精确评估。

图8为智能化历史拟合流程。首先, 对气井实测数据（如产气/水速率、井口压力等）进行整理。其次, 根据工程/地质资料设定未知参数的取值范围。利用正交实验法中的采样方法, 对参数组合进行初始迭代采样, 获得50个样本。通过模块化输入文件和EDFM模型, 同步模拟初始采样的50个气藏样本, 获得模拟结果, 随后选取两个输出参数（如井底流压、气水比或产水速率）进行全局误差计算。需要强调的是, 需选取特定时间点的数据（称之为数据点索引号）, 进行实际数据与模拟输出的比对, 如（18）式：

$\varepsilon =\frac{\sum\nolimits_{j=1}^{P}{\sum\nolimits_{i=1}^{Q}{\left| 100\frac{{{x}_{\text{model, }i\text{, }j}}-{{x}_{\text{history, }i\text{, }j}}}{{{x}_{\text{history, }i\text{, }j}}} \right|{{w}_{i, j}}}}}{\sum\nolimits_{j=1}^{P}{\sum\nolimits_{i=1}^{Q}{{{w}_{i, j}}}}}$ （18）

图8 智能化历史拟合工作流程示意图

在获得初始循环的结果后, 需运用人工智能-神经网络对50组输入变量以及输出变量建立具有关联性的代理模型。进而, 使用马尔可夫链-蒙特卡洛取样方法对输入变量进行系统性取样[25], 保证反演组合具有代表性, 且可优化全局误差值。取样完成后, 可获得数千个样本组合, 使用训练好的人工智能-神经网络代理模型对代表性样本组合进行模拟预测。再次通过（18）式计算获得数千个具有代表性样本的全局误差, 并选取25个全局误差最小的样本作为下一个迭代的输入值进行模拟运算。重复此步骤, 直至代理模型的预测具有收敛性或者达到最大迭代数量。

通过多次迭代, 此方法可确保代理模型的精度随着迭代数量的增加而提高。在运行完所有迭代后, 需要对比计算目标函数（即选定的两个输出参数）的误差来筛选历史拟合解。目标函数的计算方法类似于（18）式, 不同点在于每个拟合参数均单独计算误差, 而不对拟合参数进行求和得到目标函数误差。通过设置每个目标函数的筛选阈值, 多个符合误差要求的历史拟合解被筛选出来。最后进行统计分析, 并模拟长期、具有代表性的产能预测结果。

4 实例

以昭通地区页岩气A水平井为例, 该井区天然裂缝以近南北、北西和北东向走向为主, 倾角平均80° 。露头区裂缝延伸长度测量数据表明, 多数裂缝延伸长度小于100 m, 故实例中根据模拟需要将裂缝片长度设置为0~100 m, 平均约50 m, 裂缝片长高比为2∶ 1。该井控制范围内的天然裂缝分形参数如表4所示。该井区页岩孔隙度5.5%~6.5%, 原始基质渗透率约为（1.0~10.0）× 10-8 μ m2, 其他地质及工程参数如表5所示。该井经过分段压裂、返排后于2015年8月投产。

表4 天然裂缝基础参数
表5 数值模型基础数据

使用压裂模拟软件对该井进行水力压裂模拟, 得水力裂缝参数统计结果如图9所示, 其中水力缝长49~712 m, 平均321 m; 支撑缝长4~606 m, 平均231 m; 水力缝高2~107 m, 平均58 m; 支撑缝高0~105 m, 平均32 m。

图9 水力压裂缝模拟统计结果

由于川南地区两向水平主应力差较大, 模拟出的人工裂缝整体呈平直状, 近似于标准分段压裂水平井, 因此使用裂缝模型参数平均值设定人工裂缝几何参数。

至此, 建立了由基质、天然裂缝、人工裂缝组成的页岩气井产能评价模型。通过修正部分模型参数来拟合生产历史, 重点拟合日产气、日产水量以及井底流压等开发指标。根据生产动态数据采用自动历史拟合和反演技术, 有效开展压后效果评估, 降低核心参数的不确定性, 表6为参与拟合参数取值的上下限。

表6 地层/裂缝参数取值界限

本算例中进行12步自动迭代（每步迭代取1 000万个样本）, 从中选取325个气藏模型进行历史拟合, 总误差随迭代次数增多而降低。设定全局误差ε < 50%, 从中优选出75套历史拟合解, 拟合效果如图10所示。

图10 历史拟合解与历史生产数据对比图

参与自动历史拟合的不确定参数组合可用平行坐标图表示, 图11为全部325次自动历史拟合模拟结果。横坐标代表每个不确定参数, 纵坐标代表每个不确定参数的范围, 红线代表最优解的组合, 橘黄线代表所有优选的历史拟合解组合（75次）, 灰线代表非历史拟合解的组合（250次）。

图11 最优解、历史拟合解和非历史拟合解对应的参数对比

通过对拟合参数进行统计分析, 相应参数的P50值如表7所示。历史拟合结束后预测20年内气井累计产水量和累计产气量。根据获得的解释参数, 采用定井底流压的方式预测气井产能（井底压力为1.5 MPa）, 如图12所示。图中产能预测曲线为不同模型参数组合下的结果, 曲线簇包含了实际动态曲线。图13为20年最终采出量预测结果的统计分析, 可以看到在P10~P90的置信域内, 气井累计产气（1.05~1.22）× 108 m3, 最优解1.15× 108 m3; 气井累计产水（0.42~0.81）× 104 m3, 最优解0.72× 104 m3。联合多种产量递减模型的生产数据分析方法是页岩气井产能评价的一种重要手段, 采用该方法获得的单井累计产气量1.01× 108 m3, 累计产水量0.81× 104 m3[26], 两种方法的解释结果高度近似, 证实本方法预测结果具有较高的可信度。

表7 拟合参数统计结果的P50值

图12 不同参数组合下的气井产量拟合及预测结果

图13 不同参数组合下气井生产20年最终采出量评价结果

5 结论

基于分形理论多次信息叠加算法的三维天然裂缝生成方法可利用分形参数控制裂缝网络的整体分布, 与人工裂缝耦合, 可表征页岩压后复杂的跨尺度裂缝系统。

嵌入式离散裂缝模型具有裂缝网格数量少、运算耗时短的优点, 能够灵活表征天然裂缝及人工裂缝属性, 在有效降低计算量的同时能够精确地模拟流体在基质-裂缝中的交换过程。

嵌入式离散裂缝模型与智能化历史拟合算法相结合能够降低裂缝、储集层等未知参数计算的不确定性, 实现关键储集层、裂缝参数的有效反演, 并实现气井产能量化预测。经实际应用验证, 一体化页岩气井产能评价模型预测结果具有较高的可信度。

符号注释：

ANNC— — 连接对之间的接触面积, m2; c— — 常数; C— — 净耗热量的相关常数, 无因次; ${{C}_{2}}\left( r \right)$— — 校正函数; cag— — 吸附气孔隙体积占比, 无因次; dNNC— — 连接对间距离, m; D— — 地层深度, m; Dc— — 二维分形维数; Dc3D— — 三维分形维数, 无因次; Dl3D— — 三维裂缝长度分形维数, 无因次; Dl— — 二维裂缝长度分形维数, 无因次; F— — 0~1的随机数, 无因次; i— — 拟合参数编号; j— — 拟合数据编号; jgw— — 流体编号（jgw=g代表气相, jgw=w代表水相）; K— — 基质渗透率, 10-3 μ m2; Ki— — 原始地层压力下的渗透率, 10-3 μ m2; KNNC— — 连接渗透率, 10-3 μ m2; ${{K}_{\text{r, }{{j}_{\text{gw}}}}}$— — 流体相对渗透率; Krw, Krg— — 水、气相对渗透率; Krwo, Krgo— — 水相, 气相相对渗透率端点值, 无因次; l— — 裂缝长度, m; lratio— — 区域与子区域边长比值, 无因次, 取值为2; lmin, lmax— — 最小、最大天然裂缝长度, m; L— — 模型尺度, m; n— — 吸附层数; $n\left( l, L \right)$— — 单位尺度内长度在l~l+dl裂缝条数, 条/m; N— — 裂缝总数量, 条; Nw, Ng— — 水相, 气相相对渗透率指数, 无因次; $N\left( L \right)$— — 长度大于lmin的裂缝总数量, 条; Niter— — 迭代次数; ${{N}_{\text{p}}}\left( r \right)$— — 两条裂缝中心距离小于r的裂缝对数, 对; p— — 孔隙压力, MPa; pi— — 原始地层压力, MPa; pL— — Langmuir气体压力, MPa; po— — 气相饱和压力, MPa; P— — 拟合参数量; P10, P50, P90— — 累计概率为10%、50%、90%时的值; Q— — 拟合数据量; qg— — 气体采出速度, m3/s; qw— — 水采出速度, m3/s; r— — 天然裂缝主长度, m; Sg, Sw— — 气相, 水相饱和度, %; Sgr— — 气相端点饱和度, %; Sseg— — 裂缝面面积, m2; Swirr— — 水相端点饱和度, %; t— — 时间, s; TNNC— — 非邻近网格连接对的传导率, 10-3 μ m2· m; u— — 流体流动速率, m/s; vag— — 单位质量岩石吸附的气体体积, m3/kg; vL— — Langmuir体积, m3/kg; vm— — 单位质量岩石最大吸附气体积, m3/kg; Vb— — 微元体积, m3; Vf— — 裂缝体积, m3; wf— — 裂缝宽度, m; wi, j— — 第j个参数对应的第i个数据点的权重值, 无因次; x, y, z— — 坐标轴, m; X— — 函数自变量; xmodel— — 模拟生产指标, 无因次; xhistory— — 实际生产指标, 无因次; α — — 二维裂缝密度常数, 无因次; α 3D— — 三维裂缝密度常数, 无因次; $\Gamma \left( X \right)$— — 伽马函数; γ — — 渗透率衰减常数, MPa-1; ${{\gamma }_{{{j}_{\text{gw}}}}}$— — 流体重度, kg/（m· s）2; μ — — 流体黏度, mPa· s; u— — 流体流速, m/s; ug, uw— — 气相、水相流速, m/s; ρ g, ρ ag, ρ s, ρ w— — 气, 吸附气, 岩石密度, 水, kg/m3; ϕ — — 孔隙度, %; ϕ f— — 裂缝网格孔隙度, %; ε — — 全局误差, 无因次。

（编辑 唐俊伟）

参考文献
[1] 王军磊, 贾爱林, 位云生, 等. 基于多井模型的压裂参数-开发井距系统优化[J]. 石油勘探与开发, 2019, 46(5): 915-921.
WANG Junlei, JIA Ailin, WEI Yunsheng, et al. Optimization workflow for stimulation-well spacing design in multiwall pad[J]. Petroleum Exploration and Development, 2019, 46(5): 915-921. [本文引用:1]
[2] KUCHUK F, BIRYUKOV D. Pressure-transient behavior of continuously and discretely fractured reservoirs[J]. SPE Reservoir Evaluation & Engineering, 2014, 17(1): 82-97. [本文引用:1]
[3] 林承焰, 李辉, 马存飞, 等. 致密砂岩储层天然裂缝建模方法[J]. 中国石油大学学报(自然科学版), 2019, 43(5): 21-33.
LIN Chengyan, LI Hui, MA Cunfei, et al. Modeling method of natural fractures in tight sand stone reservoirs[J]. Journal of China University of Petroleum (Natural Science Edition), 2019, 43(5): 21-33. [本文引用:2]
[4] 郎晓玲, 郭召杰. 基于DFN离散裂缝网格模型的裂缝性储层建模方法[J]. 北京大学学报(自然科学版), 2013, 49(6): 964-972.
LANG Xiaoling, GUO Zhaojie. Fractured reservoir modeling method based on discrete fracture network model[J]. Acta Scientiarum Naturalium Universitatis Pekinensisi (Natural Science Edition), 2013, 49(6): 964-972. [本文引用:1]
[5] AJISAFE F, THACHAPARAMBIL M, LEE D, et al. Calibrated complex fracture modeling using constructed discrete fracture network from seismic data in the Avalon Shale, New Mexico[R]. SPE 179130-MS, 2016. [本文引用:1]
[6] DREUZY J, DAVY P, ERHEL J, et al. Anomalous diffusion exponents in continuous two-dimensional multifractal media[J]. Physical Review E Statistical Nonlinear & Soft Matter Physics, 2004, 70(1): 016306. [本文引用:4]
[7] 李玮, 赵欢, 李思琪, 等. 页岩裂缝网络的几何特征二维表征及连通性分析[J]. 石油钻探技术, 2017, 45(6): 70-76.
LI Wei, ZHAO Huan, LI Siqi, et al. 2D characterization of geometric features and connectivity of fracture networks in shale formations[J]. Petroleum Drilling Techniques, 2017, 45(6): 70-76. [本文引用:2]
[8] LIU R C, YU L Y, JIANG Y J, et al. Recent developments on relationships between the equivalent permeability and fractal dimension of two-dimensional rock fracture networks[J]. Journal of Natural Gas Science and Engineering, 2017, 45: 771-785. [本文引用:2]
[9] KIM T H, SCHECHETER D S. Estimation of fracture porosity of naturally fractured reservoirs with no matrix porosity using fractal discrete fracture networks[J]. SPE Reservoir Evaluation & Engineering, 2009, 12(2): 232-242. [本文引用:6]
[10] 吴明录, 丁明才. 分形离散裂缝页岩气藏多级压裂水平井数值试井模型[J]. 中国石油大学学报(自然科学版), 2020, 44(3): 98-104.
WU Minglu, DING Mingcai. A numerical well testing interpretation model for multiple fractured horizontal well in fractured shale gas reservoir based on fractal discrete fracture network[J]. Journal of China University of Petroleum (Natural Science Edition), 2020, 44(3): 98-104. [本文引用:1]
[11] SUN J L, SCHECHTER D. Pressure-transient characteristics of fractured horizontal wells in unconventional shale reservoirs with construction of data-constrained discrete-fracture network[J]. SPE Production & Operations, 2018, 33(1): 21-31. [本文引用:3]
[12] CIPOLLA C L, LOLON E P, ERDLE J C, et al. Reservoir modeling in shale-gas reservoirs[R]. SPE 125530, 2009. [本文引用:1]
[13] MIRZAEI M, CIPOLLA C L. A workflow for modeling and simulation of hydraulic fractures in unconventional gas reservoirs[R]. SPE 153022, 2012. [本文引用:1]
[14] KUCHUK F, BIRYUKOV D. Pressure-transient behavior of continuously and discretely fractured reservoirs[J]. SPE Reservoir Evaluation & Engineering, 2014, 1: 82-97. [本文引用:1]
[15] KARIMI-FARD M, DURLOFSKY L J. A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features[J]. Advances in Water Resources, 2016, 96: 354-372. [本文引用:1]
[16] 苏浩, 雷征东, 李俊超, 等. 储集层多尺度裂缝高效数值模拟模型[J]. 石油学报, 2019, 40(5): 587-594.
SU Hao, LEI Zhengdong, LI Junchao, et al. An effective numerical simulation model of multi-scale fractures in reservoir[J]. Acta Petrolei Sinica, 2019, 40(5): 587-594. [本文引用:1]
[17] LI L, LEE S H. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media[J]. SPE Reservoir Evaluation & Engineering, 2008, 11(4): 750-758. [本文引用:1]
[18] YU W, XU Y F, WEIJERMARS R, et al. A numerical model for simulating pressure response of well interference and well performance in tight oil reservoirs with complex-fracture geometries using the fast embedded-discrete-fracture-model method[J]. SPE Reservoir Evaluation & Engineering, 2017, 21(2): 489-502. [本文引用:1]
[19] 朱大伟, 胡永乐, 崔明月, 等. 局部网格加密与嵌入式离散裂缝模型耦合预测压裂改造井产能[J]. 石油勘探与开发, 2020, 47(2): 341-348.
ZHU Dawei, HU Yongle, CUI Mingyue, et al. Productivity simulation of hydraulically fractured wells based on hybrid local grid refinement and embedded discrete fracture model[J]. Petroleum Exploration and Development, 2020, 47(2): 341-348. [本文引用:2]
[20] 严侠, 黄朝琴, 姚军, 等. 基于模拟有限差分的嵌入式离散裂缝数学模型[J]. 中国科学: 技术科学, 2014, 44(12): 1333-1342.
YAN Xia, HUANG Zhaoqin, YAO Jun, et al. The embedded discrete fracture model based on mimetic finite difference method[J]. Scientia Sinica Technologica, 2014, 44(12): 1333-1342. [本文引用:1]
[21] 薛亮, 吴雨娟, 刘倩君, 等. 裂缝性油气藏数值模拟与自动历史拟合研究进展[J]. 石油科学通报, 2019, 4(4): 335-346.
XUE Liang, WU Yujuan, LIU Qianjun, et al. Advances in numerical simulation and automatic history matching of fractured reservoirs[J]. Petroleum Science Bulletin, 2019, 4(4): 335-346. [本文引用:1]
[22] DARCEL C, BOUR O, DAVY P, et al. Connectivity properties of two-dimensional fracture networks with stochastic fractal correlation[J]. Water Resources Research, 2003, 39(10): 1272. [本文引用:3]
[23] SCHERTZER D, LOVEJOY S. Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes[J]. Journal of Geophysical Research, 1987, 92(D8): 96-93. [本文引用:1]
[24] XU Y F, YU W, SEPEHRNOORI K. Modeling dynamic behaviors of complex fractures in conventional reservoir simulators[J]. SPE Reservoir Evaluation & Engineering, 2019, 22(3): 1110-1130. [本文引用:1]
[25] TRIPOPPOOM S, YU W, HUANG H Y, et al. A practical and efficient iterative history matching workflow for shale gas well coupling multiple objective functions, multiple proxy-based MCMC and EDFM[J]. Journal of Petroleum Science and Engineering, 2019, 176: 594-611. [本文引用:2]
[26] WEI Y, WANG J, JIA A, et al. Optimization of managed drawdown for a well with stress-sensitive conductivity fractures: Workflow and case study[J]. Journal of Energy Resources Technology, 2021, 143(8): 915-921. [本文引用:1]