基于自相关观测和隐马尔科夫模型的统计过程监控

张 媛，陈 震，潘尔顺，奚立峰

(上海交通大学 机械与动力工程学院，上海 200240)

0 引言

统计过程控制(Statistical Process Control, SPC)是广泛应用于现代质量管理的控制方法，休哈特控制图是其中常用的工具。传统的SPC理论基于一个基本假设：过程的观测数据相互独立。然而在现代化的制造环境下，过程数据往往具有一定的自相关性，如果仍运用常规控制图进行过程监控，则会造成大量的虚发警报，对质量管理工作产生误导[1]。

自20世纪70年代，许多学者展开了对自相关过程进行有效统计控制的研究。文献[2]使用过程原始观测值构建控制图，通过调整控制限使数据间的自相关性消失。Alwan等[3]率先提出的残差控制图，通过对自相关数据拟合自回归移动平均模型(Auto-Regressive and Moving Average models, ARMA)，得到观测值与拟合值的残差序列，由时间序列理论可知，残差序列是独立同分布的，因此能够运用常规控制图进行分析监控，这也是该领域的主流研究方向。Lu等[4]对比分析了Shewhart、CUSUM控制图对一阶自相关过程数据的残差序列具有不同的检测效果;Dawod等[5]研究了基于鲁棒性和高效率的ARMA选型问题。此外，也有学者提出，结合使用SPC和工业过程控制(Engineering Process Control, EPC)，能更有效地提高过程质量并降低成本[6]。近些年，在自相关控制图的研究中加入贝叶斯理论和神经网络等技术成为一种趋势[7-8]。何桢等[9]以神经网络特有的模式识别技术监控自相关过程中的均值突变;Pan等[10]研究了用向量自回归模型拟合多元自相关过程;崔敬巍等[11]提出以贝叶斯动态模型为基础的自相关控制图，具有建模迅速的优势。

与神经网络类似，隐马尔科夫(Hidden Markov Models, HMM)模型也是一种有效的预测算法，是语音识别、图像处理等领域内的重要工具。关于HMM改进研究的最新方向之一是基于状态过程或观测过程的变化。Shen等[12]讨论了隐状态为多个马氏链过程的HMM模型，并验证了这种模型结构能有效应用于嘈杂环境中说话者的语音识别问题;Wang等[13]利用连续观测的HMM模型对刀具状态进行监控,并预测其剩余寿命。

将HMM方法与过程控制理论相结合的研究也受到一些学者的关注[14]。Yu[15]提出以HMM模型为基础综合考虑局部信息(马氏距离)和全局信息(负对数似然概率)，能有效监控非线性多模态过程;Li等[16]建立了一种基于增距映射的HMM模型，提高了参数估计的效率，可应用于大规模动态过程的监控和故障诊断问题;Alshraideh等[17]运用标准HMM(SHMM, Standard HMM)对一阶自相关过程进行拟合[17]，但标准HMM认为各观测值之间条件独立，观测序列的概率分布仅与当前的隐含状态有关，而不考虑观测值自身在时域上的相关性。针对这一问题，本文提出一种自相关HMM(Autocorrelaion HMM， AHMM)模型，并结合残差控制图发展了基于此模型的过程监控方法，相比ARMA残差控制图，该方法具有较强的可操作性和更高的预测精度。

1 考虑自相关观测的HMM

1.1 自相关观测的模型化

HMM是在Markov链的基础上发展起来的一种不完全统计模型，它包含两个随机过程：①不能被直接观测到的Markov链，描述隐藏状态的转移；②描述观测值与隐藏状态关系的随机过程。通常，一个标准的HMM由5个基本元素组成，包括隐藏状态、状态转移概率、初始状态概率、观测序列和观测概率分布，各元素的具体描述为：

(1)有限的隐藏状态集

S={S1,S2,…,SN},

(1)

式中N为隐藏状态的个数，可以通过模型选择方法确定，如交叉验证。

(2)状态转移概率

A={aij},其中

aij=P(qt+1=Sj|qt=Si),

(2)

式中qt为时刻t时的隐藏状态。

(3)初始状态概率

π={πi},其中

(3)

(4)观测序列

O={y1,…,yT},

(4)

式中T为观测序列的长度。

(5)观测概率分布

B={bj(·)},其中

bj(yt)=P(yt|qt=Sj),1≤j≤N。

(5)

N和T隐藏在其他参数中，因此取λ=(A,B,π)作为参数集，定义标准的HMM。

