CN104615862A - 基于进化算法的高含水油田确定井位的方法 - Google Patents

Abstract

本发明提供一种基于进化算法的高含水油田确定井位的方法，包括步骤：建立预定油区的油藏模型，将油藏模型划分网格；获取油藏资料，根据油藏资料按照历史拟合得出剩余油潜力富集区；计算剩余油潜力富集区中表示剩余油分布和注入水波及体积的井位约束值，得出井位优势区域；将优势区域内的油井运用进化算法进行迭代计算适应度函数，利用克里金算法对适应度函数进行评估，进化算法根据适应度函数进行调整，克里金算法考虑多维空间、块金效应和局部权重；根据进化算法迭代计算，得到最大适应度函数值，从而确定井位。本发明以改进的进化算法和克里金算法为优化工具，构建动态潜力指数缩小进化算法的搜索空间，使优化计算收敛精度和速度大大增加。

Description

基于进化算法的高含水油田确定井位的方法

技术领域

本发明属于高含水油田开发调整井位技术领域，尤其涉及一种基于进化算法的高含水油田确定井位的方法。

背景技术

对于我国注水开发的油藏而言，大部分是已经注水开发了几十年的高含水老油田，地下储层非均质以及油水关系的复杂性，导致剩余油分布呈现“总体零散、局部集中”的复杂特征，对油田后期调整造成极大的困难，直接影响了注水开发效果和经济效益。确定合理调整井位是高含水老油田调整的重要主体技术，是提高储量动用程度和采收率的重要手段。

合理调整井位受油藏地质条件、流体性质等不确定因素及因素间复杂非线性关系的影响，常用的数值模拟法需要大量的模拟方案，而且也难以获得最优的调整井位。智能优化算法是确定最优调整井位的有效手段。

目前常用的调整井位的优化方法主要分为三类：随机优化法、梯度优化法、质量图法。通过调研分析发现，针对高含水老油田该三类方法均存在一定的不适应性：

(1)随机优化法。该类算法从多组初始解出发搜索，易获得全局最优解，但高含水期地下油水关系的复杂性及存在的诸多地质不确定性，使得该类算法搜索空间较大，计算时间长，加之优化参数选择的人为性，该类算法的性能存在较大提升空间。

(2)梯度优化法。该类算法通过计算目标函数的梯度，使井位向最大梯度方向移动，从而实现井位优化搜索。但该类算法一般需要模拟器的源代码，需要大量迭代才能收敛，且对初始解的选取非常敏感。因此，该类算法在大规模已开发的老油田中难以实施。

(3)质量图法。该类方法有两种用途：一是通过定性图件中的高值区域直观确定最优的调整井位，避免了繁琐的数值模拟；二是用作进化算法的约束条件，指导进化算法在合理区域内布井，以降低搜索空间、提高收敛速度。由于该类方法计算效率较高，且因其是基于地质模型和生产动态历史构建，计算精度也较高，因此被广泛应用于高含水期复杂地质条件下的井位优化问题。尽管目前的质量图法考虑了油藏某些静态和动态属性，但未考虑油藏驱动机理、现有井网泄流和注入水波及体积对调整井位的影响，对优化的井位造成了一定误差。

发明内容

鉴于目前调整井位优化中存在的问题，本发明提供一种基于进化算法的高含水油田确定井位的方法，其技术方案如下：

基于进化算法的高含水油田确定井位的方法，包括以下步骤：

S1：建立预定油区的油藏模型，将所述油藏模型划分网格；

S2：获取所述预定油区的油藏资料，根据所述油藏资料按照历史拟合得出剩余油潜力富集区；

S3：计算所述剩余油潜力富集区中每个网格的井位约束值，所述井位约束值表示剩余油分布和注入水波及体积，得出最大井位约束值对应的网格，所述网格表示井位优势区域；

S4：将所述优势区域内的油井运用进化算法进行迭代计算其适应度函数，并利用克里金算法对适应度函数进行评估，所述进化算法根据适应度函数进行自适应调整，所述克里金算法考虑多维空间、块金效应和局部权重；

S5：根据进化算法迭代计算，得到最大适应度函数值，从而确定井位。

如上所述的基于进化算法的高含水油田确定井位的方法，S3中确定井位优势区域的步骤是：

S31：确定质点传播时间，所述质点传播时间利用下式计算：

式中：TOF(x,y,z)为质点到达位置(x,y,z)需要的传播时间，ψ、χ称为双流函数，代表梯度；为储层孔隙度，s为沿流线的距离；

S32：在质点传播时间的基础上，计算可动油体积指数MOVI、油相流动能力指数CFI和流线动态度量DM，上述三个指数通过下式计算得到：

可动油体积指数MOVI：

MOVI＝(S_ro-S_orw)×PV×NTG

油相流动能力指数CFI：

CFI＝K×K_ro×DZ×NTG

流线动态度量DM：

DM＝(TOF_p+TOF_i)

式中：S_ro：预测开始时剩余油饱和度分布，％；

S_orw：预测结束时残余油饱和度分布，％；

PV：孔隙体积，m³；

NTG：净毛比，表示有效厚度与总厚度的比值，无量纲；

K：绝对渗透率分布，10^-3μm²；

Kro：油相相对渗透率，指在油水两相系统中，油相渗透率与绝对渗透率的比值，无量纲；

DZ：纵向网格的尺寸，m；

上述参数均可通过获取油藏资料得到；

TOF_i为以注水井为出发点计算的质点传播时间，TOF_p为以生产井为出发点计算的质点传播时间，DM为总传播时间；

S33：对S12得到的可动油体积指数MOVI、油相流动能力指数CFI和流线动态度量DM规整化，规整化的具体步骤如下：

①通过模型属性累计分布频率图上10％和90％确定属性的最大值X_max和最小值X_min；

②按照下式对各属性进行规整化：

