姜佩賀，郭 剛，吳中杰，位寅生

(1.煙臺(tái)大學(xué)光電信息科學(xué)技術(shù)學(xué)院，山東煙臺(tái) 264005；2.哈爾濱工業(yè)大學(xué)電子與信息工程學(xué)院，黑龍江哈爾濱 150001)

0 引言

通過(guò)分析振動(dòng)信號(hào)可以提取振動(dòng)源的特征，振動(dòng)信號(hào)處理被廣泛應(yīng)用在參數(shù)檢測(cè)、質(zhì)量評(píng)價(jià)、狀態(tài)檢測(cè)和故障診斷等領(lǐng)域[1-3]。獲取振動(dòng)信號(hào)是進(jìn)行信號(hào)處理的第一步，在測(cè)量機(jī)械振動(dòng)時(shí)，實(shí)時(shí)加速度信號(hào)可以使用加速度傳感器直接測(cè)量，而振動(dòng)速度、位移以及頻率的測(cè)量較困難，但更能夠體現(xiàn)振動(dòng)源的特性。例如，在汽車懸掛系統(tǒng)設(shè)計(jì)中，振動(dòng)速度和位移用于評(píng)價(jià)減震器性能；在故障檢測(cè)中，振動(dòng)頻率反映出故障的部位。因此，在得到加速度信號(hào)后，通過(guò)信號(hào)處理技術(shù)將加速器準(zhǔn)確地轉(zhuǎn)換為另外3種信號(hào)是進(jìn)行振動(dòng)分析的關(guān)鍵。

對(duì)于振動(dòng)信號(hào)處理，周小祥[4]等提出在時(shí)域使用最小二乘法去趨勢(shì)項(xiàng)的方法；顧明坤[5]等分析了3種加速度信號(hào)通過(guò)積分獲得速度和位移信號(hào)的方法，并對(duì)3種方法進(jìn)行評(píng)估；胡玉梅[6]等對(duì)加速度信號(hào)在頻域內(nèi)直接積分，利用積分精度控制方程保證積分精度；牟玉喆[7]等對(duì)振動(dòng)信號(hào)處理算法在LabVIEW中進(jìn)行實(shí)現(xiàn)。此外，還有部分研究者采用Wigner-Ville分布[8]、譜分析、小波分析[9]、盲源分離[10]、Hilbert-Huang變換[11]和高階統(tǒng)計(jì)量分析等方法對(duì)振動(dòng)信號(hào)展開(kāi)分析。

總結(jié)發(fā)現(xiàn)，現(xiàn)有研究大都側(cè)重于算法本身，算法需要借助于PC機(jī)在MATLAB或LabVIEW等大型軟件中完成，而在很多場(chǎng)合，特別是在物聯(lián)網(wǎng)時(shí)代，振動(dòng)分析經(jīng)常需要以嵌入式傳感器的形式在現(xiàn)場(chǎng)實(shí)時(shí)安裝、實(shí)時(shí)采集、實(shí)時(shí)處理，實(shí)時(shí)得到振動(dòng)源特性。在嵌入式系統(tǒng)中進(jìn)行振動(dòng)信號(hào)處理，可以很好地解決該問(wèn)題。

基于此，首先使用MATLAB評(píng)估基于去趨勢(shì)項(xiàng)的時(shí)域分析和基于快速傅里葉變換(FFT)的頻域分析兩種振動(dòng)信號(hào)方法；然后將兩種分析方法分別在嵌入式處理器STM32F103中實(shí)現(xiàn)；最后搭建硬件平臺(tái)，使用加速度傳感器芯片IIS2DH采集加速度信號(hào)，通過(guò)嵌入式系統(tǒng)的算法分析獲得振動(dòng)的速度、位移和頻率信息。

1 時(shí)域處理方法

1.1 方法概述

