基于DTW算法的參數缺失時的核動力系統故障診斷技術

趙 鑫，蔡 琦，王曉龍

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

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

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

1 DTW算法的基本原理

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

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

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

(1)

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

aij=‖xi-yj‖w

(2)

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

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

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

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

(3)

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

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

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

(4)

圖2 pd到pd+1的遞推關系Fig.2 Recursive relation between pd and pd+1

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

p=pmin/k

(5)

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

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

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

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

2.1 目標時間序列的處理

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

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

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

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

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

(6)

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

(7)

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

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

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

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

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

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

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

圖5 標準多元時間序列參數的構建Fig.5 Construction of standard multivariate time series parameter

2.3 滑動時間窗口的動態故障診斷實時模型

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

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

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

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

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

(8)

如果序列中存在監測參數缺失，則用符號χ進行填充。

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

(9)

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

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

ε=B-β

(10)

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

(11)

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

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

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

整個計算的模型結構如圖6所示。

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

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

對于待測時間序列出現監測參數缺失的情況，使原本兩條等長時間序列轉變為不等長時間序列累積距離的計算，如果標準事故序列在采集過程中出現缺失，該算法同樣可通過計算累積距離得到計算結果。

3 實例驗證

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

3.1 故障診斷模型的構建

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

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

3.2 故障診斷結果及分析

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

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

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

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

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

圖7 完整數據下的故障診斷結果Fig.7 Diagnostic result with complete data

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

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

4 結論

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