基于UKF的結構動荷載識別方法與試驗驗證

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

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

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

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

目前，國內外學者已經提出多種線性或非線性結構參數識別的系統識別方法，如最小二乘算法[2]、卡爾曼濾波器算法[3-4]、蒙特卡洛算法[5-7]、以及以損傷因子為目標函數的荷載識別方法[8]。卡爾曼濾波器算法是一種有效識別線性結構參數的均值方法，對于線性結構有著較好的估計效果。為了克服其只能應用于線性結構系統的缺點，衍生出了等效線性化的近似方法，從而發展出可用于非線性系統識別的擴展卡爾曼濾波算法。而后為了改善卡爾曼濾波算法對非線性系統的識別精度，又發展了無跡卡爾曼濾波算法等方法，該方法避免了非線性函數線性化所帶來的誤差且不用計算復雜的雅可比矩陣，具有操作簡單、濾波精度高及易收斂等優點。基于切比雪夫標準正交多項式分解的結構系統荷載識別和參數識別方法已得到眾多學者研究與驗證[9-10]，在實際運用中具有一定可行性。當前，已有一部分學者通過數值仿真分析驗證了UKF對結構體系幾何非線性、材料非線性的適用性[11]。

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

1 UKF算法流程

假定非線性離散時間系統的狀態方程與觀測方程滿足以下形式：

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

(1)

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

(2)

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

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

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

(3)

(4)

(5)

(6)

其中：

(7)

(8)

(9)

(10)

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

(11)

(12)

(13)

(14)

(15)

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

(16)

(17)

(18)

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

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

2 荷載識別方法

2.1 間接荷載識別法

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

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

(19)

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

(20)

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

(21)

式中:狀態變量k和c的作用是識別結構系統的剛度和阻尼，如果不需要進行識別，可以移除，即只留下正交展開系數w。

2.2 直接荷載識別法

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

假設有線性離散系統如下：

xk+1=Axk+Bpk+wk

(22)

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

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

pk+1=pk+ηk

(23)

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

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

(24)

由此，擴展后的狀態方程為：

(25)

矩陣Aa定義為：

(26)

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

(27)

含有未知噪聲矢量vk的經典離散時間觀測方程由于擴展狀態量變為：

(28)

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

(29)

于是，擴展狀態量的狀態方程和觀測方程表達式為

(30)

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

3 仿真算例

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

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

(31)

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

(32)

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

(33)

圖1 單自由度結構模型

3.2 仿真結果

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

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

(34)

本算例正問題采用四階龍格庫塔求解，對于非線性結構系統，通過狀態空間方程的形式可以將二階結構運動微分方程轉換為一階運動微分方程，即：

(35)

當采用直接荷載識別法時，取狀態量為：

(36)

當采用間接荷載識別法時，取狀態量為：

(37)

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

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

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

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

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

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

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

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

(38)

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

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

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

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

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

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

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

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

表2 5%噪聲時范數相對誤差

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

4.1 無線傳感系統

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

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

圖6 無線傳感系統工作原理

圖7 傳感器內部構造示意圖

圖8 集成模塊示意圖

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

4.2 模型參數識別

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

表3 待識別參數初始值

參數識別結果,如圖9所示。

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

4.3 荷載識別及改進方法

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

(a) 剛度識別結果

(b) 阻尼識別結果

(c) 位移識別結果

(d) 加速度識別結果

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

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

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

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

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

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

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

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

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

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

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

(39)

其中增廣的狀態協方差矩陣變為：

(40)

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

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

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

5 結 論

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