CN106355003A - 基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法及系统 - Google Patents

Abstract

本发明公开了一种基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法及系统，采用t分布得到初始的油藏静态参数，然后采用基于马尔科夫链蒙特卡洛的历史拟合法并调用油藏数值模拟软件对所述油藏静态参数进行迭代优化，使预测生产动态与真实值尽可能接近，得到优化的油藏静态参数和油藏数值模型。本发明先利用t分布得到初始的油藏静态参数，再采用基于马尔科夫链的蒙特卡洛法不断优化模型渗透率等油藏静态参数拟合生产实际动态，得到尽可能接近真实模型的油藏数值模型，减少拟合时间，提高历史拟合的效率和精度，使油田开发动态预测的结果更加接近实际生产。

Description

基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法及系统

技术领域

本发明涉及地球物理学中物探开发技术领域，具体涉及基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法及系统。

背景技术

在油藏数值模拟中，为了使动态预测能够尽量接近实际情况，通常需要对油藏数据进行历史拟合，根据所观测到的实际油藏动态来调整油藏模型参数，使得模型预测值与实际观测值的误差在允许范围内，为后续油藏开采服务。传统的历史拟合方法通过手工来调整模型参数，工作量大、繁琐，且效率低下。自动历史拟合方法采用优化算法自动调整油藏模型参数，大大缩短历史拟合时间，提高拟合精度。因此，研究快速的自动历史拟合方法是实现油藏历史拟合的急切需求。历史拟合问题是通过调整敏感参数(如孔隙度、渗透率等)，使得数值模拟计算的量如压力、油气比、含水率等都接近实际测量值，实质上是一个最优化问题。在求解历史拟合问题上，常见的有三种方法：梯度类方法、数据同化方法和随机类方法。

1、梯度类方法：梯度类算法中应用较多的是牛顿型方法。T.B.Tan和N.Kalogera设计了全隐式的三维三相的数值模拟器，将其应用在小型油藏模型中。但是Gauss-Newton方法不适合应用到大型的油藏模型中，因为Gauss-Newton方法在Hessian矩阵方面不便于计算。1995年孟雅杰在Gauss-Newton方法的基础上提出了改进的牛顿法，这个方法并不需要计算Hessian矩阵。2002年Razza和Reynolds又对此进行了修正，加入有限储存BFGS的策略，使得算法不再需要存储矩阵，仅仅需要计算前一步的梯度值和目标值。该方法解决了Gauss-Newton方法不适合于处理大规模的油藏历史拟合问题的弊端。2010年Razza和Reyonlds采用奇异值分解方法对算法参数进行降维，并将其应用到有限储存策略中，以此提出新的自动历史拟合的思路。梯度类算法是一种解决自动历史拟合问题的有效算法，然而因为该方法对伴随矩阵的计算的依赖性很高，并且其计算量大，不具备良好的可移植性。

2、数据同化方法：集合卡尔曼滤波方法(ENKF)是一种十分重要的数据同化方法，该方法最早主要被使用在气象学和海洋动力学中，因为集合卡尔曼滤波方法没有利用梯度类算法中伴随梯度的运算，所以在算法实现上较为方便，并且优化后的油藏模型能够体现真实油藏的不确定性。ENKF方法也存在同化循环中的滤波发散问题和在计算过程中不满秩的问题。

3、随机类算法：随机类算法是目前发展较快的一种算法，该类算法在计算过程中以随机概率和搜索策略来求解问题，它能够解决目标函数复杂和梯度求解困难的问题。2004年Tokuda和Takahashi将遗传算法应用岩心驱替的历史拟合中，实验结果表明虽然遗传算法能够有效的求解历史拟合问题，但是存在计算效率较低的问题，并且在历史拟合中可能陷入局部收敛。虽然遗传算法在计算过程中能够搜索到较优的解，但是在当油藏模型较大时计算效率较低。2009年Yasin Hajizadeh将ACO算法引入到历史拟合问题的求解中，实验结果表明该算法相对于传统的遗传算法求解效率更高，同年Yasin Hajizad将DE算法引入到历史拟合问题的求解中，该算法仅需要少量的参数就能够实现油藏自动历史拟合，但是上述两种算法在大型油藏模型中难以实现，并且存在如遗传算法易陷入早熟收敛且计算速率缓慢，模拟退火计算量大等问题。

另外，传统方法往往采用高斯分布得到渗透率等模型参数初值，但由于油藏的强非均质性，特别是在多次注水,注化学药剂的多次采油之后,储层中各物性的不确定性强，模型参数特征一般呈现尖峰厚尾的非动力学特征，并不满足高斯分布。

发明内容

本发明所要解决的技术问题是提供基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法及系统，能够采用t分布得到初始的油藏模型参数，并采用基于马尔科夫链的蒙特卡洛法不断优化油藏模型参数拟合生产实际动态，得到尽可能接近真实模型的油藏数值模型，使油田开发动态预测的结果更加接近实际生产动态。

本发明解决上述技术问题的技术方案如下：

一方面，本发明提供了基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法，其特征在于，所述方法包括：

S1、采用t分布随机初始化得到初始的油藏静态参数；

