振動分析在工程領(lǐng)域很常見，汽車中的 NVH，風(fēng)機的模態(tài)分析，塔架、葉片等機械件的振動頻率分析，發(fā)電機軸承振動，飛機發(fā)動機階次分析等等。經(jīng)常提到的知識點像信號處理，F(xiàn)FT，階次分析，頻譜估計，包絡(luò)譜分析，模態(tài)分析，頻率響應(yīng)函數(shù)估計等等。

傅里葉變換

要談信號處理就離不開傅里葉變換。我們先從最常用的傅里葉變換開始。傅里葉變換的本質(zhì)是從空間（希爾伯特空間）中找一組正交基向量 或者正交基函數(shù)，正交意思是內(nèi)積為 0，其中希爾伯特空間內(nèi)積定義為：對于向量 f， g， 有

容易證明

將時域信號投影到這組正交基上。傅里葉變換也就可以看作 f（t） 在 k 對應(yīng)的基 上的投影。對于離散傅里葉變換（DFT）

為了方便理解，我們先通過一個例子，利用 fft 函數(shù)（快速傅里葉變換，實現(xiàn)DFT的一種快速算法）將一段矩形信號（圖中深藍色信號）進行變換。通常我們得到幅值頻率圖像，可以看到對應(yīng)頻率的幅值。

為了更好的理解，我們把離散頻率對應(yīng)的正交基向量（圖淺綠色）也繪制出來，對應(yīng)圖中淺綠色。方塊綠則為幅值頻率值，此處我們首先只選了前 20 個頻率對應(yīng)的基向量進行信號重構(gòu)，得到圖中紅色信號。接下來，我們增加用于信號重建的基的數(shù)量，分別我們?nèi)∏?60，前 300，發(fā)現(xiàn)信號越來越接近原始信號。

圖表2

值得一提的是很多類似的程序會通過for循環(huán)來遍歷各個基向量來實現(xiàn)上面的程序，但 MATLAB 中的矩陣思想更應(yīng)該充分利用，既可以簡化代碼，又可以提升效率，例如程序中 baseSeries=cos（2*pi*Dfreq‘*t+Phi（1:length（Dfreq））’）計算得到的 baseSeries 變量是一個 n*m 矩陣，包含了所有的基向量n 是基向量的個數(shù)（k 的個數(shù)），m 是時域信號長度（t 的長度）。簡單理解了傅里葉變換，接下來我們就可以使用 fft 進行一個初步的轉(zhuǎn)動數(shù)據(jù)分析，我們先生成一段時長 2s，采樣頻率 600Hz 的數(shù)據(jù)，在這 2s 中，系統(tǒng)一直以 20Hz 的轉(zhuǎn)速運行。

通過 fft 我們可以得到信號在對應(yīng)轉(zhuǎn)速的 20Hz 頻率上幅值最大。這部分有很多重要的概念，包括時域/頻域分辨率，頻譜泄露，功率譜/能量譜等等。頻率分辨率簡言之就是區(qū)分兩個不同頻率的能力。假設(shè)信號的采樣頻率為 Fs， 在時域上，tn= n*Fs 反映的是時間軸，T=N*Fs 反映的是信號總時長，在頻域上，我們可以得到 N 個頻率（基向量），可以理解為 Fs 被 N 個值平分，因此頻率的分辨率為

我們從上面推導(dǎo)可以知道頻率分辨率由采樣時長決定，如果采樣時長太短，DFT 沒有能力區(qū)分兩個接近的頻率值。

這里也給大家一個小練習(xí)，例如可以合成一段信號，由 90Hz 和 100Hz 兩個頻率組成，我們只選取 0.05s 長的數(shù)據(jù)：對應(yīng) DFT 的頻率分辨率為 20Hz 》 （100Hz – 90Hz），然后進行 fft 變換，查看是否在頻域上能夠區(qū)分開這兩個信號成分。

