【论文】任文希,等：物理—数据协同驱动的页岩气井产量预测方法

页岩气开发过程中，气井产量的准确预测对于气井生产能力评价、配产优化、可采储量评估等有着重要的意义，并直接影响最终的经济效益分析^[1]。然而，页岩气井产量预测是一个高度非线性的、强耦合的问题，具有复杂性高和时变性强的特点，准确预测难度大。一方面，页岩气一般采用水平井分段多簇压裂技术开发，压裂后会形成人工裂缝和天然裂缝交织的复杂裂缝网络^[2-4]。另一方面，页岩气储层普遍发育微纳米级孔隙，气体在其中的赋存和渗流机理复杂，如吸附气和自由气共存^[5-8]、黏性流和滑脱流耦合等^[9-12]，对传统的渗流理论提出了巨大的挑战。

现阶段的页岩气井产量预测方法大致可以分为4 类：经验产量递减分析、解析/ 半解析解、数值模拟和人工智能方法。经验产量递减分析方法通常基于经验或半经验方程，通过拟合历史产量数据来确定方程参数，然后通过外推的方式预测未来产量。这种方法只需要历史产量数据，计算简单，使用方便，但一般只适用于某一流动阶段，如Arps 方法^[13] 只适用于边界流阶段，Duong 方法^[14] 则只适用于线性流阶段。部分研究者为了拓宽经验产量递减分析方法的适用范围，提出了分段递减方法，即组合式方法^[15-16]，但流动阶段的转变点难以确定。此外，经验产量递减分析方法还假设生产制度恒定， 这与页岩气井的实际生产制度不符。解析/ 半解析解方法构成了页岩气井不稳定流分析（Rate-Transient Analysis）的基石，其最大特点就是计算效率较高， 但也存在一些不足：①模型的多解性较强^[17] ；②裂缝过于简化，如平板裂缝^[18] ；③综合压缩系数等非线性项的线性化处理^[19]。数值模拟方法考虑因素全面，还可以直观地展示地下天然气的流动过程，但是复杂裂缝网络难以精确识别和表征，因此建模难度大。此外，数值模拟方法计算工作量大、耗时长。

近年来，随着人工智能的快速发展，一些研究者开始应用时间序列分析方法来分析页岩气井产量变化规律，他们大部分都选择了长短记忆神经网络（LSTM）作为产能模型^[20-22]。但是，LSTM 是纯数据驱动方法，可解释性差。此外，LSTM 采用循环策略开展产量预测，即需要利用之前的预测值来开展新的预测，因此存在误差累积的风险，可能导致预测结果失真，难以准确预测长期产量^[23-24]。

针对现有页岩气井产量预测方法的不足，笔者引入机理—数据融合建模的思路，结合连续拟稳态假设、物质平衡方程、产量递减分析方法和递推原理， 建立了物理—数据协同驱动的页岩气井产量预测方法，突破了复杂裂缝网络精确识别和表征的难题，克服了纯数据驱动方法可解释性差、预测误差易累积的不足。

1 　物理—数据协同驱动原理及实现方式

物理—数据协同驱动综合了物理和数据二者的优势，一方面可以通过数据驱动降低建模难度、简化建模过程，另一方面可以通过物理驱动提高模型泛化能力，增强模型可解释性。目前大致可以通过3 种方法实现物理—数据协同驱动：①基于输入数据（又称基于输入数据的嵌入机制）。该方法利用物理模型生成额外的特征，并与实际数据一并作为数据驱动模型的输入，或者结合领域知识构建数据集^[25]。②基于模型架构。该方法通过修改数据驱动模型的中间层结构来添加物理信息^[26]。③基于损失函数。该方法通过构建包含物理信息的损失函数来训练数据驱动模型，如Raissi 等^[27] 提出的物理信息神经网络（Physics Informed Neural Network）。

这3 类方法中，基于输入数据的方法较为灵活且有效^[26]，笔者也选用了该方法开展页岩气井产量预测研究，具体的建模思路如下（图1）：

1）引入连续拟稳态假设和物质平衡方程，简化对储层非均质性的描述，建立物理驱动产能模型， 然后通过历史拟合确定产能模型中的各项拟合参数， 进一步获得动态泄流区平均压力与累计产量的经验关系式。

