基于Welch功率譜的加速度積分改進(jìn)方法研究

李展銓， 陳太聰,2,3

(1. 華南理工大學(xué) 土木與交通學(xué)院，廣州 510641； 2. 亞熱帶建筑科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，廣州 510641； 3. 琶洲實(shí)驗(yàn)室，廣州 510335)

在工程振動實(shí)踐中，結(jié)構(gòu)的加速度響應(yīng)相對于位移和速度響應(yīng)更易測量，常用于結(jié)構(gòu)模態(tài)識別和損傷判別等分析[1-2]。而在結(jié)構(gòu)抗震性能評價、結(jié)構(gòu)阻尼評估等動力效應(yīng)評定中，位移和速度指標(biāo)更具有關(guān)鍵意義[3-5]。對加速度信號進(jìn)行時域積分獲得位移和速度是振動信號處理中的常用方法[6]，但由于儀器誤差和環(huán)境影響，測量得到的加速度信號往往帶有干擾噪聲，導(dǎo)致積分得到的速度和位移與真實(shí)響應(yīng)有較大差別，甚至完全失真。如何消除干擾噪聲以還原真實(shí)的響應(yīng)信號，一直是振動工程領(lǐng)域的重點(diǎn)關(guān)注問題。

不同于傳統(tǒng)的時域積分，徐慶華[7]嘗試頻域積分方法，通過快速傅立葉變換(fast Fourier transform，F(xiàn)FT)實(shí)現(xiàn)加速度、速度與位移間的相互轉(zhuǎn)換。Brandt等[8]通過對不同積分方法的對比研究，推薦帶低頻截止的頻域積分方法，參考加速度信號的第一階主頻，將低頻部分進(jìn)行歸零處理。方新磊等[9]在低頻截止的基礎(chǔ)上引入高頻截止濾波，以同時消除高頻干擾噪聲和低頻趨勢項(xiàng)。以上頻率截止法對截止頻率參數(shù)較為敏感，而相關(guān)參數(shù)的設(shè)置有很大的主觀性。

陳太聰?shù)萚10]針對頻率截止法的局限性提出了有效頻段法，基于加速度信號的FFT譜曲線形態(tài)擬合，進(jìn)行加速度信號的濾噪積分。該法對于中等噪聲干擾的加速度信號的積分結(jié)果精度較高，有效避免了頻率截止法的參數(shù)敏感性和主觀性等問題。但在高噪聲干擾下，積分精度隨著噪聲水平的增大而迅速變差。其次，該法的實(shí)施也需要人為給定主頻近似信息以判定分析頻率范圍，不利于工程應(yīng)用。

事實(shí)上，在信號處理領(lǐng)域，除了FFT譜曲線可以反映信號的頻域信息，功率譜曲線也具有類似的功能。傳統(tǒng)的周期圖法直接通過FFT譜的平方求解功率譜曲線，計算方便，但譜線起伏大，弱信號下的頻譜分辨率低[11]；Welch算法引入數(shù)據(jù)分段交疊和加窗函數(shù)的方法，對周期圖法進(jìn)行了改進(jìn)，得到修正功率譜，有效降低噪聲影響，所得譜線形狀平滑，易于清晰分辨主頻。如王福杰等[12]通過試驗(yàn)對比各種功率譜后發(fā)現(xiàn)，應(yīng)用Welch功率譜估計，從寬帶噪聲中檢測出窄帶信號的效果最好。如今Welch功率譜已廣泛應(yīng)用于如脈動風(fēng)譜估計[13]、鐵路軌道不平順度測試[14]、地震監(jiān)測背景噪聲分析[15]等研究和工程實(shí)踐中。

