基于多目標優化的混凝土斜拉橋靜動力有限元模型修正

彭 濤, 田仲初, 張建仁, 鄭萬泔

(1.長沙理工大學 土木工程學院，長沙 410114；2.南京安正軟件工程有限責任公司，南京 210009)

近年來，隨著人民日益增長的美好生活需要，大跨度橋梁結構的健康監測、損傷識別、安全評估以及剩余壽命預測等逐漸成為常態，這些都需要一個能滿足工程精度要求、反映結構真實性態的基準有限元模型作為基礎[1-3]。由于橋梁結構實際施工中不可避免的存在誤差、材料參數的不確定性和離散性、邊界條件的失真模擬或簡化、有限元模型固有的階次誤差和離散誤差等，依據相關圖紙和規范等建立的初始有限元模型和真實結構不可避免的會存在一定偏差，從而使得有限元分析預測的結構響應與結構表現的實際響應之間的誤差超過工程允許范圍，有必要對有限元模型進行修正[4]。有限元模型修正是在保證模態參數自身準確性和修正參數物理意義明確的前提下，利用靜、動力實測結果(位移、應變、頻率、振型等)修正有限元模型的相關參數，使模型計算結果與實測值趨于一致[5]。

靜、動力協同的有限元模型修正可以增加修正過程中的可利用信息，克服單獨運用單一測試信息進行修正帶來的不足。目前，大跨度橋梁結構從施工到運營一般會經歷施工監控、交工靜動載試驗等試驗檢測過程，這為有限元模型修正和基準模型的建立提供了豐富、可靠的初始靜動力測試數據，促進了聯合靜動力的模型修正方法的發展[6]。張啟偉等[7]利用結構實測數據修正結構剛度的子結構算法，推廣為聯合運用靜動力實測數據的模型修正和損傷識別方法，并在懸臂梁上進行了應用。宗周紅等首次采用動力模態柔度和靜力位移聯合構造目標函數，提出了一種聯合靜動力的有限元模型修正方法，并將該方法成功應用于一座加固后的剛架拱橋。姚昌榮等[8]以某斜拉橋成橋靜載實測位移和模態試驗中實測頻率作為狀態變量，以各項實測值與相應計算值組合構建目標函數，對該橋的初始有限元模型的進行了修正。方志等[9]基于一個混凝土斜拉橋模型采用參數靈敏度分析和ANSYS優化技術，對靜動力試驗獲得的靜力位移和模態測試結果進行了修正，并利用不同狀態下的實測結果檢驗了修正后的模型。魏錦輝等[10]基于結構靜載位移和振動頻率等現場實測數據，構造聯合靜動力的有限元模型修正的目標函數，提出了一種基于響應面方法的橋梁靜動力有限元模型修正技術。

有限元模型修正本質上是個優化問題，采用不同的靜、動力數據的組合方式可以形成不同的目標函數，構建合理的包含靜動力響應信息的目標函數是此種修正方法的核心問題，迄今已有不少的學者基于靜動力實測數據，針對各種橋梁結構或室內模型的有限元模型修正，提出了多種卓有成效的方法。但這些方法大都采用線性加權求和法建立優化目標函數，即給每一個目標向量賦予一個權重，將各個目標分量乘上各自相應的權重系數后組合成為一個新的目標函數，從而將其轉化成一個單目標優化問題求解。這種處理方法的不足之處是權重的賦值雖然考慮了各個目標的相對重要程度，但權重系數的取值方法目前還缺少相應的研究和理論支撐，有較大的主觀性。本文提出了一種基于多目標優化的聯合靜動力有限元模型修正方法，并將其應用到混凝土斜拉橋的有限元模型修正。

1 多目標優化方法

在工農業生產和科學研究中的優化問題大多數屬于多目標優化問題。當優化問題的目標函數為一個時稱之為單目標優化問題，目標函數有兩個或大于兩個且需要被同時處理的優化問題就被稱之為多目標優化問題(Muti-objective Optimization Problems，MOPs)。典型的多目標優化問題數學模型由n個決策變量、m個目標變量和若干約束條件組成，其數學描述如下[11]

