CN105631554A - 一种基于时间序列的油井油液含水率多模型预测方法 - Google Patents

Abstract

本发明涉及基于时间序列的油井油液含水率多模型预测方法，其特征在于，包括如下步骤：1)、利用历史数据建立油井油液含水率数据集为{x_i,i＝1,2,…,N}；2)、采用小波分析方法对油井油液含水率数据集{x_i,i＝1,2,…,N}中的数据进行预处理；3)、由近邻传播聚类算法将{x_i}^Wave进行分类；4)、将每个聚类中的数据由如下时间序列形式进行表示：5)、根据极端学习机算法建立每个聚类的时间序列模型并利用该时间序列模型得到预测值。其解决了现有油井油液含水率人工取样费时费力、影响生产监控和采油数据的实时性的问题。

Description

一种基于时间序列的油井油液含水率多模型预测方法

技术领域

本发明涉及石油生产领域，特别是一种基于时间序列的油井油液含水率多模型预测方法。

背景技术

油井油液的含水率是油田生产的一个重要指标，不仅关系到油井的开发寿命，还关系到企业的经济效益问题。因此，对其进行准确测量对计量油井产量、评价油藏的开采价值和采出程度、制定开采方案等具有极其重要的意义。目前对油井油液含水率的测量仍普遍采用人工取样再蒸馏的方法，由工人定期到油井进行人工采样，然后将油液送回技术部门进行实验分析，这种方法费时费力、影响生产监控和采油数据的实时性，不能根据油井的实际生产状况及时制定或调整合理的开采方案。

发明内容

本发明的目的是为了提供一种解决现有油井油液含水率人工取样费时费力、影响生产监控和采油数据的实时性的问题的基于时间序列的油井油液含水率多模型预测方法，利用已有的历史一段时期的油井油液含水率数据预测下一个时间点的含水率数据。

本发明的技术方案是：

一种基于时间序列的油井油液含水率多模型预测方法，其特征在于，包括如下步骤：

1)、利用历史数据建立油井油液含水率数据集为{x_i,i＝1,2,…,N}，所述油井油液含水率数据集内数据按照时间点的顺序排列，所述时间点的单位为日或月，根据油井油液含水率数据集中数据的排列顺序，记录油井油液含水率数据集中各数据的排列序号为{index_i,i＝1,2,…,N}；

2)、采用小波分析方法对油井油液含水率数据集{x_i,i＝1,2,…,N}中的数据进行预处理，由Mallat算法对{x_i,i＝1,2,…,N}中的数据进行三层小波分解，得到小波分解序列其中：表示对{x_i,i＝1,2,…,N}中的第i个数据进行小波分解后的第三层低频分量，和分别表示对{x_i,i＝1,2,…,N}中的第i个数据进行小波分解后的第三层、第二层和第一层高频分量；经过小波分解后，{x_i,i＝1,2,…,N}中的每一个数据由所对应的代替，小波分解后的数据集由{x_i}^Wave表示；

3)、由近邻传播聚类算法将{x_i}^Wave进行分类，再根据{x_i}^Wave所对应的数据的排列序号{index_i,i＝1,2,…,N}将油井油液含水率数据集中的数据划分到各个分类中，形成K个聚类，K≥1；

4)、将每个聚类中的数据按照各个数据的原始的排列序号{index_i,i＝1,2,…,N}进行排列，然后将每个聚类中的数据由如下时间序列形式进行表示：

其中，α＝1,2,…,K，t＝1,2,…,M，m为嵌入维数，τ为延迟时间，M＝N-(m-1)τ；将每个聚类的时间序列的输出表示为：

5)、根据极端学***均值最终得到预测值。

上述的基于时间序列的油井油液含水率多模型预测方法，步骤2)具体如下：

2.1、首先根据公式(1)对数据集{x_i,i＝1,2,…,N}进行分解：

其中，J＝0,1,2；H和G分别为分解低通滤波器和分解高通滤波器；

2.2、然后根据公式(2)进行重构：

$x_i^{{\mathrm{low}}_J}=H^{*}\·x_i^{{\mathrm{low}}_{J+1}}+G^{*}\·x_i^{{\mathrm{high}}_{J+1}}---\left(2\right)$

其中，H^*和G^*分别为H和G的对偶算子；

那么，油井油液含水率数据集{x_i}中的数据x_i经过三层小波分解后表示为：其中：i＝1,2,…,N，由构成的数据集表示为{x_i}^Wave。

上述的基于时间序列的油井油液含水率多模型预测方法，步骤3)利用近邻传播聚类算法将{x_i}^Wave进行分类的具体过程为：

3.1、由经过三层小波分解后的数据集{x_i}^Wave计算相似度矩阵S，{x_i}^Wave中的两个数据和的相似度s_jh由公式(3)进行计算：

$s_{jh}=-\left(\sqrt{\underset{l=1}{\overset{4}{\Σ}}{\left(x_{jl}^{Wave}-x_{hl}^{Wave}\right)}^2}\right)---\left(3\right)$

其中，和分别表示数据和的第l维；

那么，相似度矩阵S表示为：

其中，p为偏向参数，p的初始值取所有初始相似度值的均值；

3.2、设定吸引度矩阵R＝{r_jh}和归属度矩阵A＝{a_jh}，分别表示为：

其中，{r_jh}和{a_jh}的初始值为0；

3.3、设定最大迭代次数MaxLoop，吸引度矩阵R和归属度矩阵A按照如下公式进行迭代：

r_jh＝s_jh-max{a_jh+s_jh}(7)

r_kk＝p_kk-max{a_kh+s_kh}(9)

其中，k＝1,2,…,N，k为数据集{x_i}^Wave中数据的下标，表示候选聚类中心，j≠h≠k；

