姚為英， 馮高城*， 任宜偉， 尹彥君， 王中正, 李振宇

(1.中海油能源發(fā)展股份有限公司工程技術(shù)分公司， 天津 300452; 2.中國石油大學(xué)(華東)石油工程學(xué)院， 青島 266580)

基于油藏?cái)?shù)值模型的生產(chǎn)優(yōu)化是油藏開發(fā)管理的關(guān)鍵，即尋找最優(yōu)井控方案，改進(jìn)油藏開發(fā)管理策略。通過優(yōu)化每口井的控制方案，獲得最大的油藏開發(fā)經(jīng)濟(jì)效益，一般用經(jīng)濟(jì)凈現(xiàn)值 (net present value, NPV)來衡量[1-3]。傳統(tǒng)的油藏生產(chǎn)優(yōu)化算法大都使用單個精確的油藏?cái)?shù)值模型，結(jié)合優(yōu)化算法調(diào)整油藏區(qū)塊內(nèi)的注采井的工作參數(shù)，使其維持在最優(yōu)井控的工作狀態(tài)，從而最大化生產(chǎn)效益[4]。然而，此類方法忽略了油藏地質(zhì)的不確定性，由于一些未知的儲層參數(shù)，如滲透率和孔隙度，可能會顯著影響儲層性質(zhì)，因此油藏地質(zhì)具有很高的不確定性[5]。如果不考慮油藏地質(zhì)的不確定性，采用單一模型優(yōu)化得到的開發(fā)方案在現(xiàn)場實(shí)施時可能會面臨不同程度的開發(fā)風(fēng)險。因此，開發(fā)方案在實(shí)施前必須考慮油藏地質(zhì)的不確定性并進(jìn)行多目標(biāo)風(fēng)險分析。Brouwer等[6]首次提出通過采用歷史擬合得到的多個概率模型來表征地質(zhì)參數(shù)的不確定性，然后綜合考慮多個油藏模型進(jìn)行生產(chǎn)優(yōu)化，這類生產(chǎn)優(yōu)化問題被稱為魯棒生產(chǎn)優(yōu)化問題。此外，多目標(biāo)優(yōu)化方法[3]也常用于考慮地質(zhì)不確定性時進(jìn)行油藏魯棒生產(chǎn)優(yōu)化，該方法能夠同時優(yōu)化歷史擬合得到的多個概率模型的平均經(jīng)濟(jì)凈現(xiàn)值和經(jīng)濟(jì)凈現(xiàn)值的離散程度，從而得到一組權(quán)衡多個目標(biāo)的最優(yōu)生產(chǎn)方案。然而，直接將多目標(biāo)優(yōu)化算法應(yīng)用于魯棒生產(chǎn)優(yōu)化問題時，需要面臨決策變量規(guī)模大和目標(biāo)函數(shù)評估計(jì)算耗時長等問題[5]。對此，國內(nèi)外學(xué)者近年來提出了多種基于代理模型的多目標(biāo)優(yōu)化算法：Chang等[7]將多項(xiàng)式擬合引入到考慮地質(zhì)不確定性的多目標(biāo)井位魯棒優(yōu)化中，與傳統(tǒng)方法相比有效降低了計(jì)算成本；馬寧[8]提出了一種基于遺傳算法的多目標(biāo)優(yōu)化建設(shè)項(xiàng)目調(diào)度模型，大大優(yōu)化了建設(shè)項(xiàng)目的調(diào)度過程；王鵬宇等[9]使用分布式訓(xùn)練的多目標(biāo)算法，顯著提高了多目標(biāo)算法優(yōu)化的速度。然而，現(xiàn)有代理模型應(yīng)用于高維參數(shù)空間時性能會顯著下降，導(dǎo)致處理油藏多目標(biāo)生產(chǎn)優(yōu)化問題時效果欠佳[10]。為解決上述問題，現(xiàn)提出結(jié)合主成分分析降維算法和代理模型輔助的多目標(biāo)生產(chǎn)優(yōu)化新方法，將該方法應(yīng)用到標(biāo)準(zhǔn)測試模型，并將計(jì)算結(jié)果與現(xiàn)有進(jìn)化算法進(jìn)行對比分析。

