基于擴展有限元和改進海豚回聲優(yōu)化算法的損傷識別方法

師鵬杰，嚴 剛

(南京航空航天大學(xué) 航空宇航學(xué)院，機械結(jié)構(gòu)力學(xué)及控制國家重點實驗室，南京 210016)

結(jié)構(gòu)健康監(jiān)測對于保障大型和重要結(jié)構(gòu)的安全性和可靠性具有重要意義，而損傷識別是結(jié)構(gòu)健康監(jiān)測的核心環(huán)節(jié)之一。近年來在該領(lǐng)域已提出許多不同類型的損傷識別方法，其中基于模型的方法由于不僅能夠識別損傷的位置，還能識別損傷的尺寸或程度，受到了廣泛的關(guān)注和研究[1]。

通常基于模型的結(jié)構(gòu)損傷識別被歸為一類反問題來求解，需要涉及兩方面的內(nèi)容：含損傷結(jié)構(gòu)的建模和反問題的求解。在含損傷結(jié)構(gòu)的建模方面，有限元方法以其對各種形式的結(jié)構(gòu)及損傷所具有的廣泛適應(yīng)能力，在結(jié)構(gòu)損傷識別中成為最常用的一種方法[2-4]。但有限元方法在用于連續(xù)體結(jié)構(gòu)損傷識別時需要在損傷參數(shù)更新的過程中重新劃分網(wǎng)格，以適應(yīng)損傷邊界的變化，大大影響了計算效率。另一種常用的建模方法是邊界元方法，將彈性控制方程轉(zhuǎn)換為邊界積分方程，并采用格林函數(shù)進行求解[5-6]。邊界元的缺點是對于非均質(zhì)材料和任意形狀的損傷難以獲得格林函數(shù)，從而限制了它的應(yīng)用。擴展有限元是近年來在計算斷裂力學(xué)領(lǐng)域快速發(fā)展的一種數(shù)值方法，其最主要的特點是，所使用的網(wǎng)格與結(jié)構(gòu)內(nèi)部的幾何或物理界面無關(guān)，克服了在諸如裂紋尖端等高應(yīng)力和變形集中區(qū)進行高密度網(wǎng)格剖分的困難[7-10]。正是這種優(yōu)點使得擴展有限元在應(yīng)用于損傷識別等反問題時也具有很大優(yōu)勢，如RABINOVICH et al采用時間簡諧響應(yīng)和擴展有限元對薄膜中的裂紋進行識別[11]；WAISMAN et al[12]，CHATZI et al[13]進一步發(fā)展了基于擴展有限元的損傷識別方法，采用靜態(tài)響應(yīng)對平面問題的多種損傷以及多個損傷進行識別，表明了擴展有限元在結(jié)構(gòu)損傷識別中的有效性。

在反問題求解方面，各種優(yōu)化算法得到了廣泛的應(yīng)用。如早期較多采用的基于梯度的優(yōu)化算法[2]，但是由于其只能收斂到局部最優(yōu)解且求解的結(jié)果強烈依賴于初始值的選擇而限制了發(fā)展；目前廣泛應(yīng)用的各種智能算法，如遺傳算法、人工蜂群算法等[12-16]具有操作簡單，收斂速度快，全局收斂能力好等優(yōu)點而得到更多的的關(guān)注。其中KAVEH et al根據(jù)海豚運用回聲捕食獵物的行為設(shè)計了一種元啟發(fā)式算法——海豚回聲算法(dolphin echolocation optimization alogrithm，DEO)，并將其成功運用于求解剛架和桁架結(jié)構(gòu)問題[17]，算例表明DEO算法與遺傳算法、粒子群算法、蟻群算法等相比具有更好的收斂性能。

本文提出了一種搜索范圍動態(tài)變化的尋優(yōu)機制來改進海豚回聲算法，將之與擴展有限元法相結(jié)合應(yīng)用于損傷識別問題，并運用數(shù)值仿真方法來驗證所提出的改進算法以及損傷識別方法的合理性和有效性。

1 損傷識別方法

1.1 損傷識別反問題

本文將損傷識別歸為反問題來求解。如圖1所示，用一定的損傷參數(shù)來表征損傷的位置和大小，在獲得實際結(jié)構(gòu)測量響應(yīng)后，通過優(yōu)化算法最小化實際響應(yīng)與損傷結(jié)構(gòu)模型模擬響應(yīng)之差來識別優(yōu)化參數(shù)。

圖1 含損傷結(jié)構(gòu)示意圖Fig.1 A generic solid structure with void flaws

本文采用圓形孔洞來模擬損傷，并以應(yīng)變作為結(jié)構(gòu)測量響應(yīng)。對于圓形孔洞而言，其損傷參數(shù)為中心坐標(biāo)和半徑(xi,yi,ri)，i=1,2,…,n，其中n為孔洞個數(shù)。建立目標(biāo)函數(shù)：

