基于改進萬有引力搜索算法的南中環橋模型修正

秦世強， 甘耀威， 康俊濤

(武漢理工大學 土木工程與建筑學院， 武漢 430070)

建立一個準確、能夠代表實際結構行為的有限元模型對橋梁狀態評估和監測有著重要的作用。然而，按照設計圖紙建立的初始有限元模型，由于模型簡化、結構設計參數的不確定性、結構退化或損傷等因素的影響[1]，往往難以代表實際結構行為。因此，需要對初始有限元模型進行修正，以降低結構響應的有限元計算值與實測值之間的相對誤差。有限元模型修正的本質是一個優化問題，其目標函數是由結構響應有限元計算值與實測值的殘差函數構成。對于橋梁結構模型修正的目標函數，通常具有多維、高度非線性和多個局部極值點的特征；傳統搜索算法如梯度算法、牛頓迭代法在處理這類問題時，容易陷入局部最優，難以獲得目標函數的全局最優解[2]。近年來基于智能優化理論的粒子群算法[3]、遺傳算法[4]、猴群算法[5]、萬有引力搜索算法(gravitational search algorithm ，GSA)[6]等迅速發展，相比傳統搜索算法，智能進化算法在處理非線性、多極值等問題時具有一定的優勢。Alkayem等[7]對基于智能算法的有限元模型修正進行了全面的綜述，總結了各類單目標優化算法和多目標優化算法在有限元模型修正中的應用現狀；同時，作者指出提升智能算法的搜索精度和效率，仍是模型修正和損傷識別領域一個值得研究的方向。

GSA[8]是Esmat Rashedi等提出的一種元啟發式算法，是一種模擬物理學中萬有引力的新的優化算法。它通過種群中粒子間的相互作用來指導群體進行智能優化搜索；已有研究表明[9]：GSA的收斂性明顯優于粒子群算法、遺傳算法等仿生優化算法。但是GSA依然存在著容易早熟、后期迭代速度慢、不易跳出局部最優解等問題。針對GSA的上述問題，已有文獻提出了一些改進措施。谷文祥等[10]提出一種新的局部搜索算法并結合GSA求解流水線調度問題，有效的避免算法陷入局部最優值，提高了解的質量，驗證了所得算法的有效性；李超順等[11]結合粒子群算法的運動特點提出了改進的引力搜索算法，并將其應用到勵磁控制系統PID(proportional-integral-derivative)參數優化，并與傳統群體優化算法進行了充分對比試驗，驗證了其優化算法的有效性；李沛等[12]在GSA的速度更新部分引入粒子群算法中的記憶和群體信息交流功能，并提出了基于權值的粒子慣性質量更新公式，驗證了所提算法能有效規劃出無人機的最優航路；蔣建國等[13]針對基本GSA在處理優化問題時會出現發散的情況,通過限制粒子的速度同時更改算法中的參數，并將改進的GSA運用到邊坡穩定性分析中，搜索出臨界滑動面并結合極限平衡法計算出相應的最小安全系數，并驗證了其具有更好的優化性能。本文在上述研究的基礎上，結合遺傳算法提出一種改進的GSA，利用荷載試驗結果并結合Kriging代理模型，實現了一座鋼混疊合梁組合式系桿拱體系橋梁的模型修正。

1 基于Kriging模型的模型修正

有限元模型修正的本質是通過改變模型中的結構設計參數，減小有限元計算值與實測值之間的誤差。因此，模型修正的目標函數可以定義為

(1)

式中：x為待修正的設計參數；ri(x)為第i類結構響應計算值與實測值之間的殘差函數；m為模型修正中考慮的結構響應的類型種數；αi為不同殘差函數的權重系數。

