基于有限元模型修正的土木結構損傷識別方法

翁 順，朱宏平

(華中科技大學土木工程與力學學院，湖北，武漢 430074)

土木工程有限元模型廣泛應用于抗震抗風模擬、振動控制、損傷識別及優化設計等方面。有限元模型修正方法的基本思想是以結構靜動力試驗數據為基礎[1]，將試驗得到的信息與原始有限元模型分析結果進行綜合比較。通過優化約束，不斷修正模型物理參數，使理論值與試驗值基本吻合，從而獲取更精確反映實際結構特性的有限元模型。通過對比損傷前后有限元模型參數的變化，實現結構損傷識別。基于有限元模型修正的損傷識別方法，計算過程直觀、物理意義明確，能同步識別損傷位置和損傷程度，被廣泛應用于實際工程。多項標準規程中也明確規定結構安全二級評估必須實施有限元模型修正[2]。近三十年，有限元模型修正方法迅速發展[1 ? 19]，在土木工程損傷識別與健康評估中發揮著越來越重要的作用。

有限元模型修正是一個約束優化問題，屬于反問題的一種。土木工程損傷識別通常需要建立精細化的有限元模型，包含上萬乃至百萬自由度，且待修正參數多。有限元模型修正過程，需要重復計算整體結構模型的特征解，即使少數甚至一個修正參數發生改變(局部損傷)，也需要重新分析整體結構模型。大量修正參數使優化過程容易出現病態或者收斂緩慢的問題，效率極低。特別是在考慮參數不確定性、結構非線性等問題時，土木工程的有限元模型修正和損傷識別更是難以完成。例如，Xia 等[3]對西澳大利亞州的Balla Balla 橋進行了模型修正，該小型橋梁有限元模型包含4200 個自由度、1200 單元和1039 個修正參數，優化過程迭代155 次收斂，花費約420 h。對于大型土木工程，例如青馬橋[4]、廣州塔[5]包含上百萬自由度、上萬修正參數，即使使用功能強大的計算機，使用常規方法修正有限元模型也非常困難。

用于有限元模型修正的動態測量數據反映整體結構信息，但是結構損傷通常只是發生在結構的局部區域。子結構方法根據分界面上的位移協調條件和力平衡條件將整體結構分解為若干個獨立子結構，將對整體結構的模型修正轉換為對獨立子結構的模型修正。特別是在應用于損傷識別時，只需要重復分析某些特定(損傷區域)的子結構，避免重復計算整體結構特性。此外，子結構包含的待修正參數遠少于整體結構，有助于加速大規模優化問題的收斂。從而有效地減輕大型結構動力分析和有限元模型修正的計算負擔，在大型土木工程結構損傷識別中具有顯著的優勢。

本文將首先介紹基于有限元模型修正的損傷識別原理與過程。然后，詳細闡述用于土木工程損傷識別的子結構有限元模型修正方法，包括正向子結構方法和逆向子結構方法；總結基于子結構有限元模型修正的損傷識別方法在頻域、時域、不確定性分析、非線性分析中的研究進展。最后，將傳統有限元模型修正方法和子結構有限元模型方法用于一個實際超高層建筑數值模型的損傷識別中，對比分析損傷識別的精度和效率。

1 基于有限元模型修正的結構損傷識別方法

對于土木工程這一類體積龐大、參數多，且存在各類不確定性和非線性的結構，按照結構設計詳細建立的有限元模型不可避免存在各種偏差，如邊界誤差、材料參數誤差、離散化誤差等。必須采用有限元模型修正技術使得有限元模型能真實反映實際工程結構動靜力特性[5]。

現有的有限元模型修正方法可分為矩陣型修正方法和參數型修正方法。在此基礎上，針對土木工程中參數不確定性問題和非線性問題，建立了有限元模型修正的不確定性分析方法和非線性模型修正方法。

1.1 矩陣型有限元模型修正方法

矩陣型有限元模型修正方法直接修正有限元模型的剛度、質量、阻尼等系統矩陣，使得修正后的系統矩陣計算的響應與實測響應吻合。為保證系統矩陣原有的帶狀、稀疏特性，同時對矩陣優化問題施加一些約束，例如振型的正交性[6]，矩陣的對稱性或正定性[7]，使得修正后系統矩陣與實測的系統矩陣盡可能接近。為保持目標方程的一些特性，對修正矩陣的參數也施加約束，如最小范數[7]或最小秩[8]。目標方程通常是由質量、剛度、阻尼矩陣的攝動矩陣和施加約束的拉格朗日乘子的疊加[1]。楊朋超等[9]總結了各種矩陣型有限元模型修正方法，考慮對稱性、模態正交性及模態參與因子等約束，采用拉格朗日乘子算法，推導了質量矩陣和剛度矩陣的最優解。

矩陣型有限元模型修正方法具有修正結果準確和計算量小的優勢，修正后的模型能夠精確“復制”試驗測試數據[10]。其缺點是結點的連續性、修正后矩陣的稀疏性、對稱性不能得到保證；修正后的結構系統矩陣難以解釋、物理意義不明確；修正后的模型只在數學結果上與實際結構相近而不具備實際工程意義。因此，近三十年，針對矩陣型有限元模型修正方法的研究逐步減少，在結構損傷識別應用中的報道較少。

1.2 參數型有限元模型修正方法

