封周權(quán)，王文贊，華旭剛，陳政清

［1.湖南大學土木工程學院，湖南長沙 410082；2.風工程與橋梁工程湖南省重點實驗室（湖南大學），湖南長沙 410082］

在土木工程領(lǐng)域中，基于有限元技術(shù)的數(shù)值模擬方法廣泛應用于結(jié)構(gòu)靜力分析、結(jié)構(gòu)動力分析、結(jié)構(gòu)或構(gòu)件設(shè)計等各個方面［1］，通過建立有限元模型可為結(jié)構(gòu)行為預測提供重要參考.有限元模型修正實質(zhì)上是一種系統(tǒng)識別問題，通過修正模型參數(shù)來盡可能縮小模型響應數(shù)據(jù)與實測響應數(shù)據(jù)之間的差別，從而使得有限元模型更加接近實際工程結(jié)構(gòu).模態(tài)參數(shù)（頻率、振型等）是表征結(jié)構(gòu)動力特性的主要參數(shù)，因此基于模態(tài)參數(shù)的模型修正得到了廣泛的應用.求解使得模態(tài)參數(shù)模型預測值與實測值之間誤差取得最小值時的模型參數(shù)，即得到了模型參數(shù)的最優(yōu)值，因此模型修正在數(shù)學上實際是一個優(yōu)化問題.針對優(yōu)化問題的求解，現(xiàn)已有多種算法得到了應用，如模式搜索法［2］、擬牛頓法［3］、元啟發(fā)式優(yōu)化算法等.其中元啟發(fā)式優(yōu)化算法由于其優(yōu)越的全局尋優(yōu)能力和較好的準確性而受到了學者們的青睞.王家等［4］提出新的啟發(fā)式算法來求解施工現(xiàn)場設(shè)施布局的優(yōu)化問題；Feng 等［5］針對類電磁機制算法提出改進措施，并成功應用于三層剪切框架的結(jié)構(gòu)損傷識別.

螢火蟲算法（FA）作為元啟發(fā)式優(yōu)化算法的一種，具有控制參數(shù)少、物理意義簡單明了的特點，具有良好的尋優(yōu)性能［6］，但同其他元啟發(fā)式優(yōu)化算法，如遺傳算法（GA）、粒子群算法（PSO）等一樣，存在容易陷入局部最優(yōu)、早熟、不收斂的問題，因此學者針對原始螢火蟲算法的不足提出了諸多改進措施［7-12］.

本文基于模態(tài)參數(shù)識別技術(shù)，使用作者提出的改進螢火蟲算法對二維桁架數(shù)值模型和六層剪切框架實驗模型進行模型修正，并將其識別結(jié)果與原FA、GA、PSO 等算法識別結(jié)果進行對比，驗證了改進優(yōu)化算法的優(yōu)越性.

1 目標函數(shù)

對于線彈性結(jié)構(gòu)模型，頻率和振型的預測值與其質(zhì)量和剛度矩陣直接相關(guān).假設(shè)結(jié)構(gòu)的質(zhì)量矩陣M已知，結(jié)構(gòu)模型修正僅考慮剛度修正，且引入剛度修正系數(shù)θ來描述剛度變化程度.修正系數(shù)θ∈[-1，1]為一個1×N維向量，其中N為未知量個數(shù)，應為結(jié)構(gòu)的單元個數(shù).若θi>0，則表示第i個單元剛度提高，若θi<0，則表示該單元剛度降低.另假設(shè)結(jié)構(gòu)在未受損傷狀態(tài)下第i個單元的剛度矩陣為ki，那么修正狀態(tài)（損傷狀態(tài)）的剛度矩陣K(θ)可表示為：

結(jié)構(gòu)的頻率和振型預測值可以由其特征方程求解得到：

式中：ω為結(jié)構(gòu)圓頻率，與結(jié)構(gòu)物理頻率f有ω=2πf的關(guān)系；?為與之對應的特征向量，即結(jié)構(gòu)的振型.

