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

師鵬杰，嚴 剛

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

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

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

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

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

1 損傷識別方法

1.1 損傷識別反問題

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

圖1 含損傷結構示意圖Fig.1 A generic solid structure with void flaws

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

(1)

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

為了適應下文改進海豚回聲算法以最大值為優化目標的需要，將目標函數調整為：

(2)

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

1.2 擴展有限元建模

本文采用擴展有限元法對結構進行建模，計算出含有不同位置和大小的損傷結構的響應量。擴展有限元通過引入非連續的位移模式來代替傳統有限元的連續位移模式，使得對損傷結構不連續位移場的描述得到了簡化，在不用重新劃分網格的前提下可以方便地改變損傷的位置和大小。

擴展有限元中所應用的位移模式為：

(3)

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

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

(4)

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

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

(5)

(6)

2 海豚回聲算法及其改進

2.1 海豚回聲算法

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

(7)

式中：Fc1為預定義的第一個循環最優解為全局最優解的概率；pi表示第i個循環；N為總循環次數；p為收斂曲線的次數，可以通過調整收斂曲線的次數來調整收斂曲線的形狀。

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

DEO算法的求解流程為：

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

2) 計算當前循環的Fc.

3) 計算每一個可行解的目標函數值g1(Li)，其中目標函數以最大值為優化目標。

4) 計算每個變量的累計適應度：

(8)

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

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

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

(9)

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

6) 根據計算出的概率值生成下一循環的NL個可行解。

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

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

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

(10)

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

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

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

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

(11)

式中：k為縮放系數，當縮小時k>1，放大時k<1，不變時k=1.

5個搜索過程具體如下。

3) 平移搜索。當k1取值過大，導致全局最優解不在精搜索的搜索域之內，則精搜索求得的最優解會在搜索域邊緣，即當?j=1,2,…,NV使得：

(12)

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

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

(13)

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

當循環次數達到最大循環次數p時，循環結束，輸出最優解。程序框圖如圖3.

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

3 算例分析

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

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

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

3.1 單個孔洞損傷識別

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

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

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

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

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

3.2 兩個孔洞損傷識別

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

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

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

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

3.3 三個孔洞損傷識別

針對三個孔洞EDEO算法可行解個數NL=50，p1=50，最大循環次數p=200，待求解的優化參數為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算法部分循環識別結果，由圖可以看出EDEO算法經過150次迭代，孔洞1的圓心坐標為(2.185,2.461)，半徑為0.575；孔洞2的圓心坐標為(3.002,7.920)，半徑為0.592；孔洞3的圓心坐標為(7.898,7.989)，半徑為0.584，已經收斂到了不錯的結果。圖10給出了在15%噪聲下DEO算法的收斂結果，可以看出DEO算法已經無法收斂。表3給出了15%噪聲條件下EDEO算法的優化結果，結果表明對于3個孔洞EDEO算法能夠收斂到真實值。

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

圖9 15%噪聲下EDEO部分循環識別結果Fig.9 Identified process of EDEO with noise level of 15%

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

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

4 結論

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

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