參數型有限元模型修正方法選取結構自身物理參數(密度、單元剛度、材料參數等)作為修正參數，以有限元模型的動力響應特征值與實際結構試驗值的殘差為目標函數，重復迭代調整物理參數取值，使目標函數最小化，有限元模型修正的流程圖如圖1 所示。由于修正參數通常為結構物理參數，該類方法的修正結果具有明確的物理意義，修正后的參數可用于損傷識別、狀態評估。參數型有限元模型修正方法能夠克服矩陣型模型修正方法的局限性，因此得到了更加廣泛的應用。選取修正參數、構建目標函數、迭代優化算法和計算靈敏度矩陣是參數型有限元模型修正方法的關鍵步驟。

圖1 參數型有限元模型修正流程圖Fig.1 Flowchart of parameter-based finite element model updating

1.2.1 修正參數

修正參數的選擇是有限元模型修正的首要關鍵步驟[11]。選擇過多的修正參數不僅導致結構動力分析和靈敏度分析中的計算負荷增加，而且可能導致模型修正過程中的不適定或溢出問題[12]。因此，通常選擇不確定性較大、且對結構動靜力特性影響較大的部分參數進行修正。參數選擇方法主要分為三種：1) 綜合考慮物理意義和工程經驗的情況下選擇修正參數，通常選擇對結構特性影響大的參數，如材料特性(彈性模量、密度等)，幾何特性(截面面積、截面慣性矩)和邊界條件等[5]；2) 基于靈敏度分析方法，求解結構響應相對于參數的靈敏度，選取靈敏度大的參數作為待修正參數[12 ? 13]；3) 參數子集選擇方法，從各種組合的參數子集中挑選讓目標方程殘差最小的參數子集[14 ? 15]。

1.2.2 目標方程

目標函數量化了實際結構與有限元模型之間動態特性的差異，通常定義為實際結構實驗數據與有限元模型計算特征之間的殘差，例如：

式中：XE為實際結構試驗測量得到的特征；XA(r)為有限元模型的計算特征，可表達為修正參數r 的函數；r 為包含所有修正參數的向量。

根據測試的數據類型分類，用于構建目標函數的特征包括靜力數據(位移、應變、應力等)和動力數據(頻率、振型、頻響函數等)[16]。靜力數據精度高，但對大型土木工程施加較大的靜力荷載難以實施，且測試會妨礙工程正常運營。因此，基于靜力數據的模型修正方法在實時在線結構損傷識別中應用較少[10 ? 11]。大量研究集中在基于動態數據的有限元模型修正方法上[1, 10 ? 11]。

用于構建目標函數的動態數據包括頻域特征和時域特征。在頻域內，頻率和振型是最早和最廣泛應用于有限元模型修正和損傷識別的特征指標。實測頻率可以直接運用到有限元模型修正中。限于傳感器數目，實測振型的自由度要少于計算振型，需要通過模態縮聚或模態擴展技術與振型自由度相匹配[17 ? 18]。頻率和振型可以分別用于構建目標方程[19 ? 20]，也可以通過權重系數[20 ? 21]聯合構建目標函數。由于頻率和振型往往對局部小損傷不敏感，且受噪聲影響大，導致損傷識別的精度較低[19 ? 20]。Meruane[22]用從傳遞矩陣中識別反共振頻率代替振型來構建目標方程，反共振頻率的識別難度小且精度高，同時包括了頻率和振型的結構信息。一些模態參數的衍生指標，如模態應變能[23]、模態柔度[24]、應變振型[21]、小波時頻數據[25]等，比頻率、振型對結構局部變化更為敏感，也常用來構建目標方程。除模態參數外，頻響函數也用于構建目標函數[26]。Lin 和Zhu[27]通過推導頻響函數對結構參數的靈敏度，提出自動選擇有限頻率點方法，讓目標函數包含盡可能多的信息，并將基于頻響函數的有限元模型修正方法應用到有阻尼結構。Gang 等[28]通過恰當選擇測量頻率點和擴充自由度，改進有限元模型修正的精度和收斂性。

時域指標直接運用實測響應數據及其衍生量作為目標函數。最直接的時域指標為直接測量的結構加速度、速度、位移等指標[29]，以實測時域數據與模型時域數據的殘差構建目標方程。與頻域指標相比，時域指標避免了數據時頻轉換的誤差，適用于非線性結構等沒有表現模態特征的結構；并且時域數據點較多，有利于保證有限元模型修正中優化方程式的正定性[16]。

一些研究將動力數據和靜力數據聯合起來構建目標方程，增加了數據樣本點，同時也使有限元修正模型更能真實全面反映實際結構特性。張啟偉和范立礎[30]采用基于靜動力數據的有限元模型修正對一個懸臂梁進行了損傷識別。

1.2.3 優化算法

有限元模型修正是一個參數優化過程，即尋找一組最優參數使目標方程最小化。常用的優化算法，包括靈敏度算法[5]、神經網絡方法[31]、仿生優化算法[32 ? 34]等。其中，基于靈敏度分析的優化算法[5]，利用結構特征解對結構參數變化的敏感程度，根據解的搜索方向進行迭代優化，具有搜索速度快、超線性收斂等優點，被廣泛應用于實際工程的損傷識別和健康評估。其缺點是計算靈敏度矩陣耗時，且在非平滑連續處結構的靈敏度不存在。基于神經網絡的優化算法是一種類似于人類神經系統的信息處理技術，它通過學習有限元模型數據推出輸入與輸出的關系。雖然神經網絡具有超強的學習能力和非線性映射能力，但是其學習的準確度需要大量的學習樣本來保證[31]。仿生優化算法包括遺傳算法[32]、粒子群算法[33]、蟻群算法[34]等。這種優化算法依據生物智能選擇能耗最小、資源安排最高效得到啟發，仿生優化算法的發展依托處理設備計算速度與功能的提高。

1.2.4 靈敏度矩陣

