基于量子粒子群算法的風(fēng)力機(jī)葉片損傷識(shí)別

孫芳錦， 盧 琛

(1.桂林理工大學(xué) 土木與建筑工程學(xué)院 廣西 桂林 541004； 2.廣西嵌入式技術(shù)與智能系統(tǒng)重點(diǎn)實(shí)驗(yàn)室， 廣西 桂林 541006)

最近幾十年來(lái),受自然災(zāi)害和養(yǎng)護(hù)不善、爆炸、沖擊等人為因素的影響,建筑結(jié)構(gòu)發(fā)生不同程度破壞和一定程度功能喪失,結(jié)構(gòu)的健康監(jiān)測(cè)顯得格外關(guān)鍵。用于結(jié)構(gòu)損傷鑒定的重要技術(shù)指標(biāo)有:固有頻率、模態(tài)結(jié)構(gòu)振型、柔度矩陣的變化、頻響函數(shù)、神經(jīng)網(wǎng)絡(luò)等。大部分結(jié)構(gòu)破壞的本質(zhì)是結(jié)構(gòu)剛度的降低，柔度和剛度成倒數(shù)關(guān)系，可以用作結(jié)構(gòu)鑒定損傷的指標(biāo)。

近年來(lái)，柔度矩陣與智能算法相結(jié)合的結(jié)構(gòu)損傷識(shí)別技術(shù)成為國(guó)內(nèi)外科學(xué)研究的熱點(diǎn)。許多學(xué)者使用柔度矩陣差、柔度差變化率和柔度曲率差來(lái)定位結(jié)構(gòu)的損傷[1-2]，但利用柔度靈敏度[3]和廣義柔度靈敏度[4]對(duì)結(jié)構(gòu)損傷進(jìn)行定量識(shí)別時(shí)，存在穩(wěn)定性差、誤差大的缺點(diǎn)。Zhao等[5]研究頻率、振型和柔度的靈敏度,證明了柔度矩陣比頻率和振型更為靈敏。

本文結(jié)合量子粒子群優(yōu)化算法和廣義柔度矩陣提出了一種新的損傷識(shí)別方法。

1 量子粒子群算法

Eberhart與Kennedy在1995年首先提出了粒子群算法(PSO)[6]，粒子群算法的數(shù)學(xué)模型是假設(shè)N維解空間結(jié)構(gòu)中存在M個(gè)微粒。第m個(gè)粒子在t時(shí)間上位移和移動(dòng)速度分別為

Xm(t)=(Xm,1(t),...Xm,n(t),...Xm,N(t)),Vm(t)=(Vm,1(t),...Vm,n(t),...Vm,N(t)),m=1,2,...M,n=1,2,...N。

當(dāng)前時(shí)間粒子群中的局部區(qū)域最優(yōu)預(yù)測(cè)位置和全域最優(yōu)預(yù)測(cè)位置分別為

Pm(t)=(Pm,1(t),...Pm,n(t),...Pm,N(t)),G(t)=(G1(t),...Gn(t),...GN(t))。

其中

(1)

G(t)=Pg(t)=(Pg,1(t),...Pg,n(t),...Pg,N(t))

(2)

式中：f[·]表示構(gòu)建的目標(biāo)函數(shù)適應(yīng)值g=argmin{f[Pm(t)]},1≤m≤M；g為一個(gè)位于全局最優(yōu)位置粒子對(duì)應(yīng)的下標(biāo)。

隨著時(shí)間t的增長(zhǎng)，在粒子群算法中M個(gè)微粒的位置變化和速率更新，分別為：

Vm,n(t+1)=Vm,n(t)+c1u1(t)[Pm,n(t)-Xm,n(t)]+c2u2(t)[Gn(t)-Xm,n(t)]

(3)

Xm,n(t+1)=Xm,n(t)+Vm,n(t+1)

(4)

式中：c1和c2分別表示為局部最優(yōu)和全局最優(yōu)加速因子,通常取(0,2)中的一個(gè)值；u1和u2是均勻分布在區(qū)間(0,1)中的隨機(jī)數(shù)。為了防止粒子在迭代過(guò)程中飛出搜索空間，必須依次確定Vm,n(t)、Xm,n(t)搜索的上、下界。

粒子群優(yōu)化算法收斂速度慢，搜索模式單一，空間受限，易局限于局部?jī)?yōu)化。為了解決的上述缺點(diǎn)，科學(xué)家們不斷探索。Sun在量子力學(xué)的理論角度上給出了量子粒子群算法(QPSO)[7]，融入了波函數(shù)與δ勢(shì)阱的有機(jī)組合理念等，確保了搜索空間可以充分涵蓋整個(gè)區(qū)域。QPSO搜索算法把粒子在空間出現(xiàn)的機(jī)率以波函數(shù)模的平方表示，從而保證以全機(jī)率尋找到全局的最優(yōu)解，通過(guò)Monte Carlo方法得