(1)

式中：x為n維決策變量，組成決策空間；F(x)為目標函數；m為目標函數個數；gj(x)為不等式約束；hk(x)為等式約束；j，k為約束的個數。

多目標優化和單目標優化有著本質的不同，多目標優化問題各目標之間相互制約，不存在唯一的全局最優解，而是存在一個最優解的集合，其特點是至少存在一個目標優于其他所有的解，稱為Pareto最優解，可描述為：當xi優于xj時，滿足

任意k∈[1,n]，fk(xi)≤fk(xj)

(2)

存在k∈[1,n]，fk(xi)

(3)

所有Pareto最優解組成的集合成為多目標優化問題的最優解集；所有Pareto最優解對應的目標函數值構成的區域稱為Pareto前沿。以工程中常見的雙目標優化問題為例，如圖1所示，AB曲線為Pareto前沿，曲線上點的縱橫坐標分別表示兩個目標函數值，x為Pareto前沿上的一個解，其彎曲角(bend angle)定義如下

θ=θL-θR

(4)

式中：θ為解x的彎曲角；θL，θR分別為解x與左、右邊相鄰解xL，xR形成的夾角，其計算公式為

(5)

(6)

Pareto前沿曲線上彎曲角最大的點為曲線最凸點(edge-knee point)，該點可以認為是Pareto解集中的協調最優解[12-13]。

圖1 多目標優化問題的Pareto前沿Fig.1 The Pareto front of muti-objective optimization problems

2 NSGA-Ⅱ多目標優化算法

本文采用帶精英策略的快速非支配排序遺傳算法(Non-dominated Sorting Genetic Algorithm-II，NSGA-Ⅱ)來求解有限元模型修正的多目標優化問題，NSGA-Ⅱ算法是Deb等[14]于2000年在NSGA算法的基礎上提出的，是一種經典的多目標優化算法。NSGA-Ⅱ算法相對于NSGA算法具有以下優越性[15]：① 計算復雜度由O(mN3)降低為O(mN2)，其中m為目標函數個數，N為種群規模；② 采用擁擠度和擁擠度比較算子，克服了NSGA中需要人為指定共享參數的不足，保證了種群的多樣性；③ 引入精英策略，擴大采樣空間，避免優秀個體流失，改善算法的魯棒性和收斂速度。NSGA-Ⅱ算法具體步驟如下：

步驟1配置算法中所需參數，如種群規模，最大進化代數，交叉概率，變異概率等；

步驟2隨機產生種群規模為N的初始種群P0，并對P0進行非支配排序，每個個體被賦予非支配序值，令t=0；

步驟3通過二元錦標賽法從Pt選擇個體，并進行交叉和變異操作生成子代種群Qt，其種群規模也為N；

步驟4將父代種群Pt和子代種群Qt歸并形成規模為2N混合種群Rt；

步驟5對種群Rt進行快速非支配排序，根據個體的非支配解水平得到Rt中的全部2N個個體的非支配層F1,F2,F3,…；計算每一非支配層的個體局部擁擠距離，根據非支配關系以及個體的擁擠度提取N個個體作為新的父代種群Pt+1；

步驟6對種群Pt+1執行選擇、交叉、變異操作形生新子代種群Qt+1；

步驟7如果終止條件成立，則算法終止；否則跳轉到步驟4。

3 基于多目標優化的混凝土斜拉橋模型修正

3.1 工程背景

以我國某公路斜拉橋為研究對象，該橋主橋結構為(65+123+318+130)m高低塔雙索面預應力混凝土斜拉橋，采用塔梁分離、索塔下橫梁上設置豎向支座的半漂浮體系。該橋橋寬27.5 m，主梁斷面為混凝土雙肋板式截面(Π型梁)，采用C55混凝土，主梁高2.6 m，頂板厚0.28 m。橋塔采用“H”型索塔，采用C40混凝土，高塔高115.4 m，低塔高91.9 m。斜拉索采用直徑7 mm平行鋼絲斜拉索，高塔44根索，低塔32根索，采用扇形布置，除了低塔邊尾索，主梁上的索距都為8 m。高塔側邊跨設有輔助墩和過渡墩，墩身材料為C30混凝土。大橋總體布置如圖2所示，大橋成橋照片如圖3所示。