3.4、当满足最大迭代次数MaxLoop后，迭代过程终止，判断两矩阵R和A的对角线元素之和，如果满足r_kk+a_kk＞0，那么选择为聚类中心；

3.5、确定K个聚类的聚类中心后，根据公式(3)分别计算数据集{x_i}^Wave中除外的每个数据与K个聚类中心的距离，将数据集{x_i}^Wave中除外的每个数据归入与它距离最小的聚类中心所属的类中，完成数据集{x_i}^Wave中数据的分类。

上述的基于时间序列的油井油液含水率多模型预测方法，步骤5)中利用极端学习机算法建立时间序列模型的原理如下：

设有W个训练样本其中u_q为输入向量，v_q为输出向量，设定包含L个隐含层，激活函数为f(·)以及模型的训练输出表示为Q＝[c₁,c₂,…,c_q]^T，那么ELM模型由如下方程组进行描述：

其中，β_lq为第l个隐含层神经元与第q个输出神经元之间的连接权值；ω_l为隐含层神经元与输入神经元之间的连接权值；b_l为第l个隐含层神经元的偏置；

如果训练模型能以零误差逼近W个训练样本，即有那么对于式(11)有下式成立，

$\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\β}}_{lq}f({\mathrm{\ω}}_l,u_q,b_l)=v_q---\left(12\right)$

那么ELM模型的数学描述可以改写成如下矩阵形式：

Hβ＝V(13)

式(13)中，有：

H为隐含层输出矩阵，ω和b在初始化时随机给定，那么，ELM模型的训练可以转化成一个求解非线性方程最小值的问题，即：

$\left|\right|{\mathrm{H\β}}^{*}-V\left|\right|=\underset{\mathrm{\β}}{\min }\left|\right|H\mathrm{\β}-V\left|\right|---\left(15\right)$

输出权值矩阵β^*可以由下式求出，

β^*＝H⁺V(16)

其中，H⁺为隐含层输出矩阵H的Moore-Penrose广义逆；

那么，ELM的训练过程可以归纳为如下最优化问题：

$g(\mathrm{\ω},b)=min\underset{q=1}{\overset{W}{\Σ}}\left|\right|v_q-c_q\left|\right|---\left(17\right)$

其中，g(·)表示由ω和b所决定的函数，g(ω,b)表示当ω和b分别取不同值时的函数输出值，ELM训练的目的就是找到最优的β^*，使模型的训练输出值c_q和真实值v_q之间的误差最小；

对于激活函数f(·)的选择，采用高斯函数和Sigmoid函数的集成方式，f(·)的定义如下：

$f(\·)={\mathrm{\λe}}^{-\frac{{(u-1)}^2}{2{\mathrm{\σ}}^2}}+(1-\mathrm{\λ})\frac{1}{1+e^{-u}}---\left(18\right)$

其中，u表示输入向量，λ∈[0,1]为权值，σ²为高斯函数的宽度参数；

基于此，建立时间序列模型的计算步骤如下：初始化，随机生成隐含层输入权值ω，隐含层神经元的偏置b，高斯函数的宽度参数σ²和权值λ；根据公式(14)计算隐含层输出矩阵H；根据公式(16)计算输出权值矩阵β^*；根据公式(17)计算函数输出值。

上述的基于时间序列的油井油液含水率多模型预测方法，建立时间序列模型的过程中，采用改进的萤火虫算法对m、τ、ω、b、σ²和λ的取值进行优化选取，将m、τ、ω、b、σ²和λ看作一组解向量，即每个“萤火虫”为1×6维向量，其中6维变量分别为m、τ、ω、b、σ²和λ，计算步骤如下：

第一步，产生初始“萤火虫”群体，表示为：{F₁,F₂,…,F_Npop}，其中：N_pop表示初始“萤火虫”群体中“萤火虫”的数量；

第二步，计算每个“萤火虫”的亮度，表示为：{I₁,I₂,…,I_Npop}；

第三步，每个“萤火虫”向比它亮度更大的“萤火虫”移动，根据公式(22)-公式(24)进行位置更新，

Φ₁(Itenum+1)＝Φ(Itenum)+η(ε)·(Ψ_best(Itenum)-Φ(Itenum))+ξ·(rand-0.5)(22)

Φ₂(Itenum+1)＝Φ(Itenum)+η(ε)·(Ψ_secbest(Itenum)-Φ(Itenum))+ξ·(rand-0.5)(23)

$\mathrm{\Φ}(Itenum+1)=\frac{{\mathrm{\Φ}}_1(Itenum+1)+{\mathrm{\Φ}}_2(Itenum+1)}{2}---\left(24\right)$

其中，ξ表示步长，ξ∈[0,1]；rand为[0,1]上服从均匀分布的一个随机数；Ψ_best为亮度最大的“萤火虫”，Ψ_secbest为亮度大小次于Ψ_best的“萤火虫”；Itenum表示IFA算法运行过程中迭代的次数；Φ(Itenum)表示“萤火虫”在第Itenum次迭代时的位置，Φ(Itenum+1)表示“萤火虫”在第Itenum+1次迭代时的位置，Φ₁(Itenum+1)表示“萤火虫”受到亮度最大的“萤火虫”的吸引在第Itenum+1次迭代时的位置，Φ₂(Itenum+1)表示“萤火虫”受到亮度大小次于Ψ_best的“萤火虫”的吸引在第Itenum+1次迭代时的位置；