接下來我們通過一個由兩個頻率組成的正弦信號來介紹頻譜泄露。詳見 https://www.mathworks.com/help/releases/R2020b/signal/ug/amplitude-estimation-and-zero-padding.html這兩種正弦波的頻率分別為 f1 = 100Hz 和 f2 = 202.5 Hz。采樣率為 1000Hz，信號長度為 1000 個采樣點（T=1s）。我們期望的結(jié)果是頻譜幅值只在 f1 = 100Hz 和 f2 = 202.5Hz 處的結(jié)果不為零，其他位置的頻譜幅值為零。但真實結(jié)果是這樣的：

圖表4

也就是在 202.5Hz 的信號幅值被周圍的頻率值平分了，這種現(xiàn)象就是頻譜泄露。從頻譜分辨率的角度解釋通過上面的結(jié)論：我們得到DFT的頻率分辨率為 1/T = 1Hz，對于 f1 = 100Hz 正好落在分辨率的分倉邊界上，而 f2 = 202.5Hz 不能被分辨率 1Hz 取整。

從時域角度解釋：我們?nèi)×舜翱跁r間 1s 長度的數(shù)據(jù)，對應(yīng) f1 的 100 個完整周期以及 f2 的 202.5 個周期，DFT 默認(rèn)信號是采樣時間窗內(nèi)周期出現(xiàn)的無限長信號，因此當(dāng)窗重復(fù)時對于f1來說正好是整數(shù)拼接，與原信號完全一致，而 f2 在上一個 1s 采樣時間窗結(jié)束時并非對應(yīng)信號自身周期結(jié)束，如下圖：

所以在上一個周期和下一個周期之間的信號過渡時有一個突變（周期邊界由黑色虛線表示）。這個突變導(dǎo)致信號中添加了除原始頻率之外的其他頻率。這就解釋了頻譜泄露現(xiàn)象。頻譜泄露會讓我們在計算頻域幅值時不夠準(zhǔn)確，因為能量幅值被周圍的頻率平分。

了解了頻譜泄露的原因，我們可以通過調(diào)整采樣時長，保證在每個時間窗口中都有原始信號的整數(shù)個周期，例如設(shè)置采樣窗口 T = 2s（保證所有信號都是整數(shù)周期截斷）。然而通常我們可能未必能對真實數(shù)據(jù)進行延長采樣，可以考慮補零 zero-padding 的方式。

你可以用補零來對 DFT 進行插值。補零可以獲得更準(zhǔn)確的幅值估計。補零并不能提高 DFT 的譜（頻率）分辨率。根據(jù)這種方法，將 DFT 要處理的信號用 0 補長到 2000 個點。有了這個長度，DFT 分倉之間的間距是 F_s=0.5Hz，在這種情況下，202.5Hz 正弦波的能量直接落在 DFT 分倉中。獲得 DFT 并繪制振幅估計值得到更準(zhǔn)確的幅值。

時頻分析

在實際的應(yīng)用中，信號經(jīng)常是非平穩(wěn)的，他們的頻域是隨時間在變化的，因此如果只用 FFT 無法反映出在什么時間里頻域上的這種變化。為了回答在什么時刻出現(xiàn)什么頻率的信號，我們需要進行時頻分析。例如，給定一段信號，有兩個 chirp 信號組成，第一個 chirp 信號在 0.1s 到 0.68s 一直存在，頻率隨時間 t 的函數(shù)是

第二個 chirp 信號在 0.1s 到 0.75s 一直存在，頻率隨時間t的函數(shù)是

我們希望看到在不同時間軸上都各自有哪些頻率出現(xiàn)。通過傅里葉變換（圖表6），我們確實可以得到信號中包含的頻域成分，但它回答不了各頻率成分在什么時候出現(xiàn)。所以首先想到的是對原始數(shù)據(jù)加窗，每個窗口內(nèi)做傅里葉變換。這就是短時傅里葉變換 STFT。STFT 提供了信號頻率以及對應(yīng)出現(xiàn)的時間信息。但是選擇窗口大小是個問題。在短時傅里葉變換時頻分析中，根據(jù)上面傅里葉變換的結(jié)論，選擇較短的窗口大小可以得到較好的時間分辨率，但頻率分辨率差。相反，選擇較大的窗口會獲得較好的頻率分辨率，但時間分辨率差。而且一旦選擇了窗口大小，它將在整個分析中保持不變。