S2、根据贝叶斯公式构造油藏模型的目标函数；

S3、采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数；

S4、对迭代优化得到的所有优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值，并输出所述最优目标函数值及其对应的最优油藏静态参数。

本发明的有益效果：本发明提供的一种基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法，采用t分布得到初始的油藏静态参数，然后采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数，并对所有的优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值。本发明根据油藏的强非均质性，采用t分布得到初始的油藏静态参数，符合油藏模型参数特征呈现尖峰厚尾的非动力学特征，基于概率统计基本思想，采用基于马尔科夫链的蒙特卡洛法不断优化模型渗透率等油藏静态参数拟合生产实际动态，对参数空间的不确定性进行量化，使预测值与真实值尽可能接近，得到尽可能接近真实模型的油藏数值模型，使油田开发动态预测的结果更加接近实际生产。本发明自动调整油藏模型参数，以缩短拟合时间，提高历史拟合的效率和精度，研究对后期油藏开采方案的制定，以及后续生产过程优化具有十分重要的意义。

进一步的，所述S2具体包括：

S21、根据贝叶斯公式得到油藏静态参数的后验分布函数的正比公式，所述油藏静态参数的后验分布函数正比于油藏静态参数的先验t分布的概率函数与油藏生产动态数据的正态分布的似然函数的乘积；

S22、根据t分布公式和正态分布公式，得到所述油藏静态参数的后验分布函数的等式公式，具体由先验项函数和似然项函数相加得到；

S23、将所述油藏静态参数的后验分布函数的公式作为油藏模型的目标函数。

采用上述进一步方案的有益效果：可以将求解油藏静态参数的问题转化为求解使目标函数的最小值，便于求解合适的油藏静态参数。

进一步的，所述S3具体包括：

S31、设置马尔科夫链的链长和优化停止条件，将t分布得到的初始的油藏静态参数作为当前状态对应的优化油藏静态参数，并将所述当前状态放入马尔科夫链中；

S32、根据当前状态的对应的优化油藏静态参数，计算得到当前状态对应的后验分布函数值；

S33、迭代产生下一个状态对应的油藏静态参数，并根据所述下一个状态对应的油藏静态参数，计算得到所述下一个状态对应的后验分布函数值；

S34、根据所述下一个状态对应的后验分布函数值和当前状态对应的后验分布函数值，计算得到接受率R；

S35、从0～1的均匀分布中随机取一个数y，如果y≤R,则接受所述下一个状态，并替代当前状态成为新的当前状态放入马尔科夫链中，所述下一个状态对应的油藏静态参数成为新的当前状态对应的优化油藏静态参数；否则不接受所述下一个状态，依然将当前状态放入马尔科夫链中；

S36、判断是否满足所述优化停止条件，若满足则结束流程，否则返回步骤S33。

采用上述进一步方案的有益效果：采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，缩短了拟合时间，提升了拟合精度，克服了传统随机类方法运行计算量大的问题。

进一步的，根据油藏静态参数计算得到分布函数值具体包括：

根据油藏静态参数计算得到目标函数中先验项函数对应的先验项函数值；

对所述油藏静态参数采用油藏模拟器进行油藏拟合模拟计算，得到油藏生产动态数据；

根据所述油藏静态参数和所述油藏生产动态数据计算得到所述目标函数中似然项函数对应的似然项函数值；

根据所述先验项函数值和所述似然项函数值计算得到所述后验分布函数值。

采用上述进一步方案的有益效果：计算得到后验分布函数值，以便用于计算得到接受率，并且所述后验分布函数值即为目标函数值，以便对目标函数值进行最大后验估计。

另一方面，本发明提供了基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统，所述系统包括：

初始化模块，用于采用t分布随机初始化得到初始的油藏静态参数；

构造模块，用于根据贝叶斯公式构造油藏模型的目标函数；

优化模块，用于采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数；

后验估计模块，用于对迭代优化得到的所有优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值；

输出模块，用于输出所述最优目标函数值及其对应的最优油藏静态参数。

本发明的有益效果：本发明提供的一种基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统，采用t分布得到初始的油藏静态参数，然后采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数，并对所有的优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值。本发明根据油藏的强非均质性，采用t分布得到初始的油藏静态参数，符合油藏模型参数特征呈现尖峰厚尾的非动力学特征，基于概率统计基本思想，采用基于马尔科夫链的蒙特卡洛法不断优化模型渗透率等油藏静态参数拟合生产实际动态，对参数空间的不确定性进行量化，使预测值与真实值尽可能接近，得到尽可能接近真实模型的油藏数值模型，使油田开发动态预测的结果更加接近实际生产。本发明自动调整油藏模型参数，以缩短拟合时间，提高历史拟合的效率和精度，研究对后期油藏开采方案的制定，以及后续生产过程优化具有十分重要的意义。

进一步的，所述构造模块具体包括：

第一构造单元，用于根据贝叶斯公式得到油藏静态参数的后验分布函数的正比公式，所述油藏静态参数的后验分布函数正比于油藏静态参数的先验t分布的概率函数与油藏生产动态数据的正态分布的似然函数的乘积；

