采用基于KPOD的模型降阶方法提高油藏模拟速度的研究①

曹 静赵 辉 喻高明

(长江大学石油工程学院 武汉 430100)

采用基于KPOD的模型降阶方法提高油藏模拟速度的研究①

曹 静②赵 辉 喻高明

(长江大学石油工程学院 武汉 430100)

应用模型降阶技术进行了提高传统油藏模拟器运算速度的研究。考虑到该技术采用本征正交分解(POD)方法能够加速模拟器运算，但对系统的输入参数较敏感，降低模拟的精度和效率，采用了基于Krylov子空间和POD的KPOD方法，以便利用Arnoldi的矩匹配性质和POD的数据泛化能力减少POD方法对模型输入参数的依赖。油藏模拟实例表明， KPOD方法在计算速度和精度上都优于POD方法，验证了该方法的有效性和实用性。

油藏数值模拟， 模型降阶， 本征正交分解(POD)， KPOD方法

0 引 言

油藏数值模拟是现代油藏工程中不可缺少的研究工具。它普遍用于油藏开发设计、动态预测、油层参数识别、工程技术优化设计以及重大开发技术政策的研究等。传统油藏模拟器通过有限差分将偏微分方程组离散，转化成非线性代数方程组，再采用迭代法求解。然而实际油藏，数值模拟地质模型网格数可能高达几十万乃至几百万，采用传统油藏模拟器进行模拟非常耗时，严重制约了数值模拟技术在我国大油藏特别是陆相非均质严重油藏的应用，同时也造成实际油藏的生产历史拟合以及后续的生产优化过程困难重重，因此，如何大幅度提高油藏模拟运算速度是亟待解决的问题。

近年来，模型降阶(model order reduction，MOR)技术在确保精度的条件下，对减少系统的运算量显示出了一定的潜力[1-4]。该技术可将一个较大系统转化为一个近似的较小系统，以便达到降低大型复杂系统的理论分析难度和减少数据运算量的目的[5]。本征正交分解(proper orthogonal decomposition，POD)法是非线性系统模型降阶中使用最为广泛的方法，该方法已被应用于很多领域。目前，它也被应用于油藏数值模拟。Heijn等[6]提出利用POD方法获得油水两相油藏的低阶数值模型。他们的结果说明POD产生的非线性模型在比较长的时期内都是有效的，并且指出对于同一油藏采用不同的工作制度进行多次模拟时，POD方法具有提高计算效率的潜力。Van Doren等[7]应用POD方法降低前模型和伴随模型的维数来加速油藏注水的优化过程，计算时间可减少35%。Cardoso等[8]提出利用快照聚类和丢失点估计技术进一步加速POD方法产生的降阶油藏模模型，他们取得了6到10倍的加速效果。He等[9,10]将POD方法与轨迹分段线性化方法相结合，分别应用于历史拟合和组分模拟的模型降阶，计算速度获得了大幅提高。

POD方法虽然近似非线性模型具有良好的准确性，但是它对用于产生数据样本的模拟参数(如边界条件、系统输入等)较敏感。当模拟POD降阶模型的输入与产生数据样本的输入不同时，结果的准确性就会较差[11-13]。本文提出将Krylov enhanced POD (简记为KPOD)方法应用于油藏数值模拟，可减少POD方法对输入参数的依赖，并通过实例证明了该方法的有效性。

1 油藏的数学模型

地面标准状况下黑油的数学模型为[14]

(1)

方程中的未知量为po,pg,pw,so,sg,sw。辅助方程为

(2)

式中：k为油藏多孔介质的绝对渗透率；φ为孔隙度；D为深度；g为重力加速度；so,sg,sw分别为油、气、水三相的饱和度；ρo,ρg,ρw分别为油、气、水三相的密度；Rso,Rsw分别为溶解气油比和溶解气水比；Bo,Bg,Bw分别为油、气、水三相的体积系数；kro,krg,krw分别为油、气、水三相的相对渗透率；μo,μg,μw分别为油、气、水三相的粘度；po,pg,pw分别为油、气、水三相的压力；qvo,qvg,qvw分别为地面标准条件下单位时间单位体积内产出或注入油、气、水的体积；pcow为油水两相之间的毛管压力；pcgo为气油两相之间的毛管压力。

