CN108445867B - 一种基于分散式icr模型的非高斯过程监测方法 - Google Patents

Abstract

本发明公开一种基于分散式ICR模型的非高斯过程监测方法，旨在解决如何利用非高斯数据建模算法，通过数据模型将采样数据转换成误差，并以误差作为被监测对象实施非高斯过程监测的问题。具体来讲，本发明方法首先针对每个测量变量，利用独立成分回归(ICR)算法建立各变量与其他变量之间的软测量模型。然后，利用软测量模型的估计误差作为被监测对象，建立基于独立成分分析(ICA)的过程监测模型实施非高斯过程监测。可以看出本发明方法利用了分散式建模的优势，而且采用多种非高斯数据分析算法相结合的实施方式，是一种更为优选的适于非高斯过程的数据驱动的过程监测方法。

Description

一种基于分散式ICR模型的非高斯过程监测方法

技术领域

本发明涉及一种数据驱动的过程监测方法，尤其涉及一种基于分散式ICR模型的非高斯过程监测方法。

背景技术

实时地监测生产过程的运行状态是保证安全生产与维持产品质量稳定的直接途径，对于过程监测的研究一直都存在于整个工业的发展历程中，其目的在于及时而准确的发现故障。当前，故障检测的主流技术手段是数据驱动的过程监测方法。现代化工过程的大型化建设以及先进仪表与计算机技术的广泛应用，生产过程可以采集海量的数据，这为数据驱动的过程监测研究提供了坚实的数据基础。而传统基于机理模型的故障检测方法现已逐渐没落，主要是它需要使用较为精确的机理模型从而产生过程某些变量或参数的估计误差。相比之下，数据驱动的故障检测方法不需要机理模型只需要采样数据，比较适合于现代工业过程运行状态的监测。一般来讲，数据驱动的过程监测方法建模的基本出发点都是对正常工况下的采样数据进行特征挖掘，通过投影变换的方式将原始采样数据转换成潜在特征。实施在线故障检测时，将在线采样数据经同样的投影变换得到相应的潜在特征，然后依据潜在特征的统计量指标判断是否为故障样本。

由此可见，数据驱动的过程监测方法与基于机理模型的故障检测方法是存在显著差异的，前者旨在挖掘潜在特征二后者旨在生成误差。发展至今，数据驱动的过程监测研究领域涌现出了很多特征挖掘算法以及各式各样的建模思路。就特征挖掘算法而言，主成分分析(Principal Component Analysis，PCA)、独立成分分析(Independent ComponentAnalysis，ICA)、偏最小二乘(Partial Least Squares，PLS)、以及它们的非线性扩展形式是被广泛研究的算法。就建模思路而言，处理数据自相关性科采用增广矩阵法，对大规模过程对象可采用多变量块的分散式建模策略。此外，还有多阶段建模、即使学习、集成学习等等思路。可以说，数据驱动的故障检测方法已如雨后春笋般蓬勃发展。值得一提的是，ICA作为一种典型的非高斯数据分析算法，已被广泛应用于非高斯过程的数据建模。例如，基于ICA的非高斯过程监测方法、将ICA与最小二乘法结合得到非高斯过程回归的独立成分回归(Independent Component Regression，ICR)、等等。由于非高斯分布的采样数据更为常见，因此ICA或ICR的方法实用性与应用性更为广泛。

然而，数据驱动的故障检测方法似乎已于基于机理模型的故障检测方法越走越远，鲜有在数据驱动模型中考虑误差生成的策略。其实，误差变化的大小能很好地反映过程对象是否处于正常运行状态，而且误差一般都服从或近似服从高斯分布。可是精确机理模型的建立对于现代工业过程来讲可谓是困难重重，那么可否用数据模型替代机理模型后再生成过程对象中某些参数货变量的估计误差呢？由于机理模型建立的是过程输入与输出之间的关系，而数据驱动的过程监测模型旨在特征挖掘，不存在区分过程输入与输出。任何一个测量变量在故障检测中的地位是等同，也不好将某些变量作为输出而其他变量作为输出。因此，在数据驱动的故障检测模型中实施类似的误差生成方案还有待进一步商榷。可想而知，若是能实现这一目的，那么不仅能拉近数据驱动的故障检测方法与基于机理模型的故障检测方法之间的距离，而且还能重新定义现有数据驱动故障检测的方法框架。

