CN109508753A - 一种矿物浮选过程指标的在线预测方法 - Google Patents

Abstract

本发明公开一种矿物浮选过程指标的在线预测方法，首先基于关联K近邻(KNN)对样本缺失数据进行插补，再利用皮尔逊相关系数和显著性检验的方法进行特征筛选；然后建立基于时滞反馈网络(NARX)的浮选品位预测子模型，并将浮选机理和数据分析相结合，共同决定网络的时滞和反馈阶次；最后，以NARX为基学习器，支持向量回归(SVR)为次级学习器，建立基于NARX和SVR的stacking集成学习浮选品位预测模型。本发明用于铝土矿实际生产过程精矿品位和尾矿品位的预测，均方根相对误差分别为4.41％和5.19％，平均相对误差分别为3.62％和4.08％。

Description

一种矿物浮选过程指标的在线预测方法

技术领域

本发明涉及一种矿物浮选过程指标的预测方法，特别是在线进行预测的方法。

背景技术

精矿品位和尾矿品位是矿物浮选过程的指标，直接决定了选矿厂的经济效益。矿石品位可通过品位分析仪在线获得，但仪器价格昂贵，维护代价非常高，且直接测量结果也需要模型校正后其准确性才能达到要求，因此，众多的选厂仍依靠人工离线化验分析。一方面浮选流程长，物料从入口到出口需要较长时间；另一方面人工采样化验周期较长，有很大的滞后性，且劳动强度大，为降低成本选厂的化验频率极低，难以满足优化控制的需求。为此，研究浮选精矿品位和尾矿品位的在线预测，实现软测量，提前预知精矿和尾矿品位则可实时采取相应的调控措施，使过程处于最佳状态，能为工厂带来非常可观的效益，具有重要的实际意义。

目前，神经网络具有较强的学习能力和泛化能力，可以充分逼近任意复杂的非线性关系。基于泡沫视觉特征和过程操作变量的浮选指标神经网络预测模型在稳态工况下具有较好的效果。由于矿物浮选是一个连续的动态过程，基于特征的静态预测模型只能反映当前时刻的过程工况对浮选指标的影响，而实际生产中各过程之间相互影响，最终的精矿与尾矿品位也是过程工艺长时间相互作用的结果，所以当前时刻的浮选指标不仅与当前时刻的物料参数及过程操作变量相关，还应与之前时刻的物料参数及过程操作变量相关；另一方面，铝土矿浮选过程复杂且过程工况波动频繁，由于矿源的复杂性，难以保证长时间的恒矿源供矿，使得矿石性质变化频繁，过程参数也不断变化，随着时间的推移，原有的数据驱动的静态模型已经不能很好的反映当前对象，而模型的更新或重训练时间空间代价高，机理模型虽然具有很好的泛化能力，但由于做了过多简化和假设，预测精度不能满足实际生产要求；再者，由于浮选变量多、流程长、内部反应复杂且存在未知扰动，使得浮选过程具有一定的不确定性，现有软测量模型所选择的浮选特征只能在一定程度上反映浮选指标但不能完全反映浮选指标。

因此，研究如何将浮选全流程机理和数据驱动相结合，对精尾矿品位进行在线预测，是克服浮选过程的滞后性和不确定性问题的基础，对提高模型的精度和泛化性，稳定生产过程，提高生产效率，选矿厂的优化控制有重要意义。

发明内容

本发明的目的是提供一种矿物浮选过程指标的在线预测方法。

本发明将机理与数据驱动相结合，在分析浮选全流程机理的基础上，通过挖掘数据间存在的潜在关系，对工艺生产数据进行预处理，首先提出了一种关联最近邻(KNN)数据缺失值插补法，基于KNN的缺失值插补步骤如下：

S1：记预测建模数据集为data，记不含缺失值的样本集为set_data，记含有缺失值的样本集为set_missing；

记set_missing中某样本的缺失值所对应的特征为a，记该缺失值所属样本中属性值不为空的特征集合为F＝{b₁,b₂,L,b_i,L,b_l}，其中l为集合中特征的个数，b_i表示第i个特征，i＝1,2,…,l；

S2：记set_data中特征a的数据列为Y_a＝[x₁,x₂,L,x_j,L,x_m]^T，F中第i个特征对应的数据列为Y_b_i＝[y₁,y₂,L,y_j,L,y_m]^T，其中m为set_data中的样本个数，x_j，y_j分别代表数据列的第j个样本值，j＝1,2,…,m；