上述的偏微分方程、辅助方程，再加上初始条件和边界条件就构成完整的数学模型。本文中采用隐压显饱(IMPES)求解方法，求解思路如下[14]：

(1) 将油水、油气之间的毛管压力方程代入方程组(1)中，消去pw和pg，得到只含po，so，sw和sg的方程组。

(2) 将油、气、水的渗流方程组(1)乘以适当的系数后合并，以消除方程组中的so，sw和sg项，得到只含有油相压力po的综合方程，称为压力方程。

(3) 写出压力方程的差分方程，其中的传导系数、毛管压力和产量项中与时间有关的非线性项均采用显示处理方法，得到线性代数方程组，记为

ApPo=bp

(3)

(4) 求解式(3)得到包含每个网格油相压力的向量Po，然后由毛管压力方程计算出每个网格的水相和气相压力。

(5) 由水组分方程显示计算出水相饱和度，由油组分方程显示计算出油相饱和度，然后由饱和度辅助方程计算出气相饱和度。

线性方程组(3)的阶数与油藏模型的网格数相同，对于实际油藏模型的网格数是几十万甚至几百万，所以，在IMPES求解过程中，线性方程组(3)的求解非常耗时。本文提出将应用KPOD模型降阶方法减少线性方程组的阶数，进而提高油藏模拟的运算速度。

2 POD模型降阶方法

本征正交分解(POD)模型降阶方法通过利用给定的理论或者实验得到的数据样本集来构造降阶过程中所需要的变换矩阵来达到降阶的目的。关于POD方法的基本原理可参考文献[5,15]，下面给出该方法在油藏数值模拟中的具体应用过程。

本文的油藏模拟器采用IMPES求解方法，为了构造POD基矩阵，首先在给定的工作制度下运行全阶油藏模拟器(也称为训练过程)，保存每个时间步的输出向量油相压力Po(也称为快照)，构成矩阵：

(4)

Xp=USWT

(5)

实行奇异值分解，式中U和W分别为矩阵Xp的左奇异向量和右奇异向量构成的矩阵，S为矩阵Xp的奇异值构成的对角阵。可得POD的基矩阵Φ为

Φ=U=[φ1,φ2,…,φm]

(6)

进一步可得POD的降阶基矩阵Φr为

Φr=[φ1,φ2, …,φr]

(7)

降阶基矩阵Φr包含基矩阵Φ的前r个列向量，其中个数r通过如下方式确定：

(8)

式中σi为矩阵Xp的奇异值，α的值是衡量POD降阶基矩阵Φr保留矩阵Xp特征多少的标准，一般取0.9<α<1.0。

在预测的过程中，改变工作制度，运行油藏模拟器，此时可将变换矩阵Φr应用于式(3)，令Po=Φrαo，可得

ApΦrαo=bp

(9)

(10)

求解降阶的线性方程组(10)，得到降阶解αo，然后再扩展到全局解Po=Φrαo。模型降阶的优势是式(10)只需求解r个变量，然而式(3)需要求解N(N是网格的数目) 个变量，其中r<

3 KPOD模型降阶方法

由于POD方法对训练过程的输入参数较敏感，如果模拟POD降阶模型的输入与训练过程的输入不同，结果的准确性就会变差[11-13]，因此，为提高POD方法对不同输入参数的反应能力， KPOD方法被提出。该方法结合Krylov 子空间模型降阶的Arnoldi方法的矩匹配性质和POD的数据泛化能力来减少POD方法对模型输入参数的依赖。

3.1 Krylov子空间降阶方法

Krylov子空间方法是模型降阶中最基本和最重要的方法，该方法通常采用所构造的标准列正交向量基底对系统进行模型降阶。Krylov子空间方法在数学理论上相当完善，其主要优点是算法稳定，计算量较小。Arnoldi和Lanczos算法是Krylov子空间模型降阶中两种常用的方法，其中Arnoldi算法已证明比Lanczos算法具有更高的计算效率和更精确的降阶模型[16,17]。

