基于主成分分析和代理模型輔助的多目標生產優化方法

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

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

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

1 油藏魯棒生產優化問題

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

(1)

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

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

因此，在考慮實際油藏生產中地質不確定性因素的基礎上，建立以NPV最大化的水驅油藏魯棒生產模型為

(2)

式(2)中：目標函數Je(x)為平均經濟凈現值；Ne為基于概率生成的能夠反映油藏特征的模型個數。

需要注意的是，僅僅優化性能指標的數學期望并不能降低不確定性，反而會導致一個很大的方差區間。針對上述問題，魯棒生產優化具有更好的穩定性，且優化得到的開發方案更符合油田現場實際開發需求。換句話說，需要降低開發方案對地質不確定性參數的敏感性，從而提高優化結果的可靠性。

從平均經濟凈現值和最小經濟凈現值角度出發，建立降低不確定性油藏開發風險的多目標優化模型。假設井控只考慮最簡單的邊界約束，則該魯棒生產優化問題的數學模型為

(3)

2 算法機制

代理模型方法通過構建計算廉價的近似模型代替原始數值模型，可以大幅度降低原始模型的評估代價，因此能夠有效提高優化效率。基于代理模型的方法通常建立在具有強大全局搜索能力的進化算法基礎上。采用目前流行的基于分解的多目標進化算法(multi-objective evolutionary algorithm based on decomposition, MOEA/D) 作為優化求解器[12]。

2.1 基于分解的多目標進化算法

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

(4)

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

在此給出多目標優化中Pareto最優的相關概念。給定兩個可行解x1和x2，當且僅當滿足

(5)

時，稱x1支配x2，顯然，此時解x1優于x2，并記作x1

給定解x*∈X，當不存在x滿足x

基于分解的多目標優化算法首先需要在目標空間對多目標優化問題進行分解，在本文中，使用了切比雪夫的分解方法，示意圖如圖1所示。

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

假設一組均勻分布的權重向量{λ1,λ2,…,λN}，逼近真實Pareto前沿的優化問題可以通過切比雪夫分解法分解為N個標量化的子問題，其中每個子問題的表達式為

(6)

2.2 代理模型技術

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

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

g(x)=μ+ε(x)

(7)

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

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

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

(8)

(9)

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

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

(10)

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

(11)

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

(12)

(13)

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

(14)

其均方誤差為

(15)

(16)

2.3 數據預處理的降維技術

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

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

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

(17)

將原始的樣本集X=(x1,x2,…,xd)轉變為新的l維數據集，只需要將協方差矩陣XXT進行特征值分解，并排序取前l個值對應的特征向量構成W=(w1,w2,…,wl)，這就是PCA的解。

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

2.4 K-MOEA/D-PCA算法優化計算步驟

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

本文使用拉丁超立方試驗設計方法[18]從可行域中采樣生成初始訓練樣本點。然后，利用真實的油藏數值模擬器對目標函數值進行計算。這些經過真實評估的解決方案存儲在樣本庫D中。以最大的函數評估數作為進化優化過程的終止條件。

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

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

3 計算實例

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

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

圖4顯示了在優化過程中平均和最小經濟凈現值的變化情況。可以看出，未進行優化時，初始解構成的種群平均和最小經濟凈現值均較低，通過提出的K-MOEA/D-PCA算法進行優化后，能夠找到位于高經濟凈現值區域的候選解，最終得到了帶來最高回報的最優解，顯著提高了優化方案的經濟凈現值。

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

使用圖5所示的優化方案對相應的模型進行模擬開發，生產過程中累積產油量、累積產水量、累積注水量和含水率的變化情況如圖6所示。可以看出，考慮到模型的不確定性，生產3 600 d后，該模型的累積產油量可以達到(4 ～ 4.5)×105m3，含水率約為0.9。

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

為了保證解的穩定性，避免偶然性，為不同模型執行10次獨立計算，最終的優化結果如圖7所示，箱型圖中的紅線代表10次運行結果的中位數。可以看出，由于不同模型的滲透率場不同，開發過程實現的經濟凈現值也不同，經過本文提出的算法優化后的生產方案可以實現滿意的經濟效益。

圖4 優化過程經濟凈現值的變化Fig.4 The change of NPV during optimization process

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

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

圖7 不同模型的開發方案的經濟凈現值Fig.7 Net present value of development schemes for different models

4 結論

針對油藏魯棒生產優化問題的大規模與計算耗時的特性，提出了新的求解思路。由于復雜油藏注采井數多、非線性強，直接將代理模型引入進化算法會導致優化結果變差，因此提出了K-MOEA/D-PCA方法，通過主成分分析對數據進行降維預處理，然后采用先分后治的求解策略，利用權重向量將目標空間分解為多個子問題從而協同計算，并通過構建代理模型動態規劃搜索區域，實現了高效全局優化。應用K-MOEA/D-PCA方法對Egg模型進行了實例計算，優化得到的最優開發方案顯著提高了油田開發的經濟凈現值。