基于贝叶斯理论的四维地震与生产动态数据的油藏历史拟合新方法

王泽龙 刘先贵 李志勇 张 辉 张传进

（1. 中国石油国际勘探开发有限公司；2. 中国石油勘探开发研究院）

随着油气行业的智能化转型，在油藏的精细化管理和生产过程中，人们越来越重视数字化信息的实时采集，并依靠计算机对各类生产数据的分析，不断修改油藏模型，及时调整油田开发技术政策，以便实现油藏的经济、高效开发[1]。目前，油藏历史拟合技术一般只能利用井口生产动态数据[2]，通过人工手动或计算机辅助的方法修改各类参数，将油藏数值模拟得到的预测井口生产数据与井口测量的数据进行拟合，以便得到更为可靠的油藏模型。但由于井口生产数据的空间覆盖范围十分有限，油藏历史拟合缺少足够的约束条件，以此为基础修正的油藏模型往往不符合实际地质情况。

油藏历史拟合问题的数学本质是求解最优化。目前，梯度类算法和集合类算法是解决油藏历史拟合问题的两类常见的优化算法。其中，Levenberg-Marquardt（LM）算法是常用的梯度类算法，该算法用于小型反演问题效率很高，但对于大型油藏模型由于需要计算和存储超高维度的敏感度矩阵，造成计算量过于巨大而难以广泛应用[3]。集合类算法是基于贝叶斯理论的后验概率来进行最优化求解，1960年，Kalman[4]提出卡尔曼滤波，通过对平均值和均方误差的最小化问题求解，来实现对某一过程中的状态参数高效地进行循环贝叶斯概率推断；1970年，Buck[5]和Sunahara[6]将卡尔曼滤波扩展到非线性问题，提出扩展卡尔曼滤波（EKF）；1994年，Evensen[7]提出集合卡尔曼滤波（EnKF），将卡尔曼滤波与蒙特卡洛方法[8]相结合，形成了集合类历史拟合算法的雏形；2002年，Nævdal[9]等第一次将EnKF引入石油工程领域，用于更新油井附近的油藏模型，随后 Lorentzen[10]、Chen[11]、Li[12]、Jafarpour[13]等学者加入研究，使得EnKF算法成为当前油藏历史拟合最热门的研究领域；2013年，结合集合平滑算法与EnKF，Emerick[14]提出基于多次数据吸收的集合平滑算法（ES-MDA），为油藏历史拟合问题提供了一种新思路。

四维地震监测是在同一个地震采集区，间隔不同的时间部署三维地震，通过分析不同时间点的三维地震属性差异来监测油藏动态情况，可为油气田精细开发管理提供覆盖面广、信息量大的测量数据资料，据此计算出油藏压力和温度、储层生产能力和流体变化，对明确剩余油气分布、提高采收率有着重要意义。由于四维地震数据的复杂程度及与井口生产数据存在差异，当前的油藏历史拟合尚难以实现对四维地震数据的有效利用。本文研究提出了一种基于四维地震数据的历史拟合新方法，通过对四维地震数据的粗化与稀疏化处理，实现地震数据与井口数据的匹配，将由四维地震监测获得的纵波波阻抗数据与井口生产数据相结合，并通过建立起油藏状态与地震波响应的数学模型，形成油藏流体—地震属性耦合模拟程序，改进ES-MDA算法，实现四维地震数据与井口生产数据的联合油藏历史拟合。

1 常规油藏历史拟合算法

1.1 油藏历史拟合的目标函数

油藏自动历史拟合问题可以转为在约束条件下的最优化问题。这个最优化问题的目标函数可以用最小二乘法误差表示[15]，定义如下：

式中m——油藏地质模型参数矩阵；

mf——修改后的油藏地质模型参数矩阵；

f（x）——油藏流体数值模拟程序；

duc——实际测量的数据向量；

CM——油藏地质模型参数m的协方差矩阵；

CD—— 计算结果f（m）和实测数据duc的协方差矩阵。

目标函数分为两部分，每个部分的数值均大于或等于零。为使目标函数最小化，既需要模拟计算结果与实测数据相吻合，即f（m）-duc趋近于零，也需要修改后的地质模型不要偏离初始模型太远，即f（mmf）趋近于零，用来表征更新过程受到足够的地质资料约束。

