任 程，劉世洲，郜 偉，張二亮

(鄭州大學(xué) 機械與動力工程學(xué)院， 鄭州 450001)

0 引言

頻率響應(yīng)函數(shù)(frequency response function,FRF)是機械系統(tǒng)動態(tài)特性評估、振動控制和故障診斷等的有力工具。使用力錘對結(jié)構(gòu)施加沖擊載荷，是測量頻響函數(shù)的重要方式，具有操作簡單、成本低和效率高等優(yōu)點，在機械、土木和航空航天等工程領(lǐng)域得到了廣泛應(yīng)用[1]。因此，開展沖擊載荷下頻響函數(shù)的高精度辨識方法研究具有鮮明的工程應(yīng)用背景和價值。

基于力錘沖擊的FRF測量及其應(yīng)用已經(jīng)在諸多學(xué)術(shù)專著中予以詳細介紹。影響沖擊載荷下FRF測量的關(guān)鍵因素包括力錘沖擊頭材質(zhì)選取、實驗方案設(shè)計(例如預(yù)觸發(fā)確定)以及信號采集等[2]。Bediz等[3]設(shè)計了一種用于微型結(jié)構(gòu)模態(tài)實驗的沖擊激勵系統(tǒng)，具有可重復(fù)、高帶寬和沖擊力可控等優(yōu)點。Xia等[4]開展了玻璃鋼拉索的沖擊振動實驗，采用半功率點法計算實驗?zāi)B(tài)阻尼。針對大型結(jié)構(gòu)的動態(tài)特性測量，Tian等[5]提出了基于非接觸視覺測量的移動沖擊測試方法，以識別整個結(jié)構(gòu)的模態(tài)振型和柔度矩陣。

譜估計方法直接從測量的輸入輸出數(shù)據(jù)中計算FRF，是錘擊法測試中開展FRF估計的主要方法。從算法角度考慮，窗函數(shù)設(shè)計是影響FRF譜估計精度的主要因素?；谧V估計方法的計算過程，Yang等[6]建立了插值和瞬態(tài)誤差的非參數(shù)模型，分析了窗函數(shù)對兩種誤差的影響。針對譜估計窗函數(shù)帶來的插值和瞬態(tài)誤差，Shao等[7]提出了一種迭代補償方法，用于提高沖擊載荷下頻響函數(shù)的辨識精度。然而，譜估計方法完全依賴于測量數(shù)據(jù)，F(xiàn)RF的估計精度極易受到測量數(shù)據(jù)長度和質(zhì)量的影響。

高斯過程回歸(gaussian process regression,GPR)使用無限維多變量高斯分布對函數(shù)進行直接建模，已經(jīng)發(fā)展成為一種根據(jù)經(jīng)驗樣本學(xué)習(xí)映射關(guān)系的非參數(shù)高效數(shù)據(jù)驅(qū)動建模方法[8]。相對于經(jīng)典的譜估計方法，GPR不僅利用了樣本數(shù)據(jù)，還能融合函數(shù)模型的固有性質(zhì)；相比于傳統(tǒng)的參數(shù)化建模方法，GPR在模型選擇方面具有更多靈活性，能夠更好地適應(yīng)數(shù)據(jù)的變化。高斯過程的性質(zhì)與其協(xié)方差函數(shù)(即核函數(shù))有著密切聯(lián)系。通過刻畫系統(tǒng)脈沖響應(yīng)函數(shù)的固有性質(zhì)(例如光滑、穩(wěn)定)，基于最大熵原理，穩(wěn)定樣條、對角相關(guān)等核函數(shù)被相繼提出和發(fā)展，在系統(tǒng)與控制領(lǐng)域內(nèi)發(fā)展了傳遞函數(shù)的正則化估計方法[9-12]。

