基于云模型和離散卡爾曼濾波的臨時看臺結構動荷載識別

陳文博, 郭惠勇

(1. 重慶大學 土木工程學院,重慶 400045; 2. 四川電力設計咨詢有限責任公司,成都 610016)

臨時看臺結構由于其施工靈活性高、周期短、可回收循環利用等優點,被廣泛應用于各種大型文體活動中。當看臺結構上人群擁擠、搖擺以及跳躍時,其容易發生較強振動導致看臺倒塌等事故的發生[1]。例如1994年,倫敦的一場音樂會上,由于觀眾隨著音樂搖擺、跳躍,導致看臺發生倒塌,超過50人受傷。2010年9月,在巴西舉辦的一場賽車比賽中,超過500名觀眾到場觀戰。由于觀眾有節奏搖擺、跳躍,臨時看臺結構發生倒塌,現場百余人受傷;2012年2月,王菲重慶演唱會中,當王菲入場時,大部分觀眾整體起立歡呼雀躍,觀眾席隨之發生倒塌事故,造成60余人受傷,演唱會臨時取消;2021年3月22日,在印度第47屆全國青年錦標賽開幕儀式上,臨時看臺突然倒塌,超過100名觀眾受傷。2021年5月,在以色列的一座猶太教堂,由于看臺上觀眾過多導致臨時看臺坍塌,造成200余人受傷,數人死亡。可以看出,人群荷載的作用,包括人群的擁擠和跳躍、搖擺等產生的動荷載常常是看臺等臨時設施結構安全的重要影響因素[2]。而國內外學者對看臺等臨時結構荷載識別工作研究較少,因此有必要對此類結構進行人群等荷載的快速識別研究。

當前針對臨時看臺結構的研究主要集中于看臺結構優化設計[3]、安全評價[4]、看臺承載能力[5]、節點剛度[6]、風荷載分析[7]等方面,而對看臺結構荷載識別方法的研究相對較少。其中,Celik等[8]提出了一種計算機視覺方法,用于測量由于個人和人群跳躍和擺動引起的荷載時程。該方法包括使用基于光流算法跟蹤個人和人群的位移軌跡以及實驗室內的測力板測試單人和人群跳躍搖擺產生的荷載,并通過現場演示進一步驗證。因此,本文針對某比賽實際看臺的等效簡化模型進行了荷載識別研究。首先描述了云模型和離散卡爾曼濾波的基本原理。隨后,基于離散卡爾曼濾波推導出動荷載識別算法,并利用云模型技術進行優化得到基于云模型和卡爾曼濾波的動荷載識別方法。其中云模型技術有效解決了測量誤差、環境噪聲等帶來的干擾問題。最后,將該方法應用于臨時看臺結構的數值仿真,驗證該方法對臨時看臺結構動荷載識別的有效性。

1 云模型和離散卡爾曼濾波基本原理

1.1 云模型和數字特征

1995年,中國工程院院士李德毅結合模糊數學和概率論的相關知識提出云模型的概念,主要用于定量數值與定性概念之間的互相轉換,可有效解決定性和定量相結合的不確定性問題[9]。

云的數字特征主要包括期望Ex、熵En和超熵He,云的3個數字特征可以描述云模型的整體特征,從而實現定性和定量之間的相互轉化。3個數字特征的定義如下[10]。

期望Ex:期望值是云的中心位置反映在云圖橫坐標上的值,表示該點為論域中最能夠表示定性概念的點。

熵En:熵最早是一個熱力學概念,本質是一個系統的“無秩序程度”。在云模型中,熵用來表示定性概念的不確定性的度量:一方面,它反映了定性概念模糊性的裕度;另一方面,通過云滴的分布位置的稀疏和稠密反映了定性概念在此空間中出現的概率,即為該定性概念的隨機性。

超熵He:超熵可以理解為熵的熵,體現了云滴之間的離散程度,用來度量熵的不確定性。

1.2 離散卡爾曼濾波

理論上,所有基本物理時間均在連續時間內發生,但在實際研究工作中,為了方便計算機的求解,需要將連續時間系統進行離散化,即在計算機的數值計算研究中,采用的是時間離散型卡爾曼濾波器。離散卡爾曼濾波器不需要計算機存儲大量數據,并且在輸入新的觀測數據后,可以迅速得出新的狀態估計值,具有很高的時效性且容易實現[11]。

離散系統主要用以下兩個基本方程來描述。

系統狀態方程

X(k)=ΦX(k-1)+Bω(k-1)

(1)

系統觀測方程

Z(k)=HX(k)+υ(k)

(2)