发明内容

本发明所要解决的主要技术问题是：如何利用非高斯数据建模算法，通过数据模型将采样数据转换成误差，并以误差作为被监测对象实施非高斯过程监测。具体来讲，本发明方法首先针对每个测量变量，利用ICR算法建立各变量与其他变量之间的软测量模型。然后，利用软测量模型的估计误差作为被监测对象，建立基于ICA的过程监测模型实施在线非高斯过程监测。

本发明解决上述技术问题所采用的技术方案为：一种基于分散式ICR模型的非高斯过程监测方法，包括以下步骤：

(1)采集生产过程正常运行状态下的样本，组成训练数据集X∈R^n×m，并对每个测量变量进行标准化处理，得到均值为0，标准差为1的新数据矩阵

其中，n为训练样本数，m为过程测量变量数，R为实数集，R^n×m表示n×m维的实数矩阵。

(2)将标准化后的数据集

表示成

其中x_i∈R^n×1为第i个测量变量的n个数据组成的列向量，i＝1，2，…，m，并初始化i＝1。

(3)将矩阵

中第i列数据去除得到软测量模型的输入矩阵

而将x_i作为软测量模型的输出，利用独立成分回归(Independent Component Regression，ICR)算法建立输入矩阵

与输出x_i之间的软测量模型，具体的实施过程如下所示：

①设置k＝1后，初始化向量u_k＝x_i；

②根据公式

q_k＝x_i ^Ts_k/(s_k ^Ts_k)分别计算得到系数向量w_k、得分向量s_k、和系数q_k；

③根据公式u_new＝x_iq_k/q_k ²计算向量u_new；

④判断是否满足条件||u_k-u_new||＜10^-6；若否，则置u_k＝u_new后返回步骤②；若是，则执行下一步骤⑤；

⑤根据公式

计算得到第k个投影向量p_k∈R^(m-1)×1，并保留投影向量p_k、系数向量w_k、和系数q_k；

⑥根据公式

更新输入矩阵

后，判断矩阵

中的最大元素是否大于0.001；若是，则执行步骤⑦；若否，则在得到最终的投影矩阵P_i＝[p₁，p₂，…，p_k]与系数矩阵W_i＝[w₁，w₂，…，w_k]后，执行步骤⑧；

⑦判断k＜m-1；若是，则置k＝k+1后，返回步骤②；若否，则得到最终的投影矩阵P_i＝[p₁，p₂，…，p_k]与系数矩阵W_i＝[w₁，w₂，…，w_k]；

⑧根据公式

对

进行白化处理，得到白化后的矩阵Z_i∈R^n×k并初始化j＝1，其中，矩阵R_i＝W_i(P_i ^TW_i)^-1∈R^(m-1)×k、矩阵

⑨设置向量c_j为k×k维单位矩阵中的第j列；

⑩按照公式c_j←E{Z_ig(c_j ^TZ_i)}-E{h(c_j ^TZ_i)}c_j更新向量c_j，其中函数g(u)＝tanh(u)、函数h(u)＝[sech(u)]²、u＝c_j ^TZ_i表示函数自变量；

对更新后的向量c_j依次按照下式进行正交标准化处理2

c_j←c_j/||c_j|| (5)

重复步骤⑤～⑥直至向量c_j收敛，并保存向量c_j；

判断j＜k；若是，置j＝j+1后，重复步骤④～⑧；若否，则执行步骤⑨；

将得到的所有k个向量c₁，c₂，…，c_k组成矩阵C＝[c₁，c₂，…，c_k]∈R^k×k，并按照如下所示公式计算分离矩阵W₀∈R^(m-1)×k与混合矩阵A₀∈R^(m-1)×k：

A₀＝P_iD_i ^1/2C

W₀＝P_iD_i ^-1/2C (6)

根据公式

计算独立成分矩阵S₀∈R^n×k，并根据公式

计算独立成分与输出x_i之间的相关系数向量

找出行向量

中数值小于0.001的所有元素所对应的列标号，根据该列标号对应的剔除分离矩阵W₀与独立成分矩阵S₀中相应的列，从而得到新的分离矩阵W_i与独立成分矩阵S_i；

计算最小二乘回归系数向量b_i＝W_i(S_i ^TS_i)^-1S_i ^Tx_i，那么对应于输出x_i的ICR模型可表示成：

