趙 鑫，蔡 琦，王曉龍

(海軍工程大學(xué) 核科學(xué)技術(shù)學(xué)院，湖北 武漢 430033)

隨著以太網(wǎng)在核動力系統(tǒng)上的應(yīng)用，核動力系統(tǒng)運行期間的各類監(jiān)測參數(shù)可通過相關(guān)網(wǎng)絡(luò)實時獲取，因此數(shù)據(jù)驅(qū)動的核動力系統(tǒng)的故障診斷成為今后的發(fā)展趨勢，而由于設(shè)備之間的干擾和網(wǎng)絡(luò)傳輸?shù)腻e誤導(dǎo)致部分監(jiān)測參數(shù)出現(xiàn)隨機缺失，目前對監(jiān)測數(shù)據(jù)完備情況下的故障診斷相對成熟，其對應(yīng)的算法有支持向量機、時間序列、神經(jīng)網(wǎng)絡(luò)等[1-3]，對監(jiān)測參數(shù)出現(xiàn)隨機缺失的實時故障診斷研究較少：文獻[4]采用了集成學(xué)習(xí)的概念，通過對基分類器的互補優(yōu)勢提升參數(shù)缺失情況下的診斷精度；文獻[5]使用樸素貝葉斯和EM算法從缺失數(shù)據(jù)容忍和修復(fù)兩個角度對數(shù)據(jù)進行預(yù)測，提高了歷史監(jiān)測數(shù)據(jù)隨機缺失情況下的目標函數(shù)預(yù)測；文獻[6]采用了多元動態(tài)時間彎曲(DTW, dynamic time warping)算法對多維時間序列進行模式匹配，對完整的事故過程具有較好的匹配效果。

上述研究為缺失數(shù)據(jù)下的故障診斷與模式識別，主要分為兩大類，一類是對缺失數(shù)據(jù)采取了修補的措施，其又可分為均值填補、回歸填補等固定填補方式以及構(gòu)建支持向量機、EM算法等動態(tài)填補模型的方式[7]；另一類是缺失參數(shù)采取策略，通過構(gòu)建多個診斷模型進行診斷結(jié)果的互補或采用DTW算法對數(shù)據(jù)直接進行相似度比較[8]。由于核動力系統(tǒng)的監(jiān)測參數(shù)具有高維度和非線性的特點，傳統(tǒng)方式填補后的診斷效果不理想，動態(tài)填補模型如貝葉斯模型的構(gòu)建需要一定量的數(shù)據(jù)作為支撐，而核動力系統(tǒng)故障工況的小樣本特點導(dǎo)致了該方案可行性較差。為此，本文基于多元DTW[9-10]算法，提出滑動時間窗口的距離匹配模型，尋找目標模式之間的最小距離，通過選取合適的窗口大小和監(jiān)測參數(shù)減小DTW匹配模型的復(fù)雜度，提高匹配模型的計算速度，用于核動力系統(tǒng)監(jiān)測參數(shù)隨機缺失下的實時故障診斷。

1 DTW算法的基本原理

DTW算法優(yōu)化了特征參數(shù)錯位所產(chǎn)生的影響，基本原理是尋找兩個時間序列之間的最優(yōu)彎曲路徑，序列中的數(shù)據(jù)點根據(jù)坐標值去匹配另一條序列中最具相同特征的點，數(shù)據(jù)點的距離和即為最優(yōu)彎曲距離的累加和[11](圖1)。

圖1 DTW算法的基本原理Fig.1 Fundamental principle of DTW algorithm

假設(shè)有兩個長度不同的時間序列X=(x1,x2,…,xm)、Y=(y1,y2,…,yn)，其中，m、n分別為兩條時間序列的長度。根據(jù)兩條時間序列的坐標構(gòu)造成1個m×n的距離矩陣Am×n：

(1)

在距離矩陣Am×n中，元素aij是通過xi和yj坐標距離的計算得到的，其計算過程為：

aij=‖xi-yj‖w

(2)

當w=2時為2-范數(shù)，即歐式距離。

而兩條時間序列的DTW距離是通過Am×n尋找1條距離最小的彎曲路徑pmin：

pmin={p1,p2,…,pd,…,pk},

k∈{max(m,n),m+n+1}

(3)

其中，設(shè)pd為搜索至點aij時，彎曲路徑的當前累積距離。

