CN107909202B - 一种基于时间序列的油井产液量集成预测方法 - Google Patents

Abstract

本发明属于石油生产技术领域，具体提供了一种基于时间序列的油井产液量集成预测方法，利用历史数据建立油井产液量数据集，采用经验模态分解方法对油井产液量数据集{X_i}中的数据进行处理；原数据集分解成k个IMF分量数据集和一个余项数据集；采用集成预测方法，分别由基于ELM方法和基于ESN方法以并行的方式同时对k个IMF分量数据集和一个余项数据集进行预测；将每个序列由两种方法所得到的预测值求平均值，然后将所有序列的预测值求和得到最终结果。本发明提供的这种预测方法，由过去一段时间的数据去预测未来某个时间点的值，原理简单，计算复杂性低，准确度高，可以有效降低异常数据的影响。

Description

一种基于时间序列的油井产液量集成预测方法

技术领域

本发明涉及石油生产技术领域，具体涉及一种基于时间序列的油井产液量集成预测方法。

背景技术

油井产液量是油田生产的一个重要指标，不仅关系到油井的生产寿命，还与企业的经济效益密切相关。对其进行准确的预测，对合理评价油井生产状况、及时调整抽汲参数和开采方案等都具有重要的意义。油井采油过程动态变化，会受到流体特性、油藏条件、人力干扰等的影响，很难在机理上掌握产油量的变化规律。目前的油田生产只能从井下采出油液后判断油井的实际生产状况，无法提前获知产液量的变化趋势。这使得油井开采方案的制定或调整总是滞后于油井的实际生产状况，影响企业生产效益。

发明内容

本发明的目的是克服现有技术的缺陷，提供一种基于时间序列的油井产液量集成预测方法，由已有的过去一段时期的油井产液量数据预测下一个时间点的产液量数据。

本发明的技术方案如下：

一种基于时间序列的油井产液量集成预测方法，包括如下步骤：

1)利用历史数据建立油井产液量数据集为{x_i，i＝1，2，…，n}，n为数据点的个数；

2)采用经验模态分解(EMD，Empirical Mode Decomposition)方法即EMD方法对油井产液量数据集{x_i}中的数据进行处理；

3)将{x_i}中所有极大值点用三次样条插值函数拟合，得到{x_i}的上包络线；将{x_i}中所有极小值点用三次样条插值函数拟合，得到{x_i}的下包络线；{x_i}的上下包络线的平均值记作

将原数据集{x_i}减去该平均值

得到一个新的数据集{h₁(i)}，记作：

4)若{h₁(i)}满足条件，则将该{h₁(i)}看作{x_i}的第一个IMF分量，记为c₁(i)＝h₁(i)；若{h₁(i)}不满足条件，则使{h₁(i)}取代原数据集{x_i}，重复步骤3)直至计算得到一个满足条件的IMF分量，记为{c₁(i)}；

5)用{x_i}减去{c₁(i)}得到剩余值数据集{x₁(i)}，有：x₁(i)＝x_i-c₁(i)；把{x₁(i)}作为一个新的待分解数据集，重复步骤3)和步骤4)，依次提取出第2个、第3个直至第k个IMF分量，以及原数据集的余项{r_res(i)}；当满足终止条件时结束分解，终止条件为最新的数据集不能再提取IMF分量；

原数据集{x_i}分解成k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}；

6)将k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}中的数据进行相空间重构，得到时间序列模型，输入表示为：X_t＝[x_t，x_t+τ，…，x_t+(m-1)τ]，其中，t＝1，2，…，M，m为嵌入维数，τ为延迟时间，M＝n-(m-1)τ，n为数据集中数据的个数；将每个时间序列的输出表示为：Y_t＝x_t+1+(m-1)τ；

7)采用集成预测方法，分别由基于ELM方法(极限学习机的方法)和基于ESN方法(回声状态网络)以并行的方式同时对k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}进行预测；

8)对于k个IMF分量数据集时间序列和一个余项数据集时间序列，首先将每个序列由两种方法所得到的预测值求平均值，然后将所有序列的预测值求和得到最终结果；计算公式如下：

其中，

表示第i个IMF分量数据集时间序列由ELM方法得到的预测值；

表示第i个IMF分量数据集时间序列由ESN方法得到的预测值；

表示余项数据集时间序列由ELM方法得到的预测值；

表示余项数据集时间序列由ESN方法得到的预测值。

所述步骤7)中根据极限学习机(ELM)方法建立k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}的时间序列模型；每个模型的输入和输出分别为X_t和Y_t，ELM方法的基本原理如下：

设有W个训练样本

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

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

如果训练模型能以零误差逼近W个训练样本，即有