S3：计算Y_a与Y_b_i的相关性，其中i＝1,2,…,l，皮尔逊相关系数计算公式为：

其中，分别代表Y_a和Y_b_i的均值；

S4：在F中选择皮尔逊相关系数排在前h的候选特征集合，记该缺失值所属样本所对应的候选特征的值为x_miss＝[x_b₁,x_b₂,L,x_b_p,L,x_b_h]，p＝1,2,…,h，set_data中所有样本的候选特征的值为X＝[y¹,y²,L,y^j,L,y^m]^T，j＝1,2,…,m，其中

代表第j个样本的候选特征的值，按公式(2)对x_miss和X中的h个特征所对应的数据进行无量纲化：

其中，x代表某个候选特征无量纲化前的值，x'表示无量纲化后的值，x_min，x_max分别代表该候选特征列中的最小值和最大值；

S5：计算x_miss与y^j的欧氏距离d_j，公式为：

S6：选择距离缺失值所属样本最近的K个样本，记其欧氏距离分别为d₁,d₂,…,d_z,…,d_K按距离分别为K个样本分配权重w_z，z＝1,2,…,K权值计算如式(4)

式中w_z表示第z个相似样本所占权重，d_z为缺失值所属样本与第z个相似样本的欧氏距离，C为常数，C取值为0.01；

S7：根据得到的权重估计该样本中特征a的插补值：

其中x_fill为插补后特征a的值，x_za为第z个相似样本中特征a对应的值；

重复步骤S1-S7，直至所有缺失值都已插补，得到插补后不含缺失值的数据集DATA_SET；

根据公式(1)计算DATA_SET中输入特征列与输出品位列之间的皮尔逊相关系数r，构建检验统计量如式(6)：

其中r为皮尔逊相关系数，n为样本总量；

选定显著水平α_r，则置信概率P_r＝1-α_r，按照t_r服从自由度为n-2的t分布查表得到α_r和n-2所对应的置信系数其中，n为样本个数；如果表明t_r所对应的概率p_value>P_r,则认为r显著，也即该特征列与品位列显著相关；

对所有皮尔逊相关系数r进行显著性检验后，保留与品位列具有显著相关关系的特征，构成模型输入特征的数据矩阵为：

U＝[u₁,u₂,L,u_c,L,u_M] (7)

其中u_c为第c个特征的数据列，M为特征总数；

在NARX建模过程中，特征的时滞和反馈阶次可根据实际生产经验来确定，也可通过数据关联性分析来确定。工人经验具有一定的参考性，但由于浮选过程是一个非常长的流程，其内部反应复杂且不可观测，仅凭工人经验只能大致推算出滞后阶次，有一定的局限性，因此我们将工人经验和数据关联分析相结合，共同确定定模型的时滞和反馈阶数，确定模型最终的输入维度，步骤如下：

Step1：对第c个特征，记样本总数为T+W-1，选取连续T个采样点为一个滑动窗口，其构成的时间序列数据为：

u_c′＝[u_c(1),L,u_c(T-1),u_c(T)] (8)

Step2：窗口滑动步长为1，随着时间推移共得到W个窗口，其构成的数据矩阵为：

其中U_c矩阵中每一行数据代表一个滑动窗口；

Step3：记品位特征数据的第W个窗口构成的时间序列为Y，

Step4：分别对Y中最后一行和与U_c中每一行进行皮尔逊相关性分析，设置相关系数阈值为r_t，取相关系数大于r_t的个数p_c为该特征的时滞阶次；

Step5：重复步骤Step1，Step2和Step4得到所有输入特征的时滞阶次P＝[p₁,p₂,L,p_c,L,p_M]，其中p_c为第c个特征的时滞阶次；

Step6：对Y中最后一行和之前W-1行进行皮尔逊相关性分析，根据相关性大小关系确定该输出品位特征的反馈阶次q；

Step7：确定模型最终的输入维度N＝p₁+p₂+…+p_M+q；

根据所有输入输出特征的时滞阶次和反馈阶次，共得到D个训练样本，表示为输入特征矩阵的形式为：

I＝[U₁′,U₂′,L,U_c′,L,U′_M,Y′]_D×N (11)