如果可以提前知道待估計信號中我們想要的的頻率分量，則可以根據(jù)這些需求來選擇一個合適的窗口大小。兩個 chirp 信號的在初始時間點的瞬時頻率約為5Hz和15Hz。我們先選定窗口大小為 200ms（對應(yīng)的頻率分辨率為 1/0.2s = 5Hz）。瞬時頻率在信號的前段被分辨出來，但在后段就不那么好了（后段頻率變化快，時間分辨率不夠會導(dǎo)致頻帶太寬）。

我們把窗口縮短為 50ms（對應(yīng)的頻率分辨率 1/0.05s = 20Hz），雖然在后面高頻率差中可以分辨（時間分辨率更好），但在初始時刻低頻率差中分辨不開（頻率分辨率差）。

對于這種非平穩(wěn)雙曲線 chirp 信號，STFT 很難找到合適的時間窗口大小使得對整個頻域的各頻率都能分辨的很好。連續(xù)小波變換（CWT）可以解決 STFT 固有的分辨率問題。連續(xù)小波變換得時間窗口是可變的。如下圖：

圖表10如果要分析的信號主要是緩慢震蕩的低頻信號，而高頻信號只是在一些瞬態(tài)時長或突變過程出現(xiàn)，這種情況可以考慮連續(xù)小波變換。如果高頻信號持續(xù)時間很長并且是信號中的主要成分，這種情況使用連續(xù)小波變換就不合適了。

繪制 CWT 的時頻圖。時頻圖顏色對應(yīng)幅值，頻率軸用對數(shù)值，因為 CWT 中的頻率是對數(shù)的。從圖中可以清楚地看出信號中存在兩個雙曲 chirp 信號。使用連續(xù)小波變換，可以準(zhǔn)確地估計整個信號周期在某個瞬時出現(xiàn)的頻率成分，而不必手動設(shè)置窗口長度。

有人可能會問白色虛線是什么？白色虛線是 COI（cone of influence （COI）），我們可以確定白色虛線以內(nèi)的數(shù)據(jù)是準(zhǔn)確的。在陰影區(qū)域的白線之外，時頻圖中的信息應(yīng)被視為不準(zhǔn)確的信息，因為可能存在邊緣效應(yīng)。可以從時間分辨率的角度去看在低頻對應(yīng)較大尺度小波從而時間分辨率較差。

詳細解釋請查看：Boundary Effects and the Cone of Influencehttps://www.mathworks.com/help/releases/R2020b/wavelet/ug/boundary-effects-and-the-cone-of-influence.html通過 3D 可視化查看小波幅值的增長速率，同時也可以將 CWT 的結(jié)果和真實瞬時頻率結(jié)果進行對比可以發(fā)現(xiàn)：CWT 結(jié)果和真實瞬時頻率一致性較好。

MATLAB 中提供了除了上述時頻變換，還有例如 Constant-Q Gabor Transform，Empirical Mode Decomposition and Hilbert-Huang Transform 等等。歡迎大家查看文檔搜索 https://ww2.mathworks.cn/help/。

振動分析

我們再回到振動場景本身，在進行時頻分析時，我們經(jīng)常使用階次分析來確定發(fā)生在旋轉(zhuǎn)機械中的頻譜成分，從時域波形中追蹤和提取階次，計算不同階次的平均譜值，通過估計頻響函數(shù)、固有頻率、阻尼比和模態(tài)振型進行試驗?zāi)B(tài)分析，用時間同步平均法相干去除噪聲，用包絡(luò)譜分析磨損，為疲勞分析生成高周期雨流計數(shù)等等。

我們通過一個示例，對直升機在加速和減速過程機艙內(nèi)的加速度傳感器數(shù)據(jù)進行階次分析從而確定振動源并進行優(yōu)化來演示分析過程。當(dāng)設(shè)備的轉(zhuǎn)速隨時間一直變化時，階次分析可以用來計算噪聲或振動。每個階數(shù)對應(yīng)某個參考轉(zhuǎn)速的倍頻。例如，一個信號的頻率如果等于發(fā)動機轉(zhuǎn)頻的 2 倍，那階數(shù)就是 2。