(5)

式中：p是勢(shì)阱吸引子；L是δ勢(shì)阱的特征長(zhǎng)度；u是均勻分布在區(qū)間(0,1)中的隨機(jī)數(shù)，即u∈U(0,1)。

隨著時(shí)間t的增加，QPSO算法中M個(gè)粒子的位置和速度更新如下：

Xm,n(t+1)=pm,n(t)±α*|Cn(t)-Xm,n(t)|*ln(1/um,n(t))

(6)

pm,n(t)=βn(t)*Pm,n(t)+(1-βn(t))*Gn(t)

(7)

mbest=(C1(t),C2(t),...CN(t))

(8)

QPSO算法的步驟：

步驟1：初始化粒子個(gè)體位置，并設(shè)置初始粒子個(gè)體的最優(yōu)值等；

步驟2：計(jì)算各粒子的適應(yīng)度值；

步驟3：比較各粒子的適應(yīng)度值大小，然后更新個(gè)體最優(yōu)與全局最優(yōu)；

步驟4：運(yùn)用公式(5)～公式(7)更新粒子的位置和速度，得到最新的位置和速度；

步驟5：重復(fù)步驟2～步驟4，直至達(dá)到停止的條件。

2 廣義柔度矩陣的基本原理

無(wú)阻尼振動(dòng)條件下，n個(gè)自由度結(jié)構(gòu)在自由振動(dòng)時(shí)的方程為

[M]x″+[k]x=0

(9)

式中：[K]、[M]分別表示結(jié)構(gòu)的剛度矩陣和質(zhì)量矩陣；x″、x分別代表結(jié)構(gòu)的加速度響應(yīng)和位移響應(yīng)。

由(9)式的運(yùn)動(dòng)方程，得到特征方程為：

([K]-λ[M])φ=0

(10)

式中：λ、φ分別表示方程特征值、特征向量；λ表示結(jié)構(gòu)各階固有頻率的平方。

將式(10)以矩陣的形式表示為：

[K][Ψ]=[M][Ψ][Λ]

(11)

式中：[Ψ]=[φ1,φ2,...φn]為n階模態(tài)振型矩陣；[Λ]=daig(λ1,λ2,...λn)為n階對(duì)角矩陣。

由振型的質(zhì)量歸一化特性得

[Ψ]T[M][Ψ]=[I]

(12)

式中：[I]為單位矩陣。

將式(11)兩邊同時(shí)乘[Ψ]T得

[Ψ]T[K][Ψ]=[Ψ]T[M][Ψ][Λ]

(13)

式(12)代入式(11)的右側(cè)得

[Ψ]T[K][Ψ]=[Λ]

(14)

式(14)兩側(cè)同時(shí)求逆，整理得

[K]-1=[Ψ][Λ]-1[Ψ]T

(15)

其中柔度矩陣是剛度矩陣的逆矩陣，即：

(16)

從方程(16)可以看出，柔度矩陣是關(guān)于結(jié)構(gòu)振動(dòng)頻譜和強(qiáng)度振型的向量函數(shù)，與頻率倒數(shù)的平方成正比。隨著模態(tài)階數(shù)的增加，模態(tài)參數(shù)對(duì)柔度矩陣的貢獻(xiàn)越來(lái)越小。模態(tài)柔度僅通過(guò)前幾階模態(tài)就能更好地估計(jì)結(jié)構(gòu)柔度矩陣[2]，有效地克服了工程設(shè)計(jì)中僅檢測(cè)前幾階模態(tài)和高階模態(tài)不準(zhǔn)確的問(wèn)題。

為了引入廣義柔度矩陣的概念[4]，在柔度矩陣基礎(chǔ)上進(jìn)行改進(jìn)，具體函數(shù)為：

(17)

式中：Q=0,1,2,...。當(dāng)Q=0時(shí)，[G]=[F]，廣義柔度矩陣與一般柔度矩陣相等；當(dāng)Q=1時(shí)，[G]=[Ψ][Λ]-2[Ψ]T。

3 目標(biāo)函數(shù)的建立

根據(jù)式(17)建立結(jié)構(gòu)損傷前后廣義柔度矩陣差為

(18)

式中：[Gu]、[Gd]分別代表結(jié)構(gòu)損傷前后廣義柔度矩陣；φru、φrd分別代表結(jié)構(gòu)破壞損傷前后模態(tài)振型；ωru、ωrd分別代表結(jié)構(gòu)破壞損傷前后模態(tài)固有頻率；m為小于n的正整數(shù)。

(19)

式中：θi為第i個(gè)單元損傷參數(shù)因子。

令

(20)

(21)

構(gòu)造目標(biāo)函數(shù)

Z(θ)=A-B