本文將系統(tǒng)和控制領(lǐng)域的正則化方法引入實驗?zāi)B(tài)分析領(lǐng)域，較為系統(tǒng)地給出了沖擊載荷下FRF估計的GPR理論與方法。針對機械系統(tǒng)階次高及低階阻尼小的特點，豐富了FRF的復(fù)高斯過程先驗?zāi)Ｐ汀；谪惾~斯推斷框架，利用GPR模型給出了FRF的最大后驗估計及其方差，并借助QR分解改善了復(fù)高斯過程超參數(shù)優(yōu)化的數(shù)值正定性。最后，通過開展平板沖擊仿真和葉片振動實驗，驗證了本文方法的有效性和可靠性。

1 問題描述

在實際的工程測試中，采集的信號都是時域的離散信號，采用離散傅里葉變換(discrete fourier transform,DFT)將時域信號x[t],t=0,1,…,N-1變換到頻域X(k)，

(1)

式中：k=0,1,…,N/2， j2=-1。

對于線性時不變動態(tài)系統(tǒng)，沖擊激勵u[t]和系統(tǒng)響應(yīng)y[t]的DFT頻譜關(guān)系如圖1所示，可表示為

Y(k)=G(ωk)U(k)+V(k)

(2)

圖1 DFT頻譜關(guān)系

當響應(yīng)信號采集不完整時，在關(guān)系中需要引入額外的瞬態(tài)項，用于描述非周期信號傅里葉變換帶來的泄露誤差[13]。若激勵信號具備粗糙性質(zhì)，該瞬態(tài)項可以通過連同F(xiàn)RF一起辨識的方式來剔除。然而，沖擊激勵的頻譜是光滑的，無法分離瞬態(tài)項和FRF，則需要通過施加指數(shù)衰減窗函數(shù)來抑制泄露誤差。因而，沖擊載荷下的數(shù)據(jù)采集樣本長度應(yīng)盡可能足夠長。

從測量數(shù)據(jù)Z={(U(k),Y(k)),k=1,2,…,F}(直流分量和奈奎斯特頻率除外)中，使用GPR模型對FRF進行高精度建模和估計，是本文要解決的問題。

2 非參數(shù)辨識方法

2.1 復(fù)高斯過程

基于貝葉斯觀點，本文將FRF視作頻域上的復(fù)值隨機函數(shù)。進一步，假設(shè)FRF是頻域上的復(fù)高斯向量，即G服從多元復(fù)高斯分布。

一般來說，機械結(jié)構(gòu)的FRF在所關(guān)心的頻帶內(nèi)具有模型復(fù)雜度較高、低階阻尼較小得的特點，導(dǎo)致其建模難度大，這要求復(fù)高斯過程具備較強的泛化能力。為此，不同于常用的零均值假設(shè)[9-12]，本文通過引入FRF的先驗均值函數(shù)，使復(fù)高斯過程充分融合系統(tǒng)的先驗信息。復(fù)高斯向量G可表示為

G=Gp+ΔG

(3)

式中：Gp為FRF的先驗均值函數(shù)，ΔG服從零均值多元復(fù)高斯分布。

先驗均值函數(shù)Gp可以借助模態(tài)參數(shù)識別方法得到。例如，采用PolyMAX方法從測量數(shù)據(jù)Z中確定動態(tài)系統(tǒng)的某些極點pr，再將pr應(yīng)用于FRF的極點-余項模型，確定模態(tài)常數(shù)Ar和殘余項AU、AL，則

(4)

式中，上標*表示復(fù)共軛。

(5)

(6)

(7)

需要注意的是，核函數(shù)K和關(guān)系函數(shù)C滿足K(ωk,ωl)=C(ωk,ω-l)。

高斯過程的核函數(shù)刻畫了待建模系統(tǒng)的物理信息，在回歸分析中起著決定性的作用。對角相關(guān)核函數(shù)是在系統(tǒng)和控制領(lǐng)域使用最為廣泛的核函數(shù)之一，描述了脈沖響應(yīng)的光滑、穩(wěn)定等基本性質(zhì)。因此，本文采用對角相關(guān)核函數(shù)對FRF開展復(fù)高斯過程建模，其頻域表達式為