我們通過仿真得到直升機在加速和減速過程機艙內(nèi)的加速度傳感器數(shù)據(jù)。直升機有很多旋轉(zhuǎn)部件，包括發(fā)動機，齒輪箱，主旋翼和尾翼。每個部件都以一個和主發(fā)動機固定的速比運轉(zhuǎn)，每個部件都可能導(dǎo)致有害振動。示例中的信號是一個時域電壓信號 vib， 按 fs = 500Hz 的頻率進行采樣。

數(shù)據(jù)還包含渦輪發(fā)動機的轉(zhuǎn)速（rpm） vs 時間關(guān)系。通過數(shù)據(jù)可視化可以看到發(fā)動機轉(zhuǎn)速在爬升和滑行過程變化，同時振動加速度幅值也跟轉(zhuǎn)速有關(guān)。使用 RPM-Frequency 圖可視化數(shù)據(jù)為了能對振動信號進行時頻域的分析，我們可以使用 rpmfreqmap。這個函數(shù)計算信號的短時傅里葉變換生成 RPM-Frequency 圖譜。

rpmfreqmap 函數(shù)生成一幅時頻圖像同時還有對應(yīng)的轉(zhuǎn)速時間圖像，還有下面的幾個數(shù)值參數(shù)。圖中幅值默認(rèn)使用均方根（RMS）。當(dāng)然也可以通過選項設(shè)置來使用其他幅值指標(biāo)，例如峰值或功率值。瀑布圖按鈕可以生成三維的瀑布圖。

可以看到瀑布途中很多頻率對應(yīng)的幅值會隨發(fā)動機轉(zhuǎn)速變化而變化。這說明這些頻率是發(fā)動機轉(zhuǎn)頻的階次頻率。在 RPM 峰值也對應(yīng)著振動高幅值的成分，這些成分主要集中在 20-30Hz 之間。可以看到在當(dāng)前默認(rèn)的頻率分辨率（fs/128=3.906Hz，rpmfreqmap 函數(shù)默認(rèn)將采樣率做 128 等分作為分辨率）不足以分辨一些轉(zhuǎn)速峰值處的低頻成分，因此，我們嘗試將頻率分辨率調(diào)小一些，設(shè)置為 1Hz。從結(jié)果可以看出1Hz分辨率的 RPM-frequency 圖可以清晰分辨出這些成分。

從現(xiàn)在的結(jié)果來看，在 RPM 峰值處對應(yīng)的低頻成分可以被分辨出來，但同時我們也發(fā)現(xiàn)在轉(zhuǎn)速快速變化時出現(xiàn)了嚴(yán)重的頻率拖尾現(xiàn)象。這是因為在每個時間窗內(nèi)隨著發(fā)動機轉(zhuǎn)速增加或變小，振動頻率階數(shù)也在變化，覆蓋一個比較大的范圍，從而導(dǎo)致一個較大的頻譜帶寬。

分辨率越高，要求時間窗采樣時長也越長，其中包含的頻率覆蓋范圍越廣，模糊現(xiàn)象越明顯。在直升機起飛或減速滑行階段，通過提升分辨率會導(dǎo)致頻率拖尾現(xiàn)象更嚴(yán)重。在這時，我么可以考慮使用階次圖來避免這種分辨率和包含頻率帶寬大小的矛盾。

使用 RPM 階次圖可視化數(shù)據(jù)函數(shù) rpmordermap 生成階次頻譜圖與 RPM 圖，用于階次分析。由于每個階次是參考旋轉(zhuǎn)速度的固定倍數(shù)，因此階次圖包含一條條直的階次路徑，都是 RPM 的函數(shù)。函數(shù) rpmordermap 與 rpmfreqmap 使用相同的參數(shù)，對應(yīng)的頻域軸現(xiàn)在是階次，而不是頻率。使用 rpmfreqmap 可視化直升機數(shù)據(jù)的階次圖。指定階次分辨率為 0.005 倍基頻。