第二构造单元，用于根据t分布公式和正态分布公式，得到所述油藏静态参数的后验分布函数的等式公式，具体由先验项函数和似然项函数相加得到；

目标函数构造单元，用于将所述油藏静态参数的后验分布函数的公式作为油藏模型的目标函数。

采用上述进一步方案的有益效果：可以将求解油藏静态参数的问题转化为求解使目标函数的最小值，便于求解合适的油藏静态参数。

进一步的，所述优化模块具体包括：

设置单元，用于设置马尔科夫链的链长和优化停止条件，将t分布得到的初始的油藏静态参数作为当前状态对应的优化油藏静态参数，并将所述当前状态放入马尔科夫链中；

分布函数计算单元，用于根据当前状态的对应的优化油藏静态参数，计算得到当前状态对应的后验分布函数值，以及用于根据下一个状态对应的油藏静态参数，计算得到下一个状态对应的后验分布函数值；

迭代单元，用于迭代产生下一个状态对应的油藏静态参数；

接受率计算单元，用于根据所述下一个状态对应的后验分布函数值和当前状态对应的后验分布函数值，计算得到接受率R；

替换单元，用于从0～1的均匀分布中随机取一个数y，如果y≤R,则接受所述下一个状态，并替代当前状态成为新的当前状态放入马尔科夫链中，所述下一个状态对应的油藏静态参数成为新的当前状态对应的优化油藏静态参数；否则不接受所述下一个状态，依然将当前状态放入马尔科夫链中；

判断单元，用于判断是否满足所述优化停止条件，若满足则结束流程，否则转至所述迭代单元。

采用上述进一步方案的有益效果：采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，缩短了拟合时间，提升了拟合精度，克服了传统随机类方法运行计算量大的问题。

进一步的，所述分布函数计算单元具体包括：

先验项计算单元，用于根据油藏静态参数计算得到目标函数中先验项函数对应的先验项函数值；

油藏生产动态数据计算单元，对所述油藏静态参数采用油藏模拟器进行油藏拟合模拟计算，得到油藏生产动态数据；

似然项计算单元，用于根据所述油藏静态参数和所述油藏生产动态数据计算得到所述目标函数中似然项函数对应的似然项函数值；

函数值计算单元，用于根据所述先验项函数值和所述似然项函数值计算得到所述后验分布函数值。

采用上述进一步方案的有益效果：计算得到后验分布函数值，以便用于计算得到接受率，并且所述后验分布函数值即为目标函数值，以便对目标函数值进行最大后验估计。

附图说明

图1为本发明实施例1的基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法流程图；

图2为本发明实施例1的概率密度函数随自由度v的变化图；

图3为本发明实施例1的基于t分布的马尔科夫链蒙特卡洛自动历史拟合法详细流程图；

图4为本发明实施例1的PUNQS3模型的水平渗透率分布图；

图5为本发明实施例1的井1的井底压力、气油比、含水率拟合图；

图6为本发明实施例1的井4的井底压力、气油比、含水率拟合图；

图7为本发明实施例1的井5的井底压力、气油比、含水率拟合图；

图8为本发明实施例1的井11的井底压力、气油比、含水率拟合图；

图9为本发明实施例1的井12的井底压力、气油比、含水率拟合图；

图10为本发明实施例1的井15的井底压力、气油比、含水率拟合图；

图11为本发明实施例1的井1的井底压力、气油比、含水率拟合图；

图12为本发明实施例1的井4的井底压力、气油比、含水率拟合图；

图13为本发明实施例1的井5的井底压力、气油比、含水率拟合图；

图14为本发明实施例1的井11的井底压力、气油比、含水率拟合图；

图15为本发明实施例1的井12的井底压力、气油比、含水率拟合图；

图16为本发明实施例1的井15的井底压力、气油比、含水率拟合图；

图17为本发明实施例2的基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统示意图；

图18为本发明实施例2的基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统中构造模块的结构示意图；

图19为本发明实施例2的基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统中的优化模块的结构示意图；

图20为本发明实施例2的基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统中的优化模块中的分布函数计算单元的结构示意图。

具体实施方式

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

实施例1、基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法。下面结合图1至图16对本实施例提供的方法进行详细说明。

参见图1至图3，S1、采用t分布随机初始化得到初始的油藏静态参数。

具体的，油藏数值模型中需要优化的参数是各划分网格的油藏静态参数如渗透率、孔隙度等，通过某种概率分布模型随机赋初始值，传统方法往往采用高斯分布得到渗透率等油藏静态参数初始值，但由于油藏的强非均质性，特别是在多次注水、注化学药剂的多次采油之后,储层中各物性的不确定性强，模型参数特征一般呈σ现尖峰厚尾的非动力学特征，不满足高斯分布。当实际单个变量的边际分布比正常的边际分布的尾部大的时候，可使用t分布来代替正态分布。t分布曲线的形状与自由度v的大小相关，如图2所示，自由度v越小，t分布曲线越平坦，曲线中间的值越低，曲线双侧尾部越高；自由度v越大，t分布曲线越接近正态分布曲线，当自由度v→∞时，t分布曲线逐步趋近于标准正态分布曲线。因而采用t分布随机初始化得到初始的油藏静态参数。所述油藏静态参数可为个生产井的渗透率、孔隙度等参数，也可以为区块各时间片的渗透率、孔隙度等参数。