第四步，定义适应度函数fitness(m,τ,ω,b,σ²,λ)评价m、τ、ω、b、σ²和λ的取值，当m、τ、ω、b、σ²和λ取最佳值时，fitness(m,τ,ω,b,σ²,λ)函数值最小，定义变量Local_pop记录每一次迭代中m、τ、ω、b、σ²和λ的最佳值，定义变量Global_pop记录所有次迭代中m、τ、ω、b、σ²和λ的最佳值，执行完第三步后，分别更新Local_pop和Global_pop中的值；

第五步，如果达到最大迭代次数，停止迭代，分别输出Local_pop中和Global_pop中的m、τ、ω、b、σ²、λ的值；否则返回步骤第三步重新迭代。

上述的基于时间序列的油井油液含水率多模型预测方法，对于K个聚类中的每一个聚类，采用极端学习机算法建立时间序列模型的计算步骤如下：

步骤5.1初始化，设定m、τ、ω、b、σ²和λ的取值范围；

步骤5.2将每一个聚类中的数据归一化到[0,1]区间内，确定待建立模型的输入和输出分别为和

步骤5.3产生初始“萤火虫”群体，在m、τ、ω、b、σ²和λ的取值范围内随机给定m、τ、ω、b、σ²和λ的取值；

步骤5.4根据公式(14)计算隐含层输出矩阵H，根据公式(16)计算输出权值矩阵β^*，根据公式(17)计算函数输出值；

步骤5.5计算适应度函数fitness(m,τ,ω,b,σ²,λ)的值，记录Local_pop和Global_pop中的m、τ、ω、b、σ²和λ的初始值；

步骤5.6计算每个“萤火虫”的亮度，根据公式(22)-公式(24)进行位置更新；

步骤5.7重复步骤5.4，计算适应度函数fitness(m,τ,ω,b,σ²,λ)的值，更新Local_pop和Global_pop中的m、τ、ω、b、σ²、λ的值；

步骤5.8如果达到最大迭代次数，停止迭代，分别输出Local_pop中和Global_pop中的m、τ、ω、b、σ²、λ的值；否则返回步骤5.6-步骤5.7重新迭代；

步骤5.9将每一个聚类得到的最好的m、τ、ω、b、σ²和λ的值作为极端学习机算法建立时间序列模型的参数，由输入的值计算输出的值，并进行反归一化处理。

上述的基于时间序列的油井油液含水率多模型预测方法，采用单步预测，即每次预测出一个再预测下一个输出时，将加入到下一个输入的最后；每次预测出一个Λ_t，再预测下一个输出Λ_t+1时，将Λ_t加入到下一个输入Γ_t+1的最后。

上述的基于时间序列的油井油液含水率多模型预测方法，所述油井油液含水率数据集中数据的数量N一般为300-800个，以保证计算的准确性和降低计算的复杂性。

上述的基于时间序列的油井油液含水率多模型预测方法，m、τ、ω、b、σ²和λ的取值范围分别为：m∈[1,50]，τ∈[1,10]，ω∈[0,1]，b∈[0,10]，σ²∈[0.01,1000]，λ∈[0,1]。

本发明的有益效果是：

1、油井生产是一个连续的过程，一些生产参数的变化在一段时期内具有相关性，根据过去一段时间的历史数据去预测未来某个时间点的值，可以掌握油井的生产动态，及时评价油井的生产状况，并制定或调整合理的生产措施，对油田企业合理、高效的进行生产具有一定的实际意义。

2、油井油液的含水率是油田生产的一个重要指标，对原油的开采、计量、运输等都有重要的影响。油液含水率不仅关系到油井的工作寿命，还关系到油田生产企业的经济效益。因此，油液含水率的准确测量对计量油井产量、掌握油井生产动态、评价油藏的开发程度和价值、制定合理开发措施等具有十分重要的意义。

3、本发明所建立的一种基于时间序列的油井油液含水率多模型预测方法，原理简单，计算复杂性小，所需要的样本少。由聚类方法所建立的多模型结构，对于预测样本的判断，准确度高，特别是对于数据集中含有异常数据的情况具有较高的预测精度。

4、解决了现有油井油液含水率人工取样费时费力、影响生产监控和采油数据的实时性的问题，省时省力。

附图说明

图1是三层小波分解后的各维分量示意图；

图2是油井油液含水率数据集聚类后的数据分布示意图。

具体实施方式

该基于时间序列的油井油液含水率多模型预测方法，具体步骤如下：

1、利用历史数据建立油井油液含水率数据集，包含440个数据，表示为{x_i,i＝1,2,…,440}，440个数据按照时间点的顺序排列(单位为日)，根据油井油液含水率数据集中数据的排列顺序，记录油井油液含水率数据集中各数据的排列序号为{index_i,i＝1,2,…,440}。

2、采用小波分析方法对油井油液含水率数据集{x_i}中的数据进行预处理，具体为：首先由Mallat算法对{x_i}中的数据进行三层小波分解，根据公式(1)对数据集{x_i}进行分解：