那么对于式(1)有下式成立，

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

Hβ＝V (3)

式(3)中，有：

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

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

β^*＝H⁺V (6)

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

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

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

对于激活函数f(·)的选择，本发明专利采用高斯函数，定义如下：

其中，z表示输入向量，σ²为高斯函数的宽度参数。

进一步地，所述基于极限学习机方法建立时间序列模型的计算步骤如下：

初始化，随机生成隐含层输入权值ω，隐含层神经元的偏置b，高斯函数的宽度参数σ²；

根据公式(4)计算隐含层输出矩阵H；

根据公式(6)计算输出权值矩阵β^*；

根据公式(7)计算函数输出值。

进一步地，由ELM方法建立时间序列模型的过程中，m、τ、ω、b和σ²的取值决定了ELM模型的计算精度，采用一种改进的果蝇优化算法(IFOA，Improved Fruit OptimizationAlgorithm)对m、τ、ω、b和σ²的取值进行优化选取；IFOA算法的数学描述如下：

(1)设定果蝇群体的位置区间[p_low，p_up]，在位置区间内随机给定每个果蝇的初始位置Axis_x_i和Axis_y_i如下：

Axis_x_i＝p_low+rand₁·(p_up-p_low) (9)

Axis_y_i＝p_low+rand₂·(p_up-p_low) (10)

其中，i＝1，2，...，sizepop，sizepop为果蝇群体中果蝇的个数；rand₁和rand₂分别为[0，1]区间内的随机数。

(2)果蝇个体利用嗅觉搜寻食物，每个个体的位置更新公式如下：

其中，Axis_X_i和Axis_Y_i表示第i个果蝇个体移动后的位置坐标；BestAxis_x和BestAxis_y表示当前最优个体的位置坐标；Ite和Ite_max分别表示当前迭代次数和最大迭代次数；rand₃和rand₄分别为[-1，1]区间内的随机数；rand₅为[0，1]区间内的随机数。

(3)计算果蝇个体与原点的距离，如下：

(4)计算味道浓度判定值，如下：

(5)由果蝇个体的味道浓度判定值计算每个果蝇个体的味道浓度值，如下：

Smell_t＝fitness(S_t) (15)

其中，fitness为适应度函数，用以评价待优化参数的取值。

(6)找到果蝇群体中的最优个体，如下：

[bestSmell bestIndex]＝min/max(Smell) (16)

其中，bestSmell表示最优的味道浓度值，bestIndex表示最优味道浓度值对应的果蝇个体的序号；min/max表示最小或最大。

进一步地，IFOA的计算步骤如下：

步骤(1)、产生初始的与待优化参数数量相同的果蝇群体，根据待优化参数的取值范围设定每个果蝇群体的位置区间[p_low，p_up]，根据公式(9)和(10)得到每个果蝇群体中每个果蝇的初始位置；

步骤(2)、根据公式(13)和(14)分别计算每个果蝇群体中每个果蝇个体的味道浓度判定值；

步骤(3)、根据公式(15)分别计算每个果蝇群体中每个果蝇个体的味道浓度值；

步骤(4)、根据公式(16)找到每个果蝇群体中的最优果蝇个体，并记录它们的位置[BestAxis_x，BestAxis_y]，定义变量Pop_local记录每一次迭代中每个果蝇群体中的最优果蝇个体的位置，定义变量Pop_global记录所有次迭代中每个果蝇群体中的最优果蝇个体的位置；

步骤(5)、根据公式(11)和(12)更新每个果蝇群体中每个果蝇的位置，重复步骤(2)到步骤(4)；

步骤(6)、如果达到最大迭代次数，停止迭代，分别输出Pop_local和Pop_global中的果蝇个体的位置，并输出其味道浓度判定值；否则返回步骤(5)重新迭代。

进一步地，对于分别建立的k个IMF分量数据集时间序列模型和一个余项数据集时间序列模型，在采用ELM方法建立模型的同时，由IFOA方法对ELM计算过程中所涉及的m、τ、ω、b和σ²进行优化选取，计算步骤如下：

步骤7.1.1、初始化，产生初始5个果蝇群体，分别指派给m、τ、ω、b和σ²，设定m、τ、ω、b和σ²的取值范围；

步骤7.1.2、将k个IMF分量数据集和一个余项数据集中的数据分别归一化到[0，1]区间内，确定每一个待建立模型的输入和输出分别为X_t和Y_t；

步骤7.1.3、对于每一个待建立模型，由IFOA的计算步骤(1)-(6)得到最优的味道浓度判定值，分别对应m、τ、ω、b和σ²；

