基于EDM-VKF 的轉子時變復雜信號時頻表示方法

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

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

0 引 言

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

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

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

1 理論基礎

1.1 Vold-Kalman濾波器

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

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

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

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

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

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

通過設置差階s=2，上式變為：

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

其中M為N×N的矩陣。

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

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

其中C表示載波矩陣。

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

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

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

其中E表示單位矩陣。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

（1）通過基于經驗模態分解的多尺度小波閾值降噪算法處理原始信號，得到降噪后的轉子振動信號。

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

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

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

3 試驗分析

3.1 模擬信號仿真分析

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

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

x和y具有相同的瞬時頻率分量。產生兩個合成信號，采樣頻率為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

將原始信號通過經驗模態分解算法將振動信號分解為本征模態分量，以x軸信號為例分解結果如圖3所示。

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

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

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

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

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

真實諧波分量與重構諧波分量對比如圖8所示，（左側為信號x的重構諧波信號，右側為信號y的重構諧波信號）。為驗證方法的性能，通過定義特征平均絕對百分比誤差（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 所示，清晰觀察出擬合分量與真實分量擬合度較高，表明本文所提出的方法對重構諧波信號估計精度較高，噪聲魯棒性高。

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

圖9 進動方向系數圖Fig.9 Shape and directivity index

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

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

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

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

3.2 水輪機轉子不對中信號分析

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

圖11 轉子不對中信號分析圖Fig.11 Rotor misalignment signal analysis diagram

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

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

圖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

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

圖15 進動方向系數Fig.15 Shape and directivity index

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

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

4 結 論

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