大庆油田从聚合物的基本物理化学性质及其驱油机理出发，形成了一整套聚合物驱油技术，包括聚合物驱室内实验技术、筛选技术、油藏工程技术、效果评价技术和聚合物驱配制工艺、注入工艺、采油工艺等地面工程技术以及聚合物驱综合管理等技术^[1-2]。作为数值模拟重要环节之一的历史拟合，其精度直接影响着开发方案的制定与调整。对比水驱数值模拟，聚合物驱涉及到的参数更多，物理化学反应的描述方程更加复杂，传统手工历史拟合方法消耗时间较长，在文献调研中，鲜有学者进行化学驱油参数的数据同化^[3-6]。现有的化学驱油模拟软件大多不具备自动历史拟合功能。基于EN⁃ KF 方法、ESMDA 方法和大庆油田自主研制的油藏数值模拟器 CHEMEOR，对化学驱自动历史拟合技术进行了研究，有效利用油田实际开发指标，实现了油藏地质参数和化学驱参数的自动拟合。

1 聚合物驱自动历史拟合数学模型

历史拟合过程中拟合参数包含 2 种状态变量，即所有不确定的油藏模型参数（静态和动态）及生产数据，具体包括：①油藏模型的静态参数，通常不随时间变化或相对变化较小，如地层渗透率、孔隙度、油水界面的深度及表征相对渗透率曲线的参数等。②油藏模型的动态参数，在生产过程中随时间不断变化，包括压力、饱和度及针对聚合物驱的化学剂性质参数（黏浓关系、吸附量等）。③生产数据，即观测数据对应的数值模拟计算数据，包括井底流压、含水率、产油量等^[1-2]。

1.1 拟合参数

地层渗透率 渗透率决定流体流速，忽略渗透率在各方向的相关性，则达西速度为：

其中：

为避免 K_x，K_y，K_z可能出现负值，将其转换为对数坐标系下的值，即将其同化为lnK_x，lnK_y，lnK_z。

孔隙度与初始含水饱和度 孔隙度和初始含水饱和度的关系式为：

孔隙度和初始含水饱和度的取值范围均为[0， 1]，为避免参数越界，将其转换为对数坐标系下的值，即将其同化为ln（-lnϕ）与ln（-lnS_w）。

相对渗透率曲线参数 相对渗透率曲线参数包括高、低毛管数下的残余油饱和度、相对渗透率端点值和相对渗透率指数值^[7-14]，其表达式分别为：

在模拟器中，低毛管数下的相对渗透率以列表的形式给出，不再同化。在高毛管数下残余油饱和度的取值范围为[0.2，1]，为避免同化过程中残余油饱和度出现越界，将其同化为 ln[-ln（ $S_{lr}^{\text{high}}$/0.2）]，包括水油两相。相对渗透率端点值的取值范围为 [0，1]，为避免同化过程出现越界，将其同化为 ln（-ln$K_{rl}^{0_{}^{\text{high}}}$）。类似地，由于相对渗透率指数的取值范围为^[1，3]，将其同化为$\mathrm{l}\mathrm{n}\left\{-\mathrm{l}\mathrm{n}\left[\left(n_{\mathrm{r}\mathrm{l}}^{0^{\text{high}}}-1\right)/2\right]\right\}$。

聚合物黏浓曲线与剪切流变 聚合物黏浓关系式为：

高分子聚合物溶液都具有流变特征，其黏度受剪切速率影响，利用 Meter 方程表示聚合物溶液的黏度与剪切速率的函数关系为：

将A_p1，A_p2，A_p3，S_p参数进行同化。

渗透率下降系数 渗透率下降系数表达式为：

其中：

将b_rk和c_rk进行同化。

可及孔隙体积 可及孔隙体积出现在聚合物组分方程中：

由于可及孔隙体积的取值范围为[0，1]，为避免越界，将其同化为ln（-lnE）。

聚合物吸附 聚合物吸附采用 Langmuir 公式^[8]。

有效含盐量 有效含盐量的计算公式为：

同化参数为β_P。

计算模型的不确定性来源于以上拟合参数。将油水两相模型和组分模型离散后，可得拟合参数的状态向量集合，其描述为：