标准HMM，观测值被认为是相互独立的，即任意时刻t出现的观测量仅依赖于当前所处的状态，而与t时刻以前的观测值和状态无关。本文所考虑的观测值自相关情况下的AHMM，认为任意时刻t出现某观测量的概率不但依赖于系统当前所处的状态，而且依赖于t-1时刻出现的观测值。

考虑一阶自相关的观测过程。记t时刻的状态为Si，观测值yt的条件概率分布的均值为t-1时刻观测值yt-1的一阶线性函数，且yt与t-1时刻的状态无关。假设观测值的条件概率分布满足高斯分布，则对于给定的模型λ，当前观测值yt在状态Si下的条件概率分布均值为

μ(yt|yt-1,qt=Si)=ciyt-1+ci,0,

(6)

式中ci和ci,0为相关系数。构造两个二维向量，列向量x(t)=(yt-1,1)′和横向量Ci=(ci,ci,0)，则式(6)中的均值可表示为

μ(yt|yt-1,qt=Si)=Cix(t)。

(7)

bi,yt-1(yt)=P(yt|yt-1,qt=Si,λ)

(8)

1.2 参数估计

Baum-Welch算法是常见且有效的HMM参数估计方法。因为考虑了观测值之间的相关性，所以本文重新推导了基于Baum-Welch算法的参数估计表达式。对于未知参数集λ={A,C,π,σ2},给定观测序列O={y1,…,yT}，求使P(O|λ)最大的λ*即为参数估计值。

首先定义正向变量(forward variable),αt(i)为到时刻t观测到部分序列{y1,…,yt}且此时的状态是Si的概率，即

αt(i)=P(y1…yt,qt=Si|λ)。

(9)

αt(i)可通过以下递归方法计算：

(1)初始化

α1(i)=πibi,y0(y1),1≤i≤N。

(10)

式中y0为初始时刻t=0时的观测量。

(2)递归

1≤t≤T-1,1≤j≤N。

(11)

再定义一个反向变量(backward variable),βt(i)为到时刻t观测到部分序列{yt+1,…,yT}且此时的状态是Si的概率，即

βt(i)=P(yt+1…yT|qt=Si,λ)。

(12)

类似地，βt(i)也可通过递归方法计算：

(1)初始化

βT(i)=1,1≤i≤N。

(13)

(2)递归

1≤t≤T-1,1≤i≤N。

(14)

在给定观测O和参数λ时，系统在t时刻处于状态Si的概率γt(i)和系统在t时刻处于状态Si且t+1时刻处于状态Sj的概率ξt(i,j)，可分别由以下各式得到：

(15)

(16)

针对AHMM的参数估计方法，即基于Baum-Welch算法的更新算法可阐述如下：

步骤1求期望(E)。给定当前参数λ(k)，计算γt(i)和ξt(i,j)的值。

步骤2最大化(M)。由E步的计算结果可推导得到模型参数集的重新估计值为

(17)

(18)

(19)

(20)

(21)

将新的估计值记为λ(k+1)，返回E步，循环迭代直至结果收敛，得到最终的参数估计值λ*。

2 过程监控与残差控制图

下面结合改进的AHMM模型及Roberts和Wlan提出的残差控制图理论，构造基于自相关观测的过程监控方法。在合理设定隐藏状态个数N的前提下，自相关HMM模型能够较准确地拟合自相关过程的观测数据，利用所得的残差序列可以建立常规控制图对过程进行监控。具体实施步骤分为两个阶段：首先通过对历史观测数据进行模型训练，得到AHMM模型，进而实现对观测数据的预测；其次利用过程观测的预测值和实际值建立残差控制图。

2.1 预测

已知观测O和参数λ，时刻T+1的观测量为yT+1的条件概率为

qT=Si,O,λ)P(qT=Si|O,λ)

P(qT+1=Sj|qT=Si)γi(T)

(22)

(23)

式中yT+1的积分上下限由具体过程观测值的取值范围决定。

2.2 残差控制图的建立

第一阶段中确定的AHMM模型和数据预测方法是残差控制图建立的基础。

对残差序列e构建单值控制图(X控制图)，其控制限为：

(24)

延长上述X控制图的控制线，将其转化为控制用控制图。运用AHMM模型对T+1时刻的观测值进行预测，将其与实际值之差，即eT+1在控制图上描点，对后续各时刻重复此项操作。若发现异常点，需要立刻找出失控原因并采取措施加以消除，从而实现对过程质量的监控。

3 实例分析及仿真

3.1 自相关过程的控制实例

刘艳荣在《自相关过程中控制图的应用研究》里分析了一组转化率数据[18]，该组数据具有自相关性，但因为均值波动导致的非平稳性，不能使用简单的ARMA模型对其进行拟合，而需要利用差分对数据进行平稳化处理。经过分析验证，最终选定了ARIMA(0,1,1)模型，拟合结果如图1a所示。