靈敏度即是目標方程關于修正參數的偏導數，其一方面可以用于選取最優的修正參數，另一方面可給優化目標方程提供快速的搜索方向，加快優化過程收斂。

以頻域指標構建目標函數時，計算特征值(頻率)和特征向量(振型)靈敏度主要包括三類方法：1) 有限差分法，估計特征解靈敏度為設計點特征解與其附近有限步長點的特征解關于步長的斜率[35]。該方法的誤差與步長的選取有關；2) 模態法[36]，將特征向量靈敏度表達為所有特征向量的線性組合，該方法需要計算所有特征向量，大型結構通常采用模態截斷的方法取部分模態，會影響計算精度；3) Nelson 的方法[37]，是將特征向量靈敏度表達為一個齊次項和非齊次項之和，精度和效率較高。在這三種方法的基礎上，還發展出其他的計算方法，包括幾何法[38]、Lanczos 法[39]、迭代法[40]、奇異分解法[41]、攝動法[42]、子結構法[43]、模型縮聚法[44]等。頻響函數靈敏度的計算可通過模態疊加法，將頻響函數靈敏度表示為所有模態的疊加，進而通過模態截斷的方法提高計算效率，通過補充舍棄模態的貢獻以提高計算精度[45]。時域內動力響應靈敏度的計算，可對運動方程求偏導，通過Newmark 方法或者Newton Rapson 方法計算偏導方程的響應，作為時域內動力響應靈敏度[46]。

參數型有限元模型修正方法優化結構物理參數，物理意義明確，對實際工程損傷評估具有指導意義；靈敏度分析為參數優化提供快速搜索方向，是一種被廣泛應用于實際工程的有限元模型修正方法。該類方法的進步在一定程度上得益于優化算法、靈敏度計算等數學方法的發展，以及硬件設備計算能力的提高，未來將在大型土木工程損傷評估中發揮越來越重要的作用。

1.3 不確定性有限元模型修正方法

有限元模型和測試數據都存在不確定性。模型不確定包括材料制作誤差、材料本構關系建模誤差、邊界條件建模誤差、單元劃分誤差、非結構特性建模誤差等[47]。測量不確定性主要包括測量噪聲和數據處理誤差。常用三種方法來考慮模型修正過程中的不確定性：蒙特卡洛法、攝動法、貝葉斯法。

蒙特卡洛法生成大量服從特定概率密度分布的樣本點，對每個樣本點參數分別進行有限元模型修正，然后從修正結果中估計修正參數的統計特性[48]。Mares 等[49]詳細闡述了蒙特卡洛法及其在概率參數估計上的應用。蒙特卡洛法需要生成大量的分析數據來保證其修正結果的精度，計算效率低。

攝動法通過引入微小變化量考慮不確定性[50]，給有限元模型修正方程的每個變量一個微小變化量，推導一個近似的微分方程來求解函數的不確定性，參數的一階矩和二階矩常用來估計參數的均值和方差。與傳統的基于靈敏度的有限元模型修正方法一樣，攝動法要求修正參數的初值要接近真實值，當攝動量較大時會產生誤差，需要使用高階攝動方程[51]。

基于貝葉斯理論的有限元模型修正方法，根據主觀已知修正參數先驗分布信息和測試數據信息，識別修正參數的后驗概率密度分布。Beck 和Karafygiotis[52 ? 53]詳細闡述了基于貝葉斯理論的模型修正過程，該方法不僅能獲得修正參數的最優解，還能從它們聯合概率分布中評估參數不確定性水平。Lam 等[54]從環境數據中評估參數不確定性水平，并基于貝葉斯模型修正方法對IASC-ASCE基準模型進行了損傷識別。Yuen 和Katafygiotis[55]將該方法擴展到隨機模型輸入未知的情況。該方法不需要響應是平穩的，也不需要基于輸入譜密度的參數先驗信息。Cheung 和Beck[56]提出利用混合蒙特卡洛來求解貝葉斯模型修正問題，適合于求解高緯度參數空間的優化問題，且不局限于模型的類別(線性或非線性)和測試數據的類型。

有限元模型和測試數據的不確定性會影響有限元模型修正和損傷識別的結果，研究參數識別過程中的不確定性問題勢在必行。隨著統計理論與大數據分析方法的發展，為不確定性分析提供新的契機，這也將為考慮不確定性分析的有限元模型修正和損傷識別創造良好的條件。

1.4 非線性有限元模型修正方法

嚴格來講土木工程結構是非線性的，線性只是特例或一種簡化方法。例如，混凝土等材料本身的非線性引起的結構非線性、結構發生大變形時產生的幾何非線性、結構阻尼耗散的非線性、結構邊界條件及狀態的非線性等。結構的各種損傷也呈現出典型的非線性特征，如裂縫的產生和發展、節點的松動和滑移、剛度退化等。土木工程結構遭受地震、強風等較強的荷載激勵時，表現出較強的非線性特征，必須采用非線性模型。目前，線性系統的模態分析理論并不適用于非線性系統，非線性有限元模型修正方法的發展沒有線性模型修正方法成熟，非線性有限元模型修正技術尚待進一步發展[16]。

非線性模型修正方法主要包括頻域內的諧波平衡方法[57]、本構方程誤差法[58 ? 59]和時域內的恢復力面法[60]、本征正交分解[61]。頻域方法一般將非線性特征表達為線性貢獻和非線性貢獻的疊加[59]，或將非線性貢獻的多維度高階特征線性化[58]。由于線性化過程損失了一些結構信息，該方法不適用于復雜的非線性結構的有限元模型修正[62]。時域方法通過輸入輸出振動方程識別出非線性方程的恢復力[60]，或者將時間序列數據分解為非線性系統的正交基[61]，進而以恢復力或非線性正交基構建目標函數來修正非線性參數。Silva等[63]比較了頻域的非線性模型修正方法(諧波平衡法、本構方程誤差法)和時域的非線性模型修正方法(恢復力面法、本征正交分解)，認為時域方法比頻域方法計算效率高，其中恢復力面法是耗時最小且識別結果精度最高的模型修正方法。

