CN104914847A - 基于方向核偏最小二乘的工业过程故障诊断方法 - Google Patents

Abstract

本发明涉及基于方向核偏最小二乘的工业过程故障诊断方法，该方法获取易出故障的工业过程的输入变量和输出变量的历史正常数据，对该历史正常数据进行基于方向核偏最小二乘运算，计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差的控制限，采集工业过程的输入变量的采样数据，对该采样数据进行基于方向核偏最小二乘运算，得到采样数据的过程监测统计量和采样数据的平方预测误差，当采样数据的过程监测统计量控制限或者采样数据的平方预测误差超出控制限，则该采样数据中具有一种故障，获取已知故障类型的历史故障数据，对已知故障类型的历史故障数据进行基于霍特林统计量重构和基于平方预测误差重构，判断采样数据的故障类型。

Description

基于方向核偏最小二乘的工业过程故障诊断方法

技术领域

本发明属于工业过程的故障监测与诊断技术领域，具体涉及基于方向核偏最小二乘的工业过程故障诊断方法。

背景技术

随着现代科学技术的不断发展，特别是计算机科学水平与控制水平的飞速提高，现代工业过程越来越呈现出大型化、复杂化、整体化、高速化的形式。在提高工业过程生产效率和产量的同时，如何提高工业过程系统的安全性、可靠性，防止和杜绝影响系统正常运行的故障的发生和发展就成为一个重要的有待解决的问题。过程监测就是一门为了解决这类问题而发展起来的技术。其目的一方面在于通过监测系统的行为和特性，加深对系统的了解，以便更好地控制系统，提高产品的产量和质量，提高经济效益；另一方面在于避免了系统的故障、设备的损坏和人员的伤亡，减少故障所带来的经济损失，提高产业的经济和社会效益。传统的偏最小二乘(PLS)建模、监测与诊断方法虽然应用广泛，但其方法本身仍然存在一些问题，阻碍了生产过程的精确建模与准确监测。第一个问题是PLS的残差中仍然包含与输入变量有关的变异。由于输出相关变异的存在且没有被解释利用，PLS在监测和诊断与输出变量有关的故障时具有局限性，不能达到最好的效果。第二个问题是PLS残差空间中的变异量很大，不适合用平方预测误差(SPE)监测统计量对其进行监测。这是由于PLS主元不是按照方差大小排列，很多包含大方差的主元被留到了PLS残差子空间中。因此，PLS残差子空间中的变异必须被减少。传统的主元分析(PCA)故障重构方法是将故障数据空间分解成两个互相垂直的子空间，主元子空间和残差子空间。PCA是保持最主要的数据分布方向，这些方向可以有效表示数据分布特征，但是PCA模型只是研究了故障数据的内部关系，不能够有效隔离数据中的故障信息和正常信息，而且基于PCA的重构对于关注产品质量的生产过程适应性较差。另外实际工业过程中，变量之间往往呈现出非线性特征，利用传统的线性方法进行重构也是不能达到满意的效果，所以需要对传统重构方法进行改进，提高故障诊断水平。

发明内容

针对现有技术的不足，本发明提出基于方向核偏最小二乘的工业过程故障诊断方法。

本发明技术方案如下：

基于方向核偏最小二乘的工业过程故障诊断方法，包括以下步骤：

步骤1：获取易出故障的工业过程的输入变量和输出变量的历史正常数据，对该历史正常数据进行基于方向核偏最小二乘运算，令历史正常输入数据映射的高维特征空间F的PLS残差中有i个主元，求得该高维特征空间F的PLS残差中与输出变量相关的变异E_r及其主元T_r，则获得历史正常数据新的高维特征空间F的主元T_d＝[T，T_r]，计算历史正常数据的霍特林统计量和历史正常数据的平方预测误差SPE_d，并计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差SPE_d的控制限；

步骤1.1：获取易出故障的工业过程的输入变量和输出变量的历史正常数据：对m个历史正常数据进行n次采样，得到历史正常输入数据矩阵x和历史正常输出数据矩阵y；

步骤1.2：对历史正常输入数据矩阵x和历史正常输出数据矩阵y进行标准化处理，得到预处理后的历史正常输入数据X和预处理后的历史正常输出数据Y；

步骤1.3：选取非线性变换φ(X)，将预处理后的历史正常输入数据X映射到历史正常输入数据的高维特征空间F，利用径向基内积核函数，求出历史正常输入数据的初始核矩阵K₁；

步骤1.4：对预处理后的历史正常输入数据X和预处理后的历史正常输出数据Y进行KPLS运算，求得预处理后的历史正常输入数据X的PLS主元T、预处理后的历史正常输出数据Y的PLS主元U和经过α次迭代后使预处理后的历史正常输出数据Y的PLS主元U收敛的历史正常输出数据的核矩阵K_α+1；

步骤1.5：令高维特征空间F的PLS残差中有i个主元，求得高维特征空间F的PLS残差中与输出变量相关的变异E_r及其主元T_r，则获得新的高维特征空间F的主元T_d＝[T，T_r]；

步骤1.6：计算历史正常数据的霍特林统计量和历史正常数据的平方预测误差SPE_d，并计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差SPE_d的控制限。

步骤2：采集工业过程的输入变量的采样数据，对该采样数据进行基于方向核偏最小二乘运算，令采样数据映射的高维特征空间F_new的PLS残差中有i个主元，求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new及其主元T_r，new，则获得采样数据新的高维特征空间F_new的主元T_d，new＝[T_new，T_r，new]，计算采样数据的过程监测统计量和采样数据的平方预测误差SPE_d，new；

步骤2.1：采集工业过程的输入变量的采样数据：对m个采样数据进行n次采样，得到采样数据矩阵x_new，并对其进行标准化处理，得到预处理后的采样数据X_new；

步骤2.2：利用非线性变换φ(X)将预处理后的采样数据X_new映射到采样数据的高维特征空间F_new，利用径向基内积核函数，求出采样数据的初始核矩阵K_new，1；

步骤2.3：对预处理后的采样数据X_new进行KPLS运算，求得预处理后的采样数据X_new的主元T_new和经过α次迭代后的采样数据的核矩阵K_new，α+1；

步骤2.4：令采样数据的高维特征空间F_new的PLS残差中有i个主元，求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new及其主元T_r，new，获得采样数据的新的高维特征空间F_new的主元T_d，new＝[T_new，T_r，new]；

步骤2.5：计算采样数据的过程监测统计量和采样数据的平方预测误差SPE_d，new。

步骤3：当采样数据的过程监测统计量超出历史正常数据的过程监测统计量的控制限或者采样数据的平方预测误差SPE_d，new超出历史正常数据的平方预测误差SPE_d的控制限，则该采样数据中具有一种故障，执行步骤4，否则，将该采样数据视为正常数据；

步骤4：获取已知故障类型的历史故障数据，对已知故障类型的历史故障数据进行基于霍特林统计量重构和基于平方预测误差重构，判断采样数据的故障类型；

步骤4.1：获取已知故障类型的L种历史故障数据X_f，1，X_f，2，...，X_f，L；

步骤4.2：选取已知故障类型的L种历史故障数据中的第l类历史故障数据X_f，l，l＝1，2，...，L，将高维特征空间的历史正常输入数据φ(X)沿着高维特征空间的第l类历史故障数据φ(X_f，l)的故障方向进行重构，重构出高维特征空间的第l类历史故障数据φ(X_f，l)出现故障的主元方向；

