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

曹 靜趙 輝 喻高明

(長江大學石油工程學院 武漢 430100)

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

曹 靜②趙 輝 喻高明

(長江大學石油工程學院 武漢 430100)

應用模型降階技術進行了提高傳統(tǒng)油藏模擬器運算速度的研究?？紤]到該技術采用本征正交分解(POD)方法能夠加速模擬器運算，但對系統(tǒng)的輸入?yún)?shù)較敏感，降低模擬的精度和效率，采用了基于Krylov子空間和POD的KPOD方法，以便利用Arnoldi的矩匹配性質(zhì)和POD的數(shù)據(jù)泛化能力減少POD方法對模型輸入?yún)?shù)的依賴。油藏模擬實例表明， KPOD方法在計算速度和精度上都優(yōu)于POD方法，驗證了該方法的有效性和實用性。

油藏數(shù)值模擬， 模型降階， 本征正交分解(POD)， KPOD方法

0 引 言

油藏數(shù)值模擬是現(xiàn)代油藏工程中不可缺少的研究工具。它普遍用于油藏開發(fā)設計、動態(tài)預測、油層參數(shù)識別、工程技術優(yōu)化設計以及重大開發(fā)技術政策的研究等。傳統(tǒng)油藏模擬器通過有限差分將偏微分方程組離散，轉化成非線性代數(shù)方程組，再采用迭代法求解。然而實際油藏，數(shù)值模擬地質(zhì)模型網(wǎng)格數(shù)可能高達幾十萬乃至幾百萬，采用傳統(tǒng)油藏模擬器進行模擬非常耗時，嚴重制約了數(shù)值模擬技術在我國大油藏特別是陸相非均質(zhì)嚴重油藏的應用，同時也造成實際油藏的生產(chǎn)歷史擬合以及后續(xù)的生產(chǎn)優(yōu)化過程困難重重，因此，如何大幅度提高油藏模擬運算速度是亟待解決的問題。

近年來，模型降階(model order reduction，MOR)技術在確保精度的條件下，對減少系統(tǒng)的運算量顯示出了一定的潛力[1-4]。該技術可將一個較大系統(tǒng)轉化為一個近似的較小系統(tǒng)，以便達到降低大型復雜系統(tǒng)的理論分析難度和減少數(shù)據(jù)運算量的目的[5]。本征正交分解(proper orthogonal decomposition，POD)法是非線性系統(tǒng)模型降階中使用最為廣泛的方法，該方法已被應用于很多領域。目前，它也被應用于油藏數(shù)值模擬。Heijn等[6]提出利用POD方法獲得油水兩相油藏的低階數(shù)值模型。他們的結果說明POD產(chǎn)生的非線性模型在比較長的時期內(nèi)都是有效的，并且指出對于同一油藏采用不同的工作制度進行多次模擬時，POD方法具有提高計算效率的潛力。Van Doren等[7]應用POD方法降低前模型和伴隨模型的維數(shù)來加速油藏注水的優(yōu)化過程，計算時間可減少35%。Cardoso等[8]提出利用快照聚類和丟失點估計技術進一步加速POD方法產(chǎn)生的降階油藏模模型，他們?nèi)〉昧?到10倍的加速效果。He等[9,10]將POD方法與軌跡分段線性化方法相結合，分別應用于歷史擬合和組分模擬的模型降階，計算速度獲得了大幅提高。

POD方法雖然近似非線性模型具有良好的準確性，但是它對用于產(chǎn)生數(shù)據(jù)樣本的模擬參數(shù)(如邊界條件、系統(tǒng)輸入等)較敏感。當模擬POD降階模型的輸入與產(chǎn)生數(shù)據(jù)樣本的輸入不同時，結果的準確性就會較差[11-13]。本文提出將Krylov enhanced POD (簡記為KPOD)方法應用于油藏數(shù)值模擬，可減少POD方法對輸入?yún)?shù)的依賴，并通過實例證明了該方法的有效性。

1 油藏的數(shù)學模型

地面標準狀況下黑油的數(shù)學模型為[14]

(1)

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

(2)