1 油藏魯棒生產(chǎn)優(yōu)化問題

隨著油氣田開采程度的不斷加深，開采方式的經(jīng)濟(jì)效益越來越被油田生產(chǎn)管理者所重視，因此通常需要選擇最優(yōu)的工作措施和開發(fā)方案。例如，在水驅(qū)油田開發(fā)過程中，需要通過提高原油產(chǎn)量的收益去補(bǔ)償注水措施帶來的額外成本。油藏生產(chǎn)優(yōu)化的目的是最大化油藏的開發(fā)潛力，從而獲得最大的油田開發(fā)利潤。因此，采用經(jīng)濟(jì)凈現(xiàn)值作為油藏生產(chǎn)優(yōu)化問題的性能指標(biāo)，表示在整個油田生產(chǎn)周期里，項(xiàng)目預(yù)期實(shí)現(xiàn)的經(jīng)濟(jì)效益與增產(chǎn)措施所產(chǎn)生的投資成本之間的差額。此時，水驅(qū)注采優(yōu)化模型的目標(biāo)函數(shù)可被定義為

(1)

式(1)中：J(x,s)為經(jīng)濟(jì)凈現(xiàn)值，$；x為控制變量；s為油藏狀態(tài)向量；Nt為總控制時間步數(shù)；b為年折現(xiàn)率，%；pt為第n個控制步的時間尺度參數(shù)，年；Δtn為第n個調(diào)控步的時長，d；ro表示油價，$/桶(1桶=158.98 L)；Qo,n為第n個控制步的平均產(chǎn)油速度，1桶/d；rw為產(chǎn)出水的處理費(fèi)用，$/桶；Qw,n為第n個控制步的平均產(chǎn)水速度，1桶/d；rwi表示注入水費(fèi)用，$/桶；Qwi,n為第n個控制步的平均注水速度，1桶/d。

魯棒生產(chǎn)優(yōu)化問題本質(zhì)上是一個多目標(biāo)優(yōu)化問題，即整個油藏系統(tǒng)的平均性能指標(biāo)最大化，同時相關(guān)的開發(fā)風(fēng)險最小化。其中，相關(guān)的開發(fā)風(fēng)險通常采用目標(biāo)函數(shù)在所有可能模型上的方差來衡量。最后，它可以生成具有一組最優(yōu)解決方案的Pareto前沿[11]，供決策者在不同的風(fēng)險水平下分析實(shí)施各解決方案的優(yōu)劣情況。

因此，在考慮實(shí)際油藏生產(chǎn)中地質(zhì)不確定性因素的基礎(chǔ)上，建立以NPV最大化的水驅(qū)油藏魯棒生產(chǎn)模型為

(2)

式(2)中：目標(biāo)函數(shù)Je(x)為平均經(jīng)濟(jì)凈現(xiàn)值；Ne為基于概率生成的能夠反映油藏特征的模型個數(shù)。

需要注意的是，僅僅優(yōu)化性能指標(biāo)的數(shù)學(xué)期望并不能降低不確定性，反而會導(dǎo)致一個很大的方差區(qū)間。針對上述問題，魯棒生產(chǎn)優(yōu)化具有更好的穩(wěn)定性，且優(yōu)化得到的開發(fā)方案更符合油田現(xiàn)場實(shí)際開發(fā)需求。換句話說，需要降低開發(fā)方案對地質(zhì)不確定性參數(shù)的敏感性，從而提高優(yōu)化結(jié)果的可靠性。

從平均經(jīng)濟(jì)凈現(xiàn)值和最小經(jīng)濟(jì)凈現(xiàn)值角度出發(fā)，建立降低不確定性油藏開發(fā)風(fēng)險的多目標(biāo)優(yōu)化模型。假設(shè)井控只考慮最簡單的邊界約束，則該魯棒生產(chǎn)優(yōu)化問題的數(shù)學(xué)模型為

(3)

2 算法機(jī)制

