基于UKF的結(jié)構(gòu)動荷載識別方法與試驗驗證

郭麗娜， 宋開明， 張延哲， 丁 勇,3

(1.東北農(nóng)業(yè)大學 水利與土木工程學院，哈爾濱 150030；2.哈爾濱工業(yè)大學 土木工程學院，哈爾濱 150090；3.中電建冀交高速公路投資發(fā)展有限公司，石家莊 050000)

結(jié)構(gòu)動力學問題可分為兩大類：若已知結(jié)構(gòu)外界輸入，需要通過結(jié)構(gòu)試驗或數(shù)值仿真分析研究結(jié)構(gòu)響應和性能，此類問題稱為正問題；若已知結(jié)構(gòu)響應，需估計結(jié)構(gòu)外荷載、識別結(jié)構(gòu)參數(shù)，此類問題統(tǒng)稱為反問題[1]。動力荷載對結(jié)構(gòu)的破壞常為靜力荷載破壞作用的2倍～4倍。結(jié)構(gòu)外荷載常隨時間變化且具有復雜的隨機性，這使得實際工程結(jié)構(gòu)動力荷載的識別問題成為工程界亟待解決的問題之一。

通過直接測量和系統(tǒng)識別方法評估兩種方法均可獲得外荷載時程，但在實際工程中，直接測量的方法常受到測點數(shù)目和測點位置的限制難以應用。因此，通過系統(tǒng)識別的方法進行結(jié)構(gòu)外荷載識別顯得尤為重要。

目前，國內(nèi)外學者已經(jīng)提出多種線性或非線性結(jié)構(gòu)參數(shù)識別的系統(tǒng)識別方法，如最小二乘算法[2]、卡爾曼濾波器算法[3-4]、蒙特卡洛算法[5-7]、以及以損傷因子為目標函數(shù)的荷載識別方法[8]?？柭鼮V波器算法是一種有效識別線性結(jié)構(gòu)參數(shù)的均值方法，對于線性結(jié)構(gòu)有著較好的估計效果。為了克服其只能應用于線性結(jié)構(gòu)系統(tǒng)的缺點，衍生出了等效線性化的近似方法，從而發(fā)展出可用于非線性系統(tǒng)識別的擴展卡爾曼濾波算法。而后為了改善卡爾曼濾波算法對非線性系統(tǒng)的識別精度，又發(fā)展了無跡卡爾曼濾波算法等方法，該方法避免了非線性函數(shù)線性化所帶來的誤差且不用計算復雜的雅可比矩陣，具有操作簡單、濾波精度高及易收斂等優(yōu)點?；谇斜妊┓驑藴收欢囗検椒纸獾慕Y(jié)構(gòu)系統(tǒng)荷載識別和參數(shù)識別方法已得到眾多學者研究與驗證[9-10]，在實際運用中具有一定可行性。當前，已有一部分學者通過數(shù)值仿真分析驗證了UKF對結(jié)構(gòu)體系幾何非線性、材料非線性的適用性[11]。

本文利用基于UKF算法的直接荷載識別法對作用于地震模擬振動臺的外部激勵進行識別，并針對出現(xiàn)漂移現(xiàn)象的荷載識別結(jié)果提出了位移-加速度聯(lián)合觀測值和噪聲水平自適應評估兩種解決方法，改善了識別過程的收斂性，驗證了所提方法的合理性和有效性。

1 UKF算法流程

假定非線性離散時間系統(tǒng)的狀態(tài)方程與觀測方程滿足以下形式：

xk=f(xk-1,uk-1,k-1)+vk-1

(1)

yk=h(xk,uk,k)+wk

(2)

式中:k∈R為離散時間，xk∈Rn為n維的隨機狀態(tài)向量，yk∈Rm為m維觀測向量。假設非線性狀態(tài)方程f(·)和非線性觀測方程h(·)在狀態(tài)量xk處連續(xù)并且可微。vk～N(0,Qk)是過程噪聲，滿足高斯分布，wk～N(0,Rk)為觀測噪聲，同樣滿足高斯分布。

無跡卡爾曼濾波器算法的遞推過程如下：

(1) 預測步：構(gòu)造2n+1個sigma點，準則如下：

(3)

(4)

(5)

(6)

其中：

(7)

(8)

(9)

(10)

(2) 更新步：利用UT變換和觀測量可以計算得到下一步狀態(tài)量的均值與方差：

(11)

(12)

(13)

(14)

(15)

計算卡爾曼增益并更新狀態(tài)量均值與方差得：

(16)

(17)

(18)

式中：yk表示第k步的觀測值。

通過上述步驟的循環(huán)遞推運算，可以發(fā)現(xiàn)UKF算法不需要對非線性函數(shù)進行線性變換，避免了復雜的雅克比矩陣的運算。

2 荷載識別方法

2.1 間接荷載識別法

