基于加速度頻響函數(shù)小波變換的貝葉斯模型修正

王增輝， 殷 紅， 彭珍瑞， 張亞峰， 董康立

(1. 蘭州交通大學(xué) 機(jī)電工程學(xué)院，蘭州 730070;2. 浙江大學(xué) 生物醫(yī)學(xué)工程與儀器科學(xué)學(xué)院，杭州 310027)

近年來，模型修正技術(shù)在結(jié)構(gòu)動(dòng)力學(xué)領(lǐng)域逐漸成為研究熱點(diǎn)，廣泛應(yīng)用于結(jié)構(gòu)損傷識(shí)別和健康監(jiān)測(cè)等領(lǐng)域[1]。當(dāng)前絕大多數(shù)模型修正方法屬于確定性方法，但由于結(jié)構(gòu)材料、邊界條件以及試驗(yàn)過程中必然存在不確定性因素，使確定性方法的應(yīng)用受到限制[2-3]。而不確定性模型修正方法結(jié)合概率統(tǒng)計(jì)理論中的大數(shù)定理，利用結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)的統(tǒng)計(jì)特性對(duì)參數(shù)的統(tǒng)計(jì)特性進(jìn)行間接估計(jì)，可充分考慮以上不確定性因素的影響，有效克服確定性模型修正方法的不足，具有良好的應(yīng)用前景[4]。

基于貝葉斯統(tǒng)計(jì)理論的不確定性模型修正方法綜合考慮歷史數(shù)據(jù)和專家經(jīng)驗(yàn)來預(yù)設(shè)待修正參數(shù)的先驗(yàn)分布，然后結(jié)合實(shí)測(cè)響應(yīng)的統(tǒng)計(jì)信息不斷修正先驗(yàn)分布，使之靠近參數(shù)的真實(shí)后驗(yàn)概率分布，廣泛應(yīng)用于不確定性模型修正領(lǐng)域[5]。但是，當(dāng)傳統(tǒng)貝葉斯方法用于復(fù)雜結(jié)構(gòu)的模型修正時(shí)，通常存在似然函數(shù)和后驗(yàn)概率難以求解的缺點(diǎn)，因而通過直接積分的方式求取待修正參數(shù)的后驗(yàn)概率分布往往不可行。Beck等[6]引入馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo, MCMC)算法解決積分運(yùn)算問題。但傳統(tǒng)MCMC算法樣本的拒絕率高、采樣效率低，極易出現(xiàn)“采樣停滯”現(xiàn)象，即連續(xù)多次采樣的解均為同一樣本點(diǎn)，導(dǎo)致樣本集不能充分反映參數(shù)空間的分布特征。針對(duì)傳統(tǒng)MCMC算法采樣效率低的問題，Ching等[7]提出過渡MCMC(transitional MCMC，TMCMC)算法，從一系列中間目標(biāo)概率密度函數(shù)中采樣，提高了采樣效率。劉綱等[8]融合自適應(yīng)算法和相關(guān)向量機(jī)模型，提出了一種快速的不確定性模型修正方法。彭珍瑞等[9]基于最大熵值法，將布谷鳥算法和標(biāo)準(zhǔn)MCMC算法進(jìn)行融合，提高了樣本的接受率。Cheung等[10]提出混合MCMC(hybrid MCMC, HMCMC)算法，以求解維度較高時(shí)不確定參數(shù)的模型修正問題。此外，傳統(tǒng)MCMC算法在每次采樣過程中均需要調(diào)用有限元模型進(jìn)行計(jì)算，導(dǎo)致計(jì)算成本很高。因此，代理模型作為一種替代有限元模型進(jìn)行快速迭代的方法引起了研究者的關(guān)注[11]。代理模型是利用輸入(設(shè)計(jì)參數(shù))與輸出(感興趣響應(yīng))之間的關(guān)系直接構(gòu)造的一種近似模型。其中，Kriging模型是一種基于插值的代理模型，與其他代理模型不同，其利用少量樣本就可以擬合出輸入與輸出之間的關(guān)系，并且不僅能給出參數(shù)的預(yù)估值，還可給出預(yù)測(cè)的誤差，因此廣泛應(yīng)用于精密儀器設(shè)計(jì)和結(jié)構(gòu)優(yōu)化等領(lǐng)域。王巨濤等[12]將Kriging代理模型引入模型修正中進(jìn)行迭代求解，驗(yàn)證了Kriging方法的有效性。秦仙蓉等[13]融合Kriging模型和多目標(biāo)遺傳算法對(duì)岸橋結(jié)構(gòu)有限元模型進(jìn)行了修正，提高了模型修正的效率。