Krylov子空间方法的实现过程如下：对于线性方程组(3)，由矩阵Ap和向量bp生成的r维Krylov子空间为

(11)

为了构造子空间Kr(Ap;bp)的相应基底，应用Arnoldi算法对Kr(Ap;bp)进行标准的Gram-Schmidt正交化过程，最终可得由标准正交基底所构成的矩阵Vr=[v1v2…vr]。

利用降阶矩阵Vr就可对线性方程组(3)进行降阶。但是该方法仅限于线性时不变系统，即线性方程组的系数矩阵不随时间变化。

3.2 KPOD降阶方法

KPOD方法基于Arnoldi算法生成的基矩阵，再引入POD方法抽取全局基。由于Arnoldi算法的矩匹配性质，所以抽取的全局基能描述系统对不同输入参数的响应，缓解了POD方法对输入参数的依赖。该方法主要分为两步：(I)产生“局部” Krylov基；(II)POD提取全局基。

KPOD方法的具体步骤如下：

(I)产生“局部” Krylov基

(1)给定井底压力或流量，运行全阶油藏模拟器(称为“训练”过程)，记录每个时间步的矩阵Ap(ti)和向量bp(ti)，i=1,2,…,m。

(2)运行Arnoldi算法

for (i=1,2,…,m)

令s=1

while (norm>0)

令s=s+1

for (j=1,2,…,s)

endfor

end while

令V(ti)=[v1,v2,…,vs]

end for

构造矩阵集：令V=[V(t1)V(t2)…V(tm)]

(II)POD提取全局基

(1)生成全局基Φ

[U,S,W]=svd(V)

令Φ=U=[φ1φ2…φk]，其中(k=m×s)。

(2)生成截断全局基Φr

令r=1，sum=Sr

while (sum/sum(S)<ε)(ε是给定的百分比)

r=r+1，sum=sum+Sr

end while

令Φr=[φ1φ2…φr]。

将变换矩阵Φr应用于式(3)，具体降阶过程与上述POD过程相同。

4 数值实例

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

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

本文中， 笔者分别将POD方法和KPOD方法应用于该油藏，进行对比。关于POD方法在油藏模拟中的具体应用，可详见文献[8]。为了获得POD和KPOD的基矩阵，首先要运行全阶油藏模拟器(也被称为“训练”过程)。由于POD方法降阶基的预测能力直接受训练过程产生“快照”的工作制度影响，所以本文让两口生产井的井底压力在15MPa～18MPa之间任意变化，每100天(d)变化一次，工作制度如图3所示，两口注入井的井底压力不变，给定为31MPa。运行模拟器3000d，最大时间步长为10d，保存343个时间步的结果(也就是“快照”)。运用POD方法，运行结果构成的压力矩阵保留101个奇异值，获得POD基矩阵Φr∈N×r，其中r=101。这就意味着全阶模拟器要求解N=5049个未知变量，而降阶后只需求解101个变量。

为了说明KPOD比POD更有效，在训练过程中，本研究让生产井的井底压力始终保持不变，给定为18MPa，其他工作制度与POD方法的相同，保存343个时间步的矩阵Ap和向量bp。KPOD中的Arnoldi方法，取s=4，获得的基矩阵Φr∈N×r，其中r=88。式(5)的降阶模型仅仅需要求解88个未知量，它比POD方法要少。

图3 生产井井底压力的工作制度

下面采用两种不同的工作制度评估POD和KPOD的预测能力：

(1)工作制度1

在预测的过程中，给定两口生产井的井底压力为16MPa，它介于POD训练过程的15MPa～18MPa之间，注入井的井底压力与训练过程相同，最大时间步长为10d，预测2000d。

全阶油藏模拟器与应用POD方法、KPOD方法的降阶模拟器，两口生产井产油量、产水量，两口注水井的注水量计算结果对比如图4～图9。图中菱形实心点表示全阶油藏模拟器的运算结果，空心圆圈和菱形分别表示POD和KPOD降阶模拟器的运算结果。

本文中，笔者采用平均相对误差来衡量近似值的精度。如每口生产井产油量的平均相对误差定义为

(12)