為獲得式(1)所表達的目標函數在設計空間內的最小值，一般需通過迭代法求解。然而，由于結構響應與設計參數間隱函數的關系，使得每次迭代均需要調用有限元模型來計算殘差函數，這一過程十分耗時，特別是對復雜的橋梁結構有限元模型而言。基于代理模型的模型修正方法因此得以發展；常用的代理模型包括多項式響應面[14]、徑向基函數[15]、Kriging模型[16]等。有研究表明[17]：Kriging模型不僅可以給出設計空間內任意一點的預測值，同時能給出該點的預測誤差。基于Kriging模型的模型修正的主要流程可分為3個步驟:① 需要用初始有限元模型得到結構的響應和設計參數樣本，可通過試驗設計[18]完成，常用的試驗設計有中心復合設計、D最優設計、均勻設計和拉丁超立方設計等；② 利用得到的響應和設計參數樣本生成Kriging模型并做精度測試；③ 基于優化算法尋找目標函數在設計空間內的最小值。

Kriging代理模型是由確定性的回歸模型和隨機過程兩部分組成，結構響應矩陣Y(x)可以用Kriging模型表達成設計參數x的表達式，公式為

(2)

2 基本萬有引力搜索算法

GSA是一種模擬萬有引力定律的智能搜索算法。GSA中每個粒子都可以看作一個個體，每個粒子都具有自身的位置和移動速度，并且可以保留迄今為止粒子所能找到的最優位置以及所有粒子當前能找到的全局最優位置。在標準的GSA中，粒子i和粒子j之間的引力可表示為

(3)

式中：pd,i(t)為在t時刻第i個個體在d維空間中所在的位置；Mpi(t)為t時刻受力個體i的質量；Maj(t)為t時刻施力個體的質量；ε為一個很小的常數；Rij(t)為粒子i與j之間的歐幾里得距離，可表示為

(4)

G(t)為在第t時刻的萬有引力常系數，可表示為

(5)

式中：G0為常系數的初始值；α為下降系數；T為總迭代次數。因此，第i個個體在d維空間中的合力可表示為

(6)

式中：N為粒子種群的個數；rand則為0到1的隨機數。t時刻時粒子i的加速度可表示為

(7)

式中，Mii(t)為t時刻粒子i的慣性質量。粒子i在d維空間中的移動速度以及位置的更新可以表示為

Vd，i(t+1)=rand×Vd，i(t)+accd，i(t)

(8)

pd，i(t+1)=pd，i(t)+Vd，i(t+1)

(9)

式中，Vd，i(t)，pd，i(t)分別為第t次粒子移動速度和位置。每個粒子的質量是通過第t次的個體適應度計算得到，可表示為

(10)

式中：fiti(t)為第t時刻個體的適應度；worst(t)和best(t)分別為第t時刻種群中個體最差和最好的適應度。假設個體的引力質量、慣性質量和個體的質量相等，則有

Mai=Mpi=Mii=Mi,i=1,2,…,N

(11)

(12)

可知當粒子的適應度較優時，粒子的質量就較大，自身的移動速度就會較慢，會不斷吸引適應度較差的粒子向適應度更好的位置移動。

3 改進GSA

在GSA的過程中，當兩個粒子之間的距離較近時，粒子就可能陷入局部最優解并且很難跳出局部最優值，從而降低了搜索效率。本文結合遺傳算法提出一種改進的GSA。首先通過對初始種群適應度計算，按照適應度的大小對初始種群中的粒子進行排序。然后將初始種群分為優解組和劣解組兩組，對優解組進行隨機交叉的操作，對劣解組進行變異的操作；將交叉變異得到的新種群與當前迭代中的初始種群進行重新排序，取適應度值最好的前一半種群作為本次迭代后的更新種群；并且引入了概率控制，在整個算法初期，主要進行全局搜索，加大搜索范圍，在后期則進行局部挖掘，增加搜索結果的精度，通過控制進行交叉變異操作的概率，能夠有效地發揮全局搜索和局部挖掘的能力。

首先生成種群數量為N的初始種群，計算所有粒子的適應度值fitness(t)并且進行排序，然后根據適應度將適應度較好的一半粒子放入到優解組當中，再將另一半較差的粒子放入到劣解組當中，將優解組進行遺傳交叉操作，用優解組的父代種群產生優良的子代種群，并且在本改進算法中進行交叉操作的兩個粒子是隨機選取，可表示為

