曹明明，彭珍瑞，劉滿東

(蘭州交通大學(xué) 機電工程學(xué)院，蘭州 730070)

有限元模型修正技術(shù)的發(fā)展，使有限元模型與實際結(jié)構(gòu)的參數(shù)誤差縮小，提升有限元模型模擬實際結(jié)構(gòu)的能力，這對結(jié)構(gòu)健康監(jiān)測、結(jié)構(gòu)優(yōu)化設(shè)計、預(yù)測等工程應(yīng)用具有重要意義。而且有限元模型精度的提高可以減少不必要的試驗，這可以節(jié)約成本及縮短研發(fā)周期[1]。近年來，基于模態(tài)參數(shù)和基于頻響函數(shù)(FRF)的修正方法在結(jié)構(gòu)動力學(xué)模型修正中較為常見[2]，基于模態(tài)參數(shù)的模型修正方法需要模態(tài)參數(shù)識別，所以在這個過程中不可避免地會引入識別誤差。而且，在某些情況下，在工程試驗中模態(tài)參數(shù)的識別誤差比有限元模型建模誤差還大[3]。而基于頻響函數(shù)的模型修正方法避開了結(jié)構(gòu)模態(tài)分析，直接利用試驗測量得到的頻響函數(shù)來進行模型修正，從而避免了人為引入誤差導(dǎo)致最后修正效果較差，而且這種方法適合于修正模態(tài)分布相對密集的結(jié)構(gòu)[4-5]。

不管用何種方法，目前結(jié)構(gòu)有限元模型修正技術(shù)通常是經(jīng)過大量的迭代優(yōu)化搜索來實現(xiàn)的，在這個過程中需要反復(fù)的調(diào)用有限元模型進行解析計算，但是對于一些大型復(fù)雜結(jié)構(gòu)，解析計算其有限元模型不僅需要花費大量時間，而且對于計算機的性能要求比較高，甚至有些求解無法實現(xiàn)。于是，作為一種可以替代有限元模型進行快速分析的代理模型技術(shù)，引起了眾多學(xué)者的關(guān)注。代理模型通過建立一個顯式的函數(shù)來代替模型參數(shù)與結(jié)構(gòu)響應(yīng)之間復(fù)雜的隱式關(guān)系，也稱為近似模型、響應(yīng)面模型或元模型，近年來逐漸被應(yīng)用到結(jié)構(gòu)模型修正與確認中[6]。常用的代理模型有響應(yīng)面法(Response Surface Methodology,RSM)、徑向基函數(shù)(Radial Basis Function,RBF)、神經(jīng)網(wǎng)絡(luò)(Neural Network,NN)、支持向量回歸機(Support Vector Regression,SVR)、Kriging模型等[7]。Kriging模型是基于Kriging插值技術(shù)的一種等效模型，與其他代理模型不同的是，Kriging模型除了能給出對未知函數(shù)的預(yù)估值，還能得到預(yù)估值的誤差估計[8]，只用少量的樣本就可以較準確的描述結(jié)構(gòu)輸入與輸出間的關(guān)系，在航天飛行器設(shè)計和結(jié)構(gòu)優(yōu)化設(shè)計等領(lǐng)域應(yīng)用較為廣泛。文獻[9-10]將Kriging方法應(yīng)用于頻響函數(shù)模型修正中，減少了求解時間，提高了頻響函數(shù)模型修正的效率。

頻響函數(shù)包含了大量頻率點的數(shù)據(jù)信息，所以頻率點的選擇對基于頻響函數(shù)的模型修正方法至關(guān)重要，但是目前對于頻率點數(shù)據(jù)的選擇方法仍是無法統(tǒng)一。張勇等[11]將頻響函數(shù)變換為時域內(nèi)的脈沖響應(yīng)函數(shù)，通過對脈沖響應(yīng)函數(shù)的重構(gòu)相空間矩陣進行截斷奇異值分解(Singular Value Decomposition，SVD)，提取前若干較大奇異值作為原頻響函數(shù)的特征量來進行模型修正。這種方法間接使用頻響函數(shù)來進行模型修正，既保留了使用頻響函數(shù)進行模型修正的優(yōu)點，又避免了頻率點的選擇難題。文獻[12]中比較了對信號通過連續(xù)截斷和重構(gòu)吸引子構(gòu)造矩陣進行奇異值分解優(yōu)劣性，發(fā)現(xiàn)第二種構(gòu)造方法，矩陣具有良好的消噪效果，能夠有效地消除信號中的隨機噪聲，并且還具有弱信號提取能力。文獻[13]提到使用脈沖響應(yīng)函數(shù)在判斷奇異值分解時的有效秩的閾值時具有優(yōu)越性。但是使用脈沖響應(yīng)函數(shù)首先要將頻響函數(shù)進行正、逆傅里葉變換，這無疑加大了計算量。而且在不同噪聲水平下受噪聲影響的奇異值個數(shù)是變化的，恰當(dāng)?shù)钠娈愔祩€數(shù)不僅能保留住信號的主要特征信息，還能削弱隨機噪聲的干擾，所以需對保留的奇異值個數(shù)進行選擇。文獻[14]通過直接利用頻響函數(shù)構(gòu)造吸引子矩陣進行奇異值分解，并利用重構(gòu)信號中極值點數(shù)量突變點來選取有效秩，成功的對信號進行了降噪。