步骤7.1.4、对于每一个待建立模型，将得到的最优参数值分别代入相空间重构后的时间序列模型和ELM模型，由输入X_t的值计算输出Y_t的值并进行反归一化处理。

进一步地，根据回声状态网络(ESN)方法分别建立k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}的时间序列模型；每个模型的输入和输出分别为X_t和Y_t，ESN方法的基本原理如下：

ESN模型可以表示为：

其中，u(t)、x(t)和y(t)分别为输入单元、内部状态单元和输出单元在时刻t的值；M_esn、N_esn和K_esn分别为输入单元、内部状态单元和输出单元的数量；

ESN模型的学习方式可以由如下公式进行描述：

其中，x^esn(t+1)和y^esn(t+1)分别为内部状态单元和输出单元在时刻t+1的值；f(·)和f_out(·)分别为内部单元和输出单元的激活函数；

为输入层到内部单元的连接权值矩阵；W^esn为内部单元的连接权值矩阵；

为输出层反馈到内部单元的连接权值矩阵；

为内部单元到输出层的连接权值矩阵。

由如下公式进行计算：

其中，y_esn(t)为模型期望输出，y′_esn(t)为模型实际输出；w_out为权值，是

的元素。

ESN模型的求解可以转化为如下的优化求解问题：

其中，P_esn和m_esn分别表示时间点。

可以归结为如下矩阵求逆问题：

其中，B_esn为

构成的矩阵；L_esn为

构成的列矩阵。

进一步地，由ESN方法建立时间序列模型的过程中，ESN模型中有6个参数决定了模型的计算精度，分别为：m、τ、连接权值矩阵W^esn的绝对值最大的特征值SR_esn、N_esn、输入与内部单元相互连接的尺度因子IS_esn和内部单元中相互连接的单元数与总单元数的比值SD_esn；

采用IFOA方法对m、τ、SR_esn、N_esn、IS_esn和SD_esn的取值进行优化选取；对于分别建立的k个IMF分量数据集时间序列模型和一个余项数据集时间序列模型，在采用ESN方法建立模型的同时，由IFOA方法对ESN计算过程中所涉及的m、τ、SR_esn、N_esn、IS_esn和SD_esn进行优化选取，计算步骤如下：

步骤7.2.1、初始化，产生初始6个果蝇群体，分别指派给m、τ、SR_esn、N_esn、IS_esn和SD_esn，设定m、τ、SR_esn、N_esn、IS_esn和SD_esn的取值范围；

步骤7.2.2、将k个IMF分量数据集和一个余项数据集中的数据分别归一化到[0，1]区间内，确定每一个待建立模型的输入和输出分别为X_t和Y_t；

步骤7.2.3、对于每一个待建立模型，由由IFOA的计算步骤得到最优的味道浓度判定值，分别对应m、τ、SR_esn、N_esn、IS_esn和SD_esn；

步骤7.2.4、对于每一个待建立模型，将得到的最优参数值分别代入相空间重构后的时间序列模型和ESN模型，由输入X_t的值计算输出Y_t的值并进行反归一化处理。

进一步地，所采用的集成预测方法，其中所涉及的两种方法为并行计算方式，当预测开始时，两种方法同时启动运行。

进一步地，为了保证预测的有效性和降低预测的复杂性，当时间点单位是小时时，油井产液量数据集中的数据个数一般取500-1000个；当时间点单位是日时，油井产液量数据集中的数据个数一般取30-100个。

进一步地，隐含层数量L的取值一般为30-60个。

进一步地，IFOA算法中的最大迭代次数为300次。

进一步地，IFOA算法中的果蝇群体中果蝇的个数为20-30个。

进一步地，m、τ、ω、b、σ²、SR_esn、N_esn、IS_esn和SD_esn的取值范围分别为：m∈[1，30]，τ∈[1，10]，ω∈[0，1]，b∈[0，10]，σ²∈[0.01，1000]，SR_esn∈[0.1，1)，N_esn∈[10，150)，IS_esn∈[0.01，1)，SD_esn∈[0.01，1)。

本发明的有益效果：本发明提供的这种基于时间序列的油井产液量集成预测方法，由过去一段时间的数据去预测未来某个时间点的值，从数据本身的特性找到其变化趋势，可以消除流体特性变化、油藏条件变化、人力干扰等对机理分析的影响，对油田企业进行高效、稳定的生产具有一定的实际意义；原理简单，计算复杂性低，准确度高，将原始数据集进行分解得到多个子数据集，可以有效降低异常数据的影响；采用两种方法进行预测，将预测值取平均值，可以有效提高预测的稳定性；采用群智能优化方法合理选取模型参数，可以有效降低模型参数不确定性对预测结果的影响。

附图说明

图1是本发明实施例所提供的油井产液量数据集示意图；

