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

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

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

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

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

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

引 言

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

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

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

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

1.1 多分辨率分析理論

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

在該內積空間內存在嵌套式閉子空間逼近序列

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

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

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

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

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

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

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

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

1.2 提升格式與第2代小波函數和尺度函數的構造

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

構造第2代內插型尺度函數與小波函數一般可以選用懶惰小波為提升起點，經過對偶提升（預測）實現合成尺度函數與分解小波的提升，再通過提升過程實現分解尺度函數與合成小波函數的提升。

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

按下式進行提升

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

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

2 低分辨率模型更新

2.1 基本方法

設更新目標為結構的剛度曲線K（x）

其中：K（x）為結構在x處的理論剛度；κ（x）為更新剛度。

結構的質量分布不更新。根據多分辨率分析理論，在j尺度上，模型的更新剛度函數為

在j+1尺度上為

當前模型結構的模態參量為（λi，φi），i=1，2，…為剛度曲線的函數。其中：λi為第i階頻率；φi為為第i階振型；實驗模態參數為（ˉλi，ˉφi），i=1，2，…。

定義誤差向量

其中：權重矩陣為

在優化求解過程中可以強調某一階模態信息的重要性，也可以降低該信息的重要性。這可以通過調整權重矩陣對角元素的大小來實現。第i階振型的相關因數為

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

2.2 兩階段更新模式

r（α）的模可能不為零，因為即使當r（α）的維數，即方程個數大于或等于修正向量α的維數時，非線性方程r（α）=0也可能無解。因此模型更新及損傷識別過程為參數優化過程，而非定解方程求解問題。

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

圖1 算例結構簡圖（單位：mm）Fig.1 Sketch of the structure under study（unit：mm）

筆者建議根據結構自身的特性，采取適當階數的模態信息，分以下兩個階段進行模型更新和損傷識別。

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

圖2 結構兩端橫斷面（單位：mm）Fig.2 Cross section of two ends（unit：mm）

2）模型的損傷更新階段。如圖1所示，在模型中引入結構損傷，此時模型的模態參數發生變化。在前一階段剛度曲線的基礎上，重復低分辨率更新過程，擬合有損傷模型的模態參數，獲得表征結構損傷特性的尺度函數與小波函數系數以及相應的損傷剛度更新曲線。

2.3 近似誤差

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

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

3 數值算例

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

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

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

表1 模型豎向振型模態參數Tab.1 Modal parameters of the vertical modes

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

為考慮模型剛度的多分辨率表示，設梁單元內力虛功為

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

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

3.2 模型校準及單裂縫工況計算結果

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

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

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

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

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

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

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

3.3 多裂縫工況計算結果

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

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

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

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

4 結束語

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

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

［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］ 侯吉林，歐進萍.基于局部脈沖響應的約束子結構修正法［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］ 郭力，李兆霞，陳鴻天.基于子結構分析的多重子步模型修正方法［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］ 張保強，陳國平，郭勤濤.使用有效模態質量和遺傳算法的有限元模型修正［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月生，博士、教授。主要研究方向為結構動力學。曾發表《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

*河南省科技攻關資助項目（112102310453）

2013-05-08；

2013-12-10