基于共同趨勢(shì)模型的非平穩(wěn)過程在線監(jiān)控

林原靈，陳前

（南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210016）

基于共同趨勢(shì)模型的非平穩(wěn)過程在線監(jiān)控

林原靈，陳前

（南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210016）

對(duì)于非平穩(wěn)過程監(jiān)控，傳統(tǒng)的基于數(shù)據(jù)平穩(wěn)假設(shè)的多元統(tǒng)計(jì)過程控制方法是不適用的。針對(duì)上述問題，提出了一種基于共同趨勢(shì)模型的非平穩(wěn)過程監(jiān)控方法。共同趨勢(shì)模型從存在協(xié)整關(guān)系的非平穩(wěn)多元變量中辨識(shí)出共同因子，將各非平穩(wěn)過程變量分解成非平穩(wěn)的共同趨勢(shì)成分與平穩(wěn)成分之和的形式。不同于現(xiàn)有的基于協(xié)整模型的非平穩(wěn)過程監(jiān)控方法，共同趨勢(shì)模型能夠獲取各非平穩(wěn)變量中的平穩(wěn)成分，消除非平穩(wěn)共同因子的影響并體現(xiàn)變量間全部的動(dòng)態(tài)均衡關(guān)系。將對(duì)非平穩(wěn)過程的監(jiān)控變?yōu)閼?yīng)用共同趨勢(shì)模型，分解得到各非平穩(wěn)過程變量中的平穩(wěn)成分，然后應(yīng)用傳統(tǒng)的多元統(tǒng)計(jì)方法，估計(jì)平穩(wěn)成分的統(tǒng)計(jì)量及相應(yīng)的控制限進(jìn)行監(jiān)測(cè)。石油蒸餾過程監(jiān)控的實(shí)例研究結(jié)果表明，所提出的方法比基于協(xié)整新息變量的方法具有更可靠的監(jiān)控效果。

共同趨勢(shì)模型；過程監(jiān)控；非平穩(wěn)過程；協(xié)整檢驗(yàn)；過程控制；過程系統(tǒng)；系統(tǒng)工程

引 言

隨著技術(shù)的發(fā)展，現(xiàn)代工業(yè)過程系統(tǒng)的規(guī)模日益擴(kuò)大，復(fù)雜程度不斷增加，系統(tǒng)發(fā)生故障的風(fēng)險(xiǎn)隨之提高。系統(tǒng)發(fā)生故障，往往造成人員傷亡和巨大的經(jīng)濟(jì)損失。因此，必須對(duì)工業(yè)過程實(shí)時(shí)監(jiān)控以確保能夠及時(shí)發(fā)現(xiàn)并排除故障。隨著傳感器和計(jì)算機(jī)技術(shù)在過程監(jiān)控中的廣泛應(yīng)用，海量的過程數(shù)據(jù)被記錄和保存。過程監(jiān)控?cái)?shù)據(jù)通常呈現(xiàn)出高維度、結(jié)構(gòu)化分布、維度間強(qiáng)相關(guān)性的特征。為了有效分析過程監(jiān)控中產(chǎn)生的大量數(shù)據(jù)進(jìn)而揭示出故障征兆，數(shù)據(jù)驅(qū)動(dòng)方法被應(yīng)用于工業(yè)過程監(jiān)控。數(shù)據(jù)驅(qū)動(dòng)方法無需太多關(guān)于過程系統(tǒng)的先驗(yàn)知識(shí)，而是通過學(xué)習(xí)大數(shù)據(jù)樣本來獲取數(shù)據(jù)的結(jié)構(gòu)、分布規(guī)律等信息。多元統(tǒng)計(jì)過程控制（multivariate statistical process control，MSPC）作為一種數(shù)據(jù)驅(qū)動(dòng)方法，在工業(yè)生產(chǎn)中得到了最廣泛的應(yīng)用[1-2]。MSPC從表征過程系統(tǒng)正常運(yùn)行的數(shù)據(jù)中提取出數(shù)據(jù)的有效信息，應(yīng)用統(tǒng)計(jì)學(xué)理論建立統(tǒng)計(jì)監(jiān)控模型，實(shí)現(xiàn)工業(yè)過程的在線監(jiān)控。

經(jīng)典的 MSPC方法中，主元分析（principal component analysis，PCA）、偏最小二乘（partial least square，PLS）、獨(dú)立成分分析（independent component analysis，ICA）等[3-10]方法得到了廣泛應(yīng)用。進(jìn)一步，為了提取數(shù)據(jù)的非線性特征，核主元分析（kernel principal component analysis，KPCA）、核偏最小二乘（kernel partial least square，KPLS）、核獨(dú)立成分分析（kernel independent component analysis，KICA）等[11-15]方法被發(fā)展并被應(yīng)用于過程監(jiān)控。但是上述方法對(duì)數(shù)據(jù)的統(tǒng)計(jì)特性有假設(shè)限制，例如，PCA、PLS假設(shè)數(shù)據(jù)服從平穩(wěn)高斯分布。ICA則利用了數(shù)據(jù)的高階統(tǒng)計(jì)量信息，不要求數(shù)據(jù)服從高斯分布假設(shè)，但是常用的 ICA算法，如 Infomax、FastICA等，仍然假設(shè)數(shù)據(jù)分布是平穩(wěn)的[16]。實(shí)際上，受過程系統(tǒng)內(nèi)部復(fù)雜結(jié)構(gòu)的相互作用和外界環(huán)境因素的影響，工業(yè)大系統(tǒng)的隨機(jī)過程數(shù)據(jù)并非總是平穩(wěn)分布的。對(duì)于存在非平穩(wěn)變量的過程，傳統(tǒng)的基于平穩(wěn)隨機(jī)過程假設(shè)的MSPC方法是不適用的。