步骤4.3：对第l类历史故障数据X_f，l进行基于霍特林统计量重构，计算第l类历史故障数据X_f，l的新的霍特林统计量的正常部分负载向量获得第l类历史故障数据重构后的霍特林统计量的正常部分E_p，l；

步骤4.4：对第l类历史故障数据X_f，l进行基于平方预测误差重构，计算第l类历史故障数据X_f，l的新的平方预测误差的正常部分负载向量获得第l类历史故障数据重构后的平方预测误差的正常部分E_e，l；

步骤4.5：将采样数据代入第l类历史故障数据重构后的霍特林统计量的正常部分，得到采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new，将采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new进行基于方向核偏最小二乘运算，得到采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new的主元T_pd，l，new＝[T_p，l，new，T_pr，l，new]；

步骤4.6：计算相对第l类故障数据重构后的采样数据的正常部分的霍特林统计量

步骤4.7：将采样数据代入第l类故障数据重构后的平方预测误差的正常部分，得到采样数据相对第l类故障数据重构后的平方预测误差的正常部分E_e，l，new；

步骤4.8：计算相对第l类故障数据重构后的采样数据的正常部分的平方预测误差SPE_e，l，new；

步骤4.9：当相对第l类故障数据重构后的采样数据的正常部分的霍特林统计量在第l类故障数据对应的霍特林统计量的控制限以下，同时相对第l类故障数据重构后的采样数据的正常部分的平方预测误差SPE_e，l，new在第l类故障数据对应的平方预测误差控制限以下时，则该采样数据的故障类型为第l类故障，否则，该采样数据的故障类型不是第l类故障，重新选择故障类型l，返回步骤4.2。

本发明的有益效果是：

本发明提出的基于方向核偏最小二乘的工业过程故障诊断方法，解决了传统方法中残差空间中存在质量相关变异以及变异量很大的问题，提高了对于非线性数据故障的检测能力，解决了非线性系统的故障诊断问题。结果表明，对于关注产品质量的非线性过程，基于方向核偏最小二乘的监测效果更好。基于方向核偏最小二乘故障重构方法能够有效得到故障数据的故障主元方法和故障方向，以及重构恢复后的正常数据，能够使得统计量超限现象快速消除。

附图说明

图1为本发明具体实施方式中的基于方向核偏最小二乘的工业过程故障诊断方法的流程图；

图2为本发明具体实施方式中的计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差SPE_d的控制限的流程图；

图3为本发明具体实施方式中的计算采样数据的过程监测统计量和采样数据的平方预测误差SPE_d，new的流程图；

图4为传统的KPLS运算方法计算的故障数据集A的输入变量的KPLS残差与输出变量的相关值；

图5为本发明具体实施方式计算的故障数据集A的输入变量的DKPLS残差与输出变量的相关值；

图6为传统的KPLS运算方法检测故障数据集A的检测结果；

其中，(a)为传统的KPLS运算方法检测故障数据集A的T²统计量，(b)为传统的KPLS运算方法检测故障数据集A的SPE统计量；

图7为本发明具体实施方式检测故障数据集A的检测结果；

其中，(a)为本发明具体实施方式检测故障数据集A的统计量，(b)为本发明具体实施方式检测故障数据集A的SPE_d统计量；

图8为传统的KPLS运算方法检测故障数据集B的检测结果；

其中，(a)为传统的KPLS运算方法检测故障数据集B的T²统计量，(b)为传统的KPLS运算方法检测故障数据集B的SPE统计量；

图9为本发明具体实施方式检测故障数据集B的检测结果；

其中，(a)为本发明具体实施方式检测故障数据集B的统计量，(b)为本发明具体实施方式检测故障数据集B的SPE统计量；

图10为本发明具体实施方式中的检测故障数据集A的主元T_f及各主元方向上的故障幅度；

图11为本发明具体实施方式中的检测故障数据集A的主元T_fr及各主元方向上的故障幅度；

图12为本发明具体实施方式中的检测故障数据集A的主元T_fvr及各主元方向上的故障幅度；

图13为本发明具体实施方式中采样数据相对故障数据集A重构后的统计量；

其中，(a)采样数据相对于霍特林统计量重构之后的故障数据集A的监测得到的统计量，(b)采样数据对相对于平方预测误差重构之后的故障数据集A的监测得到的SPE_E统计量；

图14为本发明具体实施方式中采样数据相对故障数据集B重构后的统计量；

其中，(a)采样数据相对于霍特林统计量重构之后的故障数据集B的监测得到的统计量，(b)采样数据对相对于平方预测误差重构之后的故障数据集B的监测得到的SPE_e统计量。

具体实施方式

下面结合附图对本发明具体实施方式加以详细的说明。

针对电熔镁炉熔炼过程中容易出现的故障和不良工况，选择对电熔镁炉的温度进行监控。炉内温度值是一个重要的参数，其值由电极内的电流值和电极的位置决定，因此将三个电极中其中一个的输入电压值、三相电流值、电极相对位置三个关键变量作为电熔镁炉熔炼过程的输入变量，将电熔镁炉熔炼过程中三个电极对应的炉温值作为过程模型的输出变量。

基于方向核偏最小二乘的工业过程故障诊断方法，如图1所示，包括以下步骤：

步骤1：获取易出故障的工业过程的输入变量和输出变量的历史正常数据，对该历史正常数据进行基于方向核偏最小二乘运算，令历史正常输入数据映射的高维特征空间F的PLS残差中有i个主元，求得该高维特征空间F的PLS残差中与输出变量相关的变异E_r及其主元T_r，则获得历史正常数据新的高维特征空间F的主元T_d＝[T，T_r]，计算历史正常数据的霍特林统计量和历史正常数据的平方预测误差SPE_d，并计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差SPE_d的控制限，如图2所示。

步骤1.1：获取易出故障的工业过程的输入变量和输出变量的历史正常数据：对m个历史正常数据进行n次采样，得到历史正常输入数据矩阵x和历史正常输出数据矩阵y。

历史正常输入数据矩阵x可以表示成一个n×m＝300×3阵如式(1)所示，历史正常输出数据矩阵y也相应的表示成一个n×m＝300×3矩阵如式(2)所示：

$x={\left[\begin{array}{ccc}{x}_{\mathrm{1,1}}& x_{\mathrm{1,2}}& x_{\mathrm{1,3}}\\ x_{\mathrm{2,1}}& x_{\mathrm{2,2}}& x_{\mathrm{2,3}}\\ .& .& .\\ .& .& .\\ .& .& .\\ x_{\mathrm{300,1}}& x_{\mathrm{300,2}}& x_{\mathrm{300,3}}\end{array}\right]}_{300\×3}---\left(1\right)$

$y={\left[\begin{array}{ccc}{y}_{\mathrm{1,1}}& y_{\mathrm{1,2}}& y_{\mathrm{1,3}}\\ y_{\mathrm{2,1}}& y_{\mathrm{2,2}}& y_{\mathrm{2,3}}\\ .& .& .\\ .& .& .\\ .& .& .\\ y_{\mathrm{300,1}}& y_{\mathrm{300,2}}& y_{\mathrm{300,3}}\end{array}\right]}_{300\×3}---\left(2\right)$