图2是本发明实施例所提供采用EMD方法分解后的8个IMF分量数据集和1个余项数据集的示意图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

实施例：

本实施例提供了一种基于时间序列的油井产液量集成预测方法，包括如下步骤：

1)利用历史数据建立油井产液量数据集，数据集为{x_i，i＝1，2，...，831}，时间点单位为小时，数据集{x_i}如图1所示；

2)采用经验模态分解(EMD，Empirical Mode Decomposition)方法即EMD方法对油井产液量数据集{x_i}中的数据进行处理；

3)将{x_i}中所有极大值点用三次样条插值函数拟合，得到{x_i}的上包络线；将{x_i}中所有极小值点用三次样条插值函数拟合，得到{x_i}的下包络线；{x_i}的上下包络线的平均值记作

将原数据集{x_i}减去该平均值

得到一个新的数据集{h₁(i)}，记作：

4)若{h₁(i)}满足条件，则将该{h₁(i)}看作{x_i}的第一个IMF分量，记为c₁(i)＝h₁(i)；若{h₁(i)}不满足条件，说明它还不是一个本征模函数，则使{h₁(i)}取代原数据集{x_i}，重复步骤3)直至计算得到一个满足条件的IMF分量，记为{c₁(i)}；

5)用{x_i}减去{c₁(i)}得到剩余值数据集{x₁(i)}，有：x₁(i)＝x_i-c₁(i)；把{x₁(i)}作为一个新的待分解数据集，重复步骤3)和步骤4)，依次提取出第2个、第3个直至第k个IMF分量，以及原数据集的余项{r_res(i)}；当满足终止条件时结束分解，终止条件为最新的数据集不能再提取IMF分量；

原数据集{x_i}分解成k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}；

经过EMD分解，油井产液量数据集{x_i}分解成8个IMF分量数据集和1个余项数据集，如图2所示。

6)将8个IMF分量数据集{c_k(i)}和1个余项数据集{r_res(i)}中的数据进行相空间重构，得到时间序列模型，输入表示为：X_t＝[x_t，x_t+τ，…，x_t+(m-1)τ]，其中，t＝1，2，…，M，m为嵌入维数，τ为延迟时间，M＝n-(m-1)τ，n为数据集中数据的个数；将每个时间序列的输出表示为：Y_t＝x_t+1+(m-1)τ；

7)采用集成预测方法，分别由基于ELM方法(极限学习机的方法)和基于ESN方法(回声状态网络)以并行的方式同时对8个IMF分量数据集{c_k(i)}和1个余项数据集{r_res(i)}进行预测；

根据极限学习机(ELM)方法建立8个IMF分量数据集{c_k(i)}和1个余项数据集{r_res(i)}的时间序列模型；每个模型的输入和输出分别为X_t和Y_t，ELM方法的基本原理如下：

设有W个训练样本

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

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

如果训练模型能以零误差逼近W个训练样本，即有

那么对于式(1)有下式成立，

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

Hβ＝V (3)

式(3)中，有：

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

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

β^*＝H⁺V (6)

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

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

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

对于激活函数f(·)的选择，本发明专利采用高斯函数，定义如下：

其中，z表示输入向量，σ²为高斯函数的宽度参数。

进一步地，所述基于极限学习机方法建立时间序列模型的计算步骤如下：

初始化，随机生成隐含层输入权值ω，隐含层神经元的偏置b，高斯函数的宽度参数σ²；

根据公式(4)计算隐含层输出矩阵H；

根据公式(6)计算输出权值矩阵β^*；

根据公式(7)计算函数输出值。

在本实施例中，由ELM方法建立时间序列模型的过程中，m、τ、ω、b和σ²的取值决定了ELM模型的计算精度，本实施例采用一种改进的果蝇优化算法(IFOA，Improved FruitOptimizationAlgorithm)对m、τ、ω、b和σ²的取值进行优化选取；IFOA算法的数学描述如下：

(1)设定果蝇群体的位置区间[p_low，p_up]，在位置区间内随机给定每个果蝇的初始位置Axis_x_i和Axis_y_i如下：

Axis_x_i＝p_low+rand₁·(p_up-p_low) (9)

Axis_y_i＝p_low+rand₂·(p_up-p_low) (10)

其中，i＝1，2，...，sizepop，sizepop为果蝇群体中果蝇的个数；rand₁和rand₂分别为[0，1]区间内的随机数。

(2)果蝇个体利用嗅觉搜寻食物，每个个体的位置更新公式如下：

