CN111694327B - 一种基于混合独立成分分析算法的工业过程监测方法 - Google Patents

Abstract

本发明公开一种基于混合独立成分分析算法的工业过程监测方法，旨在提出一种混合独立成分分析算法，并利用该算法实施非高斯过程监测。具体来讲，本发明方法首先利用粒子群优化算法逐个提取训练数据中潜藏的非高斯分布的独立成分信息，其次对剩余的高斯分布信息，利用时间序列相关的主成分分析算法进一步提取时间序列独立成分，最后利用这些提取的独立成分实施非高斯过程监测。与传统方法相比，本发明方法使用粒子群优化算法优化求解非高斯分布的独立成分，保证了求解过程中的全局最优性。其次，本发明方法对服从高斯分布的独立成分进行深入挖掘，考虑到了时间序列上的二阶相关性。可以说，本发明方法是一种更为优选的非高斯过程监测方法。

Description

一种基于混合独立成分分析算法的工业过程监测方法

技术领域

本发明涉及一种数据驱动的过程监测方法，尤其涉及一种基于混合独立成分分析算法的工业过程监测方法。

背景技术

随着先进测量技术与计算技术的飞速发展，生产过程对象可以离线存储与在线测量海量的数据，现代工业过程已进入工业“大数据”时代。这些采样数据蕴含着能体现生产过程运行状态的信息，利用采样数据实施过程运行状态的监测于是乎得到了较多学者们的青睐。近年来，无论是学术界还是工业界，都投入了大量的人力与物力研究可靠的数据驱动过程监测方法。在数据驱动过程监测研究领域，统计过程监测是被研究得最多的方法，其中当以主成分分析(Principal Component Analysis，PCA)与独立成分分析(IndependentComponent Analysis，ICA)为最主流的实施技术手段。一般而言，由于ICA算法能够挖掘出数据中潜藏的非高斯成分信息，更能揭示对象的本质，因此更适合于非高斯过程对象的监测。

虽然ICA算法能通过高阶统计量挖掘出本质特征信息，但是提取的独立成分有些可能依旧服从高斯分布。而对高斯分布独立成分信息的挖掘，传统ICA算法未能考虑时间序列上的二阶统计量。因此，传统基于ICA算法的非高斯过程监测方法还有待改进与提升。此外，利用ICA算法提取独立成分时，一般都采用FastICA迭代求取。然而，FastICA算法是基于牛顿法求解优化问题的，很容易陷入局部最优，而且独立成分的非高斯性大小不是按照提取先后自动排列的。为解决这个问题，可以使用智能优化算法，比如粒子群优化算法，来保证各独立成分的非高斯性。

发明内容

本发明所要解决的主要技术问题是：提出一种混合独立成分分析算法，首先利用粒子群优化算法逐个提取训练数据中潜藏的非高斯分布的独立成分信息，其次对剩余的高斯分布信息，利用时间序列相关的主成分分析算法进一步提取时间序列独立成分，最后利用这些提取的独立成分实施非高斯过程监测。

本发明解决上述技术问题所采用的技术方案为：一种基于混合独立成分分析算法的工业过程监测方法，包括以下所示步骤：

步骤(1)：采集工业过程正常运行状态下的样本数据，组成矩阵X∈R^n×m，并计算矩阵X中各列向量的均值μ₁，μ₂，...，μ_m以及标准差δ₁，δ₂，...，δ_m，从而组成均值向量μ＝[μ₁，μ₂，...，μ_m]^T与标准差向量δ＝[δ₁，δ₂，...，δ_m]，其中，n为训练样本数，m为过程测量变量数，R为实数集，R^n×m表示n×m维的实数矩阵，上标号T表示矩阵或向量的转置。

步骤(2)：根据公式

对矩阵X实施标准化处理得到矩阵

其中，U∈R^n×m是由n个相同的均值向量μ组成的矩阵，即U＝[μ，μ，...，μ]^T，对角矩阵Φ中对角线上的元素由标准差向量δ组成。

步骤(3)：设置矩阵

确定如下所示的优化目标函数，并利用粒子群优化算法求解得到变换向量w∈R^m×1及其对应的目标函数值J_w：

上式中，E{ }表示求取向量的均值，

表示求取

的绝对值。步骤(3)中利用粒子群优化算法的求解过程如下所示：

步骤(3.1)：步骤(3.1)：初始化flag＝1，设置最大迭代次数为Imax、粒子个数为N、粒子最大飞行速度为Vmax、学习因子c₁与c₂。

步骤(3.2)：在区间[-1，1]上随机产生一个m×N维的实数矩阵P，在区间[-Vmax，Vmax]上随机产生一个m×N维的实数矩阵V，实数矩阵P中的各个列向量p₁，p₂，...，p_N即为各粒子的位置向量，实数矩阵V中的各个列向量v₁，v₂，...，v_N即为各个粒子的速度向量。