由于外部荷載的不確定性及其幅值的時變性，很難在狀態(tài)空間內(nèi)使用確定性的模型來表征外部荷載?？紤]到荷載時程是一個隨機過程，即可以通過正交分解的辦法將荷載時程分解成一系列系數(shù)和多項式正交基乘積的線性組合形式，并且這些系數(shù)之間是相互獨立的。這種分解方法可以對隨機過程進行有效的量化擬合，方便建立數(shù)學模型進行數(shù)學分析。

在分解荷載時程的過程中，發(fā)現(xiàn)基于標準正交基的隨機過程分解方法效果很好[12-13]，同時還可以避免求解復雜的Fredholm積分方程。研究發(fā)現(xiàn)，切比雪夫(Chebyshev)分解方法是適用于隨機輸入正交分解的較準確的分解方法[14-15]。切比雪夫正交多項分解如下所示：

(19)

式中:Linp是輸入長度，Nm是分解階數(shù)，基于切比雪夫標準多項式分解，結(jié)構(gòu)運動方程可以寫作：

(20)

即基于荷載正交分解的荷載識別UKF算法的狀態(tài)空間方程可以寫作如下：

(21)

式中:狀態(tài)變量k和c的作用是識別結(jié)構(gòu)系統(tǒng)的剛度和阻尼，如果不需要進行識別，可以移除，即只留下正交展開系數(shù)w。

2.2 直接荷載識別法

由于間接荷載識別法需要識別的參數(shù)很多，并且識別出來的只是正交展開系數(shù)，還需多一步系數(shù)與正交基的乘積求和才能得到待識別荷載，除此之外項數(shù)的選取也是一個問題，很難統(tǒng)一定論，所以為了減少待識別的參數(shù)個數(shù)以及簡化運算過程，可以將待識別的外部荷載直接作為狀態(tài)變量進行識別，即直接荷載識別法。

假設有線性離散系統(tǒng)如下：

xk+1=Axk+Bpk+wk

(22)

式中:A和B是系統(tǒng)矩陣，wk(wk∈Rns，0≤k≤∞)是隨機過程噪聲，ns是狀態(tài)量個數(shù)。

為了補充式(22)，假設在時間k以及k+1處有一個荷載矢量p∈Rnp：

pk+1=pk+ηk

(23)

式中:ηk是隨機過程成分，ηk∈Rnp，0≤k≤∞。通過合理選擇過程噪聲ηk的協(xié)方差矩陣，利用式(23)可以估計荷載時程。

通過結(jié)合式(22)和(23)，并且重新定義狀態(tài)量xa(上標a代表擴展)得：

(24)

由此，擴展后的狀態(tài)方程為：

(25)

矩陣Aa定義為：

(26)

噪聲矢量ζk∈Rns+np取決于模型誤差wk∈Rns以及荷載隨機增量ηk∈Rnp：

(27)

含有未知噪聲矢量vk的經(jīng)典離散時間觀測方程由于擴展狀態(tài)量變?yōu)椋?/p>

(28)

式中:d∈Rnd,代表測量數(shù)據(jù)矢量；矩陣Ga∈Rnd×(ns+np)是由輸出影響矩陣G和直接轉(zhuǎn)換矩陣J組合而成的：

(29)

于是，擴展狀態(tài)量的狀態(tài)方程和觀測方程表達式為

(30)

式中:用于識別的系統(tǒng)矩陣A、B、G和J可以根據(jù)不同方法獲得。

3 仿真算例

3.1 單自由度滯回非線性剪切模型

本部分所用模型為單自由度滯回非線性模型，如圖1所示，非線性成分采用Bouc-Wen模型模擬，在只考慮結(jié)構(gòu)單方向運動時，其運動方程可表示為：

(31)

式中，RT(x,z,t)=αkx+(1-α)kz(t)，并有

(32)

式中:z為滯回位移成分;m,c和k分別為結(jié)構(gòu)系統(tǒng)的質(zhì)量、阻尼和剛度;α為系數(shù)，用來描述結(jié)構(gòu)滯回特性與彈性程度的比值，取值范圍[0 1]。當α取1時，表示屈服后剛度與屈服前剛度相等，即結(jié)構(gòu)為純彈性；當α取0時，表示系統(tǒng)總是處于屈服的狀態(tài)，即此時認為系統(tǒng)是純滯回性的。在此算例中取α=0，即：

(33)

圖1 單自由度結(jié)構(gòu)模型

3.2 仿真結(jié)果

采用的模型參數(shù)取值為：m=0.1 t、c=0.32 kN/(m/s)和k=38 kN/m，非線性參數(shù)取值為：β=4、γ=4和n=2，外部荷載作用在該單自由度質(zhì)量上。其中外荷載考慮了兩種情況，分別為F(t)=A×sin(4×π×t)kN的周期荷載以及沖擊荷載。其中A為在一定范圍內(nèi)線性變化的常數(shù)，用來模擬時變荷載，假定荷載的作用時長是3 s，時間步長為0.001 s。