在上述背景下，為避免正、逆傅里葉變換，提高計算效率，同時避開頻率點的選擇難題，本文直接利用計算得到的頻響函數(shù)重構(gòu)吸引子矩陣，對其進行奇異值分解。根據(jù)極值點數(shù)量突變原則選擇保留主要特征信息的奇異值個數(shù)，以模型設(shè)計參數(shù)作為輸入，保留的奇異值作為結(jié)構(gòu)響應(yīng)，并通過粒子群優(yōu)化算法對Kriging模型相關(guān)系數(shù)進行尋優(yōu)來構(gòu)建Kriging代理模型。以實測的頻響函數(shù)奇異值與代理模型輸出奇異值差值的平方最小作為目標函數(shù)，通過布谷鳥算法求解模型待修正參數(shù)。數(shù)值算例表明了本文方法的有效性，特別是對于隨機噪聲的魯棒性。

1 頻響函數(shù)特征量選取

1.1 頻響函數(shù)理論

頻域內(nèi)，一個n自由度阻尼系統(tǒng)在簡諧激勵作用下的運動方程可以表示為：

(-ω2M+iωC+K)X(ω)=F(ω)

(1)

式中:M、K和C分別表示質(zhì)量矩陣、剛度矩陣和阻尼矩陣；ω為激勵頻率；F(ω)為簡諧激勵；穩(wěn)態(tài)響應(yīng)X(ω)可表示為：

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

(2)

式中，H(ω)為頻響函數(shù)。其中，頻響函數(shù)矩陣為：

H(ω)=(-ω2M+iωC+K)-1

(3)

1.2 奇異值分解理論

將一離散原始信號X=[x1,x2,…,xN]，利用相空間重構(gòu)理論重構(gòu)吸引子(Hankel)矩陣：

(4)

式中：n=N-(m-1)τ；A為m×n階的Hankel矩陣；τ為延遲步長，一般取τ=1；m為嵌入維數(shù)，當(dāng)N是偶數(shù)時，m=N/2，當(dāng)N是奇數(shù)時，m取中值[15]。

對Hankel矩陣A進行奇異值分解可得到：

A=UΣVT

(5)

式中:U與VT為正交矩陣；Σ=diag(σ1,σ2,…,σk)為對角陣；σi(i=1～k)為矩陣A的奇異值，σ1≥σ2≥…≥σk≥0，k=min[(m-1)τ+1,n]。

給定閾值λ，可以將式(5)寫成如下形式：

(6)

(7)

奇異值分解技術(shù)多用于消除測試信號中隨機噪聲的干擾，后k-r個奇異值通常較小，往往是噪聲信號。由于實測頻響函數(shù)中不可避免的會出現(xiàn)隨機噪聲干擾，所以對于實測頻響函數(shù)而言，選擇前r個奇異值作為頻響函數(shù)特征量可以有效消除噪聲的干擾。

1.3 極值點數(shù)量突變原則

由上可知，選擇奇異值的個數(shù)r是很重要的，恰當(dāng)?shù)钠娈愔祩€數(shù)不僅能保留住信號的主要特征信息，還能削弱隨機噪聲的干擾。本文利用文獻[14]所提極值點數(shù)量突變原則選擇有效奇異值個數(shù)。其基本原理為：理想信號通常是一條光滑的曲線，而當(dāng)信號受到噪聲干擾時，會有大量“毛刺”的出現(xiàn)。所以要使噪聲干擾最小，可以通過對重構(gòu)信號中“毛刺”的數(shù)量最少來確定，即信號中的極值點數(shù)量最少，而信號中“毛刺”數(shù)量隨保留奇異值數(shù)目的變化而變化，選取極值點數(shù)量突變點處奇異值數(shù)目作為有效奇異值數(shù)目。