針對上述有效頻段法的若干應(yīng)用問題，本文提出引入Welch功率譜的改進(jìn)有效頻段法，用于加速度信號濾噪積分的自動處理?；趲г爰铀俣刃盘柕腤elch功率譜，自動實(shí)現(xiàn)主頻的識別和分析頻率范圍的定義。繼而分別根據(jù)Welch功率譜曲線和Welch功率譜開方曲線進(jìn)行形態(tài)擬合，實(shí)現(xiàn)有效頻段的自動識別。最后進(jìn)行有效頻段內(nèi)的頻域積分，得到速度與位移信號。文后通過多頻激勵和隨機(jī)激勵下的數(shù)值模擬算例，開展原有效頻段法和本文改進(jìn)方法的對比分析，以檢驗(yàn)方法在不同激勵下的積分精度和抗噪性能。

1 有效頻段法的基本原理

本文方法是有效頻段法的改進(jìn)研究，在此先介紹有效頻段法的基本計算原理。

1.1 分段頻譜曲線形態(tài)假設(shè)

假設(shè)按等時距Δt測量得到的帶噪加速度信號為X(n),n=0,1,2,…,N-1，則根據(jù)離散傅里葉變化的基本原理，離散的時域加速度可轉(zhuǎn)換為頻域內(nèi)若干離散的FFT頻譜信號H(k),k=0,1,2,…,N-1，其相互轉(zhuǎn)化關(guān)系為[16]

(1)

(2)

式中：N為采樣數(shù)據(jù)量；n為時刻點(diǎn)；k為譜線序列點(diǎn)。

設(shè)加速度信號含m個峰值主頻，按從小到大的順序依次為f1,f2,…,fm，則整個頻率范圍按式(3)方式分為用于分析的m段

(3)

式中，fmax為頻譜曲線對應(yīng)的最大頻率。

對m個頻段均進(jìn)行歸一化處理，即都?xì)w一化為范圍(0,1]。假設(shè)每個頻段內(nèi)的頻譜曲線都符合相應(yīng)的高斯函數(shù)分布[17]

Gi(x)=aie-(x-bi)2/(2c2)

(4)

式中，i=1,2,…,m；x∈(0,1]為每個頻段內(nèi)的歸一化頻率；參數(shù)ai，bi和ci分別為高斯函數(shù)的幅值、中心位置和半徑。

基于式(4)定義的頻譜分布形態(tài)，相應(yīng)的每一頻段內(nèi)的頻譜累計能量分布函數(shù)可由式(5)計算得到

(5)

式中： erf(·)為誤差函數(shù)；Ei(x)在形態(tài)上表現(xiàn)為反Z型的單調(diào)遞增函數(shù)。

對累計能量分布函數(shù)進(jìn)一步進(jìn)行歸一化處理，得

(6)

1.2 分段累積信號能量計算

在第i個歸一化頻段內(nèi)，離散的加速度頻譜信號的累積能量由式(7)計算得到

(7)

對信號累計能量函數(shù)進(jìn)一步做歸一化處理，得

(8)

1.3 參數(shù)擬合與有效頻段確定

基于式(8)計算得到的離散數(shù)據(jù)，用式(6)定義的非線性連續(xù)函數(shù)進(jìn)行非線性擬合，迭代收斂后得到參數(shù)bi和ci。

根據(jù)概率統(tǒng)計理論的三倍標(biāo)準(zhǔn)差原則[18]，第i階主頻能量分布的歸一化有效頻段取為[bi-3ci,bi+3ci]，相應(yīng)的實(shí)際有效頻段為

(9)

最后經(jīng)過有效頻段截斷后得到的加速度、速度和位移分別為[19]

(10)

式中，ωk為加速度信號的頻率，ωk=2πkFs/N。

(11)

有效頻段法用于加速度的濾噪和積分的相關(guān)計算流程，如圖1所示。

圖1 有效頻段法的積分流程圖Fig.1 Integration flowchart of frequency truncation method

2 基于Welch功率譜的有效頻段法

