一种基于线性动态系统的油气田注采产量动态预测方法与流程

1.本发明涉及油气开采技术领域，更为具体地，涉及一种基于线性动态系统的油气田注采产量动态预测方法。


背景技术：

2.水驱法是油气藏开发过程中应用较为广泛的一种方法。在对低渗油气藏中水驱采气采油的过程进行动态预测时，传统的动态预测方法，如水驱特征曲线法、物质平衡方程预测法、油藏数值模拟法等在低渗储层的预测效果不是很理想，主要原因为：传统方法主要适用于高渗储层，不能表征出低渗储层中长期存在的瞬态流。目前低渗油气藏存在广泛且储量丰富，对其产量的预测是今后的研究重点。
3.在以注水为主要生产驱动方式的油气田生产动态预测的领域中，目前应用较广泛的为a.a.youself与p.gentil等人2005提出的电容模型(cm)，该模型是在alejandro albertoni与larry w.lake提出的井间连通性模型的基础上发展而来的。电容模型量化了直井之间的连通性，同时考虑了储层中岩石压缩性和流体传输性，更好地解释了注入速率与生产速率之间的变化关系。feilong liu与jerry m.mendel等根据测得的产量与注入速率数据，使用扩展的卡尔曼滤波器(ekf)预测了多个注入井和单个生产井之间的注采关系(ipr)，该研究利用连续时间单峰函数来模拟注入井和生产井之间的瞬间冲激响应，通过该方法可以表征出注入生产井之间流体的流动行为。mohammad soroush在传统电容模型的基础上考虑了生产操作对最终产量的影响。zequn zhang与heng li等人在2015扩展了电容模型，并建立了多层电容电阻模型(mlcrm)用于预测分层油藏中的水驱性能，同时应用集合卡尔曼滤波(enkf)方法估算多层电容电阻模型(mlcrm)中每一层的连通系数和井指数。现有技术中很少有将线性动态系统(lds)用于油气开采技术领域，检索公开专利仅获得一种用于油气井钻井过程的溢漏风险协同监测方法(202110160988.8)，属于钻井技术领域与本发明要应用的技术领域差别明显。
4.上述现有技术中模型普遍存在共同的缺点：参数不能通过数据学习得到，需要事先将这些参数设置为固定值，这通常会降低预测精度，尤其是在低渗油藏中。除此之外，现有技术不能对水驱油气藏中流体的瞬间流动做出较好的表征，不能较好的解决低渗水驱油气藏中产量精确预测的技术问题。


技术实现要素：