對(duì)于非平穩(wěn)過程監(jiān)控，經(jīng)典的方法是建立自回歸積分移動(dòng)平均（autoregressive integrated moving average，ARIMA）模型[17]。然而，建立ARIMA模型需對(duì)時(shí)間序列數(shù)據(jù)做差分，導(dǎo)致了數(shù)據(jù)中動(dòng)態(tài)信息的丟失，并使得用已建立的ARIMA模型對(duì)新樣本做預(yù)測(cè)成為難題。協(xié)整理論[18-20]是近年來計(jì)量經(jīng)濟(jì)學(xué)學(xué)者發(fā)展的一種非平穩(wěn)多元時(shí)間序列分析方法，已成功在催化裂化過程（FCCU）、石油蒸餾過程等[21-23]工程領(lǐng)域證明了其有效性。協(xié)整理論指出，對(duì)于存在動(dòng)態(tài)均衡關(guān)系的多變量非平穩(wěn)過程，非平穩(wěn)變量的某個(gè)線性組合是平穩(wěn)的。協(xié)整方法先用協(xié)整檢驗(yàn)方程求出協(xié)整向量，再將非平穩(wěn)變量在協(xié)整向量空間上做投影，得到平穩(wěn)的協(xié)整新息變量。協(xié)整理論將非平穩(wěn)時(shí)間序列的長期趨勢(shì)和短期動(dòng)態(tài)特征相結(jié)合，克服了ARIMA模型的不足。但是，現(xiàn)有的基于協(xié)整模型的監(jiān)控方法存在缺陷：新息變量只包含了一組協(xié)整關(guān)系，對(duì)于存在多組協(xié)整關(guān)系的變量，未能體現(xiàn)變量間的全部動(dòng)態(tài)均衡關(guān)系，因而抗干擾性不強(qiáng)。為了增強(qiáng)非平穩(wěn)過程監(jiān)控效果，本研究提出了基于共同趨勢(shì)模型的非平穩(wěn)過程監(jiān)控方法。共同趨勢(shì)模型是在協(xié)整理論基礎(chǔ)上發(fā)展的一種非平穩(wěn)時(shí)間序列模型。Stock等[24]提出了共同趨勢(shì)表達(dá)，證明了存在協(xié)整關(guān)系的多變量非平穩(wěn)過程中的各變量可以分解成非平穩(wěn)的共同趨勢(shì)成分與平穩(wěn)成分之和的形式，其中共同趨勢(shì)成分由稱為共同因子的一階單整過程驅(qū)動(dòng)生成。相較于協(xié)整模型，共同趨勢(shì)模型的最大優(yōu)點(diǎn)是能夠消除非平穩(wěn)共同因子的影響，提取出各非平穩(wěn)變量的平穩(wěn)成分，平穩(wěn)成分包含了變量間全部的動(dòng)態(tài)均衡關(guān)系。對(duì)于非平穩(wěn)過程監(jiān)控，先用共同趨勢(shì)模型提取出各非平穩(wěn)變量的平穩(wěn)成分，再計(jì)算出平穩(wěn)成分的統(tǒng)計(jì)量及其相應(yīng)的控制限，從而實(shí)現(xiàn)非平穩(wěn)工業(yè)過程的在線監(jiān)控。

1 共同趨勢(shì)模型

1.1 平穩(wěn)性檢驗(yàn)與協(xié)整檢驗(yàn)

多元隨機(jī)過程變量建立共同趨勢(shì)模型的前提條件是：①過程變量的各分量是同階單整非平穩(wěn)的；②變量間存在協(xié)整關(guān)系。因此，首先要對(duì)各變量進(jìn)行平穩(wěn)性檢驗(yàn)以確保各變量是同階單整的。對(duì)于一個(gè)非平穩(wěn)的時(shí)序變量xt，定義xt的一階差分形式為Δxt≡xt?xt?1，若 Δxt是平穩(wěn)的，則稱 xt是一階單整的，記為 I(1)。為簡單起見，下文討論的非平穩(wěn)變量均是一階單整的。目前較為通用的平穩(wěn)性檢驗(yàn)方法是 Augmented Dickey-Fuller（ADF）單位根檢驗(yàn)[25]，其核心思想是非平穩(wěn)變量可以表示為包含單位根的自回歸模型的形式。ADF檢驗(yàn)的檢驗(yàn)方程為：

式中，p為滯后階數(shù)，通常根據(jù)最小化Akaike信息準(zhǔn)則（Akaike information criterion，AIC）[26]的原則確定；Δxt為xt的差分；β1為均值；β2t為線性趨勢(shì)；π、γ1、…、γp為常系數(shù)；εt為殘差。使用最小二乘法對(duì)式(1)進(jìn)行參數(shù)估計(jì)，通過比較式(1)中的系數(shù)π的t統(tǒng)計(jì)量與給定顯著性水平下的t分布臨界值，從而判斷序列xt是否平穩(wěn)。具體的檢驗(yàn)流程可以參閱文獻(xiàn)[25]。

在確定變量的平穩(wěn)性后，接下來需要使用協(xié)整檢驗(yàn)來判定非平穩(wěn)的變量間是否存在協(xié)整關(guān)系并估計(jì)出協(xié)整向量。記Xt為n×1的多元非平穩(wěn)變量，如果存在一個(gè)非零向量β∈Rn，使得ξt=βTXt是平穩(wěn)的，則稱Xt是協(xié)整的，β是協(xié)整向量，ξt是新息變量。對(duì)于一階單整多變量的協(xié)整檢驗(yàn)，通常使用對(duì)參數(shù)進(jìn)行極大似然估計(jì)的Johansen檢驗(yàn)[27-28]方法。建立關(guān)于 n維一階單整向量 Xt的向量自回歸（vector autoregression，VAR）模型：

式中，Π1、…、Πp均為n×n的自回歸系數(shù)矩陣；殘差εt為獨(dú)立同分布的零均值協(xié)方差平穩(wěn)過程；滯后階數(shù)p通常根據(jù)最小化AIC的原則確定：

式(2)可以改寫成關(guān)于 Xt的向量誤差修正模型（vector error correction model，VECM）的形式：

Johansen檢驗(yàn)使用最小二乘法對(duì)以下兩式進(jìn)行參數(shù)估計(jì)：

得到式(5)和式(6)中殘差 ut的估計(jì)，分別記為和，滿足關(guān)系：

最后，求解式(9)的特征方程

對(duì)求解得到的特征值按照 λ1≥λ2≥…≥λn的順序排列。Johansen檢驗(yàn)根據(jù)特征值的跡統(tǒng)計(jì)量和最大特征值統(tǒng)計(jì)量檢驗(yàn)協(xié)整關(guān)系的個(gè)數(shù)，檢驗(yàn)假設(shè)為：