對(duì)加速度信號(hào)進(jìn)行一次積分可得速度，對(duì)其進(jìn)行二次積分可得位移。但由于噪聲、零漂、非整數(shù)周期采樣等原因，在積分前后要濾波和去除趨勢(shì)項(xiàng)，具體流程如圖1所示。首先，對(duì)采集到的實(shí)時(shí)加速度信號(hào)去直流和FIR濾波得到加速度信號(hào)；然后，對(duì)加速度信號(hào)一次積分后，去除一次趨勢(shì)項(xiàng)，得到速度信號(hào)；最后，對(duì)速度信號(hào)再一次積分，去除二次趨勢(shì)項(xiàng)后得到位移信號(hào)。

圖1 時(shí)域振動(dòng)信號(hào)處理流程框圖

1.2 濾波與積分

采集加速度信號(hào)會(huì)混入低頻和高頻噪聲，對(duì)噪聲的積分累加將使信號(hào)發(fā)生畸變，因而需要濾除噪聲。此外，在不同的應(yīng)用場(chǎng)合，需要的振動(dòng)頻率范圍是不同的，這也需使用濾波器濾除不需要的頻率范圍。采用32階基于窗函數(shù)的FIR濾波器，并將下限和上限截止頻率分別設(shè)置為1 Hz和100 Hz。雖然在相同性能條件下，F(xiàn)IR濾波器階數(shù)要比IIR濾波器高，但FIR濾波器總是穩(wěn)定的，且不涉及復(fù)數(shù)運(yùn)算，更適合在嵌入式系統(tǒng)中實(shí)現(xiàn)。FIR濾波器用差分方程形式可以表示為

(1)

式中：x(n)和y(n)分別為輸入和輸出時(shí)域信號(hào)序列；bk為濾波系數(shù)。

在濾波前，需要去直流，去直流也可以看作是濾波的一部分，是將不需要的直流分量濾除。去直流的方法是求出全部N個(gè)采樣點(diǎn)的平均值，然后每個(gè)采樣點(diǎn)減去該值。該過(guò)程可以表示為

(2)

式中：ai為采集到的加速度信號(hào)；ai′為去直流后的加速度信號(hào)；N為采樣點(diǎn)數(shù)；i的取值范圍為0至N-1。

常見(jiàn)的數(shù)字積分法有梯形積分法和Simpson積分法等，梯形積分法計(jì)算簡(jiǎn)單，精度能夠滿足要求，適合在嵌入式系統(tǒng)中實(shí)現(xiàn)，因此采用梯形積分法。具體公式為

(3)

(4)

式中：a為加速度；v為速度；s為位移；t1為采樣時(shí)間間隔，即采樣頻率的倒數(shù)。

1.3 去除趨勢(shì)項(xiàng)

趨勢(shì)項(xiàng)是由于零漂、噪聲、非正周期采樣等因素經(jīng)積分累加后造成的信號(hào)隨時(shí)間不斷偏離基線的現(xiàn)象，趨勢(shì)項(xiàng)影響信號(hào)的正確性，需要去除。采用最小二乘法擬合去除趨勢(shì)項(xiàng)，具體方法為：對(duì)于速度信號(hào)，對(duì)該信號(hào)進(jìn)行最小二乘法一次擬合，速度信號(hào)與擬合結(jié)果相減，即可將趨勢(shì)項(xiàng)去除。最小二乘法一次擬合的參數(shù)為

(5)

消除一次趨勢(shì)項(xiàng)的計(jì)算公式為

(6)

去除位移信號(hào)趨勢(shì)項(xiàng)的方法與此類似，但由于位移信號(hào)是經(jīng)過(guò)二次積分得到的，因此，需要使用最小二乘法對(duì)位移信號(hào)二次擬合，然后做差。在MATLAB中可以使用polyfit和polyval函數(shù)實(shí)現(xiàn)最小二乘法的參數(shù)計(jì)算與多項(xiàng)式擬合。

1.4 MATLAB仿真實(shí)驗(yàn)與結(jié)果分析

按照?qǐng)D1所示流程，在MATLAB中編寫測(cè)試程序，假設(shè)加速度信號(hào)為a=20sin(2π×25×t)，采樣頻率為400 Hz，經(jīng)本方法得到的速度和位移信號(hào)如圖2所示。