已有研究表明，有效頻段法在處理低于25%噪聲水平的加速度信號積分時有較好的效果，但在更高噪聲干擾下，由于主頻處的FFT譜混入了較大的噪聲，噪聲能量造成了擬合高斯函數(shù)的半徑變大，導(dǎo)致有效頻段的范圍擴(kuò)大，從而引入更多噪聲影響，積分精度難以保證。此外，有效頻段法需要人為指定主頻的近似取值，不能實(shí)現(xiàn)加速度積分的全程自動分析，也不便于實(shí)際工程應(yīng)用。

針對有效頻段法的以上問題，本文引入Welch功率譜對原方法進(jìn)行改進(jìn)，以實(shí)現(xiàn)自動化和高耐噪的加速度積分分析。一方面，Welch譜曲線的形態(tài)平滑，峰值清晰，易于自動識別主頻；另一方面，它作為一種平均功率譜，本身具有良好的抑噪功能，有利于消除高噪聲對擬合高斯函數(shù)參數(shù)的影響。

如圖2為文后算例中多頻激勵情況下的帶噪和不帶噪加速度信號的FFT譜曲線和Welch譜曲線對比，由圖2可知：①與FFT譜相比，Welch譜在信號主頻處的主瓣較寬、旁瓣較低、起伏性小、曲線平滑，在形態(tài)上更加符合高斯曲線的分布；②與FFT譜不同，帶噪Welch譜與不帶噪Welch譜的差別不大，且頻率影響范圍基本一致，表明Welch函數(shù)處理具有一定的消噪功能。

圖2 35%噪聲水平信號時程曲線以及其Welch功率譜圖與FFT頻譜圖Fig.2 Time-history curve, Welch power spectra and FFT spectra of 35% noise level signals

2.1 基于Welch譜的自動選峰和頻域分段

設(shè)實(shí)測加速度信號X(n)，對應(yīng)FFT頻譜為H(k)，對應(yīng)Welch功率譜為P(k)。取P(k)max×5%作為識別主頻峰值的下限，得到m階主頻，按從小到大的順序依次為f1,f2,…,fm，并選擇各個主頻間的最低波谷頻率作為頻域分段的依據(jù)，則有

(12)

2.2 分段頻譜曲線形態(tài)擬合假設(shè)

根據(jù)不同類型激勵下的結(jié)構(gòu)振動響應(yīng)特點(diǎn)，本文考察3種不同類型的加速度頻譜曲線：①FFT頻譜曲線，這也是原有效頻段法所采用的頻譜曲線；②Welch功率譜曲線；③Welch功率譜開方曲線，由于功率譜在數(shù)值上是幅值的平方關(guān)系，因此功率譜開方在數(shù)值上與FFT譜具有對應(yīng)關(guān)系。

分別假設(shè)這3種頻譜曲線符合相應(yīng)的高斯曲線分布，即可采用如式(4)所示的帶不同參數(shù)的高斯函數(shù)進(jìn)行描述。相應(yīng)地，擬合高斯函數(shù)的累計能量分布函數(shù)及其歸一化函數(shù)分別由式(5)和式(6)計算。

2.3 分段累積信號能量計算

對于實(shí)測的帶噪加速度信號，在第i個歸一化頻段內(nèi)，離散的加速度頻譜信號的累積能量，根據(jù)3種不同類型的頻譜曲線，分別按式(13)計算得到

(13)

其相應(yīng)的歸一化函數(shù)分別根據(jù)式(8)進(jìn)行計算。

2.4 參數(shù)擬合與有效頻段確定

基于式(13)和式(8)計算得到的離散數(shù)據(jù)，用式(6)定義的非線性連續(xù)函數(shù)進(jìn)行非線性擬合，迭代收斂后得到參數(shù)bi和ci。繼而由式(9)確定各有效頻段。最后采用式(10)積分計算得到濾噪后的加速度、速度和位移信號。