(13)

劣解組采用變異操作，通過變異操作產生優良的子代，充分利用每個個體；變異可表示為

(14)

對改進算法流程描述如下，改進算法的流程圖，如圖1所示。

圖1 算法流程圖

步驟1設置算法初始參數。

步驟2隨機生成初始種群。

步驟3計算初始種群的適應度，初始化個體最優值和全局最優值。

步驟4根據概率控制判斷是否進入遺傳交叉和變異，若rand>p則進入步驟5，否則進入步驟7。

步驟5進行分組策略，根據適應度優劣排列粒子順序，將適應度較好的一半粒子放入優解組進行遺傳交叉操作，將適應度較差的另一半粒子放入劣解組進行變異操作。

步驟6進行精英策略，將優解組產生的子代和劣解組產生的子代合并為一組，再與當前迭代次數下的初始種群建立一個擴大種群，將擴大種群按適應度排序，選取適應度較好的一半種群作為精英種群參加下一輪迭代。

步驟7用式(8)、式(9)更新粒子的速度與位置。

步驟8更新粒子個體最優值，全局最優值。

步驟9判斷是否滿足迭代終止標準，若算法迭代次數達到最大迭代次數T，則終止算法，否則跳回步驟3。

4 算法仿真測試

4.1 參數設置

為了驗證改進GSA的效果，基于測試函數對粒子群算法、GSA和改進GSA進行對比分析。其中粒子群(particle swarm optimization，PSO)算法是由Kennedy等[20]提出的模擬鳥群覓食行為的一種群體協作優化方法。在參數設置上，為保證各算法可比性，3個算法種群規模均取30，函數維度均取10，最大迭代次數取1 000。對于PSO，其學習因子c1和c2分別為1.5和0.5，慣性常數取0.726，上述參數是標準PSO的默認參數，已被證明適用于各類優化問題；對于GSA和本文算法，其萬有引力系數G0皆取10，下降系數α皆取20。

4.2 測試函數

選取5個函數標準測試函數，將本文改進的算法與PSO和GSA算法進行比較。其中f1為單峰函數，主要用來測試算法的精度與收斂速度；f2為多峰函數。當兩個函數維度取2時，函數圖像如圖2所示。兩個函數用來測試算法的收斂速度和精度，以及跳出局部最優解和全局搜索的能力。兩個測試函數的全局最優值均為xi=0,(i=1∶n),f(x)=0。為了消除隨機性可能對算法造成的影響，每個測試函數都用3種算法運行了100次，試驗所處環境為：MATLAB R2018a，CPU為i5-9400F@2.90 GHz。以均值和方差作為算法的評判標準。均值表示算法的精度，方差則能反映得到結果的穩定性，收斂曲線能反映算法的收斂速度。f1為單峰值球函數，所以結果越接近0就表示越精確，故設置成功率來顯示其優化的效率，當適應度小于1×10-20，就認為搜索成功；f2函數在最優值附近還有許多局部極小值，故當適應度分別小于1×10-10則認為搜索成功。

(a) Sum of different power函數圖像

(1) Sum of different power函數：x∈[-1,1]

(15)

(2) Ackley’s Path函數：x∈[-1.5,1.5]

(16)

4.3 仿真結果分析

為了比較改進GSA，PSO和GSA算法的優劣，圖3列出了兩種測試函數的平均適應度曲線圖，各算法100次搜索的統計結果，如圖3所示。從平均適應度曲線的數據可以看出，改進的優化算法在精度上都要優于PSO和GSA算法，改進算法收斂時的結果對于f1測試函數都要小于1×10-20，對于測試函數f2其收斂時的精度也明顯高于另外兩種算法。

(a) f1

