基于改進(jìn)收斂交叉映射的廠級振蕩源診斷

周學(xué)均，陳小強(qiáng)，謝 磊

（1.中電國際胡布發(fā)電有限公司,巴基斯坦 卡拉奇74200；2.中電華創(chuàng)電力技術(shù)研究有限公司，江蘇 蘇州215123；3.浙江大學(xué) 智能系統(tǒng)與控制研究所，浙江 杭州310027）

現(xiàn)代化工業(yè)過程規(guī)模龐大，有數(shù)百至數(shù)千個高度耦合的控制回路，為了使生產(chǎn)效率最大化，保證產(chǎn)品質(zhì)量，所有關(guān)鍵回路必須處于或接近最佳狀態(tài)下運(yùn)行[1]。振蕩是過程控制系統(tǒng)性能下降的常見原因[2]。Z震蕩可以在一個控制回路中產(chǎn)生并傳播到其他單元，形成廠級振蕩[3]。這些震蕩有可能造成能源和材料的浪費(fèi)，產(chǎn)品質(zhì)量波動，甚至威脅控制系統(tǒng)的安全性。因此，對工業(yè)過程回路進(jìn)行準(zhǔn)確的因果分析，快速查明故障根源并采取相應(yīng)的糾正措施，是保證安全生產(chǎn)和提高經(jīng)濟(jì)效益的關(guān)鍵。然而，由于工業(yè)控制系統(tǒng)中采集的數(shù)據(jù)往往被非平穩(wěn)因素、噪聲、周期項(xiàng)等干擾，對振蕩源進(jìn)行定位是非常有挑戰(zhàn)的任務(wù)，因此大量的因果關(guān)系分析方法被用于震蕩源的定位研究。

文獻(xiàn)[4]使用互相關(guān)分析法來進(jìn)行因果分析。方法雖然簡單，但不適用于非線性的情況。文獻(xiàn)[5]將主成分分析與多變量格蘭杰因果關(guān)系結(jié)合起來進(jìn)行廠級振蕩的定位。格蘭杰因果方法是一個廣泛應(yīng)用的方法，但對變量獨(dú)立性和線性具有較嚴(yán)苛的要求。除了這些時域的因果分析方法，譜格蘭杰、偏有向相干[6]、有向傳遞函數(shù)[7]等一些頻域的方法也被開發(fā)出來。然而，譜格蘭杰是基于傅里葉變換的方法，因此只能適用于線性平穩(wěn)過程。偏有向相關(guān)和有向傳遞函數(shù)能描述因果方向，但是不能定量衡量因果影響。

上述方法局限于分析線性數(shù)據(jù)，但工業(yè)過程往往是非線性的。貝葉斯網(wǎng)絡(luò)利用圖結(jié)構(gòu)描述變量之間的因果關(guān)系，然而它是一個有向無環(huán)圖，因此不適用于動態(tài)過程。而且貝葉斯網(wǎng)絡(luò)中的概率缺乏直觀的物理意義，難以被工程師接受。