上式中，e_i∈R^n×1为软测量估计误差向量。

(4)判断是否满足条件i＜m；若是，则置i＝i+1后返回步骤(3)；若否，则将得到的估计误差向量组成矩阵E＝[e₁，e₂，…，e_m]后继续执行下一步骤(5)。

(5)将矩阵E作为新的训练数据矩阵，对其中的每一列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

(6)利用独立成分分析(Independent Component Analysis，ICA)算法为

建立相应的ICA过程监测模型，并保留模型参数集Θ＝{G，H，Λ，D_lim，Q_lim}以备调用，其中G为分离矩阵，H为混合矩阵，Λ为独立成分的协方差矩阵，D_lim与Q_lim分别表示监测统计量D与Q的控制上限，具体的实施过程如下所示：

①计算

的协方差矩阵

②求解Ψ所有M个大于0.001的特征值η₁≥η₂≥…≥η_M所对应的特征向量β₁，β₂…，β_M；

③根据公式

对

进行白化处理，得到白化后的矩阵F∈R^n×M并初始化i＝1，其中，矩阵B＝[β₁，β₂…，β_M]∈R^m×M、对角矩阵Φ＝diag{η₁，η₂，…，η_M}；

④设置向量d_i为M×M维单位矩阵中的第i列；

⑤按照公式d_i←E{Fg(d_i ^TF)}-E{h(d_i ^TF)}d_i更新向量d_i，其中函数g(u)＝tanh(u)、函数h(u)＝[sech(u)]²、u＝d_i ^TF表示函数自变量；

⑥对更新后的向量d_i依次按照下式进行正交标准化处理：

d_i←d_i/||d_i|| (6)

⑦重复步骤⑤～⑥直至向量d_i收敛，并保存向量d_i；

⑧判断i＜M；若是，置i＝i+1后，重复步骤④～⑧；若否，则执行步骤⑨；

⑨将得到的所有M个向量c₁，c₂，…，c_M组成矩阵Ω＝[d₁，d₂，…，d_K]∈R^M×M，并按照如下所示公式计算分离矩阵G₀∈R^m×M与混合矩阵H₀∈R^m×M：

G₀＝BΦ^-1/2Ω

H₀＝BΦ^1/2Ω (7)

⑩计算H₀中每一列向量的长度，分别记为L₁，L₂，…，L_M，并将L₁，L₂，…，L_M按照数值大小进行降序排列得到l₁，l₂，…，l_M，那么保留的独立成分个数μ为满足下列条件的最小值：

将H₀中列向量长度最大的μ个列向量组成混合矩阵H∈R^m×μ，同时从G₀中取出与H对应的列向量组成分离矩阵G∈R^m×μ；

根据公式

计算出独立成分矩阵J后，计算J的协方差矩阵Λ＝(J^TJ)/(n-1)；

根据如下所示公式计算D_lim与Q_lim：

上两式中，置信水平α＝99％，

表示权重为g＝v/2a，自由度为h＝2a²/v的χ²分布，a与v分别是

对应的监测统计量Q的估计均值和估计方差。

(7)收集新采样时刻的数据样本x∈R^1×m，对其实施与步骤(1)中相同的标准化处理得到新数据向量

后，初始化i＝1。

(8)将行向量

中的第i个元素y_i取出后得到输入向量

并按照如下所示公式计算y_i的估计误差f_i：

(9)判断是否满足条件i＜m；若是，则置i＝i+1后返回步骤(8)；若否，则将得到的误差组成向量f＝[f₁，f₂，…，f_m]并继续执行下一步骤(10)。

(10)对f实施与步骤(5)中相同的标准化处理得到新向量

(11)根据如下所示公式计算监测统计量D与Q：

(12)判断D与Q的具体数值是否大于对应控制上限D_lim与Q_lim；若否，则当前样本为正常工况采样；若是，则当前采样数据来自故障工况

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

首先，本发明方法利用分散式建模的策略为每个测量变量建立其各自的非高斯软测量模型，不仅等同的对待了所有的测量变量而且还利用ICR算法应对非高斯数据建模。其次，本发明方法监测的对象为估计误差而不是原始数据或者特征成分，且进一步对误差实施基于ICA的非高斯过程监测。因此，本发明方法能显著提升传统故障检测方法的故障检测能力，是一种更为优选的故障检测方法。