(a) 立面

(b) 平面圖2 橋梁總體布置圖(單位：cm)Fig.2 The engineering drawing of the bridge(unit: cm)

圖3 橋梁成橋狀態Fig.3 The photo of the completed bridge

3.2 初始有限元模型的建立

斜拉橋屬大型復雜結構，初始有限元模型的建立以橋梁圖紙和相關規范為依據，著重對結構的剛度、質量和邊界條件的精確模擬，盡量減少模型階次誤差。節點單元劃分時綜合考慮靜載試驗加載位置、位移和應變測點位置等信息，便于荷載的施加和結果的提取。

(1)建模策略及單元選取。有限元建模策略會影響計算精度，因此需要采取合適的建模策略。斜拉橋主梁常用的模擬方式主要有：平面梁模式、空間梁模式、梁格模式、空間板殼模式、實體單元模式等。前三種模式屬于簡化處理方法，均需對主梁的質量、剛度等進行等效換算，結構較復雜時計算精度難以控制；空間板殼模式應用于箱形截面梁等薄壁結構時效果較好，但建模較復雜；實體單元模式既能精確模擬結構的空間效應還能分析主梁的局部效應，具有較高計算精度，但建模較復雜，且需要較多的計算機時。由于本大橋主梁為Π型開口截面梁，翹曲剛度對結構扭轉頻率和振型有較大的影響，為了精確模擬Π梁的空間受力行為，采用實體單元Solid65建模；索塔采用空間梁單元Beam188建模；斜拉索采用設置為只受拉模式的空間桿單元Link10模擬，采用彈性模量按Ernst公式修正的方法計入斜拉索的垂度效應影響；為了精確反映結構的真實狀態以及便于模型參數的修正，防撞墻和橋面鋪裝采用殼單元模擬，并與橋面主梁實體對應節點耦合，考慮其質量特性和空間尺寸，忽略其剛度。

(2)邊界條件。采用節點耦合方法模擬斜拉索與索塔和主梁間的連接；索塔、輔助墩、過渡墩與主梁間支座采用彈簧單元模擬，其剛度初始值按照支座剛度換算得到；橋臺處按支座剛度換算值施加節點彈性約束邊界條件；索塔、輔助墩、過渡墩底部采用固接約束邊界條件。

按照上述的建模策略、選取的單元和邊界條件，采用通用有限元軟件ANSYS建立斜拉橋的精細化初始有限元模型，模型如圖4所示，初始有限元模型的材料特性取值如表1所示。表1中的彈模和容重參數按照設計圖紙和相關規范取值；對于支座剛度參數，相同規格的支座取相同的剛度，其初始值按支座剛度換算得到，表1中列出了背景斜拉橋3類支座橫向剛度參數。

圖4 橋梁有限元模型圖Fig.4 The finite element model pictures of the bridge

表1 初始FEM材料特性取值Tab.1 Values of initial FEM material properties

3.3 現場靜動力試驗

大橋竣工后，按照在相關標準和規程指導下編制的試驗方案，對大橋進行了全橋靜載試驗和模態測試。靜載試驗主要包括主梁低塔側邊跨最大正彎矩、主梁中跨最大正彎矩、主梁輔助墩墩頂處最大負彎矩和主梁橫梁最大正彎矩等加載工況，每個加載工況都對主梁控制截面的應力、撓度，塔頂偏位，相關斜拉索的索力等進行了測試；在靜載試驗正式加載前對大橋進行了預壓，以消除結構非彈性變形的影響；根據大橋設計荷載標準，采用荷載等效的原則布置實際加載車輛，實際加載選用總重300 kN的自卸卡車(前軸重約6 t，中后軸重約24 t，前中軸距離3.6 m，中后軸距離1.4 m)加載，靜載試驗位移測點布置如圖2所示。