其中，Axis_X_i和Axis_Y_i表示第i个果蝇个体移动后的位置坐标；BestAxis_x和BestAxis_y表示当前最优个体的位置坐标；Ite和Ite_max分别表示当前迭代次数和最大迭代次数；rand₃和rand₄分别为[-1，1]区间内的随机数；rand₅为[0，1]区间内的随机数。

(3)计算果蝇个体与原点的距离，如下：

(4)计算味道浓度判定值，如下：

(5)由果蝇个体的味道浓度判定值计算每个果蝇个体的味道浓度值，如下：

Smell_i＝fitness(S_i) (15)

其中，fitness为适应度函数，用以评价待优化参数的取值。

(6)找到果蝇群体中的最优个体，如下：

[bestSmell bestIndex]＝min/max(Smell) (16)

其中，bestSmell表示最优的味道浓度值，bestIndex表示最优味道浓度值对应的果蝇个体的序号；min/max表示最小或最大。

IFOA的计算步骤如下：

步骤(1)、产生初始的与待优化参数数量相同的果蝇群体，根据待优化参数的取值范围

设定每个果蝇群体的位置区间[p_low，p_up]，根据公式(9)和(10)得到每个果蝇群体中每个果蝇的初始位置；

步骤(2)、根据公式(13)和(14)分别计算每个果蝇群体中每个果蝇个体的味道浓度判定值；

步骤(3)、根据公式(15)分别计算每个果蝇群体中每个果蝇个体的味道浓度值；

步骤(4)、根据公式(16)找到每个果蝇群体中的最优果蝇个体，并记录它们的位置[BestAxis_x，BestAxis_y]，定义变量Pop_local记录每一次迭代中每个果蝇群体中的最优果蝇个体的位置，定义变量Pop_global记录所有次迭代中每个果蝇群体中的最优果蝇个体的位置；

步骤(5)、根据公式(11)和(12)更新每个果蝇群体中每个果蝇的位置，重复步骤(2)到步骤(4)；

步骤(6)、如果达到最大迭代次数，停止迭代，分别输出Pop_local和Pop_global中的果蝇个体的位置，并输出其味道浓度判定值；否则返回步骤(5)重新迭代。

进一步地，对于分别建立的8个IMF分量数据集时间序列模型和1个余项数据集时间序列模型，在采用ELM方法建立模型的同时，由IFOA方法对ELM计算过程中所涉及的m、τ、ω、b和σ²进行优化选取，计算步骤如下：

步骤7.1.1、初始化，产生初始5个果蝇群体，分别指派给m、τ、ω、b和σ²，设定m、τ、ω、b和σ²的取值范围；m∈[1，30]，τ∈[1，10]，ω∈[0，1]，b∈[0，10]，σ²∈[0.01，1000]

步骤7.1.2、将8个IMF分量数据集和1个余项数据集中的数据分别归一化到[0，1]区间内，确定每一个待建立模型的输入和输出分别为X_t和Y_t；

步骤7.1.3、对于每一个待建立模型，由IFOA的计算步骤(1)-(6)得到最优的味道浓度判定值，分别对应m、τ、ω、b和σ²；

步骤7.1.4、对于每一个待建立模型，将得到的最优参数值分别代入相空间重构后的时间序列模型和ELM模型，由输入X_t的值计算输出Y_t的值并进行反归一化处理。

在本实施例中，根据回声状态网络(ESN)方法分别建立8个IMF分量数据集{c_k(i)}和1个余项数据集{r_res(i)}的时间序列模型；每个模型的输入和输出分别为X_t和Y_t，ESN方法的基本原理如下：

ESN模型可以表示为：

其中，u(t)、x(t)和y(t)分别为输入单元、内部状态单元和输出单元在时刻t的值；M_esn、N_esn和K_esn分别为输入单元、内部状态单元和输出单元的数量；

ESN模型的学习方式可以由如下公式进行描述：

其中，x^esn(t+1)和y^esn(t+1)分别为内部状态单元和输出单元在时刻t+1的值；f(·)和f_out(·)分别为内部单元和输出单元的激活函数；

为输入层到内部单元的连接权值矩阵；W^esn为内部单元的连接权值矩阵；

为输出层反馈到内部单元的连接权值矩阵；

为内部单元到输出层的连接权值矩阵。

由如下公式进行计算：

其中，y_esn(t)为模型期望输出，y′_esn(t)为模型实际输出；w_out为权值，是

的元素。

ESN模型的求解可以转化为如下的优化求解问题：

其中，P_esn和m_esn分别表示时间点。

可以归结为如下矩阵求逆问题：

其中，B_esn为

构成的矩阵；L_esn为

构成的列矩阵。