圖2 時(shí)域方法對(duì)正弦加速度信號(hào)的處理結(jié)果

圖2(a)～圖2(c)分別為原始加速度信號(hào)a，速度信號(hào)v和位移信號(hào)s，對(duì)比發(fā)現(xiàn)，使用本方法可以較準(zhǔn)確地恢復(fù)出速度信號(hào)和位移信號(hào)。分別選取理論值與計(jì)算值的極值點(diǎn)，經(jīng)計(jì)算，速度v平均峰值誤差為0.54%，位移s平均峰值誤差為3.99%。由于濾波器的存在，v和s存在相移，但相移并不影響實(shí)際應(yīng)用中的振動(dòng)分析。

2 頻域處理方法

2.1 方法概述

在頻域進(jìn)行振動(dòng)信號(hào)處理的流程如圖3所示。首先對(duì)加速度信號(hào)做FFT得到加速度頻譜并進(jìn)行濾波；然后，對(duì)速度頻譜依次做2次頻域積分，即可得到速度和位移頻譜；最后，使用逆快速傅里葉變換(IFFT)得到速度和位移信號(hào)。

圖3 頻域振動(dòng)信號(hào)處理流程框圖

2.2 離散傅里葉變換與頻域積分

時(shí)域信號(hào)經(jīng)過(guò)離散傅里葉變換(DFT)可以得到頻域信號(hào)，設(shè)加速度信號(hào)a(t)的頻譜為X(k)，則

(7)

(8)

式中N為采樣點(diǎn)數(shù)。

設(shè)a(n)、v(n)、s(n)分別為加速度a(t)、速度v(t)和位移s(t)信號(hào)的離散化表示。每條譜線對(duì)應(yīng)一個(gè)正弦波，有

ak(t)=X(k)ejωkt

(9)

則

(10)

(11)

推導(dǎo)得

(12)

(13)

式中：ωk=2πkΔf；Δf=Fs/n。

由上述推導(dǎo)可知，對(duì)加速度信號(hào)進(jìn)行FFT運(yùn)算后，在頻域表達(dá)式上除以jω即可實(shí)現(xiàn)頻域一次積分，得到速度的頻譜，除以(jω)2即可實(shí)現(xiàn)頻域二次積分，得到位移的頻譜。

2.3 頻域?yàn)V波

濾波的目的是去除噪聲和非需要頻率范圍的振動(dòng)信號(hào)。在頻域，只需將低于下限截止頻率和高于上限截止頻率的頻譜置0，即可完成濾波。與時(shí)域相同，設(shè)置下限和上限截止頻率分別為1 Hz和100 Hz。

2.4 MATLAB仿真實(shí)驗(yàn)與結(jié)果分析

按照?qǐng)D3所示流程編寫MATLAB程序，與時(shí)域相同，假設(shè)加速度信號(hào)為a=20sin(2π×25×t)，采樣頻率為400 Hz，采樣點(diǎn)數(shù)為2 048，即采樣時(shí)長(zhǎng)約為5.11 s，經(jīng)本方法得到的速度和位移信號(hào)頻譜，以及經(jīng)IFFT變換得到的速度和位移信號(hào)如圖4所示。為方便顯示，圖中的時(shí)域信號(hào)僅截取了0.5 s。

圖4(b)為原始加速度a，對(duì)其做FFT后得到的幅值譜見(jiàn)圖4(a)，頻域一次積分后得到速度v幅值譜見(jiàn)圖4(c)，頻域二次積分后得到位移s幅值譜見(jiàn)圖4(e)，對(duì)v和s的頻譜做IFFT變換，得到速度和位移信號(hào)，分別見(jiàn)圖4(d)和圖4(f)。在圖4(d)和圖4(f)中，理論值和計(jì)算值幾乎完全重合，為便于顯示，這里將虛線的理論值分別移位+0.02和+1×10-4。經(jīng)計(jì)算，速度v平均峰值誤差為0.39%，位移s平均峰值誤差為0.79%。位移s幅值譜的極值點(diǎn)所對(duì)應(yīng)的頻率為振動(dòng)頻率，該頻率為25 Hz，與理論值相符。

