振動(dòng)信號(hào)處理與多功能時(shí)頻分析系統(tǒng)的研究

孫苗鐘，殷 磊，談炳發(fā)，崔世海，徐元利

(天津市輕工與食品工程機(jī)械裝備集成設(shè)計(jì)與在線監(jiān)控重點(diǎn)實(shí)驗(yàn)室，天津科技大學(xué)機(jī)械工程學(xué)院，天津 300222)

振動(dòng)信號(hào)處理與分析的研究能優(yōu)化設(shè)計(jì)生產(chǎn)，實(shí)時(shí)監(jiān)測(cè)設(shè)備運(yùn)行狀態(tài)，及時(shí)發(fā)現(xiàn)故障和實(shí)施故障診斷，保障生產(chǎn)質(zhì)量，提高設(shè)備可靠性等[1]．如何高效地開發(fā)出集信號(hào)處理與分析方法于一體的振動(dòng)測(cè)試系統(tǒng)，一直是科研人員研究開發(fā)的技術(shù)難題．不少學(xué)者對(duì)振動(dòng)測(cè)試系統(tǒng)進(jìn)行研究與設(shè)計(jì)，取得了一些成果．如屈威等[2]采用 LabVIEW 開發(fā)出多通道的振動(dòng)測(cè)試系統(tǒng)，用于高鐵橋梁的振動(dòng)測(cè)試，但該系統(tǒng)不具備時(shí)頻分析功能；唐進(jìn)元等[3]采用LabVIEW生產(chǎn)者-消費(fèi)者模式開發(fā)出用于齒輪傳動(dòng)的振動(dòng)測(cè)試系統(tǒng)，分析功能比較單一；唐奕等[4]使用 LabVIEW 編程語(yǔ)言設(shè)計(jì)了振動(dòng)信號(hào)分析系統(tǒng)，時(shí)頻分析部分僅實(shí)現(xiàn)短時(shí)快速傅里葉變換(STFT)和 Helbert-Huang變換(HHT)．以上開發(fā)的振動(dòng)測(cè)試系統(tǒng)運(yùn)用LabVIEW及其他可視化編程工具，如VB、VC++等編程平臺(tái)，這些平臺(tái)對(duì)各種振動(dòng)測(cè)試信號(hào)處理和分析方法編程難度較大[5]．

MATLAB是集數(shù)值計(jì)算、時(shí)頻分析、信號(hào)處理和圖形分析等強(qiáng)大功能于一體的編程平臺(tái)，提供大量的函數(shù)庫(kù)，具有較強(qiáng)的信號(hào)分析和繪圖能力，可以生成面向?qū)ο缶幊痰慕M件對(duì)象模型(COM)組件；使用MATLAB對(duì)信號(hào)處理與分析(尤其是時(shí)頻分析)的算法設(shè)計(jì)效率高，但不能脫離 MATLAB集成環(huán)境工作，運(yùn)行速度慢，界面開發(fā)能力差[6]．Delphi是面向?qū)ο蟮目梢暬幊陶Z(yǔ)言，適合圖形界面的開發(fā)，但要實(shí)現(xiàn)對(duì)各種信號(hào)處理與分析方法的編程十分困難．如果將兩者混合編程，就可充分利用 MATLAB和Delphi各自的優(yōu)點(diǎn)，由MATLAB完成信號(hào)處理與分析過程和繪圖功能，在 Delphi語(yǔ)言開發(fā)的程序界面中嵌入 MATLAB COM 組件[7–8]．本文通過MATLAB與 Delphi混合編程，基于 COM 組件技術(shù)設(shè)計(jì)與實(shí)現(xiàn)了集32種信號(hào)處理與分析功能方法的振動(dòng)測(cè)試系統(tǒng)，可以完成信號(hào)數(shù)據(jù)回放、信號(hào)數(shù)字濾波處理、平穩(wěn)信號(hào)時(shí)域與頻域分析和非平穩(wěn)信號(hào)時(shí)頻聯(lián)合分析的方法．在時(shí)頻分析方法中，提出了時(shí)頻組合方法并實(shí)施完成，如經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)與快速傅里葉變換(FFT)組合(EMD-FFT)等，這些組合方法能更徹底地提取非平穩(wěn)信號(hào)的特征頻率成分，分析更為方便且結(jié)果表達(dá)更為深刻．該系統(tǒng)可通過擴(kuò)充 COM組件增加分析功能，且系統(tǒng)軟件可以脫離于MATLAB環(huán)境獨(dú)立運(yùn)行于Windows操作系統(tǒng)．