最終，改進(jìn)方法用于加速度的濾噪和積分的相關(guān)計算流程，如圖3所示。由圖3可知，本文改進(jìn)方法的實(shí)現(xiàn)過程與有效頻段法大致相同，不同的是利用了Welch功率譜自動進(jìn)行頻率分段，并針對基于Welch譜的兩種譜曲線進(jìn)行高斯函數(shù)形態(tài)擬合，擬合得到的高斯函數(shù)參數(shù)不會因噪聲水平而產(chǎn)生大的偏差，保證有效頻段的精度，避免高噪聲干擾下的積分效果失真情況。

圖3 基于Welch功率譜形態(tài)擬合的積分流程圖Fig.3 Integration flowchart based on morphology fitting of Welch power spectrum

本文方法在MATLAB軟件中進(jìn)行了實(shí)現(xiàn)，相關(guān)技術(shù)開發(fā)應(yīng)用了軟件內(nèi)嵌的傅里葉變換函數(shù)、傅里葉逆變換函數(shù)(inverse fast Fourier transform，IFFT)、Welch功率譜生成函數(shù)和尋峰函數(shù)等若干函數(shù)命令。

3 數(shù)值算例

為了驗(yàn)證方法效果，本文采用與原有效頻段法一致的數(shù)值算例進(jìn)行對比分析。

六自由度懸臂梁結(jié)構(gòu)，考慮豎向振動，如圖4所示。其各階模態(tài)阻尼比為0.01，質(zhì)量陣和剛度陣為

圖4 六自由度懸臂梁結(jié)構(gòu)模型Fig.4 Structural model of a 6-dofs cantilever beam

在結(jié)構(gòu)自由端分別施加豎向多頻激勵和白噪聲隨機(jī)激勵，其中多頻激勵取為

z(t)=200sin 6πt+100sin 20πt+100sin 40πt

然后按照Newmark-β法(γ=0.5，β=0.25)求出加速度、速度和位移響應(yīng)時程，作為精確參考解。在加速度解中分別加入1%，5%，10%，15%，20%，25%，30%，35%和40%的白噪聲干擾作為觀測加速度。再分別采用原有效頻段法和本文改進(jìn)方法得到修正的加速度、速度和位移信號。其中，改進(jìn)方法采用兩種頻譜曲線進(jìn)行高斯函數(shù)形態(tài)擬合。最后采用如下定義的總體誤差指標(biāo)來評價積分精度

(14)

3.1 多頻激勵情況

在多頻激勵下，以質(zhì)點(diǎn)3處的豎向響應(yīng)為例，不同噪聲水平下兩種方法的結(jié)果精度對比，如圖5所示。35%噪聲水平下兩種方法濾噪積分后所得的位移時程對比，如圖6所示。

圖5 不同噪聲水平下的結(jié)果誤差Fig.5 Result errors under different noise levels

圖6 35%噪聲下的位移時程Fig.6 Time history of displacement under 35% noise level

由圖5～圖6可知，在多頻激勵作用下：

(1) 基于Welch功率譜曲線擬合的改進(jìn)方法計算得到的位移結(jié)果，在40%噪聲水平以下均比原方法的精度更高，具有良好的抗噪穩(wěn)定性，在高噪聲情況下仍能保持結(jié)果誤差低于3%。

(2) 基于Welch功率譜開方曲線擬合的改進(jìn)方法計算得到的位移結(jié)果，在15%噪聲水平以下比原方法的精度更高，在15%～25%噪聲水平間與原方法相近，在25%噪聲水平以上具有一定的抗噪穩(wěn)定性，沒有出現(xiàn)結(jié)果失真。

綜合而言，在多頻激勵的情況下，基于兩種頻譜曲線的改進(jìn)方法均比原有效頻段法的積分精度更高，具有更好的抗噪穩(wěn)定性。其中，基于Welch功率譜曲線擬合的改進(jìn)效果最好。其原因是Welch功率譜曲線在數(shù)值上可表示為FFT頻譜曲線和Welch功率譜開方曲線的平方關(guān)系，從而放大各主頻處的真實(shí)信號與旁瓣噪聲信號之間的差異，減少了噪聲信號對形態(tài)擬合參數(shù)的影響，避免了高噪聲下有效頻段過寬的問題。因此，針對多頻激勵情況，本文建議使用基于Welch功率譜曲線擬合的改進(jìn)方法。