两口生产井产油量、产水量，两口注水井注水量的平均相对误差如表1。

图4 生产井1的产油量对比

图5 生产井1的产水量对比

图6 生产井2的产油量对比

图7 生产井2的产水量对比

图8 注水井1的注水量对比

图9 注水井2的注水量对比

表1 四口井的平均相对误差

全阶油藏模拟器，POD和KPOD降阶油藏模拟器模拟时间的比较如表2。

表2 模拟时间的比较

通过以上图表显示，在预测过程中，当生产井的井底压力介于POD训练过程的井底压力范围之内时，POD方法的运算结果接近真值，降阶模拟器4口井的平均相对误差都能控制在5%的合理范围之内。然而，对于KPOD方法，即使生产井的井底压力不在KPOD方法的训练范围内，KPOD方法的运算结果比POD更接近真值，产生的平均相对误差更小。此外，由于KPOD基向量的个数少于POD基向量的个数，所以KPOD降阶模拟器的运算速度比POD要快，与全阶模拟器相比较，提速将近2倍。结果说明，对于工作制度一，KPOD方法在计算精度和速度上都优于POD方法。

(2)工作制度2

在工作制度2中，为了说明KPOD方法比POD方法可减少对训练过程输入参数的依赖，在计算精度和速度上更高效，给定两口生产井的井底压力为14MPa，它均超出POD和KPOD训练过程的井底压力范围，注入井的井底压力与训练过程相同，预测2000d。

全阶油藏模拟器与应用POD方法、KPOD方法的降阶模拟器计算结果对比如图10～图15。图中菱形的实心点表示全阶油藏模拟器的运行结果，空心的圆圈和菱形分别表示POD和KPOD降阶模拟器的运行结果。

图10 生产井1的产油量对比

图11 生产井1的产水量对比

图12 生产井2的产油量对比

图13 生产井2的产水量对比

图14 注水井1的注水量对比

图15 注水井2的注水量对比

两口生产井产油量、产水量，两口注水井注水量的平均相对误差如表3。

全阶油藏模拟器，POD和KPOD降阶油藏模拟器模拟时间的比较如表4。

由以上图表显示，当生产井的井底压力超出POD和KPOD训练过程的井底压力范围时，POD降阶模拟器的运算结果明显偏差较大，4口井的平均相对误差都超过5%，但KPOD降阶模拟器的运算结果仍接近真值，平均相对误差都控制在5%的合理范围内，说明KPOD方法的计算精度高于POD方法，减少了对训练过程输入参数的依赖。此时，KPOD降阶模拟器的运算速度与全阶模拟器相比较，提速仍将近2倍，说明在计算速度上也优于POD方法。

表3 四口井的平均相对误差

表4 模拟时间的比较

5 结 论

本文将KPOD方法应用于油藏数值模拟。该方法首先利用Krylov子空间的Arnoldi算法产生“局部” Krylov基，然后再利用POD方法提取全局基，获得降阶基矩阵。KPOD方法结合Arnoldi算法的矩匹配性质和POD的数据泛化能力可减少POD方法对训练过程模型输入参数的依赖。

将KPOD方法和POD方法应用于一个二维油水两相，包含5049个网格和4口井的油藏。本研究采用两种不同的工作制度评估POD和KPOD的预测能力。实例结果说明，在预测过程中，当生产井的井底压力介于POD方法训练范围之内时，POD降阶模拟器四口井的平均相对误差都控制在5%的合理范围内。当生产井的井底压力超出POD方法的训练过程时，四口井的平均相对误差都超过了5%，产生的偏差较大。然而，对于KPOD方法，即使生产井的井底压力不在训练范围之内，KPOD方法的平均相对误差仍能控制在5%的合理范围内。此外，KPOD方法基向量的个数较少，KPOD降阶模拟器的运算时间比POD要少，与全阶模拟器相比较，速度提高了近2倍。KPOD方法在计算速度和精度上都优于POD方法。研究表明，KPOD方法对于提高油藏模拟运算速度具有一定的应用价值。未来的工作是将该方法进一步应用于更大、更复杂的实际油藏。