步骤1.2：对历史正常输入数据矩阵x和历史正常输出数据矩阵y进行标准化处理，得到预处理后的历史正常输入数据X和预处理后的历史正常输出数据Y。

对历史正常输入数据矩阵x和历史正常输出数据矩阵y进行标准化处理的过程如下：

定义x的第j个历史正常输入数据x_j表示为：x_j＝(x_1，j，x_2，j，…，x_300，j)^T，定义y的第j个历史正常输出数据y_j表示为：y_j＝(y_1，j，y_2，j，…，y_300，j)^T，j＝1，2，…，m＝1，2，3。

计算历史正常输入数据x_j的均值如式(3)所示，历史正常输出数据y_j的均值如式(4)所示：

${\stackrel{\‾}{x}}_j=E\left(x_j\right)=\frac{1}{n}\underset{\mathrm{\τ}=1}{\overset{n}{\mathrm{\Σ}}}{x}_{\mathrm{\τj}}=\frac{1}{300}\underset{\mathrm{\τ}=1}{\overset{10}{\mathrm{\Σ}}}{x}_{\mathrm{\τj}}---\left(3\right)$

${\stackrel{\‾}{y}}_j=E\left(y_j\right)=\frac{1}{n}\underset{\mathrm{\τ}=1}{\overset{n}{\mathrm{\Σ}}}{y}_{\mathrm{\τ},j}=\frac{1}{300}\underset{\mathrm{\τ}=1}{\overset{300}{\mathrm{\Σ}}}{y}_{\mathrm{\τ},j}---\left(4\right)$

其中，τ＝1，2，…，n＝1，2，…，300。

计算历史正常输入数据x_j的方差如式(5)所示，历史正常输出数据y_j的方差如式(6)所示：

$s_{x_j}^2=\mathrm{Var}\left(x_j\right)=\frac{1}{n-1}\underset{\mathrm{\τ}=1}{\overset{n}{\mathrm{\Σ}}}{(x_{\mathrm{\τ},j}-{\stackrel{\‾}{x}}_j)}^2---\left(5\right)$

$s_{y_j}^2=\mathrm{Var}\left(y_j\right)=\frac{1}{n-1}\underset{\mathrm{\τ}=1}{\overset{n}{\mathrm{\Σ}}}{(y_{\mathrm{\τ},j}-{\stackrel{\‾}{y}}_j)}^2---\left(6\right)$

利用历史正常输入数据x_j的方差对历史正常输入数据进行中心化处理得到X_τ，j如式(7)所示：

$X_{\mathrm{\τ},j}=\frac{x_{\mathrm{\τ},j}^{*}}{s_{x_j}}---\left(7\right)$

其中，为历史正常输入数据x_j的方差的开方。

利用历史正常输出数据y_j的方差对历史正常输出数据进行中心化处理得到Y_τ，j如式(8)所示：

$Y_{\mathrm{\τ},j}=\frac{y_{\mathrm{\τ},j}^{*}}{s_{y_j}}---\left(8\right)$

其中，为历史正常输出数据y_j的方差的开方。

因此，得到的预处理后的历史正常输入数据X和预处理后的历史正常输出数据Y如式(9)和式(10)所示：

$X={\left[\begin{array}{ccc}{X}_{\mathrm{1,1}}& X_{\mathrm{1,2}}& X_{\mathrm{1,3}}\\ X_{\mathrm{2,1}}& X_{\mathrm{2,2}}& X_{\mathrm{2,3}}\\ .& .& .\\ .& .& .\\ .& .& .\\ X_{\mathrm{300,1}}& X_{\mathrm{300,2}}& X_{\mathrm{300,3}}\end{array}\right]}_{300\×3}---\left(9\right)$

$Y={\left[\begin{array}{ccc}{Y}_{\mathrm{1,1}}& Y_{\mathrm{1,2}}& Y_{\mathrm{1,3}}\\ Y_{\mathrm{2,1}}& Y_{\mathrm{2,2}}& Y_{\mathrm{2,3}}\\ .& .& .\\ .& .& .\\ .& .& .\\ Y_{\mathrm{300,1}}& Y_{\mathrm{300,2}}& Y_{\mathrm{300,3}}\end{array}\right]}_{300\×3}---\left(10\right)$

其中，X_j＝(X_1，j，X_2，j，…，X_300，j)^T，Y_j＝(Y_1，j，Y_2，j，…，Y_300，j)^T。

步骤1.3：选取非线性变换φ(X)，将预处理后的历史正常输入数据X映射到历史正常输入数据的高维特征空间F，利用径向基内积核函数，求出历史正常输入数据的初始核矩阵K₁。

本实施方式中，根据Hilbert-Schmidt定理，核函数为k(X_τ，X_γ)＝φ(X_τ)φ(X_γ)，代表X_τ和X_γ到特征空间中的非线性映射的像的内积，其中X_τ为第τ次采样时历史正常输入数据的采样值，即X_τ＝(X_τ，1，X_τ，2，X_τ，3)，τ＝1，2，…，300，X_γ为第γ次采样时历史正常输入数据的采样值，即X_γ＝(X_γ，1，X_γ，2，X_γ，3)，γ＝1，2，…，300，利用径向基内积核函数求得历史正常输入数据的核矩阵如式(11)所示，求出历史正常输入数据的初始核矩阵K₁如式(12)所示：

$K_{1,\mathrm{\τ\γ}}^{\mathrm{raw}}=k(X_{\mathrm{\τ}},X_{\mathrm{\γ}})=\exp (-\frac{{\left|\right|X_{\mathrm{\τ}}-X_{\mathrm{\γ}}\left|\right|}^2}{c})---\left(11\right)$

$K_1=K_{1,\mathrm{\τ\γ}}^{\mathrm{raw}}-K_{1,\mathrm{\τ\γ}}^{\mathrm{raw}}I-{I_1{K}_{1,\mathrm{\τ\γ}}^{\mathrm{raw}}+I}_1{K}_{1,\mathrm{\τ\γ}}^{\mathrm{raw}}I---\left(12\right)$

其中，为历史正常输入数据的核矩阵，c为常数参数，I为单位矩阵， $I_1=\frac{1}{300}[1,...,1].$

步骤1.4：对预处理后的历史正常输入数据X和预处理后的历史正常输出数据Y进行KPLS运算，求得预处理后的历史正常输入数据X的PLS主元T、预处理后的历史正常输出数据Y的PLS主元U和经过α次迭代后使预处理后的历史正常输出数据Y的PLS主元U收敛的历史正常输出数据的核矩阵K_α+1

求得的预处理后的历史正常输入数据X的PLS主元T表示为T＝(t₁，t₂，…，t_m)，预处理后的历史正常输出数据Y的PLS主元U表示为U＝(u₁，u₂，…，u_m)，t_μ为X的第μ个得分，u_μ为Y的第μ个得分，经过α次迭代后使预处理后的历史正常输出数据Y的PLS主元U收敛的历史正常输出数据的核矩阵K_α+1如式(13)所示：

$K_{\mathrm{\α}+1}=(I-t_{{\mathrm{\μ}}_{\mathrm{\α}}}{t}_{{\mathrm{\μ}}_{\mathrm{\α}}}^T/t_{{\mathrm{\μ}}_{\mathrm{\α}}}^T{t}_{{\mathrm{\μ}}_{\mathrm{\α}}})K_{\mathrm{\α}}(I-t_{{\mathrm{\μ}}_{\mathrm{\α}}}{t}_{{\mathrm{\μ}}_{\mathrm{\α}}}^T/t_{{\mathrm{\μ}}_{\mathrm{\α}}}^T{t}_{{\mathrm{\μ}}_{\mathrm{\α}}}^T)---\left(13\right)$