S2、根据贝叶斯公式构造油藏模型的目标函数。

具体的，所述S2具体包括以下步骤：

S21、根据贝叶斯公式得到油藏静态参数的后验分布函数的正比公式，所述油藏静态参数的后验分布函数正比于油藏静态参数的先验t分布的概率函数与油藏生产动态数据的正态分布的似然函数的乘积。

具体的，传统的贝叶斯方法应用在储层数值模拟的时候，通过评估“最可能模型”来建立模型，其中先验分布主要用于描述油藏静态参数如孔隙度和渗透率等是否符合某种概率分布，后验估计在抽样以后可以得到，求解油藏静态参数m的问题可以转化为使目标函数O(m)取得最小值的问题。

由贝叶斯公式得到油藏静态参数m的后验分布函数的正比公式，如公式(1)所示：

p(m|d_obs)∝p(d_obs|m)·p(m) (1)

其中，d_obs为油藏生产动态数据，即为含水率、井底压力和气油比等参数；m为不确定的参数，即为渗透率等待优化的油藏静态参数；p(m)为m的先验t分布的概率函数；p(d_obs|m)为油藏生产动态数据的正态分布似然函数；p(m|d_obs)为m的后验分布函数。

S22、根据t分布公式和正态分布公式，得到所述油藏静态参数的后验分布函数的等式公式，具体由先验项函数和似然项函数相加得到。

具体的，根据t分布的概率密度公式和正态分布的概率密度公式，得到所述油藏静态参数m的后验分布函数的具体的等式公式，其中，正态分布的似然函数具体为公式(2)所示：

$f\left(x\right)=\frac{1}{\sqrt{\left|\mathrm{\Σ}\right|{\left(2\mathrm{\π}\right)}^d}}\exp \{-\frac{1}{2}\{{\[g\left(m\right)-d_{obs}\]}^T{C}_D^{-1}\[g\left(m\right)-d_{obs}\]\left\}\right\}---\left(2\right)$

其中，d为向量m的维度；g(m)为符合不确定性参数的先验概率分布；Σ为协方差矩阵,d_obs为油藏生产动态数据。

t分布的概率密度函数如公式(3)所示：

$f\left(x\right)=\frac{1}{|\mathrm{\Σ}{|}^{1/2}}\frac{1}{\sqrt{{\left(v\mathrm{\π}\right)}^d}}\frac{\mathrm{\Γ}\left(\right(v+d)/2)}{\mathrm{\Γ}(v/2)}{(1+\frac{m^{\′}{\mathrm{\Σ}}^{-1}m}{v})}^{-(v+d)/2}---\left(3\right)$

其中，x为向量，v为自由度,Σ为协方差矩阵,d为向量m的维度。

因而，所述油藏静态参数m的后验分布函数的具体的等式公式如公式(4)所示：

$\begin{array}{c}P\left(m\right|d_{obs})\\ =\frac{1}{|\mathrm{\Σ}{|}^{1/2}}\frac{1}{\sqrt{{\left(v\mathrm{\π}\right)}^d}}\frac{\mathrm{\Γ}\left(\right(v+d)/2)}{\mathrm{\Γ}(v/2)}{(1+\frac{m^{\′}{\mathrm{\Σ}}^{-1}m}{v})}^{-(v+d)/2}+\frac{1}{\sqrt{\left|\mathrm{\Σ}\right|{\left(2\mathrm{\π}\right)}^d}}\exp \{-\frac{1}{2}\{{\[g\left(m\right)-d_{obs}\]}^T{C}_D^{-1}\[g\left(m\right)-d_{obs}\]\left\}\right\}\end{array}---\left(4\right)$

其中，v为自由度；d为向量m的维度；g(m)为符合不确定性参数的先验概率分布；Σ为协方差矩阵。

S23、将所述油藏静态参数的后验分布函数的公式作为油藏模型的目标函数。

具体的，将所述油藏静态参数的后验分布函数的公式作为油藏模型的目标函数O(m)，具体如公式(5)所示：

$\begin{array}{c}0\left(m\right)=\\ \frac{1}{|\mathrm{\Σ}{|}^{1/2}}\frac{1}{\sqrt{{\left(v\mathrm{\π}\right)}^d}}\frac{\mathrm{\Γ}\left(\right(v+d)/2)}{\mathrm{\Γ}(v/2)}{(1+\frac{m^{\′}{\mathrm{\Σ}}^{-1}m}{v})}^{-(v+d)/2}+\frac{1}{\sqrt{\left|\mathrm{\Σ}\right|{\left(2\mathrm{\π}\right)}^d}}\exp \{-\frac{1}{2}\{{\[g\left(m\right)-d_{obs}\]}^T{C}_D^{-1}\[g\left(m\right)-d_{obs}\]\left\}\right\}\end{array}$

