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

(鄭州大學(xué) 機(jī)械與動(dòng)力工程學(xué)院， 鄭州 450001)

0 引言

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

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

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

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

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

1 問(wèn)題描述

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

(1)

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

對(duì)于線(xiàn)性時(shí)不變動(dòng)態(tài)系統(tǒng)，沖擊激勵(lì)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)系

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

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

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

2.1 復(fù)高斯過(guò)程

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

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

G=Gp+ΔG

(3)

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

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

(4)

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

(5)

(6)

(7)

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

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

(8)

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

2.2 似然函數(shù)

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

(9)

2.3 最大后驗(yàn)估計(jì)

(10)

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

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

(11)

(12)

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

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

(13)

(14)

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

(16)

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

(17)

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

(18)

(19)

(20)

最后，目標(biāo)函數(shù)J(θ)可以化簡(jiǎn)為

(21)

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

(22)

(23)

(24)

式中，Tr表示矩陣求跡運(yùn)算。

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

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

3 算例與實(shí)驗(yàn)

3.1 數(shù)值算例

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

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

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

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

(25)

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

1 dB=20lg|X/X0|

(26)

式中，默認(rèn)參考值X0=1。

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

圖4 平板的FRF估計(jì)曲線(xiàn)

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

圖5 不同SNR下的MSE均值

3.2 振動(dòng)實(shí)驗(yàn)

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

圖6 葉片振動(dòng)實(shí)驗(yàn)

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

圖7 加窗后時(shí)域信號(hào)波形

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

圖8 葉片的FRF估計(jì)曲線(xiàn)

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

4 結(jié)論

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