對于f1函數改進的優化算法得到的適應度相比其他算法其探索點的精度顯著提高，在優化的成功率上，其達到數量級1×10-15以下的成功率為100%遠遠高于其他優化算法。對于測試函數f2改進算法探索到最優位置的成功率接近100%。且測試函數f2其全局最優值相對于局部極小值的適應度相差很小且局部最優值數量較多，因此很難找到最優位置，改進GSA的成功率有98%，其他算法則很難找到最優值；這表明改進GSA顯著提升了跳出局部最優解的能力。將3種算法相對于兩種測試函數的計算效率進行了對比，列出了完成一次尋優操作所需時間具體時間，如表1所示。可以看到，相比基本GSA和PSO算法，由于改進GSA算法加入了交叉變異等操作；因此，其效率上不占優勢，但也未顯著增加計算時間。

表1 計算效率對比

5 南中環橋模型修正

5.1 工程概況

南中環橋是連接太原市東西城區的重要橋梁，如圖4所示。該橋主橋結構體系為組合式系桿拱、鋼管混凝土蝶形拱橋，跨徑布置為(60+180+60) m。主橋跨中為鋼混疊合梁段，全長92 m，主梁高為3 m；主橋邊跨60 m及中跨44 m梁段，采用的是混凝土箱梁，梁高為3.0～6.5 m。主橋標準斷面寬度53 m，跨中寬度為71 m。主跨是由4根主副拱肋組合而成的下承式鋼管混凝土拱橋，主拱肋外傾16°，矢跨比為1/4.326；副拱肋外傾26.882°；主拱肋和副拱肋之間利用拉桿連接。

圖4 南中環橋

5.2 荷載試驗概況

南中環橋荷載試驗分為靜載試驗和動載試驗兩個部分，靜載試驗主要測試橋跨結構在試驗車輛靜載作用下的位移；動載試驗主要測試結構的固有振動特性和沖擊系數。限于篇幅，僅介紹與模型修正相關的試驗工況。模型修正共選擇了兩個靜載試驗工況(A-A、B-B)，測試在車輛下主梁和主拱關鍵截面的位移。圖5為橋梁立面圖及位移測點布置圖。靜載試驗車輛為40 t三軸載重汽車，軸距為(1.8+3.3) m。動力測試部分，為獲取橋梁的固有振動特性，利用加速度傳感器測試主梁在環境脈動荷載下的三向加速度響應。在主梁縱向平均每14 m布置一個加速度測點，采樣頻率為80 Hz。利用測試得到的加速度響應，結合模態識別算法可以獲得橋梁的頻率、阻尼比和振型。靜力位移、頻率的試驗值和計算值的對比結果，如表2和表3所示。從表2和表3可知，位移最大相對誤差達28.19%，頻率的最大相對誤差達到9.51%；其中，計算值是根據Midas Civil初始有限元模型計算得到的。該模型按設計圖紙建立，全橋共有8 322個單元，6 118個節點，其中主副吊桿和斜拉桿采用桁架單元，其余部分皆用梁單元模擬。主梁采用梁格法模擬，鋼管混凝土拱采用聯合截面模擬；吊桿張拉時考慮幾何剛度貢獻，邊界條件按設計圖紙設置。

表2 兩種工況下各測試截面的位移計算值和試驗值

表3 橋面系實測自振頻率及其計算頻率

圖5 橋梁立面圖及位移測點布置圖(cm)

主梁前兩階豎向振型計算值和試驗值的對比，如圖6所示。對應的模態置信值(modal assurance criterion, MAC)分別為0.946和0.963；這表明試驗振型與計算振型吻合良好，反映試驗數據可靠。上述分析表明：初始有限元模型對結構真實剛度分布模擬不夠準確，需要通過模型修正，從而更好地代表實際結構行為。

(a) 主梁二階振型

5.3 Kriging模型建立

為了提高模型修正迭代效率，利用Kriging模型替代有限元模型預測計算值。首先選擇修正參數，然后利用拉丁超立方抽樣獲得參數在設計空間的分布，最后建立Kriging模型并評估其精度。