2）利用拟合得到的模型参数，基于建立的产能模型计算动态泄流区平均压力和物质平衡拟时间，并将其作为两个额外的特征。

3）利用第2）步生成的两个额外特征来指导数据驱动模型训练，增强数据驱动模型的可解释性。具体步骤为：①结合动态泄流区平均压力、物质平衡拟时间和历史生产数据，计算不同时刻的采气指数；②开展基于数据驱动的递减规律和特征直线分析，进一步构建采气指数与物质平衡拟时间的 关系。

4）基于第1）步获得的动态泄流区平均压力与累计产量的经验关系式和第3）步获得的采气指数与物质平衡拟时间的关系，利用递推原理开展产量预测。

2 　物理驱动产能模型

页岩气储层渗透率极低，压力波传播速度慢， 页岩气井长期处于不稳定流动状态。连续拟稳态法可以将不稳定流动问题简化为拟稳态流动问题，进而结合探测半径/ 距离方程和物质平衡方程来研究油气井的生产动态。本文亦采用连续拟稳态法，将不稳定流动近似为一系列拟稳态流动。Clarkson 等^[28] 和Yuan 等^[29-31] 提出了动态泄流区的概念，并将不同时刻的压降波及区视为动态泄流区。在此基础上，他们应用连续拟稳态法开展了针对非常规油气井的产能分析和预测。相对于传统的产能分析和预测方法， 基于连续拟稳态法的产能分析和预测方法具有以下3 大优势^[28-30] ：①对储层非均质性的描述得到了极大简化，无需预先知道储层物性参数的空间分布和裂缝网络的形状；②求解在时域中就可以完成，无需转换到频域，省去了逆变换过程；③在计算量和计算复杂度可控的条件下，考虑了真实的PVT 参数和岩石压缩性，且适用于变流压—变产量的情况。Clarkson 等^[28] 和Yuan 等^[29] 假设泄流区为矩形，裂缝为平板裂缝，建立了探测半径/ 距离方程。在此基础上，为了避免预设泄流区和裂缝形状，本文建立了压降波及区的孔隙体积计算模型，将求解探测半径/ 距离转变为求解孔隙体积，如图2 所示。物理驱动产能模型的三大核心组成如下：基于连续拟稳态法的动态泄流区孔隙体积计算、动态泄流区物质平衡方程和以累计产量误差最小为目标的历史拟合策略。

2.1 　基于连续拟稳态法的动态泄流区孔隙体积计算

井控范围通常指单井对烃类流体流动的控制范围；缝控范围是井控范围的子集，对应人工裂缝可动用的储量；动态泄流区则是缝控范围的子集，随着生产时间的推移而变化。随着生产的进行，压降波及区不断扩大，动态泄流区也不断扩大，当生产时间趋于无限时，对应的动态泄流区就是缝控范围。

动态泄流区孔隙体积计算方程的具体推导过程如下。对于任意时刻，拟稳态流动状态下，泄流区内各点的压力以相同的速度下降，此时对应的产量（q_g）为^[32]：

综合压缩系数同时考虑了孔隙体积、气体体积和吸附量的变化。采用Langmuir 方程来描述吸附量随压力的变化，因此综合压缩系数可以表示为：

TOC 可以从测井资料中获取。进一步，为了弱化产量、井底流压和甲烷物性变化带来的影响，引入拟压力（ψ）和物质平衡拟时间（t_a）。拟压力同样具有压力的量纲，其定义为：

不同时刻的v_p 主要与对应时刻的产量、物质平衡拟时间和平均压力有关。为了获得不同时刻动态泄流区的平均压力，需要求解对应时刻的物质平衡方程。

2.2 　动态泄流区物质平衡方程

页岩气藏的含水饱和度较低，因此笔者忽略了水相的影响，并假定页岩气藏为干气。弹性开采条件下，储层压力降低后，产出的气体源于气体膨胀、孔隙压缩和气体解吸，这3 项分别对应自由气含量的变化、吸附气含量的变化和孔隙体积的变化。因此， 基于物质平衡原理，累计产量（G_p）可以表示为：