代理模型方法通過構(gòu)建計(jì)算廉價的近似模型代替原始數(shù)值模型，可以大幅度降低原始模型的評估代價，因此能夠有效提高優(yōu)化效率?；诖砟Ｐ偷姆椒ㄍǔ＝⒃诰哂袕?qiáng)大全局搜索能力的進(jìn)化算法基礎(chǔ)上。采用目前流行的基于分解的多目標(biāo)進(jìn)化算法(multi-objective evolutionary algorithm based on decomposition, MOEA/D) 作為優(yōu)化求解器[12]。

2.1 基于分解的多目標(biāo)進(jìn)化算法

基于分解的多目標(biāo)進(jìn)化算法(MOEA/D)是一種流行的啟發(fā)式算法，用于解決全局優(yōu)化問題，并廣泛應(yīng)用于工程優(yōu)化問題。經(jīng)典的MOEA/D算法如下：一般來說，多目標(biāo)優(yōu)化問題 (multi-objective optimization problems, MOPs) 由多個互相沖突的目標(biāo)函數(shù)組成，一個具有m個目標(biāo)函數(shù)的多目標(biāo)優(yōu)化問題的數(shù)學(xué)模型表達(dá)式為

(4)

式(4)中：x=(x1,x2,…,xd)∈X?Rd為優(yōu)化變量，X為已知的決策空間，d為設(shè)計(jì)空間的維數(shù)；m為目標(biāo)函數(shù)的個數(shù)；F(x):Rd→Rm為目標(biāo)向量，所有F(x)構(gòu)成了未知的目標(biāo)空間；gi為p個不等式約束；hj為q個等式約束。

在此給出多目標(biāo)優(yōu)化中Pareto最優(yōu)的相關(guān)概念。給定兩個可行解x1和x2，當(dāng)且僅當(dāng)滿足

(5)

時，稱x1支配x2，顯然，此時解x1優(yōu)于x2，并記作x1

給定解x*∈X，當(dāng)不存在x滿足x

基于分解的多目標(biāo)優(yōu)化算法首先需要在目標(biāo)空間對多目標(biāo)優(yōu)化問題進(jìn)行分解，在本文中，使用了切比雪夫的分解方法，示意圖如圖1所示。

圖1 切比雪夫分解法示意圖Fig.1 Schematic diagram of Tchebycheff decomposition

假設(shè)一組均勻分布的權(quán)重向量{λ1,λ2,…,λN}，逼近真實(shí)Pareto前沿的優(yōu)化問題可以通過切比雪夫分解法分解為N個標(biāo)量化的子問題，其中每個子問題的表達(dá)式為

(6)

2.2 代理模型技術(shù)

代理模型是為了降低計(jì)算復(fù)雜度而建立的真實(shí)模型的簡化近似[13]。因此，代理模型可以使用較低的評估成本預(yù)測新解的適應(yīng)度值。如果真實(shí)的模擬值表示為f(x)，代理模型的預(yù)測值表示為f′(x)，則f′(x)=f′(x) +e(x)，其中e(x)為近似誤差。f(x)的內(nèi)部計(jì)算過程是不需要被完全了解的，只有輸入與輸出的映射關(guān)系很重要，這類統(tǒng)稱為黑盒建模。其中，廣泛應(yīng)用的代理技術(shù)有多項(xiàng)式回歸、克里金、徑向基函數(shù)以及支持向量回歸等等，主要通過數(shù)據(jù)插值或回歸的方法擬合決策變量與適應(yīng)度函數(shù)之間的映射關(guān)系。本文中采用克里金插值模型[14]來逼近目標(biāo)函數(shù)。目前為止，克里金是應(yīng)用最廣泛的代理技術(shù)，主要是因?yàn)樗軌蛱峁╊A(yù)測點(diǎn)的不確定性信息，這在代理的模型管理策略中非常有價值。

為了構(gòu)建計(jì)算廉價的代理模型y=g(x)，對于任意的個體x，克里金模型假設(shè)其目標(biāo)函數(shù)值可被估計(jì)為