在本实施例中，由ESN方法建立时间序列模型的过程中，ESN模型中有6个参数决定了模型的计算精度，分别为：m、τ、连接权值矩阵W^esn的绝对值最大的特征值SR_esn、N_esn、输入与内部单元相互连接的尺度因子IS_esn和内部单元中相互连接的单元数与总单元数的比值SD_esn；

采用IFOA方法对m、τ、SR_esn、N_esn、IS_esn和SD_esn的取值进行优化选取；对于分别建立的8个IMF分量数据集时间序列模型和1个余项数据集时间序列模型，在采用ESN方法建立模型的同时，由IFOA方法对ESN计算过程中所涉及的m、τ、SR_esn、N_esn、IS_esn和SD_esn进行优化选取，计算步骤如下：

步骤7.2.1、初始化，产生初始6个果蝇群体，分别指派给m、τ、SR_esn、N_esn、IS_esn和SD_esn，设定m、τ、SR_esn、N_esn、IS_esn和SD_esn的取值范围；m∈[1，30]，τ∈[1，10]，SR_esn∈[0.1，1)，N_esn∈[10，150)，IS_esn∈[0.01，1)，SD_esn∈[0.01，1)；

步骤7.2.2、将8个IMF分量数据集和1个余项数据集中的数据分别归一化到[0，1]区间内，确定每一个待建立模型的输入和输出分别为X_t和Y_t；

步骤7.2.3、对于每一个待建立模型，由IFOA的计算步骤得到最优的味道浓度判定值，分别对应m、τ、SR_esn、N_esn、IS_esn和SD_esn；

步骤7.2.4、对于每一个待建立模型，将得到的最优参数值分别代入相空间重构后的时间序列模型和ESN模型，由输入X_t的值计算输出Y_t的值并进行反归一化处理。

8)对于8个IMF分量数据集时间序列和1个余项数据集时间序列，首先将每个序列由两种方法所得到的预测值求平均值，然后将所有序列的预测值求和得到最终结果；计算公式如下：

其中，

表示第i个IMF分量数据集时间序列由ELM方法得到的预测值；

表示第i个IMF分量数据集时间序列由ESN方法得到的预测值；

表示余项数据集时间序列由ELM方法得到的预测值；

表示余项数据集时间序列由ESN方法得到的预测值。

本实施例中，由ELM方法得到的预测结果分别为：0.2182、0.1127、-0.0138、0.789、-0.0205、0.0227、0.0099、0.0082和0.0041；由ESN方法得到的预测结果分别为：0.2455、0.1320、0.0050、0.6685、-0.0534、-0.0031、0.1002、-0.0011和0.0095。

由公式(22)计算得到最终预测结果为：1.1168。

以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (8)

1.一种基于时间序列的油井产液量集成预测方法，其特征在于，包括如下步骤：

1)利用历史数据建立油井产液量数据集为{x_i，i＝1，2，…，n}，n为数据点的个数；

2)采用经验模态分解方法对油井产液量数据集{x_i，i＝1，2，…，n}中的数据进行处理；

3)将{x_i，i＝1，2，…，n}中所有极大值点用三次样条插值函数拟合，得到{x_i，i＝1，2，…，n}的上包络线；将{x_i，i＝1，2，…，n}中所有极小值点用三次样条插值函数拟合，得到{x_i，i＝1，2，…，n}的下包络线；{x_i，i＝1，2，…，n}的上下包络线的平均值记作

将原数据集{x_i，i＝1，2，…，n}减去该平均值

得到一个新的数据集{h₁(i)}，记作：

4)若{h₁(i)}满足条件，则将该{h₁(i)}看作{x_i，i＝1，2，…，n}的第一个IMF分量，记为c₁(i)＝h₁(i)；若{h₁(i)}不满足条件，则使{h₁(i)}取代原数据集{x_i，i＝1，2，…，n}，重复步骤3)直至计算得到一个满足条件的IMF分量，记为{c₁(i)}；

5)用{x_i，i＝1，2，…，n}减去{c₁(i)}得到剩余值数据集{x₁(i)}，有：x₁(i)＝x_i-c₁(i)；把{x₁(i)}作为一个新的待分解数据集，重复步骤3)和步骤4)，依次提取出第2个、第3个直至第k个IMF分量，以及原数据集的余项{r_res(i)}；当满足终止条件时结束分解，终止条件为最新的数据集不能再提取IMF分量；

原数据集{x_i，i＝1，2，…，n}分解成k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}；

6)将k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}中的数据进行相空间重构，得到时间序列模型，输入表示为：X_t＝[x_t，x_t+τ，…，x_t+(m-1)τ]，其中，t＝1，2，…，M，m为嵌入维数，τ为延迟时间，M＝n-(m-1)τ，n为数据集中数据的个数；将每个时间序列的输出表示为：Y_t＝x_t+1+(m-1)τ；