式（12）右边的3 项分别对应气体体积变化、吸附气含量变化和孔隙体积变化。利用历史生产数据（已知q_g 和G_p），结合式（11）、（12）可以通过迭代法求解任意时刻对应的动态泄流区平均压力。但这种计算方法无法实现外推，不能用于预测未来产量。假定动态泄流区平均压力随累计产量呈指数递减，建立一个动态泄流区平均压力与累计产量的经验关系式：

2.3 　以累计产量误差最小为目标的历史拟合策略

为了确定动态泄流区平均压力与累计产量经验关系式中的参数a₁，须开展历史拟合。常规历史拟合一般以日产量误差最小为目标，但页岩气井产量动态变化复杂，日产量拟合难度高。为克服这一困难， 以累计产量误差最小为目标开展历史拟合。这是由于相对于日产量数据，累计产量数据呈递增的单向变化趋势，不存在波动且更加光滑，有利于降低拟合难度。

具体的历史拟合策略如下：①将关井阶段的日产量数据删除；②将日产量数据进行累加，从而获得单调递增的累计产量数据；③结合最小二乘法的思想，将历史拟合的目标函数定义为累计产量残差的平方和，因此目标函数obj 可以表示为：

本文建立的产能模型采用拟压力代替压力，采用物质平衡拟时间代替时间，弱化了产量、流压和甲烷物性变化带来的影响。此外，以累计产量误差最小为目标开展历史拟合，弱化了生产制度变化带来的影响。因此，本文建立的产能模型可以自动适应流压—产量变化。

动态泄流区孔隙体积计算方程[ 式（11）]、物质平衡方程[ 式（12）] 和平均压力计算方程[ 式（13）] 构成的耦合方程组一共包含3 项拟合参数，分别为气井投产前的储层压力（p_i），地层压缩系数（c_f）和参数（a₁）。其中，p_i 无法直接通过地层压力系数计算得到，这是因为页岩气井需要进行大规模水力压裂，而压裂后压裂液的返排率较低，滞留的大量压裂液具有增能作用，导致储层压力高于原始状态（水力压裂前）。p_i 的取值范围如下：下限为原始条件下（压裂前）的储层压力，该值可以通过地层压力系数估算；上限为近水平井根端压裂段压裂后停泵压力与井筒液柱压力之和。另外一个拟合参数c_f 也很难通过岩石力学实验或是其他室内实验手段直接获得，这是因为大规模水力压裂后页岩储层会形成复杂裂缝网络， 因此c_f 反映的是基质—裂缝系统的压缩性，页岩基质孔隙压缩系数c_ϕ 一般介于(10^－4 ～ 10^－3)/MPa，而水力裂缝的压缩系数介于(10^－2 ～ 10^－1)/MPa^[36]。因此， 将c_f 的下限设置为页岩的压缩系数，对应的数量级在(10^－6 ～ 10^－5)/MPa，上限设置为水力裂缝的压缩系数，对应的数量级在(10^－2 ～ 10^－1)/MPa。

采用经典的粒子群算法来确定这3 个拟合参数， 具体流程如下：①基于给定的p_i、c_f 和a₁ 的初始值生成粒子群，进一步结合历史生产时间，利用式（13）， 计算不同时刻对应的动态泄流区平均压力；②接下来结合历史生产数据q_g，基于式（11）计算不同时刻对应的v_p ；③进一步利用动态泄流区物质平衡方程，即式（12），计算不同时刻对应的G_p ；④结合实际的G_pr 和计算得到的G_p，通过预先设定的目标函数， 即式（14），计算粒子群适应值；⑤最后基于计算得到的粒子群适应值，调整搜索方向，直到满足收敛条件（图3）。计算涉及的其他参数及其确定方法如表1 所示。

历史拟合后物理驱动产能模型可以生成两个额外特征：动态泄流区平均压力和物质平衡拟时间，这两个额外特征后续会参与数据驱动方法的 训练。

3 　物理—数据驱动产量分析