步骤(3.3)：依次以列向量p₁，p₂，...，p_N为变换向量，根据公式

计算相应的目标函数值J₁，J₂，...，J_N。

步骤(3.4)：将各个粒子在整个迭代历史中取得最大目标函数值的位置向量g₁，g₂，...，g_N记录为矩阵G＝[g₁，g₂，...，g_N]，并将整个迭代历史中取得最大目标函数值所对应的粒子位置向量记录为向量θ∈R^m×1。

步骤(3.5)：根据公式v_i＝v_i+c₁r₁(g_i-p_i)+c₂r₂(θ-p_i)更新各个粒子对应的速度向量v_i后，判断速度向量v_i中各元素是否在区间[-Vmax，Vmax]之间？若是，则保持该元素不变；若否，则将该元素变成最靠近的区间边界值。

步骤(3.6)：根据公式p_i＝p_i+v_i计算各个粒子的位置向量p_i后，判断位置向量p_i中各元素是否在区间[-1，1]之间？若是，则保持该元素不变；若否，则将该元素变成最靠近的区间边界值。

步骤(3.7)：判断是否满足条件：flag＜Imax？若是，则置flag＝flag+1后返回步骤(3.3)；若否，则得到变换向量w＝θ，并计算目标函数值

步骤(4)：根据公式

与

分别计算独立成分向量s与载荷向量a后，根据公式

更新矩阵

并判断是否满足条件：J_w＜0.1？若否，则重复步骤(3)至步骤(4)；若是，则将所有的变换向量组成变换矩阵W、将所有的独立成分向量组成独立成分矩阵S、将所有的载荷向量组成载荷矩阵A。

步骤(5)：设置自相关阶数为D后，根据公式

得到矩阵X₁，X₂，...，X_D，其中

表示将矩阵

中的第d行至第n-D+d行向量组成矩阵的操作，d＝1，2，...，D。

步骤(6)：根据公式C＝X_D ^TX₁+X_D ^TX₂+...+X_D ^TX_D-1与C＝(C+C^T)/2计算时间序列协方差矩阵C，并计算C所有特征值λ₁，λ₂，...，λ_m所对应的特征向量α₁，α₂，...，α_m，这里要求特征向量都是单位长度。

步骤(7)：将特征值λ₁，λ₂，...，λ_m取绝对值后得到|λ₁|，|λ₂|，...，|λ_m|，并根据绝对值大小选取前k个最大值所对应的特征向量，从而组成时间序列载荷矩阵U。

上述步骤最终将原数据矩阵

分解成：

其中F为误差矩阵。由此可见，S提取的是非高斯分布的独立成分，

提取的是时间序列相关的独立成分。

步骤(8)：根据公式I²＝diag(SΛ_s ^-1S^T)、T²＝diag(ΘΛ_Θ ^-1Θ^T)、和Q＝diag(FF^T)计算出监测指标向量I²、T²、和Q，并利用核密度估计法分别得到在99％的置信限条件下的控制上限

和Q_lim，其中Λ_s与Λ_Θ分别为S与Θ的协方差矩阵。

步骤(9)：收集新采样时刻的数据样本x_new∈R^m×1，并根据公式

对x_new实施标准化处理得到

步骤(10)：分别根据公式

和

计算得到非高斯独立成分向量s_new、时间序列独立成分向量u_new、和误差向量ε_new。

步骤(11)：根据公式

和Q_new＝ε_newε_new ^T计算得到当前采样时刻的监测指标

和Q_new。

步骤(12)：判断是否满足条件

且

且Q_new≤Q_lim？若是，则当前时刻生产过程对象正常运行，返回步骤(9)继续监测下一时刻的样本数据；若否，则当前时刻生产过程对象出现异常工况。

与传统方法相比，本发明方法的优势在于：

首先，本发明方法使用粒子群优化算法优化求解非高斯分布的独立成分，保证了求解过程中的全局最优性。其次，本发明方法对服从高斯分布的独立成分进行深入挖掘，考虑到了时间序列上的二阶相关性。关键的是，本发明中所涉及的混合独立成分分析算法是一种全新的非高斯过程特征提取算法。此外，具体实施案例中将会验证本发明方法的优越性，从而说明本发明方法是一种更为优选的非高斯过程监测方法。

附图说明

图1为本发明方法离线建模阶段的实施流程图。

图2为本发明方法在线监测的实施流程图。

图3为TE过程故障监测详情对比图。

具体实施方式

下面结合附图与具体的实施案例对本发明方法进行详细的说明。

