， *,

(1.大連理工大學 工程力學系 工業(yè)裝備結構分析國家重點實驗室，大連 116024；2.大連理工大學 土木工程學院，大連 116024)

1 引 言

在土木、機械和航空航天領域，長期服役或極端環(huán)境可能對結構造成損傷，這些損傷積累到一定程度便會導致結構的失效，從而造成生命和財產的重大損失。因而，采用有效的損傷識別方法對結構進行監(jiān)測和評估非常重要。

作為一種定量的全局損傷檢測方法，基于有限元模型修正的損傷識別方法在近幾十年得到廣泛的研究和應用[1]。有限元模型修正的基本思想是，結構的有限元模型與實際模型往往存在一定的偏差，可以通過修改有限元模型的物理參數(shù)，使修正后有限元模型的響應與實際結構響應盡可能一致，從而提高模型的可信度。該方法原是用作結構數(shù)值模擬分析的前置工作，以確保構建的有限元模型的準確性和有效性，其在損傷識別領域中的應用主要是通過結構的動力響應對參數(shù)進行修正，以達到全局檢測的目的，通過對比受損結構和未損結構或是受損單元與相鄰單元在剛度質量參數(shù)上的差異來判斷損傷的位置和程度，是一種同時完成損傷定位和定量的方法。