在結(jié)構(gòu)動力實測中，受制于人力、成本、器材等諸多因素，高階模態(tài)參數(shù)往往難以精確獲取，全自由度測量也經(jīng)常難以實現(xiàn)，因此假設(shè)實測僅可獲得前m階模態(tài)和nd個自由度，并重復測量ns次，即可獲得ns組m階頻率和不完整振型實測值.將模態(tài)參數(shù)實測值與預測值之間的誤差分別按頻率和振型以及模型階數(shù)進行累積，設(shè)計目標函數(shù)［5］：

式中：aj，i=為比例系數(shù)，確保預測模態(tài)振型在測量自由度處與實測模態(tài)振型最接近；分別為結(jié)構(gòu)第j階模態(tài)頻率和振型向量預測值，它們是將式（1）代入式（2）并通過模態(tài)分析計算得到的；fj，i和?j，i分別為第j階模態(tài)頻率和振型向量的第i次實測值，它們是通過實驗模態(tài)參數(shù)識別得到的.為避免振型和頻率誤差水平的不同，公式（3）以實測值（振型向量取實測值的模）為分母進行了歸一化處理.

2 原始算法與優(yōu)化算法

2.1 螢火蟲算法

螢火蟲算法最早由Yang［13］提出，是一種以自然界中螢火蟲發(fā)光行為［14］為啟發(fā)的全局優(yōu)化算法.螢火蟲算法通過建立不同螢火蟲個體之間相互吸引機制，從而達到在空間內(nèi)尋找最優(yōu)解的目的，即若某螢火蟲熒光亮度（I）越亮，則對其他個體的吸引力（β）也就越強.

螢火蟲算法做了以下3點假設(shè)：

1）螢火蟲沒有性別上的區(qū)分，即所有螢火蟲之間的相互吸引不受性別的影響.

2）螢火蟲的吸引力（β）與熒光亮度（I）成正比，且兩者都隨距離的增加而降低.

3）螢火蟲的亮度（I）取決于該螢火蟲個體所處目標函數(shù)中的位置.

在一定的導光介質(zhì)中，光線強度以高斯分布的形式隨距離的增加而衰減：

式中：r為兩相鄰螢火蟲之間的空間距離；I0為距離r=0時的熒光亮度；γ為導光介質(zhì)的吸收系數(shù)，其值的大小代表了介質(zhì)對光亮衰減的影響程度.

螢火蟲的吸引系數(shù)同樣以指數(shù)的形式隨距離的增加而衰減：

式中：β0為吸引系數(shù)初值，代表兩螢火蟲之間距離r=0時的吸引系數(shù).

兩螢火蟲之間的距離滿足笛卡兒坐標系下的空間距離計算方法：

式中：n為優(yōu)化問題的自變量總數(shù)，即解的總維數(shù)；xi，m是第i只螢火蟲第m維變量.

如前所述，熒光亮度（I）較亮的螢火蟲可吸引較暗的螢火蟲向其移動，假設(shè)第i只螢火蟲熒光亮度較第j只螢火蟲弱，則螢火蟲i的移動方式和移動量受式（7）約束：

2.2 改進螢火蟲算法

改進螢火蟲算法（modified Nelder-Mead Firefly Algorithm，m-NMFA）在原標準FA 算法的基礎(chǔ)上提出4 點改進措施，包括引入下山單純形局部優(yōu)化算法（Nelder-Mead Algorithm），并針對算法的結(jié)構(gòu)、參數(shù)以及局部搜索提出了改進措施，如引入多樣性閾值、引入邊界控制因子、對參數(shù)進行修改等.

2.2.1 參數(shù)修改

大量的研究證明，改變步長α和吸引系數(shù)β對螢火蟲算法的求解精度、收斂速度等都有較大的影響［6-7，11，20-23］.

經(jīng)過測試發(fā)現(xiàn)，在迭代過程的前期，步長偏小的往往會收斂更快，這個規(guī)律適合不同的步長公式和測試函數(shù)，但如果步長減小過快，在后期容易陷入局部最優(yōu)解而難以尋找到全局最優(yōu)解，使算法早熟.因此，步長在迭代后期逐漸趨于0 的過程中仍應保持足夠大的值，以最大限度地保持探索與開發(fā)之間的平衡［24-25］.據(jù)此設(shè)計出較為理想的步長公式：