此外，選擇合適的結(jié)構(gòu)響應(yīng)也是模型修正問題的關(guān)鍵。基于模態(tài)參數(shù)的模型修正方法需要對(duì)結(jié)構(gòu)進(jìn)行模態(tài)分析，極易引入模態(tài)識(shí)別誤差?；陬l響函數(shù)(frequency response function, FRF)的模型修正方法不需要進(jìn)行模態(tài)參數(shù)識(shí)別，可較好地避免識(shí)別誤差，同時(shí)，F(xiàn)RF包含大量頻率點(diǎn)信息[14]。但直接將FRF引入貝葉斯方法中，存在似然函數(shù)推導(dǎo)復(fù)雜且計(jì)算量大的問題。徐張明等[15]推導(dǎo)了基于FRF的靈敏度分析方程，利用不完備數(shù)據(jù)對(duì)多個(gè)參數(shù)進(jìn)行了修正。曹詩澤[16]將FRF引入貝葉斯模型修正中，推導(dǎo)了FRF的解析概率模型。

在上述理論背景下，本文提出一種基于加速度FRF小波變換的貝葉斯模型修正方法。首先，引入模態(tài)參與變異系數(shù)準(zhǔn)則[17]選取激勵(lì)點(diǎn)，模態(tài)動(dòng)能法[18]選取測(cè)點(diǎn)。其次，計(jì)算加速度FRF并進(jìn)行小波變換，提取小波總能量作為Kriging模型輸出，待修正參數(shù)作為Kriging模型輸入，采用粒子群算法尋得最優(yōu)相關(guān)系數(shù)以構(gòu)造Kriging模型。然后，以延緩拒絕策略為基礎(chǔ)，融合天牛須算法中更新天牛質(zhì)心位置的方法產(chǎn)生新的候選樣本來估計(jì)待修正參數(shù)的后驗(yàn)概率分布，以提高馬爾科夫鏈的收斂性能。最后，通過車輛三自由度系統(tǒng)和空間桁架結(jié)構(gòu)算例驗(yàn)證了本文方法的有效性。

1 FRF特征量提取

1.1 FRF理論

對(duì)于一個(gè)n自由度阻尼系統(tǒng)，其基本動(dòng)力學(xué)方程可以表示為

(1)

式中：M，C，K分別為質(zhì)量、阻尼和剛度矩陣；F為系統(tǒng)激勵(lì)；X(t)為位移向量。在簡(jiǎn)諧激勵(lì)作用下，經(jīng)傅里葉變換可得到頻域內(nèi)的穩(wěn)態(tài)響應(yīng)X(ω)

X(ω)=H(ω)F(ω)

(2)

式中：H(ω)為加速度FRF矩陣；ω為激勵(lì)頻率；則加速度FRF矩陣可表示為

(3)

1.2 小波變換

小波變換作為一種信號(hào)處理的重要研究方法，具有多分辨率分析的特點(diǎn)，小波總能量是一項(xiàng)重要特征，其對(duì)結(jié)構(gòu)的局部輕微變化十分敏感[19]。選定小波基函數(shù)，一組時(shí)域信號(hào)的連續(xù)小波變換可表示為

(4)

式中： (Wψf)(a,b)為小波系數(shù)；a和b分別為分解尺度和平移參數(shù)；ψ*(t)為小波基函數(shù)ψ(t)的共軛復(fù)數(shù)。(Wψf)(a,b)=〈f(t),ψj,k(t)〉，將a和b離散化，令a=2-j，b=2-jk，j,k∈Z，得離散小波變換

(DWψf)(j,k)=〈f(t),ψj,k(t)〉

(5)

對(duì)信號(hào)進(jìn)行離散小波變換時(shí)，首先進(jìn)行小波分解得到小波系數(shù)，每個(gè)分解尺度的小波能量(小波分量能量)為該層所有小波系數(shù)能量的總和，將各個(gè)尺度下的小波能量相加即為小波總能量。該過程可表示為

(6)

式中：Ej為第j尺度下的小波能量；Cj(k)為第j尺度下的第k個(gè)小波系數(shù)。本文對(duì)加速度FRF進(jìn)行離散小波變換，提取小波總能量替代加速度FRF進(jìn)行模型修正。

2 Kriging模型的構(gòu)造

2.1 Kriging模型理論

Kriging模型是一個(gè)基于隨機(jī)過程的代理模型，擁有對(duì)非線性函數(shù)良好的擬合預(yù)測(cè)能力以及誤差估計(jì)功能[20]。Kriging模型將未知函數(shù)看作高斯過程的具體實(shí)現(xiàn)，該過程包括線性回歸和非參數(shù)部分

(7)

式中：f(x)=[f1(x),f2(x), …,fp(x)]T為多項(xiàng)式函數(shù)；β=[β1,β2, …,βp]T為回歸模型系數(shù)；z(x)為服從正態(tài)分布N(0,σ2)的靜態(tài)隨機(jī)過程。用含兩樣本點(diǎn)空間距離的函數(shù)表示響應(yīng)z(x)之間的相關(guān)性，選取貼合工程實(shí)際的高斯函數(shù)作為相關(guān)函數(shù)，其形式為

(8)

由最大似然法訓(xùn)練模型參數(shù)，可求得

(9)

(10)

