基于独立元分析的制浆造纸废水处理过程故障检测

杨冲　宋留　刘鸿斌

摘 要：为及时、准确地做出故障诊断，本课题采用独立元分析（ICA）和主成分分析（PCA）两种常用的多元统计分析方法对制浆造纸废水处理过程中的传感器故障进行检测并对诊断效果进行对比。结果表明，对于制浆造纸废水数据中偏移和漂移两种故障，ICA模型的故障检测率分别为24%与54%，PCA模型的故障检测率分别为14%和42%，ICA模型的两种故障检测率均高于PCA模型，但是两种模型均无法达到满意的检测效果;对于完全失效故障，ICA和PCA模型的故障检测率均达到100%。

关键词：制浆造纸废水处理过程;故障检测;主成分分析;独立元分析

中图分类号：X793

文献标识码：A

随着现代工业过程的规模与复杂程度的提高，工业过程的安全运行和产品质量得到更多的关注，故障检测与分析技术也越来越不可或缺[1]。

一般来说，故障检测与分析可以分为3类：基于解析模型、基于专家知识和基于数据驱动的故障诊断方法。基于解析模型的方法可以使系统根据特定的模型得出精确的结果，但是随着系统复杂程度的提高，解析模型构造的难度与成本也随之上升[1]。基于专家知识的方法依赖于已存在的知识库与专家的经验知识以解决一些非结构化的问题，适用于过程尚未成熟的领域，但是完全依靠专家凭感觉、个人经验及简单仪表来判断和排除故障成为该方法的局限性[2]。相对于前两种故障诊断方法，基于数据驱动的方法没有以上限制，该方法在不需要知道系统精确解析模型的情况下可对系统运行数据进行采集、建模、分析，从而完成系统的故障诊断。

制浆造纸废水的处理过程具有较强的时变性、耦合性、非线性以及易受外界干扰的特性。作为制浆造纸废水处理过程中使用最广泛的方法——活性污泥法主要包括了预处理、初次沉降、曝气、二次沉降和污泥回流5个部分（如图1所示）。活性污泥法处理制浆造纸废水的原理是利用在曝气池内的微生物群体的凝聚、吸附、氧化分解等作用去除制浆造纸废水中的有机物。在每一部分处理工艺中都需要测量大量的在线参数，用来保证出水水质达标以及完成控制过程的优化，因此，该方法对测量仪表的精密性有较高的要求。测量仪表精密度的提高增加了故障出现的概率，结合制浆造纸废水处理过程的恶劣工作环境以及设备长期满负荷运行导致老化等因素，故障检测结果的时效性与准确性需进一步提高[3]。

制浆造纸废水处理复杂的过程特性使我们难以对其构建准确的数学模型与完备的专家知识库;另一方面，制浆造纸废水处理厂的长期满负荷运行使其积累了大量的过程历史数据，因此采用数据驱动的方法对制浆造纸废水处理过程的故障进行诊断逐渐受到学术界和工业界的重视。制浆造纸废水处理领域中，基于数据驱动的故障诊断方法是以机器学习与多元统计方法为主。机器学习方法的核心在于机器的自我学习、改善性能，主要包括人工神经网络、贝叶斯和支持向量机等[3]。多元统计方法的提出是区别于单变量统计方法，主要利用多变量之间的相关性及内在的统计规律进行故障检测，适用于具有强相关性的制浆造纸废水处理系统[4]。多元统计方法是将含多变量的高维数据通过多元投影的方式进行降维，在相对独立的低维空间构造统计量进行数据分析。

近年来，随着多元统计过程监测理论的发展，独立元分析法（Independent Component Analysis， ICA）获得了较大的关注。针对系统运行过程中检测变量间的相关性，ICA方法可生成一个模型，该模型假设数据变量由少数潜在变量线性组合而成，相关系数未知。这些潜在变量称为数据的独立元，相互独立且非高斯。ICA模型旨在提取潜在变量以及混合矩阵的相关系数并结合相应的统计量进行故障检测[5]。主成分分析法（Principal Component Analysis， PCA）在流程工业领域也有广泛的应用。尽管ICA模型可看作PCA模型的有效延伸，但两者迥异。PCA模型旨在降维，将多个变量映射为少数几个主成分，这些主成分可由原始变量线性组合而成，并且包含着原始变量的大部分信息，为避免主成分包含的信息相互重叠，该方法要求各主成分互不相关[6]。PCA模型在计算统计量控制限时，需要假设数据服从高斯分布，如果测量数据非高斯，那么该方法的故障诊断效果将受到影响[7]。