其中U_c′代表第c个输入特征的时滞矩阵，Y′代表输出品位特征的反馈矩阵，得到t时刻对应的品位预测模型的外部输入可表示为：

I(t)＝[U′(t),U′(t-1),L,U′(t-n_P),Y′(t-1),L,Y′(t-n_q)] (12)

式中，U′(t-n_P)表示相应模型输入参数及相应时滞n_P数值，Y′(t-n_q)表示相应模型输出参数以及相应的反馈n_q数值；

集成学习是通过构造不同的个体学习器，然后将其合并在一起来共同完成学习任务的一种方法。单学习器可能会因为误选导致模型泛化性能较差，结合多个学习器能够减少这一风险，因此集成学习是在保证模型准确度的同时也大大的提升了防止模型过拟合的能力，基于NARX和SVR建立stacking集成学习预测模型的步骤如下：

基于k折交叉验证的思想，将D个训练样本分为k个大小相等的互斥子集，则每一个子集中的样本个数为D/k，每次用k-1个子集的并集作为训练样本，剩下的一个子集作为该基学习器的测试集；每一折训练得到2β个NARX网络模型作为基学习器，其中基学习器有NARX_MLP和NARX_XGB两种，分别采用多层感知机学习算法和XGBoost学习算法，t时刻NARX_MLP第i个隐层节点数出H_i(t)为：

其中，f表示隐含层激活函数，n_P为输入特征时延长度，n_q为输出品位反馈时延长度；w_ih(t)为第t时刻第i个隐层节点与输入U′(t-h)间的权值；w_io(t)为t时刻第i个隐层节点与输出反馈量Y′(t-o)间的权值；b_i为第i个隐层节点阈值，输出层节点输出为：

式中w_i(t)为t时刻第i个隐层节点到输出层节点的权值；θ为输出层节点阈值；N为隐层节点个数；

在t时刻NARX_XGB输入和输出的关系为：

式中，f_s是表示一个具体的分类回归树CART，F表示所有可能的CART树，S是树的数量；

经过k折训练后得到k×2β个基学习器；对每一个测试样本得到2β个基学习器输出值，记第v个样本对应的基学习器的输出值向量为D个测试样本得到的基学习器输出值矩阵为其维数为D×2β；

根据小样本数据特点选择SVR作为集成学习的次级学习器，矩阵作为次级学习器的输入特征，其中每个样本输入为中的一行，对应的标签为原测试样本的标签y_v，得到次级学习器预测模型为

其中ρ为：

式中，和α_v为拉格朗日乘子，ε代表训练集中真实值和预测值的偏差，代表预测集中第g个样本对应的基学习器的输出值向量；

所得式(16)为最终的预测模型；在实际在线应用中，根据步骤Step1-Step7得到的模型输入特征及每个特征的输入和输出反馈的阶次准备模型输入，通过k×2β个NARX子模型得到k×2β个预测值，然后将这些预测值输入到SVR模型得到最终的品位预测值。

本发明通过一种关联KNN方法对缺失值进行插补，解决了因模型对缺失数据敏感而造成较大误差的问题；基于皮尔逊系数相关性分析和显著性检验探究了矿石性质和过程工况时间序列与精矿品位和尾矿品位间的关系，并筛选出精矿品位和尾矿品位的重要影响因素，解决了样本特征过多带来的冗余及误差；在此基础上，将数据关联和工人经验相结合确定模型输入时滞和输出反馈阶次，建立了基于时滞反馈神经网络的精矿品位和尾矿品位预测模型，该网络具有较强的逼近系统动态过程的能力，使得品位预测模型能够抗干扰，具有较好的稳定性和较优的动态性能；最后以时滞反馈网络为基学习器，支持向量回归为次级学习器，建立了基于stacking的集成学习模型对精尾矿品位进行预测，进一步提高了模型的精度和泛化能力，利用该方法预测得到均方根相对误差分别为4.41％和5.19％，平均相对误差分别为3.62％和4.08％。

附图说明

图1为浮选过程精尾矿品位预测方法的结构图；

图2为时滞反馈网络的结构图；

图3为本方法对精矿品位的预测结果图；

图4为本方法对尾矿品位的预测结果图；

表1为本方法对品位预测模型的性能统计结果图。

具体实施方式

