張 欣，劉 洋，高丹盈（鄭州大學土木工程學院 鄭州，450001）

基于第2代小波的有限元模型更新方法*

張 欣，劉 洋，高丹盈
（鄭州大學土木工程學院 鄭州，450001）

提出了有限元模型更新剛度曲線的低分辨率表示方法?；谟邢薜膶嶒災B(tài)信息，通過遺傳算法得到了更新剛度曲線在低分辨率表示時的尺度函數(shù)系數(shù)和小波系數(shù)。利用多分辨率分析方法實現(xiàn)了復雜精細有限元模型的幾何降維，獲得了在有限模態(tài)參數(shù)范圍內(nèi)與精細模型相似的簡單有限元模型，并在此簡單模型基礎(chǔ)上實現(xiàn)了最終的模型更新與損傷識別。證明了從多分辨率分析的角度出發(fā)，利用第2代小波變換的方法能夠減少待更新的參數(shù)數(shù)量，降低模型更新過程的奇異性。以變截面箱梁為例，驗證了提出方法對裂縫深度變化和局部裂縫數(shù)量變化的魯棒性。本研究方法也適用于多組裂縫群的辨識。

模型更新；損傷識別；多分辨率分析；小波分析

引 言

有限元模型更新與結(jié)構(gòu)損傷識別均為反分析問題，在實際工程中有廣泛的應(yīng)用。目前，使用較為頻繁的方法是通過非線性優(yōu)化［1-7］或者回歸分析［8］等數(shù)值方法尋求適合的模型參數(shù)向量，用有限元模型的理論模態(tài)擬合結(jié)構(gòu)的實驗模態(tài)，使兩者間的誤差最小，從而得到能夠反映結(jié)構(gòu)實際情況的有限元模型。此類問題具有非常強烈的奇異性［9］：a.模型更新或者損傷識別的結(jié)果對所選定的待更新參數(shù)集十分敏感；b.非線性優(yōu)化過程具有極大的不確定性，迭代過程可能會收斂于局部最小值而不是全局最小值，更新結(jié)果未必具有實際意義；c.對于比較復雜的結(jié)構(gòu)或者損傷模式，需要建立精細的有限元模型才能正確模擬，所需單元數(shù)量眾多，需要更新的參數(shù)數(shù)量也相應(yīng)增加，需要運用大規(guī)模非線性優(yōu)化算法進行分析。此時，計算結(jié)果的穩(wěn)定性更低，幾乎無法得到具有實際意義的結(jié)果；d.實驗誤差對更新結(jié)果的影響也不可忽視。因此，有效減少需要更新的參數(shù)數(shù)量［10-12］是提高模型更新和結(jié)構(gòu)損傷識別效率和可靠度的合理技術(shù)路線。

以剛度更新為例，有限元模型的實際剛度即反映結(jié)構(gòu)真實狀況的剛度，可以被視為兩種成分的疊加：一種成分為理論假設(shè)剛度；另一種成分為實驗更新剛度。前者的精細度可以任意假定，用于建立高精度的有限元模型，仿真復雜的工程結(jié)構(gòu)。后者由實驗確定，精細度無法提高。因此，模型更新的數(shù)學實質(zhì)是在一條高精度（高分辨率）理論剛度曲線上疊加低精度（低分辨率）更新剛度曲線，得到結(jié)構(gòu)的模態(tài)參數(shù)來擬合實驗數(shù)據(jù)。核心問題是如何使用最少的數(shù)據(jù)來近似更新剛度曲線。這屬于數(shù)據(jù)壓縮技術(shù)。數(shù)據(jù)壓縮技術(shù)屬于多分辨率分析的范疇，在影音圖像處理領(lǐng)域有很多的成熟應(yīng)用，使用頻率較高的應(yīng)屬小波變換。

筆者首先介紹多分辨率分析與小波變換的相關(guān)背景；然后闡述低分辨率模型更新思想。

1 多分辨率分析與小波變換

1.1 多分辨率分析理論