结合制浆造纸废水处理过程的特点，本课题采用多元统计方法中的PCA与ICA两种方法对该过程数据进行故障检测与对比分析。

1 方法原理

1.1 PCA模型

基于PCA模型可以快速、有效地处理大量过程数据的特点，从20世纪80年代起，PCA模型开始被广泛应用于数据降维、图像处理、特征提取、模式识别以及过程控制等领域[8]。

1.1.1 数据模型

假设x∈Rm表示含m个样本的向量，每个样本均有N个独立采样。由此构造数据矩阵X=[x1，x2，…，xN]T∈RN·m并进行标准化处理。PCA模型对标准化后的矩阵X进行分解[9]：

X=TPT+=TPT+T=[T][P]T≡T（1）

式中，P∈Rm·l为负载矩阵;T∈RN·l代表得分矩阵，l代表主元个数;=TPT为残差矩阵，且=[T]，=[P]。标准化后样本X的协方差矩阵估计值S可表示为：

S=1N-1XTX=T（2）

則

=1N-1T=diag{λ1，λ2，…，λm}（3）

式中，=diag{λ1，λ2，…，λm}为S的m个特征值并按照降序排列，特征值对应的特征向量为的m个列向量，这些特征向量代表的信息量也按照降序排列，为标准正交。这里采用累积贡献率来选取主元个数l，总方差中第i个主成分的方差所占的比例称为该主成分pi的贡献率，主成分pi的贡献率反映了第i个主成分包含的信息量大小。一般通过使前l（l≤m）个主成分的累计贡献率不低于某一数值（如80%）的方法来确定l的数值[6]。

PCA模型建立之后，主元子空间（Principal Component Subspace， PCS）和残差子空间（Residual Subspace， RS）可分别表示为Sp=span{P}，Sr=span{}。则任意样本向量在PCS与RS上的投影可分别表示为[9]：

=PPTx≡Cx∈Sp（4）

=Tx=（I-PPT）x=（I-C）x∈Sr（5）

因为Sp和Sr相互正交，则T=0，且x=+。

1.1.2 故障检测

一般来说，基于PCA模型的故障检测方法采用T2统计量和SPE统计量进行故障检测，分别定义为[9]：

T2=tTΛ-1t=xTPΛ-1PTx（6）

SPE≡‖‖2=‖（I-PPT）x‖2（7）

T2的控制限T2α利用F分布确定：

N（N-l）l（N2-1）T2～Fl，N-l（8）

T2α≡l（N2-1）N（N-l）Fl，N-l;α（9）

式中，F分布自由度为l和N-l，l为主成分个数，N为采样次数，α为显著性水平。SPE的控制限δ2α利用χ2分布确定：

δ2α=gχ2h;α（10）

式中，g=θ2/θ1，h=θ21/θ2，θi表示为：

θi=∑mj=l+1λij， i=1，2，3（11）

如果任意统计指标的值超出各自的控制限，则认为过程出现了故障。

1.2  ICA模型

ICA模型用来提取变量、测量值或者信号间潜在的独立元（Independent Components，ICs），这种方法最初被用来解决未知信息资源的分离问题，包括恢复与未知矩阵进行线性组合后的独立信号（如不同的声音、音乐或者噪音信号等）[10]。

1.2.1 数据模型

ICA模型旨在基于观测数据X通过线性组合寻求独立、非高斯的潜在变量S以及未知矩阵A，数据测量模型如下[5]：

X=AS+E（12）

X=[x1，x2，…，xn]∈Rd·n为观测数据矩阵，A=[a1，a2，…，am]∈Rd·m是未知的混合矩阵，S=[s1，s2，…，sn]∈Rm·n（m≤d）是独立元矩阵，E∈Rd·n是残差矩阵，n为样本的数量。一般认为X和S同维，即m=d，此时矩阵E为零矩阵，则ICA模型可转化为寻求一个混合矩阵W以满足重构矩阵各元素间最大的独立性，可表示为：