7)采用集成预测方法，分别由基于ELM方法和基于ESN方法以并行的方式同时对k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}进行预测；

根据ELM方法建立k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}的时间序列模型；每个模型的输入和输出分别为X_t和Y_t，ELM方法的基本原理如下：

设有W个训练样本

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

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

如果训练模型能以零误差逼近W个训练样本，即有

那么对于式(1)有下式成立，

那么ELM模型可以改写成如下矩阵形式：

Hβ＝V (3)

式(3)中，有：

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

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

β^*＝H⁺V (6)

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

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

其中，G(·)表示由ω和b所决定的函数，G(ω，b)表示当ω和b分别取不同值时的函数输出值；

对于激活函数f(·)的选择，采用高斯函数，定义如下：

其中，z表示输入向量，σ²为高斯函数的宽度参数；

根据ESN方法分别建立k个IMF分量数据集{c_k(i)}和一个余项数据集{r_res(i)}的时间序列模型；每个模型的输入和输出分别为X_t和Y_t，ESN方法如下：

ESN模型可以表示为：

其中，u(t)、x(t)和y(t)分别为输入单元、内部状态单元和输出单元在时刻t的值；M_esn、N_esn和K_esn分别为输入单元、内部状态单元和输出单元的数量；

ESN模型的学习方式可以由如下公式进行描述：

其中，x^esn(t+1)和y^esn(t+1)分别为内部状态单元和输出单元在时刻t+1的值；f(·)和f_out(·)分别为内部单元和输出单元的激活函数；

为输入层到内部单元的连接权值矩阵；W^esn为内部单元的连接权值矩阵；

为输出层反馈到内部单元的连接权值矩阵；

为内部单元到输出层的连接权值矩阵；

由如下公式进行计算：

其中，y_esn(t)为模型期望输出，y′_esn(t)为模型实际输出；w_out为权值，是

的元素；

ESN模型的求解可以转化为如下的优化求解问题：

其中，P_esn和m_esn分别表示时间点；

可以归结为如下矩阵求逆问题：

其中，B_esn为

构成的矩阵；L_esn为

构成的列矩阵；

8)对于k个IMF分量数据集时间序列和一个余项数据集时间序列，首先将每个序列由两种方法所得到的预测值求平均值，然后将所有序列的预测值求和得到最终结果；计算公式如下：

其中，

表示第i个IMF分量数据集时间序列由ELM方法得到的预测值；

表示第i个IMF分量数据集时间序列由ESN方法得到的预测值；

表示余项数据集时间序列由ELM方法得到的预测值；

表示余项数据集时间序列由ESN方法得到的预测值。

2.如权利要求1所述的方法，其特征在于：所述基于ELM方法建立时间序列模型的计算步骤如下：

初始化，随机生成隐含层输入权值ω，隐含层神经元的偏置b，高斯函数的宽度参数σ²；

根据公式(4)计算隐含层输出矩阵H；

根据公式(6)计算输出权值矩阵β^*；

根据公式(7)计算函数输出值。

3.如权利要求1所述的方法，其特征在于：由ELM方法建立时间序列模型的过程中，采用一种改进的果蝇优化算法对m、τ、ω、b和σ²的取值进行优化选取；改进的果蝇优化算法如下：

(1)设定果蝇群体的位置区间[p_low，p_up]，在位置区间内随机给定每个果蝇的初始位置Axis_x_i和Axis_y_i如下：

Axis_x_i＝p_low+rand₁·(p_up-p_low) (9)

Axis_y_i＝p_low+rand₂·(p_up-p_low) (10)

其中，i＝1，2，...，sizepop，sizepop为果蝇群体中果蝇的个数；rand₁和rand₂分别为[0，1]区间内的随机数；

(2)果蝇个体利用嗅觉搜寻食物，每个个体的位置更新公式如下：

其中，Axis_X_i和Axis_Y_i表示第i个果蝇个体移动后的位置坐标；BestAxis_x和BestAxis_y表示当前最优个体的位置坐标；Ite和Ite_max分别表示当前迭代次数和最大迭代次数；rand₃和rand₄分别为[-1，1]区间内的随机数；rand₅为[0，1]区间内的随机数；

(3)计算果蝇个体与原点的距离，如下：

(4)计算味道浓度判定值，如下：

(5)由果蝇个体的味道浓度判定值计算每个果蝇个体的味道浓度值，如下：

Smell_i＝fitness(S_i) (15)

其中，fitness为适应度函数，用以评价待优化参数的取值；

(6)找到果蝇群体中的最优个体，如下：

[bestSmell bestIndex]＝min/max(Smell) (16)

