加速油藏数值模拟的POD模型降阶方法∗

曹静，赵辉，喻高明，胡文明

(1.长江大学石油工程学院，湖北武汉430100；2.中原油田普光气田净化厂，四川达州636100)

油藏数值模拟是油藏开发研究的重要手段之一.它能使工程师更好地了解油藏的各种物理性质和流体在其中的流动规律，以便作出正确评价，确定合理开发方案和提高采收率.传统油藏模拟器将运动方程、状态方程和连续方程构成的偏微分方程组离散化，转化成一组非线性代数方程组，再采用迭代法求解.然而，由于实际油田建立的油藏模型可能包括几十万或几百万个网格，需要求解的方程数目非常大，所以，传统油藏模拟器的数值解决方法非常耗时.特别的，当把油藏模拟器应用于闭环油藏管理时[1−5]，其中的生产优化和历史拟合过程都需要反复运行模拟器，导致它的计算成本大大提高.因此，在保证有足够精度数值解的情况下，如何大幅度加速油藏模拟运算速度是亟待解决的问题.

模型降阶技术最早出现在自动控制和电路系统领域，它的任务是减少系统状态空间向量的维数，同时保持系统的输入、输出特性.本征正交分解法（Proper Orthogonal Decomposition，简记为POD）是非线性系统模型降阶中最为广泛使用的方法.该方法所需的计算代价较低，具有一定稳定性，同时又能保持原始系统的一些基本性质.目前POD方法已经成功地应用于计算流体力学[6]、结构力学[7]、气象学[8]及数字信号处理[9]等领域.本文将POD模型降阶方法应用于油藏模拟器，可大幅度降低油藏模型维数，从而减少计算时间，提高运算速度.

1 油藏模型的控制方程

本文利用空间上的离散将油藏模型的数学方程转换成状态空间方程的形式，以便说明POD方法的降阶过程.采用二维油水两相油藏模型，假设油、水两相不发生物质交换、过程等温、流体微可压缩，应用质量守恒方程和达西定律可得[10]

式中：K为渗透率张量；µ为流体粘度；kr为相对渗透率；p为压力；g为重力加速度；d为深度；ρ为流体密度；φ为孔隙度；S为流体饱和度；t为时间；q000为源汇相，表示单位体积的流速；上标i∈{o,w}分别表示油相和水相.等式（1）中有四个未知量po,pw,So,Sw，利用辅助方程（2）、（3），消去pw,So，使等式中只包含po,Sw两个状态变量.

式中：pc(Sw)为油水两相的毛管压力.

我们这里考虑较简单的情况，忽略重力和毛管力.利用五点块中心有限差分格式在空间上进行离散，可得非线性一阶微分方程[11].

式中：向量p和s分别为网格中心油相压力po和水相饱和度Sw；V为累积矩阵；T为传导矩阵；F为分流量矩阵；向量qwell,t为井的油水总流量.

定义状态向量x，输入向量u，输出向量y分别为

等式（4）可写成状态空间方程的形式[11]

在控制系统中，Ac称为系统矩阵，Bc称为输入矩阵，C称为输出矩阵，D称为直接传递矩阵.由于矩阵V，T，F中的元素都是状态变量x的函数，则该系统为非线性系统.

下面采用全隐式方法求解控制方程（8），用∆x/∆t近似导数dx/dt，则有

设xk表示tk时刻的状态变量值，xk−1表示tk−1时刻的状态变量值，则有

令fk=xk−1+∆tfc(uk,xk)，可得

由于等式（12）左右两端均含xk，则它为非线性方程组，一般采用牛顿迭代法求解.但当变量xk的维数很高时，求解过程非常耗时，计算代价较大，为此，我们将POD模型降阶方法应用于油藏模拟，降低xk的维数，进而提高运算速度.

2 POD模型降阶方法