=WX（13）

为方便数学表述，定义独立元具有单位方差，即E{ssT}=Im。通常求解ICA模型第一步是对观测数据X进行白化处理，使观测向量x经线性变换表示为：

z=Qx（14）

经过观测数据的白化处理，使z的各变量间线性无关且具有单位方差，即E{zzT}=Im。则可得到

z=Qx=QAs=Bs（15）

这里定义B=QA，根据式（16）得出B为正交矩阵。

E{zzT}=BE{ssT}BT=BBT=Im（16）

由于正交矩阵的限制，ICA模型问题得到简化，这里不必再寻求满秩矩阵A，转而寻求正交矩阵B。根据公式（15），重构向量可表示为：

=BTz=BTQx（17）

根据式（13）和式（17），矩阵B与W间的关系可表示为：

W=BTQ（18）

根据i=（bi）Tz，为确保第i个独立元素的最大非高斯性，B矩阵每列向量bi确定前均进行初始化处理。Hyvrinen和Oja在文献[11]中根据中心极限定理提出非高斯性代表着独立性，而峰度和负熵可用来度量非高斯性。基于对负熵的近似估计，Hyvrinen在文献[12]中提出对于矩阵B的算法。计算出矩阵B后，结合式（17）和式（18）可分别得出重构的独立元矩阵和混合矩阵W。

获取独立元后，需选取其中一部分作为主元以提高ICA模型的鲁棒性并降低分析过程的复杂性。选取主元之前，需对ICs进行排序。文献[13]中根据混合矩阵W中各行的Euclidean范数（L2）大小，将其对应的ICs进行排序，该方法简单有效。确定ICs的排序问题之后，Lee等[5]采取绘图的方式来确定主元个数。混合矩阵W中各行Euclidean范数（L2）所占的比例如图2所示。横坐标代表ICs的序号，纵坐标代表对应混合矩阵W中每行L2所占比例，该图中的混合矩阵W根据正常情况下的废水数据（2.1部分）所得。由图2可知，后3个独立元对应的L2比例明显小于其他的独立元对应的L2比例，因此，图2中ICs所建立的ICA模型包含5个主元。

1.2.2 故障检测

根据1.2.1确定主元个数（设为l）后，记Wd为W的主元部分，其由l个具有最大Euclidean范数（L2）的行向量构成，混合矩阵W中剩余部分为残差部分We，根据式（13）可得，d=Wdx，e=Wex，据式（18）得：Bd=（WdQ-1）T[5]。

这里，ICA模型的故障检测采用3种统计量作为监控指标：

I2=Tdd（19）

I2e=Tee（20）

SPE=eTe=（x-）T（x-）（21）

式中，=Q-1Bdd。以上3種统计量的控制限可通过核密度估计获取：

f（x）=1nh∑ni=1Kx-xih（22）

式中，xi为数据集中的观测值，h为平滑参数，K为核函数，n为样本数，h的取值可参考文献[14]。当x点达95%的密度函数值时，x为3种统计量的控制限。

2 仿真实验与讨论

2.1 制浆造纸废水处理过程数据

制浆造纸废水数据采取广东东莞某造纸厂的好氧段废水监测数据，该制浆造纸废水数据包括8个废水变量，每个变量包括170个测量值。制浆造纸废水处理过程数据如图3所示，左侧纵坐标分别对应进水化学需氧量（CODinf）、出水化学需氧量（CODeff）、进水悬浮固形物（SSinf）、出水悬浮固形物（SSeff）;右侧纵坐标分别对应溶解氧量（DO）、流量（Q）、温度（T）、pH值[15]。170个测量数据中，将前100个样本作为训练集，后70个样本作为测试集。

2.2 故障数据构建