式中：F為樣本點(diǎn)向量組成的矩陣；Y為響應(yīng)列向量；R為相關(guān)矩陣，其中元素Rij=R(xi,xj)(i,j=1,2,…,n)，n為試驗(yàn)點(diǎn)數(shù)；σ2和β均為θk的函數(shù)，所以Kriging模型中唯一未知數(shù)即為θk。

2.2 粒子群算法優(yōu)化Kriging模型相關(guān)系數(shù)

文獻(xiàn)[21]介紹了Kriging模型理論及影響其精度的關(guān)鍵因素，Kriging模型的相關(guān)系數(shù)θk決定著代理模型的精度。在MCMC采樣過程中，構(gòu)造滿足精度要求的Kriging模型可改善后驗(yàn)密度與目標(biāo)密度之間的匹配度，進(jìn)一步提高模型修正效率，所以相關(guān)系數(shù)θk的選取對(duì)于Kriging模型的預(yù)測(cè)精度十分重要。

本文將測(cè)試集樣本代入已訓(xùn)練的初始Kriging模型中得到均方誤差(mean square error, MSE)矩陣，然后對(duì)MSE矩陣依次按列、行求取均值，將該均值作為優(yōu)化算法的適應(yīng)度函數(shù)值。選用結(jié)構(gòu)簡(jiǎn)單、尋優(yōu)速度較快的粒子群(particle swarm optimization, PSO)算法[22]進(jìn)行相關(guān)系數(shù)尋優(yōu)。最后結(jié)合最優(yōu)相關(guān)系數(shù)和訓(xùn)練集樣本構(gòu)造Kriging模型，采用均方根誤差(root mean square error, RMSE)檢驗(yàn)其精度。

3 貝葉斯推理及天牛須搜索理論

3.1 貝葉斯理論

基于貝葉斯的模型修正結(jié)合先驗(yàn)信息(主觀信息)和測(cè)試數(shù)據(jù)(客觀信息)，采用MCMC方法推斷修正參數(shù)的后驗(yàn)概率分布。這一過程可用貝葉斯公式表達(dá)為

(11)

式中：D為觀測(cè)信息；θ為所有待修正參數(shù)構(gòu)成的向量； π(θ)為待修正參數(shù)θ的先驗(yàn)分布，通常取為廣義無偏見的均勻分布， π(θ)=1；p(D|θ)為在θ給定下的條件分布，通常稱為似然函數(shù)；c為常數(shù)因子，通常稱為正則化常數(shù)。假設(shè)進(jìn)行N次獨(dú)立試驗(yàn)，則貝葉斯公式中似然函數(shù)的表達(dá)式為

(12)

式中：yi為測(cè)試響應(yīng)；y(θ)i為模型計(jì)算響應(yīng)；covy為測(cè)試信息協(xié)方差矩陣。將式(12)代入式(11),則參數(shù)的后驗(yàn)概率分布可進(jìn)一步寫為

(13)

式中，c′為與θ無關(guān)的常數(shù)，在實(shí)際工程問題中，通常響應(yīng)y(θ)無顯式表達(dá)式，積分運(yùn)算較為困難，故采用MCMC抽樣算法計(jì)算參數(shù)的后驗(yàn)概率分布，以解決復(fù)雜積分運(yùn)算問題。

MH(Metropolis-Hastings)抽樣是目前應(yīng)用最廣的MCMC算法之一，大部分MCMC算法均為標(biāo)準(zhǔn)MH抽樣的擴(kuò)展，其最突出的特點(diǎn)是通過抽樣獲取后驗(yàn)樣本，構(gòu)造合適的馬爾科夫鏈來估計(jì)參數(shù)的統(tǒng)計(jì)特征。MH抽樣的主要步驟為：

步驟1選擇合適的初始樣本θ0，使p(θ0|D)>0，D為響應(yīng)矩陣；

步驟2選擇合適的建議分布q(θ*,θt)，利用當(dāng)前樣本值θt，從建議分布中產(chǎn)生一個(gè)新的候選樣本θ*；

步驟3由當(dāng)前樣本θt和候選樣本θ*，計(jì)算接受概率

α=min[1,p(θ*|D)/p(θt|D)]

(14)

步驟4從均勻分布U(0,1)中產(chǎn)生隨機(jī)數(shù)μ，若α>μ，則接受候選樣本θ*，即θt+1=θ*；否則拒絕候選樣本θ*，即θt+1=θt；

步驟5重復(fù)步驟2～步驟4，直至達(dá)到設(shè)定的抽樣次數(shù)時(shí)終止迭代，產(chǎn)生收斂序列，即馬爾科夫鏈；

步驟6用去掉燃燒期后的馬爾科夫鏈來估計(jì)參數(shù)的統(tǒng)計(jì)特征，從而獲得修正后的參數(shù)值。

因高斯分布在工程實(shí)際中應(yīng)用廣泛且計(jì)算簡(jiǎn)便，故本文選用高斯分布作為建議分布。

3.2 天牛須搜索算法