模擬過程中采用的噪聲為百分比噪聲，施加原則如下：

(34)

本算例正問題采用四階龍格庫塔求解，對于非線性結(jié)構(gòu)系統(tǒng)，通過狀態(tài)空間方程的形式可以將二階結(jié)構(gòu)運動微分方程轉(zhuǎn)換為一階運動微分方程，即：

(35)

當采用直接荷載識別法時，取狀態(tài)量為：

(36)

當采用間接荷載識別法時，取狀態(tài)量為：

(37)

通過數(shù)值模擬，我們得到仿真結(jié)果,如圖2和圖3所示。

圖2和圖3分別表示無噪聲和5%噪聲時，周期荷載作用下，兩種識別方法的識別結(jié)果。采用范數(shù)相對誤差對識別精度加以考慮：

(a) 直接荷載識別法-無噪聲

(b) 間接荷載識別法-無噪聲-15階

(c) 間接荷載識別法-無噪聲-20階

(a) 直接荷載識別法-5%噪聲

(b) 間接荷載識別法-5%噪聲-15階

(c) 間接荷載識別法-5%噪聲-20階

(38)

圖4和圖5給出了沖擊荷載作用下，兩種識別方法的對比情況，并考慮了噪聲的影響。

對于上述分析結(jié)果，將求得的直接荷載識別法與間接荷載識別法的范數(shù)相對誤差,見表1和表2。

從圖分析結(jié)果可以看出，間接荷載識別法精度受正交多項式分解階數(shù)的影響較為明顯，一般情況下，當正交分解階數(shù)增大時，間接荷載識別法的精度會顯著提高，與此同時，外荷載矩陣及協(xié)方差矩陣也隨之增大，因此計算效率有所降低，由此可見，對間接荷載識別法，確定恰當?shù)恼环纸怆A數(shù)是協(xié)調(diào)識別精度和計算效率的重要環(huán)節(jié)。直接荷載識別法，計算效率較高，抗干擾能力強，但識別精度不容易控制。綜上，直接荷載識別法應用于單自由度非線性模型荷載識別的合理性和有效性得以驗證。

(a) 直接荷載識別法-無噪聲

(b) 間接荷載識別法-無噪聲-40階

(a) 直接荷載識別法-5%噪聲

(b) 間接荷載識別法-5%噪聲-40階

周期荷載沖擊荷載間接識別法直接識別法間接識別法直接識別法正交分解階數(shù)1520———3040———范數(shù)相對誤差/%22.832.271.2414.869.7213.43%

表2 5%噪聲時范數(shù)相對誤差

4 基于無線傳感的振動臺荷載識別試驗驗證

4.1 無線傳感系統(tǒng)

本研究所使用的無線傳感系統(tǒng)為所在團隊自主集成研發(fā)。其中，傳感器節(jié)點有各自的物理IP，通過自身配置的WiFi模塊與無線路由器連接，PC機通過網(wǎng)線或者WiFi連接到路由器上，測試數(shù)據(jù)由傳感器節(jié)點通過WiFi在線傳輸給PC機。傳感器能量供應方式有兩種：外電源連接和電池連接。無線傳感系統(tǒng)工作原理示意圖,如圖6所示。

傳感器節(jié)點主要由加速度傳感器芯片、ARM處理器和WiFi模塊組成，如圖7所示。本系統(tǒng)最大采樣精度12bits，有三種量程±2 G、±4 G、±8 G，可根據(jù)應用場合選擇，最大采樣頻率800 Hz。該傳感器具有封裝小、功耗低、安裝方便等特點。ARM處理器選用主頻72 MHz，片內(nèi)集成20 KRAM和64K FLASH，以及多達20個GPIO，適合用于低成本、低功耗的嵌入式領域。WiFi模塊選用業(yè)界成熟的simple WiFi，支持802.11b/g標準，支持AP和STA模式，配置簡單，安全性高。電源選擇簡單、安裝方便，集成模塊示意圖如圖8所示。為了保證數(shù)據(jù)傳輸?shù)目煽啃裕瑐鞲衅鞴?jié)點和PC之間使用自定義的協(xié)議進行通信。PC機可向傳感器節(jié)點發(fā)送啟動、停止、配置工作模式等命令，傳感器節(jié)點除了響應上述命令外，還可以將采集到的數(shù)據(jù)以及節(jié)點狀態(tài)實時發(fā)送給PC。

圖6 無線傳感系統(tǒng)工作原理

圖7 傳感器內(nèi)部構(gòu)造示意圖

圖8 集成模塊示意圖

