基于組合辛普森積分的振動測試時域積分方法

許曉紅, 倪春鋒, 顏曉青

(沙洲職業(yè)工學(xué)院, 江蘇 張家港 215600)

1 前 言

振動測試是機(jī)械量測試的一個重要領(lǐng)域,關(guān)系到機(jī)器設(shè)備的工作性能和使用壽命,強(qiáng)烈的振動噪聲還會對人體的生理健康產(chǎn)生危害。振動加速度反映振動的沖擊力度,振動速度又叫振動烈度,反映振動能量的大小,振動位移直接反映振動體的位置變化量。理論上對于加速度、速度和位移的測量都有相應(yīng)的傳感器和測試方法,但是振動測量領(lǐng)域往往涉及到大型旋轉(zhuǎn)設(shè)備、運(yùn)動車輛和大型建筑結(jié)構(gòu)或橋梁,這些環(huán)境對于位移測量來說,很難找到位移傳感器安裝所需要的相對靜止結(jié)構(gòu)。而加速度傳感器不需要相對于待測結(jié)構(gòu)靜止的安裝位置,直接將加速度傳感器與待測結(jié)構(gòu)剛性連接即可[1～5]。所以實(shí)際上往往都是通過加速度傳感器獲得振動加速度,然后通過積分獲得振動速度和振動位移,這樣通過一個傳感器得到3個參量,也具有一定的經(jīng)濟(jì)性。