天牛須搜索算法是一種受天牛覓食行為啟發(fā)的智能優(yōu)化算法，其僅需一個(gè)天牛且能避開復(fù)雜的梯度計(jì)算，在一般的優(yōu)化問題中，可自適應(yīng)地實(shí)現(xiàn)尋優(yōu)[23]。算法分為搜索和檢測(cè)兩個(gè)階段。其中，搜索階段可表示為

(15)

式中：xl和xr分別為天牛左、右須的位置坐標(biāo)；di為i時(shí)刻天牛兩須之間的距離，也稱為感知距離，決定著算法的搜索空間；b為單位化的隨機(jī)方向向量。

檢測(cè)階段產(chǎn)生下一時(shí)刻天牛質(zhì)心的位置，表示為

xi=xi-1+δi·b·sign[F(xl,X)-F(xr,X)]

(16)

式中：F(·)為目標(biāo)函數(shù)，即食物氣味濃度；X為位置坐標(biāo)；δi為i時(shí)刻的搜索步長。引入指數(shù)衰減模型來不斷更新δ和d，衰減模型可表示為

(17)

式中：d0和δ0為感知距離和搜索步長的初始值；ad和aδ為衰減系數(shù)；i為搜索次數(shù)；Hd和Hδ為衰減常數(shù)。

4 基于改進(jìn)MH抽樣的模型修正

4.1 引入天牛局部隨機(jī)搜索的改進(jìn)MH抽樣

由上文可知，當(dāng)參數(shù)維度增加時(shí)，MH抽樣極易出現(xiàn)“停滯”現(xiàn)象。解決該問題較好的方法是采用延緩拒絕(delayed rejection, DR)策略，基本思想是在被拒絕的樣本點(diǎn)基礎(chǔ)上以新的建議分布產(chǎn)生一個(gè)新的候選樣本點(diǎn)。但DR策略中每個(gè)候選點(diǎn)可能需要計(jì)算多個(gè)接受率，計(jì)算量大[24]。考慮貝葉斯推理本質(zhì)是基于最大似然估計(jì)，即尋找盡可能使似然函數(shù)取得最大值的一系列解，由式(12)似然函數(shù)p(D|θ)的表達(dá)式可知，似然函數(shù)最大值對(duì)應(yīng)的解即為指數(shù)項(xiàng)[-L(θ)]最大值所對(duì)應(yīng)的解。將似然函數(shù)的指數(shù)項(xiàng)[-L(θ)]作為天牛須搜索算法的目標(biāo)函數(shù)以引導(dǎo)天牛個(gè)體搜索。考慮DR策略計(jì)算成本問題，當(dāng)候選樣本點(diǎn)被拒絕后，利用天牛個(gè)體更新質(zhì)心位置的方法，在該點(diǎn)基礎(chǔ)上進(jìn)行局部搜索以產(chǎn)生新的候選樣本點(diǎn)，并利用天牛須搜索算法左、右須的比較機(jī)制和Metropolis收斂準(zhǔn)則來保證采樣質(zhì)量。值得說明的是，雖然候選樣本點(diǎn)在Metropolis準(zhǔn)則下被拒絕，但其已經(jīng)位于較好的樣本點(diǎn)附近，因此要將天牛須搜索算法中感知距離和搜索步長的初始值設(shè)置得較小。具體采樣步驟如下：

步驟1初始化馬爾科夫鏈θ(t)(t=1,2,…,n)，選擇建議分布q(θ*,θt)，其中，t為迭代次數(shù)，n為總迭代次數(shù)；

步驟2獲取馬爾科夫鏈的初始樣本點(diǎn)θ0，即從先驗(yàn)分布中抽樣產(chǎn)生θ0；

步驟3在第t時(shí)刻，當(dāng)前樣本為θt，由MH抽樣產(chǎn)生新的候選樣本θ*，計(jì)算接受概率α并判斷是否接受候選樣本；

步驟4當(dāng)α>μ時(shí)，接受該樣本θ*，即θt+1=θ*；當(dāng)α<μ時(shí)，進(jìn)入下一步驟；

步驟5以DR策略為基礎(chǔ)，融入天牛須搜索算法中天牛個(gè)體更新質(zhì)心位置的思想：

(a) 初始化搜索參數(shù)，將步驟4的候選樣本θ*作為天牛初始的質(zhì)心位置，即x0=θ*；

(b) 以式(12)中似然函數(shù)表達(dá)式的指數(shù)項(xiàng)[-L(x)]作為搜索的目標(biāo)函數(shù)，引導(dǎo)天牛個(gè)體進(jìn)行隨機(jī)搜索；

(c) 感知距離di和搜索步長δi按式(17)隨搜索次數(shù)i的增加而遞減，保證天牛自適應(yīng)地更新其質(zhì)心位置；

(d) 達(dá)到設(shè)定的搜索次數(shù)時(shí)終止搜索并根據(jù)天牛質(zhì)心位置(即新的候選樣本)θζ，計(jì)算接受率αζ(依據(jù)MH算法判斷接受狀態(tài)θt+1=θζ，否則θt+1=θt)；