文獻(xiàn)[8]提出了一種新的因果分析方法——收斂交叉映射（Convergent Cross Mapping，CCM），用于復(fù)雜系統(tǒng)中的變量關(guān)系分析。CCM主要基于Takens引理，即如果變量X是變量Y的原因，那么X的歷史信息能夠從Y中復(fù)現(xiàn)。CCM使用單純形投影量化X和Y的狀態(tài)依賴關(guān)系，估計(jì)狀態(tài)之間的Pearson相關(guān)系數(shù)用于衡量交叉映射能力，這種能力對應(yīng)著因果影響程度。CCM已經(jīng)在各種工業(yè)過程中得到廣泛應(yīng)用，但基本CCM的收斂判據(jù)存在主觀性，因此文獻(xiàn)[9]提出了一種遲延CCM。遲延CCM認(rèn)為原因變量和結(jié)果變量之間存在時間上的滯后，最強(qiáng)的因果作用應(yīng)該出現(xiàn)在遲延處。遲延CCM克服了基本CCM的主觀性，是近些年應(yīng)用最為廣泛、效果最好的關(guān)于CCM的研究成果。但是，當(dāng)遲延CCM應(yīng)用分析控制過程中的振蕩信號時，噪聲和周期項(xiàng)會損害CCM的可靠性。因此本文將經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥mpirical Mode Decomposition，EMD）[10]與去趨勢波動分析（Detrended Fluctuation Analysis，DFA）[11]相結(jié)合實(shí)現(xiàn)消噪，并且利用奇異譜分析（Singular Spectrum Analysis，SSA）[12]消除周期項(xiàng)。將去噪去周期后的信號應(yīng)用遲滯CCM進(jìn)行分析,能準(zhǔn)確確認(rèn)回路之間的因果關(guān)系并定位故障源。最后，通過案例驗(yàn)證了所提出方法優(yōu)勢。

1 CCM簡介

1.1 基本CCM

文獻(xiàn)[13]采用CCM方法分析非線性系統(tǒng)中變量之間的因果關(guān)系。CCM假定如果變量X對變量Y具有結(jié)果影響，記作X→Y，那么Y將會包含X的信息。相應(yīng)的因果關(guān)系能夠通過X和Y的重構(gòu)流形之間的關(guān)系來度量。具體來說，來自同一個動態(tài)系統(tǒng)的兩個樣本長度為N的離散時間序列X（t）、Y（t），X（t）的時間滯后嵌入向量可表示為：

式中，E為嵌入維度；τ為遲延參數(shù)（默認(rèn)為1）；X t為X（t）的陰影流形MX。當(dāng)變量X是變量Y的原因時（X→Y），由于結(jié)果變量Y是原因變量X的特征信息，因此可以從Y中恢復(fù)關(guān)于X的狀態(tài)信息，但是原因變量X沒有結(jié)果變量Y的特征信息，因此不能從X中獲取Y的狀態(tài)信息。判斷X→Y方向的因果關(guān)系時，CCM方法先通過測量Y的歷史值，再計(jì)算X狀態(tài)的可信度來判斷因果關(guān)系。在CCM方法中，通過MY估計(jì)值X的準(zhǔn)確度來反映變量X對變量Y的因果程度，可用相關(guān)系數(shù)ρ（標(biāo)準(zhǔn)皮爾森相關(guān)系數(shù)ρ(X,?(t)|MY)）表示：

式中，N*為嵌入向量的個數(shù)。ρX?(t)|MY和X對Y的影響程度正相關(guān)。

收斂性是區(qū)分簡單相關(guān)性和因果關(guān)系的一個重要指標(biāo)，也是CCM的一個重要特性。收斂是指CCM估計(jì)的精度隨著時序長度L的增大而提高。在實(shí)際應(yīng)用中，收斂受系統(tǒng)誤差、隨機(jī)噪聲和時間序列長度等因素影響。因此，基于CCM方法分析的因果關(guān)系通過相關(guān)系數(shù)ρ來證明，理論上ρ會隨著時間序列長度的增大而收斂。

1.2 遲滯CCM

傳統(tǒng)CCM方法已廣泛應(yīng)用于低中耦合度的動態(tài)系統(tǒng)的因果分析，但其無法解決“廣義同步”現(xiàn)象，即如果X對Y的影響很大，則Y的內(nèi)在動力從屬于X，導(dǎo)致“同步”現(xiàn)象的出現(xiàn)。文獻(xiàn)[13]提出了遲延CCM以解決這個問題。遲延CCM利用MY去估計(jì)Xt＋λ，其中λ是滯后時間。實(shí)際上，可以取不同的λ來計(jì)算不同遲延下的交叉映射指標(biāo)ρ。最大ρ值即對應(yīng)最優(yōu)遲延λ*。λ*≤0對應(yīng)X→Y，因?yàn)棣?≤0意味著原因變量過去的信息包含在結(jié)果變量中；如果λ*＞0，則反映通過原因變量的未來信息能預(yù)測結(jié)果信息，這顯然是不合理的。因此，λ*＞0意味著X→Y的關(guān)系不成立。這樣遲延CCM就能自動客觀地確定因果關(guān)系。這里使用一個耦合的洛倫茲系統(tǒng)來測試基本CCM和遲延CCM的表現(xiàn)。

