基于主成分分析和代理模型辅助的多目标生产优化方法

姚为英， 冯高城*， 任宜伟， 尹彦君， 王中正, 李振宇

(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模型进行了实例计算，优化得到的最优开发方案显著提高了油田开发的经济净现值。