1.2 EnKF算法

对于油藏历史拟合，定义扩展状态向量xk为k时刻估算的状态参数向量，其中包括了静态地质模型的参数、主要油藏动态参数与实测数据[16]。

式中mk—— 静态地质模型参数（如孔隙度、渗透率）在时间k时的估值，其维度是Nm×1；

Pk—— 动态参数（如油藏压力、油气水饱和度），其维度是Np×1；

dk—— 实测数据（如油气产量、产水量、井底压力或四维地震数据），其维度是Nn×1 ；

Nm，Np，Nn—— 分别是单个地质模型中的静态参数数量、单一时刻的动态参数数量、单一时刻的实测数据量；

下角标k—— 时间序列中的第k项，即第k时间步（k为整数）。

对于油藏流体数值模拟的过程定义如下：

上角标f——油藏模型预测的数据；

上角标u——更新后的数据。

至此，完成了在k时刻的后验模型集合的估算。通过不断重复此过程，对比不同时刻的实际测量数据，并不断更新油藏地质模型集合。EnKF的算法流程如图1所示。

图1 EnKF算法流程图

EnKF算法具有如下明显的缺点：（1）应用效果十分依赖于初始模型的质量；（2）在EnKF更新方程中，模型参数和状态参数同时更新，难免造成更新后的状态参数与实际的状态参数不一致。

1.3 ES-MDA算法

在控制优化理论中，对时间序列的预测算法除滤波类算法外，还有平滑（smoother）算法，即利用过去、当前和未来的信息来估算某个时间点的状态参量。如集合平滑算法（ES）对全部时间序列数据整体吸收的方式，避免了EnKF算法更新后的状态参数与实际的状态参数不一致的问题，且对线性问题的计算效率远高于EnKF。2013年，Emerick[14]等受EnKF多次迭代的启发，将EnKF与ES相结合，通过对相同数据的重复吸收来减小修正步长，形成基于多次数据吸收的集合平滑算法（ES-MDA），并实现了对非线性问题的应用。

按照上述3个步骤，使用ES-MDA更新方程对整个集合中的所有子集模型都完成修正后，即完成了一次数据吸收。重复以上步骤Na次，即是ES-MDA算法的全部过程。ES-MDA算法的流程框图如图2所示。

图2 ES-MDA算法流程图

ES-MDA算法因提出时间较晚，仅有少量学者在地下水和油藏模拟中进行了少量理论研究，均取得了较为理想的理论性数值试验结果[17]，该算法在油藏历史拟合的适用性还需要进一步研究。对于四维地震数据结合井口生产数据的油藏历史拟合，理论上ESMDA算法契合了多维度时间序列数据的整体同化要求，但由于原始算法中的预测步骤中仅进行了一次模拟预测，无法满足地震属性—油藏流体耦合模拟的要求，并且四维地震数据与井口生产数据在时间和空间两个维度上都存在很大差异，故仍需要对ES-MDA算法进行优化改进。

2 基于四维地震数据的历史拟合

2.1 四维地震监测对于油藏历史拟合的作用

油藏历史拟合现有的方法仅通过拟合井口生产动态数据来修正油藏的模型参数。诚然，油藏状态的动态变化过程，可以反映在井口的生产动态上，但由于井口的覆盖范围有限，无法准确实现对于整个油藏状态的表征。从物理意义上讲，由于井口动态稀疏有限的井口动态生产数据必然会造成反向模型求解出现严重的欠定和多解性，所求结果往往与实际情况不符。

随着油田生产，不断有油气从油藏中采出，伴随着油藏内流体交换及状态的变化，还会造成油藏压力、温度的变化。其中，油藏压力和流体饱和度等因素的变化，在四维地震监测过程中，也会造成地震波响应的改变。其中，纵波波阻抗是常见表征油藏状态和流体变化的地震属性指标之一。

地震波阻抗反演的地震解释技术早在20世纪70年代就开始形成，经过几十年发展，已经能比较准确地获得波阻抗数据，在历史拟合中若成功应用，可以增强地质模型中流体变化特征的可信度。