該路徑的尋找過程必須滿足以下3個約束條件。1) 路徑起始點固定：搜索路徑的起點為a11，終點為amn。2) 路徑的單調(diào)性：設(shè)搜索的當前點為aij，當前累積距離為pd，pd+1=pd+ai′j′，則i′≥i，j′≥j。3) 路徑的連續(xù)性：設(shè)搜索的當前點為aij，當前累積距離為pd，pd+1=pd+ai′j′,則i′≤i+1，j′≤j+1。根據(jù)上述3個條件，第1條確定搜索路徑的起始點，第2、3條確定搜索路徑的下一個點是當前點的上方、右方或右上方中的1個(圖2)，其中當前點為pd，并假設(shè)此時正在搜索的當前點為aij，到下一個點pd+1的遞推關(guān)系為：

pd+1=pd+min[a(i+1)j,a(i+1)(j+1),ai(j+1)]

(4)

圖2 pd到pd+1的遞推關(guān)系Fig.2 Recursive relation between pd and pd+1

經(jīng)計算，得到最終的計算結(jié)果pmin，為解決序列長度不同導(dǎo)致的累積距離存在差異性，需通過式(5)對累積距離進行平均化處理：

p=pmin/k

(5)

其中，p為平均后兩條序列的累積距離。

在整個搜索過程中，根據(jù)3條約束條件，DTW算法遍歷了兩條時間序列中所有觀測點，且兩條序列中的每個點在另一條序列中找到相對應(yīng)的點，但這也反映了該算法的計算量較大。

2 滑動窗口下的多元DTW算法的模式識別模型

當通過兩個一維不等長時間序列的DTW距離去度量兩個時間序列的相似度時，由于一維時間序列在監(jiān)測過程中受到波動、噪聲等因素的影響，對DTW的準確度會產(chǎn)生一定的影響，加之目前核動力系統(tǒng)監(jiān)測數(shù)據(jù)具有高維度的特點，多元DTW算法的模式識別模型具有較強的魯棒性和更高的準確率。

2.1 目標時間序列的處理

事故引入后，改變了系統(tǒng)穩(wěn)定運行的狀態(tài)，短期內(nèi)系統(tǒng)的監(jiān)測參數(shù)會產(chǎn)生較大的波動，由于波動的大小具有一定隨機性，這會導(dǎo)致1條序列的頂點對應(yīng)另1條序列的多個點，造成兩條序列匹配點的畸形(圖3)。為優(yōu)化這種畸形匹配，需對時間序列進行平滑處理，將監(jiān)測值與該狀態(tài)之前的若干記錄值取平均處理，將波動平攤至附近的若干點上，缺點是提高了計算的復(fù)雜度。

圖3 平滑處理對DTW算法的影響Fig.3 Effect of smoothing on DTW algorithm

圖4 不同初始運行條件對算法的影響Fig.4 Effect of different initial running conditions on algorithm

系統(tǒng)的監(jiān)測參數(shù)與系統(tǒng)的運行狀態(tài)有關(guān)，如圖4所示，圖中，a、b分別為系統(tǒng)正常運行時不同運行功率系統(tǒng)監(jiān)測參數(shù)的運行變化，c、d分別為系統(tǒng)在低、高工況下發(fā)生異常時監(jiān)測參數(shù)的變化。當系統(tǒng)在高工況正常運行狀態(tài)a引入事故工況，而事故狀態(tài)下監(jiān)測參數(shù)如果具有曲線d向下波動的趨勢，經(jīng)過計算會發(fā)現(xiàn)時間序列c、d之間的累積距離大于d、b之間的累積距離，會造成d與b更接近的診斷結(jié)果，為解決該類型的誤診斷問題，對標準時間序列進行歸一化，計算公式為：

x標準=(x-xmin)/(xmax-xmin)

(6)

其中，在標準多元序列中任意一維時間序列X中，xmax為該維度中的最大值，xmin為該維度中的最小值，樣本點x經(jīng)式(6)計算后得到歸一化值x標準。而對于待測序列則采用對齊穩(wěn)態(tài)工況下監(jiān)測值的歸一化方法，計算公式為：

(7)

待測樣本點x′使用式(7)得到x待測。其中，x標穩(wěn)為標準序列穩(wěn)定時的監(jiān)測值；x待穩(wěn)為待測序列穩(wěn)定時的監(jiān)測值。

2.2 多元DTW累積距離計算模型

目前針對多元時間序列的DTW算法主要分為兩類。一類為將多維時間序列轉(zhuǎn)化成一維時間序列，文獻[12]則采用距離度量的策略進行多維時間序列的降維轉(zhuǎn)化；將多元時間序列應(yīng)用多維分段擬合，對每個擬合段求取總誤差作為DTW的1個元素，再應(yīng)用DTW算法進行累積距離的計算[13]。另一類為對每一維時間序列應(yīng)用差分、歸一化處理后，采用DTW算法計算時間序列每一維對應(yīng)序列的累積距離，按一定的方式擬合計算累積距離得到最終計算結(jié)果[14]。

