CPV(1)分段法在間歇過程中的應用

袁 杰，郭小萍，劉如有，李 元

(1.沈陽化工大學 信息工程學院，遼寧 沈陽 110142;2.中國市政工程東北設計研究總院，吉林 長春 130021)

近幾年化工產品市場品種更新迅速，需求不斷變化，間歇過程的靈活性能適應小批量、高附加值、多品種的工業要求.因此，間歇過程將成為化工生產的重要方式.動態特性是間歇過程的本質.間歇過程本身具有多工況、分階段、階段時長不一致等特點，產品質量很容易受到外部不確定因素的影響.這些特點使間歇過程的故障檢測更加復雜.

針對間歇過程的故障診斷，多向偏最小二乘法［1］(MPLS，Muti-way partial least squares)和多向主元分析［2－3］(MPCA，Muti-way principal component analysis)是最常用的方法.MPCA 最早是由Nomikos 和Macgregor 在1995 年應用到間歇過程故障診斷中的，達到了很好的檢測和診斷目的.然而，間歇過程中存在多個操作階段，各個階段內部可能具有其他階段不具有的性質.從統計學角度分析，不同階段的數據特征可能存在很大差距.在間歇過程故障診斷中，如果統一建模，某些總體的信息可能掩蓋了局部階段內的信息變化.所以，批次過程的分段故障檢測成為一個研究熱點.目前，將間歇過程合理劃分為不同時段的方法有很多，最常用方法是自動識別算法［4］，該方法通過某種算法程序能夠自動識別出過程中的各個階段.文獻［5－7］屬于自動識別算法的應用，從主元累計貢獻率的角度進行了階段劃分和故障檢測的探究.本文在將檢測數據等長處理后，進行主元分析［8］，采用第一主元貢獻率(CPV(1))的變化來判斷原數據是否分段，也屬于自動識別算法.仿真結果表明:使用第一主元貢獻率的變化能準確反映生產階段的變化.

1 間歇過程的數據處理方法

1.1 軌跡同步化

間歇過程的時變特點使得間歇過程所采集的數據呈現批次數據軌跡不等長的特點，因此，批次過程的軌跡同步化處理十分關鍵.此處理目的是將不等長的數據處理為便于操作的等長數據.目前比較常用的軌跡同步化處理方法有:最短長度法、指示變量法(indicator variable technique，IVT)等.

(1)最短長度法

最短長度法以間歇操作數據最短的批次做為基準，其余批次的數據截取長度與基準批次相同.這種方法一方面要求在公共部分的過程變量軌跡保持一致，另一方面公共時間段內要包含間歇過程的主要操作，在超出基準長度的時間位置包含的都是數據中次重要信息.

(2)指示變量法［5］

通常情況下，間歇過程變量都對應著時間而改變，但是如果以其他的變量來替代時間，就可以對不等長數據進行同步處理.把可以用來代替時間的變量命名為“指示變量”.這種方法的思想是用指示變量的進程來代替時間，將其余變量進行重新采樣，是一個數據插值過程.

文中青霉素仿真數據的軌跡同步化采用的是指示變量法，半導體數據采用的是最短長度法.

1.2 數據展開

經過軌跡同步化之后的數據是三維數據，在采用主元分析之前，需要將其展開為二維數據.傳統的展開方式有以下兩種:基于變量的展開和基于批次的展開.這兩種展開方式如圖1 所示.

圖1 批次過程的展開方式Fig.1 The unfolding means of data

兩種數據展開方式都是針對采樣時刻片進行的操作.基于變量的展開方式是將每個采樣時刻片展開為(I×K)×J 的矩陣.基于批次的展開方式是將采樣時刻片展開為I×(J ×K)的矩陣.一般的，根據數據展開后的操作區分，在建模和檢測時采用基于批次的展開方法，在數據分段時采用基于變量的展開方式.

2 基于第一主元貢獻率差值變化的批次分段

2.1 主元分析方法

主元分析(PCA)是多元統計分析中最基本的方法之一.假設x ∈Rm，代表一個包含了m 個傳感器測量樣本，每個傳感器個有n 個獨立采樣，構造如下的測量數據矩陣X=［x1，x2，…，xn］T∈Rn×m，其中X 的每一列代表一個測量變量，每一行代表一個樣本.首先對X 進行零均值單位方差的標準化處理，定義標準化后的樣本集的協方差矩陣為

對其進行特征值分解，并且按照特征值的大小順序排列.PCA 模型對X 進行如下分解:

其中，P ∈Rm×A為負載矩陣，由S 的前A 個特征向量構成.T ∈Rn×A為得分矩陣，T 的各列被稱為主元變量，A 表示主元的個數，也是得分矩陣的列數.

現使用多向主元分析的方法研究批次過程數據，在忽略次要的非線性關系的同時，假定各個過程參數分別符合高斯分布，可以將PCA 的方法用到批次展開后的數據陣中.

