基于MCMC 的页岩气藏压裂缝网不确定性反演法

贾品，封皓然，王苏冉，王志凯，牛烺昱，程林松

1 中国石油大学(北京)油气资源与工程全国重点实验室，北京 102249

2 中国石油大学(北京)石油工程学院，北京 102249

3 中海油研究总院有限责任公司，北京 100028

0 引言

我国页岩气资源分布广泛且储量丰富，可采资源量达15×1012～30×1012m3，近些年先后设立了涪陵、长宁—威远和昭通等国家级页岩气示范区，3500 m以上浅储层地质压裂工程一体化技术基本成型，实现了我国浅层页岩气的工业化生产；3500 m以下深储层页岩气也得到有效的尝试性开发，发展潜力巨大。不论对于浅层还是深层页岩气，由于页岩基质以微纳米级孔隙为主，储层渗透率极低、物性更差，能够形成三维立体复杂裂缝网络的水平井体积压裂技术成为了页岩气经济有效开发的主要手段。其中，准确反演压裂缝网参数是压裂效果评价、压裂优化设计及产能设计的关键问题，也一直是页岩气藏压裂及开发工程方面的热点问题。

页岩气藏压裂缝网的反演通常是基于一个适合于目标区的渗流数学模型，通过不断调整缝网参数对实际生产动态数据或其相关变形量进行拟合从而实现缝网关键参数的推算。目前常见的线性流分析法、图版拟合法及生产历史拟合法均属于此类方法的范畴。线性流分析法通过绘制拟压力规整化产量和物质平衡拟时间的双对数曲线，找出曲线上的线性流段，拟合线性流段的斜率和截距获得缝网参数。国内外大多数研 究 都 采 用Wattenbarger[1]、Bello和Wattenbarger[2]、Al-Ahmadi和Wattenbarger[3]、Brown[4]、Stalgorova和Mattar[5]等提出的线性流模型。Nobakht和Matter[6]、Clarkson等[7]、Qanbari和Clarkson[8]、Xu等[9]通 过 确定的原始地层压力，储层厚度，压裂段间距(通常认为段间全改造)，基质及裂缝的压缩系数，基质及裂缝的孔隙度，及气体吸附解析参数等参数对压裂缝长度、导流能力及改造区大小、渗透率进行反演。图版拟合法是基于缝网流动模型绘制裂缝关键参数影响下的产量敏感性图版，运用理论图版拟合实际生产，从而反演缝网参数。由于简化的物理模型，线性流分析法、图版拟合法能够反演的参数有限。Moghadam[10]等和Abdulal等[11]应用图版拟合法反演了页岩气井储层参数。生产历史拟合法多用于复杂缝网形态的近似表征，Li和Zhu[12]研究了水平井压裂过程的温度不稳态响应，Zhang和Zhu[13]进一步提出了裂缝长度和导流能力反演方法。结合微地震和岩心数据，Cipolla等[14]、Weng等[15]提出了一种利用微地震事件位置校准地质力学模型并生成缝网的方法，Gamboa等[16]和Sun等[17]进一步通过考虑页岩气水排采动态建立了一种压裂缝网展布形态解释方法。基于生产历史拟合的缝网参数反演存在历史拟合数据量大、耗时耗力及适用性差等不足。

确定性反演需要相对准确的输入参数，才能得到合理的反演值。由于渗流数学模型本身包含了众多因素，例如储层参数、流体性质及非线性渗流机理，这都给页岩气压裂缝网反演带来了极大的不确定性。因此，页岩气藏压裂缝网反演应该不是一个确定性过程，而是由于基础参数的不确定性引起的一个不确定过程。上述的3 类方法均是确定性反演方法，即给定了基础参数后，便反演出了一套缝网参数，人们还无法明确缝网反演参数的多组合性及其概率分布，难以客观并全面评估页岩气压裂改造效果。近几年，基于马尔可夫—蒙特卡洛方法(MCMC)的不确定反演方法已经广泛应用于岩石物理领域[18-19]，但在页岩气藏领域还仅停留在基于MCMC的不确定性产能评价，针对压裂缝网参数不确定性的反演研究鲜少，亟需提出一种适用于页岩气藏压裂缝网参数的不确定性反演方法。