$\left\{\begin{array}{c}{\mathrm{xr}}^{{\mathrm{low}}_{J+1}}=H\·x_i^{{\mathrm{low}}_J}\\ x_i^{{\mathrm{high}}_{J+1}}=G\·x_i^{{\mathrm{low}}_J}\end{array}\right.---\left(1\right)$

其中，i＝1,2,…,440；J＝0,1,2；H和G分别为分解低通滤波器和分解高通滤波器；

然后根据公式(2)进行重构：

$x_i^{{\mathrm{low}}_J}=H^{*}\·x_i^{{\mathrm{low}}_{J+1}}+G^{*}\·x_i^{{\mathrm{high}}_{J+1}}---\left(2\right)$

其中，H^*和G^*分别为H和G的对偶算子。

从而得到小波分解序列其中：表示对{x_i}中的第i个数据进行小波分解后的第三层低频分量，和分别表示对{x_i}中的第i个数据进行小波分解后的第三层、第二层和第一层高频分量；经过小波分解后，{x_i}中的每一个数据(一维)由所对应的(四维)代替，小波分解后的数据集{x_i}由{x_i}^Wave表示，i＝1,2,…,440。数据集{x_i}经过三层小波分解后的各维分量如图1所示。

3、由近邻传播聚类算法(AP)将{x_i}^Wave进行分类，再根据{x_i}^Wave所对应的数据的排列序号{index_i}将油井油液含水率数据集{x_i}中的数据划分到各分类中。具体步骤如下：

3.1、由经过三层小波分解后的数据集{x_i}^Wave计算相似度矩阵S，{x_i}^Wave中的两个数据和的相似度s_jh由如下公式进行计算：

$s_{jh}=-\left(\sqrt{\underset{l=1}{\overset{4}{\Σ}}{(x_{jl}^{Wave}-x_{hl}^{Wave})}^2}\right)---\left(3\right)$

其中，和分别表示数据和的第l维；

那么，相似度矩阵S可以表示为：

其中，N＝440，p为偏向参数，p的初始值取所有初始相似度值的均值。

3.2、设定吸引度矩阵R＝{r_jh}和归属度矩阵A＝{a_jh}，分别表示为：

其中，N＝440，{r_jh}和{a_jh}的初始值为0；

本实施例中，设定最大迭代次数MaxLoop＝500(所述最大迭代次数MaxLoop一般为300-500次)，吸引度矩阵R和归属度矩阵A按照如下公式进行迭代：

r_jh＝s_jh-max{a_jh+s_jh}(7)

$a_{jk}=\min \{0,r_{kk}+\underset{h=1}{\overset{N}{\Σ}}\{max(0,r_{hk})\left\}\right\}---\left(8\right)$

r_kk＝p_kk-max{a_kh+s_kh}(9)

$a_{kk}=\underset{h=1}{\overset{N}{\Σ}}max\{0,r_{hk}\}---\left(10\right)$

其中，k＝1,2,…,440，k为数据集{x_i}^Wave中数据的下标，表示候选聚类中心，j≠h≠k。

3.3、当满足最大迭代次数MaxLoop＝500后，迭代过程终止，判断两矩阵R和A的对角线元素之和，如果满足r_kk+a_kk＞0，那么选择为聚类中心。

3.4、确定K个聚类的聚类中心后，根据公式(3)分别计算数据集{x_i}^Wave中除外的每个数据与K个聚类中心的距离，将数据集{x_i}^Wave中除外的每个数据归入与它距离最小的聚类中心所属的类中，完成数据集{x_i}^Wave中数据的分类。

本实施例中，按照K个聚类中每个数据的排列序号index_i，将油井油液含水率数据集{x_i}中的440个数据分类到K个聚类中，分别将每一个聚类中的数据按照各数据的排列序号index_i进行排列。如图2所示，油井油液含水率数据集{x_i}中的440个数据被分类到8个聚类中，每个聚类中数据的个数分别为：82(K＝1)、57(K＝2)、25(K＝3)、101(K＝4)、43(K＝5)、69(K＝6)、22(K＝7)、41(K＝8)，分别将每一个聚类中的数据按照各数据的排列序号index_i进行排列。

4、将每一个聚类中的数据由如下时间序列形式进行表示：其中，α＝1,2,…,8，t＝1,2,…,M，m为嵌入维数，τ为延迟时间，M＝440-(m-1)τ；每一个聚类的时间序列的输出表示为：将和分别作为时间序列的输入和输出，本发明采用单步预测，即每次预测出一个再预测下一个输出时，将加入到下一个输入的最后。

5、根据极端学习机算法(ELM)建立每个聚类的时间序列模型，该模型的输入和输出分别为和

(1)、基本原理：

设有W个训练样本其中u_q为输入向量，v_q为输出向量。设定包含L个隐含层(为保证计算的准确性和降低计算的复杂性，所述隐含层数量L的取值为10-30个)，激活函数为f(·)以及模型的训练输出表示为Q＝[c₁,c₂,…,c_q]^T，那么ELM模型可以由如下方程组进行描述：

$\left\{\begin{array}{c}{c}_1=\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\β}}_{11}f({\mathrm{\ω}}_1,u_1,b_1)\\ c_2=\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\β}}_{12}f({\mathrm{\ω}}_2,u_2,b_2)\\ .\\ .\\ .\\ c_q=\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\β}}_{lq}f({\mathrm{\ω}}_l,u_q,b_l)\end{array}\right.---\left(11\right)$

其中，β_lq为第l个隐含层神经元与第q个输出神经元之间的连接权值；ω_l为隐含层神经元与输入神经元之间的连接权值；b_l为第l个隐含层神经元的偏置。

如果训练模型能以零误差逼近W个训练样本，即有那么对于式(11)有下式成立，

$\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\β}}_{lq}f({\mathrm{\ω}}_l,u_q,b_l)=v_q---\left(12\right)$

那么ELM模型的数学描述可以改写成如下矩阵形式：

Hβ＝V(13)

式(13)中，有：

H为隐含层输出矩阵，ω和b在初始化时已随机给定。那么，ELM模型的训练可以转化成一个求解非线性方程最小值的问题，即：

$\left|\right|{\mathrm{H\β}}^{*}-V\left|\right|=\underset{\mathrm{\β}}{\min }\left|\right|H\mathrm{\β}-V\left|\right|---\left(15\right)$

输出权值矩阵β^*可以由下式求出，

β^*＝H⁺V(16)

其中，H⁺为隐含层输出矩阵H的Moore-Penrose广义逆。

那么，ELM的训练过程可以归纳为如下最优化问题：