其中，$K\in {\mathbb{R}}^{3N_{\mathrm{c}}}，\varphi \in {\mathbb{R}}^{N_{\mathrm{c}}}，\alpha =\left(A_{\mathrm{p}1}，A_{\mathrm{p}2}，A_{\mathrm{p}3}，b_{\mathrm{r}\mathrm{k}}，c_{\mathrm{r}\mathrm{l}}，E，a_{\mathrm{p}}\right.$，$\left.b_{\mathrm{p}}\right)\in {\mathbb{R}}^8$。

1.2 生产数据

实际生产数据包括：①定量注入井的井底流压。②定压注入井的注入速度。③定量生产井参数，包括井底流压（每口井一个数据，每月提供一个数据）、含水率（每口井一个数据，每天都有数据）、组分浓度（软件已预留接口，但实际不测）。④定压生产井参数，包括产液量（每天每口井一个数据）、含水率（每天每口井一个数据）、组分浓度（软件已预留接口，但实际不测）。⑤矿场数据，包括矿场的产液量和含水率^[15-20]。将实际生产数据写成生产数据状态向量，其集合描述为：

历史拟合中需要调整渗透率、孔隙度和聚合物驱物性参数，将这些同化参数作为系统状态变量，渗透率和孔隙度的不确定性不随时间变化，只取决于初始变量，扩展后的模型为：

在模拟计算时，通常利用 K，ϕ，α 在 k-1 时间步的值更新状态变量 X_k，然后利用同化算法，更新 K，ϕ，α作为第k时间步的值。

2 模型求解算法

2.1 ENKF方法原理

利用随机方法对模型初始状态向量各拟合参数进行随机撒点，组成初始参数集合，可描述为：

利用离散化的参数模型计算其数学期望，其描述为：

应用数学期望来计算先验误差，其计算数学方程为：

其中：

采用集合卡尔曼滤波原理，利用 k-1 步参数集合及实际观测数据 Z 对第 k 步参数值进行更新，且假定测量误差和模型误差满足高斯分布，为保证集合成员之间的变化性，通常对观测数据加上随机扰动，即在实际观测数据（实际生产数据）Z的基础上，各个样本（集合）加入噪音作为观测值，即：

加入噪音后参数更新的数学描述方程为：

其中：

在实际模拟计算时不需要计算 R（ k|k-1）和关联矩阵 M_k，而是把模拟器当做黑盒直接计算其结果。在 Kalman Gain 的计算中，需进行矩阵求拟运算。

利用生产数据对模型计算参数进行数据同化后，计算其后验误差，R（ k|k）的计算数学方程为：

其中：

历史拟合的过程就是通过不断重复上述计算过程更新拟合参数，从而降低后验误差，最终达到历史拟合的目的。

2.2 ESMDA方法原理

在ENKF方法中，状态变量包括含水饱和度、压力、组分质量浓度，然而在实际中，这些数据并不能测量得到。而且，如果在实际生产数据中对这些变量进行校正，会造成整体系统中的物质不守恒。因此不再考虑这些状态变量，只考虑地质参数和化学剂参数。在实际拟合过程中ENKF方法处理非线性问题能力较弱，在这种情况下，通过改变观测数据的扰动强度将原同化方法转化为 ESMDA 方法，加强了其处理非线性问题的能力^[21-26]。该方法具体实现分为5步：①初始化参数，根据先验参数的先验估计产生集合：

②运行模拟器，计算模型观测值：

③扰动实际观测值：

$w_k^i\sim N\left(0，\sqrt{{\mathit{\beta }}_k}\mathit{Q}（\mathit{k}）\right)$可通过正态分布产生。

若采用 ENKF 方法，则 β_k = 1。若采用 ESMDA 算法，β_k满足：

④更新参数为：

K（ k）满足：

协方差矩阵具体计算式包括：

${\mathit{Q}}_{\mathit{y}}（\mathit{k}）$为模拟观测的协方差矩阵，其表达式为：

⑤如果需要迭代（k<N_a），则进入步骤②继续运算，否则，结束运算。

在本研究中，主要采用 ESMDA 方法处理模型非线性特征，通过不断降低系统误差从而实现参数的自动同化功能。

3 自动历史拟合程序实现与架构设计

程序模块包括：数据读入模块（同化参数的读入模块和观测数据的读入模块）、时间步长控制模块、数据同化模块、误差评估计算模块、数据的整理与输出模块。软件设计流程如图1所示。各功能模块及关系见图2。

