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

付 波，李樂(lè)全，李超順，王詩(shī)雯

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

0 引 言

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

上述方法是對(duì)單一信號(hào)進(jìn)行分析，然而同源多傳感器信號(hào)的聯(lián)合信息可以更加全面地提供相關(guān)旋轉(zhuǎn)機(jī)械健康狀況重要信息［7-9］，吳峰崎［10］等通過(guò)將全譜方法擴(kuò)展并結(jié)合級(jí)聯(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)矢量信號(hào)短時(shí)能量時(shí)頻分布的規(guī)律。綜上所述，時(shí)頻分析方法結(jié)合以全頻譜和全息譜為代表的信息融合技術(shù)可以實(shí)現(xiàn)轉(zhuǎn)子系統(tǒng)非平穩(wěn)振動(dòng)信號(hào)分析，但大多數(shù)方法都是以短時(shí)窗穩(wěn)定為假設(shè)，難以實(shí)現(xiàn)對(duì)機(jī)組非平穩(wěn)過(guò)程的準(zhǔn)確描述。

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

1 理論基礎(chǔ)

1.1 Vold-Kalman濾波器

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

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

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

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

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

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

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

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

其中M為N×N的矩陣。

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

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

其中C表示載波矩陣。

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

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

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

其中E表示單位矩陣。

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

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

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

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

轉(zhuǎn)子旋轉(zhuǎn)速度作為VKF 重要的輸入?yún)⒘?，其?zhǔn)確性直接影響VKF 的提取能力，所以在轉(zhuǎn)子多分量非平穩(wěn)信號(hào)中精確提取瞬時(shí)頻率非常重要。由此，本節(jié)提出基于EDM的時(shí)頻脊線提取方法，通過(guò)引入歐幾里得距離重構(gòu)時(shí)頻矩陣，提高時(shí)頻聚集性，解決了傳統(tǒng)方法無(wú)法高精度提取多分量復(fù)雜信號(hào)時(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)證算法有效性（為簡(jiǎn)化計(jì)算將時(shí)頻矩陣值都設(shè)置為負(fù)值），該矩陣列代表時(shí)間步長(zhǎng)，行表示頻率，設(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）對(duì)最后一列應(yīng)用懲罰產(chǎn)生最后一列的最小值tip11是8，將最后一列中的最小值與當(dāng)前索引值相加value′11為10，進(jìn)而更新第1列中其余元素的值，如下所示：

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

（5）從矩陣的最后一列開(kāi)始，找出最小值，跟蹤最小值坐標(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ù)雜信號(hào)時(shí)頻表示方法

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

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

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

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

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

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

3 試驗(yàn)分析

3.1 模擬信號(hào)仿真分析

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

式中：x表示為轉(zhuǎn)子橫軸振動(dòng)信號(hào)；y表示為轉(zhuǎn)子縱軸振動(dòng)信號(hào)；信號(hào)x基頻分量阻尼系數(shù)χx被設(shè)置為0.6；信號(hào)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è)合成信號(hào)，采樣頻率為2 048 Hz，采樣時(shí)間為1 s，在信號(hào)中人工加入高斯噪聲，信號(hào)的信噪比為2.78 dB。噪聲信號(hào)如圖1 所示。由快速傅里葉變換獲得的x和y的頻域波形如圖2（c）和（d）所示，可以看出噪聲很強(qiáng)，人們無(wú)法從頻譜中找到有用的信息。

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

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

圖3 信號(hào)x的EMD分解圖Fig.3 EMD diagram of signal x

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

圖4 降噪后信號(hào)分析圖Fig.4 Signal analysis after noise reduction

在信號(hào)時(shí)頻圖譜中通過(guò)基于EDM 的時(shí)頻脊線提取方法提取時(shí)頻脊線并獲取瞬時(shí)頻率估計(jì)值，如圖5所示（藍(lán)色線段為擬合值，紅色實(shí)線為信號(hào)真實(shí)值），將其作為轉(zhuǎn)速信號(hào)應(yīng)用于VKF提取信號(hào)諧波分量，提取結(jié)果如圖6、圖7 所示。通過(guò)希爾伯特變換求取幅值包絡(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)諧波分量對(duì)比如圖8所示，（左側(cè)為信號(hào)x的重構(gòu)諧波信號(hào)，右側(cè)為信號(hào)y的重構(gòu)諧波信號(hào)）。為驗(yàn)證方法的性能，通過(guò)定義特征平均絕對(duì)百分比誤差（Mean Absolute Percentage Error，MAPE）和均方根誤差（Root Mean Square Error，RMSE）指標(biāo)進(jìn)行評(píng)價(jià)。

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

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

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

轉(zhuǎn)子振動(dòng)信號(hào)的進(jìn)動(dòng)方向系數(shù)計(jì)算公式為SDI= sin(φxit-φyit)，其中φxit、φyit表示x、y諧波信號(hào)相位，估計(jì)結(jié)果如圖9所示，左側(cè)為含噪信號(hào)方向系數(shù)圖，右側(cè)為不含噪信號(hào)方向系數(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)表示所提出的方法對(duì)瞬時(shí)軌跡特征估計(jì)值非常接近真實(shí)值。

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

圖10 信號(hào)時(shí)頻表示對(duì)比圖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)分量，從信號(hào)的方向指向系數(shù)計(jì)算結(jié)果（如圖9所示）可以證實(shí)高分辨率時(shí)頻圖所展示進(jìn)動(dòng)分量的正確性。

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

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

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

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

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

圖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)一步通過(guò)計(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

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

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

4 結(jié) 論

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