土木工程結構實質是非線性的，非線性有限元模型修正和損傷識別更能反映土木工程的本質特征。隨著非線性動力分析的發展和設備計算能力的提高，非線性有限元模型修正將為土木工程損傷評估提供更貼切、更可靠的信息。

1.5 基于有限元模型修正的結構損傷識別方法

當結構損傷導致其剛度參數發生變化，運用有限元模型修正方法來識別損傷前后結構的剛度參數變化，從而根據剛度參數變化的位置和程度實現損傷的定位和定量。基于有限元模型修正的損傷識別過程分為兩個階段，具體流程如圖2所示。

圖2 基于有限元模型修正的損傷識別過程Fig.2 Flowchart of damage identification based on finite element model updating

第一個過程，建立結構未損狀態下的有限元模型，即基準模型。初始有限元模型尚不能準確反映真實結構的特性，需要通過對其材料參數、截面尺寸、約束條件進行修正，避免初始建模誤差對第二個過程中損傷識別結果的影響。當第一個過程的有限元模型修正完成，將修正好的模型作為精準的未損狀態結構有限元模型。在第二個過程中，由于第一個過程已經對結構質量約束等進行修正，則認為只有損傷造成剛度變化，其他因素是不變的。因此，以單元剛度作為待修正參數，利用現場實測的損傷結構的試驗數據對基準模型進行模型修正，得到損傷后的剛度參數，從而識別損傷位置和損傷程度。

Brownjohn 等[64]利用實測振型和頻率，采用基于靈敏度的有限元模型修正方法識別了實際橋梁的損傷狀況。Fritzen 等[65]利用基于頻響函數靈敏度的有限元修正方法識別實際結構損傷位置和損傷程度。Abdel[66]、Jaishi 和Ren[23]分別推導了模態曲率靈敏度和模態柔度靈敏度的計算方法，建立了基于模態曲率和模態柔度的有限元模型修正方法，完成實際橋梁的損傷識別。Lu 和Law[67]提出了基于動力響應靈敏度的模型修正方法，識別了一個小型橋梁的損傷狀況。針對實際工程的損傷識別中，為提高有限元模型修正計算效率，通常會采用簡化的有限元模型，或者從有限元模型的眾多不確定性參數中選取少數待修正參數，在不同程度上通過犧牲有限元模型的精度來提高模型修正過程的效率，修正后的模型難以準確地反映結構的一些真實損傷狀況。

基于有限元模型修正的損傷識別方法，計算過程直觀、物理意義明確，能識別結構損傷發生的位置和損傷程度，并根據模型物理參數的變化評估結構安全狀況，是了解和掌握結構安全狀況最直接、最準確的途徑。然而，土木工程結構有限元模型龐大，不確定性參數多，導致基于有限元模型修正的損傷識別過程效率低、精度差，特別是在考慮不確定性和非線性后，耗時巨大甚至是無法實施。雖然現有硬件設備具有日益強大的計算能力，為大型土木工程有限元模型修正計算創造了良好的條件，但是大型土木工程精細化有限元模型修正和損傷識別仍然難以實現[3 ? 5]，計算方法上的突破仍然是土木工程研究者與實踐者需要努力的方向。

2 基于子結構有限元模型修正的結構損傷識別方法

為保證損傷識別的精度，土木工程通常需要建立比較精細的有限元模型，包含大量單元、節點、自由度和待修正參數。有限元模型修正屬于結構動力學中的逆問題，大尺寸模型的動力分析精度極低，大量待識別參數使參數識別過程容易出現病態或者收斂緩慢的問題。

用于有限元模型修正的動態測量數據反映整體結構信息，但是結構損傷通常只是發生在結構的局部區域。子結構方法將整體結構分解為獨立的子結構，通過協調條件將對整體結構的分析轉換為對獨立子結構的分析，是解決有限元模型修正精度低、效率低的有效方法。子結構方法應用于大型土木工程結構有限元模型修正有如下優勢[68 ? 70]：

1) 各子結構之間相互獨立，采用不同的方法獨立分析，獨立存儲，并行計算；

2) 將對整體結構的分析轉換為對獨立子結構的分析，子結構模型尺寸遠小于整體結構模型，將減少模型分析時間，從而降低對分析計算設備的需求；

3) 獨立子結構修正參數的數量遠小于整體結構修正參數的數量，從而加速優化過程的收斂，獲取更精確的參數優化結果。

4) 子結構有限元模型修正方法用于損傷識別中，由于損傷往往發生在結構的局部區域，只需對少數關鍵子結構進行分析，避免了對整體結構進行反復運算。

子結構有限元模型修正方法分為兩類：第一類為正向子結構方法，通過子結構特征參數組集得到整體結構特征參數，進而修正結構參數使得整體模型特征參數與測量得到的整體特征參數相吻合；第二類為逆向子結構方法，將整體結構試驗特征參數分解為獨立子結構試驗特征參數，進而修正獨立子結構模型。

2.1 正向子結構有限元模型修正方法