其中，bestSmell表示最优的味道浓度值，bestlndex表示最优味道浓度值对应的果蝇个体的序号；min/max表示最小或最大。

4.如权利要求3所述的方法，其特征在于：改进的果蝇优化算法的计算步骤如下：

步骤(1)、产生初始的与待优化参数数量相同的果蝇群体，根据待优化参数的取值范围设定每个果蝇群体的位置区间[p_low，p_up]，根据公式(9)和(10)得到每个果蝇群体中每个果蝇的初始位置；

步骤(2)、根据公式(13)和(14)分别计算每个果蝇群体中每个果蝇个体的味道浓度判定值；

步骤(3)、根据公式(15)分别计算每个果蝇群体中每个果蝇个体的味道浓度值；

步骤(4)、根据公式(16)找到每个果蝇群体中的最优果蝇个体，并记录它们的位置[BestAxis_x，BestAxis_y]，定义变量Pop_local记录每一次迭代中每个果蝇群体中的最优果蝇个体的位置，定义变量Pop_global记录所有次迭代中每个果蝇群体中的最优果蝇个体的位置；

步骤(5)、根据公式(11)和(12)更新每个果蝇群体中每个果蝇的位置，重复步骤(2)到步骤(4)；

步骤(6)、如果达到最大迭代次数，停止迭代，分别输出Pop_local和Pop_global中的果蝇个体的位置，并输出其味道浓度判定值；否则返回步骤(5)重新迭代。

5.如权利要求3所述的方法，其特征在于：对于分别建立的k个IMF分量数据集时间序列模型和一个余项数据集时间序列模型，在采用ELM方法建立模型的同时，由改进的果蝇优化算法对ELM计算过程中所涉及的m、τ、ω、b和σ²进行优化选取，计算步骤如下：

初始化，产生初始5个果蝇群体，分别指派给m、τ、ω、b和σ²，设定m、τ、ω、b和σ²的取值范围；

将k个IMF分量数据集和一个余项数据集中的数据分别归一化到[0，1]区间内，确定每一个待建立模型的输入和输出分别为X_t和Y_t；

对于每一个待建立模型，由改进的果蝇优化算法的计算步骤(1)-(6)得到最优的味道浓度判定值，分别对应m、τ、ω、b和σ²；

对于每一个待建立模型，将得到的最优参数值分别代入相空间重构后的时间序列模型和ELM模型，由输入X_t的值计算输出Y_t的值并进行反归一化处理。

6.如权利要求3所述的方法，其特征在于：由ESN方法建立时间序列模型的过程中，ESN模型中有6个参数决定了模型的计算精度，分别为：m、τ、连接权值矩阵W^esn的绝对值最大的特征值SR_esn、N_esn、输入单元与内部状态单元相互连接的尺度因子IS_esn和内部状态单元中相互连接的单元数与总单元数的比值SD_esn；

采用改进的果蝇优化算法对m、τ、SR_esn、N_esn、IS_esn和SD_esn的取值进行优化选取；对于分别建立的k个IMF分量数据集时间序列模型和一个余项数据集时间序列模型，在采用ESN方法建立模型的同时，由改进的果蝇优化算法对ESN计算过程中所涉及的m、τ、SR_esn、N_esn、IS_esn和SD_esn进行优化选取，计算步骤如下：

初始化，产生初始6个果蝇群体，分别指派给m、τ、SR_esn、N_esn、IS_esn和SD_esn，设定m、τ、SR_esn、N_esn、IS_esn和SD_esn的取值范围；

将k个IMF分量数据集和一个余项数据集中的数据分别归一化到[0，1]区间内，确定每一个待建立模型的输入和输出分别为X_t和Y_t；

对于每一个待建立模型，由改进的果蝇优化算法的计算步骤得到最优的味道浓度判定值，分别对应m、τ、SR_esn、N_esn、IS_esn和SD_esn；

对于每一个待建立模型，将得到的最优参数值分别代入相空间重构后的时间序列模型和ESN模型，由输入X_t的值计算输出Y_t的值并进行反归一化处理。

7.如权利要求6所述的方法，其特征在于：m、τ、ω、b、σ²、SR_esn、N_esn、IS_esn和SD_esn的取值范围分别为：m∈[1，30]，τ∈[1，10]，ω∈[0，1]，b∈[0，10]，σ²∈[0.01，1000]，SR_esn∈[0.1，1)，N_esn∈[10，150)，IS_esn∈[0.01，1)，SD_esn∈[0.01，1)。

8.如权利要求1所述的方法，其特征在于：当时间点单位是小时时，油井产液量数据集中的数据个数取500-1000个；当时间点单位是日时，油井产液量数据集中的数据个数取30-100个。

Images