式中：α(t)表示第t迭代步的隨機移動的步長；C=1.0×10-10；Gen，max表示最大迭代次數(shù)；λ∈[1，3]為參考系數(shù)，其值隨著迭代次數(shù)t的增加呈非線性遞減趨勢.引入?yún)⒖枷禂?shù)λ可實現(xiàn)對步長大小變化趨勢的控制.

為使步長大小與可行域大小相適配［25］，在式（7）α?i項中引入縮放參數(shù)（scaling parameters）Sk.

式中：lk和uk分別為某螢火蟲個體的第k個未知量所對應的下界和上界.如此，式（7）可寫為：

對于吸引系數(shù)，引用Selvarasu［23］在其改進方法中給出的吸引系數(shù)計算公式，以保證算法的精確度和收斂速度：

式中：βmin=0.2 為吸引系數(shù)β的下界，βmax=1 為吸引系數(shù)β的上界，因此有β∈[βmin，βmax].

2.2.2 邊界約束處理機制

螢火蟲算法與其他元啟發(fā)式優(yōu)化算法類似，個體有可能出現(xiàn)在所定義的邊界之外［26］，這些個體將對算法的精度和收斂速度均產(chǎn)生不利影響，甚至得到完全錯誤的結(jié)果.使用有效的邊界約束處理機制可以很好地避免此類問題的發(fā)生，并提高搜索算法的優(yōu)化性能［27-28］.

在已有的邊界約束方式中，較多的學者［15，22，27-28］針對螢火蟲個體中溢出的某個未知量而非單個個體進行處理，但多維度問題的最優(yōu)解是由各個未知量共同決定的，因此提出新的邊界約束處理機制，即當某螢火蟲個體的某未知量超出邊界值時，首先在可行域范圍內(nèi)隨機生成新的螢火蟲個體，然后將該螢火蟲個體向當前最優(yōu)解移動一定距離，在移動時能夠保證移動幅度不超過該螢火蟲與當前最優(yōu)解之間距離的1/2，以最大限度地保證個體多樣性.其偽代碼如下：

其中，rand 為0～1 的隨機數(shù)，F(xiàn)∈[0，0.5]為隨機因子，xbest為當前全局最優(yōu)解.

2.2.3 下山單純形算法及多樣性閾值

下山單純形算法是一種局部優(yōu)化算法，與群體智能優(yōu)化算法通過一定規(guī)模個體隨機尋求最優(yōu)解的方法不同，它在給定的初始解附近生成N+1 個頂點的凸多面體，在此基礎(chǔ)上通過反射、擴張、收縮和壓縮等運算逐漸逼近精確解［29］.因此在螢火蟲算法迭代過程中及運行結(jié)束后，使用下山單純形算法在最優(yōu)解附近挖掘更優(yōu)解，以提高算法求解能力.

在迭代過程中，啟用單純形算法的時機將直接影響求解的精度和收斂速度，若啟用過早，則接近真實解的最優(yōu)解還未找到，單純形算法頻繁調(diào)用，將降低求解效率；若啟用過晚，通過單純形算法挖掘到的更優(yōu)解對算法整體的貢獻將會被削弱，因此引入多樣性閾值概念以解決上述矛盾.記xk，best和xk，worst分別為第k次迭代最優(yōu)解和最差解，xg，best和xg，worst分別為當前全局最優(yōu)解和最差解.多樣性閾值計算公式如下：

規(guī)定當eξ-1

經(jīng)過測試，取M=0.001 是相對合適的，此時算法已經(jīng)能夠找到接近于真實解的全局最優(yōu)解，而且經(jīng)過若干次的迭代，步長α已經(jīng)足夠小，能夠與單純形算法相適配.

3 數(shù)值方法測試

3.1 基準函數(shù)測試

為了驗證改進優(yōu)化的求解效果，選用表1 所示的4 個基準測試函數(shù)將FA、GA 和PSO 優(yōu)化算法與m-NMFA進行對比，以驗證改進算法的優(yōu)越性.測試中，種群規(guī)?；蛉旧w個數(shù)n均為30，最大迭代次數(shù)均為Gen，max=1 000.其中，m-NMFA 的步長初值為0.5，原始FA 采用Yang［15］對于多維度問題所采用的步長公式，步長初值為0.25.吸引系數(shù)βmin和βmax分別取0.2和1.0，光線吸收常數(shù)γ=1.0.對于GA，測試調(diào)用了Matlab 內(nèi)置GA 優(yōu)化算法，除染色體個數(shù)及最大迭代次數(shù)外，其他參數(shù)均采用Matlab 內(nèi)置工具箱默認值.對于PSO，學習因子取c1=c2=2，粒子最大速度為1，慣性因子采用線性變化方程：