[ 1] Heyouni M, Jbilou K. Matrix Krylov subspace methods for large scale model reduction problems.AppliedMathematicsandComputation, 2006, 181(2): 1215-1228

[ 2] Bai Z. Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems.AppliedNumericalMathematics, 2002, 43(1-2): 9-44

[ 3] Lall S, Marsden J, Glavas S. Structure-preserving model reduction for mechanical systems.PhysicaD-NonlinearPhenomena, 2003, 184(1-4): 304-318

[ 4] Awais M M, Shamail S, Ahmed N. Dimensionally reduced Krylov subspace model reduction for large scale systems.AppliedMathematicsandComputation, 2007, 191(1): 21-30

[ 5] 蒋耀林.模型降阶方法.北京：科学出版社，2010. 200-210

[ 6] Heijn T, Markovinovi R, Jansen J D. Generation of low-order reservoir models using system-theoretical concepts.SPEJournal, 2004, 9(2):202-218

[ 7] Van Doren J F M, Markovinovic R, Jansen J D. Reduced-order optimal control of water flooding using proper orthogonal decomposition.ComputationalGeosciences, 2006, 10(1):137-158

[ 8] Cardoso M A, Durlofsky L J, Sarma P. Development and application of reduced-order modeling procedures for subsurface flow simulation.InternationalJournalforNumericalMethodsinEngineering,2009, 77(9):1322-1350

[ 9] He J, Sarma P, Durlofsky L J. Reduced-order flow modeling and geological parameterization for ensemble-based data assimilation,ComputersandGeosciences, 2013, 55(55): 54-69

[10] He J, Durlofsky L J. Reduced-order modeling for compositional simulation by use of trajectory piecewise linearization.SPEJournal, 2014, 19 (5):858-872

[11] Hung E, Senturia S. Generating efficient dynamical models for microelec-tromechanical systems from a few finite element simulation runs.JournalofMicroelectromechanicalSystems, 1999, 8(3): 280-289

[12] Yvonnet J, Zahrouni H, Potier-Ferry M. A model reduction method for the post-buckling analysis of cellular microstructures.ComputerMethodsinAppliedMechanicsandEngineering, 2007, 197(1-4): 265-280

[13] Binion D, Chen X L. A Krylov enhanced proper orthogonal decomposition method for efficient nonlinear model reduction.FiniteElementsinAnalysisandDesign, 2011, 47(7): 728-738

[14] 李淑霞，谷建伟. 油藏数值模拟基础. 青岛市：中国石油大学出版社，2009.178-179

[15] Kerschen G, Golinval J C, Vakakis A F, et al. The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview.NonlinearDynamics, 2005,41(1-3): 147-169

[16] Rudnyi E B, Korvink J G. Review: automatic model reduction for transient simulation of MEMS-based devices.SensorsUpdate, 2002, 11(1): 3-33

[17] Binion D, Chen X L. Coupled electro-thermal-mechanical analysis for MEMS via model order reduction.FiniteElemAnalDesign, 2010, 46(12): 1068-1076

Research on a KPOD based model order reduction method for improving reservoir simulation speed

Cao Jing, Zhao Hui, Yu Gaoming

(College of Petroleum Engineering, Yangtze University, Wuhan 430100)

This study focused on the model order reduction technique to improve the calculation speed of a traditional reservoir simulator. Considering that when the technique uses proper orthogonal decomposition (POD) method, it can accelerate the simulator but sensitive to system input parameters, thus leading to the reduction of the simulation precision and efficiency, a Krylov enhanced POD (KPOD) based on Krylov subspace and POD, was applied to the technique, with the aim of combining the moment matching property of Arnoldi with POD’s data generalization ability to alleviate POD’s dependence on input conditions. The case study of real reservoir simulation demonstrated that the KPOD outperforms POD in the calculation speed and accuracy, and proved the soundness of the method.

reservoir simulation, model order reduction, proper orthogonal decomposition (POD), Krylov enhanced POD (KPOD)

①国家自然科学基金(51344003)资助项目。

2016-08-06)

②女，1979年生，博士生；研究方向：油藏数值模拟，最优化算法及控制系统方面的研究；联系人，E-mail: cjnew_2008@163.com