在修正參數選擇方面，因為該橋為預制構件所以截面誤差一般較小，故截面慣性矩不作為修正參數。在建立有限元模型時對于整個橋梁進行了簡化，并沒有建立節點板、螺栓等構件，所以在結構的質量分布上會存在影響，而一些細部的加勁肋則會影響到模型的整體剛度，因此將參數的選擇范圍定在各個構件的質量密度和彈性模量上。本文分別選擇了鋼管拱主拱、鋼管拱副拱、鋼混疊合梁段、混凝土段和挑臂5個構件的彈性模量和質量密度共10個參數進行敏感性分析。結合敏感性分析結果，最終選取鋼管拱主拱彈性模量E1、鋼管拱主拱密度ρ1、鋼混疊合梁段彈性模量E2、鋼混疊合梁段密度ρ2、混凝土段主梁彈性模量E3和混凝土段主梁質量ρ3這6個參數作為待修正參數，其初始值分別為2.06×108kN/m2，78.5 kN/m3，2.06×108kN/m2，78.5 kN/m3，3.45×107kN/m2，26 kN/m3。

采用拉丁超立方抽樣法對選定的修正參數進行試驗設計，對6個設計參數進行了81次抽樣，參數的抽樣范圍設定在±30%以內，再將設計得到的修正參數分別代入初始有限元模型中得到對應工況下的位移以及頻率響應。以樣本點作為輸入，節點位移值及模態頻率響應值作為輸出構建Kriging模型。選取鋼管主拱彈模E1、鋼混疊合梁段彈模E2兩個設計參數和該橋有限元模型一階頻率繪制三維Kriging響應面圖。一階頻率對于修正參數E1和E2的Kriging響應面模型以及對應的均方誤差(mean square error，MSE)，如圖7所示。其中曲面為Kriging模型響應面，散點為樣本點，可知樣本點基本位于曲面上。從對應的MSE圖可以看出，誤差基本都不超過3×10-9，說明Kriging精度較高，可替代有限元模型進行結構響應預測。

圖7 一階頻率與E1、E2之間函數擬合圖

5.4 數值驗證

考慮到實際結構損傷位置是未知的，為了證明模型修正方法的有效性，首先利用南中環橋有限元模型進行數值驗證。給參數E1,E2,E3設定數值損傷，其中：相較于初始值，參數E1,E2設定+10%的數值損傷，參數E3設定-10%的數值損傷，以未設損傷的模型代表初始有限元模型，施加損傷后的模型為真實的有限元模型。以3個設計參數作為輸入，以前6階頻率作為輸出，建立Kriging模型，然后以本文中的3種優化算法對建立的Kriging模型進行尋優，PSO算法仍取其默認參數；GSA和改進GSA算法的萬有引力系數取為10，下降系數為20；各算法的種群規模均為30，最大迭代次數為1 000。其設計參數和頻率響應的結果,如表4、表5所示。由表可知，改進算法3個設計參數的修正偏差修正前最高為11.11%，修正后誤差都降到約0.1%；頻率的響應值在修正前最大偏差達到1.78%，修正之后改進算法的最大誤差僅為-0.024%。采用本文的改進算法結合Kriging模型的模型修正方法是有效的。其中改進GSA、GSA和PSO算法的計算時間分別為0.468 s、0.353 s和0.146 s。對于橋梁模型修正而言，每次迭代中最為耗時的部分是評估目標函數，即計算設計參數修改后的目標函數值，而算法所耗時間相對較少；因此，改進GSA的計算效率可以接受。

表4 參數修正前后對比

表5 頻率修正前后對比

5.5 基于實測數據的模型修正

以位移殘差和前4階頻率殘差構建目標函數，分別利用改進算法、GSA和PSO對目標函數進行尋優，算法在參數設置上，PSO算法仍取其默認參數；GSA和改進GSA算法的萬有引力系數取為10，下降系數為20；各算法的種群規模均為30，最大迭代次數為1 000。修正后的參數、頻率和位移結果分別如表6、表7和表8所示。