式中：k為油藏多孔介質(zhì)的絕對滲透率；φ為孔隙度；D為深度；g為重力加速度；so,sg,sw分別為油、氣、水三相的飽和度；ρo,ρg,ρw分別為油、氣、水三相的密度；Rso,Rsw分別為溶解氣油比和溶解氣水比；Bo,Bg,Bw分別為油、氣、水三相的體積系數(shù)；kro,krg,krw分別為油、氣、水三相的相對滲透率；μo,μg,μw分別為油、氣、水三相的粘度；po,pg,pw分別為油、氣、水三相的壓力；qvo,qvg,qvw分別為地面標準條件下單位時間單位體積內(nèi)產(chǎn)出或注入油、氣、水的體積；pcow為油水兩相之間的毛管壓力；pcgo為氣油兩相之間的毛管壓力。

上述的偏微分方程、輔助方程，再加上初始條件和邊界條件就構成完整的數(shù)學模型。本文中采用隱壓顯飽(IMPES)求解方法，求解思路如下[14]：

(1) 將油水、油氣之間的毛管壓力方程代入方程組(1)中，消去pw和pg，得到只含po，so，sw和sg的方程組。

(2) 將油、氣、水的滲流方程組(1)乘以適當?shù)南禂?shù)后合并，以消除方程組中的so，sw和sg項，得到只含有油相壓力po的綜合方程，稱為壓力方程。

(3) 寫出壓力方程的差分方程，其中的傳導系數(shù)、毛管壓力和產(chǎn)量項中與時間有關的非線性項均采用顯示處理方法，得到線性代數(shù)方程組，記為

ApPo=bp

(3)

(4) 求解式(3)得到包含每個網(wǎng)格油相壓力的向量Po，然后由毛管壓力方程計算出每個網(wǎng)格的水相和氣相壓力。

(5) 由水組分方程顯示計算出水相飽和度，由油組分方程顯示計算出油相飽和度，然后由飽和度輔助方程計算出氣相飽和度。

線性方程組(3)的階數(shù)與油藏模型的網(wǎng)格數(shù)相同，對于實際油藏模型的網(wǎng)格數(shù)是幾十萬甚至幾百萬，所以，在IMPES求解過程中，線性方程組(3)的求解非常耗時。本文提出將應用KPOD模型降階方法減少線性方程組的階數(shù)，進而提高油藏模擬的運算速度。

2 POD模型降階方法

本征正交分解(POD)模型降階方法通過利用給定的理論或者實驗得到的數(shù)據(jù)樣本集來構造降階過程中所需要的變換矩陣來達到降階的目的。關于POD方法的基本原理可參考文獻[5,15]，下面給出該方法在油藏數(shù)值模擬中的具體應用過程。

本文的油藏模擬器采用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個列向量，其中個數(shù)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。模型降階的優(yōu)勢是式(10)只需求解r個變量，然而式(3)需要求解N(N是網(wǎng)格的數(shù)目) 個變量，其中r<

3 KPOD模型降階方法

由于POD方法對訓練過程的輸入?yún)?shù)較敏感，如果模擬POD降階模型的輸入與訓練過程的輸入不同，結果的準確性就會變差[11-13]，因此，為提高POD方法對不同輸入?yún)?shù)的反應能力， KPOD方法被提出。該方法結合Krylov 子空間模型降階的Arnoldi方法的矩匹配性質(zhì)和POD的數(shù)據(jù)泛化能力來減少POD方法對模型輸入?yún)?shù)的依賴。

3.1 Krylov子空間降階方法

Krylov子空間方法是模型降階中最基本和最重要的方法，該方法通常采用所構造的標準列正交向量基底對系統(tǒng)進行模型降階。Krylov子空間方法在數(shù)學理論上相當完善，其主要優(yōu)點是算法穩(wěn)定，計算量較小。Arnoldi和Lanczos算法是Krylov子空間模型降階中兩種常用的方法，其中Arnoldi算法已證明比Lanczos算法具有更高的計算效率和更精確的降階模型[16,17]。

Krylov子空間方法的實現(xiàn)過程如下：對于線性方程組(3)，由矩陣Ap和向量bp生成的r維Krylov子空間為

(11)