目前常用的積分方法有頻域積分和時域積分。頻域積分先對原始信號作傅里葉變換轉(zhuǎn)化為頻域信號,再對頻域信號進(jìn)行積分運(yùn)算,最后再進(jìn)行傅里葉逆變換得到時域信號。由于信號的處理經(jīng)過了變換域的正變換和逆變換,容易引起截?cái)嗾`差。時域積分對測得的加速度信號直接進(jìn)行積分,會出現(xiàn)干擾趨勢項(xiàng),需要擬合趨勢項(xiàng)以去除信號偏移。在一般的低頻振動測試中,考慮到實(shí)時分析對資源和效率的要求,采用時域積分更為準(zhǔn)確有效[6,7]。本文對時域測量中的數(shù)值積分問題進(jìn)行了研究。

2 測量原理

振動加速度可以通過加速度傳感器獲得,由于振動加速度初始值往往不為零,再有有些加速度傳感器存在著零點(diǎn)偏移,這給通過振動加速度求取振動速度和振動位移的積分過程帶來了累積誤差:

(1)

式中:a(t)為振動加速度;v(t)為振動速度;s(t)為振動位移;a0,a1,a2為積分過程帶來的累計(jì)誤差。

式(1)所示的累積誤差實(shí)際上是一種低頻信號,表示了速度v′(t)和位移s′(t)的變化趨勢,這個趨勢項(xiàng)可以通過曲線擬合的方法提取,常用的擬合方法有多項(xiàng)式擬合法和樣條函數(shù)擬合法[8,9],將積分信號v′(t)和s′(t)分別減去各自的趨勢項(xiàng),就可得到實(shí)際的振動速度和振動位移變化曲線。除此之外,加速度測量通道中的高頻噪聲也對振動速度和振動位移的積分計(jì)算產(chǎn)生影響,所以需要采用一個低通濾波器移除高頻噪聲。整個的數(shù)據(jù)處理與計(jì)算原理方框圖如圖1所示,由此得到振動速度和振動位移曲線,并由波形曲線進(jìn)行參數(shù)評定。

圖1 測量原理框圖

3 組合辛普森積分方法

在整個測量過程中,關(guān)鍵在于積分運(yùn)算。對于被積函數(shù)是確定性函數(shù),可以精確得出積分函數(shù),但實(shí)際振動加速度都是非確定性函數(shù),需要通過數(shù)值積分方法,逼近得出積分函數(shù)。常用的數(shù)值積分方法有矩形公式法和梯形公式法,它們實(shí)際上是用相鄰點(diǎn)之間的矩形和梯形逼近積分區(qū)域,具有計(jì)算簡單的特點(diǎn),但是逼近精度較差。而辛普森積分是基于牛頓-科特斯公式,它是以二次曲線逼近的方式取代矩形或梯形積分公式,具有更高的逼近精度,并且它的計(jì)算公式也相對簡單,是一種工程計(jì)算中的實(shí)用方法[10,11]。牛頓-科特斯公式如下:

(2)

(3)

(4)

對于振動速度和振動位移的測量,實(shí)際需要根據(jù)振動速度和位移的波形曲線,來求取相應(yīng)的參數(shù),例如說峰峰值、均方根值等,這就需要計(jì)算每一點(diǎn)的積分值。而辛普森計(jì)算公式只是計(jì)算了定積分,并且是3個采樣點(diǎn)之間的曲線面積,所以還不能完全直接應(yīng)用。根據(jù)牛頓-萊布尼茲公式,如式(5)所示,一個定積分式的值,就是原函數(shù)在上限的值與原函數(shù)在下限的值的差。本文在給出振動速度和位移的初始值前提下,通過對辛普森公式的改進(jìn),得出遞推算法,可求出每一點(diǎn)的振動速度與位移,進(jìn)而確定振動速度與位移曲線。而對于這兩個初始值,工程上可以取零值。

(5)

式中F為f的積分函數(shù)。

對于式(4),可以看出辛普森積分計(jì)算了相鄰2個采樣區(qū)域的定積分,涉及了3個采樣點(diǎn)。由于它不是針對單一采樣區(qū)域,且只涉及整個積分區(qū)域中首尾兩點(diǎn)的積分值,所以相鄰的辛普森積分序列不可能覆蓋所有采樣點(diǎn)。為了得到每一點(diǎn)的積分值,提出了一種組合辛普森積分方法,即對采樣點(diǎn)中奇數(shù)序列采用一套辛普森積分計(jì)算公式,對偶數(shù)序列采用另一套辛普森積分計(jì)算公式,兩者組合得到所有點(diǎn)的積分值。

這里設(shè)加速度采樣值序列為a{(0),a(1),…,a(n-2),a(n-1)},其中n為偶數(shù)。由辛普森積分原理,將奇數(shù)序列a{(1),a(3),…,a(n-3),a(n-1)}組成一個積分序列,將偶數(shù)序列a{(0),a(2),…,a(n-4),a(n-2)}組成一個積分序列。這樣可以得到如下的積分計(jì)算關(guān)系:

(6)

對于振動速度初始值v(0)可以取為0,而對于v(1)可以采用梯形積分公式計(jì)算法求得,見式(7)。2部分積分計(jì)算組合,就可得到全部采樣點(diǎn)的振動速度值,從而可以得到波形曲線,進(jìn)行參數(shù)計(jì)算。

(7)

同理,對于振動位移的積分運(yùn)算,也可在振動速度去趨勢項(xiàng)后,使用該方法處理。

4 算法校準(zhǔn)驗(yàn)證

為了對上述算法進(jìn)行驗(yàn)證,可以采用一個確定函數(shù)信號,由于確定函數(shù)通過理論計(jì)算可以求得其積分函數(shù),將這個積分函數(shù)作為標(biāo)準(zhǔn)信號與辛普森算法得到的積分函數(shù)相比較進(jìn)行驗(yàn)證。由于周期信號都可以展開傅里葉級數(shù),這里選取確定函數(shù)為一正弦函數(shù)10 sin(200 π t),其積分函數(shù)如式(8),波形如圖2所示。

(8)

圖2 標(biāo)準(zhǔn)正弦函數(shù)和其積分函數(shù)

將標(biāo)準(zhǔn)正弦函數(shù)10sin(200 π t)增加隨時間變化趨勢項(xiàng),合成得到10sin(200 π t)+50t+2,采用第3節(jié)所述組合辛普森積分算法對其積分,并對積分結(jié)果進(jìn)行趨勢項(xiàng)擬合,獲得積分變化趨勢項(xiàng),如圖3所示。得到的趨勢項(xiàng)p(t)如式(9)所示。將積分結(jié)果進(jìn)行去趨勢項(xiàng)處理,如圖4所示,可以看出組合辛普森積分去趨勢項(xiàng)后的曲線與理論計(jì)算的積分曲線高度一致,說明了本方法的處理能力。

p(t)=24.98t2+2.003t+0.0158

(9)

圖3 組合辛普森積分及趨勢項(xiàng)擬合曲線

圖4 去趨勢項(xiàng)積分曲線

辛普森積分具有3次代數(shù)精度,對于小于等于3次的積分多項(xiàng)式可以精確積分,如果f∈C4(a,b),則存在截?cái)嗾`差項(xiàng)e(f):