3.2 隨機(jī)激勵情況

在白噪聲隨機(jī)激勵下，同樣以質(zhì)點(diǎn)3處的豎向響應(yīng)為例。不同噪聲水平下兩種方法的結(jié)果精度對比，如圖7所示。35%噪聲水平下兩種方法濾噪積分后所得的位移時程對比，如圖8所示。

圖7 不同噪聲水平下的結(jié)果誤差Fig.7 Result errors under different noise levels

由圖7～圖8可知，在白噪聲隨機(jī)激勵作用下：

圖8 35%噪聲下的位移時程Fig.8 Time history of displacement under 35% noise level

(1) 基于Welch功率譜曲線擬合的改進(jìn)方法計算得到的位移結(jié)果，誤差穩(wěn)定于10%左右，在低噪聲水平時略差于原方法，在高噪聲水平時具有一定的抗噪穩(wěn)定性，沒有出現(xiàn)結(jié)果失真。

(2) 基于Welch功率譜開方曲線擬合的改進(jìn)方法計算得到的位移結(jié)果，在40%噪聲水平以下均比原方法的精度更高，具有良好的抗噪穩(wěn)定性，在高噪聲情況下仍能保持結(jié)果誤差低于10%。

綜合而言，在隨機(jī)激勵的情況下，基于Welch功率譜開方曲線擬合的改進(jìn)方法比基于Welch功率譜曲線擬合的改進(jìn)方法的效果更好，也比原有效頻段法的積分精度更高且抗噪穩(wěn)定性更強(qiáng)。其原因是白噪聲作用下的響應(yīng)信號屬于旁瓣較大的信號，Welch功率譜開方曲線相對于Welch功率譜曲線，相對弱化了信號主瓣和旁瓣之間的差異，適當(dāng)保留了旁瓣頻率的貢獻(xiàn)，避免了高噪聲下有效頻段過窄的問題。因此，針對隨機(jī)激勵情況，本文建議使用基于Welch功率譜開方曲線擬合的改進(jìn)方法。

4 結(jié) 論

針對工程振動測試中如何由帶噪加速度信號還原出真實(shí)速度和位移的問題，新近發(fā)展的有效頻段法通過對信號FFT頻譜曲線的形態(tài)擬合，實(shí)現(xiàn)了中等噪聲干擾下的高精度加速度積分。但方法仍存在高噪聲干擾下積分結(jié)果容易失真，以及需人為給定主頻近似信息等問題，不利于工程應(yīng)用。

本文引入Welch功率譜對有效頻段法進(jìn)行改進(jìn)，以實(shí)現(xiàn)自動化和高耐噪的加速度積分分析。基于帶噪加速度信號的Welch功率譜，自動實(shí)現(xiàn)主頻的識別和分析頻率范圍的定義。繼而分別根據(jù)Welch功率譜曲線和Welch功率譜開方曲線進(jìn)行形態(tài)擬合，實(shí)現(xiàn)有效頻段的自動識別。其中，針對多頻激勵情況，推薦使用基于Welch功率譜曲線擬合的改進(jìn)方法；針對隨機(jī)激勵情況，推薦使用基于Welch功率譜開方曲線擬合的改進(jìn)方法。

還需說明的是，以上方法僅適用于穩(wěn)態(tài)加速度信號的積分處理。針對非穩(wěn)態(tài)加速度信號，本文作者基于有效頻段法，綜合應(yīng)用Welch功率譜和短時傅里葉變換，獲得了良好的積分效果，后續(xù)將另文報告相關(guān)研究成果。