本文针对页岩气藏压裂缝网参数在历史拟合过程中的不确定性，基于马尔可夫—蒙特卡洛方法(MCMC)的原理，提出一种适用于页岩气藏的压裂缝网反演方法。通过随机生成反演所需的基础参数样本集，构建缝网参数解释模型，联合线性流分析法和图版拟合法，获取缝网反演参数的多组合性及其概率分布，评价由于基础参数的不确定性对压裂缝网参数解释结果的影响程度，对页岩气藏压裂效果评价及合理开发对策制定提供较全面且实际的指导。

1 压裂缝网确定性反演

如上所述，页岩气藏压裂缝网的不确定性反演要基于压裂缝网的确定性反演，即通过随机抽样方法产生若干基础参数的样本集，针对每个样本集进行对应的压裂缝网确定性反演，然后收集若干基础参数样本下的反演结果，从而得到压裂缝网参数的组合及概率分布。以下首先给出页岩气藏压裂缝网参数的确定性反演方法，然后阐述压裂缝网反演所需基础参数的不确定性，从而建立不确定性反演方法，最后进行实例分析。

1.1 页岩气藏压裂水平井物理模型

如图1 所示，页岩气藏一口水平井经体积压裂后，形成了复杂的裂缝网络系统，即压裂改造体(SRV)。整个储层可等效为由人工裂缝、压裂改造区(内区)和未改造区(外区)三个区域组成的系统。人工裂缝考虑为等长双翼缝；内区考虑压裂改造形成的复杂缝网系统，将其处理成高渗透区域；外区由于几乎没有受到压裂改造，将其处理成低渗透区域。其他基本假设为：

图1 页岩气藏压裂缝网物理模型示意图Fig.1 Schematic diagram of physical model of fracture network in shale gas reservoir

(1)储层为水平等厚，水平井位于储层中心，裂缝贯穿整个储层；

(2)考虑气体及岩石的压缩性，不考虑重力作用；

(3)各区域内气体流动都视为线性流，即外区系统流体线性流入内区系统，再由内区系统线性流入人工裂缝，最后由人工裂缝线性流入水平井筒(如图1 所示)。

1.2 压裂缝网确定性反演方法

页岩气主要以游离态、吸附态和溶解态3 种形式赋存于页岩储层中，其中游离态和吸附态占比最高。Al-Ahmadi和Wattenbarger[20]针对考虑页岩气吸附解析及非线性流动机理的三线性流模型做了大量工作，本文在此基础上构建基于产量递减分析的压裂缝网参数反演方法。对于图1 所示的物理模型，其在定压生产条件下拉普拉斯空间的无因次产量解为[21]：

式中，f(s)定义为：

其中，qD为无因次产量，yFD为无因次裂缝半长，yeD为无因次外区宽度，λIF为内区与人工裂缝窜流系数比，λIo为内区与外区窜流系数比，ηD为无因次导压系数，ω为弹性储容系数，s为拉普拉斯变量。以上所有变量均为无因次，其具体定义可参考笔者文章[22]。

大量矿场页岩气压裂井开发动态表明生产早期页岩气流动通常会出现较为明显的线性流特征，生产中后期流态以边界控制流为主[23]，直线分析方法是分析线性流及边界控制流的重要技术手段。直线分析法基于线性流或边界控制流阶段压力及产量与时间的局部解，通过绘制特定变量的坐标图，利用斜率或截距推算压裂缝网参数，如表1 所示。本文首先采用直线分析方法定量分析页岩气产量动态数据，利用线性流获取人工裂缝半长与改造区渗透率关系利用边界控制流获得改造体积大小，再以这些值为初值或约束条件利用图版拟合法反演得到人工裂缝半长、人工裂缝导流能力、改造区渗透率及改造区宽度参数。

表1 直线分析法变量及反演参数表Table 1 Table of variables and inversion parameters of linear analysis method

实际矿场生产中，页岩气压裂井往往是以变产量或变井底流压生产，生产动态分析需引入物质平衡拟时间和拟压力规整化产量分析变工作制度下生产动态。其中，物质平衡拟时间项为：

其中，Wg是累积产气量；Vpm是改造区系统的总孔隙体积，m3；Smgi为初始含气饱和度；Bgi是原始地层压力下的气体体积系数，Bg是平均压力下的气体体积系数，是压力的函数;pi为原始地层压力，MPa。平均压力的求解可以使用牛顿迭代法，是改造区系统的总孔隙体积可以表示为：