5.本发明的目的在于解决上述现有的方法无法解决低渗水驱油气藏中产量精确预测的技术问题，提供一种基于线性动态系统的油气田注采产量动态预测方法。
6.本发明是通过下述技术方案实现：一种基于线性动态系统的油气田注采产量动态预测方法，包括以下步骤：
7.步骤s1：基于水驱油气田整个储层的控制物质平衡方程，建立起线性动态系统(lds)产量预测模型；
8.步骤s2：利用em(expectation maximization，期望最大化)方法对所述线性动态系统(lds)产量预测模型中的模型参数做出估算；
9.步骤s3：动态数据的整理与划分，将数据集按照一定比例划分为训练集与测试集；
10.步骤s4：所述模型参数的初始化以及参数计算，整合预测结果并计算相对误差；
11.在本发明的一较佳实施方式中，在所述步骤s1中，在所述水驱油气田的开采过程中，所述整个储层的控制物质平衡方程为：
[0012][0013]
其中，c
t
是总可压缩性，mpa
‑1；v
p
是储层的孔隙体积，m3；是储层压力的平均变化率，mpa/t；i(t)和q(t)分别表示时间t处的总注入率和生产率，m3/t。
[0014]
在本发明的一较佳实施方式中，在所述步骤s1中，所述线性动态系统(lds)产量预测模型主要包含三部分：
[0015]
隐藏状态方程：z
n 1
＝az
n
bu
n
w,w～n(0,q)
[0016]
观测方程：x
n
＝hz
n
cu
n
v,v～n(0,r)
[0017]
模型的初始状态：
[0018]
其中，s为初始状态变量的噪声系数，无因次；w为系统噪声系数，无因次；v为观察噪声系数，无因次；q，r，v0分别为w，v，s的协方差矩阵系数，无因次。
[0019]
在本发明的一较佳实施方式中，在所述步骤s2中，分为两步对所述模型参数进行迭代更新，第一步为人为初始化所述模型参数之后得到数据的预测值；第二步为在利用得到的数据的预测值反过来估计所述模型参数；所有的计算过程通过python编程实现。
[0020]
在本发明的一较佳实施方式中，在所述步骤s3中，所述训练集：所述测试集＝7：3。
[0021]
在本发明的一较佳实施方式中，在所述步骤s4中，所述训练集与所述测试集的相对误差分别小于5％和10％即为满足精度要求；最终的月产量预测值等于10次所述预测结果的平均值。
[0022]
相比现有技术，本发明至少有如下优点：本发明针对注水驱动为主要方式的低渗油气田，基于线性动态系统(lds)建立了油气田注采产量动态预测模型，模型中隐变量的引入更好的表征了实际生产中的注采关系，可以对月产量做出准确的预测。
附图说明
[0023]
图1为本发明的基于线性动态系统的油气田注采产量动态预测方法的流程图；
[0024]
图2为本发明的实施例五点法面积井网中注入井到生产井q1的控制体积示意图；
[0025]
图3为本发明的实施例某区块月注入量与月产量变化图；
[0026]
图4为本发明的实施例某区块动态系统(lds)模型十次预测值与实际值变化对比图:(a)全部时间段，(b)测试集时间
[0025]
图5为本发明的实施例某区块动态系统(lds)模型预测平均值与实际值变化对比图。
具体实施方式
[0028]
下面结合实施例，对本发明作进一步详细说明，但本发明的实施方式不仅限于此。
[0029]
图1示出了根据本发明的基于线性动态系统的油气田注采产量动态预测方法流程。
[0030]
以某地区油气田某区块进行lds预测模型的研究。该区块自1992年6月至2009年12月布有13口注入井与11口生产井，在对月生产数据进行数据预处理之后，筛选出具有有效区块月总注入量与区块月总生产量数据共210组，采用本发明提供的一种基于线性动态系统的油气田注采产量动态预测方法，包括以下步骤：
[0031]
步骤s1：基于水驱油气田整个储层的控制物质平衡方程，建立起线性动态系统(lds)产量预测模型；
[0032]
在水驱油气田的开采过程中，整个储层的控制物质平衡方程为：
[0033][0034]
其中，c
t
是总可压缩性，mpa
‑1；v
p
是储层的孔隙体积，m3；是储层压力的平均变化率，mpa/t；i(t)和q(t)分别表示时间t处的总注入率和生产率，m3/t。
[0035]
假定储层压力是时间和位置的一般函数，的变化率为：
[0036][0037]
其中，是代表储层中特定地质位置的向量。
[0038]
设储层体积则整个储层的控制物质平衡方程将转换为：
[0039]
v
′
(t)＝i(t)
‑
q(t)
[0040]
其中，v
′
(t)是控制体积的变化率，无因次，控制体积是指特定注入井和特定生产井之间的储层空间的体积，流体会通过该储层空间从注入井流向生产井，利用注入井的流体驱动进行油气生产，通过引进控制体积的定义表征了井间的连通性，因此可以更好的对储层中的生产过程进行描述与产能预测。
[0041]
对具有i口注入井与j口生产井的油藏储层模型建立储层的控制物质平衡方程，设其中i
i
(t)表示第i口注入井的注入速率，其中q
j
(t)表示第j口生产井的生产速率，并使v
ij
表示从注入井i到生产井j的控制体积，v
′
ij
(t)表示从注入井i到生产井j的控制体积的变化率，由此得到变换之后的注入井i与生产井j的物质平衡方程：
[0042]
v
′
ij
(t)＝i
t
(t)
‑
q
j
(t)
[0043]
其中，v
′
ij
是时间和位置的函数，反应了不同地质位置的储层体积的动态变化。
[0044]
假定对于特定的生产井j，其生产率q
j
受所有注入井的影响，则使用公式表示所有i个注入井对于指定生产井j的控制体积变化率，对于指定生产井j来讲，每个注入井i对其产量的贡献率是不同的，所以引入加权贡献率的思想，一
般来讲注入井距离生产井越远，对产量的贡献率越小，所以注入井i与生产井j的物质平衡方程变为：
[0045][0046]
其中，α和β为系数列向量，向量中的每个数代表每个注入井对特定生产井的贡献率的权重。
[0047]
以图2中五点法面积井网的生产井q1为例，它受到注入井i1,
…
,i5的直接影响，由于距离每个注入井的距离不同，所以每个注入井对生产井的贡献不同，由图中控制体积(图中箭头)的大小，可以得到：注入井i1，i2，i3对生产井q1的贡献率大于注入井i4，i5的贡献率,并且得到生产井q1的物质平衡方程为：
[0048][0049]
将注入速率作为系统输入，将生产率作为系统输出；将控制体积变化率作为隐藏状态变化量，基于线性动态系统(lds)建立储层的生产预测模型。其中控制体积变化率的状态更新不仅与先前的控制体积变化率有关，也受到系统输入(即注入速率)的影响。使表示生产井j在时间t δt的控制体积变化率，则：
[0050][0051]
其中，ε(t)为用于更新的高斯噪声。
[0052]
将上式时间步长离散化得到：
[0053][0054]
其中，n为离散时间步长。
[0055]
为了简化计算过程，我们假设这里的f函数为线性递归函数，即有：
[0056][0057]
其中，的系数矩阵，b为的系数矩阵。
[0058]
同样，将上式中的物质平衡方程离散化，得到生产井j的观测方程：
[0059][0060]
其中，γ(n)为已知生产率的高斯噪声。
[0061]
综合上述两个式子可以得到生产井j的线性动态系统模型：
[0062][0063][0064]
为了简化表示方法，使向量z
n
表示控制体积变化率，表示z
n
的初始值，即
使向量u
n
表示系统输入，即注入率使标量y
n
表示系统输出，即生产井j的生产率q
j
(n)，由此得到系统方程，主要包含三部分：
[0065]
隐藏状态方程为：
[0066]
z
n 1
＝az
n
bu
n
w,w～n(0,q)
[0067]
观测方程为：
[0068]
x
n
＝hz
n
cu
n
v,v～n(0,r)
[0069]
模型的初始状态为：
[0070][0071]
其中，s为初始状态变量的噪声系数，无因次；w为系统噪声系数，无因次；v为观察噪声系数，无因次；q，r，v0分别为w，v，s的协方差矩阵系数，无因次。
[0072]
步骤s2：利用em方法对线性动态系统lds产量预测模型中的模型参数做出估算；
[0073]
分为两步对参数进行迭代更新，所有的计算过程通过python编程实现。
[0074]
第一步：人为初始化模型参数θ之后得到数据的预测值，这个过程分为两个步骤。
[0075]
(1)卡尔曼滤波过程，这属于一种前向递归过程，计算公式如下：
[0076][0077][0078][0079][0080][0081]
(2)卡尔曼平滑过程，这属于一种后向递归过程，计算公式如下：
[0082][0083][0084][0085][0086][0087][0088][0089]
第二步：在利用得到的数据的预测值反过来估计模型参数θ，主要分为两个步骤。
[0090]
(1)e步，由假设可知状态变量z
n
与观测变量x
n
的联合分布函数为：
[0091][0092]
其中，p(z1)，p(z
n
|z
n
‑1)与p(x
n
|z
n
)均满足高斯分布，且有：
[0093][0094][0095][0096]
分别计算上述三个式子的q(θ,θ
old
)＝e
z
[log p(x,z|θ)]得到：
[0097][0098][0099][0100]
(2)m步，为求极大值需对q1，q2，q3分别关于参数求导，并令其导数为零得到：
[0101]
模型参数u0,v0的估计式为：
[0102]
u0＝e[z1],
[0103]
模型参数a,b的估计值需要联立方程如下：
[0104][0105][0106]
模型参数h,c的估计值需要联立方程如下：
[0107][0108][0109]
模型参数q的估计式为：
[0110][0111]
模型参数r的估计式：
[0112][0113]
其中：a,b与h,c代入的数值为更新之后的模型参数值。
[0114]
在模型参数更新的过程中有三个重要的期望值，它们均可以根据第一步中卡尔曼滤波 平滑过程的计算结果得到，计算公式包括：
[0115][0116]
[0117][0118]
需要求解的模型参数为θ＝{a,b,h,c,q,r,v0,u0}，其中初始化参数是人为给定的，之后在利用em方法进行模型参数的更新。
[0119]
步骤s3：动态数据的整理与划分，将数据集按照一定比例划分为训练集与测试集；
[0120]
图3为该区块的月注入量与月产量数据随时间变化图，按照训练集：测试集＝7：3的比例，设定前147月的区块月总注入量与区块月总产量为lds模型训练集，后63月的注入量与产量数据作为预测集。
[0121]
步骤s4：模型参数的初始化以及参数计算，整合预测结果并计算相对误差；
[0122]
利用随机正态分布初始化模型参数，区块月总注入量作为lds模型系统输入，设定的最大迭代次数为500次，为了防止模型的过拟合，误差ε设置为0.0001，经过对lds模型总共迭代训练了10次，得到lds模型系统输出与该区块实际生产数据的关系，如图4所示。
[0123]
求得该区块10次迭代的实际产量与预测产量的相对误差，由表1可以看出模型的预测结果是有一定稳定性的，80％的模型的相对误差稳定在6％左右，将10次的预测结果进行平均计算得到最终的月产量预测值，如图5所示，并得到了训练集与测试集的相对误差分别为1.55％与5.49％，预测精度较高，足以满足生产、研究需求。
[0124]
表1区块lds模型十次预测值与实际值mape统计表
[0125][0126]
以上所述，仅是本发明的较佳实施例并非对本发明对任何形式上的限制，凡是依据本发明的技术、方法实质上以上实施例所作的任何简单修改、等同变化，均落入本发明的保护范围之内。