附图说明

图1为本发明方法基本原理示意图。

图2为TE过程故障19的故障检测详情对比图。

具体实施方式

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

如图1所示，一种基于分散式主元回归模型的故障检测方法。下面结合一个具体的工业过程的例子来说明本发明方法的具体实施过程，以及相对于现有方法的优越性。

应用对象是来自于美国田纳西-伊斯曼(TE)化工过程实验，原型是伊斯曼化工生产车间的一个实际工艺流程。目前，TE过程因其流程的复杂性，已作为一个标准实验平台被广泛用于故障检测研究。整个TE过程包括22个测量变量、12个操作变量、和19个成分测量变量。所采集的数据分为22组，其中包括1组正常工况下的数据集与21组故障数据。而在这些故障数据中，有16个是已知故障类型，如冷却水入口温度或进料成分的变化、阀门粘滞、反应动力学漂移等，还有5个故障类型是未知的。为了对该过程进行监测，选取如表1所示的33个过程变量，接下来结合该TE过程对本发明具体实施步骤进行详细的阐述。

表1：TE过程监测变量。

序号	变量描述	序号	变量描述	序号	变量描述
						1	物料A流量	12	分离器液位	23	D进料阀门位置
2	物料D流量	13	分离器压力	24	E进料阀门位置
						3	物料E流量	14	分离器塔底流量	25	A进料阀门位置
4	总进料流量	15	汽提塔等级	26	A和C进料阀门位置
						5	循环流量	16	汽提塔压力	27	压缩机循环阀门位置
6	反应器进料	17	汽提塔底部流量	28	排空阀门位置
						7	反应器压力	18	汽提塔温度	29	分离器液相阀门位置
8	反应器等级	19	汽提塔上部蒸汽	30	汽提塔液相阀门位置
						9	反应器温度	20	压缩机功率	31	汽提塔蒸汽阀门位置
10	排空速率	21	反应器冷却水出口温度	32	反应器冷凝水流量
						11	分离器温度	22	分离器冷却水出口温度	33	冷凝器冷却水流量

首先，利用TE过程正常工况下的采样数据建立故障检测模型，包括以下步骤：

(1)收集生产过程正常运行状态下的数据样本，组成训练数据集X∈R^960×33，并对每个变量进行标准化处理，得到均值为0，标准差为1的新数据矩阵

(2)将标准化后的数据集

表示成

其中x_i∈R^960×1为第i个测量变量的n数据组成的列向量，i＝1，2，…，33，并初始化i＝1。

(3)将矩阵

中第i列数据去除得到回归模型的输入矩阵

而将x_i作为回归模型的输出，利用ICR算法建立输入矩阵

与输出x_i之间的软测量模型：

(4)判断是否满足条件i＜33；若是，则置i＝i+1后返回步骤(3)；若否，则将得到的估计误差向量组成矩阵E＝[e₁，e₂，…，e₃₃]后继续执行下一步骤。

(5)将矩阵E作为新的训练数据矩阵，对其中的每一列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

(6)以

作为训练数据建立相应的ICA故障检测模型，并保留模型参数集Θ＝{G，H，Λ，D_lim，Q_lim}以备调用。

其次，为验证本发明方法在故障检测上的优越性，特以TE过程中第19类故障类型工况为例，该故障为反应器冷却水进口温度阶跃故障，针对该故障测试数据实施过程监测。

(7)收集新采样时刻的数据样本x∈R^1×33，对其实施与步骤(1)中相同的标准化处理得到新数据向量

后，初始化i＝1。

(8)将行向量

中的第i个元素y_i取出后得到输入向量

并按照如下所示公式计算y_i的估计误差

(9)判断是否满足条件i＜33？若是，则置i＝i+1后返回步骤(8)；若否，则将得到的误差组成向量f＝[f₁，f₂，…，f₃₃]并继续执行下一步骤。

(10)：对f实施与步骤(5)中相同的标准化处理得到新向量

(11)据如下所示公式计算监测统计量D与Q：

(12)判断D与Q的具体数值是否大于对应控制上限D_lim与Q_lim？若否，则当前样本为正常工况采样；若是，则当前采样数据来自故障工况。