對(duì)應(yīng)的跡統(tǒng)計(jì)量和最大特征值統(tǒng)計(jì)量的定義分別如式(10)和式(11)所示：

其中，m是樣本數(shù)，λi是解式(9)得到的第i個(gè)最大特征值。Johansen檢驗(yàn)從不存在協(xié)整關(guān)系的零假設(shè)H0(r=0)開始，如果拒絕零假設(shè)，則逐次增加零假設(shè)中的協(xié)整關(guān)系個(gè)數(shù)，檢驗(yàn)至多存在一個(gè)協(xié)整關(guān)系的零假設(shè)H0(r=1)。重復(fù)上述步驟，直到接受零假設(shè)停止檢驗(yàn)，此時(shí)零假設(shè)中r的值即變量間的協(xié)整關(guān)系個(gè)數(shù)。式(9)的前r個(gè)最大特征值λ1,…,λr所對(duì)應(yīng)的特征向量構(gòu)成了協(xié)整系數(shù)矩陣β的列向量。進(jìn)一步，調(diào)整系數(shù)矩陣α的最大似然估計(jì)表達(dá)式為：

具體的檢驗(yàn)流程及相關(guān)公式的推導(dǎo)可以參閱文獻(xiàn)[27-28]。

1.2 共同趨勢(shì)表達(dá)

Stock等[24]提出了共同趨勢(shì)表達(dá)（common trends representation）方程。共同趨勢(shì)表達(dá)是共同趨勢(shì)模型的理論基礎(chǔ)，它表明：存在協(xié)整關(guān)系的一階單整非平穩(wěn)多元變量的各個(gè)分量可以分解為一階單整的共同趨勢(shì)成分和平穩(wěn)成分之和的形式，其中共同趨勢(shì)成分由若干個(gè)被稱作共同因子的一階單整過程驅(qū)動(dòng)生成。對(duì)于多元變量Xt∈Rn，Xt的每個(gè)分量是一階單整的，則相應(yīng)的 VAR模型的特征多項(xiàng)式Π(z)可以變形為：

即多項(xiàng)式分解成單位根部分和不含單位根的平穩(wěn)部分。因此形如式(4)的誤差修正模型可以改寫為：

式(14)中，L是滯后算子，C(L)=[Π*(L)]?1，，C(0)是n×n的單位矩陣。式(14)即形為式(4)的誤差修正模型對(duì)應(yīng)的向量移動(dòng)平均（vector moving average, VMA）的表達(dá)形式。令C(L)= C(1)+(1?L)C*(L)，代入式(14)做積分，不考慮Xt中的確定性成分，得到如下的關(guān)于Xt的向量移動(dòng)平均表達(dá)式：

如果變量Xt間存在r個(gè)協(xié)整關(guān)系，記β為n×r的協(xié)整矩陣，β的每一列為Xt的一個(gè)協(xié)整向量，則βTXt是平穩(wěn)的，因此βTC(1)=0，且C(1)的秩為n?r，即存在n×(n?r)的列滿秩矩陣H1、H2，使得C(1)可以分解為C(1)=H1H2T的形式。因此，式(15)可以改寫成如下的形式：

式(16)中，H1τt稱為Xt的共同隨機(jī)趨勢(shì)成分，簡稱共同趨勢(shì)成分，τt是一個(gè)(n?r)×1的向量，其各個(gè)分量可被視作隨機(jī)游走過程，因而Xt各變量的共同趨勢(shì)成分是一階單整的非平穩(wěn)過程，St=C*(L)εt稱為Xt的平穩(wěn)成分，其各變量是協(xié)方差平穩(wěn)的。式(16)被稱為Stock-Watson共同趨勢(shì)表達(dá)。共同趨勢(shì)表達(dá)方程表明：對(duì)于存在r個(gè)協(xié)整關(guān)系的n維一階單整變量，其各變量可以分解成非平穩(wěn)的共同趨勢(shì)成分和平穩(wěn)成分之和的形式，其中共同趨勢(shì)成分由n?r個(gè)一階單整過程驅(qū)動(dòng)生成。

由式(16)可以得到共同趨勢(shì)模型的一般形式：

式(17)中，ft被稱為共同因子，其各個(gè)分量是一階單整的，是Xt非平穩(wěn)并且存在協(xié)整關(guān)系的原因；是平穩(wěn)成分，其各個(gè)分量是二階平穩(wěn)的。

1.3 Gonzalo-Granger分解

共同趨勢(shì)表達(dá)證明了存在協(xié)整關(guān)系的多元一階單整變量存在共同因子，但并沒有直接提供辨識(shí)出變量的共同趨勢(shì)成分和平穩(wěn)成分的方法。為了有效地識(shí)別出共同趨勢(shì)成分和平穩(wěn)成分，必須對(duì)共同趨勢(shì)成分或平穩(wěn)成分的性質(zhì)施加一定的限制條件?；诓煌南拗茥l件，不同的共同趨勢(shì)模型[29-30]被提出。目前，Gonzalo-Granger分解[30]是應(yīng)用最廣泛的共同趨勢(shì)模型之一。

Gonzalo等[30]提出，對(duì)于式(17)，使用以下兩個(gè)限制條件可以辨識(shí)出共同因子ft：①ft是Xt的線性組合；②式(17)可以被看成永久-瞬態(tài)（permanent-transitory, P-T）分解的形式，Aft和分別構(gòu)成了Xt的永久成分和瞬態(tài)成分。

對(duì)于各變量均是一階單整的 n維非平穩(wěn)變量Xt，Xt的各變量間存在r個(gè)協(xié)整關(guān)系，定義α, β∈Rn×r分別為調(diào)整系數(shù)矩陣和協(xié)整系數(shù)矩陣，則形如式(4)的誤差修正模型可以改寫成：

由條件①可知，共同因子ft可以表示為：

其中B是一個(gè)(n?r)×n的變換矩陣。將式(19)代入式(17)，則=(I?AB)Xt。由于是平穩(wěn)的，可以表示成=A1βTXt的形式。