具有有限能量的一個函數(shù)，如結(jié)構(gòu)剛度函數(shù)κ（x），屬于內(nèi)積空間L2（R），該空間為無窮維函數(shù)線性空間，存在無窮個線性無關(guān)向量φk（x｛｝）∞構(gòu)成的基函數(shù)族，k為采樣點的位置［13］。則

在該內(nèi)積空間內(nèi)存在嵌套式閉子空間逼近序列

其中：尺度函數(shù)φ（x-k）為Riesz基，滿足雙尺度方程

這樣構(gòu)成多分辨率逼近。計Wj=Vj+1/Vj為Vj在Vj+1中的補子空間，則Vj+1=Wj⊕Vj，⊕為子空間直和。依次類推，有

ψ（x）為W0的基函數(shù)，稱小波函數(shù)。W0∈V1，可由φ（x）線性表示

其中：φi，k（x）為低通濾波函數(shù)；ψj，k（x）為帶通濾波函數(shù)；｛hn｝為低通數(shù)字濾波器；｛gn｝為高通數(shù)字濾波器。

在j尺度層面，函數(shù)κ（x）可表示為

其中：~φj，k（x）為φj，k（x）的對偶尺度函數(shù)。

其中：~ψj，k（x）為ψj，k（x）的對偶小波函數(shù)。

這樣就形成κ（x）在j尺度上的近似κj（x）。j數(shù)值越大，就越能精確表述原函數(shù)，但是所需數(shù)據(jù)就越多；j數(shù)值越小，就會有越多的原函數(shù)信息被忽略，κj（x）就越粗糙，變化趨勢也越平緩，所需數(shù)據(jù)點數(shù)也越少?；谝陨侠碚摼涂梢詫崿F(xiàn)更新剛度曲線的低分辨率表示。然而，小波變換在無限域進行，當應(yīng)用于有限元模型更新時，必須考慮結(jié)構(gòu)只在有限范圍內(nèi)存在這一問題，需要考慮小波函數(shù)與尺度函數(shù)的邊界效應(yīng)。在這一點上，第2代小波［14］有優(yōu)勢，它具有與第1代小波相似的特性且構(gòu)造簡單方便，一般通過提升格式［15］就可以實現(xiàn)有限域上非均勻采樣間隔的小波函數(shù)和相應(yīng)尺度函數(shù)的構(gòu)造。

1.2 提升格式與第2代小波函數(shù)和尺度函數(shù)的構(gòu)造

提升格式［15］指在簡單的小波函數(shù)與尺度函數(shù)的基礎(chǔ)上，通過提升運算，層層加碼，形成復雜形式的小波變換。

構(gòu)造第2代內(nèi)插型尺度函數(shù)與小波函數(shù)一般可以選用懶惰小波為提升起點，經(jīng)過對偶提升（預測）實現(xiàn)合成尺度函數(shù)與分解小波的提升，再通過提升過程實現(xiàn)分解尺度函數(shù)與合成小波函數(shù)的提升。

定義如下多尺度分析算子：從Vj+1到Vj的運算為~H運算，Vj+1到Wj的運算為~G運算，逆運算分別為H運算和G運算。設(shè)E為偶下采樣算子，D為奇下采樣算子，懶惰小波運算為

按下式進行提升

內(nèi)插濾波運算Hintj可以采用兩點、三點、四點或更多點數(shù)的內(nèi)插濾波形式，由此可以產(chǎn)生不同形式的尺度函數(shù)與小波函數(shù)。

在以上提升格式的基礎(chǔ)上，通過迭代算法獲得每個尺度相應(yīng)采樣點上的尺度函數(shù)與小波函數(shù)。因為第2代小波能適應(yīng)非均勻采樣以及邊界的影響，在某一個尺度層內(nèi)小波函數(shù)和尺度函數(shù)可能會因點而異，各尺度層間也可能互異。

2 低分辨率模型更新

2.1 基本方法

設(shè)更新目標為結(jié)構(gòu)的剛度曲線K（x）

其中：K（x）為結(jié)構(gòu)在x處的理論剛度；κ（x）為更新剛度。