为了预测气井未来产量，需要建立产量与井底流压之间的关系，利用采气指数来定量表征这两者的关系。采气指数通过物理驱动产能模型生成的额外特征—动态泄流区平均压力，结合产量数据计算。进一步，利用物理驱动产能模型生成的另外一个额外特征—物质平衡拟时间，结合采气指数，共同指导数据驱动方法训练，从而实现物理—数据协同驱动。在双对数坐标图中绘制了采气指数—物质平衡拟时间散点图（图4）。整体来看，在双对数坐标中， 采气指数随着物质平衡拟时间的增加逐渐递减，前期递减慢，后期递减快，后期数据点可以收敛为一条斜率接近－1 的直线，即特征直线。其他井对应的采气指数—物质平衡拟时间双对数图同样存在斜率接近－1 的直线，受篇幅限制，不再一一展示。此外， 相关理论研究表明^[37]，在生产后期，页岩气井规整化产量—物质平衡时间双对数图上会出现斜率为－1 的直线，而规整化产量与采气指数的倒数接近。因此， 可以采用线性方程来表征采气指数与物质平衡拟时间的关系，即采气指数递减方程：

特征直线的确定流程如下：①找到转折点，即采气指数的最大值（图4-a）；②利用转折点之后的数据，基于最小二乘法开展线性拟合，从而确定特征直线，并通过特征直线的截距和斜率得到a₂、b₂ 值。

应用本文方法的关键就是采气指数—物质平衡拟时间双对数图中的特征直线，若双对数图中出现了特征直线，则可以开展产量预测；反之，则不可以开展产量预测。实际应用过程中，随着生产时间的增加， 可以反复应用本文方法开展产量预测，即利用新增的实际生产数据来修正产量预测结果，进一步提高产量预测结果的精度。

4 　基于递推原理的产量预测

递推法是一种重要的数学方法，它是指从已知条件出发，依据某种关系，逐次推导出所要求的各中间结果及最后结果。基于式（17），任意时刻t_j 对应的采气指数可以表示为：

结合式（19）、（20）计算t_j+1 时刻对应的动态泄流区平均压力：

重复上述过程便可以预测未来的产气量和累计产量，相应的预测流程如图1-b 所示。对于未来产量预测，井底流压取现今井底流压（对应最后一天的历史生产数据）与大气压力（0.1 MPa）的算术平均值。

5 　现场生产数据测试和对比分析

基于国内某区块页岩气井现场生产数据，对建立的物理—数据协同驱动的页岩气井产量预测方法进行测试，并与经验产量递减分析和时间序列分析方法进行了对比。选择Logistic Growth Model（ LGM）、 Stretched-Exponential Production Decline（SEPD） 和Duong 方法作为经验产量递减分析方法的代表；选择非线性自回归神经网络（NAR）作为时间序列分析方法中浅层学习的代表；选择长短记忆神经网络（LSTM）作为时间序列分析方法中深度学习的代表。

5.1 　与经验产量递减分析方法的对比

LGM 方法由Clark 等^[39] 于2011 年提出，对应的表达式如下：

将单井生产数据按照7∶3 的比例划分训练集和测试集，即利用70% 的生产数据开展历史拟合，然后利用剩余30% 的数据开展测试。以X10-5 井的生产数据为例，展示本文方法的关键步骤。X10-5 井的生产数据如图5 所示，共计1 229 d，利用前860 d 的生产数据，以累计产量误差最小为目标，采用粒子群算法开展历史拟合，确定拟合参数，对应的历史拟合结果如图6 所示。由于生产数据点较多，容易掩盖曲线变化，因此图6 只展示部分生产数据点。整体来看，本文方法给出的结果与实际值较为一致。拟合得到的p_i = 55.3 MPa，a₁ = － 6.06×10⁻⁹ m⁻³， c_f = 55.99×10⁻⁴ MPa⁻¹。拟合得到的p_i 大于原始地层压力（37.55 MPa），这是因为大规模压裂后泵入的压裂液起到了增压作用；p_i 小于该井近跟端压裂段停泵压力与液柱压力之和71.43 MPa。拟合得到的c_f 大于页岩的压缩系数，小于水力裂缝的压缩系数。因此， 拟合得到的参数具有物理意义。