伴随着油气开采的过程，一方面，溶解气在储层内脱气、气相推动油水液相，以及流体的变化造成油藏储层内孔隙中流体密度减小；另一方面，油藏内地层压力增加，则会造成孔隙膨胀，油藏岩石密度减小[19]。这两方面情形都会降低地震波反射传播速度，从而减小地震纵波波阻抗，且可以通过有色反演技术表现出来。相反，若在油气开采过程中，压力无法得到及时补充，造成油藏压力下降，导致孔隙被压缩，提高了油藏的岩石和流体整体密度，地震纵波波阻抗数值就会上升[20]。

由此可见，在油气田生产过程中发生的油藏状态变化，包括油藏压力、流体状态、饱和度等指标，都可以通过四维地震监测技术来反映。四维地震监测技术主要包括四维地震数据的采集、处理和解释，并通过岩石物理学理论将反射地震波属性与油藏状态属性联系起来，从而实现对整个油藏范围全覆盖的准动态监测。

2.2 四维地震数据的稀疏化处理

首先，按照地质—油藏模型耦合的通常处理方法，将处理后的地震纵波波阻抗数据体进行网格粗化。经过网格粗化以后，地震属性数据体与油藏流体模拟的模型尺度相一致，在获得油藏流体状态（经油藏流体数值模拟后）的基础上，通过岩石物理学模型正演可直接得到相应地震波阻抗数据。但由于油藏模型网格通常高达上百万网格量级，鉴于当前电脑硬件计算能力限制，仍需对地震波阻抗数据进一步进行稀疏化处理。

在信号处理领域，稀疏表示（sparse representation）可以对连续性的或密集的信号进行稀疏化处理，从而以较小的数据体量来表达大体量矩阵的关键性信息。对于信号s如果只有m（m＜N）个元素为非零值，其余元素均为零，那么该信号s被称为m-稀疏信号[21]。求解稀疏表示的过程被称为稀疏分解。在通常情况下，一个信号可以用K个信号原子φk的线性叠加来表示：

其中，α一般被称为信号在数据字典φ中的表示稀疏系数（representation coefficients）。字典φ是由信号原子φk为纵向组成的矩阵，通常被归一化为方差，即[22]。

对于四维地震中的地震波阻抗数据体，表示的是空间上不同油藏网格中的不同油藏状态下的纵波波阻抗数据值，虽然不具有信号那样时间上的连续性，但考虑到一定范围内地质和油藏属性的连续性，可以仿照信号数组，受信号处理领域中的稀疏表示启发，通过设计一个空间上的单位稀疏矩阵作为稀疏系数，来实现对超大型四维地震纵波波阻抗数据的稀疏化处理。四维地震监测中，分别在t1，t2，t3，...，tN时刻进行三维地震采集，经地震数据处理与解释，以及网格粗化后，每个时间点的纵波地震波阻抗数据为Ip_large（t），其维度与油藏网格的维度一致。对每个时间点的可采用相同的稀疏系数矩阵。Ip_4D为稀疏处理后的四维地震波阻抗矩阵，其过程可通过式（21）表示：

2.3 历史拟合ES-MDA算法优化

原始ES-MDA算法在油藏历史拟合中由于仅吸收了井口动态生产数据，其预测步骤为运行全时段的油藏流体模拟程序，计算出油田生产历史的完整井口生产动态预测数据。由于需要耦合油藏流体数值模拟与地震属性正向模拟，原始ES-MDA算法的预测步骤仅做一次油藏流体模拟，无法满足要求，故需要对ES-MDA算法进行优化改进。

首先，建立起油藏流体和地震波阻抗的耦合模拟与油藏地质参数的数学关系。油藏流体数值模拟不仅需要输出井口动态生产数据，为了进行地震属性的正演模拟，也需要输出油藏状态函数。油藏模拟器的计算过程可以通过式（22）表示。

地震属性的正演模拟可根据前面介绍的地震纵波波阻抗的计算方法完成，其计算可以通过式（23）简要表示。

式中I（x）——地震纵波波阻抗正演模拟程序；

Ip——模拟计算得到的地震纵波波阻抗数据。

将油藏流体模拟器输出的油藏状态参数作为地震纵波波阻抗正演模拟的输入数据，再将计算得到的地震纵波波阻抗数据与油藏流体模拟器输出的井口动态生产数据相结合，即可实现油藏流体—地震属性的耦合模拟，其数学模型可以通过式（24）表示。