$g(\mathrm{\ω},b)=min\underset{q=1}{\overset{W}{\Σ}}\left|\right|v_q-c_q\left|\right|---\left(17\right)$

其中，g(·)表示由ω和b所决定的函数，g(ω,b)表示当ω和b分别取不同值时的函数输出值。ELM训练的目的就是找到最优的β^*，使模型的训练输出值c_q和真实值v_q之间的误差最小。

对于激活函数f(·)的选择，采用高斯函数和Sigmoid函数的集成方式，f(·)的定义如下：

$f(\·)={\mathrm{\λe}}^{-\frac{{(u-1)}^2}{2{\mathrm{\σ}}^2}}+(1-\mathrm{\λ})\frac{1}{1+e^{-u}}---\left(18\right)$

其中，u表示输入向量，λ∈[0,1]为权值，σ²为高斯函数的宽度参数。

(2)、由ELM方法建立时间序列模型的过程中，m、τ、ω、b、σ²和λ的取值决定了ELM模型的计算精度，本发明采用一种改进的萤火虫算法(IFA)对m、τ、ω、b、σ²和λ的取值进行优化选取。IFA算法的数学描述如下：

Ⅰ、定义每个“萤火虫”的亮度为：

$I\left(r\right)=I_0\·e^{-{\mathrm{\γ\ϵ}}^2}---\left(19\right)$

其中，I₀表示原始亮度；γ表示光强吸收系数，为一常实数；ε表示两个“萤火虫”之间的距离，定义如下：

其中，Φ和Ψ表示任意两个“萤火虫”，φ_θ和分别表示Φ和Ψ的第θ维元素；z为Φ和Ψ的维数。

Ⅱ、定义“萤火虫”的吸引度如下公式所示：

$\mathrm{\η}\left(\mathrm{\ϵ}\right)={\mathrm{\η}}_0\·e^{-{\mathrm{\γ\ϵ}}^2}---\left(21\right)$

其中，η₀表示ε＝0处的吸引度。

Ⅲ、每个“萤火虫”受到比它亮度更大的“萤火虫”的吸引而移动，位置更新公式如下：

Φ₁(Itenum+1)＝Φ(Itenum)+η(ε)·(Ψ_best(Itenum)-Φ(Itenum))+ξ·(rand-0.5)(22)

Φ₂(Itenum+1)＝Φ(Itenum)+η(ε)·(Ψ_secbest(Itenum)-Φ(Itenum))+ξ·(rand-0.5)

(23)

$\mathrm{\Φ}(Itenum+1)=\frac{{\mathrm{\Φ}}_1(Itenum+1)+{\mathrm{\Φ}}_2(Itenum+1)}{2}---\left(24\right)$

其中，ξ表示步长，ξ∈[0,1]；rand为[0,1]上服从均匀分布的一个随机数；Ψ_best为亮度最大的“萤火虫”，Ψ_secbest为亮度大小次于Ψ_best的“萤火虫”；Itenum表示IFA算法运行过程中迭代的次数；Φ(Itenum)表示“萤火虫”在第Itenum次迭代时的位置，Φ(Itenum+1)表示“萤火虫”在第Itenum+1次迭代时的位置，Φ₁(Itenum+1)表示“萤火虫”受到亮度最大的“萤火虫”的吸引在第Itenum+1次迭代时的位置，Φ₂(Itenum+1)表示“萤火虫”受到亮度大小次于Ψ_best的“萤火虫”的吸引在第Itenum+1次迭代时的位置。

本发明在采用IFA算法对m、τ、ω、b、σ²和λ的取值进行优化选取时，将m、τ、ω、b、σ²和λ看作一组解向量，即每个“萤火虫”为1×6维向量，其中6维变量分别为m、τ、ω、b、σ²和λ。具体计算步骤如下：

第一步，产生初始“萤火虫”群体，表示为：{F₁,F₂,…,F_Npop}，其中：N_pop表示初始“萤火虫”群体中“萤火虫”的数量，取N_pop＝30；

第二步，计算每个“萤火虫”的亮度，表示为：{I₁,I₂,…,I_Npop}；

第三步，每个“萤火虫”向比它亮度更大的“萤火虫”移动，根据公式(22)-公式(24)进行位置更新；

第四步，定义如公式(25)所示的适应度函数fitness(m,τ,ω,b,σ²,λ)来评价m、τ、ω、b、σ²和λ的取值，

$fitness(m,\mathrm{\τ},\mathrm{\ω},b,{\mathrm{\σ}}^2,\mathrm{\λ})=\sqrt{\frac{1}{N}\·\underset{i=1}{\overset{N}{\Σ}}{(y_i^{real}-y_i^{predicted}(m,\mathrm{\τ},\mathrm{\ω},b,{\mathrm{\σ}}^2,\mathrm{\λ}\left)\right)}^2}---\left(25\right)$

其中，N＝440，表示第i个数据的实际输出值，表示第i个数据的模型预测值；

当m、τ、ω、b、σ²和λ取最好的值时，fitness(m,τ,ω,b,σ²,λ)函数值最小，定义变量Local_pop记录每一次迭代中m、τ、ω、b、σ²和λ的最好的值，定义变量Global_pop记录所有次迭代中m、τ、ω、b、σ²和λ的最好的值，执行完步骤第三步后，分别更新Local_pop和Global_pop中的值；

第五步，如果达到最大迭代次数MaxLoop＝500，停止迭代，分别输出Local_pop中和Global_pop中的m、τ、ω、b、σ²、λ的值；否则返回步骤4.3重新迭代。

本实施例中，对于8个聚类中的每一个聚类，在采用ELM方法建立模型的同时，由IFA方法对ELM计算过程中所涉及的m、τ、ω、b、σ²和λ进行优化选取，计算步骤如下：