結(jié)構(gòu)的質(zhì)量分布不更新。根據(jù)多分辨率分析理論，在j尺度上，模型的更新剛度函數(shù)為

在j+1尺度上為

當前模型結(jié)構(gòu)的模態(tài)參量為（λi，φi），i=1，2，…為剛度曲線的函數(shù)。其中：λi為第i階頻率；φi為為第i階振型；實驗模態(tài)參數(shù)為（ˉλi，ˉφi），i=1，2，…。

定義誤差向量

其中：權(quán)重矩陣為

在優(yōu)化求解過程中可以強調(diào)某一階模態(tài)信息的重要性，也可以降低該信息的重要性。這可以通過調(diào)整權(quán)重矩陣對角元素的大小來實現(xiàn)。第i階振型的相關(guān)因數(shù)為

其中：ˉφi，φi分別為目標和當前振型；*代表共軛轉(zhuǎn)置。模型修正的目標是使r（α）的模最小

2.2 兩階段更新模式

r（α）的?？赡懿粸榱?，因為即使當r（α）的維數(shù)，即方程個數(shù)大于或等于修正向量α的維數(shù)時，非線性方程r（α）=0也可能無解。因此模型更新及損傷識別過程為參數(shù)優(yōu)化過程，而非定解方程求解問題。

一般希望模型更新能在簡單模型上進行。例如，使用梁單元建立圖1，2中結(jié)構(gòu)的有限元模型，在此模型基礎(chǔ)上進行更新操作。由于計算模型不能夠完全擬合原始結(jié)構(gòu)的行為，因此存在模型初始誤差。另外，在結(jié)構(gòu)出現(xiàn)損傷之后，由于圣維南現(xiàn)象的影響，結(jié)構(gòu)在損傷附近的力學行為更趨復雜，用梁單元不能精確描述，這樣就進一步產(chǎn)生了模型更新誤差。此外，必然存在的實驗誤差也會使方程r（α）=0無解。此問題不能通過加入高階模態(tài)信息來解決。例如，對于四點內(nèi)插濾波，當有兩階模態(tài)信息時，理論上可以實現(xiàn)0尺度分辨率的定解更新，當有四階模態(tài)信息時，可以實現(xiàn)1尺度分辨率的定解更新。依次類推。理論上，當有足夠模態(tài)信息時，可以在足夠精細的尺度層面上實現(xiàn)定解更新。然而，隨著模態(tài)階數(shù)的不斷升高，理想梁模型的行為就越發(fā)背離實際結(jié)構(gòu)，由此將產(chǎn)生更大的模型誤差。

圖1 算例結(jié)構(gòu)簡圖（單位：mm）Fig.1 Sketch of the structure under study（unit：mm）

筆者建議根據(jù)結(jié)構(gòu)自身的特性，采取適當階數(shù)的模態(tài)信息，分以下兩個階段進行模型更新和損傷識別。

1）模型的初始校準階段。以圖1，2所示結(jié)構(gòu)為例。首先，建立結(jié)構(gòu)的精細有限元模型；然后，根據(jù)結(jié)構(gòu)幾何特征計算等效梁的剛度曲線建立梁單元模型。此時簡化模型的模態(tài)參數(shù)與精細模型的模態(tài)參數(shù)不同。依照本研究模型更新方法，得到低分辨率更新剛度曲線，疊加在幾何特征剛度曲線上，作為模型初始校準后的剛度曲線。

圖2 結(jié)構(gòu)兩端橫斷面（單位：mm）Fig.2 Cross section of two ends（unit：mm）

2）模型的損傷更新階段。如圖1所示，在模型中引入結(jié)構(gòu)損傷，此時模型的模態(tài)參數(shù)發(fā)生變化。在前一階段剛度曲線的基礎(chǔ)上，重復低分辨率更新過程，擬合有損傷模型的模態(tài)參數(shù)，獲得表征結(jié)構(gòu)損傷特性的尺度函數(shù)與小波函數(shù)系數(shù)以及相應(yīng)的損傷剛度更新曲線。

2.3 近似誤差