圖1 基本CCM和遲延CCM的因果分析效果

2 改進(jìn)CCM

2.1 噪聲與周期性對CCM的不利影響

通過構(gòu)建一個帶擾動的四輸入四輸出耦合系統(tǒng)驗(yàn)證噪聲與周期性對CCM的不利影響，其開環(huán)傳遞函數(shù)矩陣G(q－1)和擾動信號傳遞函數(shù)矩陣N(q－1)可表示為：

為了形成振蕩傳播路徑，在第三回路L3中增加閥門粘滯故障，最終導(dǎo)致四個回路的全局振蕩。CCM的周期和噪聲影響結(jié)果如圖2所示。圖2中，紅線表示錯誤的因果路徑，紅色箭頭表示與實(shí)際情況不相吻合的因果關(guān)系方向，這些錯誤的因果關(guān)系方向?qū)е聼o法判斷振蕩源到底在哪一回路。由圖2可知，噪聲和周期項(xiàng)對CCM具有不利影響。為了提高因果分析的可靠性，必須對信號進(jìn)行去噪去周期的預(yù)處理。

圖2 CCM的周期和噪聲影響結(jié)果

2.2 去噪與去周期

首先，介紹一種度量非平穩(wěn)時間序列長程相關(guān)性的方法——去趨勢波動分析（Detrended Fluctuation Analysis，DFA）方法。DFA算法的具體步驟為：

輸入：X(t),t＝1,2,…,N；

（1）對時間序列X(t)求均值并累加：

式中，Xˉ為X(i)平均值；Y(k)為累加時間序列。

（2）利用標(biāo)準(zhǔn)皮爾森相關(guān)系數(shù)對Y(k)劃分長度為n的等長區(qū)間。用最小二乘估計(jì)法得到區(qū)間序列的趨勢值Yn(k)。

（3）去除區(qū)間序列的趨勢值，即將Yn(k)從Y(k)中去除。

（4）計(jì)算均方根波動函數(shù)：

式中，N為數(shù)據(jù)長度。

（5）增加n的值，重復(fù)步驟（1）—（4），得到波動函數(shù)F(n)。如果滿足F(n)∝nα，則定義α為度量因子。設(shè)置閾值為0.5，如果0＜α＜0.5，即時序X(t)有短期記憶性；α＝0.5表示X(t)為一個隨機(jī)過程；α＞0.5則表示X(t)有長期記憶性。

輸出：α。

將經(jīng)驗(yàn)?zāi)B(tài)分解與DFA算法相結(jié)合，用經(jīng)驗(yàn)?zāi)B(tài)分解得到信號的各個模態(tài)，再計(jì)算所有模態(tài)（包含殘余分量）的度量因子α，進(jìn)而判斷其是否為噪聲。

奇異譜分析方法（Singular Spectrum Analysis，SSA）是一種強(qiáng)大的時間序列分析方法，具有以下優(yōu)點(diǎn)：

（1）適用于非平穩(wěn)非線性時間序列。

（2）分解信號具有實(shí)際物理意義，如緩慢變化的趨勢項(xiàng)、振蕩分量以及噪聲等。SSA算法的具體步驟為：

輸入：Y(t),t＝1,2,…,N；

（1）嵌入。將累加時間序列Y(t)進(jìn)行時間滯后表示為：

式中，i＝1,2,…,K，K＝N－M＋1；M為窗口長度。軌跡矩陣表示為：