在環境激勵情況下進行橋梁模態試驗，橋梁振動加速度信號由傳感器拾振并通過放大器放大再由采集儀采集，采用941B型拾振器，模態試驗測點布置如圖2所示，主梁共布置82個測點，上、下游對稱布置，采樣頻率為50 Hz，振動信號由結構模態分析軟件MaCras進行采集操作和處理，得到大橋的前13階自振頻率和振型。斜拉橋靜動力測試結果見表2、表3、表5、表6。

3.4 模型修正及驗證

3.4.1 待修正參數選取

模型修正參數選取是在充分考慮對模型有影響的不確定性因素下的所有可能存在誤差的參數，包括結構的幾何參數、材料特性參數和邊界條件參數等，參數的選取往往存在所選的參數對特征信息不敏感、漏選、誤選等問題，因此待修正參數的選取要遵循兩個原則：① 選擇的參數不能太多，過多的待修正參數易造成優化規模大求解困難，或使參數修正值不具有物理意義；② 選擇靈敏度比較大的參數，如果選擇靈敏度非常小的參數會造成靈敏度矩陣病態，導致優化不收斂或得到錯誤的求解結果。靈敏度分析可以研究系統輸出特征量對系統參數的敏感性，在參數型模型修正中，n維模型參數向量P(p1,p2,p3,…,pn)和m維結構響應向量f(f1,f2,f3,…,fm)是隱式函數關系，可表示為

fi=fi(p1,p2,p3,…,pn)

(7)

利用差分法計算各指標對各參數正則化后的相對靈敏度如下

(8)

結合背景混凝土斜拉橋的實橋情況及初始有限元模型的特點,根據工程經驗選擇初始待修正參數及其初始值如表1所示，再對初始待修正參數進行靈敏度分析，彈模和支座剛度等剛度參數對結構自振頻率的靈敏度分析如圖5所示，容重參數對結構自振頻率的靈敏度分析如圖6所示。根據靈敏度分析結果，剔除靈敏度相對較小的參數，最后篩選出靈敏度相對較大的主梁材料彈性模量、主梁材料容重、橋塔材料彈性模量、橋塔材料容重、斜拉索材料彈性模量、輔助墩和過渡墩材料彈性模量等6個參數作為待修正參數。

圖5 剛度參數靈敏度分析Fig.5 Sensitivity analysis of stiffness parameters

圖6 容重參數靈敏度分析Fig.6 Sensitivity analysis of bulk density parameters

3.4.2 目標函數

大跨度斜拉橋在靜力荷載作用下結構位移的測試方便，且測試數據的精度和可靠性能得到保證，根據試驗實測位移和模型計算位移可以構造靜力位移目標函數F1為

(9)

式中：Uaji，Utji分別為第j工況的i點處的靜力位移的有限元模型計算值和靜力試驗實測值；γj為第j工況的權重系數，本文選取背景斜拉橋主梁中跨最大正彎矩工況作為模型修正實測數據來源工況，因為在此工況下，主梁邊、中跨位移值均較大，測試誤差帶來的影響相對最小。選取主梁低塔側邊跨最大正彎矩、主梁輔助墩墩頂處最大負彎矩兩個工況為模型修正后的驗證工況，因為這兩個工況都是邊跨局部加載，中跨位移值相對較小，容易受到測量誤差的影響，而局部較大的實測數據適合用于驗證模型修正效果。

對于基于動力測試的目標函數，由于振動頻率對結構剛度變化非常敏感，可利用試驗實測振動頻率和模型計算振動頻率構造目標函數F2為

(10)

式中：fai，fti分別為振動頻率計算值和實測值；n為參與優化的實測頻率階數，為了優先保證用于修正的實測數據的豐富性，選取大橋的前11階實測動力特性為模型修正實測數據來源；選取第12和第13階頻率及振型為修正后模型驗證的數據來源，這是因為高階頻率和振型更能驗證模型計算的準確性。

聯合靜動力的有限元模型修正目標函數表示為

min{F1,F2}=

(11)