實驗時，傳感器連接無需大量的導線連接，與常規(guī)的有線連接測試相比，安裝方便，節(jié)省了大量的布線時間。在成本方面，本系統(tǒng)無需傳統(tǒng)的電荷放大器和采集儀，大大降低了實驗成本。本系統(tǒng)的使用還可避免導線振動以及連接因素等導致的測試數(shù)據(jù)的污染。

4.2 模型參數(shù)識別

由于振動臺[16]的臺面剛度很大，為了計算方便，將其計算模型簡化為單自由度結(jié)構(gòu)。對抗震實驗室振動臺進行空臺加載，通過測量振動臺臺面位移，利用UKF算法識別地震模擬振動臺的剛度和阻尼參數(shù)，其中外部激勵荷載選取El-Centro(NS, 1940)地震波，并且假定為已知力。振動臺質(zhì)量已知且假定保持不變，將待識別參數(shù)初始值,見表3。

表3 待識別參數(shù)初始值

參數(shù)識別結(jié)果,如圖9所示。

從圖9可以看出，地震模擬振動臺的剛度約為45.553 3 kN/m，阻尼約為1.540 1 kN/(m/s)。加速度實驗值與位移實驗值是通過實驗實測得到，且與識別值吻合較好，說明剛度與阻尼識別值可信度較高。

4.3 荷載識別及改進方法

在上述已識別結(jié)構(gòu)參數(shù)的基礎上，假定外部激勵即El-Centro地震波未知。僅以加速度(分別用無線傳感器和有線傳感器測得)為觀測值，利用直接荷載識別法進行外荷載識別，識別結(jié)果見圖10。

(a) 剛度識別結(jié)果

(b) 阻尼識別結(jié)果

(c) 位移識別結(jié)果

(d) 加速度識別結(jié)果

(a) 基于無線加速度的識別結(jié)果

(b) 基于有線加速度的識別結(jié)果

從圖10可以看出，僅以加速度作為觀測值進行外部荷載識別會導致識別結(jié)果漂移，這是由于加速度值積分過程將造成誤差累積，從而使識別過程不穩(wěn)定，識別結(jié)果收斂性差。

針對上述漂移現(xiàn)象，我們提出了兩種解決途徑：

(1) 加速度與位移聯(lián)合作為觀測值

由于我們直接從實驗中測得加速度值與位移值，位移值會對加速度值積分的累積誤差加以約束，從而改善了識別過程的收斂性。從識別結(jié)果我們可以發(fā)現(xiàn)，聯(lián)合位移-加速度響應下的漂移現(xiàn)象有了明顯改善；無線加速度傳感器要優(yōu)于有線加速度傳感器，因為有線加速度傳感器的接口部位容易擾動，進而產(chǎn)生較大的噪聲干擾。

(a) 基于無線加速度-位移的識別結(jié)果

(b) 基于有線加速度-位移的識別結(jié)果

Fig.11 Identified results of external force based on acceleration and displacement

(2) 噪聲水平的自適應評估

噪聲水平的自適應評估指的是將過程噪聲和觀測噪聲與原始狀態(tài)量一起組合成一個新的狀態(tài)量，因為當噪聲太大時，UKF識別的精度不高，而且不易收斂，通過將噪聲引入到狀態(tài)量中，實現(xiàn)噪聲水平的自適應評估，可以確保UKF識別過程的穩(wěn)定性和魯棒性。新的狀態(tài)量如下所示：

(39)

其中增廣的狀態(tài)協(xié)方差矩陣變?yōu)椋?/p>

(40)

除了由增廣噪聲構(gòu)建的新的狀態(tài)量和增廣狀態(tài)協(xié)方差矩陣的不同之外，其余計算流程均和上述UKF算法一致，且只使用加速度作為觀測值。從識別結(jié)果可以看出增廣噪聲法改善了識別結(jié)果的穩(wěn)定性與收斂性，且無線加速度識別優(yōu)于有線加速度識別。

(a) 使用無線加速度的識別結(jié)果

(b) 使用有線加速度的識別結(jié)果

5 結(jié) 論

本文從基于UKF算法的兩種荷載識別方法出發(fā)，首先對單自由度滯回非線性剪切模型進行仿真，分別利用直接荷載識別法和基于荷載正交分解的間接荷載識別法對該模型進行兩種外荷載識別，分別為周期荷載和沖擊荷載。驗證了直接荷載識別法應用于單自由度非線性模型荷載識別的合理性和有效性。而直接荷載識別法對于多自由度非線性模型仍是適用的，課題組人員正在從事相關領域的研究。隨后利用直接荷載法對作用于地震模擬振動臺上的未知地震波進行識別，并針對識別結(jié)果發(fā)生漂移的現(xiàn)象提出了兩種直接擴展狀態(tài)量的荷載識別方法，識別結(jié)果顯示，兩種方法對荷載識別的漂移現(xiàn)象有著明顯的改善。