正向子結構方法的核心思想是將整體結構特征參數表達為獨立子結構特征參數的疊加。該方法的過程是[70]：首先，將整體結構有限元模型劃分為若干個獨立的子結構模型，計算子結構模型的振動特性；然后，組集子結構振動特性得到整體結構的振動特性；最后，通過優化過程調整子結構物理參數，使得整體結構振動特性與實測結構試驗數據相吻合，進而實現損傷識別(如圖3 所示)。用作有限元模型修正的振動特性包括頻域特性(頻率、振型、頻響函數等)和時域特性(加速度、位移等時程響應)。

圖3 正向子結構有限元模型修正流程圖Fig.3 Flowchart of model updating based on forward substructure

1) 基于子結構的頻域特征計算方法

頻域子結構方法主要用來快速計算結構的特征解(頻率)和特征向量(振型)。Kron 子結構[68]方法通過拉格朗日乘子和虛功原理在相鄰子結構界面施加位移約束組集各個子結構特征參數。該方法精度高，且適用于復雜的分界面，但是Kron 子結構方法需要計算每個子結構所有階模態，效率極低，對較大子結構難以完成。Weng 等[69]只保留了部分模態來改進Kron 子結構方法的計算效率，同時通過剩余柔度補充舍棄模態的貢獻來保證方法的精度，該方法也被擴展到快速計算特征解靈敏度[71]。模態綜合法是將結構特征解表達為子結構模態的疊加[72 ? 75]。子結構模態分為正交模態、剛體模態、約束模態和附加模態。根據子結構界面約束類型，模態綜合法分為固定界面模態綜合法和自由界面模態綜合法。子結構方法也被用來通過子結構特征參數快速計算整體結構頻響函數[75 ? 80]。Klerk 等[75]綜述了基于子結構的頻響函數計算方法，包括阻抗耦合法[76]、導納耦合法[77]、拉格朗日乘子耦合法[78]。Lim 和Li[79]用最小二乘法和縮聚的奇異值分解法來改進子結構方法計算頻響函數的精度。Law 和Ihlenfeldt[80]將基于子結構的頻響函數計算方法拓展到計算結構多點柔度特征。

2) 基于子結構的時域特征計算方法

時域子結構方法快速計算結構時程響應[81,82]和脈沖響應函數[83]。Zhu 等[82]基于Kron 子結構方法將整體結構時程響應表達為獨立子結構模態參數的疊加，為提高計算效率只保留了部分子結構低階模態，舍棄的高階模態的貢獻通過考慮一階效應和二階效應的指標來補償，并推導了基于子結構的動態響應靈敏度計算方法，應用于超高層建筑的損傷識別。Gruber 等[84 ? 85]將模態綜合法拓展到時域計算結構位移，并轉換到狀態域計算非比例阻尼線性系統時程響應。時域子結構方法通過一定約束條件組集子結構的脈沖響應函數，也被用于快速計算整體結構脈沖響應函數。Gordis[83,86]通過對子結構脈沖響應函數進行杜哈梅積分和模態疊加法得到整體結構脈沖響應函數，用于瞬時沖擊或敲擊荷載作用下的動力分析和模型修正。基于子結構的脈沖響應函數計算方法也被應用到剛性-彈性混合節點的結構的有限元模型修正中[87]。

在正向子結構有限元模型修正過程中，當結構損傷引起局部參數發生改變，只需要重新計算一個或者幾個子結構的模態參數，其他子結構保持不變，即可快速計算整體結構振動特性，通過修正少數子結構模型實現損傷識別，極大地提高了有限元模型修正的效率。

2.2 逆向子結構有限元模型修正方法

逆向子結構方法通過探索整體結構與獨立子結構的位移、力等參數的相似性，以及子結構分界面處位移平衡、力協調等約束條件，建立整體結構振動特性與獨立子結構振動特性之間的聯系。進而，將整體結構的試驗特征參數分解為獨立子結構的試驗特征參數。將子結構完全從整體結構中分離出來，成為獨立自由的個體。最后，建立獨立子結構有限元模型，以獨立子結構模型的特征參數和子結構試驗特征參數的殘差作為優化目標，以獨立自由的子結構模型為修正對象，其他的子結構并不參與模型修正過程，基于逆向子結構的有限元模型修正過程如圖4 所示。逆向子結構方法的實質是將整體結構上測得的特征參數通過建立數學方程分解為局部區域的特征參數，進而對某一個子結構進行有限元模型修正。從整體結構試驗數據中提取的子結構特征參數包括頻域特征參數和時域特性參數。

圖4 逆向子結構有限元模型修正流程圖Fig.4 Flowchart of model updating based on inverse substructure

1) 頻域內逆向子結構法

頻域逆向子結構方法的關鍵是建立子結構特征解和整體結構特征解的關系，然后再從整體結構的模態試驗數據中提取出子結構頻域特征參數。Doebling 等[88]根據預先假定的子結構連接性和應變能量分布從整體剛度和柔度矩陣中提取出了子結構的剛度矩陣，進而從形狀方程和幾何關系中計算子結構特征向量，并通過線性最小二乘法求解子結構特征值。Alvin 和Park[89]依據力法基本原理，基于力轉換矩陣從整體柔度矩陣中提取出子結構柔度矩陣。由于轉換矩陣的計算非常耗時，這種方法主要用于簡單梁結構的計算。Park和Reich[90]綜述了從整體測量數據中提取子結構特征的柔度法，包括自由-自由子結構柔度法、基于變形的柔度法、基于應變的柔度法。Hou 等[91]根據虛擬變形方法推導了子結構隔離方法，將其他子結構對目標子結構的作用等效為一個虛擬力，利用分界面的虛擬力變形隔離出一個不受其余子結構影響的目標子結構。