步驟6增加迭代次數(shù)t，重復(fù)步驟3～步驟5，直至產(chǎn)生一個(gè)收斂序列，即馬爾科夫鏈，停止迭代。

4.2 模型修正流程

本文提取小波總能量作為FRF的響應(yīng)特征量，將FRF間接引入到貝葉斯模型修正中。為考慮馬爾科夫鏈的收斂性能以及計(jì)算成本問題，本文從全局抽樣和局部隨機(jī)游走兩個(gè)角度出發(fā)，構(gòu)建了一種貝葉斯模型修正框架。具體的模型修正流程如圖1所示。

圖1 模型修正流程圖Fig.1 Flow chart of model updating

5 激勵(lì)點(diǎn)和測(cè)點(diǎn)的選取原則

針對(duì)模態(tài)試驗(yàn)中激勵(lì)點(diǎn)和測(cè)點(diǎn)的選擇問題，本文采用模態(tài)參與變異系數(shù)準(zhǔn)則選取最優(yōu)激勵(lì)點(diǎn)，模態(tài)動(dòng)能法選取最優(yōu)測(cè)點(diǎn)。

5.1 模態(tài)參與變異系數(shù)準(zhǔn)則

模態(tài)參與度用來評(píng)價(jià)結(jié)構(gòu)各個(gè)自由度對(duì)激發(fā)感興趣模態(tài)的貢獻(xiàn)，可表示為

(18)

式中：p和q分別為輸出和輸入自由度；N為結(jié)構(gòu)的自由度總數(shù)；r為感興趣的模態(tài)階數(shù)；Apqi為留數(shù)。

假設(shè)結(jié)構(gòu)為比例阻尼，則式(3)的FRF可表示為

(19)

式中：φpi和φqi分別為模態(tài)矩陣的第(i,p)項(xiàng)和第(i,q)項(xiàng)；ωi為第i階模態(tài)頻率；ξi為第i階模態(tài)阻尼。式(18)中的留數(shù)可進(jìn)一步表示為

Apqi=φpiφqi

(20)

結(jié)構(gòu)的第q自由度對(duì)激發(fā)所有感興趣模態(tài)的貢獻(xiàn)值可表示為

(21)

Wq值越大，說明該自由度對(duì)激發(fā)感興趣模態(tài)的貢獻(xiàn)度越大。由于本文的模態(tài)試驗(yàn)為多輸入多輸出的方式，故選取貢獻(xiàn)度較大的自由度作為激勵(lì)點(diǎn)。

5.2 模態(tài)動(dòng)能法

模態(tài)動(dòng)能法是一種量化的測(cè)點(diǎn)配置方法。依據(jù)模態(tài)理論可知，結(jié)構(gòu)的自由振動(dòng)可表示為各階模態(tài)的線性疊加，因此結(jié)構(gòu)的實(shí)際振動(dòng)情況可通過模態(tài)振型進(jìn)行描述。模態(tài)振型動(dòng)能值可以反映某點(diǎn)對(duì)模態(tài)振型的貢獻(xiàn)度，如果某點(diǎn)的模態(tài)動(dòng)能值最大，則說明該點(diǎn)是振動(dòng)最敏感的點(diǎn)，即最優(yōu)測(cè)點(diǎn)。模態(tài)動(dòng)能理論可由式(22)表示

MKE=φTMφ

(22)

式中：φ為模態(tài)振型矩陣；M為結(jié)構(gòu)質(zhì)量矩陣；考慮邊界條件，M為正定矩陣，對(duì)其進(jìn)行正交分解得

M=LTL

(23)

則式(22)可進(jìn)一步寫為

MKE=φTLTLφ=ψTψ=A0

(24)

E=ψ[ψTψ]-1ψT

(25)

按照矩陣E對(duì)角元素的大小對(duì)測(cè)點(diǎn)進(jìn)行排序，每次迭代中刪除模態(tài)動(dòng)能最小值對(duì)應(yīng)的測(cè)點(diǎn)，直至保留到滿足要求的測(cè)點(diǎn)數(shù)目為止。

6 算例及結(jié)果分析

6.1 算例一：車輛三自由度系統(tǒng)模型

如圖2所示：m3為車體質(zhì)量；m2為轉(zhuǎn)向架質(zhì)量；k3和c2為二系懸掛裝置的剛度和阻尼；k2和c1為一系懸掛裝置的剛度和阻尼；k1為減震器端部剛度；c1，m1和k1構(gòu)成一個(gè)減震器模型。

圖2 車輛三自由度系統(tǒng)模型Fig.2 Model of 3DOFs vehicle

采用均勻試驗(yàn)設(shè)計(jì)方法在參數(shù)范圍內(nèi)抽取200組數(shù)據(jù)作為訓(xùn)練集，40組數(shù)據(jù)作為測(cè)試集。利用1.1節(jié)的式(3)計(jì)算加速度FRF并進(jìn)行小波變換，提取小波總能量作為Kriging模型的輸出，構(gòu)造Kriging模型。