参数先验估计数据读入模块 本模块读入拟同化参数的先验估计，包括期望值与不确定性（方差），用户需提供 CHEMEOR 数据流，其中的地质参数或化学参数为工程师认为的最可能的参数，即期望值。用户还需提供拟同化参数的不确定性，作为 AUTOM的数据流。AUTOM的数据流的输入格式与 CHEMEOR 数据流的格式完全相同，但是所表示的数值的意义不同，前者是方差，后者是期望值。AU⁃ TOM 根据方差和期望值，利用正态分布产生随机数，作为同化集合。需要注意的是，AUTOM 可以实现多套化学参数的拟合功能，要求数据格式与 CHEMEOR中化学参数的输入格式保持一致。

观测数据读入模块 本模块实现观测数据Z及其不确定性的读入。AUTOM 可以读入全部的观测数据，完成同化。由于在实际中观测数据可能很多，AUTOM 也可以分段读入观测数据。当读入一段观测数据后，AUTOM 可以实现参数的同化，当在本段的迭代次数完成后，再读入下一段观测数据，进行同化。关键字UPDATE将观测数据分段。

AUTOM 与 CHEMEOR 数据对接模块 本模块实现 AUTOM 与油藏数值模拟软件 CHEMEOR 的对接。AUTOM在CHEMEOR数据流的基础上，根据随机产生的同化参数，改变要同化的参数，生成新的 CHEMEOR 数据流，并增加关键字 AUTOM，由 CHEMEOR 实现模拟计算，并输出模拟观测数据。 CHEMEOR 模拟计算每个样本，根据数据流中的 AUTOM关键字，输出AUTOM所需的观测数据Y，并由AUTOM收集，以实现同化计算。本模块实现AU⁃ TOM 对 CHEMEOR 的并行化调用，即在 AUTOM 完成 N_E 个样本生成之后，可调用多个线程进行 CHEMEOR模拟计算，在所有线程完成模拟计算后，由 AUTOM 统一收集数据，进行同化运算。用户根据可利用的计算资源设置线程数。

参数同化模块 本模块实现同化计算，可利用的数据同化方法包括 ENKF 方法和 ESMDA 方法。其中涉及矩阵的求拟运算。在算法中，当矩阵 ${\mathit{Q}}_{\mathit{y}}（\mathit{k}）+\sqrt{{\beta }_k}\mathit{Q}（\mathit{k}）$不可逆时（矩阵中出现 0 值），利用 SVD 分解，求解该矩阵，将特征值较小的方向去掉，实现该矩阵的伪拟运算。用户需要设置舍掉的特征方向，可通过设置 Tolerence 实现，即当一个特征值的绝对值小于最大特征值的Tolerence倍时，视为该特征方向对系统的解影响很小，不再考虑^[27-28]。

误差计算与参数输出模块 在数据同化中，衡量同化效果的一个重要指标是误差。在 AUTOM 中，计算模拟观测与实际观测之间的差，衡量同化性能，并用来观测同化过程是否出现过拟合。其计算式为：

参数的后验估计可通过 AUTOM 输出数据获得，包括参数的期望值和方差。同时 AUTOM 输出每次迭代后每个样本的模拟观测值。

4 应用实例

4.1 试验区基本概况

选取 M 区块聚合物驱试验区，试验区利用目的油层聚合物驱井网进行井网重构，五点法面积井网，注采井距为175 m，共有注采井49口，其中，注入井 25 口，采出井 24 口。试验区含油面积为 1.18 km²，地质储量为 176.62×10⁴ t，孔隙体积为 305.70× 10⁴ m³，开采目的油层共发育 6 个沉积单元，平均单井射开砂岩厚度为15.0 m，有效厚度为10.8 m，平均有效渗透率为 0.522 D。数值模拟模型（图3）网格 43×40×6 共 10 320 个，dx，dy 网格尺寸分别为 31.88 和34.86 m。

4.2 聚合物驱自动历史拟合效果

在选定试验区基础上，对模型进行了聚合物驱动态计算，数值模拟计算结果（图4，图5）与实际生产动态差距较大。

在上述研究的基础上，选取了聚合物参数、初始含水饱和度等 130 个拟合参数，利用自动历史拟合软件对 120 个样本集合进行了 15 次迭代数据同化拟合计算，拟合结果如图6—图8所示。拟合结果表明，所研制的自动历史拟合软件可以降低模型计算结果系统误差（图9），全区含水率趋势基本符合历史动态，全区累积产油量误差为 3.8%，全区单井拟合结果较好，单井拟合符合率为81.2%，满足聚合物驱历史拟合要求。