(1)

式中：λ為由損傷參數(shù)(xi,yi,ri)組成的向量；ε(m)和ε(s)分別為在實際結(jié)構(gòu)中測量出的應(yīng)變響應(yīng)和損傷結(jié)構(gòu)模型計算出的應(yīng)變響應(yīng)，應(yīng)變向量由兩個方向的正應(yīng)變組成；D為應(yīng)變傳感器數(shù)量。

為了適應(yīng)下文改進海豚回聲算法以最大值為優(yōu)化目標(biāo)的需要，將目標(biāo)函數(shù)調(diào)整為：

(2)

其中：δ為篩選系數(shù)，0<δ<1，本例δ取0.85；當(dāng)g1(λ)<0時，隨機生成一個λ0替換λ，并且要求g1(λ0)≥0.

1.2 擴展有限元建模

本文采用擴展有限元法對結(jié)構(gòu)進行建模，計算出含有不同位置和大小的損傷結(jié)構(gòu)的響應(yīng)量。擴展有限元通過引入非連續(xù)的位移模式來代替?zhèn)鹘y(tǒng)有限元的連續(xù)位移模式，使得對損傷結(jié)構(gòu)不連續(xù)位移場的描述得到了簡化，在不用重新劃分網(wǎng)格的前提下可以方便地改變損傷的位置和大小。

擴展有限元中所應(yīng)用的位移模式為：

(3)

式中：NI(x)為常規(guī)有限元中的形函數(shù)；uI為節(jié)點處的未知位移；ψI(x)為擴展有限元的加強函數(shù)；aI為改進后增加的自由度。

可以通過水平集函數(shù)來表示損傷的位置和大小。對于本文中考慮的孔洞，其水平集函數(shù)為：

(4)

式中：nc為圓形孔洞數(shù)；xi和ri為第i個孔洞的中心和半徑。

針對孔洞損傷而言，其位移模式可以采用一種更簡單的表述[18]：

(5)

(6)

2 海豚回聲算法及其改進

2.1 海豚回聲算法

DEO算法是一種根據(jù)海豚回聲搜索獵物的機制所設(shè)計的一種離散型優(yōu)化算法。算法定義了收斂因數(shù)(Fc)用來表示每一次迭代的最優(yōu)解為全局最優(yōu)解的概率，一般根據(jù)一條預(yù)定的曲線來變化。本文引用文獻[16]中的收斂曲線來定義：

(7)

式中：Fc1為預(yù)定義的第一個循環(huán)最優(yōu)解為全局最優(yōu)解的概率；pi表示第i個循環(huán)；N為總循環(huán)次數(shù)；p為收斂曲線的次數(shù)，可以通過調(diào)整收斂曲線的次數(shù)來調(diào)整收斂曲線的形狀。

算法需要先將連續(xù)變量的可行域離散化然后以升序或者降序進行排列，用長度為LAj的列向量Aj表示離散化后變量lj(j=1,2,…,NV)的可取值，其中NV為變量的個數(shù)。

DEO算法的求解流程為：

1) 隨機生成NL個可行解。初始化LNL×NV=[L1,L2,…,LNL]矩陣，其中NV是變量的個數(shù)，NL為海豚所搜尋位置的數(shù)量。

2) 計算當(dāng)前循環(huán)的Fc.

3) 計算每一個可行解的目標(biāo)函數(shù)值g1(Li)，其中目標(biāo)函數(shù)以最大值為優(yōu)化目標(biāo)。

4) 計算每個變量的累計適應(yīng)度：

(8)

式中：Aij為Li(j)在Aj中的位置，F(xiàn)(Aij+k,j)是第j個變量的第(Aij+k)個可取值的累計適應(yīng)度；R是影響半徑，它表示了位置Li(j)對周圍的影響范圍，本文取LAj的1/8.

為了使分布更加均勻，將一個小量的δ加到F上，即F=F+δ，δ最好選擇比目標(biāo)函數(shù)值小的值。

5) 對于每一個變量lj，計算每一個可取值A(chǔ)j(i)(i=1,2,…,LAj)的概率值。首先找到本次循環(huán)最優(yōu)解的各變量對應(yīng)在F中的值并將其置為0，依據(jù)輪盤賭的方式計算可取值A(chǔ)j(i)的概率：

(9)

將概率Fc分配給最優(yōu)解所對應(yīng)各個變量的可行值，(1-Fc)分配給其他可行值。

6) 根據(jù)計算出的概率值生成下一循環(huán)的NL個可行解。

7)回到步驟3)，直到滿足結(jié)束條件，輸出結(jié)果。

2.2 改進海豚回聲算法(EDEO)