由于更新剛度曲線是在較低分辨率層面上實現(xiàn)的，因此可能會出現(xiàn)re（αe）殘差較大的情況。只有在較高分辨率層上實現(xiàn)更高分辨率的更新剛度曲線，才能從根本上解決這一問題。然而，更新曲線的分辨率由實驗結(jié)果決定，因此在一般情況下不太可能得到太高分辨率的更新曲線。這也是反分析問題的根本癥結(jié)之一：根據(jù)低分辨率實驗獲得高分辨率反分析結(jié)果，實現(xiàn)所謂的“超分辨率”表述是困難的。

另外，在更新剛度曲線形狀比較復雜而曲線實現(xiàn)的分辨率層又比較低時，會出現(xiàn)更新曲線線形失真的現(xiàn)象。這種現(xiàn)象與信號處理中采用較低的采樣頻率記錄較高頻率信號而產(chǎn)生的信號失真現(xiàn)象類似，是一種客觀現(xiàn)實，不能逾越，其后果不易預測。因為更新剛度曲線的低尺度表示是嵌套在結(jié)構(gòu)模態(tài)參數(shù)求解運算這一非線性運算過程中的，此時誤差討論極為困難，因此應(yīng)校驗臨近尺度層參數(shù)更新結(jié)果的收斂性，以保證更新結(jié)果的正確性。

3 數(shù)值算例

在損傷識別的過程中，如果更新過程在復雜、精細的有限元模型上進行，則顯得不夠經(jīng)濟有效，所以一般希望能在簡單有限元模型上進行。但是，對于具有復雜結(jié)構(gòu)形式或損傷方式的結(jié)構(gòu)，只有建立復雜精細的有限元模型，才能較好地仿真，因此產(chǎn)生矛盾。解決方案是先建立精細的有限元模型，獲得理論模型模態(tài)，再建立簡化模型，通過低分辨率更新方法將簡化模型的剛度進行更新，得到在有限模態(tài)范圍內(nèi)能基本反映精細模型的簡單有限元模型，在此基礎(chǔ)上再以實驗模態(tài)為目標更新并得到損傷信息。

3.1 復雜有限元模型的降維近似

以圖1，2所示的變截面箱梁為例，彈性模量為3×1010Pa，泊松比為0.2，這里僅考慮豎向振動模態(tài)。以板單元建立三維精細參照模型，單元總數(shù)為7 500。在固定端處釋放上部部分節(jié)點，仿真結(jié)構(gòu)豎向裂縫。采用ANSYS計算得到兩種工況模態(tài)參數(shù)如表1所示。梁單元每個節(jié)點僅考慮豎向和扭轉(zhuǎn)兩個自由度。單元總數(shù)為128，每個單元慣性矩由截面幾何尺寸計算得出。

表1 模型豎向振型模態(tài)參數(shù)Tab.1 Modal parameters of the vertical modes

計算數(shù)據(jù)表明，簡化模型與精細模型的模態(tài)參數(shù)存在明顯差別，必須對簡化模型的單元剛度進行更新。更新目標為無損精細有限元模型。

為考慮模型剛度的多分辨率表示，設(shè)梁單元內(nèi)力虛功為

取形函數(shù)為q（x），節(jié)點位移為u，將K（~x）的多分辨率表達式帶入得到單元剛度矩陣元素為

依照式（22）在0尺度分辨率和1尺度分辨率上分別建立有限元模型，計算模態(tài)參數(shù)更新目標函數(shù)r（α）。通過遺傳算法求得更新參數(shù)，使式（20）最小化。

3.2 模型校準及單裂縫工況計算結(jié)果

令djk=0，j=0，得到0尺度層尺度函數(shù)系數(shù)。令cjk，djk≠0，j=0，得到0尺度層尺度函數(shù)系數(shù)與小波函數(shù)系數(shù)，通過式（16）構(gòu)成1尺度分辨率層剛度函數(shù)。這兩種更新的結(jié)果不盡相同，但基本類似。圖3中更新階段A為一尺度分辨率模型初始校準階段結(jié)果。橫坐標1～5為尺度函數(shù)系數(shù)，6～9為小波函數(shù)系數(shù)。