階次圖包含每個階次的直線路徑，每一階對應(yīng)著發(fā)動機轉(zhuǎn)速的倍數(shù)對應(yīng)的振動。階次圖可以清楚的看到每個頻率成分與發(fā)動機轉(zhuǎn)速的關(guān)系。與 RPM-頻率圖相比，頻率拖尾現(xiàn)象顯著被抑制。使用平均階次頻譜確定峰值階次接下來，確定階次圖的峰值位置。尋找主旋翼和尾翼階次的整數(shù)倍的階次，對應(yīng)著這些旋翼產(chǎn)生的振動。函數(shù) rpmordermap 返回包含各階次的階次圖和與時間對應(yīng)的 RPM 值。通過分析數(shù)據(jù)以確定直升機艙內(nèi)高振幅振動的對應(yīng)階次。計算并返回階次圖譜數(shù)據(jù)。

接下來，使用 orderspectrum 計算并繪制 map 的平均譜。該函數(shù)接受 rpmordermap 生成的階次圖表作為輸入，并對時域取平均。

返回平均頻譜值，并調(diào)用 findpeaks 以返回兩個最高峰值的位置。

在圖中大約 0.05 階處可以看到兩個間隔很近的主峰。階次小于 1，因為振動頻率低于發(fā)動機轉(zhuǎn)速。分析峰值階次隨時間的變化接下來，使用 ordertrack 求峰值階次的幅值隨時間的變化。使用 map 作為輸入，通過不帶輸出參數(shù)調(diào)用 ordertrack 來繪制兩個峰值階次的振幅。

隨著發(fā)動機轉(zhuǎn)速的增加，兩個階次的振幅都會增加。接下來，使用 orderwaveform 提取每個階次對應(yīng)的時域波形。提取出來的時域波形可以直接與原始振動信號進行比較。orderwaveform 使用 Vold-Kalman 濾波器提取特定階次的時域波形。將兩個峰值階次對應(yīng)的時域波形之和與原始信號進行比較。

減少機艙振動為了確定機艙振動的來源，我們可以將振動峰值對應(yīng)的階次和直升機旋轉(zhuǎn)部件的階次配合起來看。

可以看到最大振幅分量的頻率是主旋翼頻率的四倍。而主旋翼有四個葉片，很可能是這種振動的來源，通常對于每個旋翼有 N 葉片的直升機來說，主轉(zhuǎn)速的 N 階振動是很常見的。同樣，第二大分量位于尾旋翼速度的一階次處，表明振動可能源于尾旋翼。在對主旋翼和尾旋翼進行軌跡和平衡調(diào)整后，采集新數(shù)據(jù)集。加載新數(shù)據(jù)集并比較調(diào)整前后的階次頻譜。

主峰的振幅現(xiàn)在大幅降低。

總結(jié)

此示例使用階次分析來確定直升機的主旋翼和尾翼是否為機艙內(nèi)高振幅振動的潛在來源。首先，使用了 rpmfreqmap 和 rpmordermap 對階次進行可視化。RPM-階次圖在整個 RPM 范圍內(nèi)實現(xiàn)了階次分離，且消除了在 RPM-頻率圖中出現(xiàn)的頻率拖尾。

rpmordermap 最適合可視化在發(fā)動機加速和減速期間低 RPM 時的振動分量。接下來，該示例使用 orderspectrum 確定峰值階次，使用 ordertrack 可視化峰值階次的振幅隨時間的變化情況，使用 orderwaveform 提取峰值階次的時域波形。

最大的振幅振動分量出現(xiàn)在主旋翼旋轉(zhuǎn)頻率的四倍處，表明主旋翼葉片不平衡。第二大分量出現(xiàn)在尾旋翼的旋轉(zhuǎn)頻率處。調(diào)整旋翼后，振動程度得以降低。以后，我們還將陸續(xù)為大家介紹此系列的其他兩個主題：模態(tài)分析與預(yù)測性維護：人工智能算法應(yīng)用。

編輯：jq