式中:X(k)為kΔt時刻的系統狀態向量;k為采樣時刻;Δt為數據的采樣時間間隔;Φ為狀態轉移矩陣;H為觀測矩陣;B為驅動矩陣;ω(k)和υ(k)分別為系統噪聲和測量噪聲,代表系統建模的誤差和傳感器測量的誤差。

噪聲特性假設如下。

(1)ω(k)和υ(k)均為高斯白噪聲時間序列,即

E[ω(k)]=0

(3)

E[ω(k)ωT(l)]=Q(k)δkl

(4)

E[υ(k)]=0

(5)

E[υ(k)υT(l)]=R(k)δkl

(6)

式中:Q(k)和R(k)為協方差矩陣;δkl為kronecker-δ函數,其定義如下

(7)

(2)ω(k)和υ(k)互不相關,即

E[ω(k)υT(k)]=0

(8)

(3)ω(k)和υ(k)都與初始狀態X(0)不相關,即

(9)

離散卡爾曼濾波器的工作原理[12-13]為:首先采用前一時刻的狀態估計值遞推計算當前時刻的狀態估計值,然后通過反饋當前響應的實際測量值的方法來修正當前的狀態估計值,實質上是一個不斷重復“預測-實測-修正”的順序遞推。算法流程分為時間更新和測量更新。其中,時間更新起到預測作用,包括當前時刻的狀態矩陣和協方差矩陣的估計,主要依靠系統狀態方程和噪聲參數,目的是為下一時刻的狀態估計提供先驗估計;測量更新在時間更新的基礎上,利用響應的實際測量值與當前時刻的狀態值之間的線性關系,負責反饋并校正當前估計值,同時也為下一時刻提供修正后的后驗估計。

離散卡爾曼濾波器具體公式如下所示。

時間更新

(10)

P(k/k-1)=ΦP(k-1/k-1)ΦT+BQBT

(11)

測量更新

J(k)=HP(k/k-1)HT+R

(12)

Ka(k)=P(k/k-1)HTJ-1(k)

(13)

P(k/k)=[I-Ka(k)H]P(k/k-1)

(14)

(15)

(16)

2 基于云模型和離散卡爾曼濾波的結構動荷載識別方法

2.1 基于離散卡爾曼濾波的動荷載識別方法

基于離散卡爾曼濾波的基本原理,本節從離散系統的角度研究動荷載反演。

對于一個有n個自由度的線性結構,其動態系統的運動方程為

(17)

式中:M,K,C分別為質量矩陣、剛度矩陣、阻尼矩陣;Y(t)為位移響應;F(t)為荷載列陣。將上述運動方程轉化為狀態方程,可以寫為如下形式

(18)

(19)

將狀態空間狀態方程轉化為模態空間狀態方程可得

(20)

(21)

(22)

(23)

式中:Φi,fi,Ci,Ki,ξi,ωi分別為按質量規準化的第i階振型、模態荷載、阻尼、剛度、阻尼比、自振頻率。

式(20)為矩陣形式的微分方程,其解析解為

Xi(tk+1)=Φi(k+1/k)Xi(tk)+Гi(k+1/k)fi(tk)

(24)

Φi(k+1/k)=exp[AiΔt]

(25)

(26)

式中:Г為驅動矩陣,當Δt取值一定時,Φi(k+1/k)和Гi(k+1/k)均為常量矩陣,所以式(24)可簡寫為

Xi(tk+1)=ΦiXi(tk)+Гifi(tk)

(27)

以結構的位移、速度或加速度響應為觀測值,觀測方程為如下形式

Zi(t)=HiXi(t)+Difi+εi

(28)

式中:Zi(t)為觀測值;Hi為觀測矩陣;Di為系統矩陣;εi為觀測噪聲。

當觀測響應為位移時

Hi=[10],Di=0

(29)

當觀測響應為速度時

Hi=[01],Di=0

(30)

當觀測響應為加速度時

Hi=[-Ki-Ci],Di=1

(31)

對比式(28)和離散卡爾曼濾波基本方程可得:觀測方程中多出一個未知輸入項Difi,不符合濾波基本方程的構造形式,無法直接對系統狀態響應Xi(t)進行估計,需要對濾波方程進行擴展。

系統狀態空間方程

X(k)=ΦX(k-1)+Гω(k-1)

(32)

觀測方程

Z(k)=HK(k)+Dω(k)+ε(k)

(33)

式中,ω(k)和ε(k)分別為系統噪聲和觀測噪聲,假設與1.2節相同。

令

υ(k)=Dω(k)+ε(k)

(34)