分別采用4.1節(jié)方法和標(biāo)準(zhǔn)MH算法進(jìn)行采樣，經(jīng)過20 000次迭代后得到參數(shù)樣本的馬爾科夫鏈，如圖3所示。由圖3(a)和圖3(b)對(duì)比可知，兩種方法得到的馬爾科夫鏈均能收斂到預(yù)設(shè)的參數(shù)均值，但標(biāo)準(zhǔn)MH算法得到的馬爾科夫鏈出現(xiàn)很多平滑段，說明采樣過程中出現(xiàn)“采樣停滯”現(xiàn)象，即連續(xù)多次重復(fù)同一樣本，雖然樣本協(xié)方差小，但樣本不能充分反映參數(shù)空間的分布特征；而所提算法得到的馬爾科夫鏈能夠快速收斂且一直保持波動(dòng)狀態(tài)，樣本遍歷性明顯優(yōu)于前者。同時(shí)，由式(14)計(jì)算樣本的接受率，相比標(biāo)準(zhǔn)MH算法21%的接受率，本文所提算法的接受率為62%左右，樣本的接受率明顯提高。

圖3 不同算法產(chǎn)生的馬爾科夫鏈Fig.3 Markov chains generated by different algorithms

依據(jù)概率統(tǒng)計(jì)中關(guān)于正態(tài)概率檢驗(yàn)的相關(guān)知識(shí)對(duì)圖3的后驗(yàn)樣本進(jìn)行分析：因非平穩(wěn)期樣本會(huì)對(duì)統(tǒng)計(jì)特征造成干擾，影響參數(shù)均值的估計(jì)精度，所以先剔除圖3(b)中馬爾科夫鏈前10%樣本，再對(duì)剩余樣本進(jìn)行正態(tài)概率檢驗(yàn)，結(jié)果如圖4所示。由圖4可知，去除非平穩(wěn)的樣本等劣質(zhì)樣本后，本文方法獲得的3個(gè)參數(shù)的絕大多數(shù)后驗(yàn)樣本都能落在假設(shè)的正態(tài)分布線上，所以樣本的質(zhì)量較高，正態(tài)分布假設(shè)成立。

依據(jù)后驗(yàn)樣本的統(tǒng)計(jì)特征估計(jì)修正后的參數(shù)均值，結(jié)果如表1所示。由表1可知，經(jīng)所提方法修正后的參數(shù)誤差皆可控制在0.3%以內(nèi)，且修正誤差均優(yōu)于標(biāo)準(zhǔn)MH算法，所提方法可以獲得較好的修正效果。

圖4 樣本的正態(tài)概率檢驗(yàn)圖Fig.4 Normal probability test diagram of samples

表1 車輛三自由度系統(tǒng)修正前后參數(shù)均值Tab.1 Mean values of parameters before and after updating of the vehicle 3DOFs system

6.2 算例二：空間桁架結(jié)構(gòu)模型

空間桁架結(jié)構(gòu)如圖5所示。該結(jié)構(gòu)包括28個(gè)節(jié)點(diǎn)，66個(gè)桿單元，節(jié)點(diǎn)鉸接，約束條件為4個(gè)支座固定(節(jié)點(diǎn)編號(hào)1,8,9,16)，每個(gè)節(jié)點(diǎn)只考慮Y向和Z向的平動(dòng)自由度，依次對(duì)所有節(jié)點(diǎn)的Y、Z方向的自由度進(jìn)行編號(hào)，共48個(gè)自由度，將參數(shù)的初始均值賦予該結(jié)構(gòu)作為模型修正的初始模型。利用第5.1節(jié)的模態(tài)參與變異系數(shù)準(zhǔn)則選取激勵(lì)點(diǎn)，結(jié)構(gòu)各個(gè)自由度的Wq值如圖6所示。

圖5 空間桁架結(jié)構(gòu)模型Fig.5 Structure model of the space truss

圖6 模態(tài)參與度值Fig.6 Values of modal participation

結(jié)合圖6的模態(tài)參與度和實(shí)測(cè)方便性，選取桁架第3節(jié)點(diǎn)Y方向自由度、第6節(jié)點(diǎn)Y方向自由度和第11節(jié)點(diǎn)Y方向自由度作為激勵(lì)點(diǎn)。值得說明的是，由圖6可知，對(duì)于激發(fā)前4階模態(tài)貢獻(xiàn)較大的均為Y方向自由度，這也與桁架的結(jié)構(gòu)特點(diǎn)相符。

利用5.2節(jié)的模態(tài)動(dòng)能法選取測(cè)點(diǎn)，計(jì)算各個(gè)自由度的模態(tài)動(dòng)能，選擇結(jié)構(gòu)第18節(jié)點(diǎn)Y方向自由度、第24節(jié)點(diǎn)Y方向自由度和第26節(jié)點(diǎn)Y方向自由度作為測(cè)點(diǎn)，激勵(lì)點(diǎn)和測(cè)點(diǎn)位置見圖5。

