一种油藏模拟辅助历史拟合方法、系统、存储介质及设备与流程

本发明涉及油藏开采技术领域，尤其涉及一种油藏模拟辅助历史拟合方法、系统、存储介质及设备。

背景技术：

油藏数值模拟历史拟合是油藏数值模拟中的关键环节，是正确预测油藏未来开发动态和制定油田最优生产开发方案的保证。历史拟合是典型的油藏反问题，通过反复调整油藏模型参数，使模型的计算结果与历史动态生产数据吻合或接近。传统人工历史拟合速度慢，准确性低，依赖于油藏工作者对油藏的分析认识和拟合经验，有一定盲目性和随意性，无法满足油田开发生产高效化、实时化的要求。

人工历史拟合是基于已有的油藏认知或以往经验划定一个不精确的区域。在辅助历史拟合研究中，研究者通常基于先验地质信息(如沉积相态)并结合距离截断方法得到单井敏感性区域，但距离截断存在一定的人为性，不够准确。另外，历史拟合需要反演的油藏模型参数通常数以十万计，采用如模拟退火算法、遗传算法、随机扰动近似梯度算法、集合卡尔曼滤波等主流无梯度优化算法直接计算，计算代价大。

技术实现要素：

本发明所要解决的技术问题是针对上述现有技术的不足，提供一种油藏模拟辅助历史拟合方法、系统、存储介质及设备。

本发明解决上述技术问题的技术方案如下：一种油藏模拟辅助历史拟合方法，包括如下步骤：

步骤1：采用正向最大匹配fmm算法追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域；步骤2：构建历史拟合目标函数，以所述单井敏感性区域内预设油藏模型的模型参数平均值作为历史拟合目标函数的参数变量，对所述历史拟合目标函数进行优化求解，并将所述历史拟合目标函数收敛到最小值时对应的自变量的值作为油藏模拟辅助历史拟合中模型参数的最优估计输出。

本发明的有益效果是：本发明的基于单井敏感性区域的油藏模拟辅助历史拟合方法，通过追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域，再建立历史拟合目标函数，并将预设油藏模型的模型参数进行参数化处理，将单井敏感性区域内模型参数的平均值作为历史拟合目标函数中的参数变量，并通过调整预设油藏模型的模型参数，使历史拟合效果达到最好，最大程度上提高模型准确度，相较于人工历史拟合，大大降低了变量的维数，提高了计算效率。

在上述技术方案的基础上，本发明还可以做如下改进：

进一步：所述步骤1中，所述追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间具体包括：

步骤11：将所述预设油藏模型中位于井点处的网格标记为冻结点，并作为初始起点；

步骤12：将与所述冻结点相邻的所有网格标记为相邻点，并计算压力波从所述冻结点到对应的每一个所述相邻点的飞行时间，并将所述飞行时间的最小值对应的相邻点标记为新的冻结点；

步骤13：将所述新的冻结点作为初始起点，并重复上述步骤12，直到预设油藏模型的所有网格均被标记为冻结点。

上述进一步方案的有益效果是：通过追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间，可以遍历所述预设油藏模型的所有网格，从而可以确定所述预设油藏模型中每一网格对应的飞行时间，便于后续根据所述飞行时间准确的确定对应的单井敏感性区域。

进一步：所述步骤12中，所述计算压力波从所述冻结点到对应的每一个所述相邻点的飞行时间具体包括：

将井点处油藏压力波传播方程表示为程函方程：

式中：f为速度方程，为飞行时间τ相对空间位置x的偏导数；

压力波穿过预设油藏模型中每个网格的速度定义为：

式中：k为渗透率，μ为流体粘度，为常数，ct为油藏综合压缩系数，为常数，φ为孔隙度；

考虑油藏的非均质性，基于逆风差分近似求解程函方程，计算飞行时间的方程表示为

式中：fi、fj、fk为压力波在预设油藏模型网格的x、y、z方向上的速度，δx、δy、δz为预设油藏模型网格在x、y、z方向的网格尺寸，τ、τ1、τ2、τ3的表达式分别为：

τ＝τi,j,k(4)

τ1＝min(τi-1,j,k,τi+1,j,k)(5)

τ2＝min(τi,j-1,k,τi,j+1,k)(6)

τ3＝min(τi,j,k-1,τi,j,k+1)(7)

式中：τi,j,k为预设油藏模型的任一网格(i,j,k)对应的飞行时间。