當(dāng)所要求解的問題變量取值范圍大，要求精度高時，則需要將可行域離散的更加細致，這樣會帶來計算規(guī)模加大，計算時間增長，后期收斂變慢，甚至找不到全局最優(yōu)解等問題。針對上述問題，首先考慮到當(dāng)DEO算法收斂到最優(yōu)解附近，可行解的目標(biāo)函數(shù)值已經(jīng)相差不大，所以將累計適應(yīng)度計算方法改成：

(10)

式中：Fmin為NL個可行解中目標(biāo)函數(shù)值的最小值，其他同式(8).其次本文提出了一種搜索范圍動態(tài)變化的尋優(yōu)機制來改進DEO算法。

這種尋優(yōu)機制主要包含初搜索、精搜索、平移搜索、回溯搜索和重搜索等5個過程，并且引入兩個概念：

1) 搜索域。本次循環(huán)所搜索的范圍，對于變量Dj用[mj,Mj]表示其搜索域。

2) 縮放搜索域。若Dj的搜索域為[m0j,M0j]，D0j為任意一個在搜索域內(nèi)的實數(shù)，則定義以D0j為中心的關(guān)于[m0j,M0j]的縮放可行域[mj,Mj]為：

(11)

式中：k為縮放系數(shù)，當(dāng)縮小時k>1，放大時k<1，不變時k=1.

5個搜索過程具體如下。

3) 平移搜索。當(dāng)k1取值過大，導(dǎo)致全局最優(yōu)解不在精搜索的搜索域之內(nèi)，則精搜索求得的最優(yōu)解會在搜索域邊緣，即當(dāng)?j=1,2,…,NV使得：

(12)

圖2 EDEO算法示意圖Fig.2 Diagram of EDEO Algorithm

5) 重搜索。當(dāng)調(diào)用精搜索過程過多，搜索域會變得很小，對于一般優(yōu)化問題已經(jīng)沒有了意義，為了提高算法全局收斂能力，對?j=1,2,…,NV滿足：

(13)

式中：η為容許誤差，根據(jù)求解問題精度來確定。則以可行域為搜索域，以p1為循環(huán)次數(shù)，更新Lr.

當(dāng)循環(huán)次數(shù)達到最大循環(huán)次數(shù)p時，循環(huán)結(jié)束，輸出最優(yōu)解。程序框圖如圖3.

圖3 EDEO算法流程圖Fig.3 Flowchart of EDEO Algorithm

3 算例分析

本文采用數(shù)值仿真研究驗證所提出的改進算法以及損傷識別方法的合理性和有效性。算例所采用的方形平板處于平面應(yīng)力狀態(tài)，底部受固定端約束，頂部受100 MPa的均載拉伸載荷作用，板的邊長為10 cm；平板被劃分成40 mm×40 mm的均勻網(wǎng)格；在板的左、上、右邊緣均勻分布著14個應(yīng)變傳感器，如圖4所示；材料性質(zhì)為E=70 GPa，υ=0.3；實際結(jié)構(gòu)測量響應(yīng)通過網(wǎng)格密度為50 mm×50 mm的擴展有限元模擬模型得到，考慮到噪聲的影響，向?qū)嶋H結(jié)構(gòu)響應(yīng)加不同大小的高斯白噪聲(5%，15%).待求解參數(shù)可行域為xi∈[0.5,9.5],yi∈[0.5,9.5],ri∈[0.2,2.0].

圖4 算例平板示意圖Fig.4 Diagram of plate

本算例中EDEO算法第一個循環(huán)的收斂因數(shù)Fc1=0.1，收斂曲線次數(shù)p=1.2，LAj=21，小量δ=0.05，容許誤差η=1×10-2，精搜索縮放系數(shù)k1=4，回溯搜索縮放系數(shù)k2=1/10，容許平移系數(shù)q0=2.DEO算法中其他參數(shù)和EDEO算法相同，LAj取1 801.

3.1 單個孔洞損傷識別

針對單個孔洞問題，EDEO算法可行解個數(shù)NL=20，p1=10，最大循環(huán)次數(shù)p=100，待求解的優(yōu)化參數(shù)為x=4.5,y=5.5，r=0.6.

圖5給出了15%噪聲下EDEO參數(shù)估計過程。由圖可以看出，EDEO算法經(jīng)過25次迭代，圓心坐標(biāo)為(4.47，5.62)，半徑為0.605，識別結(jié)果已經(jīng)比較精確。圖6給出了在15%噪聲條件下EDEO算法的前25個循環(huán)的識別結(jié)果和DEO算法的結(jié)果，可以看出EDEO算法在25個循環(huán)就已經(jīng)收斂到了不錯的結(jié)果，而DEO算法迭代了100循環(huán)之后仍然距離真實解較遠。表1給出了不同噪聲、不同算法的優(yōu)化結(jié)果?？梢钥闯鯡DEO算法所求得的結(jié)果非常精確，而DEO算法只能收斂到所求參數(shù)的附近而無法求得更加精確的解。

