許曉紅, 倪春鋒, 顏曉青

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

1 前 言

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

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

2 測(cè)量原理

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

(1)

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

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

圖1 測(cè)量原理框圖

3 組合辛普森積分方法

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

(2)

(3)

(4)

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

(5)

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

對(duì)于式(4),可以看出辛普森積分計(jì)算了相鄰2個(gè)采樣區(qū)域的定積分,涉及了3個(gè)采樣點(diǎn)。由于它不是針對(duì)單一采樣區(qū)域,且只涉及整個(gè)積分區(qū)域中首尾兩點(diǎn)的積分值,所以相鄰的辛普森積分序列不可能覆蓋所有采樣點(diǎn)。為了得到每一點(diǎn)的積分值,提出了一種組合辛普森積分方法,即對(duì)采樣點(diǎn)中奇數(shù)序列采用一套辛普森積分計(jì)算公式,對(duì)偶數(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)}組成一個(gè)積分序列,將偶數(shù)序列a{(0),a(2),…,a(n-4),a(n-2)}組成一個(gè)積分序列。這樣可以得到如下的積分計(jì)算關(guān)系:

(6)

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

(7)

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

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

為了對(duì)上述算法進(jìn)行驗(yàn)證,可以采用一個(gè)確定函數(shù)信號(hào),由于確定函數(shù)通過(guò)理論計(jì)算可以求得其積分函數(shù),將這個(gè)積分函數(shù)作為標(biāo)準(zhǔn)信號(hào)與辛普森算法得到的積分函數(shù)相比較進(jìn)行驗(yàn)證。由于周期信號(hào)都可以展開(kāi)傅里葉級(jí)數(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)增加隨時(shí)間變化趨勢(shì)項(xiàng),合成得到10sin(200 π t)+50t+2,采用第3節(jié)所述組合辛普森積分算法對(duì)其積分,并對(duì)積分結(jié)果進(jìn)行趨勢(shì)項(xiàng)擬合,獲得積分變化趨勢(shì)項(xiàng),如圖3所示。得到的趨勢(shì)項(xiàng)p(t)如式(9)所示。將積分結(jié)果進(jìn)行去趨勢(shì)項(xiàng)處理,如圖4所示,可以看出組合辛普森積分去趨勢(shì)項(xiàng)后的曲線(xiàn)與理論計(jì)算的積分曲線(xiàn)高度一致,說(shuō)明了本方法的處理能力。

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

(9)

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

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

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

(10)

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

(11)

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

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

對(duì)于振動(dòng)測(cè)試來(lái)說(shuō),為了保證信號(hào)采集的完整性,采樣頻率通常很高,也即采樣周期很短,所以從計(jì)算精度和計(jì)算效率的角度來(lái)考慮,辛普森積分可以滿(mǎn)足振動(dòng)測(cè)試的時(shí)域積分要求。

5 Matlab振動(dòng)信號(hào)測(cè)試

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

圖5 地震波振動(dòng)測(cè)試

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

圖6 2種算法對(duì)于振動(dòng)速度的積分對(duì)比

圖7 2種算法對(duì)于振動(dòng)位移的積分對(duì)比

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