基于EDM-VKF 的轉(zhuǎn)子時變復雜信號時頻表示方法

付 波，李樂全，李超順，王詩雯

（1. 湖北工業(yè)大學 太陽能高效利用及儲能運行控制湖北省重點實驗室，湖北 武漢 430068；2. 華中科技大學土木與水利工程學院，湖北 武漢 430074）

0 引 言

由于旋轉(zhuǎn)機械運行條件的多變性，其信號中普遍存在非平穩(wěn)特性［1，2］，而時頻分析方法是表征非平穩(wěn)信號時變特征的有效工具［3］，經(jīng)典時頻分析方法［4］包括短時傅里葉變換、小波變換和Wigner-Ville 分布。高傳昌［5］等人基于經(jīng)典時頻分析方法分析了結(jié)構(gòu)參數(shù)對壓力脈動的影響，然而受海森堡測不準原理和交叉項的限制，經(jīng)典時頻分析方法的時頻分辨率較低，難以實現(xiàn)非平穩(wěn)信號的精確描述。近年來涌現(xiàn)出一系列后處理方法［6］來提高這些經(jīng)典方法的時頻分辨率，如重分配方法、同步壓縮變換等。

上述方法是對單一信號進行分析，然而同源多傳感器信號的聯(lián)合信息可以更加全面地提供相關(guān)旋轉(zhuǎn)機械健康狀況重要信息［7-9］，吳峰崎［10］等通過將全譜方法擴展并結(jié)合級聯(lián)全譜和支持向量機實現(xiàn)了旋轉(zhuǎn)機械多種故障的分析和診斷。溫廣瑞等［11，12］將Vold-Kalman 濾波階次跟蹤技術(shù)和分數(shù)階傅里葉變換與全息譜結(jié)合實現(xiàn)轉(zhuǎn)子啟停車故障特征的提取。王麗雅等［13］將全矢譜技術(shù)與短時傅里葉變換相結(jié)合提出了短時矢譜的概念，分析了壓縮機喘振故障時振動矢量信號短時能量時頻分布的規(guī)律。綜上所述，時頻分析方法結(jié)合以全頻譜和全息譜為代表的信息融合技術(shù)可以實現(xiàn)轉(zhuǎn)子系統(tǒng)非平穩(wěn)振動信號分析，但大多數(shù)方法都是以短時窗穩(wěn)定為假設，難以實現(xiàn)對機組非平穩(wěn)過程的準確描述。

本文提出了一種轉(zhuǎn)子時變復雜信號高分辨率時頻表示方法，通過提取基于歐幾里得距離的時頻矩陣（Euclidean distance based the time-frequency matrix，EDM）的時頻脊線獲取精準的轉(zhuǎn)子旋轉(zhuǎn)速率，其次提取轉(zhuǎn)子諧波信號的瞬時軌道特征，構(gòu)建轉(zhuǎn)子時變復雜信號高分辨率時頻表示。

1 理論基礎

1.1 Vold-Kalman濾波器

Vold-Kalman 濾波器［14-17］（Vold-Kalman Filter VKF）的基本思想是設置局部約束，聲明未知相位分配的階數(shù)是平滑的，并且這些階數(shù)的和應該接近總振動信號。旋轉(zhuǎn)機械的振動信號，包括第k次諧波分量，軸轉(zhuǎn)速的第k次分量可寫為：

式中：ak(t)表示幅度包絡；fsf(t)是軸轉(zhuǎn)速。

以ΔT為時間增量的上式的離散形式可以表示為：

在真實信號中，由于機械系統(tǒng)的慣性，與載波信號相比，信號變化緩慢。VKF流程函數(shù)被描述為：

其中εk(t)是高階項，s表示微分調(diào)和。對于離散信號，上式可以修改為：

式中：εk(n)是非齊次項（在某種意義上是小的）；?表示差分算子。

通過設置差階s=2，上式變?yōu)椋?/p>

考慮到所有樣本，（5）的矩陣形式可以重寫為：

其中M為N×N的矩陣。

對于由幾個階次的諧波分量組成的信號y（n），噪聲可以表示為：

其中δ（n）是干擾噪聲和誤差。結(jié)果，上式的矩陣形式可以表示為：

其中C表示載波矩陣。

當通過時頻分析給出或估計目標階信號xk(t)的瞬時速度時，這意味著可以獲得載波信號Θk(n)。為了求解上式以獲得幅值包絡，利用最小二乘法將非齊次項和噪聲或誤差項最小化：

式中：上標?表示共軛轉(zhuǎn)置；εT表示向量ε的轉(zhuǎn)置；r是加權(quán)因子，用于通過平衡δ和ε來控制VKF帶寬。