Weng 等[92]根據位移協調條件和力平衡方程從整體模態數據提取出了子結構柔度矩陣，并建立正交投影算子剔除子結構剛體模態用于自由子結構的分析[93]。從子結構柔度矩陣中提取子結構頻率和振型，以獨立子結構為對象構建目標函數，計算靈敏度矩陣，并修正獨立子結構模型[94]。由于結構損傷通常發生在局部區域，子結構模態參數比整體結構模態參數對損傷更加敏感。該方法用于廣州塔的有限元模型修正和損傷識別[94]，局部小損傷對子結構模態參數的改變量為1.84%，而對整體結構模態參數的改變量僅為0.01%，證實了子結構模態參數比整體結構模態參數對局部損傷更加敏感。以所提取的子結構試驗模態參數修正目標子結構有限元模型，系統方程的尺寸從21 690×21 690 縮小為2736×2736，基于逆向子結構的有限元模型修正時間為整體有限元模型修正時間的10%，極大地提高了基于有限元模型修正的損傷識別的精度和效率。

2) 時域內逆向子結構法

在時域內，如果能提取目標子結構與其他子結構的界面力，那么就可以將界面力當作輸入力作用在目標子結構上，將目標子結構從整體結構中分離出來，作為一個獨立的結構進行動力分析[95]，并對該獨立子結構進行有限元模型修正和損傷識別。因此，時域內逆向子結構方法的關鍵是獲取目標子結構與其他子結構的界面力。

獲取子結構界面力最直接的方法是測量界面處的響應，通過界面處的響應來計算界面力。Koh等[95]測量每個子結構分界面結點的位移、速度、加速度，通過參數優化方法識別其他子結構對目標子結構的界面力。Yun 和Lee[96]構建界面處加速度響應與界面力的時間序列模型，通過神經網絡算法從復雜結構中提取子結構特征參數和子結構界面力[97]。這些方法需要測量界面處所有自由度的響應數據，在實驗室模型試驗中取得較好的結果。在實際工程中，測量目標子結構所有界面結點的所有響應難以實施。一些研究試圖通過數學變換來消除目標子結構振動方程中與界面結點有關的項，以避免測量子結構界面處所有結點響應。Koh 和Shankar[98]使用傳遞方程來關聯一個結點的響應和另一個結點的激勵，測量同樣激勵條件下目標子結構上多組響應數據，建立轉換矩陣將子結構界面上未測量的響應轉換為已測響應的函數。Tee 等[99]將模型縮聚和恢復方法融合，解決子結構界面處不完備測量的問題。

針對實際工程中目標子結構界面處測量不完備的難題，另一種處理方法是將子結構界面力當作未知參數，在有限元模型修正過程中與其余待識別參數同步識別。Li 和Law[100]融合響應重構方法和傳遞方程，采用基于靈敏度的有限元模型修正方法同步識別子結構界面力和損傷參數。Zhu等[101]將移動荷載和子結構界面力表示為切比雪夫多項式，推導了在狀態域內結構動力響應關于結構單元剛度參數和切比雪夫多項式因子的靈敏度，通過子結構響應重構和靈敏度分析方法同步識別移動列車荷載作用下軌道參數和列車荷載。

逆向子結構方法建立整體結構動態特征與獨立子結構動態特征的聯系，將子結構完全從整體結構中分離出來，成為獨立自由的個體，對獨立自由的子結構模型進行有限元模型修正。用于大型土木工程的損傷識別時，只需要建立局部損傷區域的有限元模型，修正獨立的局部子結構模型，而其他大部分未損區域不參與有限元建模和修正過程，是提高大型結構有限元模型修正計算效率最理想的途徑。由于整體結構動態參數包含測量噪聲等誤差，且整體結構特征參數與子結構特征參數關系復雜，在提取子結構特征參數時這些誤差會被放大和積累，因此，如何在數學方法上降低整體-子結構關系的復雜性，提高子結構試驗特征參數的精度，是目前逆向子結構方法亟待解決的難題。

2.3 考慮不確定性的子結構有限元模型修正方法

大型結構有限元模型修正過程非常耗時，考慮不確定性分析后，需要對大量統計數據進行重復分析，勢必進一步增加計算負擔。一方面，子結構方法可以對獨立子結構分別建立統計模型，對包含支座、損傷等不確定性較大的子結構進行重復統計分析，避免了反復分析整個結構。另一方面，對每個子結構分別建立統計模型，而不是對整體結構建立統一的統計模型，精度更高。

Beck 等[102]將貝葉斯方法引入子結構方法中，通過子結構參數的先驗概率估計整體結構參數的不確定性，減少了不確定性有限元模型修正的時間。Tran 等[103]將子結構模態綜合法與蒙特卡洛方法相結合，Chentouf 等[104]評估了兩種模態綜合法在基于蒙特卡洛不確定分析中的魯棒性。蒙特卡洛方法產生大量樣本點，并對大量樣本點統計數據進行重復有限元模型修正和損傷識別計算，子結構方法有效減少了大量重復計算負擔。梁鋒[105]將子結構方法與攝動法相結合，建立子結構物理參數、子結構模態參數、整體結構模態參數、整體結構統計參數的定量關系，基于攝動法推導了不確定性參數的一階和二階統計矩，提高了大型結構不確定性分析的效率。

子結構方法可以對各獨立子結構分別建立統計模型，不僅精度更高，而且將重復的統計計算限制在少數子結構內，為大型結構損傷識別及其概率估計提供了高效率、高精度的途徑。

2.4 非線性子結構有限元模型修正方法