需要指出的是，有限元模型的計算振型階次與試驗實測振型階次并不一定完全吻合，修正過程中為了保證計算模態與實測模態的振型相互匹配，采用計算模態與實測模態的模態置信準則(Modal Assurance Criterion, MAC)進行振型配對，其計算公式為

(12)

式中：φa，φt分別為計算模態與實測模態的振型向量。MAC值介于0～1，MAC值越接近于1則表示二者的相關性越好。本文計算模態與實測模態的振型配對時的MAC界限值取0.8，即當二者的MAC≥0.8時可認為二者對應同一階頻率的振型。

3.4.3 模型修正過程

本文基于多目標優化的混凝土斜拉橋靜動力有限元模型修正的本質是以構建的包含靜動力響應的目標函數為優化目標，采用帶精英策略的快速非支配排序遺傳算法(NSGA-Ⅱ)進行多目標優化，尋求非劣解。模型修正具體實施是通過在數學軟件MATLAB和有限元分析軟件ANSYS中自編程序將NSGA-Ⅱ多目標優化算法和有限元分析結合起來實現優化求解。ANSYS的批處理啟動模式(Batch模式)與MATLAB的強大運算和集成功能使得二者的聯合應用成為可能。在Batch模式下，ANSYS只在后臺啟動并運行指定的APDL命令文件，并輸出指定結果，而不用進入交互界面操作；ANSYS分析完成后，MATLAB程序再調用ANSYS的分析結果，從而實現二者的有機結合；所以模型中ANSYS的輸入、計算和輸出全部采用APDL語言編程完成。

首先在ANSYS中建立混凝土斜拉橋結構的參數化有限元模型，并對模型進行初始化，設置待修正參數的具體數值及取值范圍；然后對建立好的參數化有限元模型進行分析，包括靜力分析和模態分析，并提取優化計算需要的分析結果并寫入到輸出文件，比如結構的位移值、應力值、振動頻率值和振型向量等；隨后將以上輸入、輸出文件讀入MATLAB，并進行多目標優化的相關計算與操作，通過對優化算法的控制，判斷優化進程是否滿足終止條件，完成對混凝土斜拉橋靜動力有限元模型修正多目標優化的求解，混凝土斜拉橋靜動力有限元模型修正的流程如圖7所示。多目標優化參數設置：種群規模為60，最大進化代數為80，交叉概率為0.9，變異概率為0.1。初始種群如圖8所示，最終形成的Pareto前沿及Pareto協調最優解如圖9所示。得益于NSGA-Ⅱ算法良好的多目標尋優功能，解集的收斂性和穩定性較好，多次運行的得到的Pareto前沿基本相同，Pareto解集中的協調最優解相差在2%以內。

圖8 初始種群在解空間情況Fig.8 Initial population in objective space

圖9 Pareto最優解Fig.9 Pareto optimal solutions

3.4.4 模型修正結果

模型修正前后的靜力位移及實測值如表2所示，可以得出，模型修正后各測點的實測撓度與計算撓度更加吻合，實測與計算的最大位移均發生在11號測點，模型修正前該點的計算與實測撓度相差15.2%，修正后相差減小到2.8%。模型修正前后的動力特性及實測自振頻率如表3所示，可以看出，修正后的各階計算自振頻率與實測值偏差明顯減小，修正后的豎向彎曲頻率計算值與實測值相差最大為二階豎彎，相差1.38%；修正后的橫向彎曲頻率計算值與實測值相差最大為三階橫彎，相差2.83%；修正后的一階扭轉頻率計算值與實測值相差1.71%；可見，修正后豎向、橫向及扭轉振動頻率都與實測值均較為吻合，各階的MAC值也有一定的改善，修正后的MAC>0.85。綜上可得，基于多目標優化的靜動力有限元模型修正應用于混凝土斜拉橋實橋的有限元模型修正能夠取得較好的效果。