根據(jù)條件②，Gonzalo等[30]定義了一階單整序列Xt的P-T分解形式：①Xt=Pt+Tt，其中Pt是永久成分，Tt是瞬態(tài)成分；②Pt的各個(gè)分量是一階單整的，Tt的各分量是平穩(wěn)的；③令

為(ΔPt,Tt)的自回歸表達(dá)，并且uPt和uTt是不相關(guān)的，則

其中Et表示對(duì)過去歷史的條件期望。該限制條件表明只有永久成分Pt的新息項(xiàng)uPt對(duì)Xt有長期影響。結(jié)合式(18)，由于ft為Xt的線性組合，要使得對(duì)Xt沒有長期影響，即ΔXt的線性組合不包含 βTXt?1的成分，則唯一的線性組合為，其中α⊥是α的正交補(bǔ)矩陣。對(duì)矩陣(α⊥, β)T求逆，得到關(guān)于 Xt的P-T分解表達(dá)式：

式(20)中，等號(hào)右邊第1項(xiàng)表示Xt的永久成分，即共同趨勢(shì)模型中的共同趨勢(shì)成分，等式右邊第2項(xiàng)表示 Xt的瞬態(tài)成分，即共同趨勢(shì)模型中的平穩(wěn)成分。式(20)即Gonzalo-Granger分解共同趨勢(shì)模型。

在上述共同趨勢(shì)模型研究工作的基礎(chǔ)上，Escribano等[31]進(jìn)一步系統(tǒng)地論證了協(xié)整模型、Stock-Watson共同趨勢(shì)表達(dá)和Gonzalo-Granger分解之間的等價(jià)關(guān)系，表明①共同因子是非平穩(wěn)變量間存在協(xié)整關(guān)系的原因，②平穩(wěn)成分體現(xiàn)了變量間的全部動(dòng)態(tài)均衡關(guān)系。相關(guān)結(jié)論和具體證明可以參閱文獻(xiàn)[31]。

1.4 數(shù)值仿真

為了驗(yàn)證共同趨勢(shì)模型的有效性，應(yīng)用Gonzalo-Granger分解來辨識(shí)仿真生成的多元一階單整過程數(shù)據(jù)的平穩(wěn)成分。根據(jù)共同趨勢(shì)模型理論，共同趨勢(shì)成分是一階單整的，其2階矩是隨時(shí)間變化的，而平穩(wěn)成分是協(xié)方差平穩(wěn)的，其2階矩是時(shí)不變的，因此可以通過比較平穩(wěn)成分各分量的實(shí)際標(biāo)準(zhǔn)差和由共同趨勢(shì)模型辨識(shí)出的平穩(wěn)成分各分量的標(biāo)準(zhǔn)差來衡量共同趨勢(shì)模型的分解效果。

仿真數(shù)據(jù)為包含 3個(gè)變量的多變量非平穩(wěn)過程，每個(gè)變量均受到同一個(gè)隨機(jī)游走過程影響。該過程如下：

其中εt,,均服從均值為0，標(biāo)準(zhǔn)差為0.5的正態(tài)分布， ft是隨機(jī)游走過程，和均為平穩(wěn)的AR(1)過程。

由于仿真數(shù)據(jù)的隨機(jī)性，為了保證結(jié)果的可靠性，仿真數(shù)據(jù)的樣本長度為5000，仿真10000次取平均值。每次仿真時(shí)，先對(duì)Xt進(jìn)行Johansen檢驗(yàn)，估計(jì)出協(xié)整系數(shù)矩陣β和調(diào)整系數(shù)矩陣α，再將β和α代入Gonzalo-Granger分解的表達(dá)式(20)，得到各變量的共同趨勢(shì)成分和平穩(wěn)成分。

表1 仿真數(shù)據(jù)各變量的平穩(wěn)成分的標(biāo)準(zhǔn)差Table 1 Standard derivations of simulated stationary components

由表1可知，仿真實(shí)驗(yàn)中，Gonzalo-Granger分解辨識(shí)出的各變量平穩(wěn)成分的標(biāo)準(zhǔn)差非常接近實(shí)際值，說明該模型能夠有效地識(shí)別出平穩(wěn)成分。因此，選擇使用 Gonzalo-Granger分解模型用于非平穩(wěn)過程監(jiān)控。

2 基于共同趨勢(shì)模型的非平穩(wěn)過程在線監(jiān)控方法

根據(jù)共同趨勢(shì)模型的相關(guān)理論，存在協(xié)整關(guān)系的多元非平穩(wěn)過程變量可以分解成非平穩(wěn)的共同趨勢(shì)成分和平穩(wěn)成分，其中平穩(wěn)成分體現(xiàn)了變量間全部的動(dòng)態(tài)均衡關(guān)系。由于平穩(wěn)成分是協(xié)方差平穩(wěn)的，因此可以將平穩(wěn)成分提取出來，消除非平穩(wěn)共同趨勢(shì)成分的影響，應(yīng)用經(jīng)典的多元統(tǒng)計(jì)方法對(duì)其進(jìn)行監(jiān)控。對(duì)于存在協(xié)整關(guān)系的非平穩(wěn)工業(yè)過程，先使用共同趨勢(shì)模型提取出正常工況下各非平穩(wěn)過程變量的平穩(wěn)成分。當(dāng)工業(yè)過程系統(tǒng)發(fā)生故障，變量間的動(dòng)態(tài)均衡關(guān)系被破壞，此時(shí)，依照正常模型從故障數(shù)據(jù)提取出的“平穩(wěn)成分”不再服從正常工況時(shí)平穩(wěn)成分的統(tǒng)計(jì)規(guī)律?；谏鲜鏊枷?，提出了基于共同趨勢(shì)模型的非平穩(wěn)過程監(jiān)控方法。與文獻(xiàn)[22]中基于協(xié)整新息變量的監(jiān)控方法的區(qū)別在于，本方法是對(duì)非平穩(wěn)變量的平穩(wěn)成分進(jìn)行監(jiān)控，而后者是對(duì)非平穩(wěn)變量協(xié)整模型的新息變量進(jìn)行監(jiān)控，但是新息變量不能全面地反映變量間的動(dòng)態(tài)均衡關(guān)系。相似之處在于新息變量與平穩(wěn)成分都是二階平穩(wěn)的，兩者可以使用相同的統(tǒng)計(jì)指標(biāo)。本文使用的監(jiān)控指標(biāo)是Hotelling的T2統(tǒng)計(jì)量[32]，與文獻(xiàn)[22]相同。