5 结论

多参数聚合物驱数值模拟数学模型可以较好地反映油藏中聚合物驱替的基本规律，通过综合调整聚合物的吸附量、黏浓关系、残余阻力系数、不可及孔隙体积等参数，与油藏实际动态相符合，达到历史拟合目的。

所建立的自动历史拟合数学模型满足算法要求，通过改变观测数据的扰动强度加入符合正态分布的扰动系数将原算法转化为 ESMDA 方法，加强了历史拟合过程中处理非线性问题的能力。

使用改进后的 ESMDA 自动历史拟合方法所研制的自动历史拟合软件具有聚合物驱参数的自动调整功能，实际区块测试结果表明该软件可以降低模型计算结果系统误差，全区含水率趋势基本符合历史动态，全区累积产油量误差为 3.8%，全区单井拟合结果较好，单井拟合符合率为81.2%，满足聚合物驱历史拟合要求。

符号解释

a_p，b_p——动吸附所需拟合参数；

A_p1，A_p2，A_p3，S_p——由实验资料确定的常数；

b_rk，c_rk——输入参数，由实验数据确定；

${\tilde{C}}_i$——第i种物质组分的总体积分数；

C_il ——液相中第i种物质组分的质量浓度，mg/L；

C_p——溶液中聚合物的质量浓度，g/L；

C_SEP——聚合物有效含盐量当量浓度；

C₃，C₄，C₅——CHEMEOR软件组分卡中第3、第4、第5组分，分别代表聚合物质量浓度、阴离子当量浓度、阳离子当量浓度；

C₃（0|0），C₄（0|0），C₅（0|0）——随机产生的聚合物质量浓度、阴离子当量浓度和阳离子当量浓度；

C₄₀，C₅₀——阴离子、阳离子在水相中的当量浓度；

${\tilde{D}}_{iw}$——弥散流量；

E——可及体积分数；

g——重力加速度，m/s²；

i——拟合样本序号；

k——时间步；

K——渗透率，mD；

K_k——k时间步渗透率，mD；

K_rl ——l相相对渗透率；

Kⁱ（ 0|0）——随机产生的渗透率场；

$K_{\mathrm{r}\mathrm{l}}^0$——l相相对渗透率端点值；

$K_{\mathrm{r}\mathrm{l}}^{0^{\text{high}}}$——l相高毛管数下相对渗透率端点值；

$K_{\mathrm{r}\mathrm{l}}^{0^{\text{low}}}$ ——l相低毛管数下相对渗透率端点值；

K_rα——α相相对渗透率；

K_x，K_y，K_z——x，y，z方向的渗透率，mD；

K（k）——卡曼根即最优权重值；

下标l——液相，包括油和水，即o和w；

L（ k|k）——k时间步每个参数样本与期望值之差的向量集合；

L _x（ k|k-1）——k-1时间步每个参数样本与期望值之差的向量集合；

L_y（ k|k-1）——k-1时间步每个数值模拟计算值与期望值之差的向量集合；

M_k——状态与观测的关联矩阵，${\mathit{M}}_{\mathit{k}}\in {\mathbb{R}}^{N_{\mathrm{m}}\times N_{\mathrm{T}}}$；

MSE——均方根误差；

n_l ——l相相对渗透率指数；

$n_{\mathrm{r}l}^0$——相对渗透率指数；

$n_{\mathrm{r}\mathrm{l}}^{0^{\text{high}}}$ ——l相高毛管数下相对渗透率指数；

$n_{\mathrm{r}\mathrm{l}}^{0^{\text{low}}}$——l相低毛管数下相对渗透率指数；

N_a——迭代步数；

N_E——样本数；

N_G——网格数；

N_I ——注入井个数；

N_m——观测变量的维数，满足N_m=N_I +2N_P；

N_obs——观测的个数；

N_P——生产井个数；

N_T——总矩阵维数，根据所需拟合参数计算得出；

$N_{T_l}$ ——毛管数；

O（ X_k）——均方根误差；

p_α——α相压力，MPa；

q_a——由实验数据确定的系数；

P_bh——井底流压，MPa；

P（0|0）——随机产生的压力集合；

q_l ——注入量或者开采量，m³ /d；

Q——产液量，m³；

Q（ k）——实际生产数据的协方差矩阵，$\mathit{Q}（\mathit{k}）\in {\mathbb{R}}^{N_{\mathrm{m}}\times N_{\mathrm{m}}}$；