時(shí)域方法與頻域方法計(jì)算結(jié)果對(duì)比如表1所示，在速度方面計(jì)算精度相差不大，但在位移方面，頻域方法的計(jì)算精度顯著高于時(shí)域方法，即使用頻域方法可以獲得更高的精度，而且計(jì)算量更小。若要進(jìn)一步提高計(jì)算精度，可以增加采樣點(diǎn)個(gè)數(shù)或修改窗函數(shù)。

表1 時(shí)域和頻域方法計(jì)算結(jié)果對(duì)比 %

3 算法在嵌入式系統(tǒng)中的軟硬件實(shí)現(xiàn)及實(shí)驗(yàn)驗(yàn)證

3.1 硬件平臺(tái)

硬件平臺(tái)框圖如圖5所示，加速度傳感器與嵌入式處理器共同構(gòu)成振動(dòng)傳感器。使用加速度傳感器IIS2DH實(shí)時(shí)采集加速度信號(hào)，IIS2DH是一款低功耗高性能三軸加速度傳感器，量程可配置為±2g、±4g、±6g或±16g，采樣頻率可配置為1 Hz到5.3 kHz，對(duì)外接口為SPI和I2C。嵌入式處理器STM32F103采用Cortex-M3內(nèi)核，通過(guò)SPI接口與加速度傳感器IIS2DH通信，在完成信號(hào)處理后，通過(guò)串口發(fā)送處理結(jié)果。在不同的應(yīng)用場(chǎng)合，可以通過(guò)修改硬件和軟件程序采用USB、LoRa、Wi-Fi、藍(lán)牙等不同類型的有線或無(wú)線接口發(fā)送處理結(jié)果。

圖5 振動(dòng)傳感器與測(cè)試平臺(tái)框圖

為方便測(cè)試，搭建簡(jiǎn)易測(cè)試平臺(tái)，振動(dòng)傳感器固定在測(cè)試平臺(tái)上。測(cè)試平臺(tái)的“振動(dòng)”由舵機(jī)產(chǎn)生，舵機(jī)型號(hào)為SG90，扭矩1.6 kg/cm，舵機(jī)在水平方向以圓弧形式往復(fù)運(yùn)動(dòng)，這樣振動(dòng)傳感器IIS2DH即在x和y方向發(fā)生“振動(dòng)”。

3.2 時(shí)域和頻域方法在STM32F103中的實(shí)現(xiàn)

STM32F103內(nèi)部不含F(xiàn)PU等硬件信號(hào)處理單元，所有算法需要用軟件實(shí)現(xiàn)。通過(guò)移植DSP庫(kù)，可以比較容易地實(shí)現(xiàn)時(shí)域和頻域分析算法。

時(shí)域與頻域分析的嵌入式實(shí)現(xiàn)流程與圖1和圖3相同。在時(shí)域分析算法中，濾波器的設(shè)計(jì)較關(guān)鍵，利用MATLAB的fdatool工具生成濾波系數(shù)，調(diào)用arm_fir_f32函數(shù)實(shí)現(xiàn)濾波器的設(shè)計(jì)。面向不同的上限和下限截止頻率，可以生成多組參數(shù)，將參數(shù)以數(shù)組的形式保存在嵌入式處理器中，以方便在不同應(yīng)用場(chǎng)合通過(guò)軟件為系統(tǒng)配置不同的截止頻率；在頻域處理算法中，調(diào)用arm_rfft_fast_f32函數(shù)和arm_rfft_fast_f32完成FFT和IFFT運(yùn)算。

3.3 實(shí)驗(yàn)測(cè)試與結(jié)果分析

