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

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

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

0 引 言

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

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

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

1 理論基礎(chǔ)

1.1 Vold-Kalman濾波器

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

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

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

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

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

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

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

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

其中M為N×N的矩陣。

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

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

其中C表示載波矩陣。

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

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

基于上式，給出目標(biāo)分量的幅值包絡(luò)矩陣A：

其中E表示單位矩陣。

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

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

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

2.1 基于EDM的時(shí)頻脊線提取方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 試驗(yàn)分析

3.1 模擬信號仿真分析

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

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

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

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

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

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

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

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

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

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

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

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

式中：N為信號長度；下標(biāo)r表示實(shí)際值；下標(biāo)e為估計(jì)值。計(jì)算含噪各分量的各項(xiàng)指標(biāo)，如表1所示。

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

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

圖9 進(jìn)動(dòng)方向系數(shù)圖Fig.9 Shape and directivity index

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

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

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

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

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

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

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

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

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

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

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

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

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

圖15 進(jìn)動(dòng)方向系數(shù)Fig.15 Shape and directivity index

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

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

4 結(jié) 論

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