2.1 監(jiān)控統(tǒng)計(jì)量

對(duì)于多元變量監(jiān)控，Hotelling的T2統(tǒng)計(jì)量是最常用的監(jiān)控指標(biāo)之一。對(duì)于數(shù)據(jù)集 X∈Rn×m，其中n表示變量維度，m表示樣本個(gè)數(shù)，為了消除變量尺度不同帶來的影響，先對(duì)X做歸一化處理，即每個(gè)樣本分別減去變量的樣本均值再除以相應(yīng)的樣本標(biāo)準(zhǔn)差。記Z∈Rn×m為標(biāo)準(zhǔn)化的數(shù)據(jù)集，計(jì)算Z的樣本協(xié)方差矩陣：C=ZZT/(m?1)，接著對(duì)協(xié)方差矩陣進(jìn)行特征值分解，并按照特征值的大小進(jìn)行降序排列：C=PΛPT，其中Λ表示k個(gè)最大特征值構(gòu)成的對(duì)角矩陣，即Λ=diag(λ1, λ2,…,λk)，P為由k個(gè)特征值對(duì)應(yīng)的特征向量構(gòu)成的載荷矩陣。在確定k的大小時(shí)，通常采用累積方差貢獻(xiàn)率方法。累積貢獻(xiàn)率定義為，通常取累積方差貢獻(xiàn)率不小于 85%即認(rèn)為包含了原始數(shù)據(jù)足夠多的信息。對(duì)于每個(gè)樣本z，按下式計(jì)算Hotelling的T2統(tǒng)計(jì)量：

在數(shù)據(jù)服從協(xié)方差平穩(wěn)分布的假設(shè)下，Hotelling的T2統(tǒng)計(jì)量的控制限由F分布估計(jì)：

式(22)中，F(xiàn)α(k, m?k)表示在置信水平α上自由度為k和(m?k)的F分布。

2.2 監(jiān)控流程

基于共同趨勢(shì)模型的非平穩(wěn)過程監(jiān)控流程如圖1所示。

圖1 基于共同趨勢(shì)模型的非平穩(wěn)過程監(jiān)測(cè)方法Fig.1 Non-stationary process monitoring using common trends models

2.2.1 建立正常工況下的模型

（1）對(duì)正常數(shù)據(jù)集樣本的各個(gè)變量及其一階差分變量進(jìn)行 ADF單位根檢驗(yàn)，選取出一階單整變量Xt=(x1t, x2t, …, xnt)T，將它們作為建模變量。

（2）由式(3)的AIC準(zhǔn)則確定Xt的VAR模型的滯后階數(shù)，再依據(jù)式(5)～式(12)及相應(yīng)流程，對(duì)Xt進(jìn)行Johansen檢驗(yàn)，確定Xt各變量間的協(xié)整關(guān)系個(gè)數(shù)，并計(jì)算出協(xié)整系數(shù)矩陣β和調(diào)整系數(shù)矩陣α。

（3）根據(jù) Gonzalo-Granger分解表達(dá)式(20)，得到 Xt中的平穩(wěn)成分 St=WXt，其中投影矩陣W=α(βTα)?1βT。

（4）對(duì)平穩(wěn)成分St進(jìn)行歸一化處理，使得歸一化后的平穩(wěn)成分的各分量為零均值和單位方差。

（5）選擇合適的方差累積貢獻(xiàn)率，計(jì)算歸一化的平穩(wěn)成分的載荷矩陣P和對(duì)角線矩陣Λ，并由式(22)計(jì)算出監(jiān)控指標(biāo)Hotelling的T2統(tǒng)計(jì)量的控制限。

2.2.2 在線監(jiān)控

（1）獲取新的實(shí)時(shí)監(jiān)控?cái)?shù)據(jù)，根據(jù)建模階段步驟(1)的結(jié)果選取出數(shù)據(jù)中的非平穩(wěn)過程變量xt。

（2）根據(jù)建模階段步驟(3)獲得的投影矩陣求得xt中的平穩(wěn)成分st=Wxt。

（3）用建模階段步驟(4)得到的均值和方差對(duì)平穩(wěn)成分st進(jìn)行歸一化處理。

（4）將建模階段步驟(5)得到的載荷矩陣P和對(duì)角線矩陣Λ代入式(21)，計(jì)算出歸一化平穩(wěn)成分的Hotelling的T2統(tǒng)計(jì)量。

（5）監(jiān)控Hotelling的T2統(tǒng)計(jì)量是否超限。

3 案例研究

3.1 石油蒸餾過程

為了驗(yàn)證本文提出的監(jiān)控方法的有效性，將該方法應(yīng)用于某個(gè)非平穩(wěn)的石油蒸餾過程的在線監(jiān)控。該蒸餾過程用于精餾烴的混合物，出于產(chǎn)品質(zhì)量和安全因素考慮，共有16個(gè)過程變量被監(jiān)測(cè)。關(guān)于該過程的詳細(xì)信息可以參閱文獻(xiàn)[22]。該蒸餾過程的工藝流程圖如圖2所示。

圖2 石油蒸餾過程的工藝流程Fig.2 Process flow sheet of distillation process

本文監(jiān)控的故障是新鮮原料的異常流動(dòng)，與文獻(xiàn)[22]相同。在該過程中，蒸餾塔頂部和底部的產(chǎn)品濃度不受反饋控制，而新鮮原料的流速一直在變化。因此，原料或原料成分的任何異常減少或增加都會(huì)導(dǎo)致產(chǎn)品濃度發(fā)生變化，降低產(chǎn)品質(zhì)量。因此必須及時(shí)監(jiān)測(cè)到新鮮原料流動(dòng)異常變化的發(fā)生。

樣本數(shù)據(jù)的采樣間隔是 30 s，每個(gè)樣本由 16個(gè)過程變量組成。正常數(shù)據(jù)集包含8000個(gè)樣本，作為訓(xùn)練集用來建立監(jiān)控模型。故障數(shù)據(jù)集包含大約4700個(gè)樣本，用來驗(yàn)證本文提出方法的有效性。