(8)

式中：比例因子γ∈R+和核參數(shù)α,β∈R+稱為超參數(shù)。

2.2 似然函數(shù)

根據(jù)中心極限定理和輸出噪聲的獨立同分布假設(shè)，當數(shù)據(jù)量F→∞時，干擾噪聲的DFT頻譜是圓形復(fù)高斯分布，其關(guān)系矩陣為零矩陣。為了與式(5)一致，F(xiàn)RF的似然函數(shù)寫成如下形式

(9)

2.3 最大后驗估計

(10)

式中：∝表示正比于，概率密度函數(shù)p(θ)通常假設(shè)為均勻分布。

最大后驗概率估計常用于參數(shù)的推演。聯(lián)合式(5)(9)和(10)，F(xiàn)RF的后驗分布仍為復(fù)高斯分布，其最大后驗概率估計為

(11)

(12)

2.4 超參數(shù)優(yōu)化

超參數(shù)θ是未知的，需要從測量數(shù)據(jù)Z中學(xué)習(xí)獲得，通常采用極大邊際似然估計方法。邊際似然函數(shù)可視作似然函數(shù)在先驗分布上的期望分布，它能夠自動權(quán)衡數(shù)據(jù)擬合精度與模型復(fù)雜程度。超參數(shù)θ的極大似然估計為

(13)

(14)

式中：L是下三角矩陣。根據(jù)Woodbury矩陣求逆公式和Sylvester行列式定理可得

(16)

受到式(15)(16)啟發(fā)，構(gòu)造如下矩陣，并將其改寫成QR分解式，

(17)

式中：Q為正交矩陣，R1為2F維上三角方陣，R2為2F維列向量，r為標量?；谑?17)可以得到

(18)

(19)

(20)

最后，目標函數(shù)J(θ)可以化簡為

(21)

采用基于梯度的優(yōu)化算法對目標函數(shù)進行最小化，目標函數(shù)J(θ)關(guān)于超參數(shù)θ的偏導(dǎo)數(shù)為

(22)

(23)

(24)

式中，Tr表示矩陣求跡運算。

基于式(21)和(22)，使用Matlab的優(yōu)化工具箱中fmincon函數(shù)進行尋優(yōu)計算，便可獲得最優(yōu)的超參數(shù)估計。該求解器可用于求解約束非線性多元函數(shù)最小值，包含內(nèi)點優(yōu)化、SQP優(yōu)化和信賴域反射優(yōu)化等不同類型的算法。

最后，根據(jù)超參數(shù)的最優(yōu)估計值，結(jié)合式(11)和(12)即可獲得FRF的估計及其方差，完整的算法流程如圖2所示。

3 算例與實驗

3.1 數(shù)值算例

為了驗證本文方法的有效性和可靠性，開展沖擊載荷下簡支平板的FRF估計。平板的彈性模量為200 GPa，尺寸為1 m×1 m×0.002 m，密度為7 850 kg/m3，泊松比為0.3。在中心位置施加沖擊載荷，采用模態(tài)疊加法生成系統(tǒng)響應(yīng)(本算例以速度作為響應(yīng))，添加高斯白噪聲來模擬測量數(shù)據(jù)，如圖3所示。使用的采樣頻率為1 024 Hz，感興趣的頻帶取1～300 Hz，數(shù)據(jù)樣本量F=300。

圖2 沖擊載荷下FRF的GPR方法流程框圖

圖3 平板1/4有限元模型 (a)和無噪聲響應(yīng)波形(b)

為了提高仿真運算結(jié)果的可靠度，以及對仿真得到的結(jié)果進行評價，可進行多次的蒙特卡洛模擬?；诿商乜迥M，計算均方誤差(mean-square error,MSE)對FRF的估計結(jié)果進行評價，

(25)

為了顯示更寬的數(shù)值范圍，本文使用分貝(dB)作為FRF幅值和MSE的單位。dB是振動、聲學(xué)和電信等研究領(lǐng)域常用的無量綱量，可定義為2個數(shù)值(測量值X和參考值X0)的對數(shù)比率，