针对制浆造纸废水处理所面临的复杂工作环境并结合传感器在数据监测中起到的重要作用，本课题考虑构建以下3种传感器故障：偏移故障、漂移故障和完全失效故障[16]。3种故障数据的构建如表1所示。在表1中，对变量CODinf加入均值的20%作为偏移故障，对变量pH值加入以0.05为系数的漂移故障，将SSeff的数据改至30，作为完全失效故障。如图4所示，在测试集后50组数据中加入故障信息。

2.3 故障检测

PCA和ICA两种方法模型对3种故障类型的检测结果如表2和表3及图5～图7所示。所有图中的虚线代表95%的置信度控制限，超出控制限的统计量代表出现了故障，根据测试集后50组数据中超出样本控制限的检测点个数所占的比例计算故障检测率，对于每种故障对应的多个统计量，选择故障检测率最高的统计量对应结果作为该故障的检测率。

由表2和表3及圖5～图7可知，对于偏移故障，PCA模型的SPE统计量与ICA模型的I2统计量从测试集样本点21处开始均有明显的上升，但是两者的检测效果不佳;PCA模型的偏移故障检测率为14%，ICA模型的偏移故障检测率为24%。对于漂移故障，PCA模型的T2统计量与ICA模型的I2e、SPE统计量均表现出明显的线性增长趋势;相对而言，ICA模型的漂移故障检测率达到54%，检测效果优于PCA模型。对于完全失效故障，PCA模型与ICA模型的各个统计量在故障出现后的第1个样本开始均有明显变化，且PCA的SPE统计量和ICA的I2和SPE统计量对应的完全失效故障检测率均达到100%。

为了优化故障检测的效率，结合变量贡献图与统计量的值进行故障检测可高效地检测出故障所属的变量类型。以ICA模型对于图7（c）所检测出的完全失效故障点为例，对位于测试集第70个样本点进行故障的变量贡献值检测，图8（a）～图8（c）分别显示了统计量I2、Ie2和SPE对于故障点的变量贡献值，结果显示，变量SSinf和变量SSeff对完全失效故障的贡献值最高，变量CODinf、变量CODeff以及变量DO的贡献值次之，其余变量贡献值均很低。由此，可判断变量SSinf和变量SSeff与完全失效故障有较强的相关性，根据以上信息可有效缩小故障排查范围，提高故障诊断效率。

3 结 论

本课题分别采用主成分分析（PCA）模型与独立元分析（ICA）模型对制浆造纸废水处理过程进行故障诊断。考虑到PCA模型控制限的设置对于测量数据高斯性的要求，结合制浆造纸废水处理过程中数据呈非高斯分布的特点，PCA模型的限制性相比于ICA模型较为明显。相比于PCA模型采用的二阶统计信息，ICA模型的优势在于采用I2、I2e、SPE 3种统计量进行故障检测。其中统计量I2替代了PCA模型中T2统计量，ICA模型在确定主元个数时所排除掉的独立元矩阵可用于构建统计量I2e，以此来弥补主元个数选取不当引起的故障检测能力下降的缺点。此外，ICA模型在故障诊断中要求潜变量之间统计独立，没有对潜变量正交性与高斯性的限制。

制浆造纸废水处理过程的故障检测结果表明，ICA模型的整体故障检测率要高于PCA模型，但优势并不明显。制浆造纸废水处理过程有着较强时变性，随着时间的改变，变量的均值与方差会随之改变。而常规的多元统计过程监测方法（如PCA和ICA方法）属于静态的统计监测方法，因此在废水数据故障检测中可能造成较大的误报与漏报现象，造成故障检测效果不佳。鉴于此，未来考虑将动态ICA方法应用到制浆造纸废水处理过程监测领域。

参 考 文 献

[1] Ge Z， Song Z， Gao F. Review of recent research on data-based process monitoring[J]. Industrial & Engineering Chemistry Research， 2013， 52（10）： 3543.

[2] Liu T L， Shen W H. A review of application of fault diagnostic expert system in wastewater treatment[J]. Paper Science & Technology， 2011， 30（2）： 75.

刘天龙， 沈文浩. 污水处理过程中故障诊断专家系统的应用[J]. 造纸科学与技术， 2011， 30（2）： 75.