如图1所示，本发明基于关联K近邻(KNN)及Pearson相关系数和显著性检验的特征筛选方法对数据进行预处理，结合浮选粗选、扫选、精选流程确定时滞反馈网络框架，然后通过大量的数据相关性分析来选取时滞反馈网络的时滞和反馈阶次，建立了基于时滞反馈网络(NARX)的浮选品位预测模型，并以多层感知机模型(MLP)和XGBoost模型为基学习器，支持向量回归模型(SVR)为次级学习器，基于stacking集成学习方法(图1虚线框中内容)对学习器进行模型融合，得到精矿品位预测值y_1和尾矿品位预测值y_2。本实例采用某浮选厂工业生产数据共450个，采样间隔为8小时，其中360个作为训练集(验证集从训练集中划分)，未来30天数据(90个)作为测试集。

1.对工艺生产数据进行预处理

对缺失值数据进行插补，步骤如下：

S1：记set_missing中某样本的缺失值所对应的特征为a，记该缺失值所属样本中属性值不为空的特征集合为F＝{b₁,b₂,L,b_i,L,b_l}，其中l为集合中特征的个数，b_i表示第i个特征，i＝1,2,…,l；

S2：记set_data中特征a的数据列为Y_a＝[x₁,x₂,L,x_j,L,x_m]^T，F中第i个特征对应的数据列为Y_b_i＝[y₁,y₂,L,y_j,L,y_m]^T，其中m为set_data中的样本个数，x_j，y_j分别代表数据列的第j个样本值，j＝1,2,…,m；

S3：计算Y_a与Y_b_i的相关性，其中i＝1,2,…,l，皮尔逊相关系数计算公式为：

其中，分别代表Y_a和Y_b_i的均值；

S4：在F中选择皮尔逊相关系数排在前5的候选特征集合，记该缺失值所属样本所对应的候选特征的值为x_miss＝[x_b₁,x_b₂,L,x_b_p,L,x_b₅]，p＝1,2,…,5，set_data中所有样本的候选特征的值为X＝[y¹,y²,L,y^j,L,y^m]^T，j＝1,2,…,m，其中代表第j个样本的候选特征的值，按公式(2)对x_miss和X中的5个特征所对应的数据进行无量纲化：

其中，x代表某个候选特征无量纲化前的值，x'表示无量纲化后的值，x_min，x_max分别代表该候选特征列中的最小值和最大值；

S5：计算x_miss与y^j的欧氏距离d_j，公式为：

S6：选择距离缺失值所属样本最近的3个样本，记其欧氏距离分别为d₁,d₂,…,d_z,…,d₃按距离分别为3个样本分配权重w_z，z＝1,2,3权值计算如式(4)

式中w_z表示第z个相似样本所占权重，d_z为缺失值所属样本与第z个相似样本的欧氏距离，C为常数，C取值为0.01；

S7：根据得到的权重估计该样本中特征a的插补值：

其中x_fill为插补后特征a的值，x_za为第z个相似样本中特征a对应的值；

重复步骤S1-S7，直至所有缺失值都已插补，得到插补后不含缺失值的数据集DATA_SET；

根据公式(1)计算DATA_SET中输入特征列与输出品位列之间的皮尔逊相关系数r，构建检验统计量如式(6)：

其中r为皮尔逊相关系数，n为样本总量；

选定显著水平为0.05，则置信概率为0.95，按照t_r服从自由度为n-2的t分布查表得到α_r和n-2所对应的置信系数其中，n为样本个数；如果表明t_r所对应的概率p_value>0.95,则认为r显著，也即该特征列与品位列显著相关。

对所有皮尔逊相关系数r进行显著性检验，保留与品位列具有显著相关关系的特征，并结合人工经验保留非显著相关的特征，非显著性相关的特征如入矿品位和粗选I段的泡沫层厚度经工人经验分析后发现与品位存在明显相关关系，也予以保留。最终确定与尾矿品位紧密相关的特征为下料量、粗选和扫选过程捕收剂的添加量、碳酸钠加入量、泡沫层厚度以及风机频率和分散剂加入量等20个特征；与精矿品位紧密相关的是风机频率、矿浆浓度、粗选I段及精选I段的捕收剂添加量、碳酸钠的加入量、下料量及粗选过程的分散剂加入量和磨矿细度等19个特征。

根据上述特征构成品位预测模型的输入特征矩阵：