由于核動力系統(tǒng)的監(jiān)測數(shù)據(jù)是隨機缺失的，同一時刻的某一監(jiān)測數(shù)據(jù)可能存在缺失情況，將多維度的監(jiān)測數(shù)據(jù)轉(zhuǎn)化為一維時間序列會產(chǎn)生對應(yīng)監(jiān)測參數(shù)錯位的情況，加之計算每一維對應(yīng)時間序列的累積誤差更能說明兩條時間序列的相似性，故本文采用第2類的多元DTW算法進行相似度擬合。

考慮到每種事故工況的監(jiān)測參數(shù)有相對固定的變化趨勢，設(shè)定事故發(fā)生后有明顯變化的監(jiān)測參數(shù)為該事故下的敏感參數(shù)。本文采用多元DTW算法進行累積距離計算時，選用事故發(fā)生后波動較大的幾類監(jiān)測參數(shù)構(gòu)成多元時間序列。為此本文多元時間序列采用了每種事故工況參數(shù)變化最明顯的前幾類監(jiān)測數(shù)據(jù)組合在一起構(gòu)成標準多元事故序列，構(gòu)建形式如圖5所示。

最終篩選出由上述監(jiān)測參數(shù)及差分數(shù)據(jù)組成的多元時間序列，計算兩個多元時間序列中對應(yīng)的每一維的DTW累積距離，每個一維時間序列的相似性越高，DTW的累積距離越小，對每一維時間序列累積距離求和即得到多元時間序列的累積距離。

在每個事故工況下，將上述敏感參數(shù)按照時間軸構(gòu)成多元時間序列，作為模式識別的標準序列。將未知工況的多元時間序列保留上述監(jiān)測參數(shù)，并計算差分數(shù)據(jù)，整理生成待測序列。通過待測多元時間序列依次與各類模式識別的標準序列計算得到DTW累積距離。由于累積距離越小則兩條時間序列越相似，待測序列與累積距離最小的多元時間序列模式一致。

圖5 標準多元時間序列參數(shù)的構(gòu)建Fig.5 Construction of standard multivariate time series parameter

2.3 滑動時間窗口的動態(tài)故障診斷實時模型

本文假設(shè)原本完整的多元事故序列由于傳輸過程不穩(wěn)定性造成了數(shù)據(jù)的隨機缺失，導(dǎo)致原本等長的時間序列出現(xiàn)部分監(jiān)測數(shù)值的缺失，通過DTW算法對不等長的序列進行模式識別。多元時間序列的模式識別模型對于進行實時的故障診斷存在兩個問題：1) 實時的多元時間序列不能準確定位事故的起點，系統(tǒng)運行時監(jiān)測數(shù)據(jù)是實時生成的，事故的引入時間點也是未知的，而事故的發(fā)生起點決定了待測多元序列的起點；2) 隨運行時間的增長，多元待測序列的長度越來越長，隨之而來的DTW算法的計算量不斷增大，實時性診斷難以實現(xiàn)。如何定位事故的起點、減少計算的時間成為診斷模型的關(guān)鍵。

無線通信網(wǎng)絡(luò)中，移動站和移動站之間、移動站和網(wǎng)絡(luò)控制中心之間的信息傳遞都是通過無線信道來進行的。但是無線信道的開放性使得用戶身份信息完全曝光在信道中，攻擊者可以通過信道獲得合法用戶的身份信息，然后假冒合法用戶身份進入網(wǎng)絡(luò)，并假冒合法身份進行網(wǎng)絡(luò)資源訪問、使用通信服務(wù)，或者假冒網(wǎng)絡(luò)端基站欺騙其他移動用戶。

為解決DTW算法計算量過大的情況，本文通過構(gòu)造滑動窗口DTW模型來簡化運算量，具體步驟如下。

1) 實時獲取核動力系統(tǒng)監(jiān)測參數(shù)的運行數(shù)據(jù)，按照圖5確定的敏感參數(shù)對獲取的監(jiān)測數(shù)據(jù)進行篩選，將t時刻的參數(shù)記為αt=(αt1,αt2,…,αtp)，其中p為監(jiān)測敏感參數(shù)的數(shù)量。

2) 設(shè)定多元時間序列的長度為β，則記錄t、t-1時刻，直至t-(β-1)時刻，并組成多元待測時間序列αtβ(其中p為多元時間序列的維度)：