1 系統(tǒng)編程框圖與功能設(shè)計(jì)

系統(tǒng)的編程框圖如圖 1所示．“振動(dòng)信號(hào)處理與多功能時(shí)頻分析”系統(tǒng)由信號(hào)分析設(shè)置、信號(hào)數(shù)字濾波處理、數(shù)據(jù)回放、平穩(wěn)信號(hào)分析和非平穩(wěn)信號(hào)分析組成．

圖1 系統(tǒng)的編程框圖Fig. 1 Flow chart of system programming

1.1 信號(hào)分析設(shè)置

信號(hào)分析設(shè)置分為4類分析：在線采集單通道分析、在線采集雙通道分析、非實(shí)時(shí)采集單通道分析及非實(shí)時(shí)采集雙通道分析．在線采集分析指的是對(duì)實(shí)時(shí)采集后的信號(hào)分析．非實(shí)時(shí)采集分析指的是對(duì)以前數(shù)據(jù)的信號(hào)分析(包括分析后保存的數(shù)據(jù)文件)．信號(hào)分析設(shè)置中要求輸入分析通道數(shù)、采樣頻率、分析點(diǎn)數(shù)、輸入文件名、輸出文件名(要求輸出)．

1.2 信號(hào)數(shù)字濾波處理

采集后的信號(hào)可進(jìn)行選擇性頻率的數(shù)字濾波，增加信號(hào)的預(yù)處理能力．運(yùn)用 MATLAB平臺(tái)中的頻域?yàn)V波方法和時(shí)域?yàn)V波方法，能進(jìn)行高通、低通、帶阻和帶通 4種方式濾波，時(shí)域?yàn)V波可選擇巴特沃斯、切比雪夫Ⅰ和切比雪夫Ⅱ 3種方法．需要輸入不同的截止頻率、通帶波動(dòng)系數(shù)和阻帶波動(dòng)系數(shù)．

1.3 數(shù)據(jù)回放

要對(duì)已存測(cè)試信號(hào)進(jìn)行不丟點(diǎn)和準(zhǔn)確性的驗(yàn)證，就需要進(jìn)行數(shù)據(jù)重新回放，觀察信號(hào)整個(gè)時(shí)間歷程．圖 2中的波形回放主窗口內(nèi)是采集到的 8個(gè)通道波形回放．即是一臺(tái)雙通道 SF2-DDS函數(shù)信號(hào)發(fā)生器分別發(fā)生 A路 100Hz、峰值電壓 1V的正弦波和B路100Hz、峰值電壓1.5V的方波，A路正弦波同時(shí)接數(shù)據(jù)采集卡的 0、2、4、6通道、B 路方波同時(shí)接數(shù)據(jù)采集卡的 1、3、5、7通道．通過回放各個(gè)波形完好無(wú)缺．下面是數(shù)據(jù)回放窗口，可檢查每個(gè)采集點(diǎn)的數(shù)據(jù)大小．其左邊是信號(hào)回放的信息和采樣信息，包括回放速度、文件長(zhǎng)度、通道點(diǎn)數(shù)、屏幕位置、量程、通道頻率和通道總數(shù)．本次選擇回放速度點(diǎn)數(shù)為20．總數(shù)據(jù)文件長(zhǎng)度達(dá) 139420點(diǎn)，每個(gè)通道數(shù)據(jù)為 17427點(diǎn)．

圖2 振動(dòng)信號(hào)處理及多功能時(shí)頻分析系統(tǒng)界面Fig. 2 System interface of vibration signal processing and multi-function time-frequency analysis