3.2 監(jiān)控模型的建立

首先對(duì)正常數(shù)據(jù)集的16個(gè)過程變量按式(1)的檢驗(yàn)方程進(jìn)行 ADF檢驗(yàn)以確定變量的平穩(wěn)性，滯后階數(shù)p由AIC準(zhǔn)則確定，相應(yīng)的檢驗(yàn)結(jié)果見表2。

表2 正常數(shù)據(jù)集過程變量的ADF檢驗(yàn)結(jié)果Table 2 ADF test results of normal process variables

由ADF檢驗(yàn)結(jié)果可知，在10%的置信水平上，存在4個(gè)一階單整的非平穩(wěn)變量，這與文獻(xiàn)[22]給出的變量平穩(wěn)性檢驗(yàn)結(jié)果一致。接著，對(duì)這4個(gè)非平穩(wěn)變量進(jìn)行Johansen檢驗(yàn)，以確定變量間存在的協(xié)整關(guān)系個(gè)數(shù)，并估計(jì)出協(xié)整系數(shù)矩陣和調(diào)整系數(shù)矩陣。對(duì)上述非平穩(wěn)變量建立形如式(2)的VAR模型，由式(3)的AIC準(zhǔn)則確定其滯后階數(shù)為13，不同滯后階數(shù)對(duì)應(yīng)的AIC見圖3。Johansen跡檢驗(yàn)和最大特征值檢驗(yàn)的結(jié)果分別見表3和表4。

圖3 由AIC確定的VAR模型的最佳滯后階數(shù)Fig.3 Optimal VAR lag order determined by AIC

Johansen檢驗(yàn)的結(jié)果說明，非平穩(wěn)變量間存在2個(gè)協(xié)整關(guān)系。在按照式(5)～式(12)計(jì)算出相應(yīng)的協(xié)整系數(shù)矩陣β和調(diào)整系數(shù)矩陣α之后，根據(jù)式(20)的Gonzalo-Granger分解，提取出各個(gè)非平穩(wěn)變量的平穩(wěn)成分St和共同趨勢(shì)成分Nt。為了驗(yàn)證提取出的共同趨勢(shì)成分和平穩(wěn)成分的平穩(wěn)性，對(duì)St和Nt的各分量進(jìn)行ADF檢驗(yàn)。St各分量的檢驗(yàn)結(jié)果見表5，Nt各分量的檢驗(yàn)結(jié)果見表6。檢驗(yàn)結(jié)果顯示，St的各個(gè)分量均是平穩(wěn)的，Nt的各個(gè)分量均是一階單整的，表明共同趨勢(shì)模型用于非平穩(wěn)工業(yè)過程數(shù)據(jù)建模是有效的。

表3 Johansen跡檢驗(yàn)結(jié)果Table 3 Johansen trace test results

表4 Johansen最大特征值檢驗(yàn)結(jié)果Table 4 Johansen max-eigenvalue test results

表5 平穩(wěn)成分的ADF檢驗(yàn)結(jié)果Table 5 ADF test results of stationary components

接著對(duì)從非平穩(wěn)過程中提取出的平穩(wěn)成分做歸一化處理。計(jì)算歸一化后的平穩(wěn)成分的載荷矩陣P和對(duì)角線矩陣Λ，由于累積方差貢獻(xiàn)率的選取應(yīng)當(dāng)不小于85%，因此相應(yīng)的k值選取為2。再根據(jù)式(22)計(jì)算出監(jiān)控指標(biāo)的控制限，控制限的置信水平為 99%，與文獻(xiàn)[22]一致。至此，正常工況下基于共同趨勢(shì)模型的非平穩(wěn)過程監(jiān)控模型被建立。

表6 共同趨勢(shì)成分的ADF檢驗(yàn)結(jié)果Table 6 ADF test results of common trends components

3.3 故障監(jiān)控分析與討論

故障數(shù)據(jù)表示的是在第1800和第2900個(gè)樣本點(diǎn)處發(fā)生了兩次新鮮原料減少。在第1次原料減少發(fā)生后，設(shè)備操作人員注意到了此次異常的發(fā)生，通過降低回流流速抵消了此次異常的影響。但是，當(dāng)?shù)?次原料減少發(fā)生后，設(shè)備操作人員未能及時(shí)注意到此次異常的發(fā)生，直到在大約第4500個(gè)樣本點(diǎn)處，此時(shí)丁烷產(chǎn)品中的雜質(zhì)含量有了顯著增加，設(shè)備操作人員才再一次降低了回流流速，抵消了第2次原料減少造成的影響[22]。

將正常模型的協(xié)整系數(shù)矩陣β和調(diào)整系數(shù)矩陣α代入式(20)，提取出故障狀態(tài)下非平穩(wěn)過程變量的“平穩(wěn)成分”。接著將正常模型的載荷矩陣 P和對(duì)角線矩陣Λ代入式(21)，計(jì)算出從歸一化后的“平穩(wěn)成分”的Hotelling的T2統(tǒng)計(jì)量，如圖4所示。圖4中，實(shí)橫線表示控制限。分析監(jiān)控結(jié)果，在第1800個(gè)樣本點(diǎn)處第1次原料減少發(fā)生后，監(jiān)控指標(biāo)超出了控制限，及時(shí)做出了預(yù)警，由于設(shè)備操作人員降低了回流流量，抵消了這次減少造成的影響，監(jiān)控指標(biāo)回落到正常水平。在第2900個(gè)樣本點(diǎn)處新鮮原料再次減少，這次減少未引起操作人員的注意，但監(jiān)控指標(biāo)再次超過了控制限，及時(shí)做出了預(yù)警，證明了本文提出的監(jiān)控方法的有效性。

圖4 基于共同趨勢(shì)模型的監(jiān)測(cè)結(jié)果Fig.4 Monitoring results for fresh feed drops using common trends model

圖5 基于文獻(xiàn)[22]方法的監(jiān)測(cè)結(jié)果Fig.5 Monitoring results for fresh feed drops using methods in Ref.[22]