綜上，本文考慮直接利用計算得到的頻響函數(shù)重構(gòu)吸引子矩陣，對其進行奇異值分解。根據(jù)極值點數(shù)量突變原則選擇保留主要特征信息的奇異值，用其表征原頻響函數(shù)來進行模型修正。

2 基于Kriging模型的模型修正理論

2.1 Kriging模型的建立

Kriging模型是一個基于隨機過程的代理模型，包含線性回歸部分和非參數(shù)部分：

(8)

式中:β為回歸模型系數(shù)；f(x)為樣本點x的多項式函數(shù)；z(x)為服從正態(tài)分布N(0,σ2)的隨機過程。任意兩樣本點的z(x)之間的相關(guān)性可以表示為其空間距離的函數(shù)，選計算效率較好的高斯函數(shù)為相關(guān)函數(shù)，其形式如下：

(9)

根據(jù)極大似然函數(shù)法，可以求得：

(10)

(11)

式中:F為樣本點的向量組成的矩陣；Y為樣本點響應(yīng)組成的列向量；R為空間相關(guān)矩陣，其中元素Rij=R(xi,xj)(i,j=1,2,…,n)，n為試驗點數(shù)；σ2和β均為θk的函數(shù)，那么Kriging模型中唯一未知數(shù)即為相關(guān)系數(shù)θk，其決定著預(yù)測響應(yīng)精度，本文采用粒子群算法對相關(guān)系數(shù)θk進行尋優(yōu)。

(12)

式中，rT(x0)=[R(x0,x1),R(x0,x2),…,R(x0,xn)]

淺談道路橋梁施工的常見問題及質(zhì)量檢測技術(shù)的應(yīng)用…………………………………………………… 田文澤（11-100）

2.2 粒子群算法優(yōu)化Kriging模型

Kriging模型中的相關(guān)系數(shù)θk決定著預(yù)測響應(yīng)值精度，只有當(dāng)所建立Kriging模型精度足夠高時，代理模型對修正結(jié)果的誤差影響才會最小。所以相關(guān)系數(shù)的確定對于構(gòu)建代理模型至關(guān)重要。游海龍等[16]討論了利用傳統(tǒng)數(shù)值優(yōu)化方法確定Kriging模型相關(guān)系數(shù)存在依賴搜索起始點等缺點，并利用遺傳算法對相關(guān)系數(shù)進行了尋優(yōu)，解決了對于起始點的依賴問題。粒子群算法(Partical Swarm Optimization,PSO)[17]由Kennedy和Eberhart提出,屬于進化算法的一種,具有算法結(jié)構(gòu)簡單、尋找最優(yōu)解速度快的優(yōu)點。根據(jù)以上分析，本文打算利用PSO優(yōu)化Kriging模型相關(guān)系數(shù)構(gòu)建滿足修正要求精度的Kriging模型。首先采用拉丁超立方抽樣方法抽取一定數(shù)量的樣本點，然后把抽取的樣本點分成訓(xùn)練集和測試集。以訓(xùn)練集作為Kriging模型的輸入變量，以所保留的r個奇異值作為響應(yīng)值來構(gòu)建Kriging模型，以測試集Kriging模型的均方誤差(Mean Squared Error,MSE)均值平均值作為粒子群算法的適應(yīng)度函數(shù)，尋得最優(yōu)相關(guān)系數(shù)。最后以訓(xùn)練集作為樣本點建立Kriging模型。

2.3 模型修正流程

當(dāng)Kriging模型構(gòu)造完成時，模型修正歸根結(jié)底是一優(yōu)化問題，在目標函數(shù)最小的條件下，通過求解該優(yōu)化問題，得到設(shè)計參數(shù)的修正值。目標函數(shù)為：

(13)

在優(yōu)化求解時，除了靈敏度方法外，智能算法也是很簡便的。布谷鳥算法[18]由于特有的萊維性能，使其具有較強的全局搜索能力，而且它參數(shù)少，操作簡單，易于實現(xiàn)，尋優(yōu)性能佳。故本文選用布谷鳥算法進行優(yōu)化求解。模型修正流程如圖1所示。

圖1 模型修正流程圖

3 數(shù)值算例

3.1 算例一

通過一個車輛三自由度垂向頻率響應(yīng)模型來驗證本文方法，該模型如圖2所示。

圖2 車輛三自由度垂向頻率響應(yīng)模型

圖中k3和c2為二系懸掛剛度和阻尼，k2和c1為一系懸掛剛度和阻尼矩陣，k1為減震器端部剛度(串聯(lián)剛度)，c1、m1和k1組成了一個完整的減震器模型，m2為轉(zhuǎn)向架質(zhì)量，m3為車體質(zhì)量。