將桁架分為上弦桿、中部桿(直腹桿和斜腹桿)和下弦桿三部分。因大多數(shù)情況下發(fā)生的損傷會(huì)導(dǎo)致結(jié)構(gòu)剛度發(fā)生顯著變化，主要表現(xiàn)為彈性模量存在不確定性。選取三組桿件的彈性模量作為待修正參數(shù)，分別為θ1=E1/E0，θ2=E2/E0，θ3=E3/E0，其中初始彈性模量E0為190 GPa，假設(shè)各參數(shù)實(shí)測(cè)均值服從正態(tài)分布矩陣[1.5，1.2，0.8]，標(biāo)準(zhǔn)差都為0.1。設(shè)置參數(shù)均值的取值范圍分別為[0.9，2.1]，[0.6，1.8]和[0.3，1.3]。

利用式(3)計(jì)算對(duì)應(yīng)激勵(lì)點(diǎn)和測(cè)點(diǎn)的加速度FRF并進(jìn)行小波變換，選取小波基函數(shù)為“db1”，分解尺度數(shù)為3，提取小波總能量作為模型修正的響應(yīng)特征量。采用2.2節(jié)方法構(gòu)造滿足精度要求的Kriging模型，選用4階小波總能量的RMSE值來檢驗(yàn)所構(gòu)造的Kriging模型的精度，將程序運(yùn)行10次，統(tǒng)計(jì)結(jié)果如表2所示。由表2可知，構(gòu)造的Kriging模型預(yù)測(cè)精度很高，可替代有限元模型進(jìn)行迭代計(jì)算。分別采用本文4.1節(jié)方法和標(biāo)準(zhǔn)MH算法進(jìn)行采樣，采用4.1節(jié)方法所抽取樣本的分布云圖以及相應(yīng)的置信橢圓如圖7所示，橢圓內(nèi)的區(qū)域表示含95%樣本的區(qū)域。由圖7可知，本文方法所抽取的絕大多數(shù)樣本點(diǎn)都能位于95%的高概率區(qū)域，后驗(yàn)樣本的質(zhì)量較高。

表2 Kriging代理模型精度檢驗(yàn)Tab.2 Accuracy test of Kriging surrogate model

圖7 分布云圖及置信橢圓Fig.7 Distribution cloud maps and confidence ellipse

采用標(biāo)準(zhǔn)MH算法和本文所提方法采樣20 000次獲得的馬爾科夫鏈，如圖8所示。由圖8(a)和圖8(b)對(duì)比可知，所提算法得到的馬爾科夫鏈能夠一直保持波動(dòng)狀態(tài)，很少出現(xiàn)標(biāo)準(zhǔn)MH算法的“采樣停滯”現(xiàn)象，樣本的遍歷性明顯優(yōu)于標(biāo)準(zhǔn)MH算法。由式(14)分別計(jì)算兩種方法所抽取樣本的接受率，標(biāo)準(zhǔn)MH算法樣本的接受率約為21%，本文所提方法約為52%，樣本接受率明顯提高。

圖8 不同算法獲得的馬爾科夫鏈Fig.8 Markov chains generated by different algorithms

由馬爾科夫鏈估計(jì)參數(shù)均值，參數(shù)均值的修正結(jié)果如表3所示。由表3可知，參數(shù)均值的修正誤差皆小于0.8%，且3個(gè)參數(shù)的修正誤差均優(yōu)于標(biāo)準(zhǔn)MH算法，表明所提方法進(jìn)行不確定性模型修正取得了很好的效果。為進(jìn)一步驗(yàn)證所提方法的修正效果，使用表3中的參數(shù)均值分別計(jì)算試驗(yàn)、修正后以及初始模型的加速度FRF，其對(duì)應(yīng)的FRF曲線，如圖9所示。由圖9可知，修正后均值模型加速度FRF曲線與試驗(yàn)均值模型加速度FRF曲線充分接近，基本重合。同時(shí)，在處理器為AMD R7 4800U的計(jì)算機(jī)上運(yùn)行DR算法，DR算法修正時(shí)間為48.6 s，改進(jìn)算法修正時(shí)間為42.3 s，本文方法修正時(shí)間小于DR算法。

表3 空間桁架結(jié)構(gòu)修正前后參數(shù)均值Tab.3 Mean values of parameters before and after updating of the space truss structure

圖9 頻響函數(shù)曲線Fig.9 Frequency response function curves

7 模型修正影響因素分析

7.1 考慮噪聲影響的模型修正

由于實(shí)際試驗(yàn)過程中，考慮測(cè)得的FRF信號(hào)將不可避免地受到噪聲的干擾。因此，為了驗(yàn)證本文方法的抗噪性，對(duì)算例二中的加速度FRF信號(hào)施加一定的高斯白噪聲來模擬實(shí)測(cè)數(shù)據(jù)的噪聲污染，信噪比分別為30 dB，20 dB，10 dB和5 dB。將程序運(yùn)行10次取均值，修正后各參數(shù)的后驗(yàn)均值誤差如圖10所示。由圖10可知，隨著噪聲增大，雖然各參數(shù)的修正誤差也逐漸增大，但基本都能保持在5%以內(nèi)，表明所提方法對(duì)噪聲具有良好的魯棒性。