表1 基準函數(shù)Tab.1 Benchmark function

式中：ωi為第i次迭代的慣性因子，ω1取0.4，ω2取0.9.當某粒子某一維度的速度或位置達到或超出邊界時，該粒子全部維度的速度或位置將會全被重置為邊界值.

為對比各優(yōu)化算法的穩(wěn)定性，每個測試函數(shù)單獨運算100 次，得到優(yōu)化結(jié)果的平均值、標準差及最小值如表2所示.

從表2可以看出，m-NMFA 除了在函數(shù)4中的平均值和標準差略差于GA 外，其余計算結(jié)果均優(yōu)于其他優(yōu)化算法，尤其是在函數(shù)2 和函數(shù)3 中的平均值、函數(shù)2～4 中最小值、函數(shù)2 中的標準差精度較其他算法均大幅提高.由此看出，m-NMFA 可以使求解的精度大大提高，并且有足夠的穩(wěn)定性.圖1為4 種優(yōu)化算法收斂情況的對比結(jié)果，從圖1 可以看出，m-NMFA 算法在收斂速度和運算精度上均明顯優(yōu)于其他3 種算法，并且存在第二、第三次深度求解的能力.

表2 優(yōu)化算法解的平均值、標準差和最小值Tab.2 The mean，standard deviation and minimum of the solutions of the four optimization algorithms

圖1 4種優(yōu)化算法的收斂圖Fig.1 Convergence diagrams of four optimization algorithms

3.2 桁架模型數(shù)值模擬

數(shù)值模擬使用如圖2 所示的兩端簡支桁架模型進行測試.圖2 中桁架長10 m，材料彈性模量為6.88×1010Pa，密度為2 780 kg/m3，桿的橫截面積為0.002 5 m2.數(shù)值模擬初始工況有如下假設(shè)：①測試值僅為前6階模態(tài)以及圖中箭頭所指的9個自由度；②受環(huán)境影響，頻率和振型的測試結(jié)果中均存在1%的白噪聲；③結(jié)構(gòu)的損傷僅考慮剛度損傷，第4、5、6 號單元分別受到了20%、30%和20%的剛度折減.

圖2 簡支桁架結(jié)構(gòu)模型（單位：m）Fig.2 Schematic of the 17-bar planar truss（unit：m）

本次模擬同樣采用FA、GA、PSO、m-NMFA 這4種優(yōu)化算法做對比，除自變量個數(shù)做了調(diào)整外，其余參數(shù)均與基準函數(shù)測試保持一致.每種算法各測試100 次，得到各單元損傷系數(shù)θn的平均值和標準差，結(jié)果如圖3 所示.此外，模擬試驗采用兩種對比工況與原工況作對比，以驗證改進優(yōu)化算法計算模型修正問題結(jié)果的準確性與穩(wěn)定性.對比工況Ⅰ：在初始工況假設(shè)的基礎(chǔ)上，將可測得的模態(tài)數(shù)從6 階減少為4 階；對比工況Ⅱ：在初始工況假設(shè)的基礎(chǔ)上，將白噪聲由1%增加為5%.

圖3 原工況下不同算法計算得到的桁架模型損傷系數(shù)θn的平均值及標準差Fig.3 The mean and standard deviation of damage coefficient of the truss model of the original condition run by different algorithms

從圖3 可以看出，F(xiàn)A、GA、PSO 三種優(yōu)化算法的最優(yōu)解與目標值存在較大偏差，離散程度較大.通過改進算法得到的受損處單元損傷系數(shù)平均值分別為20.05%、29.86%和20.03%，非受損處單元損傷系數(shù)誤差均小于5%，與假設(shè)情況符合良好；最優(yōu)解的標準差均小于10-3，離散程度低.因此可以認為，使用m-NMFA算法尋找模態(tài)參數(shù)識別問題的最優(yōu)解是可靠的.