利用拟合得到的参数，基于式（13）计算了动态泄流区平均压力，如图7 所示，将动态泄流区平均压力作为一个额外特征，进一步结合历史生产数据中的井底流压和日产量，基于式（16）计算了采气指数。此外，利用动态泄流区平均压力和日产量， 基于式（8）计算了物质平衡拟时间，并将其作为另一个额外特征。在双对数坐标系中绘制了采气指数和物质平衡拟时间的散点图（图8），从图8 中可以看出，数据点收敛为一条斜率为负数的直线，通过该直线的斜率和截距可以确定式（15），从而获得了采气指数递减方程。基于该递减方程，结合递推原理， 可以预测未来产量，预测结果如图9 所示。整体来看，本文方法预测结果与实际值较为吻合。进一步，开展了对比分析，对于X10-5 井，本文方法、LGM、SEPD 和Duong 方法的拟合结果如图6 所示，经验产量递减分析方法对应的拟合参数见表2 所示。从图6 中可以看出，本文方法的拟合效果最好，对应均方根误差（RMSE） 为15.25×10⁴ m³ ；LGM 和SEPD 方法的拟合效果次之，对应的RMSE 分别为71.99×10⁴ m³ 和81.04×10⁴ m³ ；Duong 方法的拟合效果最差，对应的RMSE 达到了202.99×10⁴ m³。进一步，预测了后369 天的累计产量，并与实际的累计产量进行了对比。从图9-a 中可以看出，本文方法的预测效果最好，预测结果与实际值较为接近，对应的RMSE 为16.73×10⁴ m³， 对应的相对误差（RE） 为0.46% ；LGM 和Duong 方法的预测结果偏高，对应的RMSE 分别为151.3×10⁴ m³ 和677.31×10⁴ m³，对应的RE 分别为5.83%和26.12%；SEPD 方法的预测效果最差， 对应的RMSE 高达2.4×10⁷ m³，RE 为97.96%。图9-b 展示了各方法对应的长期预测结果，共20 年，7 300 d。从图9-b 中可以看出，20 年后Duong 方法给出的预测结果为4.32×10⁸ m³，不符合常规认识。Duong 方法假设气井长期处于线性流状态，若后期气井流动状态发生改变，再用以前的生产数据来预测产量，将会得到比实际值大得多的结果^[41]。20 年后本文方法的预测结果为0.3×10⁸ m³，LGM 和SEPD 方法给出的预测结果分别为0.56×10⁸ m³ 和0.24×10⁸ m³。

5.2 　与时间序列分析方法的对比

接下来，将本文方法与时间序列分析方法进行对比。选择了两种时间序列分析方法，即NAR 和LSTM。NAR 具有自适应性和自组织特性，是一种常用的动态神经网络，属于浅层学习的范畴，其在传统的神经网络基础上，添加了反馈和延迟单元，因此具备了“记忆能力”^[42]。NAR 通过反馈和延迟单元保留历史数据，并将其应用于预测，一般由输入层、滞后层、隐藏层和输出层组成，其对应的超参数为：隐藏层中的神经元数量（HN_N）和延迟阶数（lag_N）。

LSTM 是一种改进的循环神经网络，其解决了循环神经网络梯度“爆炸”和“消失”的问题，可以学习到时间跨度较大的依赖关系^[43]。LSTM 的隐藏层由1 个记忆单元和3 个控制单元（“门”结构）组成：记忆单元“C”、输入门“i”、遗忘门“f”和输出门“o”。记忆单元用于传递信息，“门”结构则用于添加和移除信息，从而实现长期记忆。LSTM 采用经典的循环策略开展预测，对于产量预测，其具体过程如下：①采用过去几天的累计产量数据预测未来第一天的累计产量；②采用过去几天的累计产量和未来第一天的累计产量（上一步的预测值）预测未来第二天的累计产量；③基于上述循环策略预测未来的累计产量。LSTM 涉及2 个超参数：延迟天数（lag_L）和隐藏层单元数目（HN_L）。