基于上式，給出目標分量的幅值包絡矩陣A：

其中E表示單位矩陣。

在已知幅值包絡矩陣A和載波矩陣C的情況下，目標諧波分量可以表示為：

1.2 基于全頻譜的軸心軌跡瞬時特征參數(shù)

2 基于EDM-VKF 的轉(zhuǎn)子時變復雜信號時頻表示方法

2.1 基于EDM的時頻脊線提取方法

轉(zhuǎn)子旋轉(zhuǎn)速度作為VKF 重要的輸入?yún)⒘?，其準確性直接影響VKF 的提取能力，所以在轉(zhuǎn)子多分量非平穩(wěn)信號中精確提取瞬時頻率非常重要。由此，本節(jié)提出基于EDM的時頻脊線提取方法，通過引入歐幾里得距離重構(gòu)時頻矩陣，提高時頻聚集性，解決了傳統(tǒng)方法無法高精度提取多分量復雜信號時頻脊線的弊端。

此方法以時頻矩陣第一列為起點，逐列計算更新矩陣列向量，然后根據(jù)當前列更新值，繼續(xù)計算后續(xù)列，直到計算時頻矩陣完成。然后從新矩陣每一列，找出最小值其坐標就是時頻脊線坐標，進而在時頻矩陣中提取最大能量脊。以三階正弦波矩陣為例詳細闡述基于EDM 的時頻脊線提取方法，算法步驟如下：

（1）構(gòu)建以3×3矩陣驗證算法有效性（為簡化計算將時頻矩陣值都設置為負值），該矩陣列代表時間步長，行表示頻率，設時頻譜圖矩陣TFM=-ST如下所示。

（2）以更新ST11為例時，保持第一個時間點的值不變。根據(jù)第一列中的值與根據(jù)歐式距離計算的懲罰值重構(gòu)坐標值，以達到尋優(yōu)和限制頻率跳躍的目的。

式中：D代表兩個點之間歐氏距離；tip表示當前索引最小值；value表示更新前坐標值；value′表示更新后坐標值；p為懲罰系數(shù)；本次以懲罰系數(shù)2為例計算。

（3）對最后一列應用懲罰產(chǎn)生最后一列的最小值tip11是8，將最后一列中的最小值與當前索引值相加value′11為10，進而更新第1列中其余元素的值，如下所示：

（4）使用步驟上述中相同的過程，用懲罰因子重新計算前一列的值，獲取后序列的新值。最終矩陣為：

（5）從矩陣的最后一列開始，找出最小值，跟蹤最小值坐標索引，它們構(gòu)成了組成時頻脊線的路徑fi(t)。

在本例中，重構(gòu)時頻矩陣的列最小值坐標為（1，2），（2，2），（3，2），這與所構(gòu)成矩陣的第2 行的正弦波能量路徑相匹配，驗證了方案的可行性。

2.2 基于EDM-VKF 的轉(zhuǎn)子時變復雜信號時頻表示方法

基于EDM-VKF 的轉(zhuǎn)子時變復雜信號時頻表示的構(gòu)建由以下步驟組成，流程圖如圖1所示。

圖1 基于EDM-VKF的轉(zhuǎn)子時變復雜信號時頻表示方法流程圖Fig.1 Flowchart of EDM-VKF based time-frequency representation of rotor time-varying complex signals

（1）通過基于經(jīng)驗模態(tài)分解的多尺度小波閾值降噪算法處理原始信號，得到降噪后的轉(zhuǎn)子振動信號。

（2）對于基頻分量及其諧波分量，通過基于EDM 的時頻脊線提取方法從時頻矩陣中提取時頻脊線并獲得其瞬時頻率估計值，并使用兩個基頻分量均值作為轉(zhuǎn)速信號的估計。

（3）在已知諧波分量階次和轉(zhuǎn)速信號前提下，通過VKF 將基頻和諧波階次分量從原始信號中分離。通過希爾伯特變換得到所有單分量的幅值包絡、瞬時相位和瞬時頻率。進一步計算各階諧波分量的瞬時軌道的長軸、短軸、瞬時軌道的正反向分量以及形狀指向系數(shù)。

（4）將諧波分量軸軌道的正向分量和反向分量幅值投影到時頻圖三維坐標中，得到諧波分量高分辨率時頻表示。所有諧波時頻表示的疊加是原始信號的轉(zhuǎn)子時頻表示。

3 試驗分析

3.1 模擬信號仿真分析

在仿真分析中，對仿真信號的參數(shù)進行了精確獲知和精確控制，實現(xiàn)本方法揭示的信號特征與預定義的理論信號特征進行比較，檢驗該方法的性能。在本節(jié)中，將通過數(shù)值模擬信號分析來說明和測試所提出的方法。第一組模擬信號包含兩個非平穩(wěn)分量，它們模擬停機過程中轉(zhuǎn)子的振動響應，用數(shù)學表達式為：

