李天之，劉 綱,b，張亮亮,b

(重慶大學(xué) a土木工程學(xué)院；b.山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，重慶 400045)

結(jié)構(gòu)參數(shù)和荷載的識別對結(jié)構(gòu)健康監(jiān)測、有限元模型更新等具有重要作用。在假設(shè)荷載已知的情況下，基于最小二乘估計(jì)[1-2]、擴(kuò)展卡爾曼濾波(EKF)[3-4]、無味卡爾曼濾波(UKF)[5-7]、粒子濾波(PF)[5,8-9]等眾多時(shí)域方法被廣泛應(yīng)用于非線性系統(tǒng)的參數(shù)識別。當(dāng)已知結(jié)構(gòu)參數(shù)時(shí)，可采用卡爾曼濾波[10-13]等對荷載激勵(lì)進(jìn)行辨識。在土木結(jié)構(gòu)服役過程中，結(jié)構(gòu)參數(shù)的變化和荷載激勵(lì)往往都是未知的，將結(jié)構(gòu)參數(shù)識別和荷載辨識分開進(jìn)行研究較難滿足土木工程的實(shí)際需求。為解決這一難題，各國學(xué)者逐步發(fā)展了結(jié)構(gòu)參數(shù)與荷載同步識別方法，丁勇等[14]運(yùn)用正交分解法建立荷載方程，帶入狀態(tài)空間方程，再運(yùn)用改進(jìn)UKF算法識別參數(shù)和荷載，在三層非線性框架的數(shù)值模擬中，可成功識別剛度、阻尼等結(jié)構(gòu)參數(shù)和隨機(jī)荷載；萬志敏等[15]在運(yùn)用改進(jìn)粒子濾波算法識別結(jié)構(gòu)參數(shù)的同時(shí)，運(yùn)用最小二乘估計(jì)求得荷載，數(shù)值模擬表明，該法同時(shí)識別參數(shù)和荷載的效果較好。以上解決方案都采用了2種不同的算法，流程相對復(fù)雜。同步識別結(jié)構(gòu)參數(shù)與荷載還面臨另一個(gè)難題，當(dāng)僅運(yùn)用加速度進(jìn)行長時(shí)間荷載識別時(shí)，其結(jié)果隨著時(shí)間的推移逐漸發(fā)生偏移，這是基于時(shí)域方法的荷載識別過程中難以避免的現(xiàn)象[15]。雖然可以將速度或位移作為觀測值來解決這一問題[5]，但在結(jié)構(gòu)健康監(jiān)測系統(tǒng)的實(shí)際應(yīng)用中，由于較難測量結(jié)構(gòu)的動(dòng)態(tài)速度或位移，一般僅將加速度作為觀測值；另一種可行的方法是通過額外的步驟，如虛擬測量技術(shù)對位移進(jìn)行調(diào)整[11]。

文中將結(jié)構(gòu)參數(shù)、荷載、速度和位移同時(shí)作為狀態(tài)值[16]，構(gòu)造狀態(tài)空間方程，采用粒子濾波算法直接對結(jié)構(gòu)參數(shù)與荷載進(jìn)行識別，步驟相對簡單。同時(shí)，針對荷載辨識結(jié)果漂移的問題，將加速度數(shù)據(jù)進(jìn)行梯形積分和零相位高通濾波，求得速度和位移的計(jì)算值，利用計(jì)算值來更新由粒子濾波求得的速度和位移狀態(tài)值，長時(shí)間確保荷載識別的準(zhǔn)確性。

1 狀態(tài)空間方程

一般來說，土木結(jié)構(gòu)的狀態(tài)空間方程為

(1)

其中，xk、uk和yk分別是第k時(shí)間步的狀態(tài)值(包括位移和速度)、荷載激勵(lì)和結(jié)構(gòu)響應(yīng)；ωk和νk分別是第k時(shí)間步的過程噪聲和觀測噪聲；f(·)和h(·)是非線性函數(shù)，描述了狀態(tài)值和結(jié)構(gòu)響應(yīng)隨時(shí)間變化的關(guān)系。

為實(shí)現(xiàn)結(jié)構(gòu)參數(shù)和荷載的同步識別，可將結(jié)構(gòu)參數(shù)mk、荷載激勵(lì)uk與原有的狀態(tài)值xk組成新的狀態(tài)向量zk[16]。

(2)

新的狀態(tài)向量zk與結(jié)構(gòu)響應(yīng)yk組成新的狀態(tài)空間方程：

(3)

式中，θk和ηk分別是結(jié)構(gòu)參數(shù)與荷載激勵(lì)的過程噪聲。文中所有噪聲均假設(shè)為零均值高斯噪聲。