POD方法也被称为经验正交函数法和Karhunen–Lo`eve分解法.该方法利用给定的理论或者实验得到的数据来构造一组最优正交基，达到降阶目的.下面我们给出该方法的基本原理及在油藏数值模拟中的具体应用，关于POD方法更详细的介绍可参考文献[12-15].

2.1 POD方法的基本原理[16]

假设φ1,φ2,…,φn是Rn中数据集{x1,x2,…,xm}的标准正交基向量，则对任意样本向量xi∈Rn，就可以由基底按如下形式线性表出

其中cji=.若选择前k个基向量来逼近样本向量xi，则有

对于获得的数据集{x1,x2,…,xm}，POD方法的目标是找到一组最优基向量{φ1,φ2,…,φn}，以使误差ε2达到最小.引入Lagrange因子uij(i,j=k+1,k+2,…,n)，并且构造Lagrange函数

（16）式两端对φj求偏导，得到

其中Un−k=[uk+1uk+2…un].为求最优解，让（18）式左端为0，则有

其中Λ是对角矩阵.对（19）式两端右乘矩阵P，并结合（20）式，就有

由（21）式可知，矩阵Λ的对角元素由矩阵XXT的特征值λi构成，矩阵Φn−kP由矩阵XXT相应于λi的特征向量构成.实际上，矩阵Λ的对角元素是矩阵XT的奇异值σi的平方.考虑到正交矩阵在Frobenius范数k•kF下是保范的，因此得到

为了获得最小误差，矩阵Λ的对角元素需要选择矩阵XT的最后n−k个奇异值.这样，就有minε2=

由以上讨论可知，当选择矩阵XT的右奇异值向优性条件，而且也获得了最小误差，其最小误差为矩阵XT的最后n−k个奇异值的平方和.

通常在输入、输出系统中，较大的特征值对应系统的主要特征.为此，我们选择的k个POD基向要能代表原始向量较多的特征相关矩阵按降序排列的特征值.

假如I(k)≥d%，则k个POD基向量保持原始数据集的d%特征信息，其中k满足

2.2 POD方法在油藏模拟中的应用

为了构造POD基向量，首先要运行全阶油藏模拟器（也称为训练过程），并保存每个时间步的状态向量x（也称为快照，包括所有网格的油相压力po和水相饱和度Sw），由于压力和饱和度具有不同的物理性质，我们分别用矩阵Xp，XS保存po，Sw（以下简记为p和S）.

假设油藏模型的网格数为N，则矩阵Xp和XS中的每个向量（上标i表示快照的数目）都是N维，而系统的状态向量x的维数为n=2N.快照获得后，需要计算快照的均值

并且数据矩阵Xp和XS中的每个快照都减去均值

将式（27）去掉下标l代入（12）中，可得

等式两端同时左乘ΦT

这就实现了对原系统降阶的目的，状态变量的维数从n=2N减少为l=lp+ls.再采用牛顿迭代法求解非线性方程组（29），由于状态变量的维数较低，可大大减少计算量.

3 算例验证

本文采用文献[11]中的算例.该算例描述的是一个二维油水两相各向异性的方形油藏，网格划分为21×21，渗透率和孔隙度的分布如图1、2所示.油藏模型相关参数：油层厚度h=2 m，网格的长、宽∆x=∆y=33.33m，原油粘度µo=5.0×10−4Pa·s，地层水粘度µw=1.0×10−3Pa·s，综合压缩系数ct=3.0×10−9Pa−1，原始地层压力pi=30×106Pa，井孔半径rwell=0.114 m，油相端点相对渗透率=0.9，水相端点相对渗透率=0.6，油相Corey指数no=2.0，水相Corey指数nw=2.0，残余油饱和度Sor=0.2，束缚水饱和度Swc=0.2.采用反五点法井网生产，中心一口注水井，四角为四口生产井，忽略重力和毛管力影响.

图1 油藏模型的渗透率分布

图2 油藏模型的孔隙度分布

我们采用Matlab编辑的全隐式处理的数模程序模拟该算例，修改源代码实现POD模型降阶过程，验证该方法的有效性.

3.1 训练过程

生产井给定井底压力为295×105Pa，注入井给定井底流量为0.0015m3/s，运行全阶模拟器模拟1400天，保存66个时间步的结果.压力矩阵保留37个奇异值，饱和度矩阵保留38个奇异值，最后获得POD的基矩阵Φl的维数为2N×l，其中l=37+38=75.这就意味着全阶模拟器要求解2N=882个未知变量，而降阶后只需求解75个变量.

训练过程中，全阶油藏模拟器与应用POD方法的降阶模拟器计算结果对比见图3，4.

在文中我们采用平均相对误差来衡量近似值的精度.如每口生产井产油量的平均相对误差定义为

式中：i表示时间步；nt表示时间步总数；表示生产井m第i时间步全阶模拟器的产油量；表示生产井m第i时间步POD降阶模拟器的产油量.同理也可定义每口生产井产水量的平均相对误差.

图3 四口生产井的产油量对比（训练过程）

图4 四口生产井的产水量对比（训练过程）

训练过程，四口生产井产油量、产水量的平均相对误差如表1.

表1 四口生产井产油量、产水量的平均相对误差（训练过程）

以上结果显示，在训练过程中，降阶模拟器与全阶模拟器四口生产井的产油量和产水量基本完全相同，产生的平均相对误差都非常小，但模拟时间提高了近3倍，全阶模拟器的运行时间为35.527 s，而降阶模拟器的运行时间为12.019 s.

3.2 预测过程

下面我们用训练过程获得的基矩阵Φl来验证POD降阶模型的预测能力.此时生产井的井底压力改为285×105Pa，注入井的井底流量保持不变，全阶模拟器与降阶模拟器的计算结果对比见图5、6.

图5 四口生产井的产油量对比（预测过程）

图6 四口生产井的产水量对比（预测过程）

预测过程，四口生产井产油量、产水量的平均相对误差见表2.

表2 四口生产井产油量、产水量的平均相对误差（预测过程）

上面结果显示，当预测过程生产制度与训练过程生产制度不同时，降阶与全阶模拟器输出结果的平均相对误差有所提高，但仍控制在5%的合理范围内.此时，模拟时间也提高了将近3倍，全阶模拟器的运行时间为35.851 s，降阶模拟器的运行时间为12.404 s.

4 结论

（1）将油水两相渗流的偏微分方程，利用五点块中心有限差分格式在空间上进行离散，转化成状态空间方程的形式；

（2）POD模型降阶方法应用于油藏模拟器，可大幅度减少油藏模型维数，实例显示在保证一定精度的条件下，模拟器的运算速度提高近3倍，验证了该方法的有效性；

（3）当训练与检测过程的生产制度不同时，降阶模拟器的计算结果平均相对误差有所提高，并且两个过程的生产制度差别越大，产生误差也越大，所以POD模型降阶方法还需进一步改进.