從圖4 可以看出，對于工況Ⅰ，模態(tài)階數(shù)減少至4階后，改進算法仍能在損傷部位找到較準確的損傷系數(shù)，且第5號單元損傷系數(shù)較原工況的損傷系數(shù)偏差僅有0.22%；工況Ⅱ中受損處單元損傷系數(shù)相較于原工況分別存在3.33%、1.83%和3.57%的增加，對于土木工程結(jié)構(gòu)而言，這將是偏保守和偏安全的結(jié)果.

圖4 不同工況下m-NMFA計算得到的桁架模型損傷系數(shù)θn的平均值及標準差Fig.4 The mean and standard deviation of damage coefficient of the truss model of different working conditions run by m-NMFA

4 剪切框架試驗

選用6 層鋼結(jié)構(gòu)剪切框架，通過測試模態(tài)參數(shù)并采用本文提出的改進螢火蟲算法對模型剛度矩陣進行修正從而實現(xiàn)損傷識別.鋼結(jié)構(gòu)框架如圖5（a）所示，結(jié)構(gòu)由地腳螺栓固定在地面上，層間板長、寬、厚分別為450 mm、450 mm、10 mm，各板在厚度中點處距離均為250 mm，各個層間板之間的剛度由分布在層間板4 個頂點處的主要柱和側(cè)邊中點處的附加柱提供.剪切框架的有限元模型采用如圖5（b）所示的6 自由度集中質(zhì)量模型.各層間板中心位置處設(shè)置一塊豎直擋板，并使用采樣頻率為200 Hz 的IL-300 激光位移計來測量擋板底部的位移.試驗測試了5 層附加柱缺失（測試Ⅰ）、第三和第五層附加柱缺失（測試Ⅱ）兩種工況，為保證測試結(jié)果的準確性與穩(wěn)定性，每種工況均重復測量50 次，每次持續(xù)時長5 min.剪切框架前6 階固有頻率和振型如表3 和圖6 所示，其中表3 所列頻率為50 次測量的平均值，振型均為位移歸一化后的結(jié)果.

圖5 鋼結(jié)構(gòu)剪切框架及其集中質(zhì)量有限元模型Fig.5 Steel structure shear frame and its lumped mass finite element model

圖6 剪切框架試驗前6階振型圖Fig.6 The first 6 mode shapes of the shear frame tests

表3 剪切框架試驗前6階頻率Tab.3 The first 6 natural frequencies of the shear frame tests Hz

試驗測得的結(jié)果如圖7和圖8所示.在測試Ⅰ中，第五層剛度減少21.36%；在測試Ⅱ中，第三層和第五層結(jié)構(gòu)剛度分別減少21.76%和22.81%.兩組測試中，未損傷層的結(jié)構(gòu)剛度平均值誤差均在5%以內(nèi)，因此可以認為，改進的優(yōu)化算法能夠準確識別結(jié)構(gòu)損傷的位置，并且能夠得到符合實際期望的最優(yōu)解.

圖7 測試Ⅰ結(jié)構(gòu)損傷系數(shù)平均值及標準差Fig.7 Mean and standard deviation of structural damage coefficient of test Ⅰ

圖8 測試Ⅱ結(jié)構(gòu)損傷系數(shù)平均值及標準差Fig.8 Mean and standard deviation of structural damage coefficient of test Ⅱ

5 結(jié)論

基于結(jié)構(gòu)的模態(tài)參數(shù)數(shù)據(jù)，通過m-NMFA 改進螢火蟲算法尋找模型參數(shù)最優(yōu)解，成功地識別出受損結(jié)構(gòu)的剛度變化，主要結(jié)論如下：

1）改進螢火蟲算法可在保持較強全局搜索能力的同時，引入4點改進措施，極大地提高了局部搜索能力，其最優(yōu)解相比于其他優(yōu)化算法更精確、更穩(wěn)定.

2）基于結(jié)構(gòu)的頻率、振型等模態(tài)參數(shù)，采用改進算法可利用較少的模態(tài)階數(shù)和自由度，在一定白噪聲影響下準確地找到結(jié)構(gòu)損傷的位置，并求得有較高可信度的剛度損傷，為有限元模型修正準確度提供了有力保證，具有很好的實用價值.