本发明公开一种基于混合独立成分分析算法的工业过程监测方法，下面结合一个具体的案例来说明本发明方法的具体实施过程，以及相对于传统方法的优越性。

根据如下所示公式②仿真一非高斯过程对象：

上式中，f_s＝1000，g_1，t+1与g_2，t+1为高斯分布的白噪声。按照上式②产生1000个样本数据作为训练数据矩阵X∈R^100×3，离线建模阶段的实施流程如图1所示，具体包括如下所示步骤：

(1)采集工业过程正常运行状态下的样本数据，组成矩阵X∈R^1000×3，并计算矩阵X中各列向量的均值μ₁，μ₂，μ₃以及标准差δ₁，δ₂，δ₃，即可组成均值向量μ＝[μ₁，μ₂，μ₃]^T与标准差向量δ＝[δ₁，δ₂，δ₃]。

(2)根据公式

对矩阵X实施标准化处理得到

其中U∈R^n×m是由n个相同的均值向量μ组成的矩阵，对角矩阵Φ中对角线上的元素由标准差向量δ组成。

(3)设置矩阵

后，确定如前述公式①所示的优化目标函数，并利用粒子群优化算法求解得到变换向量w∈R^m×1及其对应的目标函数值J_w。

(4)根据公式

与

分别计算独立成分向量s与载荷向量a后，根据公式

更新矩阵

并判断是否满足条件：J_w＜0.1？若否，则重复步骤(3)至步骤(4)；若是，则将所有的变换向量组成变换矩阵W、将所有的独立成分向量组成独立成分矩阵S、将所有的载荷向量组成载荷矩阵A。

(5)设置自相关阶数为D＝2后，根据公式

得到矩阵X₁，X₂。

(6)根据公式C＝X₂ ^TX₁与C＝(C+C^T)/2计算矩阵

的时间序列协方差矩阵C，并计算C所有特征值λ₁，λ₂，...，λ₃₃所对应的特征向量α₁，α₂，...，α₃₃，这里要求特征向量都是单位长度。

(7)将特征值λ₁，λ₂，...，λ₃₃取绝对值后得到|λ₁|，|λ₂|，...，|λ₃₃|，并根据绝对值大小选取前k＝2个最大值所对应的特征向量，从而组成时间序列载荷矩阵U。

(8)根据公式I²＝diag(SΛ_s ^-1S^T)、T²＝diag(ΘΛ_Θ ^-1Θ^T)、和Q＝diag(FF^T)计算出监测指标向量I²、T²、和Q，并利用核密度估计法分别得到在99％的置信限条件下的控制上限

和Q_lim。

利用公式②再次生成800个样本数据，并且从第301个样本数据采样点开始，给s_1，t+1添加幅值为+4的阶跃跳变，并按照如图2所示的在线过程监测流程实施在线故障检测，具体包括以下实施步骤。

(9)收集新采样时刻的数据样本x_new∈R^3×1，并根据公式

对x_new实施标准化处理得到

(10)分别根据公式

和

计算得到非高斯独立成分向量s_new、时间序列独立成分向量u_new、和误差向量ε_new。

步骤(11)：根据公式

和Q_new＝ε_newε_new ^T计算得到当前采样时刻的监测指标

和Q_new；

步骤(12)：判断是否满足条件

且

且Q_new≤Q_lim？若是，则当前时刻生产过程对象正常运行，返回步骤(9)继续监测下一时刻的样本数据；若否，则当前时刻生产过程对象出现异常工况。

如图3所示，本发明方法与基于PCA的传统过程监测方法在监测故障数据时的监测详情。从图3中可以很明显地发现，本发明方法在故障发生后漏报情况明显优越于传统。因此，可以说本发明方法具有更可靠的过程监测性能。

上述实施案例只用来解释说明本发明的具体实施，而不是对本发明进行限制。在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改，都落入本发明的保护范围。

Claims (2)

1.一种基于混合独立成分分析算法的工业过程监测方法，其特征在于，包括以下步骤：

步骤(1)：采集工业过程正常运行状态下的样本数据，组成矩阵X∈R^n×m，并计算矩阵X中各列向量的均值μ₁，μ₂，…，μ_m以及标准差δ₁，δ₂，…，δ_m，从而组成均值向量μ＝[μ₁，μ₂，…，μ_m]^T与标准差向量δ＝[δ₁，δ₂，…，δ_m]，其中，n为训练样本数，m为过程测量变量数，R为实数集，R^{n ×m}表示n×m维的实数矩阵，上标号T表示矩阵或向量的转置；

步骤(2)：根据公式

对矩阵X实施标准化处理得到矩阵

其中，U∈R^n×m是由n个相同的均值向量μ组成的矩阵，即U＝[μ，μ，…，μ]^T，对角矩阵Φ中对角线上的元素由标准差向量δ组成；