其中，μ为先验值，为先验项函数；为似然项函数。

S3、采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数。

具体的，在连续的历史拟合过程中，采用基于马尔科夫链的蒙特卡洛法更新模型参数。其原理是采用先验t分布随机采样得到油藏静态参数的初始状态，并循环进行状态转移，当随机取的数小于等于接受率时，接收下一状态，否则舍弃，将当前状态放入链中；重复上述操作以得到优化的油藏静态参数。

对于每一条马尔科夫链需要考虑四个参数：

(1)、马尔科夫链的初始状态initial，表征马尔科夫链随机取样的起点；

(2)、先验项logprior，表征计算先验项；

(3)、后验项loglikelihood，表征计算后验项；

(4)、马尔科夫链链长mccount，表征马尔科夫链转移状态的长度。

基于先验t分布随机生成初始状态，所述初始状态对应t分布随机初始化得到的油藏静态参数的初始值，根据马尔科夫链生成下一个状态，计算得到接受率，从均匀分布中随机取一个数，当随机取的数小于等于接受率时，接收下一状态，否则舍弃下一状态，将当前状态放入链中。从而进行不断地循环，修改油藏参数，求解得到与生产历史拟合匹配的最优解。

主要输入的数据包括各类静动态数据，如生产动态数据包括各油井的产液量、日产油、含水率等，网格数据，相对渗透率、毛管压力以及油藏流体的PVT属性数据，油、水、气的地面密度，岩石压缩系数等物性参数。

具体的，所述S3具体包括以下步骤：

S31、设置马尔科夫链的链长和优化停止条件，将t分布得到的初始的油藏静态参数作为当前状态对应的优化油藏静态参数，并将所述当前状态放入马尔科夫链中。所述优化停止条件具体为到达马尔科夫链的链长，也可根据情况设置其他停止条件。

S32、根据当前状态的对应的优化油藏静态参数，计算得到当前状态对应的后验分布函数值。

S33、迭代产生下一个状态对应的油藏静态参数，并根据所述下一个状态对应的油藏静态参数，计算得到所述下一个状态对应的后验分布函数值。

S34、根据所述下一个状态对应的后验分布函数值和当前状态对应的后验分布函数值，计算得到接受率R。

具体的，根据所述下一个状态对应的后验分布函数值和当前状态对应的后验分布函数值，计算得到接受率R，所述接收率R的计算公式如公式(6)所示：

$R=\frac{P\left(m^{t_{i+1}}\right|d_{obs}^{t_{i+1}})}{P\left(m^{t_i}\right|d_{obs}^{t_i})}---\left(6\right)$

其中，所述为下一个状态对应的后验分布函数值，为当前状态对应的后验分布函数值。

S35、从0～1的均匀分布中随机取一个数y，如果y≤R,则接受所述下一个状态，并替代当前状态成为新的当前状态放入马尔科夫链中，所述下一个状态对应的油藏静态参数成为新的当前状态对应的优化油藏静态参数；否则不接受所述下一个状态，依然将当前状态放入马尔科夫链中。

具体的，只有在一个状态被接受后，该状态对应的油藏静态参数才会成为迭代得到的优化油藏静态参数；如果一个状态没有被接受，则该状态对应的油藏静态参数不会成为迭代得到的优化油藏静态参数，该油藏静态参数会被舍弃。

S36、判断是否满足所述优化停止条件，若满足则结束流程，否则返回步骤S33。

具体的，根据油藏静态参数计算得到分布函数值具体包括以下步骤：

a、根据油藏静态参数计算得到目标函数中先验项函数对应的先验项函数值。

b、对所述油藏静态参数进行油藏拟合模拟计算，得到油藏生产动态数据。具体为调用油藏数值模拟器进行油藏拟合模拟计算得到油藏生产动态数据。

c、根据所述油藏静态参数和所述油藏生产动态数据计算得到所述目标函数中似然项函数对应的似然项函数值。

d、根据所述先验项函数值和所述似然项函数值计算得到所述后验分布函数值。所述后验分布函数值与目标函数值相同。

无论是当前状态还是下一个状态对应的后验分布函数值均用上述方法进行计算。

算法马尔科夫链的蒙特卡洛法具体包括以下步骤：

首先设置输入输出参数：输入参数为：链长t为N(N正整数)、先验值μ、先验项的计算函数、似然项的计算函数以及输出结果设置；输出参数为：马尔科夫链。

步骤1、初始化马尔科夫链的初始状态即链长t＝1时的初始状态，即第1个状态。

步骤2、对链长t＝2，3，4,…,N；i＝1，2，3，…，N,循环以下过程进行采样：

2.1、从第i个状态到第i+1个状态，计算下一个状态的值；

2.2、计算接收率R；

2.3、从0～1的均匀分布中随机取一个数y，如果y≤R,则接受下一状态，并放入马尔科夫链中；否则又将当前状态放入马尔科夫链中。

基于t分布的马尔科夫链蒙特卡洛油藏参数自动历史拟合法具体包括以下步骤：

首先设置输入输出参数：输入参数为：油藏模型数据文件、链长t为N(N正整数)、先验值μ、输出结果位置和格式；输出参数为：渗透率优化参数和目标函数值。