利用本文提出的AHMM模型对该组转化率数据进行拟合得到的结果如图1b所示。对比图1b和图1a可以看出，除去3个峰值点，AHMM模型能更贴切地反映出数值的变化趋势。

另有某钢厂热轧过程的一组温度数据[19]，由自相关和偏自相关函数图判断该过程符合AR(1)模型(即ARIMA(1,0,0)模型)，拟合结果如图2a所示。利用AHMM模型对该组温度数据进行拟合，结果如图2b所示。

以上两个案例都反映，尽管AHMM模型在预测初期存在一定偏差，但是随着数据资料的累积，预测误差迅速减小，从而体现了其强大的学习能力。整体来看，基于AHMM预测的精确度明显高于ARIMA模型。表1列出的预测误差方差、平均绝对误差和平均相对误差绝对值3项参数的对比，也有力地说明了这一优势。此外，对于案例一中的非平稳时间序列，AHMM模型无需对数据采取平稳化处理，在实际应用中减少了数据处理的工作量。

表1 两种模型预测精度相关参数表

基于AHMM的预测结果，进一步建立两个案例的残差控制图(如图3和图4)。图3显示转化率数据对应的第38,69,83个点落在了控制限之外，为异常点,这一结果与文献[18]中的分析一致，这些数据对应的时刻确实有“异常原因”存在。同时也说明图1b中的预测值在3个峰值点的偏离是合理的，因其不符合数据正常的波动趋势。

图4中所有的残差数据均落入控制图上下控制限范围内，说明热轧温度比较平稳，以此为指标监控热轧过程可判断该过程正常，这一结论也与实际相符。

3.2 基于平均运行链长的仿真分析

平均运行链长(Average Run Length, ARL)是评价控制图性能的一个重要指标。本文采用蒙特卡洛模拟法，通过计算监视过程均值偏移的ARL，比较ARMA,SHMM,AHMM 3种模型在自相关过程中捕捉均值偏移异常的灵敏度。仿真数据由最简单的自相关过程AR(1)模型产生。

AR(1)模型可以表示为yt=μ+φyt-1+εt，设μ=100，εt～N(0,1)，自相关系数φ分别取0.5和0.7，通过MATLAB仿真产生400个数据，依次训练3种模型。对于任一模型，通过拟合得到400个训练数据的预测值，进而与真值相减得到残差，利用残差序列计算出控制图的控制限，用于后序ARL的计算。过程失控时，均值偏移量分别取1σ,2σ，针对3种模型分别进行了800组仿真，得到其对应的控制图ARL如表2所示。

表2 基于AR(1)、标准HMM、AHMM模型的残差控制图ARL

由表中数据可以看出，与AR(1)模型相比，基于HMM模型的残差控制图对异常的捕捉更灵敏，而AHMM模型的表现最为优异，这与其对数据拟合的高精度有着必然的关系。值得注意的是，仿真数据根据AR(1)模型的表达式生成，可以认为数据的变化趋势与AR(1)模型完全契合，对于一般性的自相关过程，AR(1)模型的表现会更差。

对比两种HMM模型的残差控制图ARL，考虑到AHMM能更准确地描绘数据间的相关关系，当数据的相关系数增大时，它对均值偏移异常的捕捉速度有显著的提高，而标准HMM模型只是通过状态过程的变化来预测数据的波动趋势，因此其表现并没有随着数据相关性的增加而明显提升。

正如文献[17]指出的，利用ARMA模型对数据进行拟合时，需要经过多次尝试来确定自回归项和移动平均项的系数，对于非平稳数据还要确定差分的阶数，这一过程中往往需要以经验为依据进行主观判断，因此选型过程存在一定的误差概率。而HMM模型的确定更有据可循，只要选定合适的状态数即可，简化了选型过程，在生产实际中减少了一定的时间和人力成本。

4 结束语

在连续生产过程中，采集到的数据往往具有相关性，经典控制理论无法发挥作用。残差控制图是解决自相关过程监控问题的有力工具，通常要以ARMA模型的数据预测为基础，但其存在选型过程复杂、对异常识别不够迅速等问题。本文提出的基于自相关HMM模型建立残差控制图，参数确定简便，并且能准确地把握数据的变动趋势，提高了生产实际中对自相关过程进行质量监控的效率，具有良好的应用前景。

本文对HMM进行优化建模时，利用一阶线性函数来描述观测的自相关性。未来可以进一步研究高阶非线性函数在HMM 模型中对观测序列自相关性的表达，从而更准确地捕捉观测的变化趋势。