步骤(3)：设置矩阵

后，确定如下所示的优化目标函数，并利用粒子群优化算法求解得到变换向量w∈R^m×1及其对应的目标函数值J_w：

上式中，E{}表示求取向量的均值；

步骤(4)：根据公式

与

分别计算独立成分向量s与载荷向量a后，根据公式

更新矩阵

并判断是否满足条件：J_w＜0.1；若否，则重复步骤(3)至步骤(4)；若是，则将所有的变换向量组成变换矩阵W、将所有的独立成分向量组成独立成分矩阵S、将所有的载荷向量组成载荷矩阵A；

步骤(5)：设置自相关阶数为D后，根据公式

得到矩阵X₁，X₂，…，X_D，其中

表示将矩阵

中的第d行至第n-D+d行向量组成矩阵的操作，d＝1，2，…，D；

步骤(6)：根据公式C＝X_D ^TX₁+X_D ^TX₂+…+X_D ^TX_D-1与C＝(C+C^T)/2计算时间序列协方差矩阵C，并计算C所有特征值λ₁，λ₂，…，λ_m所对应的特征向量α₁，α₂，…，α_m，这里要求特征向量都是单位长度；

步骤(7)：将特征值λ₁，λ₂，…，λ_m取绝对值后得到|λ₁|，|λ₂|，…，|λ_m|，并根据绝对值大小选取前k个最大值所对应的特征向量，从而组成时间序列载荷矩阵U；

上述步骤最终将原数据矩阵

分解成：

其中F为误差矩阵，

是时间序列相关的独立成分矩阵；

步骤(8)：根据公式I²＝diag(SΛ_s ^-1S^T)、T²＝diag(ΘΛ_Θ ^-1Θ^T)、和Q＝diag(FF^T)计算出监测指标向量I²、T²、和Q，并利用核密度估计法分别得到在99％的置信限条件下的控制上限

和Q_lim，其中Λ_s与Λ_Θ分别为S与Θ的协方差矩阵；

步骤(9)：收集新采样时刻的数据样本x_new∈R^m×1，并根据公式

对x_new实施标准化处理得到

步骤(10)：根据公式

和

分别计算得到非高斯独立成分向量s_new、时间序列独立成分向量u_new、和误差向量ε_new；

步骤(11)：根据公式

和Q_new＝ε_newε_new ^T计算得到当前采样时刻的监测指标

和Q_new；

步骤(12)：判断是否满足条件

且

且Q_new≤Q_lim；若是，则当前时刻生产过程对象正常运行，返回步骤(9)继续监测下一时刻的样本数据；若否，则当前时刻生产过程对象出现异常工况。

2.根据权利要求1所述的一种基于混合独立成分分析算法的工业过程监测方法，其特征在于，所述步骤(3)的具体实施过程如下所示：

步骤(3.1)：初始化flag＝1，设置最大迭代次数为Imax、粒子个数为N、粒子最大飞行速度为Vmax、学习因子c₁与c₂；

步骤(3.2)：在区间[-1，1]上随机产生一个m×N维的实数矩阵P，在区间[-Vmax，Vmax]上随机产生一个m×N维的实数矩阵V，实数矩阵P中的各个列向量p₁，p₂，…，p_N即为各粒子的位置向量，实数矩阵V中的各个列向量v₁，v₂，…，v_N即为各个粒子的速度向量；

步骤(3.3)：依次以列向量p₁，p₂，…，p_N为变换向量，根据公式

计算相应的目标函数值J₁，J₂，…，J_N；

步骤(3.4)：将各个粒子在整个迭代历史中取得最大目标函数值的位置向量g₁，g₂，…，g_N记录为矩阵G＝[g₁，g₂，…，g_N]，并将整个迭代历史中取得最大目标函数值所对应的粒子位置向量记录为向量θ∈R^m×1；

步骤(3.5)：根据公式v_i＝v_i+c₁r₁(g_i-p_i)+c₂r₂(θ-p_i)更新各个粒子对应的速度向量v_i后，判断速度向量v_i中各元素是否在区间[-Vmax，Vmax]之间；若是，则保持该元素不变；若否，则将该元素变成最靠近的区间边界值；

步骤(3.6)：根据公式p_i＝p_i+v_i计算各个粒子的位置向量p_i后，判断位置向量p_i中各元素是否在区间[-1，1]之间；若是，则保持该元素不变；若否，则将该元素变成最靠近的区间边界值；

步骤(3.7)：判断是否满足条件：flag＜Imax；若是，则置flag＝flag+1后返回步骤(3.3)；若否，则得到变换向量w＝θ，并计算目标函数值

Images