U＝[u₁,u₂,L,u_c,L,u_M] (7)

其中u_c为第c个特征的数据列，M为特征总数。

2.基于生产经验和数据相关性分析选取时滞及反馈阶次

在NARX建模过程中，特征的时滞和反馈阶次可根据实际生产经验来确定，也可通过数据关联性分析来确定。工人经验具有一定的参考性，但由于浮选过程是一个非常长的流程，其内部反应复杂且不可观测，仅凭工人经验只能大致推算出滞后阶次，有一定的局限性，因此我们将工人经验和数据关联分析相结合，共同决定模型的时滞和反馈阶数，确定模型最终的输入维度，步骤如下：

Step1：对第c个特征，记样本总数为W+T+1，选取连续T个采样点为一个滑动窗口，其构成的时间序列数据为：

u_c′＝[u_c(1),L,u_c(T-1),u_c(T)] (8)

Step2：窗口滑动步长为1，随着时间推移共得到W(W＝12)个窗口，其构成的数据矩阵为：

其中U_c矩阵中每一行数据代表一个滑动窗口；

Step3：记品位特征数据的第W个窗口构成的时间序列为Y，

Step4：分别对Y中最后一行和与U_c中每一行进行皮尔逊相关性分析，设置相关系数阈值为r_t，取相关系数大于r_t的个数p_c为该特征的时滞阶次；

Step5：重复步骤Step1，Step2和Step4得到所有输入特征的时滞阶次P＝[p₁,p₂,L,p_c,L,p_M]，其中p_c为第c个特征的时滞阶次，精矿品位预测模型特征阶次为P_精＝[2，3，0，3，1，3，3，3，3，2，4，2，4，2，11，11，11，11，6]，尾矿品位预测模型特征阶次为P_尾＝[2，5，0，2，5，1，3，5，5，4，5，3，4，3，11，11，11，11，11，5]；

Step6：对Y中最后一行和之前11行进行皮尔逊相关性分析，根据相关性大小关系确定该输出品位特征的反馈阶次q_精＝q_尾＝4；

Step7：确定模型最终的输入维度N＝p₁+p₂+…+p_M+q，精矿品位预测模型输入维度为85，尾矿品位预测模型输入维度为111；

根据所有输入输出特征的时滞阶次和反馈阶次，共得到D(D＝360)个训练样本，表示为输入特征矩阵的形式为：

I＝[U₁′,U₂′,L,U_c′,L,U′_M,Y′]_D×N (11)

其中U_c′代表第c个输入特征的时滞矩阵，Y′代表输出品位特征的反馈矩阵。

3.基于时滞反馈网络NARX和支持向量回归SVR建立stacking集成学习预测模型

根据上一步确定的各输入特征的时滞和输出特征的反馈阶次，t时刻对应的NARX网络外部输入可表示为：

I(t)＝[U′(t),U′(t-1),L,U′(t-n_P),Y′(t-1),L,Y′(t-n_q)] (12)

式中，U′(t-n_P)表示相应模型输入参数及相应时滞n_P数值，Y′(t-n_q)表示相应模型输出参数以及相应的反馈n_q数值，NARX网络结构如图2所示；

本实施例中，选择多层感知机算法和Xgboost算法作为基学习器的模型算法，分别建立NARX_MLP和NARX_XGB两种基学习器作为预测子模型。在t时刻NARX_MLP第i个隐层节点数出H_i(t)为：

其中，f表示隐含层激活函数，n_P为输入特征时延长度，n_q为输出品位反馈时延长度；w_ih(t)为第t时刻第i个隐层节点与输入U′(t-h)间的权值；w_io(t)为t时刻第i个隐层节点与输出反馈量Y′(t-o)间的权值；b_i为第i个隐层节点阈值。输出层节点输出为：

式中w_i(t)为t时刻第i个隐层节点到输出层节点的权值；θ为输出层节点阈值；N为隐层节点个数；

在t时刻NARX_XGB输入和输出的关系为：

式中，f_s是表示一个具体的分类回归树CART，F表示所有可能的CART树，S是树的数量；