1.4 平穩(wěn)信號(hào)分析

平穩(wěn)信號(hào)分析又分為單通道分析和雙通道分析，單通道分析的方法有實(shí)時(shí)信號(hào)波形顯示、快速傅里葉變換(FFT)、自相關(guān)函數(shù)、自功率譜密度、倒譜和概率密度函數(shù)；雙通道分析的方法有互相關(guān)函數(shù)、互功率譜密度、頻率響應(yīng)函數(shù)和相干函數(shù)．

1.5 非平穩(wěn)信號(hào)分析

該部分為系統(tǒng)分析的核心，分為變換與分布和組合分析兩部分，變換與分布的方法有短時(shí)快速傅里葉變換(STFT)、Helbert-Huang變換(HHT)、魏格爾分布(WVD)、偽魏格爾分布(PWVD)、平滑偽魏格爾分布(SPWVD)和 Choi-Williams分布(CWD)．組合分析的方法有小波分解與 FFT、小波降噪與 FFT、小波包分解、小波包降噪與FFT、EMD-FFT、EEMD-FFT、EMD-WVD和EMD-WVD的疊加．

本系統(tǒng)還具有保存或打印分析波形的功能．

2 混合編程過程

2.1 COM組件技術(shù)

組件對(duì)象模型(COM)為Windows提供了一種面向?qū)ο蟮?、可擴(kuò)充的通信協(xié)議，通過此協(xié)議，任何編程語(yǔ)言、不同平臺(tái)和彼此獨(dú)立的對(duì)象都可以實(shí)現(xiàn)連通．MATLAB通過調(diào)用外部編譯器對(duì)函數(shù)文件編譯，使用COM Builder生成COM組件的動(dòng)態(tài)鏈接庫(kù)(DLL)文件，Delphi通過安裝、調(diào)用該 DLL文件，實(shí)現(xiàn)MATLAB和Delphi之間的數(shù)據(jù)通信，從而實(shí)現(xiàn)混合編程[8]．

2.2 生成DLL文件并注冊(cè)

MATLAB平臺(tái)擁有各種信號(hào)處理與分析函數(shù)和強(qiáng)大的繪圖函數(shù)．可編寫各種函數(shù) M 文件(子程序)，該文件可以接受和返回參數(shù)．在MATLAB命令窗口中輸入 mex-setup，選擇外部編譯器 Microsoft Visual C++ 6.0．再輸入命令 comtool，彈出 COM 組件創(chuàng)建窗口，輸入組件名字，點(diǎn)擊菜單中的 Build COM Object進(jìn)行編譯，編譯成功后形成相應(yīng)的動(dòng)態(tài)鏈接庫(kù)(DLL)文件，該文件會(huì)在 Windows系統(tǒng)中完成注冊(cè)[9]．

2.3 安裝DLL文件

在 Delphi平臺(tái)中，點(diǎn)擊菜單選項(xiàng)中的 Project Import Type Library，彈出窗會(huì)顯示上面注冊(cè)后的各種類型庫(kù)(DLL)文件，選擇所需DLL文件，安裝完畢后就會(huì)在ActiveX類中出現(xiàn)相應(yīng)的組件．雙擊組件會(huì)在 Delphi的編程窗口中出現(xiàn)相應(yīng)的函數(shù)類，如圖 3為系統(tǒng)主界面窗口調(diào)用各函數(shù)類情況，在 Delphi的子程序中調(diào)用相應(yīng)的分析方法函數(shù)類，實(shí)現(xiàn)該方法與MATLAB數(shù)據(jù)接口的連接與通信．

圖3 系統(tǒng)編程主界面窗口調(diào)用各函數(shù)類Fig. 3 Function classes in the main system programming interface

混合編程的流程框圖如圖4所示．

圖4 混合編程的流程框圖Fig. 4 Flow chart of mixed programming

3 時(shí)頻分析的組合方法實(shí)施