圖3 模型更新尺度函數(shù)與小波函數(shù)系數(shù)Fig.3 Updating coefficients of the scale function and wavelet function

由小波函數(shù)和尺度函數(shù)系數(shù)合成的更新剛度曲線如圖4所示。圖4結(jié)果表明，三者吻合程度較好，0尺度結(jié)果基本收斂于1尺度結(jié)果。

在經(jīng)過初始校準的模型上，針對有損模型的更新結(jié)果也繪制在圖3，4中，用更新階段B表示。單元剛度折減曲線如圖5所示。

需要說明的是，筆者采用連續(xù)函數(shù)來表示結(jié)構(gòu)剛度的方法與目前流行的做法不同。目前流行的局部損傷處理方法是采用集中剛度折減法，即將受損單元的剛度降低而相鄰單元的剛度保持不變，相當于用方波函數(shù)來表示結(jié)構(gòu)剛度的空間變化。這樣的做法在單元劃分比較粗糙時適用。如果需要精細劃分單元則會出現(xiàn)問題，因為從圖形分辨率角度來說，由于存在階躍，理想的方波函數(shù)相當于無限分辨率圖形，在單元精細劃分的情況下需要很多實驗信息量才能有效識別。基于有限實驗信息，采用低分辨率函數(shù)表述剛度折減，就如同對方波函數(shù)采取小波分解，只保留緩慢變化的那一部分低分辨率趨勢線，從而降低了對實驗信息量的要求，但是必然會有信息擴散效應(yīng)。如圖5所示，剛度折減已經(jīng)擴散到較多單元，這正是信息擴散效應(yīng)所致。結(jié)構(gòu)損傷的大致位置可以通過剛度折減曲線的谷值位置來判斷。

圖4 梁模型更新結(jié)果Fig.4 Updating results of the beam

圖5 損傷引起的剛度折減Fig.5 Reduction of stiffness

3.3 多裂縫工況計算結(jié)果

多裂縫工況由圖1中Cr1～Cr4裂縫組合而成。具體組合及相應(yīng)的前3階豎向頻率如表2所示。表中Y代表采用該項裂縫組合。

表2 損傷狀況組合Tab.2 Combination of damages

利用上述方法得到3種損傷組合工況下剛度折減曲線，如圖6所示?？梢钥闯觯篴.因為圖6中工況1的裂縫預設(shè)深度大于圖5，所以圖6工況1的剛度折減幅度比圖5大，但是折減函數(shù)曲線的形狀與圖5大致相同，說明本研究方法對局部裂縫深度變化具有魯棒性；b.工況2與工況1相比，由于存在較多的裂縫，因此產(chǎn)生了較大的剛度折減，但是折減函數(shù)曲線的形狀與工況1大致相同，說明本研究方法對局部裂縫數(shù)量變化具有魯棒性；c.工況3與工況2在梁的左端具有相同的損傷組合，圖6表明，兩種工況在此處的剛度折減也相仿，說明本研究方法對裂縫群的數(shù)量變化具有魯棒性；d.工況3在80號單元附近有正彎矩區(qū)裂縫，表明在該位置有明顯的剛度折減，說明本研究方法適用于多組裂縫群的識別。

圖6 組合損傷工況下剛度折減Fig.6 Stiffness reduction of various damage conditions

4 結(jié)束語

基于多分辨率分析的有限元模型更新與損傷識別方法可以在保持結(jié)構(gòu)更新剛度曲線基本特征的前提下，有效減少需要更新的模型參數(shù)，降低了非線性優(yōu)化過程的復雜程度和不確定性。第2代小波是有效的多分辨率分析工具，可應(yīng)用于有限尺寸結(jié)構(gòu)，符合實際應(yīng)用條件下結(jié)構(gòu)損傷識別的要求。

在本研究方法的框架內(nèi)，結(jié)構(gòu)的損傷識別過程可以分解為兩個步驟來進行：a.復雜精細有限元模型的幾何降維，獲得簡化模型；b.在簡化模型的基礎(chǔ)上實現(xiàn)結(jié)構(gòu)損傷識別。以上步驟均可通過多分辨率模型更新方法結(jié)合遺傳算法實現(xiàn)。該方法具有對裂縫深度變化和局部裂縫數(shù)量變化的魯棒性，并能夠適應(yīng)多組裂縫群的識別。