其中，为α次迭代后的历史正常输入数据X的第μ个得分，K_α为经过α-1次迭代后的历史正常输入数据的核矩阵。

步骤1.5：令高维特征空间F的PLS残差中有i个主元，求得高维特征空间F的PLS残差中与输出变量相关的变异E_r及其主元T_r，则获得新的高维特征空间F的主元T_d＝[T，T_r]。

令高维特征空间F的PLS残差中有i个主元，求得高维特征空间F的PLS残差E中与输出变量相关的变异E_r如式(14)所示：

E_r＝Cφ_i(X)    (14)

其中，C＝Y(Y^TY)^-1Y^T，φ_i(X)为历史正常输入数据的高维特征空间φ(X)中具有i个主元的数据。

计算高维特征空间F的PLS残差中与输出变量相关的变异E_r的协方差矩阵S如式(15)所示：

S＝(1/n)φ_i(X)^TC^TCφ_i(X)    (15)

设P_r和λ是协方差矩阵S的特征向量和特征值，则有公式(16)：

SP_r＝λP_r    (16)

当λ≠0时，P_r可以看作是Cφ_i(X)的线性组合，则有P_r＝φ_i(X)^TC^TA，其中，A为(1/n)CK_α+1C^T的特征值。

因此高维特征空间F的残差中与输出变量相关的变异E_r的主元如式(17)所示：

T_r＝Cφ_i(X)P_r＝Cφ_i(X)φ_i(X)^TC^TA＝CK_α+1C^TA    (17)

获得新的高维特征空间F的主元为T_d＝[T，T_r]。

步骤1.6：计算历史正常数据的霍特林统计量和历史正常数据的平方预测误差SPE_d，并计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差PSE_d的控制限。

计算历史正常数据的霍特林统计量如式(18)所示：

$T_d^2=T_d{\mathrm{\Λ}}^{-1}{T}_d^T---\left(18\right)$

其中， $\mathrm{\Λ}=T_d^T{T}_d.$

计算历史正常数据的平方预测误差SPE_d如式(19)所示：

SPE_d＝||φ_i(X)-Cφ_i(X)||²＝||(I-C)φ_i(X)||²＝θ²K_α+1    (19)

其中，θ是(I-C)的特征值。

步骤2：采集工业过程的输入变量的采样数据，对该采样数据进行基于方向核偏最小二乘运算，令采样数据映射的高维特征空间F_new的PLS残差中有i个主元，求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new及其主元T_r，new，则获得采样数据新的高维特征空间F_new的主元T_d，new＝[T_new，T_r，new，计算采样数据的过程监测统计量和采样数据的平方预测误差SPE_d，new，如图3所示。

步骤2.1：采集工业过程的输入变量的采样数据：对m个采样数据进行n次采样，得到采样数据矩阵x_new，并对其进行标准化处理，得到预处理后的采样数据X_new。

得到预处理后的采样数据X_new如式(20)所示：

$X_{\mathrm{new}}={\left[\begin{array}{ccc}{X}_{\mathrm{new},\mathrm{1,1}}& X_{\mathrm{new},\mathrm{1,2}}& X_{\mathrm{new},\mathrm{1,3}}\\ X_{\mathrm{new},\mathrm{2,1}}& X_{\mathrm{new},\mathrm{2,2}}& X_{\mathrm{new},\mathrm{2,3}}\\ .& .& .\\ .& .& .\\ .& .& .\\ X_{\mathrm{new},\mathrm{300,1}}& X_{\mathrm{new},\mathrm{300,2}}& X_{\mathrm{new},\mathrm{300,3}}\end{array}\right]}_{300\×3}---\left(20\right)$

步骤2.2：利用非线性变换φ(X)将预处理后的采样数据X_new映射到采样数据的高维特征空间F_new，利用径向基内积核函数，求出采样数据的初始核矩阵K_new，1。

利用径向基内积核函数求得采样数据的核矩阵如式(21)所示，求出的采样数据的初始核矩阵K_new，1如式(22)所示：

$K_{\mathrm{new},\mathrm{\τ\γ}}^{\mathrm{raw}}=k(X_{\mathrm{new},\mathrm{\τ}},X_{\mathrm{new},\mathrm{\γ}})=\exp (-\frac{{\left|\right|X_{\mathrm{new},\mathrm{\τ}}-X_{\mathrm{new},\mathrm{\γ}}\left|\right|}^2}{c})---\left(21\right)$

$K_{\mathrm{new},1}=K_{\mathrm{new},\mathrm{\τ\γ}}^{\mathrm{raw}}-K_{\mathrm{new},\mathrm{\τ\γ}}^{\mathrm{raw}}I-{I_1{K}_{\mathrm{new},\mathrm{\τ\γ}}^{\mathrm{raw}}+I}_1{K}_{\mathrm{new},\mathrm{\τ\γ}}^{\mathrm{raw}}I---\left(22\right)$

其中，为采样数据的核矩阵，X_new，τ＝(X_new，τ，1，X_new，τ，2，X_new，τ，3)，X_new，γ＝(X_new，γ，1，X_new，γ，2，X_new，γ，3)^T。

步骤2.3：对预处理后的采样数据X_new进行KPLS运算，求得预处理后的采样数据X_new的主元T_new和经过α次迭代后的采样数据的核矩阵K_new，α+1。

求得的预处理后的采样数据X_new的主元T_new如式(23)所示：

T_new＝K_new，1U    (23)

经过α次迭代后的采样数据的核矩阵K_new，α+1如式(24)所示：

$K_{\mathrm{new},\mathrm{\α}+1}=(I-\frac{t_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}{t}_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}^T}{t_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}^T{t}_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}})K_{\mathrm{new},\mathrm{\α}}(I-\frac{t_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}{t}_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}^T}{t_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}^T{t}_{\mathrm{new},{\mathrm{\μ}}_{\mathrm{\α}}}})---\left(24\right)$

其中，为α次迭代后的预处理后的采样数据X_new的第μ个得分，K_new，α为经过α-1次迭代后的采样数据的核矩阵。

步骤2.4：令采样数据的高维特征空间F_new的PLS残差中有i个主元，求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new及其主元T_r，new，则获得采样数据的新的高维特征空间F_new的主元T_d，new＝[T_new，T_r，new]。

令采样数据的高维特征空间F_new的PLS残差中有i个主元，求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new如式(25)所示：

E_r，new＝Cφ_i(X_new)    (25)

其中，φ_i(X_new)为采样数据的高维特征空间φ(X_new)中具有i个主元的数据。

求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new的主元T_r，new如式(26)所示：

T_r，new＝CK_new，α+1C^TA′    (26)

其中，A′为(1/n)CK_new，α+1C^T的特征值。

获得采样数据的新的高维特征空间F_new的主元T_d，new＝[T_new，T_r，new]。

步骤2.5：计算采样数据的过程监测统计量和采样数据的平方预测误差SPE_d，new。

计算采样数据的过程监测统计量如式(27)所示：