GI(m)——油藏流体—地震属性耦合模拟程序，其程序架构图如图3所示。

图3 油藏流体—地震属性耦合模拟程序架构图

使用地震纵波波阻抗数据和井口生产数据进行历史拟合，还需要对这两类动态数据在空间和时间两个维度上进行匹配。一方面，由于四维地震数据采集的时间间隔往往是1年以上，远远高于井口生产数据采集的时间间隔；另一方面，地震纵波波阻抗数据的空间密度要远远高于井口的生产数据。因此，需要对四维地震数据与井口生产数据进行数据匹配。

基于减小数据体量、节省计算时间、保证运算数据稳定性这3个原则，数据匹配工作需要实现以下3个目标：（1）将地震波阻抗的空间维度降低。通过在空间上设计合理的稀疏采样的方法，将地震纵波波阻抗数据的空间体量降低2～3个数量级，从而达到减小数据体量的效果。（2）保持井口生产数据的时间维度。将地震纵波波阻抗数据在时间维度上插入到井口生产数据上，以保证ES-MDA算法所吸收的数据有足够反映油藏动态变化的时间密度，以及保证后验概率更新的可靠性。（3）不同类型实际测量数据的数值大小不同，数据吸收算法会把数值大的数据给予更高的影响权重，且如果数值的数量级差距过大，也会造成算法的不稳定。为此，需要设计实测数据矩阵，使得不同类型的实际测量数据的数量级一致，如式（26）所示。

式中Dobs—— 匹配后的包含井口生产动态和四维地震纵波波阻抗的实际测量数据；

标量fd，fI—— 分别是井口生产数据和稀疏化处理后的四维地震纵波波阻抗数据的匹配系数；

dobs—— 包含原油产量、产水量、注水井井底流压等井口生产动态的实际测量数据。

根据式（26），数据匹配的关键在于匹配系数设计。良好的匹配系数设计会使得需要拟合的实际测量数据dobs中的每个非零元素数值的数量级相当，从而实现在历史拟合过程中，每种类型的实际测量数据都能得到充分利用和拟合，也保证了程序的数值稳定性。之后，对ES-MDA算法中的调整实际测量数据步骤进行优化改进，结合四维地震与井口动态的数据匹配方法，根据式（27）对匹配后的实际测量数据加上白色噪声。

式中Duc—— 加误差后的包含井口生产动态和四维地震波阻抗的数据。

通过上述对预测步骤和调整实测数据步骤的改进，形成了ES-MDA算法的改进型，其算法流程图如图4所示。

图4 ES-MDA算法改进型的流程框图

2.4 新方法流程

综合上文所述，笔者提出了一套基于四维地震波阻抗数据与生产数据的油藏自动历史拟合新方法体系，主要步骤如下：

（1）储层随机建模形成初始油藏地质模型集合。利用测井数据、岩心数据、地震数据等所有可利用的资料，根据储层类型特点，选取合理的储层随机建模方法，建立一整套初始油藏地质模型集合；

（2）井口生产数据预处理。分析所有的油井生产历史数据，去除掉明显异常值，尽量减少产量的高频波动，选取适当的时间间隔。

（3）地震数据与井口生产数据匹配。原始地震数据体的空间密度往往非常高，而时间密度却远低于井口生产数据。按照前面介绍的方法，在空间维度上降低地震数据体量，时间维度上保持井口生产数据的密度，并形成实际测量数据向量dobs。

（4）运行ES-MDA自动历史拟合程序。设定好数据吸收次数Na和加权吸收α，输入初始地质模型集合及数据配合后的实际测量数据dobs，运行ESMDA自动历史拟合程序，不断调用油藏流体模拟程序，修正得到后验地质模型集合。

（5）地质模型验证与预测数据分析。ES-MDA油藏自动历史拟合得到后验地质模型，需要仔细对照地质模型是否符合所在储层类型成藏过程的地质形态，再对后验地质模型运行油藏流体模拟程序，查看预测的生产数据是否与实际测量的动态数据相符合。

如果历史拟合结果不符合预期，则需要重回第二步和第三步，重新进行数据匹配和ES-MDA自动历史拟合程序的参数设计，直到得到较好的结果。基于四维地震数据的油藏自动历史拟合程序架构图见图5。