步骤5.1初始化，设定m、τ、ω、b、σ²和λ的取值范围，m∈[1,50]，τ∈[1,10]，ω∈[0,1]，b∈[0,10]，σ²∈[0.01,1000]，λ∈[0,1]；

步骤5.2将每一个聚类中的数据归一化到[0,1]区间内，确定待建立模型的输入和输出分别为和

步骤5.3产生初始“萤火虫”群体，在m、τ、ω、b、σ²和λ的取值范围内随机给定m、τ、ω、b、σ²和λ的取值；

步骤5.4根据公式(14)计算隐含层输出矩阵H，根据公式(16)计算输出权值矩阵β^*，根据公式(17)计算函数输出值；

步骤5.5根据公式(25)计算适应度函数fitness(m,τ,ω,b,σ²,λ)的值，记录Local_pop和Global_pop中的m、τ、ω、b、σ²和λ的初始值；

步骤5.6计算每个“萤火虫”的亮度，根据公式(22)-公式(24)进行位置更新；

步骤5.7重复步骤5.4，根据公式(25)计算适应度函数fitness(m,τ,ω,b,σ²,λ)的值，更新Local_pop和Global_pop中的m、τ、ω、b、σ²、λ的值；

步骤5.8如果达到最大迭代次数MaxLoop＝500，停止迭代，分别输出Local_pop中和Global_pop中的m、τ、ω、b、σ²、λ的值；否则返回步骤5.6-步骤5.7重新迭代。

步骤5.9将每一个聚类得到的最好的m、τ、ω、b、σ²和λ的值作为ELM方法建立时间序列模型的参数，由输入的值计算输出的值并进行反归一化处理。

6、对于8个聚类，共有8个ELM模型，那么由输入共得到8个输出值本实施例中，分别为：60、79、77、75、73、67、76和74。认为这8个输出值是由每一个聚类中的数据所预测出的下一个值，将其排列序号分别标记为{index^α _max+1}，其中，index^α _max表示第α个聚类中最大的排列序号。根据ELM模型的计算结果，这8个输出值的排列序号分别标记为：440、437、418、441、413、439、420和426。

再将排列序号标记为{index^α _max+1}(index^α _max≠440)的7个输出值按照排列序号的大小顺序由小到大重新排列，构成新时间序列的数据集{Δ_d}(d＝1,2,…,7)。重新排列后的7个输出值的排列序号为：413、418、420、426、437、439、440，所构成新时间序列的数据集{Δ_d}＝{73,77,76,74,79,67,60}，其中第4个聚类(α＝4)中的index⁴ _max＝440。

然后将其表示为新时间序列的输入：Γ_t＝[Δ_t,Δ_t+τ,…,Δ_t+(m-1)τ]，其中，t＝1,2,…,M，M＝440-(m-1)τ，新时间序列的输出表示为：Λ_t＝Δ_t+1+(m-1)τ。将Γ_t和Λ_t分别作为新时间序列的输入和输出，根据步骤5建立新时间序列模型并计算输出值，采用单步预测，每次预测出一个Λ_t，再预测下一个输出Λ_t+1时，将Λ_t加入到下一个输入Γ_t+1的最后。本实施例中，得到新时间序列的输出值Λ_t＝72。

7、将新时间序列的输出值Λ_t＝72与排列序号标记为{index⁴ _max+1}(index⁴ _max＝440)的输出值Y_t ⁴＝75取平均值得到最终输出为73.5。

对于聚类数目K取值为1，即只有一个聚类时，由油井油液含水率数据集中的N个数据训练模型进行计算。

为了进一步说明本发明所提出方法的有效性，采用国内某油田作业区十口生产井的数据进行验证，结果如表1所示。

表1十口生产井油液含水率的预测结果

Claims (9)

1.一种基于时间序列的油井油液含水率多模型预测方法，其特征在于，包括如下步骤：

1)、利用历史数据建立油井油液含水率数据集为{x_i,i＝1,2,…,N}，所述油井油液含水率数据集内数据按照时间点的顺序排列，所述时间点的单位为日或月，根据油井油液含水率数据集中数据的排列顺序，记录油井油液含水率数据集中各数据的排列序号为{index_i,i＝1,2,…,N}；

2)、采用小波分析方法对油井油液含水率数据集{x_i,i＝1,2,…,N}中的数据进行预处理，由Mallat算法对{x_i,i＝1,2,…,N}中的数据进行三层小波分解，得到小波分解序列其中：表示对{x_i,i＝1,2,…,N}中的第i个数据进行小波分解后的第三层低频分量，和分别表示对{x_i,i＝1,2,…,N}中的第i个数据进行小波分解后的第三层、第二层和第一层高频分量；经过小波分解后，{x_i,i＝1,2,…,N}中的每一个数据由所对应的代替，小波分解后的数据集由{x_i}^Wave表示；

3)、由近邻传播聚类算法将{x_i}^Wave进行分类，再根据{x_i}^Wave所对应的数据的排列序号{index_i,i＝1,2,…,N}将油井油液含水率数据集中的数据划分到各个分类中，形成K个聚类，K≥1；

4)、将每个聚类中的数据按照各个数据的原始的排列序号{index_i,i＝1,2,…,N}进行排列，然后将每个聚类中的数据由如下时间序列形式进行表示：其中，α＝1,2,…,K，t＝1,2,…,M，m为嵌入维数，τ为延迟时间，M＝N-(m-1)τ；将每个聚类的时间序列的输出表示为： $Y_t^{\mathrm{\α}}=x_{t+1+(m-1)\mathrm{\τ}}^{\mathrm{\α}};$

5)、根据极端学***均值最终得到预测值。