(10)

對于組合辛普森積分,需要將截?cái)嗾`差進(jìn)行累積,設(shè)M=max|f(4)(x)|,則總的截?cái)嗾`差E(f)絕對值為:

(11)

式中f(4)(ξ)為函數(shù)f的4階導(dǎo)數(shù)。

從式(11)可以看出,總的截?cái)嗾`差E(f)為h4的高階無窮小,即o(h4),當(dāng)采樣步長h足夠小的情況下,總的截?cái)嗾`差也相當(dāng)小。例如對于上面的標(biāo)準(zhǔn)正弦信號,采樣步長為0.000 1,所得總的截?cái)嗾`差小于0.173 2×10-6,可以忽略。

對于振動測試來說,為了保證信號采集的完整性,采樣頻率通常很高,也即采樣周期很短,所以從計(jì)算精度和計(jì)算效率的角度來考慮,辛普森積分可以滿足振動測試的時域積分要求。

5 Matlab振動信號測試

Matlab中包含了一個可視化地震波信號quake,它是一個開放的數(shù)據(jù)源,記錄了1989年美國加利福尼亞州洛馬·普雷塔大地震中200 Hz采樣頻率的地震波信號,它通過加速度傳感器記錄了東-西、南-北和垂直3個方向的地震波信號[12],實(shí)驗(yàn)通過這些數(shù)據(jù)可以得到3個方向的振動速度和振動位移,圖5分別給出了采用本文所述的組合辛普森積分方法得到的3個方向8～15 s間的振動速度和振動位移波形,采樣間隔為5 ms,每個方向共 1 401 個采樣點(diǎn)。

圖5 地震波振動測試

實(shí)際上Matlab提供的演示程序quake,也對地震波信號進(jìn)行了積分處理,得到了3個方向的振動速度和振動位移信號,然而Matlab中采用的積分算法為cumsum,它實(shí)際上是矩形積分算法,就是以采樣間隔內(nèi)的矩形面積近似代替這一區(qū)域的積分值。圖6(a)和圖7(a)是Matlab中矩形積分算法得到的振動速度和振動位移波形曲線,圖6(b)和圖7(b)是本文論述的組合辛普森積分算法得到的振動速度和振動位移波形曲線。從圖6(a)和圖6(b)、圖7(a)和圖7(b)的圖形曲線對比中可以看出兩種曲線形狀基本一致,這也說明了本算法的有效性。

圖6 2種算法對于振動速度的積分對比

圖7 2種算法對于振動位移的積分對比

對于辛普森積分,依據(jù)式(11),其截?cái)嗾`差為:

24.31×10-12f(4)(c)

(12)

而矩形積分,其截?cái)嗾`差為:

(13)

從式(12)和式(13)可以看出,組合辛普森積分總的截?cái)嗾`差為o(h4),而矩形積分總的截?cái)嗾`差為o(h2)。所以就算法本身而言,辛普森積分具有更高的計(jì)算精度,更小的截止誤差。

5 結(jié) 論

本文根據(jù)振動測試要求,由辛普森積分公式和牛頓-萊布尼茲公式,將原始信號分成奇數(shù)序列和偶數(shù)序列,推導(dǎo)得出了組合辛普森積分算法。結(jié)合對積分初始值的有效處理,得出了相應(yīng)的遞推計(jì)算公式。應(yīng)用此算法可以得到測試數(shù)據(jù)序列的完整積分波形,對地震波數(shù)據(jù)的實(shí)驗(yàn),驗(yàn)證了本方法的有效性。并且該算法具有更高的計(jì)算精度,非常適于振動測試的時域積分應(yīng)用。