其中h为储层厚度，m；xf为人工裂缝半长，m；这里人工假设半缝长初值用于计算物质平衡拟时间，与解释出的结果对比迭代。Ls为改造区宽度，m；通过统计水平井长度和压裂段数可折算得到。

拟压力归整化产量项为：

式中，ψi为原始地层压力下对应的拟压力，MPa2/mPa·s；ψw为井底流压对应的拟压力，MPa2/mPa·s。以一口体积压裂水平井数值模拟模型验证上述线性流分析法和图版拟合法联合解释压裂缝网参数的流程，模型示意如图2 所示，水平段长度为1500 m，压裂25 段，段间完全改造。人工裂缝渗透率为5000 mD，人工裂缝半长为102.5 m，改造区渗透率为5×10-4mD，改造区宽度为30 m，压裂液仅赋存于人工主缝中，基质内不含水。

图2 数值模型示意图Fig.2 Schematic diagram of numerical model

基于线性流分析法和图版法联合解释压裂缝网参数，具体解释流程为：(1)通过实际产气量及井底流压动态数据计算物质平衡拟时间和归一化拟压力，绘制两者的双对图版曲线并识别流动阶段，如图3 所示。图3 中可以明显看出线性流，边界控制流两个流动阶段，在边界流之前，由于数据点时间间隔过短，没有明显的双线性流阶段，通过流态的划分可以进一步验证三线性流的流动阶段划分合理；(2)提取双对数图版线性流阶段曲线特征，拟合直线段斜率mL(图4)，利用表1 中的关系式确定人工裂缝半长与改造区渗透率关系，进一步提取双对数图版中的边界控制流阶段曲线特征，拟合直线与x轴的截距(图5)，利用表1中的关系式确定改造体大小；(3)以线性流拟合关系式及边界控制流得到的改造体大小为约束，不断调整人工裂缝半长、人工裂缝渗透率、改造区渗透率及改造区宽度参数，拟合全生产历史动态下的物质平衡拟时间和规整化产量图版，反演缝网参数，如图6 所示。反演参数见表2，可以看出压裂缝网的4 个关键参数反演值与数值模型值几乎一致，吻合度均在94.4%以内，验证了确定性反演方法的准确性和可靠性。

表2 数值模型参数反演参数表Table 2 Numerical model parameter inversion parameter table

图3 生产不同阶段流态划分Fig.3 Flow regime division in different production stages

图4 线性流直线分析图Fig.4 Linear flow line analysis diagram

图5 边界控制流直线分析图Fig.5 Boundary -dominated flow analysis diagram

图6 产量递减图版拟合结果Fig.6 Production decline type curve fitting

2 基础参数不确定性分析

压裂缝网反演过程中所需基础参数的不确定性影响着线性流分析和图版拟合的参数解释结果。页岩气藏压裂缝网参数反演的不确定性主要来源于以下两方面：一方面由于输入参数多且适用范围有限，如天然岩心及测井资料获取的储层参数本质上仅代表了近井筒附近储层物性，而无法准确反映储层整体性质；另一方面复杂缝网多参数反演易出现局部最优解，图版拟合的多解性强。以下给出了页岩气藏压裂缝网反演中的不确定性参数分析。

页岩储层参数不确定性：页岩储层参数主要包括储层有效厚度、原始地层压力、原始地层温度、含气饱和度、基质孔隙度。页岩纵向不同埋深及平面井间参数存在差异，基于目标区储层物性评价资料可确定该类参数的分布及范围。

页岩吸附解析及流动参数不确定性：页岩不同井区的干酪根及初始含气量存在差异，且含气量随着储层埋深增加而增加。另外，由于不同区域储层页岩基质孔隙直径变化大，气体流动参数在不同位置处也不同。因此，页岩吸附解析及流动参数不确定性分别体现在兰氏体积、兰氏压力和扩散系数中。