模型修正參數的初始值與修正值見表4，其中混凝土主梁和主塔的彈性模量值的修正幅度較大，主要是因為混凝土材料的彈性模量本身具有較大離散性，當選材品質好、施工質量較好時，容易出現彈性模量值大于設計值的情況。此外，對于主梁，由于初始模型未考慮梁內配有較多的普通鋼筋和預應力鋼筋，以及未計入橋面鋪裝和混凝土防撞墻對主梁剛度的提高作用，這也是主梁的彈性模量修正幅度較大的重要原因；對于混凝土主塔，混凝土內部配置強大鋼結構勁性骨架和大量的普通鋼筋，而初始模型未考慮這些因素影響，這也是其彈性模量值修正幅度較大的另一個重要原因。

表2 模型修正前后靜力位移值比較Tab.2 Comparisons of static displacements before and after model updating

3.4.5 模型修正效果驗證

通過上述的模型修正，使得大橋參與修正的計算靜力位移和模態結果與實測結果都已吻合較好，如果修正后的模型能更加真實的模擬實際結構，那么未參與修正的靜動力實測值與相應修正后模型的理論計算值會更加接近；若兩者差別沒有減小，則說明修正后的有限元模型還不能全面真實模擬實際結構，還有其他產生誤差的因素未考慮到或修正方法還存在問題。采用靜力試驗的主梁低塔側邊跨最大正彎矩工況(工況二)和主梁輔助墩墩頂處最大負彎矩工況(工況三)兩個工況為驗證工況，修正前后的靜力位移比較如表5所示，可得模型修正后兩個工況的靜力位移計算值都與實測值更為吻合。采用大橋的第12和第13階自振頻率對動力修正效果進行驗證，修正前后的動力結果比較如表6所示，可見修正后的這兩階頻率的計算值與實測值更為接近。驗證表明，未參與修正的靜力位移和動力計算結果也與實測值更加接近，修正后的有限元模型能夠較準確全面的反映實際結構的靜動力特性。

表3 模型修正前后動力特性比較Tab.3 Comparisons of dynamic properties before and after model updating

表4 模型修正前后參數值比較Tab.4 Comparisons of parameter values before and after model updating

表5 模型驗證的靜力位移比較Tab.5 Comparisons of static displacements for model verification

表6 模型驗證的動力特性比較Tab.6 Comparisons of dynamic properties for model verification

4 結 論

針對混凝土斜拉橋結構，提出一種基于多目標優化的靜動力有限元模型修正方法，并將這種方法運用于國內某混凝土斜拉橋的有限元模型修正，主要得到以下幾點結論：

(1) 本文以靜動力實測數據構建目標函數，運用了較全面的實測信息，獲得的修正結果能更全面的反映混凝土斜拉橋的靜動力特性，比單一運用橋梁的靜力試驗或者動力試驗響應獲得的修正結果更為可靠。

(2) 采用三維實體單元模擬背景斜拉橋的Π形截面混凝土主梁，最大進化代數取為80代即可得到較好的多目標優化結果，修正后主梁的一階扭轉頻率計算值與實測值相差1.71%，二階扭轉頻率計算值與實測值相差1.42%，表明采用三維實體單元模擬混凝土斜拉橋Π形主梁能較真實地反映實際結構的空間效應，減少模型固有的離散誤差，取得更好的修正效果。

(3) 采用帶精英策略的快速非支配排序遺傳算法(NSGA-Ⅱ)對國內某混凝土斜拉橋有限元模型進行了多目標優化修正，獲得了分布較均勻的Pareto非劣解集，在此基礎上，利用最大彎曲角法得到了Pareto解集中的協調最優解，最后利用未參與修正的靜動力實測數據對修正后的有限元模型進行了驗證。結果表明：基于多目標優化的靜動力有限元模型修正方法不但能避免現有的聯合靜動力模型修正方法中將多個子目標函數組合轉化為單目標函數時人為確定權重系數給模型修正帶來的主觀性影響，而且模型修正效果也較為理想。模型修正后不僅參與優化的靜力位移和動力計算結果與實測值更加吻合，而且未參與修正的靜力位移和動力計算結果也與實測值更加接近，表明修正后的有限元模型能夠較準確全面的反映實際結構的靜動力特性，可作為該橋的基準有限元模型。