g(x)=μ+ε(x)

(7)

式(7)中：ε(x)～N(0,σ2)；μ和σ為兩個與x無關(guān)的常數(shù)。換句話說，g(x)的先驗(yàn)分布為高斯分布，且均值為μ，方差為σ2。

此外，對于任意x，x′∈Xd，ε(x)與ε(x′)之間的關(guān)聯(lián)式c(x,x′)，取決于|x-x′|。更準(zhǔn)確地表達(dá)為

c(x,x′)=exp[-d(x,x′)]

(8)

(9)

式(9)中：θi>0 且 1≤pi≤2；當(dāng)d(x,x′)→0時，c(x,x′)→1，當(dāng)d(x,x′)→+∞時，c(x,x′)→0，這種特征在為連續(xù)函數(shù)f(x)建模時是符合要求的；pi與g(x)在xi處的平滑度有關(guān)；θi表示xi對g(x)的重要性。更多關(guān)于克里金建模的細(xì)節(jié)見文獻(xiàn)[14]。

給定K個點(diǎn)x1,x2,…,xK∈Xd，與其對應(yīng)的函數(shù)值y1,y2,…,yK，超參數(shù)μ,θ1,θ2,…,θd和p1,p2,…,pd可以通過g(x)=yi在x=xi(i=1,2,…,K) 處的最大似然估計(jì)得到，即

(10)

式(10)中：y=(y1,y2,…,yK)T；I為K維列向量；C為K×K的矩陣，表達(dá)式為

(11)

為了最大化式(10)，μ和σ2的值分別取為

(12)

(13)

將式(12)和式(13)代入式(10)中，消除了式(10)中的未知參數(shù)μ和σ2。因此，似然函數(shù)僅依賴于θi和pi。然后，可以使用優(yōu)化方法最大化式(10)來估計(jì)θi和pi，再通過式(12)和式(13)可以很容易地估算出μ和σ2。

(14)

其均方誤差為

(15)

(16)

2.3 數(shù)據(jù)預(yù)處理的降維技術(shù)

盡管克里金是一個很流行的能夠提供不確定性信息的代理模型，但由于訓(xùn)練模型的計(jì)算復(fù)雜度，它也存在潛在的嚴(yán)重缺陷。研究發(fā)現(xiàn)，大規(guī)模優(yōu)化問題的克里金模型訓(xùn)練時間是無法忽略的，因此，訓(xùn)練時間可能會比真實(shí)評估目標(biāo)函數(shù)的時間長，從而破壞了降低計(jì)算成本的整體目標(biāo)。

由于油藏生產(chǎn)優(yōu)化方案設(shè)計(jì)的要求不斷提高，從決策變量空間來講，設(shè)計(jì)的變量愈加復(fù)雜，由此帶來變量數(shù)目的劇增。“維數(shù)災(zāi)難”問題的限制導(dǎo)致直接使用克里金插值構(gòu)建代理模型近似誤差大，且近似模型訓(xùn)練參數(shù)時間長，失去了代理輔助的意義，因此本文中采用主成分分析方法[17](principal component analysis, PCA) 對優(yōu)化問題進(jìn)行降維處理。

PCA的目標(biāo)就是在盡力保證信息量不損失的前提下，對原始的數(shù)據(jù)特征進(jìn)行降維，從而在較小的維度下展示數(shù)據(jù)的特征。PCA的原理是利用正交變換將潛在線性相關(guān)的樣本轉(zhuǎn)變?yōu)橐唤M線性不相關(guān)的低維變量，轉(zhuǎn)變完之后的這組變量稱作主成分。假設(shè)原始樣本點(diǎn)xi在新的低維空間上的投影為WTxi，根據(jù)最大可分性，所有的樣本點(diǎn)的投影需盡量分開，即最大化樣本點(diǎn)投影后的方差，因此優(yōu)化的目標(biāo)函數(shù)為

(17)