圖5 15%噪聲下EDEO參數(shù)估計過程Fig.5 Process of parameter estimation with noise level of 15%

圖6 15%噪聲下EDEO部分循環(huán)識別結(jié)果和DEO最終結(jié)果Fig.6 Identified process of EDEO and DEO with noise level of 15%

表1 單個孔洞缺陷求解結(jié)果Table 1 Identified results of a single void

3.2 兩個孔洞損傷識別

針對兩個孔洞問題，EDEO算法可行解個數(shù)NL=30，p1=30，最大循環(huán)次數(shù)p=200，待求解的優(yōu)化參數(shù)為x1=8，y1=8，r1=0.6，x2=2.5，y2=2.5，r2=0.6.

圖7為15%噪聲下EDEO參數(shù)估計過程，由圖可以看出在經(jīng)過40次迭代之后，孔洞1圓心坐標(biāo)(7.81，7.81)，半徑0.67，孔洞2圓心坐標(biāo)(2.64,2.19)，半徑為0.72，已經(jīng)能夠得到較為不錯的結(jié)果。圖8給出了在15%噪聲條件下EDEO算法的部分循環(huán)的識別結(jié)果和DEO算法的結(jié)果，由圖可以看出當(dāng)EDEO算法迭代40次之后所得到的結(jié)果比DEO算法最終結(jié)果還要精確，DEO算法已經(jīng)陷入局部最優(yōu)無法收斂。表2給出了不同噪聲，不同算法的優(yōu)化結(jié)果，可以看出，在不同噪聲下EDEO算法比DEO算法求解得到的結(jié)果都要精確很多。

圖7 15%噪聲下EDEO參數(shù)估計過程Fig.7 Process of parameter estimation with noise level of 15%

表2 兩個孔洞缺陷求解結(jié)果Table 2 Identified results of two voids

3.3 三個孔洞損傷識別

針對三個孔洞EDEO算法可行解個數(shù)NL=50，p1=50，最大循環(huán)次數(shù)p=200，待求解的優(yōu)化參數(shù)為x1=2.5，y1=2.5，r1=0.6，x2=3.0，y2=8.0，r2=0.6，x3=8.0，y3=8.0，r3=0.6.

圖9給出了在15%噪聲下EDEO算法部分循環(huán)識別結(jié)果，由圖可以看出EDEO算法經(jīng)過150次迭代，孔洞1的圓心坐標(biāo)為(2.185,2.461)，半徑為0.575；孔洞2的圓心坐標(biāo)為(3.002,7.920)，半徑為0.592；孔洞3的圓心坐標(biāo)為(7.898,7.989)，半徑為0.584，已經(jīng)收斂到了不錯的結(jié)果。圖10給出了在15%噪聲下DEO算法的收斂結(jié)果，可以看出DEO算法已經(jīng)無法收斂。表3給出了15%噪聲條件下EDEO算法的優(yōu)化結(jié)果，結(jié)果表明對于3個孔洞EDEO算法能夠收斂到真實值。

圖8 15%噪聲下EDEO部分循環(huán)識別結(jié)果和DEO最終結(jié)果Fig.8 Identified process of EDEO and DEO with noise level of 15%

圖9 15%噪聲下EDEO部分循環(huán)識別結(jié)果Fig.9 Identified process of EDEO with noise level of 15%

圖10 15%噪聲下DEO最終識別結(jié)果Fig.10 Identified result of DEO with noise level of 15%

表3 3個孔洞EDEO求解結(jié)果Table 3 EDEO identified results of three voids

4 結(jié)論

本文結(jié)合擴展有限元與改進海豚回聲算法建立了損傷識別模型，并重點介紹了海豚回聲算法的建立及改進過程。海豚回聲算法運用收斂因數(shù)這一參數(shù)控制迭代過程中最優(yōu)解為全局最優(yōu)解的概率。為了提高收斂性能，提出了一種搜索范圍動態(tài)變化的尋優(yōu)機制，根據(jù)其搜索域大小和范圍的不同分為初搜索、精搜索、平移搜索、回溯搜索、重搜索5個過程。通過數(shù)值仿真研究結(jié)果表明，改進海豚回聲算法與原算法相比有較高的效率和較好的收斂性能，算法在存在模型誤差和隨機測量誤差的情況下仍然能夠得到較好的識別結(jié)果，驗證了改進算法的有效性和損傷識別方法的合理性。

本文未討論在未知損傷數(shù)量條件下如何識別損傷參數(shù)。尋求一種高效的判斷損傷數(shù)量并將其引入優(yōu)化過程的方法是本文之后研究的重點。