基于Welch功率譜的加速度積分改進方法研究

李展銓， 陳太聰,2,3

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

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

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

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

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

針對上述有效頻段法的若干應用問題，本文提出引入Welch功率譜的改進有效頻段法，用于加速度信號濾噪積分的自動處理。基于帶噪加速度信號的Welch功率譜，自動實現主頻的識別和分析頻率范圍的定義。繼而分別根據Welch功率譜曲線和Welch功率譜開方曲線進行形態擬合，實現有效頻段的自動識別。最后進行有效頻段內的頻域積分，得到速度與位移信號。文后通過多頻激勵和隨機激勵下的數值模擬算例，開展原有效頻段法和本文改進方法的對比分析，以檢驗方法在不同激勵下的積分精度和抗噪性能。

1 有效頻段法的基本原理

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

1.1 分段頻譜曲線形態假設

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

(1)

(2)

式中：N為采樣數據量；n為時刻點；k為譜線序列點。

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

(3)

式中，fmax為頻譜曲線對應的最大頻率。

對m個頻段均進行歸一化處理，即都歸一化為范圍(0,1]。假設每個頻段內的頻譜曲線都符合相應的高斯函數分布[17]

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

(4)

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

基于式(4)定義的頻譜分布形態，相應的每一頻段內的頻譜累計能量分布函數可由式(5)計算得到

(5)

式中： erf(·)為誤差函數；Ei(x)在形態上表現為反Z型的單調遞增函數。

對累計能量分布函數進一步進行歸一化處理，得

(6)

1.2 分段累積信號能量計算

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

(7)

對信號累計能量函數進一步做歸一化處理，得

(8)

1.3 參數擬合與有效頻段確定

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

根據概率統計理論的三倍標準差原則[18]，第i階主頻能量分布的歸一化有效頻段取為[bi-3ci,bi+3ci]，相應的實際有效頻段為

(9)

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

(10)

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

(11)

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

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

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

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

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

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

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

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

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

(12)

2.2 分段頻譜曲線形態擬合假設

根據不同類型激勵下的結構振動響應特點，本文考察3種不同類型的加速度頻譜曲線：①FFT頻譜曲線，這也是原有效頻段法所采用的頻譜曲線；②Welch功率譜曲線；③Welch功率譜開方曲線，由于功率譜在數值上是幅值的平方關系，因此功率譜開方在數值上與FFT譜具有對應關系。

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

2.3 分段累積信號能量計算

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

(13)

其相應的歸一化函數分別根據式(8)進行計算。

2.4 參數擬合與有效頻段確定

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

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

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

本文方法在MATLAB軟件中進行了實現，相關技術開發應用了軟件內嵌的傅里葉變換函數、傅里葉逆變換函數(inverse fast Fourier transform，IFFT)、Welch功率譜生成函數和尋峰函數等若干函數命令。

3 數值算例

為了驗證方法效果，本文采用與原有效頻段法一致的數值算例進行對比分析。

六自由度懸臂梁結構，考慮豎向振動，如圖4所示。其各階模態阻尼比為0.01，質量陣和剛度陣為

圖4 六自由度懸臂梁結構模型Fig.4 Structural model of a 6-dofs cantilever beam

在結構自由端分別施加豎向多頻激勵和白噪聲隨機激勵，其中多頻激勵取為

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

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

(14)

3.1 多頻激勵情況

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

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

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

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

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

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

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

3.2 隨機激勵情況

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

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

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

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

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

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

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

4 結 論

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

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

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