上述进一步方案的有益效果是：根据上述公式可以准确的计算出压力波从所述冻结点到对应的每一个所述相邻点的飞行时间，这样既可确定所述预设油藏模型中每一网格对应的飞行时间，便于后续根据每一网格对应所有井的最小飞行时间准确的确定对应的单井敏感性区域，提高历史拟合效果。

进一步：所述步骤1中，所述根据所述飞行时间划定三维空间中每口井的单井敏感性区域的具体实现为：

当油藏存在多口井时，获取预设油藏模型中每一个网格对应每一口井的飞行时间的最小值，并将该网格划定为所述最小值对应的井的敏感性区域。

上述进一步方案的有益效果是：根据预设油藏模型中每一个网格对应每一口井的飞行时间的最小值将该网格划定为所述最小值对应的井的敏感性区域，从而精确的得到每口井对应的敏感性区域，从而根据单井敏感性区域得到更加精确的历史拟合结果。

进一步：所述步骤2的具体实现为：

步骤21：构建历史拟合目标函数，获取预设油藏模型中每一层所述单井敏感性区域的模型参数，并计算预设油藏模型中每一层所述模型参数的初始平均值；

其中，所述历史拟合目标函数定义如下：

式中：m为预设油藏模型参数组成的nm维向量；dobs为nd维向量，其包含实际测量的动态生产数据；g为油藏系统，cd为油藏历史动态生产数据测量误差的协方差矩阵；

步骤22：对所述预设油藏模型计算的生产数据预测值与历史动态生产数据进行拟合，将模型参数的初始平均值作为所述历史拟合目标函数的参数变量计算历史拟合目标函数初始值，并采用近似扰动梯度升级算法调整所述预设油藏模型中每一层模型参数的平均值，获取所述历史拟合目标函数的最小值；

步骤23：将所述历史拟合目标函数收敛到最小值时对应的自变量的值作为油藏模拟辅助历史拟合中模型参数的最优估计输出。

上述进一步方案的有益效果是：通过将预设油藏模型中每一层所述单井敏感性区域的模型参数的平均值作为历史拟合目标函数的变量，可以达到降低了变量的维数，简化计算，提高了计算效率，并且提高了预设油藏模型的精确度。

进一步：所述步骤22中，采用近似扰动梯度升级算法调整所述预设油藏模型中每一层模型参数的平均值，获取所述历史拟合目标函数的最小值具体为：

对于目标函数o(m)，考虑在第l个迭代步的模型参数的平均值为ml，在其周围进行扰动生成n个模型变量：

ml,i＝ml+γδl,i,i＝1,2,…,n(13)

式中：初始状态时，m0取所述模型参数的初始平均值；ml,i为第l步第i个模型扰动变量；δl,i为扰动向量；γ为扰动步长；

构造近似扰动梯度一般式为：

式中：δ为n个扰动向量组成的矩阵；δo为n维列向量，其中的元素为模型变量实现的目标函数值与当前最优目标函数值的差值；c为常数；l为优化后的下三角矩阵；

根据拟合前后模型参数的平均值对所述预设油藏模型进行更新，具体表达式如下：

式中，αl为迭代步长；

当相邻两次迭代对应的所述历史拟合目标函数值之间的差值为零，则将对应的所述历史拟合目标函数值作为其最小值。

上述进一步方案的有益效果是：通过对所述最优模型变量进行扰动，并生成模型扰动变量，并根据拟合前后模型变量和近似扰动梯度对所述预设油藏模型进行更新，直到得到所述历史拟合目标函数的最小值，从而得到油藏模拟辅助历史拟合中模型参数的最优估计并输出。

进一步：所述下三角矩阵l的优化过程如下：

近似梯度与真实梯度g的向量积为：

近似扰动梯度和真实梯度g夹角的余弦值为

下三角矩阵l的优化目标函数为：

采用优化算法对式(17)进行迭代，对所述下三角矩阵l进行优化，得到优化后的下三角矩阵l。

上述进一步方案的有益效果是：通过对所述下三角矩阵l进行优化，可以使得所述近似梯度更接近真实梯度，从而使得历史拟合的结果更加精确。

本发明还提供了一种基于单井敏感性区域的油藏模拟辅助历史拟合系统，跟踪划分模块和计算拟合模块；

所述跟踪划定模块，用于采用正向最大匹配fmm算法追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域；

所述计算拟合模块，用于构建历史拟合目标函数，以所述单井敏感性区域内预设油藏模型的模型参数平均值作为历史拟合目标函数的参数变量对所述历史拟合目标函数进行优化求解，并将所述历史拟合目标函数收敛到最小值时对应的自变量的值作为油藏模拟辅助历史拟合中模型参数的最优估计输出。