設(shè)置振動(dòng)測(cè)試平臺(tái)舵機(jī)的轉(zhuǎn)動(dòng)角為90°，轉(zhuǎn)動(dòng)周期為640 ms。設(shè)置加速度傳感器IIS2DH的量程為±2g，采樣頻率為400 Hz，設(shè)置濾波的下限和上限截止頻率分別為1 Hz和100 Hz。將信號(hào)處理的全部結(jié)果通過(guò)串口發(fā)送至上位機(jī)，在實(shí)際應(yīng)用中，可以根據(jù)場(chǎng)景需求有選擇的對(duì)信號(hào)進(jìn)行取舍。

時(shí)域處理結(jié)果如圖6所示，圖6(a)為IIS2DH采集到的加速度信號(hào)a，圖6(b)為時(shí)域處理得到的速度信號(hào)v，圖6(c)為時(shí)域處理得到的位移信號(hào)s，為實(shí)施對(duì)比驗(yàn)證，將采集到的加速度信號(hào)在MATLAB進(jìn)行數(shù)據(jù)處理，安裝MATLAB的PC機(jī)配置為Intel i5-3320M CPU，8GB內(nèi)存，Windows 7 64 bit操作系統(tǒng)。圖6(b)和圖6(c)中，為方便顯示，將MATLAB處理結(jié)果在縱向進(jìn)行了平移。在處理時(shí)長(zhǎng)方面，MATLAB進(jìn)行一次數(shù)據(jù)處理的時(shí)長(zhǎng)約為5 ms，STM32F103的處理時(shí)長(zhǎng)約為420 ms。

圖6 STM32F103中振動(dòng)加速度時(shí)域處理結(jié)果

頻域處理結(jié)果如圖7所示，圖7(b)為采集到的加速度信號(hào)，對(duì)其進(jìn)行FFT變換得到的幅值譜見(jiàn)圖7(a)，在頻域積分得到的速度和位移幅值譜分別如圖7(c)和圖7(e)所示，進(jìn)行IFFT得到的速度和位移信號(hào)分別見(jiàn)圖7(d)和圖7(f)，為方便顯示，對(duì)MATLAB處理結(jié)果在縱向進(jìn)行了平移。在處理時(shí)長(zhǎng)方面，MATLAB進(jìn)行一次數(shù)據(jù)處理的時(shí)長(zhǎng)約0.4 ms，STM32F103的處理時(shí)長(zhǎng)約150 ms。

圖7 STM32F103中振動(dòng)加速度頻域處理結(jié)果

對(duì)比STM32F103和MATLAB的處理結(jié)果可知，無(wú)論是時(shí)域還是頻域，微小的差別源于串口發(fā)送時(shí)對(duì)浮點(diǎn)數(shù)的截?cái)?，這驗(yàn)證了算法在嵌入式系統(tǒng)中實(shí)現(xiàn)的正確性。嵌入式系統(tǒng)的數(shù)據(jù)處理能力顯著低于PC機(jī)，因而算法在STM32F103的處理時(shí)長(zhǎng)明顯大于MATLAB的處理時(shí)長(zhǎng)，但STM32F103的處理時(shí)長(zhǎng)仍在可接受的范圍內(nèi)，且在實(shí)際振動(dòng)分析中，通常是先采集信號(hào)再進(jìn)行分析，因而對(duì)系統(tǒng)實(shí)時(shí)性要求不高。

4 結(jié)論

小型化、現(xiàn)場(chǎng)化和實(shí)時(shí)化是振動(dòng)分析系統(tǒng)的發(fā)展方向。本研究采用加速度傳感器IIS2DH采集加速度信號(hào)，在STM32F103中分別實(shí)現(xiàn)了基于去趨勢(shì)項(xiàng)的時(shí)域分析方法和基于IFFT的頻域分析方法。測(cè)試結(jié)果表明嵌入式中振動(dòng)分析結(jié)果與MATLAB計(jì)算結(jié)果吻合，算法實(shí)現(xiàn)正確。利用所提出的方法可以將振動(dòng)分析系統(tǒng)以傳感器的形式安裝在測(cè)試現(xiàn)場(chǎng)，實(shí)現(xiàn)實(shí)時(shí)安裝、實(shí)時(shí)采集、實(shí)時(shí)處理、實(shí)時(shí)得到振動(dòng)源特性。