式中：x表示為轉(zhuǎn)子橫軸振動信號；y表示為轉(zhuǎn)子縱軸振動信號；信號x基頻分量阻尼系數(shù)χx被設置為0.6；信號y基頻阻尼系數(shù)χy被設置為0.5；f1(t) = 160 - 220t+ 110t2、f2(t) = 1.5f1(t)，分別表示兩個分量的瞬時頻率；x和y的相位差由φi表示，分別設置φ1= π/3、φ2= 4π/3。

x和y具有相同的瞬時頻率分量。產(chǎn)生兩個合成信號，采樣頻率為2 048 Hz，采樣時間為1 s，在信號中人工加入高斯噪聲，信號的信噪比為2.78 dB。噪聲信號如圖1 所示。由快速傅里葉變換獲得的x和y的頻域波形如圖2（c）和（d）所示，可以看出噪聲很強，人們無法從頻譜中找到有用的信息。

圖2 原始信號時域波形及頻域波形Fig.2 Original signal time domain waveform and frequency domain waveform

將原始信號通過經(jīng)驗模態(tài)分解算法將振動信號分解為本征模態(tài)分量，以x軸信號為例分解結(jié)果如圖3所示。

圖3 信號x的EMD分解圖Fig.3 EMD diagram of signal x

通過選取“sym6”的小波基函數(shù)和三層小波分解層數(shù)，對本征模態(tài)分量進行小波分解，采用軟閾值方法對信號分量進行處理，得到降噪后的信號，如圖4降噪后信號分析圖所示。由短時傅里葉變換獲得的x1和y1的頻譜如圖4（c）和（d）所示，可以直觀看出經(jīng)過經(jīng)驗模態(tài)分解結(jié)合小波閾值降噪，有效濾除了部分高頻噪聲，信號信噪比提高了一倍同時保留了信號特征，根據(jù)皮爾遜系數(shù)計算，得到重構(gòu)信號與無噪聲信號相關(guān)性系數(shù)大于0.91，從時頻譜圖可以看出，信號由兩個分量組成，然而，兩個分量從0.4 s 開始出現(xiàn)時頻線混疊跡象，難以直觀區(qū)分，從原始轉(zhuǎn)子軌跡圖和全譜也無法獲得更多信息。

圖4 降噪后信號分析圖Fig.4 Signal analysis after noise reduction

在信號時頻圖譜中通過基于EDM 的時頻脊線提取方法提取時頻脊線并獲取瞬時頻率估計值，如圖5所示（藍色線段為擬合值，紅色實線為信號真實值），將其作為轉(zhuǎn)速信號應用于VKF提取信號諧波分量，提取結(jié)果如圖6、圖7 所示。通過希爾伯特變換求取幅值包絡、瞬時相位和瞬時頻率，并構(gòu)建轉(zhuǎn)子諧波分量的瞬時軸心軌跡。

圖7 1.5倍諧波分量時域分析圖Fig.7 Time domain analysis diagram of 1.5X harmonic component

真實諧波分量與重構(gòu)諧波分量對比如圖8所示，（左側(cè)為信號x的重構(gòu)諧波信號，右側(cè)為信號y的重構(gòu)諧波信號）。為驗證方法的性能，通過定義特征平均絕對百分比誤差（Mean Absolute Percentage Error，MAPE）和均方根誤差（Root Mean Square Error，RMSE）指標進行評價。

圖8 真實諧波分量與擬合分量對比圖Fig.8 Comparison of the real harmonic component and the fitted component

式中：N為信號長度；下標r表示實際值；下標e為估計值。計算含噪各分量的各項指標，如表1所示。

通過圖8 所示，清晰觀察出擬合分量與真實分量擬合度較高，表明本文所提出的方法對重構(gòu)諧波信號估計精度較高，噪聲魯棒性高。

轉(zhuǎn)子振動信號的進動方向系數(shù)計算公式為SDI= sin(φxit-φyit)，其中φxit、φyit表示x、y諧波信號相位，估計結(jié)果如圖9所示，左側(cè)為含噪信號方向系數(shù)圖，右側(cè)為不含噪信號方向系數(shù)圖。計算量化指標如表1所示。

圖9 進動方向系數(shù)圖Fig.9 Shape and directivity index

進動方向系數(shù)實現(xiàn)了瞬時軌道方向性信息的定量表達，量化指標表示所提出的方法對瞬時軌跡特征估計值非常接近真實值。