2.2 第一主元貢獻率

原始矩陣X ∈Rn×m所對應的協方差矩陣是S ∈Rm×m，協方差矩陣S 對應的特征向量矩陣是Q ∈Rm×m，其中前A 個組成的矩陣就是負載矩陣，根據特征值特征向量公式:SQ=λQ 可以計算出協方差矩陣的m 個特征值:λ1，λ2，…，λm.

累積方差百分比表征的是前k 個特征值的和占總體方差總和的比率.它的計算方法如下:

其中，第一主元的累積方差百分比可以表示為:

CPV(1)表示第一主元所攜帶的信息占原始數據總信息的比例，能反映原始數據的主要信息.CPV(1)在主元分析中的作用一方面表現了第一主元提取原始數據信息量的多少，另一方面也可反映數據變量之間的相關性信息和動態關系.如果數據間的動態關系保持不變，那么CPV(1)的值將基本恒定［6］.用相鄰兩個采樣時刻對應的數據段的CPV(1)的差值表示2 個數據段所含信息量的差異.當差值最大時，表明數據段的主要信息發生了很大變化，段落由此劃分.因此，以CPV(1)差值的變化為指標來確定批次內部不同階段間動態特性的轉變.

2.3 基于第一主元貢獻率差值的分段方法

批次過程中的變量并不是全都可以用于分段算法.在批次分段過程中引入某些恒值、多變性的變量，會使分段效果不明顯，甚至使過程不能分段.所以，在對批次過程進行分段的過程中，首先要分析原始數據，篩選出適合分段的變量.篩選變量的基本原則是選取變量軌跡具有分段特點的變量，不選取波動不大和波動過于劇烈的變量.進行分段算法的具體流程如圖2 所示.

基于第一主元貢獻率變化分階段算法如下:

(1)篩選出合適變量的數據集合(I ×K ×J)，在此基礎上，得到K－1 個子階段，每個階段的數據長度為［t，t+1，…，K］(t=1，2，3，…，K－1)，其中t 為采樣時刻;

(2)對于每個子階段應用基于變量的展開方式展開為X(I(K－t+1)×J)，建立PCA 模型;

(3)每個子階段，計算CPV (1)1×t={CPV(1)t}，當t ＞2 時，計算每個子階段與前一子階段CPV(1)的變化，即d=CPV(1)t－CPV(1)t－1，得到相鄰兩個子段第一主元累計貢獻率差值集合{d(t)};

(4)當{d(t)}達到最大值時，若其所對應的采樣時刻t 滿足K－t ＞l 且t－1 ＞l，其中l 為根據經驗確定的邊緣效應值，不同生產過程l 的數值一般不同，則第一次分階段的時刻即為t 時刻，并將整個過程分為2 個階段［1，2，…，t－1］和［t，t+1，…，K］;若不滿足上述條件，則不分段.

(5)對分好的各個段重復步驟(2)到(4)，直到每個階段不可以繼續分段，結束［7］.

圖2 第一主元貢獻率差值分段流程Fig.2 Flow chart of phase-dividing using CPV(1)method

3 仿真實驗

3.1 青霉素仿真數據的分段實驗

青霉素發酵是一個典型的補料分批發酵過程，根據操作主要分為兩個部分:第一部分為發酵前期，是間歇式操作.在這個過程中菌體會消耗掉大量的葡萄糖(反應底物)，當葡萄糖的溶度過低時反應進行到第二階段，是半間歇補料操作.為了保持青霉素有較高的產率，細胞生長速率必須保持在某一最小值以上，因此，在這一操作階段葡萄糖要連續地添加到發酵反應罐中.這一操作過程又包括青霉素合成期和菌體死亡期.從青霉素菌體的生長過程看，該過程可以分為3個階段.

青霉素仿真軟件產生的是一個18 個變量的數據集.仿真實驗中，采用文獻［5］中選用和處理變量的方法.選擇這18 個變量中的第4、8、9、10、14、17 變量進行分段實驗.這幾個變量的物理含義分別為:底物供給速度、生物總量、青霉素溶量、反應體積、產熱量、冷水流量.

用pensim 軟件產生50 批次的青霉素數據.為體現間歇過程特點，所得到的50 批次數據是不等長的，即整體反應時間和每個批次內部階段工時都不等長.在間歇階段中選取變量10 為指示變量，通過體積的減小來反映發酵過程的進程.通過底物加入速度(變量4)可以計算出底物加入量占總體加入量的百分比，是遞增變量，選作半間歇過程指示變量，并假設當14 L 的葡萄糖加入到反應器時整個發酵過程結束.在間歇操作階段，隨著體積每減少1 %，其他的變量通過線性插值的方法進行重新采樣;在半間歇操作階段，葡萄糖每加入總量的0.2 %，其他變量相應地重新進行采樣.通過批次軌跡同步化后達到所有批次的采樣個數及數據長度相同，這兩個階段的數據長度分別為101 個和501 個，每一批次的數據長度都是602 個.軌跡同步化后的訓練集X變為(50 ×602 ×6)的三維矩陣.基于變量展開為二維矩陣(50 ×602)×6.