本发明的基于单井敏感性区域的油藏模拟辅助历史拟合系统，通过跟踪划定模块追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域，再通过计算拟合模块建立历史拟合目标函数，并将预设油藏模型的模型参数进行参数化处理，将单井敏感性区域内模型参数的平均值作为历史拟合目标函数中的参数变量，并通过调整预设油藏模型的模型参数，使历史拟合效果达到最好，最大程度上提高模型准确度，相较于人工历史拟合，大大降低了变量的维数，提高了计算效率。

本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，实现所述的基于单井敏感性区域的油藏模拟辅助历史拟合方法。

本发明还提供了一种基于单井敏感性区域的油藏模拟辅助历史拟合设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述的基于单井敏感性区域的油藏模拟辅助历史拟合方法的步骤。

附图说明

图1为本发明的基于单井敏感性区域的油藏模拟辅助历史拟合方法流程示意图；

图2为本发明的飞行时间追踪计算过程示意图

图3为本发明一实施例历史拟合目标函数的优化过程示意图；

图4为本发明一实施例的区块累产油量示意图；

图5为本发明一实施例的一井日产油量示意图；

图6为本发明一实施例的另一井日产油量示意图；

图7为本发明一实施例的又一井日产油量示意图；

图8为本发明的基于单井敏感性区域的油藏模拟辅助历史拟合系统结构示意图。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

如图1所示，一种基于单井敏感性区域的油藏模拟辅助历史拟合方法，包括如下步骤：

步骤1：采用正向最大匹配fmm算法追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域；

步骤2：构建历史拟合目标函数，以所述单井敏感性区域内预设油藏模型的模型参数平均值作为历史拟合目标函数的参数变量对所述历史拟合目标函数进行优化求解，并将所述历史拟合目标函数收敛到最小值时对应的自变量的值作为油藏模拟辅助历史拟合中模型参数的最优估计输出。

本发明的基于单井敏感性区域的油藏模拟辅助历史拟合方法，通过追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域，再建立历史拟合目标函数，并将预设油藏模型的模型参数进行参数化处理，将单井敏感性区域内模型参数的平均值作为历史拟合目标函数中的参数变量，并通过调整预设油藏模型的模型参数，使历史拟合效果达到最好，最大程度上提高模型准确度，相较于人工历史拟合，大大降低了变量的维数，提高了计算效率。

本发明的实施例中，所述步骤1中，所述追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间具体包括：

步骤11：将所述预设油藏模型中位于井点处的网格标记为冻结点，并作为初始起点；

步骤12：将与所述冻结点相邻的所有网格标记为相邻点，并计算压力波从所述冻结点到对应的每一个所述相邻点的飞行时间，并将所述飞行时间的最小值对应的相邻点标记为新的冻结点；

步骤13：将所述新的冻结点作为初始起点，并重复上述步骤12，直到预设油藏模型的所有网格均被标记为冻结点。

通过追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间，可以遍历所述预设油藏模型的所有网格，从而可以确定所述预设油藏模型中每一网格对应的飞行时间，便于后续根据所述飞行时间准确的确定对应的单井敏感性区域。

如图2所示，将某一确定网格作为初始起点(τ＝0)，并标记为冻结点，如图2(a)中网格13；查找与冻结点相连接的网格，标记为相邻点，如图2(b)中网格8、12、14、18；利用式(3)计算压力波从冻结点到达相邻点的飞行时间，选出飞行时间最小的相邻点，标记为冻结点，如图2(c)中网格12；将网格12作为新的初始起点，并找到网格12的相邻点(不包括已经计算过飞行时间的网格)，如图2(d)中网格7、11、17；计算压力波从网格12到达网格7、11、17的飞行时间，并找到所有相邻点中飞行时间最小的点，标记为冻结点，如图2(e)中网格18，如此重复，直到所有的网格被标记为冻结点。

需要指出的是，在计算整个油藏网格的飞行时间过程中，将井点处网格作为初始起点，考虑的是离井点越近，飞行时间越小；离井点越远，飞行时间越大。

优选地，本发明的实施例中，所述步骤12中，所述计算压力波从所述冻结点到对应的每一个所述相邻点的飞行时间具体包括：

将井点处油藏压力波传播方程表示为程函方程：

式中：f为速度方程，为飞行时间τ相对空间位置x的偏导数；

压力波穿过预设油藏模型中每个网格的速度定义为：

式中：k为渗透率，μ为流体粘度，为常数，ct为油藏综合压缩系数，为常数，φ为孔隙度；

考虑油藏的非均质性，基于逆风差分近似求解程函方程，计算飞行时间的方程表示为