在非平穩(wěn)信號(hào)的時(shí)頻分析中，為了更好地提取非平穩(wěn)信號(hào)的頻率特征信息，在實(shí)施小波變換(WT)、EMD和EEMD等方法時(shí)采用2種或3種方法的組合(組合方法)分析，如 WT-FFT、EMD-WVD 和EEMD-FFT等，這可使結(jié)果分析更為深刻和清晰[10].以經(jīng)驗(yàn)?zāi)Ｊ椒纸馀c魏格爾分布(EMD-WVD)組合方法為例說(shuō)明其實(shí)施過程，該方法采用了EMD、相關(guān)系數(shù)法和WVD 3種方法進(jìn)行組合對(duì)非平穩(wěn)信號(hào)分析.

經(jīng)驗(yàn)?zāi)Ｊ椒纸庥捎诜N種原因存在過分解現(xiàn)象，過分解產(chǎn)生的原因主要有以下幾方面：局部均值的數(shù)值計(jì)算方法的插值誤差、邊界效應(yīng)的影響、終止篩選標(biāo)準(zhǔn)的不嚴(yán)格等．過分解使得到的分量個(gè)數(shù)比原信號(hào)組成分量多，把非原信號(hào)組成成分的分量稱為“偽分量”[10]．用相關(guān)系數(shù)法來(lái)剔除“偽分量”，采用 EMD分解后的各個(gè)內(nèi)稟模式函數(shù)(IMF)與原信號(hào)進(jìn)行相關(guān)得到各個(gè)相關(guān)系數(shù)，當(dāng)相關(guān)系數(shù)很小時(shí)，就可判斷該分量為偽分量，可能剔除．由于 WVD在非平穩(wěn)信號(hào)分析具有較好的時(shí)頻聚集性和分辨率[11]，所以再用WVD方法對(duì)剩余的IMF進(jìn)行分析，得到各WVD分布的時(shí)頻譜圖，分析結(jié)果更為清晰，信號(hào)被分解后的各頻率成分更清晰地表達(dá)出來(lái)．

4 仿真信號(hào)驗(yàn)證系統(tǒng)

為了驗(yàn)證系統(tǒng)的分析效果與可靠性，運(yùn)用系統(tǒng)中多種分析方法對(duì)一仿真信號(hào)進(jìn)行分析比較．

設(shè)仿真信號(hào)的解析表達(dá)式為

該信號(hào)是由一個(gè)線性趨勢(shì)項(xiàng) t，兩個(gè)基頻分別為20Hz、200Hz的正弦波和一個(gè)基頻100Hz并被基頻20Hz正弦波調(diào)頻的余弦波疊加組成的非平穩(wěn)信號(hào)．對(duì)其調(diào)頻部分的頻率分析，得角頻率：

由式(3)可得其頻率變化范圍為 80～120Hz，且以余弦波規(guī)律變化．

仿真信號(hào)的波形和 FFT分析后的頻譜圖(輸出打印得到)如圖 5所示．輸入采樣點(diǎn)數(shù) 1000，采樣頻率 1000Hz．

從頻譜圖可知 20、100、200Hz基頻成分出現(xiàn)分明，這是符合原信號(hào)成分的．由于在FFT計(jì)算中對(duì)信號(hào)進(jìn)行截?cái)?，并采用了窗函?shù)(本系統(tǒng)采用漢寧窗)，會(huì)產(chǎn)生截?cái)嗪蜐B漏效應(yīng)，明顯地在 100Hz處 20Hz正弦波調(diào)頻成分表達(dá)不清楚，在頻率為零附近的低頻成分是不符合趨勢(shì)項(xiàng)分析要求的，說(shuō)明 FFT對(duì)非平穩(wěn)信號(hào)分析存在缺陷．

