基于改進螢火蟲算法的有限元模型修正

劉 綱，陳 奇，雷振博，熊 軍

(1. 山地城鎮建設與新技術教育部重點實驗室(重慶大學)，重慶 400045；2. 重慶大學土木工程學院，重慶 400045；3. 中鐵十一局集團第五工程有限公司，重慶 400037)

有限元分析技術具有應用成本低、計算過程直觀等顯著優點，已廣泛應用于各種復雜機械、航空及土木工程結構的設計、性能分析與損傷識別[1?2]。但實際結構的復雜性、邊界條件等的不確定性，使得有限元計算結果與實際結構響應往往存在一定差異[3]。為使有限元模型能夠更為準確地反映結構實際狀態，有限元模型修正技術得到了大力發展，已成為當前重點研究方向和前沿[4]。

目前，有限元模型修正技術多采用參數型修正法，即對有限元模型中的材料特性、截面尺寸等參數進行修正，具有簡單直觀、物理意義明確等特點[5]。修正步驟主要包括目標函數建立、修正參數確定及優化算法選取等，其中，合理算法的選取是保證修正精度的關鍵因素[6]。因為傳統的牛頓法、一階優化法等確定性搜索優化算法對初始值較為敏感，且易陷入局部最優[7]，故而近年來各類隨機搜索優化算法得到了更多的發展及應用。

李成等[8]提出了基于人工魚群算法的結構模型修正方法，以每層層間剛度為修正參數，完成了實驗室3層框架的模型修正；杜大華等[9]提出一種基于模擬退火算法的有限元模型修正方法，并以火箭發動機噴管模型為例實現了7參數的有限元模型修正；蘇越[10]提出了基于蟻群優化算法的有限元模型修正方法，并完成了虎門大橋的有限元模型修正；夏志遠等[11]提出了基于高斯白噪聲擾動的粒子群優化有限元模型修正方法，某在役橋梁結構模型修正結果表明，各階頻率誤差降低至8.86%以下。此外，差分進化算法、蜂群算法、遺傳算法等[12?16]也被大量應用到結構動力學模型修正中，形成了一系列基于仿生學的有限元模型修正方法，并在土木、機械等結構中顯示出良好的修正效果。

螢火蟲算法(firefly algorithm，FA)作為一種新興啟發式智能仿生算法，具有參數少、易實現、全局尋優能力強等特點，在狀態預測、傳感器優化布置、信號降噪處理等領域得到初步應用[17?21]。但標準FA算法存在著后期收斂速度慢、收斂精度低等不足，嚴重制約了FA在實際工程中的應用。針對這一問題，本文通過對隨機吸引度因子及自適應步長因子的改進，提升算法收斂效率與精度。并通過函數、簡支梁數值模擬和某實際連續剛構橋的有限元模型修正，驗證改進FA算法有效性，從而為有限元模型修正提供優化方法支撐。

1 螢火蟲算法及其改進

1.1 螢火蟲算法基本原理

螢火蟲算法是一種隨機搜索算法，其基本思路是借鑒螢火蟲的自然行為，將優化問題中的迭代搜索看作是螢火蟲相互吸引而自發聚集的過程，從而通過模擬螢火蟲群體行為實現最優化問題求解。在自然界中，螢火蟲個體的吸引力與發光強度正相關，與它們之間距離的平方負相關；亮度小的螢火蟲個體向亮度大的螢火蟲靠近，從而完成螢火蟲種群的自發式移動和聚集[17]，如圖1所示。圖1中，每個黑色原點為一只螢火蟲個體，代表優化問題的一個可能解；螢火蟲周圍的光暈表示其發光強度，代表該只螢火蟲的適應度函數值大小，適應度值越優，螢火蟲的發光強度越強、光暈范圍越大；兩只螢火蟲之間的虛線代表吸引力，虛線越粗，代表吸引力越大。

圖1 螢火蟲發光強度示意圖Fig. 1 Luminous intensity of firefly

螢火蟲的移動過程就是可能解不斷向最優解的迭代過程，下面給出該過程的數學描述。假設Xi=(x1i,x2i,···,xdi)是群體中的第i只螢火蟲，其中i =1, 2, 3, ···, N。N與d分別表示種群中螢火蟲的個數與修正參數的維度。對任意兩只螢火蟲Xi與Xj(i≠j)，定義兩者之間的歐氏距離rij為：

則兩者之間的吸引度因子為：

式中：β0為兩只螢火蟲距離rij=0時的吸引力，一般取為1；γ為光在空氣中的衰減系數，一般取為0.01～100。假設螢火蟲Xj(j≠i)的適應度值優于Xi，則螢火蟲Xi會被螢火蟲Xj吸引而更新，有：