(22)

式中：θ=(θ1,θ2,...θn)，θ的取值范圍為0≤θ≤1?？紤]到QPSO算法粒子搜索的最佳取值范圍和結(jié)構(gòu)損傷大于90%時(shí)將喪失使用功能，將θ的取值范圍調(diào)整為0.1≤θ≤0.9。

變剛度法是將損傷部位的各單元的剛度減小，再與其他常規(guī)單元共同進(jìn)行計(jì)算。根據(jù)圣維南定律，這種損傷只能影響到周?chē)膮^(qū)域，而不會(huì)影響到更遠(yuǎn)的地方。如果單元選擇得足夠大，則損傷部位的單元?jiǎng)偠葍H會(huì)影響到損傷部位的單元?jiǎng)偠龋渌麊卧膭偠炔粫?huì)受到影響。用變剛度法模擬了葉片的損傷，在損傷單元的縱向、橫向和剪切模量均減少的情況下，利用各部位的剛度變化來(lái)模擬損傷部位的損傷，以同一部位的剛度變化來(lái)反映損傷的嚴(yán)重性。假定E和E*分別為葉片完好和損傷兩種狀態(tài)下的彈性模量，則反映損傷程度的單元?jiǎng)偠日蹨p系數(shù)x可表示為：

(23)

4 基于QPSO算法對(duì)風(fēng)力機(jī)葉片的損傷識(shí)別

葉片是風(fēng)力機(jī)的主要部件之一，很容易惡化和磨損，影響使用壽命。進(jìn)行定期檢查，在早期階段發(fā)現(xiàn)故障，能減輕這些損失。這里采用一種以損傷前后結(jié)構(gòu)的廣義柔度矩陣差構(gòu)造目標(biāo)函數(shù)來(lái)識(shí)別結(jié)構(gòu)損傷的方法。

選取20 kW風(fēng)力機(jī)單獨(dú)的一個(gè)5.532 m葉片如圖1所示,將該葉片劃分為92個(gè)單元如圖2所示，葉片材料屬性與參數(shù)如表1所示，實(shí)驗(yàn)數(shù)據(jù)[8]如表2所示。

圖1 20 kW風(fēng)力機(jī)葉片模型

圖2 20 kW風(fēng)力機(jī)葉片網(wǎng)格劃分示意圖

表1 20 kW風(fēng)力機(jī)葉片材料屬性與參數(shù)

表2 20 kW風(fēng)力機(jī)葉片實(shí)驗(yàn)數(shù)據(jù)

迭代次數(shù)=評(píng)價(jià)次數(shù)/種群數(shù)量[9-11]。在相同的種群數(shù)量下，評(píng)價(jià)次數(shù)與迭代次數(shù)成正比，達(dá)到相同的適應(yīng)度值，看所需的評(píng)價(jià)次數(shù)來(lái)評(píng)估算法的計(jì)算效率。即評(píng)價(jià)次數(shù)越少，迭代次數(shù)越少，耗時(shí)越少，計(jì)算速度越快。為了證明研究的實(shí)效性以及計(jì)算的優(yōu)越性，分別按變剛度法折減剛度為30%、50%、70%設(shè)置損傷工況，利用ABAQUS軟件提取單元?jiǎng)偠染仃?、整體剛度矩陣、質(zhì)量矩陣代入構(gòu)造的目標(biāo)函數(shù)計(jì)算得到相應(yīng)損傷工況的目標(biāo)函數(shù)，再用MATLAB得到每個(gè)損傷工況的單元損傷參數(shù)、評(píng)價(jià)次數(shù)與適應(yīng)度值的折線圖。分別用PSO、QPSO算法來(lái)進(jìn)行3種工況的損傷識(shí)別，如圖3～圖5所示。

(a)PSO算法 (b)QPSO算法圖3 工況1下20 kW風(fēng)力機(jī)葉片在不同算法下的收斂折線圖

(a)PSO算法 (b)QPSO算法圖4 工況2下20 kW風(fēng)力機(jī)葉片在不同算法下的收斂折線圖

(a)PSO算法 (b)QPSO算法圖5 工況3下20 kW風(fēng)力機(jī)葉片在不同算法下的收斂折線圖

選取5 MW風(fēng)力機(jī)單獨(dú)的一個(gè)61.5 m葉片如圖6所示，將葉片劃分為55 366個(gè)單元如圖7所示，結(jié)構(gòu)物理參數(shù)如表3所示。

圖6 5 MW風(fēng)力機(jī)葉片模型

圖7 5 MW風(fēng)力機(jī)葉片網(wǎng)格劃分示意圖

表3 5 MW風(fēng)力機(jī)葉片材料屬性與參數(shù)

同20 kW風(fēng)力機(jī)葉片的計(jì)算方法，得到識(shí)別結(jié)果如圖8～圖10所示。