系統(tǒng)運動方程為：

(14)

以上文提到的數(shù)據(jù)為理論模型參數(shù)，試驗驗證參數(shù)將減震器端部剛度k1減小20%，將二系懸掛剛度k3增加20%。以這兩個參數(shù)為修正設(shè)計參數(shù)，在其理論模型初始值的±40%變化區(qū)間里，采用拉丁超立方抽樣抽取80個樣本點。用前40個數(shù)據(jù)作為訓(xùn)練集，后40個數(shù)據(jù)作為測試集。采用前文所述方法進行模型修正，修正結(jié)果如表1所示，其第1個奇異值得到的擬合響應(yīng)曲面，如圖3所示。

表1 車輛三自由度系統(tǒng)的模型修正結(jié)果

圖3 車輛三自由度系統(tǒng)的Kriging響應(yīng)面

由圖3可知Kriging模型中樣本響應(yīng)和預(yù)測響應(yīng)重合度很好，由表1可知用本文方法進行模型修正所得修正誤差很小，證明用奇異值來表征頻響函數(shù)作為結(jié)構(gòu)響應(yīng)進行模型修正是有效的，并且也達到了很高的修正精度。

3.2 算例二

實際試驗中，測試數(shù)據(jù)不可避免的會受到噪聲的干擾。所以為進一步驗證本文方法對噪聲的抗干擾性，選擇一三維桁架結(jié)構(gòu)來驗證。該結(jié)構(gòu)包括28個節(jié)點，66個單元和48個自由度。桁架約束條件為4個支座固定(節(jié)點編號1,8,9,16)，節(jié)點鉸接，每個節(jié)點只考慮豎向Y和橫向Z的平動自由度。再運用動力學(xué)計算其質(zhì)量矩陣、剛度矩陣，將其代入式(3)中計算出其頻響函數(shù)矩陣。由于本文所用方法避開了頻率點選擇問題，只需保證頻響函數(shù)在感興趣頻帶內(nèi)有足夠數(shù)量的共振峰，所選結(jié)構(gòu)的激勵點和響應(yīng)點如圖4所示，桁架結(jié)構(gòu)參數(shù),如表2所示。

表2 三維桁架結(jié)構(gòu)參數(shù)表

圖4 三維桁架模型結(jié)構(gòu)圖

用表2中的模型參數(shù)作為“試驗?zāi)Ｐ汀眳?shù)，以所有桿的彈性模量和密度作為修正變量，然后對“試驗?zāi)Ｐ汀眳?shù)中修正參數(shù)偏移得到“有限元模型”。彈性模量增加10%，材料密度減小10%，得到的對應(yīng)的有限元模型參數(shù)值,如表3所示。

表3 試驗?zāi)Ｐ秃陀邢拊Ｐ蛥?shù)

采用拉丁超立方抽樣在上述有限元模型兩個參數(shù)(彈性模量、材料密度)的±20%區(qū)間內(nèi)抽取80個樣本進行模型修正。圖5為用粒子群算法進行Kriging模型相關(guān)模型參數(shù)尋優(yōu)的適應(yīng)度曲線，由圖5可知當(dāng)進化代數(shù)為31次時，適應(yīng)度曲線達到收斂，尋得最優(yōu)值。圖6為Kriging模型在第1個奇異值處得到的擬合響應(yīng)面和MSE曲面。

圖5 粒子群優(yōu)化適應(yīng)度曲線

由圖6(a)可以看出Kriging模型中樣本點都落在了所構(gòu)建的響應(yīng)面上，其與預(yù)測響應(yīng)重合度較高。而且MSE值最大也小于3×10-13，可見所構(gòu)建的代理模型精度很高。

(a)奇異值響應(yīng)面

Kriging代理模型構(gòu)造完成之后，便可以代替有限元模型進行迭代尋優(yōu)。假定試驗參數(shù)值在有限元參數(shù)值的±25%區(qū)間內(nèi)，使用布谷鳥算法進行迭代尋優(yōu)。鳥巢數(shù)為20，最大迭代次數(shù)為100。為了證明算法的穩(wěn)定性，算法運行30次，取其中的誤差最優(yōu)的一次、最差的一次和將30次取平均，迭代收斂曲線,如圖7所示。模型修正結(jié)果,如表4所示。

由圖7可以看出，布谷鳥算法尋優(yōu)能力穩(wěn)定，迭代次數(shù)為40之前就可以收斂，30次尋優(yōu)結(jié)果中最優(yōu)值和最差值差別也很小。由表4可以看出使用本文所提方法在沒有噪聲干擾的情況下達到了很好的修正效果。使用所得修正參數(shù)的修正值，Kriging模型、試驗?zāi)Ｐ图坝邢拊Ｐ皖l響函數(shù)曲線,如圖8所示。