页岩压裂缝网基础参数不确定性：页岩储层经过大型水力压裂后，人工裂缝半长、人工裂缝渗透率、改造区渗透率及改造区宽度参数需要反演，但裂缝压缩系数、改造区孔隙度、改造区压缩系数需要作为反演的基础参数提前提供，该类参数具有不确定性。通过统计目标区块的各井的原始地层压力，储层有效厚度，朗格缪尔体积，朗格缪尔压力，水平井长度，基质孔隙度，基质压缩系数，裂缝压缩系数，选择合适的概率密度函数(表3)进行拟合，得到针对该区块的不同储层参数的概率密度函数。

表3 常用概率密度函数Table 3 Common probability density function

3 基于MCMC 的压裂缝网不确定反演

本文基于马尔可夫-蒙特卡洛方法(MCMC)对页岩气藏水平井压裂缝网参数进行不确定性反演。基于页岩气藏压裂水平井基础参数，采用MCMC抽样方法获取储层参数、吸附解析及流动参数、压裂缝网基础参数的概率分布，结合线性流分析法和图版拟合法识别典型流动阶段，拟合生产动态，获取各反演参数的概率分布。MCMC抽样方法的基本思想是采用Metropolis抽样方法，构造一个遍历的马尔可夫链，从任意概率向量出发，通过多次迭代计算，获得不变分布，使其成为所需要的抽样分布。抽样次数越大，抽样分布越趋于一个稳定分布。MCMC抽样基于贝叶斯理论，其表达式为：

① 生成基础参数：设置储层、流动及裂缝基础参数的初始样本θs，通常取各参数分布的中值；

② 基础参数分布抽样：根据区块基础参数分布概率(统计获得)，产生新样本θp；

③ 计算新样本θp接收概率p：接受概率p的计算方法：

④ 判断是否更新样本：产生一个(0，1)随机数r，若r＜p，接受新样本θs=θp，反之θs=θs-1；

⑤ 判断是否结束：若达到最大迭代步数，则抽样结果统计，否则返回步骤②。具体的不确定性反演方法流程如图7 所示。

图7 页岩气藏压裂缝网参数不确定性反演流程图Fig.7 Flow chart of parameter uncertainty inversion of pressure fracture network in shale gas reservoir

4 实例应用

4.1 基础参数不确定性分析

以我国西南地区某页岩气藏一口体积压裂水平井为例进行分析，该气藏埋深约3500 m，储层有效厚度约30 m，储层中深压力为67.3 MPa，地层温度375 K，初始含水饱和度约为0.45。岩心分析结果表明，页岩基质孔隙度约为0.04，渗透率跨度较大(10-5mD～10-3mD)。该水平井分为26 段压裂，每段压裂3 簇，压裂总段长1430 m，入地总液量53831 m3，入地总砂量2442.4 t。该页岩气藏压裂水平井的产气动态及井底流压如图8 所示。根据该地区临井基础资料分析，获取了储层参数、流动参数及压裂缝网基础参数的范围及概率分布，如表4 所示。

表4 各相关参数取值Table 4 Values of relevant parameters

图8 生产动态数据图Fig.8 Production dynamic data fitting

4.2 压裂缝网参数解释

基于表4 中基础参数的概率分布，利用MCMC随机生成至少1000 个基础参数样本集合，再利用不同的基础参数样本结合等式(3)和(4)计算物质平衡拟时间及归一化拟压力，绘制出1000 个页岩气流动特征双对数曲线。为了便于展示，图9(a)给出了其中基于6 个基础参数样本的曲线，可以看出生产早期气体流动表现出斜率为1/2 的线性流特征，表明改造区得到有效动用，改造区内气体线性流入人工裂缝中；随后出现斜率为1 的边界控制流，改造区压力不断衰竭。如图9(b)和9(c)所示通过拟合线性流阶段和边界控制流阶段的直线段，分别可获得人工裂缝半长与改造区渗透率关系以及改造体大小，进一步通过拟合物质平衡拟时间和归一化产量图版，如图9(d)所示，进而反演获得压裂缝网在每个基础参数样本下的参数。

图9 实际井压裂缝网参数不确定性反演结果Fig.9 Uncertainty inversion results of actual well pressure fracture network parameters.(a) Flow regime division in different production stages (b) Linear flow line analysis diagram (c) Boundary-dominated flow analysis diagram (d) Production decline type curve fitting

4.3 压裂缝网参数不确定性分析

基于图9 所展示的压裂缝网反演流程，获得了人工裂缝半长、人工裂缝导流能力、改造区渗透率及改造区宽度参数组合下的不确定分布，其概率分布以及累积概率分布如图10～13 所示。同时也给出了各反演参数的拟合值及P10、P50、P90 下的反演结果，见表5。