如图2所示，本发明方法以及传统基于PCA的故障检测方法在故障19上的监测详情对比图，可以很明显地发现本发明方法在该故障发生后立即检测出该故障且一直持续不断的触发故障警报上限，而传统PCA与ICA方法存在很明显的故障漏报现象。因此，本发明方法的监测效果要优越于传统PCA与ICA方法。

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

Claims (3)

1.一种基于分散式ICR模型的非高斯过程监测方法，其特征在于，包括以下步骤：

离线建模阶段的实施过程如下所示：

步骤(1)：采集生产过程正常运行状态下的样本，组成训练数据集X∈R^n×m，并对每个测量变量进行标准化处理，得到均值为0，标准差为1的新数据矩阵

其中，n为训练样本数，m为过程测量变量数，R为实数集，R^n×m表示n×m维的实数矩阵；

步骤(2)：将将标准化后的数据集

表示成

其中x_i∈R^n×1为第i个测量变量的n个数据组成的列向量，i＝1，2，…，m，并初始化i＝1；

步骤(3)：将矩阵

中第i列数据去除得到软测量模型的输入矩阵

而将x_i作为软测量模型的输出，利用独立成分回归(Independent Component Regression，ICR)算法建立输入矩阵

与输出x_i之间的软测量模型：

上式中，e_i∈R^n×1为软测量估计误差向量；

步骤(4)：判断是否满足条件i＜m；若是，则置i＝i+1后返回步骤(3)；若否，则将得到的估计误差向量组成矩阵E＝[e₁，e₂，…，e_m]后继续执行下一步骤(5)；

步骤(5)：将矩阵E作为新的训练数据矩阵，对其中的每一列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

步骤(6)：利用独立成分分析(Independent Component Analysis，ICA)算法为

建立相应的ICA过程监测模型，并保留模型参数集Θ＝{G，H，Λ，D_lim，Q_lim}以备调用，其中G为分离矩阵，H为混合矩阵，Λ为独立成分的协方差矩阵，D_lim与Q_lim分别表示监测统计量D与Q的控制上限；

步骤(7)：收集新采样时刻的数据样本x∈R^1×m，对其实施与步骤(1)中相同的标准化处理得到新数据向量

后，初始化i＝1；

步骤(8)：将行向量

中的第i个元素y_i取出后得到输入向量

并按照如下所示公式计算y_i的估计误差f_i：

步骤(9)：判断是否满足条件i＜m；若是，则置i＝i+1后返回步骤(8)；若否，则将得到的误差组成向量f＝[f₁，f₂，…，f_m]并继续执行下一步骤；

步骤(10)：对f实施与步骤(5)中相同的标准化处理得到新向量

步骤(11)：根据如下所示公式计算监测统计量D与Q：

步骤(12)：判断D与Q的具体数值是否大于对应控制上限D_lim与Q_lim；若否，则当前样本为正常工况采样；若是，则当前采样数据来自故障工况。

2.根据权利要求1所述的一种基于分散式ICR模型的非高斯过程监测方法，其特征在于，所述步骤(3)中建立软测量模型的具体实施过程如下所示：

①设置k＝1后，初始化向量u_k＝x_i；

②根据公式

q_k＝x_i ^Ts_k/(s_k ^Ts_k)分别计算得到系数向量w_k、得分向量s_k、和系数q_k；

③根据公式u_new＝x_iq_k/q_k ²计算向量u_new；

④判断是否满足条件||u_k-u_new||＜10^-6；若否，则置u_k＝u_new后返回步骤②；若是，则执行下一步骤⑤；

⑤根据公式

计算得到第k个投影向量p_k∈R^(m-1)×1，并保留投影向量p_k、系数向量w_k、和系数q_k；

⑥根据公式

更新输入矩阵

后，判断矩阵

中的最大元素是否大于0.001；若是，则执行步骤⑦；若否，则在得到最终的投影矩阵P_i＝[p₁，p₂，…，p_k]与系数矩阵W_i＝[w₁，w₂，…，w_k]后，执行步骤⑧；

⑦判断k＜m-1；若是，则置k＝k+1后，返回步骤②；若否，则得到最终的投影矩阵P_i＝[p₁，p₂，…，p_k]与系数矩阵W_i＝[w₁，w₂，…，w_k]；

⑧根据公式

对