將原始的樣本集X=(x1,x2,…,xd)轉(zhuǎn)變?yōu)樾碌膌維數(shù)據(jù)集，只需要將協(xié)方差矩陣XXT進(jìn)行特征值分解，并排序取前l(fā)個值對應(yīng)的特征向量構(gòu)成W=(w1,w2,…,wl)，這就是PCA的解。

通過PCA對高維的設(shè)計(jì)空間進(jìn)行降維處理，把較高維數(shù)的油藏設(shè)計(jì)變量在低維空間表征，從而提高代理模型的近似精度以及優(yōu)化效率。

2.4 K-MOEA/D-PCA算法優(yōu)化計(jì)算步驟

K-MOEA/D-PCA算法包含兩個循環(huán)，如圖2所示，主循環(huán)和次循環(huán)分別是使用油藏?cái)?shù)值模擬器和代理模型的優(yōu)化過程。該算法有3個步驟：①初始化；②克里金代理模型輔助MOEA/D進(jìn)化搜索；③更新代理模型。

本文使用拉丁超立方試驗(yàn)設(shè)計(jì)方法[18]從可行域中采樣生成初始訓(xùn)練樣本點(diǎn)。然后，利用真實(shí)的油藏?cái)?shù)值模擬器對目標(biāo)函數(shù)值進(jìn)行計(jì)算。這些經(jīng)過真實(shí)評估的解決方案存儲在樣本庫D中。以最大的函數(shù)評估數(shù)作為進(jìn)化優(yōu)化過程的終止條件。

隨后，為每個目標(biāo)函數(shù)訓(xùn)練克里金代理模型，以減少計(jì)算成本。由于克里金模型能夠在誤差范圍內(nèi)進(jìn)行預(yù)測，因此得到了廣泛應(yīng)用。另一方面，處理高維優(yōu)化問題時，構(gòu)造克里金模型的計(jì)算時間過長又導(dǎo)致該方法使用受限。因此，在構(gòu)建代理模型之前，采用主成分分析克里金將樣本數(shù)據(jù)從高維映射到低維。因此，低維數(shù)據(jù)可以提高克里金模型的精度。代理模型的管理策略對代理輔助進(jìn)化計(jì)算的性能是至關(guān)重要的，即選擇個體使用數(shù)值模擬器重新評估的策略，然后使用這些重新評估的解決方案來更新代理模型，并將其添加到樣本庫D中。因此，需要綜合考慮這些解決方案的收斂性和多樣性。

圖2 K-MOEA/D-PCA算法流程Fig.2 Flow diagram of K-MOEA/D-PCA

3 計(jì)算實(shí)例

利用三維油藏模型測試K-MOEA/D-PCA算法在處理考慮地質(zhì)不確定性條件下魯棒生產(chǎn)優(yōu)化問題時的有效性和收斂速度。采用基準(zhǔn)測試模型(Egg模型)的集合版本(包含100個不確定模型)，該版本是公開的，詳細(xì)地質(zhì)統(tǒng)計(jì)學(xué)參數(shù)可見文獻(xiàn)[19]。由于模型之間具有相似性，基于100個數(shù)值模型的評估耗時長且浪費(fèi)計(jì)算資源，因此在魯棒生產(chǎn)優(yōu)化過程前，采用多維尺度變換和K-means聚類算法來減少計(jì)算量，但不影響結(jié)果的準(zhǔn)確性。利用聚類算法，將隨機(jī)模型按滲透率場的相似度分為6類，并從每類中選擇一個具有代表性的模型。如圖3所示，每個模型共有12口直井，其中8口注水井，4口生產(chǎn)井。每個模型包含60 × 60 × 7個網(wǎng)格，其中18 553個為活網(wǎng)格。

生產(chǎn)井井底壓力恒為395 bar(1 bar=0.1 MPa)，因此控制輸入變量為每口注水井的流量，其上下限分別設(shè)置為0、100 桶/d。經(jīng)濟(jì)參數(shù)中的油價、注水成本、生產(chǎn)成本分別設(shè)置為45、5、5 $/桶，假設(shè)折現(xiàn)系數(shù)為0。每個油藏模型的生產(chǎn)周期為3 600 d，時間步長為180 d，因此，優(yōu)化決策變量的個數(shù)為8 × 20 = 160。注意，每個解決方案的目標(biāo)函數(shù)評估是基于六個油藏?cái)?shù)值模型的。