步骤1、基于t分布初始化马尔科夫链的初始状态并计算得到初始状态对应的后验分布函数值；

步骤2、对t＝2，3，4，…，N；i＝1，2，3，…，N,直到达到马尔科夫链的链长，循环以下过程进行渗透率的参数历史拟合；

2.1、产生可能的下一个状态，并计算下一个状态对应的先验项函数值；

2.2、运行油藏数值模拟程序ECLI PSE，计算得到油藏生产动态数据；

2.3、计算所述下一个状态对应的似然函数值；

2.4、计算下一个状态对应的后验分布函数值；

2.5、使用马尔科夫链的蒙特卡洛法判断是否跳转到下一个状态。

步骤3、对目标函数进行最大后验估计，得到渗透率分布。

S4、对迭代优化得到的所有优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值，并输出所述最优目标函数值及其对应的最优油藏静态参数。

具体的，在采用所述马尔科夫链蒙特卡洛法自动历史拟合迭代更新结束后，会将在迭代更新过程中所有的优化油藏静态参数值以及其对应的所有的目标函数值进行输出，所述目标函数值即所述后验分布函数值，然后对所有的目标函数值进行最大后验估计，得到最优目标函数值，并输出所述最优目标函数值及其对应的最优油藏静态参数，所述最优油藏静态参数即渗透率分布。所有的优化油藏静态参数值具体指在迭代优化过程中被接受的状态的对应的油藏静态参数，即被放入马尔科夫链中的状态对应的油藏静态参数。

另外，也可以在每次迭代更新得到一个优化油藏静态参数，就对所述优化油藏静态参数对应的目标函数值进行最大后验估计，然后在迭代结束之后直接输出得到的最优目标函数值及其对应的最优油藏静态参数。

综上所述，采用基于t分布的马尔科夫链蒙特卡洛油藏自动历史拟合方法，通过调用油藏数值模拟软件计算，使预测值与真实值尽可能接近，得到和真实油藏模型较为一致的数值模型。

具体实例：

1、主要通过对基于t分布的马尔科夫链蒙特卡洛油藏参数自动历史拟合法进行实验以检验其效果。实验采用的是PUNQ-S3油藏数据模型，PUNQ-S3油藏数据模型是一个三维三相的油藏工程模型，由19*28*25个网格块构成，共分为五层，每层为2660个网格块，每个网格块的大小一致，其中包含1761个有效网格块。如图2所示，空白部分表示的是无效网格，线段特点的网格表示的是不同数值的水平渗透率，对于模型的每一层可以将水平渗透率分为不同的区块，综上所述，可将PUNQS3油藏模型的1761个网格，分为5*9共45个区块，达到历史拟合分区分块的目的。PUNQS3模型每层的水平渗透率分布图如图4所示。

2、油藏单井历史拟合情况比较

根据基于t分布的马尔科夫链蒙特卡洛油藏自动历史拟合法的实验结果，对历史拟合进行了相关的实验和分析，其中链长设为500。目标函数值越小说明拟合值与实际测量值之间的差异程度越小，即拟合效果越好，效果越优。

为进一步比较说明基于正态分布和t分布的马尔科夫链蒙特卡洛自动历史拟合方法的效果，将计算出的单口井的含水率(WWCT)、井底压力(WBHP)和气油比(WGOR)等参数同模型真实值分别进行对比，如图5到图10所示。其中，图5-10分别是生产井1、4、5、11、12和15六口井的参数拟合图，其中节点为圆圈的曲线为基于正态分布的马尔科夫链的蒙特卡洛自动历史拟合方法的拟合曲线，节点为三角的曲线为基于t分布的马尔科夫链的蒙特卡洛自动历史拟合方法的拟合曲线，节点为正方形的曲线为模型真实值。

为进一步说明拟合效果，对拟合结果采用均方根误差(RE)及整体误差(EE)进行了计算和统计，计算公式如公式(7)和(8)所示：

$RE=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{\[D_i\]}^2}---\left(7\right)$

$EE=\underset{i=1}{\overset{N}{\Σ}}\[|N_i-{N_i}^{\′}|\]---\left(8\right)$

其中，N为维度，Di＝N_i-N′_i，N_i为真实值，N′_i为拟合值。

基于正态分布马尔科夫链的蒙特卡洛自动历史拟合方法的拟合误差统计结果如表1所示，基于t分布马尔科夫链的蒙特卡洛自动历史拟合方法的拟合误差统计结果如表2所示：

表1基于正态分布马尔科夫链的蒙特卡洛自动历史拟合方法的误差表

表2基于t分布马尔科夫链的蒙特卡洛自动历史拟合方法的误差表

3、预测

为进一步说明马尔科夫链的蒙特卡洛自动历史拟合方法生产预测的效果，取每口井前1460天的生产数据作为训练集，训练得到模型，然后采用训练模型预测后1540天的生产数据并与真实数据进行对比，如图11-16所示，其中，图11-16分别是生产井1、4、5、11、12和15六口井的含水率(WWCT)、井底压力(WBHP)和气油比(WGOR)的预测结果。其中节点为圆圈的曲线为真实值，节点为三角的曲线为基于正态分布的马尔科夫链的蒙特卡洛法的拟合曲线，节点为正方形的曲线为基于t分布的马尔科夫链的蒙特卡洛法的拟合曲线。