图5 基于四维地震数据的油藏自动历史拟合程序架构图

3 应用实例

3.1 历史拟合案例设置

本文研究区以布伦特油田某海相砂岩油藏为原型，建立油藏模型作为真实油藏地质情况的参照。该油藏为简单背斜构造，边界由断层和地质层位构成，内部无明显断层，油藏边界为密闭断层，不含活跃边底水，缺少外界能量供给。对于海相沉积环境，通常在海平面上升时期进行沉积，在海平面下降时则没有沉积物沉降，但海的作用影响到海岸，从而会造成沉积边界形成[23]。油藏具有较强的非均质性，存在未知的流体高渗通道。

参照油藏模型由40×120×20个正交立方网格组成，每个网格尺寸是75m×75m×4m，其中78720个网格为有效储层网格，其余为非储层死网格。其中，孔隙呈高斯分布；海相砂岩储层渗透率基本呈对数高斯分布（图6）。

图6 对数渗透率参照模型与频率分布直方图

通过序贯高斯模拟建立100套孔隙度和对数渗透率的初始地质模型。对数渗透率初始模型的平均值和其中4个初始模型如图7所示。

图7 对数渗透率初始地质模型

由于序贯高斯模拟使用的是随机路径条件下根据变差函数的克里金的估算，故初始地质模型的不确定性主要表现在估算点的空间位置上。初始地质模型的平均地质参数在空间分布上较为均匀，难以判断油藏的地质特征；而对单独的初始地质模型则有明显的地质属性空间分布规律。这说明未经历史拟合的初始地质模型在空间上存在很大的不确定性。为减小地质模型的不确定性，需要纳入油藏动态数据，在此基础上进行贝叶斯后验概率估计，从而提高地质模型的置信度。

为保持油藏压力，在油田生产初始阶段即采用注水采油方式生产，有5口注水井和5口采油井，以反五点法方式注水，P1、P2、P3、P4和P5为采油井，Ip1、Ip2、Ip3、Ip4和Ip5为注水井。已知该油藏18年的井口生产动态数据，如生产井产油量、产水量，以及注水井的井底流压。对该油田海底设置永久地震检波器，每2年进行一次四维地震采集，经过处理与解释，获得纵波波阻抗数据体。

下面通过只使用井口生产数据和综合使用四维地震数据与井口生产数据进行油藏自动历史的研究来检验基于四维地震数据的油藏历史拟合方法的可靠性，探讨四维地震数据对于修正油藏模型的重要意义。

3.2 方法1：井口生产数据的油藏历史拟合

这里仅对井口生产数据进行拟合，用以验证ESMDA历史拟合算法程序的可靠性。拟合的井口生产数据有3种：（1）油井产油量；（2）油井产水量；（3）注水井井底流压。经大量试验确定[24]，ES-MDA历史拟合中，设置数据吸收次数Na=4，每次数据吸收的α分别为α1=α2=α3=α4=4。

井口生产数据历史拟合结果如图8所示。图中绿色细实线为初始地质模型通过油藏流体模拟预测的油田生产数据，包括油田日产油量、油井P1的日产油量和日产水量、注水井I1的井底流压。图中，红色粗虚线为实际测量的生产数据，黑色细实线为经过历史拟合修正地质模型后，在油藏流体模拟器上预测得到的生产数据结果。

图8 方法1：井口生产数据历史拟合结果

图9为修正后的对数渗透率模型。可以看出，对于生产井的产油量和产水量拟合结果良好，黑色细实线基本围绕在红色粗虚线附近，且分布远比绿色细实线更集中，说明经ES-MDA算法修正地质模型后，对后验地质模型进行油藏流体模拟预测，基本上拟合了实际测量的井口生产数据。

图9 方法1：修正后的对数渗透率模型

3.3 方法2：四维地震数据与井口生产数据的油藏历史拟合

在前面针对井口生产数据的油藏历史拟合基础上，将通过智能地震处理解释技术[25]获得的四维地震的纵波波阻抗数据也纳入历史拟合，真正实现了四维地震数据与井口生产数据的油藏自动历史拟合研究[16]。