2.根据权利要求1所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：步骤2)具体如下：

2.1、首先根据公式(1)对数据集{x_i,i＝1,2,…,N}进行分解：

$\left\{\begin{array}{c}{x}_i^{{\mathrm{low}}_{J+1}}=H\·x_i^{{\mathrm{low}}_J}\\ x_i^{{\mathrm{high}}_{J+1}}=G\·x_i^{{\mathrm{low}}_J}\end{array}\right.---\left(1\right)$

其中，J＝0,1,2；H和G分别为分解低通滤波器和分解高通滤波器；

2.2、然后根据公式(2)进行重构：

$x_i^{{\mathrm{low}}_J}=H^{*}\·x_i^{{\mathrm{low}}_{J+1}}+G^{*}\·x_i^{{\mathrm{high}}_{J+1}}---\left(2\right)$

其中，H^*和G^*分别为H和G的对偶算子；

那么，油井油液含水率数据集{x_i}中的数据x_i经过三层小波分解后表示为：其中：i＝1,2,…,N，由构成的数据集表示为{x_i}^Wave。

3.根据权利要求1所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：步骤3)利用近邻传播聚类算法将{x_i}^Wave进行分类的具体过程为：

3.1、由经过三层小波分解后的数据集{x_i}^Wave计算相似度矩阵S，{x_i}^Wave中的两个数据和(j≠h)的相似度s_jh由公式(3)进行计算：

$s_{jh}=-\left(\sqrt{\underset{l=1}{\overset{4}{\Σ}}{(x_{jl}^{Wave}-x_{hl}^{Wave})}^2}\right)---\left(3\right)$

其中，和分别表示数据和的第l维；

那么，相似度矩阵S表示为：

其中，p为偏向参数，p的初始值取所有初始相似度值的均值；

3.2、设定吸引度矩阵R＝{r_jh}和归属度矩阵A＝{a_jh}，分别表示为：

其中，{r_jh}和{a_jh}的初始值为0；

3.3、设定最大迭代次数MaxLoop，吸引度矩阵R和归属度矩阵A按照如下公式进行迭代：

r_jh＝s_jh-max{a_jh+s_jh}(7)

$a_{jk}=min\{0,r_{kk}+\underset{h=1}{\overset{N}{\Σ}}\{max(0,r_{hk})\left\}\right\}---\left(8\right)$

r_kk＝p_kk-max{a_kh+s_kh}(9)

$a_{kk}=\underset{h=1}{\overset{N}{\Σ}}max\{0,r_{hk}\}---\left(10\right)$

其中，k＝1,2,…,N，k为数据集{x_i}^Wave中数据的下标，表示候选聚类中心，j≠h≠k；

3.4、当满足最大迭代次数MaxLoop后，迭代过程终止，判断两矩阵R和A的对角线元素之和，如果满足r_kk+a_kk＞0，那么选择为聚类中心；

3.5、确定K个聚类的聚类中心后，根据公式(3)分别计算数据集{x_i}^Wave中除外的每个数据与K个聚类中心的距离，将数据集{x_i}^Wave中除外的每个数据归入与它距离最小的聚类中心所属的类中，完成数据集{x_i}^Wave中数据的分类。

4.根据权利要求1所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：步骤5)中利用极端学习机算法建立时间序列模型的原理如下：

设有W个训练样本其中u_q为输入向量，v_q为输出向量，设定包含L个隐含层，激活函数为f(·)以及模型的训练输出表示为Q＝[c₁,c₂,…,c_q]^T，那么ELM模型由如下方程组进行描述：

$\left\{\begin{array}{c}{c}_1=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_{11}f({\mathrm{\ω}}_1,u_1,b_1)\\ c_2=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_{12}f({\mathrm{\ω}}_2,u_2,b_2)\\ \·\\ \·\\ \·\\ c_q=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_{lq}f({\mathrm{\ω}}_l,u_q,b_l)\end{array}\right.---\left(11\right)$

其中，β_lq为第l个隐含层神经元与第q个输出神经元之间的连接权值；ω_l为隐含层神经元与输入神经元之间的连接权值；b_l为第l个隐含层神经元的偏置；

如果训练模型能以零误差逼近W个训练样本，即有那么对于式(11)有下式成立，

$\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\β}}_{lq}f({\mathrm{\ω}}_l,u_q,b_l)=v_q---\left(12\right)$

那么ELM模型的数学描述可以改写成如下矩阵形式：

Hβ＝V(13)

式(13)中，有：

H为隐含层输出矩阵，ω和b在初始化时随机给定，那么，ELM模型的训练可以转化成一个求解非线性方程最小值的问题，即：

$\left|\right|{\mathrm{H\β}}^{*}-V\left|\right|=\underset{\mathrm{\β}}{min}\left|\right|H\mathrm{\β}-V\left|\right|---\left(15\right)$

输出权值矩阵β^*可以由下式求出，

β^*＝H⁺V(16)

其中，H⁺为隐含层输出矩阵H的Moore-Penrose广义逆；

那么，ELM的训练过程可以归纳为如下最优化问题：

$g(\mathrm{\ω},b)=min\underset{q=1}{\overset{W}{\Σ}}\left|\right|v_q-c_q\left|\right|---\left(17\right)$

其中，g(·)表示由ω和b所决定的函数，g(ω,b)表示当ω和b分别取不同值时的函数输出值，ELM训练的目的就是找到最优的β^*，使模型的训练输出值c_q和真实值v_q之间的误差最小；

对于激活函数f(·)的选择，采用高斯函数和Sigmoid函数的集成方式，f(·)的定义如下：