式中：j=1, 2,···, i?1, i+1,···, N，上標t為當前迭代步數；α為步長因子，取值范圍為0～1；ε為d維隨機向量，服從正態分布N(0,1)。若Xi為當前最佳，則對Xi進行隨機移動：

螢火蟲Xi在第t步的具體更新過程如圖2所示。

圖2 第i只螢火蟲更新步驟Fig. 2 Updating steps of the ith firefly

1.2 螢火蟲算法改進

在式(3)中，螢火蟲Xi的更新過程可細分為跳躍操作dxβ和隨機游走操作dxα：

因此，該算法的優化性能主要取決于吸引度因子βij和步長因子α，若這兩者設置不當，會導致算法精度下降甚至難以收斂。針對這一問題，本文提出改進自適應參數控制算法，以提升計算效率及精度。

1) 隨機吸引度因子

吸引度因子βij在不同γ取值下隨著歐氏距離rij的變化規律如圖3所示。

由圖3知，γ較大時，吸引度因子βij及跳躍距離dxβ均會隨著距離增加快速趨近于0，致使搜索停滯；γ較小時，βij往往過早趨近于1，此時，跳躍操作可近似等效為：

圖3 βij迭代過程Fig. 3 The iterative process of βij

此時，亮度較低的螢火中直接跳轉到高亮度的螢火蟲附近，難以對狀態空間內其他位置進行探索、產生有效跳躍，導致發現率低、易陷入局部最優。因此，衰減系數γ大小決定了螢火蟲的可視范圍，在原螢火蟲算法中，衰減系數γ為定值，不會隨著迭代過程而改變，故無論γ取值如何，不可能同時兼顧迭代初期大范圍快速搜索、迭代后期小范圍精細搜索的需求。

針對這一現象，本文取γ=0.01，使各螢火蟲的βij大于0.8，擁有較大的初始搜索空間，避免計算停滯；同時，引入隔代隨機吸引度因子，如式(7)所示，使βij有機會取較小值，解決因γ取值過小導致βij過早趨于1而致使解的同質化現象，以提升算法遍歷性。

式中：rand表示從均勻分布U(0,1)中抽取隨機數；mod(t,2)=0表示迭代步數t為偶數。通過設置隔代隨機吸引度因子改善算法搜索空間。以某二維函數優化過程為例，跳躍操作改進前后搜索變化如圖4所示。由圖4可知，改進前由于βij過早趨近于1，跳躍距離dxβ沿固定線性路徑分布，搜索范圍單一。而改進吸引度因子有效增加算法搜索范圍，擴大了算法的搜索路徑，提升了算法遍歷性，提高了計算精度。

圖4 跳躍操作產生的距離dxβFig. 4 The distance dxβ produced by jump operation

2) 自適應步長因子

在螢火蟲算法迭代初期，隨機選取的初始值往往離真值較遠，則此時步長因子α應取較大值以保證獲得較廣的探索空間，而隨著算法不斷迭代優化，種群中所有個體將逐漸聚集，最終收斂到真值，即：

此時，α應趨近于0。否則，算法將在最優解附近大幅振蕩，影響算法的收斂性能。從以上分析可知，α值應隨著迭代演化而逐漸減小，但縮減的速度大小會影響最終計算效率與精度。縮減速度過大，會導致dxα過早趨于0，迭代過早停滯，無法收斂到最優值；縮減速度過小，會導致迭代后期的dxα過大，使結果在最優解附近振蕩。因此，本文提出參數α自適應動態更新策略：

式中：α0為步長因子初始值，建議取為1；T為最大迭代步數；c為常數縮減因子，可根據實際情況取值以獲取適宜的步長因子自適應縮減速度。根據算例實際計算結果分析，建議自適應步長因子α的縮減系數c取為15，此時，α及dxα變化規律如圖5所示。

圖5 α及隨機游走距離dxα變化曲線Fig. 5 Curves of α and random walk distance dxα

改進后的螢火蟲算法如圖6所示，其中，虛線框內內容為算法改進部分。

圖6 算法流程圖Fig. 6 Algorithm flow chart

2 算法性能測試

2.1 試驗設置

采用3個常用的基準測試函數檢驗改進FA算法的效果(如圖7所示)，算式如下：式中：d 為參數的維數，所有函數均取為2； f1為單峰函數，用于考察算法的收斂速度； f2雖是單峰函數，但其全局最優點隱藏于一條狹長的通道中不易獲得； f3是典型的非線性多模態函數，具有大量局部極值，可有效檢驗算法的全局搜索性能。f1和 f2的參數取值范圍為(?10, 10)，f3的參數取值范圍為(?32.768, 32.768)，各函數圖像如圖7所示。改進及標準FA算法的螢火蟲總數均取為60，最大迭代步數取為1500。