應(yīng)用時(shí)頻組合方法(EMD-WVD)對(duì)仿真信號(hào)EMD分解后各IMF波形以及對(duì)各IMF用WVD分析后的時(shí)頻譜圖如圖 6所示．在圖 6(b)中標(biāo)出了各IMF與原信號(hào)的相關(guān)系數(shù)，分別為：0.6361、0.5489、0.5544和 0.0327．顯然，IMF1、IMF2和 IMF3與原信號(hào)相關(guān)，IMF4與原信號(hào)幾乎無(wú)相關(guān)．從 WVD分析可知，IMF1、IMF3分別是 200Hz、20Hz正弦波，IMF2是基頻 100Hz被基頻 20Hz正弦波調(diào)頻的余弦波，IMF4為偽 IMF，應(yīng)剔除，Res是線性趨勢(shì)項(xiàng)t．可見采用這種組合方法來(lái)分析非平穩(wěn)仿真信號(hào)顯示其優(yōu)點(diǎn)．

圖5 仿真信號(hào)的波形及FFT分析后的頻譜Fig. 5 Simulation signal waveform and its spectrum

圖6 仿真信號(hào)用EMD-WVD分析后的結(jié)果圖Fig. 6 Results of simulation signal analyzed with EMD-WVD

5 實(shí)例分析

一臺(tái) S3SL型砂輪機(jī)轉(zhuǎn)速為 2850r/min，即工頻f0為 47.5Hz，振動(dòng)劇烈噪聲大，達(dá)到 96.2dB．測(cè)試系統(tǒng)的硬件組成框圖如圖7所示．

用傳感器提取砂輪機(jī)在工頻轉(zhuǎn)速下其機(jī)殼上的垂直振動(dòng)信號(hào)，選擇采樣頻率為 1000Hz，取數(shù)據(jù)長(zhǎng)度1000點(diǎn)，運(yùn)用本系統(tǒng)分析．砂輪機(jī)原信號(hào)、頻域方法濾波后的信號(hào)及其頻譜圖如圖 8所示．頻域方法濾波是利用 FFT快速算法對(duì)輸入信號(hào)采樣數(shù)據(jù)進(jìn)行離散傅里葉變換分析其頻譜，根據(jù)濾波要求，將需要濾除的頻率成分直接設(shè)置成零或加漸變過渡帶后再設(shè)置成零．濾波時(shí)采用低通濾波方式，截止頻率取300Hz，由于其頻率特性，對(duì)頻域數(shù)據(jù)的突然截?cái)嘣斐傻淖V泄會(huì)造成濾波后的時(shí)域信號(hào)出現(xiàn)失真變形，但頻域方法具有較好的頻率選擇性和靈活性．砂輪機(jī)原信號(hào)、小波濾波后的信號(hào)及其頻譜圖如圖 9．小波分解時(shí)，對(duì)于不同的信號(hào)，應(yīng)選擇不同的小波基．小波變換相當(dāng)于通過小波的尺度因子和時(shí)移因子的變化去觀察信號(hào)，尺寸因子大可觀察信號(hào)的總體，尺寸因子小可觀察信號(hào)的細(xì)節(jié)．本系統(tǒng)選擇廣泛使用的Daubechies(db)小波基，適中選擇 5層小波，適合大多數(shù)信號(hào)．分解后取對(duì)第一層重構(gòu)系數(shù) a1(信號(hào)的低頻成分)，對(duì)上面兩圖中信號(hào)濾波后的頻譜圖進(jìn)行比較，結(jié)果吻合很好，證實(shí)了兩種濾波方法的可靠性．從濾波后的信號(hào)頻譜圖中反映了砂輪機(jī)振動(dòng)信號(hào)低頻區(qū)域內(nèi)各倍頻的頻率特征成分信息，即低頻處的1～5倍頻峰值十分明顯，尤其是1倍頻(1f0)．

圖7 砂輪機(jī)測(cè)試系統(tǒng)的硬件組成框圖Fig. 7 Hardware framework of the testing system of the grinding machine

圖9 砂輪機(jī)原信號(hào)、小波濾波后信號(hào)及其頻譜圖Fig. 9 Grinding machine original signal，wavelet denoise signal and their spectrograms