$T_{d,\mathrm{new}}^2=T_{d,\mathrm{new}}{\mathrm{\Λ}}^{\′-1}{T}_{d,\mathrm{new}}^T---\left(27\right)$

其中， ${\mathrm{\Λ}}^{\′}=T_{d,\mathrm{new}}^T{T}_{d,\mathrm{new}}.$

计算采样数据的平方预测误差SPE_d，new如式(28)所示：

SPE_d，new＝||φ_i(X_new)-Cφ_i(X_new)||²＝||(I-C)φ_i(X_new)||²＝θ²K_new，α+1(X_new，X_new)    (28)

步骤3：当采样数据的过程监测统计量超出历史正常数据的过程监测统计量的控制限或者采样数据的平方预测误差SPE_d，new超出历史正常数据的平方预测误差SPE_d的控制限，则该采样数据中具有一种故障，执行步骤4，否则，将该采样数据视为正常数据。

步骤4：获取已知故障类型的历史故障数据，对已知故障类型的历史故障数据进行基于霍特林统计量重构和基于平方预测误差重构，判断采样数据的故障类型。

步骤4.1：获取已知故障类型的L种历史故障数据X_f，1，X_f，2，...，X_f，L。

步骤4.2：选取已知故障类型的L种历史故障数据中的第l类历史故障数据X_f，l，l＝1，2，...，L，将高维特征空间的历史正常输入数据φ(X)沿着高维特征空间的第l类历史故障数据φ(X_f，l)的故障方向进行重构，重构出高维特征空间的第l类历史故障数据φ(X_f，l)出现故障的主元方向。

将第l类历史故障数据X_f，l在其高维特征空间的映射φ(X_f，l)表示如式(29)所示：

φ(X_f，l)＝φ(X_f，l)^*+Δφ(X_f，l)＝φ(X_f，l)^*+∑_lf_l    (29)

其中，φ(X_f，l)^*表示第l类历史故障数据重构后的正常数据，Δφ(X_f，l)表示第l类历史故障数据与其重构后的正常数据的偏差，∑_l表示第l类历史故障数据出现故障的变量方向，f_l表示第l类历史故障数据中的故障幅值。

将高维特征空间的历史正常输入数据φ(X)沿着高维特征空间的第l类历史故障数据φ(X_f，l)的故障方向∑_l进行重构，则φ(X)可表示如式(30)所示：

φ(X)＝φ(X_f，l)^*+Δφ_z(X_f，l)＝φ(X_f，l)^*+∑_lf_zl    (30)

其中，Δφ_z(X_f，l)表示历史正常输入数据与第l类历史故障数据重构后的正常数据的偏差，f_zl表示历史正常输入数据沿φ(X_f，l)的故障方向∑_l的故障幅值。

将公式(29)和公式(30)左右两边都乘以高维特征空间的历史正常输入数据φ(X)的负载向量P，推导出公式(31)和公式(32)：

φ(X_f，l)P＝φ(X_f，l)^*P+Δφ(X_f，l)P＝T_f，l    (31)

φ(X)P＝φ(X_f，l)^*P+Δφ_z(X_f，l)P＝T    (32)

其中，T_f，l为高维特征空间的第l类历史故障数据的主元。

重构出高维特征空间的第l类历史故障数据φ(X_f，l)出现故障的主元方向，如式(33)所示：

T_f，l＝T+Δφ(X_f，1)P-Δφ_z(X_f，l)P＝T+∑_T，lf_T，l    (33)

其中，∑_T，l为高维特征空间的第l类历史故障数据φ(X_f，l)出现故障的主元方向，f_T，l为高维特征空间的第l类历史故障数据φ(X_f，l)在出现故障的主元方向上的故障幅值。

步骤4.3：对第l类历史故障数据X_f，l进行基于霍特林统计量重构，计算第l类历史故障数据X_f，l的新的霍特林统计量的正常部分负载向量获得第l类历史故障数据重构后的霍特林统计量的正常部分E_p，l。

根据基于方向核偏最小二乘运算方法，高维特征空间的历史正常输入数据φ(X)可以表示如式(34)所示：

$\mathrm{\φ}\left(X\right)=T_d{P}_d^T+E_{\mathrm{vr}}=[T,T_r]{[P,P_r]}^T+E_{\mathrm{vr}}---\left(34\right)$

其中，E_vr＝E-E_r。

将高维特征空间的第l类历史故障数据φ(X_f，l)利用高维特征空间的历史正常输入数据φ(X)的负载向量P及E_r的协方差矩阵的负载向量P_r投影到历史正常输入数据的空间，如式(35)所示：