$f(\·)={\mathrm{\λe}}^{-\frac{{(u-1)}^2}{2{\mathrm{\σ}}^2}}+(1-\mathrm{\λ})\frac{1}{1+e^{-u}}---\left(18\right)$

其中，u表示输入向量，λ∈[0,1]为权值，σ²为高斯函数的宽度参数；

基于此，建立时间序列模型的计算步骤如下：初始化，随机生成隐含层输入权值ω，隐含层神经元的偏置b，高斯函数的宽度参数σ²和权值λ；根据公式(14)计算隐含层输出矩阵H；根据公式(16)计算输出权值矩阵β^*；根据公式(17)计算函数输出值。

5.根据权利要求4所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：建立时间序列模型的过程中，采用改进的萤火虫算法对m、τ、ω、b、σ²和λ的取值进行优化选取，将m、τ、ω、b、σ²和λ看作一组解向量，即每个“萤火虫”为1×6维向量，其中6维变量分别为m、τ、ω、b、σ²和λ，计算步骤如下：

第一步，产生初始“萤火虫”群体，表示为：{F₁,F₂,…,F_Npop}，其中：N_pop表示初始“萤火虫”群体中“萤火虫”的数量；

第二步，计算每个“萤火虫”的亮度，表示为：{I₁,I₂,…,I_Npop}；

第三步，每个“萤火虫”向比它亮度更大的“萤火虫”移动，根据公式(22)-公式(24)进行位置更新，

Φ₁(Itenum+1)＝Φ(Itenum)+η(ε)·(Ψ_best(Itenum)-Φ(Itenum))+ξ·(rand-0.5)(22)

Φ₂(Itenum+1)＝Φ(Itenum)+η(ε)·(Ψ_secbest(Itenum)-Φ(Itenum))+ξ·(rand-0.5)(23)

$\mathrm{\Φ}(Itenum+1)=\frac{{\mathrm{\Φ}}_1(Itenum+1)+{\mathrm{\Φ}}_2(Itenum+1)}{2}---\left(24\right)$

其中，ξ表示步长，ξ∈[0,1]；rand为[0,1]上服从均匀分布的一个随机数；Ψ_best为亮度最大的“萤火虫”，Ψ_secbest为亮度大小次于Ψ_best的“萤火虫”；Itenum表示IFA算法运行过程中迭代的次数；Φ(Itenum)表示“萤火虫”在第Itenum次迭代时的位置，Φ(Itenum+1)表示“萤火虫”在第Itenum+1次迭代时的位置，Φ₁(Itenum+1)表示“萤火虫”受到亮度最大的“萤火虫”的吸引在第Itenum+1次迭代时的位置，Φ₂(Itenum+1)表示“萤火虫”受到亮度大小次于Ψ_best的“萤火虫”的吸引在第Itenum+1次迭代时的位置；

第四步，定义适应度函数fitness(m,τ,ω,b,σ²,λ)评价m、τ、ω、b、σ²和λ的取值，当m、τ、ω、b、σ²和λ取最佳值时，fitness(m,τ,ω,b,σ²,λ)函数值最小，定义变量Local_pop记录每一次迭代中m、τ、ω、b、σ²和λ的最佳值，定义变量Global_pop记录所有次迭代中m、τ、ω、b、σ²和λ的最佳值，执行完第三步后，分别更新Local_pop和Global_pop中的值；

第五步，如果达到最大迭代次数，停止迭代，分别输出Local_pop中和Global_pop中的m、τ、ω、b、σ²、λ的值；否则返回步骤第三步重新迭代。

6.根据权利要求5所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：对于K个聚类中的每一个聚类，采用极端学习机算法建立时间序列模型的计算步骤如下：

步骤5.1初始化，设定m、τ、ω、b、σ²和λ的取值范围；

步骤5.2将每一个聚类中的数据归一化到[0,1]区间内，确定待建立模型的输入和输出分别为和

步骤5.3产生初始“萤火虫”群体，在m、τ、ω、b、σ²和λ的取值范围内随机给定m、τ、ω、b、σ²和λ的取值；

步骤5.4根据公式(14)计算隐含层输出矩阵H，根据公式(16)计算输出权值矩阵β^*，根据公式(17)计算函数输出值；

步骤5.5计算适应度函数fitness(m,τ,ω,b,σ²,λ)的值，记录Local_pop和Global_pop中的m、τ、ω、b、σ²和λ的初始值；

步骤5.6计算每个“萤火虫”的亮度，根据公式(22)-公式(24)进行位置更新；

步骤5.7重复步骤5.4，计算适应度函数fitness(m,τ,ω,b,σ²,λ)的值，更新Local_pop和Global_pop中的m、τ、ω、b、σ²、λ的值；

步骤5.8如果达到最大迭代次数，停止迭代，分别输出Local_pop中和Global_pop中的m、τ、ω、b、σ²、λ的值；否则返回步骤5.6-步骤5.7重新迭代；

步骤5.9将每一个聚类得到的最好的m、τ、ω、b、σ²和λ的值作为极端学习机算法建立时间序列模型的参数，由输入的值计算输出的值，并进行反归一化处理。

7.根据权利要求1所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：采用单步预测，即每次预测出一个再预测下一个输出时，将加入到下一个输入的最后；每次预测出一个Λ_t，再预测下一个输出Λ_t+1时，将Λ_t加入到下一个输入Γ_t+1的最后。

8.根据权利要求1所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：所述油井油液含水率数据集中数据的数量N为300-800个。

9.根据权利要求6所述的基于时间序列的油井油液含水率多模型预测方法，其特征在于：m、τ、ω、b、σ²和λ的取值范围分别为：m∈[1,50]，τ∈[1,10]，ω∈[0,1]，b∈[0,10]，σ²∈[0.01,1000]，λ∈[0,1]。