采用均方根误差(RE)及整体误差(EE)进一步计算和统计，基于正态分布马尔科夫链的蒙特卡洛自动历史拟合方法的拟合误差统计结果如表3所示，基于t分布马尔科夫链的蒙特卡洛自动历史拟合方法的拟合误差统计结果如表4所示：

表3基于正态分布马尔科夫链的蒙特卡洛自动历史拟合方法预测误差表

表4基于t分布马尔科夫链的蒙特卡洛自动历史拟合方法预测误差表

综上所述，基于t分布的马尔科夫链蒙特卡洛油藏自动历史拟合法与基于正态分布的马尔科夫链蒙特卡洛自动历史拟合法相比，大部分井参数预测值与真实值间的误差更小，减小了油藏模型认识的不确定性，提高模型的预测能力，效果更优。

实施例2、基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统。下面结合图17至图20对本实施例提供的系统进行详细说明。

参见图17至图20，基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统，其特征在于，所述系统包括初始化模块、构造模块、优化模块、后验估计模块以及输出模块。

初始化模块，用于采用t分布随机初始化得到初始的油藏静态参数。

构造模块，用于根据贝叶斯公式构造油藏模型的目标函数。

具体的，所述构造模块具体包括第一构造单元、第二构造单元以及目标函数构造单元。

第一构造单元，用于根据贝叶斯公式得到油藏静态参数的后验分布函数的正比公式，所述油藏静态参数的后验分布函数正比于油藏静态参数的先验t分布的概率函数与油藏生产动态数据的正态分布的似然函数的乘积。

第二构造单元，用于根据t分布公式和正态分布公式，得到所述油藏静态参数的后验分布函数的等式公式，具体由先验项函数和似然项函数相加得到。

目标函数构造单元，用于将所述油藏静态参数的后验分布函数的公式作为油藏模型的目标函数。

优化模块，用于采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数。

所述优化模块具体包括设置单元、分布函数计算单元、迭代单元、接受率计算单元、替换单元以及判断单元。

设置单元，用于设置马尔科夫链的链长和优化停止条件，将t分布得到的初始的油藏静态参数作为当前状态对应的优化油藏静态参数，并将所述当前状态放入马尔科夫链中。

分布函数计算单元，用于根据当前状态的对应的优化油藏静态参数，计算得到当前状态对应的后验分布函数值，以及用于根据下一个状态对应的油藏静态参数，计算得到下一个状态对应的后验分布函数值。

所述分布函数计算单元具体包括先验项计算单元、油藏生产动态数据计算单元、似然项计算单元以及函数值计算单元。

先验项计算单元，用于根据油藏静态参数计算得到目标函数中先验项函数对应的先验项函数值。

油藏生产动态数据计算单元，对所述油藏静态参数采用油藏模拟器进行油藏拟合模拟计算，得到油藏生产动态数据。

似然项计算单元，用于根据所述油藏静态参数和所述油藏生产动态数据计算得到所述目标函数中似然项函数对应的似然项函数值。

函数值计算单元，用于根据所述先验项函数值和所述似然项函数值计算得到所述后验分布函数值。

所述迭代单元，用于迭代产生下一个状态对应的油藏静态参数。

所述接受率计算单元，用于根据所述下一个状态对应的后验分布函数值和当前状态对应的后验分布函数值，计算得到接受率R。

所述替换单元，用于从0～1的均匀分布中随机取一个数y，如果y≤R,则接受所述下一个状态，并替代当前状态成为新的当前状态放入马尔科夫链中，所述下一个状态对应的油藏静态参数成为新的当前状态对应的优化油藏静态参数；否则不接受所述下一个状态，依然将当前状态放入马尔科夫链中。具体的，若所述下一个状态不被接受，该状态对应的油藏静态参数无法成为优化油藏静态参数，而是被直接舍弃掉。

所述判断单元，用于判断是否满足所述优化停止条件，若满足则结束流程，否则转至所述迭代单元。

后验估计模块，用于对迭代优化得到的所有优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值。

具体的，所述目标函数值即为计算得到的后验分布函数值，二者的的值相同，另外，只有优化油藏静态参数对应的目标函数值才会进行最大后验估计。

输出模块，用于输出所述最优目标函数值及其对应的最优油藏静态参数。

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

Claims (8)

1.基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法，其特征在于，所述方法包括：

S1、采用t分布随机初始化得到初始的油藏静态参数；

S2、根据贝叶斯公式构造油藏模型的目标函数；

S3、采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数；

S4、对迭代优化得到的所有优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值，并输出所述最优目标函数值及其对应的最优油藏静态参数。

2.如权利要求1所述的基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法，其特征在于，所述S2具体包括：

S21、根据贝叶斯公式得到油藏静态参数的后验分布函数的正比公式，所述油藏静态参数的后验分布函数正比于油藏静态参数的先验t分布的概率函数与油藏生产动态数据的正态分布的似然函数的乘积；