为了得到充足的样本来训练次级学习器，同时又避免使用训练过基学习器的样本，采用10折交叉验证法来训练集成模型，将360个训练样本分为10个大小相等的互斥子集，则每一个子集中的样本个数为36，每次用9个子集的并集作为训练样本，剩下的一个子集作为该基学习器的测试集；每一折训练得到4个NARX网络模型作为基学习器；经过10折训练后得到40个基学习器；对每一个测试样本得到4个基学习器输出值，记第v个样本对应的基学习器的输出值向量为360个测试样本得到的基学习器输出值矩阵为其维数为360 14；

根据小样本数据特点选择SVR作为集成学习的次级学习器，矩阵作为次级学习器的输入特征，其中每个样本输入为中的一行，对应的标签为原测试样本的标签y_v，得到次级学习器预测模型为：

其中ρ为：

式中，和α_v为拉格朗日乘子，ε代表训练集中真实值和预测值的偏差，代表预测集中第g个样本对应的基学习器的输出值向量；

所得式(16)为最终的预测模型；在实际应用中，根据步骤Step1-Step7得到所选输入特征及其阶次准备模型输入，首先通过40个NARX子模型得到40个预测值，然后将这些预测值输入到SVR模型得到最终的品位输出值。

对工厂连续生产30天的精矿品位和尾矿品位预测结果如图3和图4所示，模型性能统计结果见表1，均方根相对误差分别为4.41％和5.19％，平均相对误差分别为3.62％和4.08％。，表明预测模型具有较强的逼近系统动态过程的能力，使得品位预测模型能够抗干扰，具有较好的泛化性能和预测精度。

表1

Claims (1)

1.一种矿物浮选过程指标的在线预测方法，其特征在于包括以下步骤：

(1)对工艺生产数据进行预处理

记预测建模数据集为data，记不含缺失值的样本集为set_data，记含有缺失值的样本集为set_missing，处理步骤如下：

S1：记set_missing中某样本的缺失值所对应的特征为a，记该缺失值所属样本中属性值不为空的特征集合为F＝{b₁,b₂,L,b_i,L,b_l}，其中l为集合中特征的个数，b_i表示第i个特征，i＝1,2,…,l；

S2：记set_data中特征a的数据列为Y_a＝[x₁,x₂,L,x_j,L,x_m]^T，F中第i个特征对应的数据列为Y_b_i＝[y₁,y₂,L,y_j,L,y_m]^T，其中m为set_data中的样本个数，x_j，y_j分别代表数据列Y_a和数据列Y_b_i的第j个样本值，j＝1,2,…,m；

S3：计算Y_a与Y_b_i的相关性，其中i＝1,2,…,l，皮尔逊相关系数计算公式为：

其中，分别代表Y_a和Y_b_i的均值；

S4：在F中选择皮尔逊相关系数排在前h的候选特征集合，记该缺失值所属样本所对应的候选特征的值为x_miss＝[x_b₁,x_b₂,L,x_b_p,L,x_b_h]，p＝1,2,…,h，set_data中所有样本的候选特征的值为X＝[y¹,y²,L,y^j,L,y^m]^T，j＝1,2,…,m，其中p＝1,2,…,h，代表第j个样本的候选特征的值，按公式(2)对x_miss和X中的h个特征所对应的数据进行无量纲化：

其中，x代表某个候选特征无量纲化前的值，x'表示无量纲化后的值，x_min，x_max分别代表该候选特征列中的最小值和最大值；

S5：计算x_miss与y^j的欧氏距离d_j，公式为：

S6：选择距离缺失值所属样本最近的K个样本，记其欧氏距离分别为d₁,d₂,…,d_z,…,d_K按距离分别为K个样本分配权重w_z，z＝1,2,…,K权值计算如式(4)

式中w_z表示第z个相似样本所占权重，d_z为缺失值所属样本与第z个相似样本的欧氏距离，C为常数，C取值为0.01；

S7：根据得到的权重估计该样本中特征a的插补值：

其中x_fill为插补后特征a的值，x_za为第z个相似样本中特征a对应的值；

重复步骤S1-S7，直至所有缺失值都已插补，得到插补后不含缺失值的数据集DATA_SET；

(2)根据预处理后的工艺数据进行输入特征筛选

计算输入特征列与输出品位列之间的皮尔逊相关系数r，构建检验统计量如式(6)

其中r为皮尔逊相关系数，n为样本总量；