在該算法中，選取青霉素工藝過程對應的邊緣效應為25.在圖3(a)中，d 的最大值在第101采樣時刻處.在已經分好的2 個階段內部重新采用分段算法，前101 個采樣時刻所對應的第一主元貢獻率差值的變化如圖3(b)所示，第一主元貢獻率差值最大值在邊緣效應范圍內，不可以再分段.對102－602 采樣時刻進行分段，第一主元貢獻率的差值如圖3(c)所示，在第405 時刻達到最大值，可以對這一部分繼續分段.從這一階段又分出的兩部分應用該算法所得到的第一主元貢獻率差值如圖3(d)和圖3(e)所示，均不可以再分段.

圖3 青霉素數據第一主元貢獻率差值Fig.3 The differential chart of CPV(1)of pensim data

通過以上操作，青霉素602 個采樣時刻被分為3 個階段:第1 時刻到第101 時刻對應第1 階段發酵前期;第102 時刻到第405 時刻對應第2階段青霉素合成期;第406 時刻到602 時刻對應第3 階段菌體死亡期.分段結果與青霉素菌體生長的實際生物過程一致.

3.2 半導體生產過程分段實驗

半導體數據是從鋁板蝕刻機的工業生產過程中采集的［9］.共采集了107 批的正常批次數據.這些數據記錄了3 個不同模態下過程變量變化.該半導體生產過程實際可分為3 個階段.

圖4 給出了第7 和第20 個變量數值波動情況.從這兩個變量的波動軌跡可以直觀地分辨出數據的分段位置.該過程產生的數據含有21 個變量.選取利于分段的第7、9、12、15、20 變量用于半導體數據的分段.考慮到傳感器初始化，不使用前6 個采樣時刻數據.從第7 個采樣數據開始，采用最短長度法，每個批次保留85 個采樣時刻，然后進行分段.首先根據經驗選擇邊緣效應長度為5.不同采樣時刻的第一主元貢獻率差值結果如圖5(a)～(e)所示.

圖4 半導體數據變量Fig.4 The variables of the semi-conductor data

圖5 半導體數據第一主元貢獻率差值Fig.5 The differential chart of CPV(1)of semi-conductor data

由圖5 可知:半導體數據可分為3 段:第1時刻到第43 時刻對應第1 階段，第44 時刻到第61 時刻對應第2 階段，第62 時刻到第85 時刻對應第3 階段.考慮到去除前6 個采樣時刻，實際的分段時刻為49 采樣時刻和70 采樣時刻兩處，很好地反映出實際的生產階段.

4 結論

對青霉素批次數據訓練集進行主元分析后，用第一主元貢獻率差值的方法尋求批次過程分段.實驗表明:這種算法能較好地完成批次過程分段，并且與實際的青霉素生長階段一致.然后將這種算法應用到半導體實際生產數據中，分段效果同樣明顯，驗證了分段理論的有效性.批次過程分段理論能將批次故障范圍縮小，對于故障檢測具有重要意義.

［1］Nomikos Paul，Macgregor John F.Muti-way Partial Least Squares in Monitoring Batch Processes［J］.Chemometries and Intelligent Laboratory Systems，1995，30(1):97－108.

［2］Nomikos Paul，Macgregor John F.Monitoring Batch Processes Using Multi-way Principal Component Analysis［J］.AIChE Journal，1994，40(8):1361－1375.

［3］李元.MPCA 在間歇反應過程故障診斷中的應用［J］.化工自動化及儀表，2003，30(4):10－12.

［4］趙春暉，王福利，姚遠，等.基于時段的間歇過程統計建模、在線監測及質量預報［J］.自動化學報，2010，36(3):366－372.

［5］孟艷.多階段間歇過程監測中動態特性轉變的識別方法的研究［D］.北京:北京化工大學，2009:15－18.

［6］張佳.基于多階段MPCA 方法的間歇過程監測研究［D］.北京:北京化工大學，2009:43－48.

［7］Sun Wei，Meng Yan，Palazoglu Ahmet，et al.A Method for Multiphase Batch Process Monitoring Based on Auto Phase Identification［J］.Journal of Process Control，2011，21(4):627－638.

［8］周東華，李鋼，李元.數據驅動的工業過程故障診斷技術［M］.北京:科學出版社，2011:23－29.

［9］Wise B M，Gallagher N B，Butler S W，et al.A Comparison of Principal Component Analysis，Multiway Prin-cipal Component Analysis，Trilinear Decomposition and Parallel Factor Analysis for Fault Detection in a Semiconductor Etch Process［J］.Chemo-motrics，1999，13:379－396.