2 粒子濾波與加速度更新方法

2.1 粒子濾波算法

以式(3)為例，粒子濾波的流程如下：

(4)

(5)

(6)

計(jì)算有效粒子數(shù)Neff：

(7)

(8)

2.2 加速度更新策略

土木結(jié)構(gòu)通常僅能得到加速度響應(yīng){ak:k=1,…,n}，采用梯形積分得到k時(shí)間步的速度vk：

(9)

其中，T是采樣周期。

3 數(shù)值模擬

3.1 三層框架結(jié)構(gòu)

采用如圖1所示的三層框架對文中方法進(jìn)行驗(yàn)證，每層的剛度s1=s2=s3=500 N/m，每層的質(zhì)量m1=m2=m3=20 kg。通過集中質(zhì)量矩陣和剛度矩陣計(jì)算得到前三階的自振頻率：f1=2.225 rad/s，f2=6.235 rad/s，f3=9.010 rad/s，假設(shè)第一階阻尼比ξ1、第三階阻尼比ξ3均為0.05，結(jié)構(gòu)阻尼為瑞利阻尼，可通過式(10)求得瑞利阻尼參數(shù)為α=0.178，β=0.009。

圖1 三層框架數(shù)值模擬

(10)

以NS方向的El-centro地震波生成外部荷載F3，水平作用在第三層上，最大值為34.87 N。

采樣頻率設(shè)置為200 Hz，通過Newmark-β法計(jì)算得到結(jié)構(gòu)前30 s的加速度數(shù)據(jù)。假定觀測噪聲是零均值高斯噪聲，其標(biāo)準(zhǔn)差為加速度響應(yīng)數(shù)據(jù)的5%均方根值。計(jì)算所用的荷載及各層加速度如圖2所示。

圖2 加速度與外部荷載

3.2 框架結(jié)構(gòu)的狀態(tài)空間方程

圖1所示三層框架結(jié)構(gòu)的動(dòng)力學(xué)方程為

(11)

當(dāng)對結(jié)構(gòu)剛度、荷載同時(shí)進(jìn)行識別時(shí)，公式(11)可轉(zhuǎn)換為離散化的狀態(tài)空間方程[5]：

(12)

其中，{ωi:i=1,…,6}、{θi:i=1,2,3}和η是過程噪聲，{υi:i=1,2,3}為觀測噪聲，T=0.005 s。圖1中的三層框架是線性結(jié)構(gòu)，由于式(12)中存在未知參數(shù)與未知速度、位移相乘的情況，該公式是非線性方程[5,9]，故采用粒子濾波進(jìn)行參數(shù)識別。

運(yùn)用梯形積分和截止頻率為0.1 Hz的零相位高通濾波，可通過加速度求得速度、位移計(jì)算值。在此數(shù)值模擬中，取Δt=1 s，即每隔1 s對速度、位移狀態(tài)值進(jìn)行調(diào)整。在粒子濾波算法中，重采樣的閾值NT取0.5，粒子個(gè)數(shù)N=16 000。初始剛度在400～800 N/m的范圍內(nèi)隨機(jī)取值，荷載、位移和速度的初始值均取0。假設(shè)各過程噪聲均為獨(dú)立的零均值高斯噪聲，根據(jù)試算，速度、位移的過程噪聲的標(biāo)準(zhǔn)差分別取速度、位移計(jì)算值的0.2%均方根值，剛度、荷載的過程噪聲的標(biāo)準(zhǔn)差分別取0.1、2.5。

3.3 識別結(jié)果分析

分別將不調(diào)整速度、位移狀態(tài)值以及按文中方法進(jìn)行調(diào)整后的結(jié)果繪制于圖3中，其縱坐標(biāo)代表在16 000個(gè)粒子中剛度或荷載的均值。

圖3 結(jié)構(gòu)參數(shù)與荷載辨識過程

將最后一步的均值作為參數(shù)識別結(jié)果，并用相對誤差(RE)來評價(jià)各參數(shù)識別結(jié)果的精度：

(13)

其中，se表示剛度識別結(jié)果；sr表示剛度真實(shí)值。

用均方誤差(MSE)來評價(jià)荷載識別的精度：

(14)

其中，F(xiàn)e(k)表示第k步的荷載識別值；Fr(k)表示第k步的荷載真實(shí)值；M是時(shí)間步長。