(8)

如果序列中存在監(jiān)測參數(shù)缺失，則用符號χ進行填充。

3) 對每個監(jiān)測參數(shù)進行δ范圍內(nèi)平滑處理。對αtβ矩陣中每個元素的值由其本身及其前δ-1個時刻的原始監(jiān)測數(shù)據(jù)取平均值得到，如果前δ-1個時刻監(jiān)測值存在缺失，則將缺失值χ取0，缺失數(shù)據(jù)的個數(shù)記為η，計算公式為：

(9)

4) 將已知事故工況下的多元時間序列按照步驟1～3得到標準DTW多元事故序列，設(shè)該標準時間長度為B，與待測時間序列不同，B?β。

5) 對于1條已知的事故標準序列，從事故起點依次滑動選取長度為β多元序列，與待測序列αtβ應(yīng)用多元DTW算法計算累積距離，如果出現(xiàn)監(jiān)測數(shù)據(jù)缺失的符號χ，則直接跳過該符號，比較兩條不等長的時間序列，該事故工況下計算的次數(shù)ε為：

ε=B-β

(10)

φt=(φ1,φ2,…,φε)

(11)

其中：待測序列φε為第ε次的累積距離；φt中最小值φtmin即為待測多元序列在t時刻的最小累積距離。

6) 應(yīng)用步驟5依次計算待測時間序列與其他事故標準多元時間序列的最小累積距離。

7) 比較各事故得到的最小累積距離，待測時間序列應(yīng)與產(chǎn)生最小累積距離值最小的事故標準事故序列中的某一段最為相似，與對應(yīng)的標準事故模式相一致。

整個計算的模型結(jié)構(gòu)如圖6所示。

由于待測多元時間序列的長度β不需過大，事故發(fā)生后監(jiān)測參數(shù)的變化會在待測序列中直接反映，通過滑動窗口進行距離匹配，在標準事故序列中不需要定位事故的起點。β是固定的，不會隨著時間的推移增長，而不同的事故序列DTW的計算過程并不沖突，故步驟5、6計算采用并行計算的形式，對于單個事故序列長度B較長的情況，可對事故序列選取多個斷點，同時計算累積距離，以滿足實時性的要求。

圖6 滑動時間窗口的動態(tài)實時故障診斷模型Fig.6 Dynamic real-time fault diagnosis model of sliding time window

對于待測時間序列出現(xiàn)監(jiān)測參數(shù)缺失的情況，使原本兩條等長時間序列轉(zhuǎn)變?yōu)椴坏乳L時間序列累積距離的計算，如果標準事故序列在采集過程中出現(xiàn)缺失，該算法同樣可通過計算累積距離得到計算結(jié)果。

3 實例驗證

本文為檢驗滑動時間窗口動態(tài)故障診斷實時模型辨識的準確性，選用了核動力系統(tǒng)一回路失水事故(右環(huán)路U型管破口百分比53%)、一回路失流事故(單環(huán)路兩個主閘閥開度百分比50%)、二回路主蒸汽管道破口事故(右主蒸汽管道破口百分比11%)與正常工況(與上述事故工況發(fā)生前的運行狀態(tài)相一致)進行實時在線監(jiān)測，事故工況選取了中等嚴重程度。

3.1 故障診斷模型的構(gòu)建

經(jīng)過多次計算擬合，根據(jù)3類事故工況發(fā)生的基本原理，模型最終選用的敏感參數(shù)為：堆平均溫度、穩(wěn)壓器壓力、穩(wěn)壓器水位、穩(wěn)壓器溫度、左右環(huán)路冷卻劑流量、左右蒸汽發(fā)生器蒸汽流量、左右蒸汽發(fā)生器壓力、左右蒸汽發(fā)生器水位，構(gòu)建標準事故序列[15-16]。

與模型相關(guān)的其他參數(shù)值為：標準失水事故序列的長度為300，標準序列的失流事故長度為300，標準序列的主蒸汽管道破口長度為250，標準序列的正常工況長度為100，待測序列滑動窗口的長度為10，平滑處理參數(shù)δ為5。

3.2 故障診斷結(jié)果及分析

對模型的測試過程是在系統(tǒng)正常運行時投入事故工況，得到相關(guān)的監(jiān)測結(jié)果。本文主要從兩個角度對模型進行測試，一種為標準事故序列數(shù)據(jù)完整而待測序列數(shù)據(jù)出現(xiàn)缺失，另一種為標準事故序列數(shù)據(jù)出現(xiàn)缺失而待測序列數(shù)據(jù)同樣出現(xiàn)缺失。通過實例進行測試，并對照了監(jiān)測參數(shù)無缺失的情況。