Q_xy——协方差矩阵；

Q_xy（ k）——参数 X（ k|k-1）与 Y（ k|k-1）的协方差矩阵；

Q_y（ k）——模拟观测的协方差矩阵；

R（k|k-1）——先验误差；

R_i ——第i组分源汇项^[10]；

R_K——渗透率下降系数；

R（ k|k）——后验误差；

R_KMAX——最大渗透率下降系数；

R_α——α相溶解汽油比，m³ /m³；

S_l ——l相饱和度；

S_lr ——l相残余油饱和度；

$S_{l\mathrm{r}}^{\text{high}}$——l相高毛管数下残余油饱和度；

$S_{l\mathrm{r}}^{\text{low}}$——l相低毛管数下残余油饱和度；

S（0|0）——随机产生的饱和度集合；

t——组分运动时间，s；

T_l——实验测定参数；

$w_k^i$——观测噪音，可通过高斯分布产生随机数生成；

W_c——采出液含水率，%；

u_l ——l相达西速度，m/d；

u_α——α相达西速度，m/d；

$\widehat{X}$（ k|k）——数学期望；

$\widehat{X}$（ k|k-1）——离散化参数的数学期望；

$X^i$（ k|k）——k时间步更新后的拟合参数集合；

$X^i$（ k|k-1）——k-1时间步拟合参数集合；

$X^i$（ 0|0）——拟合参数集合；

X_k——拟合参数状态向量，包含压力、饱和度、聚合物质量浓度、阴离子当量浓度和阳离子当量浓度；

X^* _k—— 扩展后拟合参数状态向量，包括 X_k，ln（-lnϕ_k），lnK_k，α_k；

X_k-1——k-1时间步拟合参数状态向量；

X^* _k-1——k-1时间步扩展后的参数状态向量；

$X^{N_{\mathrm{E}}}（k\mid k-1）$——第N_E个拟合参数集合；

Y——数值模拟计算结果状态向量；

Y_k——数值模拟计算结果集合；

$\widehat{Y}$（ k|k-1）——数值模拟计算值期望；

$Y^i$（ k|k-1）——k-1时间步第i个样本数值模拟计算值；

$Y^{N_{\mathrm{E}}}（k\mid k-1）$——第N_E个样本数值模拟计算值；

Y（ k|k-1）——k-1时间步数值模拟计算值；

z——深度，m；

Z——实际观测数据；

Z（ k|k）——实际生产数据；

Z^* （ k|k）——加入扰动的实际生产数据；

Z_k——实际生产数据（产液量、含水率等）；

Z^* _k——加入扰动后的实际生产数据；

μ_α——α相黏度，mPa·s；

下标α——相，α为o，g，w，分别表示油、气、水；

α——需要同化的8个聚合物物性参数；

αⁱ（ 0|0）——随机产生的聚合物物性参数；

α_k——k时间步聚合物物性参数；

β_k——相关系数；

β_P——由实验确定的参数；

$\dot{\gamma }$——多孔介质中流体的等效剪切速率，s^-1 ；

γ_ref ——参考剪切速率，s^-1；

$\tilde{\mu }$——聚合物溶液的本征黏度，mPa·s；

μ_o，μ_w——油相、水相黏度，mPa·s；

μ_p——聚合物溶液的黏度，mPa·s；

${\mu }_{\mathrm{p}}^0$——零剪切速率下聚合物溶液的黏度，mPa·s；

${\xi }^i$（ 0|0）——初始状态向量集合；

${\xi }^i$（ k|k）——第i个拟合参数的状态向量；

${\xi }^{N_{\mathrm{E}}}$（ k|k）——第N_E个拟合参数的状态向量；

${\xi }^{N_{\mathrm{E}}}$（ k|k-1）——k-1 时间步第 N_E个拟合参数的状态向量；

ρ_i ——i组分的密度，g/L；

ρ_l ——l相密度，g/cm³；

ρ_α——α相密度，g/cm³；

ϕ——孔隙度；

ϕⁱ（ 0|0）——随机产生的孔隙度场；

ϕ_k——k时间步孔隙度；

ϕ_p——可及孔隙体积，m³；

$\varnothing$——孔隙体积，m³；

Φ，Φ^* ——相应油藏模拟求解算子；

ω_k——观测噪音，可通过高斯分布产生随机数生成。

参考文献