（2）奇異值分解。計(jì)算S＝XXT的特征值并按幅值降序排列（λ1≥λ2≥…≥λM≥0），設(shè)置d＝max(i)(λi＞0)＝rank(X)，則X的SVD結(jié)果為：

即計(jì)算逆對角線上元素的平均值。將XIk轉(zhuǎn)化為與Y(t)長度相同的時間序列，則Y(t)就可以由m個長度為N的時間序列重構(gòu)為：

SSA方法基于時間序列協(xié)方差矩陣進(jìn)行分析，該方法對時間序列的長度沒有要求，并且也適用于短時間序列。

對信號進(jìn)行SSA分解可以獲得一系列特征值以及相應(yīng)的特征向量和重構(gòu)時間序列。特征向量及其相應(yīng)的重構(gòu)時間序列具有相同的屬性。通常情況下，占比較大的特征值對應(yīng)趨勢項(xiàng)，兩個相近的特征值對應(yīng)周期項(xiàng)，遞減緩慢的特征值序列對應(yīng)噪聲項(xiàng)。去噪去周期的CCM流程如圖3所示。通過去噪去周期的步驟，改進(jìn)后的CCM的因果分析效果更加可靠準(zhǔn)確。

圖3 去噪去周期的CCM流程

3 仿真實(shí)驗(yàn)

為驗(yàn)證所提方法的有效性，將該方法應(yīng)用到四回路仿真控制系統(tǒng)上。去噪去周期后的CCM的結(jié)果如圖4所示。

圖4 去噪去周期后的CCM的結(jié)果

從圖4可以看出，沒有錯誤的因果關(guān)系方向，所提方法的因果分析結(jié)果與預(yù)設(shè)的系統(tǒng)結(jié)構(gòu)完全一致，也與已有研究結(jié)果相吻合[6]，因此所提方法有效且明顯改進(jìn)了標(biāo)準(zhǔn)CCM的分析效果，對于準(zhǔn)確分析因果關(guān)系是非常有必要的。

4 工業(yè)案例

為了驗(yàn)證所提方法在工業(yè)控制系統(tǒng)中的實(shí)用效果，采用田納西－伊斯曼（TE）仿真平臺作為實(shí)驗(yàn)對象。本文研究TE過程其中的17個主要回路。TE過程采樣間隔為0.1 h，共運(yùn)行72.0 h。系統(tǒng)穩(wěn)定運(yùn)行20.0 h后，在第16回路加入振蕩擾動，最終引發(fā)17個回路發(fā)生全局振蕩。按照文獻(xiàn)[13]的方法將因果結(jié)構(gòu)圖轉(zhuǎn)化為基于改進(jìn)CCM的因果流圖如圖5所示。從圖5中可知，第16和第17個變量所對應(yīng)的回路是潛在的振蕩源。因?yàn)榈?6和第17個變量都在同一個回路中，且數(shù)值越大表明越有可能是振蕩源，所以將第16或者第17個變量識別為震蕩源是準(zhǔn)確的。通過所提出的改進(jìn)CCM方法準(zhǔn)確地找到了TE工業(yè)過程的振蕩源，這說明該方法具有工業(yè)實(shí)際有效性。

圖5 基于改進(jìn)CCM的因果流圖

5 結(jié) 論

提出了一種基于去噪和去周期的CCM因果分析算法框架。首先，研究了噪聲和周期性的不利影響。然后，將EMD和FDA結(jié)合起來實(shí)現(xiàn)去噪。用SSA去除周期性。仿真結(jié)果表明，所提方法能夠有效地改善CCM的性能，可以減少對因果關(guān)系的誤判。最后，將所提出方法應(yīng)用于TE過程，以找出引起全局振蕩的振蕩源。應(yīng)用結(jié)果表明，該方法對過程控制系統(tǒng)的振蕩源定位有效。