使用短時(shí)快速傅里葉變換(STFT)和魏格爾分布(WVD)兩種時(shí)頻方法對(duì)砂輪機(jī)的振動(dòng)信號(hào)進(jìn)行分析，這兩種方法分析結(jié)果如圖10所示．

圖10 砂輪機(jī)振動(dòng)信號(hào)STFT、WVD時(shí)頻譜圖Fig. 10 STFT and WVD time-frequency spectrograms of the grinding machine

經(jīng)STFT分析后，得到1倍頻約50Hz頻率成分最為突出，而3倍頻約150Hz的頻率成分也較明顯，2、4、5倍頻的頻率成分隱約可見，在300～400Hz范圍內(nèi)存在比較多頻率成分．由 WVD分析結(jié)果可見，頻率分辨率明顯增加，但出現(xiàn)不可知的頻率成分交叉項(xiàng)，1倍頻的頻率成分也最為突出．

上述幾種方法的分析得到砂輪機(jī)振動(dòng)信號(hào) 1、3倍頻的特征頻率成分突出，根據(jù)旋轉(zhuǎn)機(jī)械故障診斷理論[12]可知，砂輪機(jī)運(yùn)行時(shí)存在嚴(yán)重偏心而造成運(yùn)轉(zhuǎn)不平衡，導(dǎo)致振動(dòng)劇烈．測(cè)后更換砂輪，振動(dòng)與噪聲恢復(fù)正常．

6 結(jié) 語(yǔ)

應(yīng)用 COM 組件技術(shù)，可以將 MATLAB強(qiáng)大的信號(hào)處理與分析功能、圖形處理能力和 Delphi靈活高效的程序設(shè)計(jì)能力相結(jié)合，實(shí)現(xiàn)混合編程．開發(fā)集信號(hào)處理與多功能分析于一體的振動(dòng)測(cè)試系統(tǒng)可以脫離于 MATLAB環(huán)境，以執(zhí)行文件的形式獨(dú)立運(yùn)行于 Windows環(huán)境．本系統(tǒng)為更好提取非平穩(wěn)信號(hào)頻率特征成分信息提出了多種新的時(shí)頻組合方法，是對(duì)時(shí)頻分析進(jìn)行補(bǔ)充和完善．通過仿真信號(hào)分析和實(shí)例測(cè)試，驗(yàn)證了系統(tǒng)的可靠性，可應(yīng)用于振動(dòng)信號(hào)處理與分析，尤其是非平穩(wěn)信號(hào)的分析．

該系統(tǒng)主要優(yōu)點(diǎn)：實(shí)現(xiàn) 32通道信號(hào)同時(shí)大容量數(shù)據(jù)采集顯示與存盤；對(duì)采集后的數(shù)據(jù)可以回放分析；根據(jù)信號(hào)分析要求選擇不同的數(shù)字濾波處理；根據(jù)非平穩(wěn)信號(hào)的特點(diǎn)，選擇合適的組合時(shí)頻分析，使分析結(jié)果更為清晰．該系統(tǒng)存在的缺點(diǎn)：缺少邊采集信號(hào)邊實(shí)時(shí)信號(hào)分析功能；對(duì)于有關(guān)時(shí)頻分析方法中，輸入控制參數(shù)不夠全面，分析有局限性，如小波分析的基函數(shù)、層次數(shù)，EMD 分析中的標(biāo)準(zhǔn)差、循環(huán)次數(shù)等．

本系統(tǒng)具有在線幫助功能，對(duì) MATLAB平臺(tái)的每個(gè)函數(shù)應(yīng)用都有較好的解釋和理論指導(dǎo)，對(duì)各個(gè)分析函數(shù)使用時(shí)作了詳細(xì)的說(shuō)明．

在使用本系統(tǒng)時(shí)，強(qiáng)調(diào)說(shuō)明對(duì)于不同振動(dòng)信號(hào)進(jìn)行預(yù)處理時(shí)，分析要求不同，處理是不同的．如對(duì)軸承的振動(dòng)信號(hào)不應(yīng)進(jìn)行濾波，要保存其信號(hào)的內(nèi)部頻率成分．