圖7 各測試函數3維圖Fig. 7 Graphical model of each test function

2.2 試驗結果

因螢火蟲算法是一種隨機算法，為確定統計意義上的優化效果，式(10)中每個函數分別優化20次，統計最終計算結果及計算精度達到0.01[11]所需平均迭代步數。計算結果如表1所示。限于篇幅，圖8僅給出f1與f3函數20次平均優化迭代過程函數值與迭代步數的關系。

表1 測試函數計算結果Table 1 Computed results of test functions

圖8 f1、f3函數20次平均優化迭代過程Fig. 8 The 20 average iterative processes of f1、f3

由表1知，對于所有基準測試函數，改進FA算法在優化平均值方面均有較明顯提升，且從標準差來看，改進FA算法具有更強的算法穩定性；由圖8知，在優化計算過程中，標準與改進FA算法在計算初期收斂速率相當，但迭代50次左右以后，由于參數沒有適應性功能，標準FA的優化速率急劇降低，陷入局部最優解，在多峰函數計算時，未能收斂。而改進FA算法，無論處理單峰函數還是復雜多峰函數，在900步前均能保持較好的收斂趨勢，有效提升了算法中、后期收斂速度，且計算精度遠超設定閾值，擁有較高的收斂精度。

3 基于改進FA的有限元模型修正

3.1 基于改進FA算法的有限元模型修正方法

有限元模型修正的3個主要步驟是：確定目標函數、修正參數及修正算法。目標函數往往根據可獲得的實測數據確定。在土木工程中，結構動態測試在隨機激勵下容易實現，故選取識別的固有頻率與模態振型，建立如下目標函數：

式中：n、m分別表示頻率ω與振型φ的階數；下標e、r分別表示有限元模型計算值與實測值；CV表示測量變異系數。通過變異系數加權，給予測量誤差較小的項以較大的權重，減小測量誤差對修正結果的影響。

其次，選取彈性模量、密度、幾何尺寸等作為修正參數[5]，并用向量X=[x1,x2,···,xd]表示。

最后，使用改進FA算法進行有限元模型修正。通過抽樣形成初始螢火蟲種群，每個螢火蟲代表一個參數向量X，將參數向量代入有限元模型計算結構響應，進而通過式(11)確定每個螢火蟲的適應度值，再按照圖5中的流程更新螢火蟲位置，直到達到迭代終止條件。

3.2 簡支梁模型修正

對一跨度為10 m的簡支梁進行模型修正。該梁高H=0.6 m，寬B=0.4 m，材料密度為2500 kg/m3，沿長度方向將梁均分為10個單元，記為編號1～10，并設置9個振型測點，記為編號①～⑨，如圖9所示。

圖9 簡支梁模型Fig. 9 The model of the simple-supported beam

選取10個單元的彈性模量為修正參數Ei(i=1,2, ···, 10)。為方便修正，對修正參數值進行歸一化無量綱處理，即θi=Ei/E0，其中，E0=3.3×104MPa，則修正參數為無量綱量θi。考慮環境、施工工藝等因素影響，假設混凝土實際彈性模量存在一定離散性，隨機設置各單元的彈性模量“真值”，如表2所示。

表2 各單元彈性模量真值Table 2 The actual elastic modulus of each element

在實際問題中，由于測量誤差，每次測量結果往往不盡相同，故一般使用多次測量結果以減少誤差影響。為模擬這一過程，首先將預設真值代入有限元模型計算獲取前6階頻率和前3階豎向振型；然后對各階頻率及振型分別添加1%噪聲以模擬測量誤差，最終“實測”頻率信息如表3所示。

表3 實測頻率信息Table 3 The information of measured frequency

基于上述“實測”值，分別使用原始FA及改進FA對簡支梁進行有限元模型修正。修正時，各螢火蟲參數初始值從均勻分布U(0, 2)中隨機抽取，種群中螢火蟲個數N=20，每只螢火蟲包含變量維數d=10，最大迭代步T=3000。本文僅列出某只螢火蟲迭代收斂情況，如圖10所示。參數修正結果如圖11所示。

圖10 某只螢火蟲各維參數收斂情況Fig. 10 Convergence of parameters of one firefly

圖11 修正前、后參數及誤差對比Fig. 11 Comparisons of parameters and errors before and after updating