由圖3(a)和表1可知，當(dāng)不調(diào)整速度、位移狀態(tài)值時(shí)，1.5 s左右結(jié)構(gòu)參數(shù)收斂，識別結(jié)果較好；但在大約10 s后，荷載辨識值發(fā)生漂移，程度越來越嚴(yán)重。此時(shí)，剛度識別值已接近真實(shí)值，故剛度值的辨識誤差不是荷載出現(xiàn)漂移的原因。正如文獻(xiàn)[11]所述，這是僅用加速度數(shù)據(jù)進(jìn)行長時(shí)間荷載辨識過程中難以避免的現(xiàn)象。從圖3(b)和表1可知，對結(jié)構(gòu)參數(shù)而言，文中方法的參數(shù)識別收斂時(shí)間、精度和傳統(tǒng)方法基本一致；但對于荷載識別而言，全過程中未出現(xiàn)漂移現(xiàn)象，有效解決了傳統(tǒng)方法的不足。

表1 參數(shù)識別和荷載識別結(jié)果

由于初始值、噪聲的影響，識別結(jié)果具有一定的隨機(jī)性，在同樣的情況下再對同一組數(shù)據(jù)進(jìn)行9次識別，計(jì)算結(jié)果如表2所示。從表2可知，結(jié)構(gòu)參數(shù)識別的最大誤差為-14.60%，最小誤差為0.26%；荷載激勵(lì)識別的最大均方誤差為1.88，最小均方誤差為0.79。

表2 結(jié)構(gòu)參數(shù)與荷載10次同步識別的結(jié)果

3.4 調(diào)整時(shí)間間隔對識別的影響

為研究調(diào)整時(shí)間間隔Δt對結(jié)構(gòu)參數(shù)、荷載識別的影響，選取{Δt=0.1, 0.3, 0.5, 1, 1.5, 2, 3}共7種情況各進(jìn)行10次運(yùn)算，不同Δt時(shí)的參數(shù)相對誤差的絕對值的均值、10次荷載識別均方誤差的均值如圖4所示。

圖4 不同時(shí)間間隔Δt的識別情況

由圖4可知，當(dāng)Δt較短(如0.1 s、0.3 s)時(shí)，第3層的剛度s3的誤差和荷載識別的均方誤差相對較大。因?yàn)樵谠摃r(shí)段內(nèi)結(jié)構(gòu)參數(shù)識別尚未收斂，此時(shí)用帶有誤差的速度、位移計(jì)算值來調(diào)整狀態(tài)值，參數(shù)識別受到的影響相對較大；同時(shí)，誤差相對較大的s3會對荷載識別造成影響。隨著Δt逐漸增大，各剛度識別的誤差基本保持不變，但荷載識別的均方誤差會逐漸增大。各參數(shù)趨于收斂或已收斂，計(jì)算值的誤差對參數(shù)識別的影響相對較?。煌瑫r(shí)，由于不能及時(shí)修正速度和位移，荷載識別值會逐漸漂移，故其均方誤差會隨時(shí)間間隔的增大而增大。綜上所述，為同時(shí)保持結(jié)構(gòu)參數(shù)和荷載識別的精確性，應(yīng)合理設(shè)置Δt，不宜過長或過短，建議以參數(shù)收斂所需的時(shí)間(數(shù)值模擬算例的1.5 s)來設(shè)置間隔Δt。

4 結(jié)束語

為了土木工程結(jié)構(gòu)不僅能測試結(jié)構(gòu)的加速度響應(yīng)，同時(shí)還能識別結(jié)構(gòu)參數(shù)和輸入荷載，文中將結(jié)構(gòu)參數(shù)和荷載擴(kuò)充為狀態(tài)值并建立狀態(tài)空間方程，基于粒子濾波算法進(jìn)行結(jié)構(gòu)參數(shù)和荷載的同步識別。同時(shí)，利用測試的加速度響應(yīng)，通過梯形積分和零相位高通濾波計(jì)算各時(shí)間步的速度、位移，然后將該值以一定時(shí)間間隔Δt用于更新對應(yīng)狀態(tài)識別值。三層框架結(jié)構(gòu)的數(shù)值算例結(jié)果表明：

1)文中方法可同時(shí)識別結(jié)構(gòu)參數(shù)與荷載，識別結(jié)果較為準(zhǔn)確、穩(wěn)定。

2)采用積分和零相位高通濾波所得的速度、位移計(jì)算值來調(diào)整速度、位移狀態(tài)值可克服荷載辨識值出現(xiàn)偏移這一問題。

3)應(yīng)合理設(shè)置速度、位移的時(shí)間間隔Δt，不宜過長或過短，建議根據(jù)結(jié)構(gòu)參數(shù)識別收斂所需的時(shí)間來設(shè)置Δt的大小。