最后，將本文提出方法的監(jiān)控結(jié)果與文獻(xiàn)[22]方法的監(jiān)控結(jié)果進(jìn)行比較。文獻(xiàn)[22]提出了一種基于協(xié)整新息變量的非平穩(wěn)過程監(jiān)控方法，對(duì)非平穩(wěn)過程變量建立協(xié)整模型，估計(jì)出一組協(xié)整向量，得到新息變量，對(duì)新息變量進(jìn)行逆自回歸濾波，得到新息變量的自回歸殘差。圖5展示了文獻(xiàn)[22]的監(jiān)控結(jié)果，自上至下依次是原始新息變量自回歸殘差、新息變量的自回歸殘差及補(bǔ)償、新息變量去除補(bǔ)償后的殘差的Hotelling的T2統(tǒng)計(jì)量控制圖。對(duì)于第1800個(gè)樣本點(diǎn)處的第1次原料減少，文獻(xiàn)[22]和本文提出的方法均監(jiān)測(cè)到了此次異常的發(fā)生，及時(shí)地做出了預(yù)警。對(duì)于第2900個(gè)樣本點(diǎn)處的第2次原料減少，雖然文獻(xiàn)[22]的3組監(jiān)測(cè)指標(biāo)均監(jiān)測(cè)到了此次異常的發(fā)生，但是只有圖5(b)的監(jiān)測(cè)指標(biāo)及時(shí)做出了預(yù)警，其他兩組監(jiān)測(cè)指標(biāo)均有明顯的滯后性，直到第3400個(gè)樣本點(diǎn)附近才做出預(yù)警，而本文提出的方法則及時(shí)地做出了預(yù)警。此外，文獻(xiàn)[22]的 3組監(jiān)控指標(biāo)在抗干擾方面存在不足，3組監(jiān)控指標(biāo)在第1000個(gè)樣本處均做出了明顯的誤報(bào)，在第2200個(gè)樣本點(diǎn)附近，圖5(b)的監(jiān)控指標(biāo)也做出了明顯的誤報(bào)[22]，相比之下，本文提出的方法沒有出現(xiàn)明顯的誤報(bào)，具有更好的抗干擾性。

4 結(jié) 論

針對(duì)傳統(tǒng)的基于平穩(wěn)假設(shè)的多元統(tǒng)計(jì)方法不適用于非平穩(wěn)過程監(jiān)控的困境，提出了基于共同趨勢(shì)模型的非平穩(wěn)過程監(jiān)控方法。該方法首先應(yīng)用共同趨勢(shì)模型提取出各非平穩(wěn)變量中的平穩(wěn)成分，然后再對(duì)平穩(wěn)成分應(yīng)用Hotelling的T2統(tǒng)計(jì)量進(jìn)行狀態(tài)監(jiān)測(cè)。與現(xiàn)有的基于協(xié)整新息變量的非平穩(wěn)過程監(jiān)控方法相比，共同趨勢(shì)模型可以有效地提取出非平穩(wěn)變量中的平穩(wěn)成分，消除非平穩(wěn)的共同趨勢(shì)成分的影響，平穩(wěn)成分包含了反映系統(tǒng)非平穩(wěn)變量間全部的動(dòng)態(tài)均衡關(guān)系的信息。對(duì)石油蒸餾過程監(jiān)控的結(jié)果表明，基于共同趨勢(shì)模型的監(jiān)控方法好于基于協(xié)整新息變量的監(jiān)控方法。

然而，協(xié)整及共同趨勢(shì)模型的相關(guān)理論近年來仍在不斷發(fā)展。進(jìn)一步的研究包括：①改進(jìn)協(xié)整系數(shù)矩陣的計(jì)算方法；②非平穩(wěn)過程變量與平穩(wěn)過程變量的聯(lián)合監(jiān)控。

[1]陳前. 關(guān)于工程大系統(tǒng)的狀態(tài)監(jiān)測(cè)與故障診斷[J]. 振動(dòng)、測(cè)試與診斷, 2002, 22(3): 163-230. CHEN Q. Condition monitoring and fault diagnosis of large-scale engineering system[J]. Journal of Vibration Measurement and Diagnosis, 2002, 22(3): 163-230.

[2]GE Z, SONG Z, GAO F. Review of recent research on data-based process monitoring [J]. Industrial & Engineering Chemistry Research, 2013, 52(10): 3543-3562.

[3]CHEN Q, WYNNE R J, GOULDING P, et al. The application of principal component analysis and kernel density estimation to enhance process monitoring [J]. Control Engineering Practice, 2000, 8(5): 531-543.

[4]童楚東, 史旭華. 基于互信息的PCA方法及其在過程監(jiān)測(cè)中的應(yīng)用[J]. 化工學(xué)報(bào), 2015, 66(10): 4101-4106. TONG C D, SHI X H. Mutual information based PCA algorithm with application in process monitoring [J]. CIESC Journal, 2015, 66(10): 4101-4106.

[5]KANO M, TANAKA S, HASEBE S, et al. Monitoring independent components for fault detection [J]. AIChE Journal, 2003, 49(4): 969-976.

[6]GE Z, SONG Z. Performance-driven ensemble learning ICA model for improved non-Gaussian process monitoring [J]. Chemometrics and Intelligent Laboratory Systems, 2013, 123: 1-8.

[7]HUANG J, YAN X. Gaussian and non-Gaussian double subspace statistical process monitoring based on principal component analysis and independent component analysis [J]. Industrial & Engineering Chemistry Research, 2015, 54(3): 1015-1027.

[8]衷路生, 何東, 龔錦紅, 等. 基于分布式 ICA-PCA模型的工業(yè)過程故障監(jiān)測(cè)[J]. 化工學(xué)報(bào), 2015, 66(11): 4546-4554. ZHONG L S, HE D, GONG J H, et al. Fault monitoring of industrial process based on distributed ICA-PCA model [J]. CIESC Journal, 2015, 66(11): 4546-4554.

[9]GUNTHER J C, CONNER J S, SEBORG D E. Process monitoring and quality variable prediction utilizing PLS in industrial fed-batch cell culture [J]. Journal of Process Control, 2009, 19(5): 914-921.