1 dB=20lg|X/X0|

(26)

式中，默認參考值X0=1。

考慮響應(yīng)信號的信噪比(signal-to-noise ratio,SNR)為10 dB的情形，分別使用譜估計方法和本文方法對平板的FRF進行估計，結(jié)果如圖4所示。與譜估計方法相比，由于借助復(fù)高斯過程對FRF的固有性質(zhì)(光滑、穩(wěn)定)進行刻畫和融合，本文方法可以獲得更小的MSE值，F(xiàn)RF的估計精度更高，證明了其有效性。

圖4 平板的FRF估計曲線

令式(25)在所有關(guān)心頻率上取平均，進一步考慮不同SNR情形，MSE均值的計算結(jié)果如圖5所示。相比于譜估計方法，本文方法在不同的噪聲水平下均能夠獲得更高精度的FRF估計結(jié)果，且在低SNR(例如-10 dB)的嚴苛環(huán)境下效果更為顯著，證明了其可靠性。

圖5 不同SNR下的MSE均值

3.2 振動實驗

以長約44 cm的無人機螺旋槳葉片為例，開展錘擊振動實驗，驗證本文方法的實際應(yīng)用效果。使用PCB 086C03模態(tài)力錘激振葉片，利用數(shù)字圖像相關(guān)VIC-3D高速攝像系統(tǒng)測量葉片的振動響應(yīng)，如圖6所示。高速攝像機型號為Photron FastCam Mini UX100，配備了Nikon 60 mm f/2.8 D廣角鏡頭，板載內(nèi)存容量為8 GB，分辨率為1 024像素×576像素，曝光時間為0.25 ms，并以2 000 fps的幀率進行拍照記錄。VIC-3D高速攝像系統(tǒng)的位移測量精度約為6 μm，本實驗直接采用位移作為系統(tǒng)響應(yīng)。

圖6 葉片振動實驗

為了進一步抑制噪聲干擾和減少頻譜泄露，本實驗對時域信號進行加窗處理：沖擊信號施加單位增益窗，響應(yīng)信號施加指數(shù)窗。經(jīng)過加窗處理的時域信號波形如圖7所示。

圖7 加窗后時域信號波形

選取感興趣的頻帶為1～400 Hz，F(xiàn)=400。分別采用譜估計方法和本文方法對螺旋槳葉片的FRF進行估計，結(jié)果如圖8所示。2種方法估計的FRF基本吻合，但相較于譜估計方法，由于借助先驗信息和全局建模的優(yōu)勢，本文方法估計的FRF波動更小，即隨機擾動更小。此外，F(xiàn)RF的譜估計在145 Hz附近存在異常凸起，這極可能是高速相機風(fēng)扇運轉(zhuǎn)產(chǎn)生的噪聲峰值[16]，本文方法則能獲得更加準確的FRF估計。

圖8 葉片的FRF估計曲線

根據(jù)式(12)預(yù)測了FRF估計的標準偏差，反映了FRF估計的精度，同時刻畫了FRF估計的不確定度在頻域的分布情況，在120～150 Hz附近不確定度最大，與該頻段存在異常噪聲的事實相切合。

4 結(jié)論

為了提高沖擊載荷下實驗?zāi)B(tài)分析的精度，本文發(fā)展了一種基于GPR方法的FRF非參數(shù)辨識方法。該方法利用貝葉斯學(xué)習(xí)技術(shù)融合了測試系統(tǒng)的樣本信息和模型信息，給出了FRF的最大后驗估計及其方差，其有效性和可靠性在平板沖擊仿真和葉片振動實驗中得到了驗證。相比于經(jīng)典的譜估計方法，本文提出的方法具有更高的FRF估計精度，同時還可給出FRF估計的不確定度描述，為模態(tài)測試等領(lǐng)域提供了有力的工具。在下一步工作中，將致力把該方法擴展至有色噪聲情形。