土木工程非線性通常存在于結構的局部位置，例如非線性支撐或非線性連接，結構損傷也是一種局部非線性發展過程。局部非線性特征使整個結構呈現非線性行為，必須對整體結構進行非線性動力分析。常用直接積分法、諧波平衡法求解非線性方程，大型結構的非線性計算非常耗時。子結構方法考慮局部非線性特征分別建立非線性子結構和線性子結構，對線性子結構采用線性疊加、動態縮聚等高效的動力分析方法，將耗時的非線性動力分析限制在局部子結構內，極大地提高非線性動力分析及有限元模型修正的效率。

Praveen 和Padmanabhan[106]依據整體結構的非線性特征，將剛度矩陣和質量矩陣分為主自由度和從自由度，與非線性區域和外力自由度相關的部分作為主自由度，其余的線性區域作為從自由度。對主自由度進行耗時的非線性動力分析，非線性修正力施加到線性部分作為非線性效應影響，以提高計算效率。Weng 等[107]將整體非線性結構分解為線性子結構和非線性子結構，將線性子結構表達為低階模態的疊加從而縮減線性子結構方程尺寸，非線性結構保持原狀。由于土木工程中大部分區域為線性，耗時的非線性動力分析和參數識別被限制在少數非線性子結構內，所計算的地震荷載作用下的動力響應被用于非線性參數識別中。Apiwattanalunggarn 等[108]引入子結構非線性正交模態，通過模態綜合法縮減線性子結構的運動方程，通過非線性正交模態縮減非線性子結構的運動方程，非線性正交模態是線性正交模態的拓展，它們通過線性約束模態耦合起來[109]。

非線性子結構方法將耗時的非線性動力分析限制在局部子結構內，并極大地縮聚線性子結構尺寸，為土木工程局部非線性動力分析、參數識別和損傷識別開辟了一個高效的途徑。

2.5 基于子結構有限元模型修正的結構損傷識別方法

由于土木工程損傷通常發生在局部區域，在判定結構損傷區域后，可只對發生損傷的子結構進行有限元模型修正，而其他大部分未損結構保持不變，將提高損傷識別的精度和效率。

Weng 等[70]推導了基于正向子結構的特征解和特征解靈敏度求解方法，建立了基于子結構有限元模型修正的損傷識別方法，完成澳大利亞Balla Balla 橋的損傷識別；并進一步提取獨立子結構試驗柔度矩陣，建立了逆向子結構有限元模型修正方法，應用于廣州塔的損傷識別[92,94]。Zhu 等[83]利用子結構方法推導結構位移和位移靈敏度的快速計算方法，通過子結構有限元模型修正過程對武漢長江航運中心進行損傷識別。Yu 等[110]推導了固定界面模態綜合法、自由界面模態綜合法[111]計算結構特征靈敏度，并基于特征解和特征靈敏度對一框架建筑進行了損傷識別。Lam 和Yang[112]綜合貝葉斯模型修正和子結構方法，依據試驗模態數據實現了鋼塔的損傷識別。Jensen 等[113]在時域內建立了貝葉斯模型修正和損傷識別方法，其中子結構模態綜合法用來縮小模型的尺寸。Xu 等[29]提出了多尺度子結構模型修正方法，首先在子結構層面識別出損傷的子結構，然后在單元層面識別出損傷發生的單元和程度。

土木工程結構龐大而損傷往往只是發生在局部區域，子結構方法將對整體結構的分析轉換為對獨立子結構的分析，將耗時的優化計算、不確定性分析和非線性計算限制在局部子結構內，避免對整體結構大尺寸方程實施各類耗時的運算，在保證損傷識別高精度的同時提高損傷識別的效率，為土木工程基于有限元模型修正的損傷識別開辟了新的途徑。

3 工程應用

本節將上述有限元模型修正方法和子結構有限元模型修正方法用于超高層建筑武漢長江航運中心模型的損傷識別。武漢長江航運中心大樓地上66 層，建筑高度335 m，如圖5 所示。主體結構為第1 層～第64 層，高度306 m，采用外框架-核心筒結構體系。外框架由4 根鋼管混凝土柱和16 根型鋼混凝土柱構成，截面尺寸50 m×50 m。核心筒為鋼筋混凝土剪力墻，截面尺寸30 m×30 m。

依據施工圖紙建立該實際工程有限元模型，包含9112 個單元、3950 個結點、23 364 個自由度。首先，采用現場模態試驗測得的頻率和振型修正有限元模型的剛度參數，完成第一階段有限元模型修正。然后，以修正后的有限元模型為基礎，進行第二階段有限元模型修正，即損傷識別。

圖5 武漢長江航運中心及其有限元模型Fig.5 Wuhan Yangtze River Navigation Center and its finite element model

在第二階段有限元模型修正過程中，分別采用傳統的有限元模型修正和子結構有限元模型修正方法作對比分析。在采用子結構方法時，將有限元模型劃分為9 個子結構。由于該實際工程目前沒有發生損傷，因此數值模擬結構損傷。對局部區域(子結構3)中單元21 和單元77(剪力墻單元)的單元剛度分別折減20%和30%來模擬損傷，對整體有限元模型施加如圖6 所示水平地震荷載作用，地震荷載持續30 s，采樣頻率是0.03 s，總的時間步為1001，采用Newmark 法計算損傷后模型的動位移作為現場試驗數據。有限元模型修正過程中，目標方程為試驗動位移和有限元模型計算位移的殘差，當目標方程小于2 × 10?15時，模型修正過程停止。有限元模型修正的待修正參數選為第三個子結構對應所有單元的剛度，共336 個剛度參數。傳統以整體結構為對象進行有限元模型修正時，一般取所有單元剛度參數為待修正參數，即9112 個單元參數，現有計算設備難以完成對9112 個單元參數的優化。因此，本算例中兩種有限元模型修正方法均采用336 個修正參數，即以整體結構為對象進行損傷識別時只取局部336 個參數進行優化。