為了構造子空間Kr(Ap;bp)的相應基底，應用Arnoldi算法對Kr(Ap;bp)進行標準的Gram-Schmidt正交化過程，最終可得由標準正交基底所構成的矩陣Vr=[v1v2…vr]。

利用降階矩陣Vr就可對線性方程組(3)進行降階。但是該方法僅限于線性時不變系統(tǒng)，即線性方程組的系數(shù)矩陣不隨時間變化。

3.2 KPOD降階方法

KPOD方法基于Arnoldi算法生成的基矩陣，再引入POD方法抽取全局基。由于Arnoldi算法的矩匹配性質(zhì)，所以抽取的全局基能描述系統(tǒng)對不同輸入?yún)?shù)的響應，緩解了POD方法對輸入?yún)?shù)的依賴。該方法主要分為兩步：(I)產(chǎn)生“局部” Krylov基；(II)POD提取全局基。

KPOD方法的具體步驟如下：

(I)產(chǎn)生“局部” 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 數(shù)值實例

圖1 油藏模型的滲透率分布

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

本文中， 筆者分別將POD方法和KPOD方法應用于該油藏，進行對比。關于POD方法在油藏模擬中的具體應用，可詳見文獻[8]。為了獲得POD和KPOD的基矩陣，首先要運行全階油藏模擬器(也被稱為“訓練”過程)。由于POD方法降階基的預測能力直接受訓練過程產(chǎn)生“快照”的工作制度影響，所以本文讓兩口生產(chǎn)井的井底壓力在15MPa～18MPa之間任意變化，每100天(d)變化一次，工作制度如圖3所示，兩口注入井的井底壓力不變，給定為31MPa。運行模擬器3000d，最大時間步長為10d，保存343個時間步的結果(也就是“快照”)。運用POD方法，運行結果構成的壓力矩陣保留101個奇異值，獲得POD基矩陣Φr∈N×r，其中r=101。這就意味著全階模擬器要求解N=5049個未知變量，而降階后只需求解101個變量。

為了說明KPOD比POD更有效，在訓練過程中，本研究讓生產(chǎn)井的井底壓力始終保持不變，給定為18MPa，其他工作制度與POD方法的相同，保存343個時間步的矩陣Ap和向量bp。KPOD中的Arnoldi方法，取s=4，獲得的基矩陣Φr∈N×r，其中r=88。式(5)的降階模型僅僅需要求解88個未知量，它比POD方法要少。

圖3 生產(chǎn)井井底壓力的工作制度

下面采用兩種不同的工作制度評估POD和KPOD的預測能力：

(1)工作制度1

在預測的過程中，給定兩口生產(chǎn)井的井底壓力為16MPa，它介于POD訓練過程的15MPa～18MPa之間，注入井的井底壓力與訓練過程相同，最大時間步長為10d，預測2000d。

全階油藏模擬器與應用POD方法、KPOD方法的降階模擬器，兩口生產(chǎn)井產(chǎn)油量、產(chǎn)水量，兩口注水井的注水量計算結果對比如圖4～圖9。圖中菱形實心點表示全階油藏模擬器的運算結果，空心圓圈和菱形分別表示POD和KPOD降階模擬器的運算結果。

本文中，筆者采用平均相對誤差來衡量近似值的精度。如每口生產(chǎn)井產(chǎn)油量的平均相對誤差定義為

(12)

兩口生產(chǎn)井產(chǎn)油量、產(chǎn)水量，兩口注水井注水量的平均相對誤差如表1。

圖4 生產(chǎn)井1的產(chǎn)油量對比

圖5 生產(chǎn)井1的產(chǎn)水量對比

圖6 生產(chǎn)井2的產(chǎn)油量對比

圖7 生產(chǎn)井2的產(chǎn)水量對比

圖8 注水井1的注水量對比

圖9 注水井2的注水量對比

表1 四口井的平均相對誤差

全階油藏模擬器，POD和KPOD降階油藏模擬器模擬時間的比較如表2。

表2 模擬時間的比較