S22、根据t分布公式和正态分布公式，得到所述油藏静态参数的后验分布函数的等式公式，具体由先验项函数和似然项函数相加得到；

S23、将所述油藏静态参数的后验分布函数的公式作为油藏模型的目标函数。

3.如权利要求2所述的基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法，其特征在于，所述S3具体包括；

S31、设置马尔科夫链的链长和优化停止条件，将t分布得到的初始的油藏静态参数作为当前状态对应的优化油藏静态参数，并将所述当前状态放入马尔科夫链中；

S32、根据当前状态的对应的优化油藏静态参数，计算得到当前状态对应的后验分布函数值；

S33、迭代产生下一个状态对应的油藏静态参数，并根据所述下一个状态对应的油藏静态参数，计算得到所述下一个状态对应的后验分布函数值；

S34、根据所述下一个状态对应的后验分布函数值和当前状态对应的后验分布函数值，计算得到接受率R；

S35、从0～1的均匀分布中随机取一个数y，如果y≤R,则接受所述下一个状态，并替代当前状态成为新的当前状态放入马尔科夫链中，所述下一个状态对应的油藏静态参数成为新的当前状态对应的优化油藏静态参数；否则不接受所述下一个状态，依然将当前状态放入马尔科夫链中；

S36、判断是否满足所述优化停止条件，若满足则结束流程，否则返回步骤S33。

4.如权利要求3所述的基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法，其特征在于，根据油藏静态参数计算得到分布函数值具体包括：

根据油藏静态参数计算得到目标函数中先验项函数对应的先验项函数值；

对所述油藏静态参数采用油藏模拟器进行油藏拟合模拟计算，得到油藏生产动态数据；

根据所述油藏静态参数和所述油藏生产动态数据计算得到所述目标函数中似然项函数对应的似然项函数值；

根据所述先验项函数值和所述似然项函数值计算得到所述后验分布函数值。

5.基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统，其特征在于，所述系统包括：

初始化模块，用于采用t分布随机初始化得到初始的油藏静态参数；

构造模块，用于根据贝叶斯公式构造油藏模型的后验分布函数；

优化模块，用于采用马尔科夫链蒙特卡洛历史拟合法对所述油藏静态参数进行迭代优化，得到优化油藏静态参数；

后验估计模块，用于对迭代优化得到的所有优化油藏静态参数对应的所有目标函数值进行最大后验估计，得到最优目标函数值；

输出模块，用于输出所述最优目标函数值及其对应的最优油藏静态参数。

6.如权利要求5所述的基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统，其特征在于，所述构造模块具体包括：

第一构造单元，用于根据贝叶斯公式得到油藏静态参数的后验分布函数的正比公式，所述油藏静态参数的后验分布函数正比于油藏静态参数的先验t分布的概率函数与油藏生产动态数据的正态分布的似然函数的乘积；

第二构造单元，用于根据t分布公式和正态分布公式，得到所述油藏静态参数的后验分布函数的等式公式，具体由先验项函数和似然项函数相加得到；

目标函数构造单元，用于将所述油藏静态参数的后验分布函数的公式作为油藏模型的目标函数。

7.如权利要求6所述的基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统，其特征在于，所述优化模块具体包括：

设置单元，用于设置马尔科夫链的链长和优化停止条件，将t分布得到的初始的油藏静态参数作为当前状态对应的优化油藏静态参数，并将所述当前状态放入马尔科夫链中；

分布函数计算单元，用于根据当前状态的对应的优化油藏静态参数，计算得到当前状态对应的后验分布函数值，以及用于根据下一个状态对应的油藏静态参数，计算得到下一个状态对应的后验分布函数值；

迭代单元，用于迭代产生下一个状态对应的油藏静态参数；

接受率计算单元，用于根据所述下一个状态对应的后验分布函数值和当前状态对应的后验分布函数值，计算得到接受率R；

替换单元，用于从0～1的均匀分布中随机取一个数y，如果y≤R,则接受所述下一个状态，并替代当前状态成为新的当前状态放入马尔科夫链中，所述下一个状态对应的油藏静态参数成为新的当前状态对应的优化油藏静态参数；否则不接受所述下一个状态，依然将当前状态放入马尔科夫链中；

判断单元，用于判断是否满足所述优化停止条件，若满足则结束流程，否则转至所述迭代单元。

8.如权利要求7所述的基于t分布的马尔科夫链蒙特卡洛自动历史拟合系统，其特征在于，所述分布函数计算单元具体包括：

先验项计算单元，用于根据油藏静态参数计算得到目标函数中先验项函数对应的先验项函数值；

油藏生产动态数据计算单元，对所述油藏静态参数采用油藏模拟器进行油藏拟合模拟计算，得到油藏生产动态数据；

似然项计算单元，用于根据所述油藏静态参数和所述油藏生产动态数据计算得到所述目标函数中似然项函数对应的似然项函数值；

函数值计算单元，用于根据所述先验项函数值和所述似然项函数值计算得到所述后验分布函数值。