最后，將各個諧波分量軸軌道的正向分量和反向分量振幅投影到時頻圖三維坐標中，得到諧波分量軸軌道時頻表示。原始信號的轉(zhuǎn)子信號高分辨率時頻表示由兩個諧波分量的高分辨率時頻表示疊加得到。為驗證方法的準確性和有效性，本文將短時傅里葉變換如圖10（a）、同步壓縮變換的時頻表示如圖10（b）、理論上的時頻表示如圖10（c）所示和本文所提出方法的結(jié)果如圖10（d）進行對比，它具有良好的時頻分辨率，C1 和C2分量的時頻圖可以清晰地分辨出來，與理論時頻表示非常接近，并且從時頻譜圖上可以清晰地觀察出各諧波分量的進動方向信息。

圖10 信號時頻表示對比圖Fig.10 Signal time and frequency representation comparison chart

由本方法提出的時頻表示結(jié)果可以觀察出，C1分量反進動分量大于正進動分量，C2 分量正進動分量大于反進動分量，從信號的方向指向系數(shù)計算結(jié)果（如圖9所示）可以證實高分辨率時頻圖所展示進動分量的正確性。

3.2 水輪機轉(zhuǎn)子不對中信號分析

為進一步驗證本文所提出方法的有效性，選取了某旋轉(zhuǎn)機械發(fā)生轉(zhuǎn)子不對中時轉(zhuǎn)子振動信號進行分析說明，該信號采集時轉(zhuǎn)子轉(zhuǎn)速為3 180 r/min（53 Hz），信號分析圖如圖11所示。

圖11 轉(zhuǎn)子不對中信號分析圖Fig.11 Rotor misalignment signal analysis diagram

從信號頻譜分析，信號主要以一倍頻和二倍頻分量為主，在x的頻譜上二倍頻分量幅值大于一倍頻。從全頻譜直接看出信號一倍頻分量是反進動，三倍頻分量是正進動除此之外無法直觀獲取更多有價值信息。

接下來，使用本文所提出的算法對上述信號進行分析，所提出的一倍頻、二倍頻和三倍頻諧波分量如圖12、圖13 和圖14所示。從圖中可以觀察到，諧波分量可以很好的分離出來，三倍頻諧波分量是橢圓形，并且信號二倍頻諧波分量的軸心軌跡符合轉(zhuǎn)子不對中特征。

圖12 一倍頻諧波分量時域分析圖Fig.12 Time domain analysis diagram of 1X harmonic component

圖13 二倍頻諧波分量時域分析圖Fig.13 Time domain analysis diagram of 2X harmonic component

圖14 三倍頻諧波分量時域分析圖Fig.14 Time domain analysis diagram of 3X harmonic component

計算諧波分量瞬時特征參數(shù)，進一步通過計算進動方向系數(shù)（如圖15 所示），獲取諧波分量的進動方向系數(shù)，C3 始終保持正進動，C1始終保持負進動。

圖15 進動方向系數(shù)Fig.15 Shape and directivity index

接下來，構(gòu)建高分辨率時頻譜圖［如圖16（d）所示］，并給出傳統(tǒng)時頻表示方法作相比，此圖清晰地展示了信號時頻譜圖，并且提供了基頻及諧波分量的進動方向及信號時頻域的變化。

圖16 信號時頻表示對比圖Fig.16 Signal time and frequency representation comparison chart

4 結(jié) 論

本文提出了引入歐幾里得距離的重構(gòu)時頻矩陣的時頻脊線提取方法，提高了時頻聚集性和時頻脊線連貫性，為VKF 提取諧波信號分量提供精準的轉(zhuǎn)子速度信號，并基于諧波分量瞬時特征參數(shù)構(gòu)建轉(zhuǎn)子軸心高分辨率時頻表示。該方法實現(xiàn)了轉(zhuǎn)子時變復雜信號的諧波分量提取，同時計算轉(zhuǎn)子諧波信號的瞬時特征參數(shù)，可用于描述非平穩(wěn)過程中轉(zhuǎn)子諧波分量瞬時振動狀態(tài)。通過仿真實驗與轉(zhuǎn)子不對中數(shù)據(jù)分析，基于EDMVKF 的轉(zhuǎn)子時變復雜信號時頻表示方法對信號諧波分量重構(gòu)精度高，同時高分辨率的時頻圖提供了轉(zhuǎn)子諧波分量的進動方向信息，清晰地展現(xiàn)轉(zhuǎn)子不對中二倍頻信號的特征信息，為旋轉(zhuǎn)機械振動信號分析提供更直觀的特征，有利于非平穩(wěn)過程中旋轉(zhuǎn)機械轉(zhuǎn)子的狀態(tài)監(jiān)測。