則觀測方程為

Z(k)=HX(k)+υ(k)

(35)

將以上結果代入式(10)～式(16),得到

(36)

G(k-1)=ΓQ(k-1)DT+
[DQ(k-1)DT+R(k-1)]-1

(37)

P(k/k-1)=[Φ-G(k-1)H]P(k-1/k-1)·

[Φ-G(k-1)H]T+ΓQΓT-G(k-1)DQΓT

(38)

J(k)=HP(k-1/k-1)HT+DQDT+R

(39)

Ka(k)=P(k/k-1)HTJ-1(k)

(40)

P(k/k)=[I-KaH]P(k/k-1)

(41)

(42)

(43)

(44)

求得

(45)

(46)

(47)

以上就是基于卡爾曼濾波算法的載荷識別具體過程。

繪制出基于卡爾曼濾波算法的動荷載識別流程圖如圖1所示。

圖1 基于卡爾曼濾波器的荷載識別原理Fig.1 Load identification principle based on the Kalman filter

2.2 云模型優化的離散卡爾曼濾波動荷載識別

在實際測量中,結構的動力測量響應會受到噪聲和誤差等問題的干擾,進而導致荷載識別結果產生較大偏差,因此僅僅利用該方法進行荷載識別難以處理噪聲和誤差帶來的干擾問題。因此,本文采用云模型理論來解決此類干擾問題。云的生成的具體算法稱為云發生器[14],可以分為正向云發生器和逆向云發生器,其中,逆向云發生器通過輸入一定數量的云滴求出云模型的3個數字特征。本文利用多次荷載識別結果生成云滴,借鑒逆向云發生器的原理產生云的3個數字特征,解決測量誤差及噪聲帶來的干擾問題。具體過程如下:

(1)根據每個時刻n次采樣數據的識別結果計算該組數據的樣本均值為

(48)

一階樣本絕對中心距

(49)

樣本方差

(50)

(51)

(52)

(53)

(3) 利用數字特征計算確定度μ

(54)

(4) 利用確定度確定識別值

由云模型性質可知,該點確定度越高,則越接近云的期望值,而期望值是云的中心位置反映在云圖橫坐標上的值,表示該點在論域中最能夠代表定性概念,所以本文計算每一時刻每次采樣數據的識別結果對應的確定度μ,此時求得n個μ,然后選取其中的最大值對應的識別值,表示該識別值最接近該時刻真實的荷載值,其余時刻以此類推。

通過圖2可以進一步理解如何利用確定度確定識別值。圖2中:縱軸為某時刻識別值的確定度大小;橫軸為該確定度對應的荷載大小,確定度最大的點對應的荷載值即為本次識別最終的識別值。

圖2 某時刻識別值的確定度Fig.2 Determination of the recognition value at some time

綜合本節內容,建立基于云模型優化的卡爾曼濾波的荷載識別處理策略如下:

(1)進行多次測量,獲取多次含測量噪聲的結構動力響應數據。

(2)將動力響應數據代入基于卡爾曼濾波的荷載識別算法程序,多次采樣數據得到多個識別結果,從而構造出包含多個云滴的云模型。

(3)計算所有云滴的數字特征,即期望、熵和超熵,并計算每個時刻識別值的確定度。

(4)根據每個時刻識別值確定度的大小確定最終識別值,確定度最大的點對應的荷載即為本次的荷載識別值。

3 數值仿真

3.1 有限元建模

為了驗證本文提出的基于云模型和卡爾曼濾波的動荷載識別方法的有效性,本文針對具體實際工程臨時看臺結構建立有限元模型,設置多種荷載類型和噪聲水平進行數值仿真。

某比賽賽場臨時看臺結構如圖3所示,看臺支撐結構為承插型盤扣式鋼骨架,立桿縱向距離為2.0 m,橫向距離為2.4 m,看臺首排高度為1 m,后排高度為8.0 m。看臺主體由承插型盤扣式鋼骨架、看臺板及其下部的三角形桁架、欄桿和通道體系組成,立桿底座采用槽鋼固定至樓板。本文僅針對承插型盤扣式鋼管骨架進行研究。

圖3 某比賽賽場臨時看臺結構[15](mm)Fig.3 Structure of temporary stand for a competition field (mm)

臨時看臺結構所用鋼材為Q345B,彈性模量為E=2.06×1011Pa,泊松比為0.3,密度7 850 kg/m3。桿件的截面屬性如表1所示。表1中:P和B分別代表圓鋼管和方鋼管。

表1 桿件截面屬性