圖6 地震荷載激勵Fig.6 Earthquake excitation

首先采用基于正向子結構有限元模型修正的損傷識別方法，每個子結構只保留了少量的低階模態來組集得到縮聚的振動方程，進而計算出整體結構的動位移，并對縮聚的振動方程求偏導計算動位移靈敏度。在振動方程中補充一個額外項來考慮舍棄的高階模態的貢獻。保留的模態數量會影響子結構方法的精度，保留的模態越多，精度越高。但是，保留模態的數目越多，計算效率會大大降低。模型修正越接近最優解，對位移和位移靈敏度的精度要求越高。因此，為了同時保證子結構模型修正的精度和效率，在不同模型修正階段保留不同數量的模態。最初的前6 階迭代中，保留30 階模態。隨著有限元模型修正過程逐漸接近收斂，保留模態的數目不斷增加。在最后的3 次迭代中，保留130 階模態。有限元模型修正過程在滿足目標方程小于2 × 10?15時完成運算。基于子結構有限元模型修正的損傷識別結果如圖7所示，識別出單元21 和77 的剛度折減值分別為0.2 和0.3，而其他單元識別的剛度折減值均為0，該方法高精度地識別出預先設定的損傷位置和損傷程度。同樣也應用傳統的基于整體結構有限元模型修正的損傷識別方法，采用相同的目標函數、靈敏度計算和優化過程。由于本算例數值模擬結構損傷，且未考慮不確定性和噪聲影響，傳統的基于整體結構有限元模型修正的損傷識別方法同樣能高精度地識別結構損傷位置和損傷程度。

圖7 基于有限元模型修正的損傷識別結果Fig.7 Damage identification by finite element model updating

表1 比較了兩種方法的迭代過程和計算時間，圖8 對比了兩種方法的收斂曲線。基于子結構有限元模型修正的損傷識別方法，當模態保留30 階時，振動方程和靈敏度方程的尺寸被縮聚為270×270，每次迭代花費時間為0.18 h。在最后3 次迭代，振動方程/靈敏度方程的尺寸為1170×1170，每次迭代花費時間為0.90 h。子結構方法總共需要23 個迭代步和9.78 h 收斂。整體結構方法計算整體模型的動力響應與靈敏度，振動方程/靈敏度方程尺寸為23 364×23 364，每個迭代步花費14.17 h，一共花費了14 個迭代步和198.38 h 達到收斂。綜上所述，在相同收斂條件、相同計算精度的前提下，基于子結構有限元模型修正的損傷識別方法花費的時間是傳統整體模型修正方法的4.93%。

表1 有限元模型修正計算時間和迭代步數的比較Table 1 Comparison of computational time and iterations by finite element model updating

子結構有限元模型修正能極大地提高損傷識別的效率，原因主要在如下三個方面：1) 子結構方法極大地縮減了模型尺寸，當子結構模態保留30 階時，振動方程和靈敏度方程的尺寸由23 364×23 364 縮減為270×270；2) 計算靈敏度是一個比較耗時的過程，在采用子結構方法計算靈敏度時，各個子結構是相互獨立的，通過計算第3 個子結構的靈敏度矩陣來組集得到整體結構的靈敏度，而其余8 個子結構的靈敏度矩陣為零，對獨立子結構計算靈敏度比對整體結構計算靈敏度效率高；3) 局部子結構的待修正參數遠小于整體結構，本算例中對子結構方法和整體結構方法均取局部336 個修正參數，整體結構有限元模型修正消耗198.38 h。在實際結構損傷識別中，以整體結構為對象時，一般取所有單元剛度參數為待修正參數，也即9112 個單元參數。可以預測，同時優化9112 個單元參數的耗時將成級數倍增加，現有的計算設備難以完成。

圖8 整體結構有限元模型修正和子結構有限元模型修正的收斂曲線Fig.8 Convergence of global-based and substructurebased finite element model updating process

4 結論

本文詳細介紹了有限元模型修正方法及其在土木工程損傷識別中的應用，并總結了不確定性和非線性有限元模型修正方法。土木工程結構龐大而損傷往往發生在局部區域，針對大型土木工程有限元模型修正效率低的難題，詳細闡述了子結構有限元模型修正方法，并介紹了適合于土木結構損傷識別的正向子結構方法和逆向子結構方法。主要結論如下：

(1)基于有限元模型修正的損傷識別方法計算過程直觀、物理意義明確，能同步識別結構損傷位置和損傷程度，是一種直接、有效的土木工程損傷識別技術。然而，土木工程有限元模型包含大量節點、單元和待修正參數，大型結構有限元模型修正效率極低，甚至無法完成。

(2)土木工程尺寸龐大而損傷通常發生在局部區域，子結構方法將對整體結構的分析轉換為對獨立子結構的分析，將耗時的優化計算、不確定性分析和非線性計算限制在局部子結構內，通過局部子結構有限元模型修正實現損傷識別，有效提高了損傷識別的精度和效率，為傳統基于有限元模型修正的損傷識別技術開辟了新的途徑。

基于有限元模型修正的結構損傷識別方法在結構健康監測領域已經取得了較好的理論成果和工程應用，多個結構健康評估規范中均規定結構二次評估需要實施基于有限元模型修正的損傷識別，該方法正在為結構安全評估提供越來越重要的信息。基于有限元模型修正的損傷識別方法與子結構方法、貝葉斯方法、人工智能與大數據方法等先進方法的深入融合，以及計算設備的高速化、智能化，將推動土木工程安全評估進一步向高精高效、智能化、信息化的可持續方向發展。