本发明公开了一种适用于低渗透裂缝性油藏的快速模拟及优化方法，属于非常规油藏数值模拟领域，包括如下步骤：步骤1、将油藏系统抽象为基质网格和裂缝网格，根据基质和裂缝的渗流特征构建双重介质模型；步骤2、利用ES‑MDA历史拟合算法，结合约束条件，建立自动历史拟合数学模型；步骤3、将经济净现值作为生产优化的目标函数，通过对油藏地质储量、各井压力以及注采液量上下限进行约束，建立油藏注采优化模型，采用差分进化算法得到使油藏经济净现值最大化的生产调控制度。本发明有助于更好地了解
1.一种适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，包括如下步骤：步骤1、将油藏系统抽象为基质网格和裂缝网格，根据基质和裂缝的渗流特征构建双重介质模型；具体过程如下：步骤1.1、按照油藏系统中井点的分布将储层离散为一系列的节点，对所有节点两两相连形成初始井网模型；
步骤1.2、引入判定条件对初始井网模型进行筛选，生成两两相连的井间连接管道，并将储层体积划分到每条连接管道中；
步骤1.3、在每条连接管道上划分网格，分为基质网格和裂缝网格，所有的基质网格构成基质系统，所有的裂缝网格构成裂缝系统；网格间流体流动能力由网格间传导率大小决定；
步骤1.4、在基质系统中考虑非达西渗流，在裂缝系统中考虑岩石的应力敏感效应；
步骤1.5、结合基质和裂缝内流体的渗流规律，构建低渗透裂缝性油藏双重介质模型；
双重介质模型中包括基质系统的连续性方程和裂缝系统的连续性方程；
其中， 为油相密度； 为基质的孔隙度； 为基质中油相饱和度； 为基质网格体积；为时间步序号； 表示散度； 为油相在基质的渗流速度； 为基质网格内单位时间单位体积内油相流体体积流量的变化值； 为水相密度； 为水相在基质内的饱和度； 为水相在基质的渗流速度； 为基质网格内单位时间单位体积内水相流体体积流量的变化值；
其中， 为裂缝的孔隙度； 为油相在裂缝内的饱和度； 为油相在裂缝内的渗流速度； 为裂缝网格内单位时间单位体积内油相流体体积流量的变化值； 为水相在裂缝内的饱和度； 为水相在裂缝内的渗流速度； 为裂缝网格单位时间单位体积内水相流体体积流量的变化值；
步骤1.6、对连续性方程进行线性插值近似，利用Newton‑Raphson方法全隐式求解裂缝及基质系统中的参数；
步骤2、利用ES‑MDA历史拟合算法，结合约束条件，建立自动历史拟合数学模型；
步骤3、将经济净现值作为生产优化的目标函数，通过对油藏地质储量、各井压力以及注采液量上下限进行约束，建立油藏注采优化模型，采用差分进化算法得到使油藏经济净现值最大化的生产调控制度。
2.根据权利要求1所述适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，所述步骤1.2中，判定条件包括距离条件、角度条件和最大连接点数条件；其中，距离条件的阈值为最小井距的3倍，超过阈值的井视为不连接；角度条件具体设置为：任意三条连接管道组成的三角形中，最小内角不能小于设定的最小角度；最大连接点数条件具体设置为：每个井点允许连接的井点数量为4‑8个。
3.根据权利要求2所述适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，所述步骤1.3中，划分网格是在每条连接管道中平均分为5‑10个网格，当指定井点网格体积时，井点网格体积默认值的计算公式如下： (1)；
其中， 为井点 所在网格体积； 为总连接体积； 为连接管道数量； 为每根连接管道中网格数量； 为井点数量；
确定好井点网格体积后，连接管道中网格体积的计算公式如下： (2)；
其中，、为位于连接管道两端的不同井点；表示连接管道中第 个网格； 表示 、井点的连接管道中第 个网格的体积； 为 、井点间连接管道的体积； 和 分别为井点和井点 所在网格体积； 为井点 连接的管道数量； 为井点 连接的管道数量；
基质网格和裂缝网格中，裂缝网格体积 、基质网格体积 计算如下： (3)；
其中， 表示粗网格内裂缝间距单元数； 表示微裂缝占整个单元的体积分数； 为油藏中微裂缝间距； 为油藏系统中微裂缝开度；
其中， 为网格 和网格 间传导率； 为网格 和网格 的接触面积； 为网格 的中心到网格 和网格 的接触面的距离； 为网格 的中心到网格 和网格 的接触面的距离；
基质与裂缝间存在流体流动，流体的窜流计算如下： (8)；
其中，为基质裂缝间的窜流量； 为形状因子； 为流体密度； 为基质网格的渗透率；为流体粘度； 为基质系统压力； 为裂缝系统压力； 为基质网格体积。
4.根据权利要求3所述适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，所述步骤1.4中，考虑非达西渗流后，基质系统的运动方程为： (9)；
其中， 为基质渗流速度； 为基质网格的渗透率； 为流体粘度； 为流动压差； 和为两个不同的非线性渗流系数；
考虑岩石的应力敏感效应后，裂缝系统的运动方程为： (10)；
其中， 为裂缝渗流速度； 为裂缝网格的渗透率；为常数； 为流体粘度； 为应力敏感系数； 为裂缝网格压力； 为裂缝网格初始压力； 为流动压差；
将公式（9）和公式（10）统一表述为如下运动方程的形式： (11)；
其中， 为某相流体； ，表示 为油相 或水相 ； 为 的渗流速度； 为 的拟渗透率； 为 的流体粘度； 为流动压差；
考虑非达西渗流和应力敏感效应后，基质和裂缝的拟渗透率表示如下： (12)；
其中， 和 分别表示基质和裂缝； 为网格 和网格 间考虑非达西渗流后基质的拟渗透率； 分别为网格 和网格 间考虑应力敏感后裂缝 内的拟渗透率； 为网格 和网格 间基质 的初始渗透率； 表示网格 和网格 间裂缝 的初始渗透率； 为流体的粘度； 为基质内网格 和网格 两节点的压力梯度； 和 为两个不同的非线性渗流系数； 为应力敏感系数； 为网格 和网格 间裂缝的当前压力； 为网格 和网格 间裂缝的初始压力。
5.根据权利要求4所述适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，所述步骤1.6的具体过程为：步骤1.6.1、定义流度系数如下：
其中， 为某相流体； ，表示 为油相 或水相 ； 为某相流体流度； 为某相流体的相对渗透率； 为某相流体粘度；
步骤1.6.2、对基质系统和裂缝系统采用同样的格式进行线性插值近似，线性插值近似采用如下的数值离散方程： (23)；
式中， 为某相流体的密度；为孔隙度； 为某相流体的饱和度； 为网格体积；为模拟迭代过程中的上一时刻； 为模拟迭代过程中的当前时刻； 为时间步长； 为网格相邻的单元集合；下标 表示两个相邻网格 和 界面上某一参数的平均值； 为网格 和网格 间传导率； 为某相流体在 时刻网格 的压力； 为某相流体在 时刻网格 的压力； 为网格 内某相流体在 时刻的注采速率；
步骤1.6.3、将上述数值离散方程写为如下的残差形式： (24)；
采用Newton‑Paphson迭代求解方法对裂缝及基质系统中的压力和含水饱和度参数进行求解，求解过程为：从初始点开始，通过不断的迭代计算，得到一个接近函数零点的新点，这个新点是通过在当前点处进行泰勒展开计算得到的，然后将新点作为下次迭代的初始点继续进行迭代计算，直到找到满足要求的解。
6.根据权利要求5所述适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，所述步骤2的具体过程如下：步骤2.1、构建自动历史拟合数学模型的目标函数，如下： (25)；
式中， 为目标函数值； 表示模型计算结果； 表示实际观测数据； 表示观测误差的协方差矩阵；
步骤2.2、在自动历史拟合数学模型的历史拟合过程中，仅需要调整每根管道的体积、网格间传导率以及幂律相渗曲线参数；
式中， 和 分别表示油相和水相的相对渗透率； 表示无因次饱和度； 和 分别决定水相和油相的相对渗透率曲线的弯曲程度； 表示束缚水饱和度； 表示残余油饱和度； 表示在含水饱和度为束缚水饱和度时的水相相对渗透率值； 为含水饱和度值；
其中， 为网格 和网格 间传导率，包括基质与基质网格间、裂缝与裂缝网格间、基质与裂缝网格间； 为 、井点的连通管道的基质传导体积， 为 、井点的连通管道的裂缝传导体积； 、 、 为 、井点的连通管道相渗曲线的不同参数；和 为两个不同的非线性渗流系数；
当所有井对间共用一组相渗数据时，历史拟合参数 表示如下： (29)；
其中， 、 、 为所有井对间共用的三个相渗曲线参数；
步骤2.4、在历史拟合过程中，对模型参数增加约束；自动历史拟合数学模型的约束条件如下： (30)；
其中， 为网格 和网格 的基质体积； 为基质系统的总体积； 为网格 和网格 的裂缝体积； 为裂缝系统的总体积； 为总的井数量。
7.根据权利要求6所述适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，所述步骤3的具体过程如下：步骤3.1、采用经济净现值作为生产优化的目标函数，数学表达式如下： (31)；
其中，为经济净现值；为控制变量；为时间步序号，表示总时间步数；表示生产井序号， 表示生产井数量； 表示注水井序号， 表示注水井数量；为每吨原油的价格； 为生产井 在 时间步的采油量； 为处理每吨污水的成本； 表示生产井 在 时间步的采水量； 表示每吨水的注水成本； 表示注水井 在 时间步的注水量；
步骤3.2、在生产优化过程中，采用边界约束条件对油藏参数进行控制，油藏参数的约束条件如下所示： (33)；
其中， 为注水井 在 时间步的注水量； 为生产井 在 时间步的采液量； 为注采井 在 时间步注采液量 的下限； 为注采井 在 时间步注采液量 的上限； 为总的地质储量； 和 分别为注采井 控制压力 的下限和上限；
步骤3.3、采用差分进化算法得到使油藏经济净现值最大化的生产调控制度。
8.根据权利要求7所述适用于低渗透裂缝性油藏的快速模拟及优化方法，其特征在于，所述步骤3.3的具体过程为：首先初始化形成一系列的控制变量 ，对这一系列控制变量通过以下方式进行变异，变异后的控制变量记为 ；
式中， 代表第 次迭代； 表示变异算子； 代表三个不同的个体序号； 表示第 次迭代变异后的序号为 的控制变量； 表示第 次迭代的序号为 的初始控制变量； 表示第 次迭代的序号为 的初始控制变量； 表示第 次迭代的序号为 的初始控制变量；
接下来进行交叉，交叉后的变量记为 ，对于每个 ，有： (35)；
其中， 表示0‑1的随机数； 表示交叉算子； 表示第 次迭代交叉后序号为 的控制变量的部分片段； 表示第 次迭代变异后序号为 的控制变量的部分片段； 表示第 次迭代序号为 的初始控制变量的部分片段；
经过变异交叉后的变量 ，对 中的每个子变量进行求目标函数操作，然后将 中第 个子变量 与第 次迭代中的第 个子变量 的函数值进行比较，如果 的函数值比 的大，则将 作为第 次迭代的 。
一种适用于低渗透裂缝性油藏的快速模拟及优化方法
技术领域
[0001] 本发明属于非常规油藏数值模拟领域，具体涉及一种适用于低渗透裂缝性油藏的快速模拟及优化方法。
背景技术
[0002] 随着常规油气资源的逐渐枯竭，低渗透裂缝性油藏的开发地位日益提高，是未来增储上产的主要油藏类型，因此如何对其高效开发是我国应对能源危机的重中之重。与常规砂岩油藏不同，低渗透裂缝性油藏空隙吼道小、渗流阻力大，且流体流动过程中存在启动压力梯度现象，因此传统的达西渗流方程不再适应，加上该类油藏大多采用压裂的方式增加其渗流通道，所以油藏中裂缝充分发育，使用单一的介质模型进行模拟并不合适，目前比较常见的模型有双重介质模型和离散裂缝模型，两种方式都能够很好地对裂缝系统进行描述，但是离散裂缝模型需要提供裂缝的具体参数，包括延展性、裂缝开度、裂缝间距等，这些参数相对来说较难获取，因此本发明采用的是双重介质模型对低渗透裂缝性油藏进行建模并进行后续生产优化。
发明内容
[0003] 为了解决上述问题，本发明提出了一种适用于低渗透裂缝性油藏的快速模拟及优化方法，首先根据油藏的布井方式将其划分为一系列井与井相连的连通管道，并在每根管道上划分网格，将划分好的网格矩阵抽象为基质网格和裂缝网格，每个网格有自己的渗透率、有效厚度、含水饱和度等参数，从而构建双孔双渗模型，有效实现低渗透裂缝性油藏中基质与裂缝内流体流动数值模拟的精确描述和刻画。同时利用ES‑MDA算法对油田生产历史进行拟合，对所建模型进行实时更新，并结合最优化理论，以经济净现值或累积采油量为目标函数，以各井及区块压力及地质储量为约束，根据各井的生产历史进行后续的注采调控，优化结果对油田增储上产、控水增油具有指导意义。
[0004] 本发明的技术方案如下：
[0005] 一种适用于低渗透裂缝性油藏的快速模拟及优化方法，包括如下步骤：
[0006] 步骤1、将油藏系统抽象为基质网格和裂缝网格，根据基质和裂缝的渗流特征构建双重介质模型；
[0007] 步骤2、利用ES‑MDA历史拟合算法，结合约束条件，建立自动历史拟合数学模型；
[0008] 步骤3、将经济净现值作为生产优化的目标函数，通过对油藏地质储量、各井压力以及注采液量上下限进行约束，建立油藏注采优化模型，采用差分进化算法得到使油藏经济净现值最大化的生产调控制度。
[0009] 进一步地，步骤1的具体过程如下：
[0010] 步骤1.1、按照油藏系统中井点的分布将储层离散为一系列的节点，对所有节点两两相连形成初始井网模型；
[0011] 步骤1.2、引入判定条件对初始井网模型进行筛选，生成两两相连的井间连接管道，并将储层体积划分到每条连接管道中；
[0012] 步骤1.3、在每条连接管道上划分网格，分为基质网格和裂缝网格，所有的基质网格构成基质系统，所有的裂缝网格构成裂缝系统；网格间流体流动能力由网格间传导率大小决定；
[0013] 步骤1.4、在基质系统中考虑非达西渗流，在裂缝系统中考虑岩石的应力敏感效应；
[0014] 步骤1.5、结合基质和裂缝内流体的渗流规律，构建低渗透裂缝性油藏双重介质模型；双重介质模型中包括基质系统的连续性方程和裂缝系统的连续性方程；
[0015] 步骤1.6、对连续性方程进行线性插值近似，利用Newton‑Raphson方法全隐式求解裂缝及基质系统中的参数。
[0016] 进一步地，步骤1.2中，判定条件包括距离条件、角度条件和最大连接点数条件；其中，距离条件的阈值为最小井距的3倍，超过阈值的井视为不连接；角度条件具体设置为：任意三条连接管道组成的三角形中，最小内角不能小于设定的最小角度；最大连接点数条件具体设置为：每个井点允许连接的井点数量为4‑8个。
[0017] 进一步地，步骤1.3中，划分网格是在每条连接管道中平均分为5‑10个网格，当指定井点网格体积时，井点网格体积默认值的计算公式如下：
[0018] (1)；
[0019] 其中， 为井点 所在网格体积； 为总连接体积； 为连接管道数量；为每根连接管道中网格数量； 为井点数量；
[0020] 确定好井点网格体积后，连接管道中网格体积的计算公式如下：
[0021] (2)；
[0022] 其中，、为位于连接管道两端的不同井点； 表示连接管道中第 个网格；表示 、 井点的连接管道中第 个网格的体积； 为 、 井点间连接管道的体积； 和分别为井点 和井点 所在网格体积； 为井点 连接的管道数量； 为井点
连接的管道数量；
[0023] 基质网格和裂缝网格中，裂缝网格体积 、基质网格体积 计算如下：
[0024] (3)；
[0025] (4)；
[0026] (5)；
[0027] (6)；
[0028] 其中， 表示粗网格内裂缝间距单元数； 表示微裂缝占整个单元的体积分数； 为油藏中微裂缝间距； 为油藏系统中微裂缝开度；
[0029] 网格间传导率计算公式如下：
[0030] (7)；
[0031] 其中， 为网格 和网格 间传导率； 为网格 和网格 的接触面积； 为网格 的中心到网格 和网格 的接触面的距离； 为网格 的中心到网格 和网格的接触面的距离； 为网格 拟渗透率 和网格 拟渗透率 的平均值；
[0032] 基质与裂缝间存在流体流动，流体的窜流计算如下：
[0033] (8)；
[0034] 其中， 为基质裂缝间的窜流量； 为形状因子； 为流体密度； 为基质网格的渗透率； 为流体粘度； 为基质系统压力； 为裂缝系统压力； 为基质网格体积。
[0035] 进一步地，步骤1.4中，考虑非达西渗流后，基质系统的运动方程为：
[0036] (9)；
[0037] 其中， 为基质渗流速度； 为基质网格的渗透率； 为流体粘度； 为流动压差； 和 为两个不同的非线性渗流系数；
[0038] 考虑岩石的应力敏感效应后，裂缝系统的运动方程为：
[0039] (10)；
[0040] 其中， 为裂缝渗流速度； 为裂缝网格的渗透率； 为常数； 为流体粘度；为应力敏感系数； 为裂缝网格压力； 为裂缝网格初始压力； 为流动压差；
[0041] 将公式（9）和公式（10）统一表述为如下运动方程的形式：
[0042] (11)；
[0043] 其中， 为某相流体； ，表示 为油相 或水相 ； 为 的渗流速度； 为 的拟渗透率； 为 的流体粘度； 为流动压差；
[0044] 考虑非达西渗流和应力敏感效应后，基质和裂缝的拟渗透率表示如下：
[0045] (12)；
[0046] (13)；
[0047] 其中， 和 分别表示基质和裂缝； 为网格 和网格 间考虑非达西渗流后基质 的拟渗透率； 分别为网格 和网格 间考虑应力敏感后裂缝 内的拟渗透率； 为网格 和网格 间基质 的初始渗透率； 表示网格 和网格 间裂缝的初始渗透率； 为流体的粘度； 为基质内网格 和网格 两节点的压力梯度；
和 为两个不同的非线性渗流系数； 为应力敏感系数； 为网格 和网格 间裂缝的当前压力； 为网格 和网格 间裂缝的初始压力。
[0048] 进一步地，步骤1.5中，基质系统的连续性方程为：
[0049] (18)；
[0050] (19)；
[0051] 其中， 为油相密度； 为基质的孔隙度； 为基质中油相饱和度； 为基质网格体积；为时间步序号； 表示散度； 为油相在基质的渗流速度； 为基质网格内单位时间单位体积内油相流体体积流量的变化值； 为水相密度； 为水相在基质内的饱和度； 为水相在基质的渗流速度； 为基质网格内单位时间单位体积内水相流体体积流量的变化值；
[0052] 裂缝系统的连续性方程为：
[0053] (20)；
[0054] (21)；
[0055] 其中， 为裂缝的孔隙度； 为油相在裂缝内的饱和度； 为油相在裂缝内的渗流速度； 为裂缝网格内单位时间单位体积内油相流体体积流量的变化值； 为水相在裂缝内的饱和度； 为水相在裂缝内的渗流速度； 为裂缝网格单位时间单位体积内水相流体体积流量的变化值。
[0056] 进一步地，步骤1.6的具体过程为：
[0057] 步骤1.6.1、定义流度系数如下：
[0058] (22)；
[0059] 其中， 为某相流体； ，表示 为油相 或水相 ； 为某相流体流度； 为某相流体的相对渗透率； 为某相流体粘度；
[0060] 步骤1.6.2、对基质系统和裂缝系统采用同样的格式进行线性插值近似，线性插值近似采用如下的数值离散方程：
[0061] (23)；
[0062] 式中， 为某相流体的密度； 为孔隙度； 为某相流体的饱和度； 为网格体积； 为模拟迭代过程中的上一时刻； 为模拟迭代过程中的当前时刻； 为时间步长； 为网格 相邻的单元集合；下标 表示两个相邻网格 和 界面上某一参数的平均值； 为网格 和网格 间传导率； 为某相流体在 时刻网格 的压力；为某相流体在 时刻网格 的压力； 为网格 内某相流体在 时刻的注
采速率；
[0063] 步骤1.6.3、将上述数值离散方程写为如下的残差形式：
[0064] (24)；
[0065] 式中， 为方程残差；
[0066] 采用Newton‑Paphson迭代求解方法对裂缝及基质系统中的压力和含水饱和度参数进行求解，求解过程为：从初始点开始，通过不断的迭代计算，得到一个接近函数零点的新点，这个新点是通过在当前点处进行泰勒展开计算得到的，然后将新点作为下次迭代的初始点继续进行迭代计算，直到找到满足要求的解。
[0067] 进一步地，步骤2的具体过程如下：
[0068] 步骤2.1、构建自动历史拟合数学模型的目标函数，如下：
[0069] (25)；
[0070] 式中， 为目标函数值； 表示模型计算结果； 表示实际观测数据；表示观测误差的协方差矩阵；
[0071] 步骤2.2、在自动历史拟合数学模型的历史拟合过程中，仅需要调整每根管道的体积、网格间传导率以及幂律相渗曲线参数；
[0072] 幂律相渗曲线参数由如下的幂律相渗模型调整：
[0073] (26)；
[0074] 其中，
[0075] (27)；
[0076] 式中， 和 分别表示油相和水相的相对渗透率； 表示无因次饱和度； 和分别决定水相和油相的相对渗透率曲线的弯曲程度； 表示束缚水饱和度； 表示残余油饱和度； 表示在含水饱和度为束缚水饱和度时的水相相对渗透率值； 为含水饱和度值；
[0077] 步骤2.3、进行历史拟合调整的参数 如下：
[0078] (28)；
[0079] 其中， 为网格 和网格 间传导率，包括基质与基质网格间、裂缝与裂缝网格间、基质与裂缝网格间； 为 、 井点的连通管道的基质传导体积， 为 、 井点的连通管道的裂缝传导体积； 、 、 为 、 井点的连通管道相渗曲线的不同参数；和 为两个不同的非线性渗流系数；
[0080] 当所有井对间共用一组相渗数据时，历史拟合参数 表示如下：
[0081] (29)；
[0082] 其中， 、 、 为所有井对间共用的三个相渗曲线参数；
[0083] 步骤2.4、在历史拟合过程中，对模型参数增加约束；自动历史拟合数学模型的约束条件如下：
[0084] (30)；
[0085] 其中， 为网格 和网格 的基质体积； 为基质系统的总体积； 为网格 和网格 的裂缝体积； 为裂缝系统的总体积； 为总的井数量。
[0086] 进一步地，步骤3的具体过程如下：
[0087] 步骤3.1、采用经济净现值作为生产优化的目标函数，数学表达式如下：
[0088] (31)；
[0089] 其中， 为经济净现值； 为控制变量；为时间步序号， 表示总时间步数； 表示生产井序号， 表示生产井数量； 表示注水井序号， 表示注水井数量； 为每吨原油的价格； 为生产井 在 时间步的采油量； 为处理每吨污水的成本；表示生产井 在 时间步的采水量； 表示每吨水的注水成本； 表示注水井 在 时间步的注水量；
[0090] 步骤3.2、在生产优化过程中，采用边界约束条件对油藏参数进行控制，油藏参数的约束条件如下所示：
[0091] (33)；
[0092] 其中， 为注水井 在 时间步的注水量； 为生产井 在 时间步的采液量； 为注采井 在 时间步注采液量 的下限； 为注采井 在 时间步注采液量 的上限； 为总的地质储量； 和 分别为注采井 控制压力 的下限和上限；
[0093] 步骤3.3、采用差分进化算法得到使油藏经济净现值最大化的生产调控制度。
[0094] 进一步地，步骤3.3的具体过程为：首先初始化形成一系列的控制变量 ，对这一系列控制变量通过以下方式进行变异，变异后的控制变量记为 ；
[0095] (34)；
[0096] 式中， 代表第 次迭代； 表示变异算子； 代表三个不同的个体序号；表示第 次迭代变异后的序号为 的控制变量； 表示第 次迭代的序号为
的初始控制变量； 表示第 次迭代的序号为 的初始控制变量； 表示第 次迭代的序号为 的初始控制变量；
[0097] 接下来进行交叉，交叉后的变量记为 ，对于每个 ，有：
[0098] (35)；
[0099] 其中， 表示0‑1的随机数； 表示交叉算子； 表示第 次迭代交叉后序号为 的控制变量的部分片段； 表示第 次迭代变异后序号为
的控制变量的部分片段； 表示第 次迭代序号为 的初始控制变量的部分片段；
[0100] 经过变异交叉后的变量 ，对 中的每个子变量进行求目标函数操作，然后将 中第 个子变量 与第 次迭代中的第 个子变量 的函数值进行比较，如果 的函数值比 的大，则将 作为第 次迭代的 。
[0101] 本发明所带来的有益技术效果：通过对低渗透裂缝性油藏进行简化，构建连通管道的双重介质模型，可以在保证模拟精度的前提下对油藏进行快速建模，本发明还提供了自动历史拟合和自动生产优化模型的建立，使决策更加高效，有助于更好地了解低渗透油藏的地质特征和裂缝性质，提高效率、降低成本，从而增加油田开发的经济效益。
附图说明
[0102] 图1为本发明适用于低渗透裂缝性油藏的快速模拟及优化方法的流程图。
[0103] 图2为本发明实施例历史拟合时第一次迭代的日产油曲线对比图。
[0104] 图3为本发明实施例历史拟合时第二次迭代的日产油曲线对比图。
[0105] 图4为本发明实施例历史拟合时第三次迭代的日产油曲线对比图。
[0106] 图5为本发明实施例历史拟合后的生产井1的日产油曲线与数值模拟器日产油曲线的对比图。
[0107] 图6为本发明实施例四个生产井的生产优化制度图。
[0108] 图7为本发明实施例的生产优化结果对比图。
具体实施方式
[0109] 下面结合附图以及具体实施方式对本发明作进一步详细说明：
[0110] 本发明提出了一种适用于低渗透裂缝性油藏的快速模拟及优化方法，主要包含如下内容：
[0111] 本发明以油藏系统中井点相对坐标为基础，通过井网剖分的限制条件将油藏划分为井井相连的连通管道，并在管道中均匀划分网格，将划分好的网格矩阵抽象为基质网格和裂缝网格，每个网格具有自己的地质参数，在基质系统中考虑低渗透裂缝性油藏的非达西渗流，在裂缝系统中考虑岩石的应力敏感效应，结合物质平衡方程，从而构建双重介质模型，对每个网格的压力及含水饱和度等参数进行全隐式求解；
[0112] 本发明利用ES‑MDA历史拟合算法，结合约束条件，建立自动历史拟合数学模型。构建历史拟合目标函数，通过对模型参数进行调整，使模型模拟结果可以与实际观测到的油藏生产动态相吻合，得到能够指导后续生产优化的油藏模型，在历史拟合过程中，为降低模型的随机性，保证模型与实际地质认识相符，对模型参数增加约束，建立自动历史拟合数学模型。利用ES‑MDA历史拟合算法进行求解，并对模型参数向量进行更新，直至完成整个迭代过程或者拟合精度满足工程需求；
[0113] 本发明将经济净现值作为生产优化的目标函数，通过对油藏地质储量、各井压力以及注采液量上下限进行约束，建立油藏注采优化模型，采用差分进化算法，通过选择、交叉、变异得到使油藏经济净现值最大化的生产调控制度。
[0114] 如图1所示，一种适用于低渗透裂缝性油藏的快速模拟及优化方法，具体包括以下步骤：
[0115] 步骤1、将油藏系统抽象为基质网格和裂缝网格，根据基质和裂缝的渗流特征构建双重介质模型。
[0116] 以油藏系统中相对井点坐标为基础，通过井网剖分的限制条件将油藏划分为井井相连的连通管道，并在管道中均匀划分网格，将划分好的网格矩阵抽象为基质网格和裂缝网格，每个网格具有自己的地质参数，在基质系统中考虑低渗透裂缝性油藏的非达西渗流，在裂缝系统中考虑岩石的应力敏感效应，结合物质平衡方程，从而构建双重介质模型，对每个网格的压力及含水饱和度等参数进行全隐式求解。
[0117] 步骤1的具体过程如下：
[0118] 步骤1.1、按照油藏系统中井点的分布将储层离散为一系列的节点，对所有节点两两相连形成初始井网模型；所述初始井网模型，即将油藏中所有井点进行两两相连所得到的井网模型；
[0119] 步骤1.2、并不是所有连接对应的两口井之间都有流体交换，因此需要对初始井网模型进行优化筛选，因此通过引入判定条件对初始井网模型进行筛选，从而生成一系列符合条件的两两相连的井间连接管道，并将储层体积划分到每条连接管道中，每条连接管道的体积取其井点与两侧连通三角形重心组成四边形的面积乘以有效厚度。
[0120] 其中，引入的判定条件包括距离条件、角度条件和最大连接点数条件；
[0121] 每口井都有自己作用范围，注水井不可能作用到距离很远的生产井，因此需要对初始井网模型进行距离筛选，距离条件的阈值设定为：一般取值为最小井距的2至3倍，超过该阈值的井视为不连接。
[0122] 所述角度条件具体设置为：任意三条连接管道组成的三角形中，其最小内角不能小于所设定的最小角度，否则视该三口井是在同一条路径中，注水井通过中间井流向另一口，而不是直接见效。
[0123] 所述最大连接点数条件具体设置为：每个井点允许连接的井点数量为4‑8个。
[0124] 步骤1.3、在每条连接管道上划分网格，分为基质网格和裂缝网格，所有的基质网格构成基质系统，所有的裂缝网格构成裂缝系统；每个网格都有自己的渗透率、孔隙度、体积、压力及含水饱和度等参数，其中基质网格为主要的储集单元，裂缝网格为主要的流动单元，网格间流体流动能力由网格间传导率大小决定，基质与基质间、裂缝与裂缝间、基质与裂缝间都有流体流动，基质与裂缝间的流体流量取决于基质系统、裂缝系统两个系统的压力，从而将初始井网模型的管道连接转换为双重介质模型中网格与网格间的连接。
[0125] 所述划分网格是在每条连接管道中平均分为5‑10个网格，井点所在网格体积可以指定（比如 ），如果指定井点网格体积，则网格体积默认值计算如下：
[0126] (1)；
[0127] 其中， 为井点 所在网格体积； 为总连接体积； 为连接管道数量；为每根连接管道中网格数量； 为井点数量。公式（1）的含义为井点网格的体积=油藏总体积/（网格数量+井点数量）。
[0128] 确定好井点网格体积后，连接管道中网格体积的计算公式如下：
[0129] (2)；
[0130] 其中，、为位于连接管道两端的不同井点； 表示连接管道中第 个网格；表示 、 井点的连接管道中第 个网格的体积； 为 、 井点间连接管道的体积； 和分别为井点 和井点 所在网格体积； 为井点 连接的管道数量； 为井点
连接的管道数量。
[0131] 所述基质网格和裂缝网格中，裂缝网格体积 、基质网格体积 计算如下：
[0132] (3)；
[0133] (4)；
[0134] (5)；
[0135] (6)；
[0136] 其中， 表示粗网格内裂缝间距单元数； 表示微裂缝占整个单元的体积分数，1； 为油藏中微裂缝间距，m； 为油藏系统中微裂缝开度，m。
[0137] 所述每个网格都有自己的渗透率、孔隙度、压力及含水饱和度等参数，其值由两井点参数值差值取得，以孔隙度为例，比如两井点的孔隙度分别为 和 ，在其间划分 个网格，则第 个网格的孔隙度为 ，其余参数以此类推。
[0138] 所述网格间传导率计算方法如下：
[0139] (7)；2
[0140] 其中， 为网格 和网格 间传导率； 为网格 和网格 的接触面积，m；为网格 的中心到网格 和网格 的接触面的距离； 为网格 的中心到网格 和网格的接触面的距离； 为网格 拟渗透率 和网格 拟渗透率 的平均值，其拟渗透率在步骤1.4中提到。
[0141] 所述基质与裂缝间的流体流动取决于基质系统、裂缝系统两个系统的压力，基质与裂缝间流体的窜流计算如下：
[0142] (8)；2
[0143] 其中， 为基质裂缝间的窜流量，kg/s； 为形状因子，1/m； 为流体密度，kg/3
m； 为基质网格的渗透率， ； 为流体粘度， ； 为基质系统压力， ； 为
3
裂缝系统压力， ； 为基质网格体积，m。
[0144] 步骤1.4、低渗透裂缝性油藏需要在基质系统中考虑非达西渗流，在裂缝系统中考虑岩石的应力敏感效应；
[0145] 考虑非达西渗流后，基质系统的运动方程表示如下：
[0146] (9)；
[0147] 其中， 为基质渗流速度， ； 为基质网格的渗透率， ； 为流体粘度，； 为流动压差， ； 和 为两个不同的非线性渗流系数。
[0148] 考虑岩石的应力敏感效应后，裂缝系统的运动方程为：
[0149] (10)；
[0150] 其中， 为裂缝渗流速度， ； 为裂缝网格的渗透率， ； 为常数，； 为流体粘度， ； 为应力敏感系数， ； 为裂缝网格压力，； 为裂缝网格初始压力， ； 为流动压差， 。
[0151] 将公式（9）和公式（10）统一表述为如下运动方程的形式：
[0152] (11)；
[0153] 其中， 为某相流体； ，表示 为油相 或水相 ； 为 的渗流速度， ； 为 的拟渗透率， ； 为 的流体粘度， ； 为流动压差，
。
[0154] 考虑非达西渗流和应力敏感效应后，渗透率可视为压力的函数，那么基质和裂缝的拟渗透率表示如下：
[0155] (12)；
[0156] (13)；
[0157] 其中， 和 分别表示基质和裂缝； 为网格 和网格 间考虑非达西渗流后基质 的拟渗透率， ； 分别为网格 和网格 间考虑应力敏感后裂缝 内的拟渗透率， ； 为网格 和网格 间基质 的初始渗透率， ； 表示网格 和网格 间裂缝 的初始渗透率， ； 为流体的粘度， ； 为基质内网格 和
网格 两节点的压力梯度， ； 和 为两个不同的非线性渗流系数，其值由实验测定； 为应力敏感系数， ； 为网格 和网格 间裂缝的当前压力， ；
为网格 和网格 间裂缝的初始压力， 。
[0158] 通过这种处理，能够使基质系统和裂缝系统的运动方程呈现出一致的形式，这将为进一步的数学模型构建以及程序编写工作带来便利。
[0159] 步骤1.5、结合基质和裂缝内流体的渗流规律，根据公式（11）构建低渗透裂缝性油藏双重介质模型。
[0160] 构建的双重介质模型应符合以下假设：渗流过程是等温的，油藏内部的渗流行为不受温度变化影响；模型考虑了同时存在的油相和水相，两相不相混溶，且都遵循相同的渗透规律；可考虑毛管力的影响，暂未考虑重力的影响；地层中流体和岩石均具有可压缩性。
[0161] 裂缝系统和基质系统为抽象出来的两套网格系统，其在空间上同时存在、相互叠加，且基质系统内流体渗流为非线性渗流，裂缝系统内流体流动为达西渗流，但考虑应力敏感效应。
[0162] 双重介质模型内包括基质系统的连续性方程和裂缝系统的连续性方程；
[0163] 其中，基质系统的连续性方程如下：
[0164] (14)；
[0165] (15)；
[0166] 其中， 为油相密度， ； 为基质的孔隙度，1； 为基质中油相饱和度，1； 为基质网格体积， ；为时间，s； 为油相在基质的渗流速度， ；表示散度； 为形状因子， ； 为基质内的拟渗透率， ； 为油相粘度， ；
为基质内油相压力， ； 为裂缝内油相压力， ； 为基质内油相的注采速
率，注入为正，采出为负， ； 为水相密度， ； 为水相在基质内的饱和
度，1； 为水相在基质的渗流速度， ； 为水相粘度， ； 为基质内水相
压力， ； 裂缝内水相压力， ； 为基质内水相的注采速率，注入为正，采出为负， ；
[0167] 裂缝系统的连续性方程如下：
[0168] (16)；
[0169] (17)；
[0170] 其中， 为裂缝的孔隙度，1； 为油相在裂缝内的饱和度，1； 为水相在裂缝内的饱和度，1； 为裂缝网格体积， ； 为油相在裂缝内的渗流速度， ；为水相在裂缝内的渗流速度， ； 为裂缝内油相的注采速率，注入为正，采出为负，； 为裂缝内水相的注采速率，注入为正，采出为负， ；
[0171] 令 表示基质网格内单位时间单位体积内相流体 体积流量的变化值，表示裂缝网格内单位时间单位体积内相流体 体积流量的变化值，即， ，其中，
为相流体； ，表示 为油相 或水相 ； 为基质内相流体 的注采速率，
注入为正，采出为负， ； 为相流体 的密度， ； 为基质内相流体 的
压力， ； 为裂缝内相流体 的压力， ； 为裂缝内相流体 的注采速率，注入为正，采出为负， ；此时，基质系统的连续性方程可变为：
[0172] (18)；
[0173] (19)；
[0174] 其中， 为基质网格内单位时间单位体积内油相流体体积流量的变化值，； 为基质网格内单位时间单位体积内水相流体体积流量的变化值， 。
[0175] 裂缝系统的连续性方程可变为：
[0176] (20)；
[0177] (21)；
[0178] 其中， 为裂缝网格内单位时间单位体积内油相流体体积流量的变化值，； 为裂缝网格单位时间单位体积内水相流体体积流量的变化值， 。
[0179] 步骤1.6、对所建立的连续性方程进行线性插值近似，连续性方程的时间项采用向后一阶差分，空间项采用中心差分格式，将数值离散结果写为残差的形式，利用Newton‑Raphson方法（牛顿‑拉裴森方法）全隐式求解裂缝及基质系统中压力和含水饱和度等参数。具体过程如下：
[0180] 步骤1.6.1、定义流度系数如下：
[0181] (22)；
[0182] 其中， 为某相流体； ，表示 为油相 或水相 ； 为某相流体流度， ； 为某相流体的相对渗透率， ； 为某相流体粘度， 。
[0183] 步骤1.6.2、因为基质系统和裂缝系统的渗流特征和互相之间的窜流可以综合考虑在其网格间的传导率中，这样，基质网格与裂缝网格除本身的参数不同外，没有本质的区别，因此，可以将两套系统用统一的格式进行线性插值近似，这样可以大大降低后续程序编制的复杂性。
[0184] 线性插值近似采用如下的数值离散方程：
[0185] (23)；
[0186] 式中， 为某相流体； ，表示 为油相 或水相 ； 为某相流体的密度， ； 为孔隙度； 为某相流体的饱和度； 为网格体积， ； 为模拟迭代过程中的上一时刻； 为模拟迭代过程中的当前时刻； 为时间步长，； 为网格 相邻的单元集合，在该模型中，以基质网格为例，其相邻的单元只有旁边的同一管道内的网格以及其对应的裂缝网格；下标 表示两个相邻网格 和 界面上某一参数的平均值；为网格 和网格 间传导率； 为某相流体在 时刻网格 的压力， ； 为
某相流体在 时刻网格 的压力， ； 为网格 内某相流体在 时刻的注采
速率，注入为正，采出为负， 。
[0187] 步骤1.6.3、将上述数值离散方程写为如下的残差形式：
[0188] (24)；
[0189] 式中， 为方程残差； ， 为网格总数，此时间步选项为 ，因此是全隐式求解，利用Newton‑Paphson迭代求解方法可以对裂缝及基质内的压力和含水饱和度进行求解。
[0190] 所述Newton‑Paphson迭代求解方法，即从初始点开始，通过不断的迭代计算，得到一个更接近函数零点的新点，这个新点是通过在当前点处进行泰勒展开计算得到的，然后将新点作为下次迭代的初始点继续进行迭代计算，直到找到满足要求的解。
[0191] 步骤2、利用ES‑MDA（集合平滑多重数据同化）历史拟合算法，结合约束条件，建立自动历史拟合数学模型。
[0192] 构建历史拟合目标函数，通过对模型参数进行调整，使模型模拟结果可以与实际观测到的油藏生产动态相吻合，得到能够指导后续生产优化的油藏模型，在历史拟合过程中，为降低模型的随机性，保证模型与实际地质认识相符，对模型参数增加约束，建立自动历史拟合数学模型。利用ES‑MDA历史拟合算法进行求解，并对模型参数向量进行更新，直至完成整个迭代过程或者拟合精度满足工程需求。
[0193] 所述步骤2具体包括以下步骤：
[0194] 步骤2.1、自动历史拟合数学模型目标函数的构建在本质上其实是最小化模型计算结果与实际观测数据的误差，因此，可以表示为：
[0195] (25)；
[0196] 式中， 为目标函数值，即计算结果与观测数据的差异； 表示模型计算结果，维度为 ，其中 表示观测数据的数量； 表示实际观测数据，维度为； 表示观测误差的协方差矩阵，维度为 。
[0197] 步骤2.2、该自动历史拟合数学模型通过简化方法大幅减少了油藏模型的特征参数，因此在历史拟合过程中，仅需要调整每根管道的体积、网格间传导率以及幂律相渗曲线参数。
[0198] 幂律相渗曲线参数由如下的幂律相渗模型调整：
[0199] (26)；
[0200] 其中，
[0201] (27)；
[0202] 式中， 和 分别表示油相和水相的相对渗透率； 表示无因次饱和度； 和分别决定水相和油相的相对渗透率曲线的弯曲程度； 表示束缚水饱和度； 表示残余油饱和度； 表示在含水饱和度为束缚水饱和度时的水相相对渗透率值； 为含水饱和度值。这些参数和值的变化关系共同描述了多相流体在多孔介质中渗透特性的变化情况。
[0203] 步骤2.3、进行历史拟合调整的参数 如下：
[0204] (28)；
[0205] 其中， 为网格 和网格 间传导率，包括基质与基质网格间、裂缝与裂缝网格间、基质与裂缝网格间； 为 、 井点的连通管道的基质传导体积， 为 、 井点的连通管道的裂缝传导体积； 、 、 为 、 井点的连通管道相渗曲线的不同参数；和 为两个不同的非线性渗流系数；
[0206] 当所有井对间共用一组相渗数据时，历史拟合参数 表示如下：
[0207] (29)；
[0208] 其中， 、 、 为所有井对间共用的三个相渗曲线参数。
[0209] 步骤2.4、在历史拟合过程中，对模型参数增加约束；自动历史拟合数学模型的约束条件如下：
[0210] 步骤2.4、在历史拟合过程中，为降低模型的随机性，保证模型与实际地质认识相符，需要对模型参数增加约束，从而建立自动历史拟合数学模型。
[0211] 根据传导率和体积的定义可知，传导率和体积的值是非负的，并且每根管道的体积要小于总传导体积，并且，所有管道体积之和应该等于总传导体积，相渗模型中的参数也应该服从以下约束，非线性渗流曲线的特征参数应该大于0，综合压缩系数取值一般在0.00001 0.1之间，因此，自动历史拟合数学模型的约束条件如下：
~
[0212] (30)；
[0213] 其中， 为网格 和网格 的基质体积； 为基质系统的总体积； 为网格 和网格 的裂缝体积； 为裂缝系统的总体积； 为总的井数量。
[0214] 步骤3、将经济净现值作为生产优化的目标函数，通过对油藏地质储量、各井压力以及注采液量上下限进行约束，建立油藏注采优化模型，采用差分进化算法得到使油藏经济净现值最大化的生产调控制度。
[0215] 所述步骤3具体包括以下步骤：
[0216] 步骤3.1、采用经济净现值作为生产优化的目标函数，数学表达式如下：
[0217] (31)；
[0218] 其中， 为经济净现值； 为控制变量；为时间步序号， 表示总时间步数； 表示生产井序号， 表示生产井数量； 表示注水井序号， 表示注水井数量； 为每吨原油的价格； 为生产井 在 时间步的采油量； 为处理每吨污水的成本；表示生产井 在 时间步的采水量； 表示每吨水的注水成本； 表示注水井 在 时间步的注水量；
[0219] 步骤3.2、在生产优化过程中，采用边界约束条件对油藏参数进行控制，油藏参数的约束条件如下所示：
[0220] (33)；
[0221] 其中， 为注水井 在 时间步的注水量； 为生产井 在 时间步的采液量；为注采井 在 时间步注采液量 的下限，注采井包括注水井和生产井； 为注采井 在 时间步注采液量 的上限； 为总的地质储量； 和 分别为注采井
控制压力 的下限和上限；
[0222] 步骤3.3、采用差分进化算法得到使油藏经济净现值最大化的生产调控制度。
[0223] 所述差分进化算法，即首先初始化形成一系列的控制变量 ，对这一系列控制变量通过以下方式进行变异，变异后的控制变量记为 ；
[0224] (34)；
[0225] 式中， 代表第 次迭代； 表示变异算子，可取0‑2之间的任意数； 代表三个不同的个体序号； 表示第 次迭代变异后的序号为 的控制变量； 表示第 次迭代的序号为 的初始控制变量； 表示第 次迭代的序号为 的初始控制变量； 表示第 次迭代的序号为 的初始控制变量；也就是说在初始控制变量中选择三个个体进行上述操作，得到下一次迭代的序号为 的控制变量。
[0226] 接下来进行交叉，交叉后的变量记为 ，对于每个 ，有：
[0227] (35)；
[0228] 其中， 表示0‑1的随机数； 表示交叉算子； 表示第 次迭代交叉后序号为 的控制变量的部分片段； 表示第 次迭代变异后序号为
的控制变量的部分片段； 表示第 次迭代序号为 的初始控制变量的部分片段；如果随机数大于交叉算子，则交叉后的变量取变异后的 ，如果随机数小于交叉算子，则取原先 次迭代的初始种群。
[0229] 经过变异交叉后的变量 ，对 中的每个子变量进行求目标函数操作，然后将 中第 个子变量 与第 次迭代中的第 个子变量 的函数值进行比较，如果 的函数值比 的大，则将 作为第 次迭代的 。
[0230] 为了证明本发明的可行性与优越性，给出如下实施例。
[0231] 首先，利用有限体积法油藏数值模拟器建立一个双重介质模型，其网格规模：25*25*1，模型大小：100m*100m*10m，油藏基质孔隙度为0.2；裂缝孔隙度为0.02，微裂缝开度为‑5
1*10 m，为裂缝间距为5m，油藏的初始含水饱和度为0.2，油水粘度分别为4 mPa*s、1 mPa*s，油藏中有1口注水井，4口生产井，注水井在油藏中央，四口生产井在油藏四个角点，采用
3 3
油藏数值模拟器模拟生产1500天，注水井注水速率为0.8m/d，生产井产液速率均为0.2m /d，区块整体注采平衡。
[0232] 然后使用本发明所述方法，使用相同的油藏参数建模，其单次正向运行时间为15s，相比于有限体积油藏数值模拟器耗时70.3s，计算速度提升接近5倍，并采用ESMDA算法进行历史拟合，种群数量选100个，迭代5次，其历史拟合中生产井1的产油速度变化如图2、图3、图4所示，图2、图3、图4分别为历史拟合第一次迭代、第二次迭代、第三次迭代时的产油速度对比图，其中横坐标为时间，纵坐标为日产油速度，灰色细曲线为每个控制变量得到的日产油曲线，反映每个控制变量的日产油速度变化，黑色粗曲线为油藏数值模拟器的日产油曲线，反映数值模拟器的日产油速度变化。图5为本发明经历史拟合后得到的生产井1的日产油曲线与常规油藏数值模拟器计算得到的产油速度对比曲线，可以看出，经过历史拟合后，生产井1的日产油曲线与数值模拟器的模拟真实日产油曲线基本相符。最后，利用历史拟合后的模型进行生产优化，每个生产井的生产优化制度均包括初始调控和第一、二、三次调控，调控后的产液速度如图6所示。本实施例优化三个时间步，每个时间步为100天，总共优化300天，优化结果包括区块累产油和区块含水率，具体数据如图7所示，优化后的累产油明显提高，含水率相对稳定，因此验证了本发明方法在低渗透裂缝性油藏模拟及优化计算的可靠性。
[0233] 当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。