由JGJ 231—2010《建筑施工承插型盤扣式鋼管支架安全技術規程》[16]可知承插型盤扣式鋼管骨架結構本質上是一種半剛性空間框架鋼結構,水平桿與立桿之間介于鉸接與剛接之間的一種連接形式,設計時可以偏于保守的將骨架結構考慮為鋼桁架形式,即可以將半剛性節點按鉸接建模。因此,將看臺骨架結構簡化為節點鉸接,底部固接的結構。為了簡化計算,本文僅選取部分看臺支撐結構作為研究對象。本文利用ABAQUS軟件中的T3D2H桁架單元創建臨時看臺支撐結構模型如圖4所示。

圖4 臨時看臺支撐結構簡化模型Fig.4 Temporary stand support structure simplified model

3.2 參數和工況設置

定義時程計算工況為:分析時長為10 s,分析時間步長0.001 s,輸出時間步長1 000步;定義荷載工況依次為正弦荷載、方波荷載、鋸齒荷載以及多個正弦疊加荷載作用在臨時看臺上表面,假設荷載作用方向豎直向下,并且每個節點均勻分擔荷載;測點布置如圖5所示,即本文數值模擬是通過該16個節點的豎向位移響應進行看臺動荷載識別,看臺節點豎向位移時程如圖6所示。

圖5 測點布置圖Fig.5 Measuring point layout

圖6 節點豎向位移時程圖Fig.6 Time history diagram of node vertical displacement

另外,為了驗證該方法在位移響應的魯棒性,將在觀測響應中混入高斯白噪聲,噪聲施加方式如下

(55)

3.3 荷載識別結果分析

為驗證本文提出的云模型優化卡爾曼濾波的荷載識別方法在荷載識別中具有更高的精度,選擇與基本離散卡爾曼濾波荷載識別方法進行對比。以下將從荷載時程和標準均方根誤差兩個角度評估荷載識別精度。由于篇幅限制,僅選取節點5的識別結果進行展示,如圖7所示。其余節點識別效果與節點5類似。

圖7 荷載識別結果時程圖Fig.7 Load identification result time course chart

為了評估荷載識別的精度,本文將采用標準均方根誤差(normalized root mean square error,NRMSE)來量化誤差,其表達式如下

(56)

由于噪聲的生成具有隨機性,為了更好地評價該方法的荷載識別性能,對每個噪聲工況獨立運行5次并求均值。該過程采用MATLAB編程實現,結果如表2所示。

表2 不同噪聲水平下荷載識別的NRMSE值

本文提出的云模型優化卡爾曼濾波荷載識別算法識別精度高,且該方法對計算硬件要求較低,計算速度與基本卡爾曼濾波方法接近。通過對結構位移響應信號施加不同水平的噪聲及其荷載識別效果圖可以看出該方法具有較好的穩定性及誤差處理能力。對看臺結構模型依次施加正弦、方波、鋸齒波以及正弦疊加荷載進行數值仿真,結果表明:正弦荷載識別效果最優,而后依次為正弦疊加荷載、鋸齒波和方波荷載。由圖6可知方波和鋸齒波荷載在曲線轉折處識別誤差較大,尤其是方波的初始點和轉折處誤差很大,導致NRMSE值過大,但平緩荷載處識別效果明顯較好,這是由于荷載突變處對識別算法的實時追蹤能力要求較高。

由表2的NRMSE值可知,本文所建議的云模型優化卡爾曼濾波方法的識別誤差普遍小于基本卡爾曼濾波的荷載識別算法。隨著噪聲水平的升高,本文方法的NRMSE值變化較小,仍舊保持良好的識別性能,說明該方法具有良好的抗噪性能。

4 結 論

本文提出了一種基于云模型和離散卡爾曼濾波的動荷載識別方法。為了驗證該方法對看臺結構荷載識別的有效性,本文采用了簡化的臨時看臺結構模型進行數值仿真,并與基本離散卡爾曼濾波的荷載識別方法進行了對比。通過不同類型荷載識別的對比結果可看出所建議方法具有良好的穩定性和較好的處理能力;另外,本文設置了2%和5%兩種噪聲水平,可以看出隨著噪聲水平的升高,基于卡爾曼濾波方法的識別結果出現明顯偏差,荷載時程曲線擬合程度降低,標準均方根誤差值也隨即增大,而本文提出的方法受噪聲影響較小,仍保持著良好的識別性能。因此,本文提出的基于云模型和離散卡爾曼濾波的動荷載識別方法識別精度較高,抗噪性良好,而且和基本卡爾曼濾波方法類似,具有一定的在線監測荷載的功能,故具有較好的工程應用前景。