在早期的模型修正方法中，基于靈敏度的模型修正方法取得了較大成功[2]。在具體應用方面，根據所依據的不同動力響應，該方法可以分為基于頻率、振型以及頻響函數(shù)三大類別。其中，振型響應主要作為頻率響應的補充，用以應對修正參數(shù)較多的情形[3,4]。而頻響函數(shù)自身即可提供足夠多的響應信息[5,6]，但對測試環(huán)境要求較高，一定程度上限制了其應用?；陟`敏度的模型修正主要通過構建靈敏度矩陣來反映響應殘差與修正參數(shù)的關系，靈敏度矩陣通過殘差對參數(shù)的一階導數(shù)求得，進而通過迭代的方式不斷最小化殘差，直到結果收斂，以達到模型修正的目的。模型修正方法可以看作一種優(yōu)化問題，將殘差作為需要最小化的目標函數(shù)，靈敏度法可以看作無約束最速下降法求解。此外，對于優(yōu)化問題，直接使用人工智能算法，如遺傳算法(GA)[7]，也是有效求解手段，其優(yōu)點在于不需要預知參數(shù)與殘差的關系，同時可以對參數(shù)的變化范圍加以限制，避免出現(xiàn)無意義的修正結果。

模型修正方法在工程應用中的主要限制之一是計算效率問題。理論上，有限元網絡劃分得越密，用以預測結構響應就越準確，模型修正的可靠性也越高。然而，更加密集的有限元網絡將導致單次理論模型響應的計算時間急劇增多，如果模型修正程序需要足夠多的迭代次數(shù)才能收斂，其所需的計算時間將放大到不可接受。

基于此，代理模型(或稱為響應面法)作為一種快速的有限元分析替代手段，引起了眾多研究者的關注。代理模型的核心思想是通過計算機仿真模型代替有限元分析模型構建參數(shù)與響應之間的輸入輸出關系。在確保了代理模型的準確度后，于迭代過程中直接通過代理模型計算理論模型的響應輸出，以達到快速求解的目的。

目前，代理模型在有限元模型修正中的應用已經有不少成果。文獻[8，9]使用響應面法進行模型修正和損傷識別，使用特定階次的多項式擬合參數(shù)與響應的關系，通過合理的初始采樣(使用D最優(yōu)設計或中心復合設計)，用最小二乘法確認多項式的待定系數(shù)，可以有效地利用少量采樣點預測結構在不同參數(shù)下的響應。此外，為了增加響應預測點附近采樣點的影響， Chakraborty等[10]使用移動最小二乘法求解響應面多項式的系數(shù)，并證明其對于響應的預測更準確。然而，該方法也導致響應面系數(shù)變?yōu)轭A測點位置的函數(shù)，在每一次預測中都需要重新計算，增加了收斂所需的計算量。Zhou等[11]比較了使用不同的徑向基函數(shù)構建響應面的效果，并對一個大跨斜拉橋結構進行了模型修正。

此外，作為與高斯基函數(shù)相類似的一種特殊的響應面構建方式，Kriging模型在近年來受到越來越多的關注。Khodaparast等[12]使用Kriging模型建立了一個框架結構中兩個梁的位置參數(shù)與頻率的關系，用于快速進行區(qū)間模型修正，以界定參數(shù)的不確定性范圍。胡俊亮等[13]使用Kriging模型對一座大跨度鋼管混凝土連續(xù)梁拱組合體系橋進行了模型修正。Jin等[14]提出了一種順序代理建模法，在一定量的初始采樣點的基礎上，對不同響應分別進行基于特定加點準則的響應面構筑。

上述方法大多選擇了頻率作為代理模型的輸出響應，這主要是因為頻率相對于不同參數(shù)的變化較為平緩，尤其是對多項式而言，其對于參數(shù)交叉項的影響并不敏感，因此通過少量的采樣點即可得到較為精確的響應面。相較于頻率而言，頻響函數(shù)同樣是模型修正的重要組成部分之一，可以提供更多的結構振動信息，但是其在代理模型中的應用卻很少。因此，本文對使用頻響函數(shù)作為代理模型輸出響應的模型修正方法進行了研究，利用對非線性響應面模擬效果較好的Kriging模型建立其與參數(shù)的輸入輸出關系。為避免由于輸入參數(shù)增加而導致所需采樣點過多的情況，引入高效全局優(yōu)化算法(EGO)[15]來實現(xiàn)模型修正的快速收斂，并將本文方法在兩個數(shù)值模擬算例中進行應用，以驗證其用于損傷識別的可行性及有效性。

2 Kriging模型

Kriging代理模型法是一種在數(shù)學和地質學中廣泛應用的基于隨機過程的統(tǒng)計預測方法，其構建的描述輸入輸出關系的函數(shù)主要用于計算機仿真模型的快速計算。對于一組給定的m個n維樣本點輸入X={x1,x2,…,xm}T(X∈Rm×n)，和其對應的響應輸出Y={y1,y2,…,ym}T(Y∈Rm×1)，反應它們之間關系的Kriging模型的數(shù)學表達式為

(1)

式中 第一部分是數(shù)據的線性回歸，通過p個多項式模擬響應面的全局近似；第二部分為服從正態(tài)分布N(0,σ2) 的隨機過程。

需要說明的是，Kriging模型的基本假設是同樣的輸入會得到固定的輸出。因此，對于響應輸出與回歸部分之間的偏差只因模型誤差本身產生，無需考慮測量誤差等隨機因素的影響。該方法并不依賴多項式部分對響應面的模擬精度，而是著重于通過隨機過程部分的有效填補來構建合適的代理模型，因而更適合處理非線性較強的問題，其他文獻對Kriging模型的多項式部分經常取作常數(shù)。

對于隨機過程z(x)，Kriging假設所構建函數(shù)的真實響應面是連續(xù)的，任意兩點隨著距離的接近將趨向擁有相同的函數(shù)值，兩點的z(x)同理。因此，任意兩個樣本點的z(x)之間的相關性可以表示為其空間距離的函數(shù)。這里采用應用最廣泛的高斯相關模型：

(2)

(3)

在此數(shù)學模型下，樣本點真實響應出現(xiàn)的似然函數(shù)為

(4)

式中F為每個樣本點處f(x)向量值組成的矩陣。

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

(5)

(6)

在此基礎上，極大似然函數(shù)的對數(shù)形式為

(7)

通過遺傳算法求解該函數(shù)的最大值，即可確定不同維度上的衰減速率θk。

至此，一個聯(lián)系樣本點輸入輸出關系的Kriging模型構建完成。下一步則是對新點的預測。對于任意一點x0，遵循該點的預測值繼續(xù)使樣本點和預測點的似然函數(shù)極大化原則，可以得到其響應預測值為

(8)

式中rT(x0)為待求點與每個樣本點之間的相關函數(shù)向量,

rT(x0)=[R(x0,x1),…,R(x0,xn)]

(9)

值得注意的是，對于第i個樣本點響應的預測，由于rT(xi)R-1為第i階單位向量，

(10)

因此，可以認為Kriging模型是一種插值技術，對于更詳細的原理推導，可以參考文獻[15]。

3 基于Kriging的模型修正方法

傳統(tǒng)的基于靈敏度的模型修正方法是利用靈敏度矩陣將修正參數(shù)與響應的殘差聯(lián)系起來，靈敏度矩陣為響應殘差相對于修正參數(shù)的一階導數(shù)，然后通過使響應殘差為0的方式迭代求解。但是，這種方法通常需要有限元模型全部自由度處的振型或頻響函數(shù)測量值，對于不能滿足條件的測點，則需要使用理論模型縮聚或是實驗模型擴展來進行替代，在測點數(shù)遠小于模型自由度數(shù)的情況下，將由于模型誤差過大而失效。

由于基于靈敏度的模型修正方法可以看作是一種使用最速下降法求解的優(yōu)化問題，因此，可以直接使用人工智能算法求解，該優(yōu)化問題可以描述為

lbj

(11)

式中wi為響應權重，n為響應個數(shù)，resi(x)為理論響應值，res*i為實測響應值，lbj和ubj為相應修正參數(shù)的上下限。在具體應用中，所有數(shù)值均應轉換為無量綱量求解。

使用應用廣泛的遺傳算法求解上述問題，可方便地得到參數(shù)的修正結果。其優(yōu)點在于,修正參數(shù)的范圍變大，由于不需要計算靈敏度矩陣，許多靈敏度計算復雜或是無法計算的參數(shù)變得可以修正；其次是規(guī)定了參數(shù)的變化區(qū)間，可以避免參數(shù)收斂到無意義的結果上。

由于模型修正的參數(shù)間存在較強的互補性，導致模型修正的優(yōu)化是一個典型的多峰值問題，如果修正參數(shù)較多，算法很容易收斂到一個局部解上。此外，該方法的核心問題在于，需要進行大量理論響應的計算來完成全局搜索。因此，利用Kriging模型代替原始有限元模型進行響應計算，可以在保留此方法優(yōu)點的同時，顯著提升其計算效率。

對于頻響函數(shù)，為了構建平滑的響應面，將測點頻響函數(shù)模值的對數(shù)作為輸出響應，同時應避免使用參數(shù)變化空間內產生的共振頻率作為頻率點。但頻響函數(shù)的響應面較為復雜，無法僅憑多項式部分進行精確擬合，因此只能通過足夠的采樣點來提升隨機過程的預測精度。而足夠精度的響應面所需的采樣點個數(shù)將隨修正參數(shù)數(shù)目的增多而快速增加至不可接受的程度。因而，直接構建頻響函數(shù)響應面只適用于少量參數(shù)的情形。

由于只有優(yōu)化問題的結果值得關注，而遠離最優(yōu)解區(qū)域的代理模型準確程度并不重要，因此，本文將優(yōu)化問題的目標函數(shù)，即理論頻響與實際頻響殘差的加權和f(x)直接作為響應輸出，構建修正參數(shù)與該響應輸出之間的Kriging模型，再使用高效全局優(yōu)化算法EGO[15]來實現(xiàn)快速建模求解f(x)的最小值，其對應的輸入參數(shù)即為參數(shù)的修正結果。

對于少量初始樣本點構建的Kriging模型，可以通過增加樣本點來提高其對于遠離已有樣本點位置的未知點的預測精度。由于本文的優(yōu)化問題只關心最小值的位置，對于一定初始樣本點構建的Kriging模型，可以僅在可能產生最小值的位置進行加點，即只提升有用區(qū)域的模型精確度，達到只使用少量采樣點求解優(yōu)化問題的目的。在具體的加點過程中，使用EGO算法具有較好效果，可以很好地平衡局部和全局搜索，避免優(yōu)化結果收斂到局部極小值上，其原理如下。

對于Kriging模型，除了可以對未知點的輸出進行預測，還可以對該預測值的準確程度進行評估。

(12)

在得到了未知點的預測值和均方差之后，EGO算法將該點的取值理解為一個均值為預測值、方差為均方差的高斯概率密度函數(shù)。該算法選擇的加點位置并不是單純的Kriging模型的最小值，而是計算Kriging模型中引入未知點對于已有樣本點的響應輸出最小值的提升的預期，并尋找該預期最大值的位置。

I(x)=ymin-y(x)

(13)

式中 erf(·)為累積分布函數(shù)，ymin為樣本點響應的最小值，I(x)為未知點的加入對最小值的提升，E[I(x)]為該提升的期望，通常未知點的預測值越小、方差越大，則E[I(x)]的值越大。因此EGO算法不僅會搜索Kriging模型中響應預測值最小的點，即求解局部極小值，還會加入響應預測值較小且預測準確度較低(均方差較大)的點，即進行全局搜索。本文對于最大期望所在位置的求解主要通過遺傳算法來完成。

基于Kriging模型用于頻響函數(shù)的模型修正的主要步驟可以總結如下。

(1) 根據拉丁超立方采樣，生成一定數(shù)量的關于修正參數(shù)的初始樣本點集。

(2) 運行有限元模擬程序，計算生成樣本點的目標函數(shù)輸出向量，并建立初始的Kriging模型。

(3) 使用遺傳算法，計算當前Kriging模型中改進期望最大值所在的位置，并將其作為新加入的樣本點。

(4) 計算新樣本點的真實響應輸出，并將其加入已有的樣本點中，更新Kriging模型。

(5) 檢驗是否達到終止判據，達到則停止，并將輸出向量中的最小值所對應的樣本點作為修正結果；否則，回到步驟(3)，繼續(xù)加點運算。

Kriging模型的加點終止準則很多，本文選擇改進期望的最小值作為判定依據，即Kriging模型將不斷向未知區(qū)域進行探索，直到所有位置的對當前最小值的最大改進期望低于一個值。

基于Kriging模型的模型修正不僅繼承了遺傳算法不需要計算敏感度矩陣以及布置大量測點的特點，其優(yōu)勢還在于，首先算法會率先尋找當前代理模型下的最小值；然后，會受迫向可能產生更小值的區(qū)域進行探索，降低算法收斂到局部極小值的可能，且算法求解的精確程度易控。最重要的是，其對于有限元模型計算次數(shù)的需求遠少于直接使用遺傳算法。

4 數(shù)值仿真算例

4.1 算例1

首先通過一個3自由度的質量-彈簧系統(tǒng)來驗證本文方法，如圖1所示，理論模型中各參數(shù)的取值為

m1=1 kg,m2=3 kg, m3=1 kg

k1=k2=k3=k4=k5=10 N/m,k6=30 N/m

假設實際系統(tǒng)中存在損傷，實驗模型與理論模型的區(qū)別在于k2=8 N/m,k4=9 N/m，首先對兩者進行特征值分析，得到其前三階的頻率,列入表1。接著以k2和k4的系數(shù)作為修正參數(shù)，其變化區(qū)間設為[0.5,1.5]，以頻響函數(shù)矩陣的H(1,3)分量作為參考響應，激振頻率分別取0.5 Hz,1 Hz和1.5 Hz，其數(shù)目略多于修正參數(shù)的數(shù)量，以便驗證兩種方法在實驗數(shù)據有限時的修正結果。需要注意的是，頻響函數(shù)的激振頻率應避免出現(xiàn)在自然頻率隨修正參數(shù)改變的變化范圍內，否則，將會造成響應面的突變，增加代理模型的構建難度。

表1 3自由度質量-彈簧系統(tǒng)理論和實驗模型的自然頻率(單位：Hz)

Tab.1 Natural frequencies of the analytical and experimental model of the three degree of freedom mass-spring system(unit:Hz)

模態(tài)階數(shù)理論模型實驗模型10.37070.351720.78900.773331.42351.4191

圖1 3自由度質量-彈簧系統(tǒng)

Fig.1 Three degree of freedom mass-spring system

運用Kriging模型修正方法對參數(shù)進行修正求解，響應的權重均取為1，同時使用基于遺傳算法的模型修正作為參考，對比兩種方法修正結果的優(yōu)劣。本文所有用到的遺傳算法均使用Matlab自帶的GA函數(shù)，沒有加入適用于模型修正的改進，以驗證Kriging模型修正在原始遺傳算法上的效果。

在修正參數(shù)的變化區(qū)間里，采用拉丁超立方采樣獲取20個樣本點作為初始樣本。依次進行Kriging模型的修正步驟，終止標準設定為改進期望的最大值E[I(x)]<10-10時停止。僅經過8次迭代，Kriging模型便構建完成，如圖2所示。

由于所需修正參數(shù)數(shù)量較少，僅適用于少量樣本點即可構建足夠精度的響應面。本算例中，Kriging模型構建的響應面與結構真實的目標函數(shù)響應面幾乎完全一致。兩種方法的修正結果列入表2，可以看出，兩種方法均可得到較準確的修正結果，而Kriging模型的結果精確度較高。

對于這種參數(shù)較少，可以直接構建出精確響應面的問題，還可以直接將頻響函數(shù)作為輸出，分別構建每個頻響函數(shù)的響應面，再通過遺傳算法進行求解，同樣可以大幅減少理論模型本身的計算次數(shù)，所構建響應面還具有可以重復利用的優(yōu)點。然而，隨著修正參數(shù)數(shù)量的增多，構建精確響應面所需樣本點的個數(shù)也急劇增多。樣本點個數(shù)是制約Kriging模型計算效率的關鍵，因為對每一個新點的預測都需要計算該點與所有樣本點的相關性。因此，為保證計算效率，Kriging模型的樣本點應盡量控制在500以下。由于頻響函數(shù)隨參數(shù)變化復雜，直接構建精確響應面可接受的修正參數(shù)數(shù)量最多不超過5，否則就必須依賴EGO算法來進行求解。

表2 3自由度系統(tǒng)的模型修正結果對比

Tab.2 Comparison of model updating results for 3 degree of freedom system

修正參數(shù)實際值Kriging誤差/%GA誤差/%k20.80.83534.410.85807.25k40.90.89750.280.89340.73

圖2 3自由度系統(tǒng)的Kriging響應面

Fig.2 Kriging response surface of 3 degree of freedom system

4.2 算例2

為進一步驗證本文方法在修正參數(shù)較多時的應用效果，使用一個懸臂梁有限元模型來進行修正。該模型長為0.9 m，截面慣性矩為9.135×10-10m4，材料的彈性模量為200 GPa，密度為7670 kg/m3，將其劃分為9個單元，如圖3所示。

本算例中，梁單元采用伯努利-歐拉假定，每個節(jié)點上有橫向位移和轉角兩個自由度，假定結構存在損傷，造成了單元剛度的縮減，單元2減少20%，單元4減少30%，單元5減少10%。理論模型和損傷模型的自然頻率列入表3。

以所有單元的剛度系數(shù)作為修正參數(shù)，其變化區(qū)間設為[0.7,1]。由于本算例要修正的參數(shù)多達9個，且參數(shù)之間的相互補償狀況較強，所以將參數(shù)的變化區(qū)間適當縮小，以降低響應面的構建成本，提高算法的收斂效率。

參考響應選取為頻響函數(shù)矩陣的H(3,17)分量，激振頻率分別取為10 Hz,45 Hz,115 Hz,220 Hz,360 Hz,600 Hz,760 Hz,1010 Hz和1300 Hz。這些激振頻率基本都選在理論自然頻率稍微靠右的位置，首先是因為這些位置的頻響函數(shù)的變化較為劇烈，可以有效避免實際測試中模型誤差和噪音誤差的影響；其次構建響應面時，勢必會考慮最極端的情況，即所有單元都達到理論的最大損傷，由于剛度降低，此時模型的自然頻率較初始模型會產生較大的左移，為了避免激振頻率與自然頻率重合，必須將激振頻率選定在初始模型兩個自然頻率之間靠右的位置。

表3 懸臂梁理論和損傷模型的自然頻率(單位：Hz)

Tab.3 Natural frequencies of the analytical and damaged model of the cantilever(unit:Hz)

模態(tài)階數(shù)理論模型實驗模型16.115.79238.2936.663107.26104.54210.41203.495348.63330.786522.95504.747734.98709.418985.66949.1491261.51210.95

圖3 懸臂梁有限元模型

Fig.3 Finite element model of a cantilever

采用拉丁超立方采樣獲取90個樣本作為初始樣本，進行Kriging模型修正，終止標準不變。本算例中，Kriging模型經過179次迭代構建完成。由于參數(shù)的維數(shù)較高，僅將單元1和單元2的剛度系數(shù)作為變量，其他單元的剛度系數(shù)設為1不變，Kriging 模型預測的響應面與真實響應面的比較如圖4所示?？梢钥闯?，當修正參數(shù)較多時，Kriging模型對于遠離最小值附近的目標函數(shù)的預測誤差非常大。但是，這并不影響Kriging模型給出滿意的修正結果。兩種方法的修正結果對比列入 表4?？梢钥闯?，Kriging方法修正的效果較好，而GA法則出現(xiàn)了明顯的誤判。

此外，本算例還比較了兩種方法所需的理論模型有限元程序計算次數(shù)。分別使用前6,7,8,9個單元的剛度系數(shù)作為修正參數(shù)，兩種方法所需的有限元程序計算次數(shù)列入表5。可以看出，基于Kriging模型的有限元模型修正需要的計算次數(shù)遠比直接使用GA方法少得多。隨著結構有限元模型的精細化和復雜化，直接使用GA方法所耗費的時間將變得不可接受，這也正是使用代理模型進行模型修正的意義所在。

表4 懸臂梁的模型修正結果對比

Tab.4 Comparison of model updating results for the cantilever

修正參數(shù)實際值Kriging誤差/%GA誤差/%E110.99570.430.99540.46E20.80.79540.580.78202.30E310.99940.060.98911.0E40.70.70991.390.810213.6E50.90.92182.360.92512.71E610.99900.10.904110.61E710.99080.9310E810.99140.870.96453.68E910.98931.080.96613.51

表5 兩種方法所需的有限元模型計算次數(shù)

Tab.5 Number of finite element models required for the two methods

參數(shù)個數(shù)KrigingGA612810800713711000828811400926912000

圖4 Kriging模型響應面與真實響應面的對比

Fig.4 Comparison between Kriging model response surface and real response surface

5 結 論

本文提出一種新的使用頻響函數(shù)作為測量響應的模型修正方法。該方法主要通過構建Kriging代理模型來建立修正參數(shù)輸入與響應殘差之和的目標函數(shù)輸出之間的聯(lián)系，再通過全局優(yōu)化算法EGO快速求解。該方法主要具有以下優(yōu)點。

(1) 不需要已知修正參數(shù)與響應殘差之間的關系，對于測點的位置分布同樣要求較低，在算例中，只需要頻響函數(shù)矩陣的一個分量即可實現(xiàn)修正目的，其應用的可行性和方便性優(yōu)于傳統(tǒng)的靈敏度方法。

(2) 可以通過控制終止標準來提高算法搜索的廣度和深度，即便是使用未加處理的遺傳算法，也可以顯著提升修正結果的精確度。

(3) 由于使用了代理模型代替有限元模型進行優(yōu)化問題求解，對有限元模型的計算需求大大降低，對于單次計算時間較長的有限元模型，計算效率顯著提高。