[10]韓敏, 張占奎. 基于改進(jìn)核主成分分析的故障檢測(cè)與診斷方法[J].化工學(xué)報(bào), 2015, 66(6): 2139-2149. HAN M, ZHANG Z K. Fault detection and diagnosis method based on modified kernel principal component analysis [J]. CIESC Journal, 2015, 66(6): 2139-2149.

[11]LEE J M, YOO C K, CHOI S W, et al. Nonlinear process monitoring using kernel principal component analysis [J]. Chemical Engineering Science, 2004, 59(1): 223-234.

[12]LEE J M, QIN S J, LEE I B. Fault detection of non-linear processes using kernel independent component analysis [J]. The Canadian Journal of Chemical Engineering, 2007, 85(4): 526-536.

[13]ZHANG Y. Enhanced statistical analysis of nonlinear processes using KPCA, KICA and SVM [J]. Chemical Engineering Science, 2009, 64(5): 801-811.

[14]胡益, 王麗, 馬賀賀, 等. 基于核PLS方法的非線性過程在線監(jiān)控[J]. 化工學(xué)報(bào), 2011, 62(9): 2555-2561. HU Y, WANG L, MA H H, et al. Online nonlinear process monitoring using kernel partial least squares [J]. CIESC Journal, 2011, 62(9): 2555-2561.

[15]GODOY J L, ZUMOFFEN D A, VEGA J R, et al. New contributions to non-linear process monitoring through kernel partial least squares [J]. Chemometrics and Intelligent Laboratory Systems, 2014, 135: 76-89.

[16]COMON P, JUTTEN C. Handbook of Blind Source Separation: Independent Component Analysis and Applications [M]. Academic Press, 2010.

[17]BERTHOUEX P M, BOX G E. Time series models for forecasting wastewater treatment plant performance [J]. Water Research, 1996, 30(8): 1865-1875.

[18]ENGLE R F, GRANGER C W J. Co-integration and error correction: representation, estimation, and testing [J]. Econometrica, 1987, 55(2): 251-276.

[19]PFAFF B. Analysis of Integrated and Cointegrated Time Series with R [M]. Springer Science & Business Media, 2008.

[20]ENDERS W. Applied Econometric Time Series [M]. John Wiley & Sons, 2008.

[21]陳前, 潘昱昱. 協(xié)整理論應(yīng)用于非平穩(wěn) FCCU系統(tǒng)的狀態(tài)監(jiān)測(cè)與故障診斷[J]. 石油學(xué)報(bào)（石油加工）, 2007, 23(1): 69-76. CHEN Q, PAN Y Y. Application of cointegration testing method to condition monitoring and fault diagnosis of nonstationary FCCU system [J]. Acta Petrolet Sinica (Petroleum Processing Section), 2007, 23(1): 69-76.

[22]CHEN Q, KRUGER U, LEUNG A Y T. Cointegration testing method for monitoring nonstationary processes [J]. Industrial & Engineering Chemistry Research, 2009, 48(7): 3533-3543.

[23]石海忱, 陳前, 林原靈. 協(xié)整系數(shù)矩陣的非平穩(wěn)工程系統(tǒng)故障診斷應(yīng)用研究[J]. 振動(dòng)與沖擊, 2015, 34(1): 146-150. SHI H C, CHEN Q, LIN Y L. Fault diagnosis of non-stationary engineering system using cointegration cofficients matrix [J]. Journal of Vibration and Shock, 2015, 34(1): 146-150.

[24]STOCK J H, WATSON M W. Testing for common trends [J]. Journal of the American statistical Association, 1988, 83(404): 1097-1107.

[25]DICKEY D A, FULLER W A. Distribution of the estimators for autoregressive time series with a unit root [J]. Journal of the American statistical association, 1979, 74(366a): 427-431.

[26]AKAIKE H. A new look at the statistical model identification [J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716-723.

[27]JOHANSEN S, JUSELIUS K. Maximum likelihood estimation and inference on cointegration - with applications to the demand for money [J]. Oxford Bulletin of Economics and statistics, 1990, 52(2): 169-210.

[28]JOHANSEN S. Estimation and hypothesis testing of cointegration vectors in Gaussian vector autoregressive models [J]. Econometrica, 1991, 59(59): 1551-1580.

[29]KASA K. Common stochastic trends in international stock markets [J]. Journal of monetary Economics, 1992, 29(1): 95-124.

[30]GONZALO J, GRANGER C. Estimation of common long-memory components in cointegrated systems [J]. Journal of Business and Economic Statistics, 1995, 13(1): 27-35.

[31]ESCRIBANO A, PE?A D. Cointegration and common factors [J]. Journal of Time Series Analysis, 1994, 15(6): 577-586.

[32]HOTELLING H. The Generalization of Student's Ratio [M]. New York: Springer, 1992.

Online non-stationary process monitoring by common trends model

LIN Yuanling, CHEN Qian
(State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, Jiangsu, China)

Non-stationary process monitoring based on common trends model was proposed, because conventional multivariate statistical process control methods with stationary data assumption were inapplicable to non-stationary process monitoring. The new common trends model was capable of identifying common factors from co-integrated non-stationary multiple variables and decomposing each non-stationary process variable into summation of a non-stationary common trends component and a stationary counterpart. Contrary to existing non-stationary process monitoring technique from cointegration model, the common trends model captured the stationary component of each non-stationary process variable, eliminated effects of non-stationary common factors and unveiled overall dynamic equilibrium relationships among variables. Hence, non-stationary process monitoring was transformed to an application of common trends model, which involved obtaining stationary component of each process variable, creating estimation for the stationary components by conventional multivariate statistical methods and setting up monitoring on corresponding control limits. A case study of monitoring petroleum distillation process showed that the proposed approach possessed more reliable process monitoring performance than the method of cointegration model.

common trends model; process monitoring; non-stationary process; cointegration testing; process control; process systems; system engineering

Prof. CHEN Qian, q.chen@nuaa.edu.cn

TP 277

：A

：0438—1157（2017）01—0178—10

10.11949/j.issn.0438-1157.20160670

2016-05-16收到初稿，2016-09-25收到修改稿。

聯(lián)系人：陳前。

：林原靈（1989—），男，博士研究生。

江蘇高校優(yōu)勢(shì)學(xué)科建設(shè)工程資助項(xiàng)目。

Received date: 2016-05-16.

Foundation item: supported by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).