表5 W1 井裂缝参数不确定性反演小结Table 5 Summary of uncertainty inversion of fracture parameters in well W1

图10 人工裂缝长度分布图Fig.10 Distribution of hydraulic fracture length.(a) Uncertainty inversion probability distribution diagram (b) Cumulative probability distribution

图11 人工裂缝导流能力不确定性分析结果Fig.11 Uncertainty analysis results of artificial fracture conductivity.(a) Uncertainty inversion probability distribution diagram (b) Cumulative probability distribution

图12 改造区渗透率不确定性分析结果Fig.12 Analysis results of permeability uncertainty in the reconstruction area.(a) Uncertainty inversion probability distribution diagram (b) Cumulative probability distribution

图13 改造区宽度不确定性分析结果Fig.13 Uncertainty analysis results of width of reconstruction area.(a) Uncertainty inversion probability distribution diagram(b) Cumulative probability distribution

从人工裂缝半长不确定性反演分析图可以看出，反演得到的人工裂缝半长基本符合正态分布且跨度较大，同时累计频率分布P10、P50、P90 对应的参数数值大小较为分散，表明其不确定性较强。

从人工裂缝导流能力不确定性反演分析图可以看出，反演得到的人工裂缝导流能力基本符合正态分布且跨度较大，同时累计频率分布P10、P50、P90 对应的参数数值大小较为分散，表明其不确定性较强。

从改造区渗透率不确定性反演分析图中可以看出，反演得到的改造区渗透率也基本符合正态分布且跨度较大但整体范围较小，同时累计频率分布P10、P50、P90 对应的参数数值大小较为分散，表明其不确定性一般。

由改造区宽度不确定性反演分析图中可以看出，反演得到的改造区宽度的概率分布较为集中且主要集中在45 m左右，且累计频率分布P10、P50、P90 对应的参数数值大小也较为相近，表明其不确定性较弱。

综合分析基于随机生成的基础参数样本反演得到的缝网参数不确定性可得：该地区页岩气藏体积压裂后人工裂缝半长和人工裂缝导流能力的不确定性较强，改造区渗透率不确定性一般，改造区宽度不确定性较弱。通过对压裂裂缝参数不确定性的分析，即可为实际区块中压裂裂缝参数多解性提供参考。

5 结论

本文基于页岩气单井产能模型，首先构建了压裂缝网参数确定性反演方法，然后结合马尔可夫—蒙特卡洛方法(MCMC)形成了压裂缝网参数不确定性反演方法，最后对我国西南地区一口页岩气井压裂缝网的4 个关键参数进行了不确定性反演，得到了不同基础参数下的压裂裂缝反演参数的概率分布及累积概率分布，得到结论如下：

(1)构建了页岩气压裂井确定参数下的储层参数反演方法，结合直线分析法和图版拟合法对人工裂缝半长、人工裂缝渗透率、改造区渗透率和改造区宽度进行确定性解释，并结合数值模型论证了方法的可靠性；

(2)页岩储层缝网参数反演参数多、不确定性大，基于马尔可夫-蒙特卡洛方法(MCMC)的页岩气藏压裂缝网不确定性反演方法主要包括蒙特卡洛抽样，典型流动阶段分析及生产动态拟合等步骤，能快速获取人工裂缝长度、导流能力、改造区渗透率及改造区裂缝宽度等关键缝网参数及其概率分布；

(3)针对我国西南地区某页岩气区块一口体积压裂水平井生产动态，对其进行不确定性反演，明确了压裂缝网四个关键参数的不确定性排序，人工裂缝半长的拟合值为75 m，累计概率分布P10 为40.5 m，P50为69.3 m，P90 为104.3 m。人工裂缝导流能力的拟合值为125 mD，累计概率分布P10 为84.1 mD，P50 为166.9 mD，P90 为270.3 mD。改造区渗透率的拟合值为5.0×10-4mD，累计概率分布P10 为4.4×10-4mD，P50 为5.8×10-4mD，P90 为7.0×10-4mD。改造区宽度的拟合值为45.3 m，累计概率分布P10 为40.2 m，P50 为45.3 m，P90 为47 m。