首先，分析各类型动态数据进行数据匹配，设置实际测量数据匹配矩阵，如式（28）所示。将油井产油量fQ、油井产水量fW、注入井井底流压fB、地震纵波波阻抗的匹配系数fI分别设为1、10、0.005和0.1，从而得到数据匹配后的实际测量数据向量Dobs。

运行基于四维地震数据的油藏自动历史拟合程序，井口动态生产数据的拟合结果如图10所示。对比图8，图10中的各项生产数据拟合结果要明显更好。图中黑色细实线不仅都围绕在红色粗虚线附近，而且其分布更集中于实际测量数值。说明同时拟合四维地震数据与井口生产数据，要比仅对井口生产数据进行历史拟合数值收敛性增强很多，对后验地质模型进行油藏流体数值模拟的预测结果更接近于真实测量数据，间接表明了后验地质模型的置信度的提高。修正后的对数渗透率模型如图11所示。

图10 方法2：井口生产数据历史拟合结果

图11 方法2：修正后的对数渗透率模型

将后验地质模型与方法1的结果对比，不难发现方法2的后验地质模型集合的平均地质参数分布拥有更明显的特征，更接近于真实的油藏情况。且单个集合中的子集地质模型相互之间的差异也减小更多。说明加入地震数据进行历史拟合的后验地质模型的置信度更高的同时，也证实4次数据吸收的次数已足够充分。增加地震数据作为修正条件之后，ES-MDA算法的数据吸收次数明显减少。

3.4 不同方法后验地质模型置信度定量分析

为了定量分析对比不同油藏历史拟合方法的后验地质模型置信度，设计了两个统计学指标，用于对比这两种方法的后验地质模型进行油藏流体数值模拟预测结果的准确性[26]。

后验地质模型集合的平均地质参数模型的油藏流体数值模拟的预测相对误差e，定义如下：

式中，预测相对误差可以通过e的下角标来表示。油田生产的最后一天，对于油井产油量、产水量的预测相对误差，以eoPi和ewPi表示；对于注水井井底流压的预测相对误差以ebhpIi来表示。其中下角标中的i表示模型编号，P和I则分别表示采油井与注水井。

（2）后验地质模型集合的平均地质参数模型的油藏流体数值模拟的预测方差，定义如下：

同样，式中预测方差可以通过σ2的下角标来表示，下角标的标准准则与预测相对误差e相同。

表1和表2分表列出对于预测相对误差e和预测方差σ2的值。对比表中两种方法的预测误差和方差结果可知，方法2得到的后验地质模型运行油藏流体模拟程序之后的各项预测指标，明显好于方法1的预测，这也验证了四维地震数据对于油藏历史拟合结果的积极意义。

表1 两个方法的预测相对误差指标统计表 （单位：%）

表2 两个方法的预测方差指标统计表

4 结论

为了充分利用四维地震监测技术提供的覆盖全油藏的状态信息数据，解决现有油藏自动历史拟合方法仅能处理井口生产动态数据造成的数据利用程度低、拟合效果差、油藏模型更新可靠性低等问题，本文基于贝叶斯理论方法，优化改进了ES-MDA算法，提出一种基于四维地质历史拟合的新方法，大幅提高了油藏自动历史拟合的质量，并在某海相砂岩油藏的历史拟合应用中取得了良好效果。

本文通过理论推导和应用研究，得到以下认识：

（1）针对四维地震数据的空间密度要远高于井口生产动态数据，而时间密度远低于井口生产动态数据的问题，提出了四维地震数据的粗化与稀疏化处理方法，并通过设计匹配系数形成实际测试数据矩阵，实现了四维地震数据和井口生产动态数据的匹配。

（2）通过建立起油藏流体和地震波阻抗的耦合模拟与油藏地质参数的数学关系，形成油藏流体—地震属性耦合模拟程序，并结合数据匹配方法，改进优化ES-MDA算法，建立了基于四维地震数据的历史拟合新方法。

（3）同时拟合四维地震数据与井口生产动态数据，比仅对井口生产数据进行历史拟合所得到的数据拟合结果更好，后验地质模型的置信度更高，且增加了四维地震数据作为修正条件之后，也能明显减少ES-MDA算法的数据吸收次数，体现了四维地震监测的价值。四维地震与井口生产动态数据的油藏历史拟合具有实际意义和应用前景。