选定显著水平α_r，则置信概率P_r＝1-α_r，按照t_r服从自由度为n-2的t分布查表得到α_r和n-2所对应的置信系数其中，n为样本个数；如果表明t_r所对应的概率p_value>P_r,则认为r显著，也即该特征列与品位列显著相关；

对所有皮尔逊相关系数r进行显著性检验后，保留与品位列具有显著相关关系的特征，构成模型输入特征的数据矩阵为：

U＝[u₁,u₂,L,u_c,L,u_M] (7)

其中u_c为第c个特征的数据列，M为特征总数；

(3)基于生产经验和数据相关性分析选取时滞及反馈阶次

时滞反馈阶次确定步骤如下Step1到Step7：

Step1：对第c个特征，记样本总数为T+W-1，选取连续T个采样点为一个滑动窗口，其构成的时间序列数据为：

u′_c＝[u_c(1),L,u_c(T-1),u_c(T)] (8)

Step2：窗口滑动步长为1，随着时间推移共得到W个窗口，其构成的数据矩阵为：

其中U_c矩阵中每一行数据代表一个滑动窗口；

Step3：记品位特征数据的第W个窗口构成的时间序列为Y，

Step4：分别对Y中最后一行和与U_c中每一行进行皮尔逊相关性分析，设置相关系数阈值为r_t，取相关系数大于r_t的个数p_c为该特征的时滞阶次；

Step5：重复步骤Step1，Step2和Step4得到所有输入特征的时滞阶次P＝[p₁,p₂,L,p_c,L,p_M]，其中p_c为第c个特征的时滞阶次；

Step6：对Y中最后一行和之前W-1行进行皮尔逊相关性分析，根据相关性大小关系确定该输出品位特征的反馈阶次q；

Step7：确定模型最终的输入维度N＝p₁+p₂+…+p_M+q；

根据所有输入输出特征的时滞阶次和反馈阶次，共得到D个训练样本，表示为输入特征矩阵的形式为：

I＝[U′₁,U′₂,L,U′_c,L,U′_M,Y′]_D×N (11)

其中U′_c代表第c个输入特征的时滞矩阵，Y′代表输出品位特征的反馈矩阵，得到t时刻对应的品位预测模型的外部输入，表示为：

I(t)＝[U′(t),U′(t-1),L,U′(t-n_P),Y′(t-1),L,Y′(t-n_q)](12)

式中，U′(t-n_P)表示相应模型输入参数及相应时滞n_P数值，

Y′(t-n_q)表示相应模型输出参数以及相应的反馈n_q数值；

(4)基于时滞反馈网络NARX和支持向量回归SVR建立stacking集成学习预测模型

基于k折交叉验证的思想，将D个训练样本分为k个大小相等的互斥子集,则每一个子集中的样本个数为D/k，每次用k-1个子集的并集作为训练样本，剩下的一个子集作为该基学习器的测试集；每一折训练得到2β个NARX网络模型作为基学习器，其中基学习器有NARX_MLP和NARX_XGB两种，分别采用多层感知机学习算法和XGBoost学习算法，在t时刻NARX_MLP第i个隐层节点输出H_i(t)为：

其中，f表示隐含层激活函数，n_P为输入特征时延长度，n_q为输出品位反馈时延长度；w_ih(t)为t时刻第i个隐层节点与输入U′(t-h)间的权值；w_io(t)为t时刻第i个隐层节点与输出反馈量Y′(t-o)间的权值；b_i为第i个隐层节点阈值；输出层节点输出为：

式中w_i(t)为t时刻第i个隐层节点到输出层节点的权值；θ为输出层节点阈值；N为隐层节点个数；

在t时刻NARX_XGB输入和输出的关系为：

式中，f_s是表示一个具体的分类回归树CART，F表示所有可能的CART树，S是树的数量；

经过k折训练后得到k×2β个基学习器；对每一个测试样本得到2β个基学习器输出值，记第v个样本对应的基学习器的输出值向量为D个测试样本得到的基学习器输出值矩阵为其维数为D×2β；

根据小样本数据特点选择SVR作为集成学习的次级学习器，矩阵作为次级学习器的输入特征，其中每个样本输入为中的一行，对应的标签为原测试样本的标签y_v，得到次级学习器预测模型为

其中ρ为：

式中，和α_v为拉格朗日乘子，ε代表训练集中真实值和预测值的偏差，代表预测集中第g个样本对应的基学习器的输出值向量。