圖7 迭代曲線

表4 模型修正參數(shù)及誤差

圖8 Kriging模型、試驗?zāi)Ｐ图坝邢拊Ｐ皖l響函數(shù)曲線

由圖8可以看出修正過后的結(jié)構(gòu)頻響函數(shù)與試驗頻響函數(shù)幾乎重合，與有限元模型頻響函數(shù)形狀相似，只是沿坐標軸有一定的移動。

為進一步檢驗本文修正方法的修正效果，比較頻響函數(shù)修正前后實、虛部曲線，如圖9所示。由圖9可以看出修正后的有限元模型無論是實部還是虛部，都能很好的與試驗?zāi)Ｐ偷念l響函數(shù)曲線重合。這進一步證明了本文所提方法的有效性。

(a)實部曲線

由于試驗測得的頻響函數(shù)不可避免的會受到噪聲的干擾。所以為了驗證本文在抗噪性方面的性能，在目標頻響函數(shù)中加入高斯白噪聲，信噪比分別為30 dB、20 dB、10 dB和5 dB，如圖10所示。

(a)信噪比30 dB

由前文理論可知，對于受到噪聲干擾時頻響函數(shù)的較小奇異值往往是噪聲信號，所以本文選擇奇異值來表征頻響函數(shù)時，應(yīng)當(dāng)選擇合適的奇異值數(shù)目，將噪聲信號盡可能的消除，保留有效的奇異值，這樣模型修正的誤差也會盡可能的減小，所以選擇合理的閾值相當(dāng)重要。根據(jù)前文所提極值點數(shù)量突變點原則來選擇合理的閾值，不同信噪比下極值點數(shù)量突變點如圖11所示，在各信噪比情況下，程序運行10次的模型修正結(jié)果均值,如表5所示。

表5 不同信噪比下修正結(jié)果

圖11為在不同信噪比下保留不同奇異值進行信號重構(gòu)以后，信號中所對應(yīng)的極值點數(shù)量。然后根據(jù)算出的極值點數(shù)量的突變點來確定有效奇異值的個數(shù)，由圖11可以看出在信噪比分別為30 dB、20 dB、10 dB和5 dB情況下選出的奇異值個數(shù)分別為55、37、20、12。圖中極值點數(shù)量在不同信噪比情況下，前面和后面變化平穩(wěn)，前面極值點數(shù)量較少，后面極值點數(shù)量多，這說明前面奇異值是受噪聲影響較小的有效奇異值，后面是受噪聲影響較大的奇異值。并且隨著信噪比的降低有效奇異值個數(shù)逐漸減小，這也符合前文所說當(dāng)信號隨著噪聲污染加大時，能表征頻響函數(shù)的有效奇異值數(shù)目逐漸減少。

(a)信噪比30 dB極值點數(shù)量突變點

由表5可以看出在相同信噪比情況下，選擇奇異值數(shù)目比不選奇異值數(shù)目修正誤差小，尤其是當(dāng)信噪比逐漸減小時，誤差尤為明顯。而且隨著噪聲干擾加大，選秩的修正誤差變化不大，而不選秩的修正誤差變化較大。當(dāng)信噪比較高時，不選秩所得修正誤差不是那么大，這表明用奇異值表征頻響函數(shù)本身就有一定的抗噪性。當(dāng)信噪比為5 dB時選秩的修正誤差也低于3%，由此也說明了本文方法對于噪聲干擾大時的有效性。

4 結(jié) 論

本文提出了一種基于頻響函數(shù)奇異值的模型修正方法，通過仿真算例證明該方法修正效果較好。結(jié)論如下：

(1)修正方法間接使用頻響函數(shù)，無需進行模態(tài)分析與振型匹配，避免了模態(tài)分析誤差，回避了傳統(tǒng)頻響函數(shù)中頻率點的選擇難題。

(2)使用頻響函數(shù)的奇異值表征頻響函數(shù)作為Kriging模型的輸出，修正精度較高。

(3)當(dāng)信號受到噪聲干擾時，使用本文方法仍能得到較滿意的修正效果。算例表明，當(dāng)信號受到弱噪聲干擾時，即使不選有效奇異值也具有一定的抗噪性；當(dāng)信號受到強噪聲干擾時，根據(jù)極值點數(shù)量突變原則選出的有效奇異值進行模型修正，修正誤差也低于3%，證明本文方法具有較強的抗噪性。