式中：fi、fj、fk为压力波在预设油藏模型网格的x、y、z方向上的速度，δx、δy、δz为预设油藏模型网格在x、y、z方向的网格尺寸，τ、τ1、τ2、τ3的表达式分别为：

τ＝τi,j,k(4)

τ1＝min(τi-1,j,k,τi+1,j,k)(5)

τ2＝min(τi,j-1,k,τi,j+1,k)(6)

τ3＝min(τi,j,k-1,τi,j,k+1)(7)

式中：τi,j,k为预设油藏模型的任一网格(i,j,k)对应的飞行时间。

根据上述公式可以准确的计算出压力波从所述冻结点到对应的每一个所述相邻点的飞行时间，这样既可确定所述预设油藏模型中每一网格对应的飞行时间，便于后续根据每一网格对应所有井的最小飞行时间准确的确定对应的单井敏感性区域，提高历史拟合效果。

优选地，本发明的实施例中，所述步骤1中，所述根据所述飞行时间划定三维空间中每口井的单井敏感性区域的具体实现为：

当油藏存在多口井时，获取预设油藏模型中每一个网格对应每一口井的飞行时间的最小值，并将该网格划定为所述最小值对应的井的敏感性区域。

根据预设油藏模型中每一个网格对应每一口井的飞行时间的最小值将该网格划定为所述最小值对应的井的敏感性区域，从而精确的得到每个井对应的敏感性区域，从而根据单井敏感性区域得到更加精确的历史拟合结果。

本发明的实施例中，所述步骤2的具体实现为：

步骤21：构建历史拟合目标函数，获取预设油藏模型中每一层所述单井敏感性区域的模型参数，并计算预设油藏模型中每一层所述模型参数的初始平均值；

其中，所述历史拟合目标函数定义如下：

式中：m为预设油藏模型参数组成的nm维向量；dobs为nd维向量，其包含实际测量的动态生产数据；g为油藏系统，cd为油藏历史动态生产数据测量误差的协方差矩阵；

步骤22：对所述预设油藏模型计算的生产数据预测值与历史动态生产数据进行拟合，将模型参数的初始平均值作为所述历史拟合目标函数的参数变量计算历史拟合目标函数值，并采用近似扰动梯度升级算法调整所述预设油藏模型中每一层模型参数的平均值，获取所述历史拟合目标函数的最小值；

步骤23：将所述历史拟合目标函数收敛到最小值时对应的自变量的值作为油藏模拟辅助历史拟合中模型参数的最优估计输出。

通过将预设油藏模型中每一层所述单井敏感性区域的模型参数的平均值作为历史拟合目标函数的变量，可以达到降低了变量的维数，简化计算，提高了计算效率，并且提高了预设油藏模型的精确度。

由于预设油藏模型的模型参数在每一层存在差异，因此，模型参数的平均值是对每一层的敏感性区域进行计算的。最终变量m的维数仅与油藏井数、油藏模型的层数和反演参数的类型个数有关，与油藏模型的网格数量没有关系，可以大幅度的降低优化变量的维数，提高优化算法的计算效率。

以一个二维油藏模型为例，其含有4口井，模型参数的平均值计算公式如式(9)所示(以孔隙度为例)：

式中：φavei表示第i井敏感性区域内网格孔隙度的平均值；φi,j表示第i井敏感性区域内第j个网格的孔隙度；ni表示第i井敏感性区域内网格总数；nw表示总井数。

求得模型参数平均值后，历史拟合目标函数中变量m如式(10)所示，拟合计算后变量m'如式(11)所示，根据拟合前后的变量的倍数变化关系对油藏模型进行更新(以孔隙度为例)，如式(12)所示：

m＝[φave1,…,φave4,kave1,…,kave4]^t,i＝1,…,nw(10)

m'＝[φ′ave1,…,φ′ave4,k'ave1,…,k'ave4]^t,i＝1,…,nw(11)

式中：φavei表示优化前第i井敏感性区域内网格的平均孔隙度，φ′avei表示优化后第i井敏感性区域内网格的平均孔隙度，kavei表示优化前第i井敏感性区域内网格的平均渗透率，k'avei表示优化后第i井敏感性区域内网格的平均渗透率，φ′i，j表示第i井敏感性区域内第j个网格更新后的孔隙度。

优选地，在上述实施例的基础上，所述步骤22中，采用近似扰动梯度升级算法调整所述预设油藏模型中每一层模型参数的平均值，获取所述历史拟合目标函数的最小值具体为：