(a)PSO算法 (b)QPSO算法圖8 工況1下5 MW風(fēng)力機(jī)葉片在不同算法下的收斂折線圖

(a)PSO算法 (b)QPSO算法圖9 工況2下5 MW風(fēng)力機(jī)葉片在不同算法下的收斂折線圖

(a)PSO算法 (b)QPSO算法圖10 工況3下5 MW風(fēng)力機(jī)葉片在不同算法下的收斂折線圖

5 結(jié)果分析

20 kW和5 MW風(fēng)力機(jī)葉片在不同工況下兩種算法的收斂坐標(biāo)如表4和表5所示。

表4 20 kW風(fēng)力機(jī)葉片不同算法在不同工況下的收斂坐標(biāo)

表5 5 MW風(fēng)力機(jī)葉片不同算法在不同工況下的收斂坐標(biāo)

在MATLAB中查看結(jié)果，每個(gè)工況每個(gè)單元損傷情況一致，只需設(shè)置一個(gè)單元損傷參數(shù)，對(duì)20 kW和5 MW風(fēng)力機(jī)葉片由兩種算法在不同工況下單元損傷參數(shù)如表6和表7所示。

表6 20 kW風(fēng)力機(jī)葉片不同算法在不同工況下的單元損傷參數(shù)

表7 5 MW風(fēng)力機(jī)葉片不同算法在不同工況下的單元損傷參數(shù)

由圖3～圖5、圖8～圖10和表4～表7分析得：

(1)QPSO算法比PSO算法識(shí)別損傷程度更精確，精度提高接近2%。

(2)由圖3可知?jiǎng)傞_(kāi)始適應(yīng)度值達(dá)到0.01時(shí)，2種算法的評(píng)價(jià)次數(shù)相差不大，隨著適應(yīng)度值達(dá)到0.001時(shí)，PSO算法評(píng)價(jià)次數(shù)達(dá)到QPSO算法的7倍，即迭代次數(shù)達(dá)到QPSO算法的7倍。收斂時(shí)的適應(yīng)度值大概為0.000 1左右，PSO算法評(píng)價(jià)次數(shù)達(dá)到QPSO算法的6倍。在圖8中適應(yīng)度值達(dá)到0.001時(shí)，兩種算法的評(píng)價(jià)次數(shù)相差不大，適應(yīng)度值達(dá)到0.000 1時(shí)，PSO算法評(píng)價(jià)次數(shù)達(dá)到QPSO算法的10倍。說(shuō)明當(dāng)達(dá)到相同的適應(yīng)度值時(shí)，QPSO算法比PSO算法收斂更快，計(jì)算效率更快。

(3)由圖4可知隨著風(fēng)力機(jī)葉片損傷程度加大后，當(dāng)適應(yīng)度值為0.001左右時(shí)，PSO算法評(píng)價(jià)次數(shù)達(dá)到QPSO算法的3倍，即迭代次數(shù)達(dá)到QPSO算法的3倍。在圖9中當(dāng)適應(yīng)度值為0.001左右時(shí)，PSO算法評(píng)價(jià)次數(shù)達(dá)到QPSO算法的2倍，即迭代次數(shù)達(dá)到QPSO算法的2倍，可以看出QPSO算法計(jì)算效率仍然大于PSO算法。

(4)隨著葉片損傷程度達(dá)到70%，當(dāng)適應(yīng)度值為0.01左右時(shí)，圖5中PSO算法評(píng)價(jià)次數(shù)達(dá)到QPSO算法的1倍左右，即迭代次數(shù)達(dá)到QPSO算法的1倍左右。圖10中當(dāng)適應(yīng)度值為0.01左右時(shí)，PSO算法評(píng)價(jià)次數(shù)達(dá)到QPSO算法的11倍左右，即迭代次數(shù)達(dá)到QPSO算法的11倍左右。說(shuō)明隨著葉片損傷程度增加，QPSO算法對(duì)于大型風(fēng)力機(jī)葉片的計(jì)算效率遠(yuǎn)大于PSO算法。

6 結(jié) 語(yǔ)

通過(guò)算法檢測(cè)不同風(fēng)機(jī)損傷程度工況的同型風(fēng)力機(jī)不同葉片,且雖然風(fēng)力機(jī)不同葉片分別劃分的檢測(cè)單元多與雜,具有一定的不準(zhǔn)確性,但采用QPSO檢測(cè)算法與PSO檢測(cè)算法都仍然能得到較好的準(zhǔn)確識(shí)別,表明采用QPSO檢測(cè)算法目前具有良好的魯棒性,QPSO檢測(cè)算法未來(lái)在實(shí)際應(yīng)用工程風(fēng)機(jī)損傷工況檢測(cè)領(lǐng)域具有廣闊的技術(shù)應(yīng)用性和前景。