[3] Huang D P， Qiu Y， Liu Y Q， et al. Review of data-driven fault diagnosis and prognosis for wastewater treatment[J]. Journal of South China University of Technology （Natural Science Edition）， 2015， 43（3）： 111.

黄道平， 邱 禹， 刘乙奇， 等. 面向污水处理的数据驱动故障诊断及预测方法综述[J]. 华南理工大学学报（自然科学版）， 2015， 43（3）： 111.

[4] Zhou D H， Hu Y Y. Fault diagnosis techniques for dynamic systems[J]. Acta Automatic Sinica， 2009， 35（6）： 748.

周东华， 胡艳艳. 动态系统的故障诊断技术[J]. 自动化学报， 2009， 35（6）： 748.

[5] Lee J M， Yoo C K， Lee I B. Statistical process monitoring with independent component analysis[J]. Journal of Process Control， 2004， 14（5）： 467.

[6] Xie Z H. Matlab statistical analysis and application of 40 case studies[M]. Beijing： Beihang University Press， 2010.

謝中华. MATLAB统计分析与应用： 40个案例分析[M]. 北京： 北京航空航天大学出版社， 2010.

[7] Qin S J. Survey on data-driven industrial process monitoring and diagnosis[J]. Annual Reviews in Control， 2012， 36（2）： 220.

[8] Shen Y， Ding S X， Haghani A， et al. A comparison study of basic data-driven fault diagnosis and process monitoring methods on the benchmark Tennessee Eastman process[J]. Journal of Process Control， 2012， 22（9）： 1567.

[9] Qin S J. Statistical process monitoring： basics and beyond[J]. Journal of Chemometrics， 2003， 17（8/9）： 480.

[10] Olsson G. ICA and me—A subjective review[J]. Water Research， 2012， 46（6）： 1585.

[11] Hyvrinen A， Oja E. Independent component analysis： algorithms and applications[J]. Neural Networks， 2000， 13（5）： 411.

[12] Hyvrinen A. Fast and robust fixed-point algorithms for independent component analysis[J]. IEEE Transactions on Neural Networks， 1999， 10（3）： 626.

[13] Cardoso F， Souloumiac A. Blind beamforming for non Gaussian signals[J]. Radar & Signal Processing IEE Proceedings F， 1993， 140（6）： 362.

[14] Silverman B W. Density estimation for statistics and data analysis[M]. London： Chapman & Hall， 1986.

[15] Liu Y Y， Li X Y， Zhang G R， et al. Sensor fault detection of papermaking wastewater treatment processses based on multivariate statistical analysis[J]. China Pulp & Paper Industry， 2017， 38（8）： 41.

刘耀瑶， 李祥宇， 张光锐， 等. 基于多元统计分析的造纸废水处理过程传感器故障检测[J]. 中华纸业， 2017， 38（8）： 41.

[16] Dunia R， Qin S J， Edgar T F， et al. Identification of faulty sensors using principal component analysis[J]. AICHE Journal， 1996， 42（10）： 2797.

Fault Detection of Papermaking Wastewater Treatment Process Based on Independent Component Analysis

YANG Chong1 SONG Liu1 LIU Hongbin1，2，*

（1Co-Innovation Center of Efficient Processing and Utilization of Forest Resources， Nanjing Forestry University，

Nanjing， Jiangsu Province， 210037; 2State Key Laboratory of Pulp and Paper Engineering，

South China University of Technology， Guangzhou， Guangdong Province， 510640）

（*E-mail：hongbinliu@njfu.edu.cn）[JZ）]

Abstract：To monitor and control papermaking wastewater treatment process（WWTP） effectively， two common methods of multivariate statistical analysis named independent component analysis （ICA） and principal component analysis （PCA） were used to detect the sensor faults in a papermaking WWTPThe results showed that the detection rates of the bias and drifting faults using ICA were 24% and 54%， respectivelyMeanwhile， the bias and drifting faults detection rates of PCA were 14% and 42%The fault detection rates of ICA were higher than those of PCA， but neither of the two methods achieved satisfactory result of detecting the bias and drifting faultsConcerning the complete failure fault， both the fault detection rates of the two methods were 100%

Keywords：papermaking wastewater treatment process; fault detection; principal component analysis; independent component analysis