对于目标函数o(m)，考虑在第l个迭代步的模型变量为ml，在其周围进行扰动生成n个模型变量：

ml,i＝ml+γδl,i,i＝1,2,…,n(13)

式中：初始状态时，m0取所述模型参数的初始平均值；ml,i为第l步第i个模型扰动变量；δl,i为扰动向量；γ为扰动步长；

构造近似扰动梯度一般式为：

式中：δ为n个扰动向量组成的矩阵；δo为n维列向量，其中的元素为模型变量实现的目标函数值与当前最优目标函数值的差值；c为常数；l为优化后的下三角矩阵，其为单位阵时即是传统的spsa梯度；

根据拟合前后模型变量对所述预设油藏模型进行更新，具体表达式如下：

式中，αl为迭代步长；

当相邻两次迭代对应的所述历史拟合目标函数值之间的差值为零，则将对应的所述历史拟合目标函数值作为其最小值。

通过对所述最优模型变量进行扰动，并生成模型扰动变量，并根据拟合前后模型变量和近似扰动梯度对所述预设油藏模型进行更新，直到得到所述历史拟合目标函数的最小值，从而得到油藏模拟辅助历史拟合中模型参数的最优估计并输出。

具体地，所述下三角矩阵l的优化过程如下：

近似梯度与真实梯度g的向量积为：

近似扰动梯度和真实梯度g夹角的余弦值为

下三角矩阵l的优化目标函数为：

采用优化算法对式(17)进行迭代，对所述下三角矩阵l进行优化，得到优化后的下三角矩阵l。

通过对所述下三角矩阵l进行优化，可以使得所述近似梯度更接近真实梯度，从而使得历史拟合的结果更加精确。

使用本发明提出的历史拟合方法对某实际碳酸盐岩油藏进行了历史拟合测试，该油藏非均质性极强，油藏模型网格划分为76×115×20，总网格数为174800。共有油水井29口，选择对该油藏4564天生产历史的区块累产油量和单井日产油量进行拟合，反演的油藏模型参数渗透率和孔隙度，参数化后变量的维数为1160。

历史拟合目标函数优化结果如图3所示，在经过约35次迭代步后目标函数逐渐收敛。油藏的区块累产油量和单井日产油量拟合结果如图4-图7所示，拟合后模型的生产数据较于初始模型更接近实际值，拟合效果好。

本发明针对历史拟合中单井敏感性区域的确定，基于油藏中压力传播方程，采用fmm方法求解，提出了一种可定量表征油藏中单井敏感性区域的计算方法，该方法可以快速而准确的确定单井在三维空间中敏感性区域，具有稳定性好、计算效率高的优点，适用于大型实际油藏。

本发明针对实际油藏应用中历史拟合目标函数的参数变量维数大的问题，提出以油藏各层中单井敏感性区域内模型参数的平均值作为历史拟合目标函数中的变量，大幅度降低了变量的维数，提高了算法的计算效率。

实验证明，本发明的基于单井敏感性区域的油藏模拟辅助历史拟合方法，对油藏区块和单井动态生产指标的拟合取得了较好的效果，可靠性较高。

如图8所示，本发明还提供了一种基于单井敏感性区域的油藏模拟辅助历史拟合系统，跟踪划分模块和计算拟合模块；

所述跟踪划定模块，用于采用正向最大匹配fmm算法追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域；

所述计算拟合模块，用于构建历史拟合目标函数，以所述单井敏感性区域内预设油藏模型的模型参数平均值作为历史拟合目标函数的参数变量对所述历史拟合目标函数进行优化求解，并将所述历史拟合目标函数收敛到最小值时对应的自变量的值作为油藏模拟辅助历史拟合中模型参数的最优估计输出。

本发明的基于单井敏感性区域的油藏模拟辅助历史拟合系统，通过跟踪划定模块追踪压力波从井点穿过预设油藏模型中每个网格的飞行时间,并根据所述飞行时间划定三维空间中每口井的单井敏感性区域，再通过计算拟合模块建立历史拟合目标函数，并将预设油藏模型的模型参数进行参数化处理，将单井敏感性区域内模型参数的平均值作为历史拟合目标函数中的参数变量，并通过调整预设油藏模型的模型参数，使历史拟合效果达到最好，最大程度上提高模型准确度，相较于人工历史拟合，大大降低了变量的维数，提高了计算效率。

本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，实现所述的基于单井敏感性区域的油藏模拟辅助历史拟合方法。

本发明还提供了一种基于单井敏感性区域的油藏模拟辅助历史拟合设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述的基于单井敏感性区域的油藏模拟辅助历史拟合方法的步骤。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。