基于頻響函數奇異值的模型修正方法

曹明明，彭珍瑞，劉滿東

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

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

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

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

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

1 頻響函數特征量選取

1.1 頻響函數理論

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

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

(1)

式中:M、K和C分別表示質量矩陣、剛度矩陣和阻尼矩陣；ω為激勵頻率；F(ω)為簡諧激勵；穩態響應X(ω)可表示為：

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

(2)

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

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

(3)

1.2 奇異值分解理論

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

(4)

式中：n=N-(m-1)τ；A為m×n階的Hankel矩陣；τ為延遲步長，一般取τ=1；m為嵌入維數，當N是偶數時，m=N/2，當N是奇數時，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)

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

1.3 極值點數量突變原則

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

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

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

2.1 Kriging模型的建立

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

(8)

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

(9)

根據極大似然函數法，可以求得：

(10)

(11)

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

(12)

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

淺談道路橋梁施工的常見問題及質量檢測技術的應用…………………………………………………… 田文澤（11-100）

2.2 粒子群算法優化Kriging模型

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

2.3 模型修正流程

當Kriging模型構造完成時，模型修正歸根結底是一優化問題，在目標函數最小的條件下，通過求解該優化問題，得到設計參數的修正值。目標函數為：

(13)

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

圖1 模型修正流程圖

3 數值算例

3.1 算例一

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

圖2 車輛三自由度垂向頻率響應模型

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

系統運動方程為：

(14)

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

表1 車輛三自由度系統的模型修正結果

圖3 車輛三自由度系統的Kriging響應面

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

3.2 算例二

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

表2 三維桁架結構參數表

圖4 三維桁架模型結構圖

用表2中的模型參數作為“試驗模型”參數，以所有桿的彈性模量和密度作為修正變量，然后對“試驗模型”參數中修正參數偏移得到“有限元模型”。彈性模量增加10%，材料密度減小10%，得到的對應的有限元模型參數值,如表3所示。

表3 試驗模型和有限元模型參數

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

圖5 粒子群優化適應度曲線

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

(a)奇異值響應面

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

由圖7可以看出，布谷鳥算法尋優能力穩定，迭代次數為40之前就可以收斂，30次尋優結果中最優值和最差值差別也很小。由表4可以看出使用本文所提方法在沒有噪聲干擾的情況下達到了很好的修正效果。使用所得修正參數的修正值，Kriging模型、試驗模型及有限元模型頻響函數曲線,如圖8所示。

圖7 迭代曲線

表4 模型修正參數及誤差

圖8 Kriging模型、試驗模型及有限元模型頻響函數曲線

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

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

(a)實部曲線

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

(a)信噪比30 dB

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

表5 不同信噪比下修正結果

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

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

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

4 結 論

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

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

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

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