通過以上圖表顯示，在預測過程中，當生產(chǎn)井的井底壓力介于POD訓練過程的井底壓力范圍之內(nèi)時，POD方法的運算結果接近真值，降階模擬器4口井的平均相對誤差都能控制在5%的合理范圍之內(nèi)。然而，對于KPOD方法，即使生產(chǎn)井的井底壓力不在KPOD方法的訓練范圍內(nèi)，KPOD方法的運算結果比POD更接近真值，產(chǎn)生的平均相對誤差更小。此外，由于KPOD基向量的個數(shù)少于POD基向量的個數(shù)，所以KPOD降階模擬器的運算速度比POD要快，與全階模擬器相比較，提速將近2倍。結果說明，對于工作制度一，KPOD方法在計算精度和速度上都優(yōu)于POD方法。

(2)工作制度2

在工作制度2中，為了說明KPOD方法比POD方法可減少對訓練過程輸入?yún)?shù)的依賴，在計算精度和速度上更高效，給定兩口生產(chǎn)井的井底壓力為14MPa，它均超出POD和KPOD訓練過程的井底壓力范圍，注入井的井底壓力與訓練過程相同，預測2000d。

全階油藏模擬器與應用POD方法、KPOD方法的降階模擬器計算結果對比如圖10～圖15。圖中菱形的實心點表示全階油藏模擬器的運行結果，空心的圓圈和菱形分別表示POD和KPOD降階模擬器的運行結果。

圖10 生產(chǎn)井1的產(chǎn)油量對比

圖11 生產(chǎn)井1的產(chǎn)水量對比

圖12 生產(chǎn)井2的產(chǎn)油量對比

圖13 生產(chǎn)井2的產(chǎn)水量對比

圖14 注水井1的注水量對比

圖15 注水井2的注水量對比

兩口生產(chǎn)井產(chǎn)油量、產(chǎn)水量，兩口注水井注水量的平均相對誤差如表3。

全階油藏模擬器，POD和KPOD降階油藏模擬器模擬時間的比較如表4。

由以上圖表顯示，當生產(chǎn)井的井底壓力超出POD和KPOD訓練過程的井底壓力范圍時，POD降階模擬器的運算結果明顯偏差較大，4口井的平均相對誤差都超過5%，但KPOD降階模擬器的運算結果仍接近真值，平均相對誤差都控制在5%的合理范圍內(nèi)，說明KPOD方法的計算精度高于POD方法，減少了對訓練過程輸入?yún)?shù)的依賴。此時，KPOD降階模擬器的運算速度與全階模擬器相比較，提速仍將近2倍，說明在計算速度上也優(yōu)于POD方法。

表3 四口井的平均相對誤差

表4 模擬時間的比較

5 結 論

本文將KPOD方法應用于油藏數(shù)值模擬。該方法首先利用Krylov子空間的Arnoldi算法產(chǎn)生“局部” Krylov基，然后再利用POD方法提取全局基，獲得降階基矩陣。KPOD方法結合Arnoldi算法的矩匹配性質(zhì)和POD的數(shù)據(jù)泛化能力可減少POD方法對訓練過程模型輸入?yún)?shù)的依賴。

將KPOD方法和POD方法應用于一個二維油水兩相，包含5049個網(wǎng)格和4口井的油藏。本研究采用兩種不同的工作制度評估POD和KPOD的預測能力。實例結果說明，在預測過程中，當生產(chǎn)井的井底壓力介于POD方法訓練范圍之內(nèi)時，POD降階模擬器四口井的平均相對誤差都控制在5%的合理范圍內(nèi)。當生產(chǎn)井的井底壓力超出POD方法的訓練過程時，四口井的平均相對誤差都超過了5%，產(chǎn)生的偏差較大。然而，對于KPOD方法，即使生產(chǎn)井的井底壓力不在訓練范圍之內(nèi)，KPOD方法的平均相對誤差仍能控制在5%的合理范圍內(nèi)。此外，KPOD方法基向量的個數(shù)較少，KPOD降階模擬器的運算時間比POD要少，與全階模擬器相比較，速度提高了近2倍。KPOD方法在計算速度和精度上都優(yōu)于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] 李淑霞，谷建偉. 油藏數(shù)值模擬基礎. 青島市：中國石油大學出版社，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年生，博士生；研究方向：油藏數(shù)值模擬，最優(yōu)化算法及控制系統(tǒng)方面的研究；聯(lián)系人，E-mail: cjnew_2008@163.com