进行白化处理，得到白化后的矩阵Z_i∈R^n×k并初始化j＝1，其中，矩阵R_i＝W_i(P_i ^TW_i)^-1∈R^(m-1)×k、矩阵

⑨设置向量c_j为k×k维单位矩阵中的第j列；

⑩按照公式c_j←E{Z_ig(c_j ^TZ_i)}-E{h(c_j ^TZ_i)}c_j更新向量c_j，其中函数g(u)＝tanh(u)、函数h(u)＝[sech(u)]²、u＝c_j ^TZ_i表示函数自变量；

对更新后的向量c_j依次按照下式进行正交标准化处理：

c_j←c_j/||c_j|| (5)

重复步骤⑤～⑥直至向量c_j收敛，并保存向量c_j；

判断j＜k；若是，置j＝j+1后，重复步骤④～⑧；若否，则执行步骤⑨；

将得到的所有k个向量c₁，c₂，…，c_k组成矩阵C＝[c₁，c₂，…，c_k]∈R^k×k，并按照如下所示公式计算分离矩阵W₀∈R^(m-1)×k与混合矩阵A₀∈R^(m-1)×k：

A₀＝P_iD_i ^1/2C

W₀＝P_iD_i ^-1/2C (6)

根据公式

计算独立成分矩阵S₀∈R^n×k，并根据公式

计算独立成分与输出x_i之间的相关系数向量

找出行向量

中数值小于0.001的所有元素所对应的列标号，根据该列标号对应的剔除分离矩阵W₀与独立成分矩阵S₀中相应的列，从而得到新的分离矩阵W_i与独立成分矩阵S_i；

计算最小二乘回归系数向量b_i＝W_i(S_i ^TS_i)^-1S_i ^Tx_i，那么对应于输出x_i的ICR模型为：

上式中，e_i∈R^n×1为软测量估计误差向量。

3.根据权利要求1所述的一种基于分散式ICR模型的非高斯过程监测方法，其特征在于，所述步骤(6)中建立ICA过程监测模型的具体实施过程如下所示：

①计算

的协方差矩阵

②求解Ψ所有M个大于0.001的特征值η₁≥η₂≥…≥η_M所对应的特征向量β₁，β₂…，β_M；

③根据公式

对

进行白化处理，得到白化后的矩阵F∈R^n×M并初始化i＝1，其中，矩阵B＝[β₁，β₂…，β_M]∈R^m×M、对角矩阵Φ＝diag{η₁，η₂，…，η_M}；

④设置向量d_i为M×M维单位矩阵中的第i列；

⑤按照公式d_i←E{Fg(d_i ^TF)}-E{h(d_i ^TF)}d_i更新向量d_i，其中函数g(u)＝tanh(u)、函数h(u)＝[sech(u)]²、u＝d_i ^TF表示函数自变量；

⑥对更新后的向量d_i依次按照下式进行正交标准化处理：

d_i←d_i/||d_i|| (9)

⑦重复步骤⑤～⑥直至向量d_i收敛，并保存向量d_i；

⑧判断i＜M；若是，置i＝i+1后，重复步骤④～⑧；若否，则执行步骤⑨；

⑨将得到的所有M个向量c₁，c₂，…，c_M组成矩阵Ω＝[d₁，d₂，…，d_K]∈R^M×M，并按照如下所示公式计算分离矩阵G₀∈R^m×M与混合矩阵H₀∈R^m×M：

G₀＝BΦ^-1/2Ω

H₀＝BΦ^1/2Ω (10)

⑩计算H₀中每一列向量的长度，分别记为L₁，L₂，…，L_M，并将L₁，L₂，…，L_M按照数值大小进行降序排列得到l₁，l₂，…，l_M，那么保留的独立成分个数μ为满足下列条件的最小值：

将H₀中列向量长度最大的μ个列向量组成混合矩阵H∈R^m×μ，同时从G₀中取出与H对应的列向量组成分离矩阵G∈R^m×μ；

根据公式

计算出独立成分矩阵J后，计算J的协方差矩阵Λ＝(J^TJ)/(n-1)；

根据如下所示公式计算D_lim与Q_lim：

上两式中，置信水平α＝99％，

表示权重为g＝v/2a，自由度为h＝2a²/v的χ²分布，a与v分别是

对应的监测统计量Q的估计均值和估计方差。

Images