圖4顯示了在優(yōu)化過程中平均和最小經(jīng)濟(jì)凈現(xiàn)值的變化情況?？梢钥闯?，未進(jìn)行優(yōu)化時，初始解構(gòu)成的種群平均和最小經(jīng)濟(jì)凈現(xiàn)值均較低，通過提出的K-MOEA/D-PCA算法進(jìn)行優(yōu)化后，能夠找到位于高經(jīng)濟(jì)凈現(xiàn)值區(qū)域的候選解，最終得到了帶來最高回報的最優(yōu)解，顯著提高了優(yōu)化方案的經(jīng)濟(jì)凈現(xiàn)值。

圖5顯示了經(jīng)過本文提出的K-MOEA/D-PCA算法優(yōu)化后，6個不同模型的最優(yōu)注水開發(fā)方案。在每個圖中，橫坐標(biāo)為調(diào)控的時間步，縱坐標(biāo)為不同注水井的注水量。所有模型不同時間步的注水量被限制在0～100 桶/d。

使用圖5所示的優(yōu)化方案對相應(yīng)的模型進(jìn)行模擬開發(fā)，生產(chǎn)過程中累積產(chǎn)油量、累積產(chǎn)水量、累積注水量和含水率的變化情況如圖6所示?？梢钥闯?，考慮到模型的不確定性，生產(chǎn)3 600 d后，該模型的累積產(chǎn)油量可以達(dá)到(4 ～ 4.5)×105m3，含水率約為0.9。

圖3 Egg模型對數(shù)滲透率分布Fig.3 Log-permeability distribution of Egg models

為了保證解的穩(wěn)定性，避免偶然性，為不同模型執(zhí)行10次獨(dú)立計(jì)算，最終的優(yōu)化結(jié)果如圖7所示，箱型圖中的紅線代表10次運(yùn)行結(jié)果的中位數(shù)?？梢钥闯?，由于不同模型的滲透率場不同，開發(fā)過程實(shí)現(xiàn)的經(jīng)濟(jì)凈現(xiàn)值也不同，經(jīng)過本文提出的算法優(yōu)化后的生產(chǎn)方案可以實(shí)現(xiàn)滿意的經(jīng)濟(jì)效益。

圖4 優(yōu)化過程經(jīng)濟(jì)凈現(xiàn)值的變化Fig.4 The change of NPV during optimization process

圖5 不同概率模型中各注水井的最優(yōu)注水方案Fig.5 The optimal water-injection scheme for each injection well within different models

圖6 不同模型開發(fā)方案下各指標(biāo)的變化Fig.6 The change of each indicator for the model with different development schemes

圖7 不同模型的開發(fā)方案的經(jīng)濟(jì)凈現(xiàn)值Fig.7 Net present value of development schemes for different models

4 結(jié)論

針對油藏魯棒生產(chǎn)優(yōu)化問題的大規(guī)模與計(jì)算耗時的特性，提出了新的求解思路。由于復(fù)雜油藏注采井?dāng)?shù)多、非線性強(qiáng)，直接將代理模型引入進(jìn)化算法會導(dǎo)致優(yōu)化結(jié)果變差，因此提出了K-MOEA/D-PCA方法，通過主成分分析對數(shù)據(jù)進(jìn)行降維預(yù)處理，然后采用先分后治的求解策略，利用權(quán)重向量將目標(biāo)空間分解為多個子問題從而協(xié)同計(jì)算，并通過構(gòu)建代理模型動態(tài)規(guī)劃搜索區(qū)域，實(shí)現(xiàn)了高效全局優(yōu)化。應(yīng)用K-MOEA/D-PCA方法對Egg模型進(jìn)行了實(shí)例計(jì)算，優(yōu)化得到的最優(yōu)開發(fā)方案顯著提高了油田開發(fā)的經(jīng)濟(jì)凈現(xiàn)值。