［1］ Zapico-Valle J L，Alonso-Camblor R，González-Martínez M P，et al.A new method for finite elementmodel updating in structural dynamics［J］.Mechanical Systems and Signal Processing，2010，24（7）：2137-2159.

［2］ Friswell M I，Mottershead J E.Finite element model updating in structural dynamics［M］.Dordrecht：Kluwer Academic Publishers，1999：1-20.

［3］ Teughels A，DeRoeck G.Structural damage identification of the highway bridge Z24 by FE model updating ［J］.Engineering Structures，2003，278：589-610.

［4］ 侯吉林，歐進萍.基于局部脈沖響應(yīng)的約束子結(jié)構(gòu)修正法［J］.工程力學，2009，26（11）：23-30.Hou Jilin，Ou Jinping.Isolated substructure model updating method based on local impulse response［J］.Engineering Mechanics，2009，26（11）：23-30.（in Chinese）

［5］ 郭力，李兆霞，陳鴻天.基于子結(jié)構(gòu)分析的多重子步模型修正方法［J］.中國工程科學，2006，8（9）：42-48.Guo Li，Li Zhaoxia，Chen Hongtian.Multi-stage model updating method via substructure analysis［J］.Engineering Science，2006，8（9）：42-48.（in Chinese）

［6］ 張保強，陳國平，郭勤濤.使用有效模態(tài)質(zhì)量和遺傳算法的有限元模型修正［J］.振動、測試與診斷，2012，32（4）：577-580.Zhang Baoqiang，Chen Guoping，Guo Qintao.Finite element model updating using effective modal mass with genetic algorithm［J］.Journal of Vibration，Measurement＆amp；Diagnosis，2012，32（4）：577-580.（in Chinese）

［7］ Koh C G，Chen Y F，Liaw C Y.A hybrid computational strategy for identification of structural parameters［J］.Computers and Structures，2003，81：107-117.

［8］ Ren Weixin，Chen Huabing.Finite element model updating in structural dynamics by using the response surface method［J］.Engineering Structures，2010，32：2455-2465.

［9］ Vestroni F，Capecchi D.Damage detection in beam structures based on frequency measurements［J］.Journal of Engineering Mechanics，ASCE，2000，126（7）：761-768.

［10］Kim G，Park Y.An automated parameter selection procedure for finite-element model updating and its applications［J］.Journal of Sound and Vibration，2008，309：778-793.

［11］Kim G，Park Y.An improved updating parameter selection method and finite element model update using multi-objective optimization technique［J］.Mechanical Systems and Signal Processing，2004，18：59-78.

［12］Fang S，Perera R，Roeck G.Damage identification of a reinforced concrete frame by finite element model updating using damage parameterization［J］.Journal of Sound and Vibration，2008，313：544-559.

［13］Mallat S.A wavelet tour of signal processing：the sparse way［M］.3rd ed.Singapore：Elsevier，2009：328-338.

［14］Sweldens W.The lifting scheme：a construction of second generation wavelets［J］.Siam Journal on Mathematical Analysis，1997，29：511-546.

［15］Sweldens W.The lifting scheme：a custom-design construction of bi-orthogonal wavelets［J］.Applied and Computational Harmonic Analysis，1996（3）：186-200.

TU311.3；U441；TH113

10.16450/j.cnki.issn.1004-6801.2015.04.010

張欣，男，1971年1月生，博士、教授。主要研究方向為結(jié)構(gòu)動力學。曾發(fā)表《Frequency modulated empirical mode decomposition method for the identification of instantaneous modal parameters of aeroelastic systems》（《Journal of Wind Engineering and Industrial Aerodynamics》2012，Vol.101）等論文。

E-mail：zhangxin@zzu.edu.cn

*河南省科技攻關(guān)資助項目（112102310453）

2013-05-08；

2013-12-10