待測多元時間序列選擇系統(tǒng)正常運行時突發(fā)失流事故(單環(huán)路兩個主閘閥開度為30%)，由于篇幅限制，部分監(jiān)測參數(shù)的變化趨勢如圖7所示，并給出了監(jiān)測參數(shù)未缺失情況下的診斷結(jié)果。

從故障診斷結(jié)果的角度進行分析，事故在25 s時引入，區(qū)域A為事故發(fā)生前，由于監(jiān)測參數(shù)在事故發(fā)生前處于穩(wěn)定狀態(tài)，故無論是累積距離(圖7b)還是最小累積距離出現(xiàn)的位置(圖7c)的曲線都比較穩(wěn)定，診斷結(jié)果為正常運行工況；在區(qū)域B，由于事故的引入，使系統(tǒng)運行狀態(tài)從正常轉(zhuǎn)為失流事故，故表征正常工況的最小累積距離不斷增大，失流事故的最小累積距離不斷減小，加上其他故障曲線的干擾，在25～35 s的診斷結(jié)果出現(xiàn)錯誤波動，35 s以后的診斷基本正確；在區(qū)域C中，由于事故的監(jiān)測參數(shù)波動更符合失流事故的波動形式，故失流事故最小累積距離與其他事故的相比最小，故診斷結(jié)果為失流事故，而最小累積距離出現(xiàn)的位置也符合失流序列的發(fā)展過程，故圖7c的失流事故曲線大致呈現(xiàn)直線增長的趨勢，其斜率應(yīng)與事故的嚴重程度有關(guān)，由于標準事故為50%開度，測試事故為30%開度，故斜率小于1；在區(qū)域D，由于超出了標準序列的測試長度，故在該區(qū)域的后期出現(xiàn)預(yù)測偏差。

在監(jiān)測參數(shù)缺失率為10%的情況下，模型計算的累積距離及最小累積距離出現(xiàn)的位置曲線如圖8a、b所示。在監(jiān)測參數(shù)缺失率為10%、標準序列缺失率為10%的情況下，模型計算的累積距離和最小累積距離出現(xiàn)的位置如圖8c、d所示，最終故障診斷結(jié)果如圖8e、f所示。

從故障診斷結(jié)果的角度進行分析，事故在25 s時引入，各曲線的最小累積距離出現(xiàn)的位置的變化趨勢大致與圖7c相同，但隨缺失數(shù)據(jù)的增加，診斷曲線的波動不斷增大，原因在于DTW算法在比較不等長的時間序列時，會通過相鄰點尋找最小的累積距離，此時的計算距離一般較參數(shù)完整情況下的大，進而導(dǎo)致了曲線的波動，但從計算的辨識結(jié)果看，曲線的波動對累積距離大小之間的比較影響較小，辨識結(jié)果并未產(chǎn)生錯誤。正常工況的最小累積距離曲線在事故發(fā)生后一直呈現(xiàn)不斷增大的情況，這反映了事故工況下系統(tǒng)偏離正常狀態(tài)逐漸增大。

圖7 完整數(shù)據(jù)下的故障診斷結(jié)果Fig.7 Diagnostic result with complete data

a、b、e——標準序列完整，待測序列缺失率10%；c、d、f——標準序列缺失率10%，待測序列缺失率10%圖8 缺失數(shù)據(jù)下的故障診斷結(jié)果Fig.8 Diagnostic result with missing data

經(jīng)過上述測試，該模型的診斷準確率較高，并且可對故障工況下的事故序列進行復(fù)盤和溯源分析。在計算量允許的范圍內(nèi)可添加其他事故的標準序列來拓展模型的辨識工況。

4 結(jié)論

通過計算和分析，證明本文提出的滑動時間窗口的動態(tài)故障診斷實時模型能以較高的準確率識別出事故工況的發(fā)生，對存在缺失的監(jiān)測數(shù)據(jù)有較強的冗余能力，同時根據(jù)模式識別結(jié)果中最小累積距離出現(xiàn)的位置可進行多元時間序列的溯源，在標準事故序列中尋找到與當前系統(tǒng)最為接近的狀態(tài)。該方法較好地平衡了DTW算法計算時間過長，發(fā)揮了算法在比較不等長多元時間序列的優(yōu)勢。對于模型中未考慮到的故障類型，可通過添加事故標準序列進行識別故障類型的拓展，為核動力系統(tǒng)在線運行支持提供了新的手段。