圖10 不同噪聲水平下的參數(shù)修正誤差Fig.10 Parameter updating errors at different noise levels

7.2 考慮小波基函數(shù)種類影響的模型修正

小波基函數(shù)的選取對(duì)于離散小波變換十分重要，采用不同的小波基函數(shù)對(duì)同一分析信號(hào)進(jìn)行變換得到的小波系數(shù)不同。為驗(yàn)證本文方法在不同小波基函數(shù)下的修正效果，本文以算例二的加速度FRF為分析信號(hào)，選取工程中常用的db4小波、haar小波、coif3小波和sym6小波作為小波變換的基函數(shù)，按照?qǐng)D1的流程進(jìn)行模型修正。修正后得到的馬爾科夫鏈如圖11所示。圖11中的小波基函數(shù)1～4分別對(duì)應(yīng)于上述4種小波。將程序運(yùn)行10次后取均值，修正結(jié)果如表4所示。

圖11 不同小波基函數(shù)下的馬爾科夫鏈Fig.11 Markov chains of different wavelet basis functions

表4 不同種類小波基函數(shù)下的修正誤差Tab.4 Updating errors of different wavelet basis functions %

由圖11和表4可知，在以上4種小波基函數(shù)的情況下，本文方法所獲得的馬爾科夫鏈均能較好地收斂到參數(shù)的預(yù)設(shè)均值，且參數(shù)的修正誤差基本保持在1.5%以內(nèi)，說明所提方法對(duì)于小波基函數(shù)種類的選取具有魯棒性。值得說明的是，參數(shù)θ1的馬氏鏈較另外兩個(gè)參數(shù)的馬氏鏈波動(dòng)較大，反映出參數(shù)θ1的靈敏度略差，但不影響對(duì)小波基函數(shù)魯棒性的判斷。

7.3 考慮分解尺度數(shù)影響的模型修正

為研究小波分解尺度對(duì)于修正效果的影響，以算例二的加速度FRF為分析信號(hào)，選取‘db1’為小波基函數(shù)，分別在不同分解尺度數(shù)下按照?qǐng)D1的流程進(jìn)行模型修正。限于篇幅，僅給出分解尺度數(shù)為2層、3層、4層和5層的情況下所得到的馬爾科夫鏈，如圖12所示。將程序運(yùn)行10次后取均值，修正誤差如表5所示。

圖12 不同小波分解尺度下的馬爾科夫鏈Fig.12 Markov chains of different wavelet decomposition scales

表5 不同小波分解尺度下的修正誤差Tab.5 Updating errors of different wavelet decomposition scales %

綜合圖12和表5可知，分解尺度數(shù)為4和5時(shí)的馬氏鏈波動(dòng)略大，其與不同分解尺度數(shù)下構(gòu)造的Kriging模型的精度不同有關(guān)，故在模型修正過程中要盡可能地保證代理模型的精度。但從整體來看，本文方法在不同分解尺度下獲得的馬爾科夫鏈均能收斂到預(yù)設(shè)的參數(shù)均值，且修正精度基本保持在2%以內(nèi)，小波分解尺度對(duì)修正結(jié)果基本沒有影響。分析原因是：不同的分解尺度會(huì)改變小波系數(shù)在特定尺度下的分布情況，但不會(huì)改變小波總能量值，所以提出的模型修正方法對(duì)于小波分解尺度具有良好的包容性。

8 結(jié) 論

(1) 修正方法將FRF間接引入到貝葉斯模型修正中，無需進(jìn)行模態(tài)分析與振型匹配，避開了傳統(tǒng)方法中頻率點(diǎn)選擇困難的問題，回避了基于FRF的貝葉斯方法中似然函數(shù)和后驗(yàn)概率分布的復(fù)雜推導(dǎo)過程，有效減輕了計(jì)算負(fù)擔(dān)。

(2) 在MCMC框架下，結(jié)合全局建議方案的快速搜索特性和局部隨機(jī)游走的搜索魯棒性，在DR策略的基礎(chǔ)上，引入更新天牛質(zhì)心位置的方法來估計(jì)參數(shù)的后驗(yàn)概率分布，可以在保證較優(yōu)修正精度的同時(shí)提高馬爾科夫鏈樣本的接受率。

(3) 在噪聲污染、不同小波基函數(shù)以及不同小波分解尺度的情況下，基于小波總能量的模型修正方法仍能保持較好的修正效果，對(duì)于工程實(shí)際具有一定的參考價(jià)值。

本文只研究了改進(jìn)MCMC方法和標(biāo)準(zhǔn)MH算法在三維度情況下的修正效果。為進(jìn)一步證實(shí)該方法的使用效果，應(yīng)采用其他更加復(fù)雜的結(jié)構(gòu)，并在激勵(lì)和響應(yīng)上添加噪聲進(jìn)行研究。作者已開展了相關(guān)初步研究并取得了較好的結(jié)果。