$X=\left\{\begin{array}{cc}0,& X<X_{\min }\\ \frac{X-X_{\max }}{X_{\max }-X_{\min }},& X_{\min }\&le;X<X_{\max }\\ 1,& X\&GreaterEqual;X_{\max }\end{array}\right.$

S34：基于上述规整化的三个指数，利用下式计算动态潜力指数DOI：

$\mathrm{DOI}=\sqrt[3]{\mathrm{MVOI}\&times;\mathrm{CFI}\&times;\mathrm{DM}}$

动态潜力指数DOI是剩余油分布状况和井网波及状况的综合反映，DOI越接近1的区域，表明该区域剩余油储量越大、原油流动性越好，而且尚未被注水波及，因而这些区域成为钻调整井的优势区域。

如上所述的基于进化算法的高含水油田确定井位的方法，S4中考虑块金效应和局部权重的克里金算法的步骤是：

S41：考虑多维空间，利用下式计算经验协方差：

$\mathrm{Cov}{\left(h\right)}_1=\frac{\mathrm{\&Sigma;}\left[z\left(x_i\right)z(x_i+h)\right]}{N\left(h\right)-1}-\frac{\mathrm{\&Sigma;z}\left(x_i\right)\mathrm{\&Sigma;z}(x_i+h)}{N\left(h\right)[N\left(h\right)-1]}$

其中：Cov(h)₁表示距离h以内样本点的经验协方差函数值，z(x_i)为与样本点的自适应度函数有关函数值，N(h)表示距离h以内的总样本点个数；

S42：考虑块金效应和局部权重，利用下式计算理论协方差：

$\mathrm{Cov}{\left(h\right)}_2=\mathrm{\&rho;}\left(h\right){\widehat{\mathrm{\&sigma;}}}^2$

$\mathrm{\&rho;}\left(h\right)=2(1-\mathrm{nu}){\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right){\left(\frac{\mathrm{bh}}{2}\right)}^{s-1}{K}_{s-1}\left(\mathrm{bh}\right)\mathrm{\&Gamma;}{(s-1)}^{-1}$

其中：Cov(h)₂为理论协方差，ρ(h)为相距h任意两点的空间相关函数，为克里金估计方差，nu为块金值，K为第二阶修正贝塞尔函数，Г为伽马函数，参数b、s为表征相关性函数形状的实数，通过求解如下的无约束优化问题可得到参数b、s：

$\underset{b,s}{\min }\mathrm{\&Sigma;}{[2{\widehat{\mathrm{\&sigma;}}}^2(1-\mathrm{nu}){\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right){\left(\frac{\mathrm{bh}}{2}\right)}^{s-1}{K}_{s-1}\left(\mathrm{bh}\right)\mathrm{\&Gamma;}{(s-1)}^{-1}-\mathrm{Cov}\left(h\right)]}^2$

ρ_θ(h)为局部权重函数，利用下式计算：

${\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right)=\left\{\begin{array}{cc}(1-\frac{h}{\mathrm{\&theta;}})\cos \left(\frac{\mathrm{\&pi;h}}{\mathrm{\&theta;}}\right)+\frac{\sin \left(\frac{\mathrm{\&pi;h}}{\mathrm{\&theta;}}\right)}{\mathrm{\&pi;}},& h<\mathrm{\&theta;}\\ 0,& h\&GreaterEqual;\mathrm{\&theta;}\end{array}\right.$

其中：θ为支撑度；

S43：根据计算得到的经验协方差和理论协方差，利用最小二乘法拟合协方差；

S44：根据协方差与变差函数的线性关系，得出变差函数；

S45：采用全主元消去法求解线性方程组，求解得到权重系数，完成克里金算法的建立。

如上所述的基于进化算法的高含水油田确定井位的方法，S5中所述的进化算法包括遗传算法，所述遗传算法的交叉概率和变异概率根据自适应函数进行自动调节：

利用下式使交叉概率根据适应度值变化自动调节：

$p_x=\left\{\begin{array}{cc}\frac{a(F_{\max }-F_x)}{F_{\max }-F_{\mathrm{avg}}},& F_x\&GreaterEqual;F_{\mathrm{avg}}\\ b,& F_x<F_{\mathrm{avg}}\end{array}\right.$

式中：F_max为种群中最大适应值，F_avg为种群平均适应值，F_x为要交叉两个个体中较大的适应度值，a＝0.8，b＝1；

利用下式使变异概率根据适应度值变化自动调节：

$p_m=\left\{\begin{array}{cc}\frac{c(F_{\max }-F_m)}{F_{\max }-F_{\mathrm{avg}}},& F_m\&GreaterEqual;F_{\mathrm{avg}}\\ d,& F_m<F_{\mathrm{avg}}\end{array}\right.$

式中：F_max为种群中最大适应值，F_avg为种群平均适应值，F_m为要变异的个体的适应度值，c＝0.01，d＝0.05。

如上所述的基于进化算法的高含水油田确定井位的方法，S5中所述的进化算法包括粒子群算法，所述粒子群算法的粒子位移更新和飞行速度根据自适应函数进行自动调节：

$x_{i,j}^{k+1}=x_{i,j}^k+v_{i,j}^{k+1}\mathrm{\&Delta;t}$

式中：分别为进化到第k+1个时间步时第i个粒子的位置和速度的第j个元素；为进化到第k个时间步时第i个粒子本身找到的最优解的第j个元素；为进化到第k个时间步时整个种群目前找到的全局最优解的第j个元素；c₁、c₂为正的学习因子；r₁、r₂为[0，1]之间均匀分布的随机数；Δt为时间增量，为压缩系数：

ω为惯性权重：

$w=\left\{\begin{array}{cc}{w}_{\min }-\frac{(w_{\max }-w_{\min })(F-F_{\min })}{(F_{\mathrm{avg}}-F_{\min })},& F\&le;F_{\mathrm{avg}}\\ w_{\max }& F>F_{\mathrm{avg}}\end{array}\right.$

式中：ω_min、ω_max分别为惯性权重的最大值和最小值，F为粒子当前的适应值；F_avg、F_min为当前所有粒子的平均目标适应值和最小适应值。

借由以上的技术方案，本发明的有益效果在于：本发明的基于进化算法的高含水油田确定井位的方法以兼顾精度和求解速度的进化算法和克里金算法为主要优化工具，以累计产油量为目标函数，并构建了综合考虑老井干扰、断层边界、无效储层、井距限制和体现剩余油分布和注入水波及体积的动态潜力指数缩小进化算法的搜索空间，从而使优化计算过程快捷、省时，收敛精度和速度大大增加。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施方式的基于进化算法的高含水油田确定井位的方法的流程图；

图2a为待调整井位的油区的剩余油栅状图；

图2b为待调整井位的油区的剩余油剖面分布图；

图2c为待调整井位的主力油区一的剩余油潜力富集区；

图2d为待调整井位的主力油区二的剩余油潜力富集区；

图3a为质点传播时间示意图；

图3b为质点传播时间平面分布；

图4a为可动油体积指数MVOI的模型属性累计分布频率图；

图4b为油相流动能力指数CFI的模型属性累计分布频率图；

图4c为流线动态度量DM的模型属性累计分布频率图；

图5a为规整化后的可动油体积指数MVOI的栅状图；

图5b为规整化后的油相流动能力指数CFI的栅状图；

图5c为规整化后的流线动态度量DM的栅状图；

图6a为规整化后的动态潜力指数DOI分布图；

图6b为规整化后规整化后动态潜力指数DOI纵向叠加平面分布图；

图7a为动态潜力指数DOI>0.2确定的井位优势区域；

图7b为动态潜力指数DOI>0.4确定的井位优势区域；

图7c为动态潜力指数DOI>0.6确定的井位优势区域；

图7d为动态潜力指数DOI>0.8确定的井位优势区域；

图8a为现有的GA算法得到的累计产油量进化曲线；

图8b为现有的GA算法得到的最优井位坐标；

图9a为现有的PSO算法得到的累计产油量进化曲线；

图9b为现有的PSO算法得到的最优井位坐标；

图10a为本发明实施方式的HGA算法得到的累计产油量进化曲线；

图10b为本发明实施方式的HGA算法得到的最优井位坐标；

图11a为本发明实施方式的HPSO算法得到的累计产油量进化曲线；

图11b为本发明实施方式的HPSO算法得到的最优井位坐标；

图12为DOI>0.6约束条件下不同井位优化算法优化效果对比。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参见图1。基于进化算法的高含水油田确定井位的方法，包括：

S1：建立预定油区的油藏模型，将所述油藏模型划分网格；

S2：获取所述预定油区的油藏资料，根据所述油藏资料按照历史拟合得出剩余油潜力富集区；

S3：计算所述剩余油潜力富集区中每个网格的井位约束值，所述井位约束值表示剩余油分布和注入水波及体积，得出最大井位约束值对应的网格，所述网格表示井位优势区域；

S4：将所述优势区域内的油井运用进化算法进行迭代计算其适应度函数，并利用克里金算法对适应度函数进行评估，所述进化算法根据适应度函数进行自适应调整，所述克里金算法考虑多维空间、块金效应和局部权重；

S5：根据进化算法迭代计算，得到最大适应度函数值，从而确定井位。

借由以上的技术方案，本发明实施方式的有益效果在于：本发明的基于进化算法的高含水油田确定井位的方法以兼顾精度和求解速度的进化算法和克里金算法为主要优化工具，以累计产油量为目标函数，并构建了综合考虑老井干扰、断层边界、无效储层、井距限制和体现剩余油分布和注入水波及体积的动态潜力指数缩小进化算法的搜索空间，从而使优化计算过程快捷、省时，收敛精度和速度大大增加。

下面以本发明在某高含水油田某区块的应用为例，进一步说明本发明实施例的基于进化算法的高含水油田确定井位的方法。

目前该油田这一区块已水驱开发20年，地下油水关系较为复杂。油田注采井网进行了多次的加密调整后，水驱控制程度(仅55％)和动用程度(20％～35％)都依然较低，注入水单层突进现象仍比较严重，采油速度低，产油量递减快。需要结合精细地质研究结果，进一步完善井网和注采关系，改善高含水后期的开发效果。采用多种方法测算，该区块合理井距为150-200m，目前平均井距为220m，因此目前仍然存在可观的井位调整潜力。

本发明针对上述油田区块基于进化算法调整井位的方法步骤如下：

S1：建立预定油区的油藏模型，将所述油藏模型划分网格；

本发明方法选取中部物性较好、地层倾角约16°的主力开发单元进行研究。对选取的区块建立油藏模型，并对其进行网格划分，网格为77×49×15，总节点数为56595个。

S2：获取所述预定油区的油藏资料，根据所述油藏资料按照历史拟合得出剩余油潜力富集区；

采用历史拟合辅助以手动修改参数，得到该区块剩余油的整体分布状况如图2a、2b所示，以及主力层位的剩余油分布状况如图2c、2d所示。

油田开采进入中高含水期后，剩余油分布零散的原因是由于存在大量小断层，油水关系极为复杂，加上砂体规模小、变化快的特征，受断层切割后，决定了储层平面上连通程度不高，注采井网不规则，开发后期剩余油分布较零散。这些特征为后期调整井位的优化带来了困难。

S3：计算所述剩余油潜力富集区中每个网格的井位约束值，所述井位约束值表示剩余油分布和注入水波及体积，得出最大井位约束值对应的网格，所述网格表示井位优势区域；

根据上述区块剩余油潜力富集区分布状况，计算研究区动态潜力指数DOI，根据动态潜力指数DOI确定井位优势区域。具体方法和步骤如下：

S31、确定质点传播时间

流线模拟是水驱开发管理和高含水后期精细油藏数值模拟的有利工具，其方法是将三维模拟还原为一系列的一维流线模型，在还原过程中计算出流体流动数据，具有处理更大数量级数据的优势。由于在驱替过程中具有保持明显的驱替前缘和减少网格方位影响的特点，使流动计算不受地质网格的限制，有效避免了有限差分网格造成的数值弥散和网格取向效应的影响，流线模拟对于高含水老油田复杂油水关系的刻画更加精确。

流线模拟中的一个重要参数为质点传播时间(Time Of Flight，TOF)，表示某质点从生产井或注水井出发到达空间中某点(坐标x、y和z的函数)所需要的时间，是一个与空间有关的物理量。

根据流线模拟理论，质点传播时间可定义为：

式中：TOF(x,y,z)为质点到达位置(x,y,z)需要的传播时间，单位是天，ψ、χ称为双流函数，代表梯度；为储层孔隙度；s为沿流线的距离，即质点移动形成的轨迹线。通常采用流线模拟或有限差分模拟器计算的三维速度场还原计算一维流线模型质点传播时间。

质点传播时间可以刻画地下油藏不同流体相的流动快慢，并可视化水驱前缘。本发明通过对有限差分网格内所有流线TOF值求算术平均，并将平均TOF值作为该网格的TOF值，即可得到TOF分布如图3a、3b所述。TOF分布是反映井网下注水井波及状况和生产井泄油状况的良好指数。

S32、计算可动油体积指数MOVI、油相流动能力指数CFI和流线动态度量DM

在质点传播时间的基础上，引入表征可动剩余油储量的可动油体积指数MOVI、表征油相流动能力的油相流动能力指数CFI和表征注入水波及、生产井泄流状况的流线动态度量DM。上述三个指数通过下式计算得到：

(1)可动油体积指数MOVI

MOVI＝(S_ro-S_orw)×PV×NTG   (2)

(2)油相流动能力指数CFI

CFI＝K×K_ro×DZ×NTG   (3)

(3)流线动态度量DM

DM＝(TOF_p+TOF_i)   (4)

式中：S_ro：预测开始时剩余油饱和度分布，％；

S_orw：预测结束时残余油饱和度分布，％；

PV：孔隙体积，m³；

NTG：净毛比，表示有效厚度与总厚度的比值，无量纲；

K：绝对渗透率分布，10^-3μm²；

Kro：油相相对渗透率，指在油水两相系统中，油相渗透率与绝对渗透率的比值，无量纲；

DZ：纵向网格的尺寸，m。

上述参数均可通过获取油藏资料得到。

剩余油中有可动和不可动两部分，MOVI指可动油部分的分布情况；而CFI表征可动油部分的流动条件；TOF_i为以注水井为出发点计算的质点传播时间，体现了注入水的波及体积；TOF_p为以生产井为出发点计算的质点传播时间，反映了生产井的泄流状况；DM为总传播时间，其体现了生产井泄流体积和注水井波及体积的综合影响。

S33、对上述的三个指数规整化

由于可动油体积指数MOVI、油相流动能力指数CFI和流线动态度量的量纲DM不同，数值差别可能大到1～2个数量级。为了使各个指数具有可对比性，需要根据各指数的最大最小值将各指数进行规整化。规整化的具体步骤如下：

①通过模型属性累计分布频率图上10％和90％确定属性的最大值X_max和最小值X_min；

②然后按照下式对各属性进行规整化：

$X=\left\{\begin{array}{cc}0,& X<X_{\min }\\ \frac{X_{\max }-X}{X_{\max }-X_{\min }},& X_{\min }\&le;X<X_{\max }\\ 1,& X\&GreaterEqual;X_{\max }\end{array}\right.---\left(5\right)$

本发明实施方式中，根据主力断块孔隙度、渗透率、储层厚度分布及历史拟合后的剩余油分布，计算出可动油体积指数MVOI、油相流动能力指数CFI以及流线动态度量DM三个指数，得出可动油体积指数MVOI、油相流动能力指数CFI以及流线动态度量DM的模型属性累计分布频率图图4a、4b、4c所示。统计了模型属性累计分布频率图上10％和90％确定属性的最大值X_max和最小值X_min如表1所示。

表1 各指数10％和90％累计概率对应的属性值

规整化后三个指数栅状图如图5a、5b、5c所示。

S34、计算动态潜力指数DOI

基于上述规整化的三个指数，利用下式计算反映井位优势区域的动态潜力指数DOI：

$\mathrm{DOI}=\sqrt[3]{\mathrm{MVOI}\&times;\mathrm{CFI}\&times;\mathrm{DM}}---\left(6\right)$

动态潜力指数DOI是剩余油分布状况和井网波及状况的综合反映，DOI越接近1的区域，表明该区域剩余油储量越大、原油流动性越好，而且尚未被注水波及，因而这些区域成为钻调整井的优势区域。以动态潜力指数DOI作为井位优化的优势区域具有如下优点：高潜力区识别精度更高；考虑注入水波及体积的新指数对可采剩余油的识别更加合理。

根据规整化后三个指数，计算得到动态潜力指数DOI分布图如图6a所示，对动态潜力指数DOI进行厚度加权处理，得到纵向叠加平面分布图如图6b所示。

动态潜力指数分布能综合体现剩余油潜力、流动能力和井网波及状况，对井位优势区域的高潜力区识别精度更高。

根据计算式(6)，动态潜力指数DOI越大的区域，越适宜钻调整井。但实际操作中难以确定井位优势区域与动态潜力指数DOI的定量关系，因此本发明分别做出DOI>0.2、0.4、0.6、0.8四个级别动态潜力指数分布图作为调整井位优化的新约束条件。根据动态潜力指数分布图，可以确定不同级别的井位优势区域，如图7a、7b、7c、7d所示。

不同级别动态潜力指数分布图对井位优化算法的搜索空间和最优井位均有影响。DOI值较小的质量图，搜索空间较大，井位优化算法的不确定性较高，甚至有可能使得算法收敛到局部解而出现“早熟”。通过上述约束条件的施加，大大降低了算法的搜索空间。

S4：将所述优势区域内的油井运用进化算法进行迭代计算其适应度函数，并利用克里金算法对适应度函数进行评估，所述进化算法根据适应度函数进行自适应调整，所述克里金算法考虑多维空间、块金效应和局部权重；

为了解决进化算法调用数值模拟器次数多、计算时间长的问题，本发明将地质统计学中的普通克里金算法推广至多维，并引入块金效应、局部权重对其进行改进。改进的克里金算法以拟合协方差代替现有克里金算法中的拟合变差函数的过程。

地质统计学中的普通克里金方程可用如下矩阵形式表示：

[K][λ]＝[D]   (7)

其中：

$\left[K\right]=\left[\begin{array}{ccccc}{\mathrm{\&gamma;}}_{11}& {\mathrm{\&gamma;}}_{12}& ...& {\mathrm{\&gamma;}}_{1n}& 1\\ {\mathrm{\&gamma;}}_{21}& {\mathrm{\&gamma;}}_{22}& ...& {\mathrm{\&gamma;}}_{2n}& 1\\ .& .& & .& .\\ .& .& & .& .\\ .& .& & .& .\\ {\mathrm{\&gamma;}}_{n1}& {\mathrm{\&gamma;}}_{n2}& ...& {\mathrm{\&gamma;}}_{\mathrm{nn}}& 1\\ 1& 1& ...& 1& 0\end{array}\right],\left[\mathrm{\&lambda;}\right]=\left[\begin{array}{c}{\mathrm{\&lambda;}}_1\\ {\mathrm{\&lambda;}}_2\\ .\\ .\\ .\\ {\mathrm{\&lambda;}}_u\\ \mathrm{\&mu;}\end{array}\right],\left[D\right]=\left[\begin{array}{c}\mathrm{\&gamma;}(x_1,x)\\ \mathrm{\&gamma;}(x_2,x)\\ .\\ .\\ .\\ \mathrm{\&gamma;}(x_n,x)\\ 1\end{array}\right]---\left(8\right)$

采用全主元消去法求解线性方程组，求解得到克里金系数λ＝K^-1D和拉格朗日系数μ，根据权重系数λ可计算得到普通克里金计算的估计值。进而可以求得克里金估计方差为：

${\widehat{\mathrm{\&sigma;}}}^2\left(x\right)=D^T\mathrm{\&lambda;}---\left(9\right)$

普通克里金方程中的系数矩阵隐含了变差函数。因此，构建克里金算法需要选取合理的理论变差函数，通过理论与经验变差函数拟合，确定出系数矩阵中的变差函数。

根据地质统计学理论，协方差函数与变差函数呈线性关系。

S41、考虑多维空间，计算经验协方差。

在井位优化中，构建克里金的样本为所要设计的井位的坐标及对应的累计产油量。因此，样本点为多维空间的点，其维数取决于优化的井数，一般可表示为：样本点的维数＝优化井数×2+1。例如，优化3口井时，每个个体的样本点为7×n维，如式：

$M=\left(\begin{array}{ccccccc}{X}_{11}& Y_{11}& X_{21}& Y_{21}& X_{31}& Y_{31}& Q_{c1}\\ X_{12}& Y_{12}& X_{22}& Y_{22}& X_{32}& Y_{32}& Q_{c2}\\ ...& ...& ...& ...& ...& ...& ...\\ X_{1n}& Y_{1n}& X_{2n}& Y_{2n}& X_{3n}& Y_{3n}& Q_{\mathrm{cn}}\end{array}\right)---\left(10\right)$

其中：X_ij、Y_ij分别表示第i个个体中第j口井的x、y坐标，i＝1,2...，j＝1,2,3；Q_ci为第i个个体的累计产油量，m³。

鉴于井位优化的多维性，本发明引入线性代数中欧几里得距离表示任意维空间中两点之间的真实距离。假设最大优化井数为m，即克里金算法需要评估n＝2m维空间中两个点之间的距离。根据欧式距离的公式，得到n维空间中的点A(X_1i,Y_1i,X_2i,Y_2i,…,X_mi,Y_mi)到点B(X_1i+1,Y_1i+1,X_2i+1,Y_2i+1,…,X_mi+1,Y_mi+1)之间的距离为：

$h(A,B)=\sqrt{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}[{(x_{i1}-x_{i2})}^2+{(y_{i1}-y_{i2})}^2]}---\left(11\right)$

其中：h(A,B)表示AB两点的欧氏距离；X_ji表示第i个个体中第j口井的x坐标；Y_ji表示第i个个体中第j口井的y坐标。

经验协方差函数用下式计算：

$\mathrm{Cov}{\left(h\right)}_1=\frac{\mathrm{\&Sigma;}\left[z\left(x_i\right)z(x_i+h)\right]}{N\left(h\right)-1}-\frac{\mathrm{\&Sigma;z}\left(x_i\right)\mathrm{\&Sigma;z}(x_i+h)}{N\left(h\right)[N\left(h\right)-1]}---\left(12\right)$

其中：Cov(h)₁、γ(h)表示距离h以内样本点的经验协方差函数值、经验变差函数值z(x_i)为样本点的函数值，在本发明中与自适应度函数有关N(h)表示距离h以内的总样本点个数，由于一般油藏模型具有不规则的形状，因此随着距离h变化。

S42、考虑块金效应和局部权重，计算理论协方差。

理论变差函数的选取对克里金算法预测的精度非常重要。常用的理论变差函数有：线性、球状、高斯和指数模型，在复杂井位优化问题中，采用这些函数难以达到对经验变差函数的拟合。为了解决这个问题，本发明引入考虑块金值效应、局部权重对克里金算法进行改进。基本思路如下：

$\mathrm{Cov}{\left(h\right)}_2=\mathrm{\&rho;}\left(h\right){\widehat{\mathrm{\&sigma;}}}^2---\left(13\right)$

$\mathrm{\&rho;}\left(h\right)=2(1-\mathrm{nu}){\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right){\left(\frac{\mathrm{bh}}{2}\right)}^{s-1}{K}_{s-1}\left(\mathrm{bh}\right)\mathrm{\&Gamma;}{(s-1)}^{-1}---\left(14\right)$

其中：Cov(h)₂为理论协方差，ρ(h)为相距h任意两点的空间相关函数，nu为块金值，K为第二阶修正贝塞尔函数，Г为伽马函数，参数b、s为表征相关性函数形状的实数，通过求解如下的无约束优化问题可得到参数b、s：

$\underset{b,s}{\min }\mathrm{\&Sigma;}{[2{\widehat{\mathrm{\&sigma;}}}^2(1-\mathrm{nu}){\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right){\left(\frac{\mathrm{bh}}{2}\right)}^{s-1}{K}_{s-1}\left(\mathrm{bh}\right)\mathrm{\&Gamma;}{(s-1)}^{-1}-\mathrm{Cov}\left(h\right)]}^2---\left(15\right)$

ρ_θ(h)为局部权重函数(或紧支撑相关函数)，可用下式表示：

${\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right)=\left\{\begin{array}{cc}(1-\frac{h}{\mathrm{\&theta;}})\cos \left(\frac{\mathrm{\&pi;h}}{\mathrm{\&theta;}}\right)+\frac{\sin \left(\frac{\mathrm{\&pi;h}}{\mathrm{\&theta;}}\right)}{\mathrm{\&pi;}},& h<\mathrm{\&theta;}\\ 0,& h\&GreaterEqual;\mathrm{\&theta;}\end{array}\right.---\left(16\right)$

其中：θ为支撑度。

S43、根据计算得到的经验协方差和理论协方差，利用最小二乘法拟合协方差。

S44、根据协方差与变差函数的线性关系，得出变差函数。

S45、采用全主元消去法求解线性方程组，求解得到权重系数，得到考虑多维空间、块金效应和局部权重的克里金算法。

进化算法是以达尔文的进化论思想为基础，通过模拟生物进化过程与机制的求解问题的自组织、自适应的人工智能技术。如同生物进化中的繁殖、变异、竞争和选择现象，进化算法主要通过选择、重组和变异这三种操作实现优化问题的求解。常用的进化算法有遗传算法、粒子群算法。

本发明实施方式对遗传算法和粒子群算法进行改进，利用改进后的遗传算法和粒子群算法对确定的新的剩余油潜力富集区进行迭代计算自适应函数。

本发明选用累计产油量作为井位优化的自适应函数，其函数关系式是：

F＝f(I,J,T,S,N)   (17)

式中：F为累计产油量，为待求井参数的函数；I、J分别为直井X、Y方向的网格指标；T为井的类别，注水井或生产井；S为井的状态，开启或关闭；N为总井数。井位优化问题即为采用一种稳健的优化算法确定使得F最大的参数组合(I,J,T,S,N)。

改进的遗传算法

遗传算法不是采用确定性规则，而是采用概率的变迁规则来指导他的搜索方向。

遗传算法(Genetic Algorithm，GA)是一种简单通用、鲁棒性强、应用范围广的全局随机优化算法，模仿了生物遗传学上适者生存的进化过程。编码规则、初始种群、交叉概率和变异概率是影响遗传算法收敛性和鲁棒性的重要参数，本发明提出了综合考虑可优化多变量的编码规则、基于油藏工程经验的初始种群和自适应交叉概率和变异概率的改进遗传算法。

(1)个体编码

遗传算法的运算对象是表示个体的符号串，所以必须把待求井参数编码为一种符号串。本发明采用简单易懂的二进制编码，在设置最大优化井数后，使用二进制编码同时对待优化的井的井别和井位坐标进行编码。

对于不同的网格系统，每口井的编码位数不同。以优化一口直井为例进行说明，假设每口井需要10位数编码，前2位分别表示井的状态(开井/关井)和井的类别(生产井/注水井)，后8位表示该井的网格坐标。例如：某代种群中某个体的二进制编码为1111011111，解释如下：

1：状态位，1表示开井，0表示关井，因此该井为开启状态。

1：井别位，1表示生产井，0表示注水井，因此该井为生产井。

1101：根据二进制与十进制换算，I坐标网格数I＝13。

1111：根据二进制与十进制换算，J坐标网格数J＝15。

则基因型为1111011111的个体，其所对应的表现型为位于(13,15)网格处1口开启状态的生产井。

(2)确定初始种群

初始种群对遗传算法收敛速度和准确性有重要的影响。初始种群质量不好，可能会大大增加算法的进化代数，导致计算效率降低，甚至有可能过早收敛到局部最优解导致“早熟”。本发明确定初始种群时，根据油藏工程经验给出初始种群的10％的个体，其余个体通过随机方法产生。这样的初始化方法既确保了算法的计算效率，又提高了优化井位的合理性。

(3)适应度汁算

遗传算法中以个体适应度的大小来评定各个个体的优劣程度，从而决定其遗传机会的大小。本发明实施方式中，利用克里金算法对个体的适应度进行评估。

(4)选择运算

选择运算又称为复制运算，把当前群体中适应度较高的个体按某种规则或模型遗传到下一代群体中。一般要求适应度较高的个体将有更多的机会遗传到下一代群体中。

(5)交叉运算

交叉运算是遗传算法中产生新个体的主要操作过程，以选择运算后的种群中的个体为运算基础，它以某一交叉概率相互交换某两个个体之间的部分染色体。

交叉概率是指染色体中某个基因被选择作为交叉点的概率，较高的交叉概率使得进化更加充分，算法的局部优化能力加强，但收敛速度变慢。本发明采用如下的自适应交叉概率，根据适应度值变化自动调节交叉概率：

$p_x=\left\{\begin{array}{cc}\frac{a(F_{\max }-F_x)}{F_{\max }-F_{\mathrm{avg}}},& F_x\&GreaterEqual;F_{\mathrm{avg}}\\ b,& F_x<F_{\mathrm{avg}}\end{array}\right.---\left(18\right)$

式中：F_max为种群中最大适应值，本发明为累计产油量值；F_avg为种群平均适应值；F_x为要交叉两个个体中较大的适应度值；a、b为常数，本发明中a＝0.8，b＝1。

(6)变异运算

变异运算是对个体的某一个或某一些基因座上的基因值按某一较小的概率进行改变，以选择交叉后的种群中的个体为运算基础，它也是产生新个体的一种操作方法。

高变异概率对种群进化的多样性有积极作用，但可能破坏已有的种群。在进化早期，应采用低变异概率；在进化晚期，应选用高变异概率。这样既避免因为高变异概率掩盖早期个体的交叉作用，又会在进化陷入提前收敛而对解不再有所改善，高变异概率使种群更具进化能力和多样性。

考虑到最优变异概率随适应度的变化，本发明引入自适应变异概率函数，使得变异概率根据适应度自动设定。本发明在优化早期选用变异概率0.01，当进化到一定代数，低变异概率不能改善解时，选择高变异概率0.05。自适应变异概率通过下式计算：

$p_m=\left\{\begin{array}{cc}\frac{c(F_{\max }-F_m)}{F_{\max }-F_{\mathrm{avg}}},& F_m\&GreaterEqual;F_{\mathrm{avg}}\\ d,& F_m<F_{\mathrm{avg}}\end{array}\right.---\left(19\right)$

式中：F_m为要变异的个体的适应度值；c、d为常数，本发明中c＝0.01，d＝0.05；其余参数同前。

本发明实施方式，采用适应度函数对交叉概率和变异概率进行自适应，得到可根据适应度函数对进化算法进行自适应的改进进化算法。运用改进的进化算法进行迭代计算适应度函数，融合建立的克里金算法对适应度函数进行评估，形成了混合遗传算法(HGA)。

改进的粒子群算法

与成熟的遗传算法相比，粒子群算法(Particle Swarm Optimization，PSO)是近年来发展起来的一种新的进化算法。不同的是它没有遗传算法的交叉和变异，而是粒子在解空间追随最优粒子进行搜索，其优势在于简单容易实现并且没有许多参数需要调整。为了增强PSO的全局探索能力和局部开发能力，本发明同时引入压缩因子和非线性惯性权重系数对现有粒子群算法进行改进。

粒子群算法中粒子位移更新公式为：

$x_{i,j}^{k+1}=x_{i,j}^k+v_{i,j}^{k+1}\mathrm{\&Delta;t}---\left(20\right)$

为有效控制粒子的飞行速度，本发明引入考虑压缩因子的粒子速度更新公式：

式中：分别为进化到第k+1个时间步时第i个粒子的位置和速度的第j个元素；为进化到第k个时间步时第i个粒子本身找到的最优解的第j个元素；为进化到第k个时间步时整个种群目前找到的全局最优解的第j个元素；c₁、c₂为正的学习因子；r₁、r₂为[0，1]之间均匀分布的随机数；Δt为时间增量，为保证算法稳定，Δt设置为1。为压缩系数，定义为：

一般情况下，设置k为1，学习因子c₁＝c₂＝2.05，此时C为4.1，计算压缩系数为0.7298。

为了均衡粒子的全局探索能力和局部开发能力，本发明引入非线性的惯性权重系数公式：

$w=\left\{\begin{array}{cc}{w}_{\min }-\frac{(w_{\max }-w_{\min })(F-F_{\min })}{(F_{\mathrm{avg}}-F_{\min })},& F\&le;F_{\mathrm{avg}}\\ w_{\max }& F>F_{\mathrm{avg}}\end{array}\right.---\left(23\right)$

式中：ω为惯性权重；ω_min、ω_max分别为惯性权重的最大值和最小值，本发明设置ω_min＝0.4、ω_max＝0.9；F为粒子当前的目标函数值；F_avg、F_min为当前所有粒子的平均目标函数值和最小目标函数值。本发明实施方式中，一个粒子代表一组布井方案，具体为井数、每口井的井位。

本发明实施方式，采用适应度函数对惯性权重进行自适应，得到可根据适应度函数对粒子群算法进行自适应的改进粒子群算法。运用改进的粒子群算法进行迭代计算适应度函数，融合建立的克里金算法对适应度函数进行评估，形成了混合遗传算法(HPSO)。

S5：根据进化算法迭代计算，得到最大适应度函数值，从而确定井位。

以求得的动态潜力指数质量图作为约束条件，以油区经过历史拟合校正的数值模型为基础，采用两种现有井位优化算法-遗传算法(GA)和粒子群算法(PSO)与本发明得到的混合遗传算法(HGA)和混合遗传算法(HPSO)对主力断块高含水后期最优调整井数、井别及井位进行优化。设置最大进化代数50，初始种群的个体规模为50，约束条件采用四个级别的动态潜力指数约束，根据油藏工程分析设置最大井数为3口。新钻生产井定产量50m³/d生产，最小井底流压设为15MPa；新钻注水井定注入量30m³/d注水，最大注水压力为50MPa。优化时间为5年(从历史拟合结束后的2012年3月计算到2017年3月)。表2为不同算法在不同约束条件，即DOI>0.2、0.4、0.6、0.8的优化结果。

表2 不同井位约束条件下不同算法的优化3口井的计算结果

从表2中可以看出，以本发明提出的动态潜力指数DOI作为井位优化约束条件，提高了井位优化算法的搜索效率。不同动态潜力指数约束条件，计算时间和收敛速度略有不同。由于在约束条件DOI>0.2下，算法搜索空间更大，导致进化过程中随机搜索到的无效井位增多，因此计算效率较低、计算时间较长。

图8～11分别为不同井位约束条件下各算法优化井位的累积产油量进化曲线及优化井位的坐标。

如图8a和8b所示，现有的GA算法最优解在DOI>0.6约束下得到，最优井位为(21,18)和(8,9)，最大累计产油量为15.54×10⁴m³；如图9a和9b所示，现有的PSO算法最优解也在DOI>0.6约束下得到，最优井位为(64,12)和(6,28)，最大累计产油量为15.69×10⁴m³。

如图10a和10b所示，本发明的HGA算法最优解在DOI>0.6约束下在第28代达到，总计算时间为5小时，搜索到最优井数2口，均为生产井，最优井位为(62,16)和(6,9)。

如图11a和11b所示，本发明的HPSO算法最优解同样在DOI>0.6时达到，但仅需进化到第25代算法即可收敛，总计算时间为5小时，搜索到最优井数2口，均为生产井，最优井位为(60,14)和(6,7)。

相对于现有两种井位优化算法，本发明提出的两种算法计算效率更高，计算井位的累计产油量更高，如图12所示，5年累计产油量高2.04×10⁴m³，表明本发明方法是可靠的。本发明的HPSO算法和HGA算法计算时间相同，累计产油量差值在计算误差以内，HPSO略高，最优井位相距不远。因此在大致范围内打井均是合理的，表明本发明提出的两种算法均为计算速度快、结果可靠的井位优化方法，鉴于本发明的HPSO算法累计产油量较高，因此优选为最终井位优化算法。

综合分析确定，研究区最优调整井数为2口，均为生产井，最优井位在地质模型中的网格分别为(60,14)和(6,7)，5年后累计产油量可达17.73×10⁴m³，地质储量采出程度达到26.1％。

本发明理论基础较强，针对高含水老油田复杂的油水关系等状况提供了一种确定高含水老油田确定井位的方法，与现有井位优化方法相比，该方法能获得更高的累计产油量、水驱波及体积更大，且优化计算过程快捷、省时，收敛速度大大增加。通过实例分析表明，本发明确定调整井位合理有效，对高含水老油田调整井数井位的合理布署具有应用意义，为老油田高含水阶段进一步挖掘剩余潜力提供了科学依据和计算方法。

以上所述仅为本发明的几个实施例，本领域的技术人员依据申请文件公开的内容可以对本发明实施例进行各种改动或变型而不脱离本发明的精神和范围。

Claims (5)

1.基于进化算法的高含水油田确定井位的方法，其特征在于，包括以下步骤：

S1：建立预定油区的油藏模型，将所述油藏模型划分网格；

S2：获取所述预定油区的油藏资料，根据所述油藏资料按照历史拟合得出剩余油潜力富集区；

S3：计算所述剩余油潜力富集区中每个网格的井位约束值，所述井位约束值表示剩余油分布和注入水波及体积，得出最大井位约束值对应的网格，所述网格表示井位优势区域；

S4：将所述优势区域内的油井运用进化算法进行迭代计算其适应度函数，并利用克里金算法对适应度函数进行评估，所述进化算法根据适应度函数进行自适应调整，所述克里金算法考虑多维空间、块金效应和局部权重；

S5：根据进化算法迭代计算，得到最大适应度函数值，从而确定井位。

2.如权利要求1所述的基于进化算法的高含水油田确定井位的方法，其特征在于：S3中确定井位优势区域的步骤是：

S31：确定质点传播时间，所述质点传播时间利用下式计算：

式中：TOF(x,y,z)为质点到达位置(x,y,z)需要的传播时间，单位是天，ψ、χ称为双流函数，代表梯度；为储层孔隙度，s为沿流线的距离；

S32：在质点传播时间的基础上，计算可动油体积指数MOVI、油相流动能力指数CFI和流线动态度量DM，上述三个指数通过下式计算得到：

可动油体积指数MOVI：

MOVI＝(S_ro-S_orw)×PV×NTG

油相流动能力指数CFI：

CFI＝K×K_ro×DZ×NTG

流线动态度量DM：

DM＝(TOF_p+TOF_i)

式中：S_ro：预测开始时剩余油饱和度分布，％；

S_orw：预测结束时残余油饱和度分布，％；

PV：孔隙体积，m³；

NTG：净毛比，表示有效厚度与总厚度的比值，无量纲；

K：绝对渗透率分布，10^-3μm²；

Kro：油相相对渗透率，指在油水两相系统中，油相渗透率与绝对渗透率的比值，无量纲；

DZ：纵向网格的尺寸，m；

上述参数均可通过获取油藏资料得到；

TOF_i为以注水井为出发点计算的质点传播时间，TOF_p为以生产井为出发点计算的质点传播时间，DM为总传播时间；

S33：对S12得到的可动油体积指数MOVI、油相流动能力指数CFI和流线动态度量DM规整化，规整化的具体步骤如下：

①通过模型属性累计分布频率图上10％和90％确定属性的最大值X_max和最小值X_min；

②按照下式对各属性进行规整化：

$X=\left\{\begin{array}{cc}0,& X<X_{\min }\\ \frac{X-X_{\max }}{X_{\max }-X_{\min }},& X_{\min }\&le;X<X_{\max }\\ 1,& X\&GreaterEqual;X_{\max }\end{array}\right.$

S34：基于上述规整化的三个指数，利用下式计算动态潜力指数DOI：

$\mathrm{DOI}=\sqrt[3]{\mathrm{MVOI}\&times;\mathrm{CFI}\&times;\mathrm{DM}}$

动态潜力指数DOI是剩余油分布状况和井网波及状况的综合反映，DOI越接近1的区域，表明该区域剩余油储量越大、原油流动性越好，而且尚未被注水波及，因而这些区域成为钻调整井的优势区域。

3.如权利要求1所述的基于进化算法的高含水油田确定井位的方法，其特征在于：S4中考虑块金效应和局部权重的克里金算法的步骤是：

S41：考虑多维空间，利用下式计算经验协方差：

$\mathrm{Cov}{\left(h\right)}_1=\frac{\mathrm{\&Sigma;}\left[z\left(x_i\right)z(x_i+h)\right]}{N\left(h\right)-1}-\frac{\mathrm{\&Sigma;z}\left(x_i\right)\mathrm{\&Sigma;z}(x_i+h)}{N\left(h\right)[N\left(h\right)-1]}$

其中：Cov(h)₁表示距离h以内样本点的经验协方差函数值，z(x_i)为与样本点的自适应度函数有关函数值，N(h)表示距离h以内的总样本点个数；

S42：考虑块金效应和局部权重，利用下式计算理论协方差：

$\mathrm{Cov}{\left(h\right)}_2=\mathrm{\&rho;}\left(h\right){\widehat{\mathrm{\&sigma;}}}^2$

$\mathrm{\&rho;}\left(h\right)=2(1-\mathrm{nu}){\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right){\left(\frac{\mathrm{bh}}{2}\right)}^{s-1}{K}_{s-1}\left(\mathrm{bh}\right)\mathrm{\&Gamma;}{(s-1)}^{-1}$

其中：Cov(h)₂为理论协方差，ρ(h)为相距h任意两点的空间相关函数，为克里金估计方差，nu为块金值，K为第二阶修正贝塞尔函数，Г为伽马函数，参数b、s为表征相关性函数形状的实数，通过求解如下的无约束优化问题可得到参数b、s：

$\underset{b,s}{\min }\mathrm{\&Sigma;}{[2{\widehat{\mathrm{\&sigma;}}}^2(1-\mathrm{nu}){\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right){\left(\frac{\mathrm{bh}}{2}\right)}^{s-1}{K}_{s-1}\left(\mathrm{bh}\right)\mathrm{\&Gamma;}{(s-1)}^{-1}-\mathrm{Cov}\left(h\right)]}^2$

ρ_θ(h)为局部权重函数，利用下式计算：

${\mathrm{\&rho;}}_{\mathrm{\&theta;}}\left(h\right)=\left\{\begin{array}{cc}(1-\frac{h}{\mathrm{\&theta;}})\cos \left(\frac{\mathrm{\&pi;h}}{\mathrm{\&theta;}}\right)+\frac{\sin \left(\frac{\mathrm{\&pi;h}}{\mathrm{\&theta;}}\right)}{\mathrm{\&pi;}},& h<\mathrm{\&theta;}\\ 0,& h\&GreaterEqual;\mathrm{\&theta;}\end{array}\right.$

其中：θ为支撑度；

S43：根据计算得到的经验协方差和理论协方差，利用最小二乘法拟合协方差；

S44：根据协方差与变差函数的线性关系，得出变差函数；

S45：采用全主元消去法求解线性方程组，求解得到权重系数，完成克里金算法的建立。

4.如权利要求1所述的基于进化算法的高含水油田确定井位的方法，其特征在于：S5中所述的进化算法包括遗传算法，所述遗传算法的交叉概率和变异概率根据自适应函数进行自动调节：

利用下式使交叉概率根据适应度值变化自动调节：

$p_x=\left\{\begin{array}{cc}\frac{a(F_{\max }-F_x)}{F_{\max }-F_{\mathrm{avg}}},& F_x\&GreaterEqual;F_{\mathrm{avg}}\\ b,& F_x<F_{\mathrm{avg}}\end{array}\right.$

式中：F_max为种群中最大适应值，F_avg为种群平均适应值，F_x为要交叉两个个体中较大的适应度值，a＝0.8，b＝1；

利用下式使变异概率根据适应度值变化自动调节：

$p_m=\left\{\begin{array}{cc}\frac{c(F_{\max }-F_m)}{F_{\max }-F_{\mathrm{avg}}},& F_m\&GreaterEqual;F_{\mathrm{avg}}\\ d,& F_m<F_{\mathrm{avg}}\end{array}\right.$

式中：F_max为种群中最大适应值，F_avg为种群平均适应值，F_m为要变异的个体的适应度值，c＝0.01，d＝0.05。

5.如权利要求1所述的基于进化算法的高含水油田确定井位的方法，其特征在于：S5中所述的进化算法包括粒子群算法，所述粒子群算法的粒子位移更新和飞行速度根据自适应函数进行自动调节：

$x_{i,j}^{k+1}=x_{i,j}^k+v_{i,j}^{k+1}\mathrm{\&Delta;t}$

式中：分别为进化到第k+1个时间步时第i个粒子的位置和速度的第j个元素；为进化到第k个时间步时第i个粒子本身找到的最优解的第j个元素；为进化到第k个时间步时整个种群目前找到的全局最优解的第j个元素；c₁、c₂为正的学习因子；r₁、r₂为[0，1]之间均匀分布的随机数；Δt为时间增量，为压缩系数：

ω为惯性权重：

$w=\left\{\begin{array}{cc}{w}_{\min }-\frac{(w_{\max }-w_{\min })(F-F_{\min })}{(F_{\mathrm{avg}}-F_{\min })},& F\&le;F_{\mathrm{avg}}\\ w_{\max },& F>F_{\mathrm{avg}}\end{array}\right.$

式中：ω_min、ω_max分别为惯性权重的最大值和最小值，F为粒子当前的适应值；F_avg、F_min为当前所有粒子的平均目标适应值和最小适应值。