同样利用X10-5 井的生产数据开展测试，将生产数据按照7∶3 的比例划分训练集和测试集。首先开展历史拟合，本文方法、NAR 和LSTM 方法的拟合结果如图10-a 所示，对应的拟合参数见表3。整体来看，对于训练集，3 种方法的拟合效果较好， 对应的结果与实际值较为一致，本文方法、NAR 和LSTM 方法对应的RMSE 分别为15.25×10⁴ m³、 9.6×10³ m³、3.24×10⁴ m³。但是，对于测试集，本文方法预测效果最好，预测结果与实际值较为一致（图10-b），对应的RMSE 为16.73×10⁴ m³，对应的RE 为0.46%。而NAR 和LSTM 方法的预测结果与实际值存在差异，对应的RMSE 分别为67.41×10⁴ m³ 和82.19×10⁴ m³，对应的RE 分别为2.32% 和2.75%。随着生产时间的延长，NAR 和LSTM 方法的预测结果与实际值的差异逐渐增大，LSTM 方法的预测结果还出现了锯齿状波动。图10-c 展示了各方法的长期预测结果，NAR 和LSTM 方法预测结果近似呈直线，与实际情况不符。这是因为：① NAR 和LSTM 方法是纯数据驱动方法，开展外推预测时需要用到之前的预测值，若之前的预测值存在误差，那么在物理机制缺失的情况下，随着预测步数的增加，误差会逐渐累积且不受控制，从而导致长期预测结果失真^[24] ；② NAR 和LSTM 方法仅考虑了预测变量（累计产量）与自身历史变化之间的关系，忽略了井底流压变化的影响，导致长期预测结果的准确性较差。相对于NAR 和LSTM 方法，本文方法预测的累计产量随时间递增，当t = 7 300 d，预测的累计产量为0.3×10⁸ m³。

5.3 　不同方法的综合对比分析

进一步结合研究区块其他5 口井的生产数据， 开展了测试。图11 展示了各方法对应的历史拟合、短期和长期预测结果，各方法对应的RMSE 如图12 所示。整体来看，对于训练集，时间序列分析方法的非线性拟合能力强，历史拟合效果最好，对应的RMSE ＜ 10×10⁴ m³ ；本文方法拟合效果次之，经验产量递减分析方法的拟合效果较差。

从图12 中可以看出， 对于NAR 方法， 测试集和训练集对应RMSE 的差异最大，具体的，对于X3-2 井，训练集对应的RMSE 为5.18×10³ m³，测试集对应的RMSE 高达355.49×10⁴ m³。这说明NAR 方法的拟合效果好，但是预测效果差，即泛化能力不强。同样，对于LSTM 方法，测试集和训练集对应RMSE 的差异较大。LSTM 和NAR 方法同属于纯数据驱动方法，其非线性拟合能力强，但可解释性差， 还存在误差易累积且不可控的问题。因此，NAR 和LSTM 方法不适用于长期预测。从图12 中还可以看出，对于本文方法，测试集和训练集对应RMSE 的差异较小，说明本文方法的泛化能力较强。

对于长期预测，Duong 方法出现了“翘头”的现 象，即累计产量随时间的增加快速增加，如图11-a、c、e，这与实际情况不符。因此，Duong 方法会高估气井产量。对于部分井，LSTM 和NAR 方法的长期预测结果近似为水平直线（图11-c ～ e），同样不符合实际。SEPD 方法的预测结果一般偏保守^[44]，本文方法给出的长期预测结果一般大于SEPD 方法给出的结果，这也间接验证了本文方法的正确性。此外，本文方法的长期预测结果未出现波动，也未出现预测结果为水平直线的情况。因此，本文方法的长期预测能力优于其他方法。

6 　结论

1）基于机理—数据融合建模的思路，利用动态泄流区孔隙体积计算方程和物质平衡方程，生成了动态泄流区平均压力和物质平衡拟时间两个额外特征， 进一步结合递减规律、特征直线分析和递推原理，建立了物理—数据协同驱动的页岩气井产量预测方法。

2）物理—数据协同驱动的页岩气井产量预测方法简化了对储层非均质性的描述，避开了裂缝网络精确识别和表征的难题，克服了数值模拟计算量大的不足，解决了纯数据方法物理意义不明确、可解释性差的问题。

3）物理—数据协同驱动的页岩气井产量预测方法自动适应变流压—变产量的情况，精度高且泛化能力强。生产数据测试表明：本文方法精度高，长期预测结果稳定，优于LGM、Duong 和SEPD 经验产量递减分析方法，也优于NAR、LSTM 时间序列分析方法。

符　号　说　明