由表6可知，參數鋼管拱主拱彈性模量E1、鋼混疊合梁段彈性模量E2、混凝土段彈性模量E3利用改進萬有引力搜索算法的修正結果相較于初始值分別提高了13.27%，21.72%，24.37%，主要原因為在建立初始模型時對橋梁整體進行了簡化，在建立有限元模型時并沒有考慮局部加勁肋和橫隔板，并且鋼管拱以及鋼混疊合梁段在實際橋梁中數量較多，所以有限元模型比實橋的剛度較低；鋼管拱主拱密度ρ1的修正值都與實際值差別不大，都略微偏小；而鋼混疊合梁段密度ρ2的修正值都比實際值偏大，混凝土段質量ρ3的修正值則比實際值都較為偏小，原因是在混凝土施工質量存在一定的離散性，使分布不均勻，并在模型修正中體現出來。需要指出的是，實際橋梁參數修正前后的變化參數修正前后的變化，僅表示模型整體剛度和質量分布與實際情況的差異，并非材料參數發生了變化，即修正后的參數可理解為等效彈性模量和等效質量密度。并且通過查看設計圖紙可以得知鋼混疊合梁段鋼梁節段一的總質量為107 434.7 kg，而通過Midas Civil查詢單元質量可知初始有限元模型鋼梁節段一的總質量為103 025.2 kg，即有限元模型中的質量低估了約4.28%，這與鋼混疊合梁段密度ρ2的修正值3.77%較為接近。

表6 參數修正前后對比

由表7可知，利用改進GSA、GSA和PSO修正后頻率值的對比，其相對誤差皆為修正前后計算值相對于實測值的誤差。從整體看南中環橋的實測頻率比計算頻率都偏大，說明有限元模型的整體剛度都偏低，與參數的修正結果符合，修正后的頻率值都有所提升，說明修正結果可作參考。相比初始模型，各算法修正后的頻率相對誤差均有不同程度降低，改進GSA取得了相對于PSO和GSA更好的修正結果。

表7 頻率修正前后對比

由表8可知，3個優化算法修正后位移值的對比，除了節點1和節點3在修正前其計算值和位移值差距不大導致修正空間較小，其余7個節點都相對于修正前精度有明顯的提高。整體上有限元模型計算值與實測值相比都偏大，也證明了有限元模型在整體剛度上都偏小，修正后結果都比修正前更加接近實測值。修正前，實測值與計算值的最大誤差達到了28.19%，在修正后，改進算法所得結果的相對誤差為13.62%，在精度上得到明顯的提升，且GSA和PSO的修正結果相對誤差分別為16.39%和15.90%，改進算法也比這兩種算法修正精度要高。上述結果表明：經過模型修正，模型的位移值計算值和測試值更為接近，其相對誤差減小，表明了修正后更能反映實際情況。相對于GSA和PSO，改進GSA的修正精度更高。

表8 位移修正前后對比

6 結 論

通過試驗與有限元模擬，基于Kriging模型對南中環橋進行了有限元模型修正。為提高修正的精度，提出了一種改進的優化算法，得出了以下結論：

(1) 對于測試函數的修正結果，改進的優化算法對于普通GSA和PSO算法都有著較高的成功率和穩定性，從而提高了優化結果合理解的可能性；對于南中環橋的修正結果，改進的優化算法獲得更小的目標函數，且修正后頻率和位移整體相對誤差較小，剛度分布更為合理。

(2) 經過修正后位移和頻率的計算值與試驗值之間的誤差更小，其位移頻率除個別值外其誤差都有顯著降低，其中位移較大的兩個節點誤差分別從23.92%和28.19%降低到14.35%和13.62%，表明修正后的模型能更好代表實際結構。

(3) 因為實際結構的復雜性和試驗結果的不確定性，故模型修正的目標函數在修正參數的空間分布較復雜，本文提出的優化算法對修正結果的精度和穩定性有所提升。