$\mathrm{\φ}\left(X_{f,l}\right)=T_{\mathrm{fd},l}{P}_{\mathrm{fd},l}^T+E_{\mathrm{fvr},l}=[T_{f,l},T_{\mathrm{fr},l}{[P,P_r]}^T+E_{\mathrm{fvr},l}=T_{f,l}{P}^T+T_{\mathrm{fr},l}{P}_r^T+E_{\mathrm{fvr},l}---\left(35\right)$

其中，E_fvr，l为高维特征空间的第l类历史故障数据的残差中与输出变量相关的变异，T_fr，l为高维特征空间的第l类历史故障数据的残差中与输出变量相关的变异的主元，T_fd，l＝[T_f，l，T_fr，l]，P_fd，l为E_fvr，l的协方差矩阵的负载向量。

高维特征空间的第l类历史故障数据中能够使第l类历史故障数据的新的霍特林统计量超限部分φ(X_f，l)^*如式(36)所示：

$\mathrm{\φ}{\left(X_{f,l}\right)}^{*}=T_{f,l}^{*}{P}^{*T}+T_{\mathrm{fr},l}^{*}{P}_r^{*T}=T_{f,l}^{*}\frac{T^{*T}}{T^{*T}{T}^{*}}\mathrm{\φ}\left(X\right)+T_{\mathrm{fr},l}^{*}{A}^{*T}\mathrm{C\φ}\left(X\right)=B_l\mathrm{\φ}\left(X\right)---\left(36\right)$

其中，A^*为A对应的方向，为T_f，l中的故障主元方向，P^*为P对应的方向，为T_fr，l中的故障主元方向，为P_r对应的方向，T^*为T中的故障主元方向。

对高维特征空间的第l类历史故障数据中能够使第l类历史故障数据的新的霍特林统计量超限部分φ(X_f，l)^*进行KPCA运算，得到φ(X_f，l)^*的协方差S_f，l如式(37)所示：

S_f，l＝(1/n)φ(X)^TB_l ^TB_lφ(X)    (37)

则φ(X_f，l)*的主元T_p，l如式(38)所示：

T_p，l＝φ(X_f，l)^*P_p，l＝B_lφ(X)φ(X)^TB_l ^TA_f，l    (38)

其中，P_p，l＝φ(X)^TB_l ^TA_f，l为φ(X_f，l)^*的负载向量，A_f，l为的特征值。

令φ(X_f，l)^*的主元子空间的主元为φ(X_f，l)^*的主元子空间的主元的负载向量为φ(X_f，l)^*残差子空间的主元为φ(X_f，l)^*残差子空间的主元的负载向量为即为l类故障数据X_f，l的霍特林统计量的正常部分负载向量。

所以，第l类历史故障数据的霍特林统计量的正常部分负载向量如式(39)所示：

${\tilde{P}}_{p,l}=\mathrm{\φ}{\left(X_{f,l}\right)}^T{B}_l^T{\tilde{A}}_{f,l}---\left(39\right)$

其中，为的特征值。

高维特征空间的第l类历史故障数据中能够使第l类历史故障数据的新的霍特林统计量超限部分φ(X_f，l)^*可以表示如式(40)所示：

所以，根据公式(40)可得第l类历史故障数据重构后的霍特林统计量的正常部分E_p，l如式(41)所示：

$E_{p,l}=K_{f,\mathrm{\α}+1}{B}_l^T{\tilde{A}}_{f,l}{\tilde{A}}_{f,l}^T{B}_l\mathrm{\φ}\left(X\right)---\left(41\right)$

其中，K_f，α+1为经过α次迭代后的历史故障数据的核矩阵。

步骤4.4：对第l类历史故障数据X_f，l进行基于平方预测误差重构，计算第l类历史故障数据X_f，l的新的平方预测误差的正常部分负载向量获得第l类历史故障数据重构后的平方预测误差的正常部分E_e，l，

根据公式(34)将高维特征空间的第l类历史故障数据φ(X_f，l)利用高维特征空间的历史正常输入数据φ(X)的负载向量P及E_r的协方差矩阵的负载向量P_r投影到历史正常输入数据的空间，得到公式(42)：

$E_{\mathrm{fvr},l}=\mathrm{\φ}\left(X_{f,l}\right)-T_{f,l}{P}^T-T_{\mathrm{fr},l}{P}_r^T=\mathrm{\φ}\left(X_{f,l}\right)-{T_{f,l}}^T\frac{{T_{f,l}}^T}{{T_{f,l}}^T{T}_{f,l}}\mathrm{\φ}\left(X_{f,l}\right)-T_{\mathrm{fr},l}{A}^T\mathrm{C\φ}\left(X_{f,l}\right)=B_{\mathrm{fvr},l}\mathrm{\φ}\left(X_{f,l}\right)---\left(42\right)$

其中， $B_{\mathrm{fvr},l}=I-{T_{f,l}}^T\frac{{T_{f,l}}^T}{{T_{f,l}}^T{T}_{f,l}}-T_{\mathrm{fr},l}{A}^{T}C.$

高维特征空间的第l类历史故障数据中能够使第l类历史故障数据的新的平方预测误差SPE_fd，l超限部分如式(43)所示：

$E_{\mathrm{fvr},l}^{*}=T_{\mathrm{fvr},l}^{*}{P}_{\mathrm{fvr},l}^{*T}=E_{f,l}{K}_{f,l}{E}^T{A}_{\mathrm{vr}}^{*}{A}_{\mathrm{vr}}^{*T}\mathrm{E\φ}\left(X\right)=B_{f,l}\mathrm{\φ}\left(X\right)---\left(43\right)$

其中， $B_{f,l}=B_{\mathrm{fvr},l}{K}_{f,l}{B}^T{A}_{\mathrm{vr}}^{*}{A}_{\mathrm{vr}}^{*T}B,B=I-T^T\frac{T^T}{T^{T}T}-T_r{A}^{T}C,$ 为高维特征空间的第l类历史故障数据φ(X_f，l)的E_fvr，l的主元T_fvr，l的故障主元方向，为高维特征空间的第l类历史故障数据φ(X_f，l)的E_fvr，l的主元负载向量P_fvr，l的方向，A_vr为的特征值，为A_vr的方向。

对高维特征空间的第l类历史故障数据中能够使第l类历史故障数据的新的平方预测误差SPE_fd，l超限部分进行KPCA计算，得到的主元T_e，l如式(44)所示：

$T_{e,l}=E_{\mathrm{fvr},l}^{*}{P}_{e,l}=B_{f,l}{K}_{\mathrm{\α}+1}{B}_{f,l}^T{A}_{\mathrm{ef},l}---\left(44\right)$

其中， $P_{e,l}=\mathrm{\φ}{\left(X\right)}^T{B}_{f,l}^T{A}_{\mathrm{ef},1}$ 为负载向量，A_ef，l为的特征值。

令的主元子空间的主元为的主元子空间的主元的负载向量为的残差子空间的主元为残差子空间的主元的的负载向量为即为第l类历史故障数据X_f，l的新的平方预测误差的正常部分负载向量。

所以，第l类历史故障数据X_f，l的新的平方预测误差的正常部分负载向量如式(45)所示：

${\tilde{P}}_{e,l}=\mathrm{\φ}{\left(X\right)}^T{B}_{f,l}^T{\tilde{A}}_{\mathrm{ef},l}---\left(45\right)$

其中，为的特征值。

因此，高维特征空间的第l类历史故障数据中能够使第l类历史故障数据的新的平方预测误差SPE_fd，l超限部分可以表示如式(46)所示：

根据公式(46)可得第l类历史故障数据重构后的平方预测误差的正常部分E_e，l如式(47)所示：

$E_{e,l}=K_{f,\mathrm{\α}+1}{B}_{f,l}^T{\tilde{A}}_{\mathrm{ef},l}{\tilde{A}}_{\mathrm{ef},l}^T{B}_{f,l}\mathrm{\φ}\left(X\right)---\left(47\right)$

步骤4.5：将采样数据代入第l类历史故障数据重构后的霍特林统计量的正常部分，得到采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new，将采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new进行基于方向核偏最小二乘运算，得到采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new的主元T_pd，l，new＝[T_p，l，new，T_pr，l，new]。

得到的采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new如式(48)所示：

$E_{p,l,\mathrm{new}}=\mathrm{\φ}\left(X_{\mathrm{new}}\right){\tilde{P}}_{p,l}{\tilde{P}}_{p,l}^T=\mathrm{\φ}\left(X_{\mathrm{new}}\right)\mathrm{\φ}{\left(X\right)}^T{B}^T{\tilde{A}}_{f,l}{\tilde{A}}_{f,l}^T\mathrm{B\φ}\left(X\right)=K_{\mathrm{new},\mathrm{\α}+1}{B}^T{\tilde{A}}_{f,l}{\tilde{A}}_{f,l}^T\mathrm{B\φ}\left(X\right)---\left(48\right)$

得到采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new的主元如式(49)、式(50)和式(51)所示：

T_pd，l，new＝[T_p，l，new，T_pr，new]    (49)

$T_{p,l,\mathrm{new}}=E_{p,l,\mathrm{new}}\mathrm{\φ}{\left(X\right)}^{T}U=K_{\mathrm{new},\mathrm{\α}+1}{B}_l^T{\tilde{A}}_{f,1}{\tilde{A}}_{f,1}^T{B}_l\mathrm{\φ}\left(X\right)\mathrm{\φ}{\left(X\right)}^{T}U=K_{\mathrm{new},\mathrm{\α}+1}{B}_l^T{\tilde{A}}_{f,1}{\tilde{A}}_{f,1}^T{B}_l\mathrm{KU}---\left(50\right)$

$T_{p,l,\mathrm{new}}=E_{p,l,\mathrm{new}}{P}_r=K_{\mathrm{new},\mathrm{\α}+1}{B}_l^T{\tilde{A}}_{f,1}{\tilde{A}}_{f,l}^T{B}_l\mathrm{\φ}\left(X\right)\mathrm{\φ}{\left(X\right)}^T{C}^{T}A=K_{\mathrm{new},\mathrm{\α}+1}{B}_l^T{\tilde{A}}_{f,1}{\tilde{A}}_{f,1}^T{B}_l{K}_{\mathrm{\α}+1}{C}^{T}A---\left(51\right)$

步骤4.6：计算相对第l类故障数据重构后的采样数据的正常部分的霍特林统计量

计算相对第l类故障数据重构后的采样数据的正常部分的霍特林统计量如式(52)所示：

步骤4.7：将采样数据代入第l类故障数据重构后的平方预测误差的正常部分，得到采样数据相对第l类故障数据重构后的平方预测误差的正常部分E_e，l，new。

采样数据相对第l类故障数据重构后的平方预测误差的正常部分E_e，l，new如式(53)所示：

$E_{e,l,\mathrm{new}}=\mathrm{\φ}\left(X_{\mathrm{new}}\right){\tilde{P}}_{e,l}{P}_{e,l}^T=\mathrm{\φ}\left(X_{\mathrm{new}}\right)\mathrm{\φ}{\left(X\right)}^T{B}_{f,l}^T{\tilde{A}}_{\mathrm{ef},l}{\tilde{A}}_{\mathrm{ef},l}^T{B}_{f,l}\mathrm{\φ}\left(X\right)=K_{\mathrm{new},\mathrm{\α}+1}{B}_{f,l}^T{\tilde{A}}_{\mathrm{ef},l}{\tilde{A}}_{\mathrm{ef},l}^T{B}_{f,l}\mathrm{\φ}\left(X\right)---\left(53\right)$

步骤4.8：计算相对第l类故障数据重构后的采样数据的正常部分的平方预测误差SPE_e，l，new。

相对第l类故障数据重构后的采样数据的正常部分的平方预测误差SPE_e，l，new如式(54)所示：

${\mathrm{SPE}}_{e,l,\mathrm{new}}={\left|\right|E_{e,l,\mathrm{new}}\left|\right|}^2=K_{\mathrm{new},\mathrm{\α}+1}{B}_{f,l}^T{\tilde{A}}_{\mathrm{ef},l}{\tilde{A}}_{\mathrm{ef},l}^T{B}_{f,l}{K}_{\mathrm{\α}+1}{B}_{f,l}^T{\tilde{A}}_{\mathrm{ef},l}{\tilde{A}}_{\mathrm{ef},l}^T{B}_{f,l}{K}_{\mathrm{new},\mathrm{\α}+1}^T---\left(54\right)$

步骤4.9：当相对第l类故障数据重构后的采样数据的正常部分的霍特林统计量在第l类故障对应的霍特林统计量的控制限以下，同时相对第l类故障数据重构后的采样数据的正常部分的平方预测误差SPE_e，l，new在第l类故障数据对应的平方预测误差控制限以下时，则该采样数据的故障类型为第l类故障，否则，该采样数据的故障类型不是第l类故障，重新选择故障类型l，返回步骤4.2。

本实施方式中，选取已知故障类型为A的300个故障数据集和已知故障类型为B的300个故障数据集进行上述步骤的工业过程故障监测与诊断，分别对比了传统的KPLS运算方法和基于方向核偏最小二乘的工业过程故障诊断方法的故障诊断过程。

已知故障类型为A的300个故障数据集大约从第50个采样数据开始发生，到大约第150个采样数据时结束，故障是由于变压器异常导致电熔镁炉电流大幅下降，出现炉温异常。

已知故障类型为B的300个故障数据集大约从第150个采样数据开始发生，到大约第300个采样数据时结束，故障是由于电极执行器异常导致电熔镁炉电流缓慢变化，造成炉温异常。

由图4和图5可以看出，对已知故障类型为A的300个故障数据集进行不同方法的测试，利用本发明基于方向核偏最小二乘的工业过程故障诊断方法计算的输入数据的DKPLS(基于方向核偏最小二乘)残差与输出数据的相关值小于利用传统的KPLS运算方法计算的输入数据的KPLS残差与输出数据的相关值，表明本发明已经提取出了KPLS残差中的输出相关变异。

传统的KPLS运算方法检测故障数据集A的检测结果，如图6所示，图6(a)中，传统的KPLS运算方法的T²统计量并未显示过程出现故障，大部分T²统计量数值都在控制限之下，只有第140，205，245和265个采样处的T²统计量数值超出了其控制限，可以判断这四个超限T²统计量为误报，这表明了传统的KPLS运算方法的T²统计量不能够检测出故障数据集A中添加的故障。图6(b)中，传统的KPLS运算方法检测故障数据集A的SPE统计量检测到了过程内部产生的故障。传统的KPLS运算方法的SPE统计量数值从大约第51个采样点开始超过控制限，一直持续到第150个采样，这与故障的时间相符，表明传统的KPLS运算方法的SPE统计量能够检测到故障数据集A中添加的故障。除此之外，传统的KPLS运算方法的SPE统计量也在第205，245和265个采样处出现超限现象，与传统的KPLS运算方法的T²统计量超限处一致，表明是过程内部存在的特性变化造成了传统KPLS的T²和SPE统计量的误报。

本发明具体实施方式检测故障数据集A的检测结果，如图7所示，图7(a)中，本发明具体实施方式的统计量显示出了过程内部出现的故障，从大约第50个采样开始，持续到大约第150个采样，其统计量一直大于控制限。统计量检测得到的故障与实际加入的故障时间相符，表明本发明具体实施方式的统计量能够检测出故障数据集A中添加的故障。图7(b)中，本发明具体实施方式的SPE_d统计量也检测到了过程内部产生的故障。本发明具体实施方式的SPE_d统计量数值从大约第51个采样点开始超过控制限，一直持续到第150个采样，表明DKPLS的SPE_d统计量能够检测到故障数据集A中添加的故障。对于故障数据集A中的故障来说，本发明的效果优于传统的KPLS方法。传统的KPLS方法中，只有SPE统计量能够检测到故障，而本发明的方法的和SPE_d均能够检测到故障。

传统的KPLS运算方法检测故障数据集B的检测结果，以及本发明具体实施方式检测故障数据集B的检测结果，如图8、图9所示，对于故障数据集B中的故障来说，本发明的效果优于传统的KPLS方法，与故障数据集A的检测结果相似，使用传统的KPLS方法，只有SPE统计量能够检测到故障，而使用本发明的方法，其和SPE_d统计量均能够检测到故障。

对故障数据集A的故障方向进行重构，提取故障方向，恢复出相对于监测统计量的正常数据。故障数据集A的数据向正常数据负载投影得到主元子空间的主元T_f和T_fr中，各主元方向的故障幅度如图10和图11所示。图10中，T_f的第三主元方向的故障幅度远远大于其他方向，因此选择T_f的第三主元方向作为其故障主元方向。图11中，T_fr的第三主元方向的故障幅度比其他主元方向的故障幅度大很多，因此选择T_fr的第三主元方向作为其故障主元方向。故障数据集A的数据向正常数据残差子空间负载投影得到残差空间的主元T_fvr中，各主元方向的故障幅度如图12所示。图12中，T_fvr的第三主元方向的故障幅度最大，第四主元方向的故障幅度次之，两个主元方向的故障幅度都要远远大于其他主元方向的故障幅度，因此选取T_fvr的第三和第四主元方向作为其故障主元方向。

图13(a)中，是采样数据相对于霍特林统计量重构之后的故障数据集A的监测得到的统计量，从图中可以看到，重构之后统计量消除了超限现象，全部处于控制限以下。图13(b)中，是采样数据对相对于平方预测误差重构之后的故障数据集A的监测得到的SPE_e统计量，与统计量类似，其超限现象也被消除，位于控制限以下。可以看出采样数据重构后的统计量全部位于其控制限以下，这样便可以将监测得到的采样数据的故障诊断为故障类型A。

图14(a)中，是采样数据相对于霍特林统计量重构之后的故障数据集B的监测得到的统计量，从图中可以看到，重构之后统计量与重构之前的统计量图形基本一致，超限现象并未消除，图14(b)中，是采样数据对相对于平方预测误差重构之后的故障数据集A的监测得到的SPE_e统计量，与重构之前的SPE_p统计量相似，其超限现象也没有消除。可以看出重构后的统计量超限现象并未消除，表明采样数据并非为故障类型B。

Claims (3)

1.基于方向核偏最小二乘的工业过程故障诊断方法，其特征在于，包括以下步骤：

步骤1：获取易出故障的工业过程的输入变量和输出变量的历史正常数据，对该历史正常数据进行基于方向核偏最小二乘运算，令历史正常输入数据映射的高维特征空间F的PLS残差中有i个主元，求得该高维特征空间F的PLS残差中与输出变量相关的变异E_r及其主元T_r，则获得历史正常数据新的高维特征空间F的主元T_d＝[T，T_r]，计算历史正常数据的霍特林统计量和历史正常数据的平方预测误差SPE_d，并计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差SPE_d的控制限；

步骤2：采集工业过程的输入变量的采样数据，对该采样数据进行基于方向核偏最小二乘运算，令采样数据映射的高维特征空间F_new的PLS残差中有i个主元，求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new及其主元T_r，new，则获得采样数据新的高维特征空间F_new的主元T_d，new＝[T_new，T_r，new]，计算采样数据的过程监测统计量和采样数据的平方预测误差SPE_d，new；

步骤3：当采样数据的过程监测统计量超出历史正常数据的过程监测统计量的控制限或者采样数据的平方预测误差SPE_d，new超出历史正常数据的平方预测误差SPE_d的控制限，则该采样数据中具有一种故障，执行步骤4，否则，将该采样数据视为正常数据；

步骤4：获取已知故障类型的历史故障数据，对已知故障类型的历史故障数据进行基于霍特林统计量重构和基于平方预测误差重构，判断采样数据的故障类型；

步骤4.1：获取已知故障类型的L种历史故障数据X_f，1，X_f，2，…，X_f，L；

步骤4.2：选取已知故障类型的L种历史故障数据中的第l类历史故障数据X_f，l，l＝1，2，…，L主将高维特征空间的历史正常输入数据φ(X)沿着高维特征空间的第l类历史故障数据φ(X_f，l)的故障方向进行重构，重构出高维特征空间的第l类历史故障数据φ(X_f，l)出现故障的主元方向；

步骤4.3：对第l类历史故障数据X_f，l进行基于霍特林统计量重构，计算第l类历史故障数据X_f，l的新的霍特林统计量的正常部分负载向量获得第l类历史故障数据重构后的霍特林统计量的正常部分E_p，l；

步骤4.4：对第l类历史故障数据X_f，l进行基于平方预测误差重构，计算第l类历史故障数据X_f，l的新的平方预测误差的正常部分负载向量获得第l类历史故障数据重构后的平方预测误差的正常部分E_e，l；

步骤4.5：将采样数据代入第l类历史故障数据重构后的霍特林统计量的正常部分，得到采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new，将采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new进行基于方向核偏最小二乘运算，得到采样数据相对第l类故障数据重构后的霍特林统计量的正常部分E_p，l，new的主元T_pd，l，new＝[T_p，l，new，T_pr，l，new]；

步骤4.6：计算相对第l类故障数据重构后的采样数据的正常部分的霍特林统计量

步骤4.7：将采样数据代入第l类故障数据重构后的平方预测误差的正常部分，得到采样数据相对第l类故障数据重构后的平方预测误差的正常部分E_e，l，new；

步骤4.8：计算相对第l类故障数据重构后的采样数据的正常部分的平方预测误差SPE_e，l，new；

步骤4.9：当相对第l类故障数据重构后的采样数据的正常部分的霍特林统计量在第l类故障数据对应的霍特林统计量的控制限以下，同时相对第l类故障数据重构后的采样数据的正常部分的平方预测误差SPE_e，l，new在第l类故障数据对应的平方预测误差控制限以下时，则该采样数据的故障类型为第l类故障，否则，该采样数据的故障类型不是第l类故障，重新选择故障类型l，返回步骤4.2。

2.根据权利要求1所述的基于方向核偏最小二乘的工业过程故障诊断方法，其特征在于，所述的步骤1包括以下步骤：

步骤1.1：获取易出故障的工业过程的输入变量和输出变量的历史正常数据：对m个历史正常数据进行n次采样，得到历史正常输入数据矩阵x和历史正常输出数据矩阵y；

步骤1.2：对历史正常输入数据矩阵x和历史正常输出数据矩阵y进行标准化处理，得到预处理后的历史正常输入数据X和预处理后的历史正常输出数据Y；

步骤1.3：选取非线性变换φ(X)，将预处理后的历史正常输入数据X映射到历史正常输入数据的高维特征空间F，利用径向基内积核函数，求出历史正常输入数据的初始核矩阵K₁；

步骤1.4：对预处理后的历史正常输入数据X和预处理后的历史正常输出数据Y进行KPLS运算，求得预处理后的历史正常输入数据X的PLS主元T、预处理后的历史正常输出数据Y的PLS主元U和经过α次迭代后使预处理后的历史正常输出数据Y的PLS主元U收敛的历史正常输出数据的核矩阵K_α+1；

步骤1.5：令高维特征空间F的PLS残差中有i个主元，求得高维特征空间F的PLS残差中与输出变量相关的变异E_r及其主元T_r，则获得新的高维特征空间F的主元T_d＝[T，T_r]；

步骤1.6：计算历史正常数据的霍特林统计量和历史正常数据的平方预测误差SPE_d，并计算历史正常数据的霍特林统计量的控制限和历史正常数据的平方预测误差SPE_d的控制限。

3.根据权利要求1所述的基于方向核偏最小二乘的工业过程故障诊断方法，其特征在于，所述的步骤2包括以下步骤：

步骤2.1：采集工业过程的输入变量的采样数据：对m个采样数据进行n次采样，得到采样数据矩阵x_new主并对其进行标准化处理，得到预处理后的采样数据X_new；

步骤2.2：利用非线性变换φ(X)将预处理后的采样数据X_new映射到采样数据的高维特征空间F_new主利用径向基内积核函数，求出采样数据的初始核矩阵K_new，1；

步骤2.3：对预处理后的采样数据X_new进行KPLS运算，求得预处理后的采样数据X_new的主元T_new和经过α次迭代后的采样数据的核矩阵K_new，α+1；

步骤2.4：令采样数据的高维特征空间F_new的PLS残差中有i个主元，求得该高维特征空间F_new的PLS残差中与输出变量相关的变异E_r，new及其主元T_r，new，获得采样数据的新的高维特征空间F_new的主元T_d，new＝[T_new，T_r，new]；

步骤2.5：计算采样数据的过程监测统计量和采样数据的平方预测误差SPE_d，new。