從圖10可看出，標準FA算法由于參數不可適應性，迭代過程易陷入停滯，修正效率低，且最終計算結果未能全部收斂到真值；所提改進算法在1000步左右即可收斂，且均能收斂到真值。從圖11可知，修正前各參數誤差均大于9%，其中單元4的彈性模量誤差最大，高達66.7%；標準FA修正后模型的個別參數計算誤差不減反增，修正效果較差；改進算法修正后模型各參數誤差均降低至1.08%以下，其中單元1參數誤差僅為0.08%。故所提改進算法極大地提升了原始算法的計算性能，具有良好的修正效果。

3.3 某連續剛構橋修正

某三跨預應力混凝土連續剛構橋跨徑組合為140 m+240 m+140 m，根部梁高為15 m，跨中及邊跨直線段梁高均為4.5 m，主梁采用單箱雙室截面，箱梁腹板采用直腹板，主梁材料采用C50混凝土；橋墩主墩采用雙薄壁墩，墩身薄壁截面尺寸為16 m×25 m，兩薄壁間凈距為7 m，墩高為55 m，墩身材料采用C40混凝土。C40、C50混凝土初始彈性模量分別取為3.25×104MPa、3.45×104MPa，初始密度均取為2600 kg/m3。

利用商業有限元軟件建立該橋模型，所有構件均采用歐拉三維梁單元模擬。在橋墩底端，樁基礎與承臺基礎具有很大的剛度且與基巖嵌固良好，故將其視作固端約束；主梁與墩頂采用固接；橋梁兩端為活動支座，在有限元中各設一個單向支座，提供X、Z(垂直紙面)方向的平移約束與繞Z方向的轉動約束。初始有限元模型如圖12(a)所示，模型共包含547個節點，272個單元。

圖12 初始有限元模型、傳感器布設及參數選取示意圖Fig. 12 The initial finite element, the sensor arrangement and the selection of the parameters

傳感器布設位置如圖12(b)所示，共采用13個加速度傳感器進行實橋環境隨機激勵測試，采樣頻率設為25 Hz，分別進行5次采樣以減少隨機誤差的影響，每次采樣時間持續30 min。通過頻域分解法(frequency domain decomposition)處理測試數據獲得的前4階豎向自振頻率、豎向振型與初始模型計算結果對比如表4所示，除第4階頻率誤差較小之外，其余階頻率誤差均達到6%以上。

表4 模型修正前后計算結果對比Table 4 Comparisons of the calculation results before and after model updating

該橋主梁施工采用懸臂澆筑工藝，考慮到相鄰節段梁施工間隔較近、差異較小，為減少計算參數，將相鄰節段梁的彈性模量、密度分別視為一個參數，參數選取如圖12所示，其中1～10為選取的合并后節段編號，故全橋共選取20個修正參數。修正時，螢火蟲個數N=20，最大迭代步數T=2000。

由表4知，改進FA算法修正后有限元模型的各階頻率誤差均明顯下降，均低于3.25%，修正效果較好；而原始FA算法修正后模型的第3階、4階頻率誤差反而略微增大，最大誤差出現在第3階，達到7.34%，修正效果較差；根據表中MAC數值及圖13(限于篇幅，只列出前2階振型)可知，兩種算法修正后各階振型與實測振型吻合度均有提升，但總體上改進FA計算精度優于原始FA。

圖13 振型對比圖Fig. 13 Comparisons of mode shapes

綜上所述，改進FA算法相較于原始FA算法在實橋有限元模型修正中表現出更加良好的修正效果。

4 結論

針對標準FA算法缺陷，提出了參數自適應改進方法，通過多種單、多峰函數以及簡支梁數值算例與某大橋有限元模型修正，驗證了所提方法的有效性。所得結論如下：

(1) 根據FA算法更新原理，分析吸引度因子βij與步長因子α的變化規律，提出參數自適應方法，并通過4個具有不同形式特點的單、多峰測試函數驗證了改進FA算法相較于原始FA算法在尋優速率、精度方面的優越性。

(2) 分別將改進前后的FA算法應用于10維簡支梁模型修正，結果表明，原始算法在修正中易陷入停滯，修正精度較低，而改進算法擁有較好的修正效果。修正前后，簡支梁參數最大誤差由初始模型的66.7%降低至修正后的1.08%，且修正后最小誤差降低至0.08%。驗證了改進FA算法在結構有限元模型修正中的應用可行性。

(3) 基于環境隨機激勵測試與頻域分解法，識別某連續剛構橋的豎向前4階頻率、振型。結果表明，修正前頻率最大誤差為14.47%，改進算法修正后各階頻率相對誤差明顯下降，均控制在3.25%以內。此外，修正后的各階振型與實測振型吻合度均有提升，表明所提方法擁有良好的修正效果。