張福儉， 李 惠， 毛晨曦,2

(1.哈爾濱工業(yè)大學(xué) 土木工程學(xué)院， 黑龍江 哈爾濱 150090； 2.中國地震局工程力學(xué)所， 黑龍江 哈爾濱 150080)

斜拉索是斜拉橋上最重要的結(jié)構(gòu)構(gòu)件之一，索力的往復(fù)變化會造成索的疲勞損傷。索力的識別對于斜拉橋的維護與安全評定具有重要的意義，國內(nèi)外有很多關(guān)于索力識別的文獻報道?；谡駝臃ǖ乃髁ψR別方法由于其低成本和便捷性，在健康監(jiān)測實踐中廣泛采用。

Casas[1]、Hiroshi[2]各自采用不同的索力與基頻之間的關(guān)系，建立了基于基頻的索力識別方法；Ressel和Lardner在他們的索力識別中重點強調(diào)了模態(tài)階躍現(xiàn)象[3]；Ren[4]、Geier[5]采用了基于模態(tài)分析的索力識別方法；Park等[6]、Kim和Park[7]采用了基于靈敏度分析的方法進行了索力識別。Kim等[8]還對上述各種識別方法進行了對比研究，他們建議采用系統(tǒng)辨識方法和張緊弦理論方法。Ceballo和Preto[9]、Ren[10]、Fang和Wang[11]采用曲線擬合的方法建立了兩端固支索的索力與基頻之間的關(guān)系，并基于該關(guān)系提出了索力識別的方法。

國內(nèi)外與此相關(guān)的研究還有許多，但是以往振動法一般只能識別在特定時間內(nèi)的索力的平均值，無法識別時變的索力。因此，所得到的索力識別結(jié)果無法用于索的疲勞損傷評估，這是非常遺憾的。本文提出了一種結(jié)合擴展卡爾曼濾波和小波級數(shù)的時變索力識別方法。

1 理論推導(dǎo)

1.1 斜拉索的運動方程

如圖1所示，斜拉索的支座為A、B。B點是固定支座，A點是沿著AB方向(記為ξ方向)的滑動支座。索承受兩個方向的外力：垂直索長方向的(如風(fēng)荷載等)和由車輛引起的索支座移動。如前所述，車輛的通過會引起索力的變化。另外，采用如下假定：(1)斜拉索的自重遠(yuǎn)小于張力，拉索處于小垂度狀態(tài)，因此自重和抗彎剛度引起的位移可以忽略，因為它們都是二階效應(yīng)[12]；(2)由風(fēng)荷載引起的橫向位移獨立于由車輛荷載引起的軸向變形。這樣η方向的平衡方程如下(受力分析如圖2)：

圖1 斜拉索模型

圖2 微元體受力分析

(1)

式中，T(t)是拉索張力，它隨著車輛荷載的通過是時變的；F(ξ,t)是垂直于索長方向的外部激勵(這里主要考慮橫向風(fēng)荷載)；v是沿y方向的橫向撓度；m是單位長度質(zhì)量；c是單位長度阻尼系數(shù)；s代表微元長度。

沿ξ方向的微元平衡方程為

(2)

方程(2)可以直接積分得到

(3)

式中H(t)代表索力的水平分量，可認(rèn)為沿索長不變。將方程(3)帶入(1)可得

(4)

(5)

邊界條件：

v(0,t)=0,v(L-u(t),t)=0

(6)

式中，L是索的長度；u(t)代表索端沿ξ方向的運動。.考慮到沿索長方向的索端移動遠(yuǎn)小于索長，可忽略索的軸向慣性力。邊界條件(6)可以簡化為

v(0,t)=0,v(L,t)=0

(7)

1.2 索運動方程的近似級數(shù)解

索的動撓度可以用有限級數(shù)表示為

(8)

式中，qj(t)是第j階模態(tài)坐標(biāo)；φj(ξ)是滿足幾何邊界條件的連續(xù)的形函數(shù)，并且φj(0,t)=0,φj(L,t)=0滿足正交條件[13]：

(9)

將方程(8)帶入(5)，可以得到

H(t)qj(t)φj″(ξ)]=F(ξ,t)

(10)

式中，r是待研究的模態(tài)階數(shù)。

在方程(10)兩側(cè)同時乘以φi(ξ)并沿索長積分即可得到矩陣形式下的模態(tài)域振動方程為

(11)

式中，質(zhì)量矩陣M=[mij]； 阻尼矩陣C=[cij]； 剛度矩陣K=[kij]；模態(tài)激勵Fq=[Fq1Fq2…Fqr]T，各矩陣元素可由如下方程確定：

(12)

(13)

(14)

(15)

1.3 索力的識別

取任一階模態(tài)的振動方程

i=1,2,…,r

(16)

其中的時變剛度ki可以用Vj子空間的小波尺度函數(shù)表示為

(17)

假定阻尼已知，這樣方程(16)可以表示為

AP=Q

(18)

其中

An×(l+1)=

(19)

(20)

(21)

矩陣A可以由小波尺度函數(shù)及模態(tài)形函數(shù)的離散點采樣值得到。垂直索力方向的外荷載可以由風(fēng)速儀測量得到，這樣Q也可以得到。由于廣義坐標(biāo)是不能直接測量得到的，它們需要由測量到的位移、速度或加速度變換得到。

如果傳感器數(shù)目與研究模態(tài)的階數(shù)相同，那么廣義坐標(biāo)可以由如下方程得到：

(22)

1.4 由有限觀測進行索力識別

通常在索上只能安裝有限數(shù)量的傳感器。因此傳感器的數(shù)量是小于控制模態(tài)的數(shù)量的，這時可以采用如下方法確定廣義坐標(biāo)。

當(dāng)只有有限觀測時，如{v(ξ1,t)，…，v(ξi,t)}T，而{v(ξi+1,t)，…，v(ξr,t)}T沒有被觀測到。我們采用擴展卡爾曼濾波來根據(jù)已觀測的部分對未觀測的部分進行估計。令v0={v(ξ1,t)，…，v(ξi,t),v(ξi+1,t)，…，v(ξr,t)}T，則根據(jù)公式(22)

q0=Φ-1v0

(23)

將公式(23)帶入(11)，得到如下在物理坐標(biāo)下的運動方程：

(24)

(25)

(26)

p(v*)～N(0,R)

(27)

式中，R為矩陣v*(t)的協(xié)方差矩陣。

將方程(25)和方程(26)變換為離散形式[14]

(28)

(29)

其中

(30)

(31)

(32)

式中，Δt是采樣周期。

(2)根據(jù)過程方程，有對tn+1時的先驗估計

(33)

(34)

(3) 計算卡爾曼濾波增益κn+1

(35)

(4)得到狀態(tài)變量和方差矩陣的最優(yōu)估計為

(36)

(37)

圖3 小波-擴展卡爾曼濾波識別數(shù)據(jù)流程

2 數(shù)值試驗

研究方法如下：首先，由ANSYS有限元程序?qū)蛄涸谲囕v荷載下的索力響應(yīng)進行模擬；然后，根據(jù)所得到的時變索力采用Newmark-β法求解方程(11)，得到索在環(huán)境激勵和軸向激勵共同作用下的響應(yīng)；最后，選取部分索上的響應(yīng)作為觀測，進行索力的識別。

2.1 背景資料

本文選取濱州黃河公路大橋上最長的索N26作為研究對象，其基本信息如表1所示。

表1 N26索的基本信息

軸向激勵，由ANSYS模擬生成，模擬工況為50 t卡車以20 m/s的速度通過大橋，結(jié)果如圖4(a)。索的橫向荷載采用人工風(fēng)，模擬風(fēng)速為20 m/s如圖4(b)所示。這樣，模態(tài)位移、模態(tài)速度、模態(tài)加速度都可以通過求解方程(11)得到。

(a) 由自重和車載引起的索應(yīng)力 (b) 人工風(fēng)時程圖4 索的軸向及橫向激勵

2.2 通過對一個模態(tài)所有狀態(tài)變量觀測而得到的索力識別結(jié)果

這里采用Db16小波母函數(shù)進行索力識別。識別結(jié)果(如圖5)顯示此方法具有很高的精度，但是存在比較明顯的邊端效應(yīng)。

(a) 識別到的索力變化系數(shù) (b) 放大的索力變化系數(shù) 圖5 全觀測下索力識別結(jié)果

2.3 通過有限觀測得到的索力識別結(jié)果

沿索長觀測所有的位移、速度和加速度是不現(xiàn)實的，因此，這里采用擴展卡爾曼濾波來進行有限觀測條件下的索力識別。在下面的識別中，只觀測到索上兩點的位移。

(a) 軸向輸入 (b) 橫向人工風(fēng)輸入圖6 有限觀測下工況的軸向和橫向輸入

(a) 識別的剛度變化系數(shù) (b) 放大的剛度變化系數(shù)圖7 有限觀測下的剛度變化識別結(jié)果

為了驗證所提出方法的有效性，這里采用了一個稍微復(fù)雜一些的工況：模擬兩輛50 t重車先后分別過橋，第一輛車的速度為20 m/s，第二輛車的速度為30 m/s。這樣的工況下(圖6)，識別結(jié)果如圖7所示。

3 結(jié) 論

本文提出了一種基于小波級數(shù)結(jié)合卡爾曼濾波的時變索力識別方法。以濱州黃河公路大橋為背景，進行了數(shù)值模擬研究，結(jié)果表明：(1)在觀測全局信息時，本文所提出的方法對于較小的索力變化亦有較高的精度；(2)在有限觀測條件下本文方法依然可行，但識別結(jié)果精度降低。

[1] Casas J R. A combined method for measuring cable forces: the cable-stayed Alamillo bridge, Spain [J]. Structural Engineering International, 1994, 4(4) : 235-240.

[2] Zui H, Shinke T, Namita Y. Practical formulas for estimation of cable tension by vibration method [J]. Journal of Structural Engineering, 1996, 122(6): 651-656.

[3] Russell J C, Lardner T J. Experimental determination of frequencies and tension for elastic cables [J]. Journal of Engineering Mechanics, 1998, 124(10): 1067-1072.

[4] Ren W X, Chen G , Hu W H. Empirical formulas to estimate cable tension by cable fundamental frequency [J]. Structural Engineering and Mechanics, 2005, 20(3): 363-380.

[5] Geier R, De Roeck G , Flesch R. Accurate cable force determination using ambient vibration measurements [J]. Structure and Infrastructure Engineering, 2006, 2(1): 43-52.

[6] Park S, Choi S, Oh S T,et al. Identification of the tensile force in high-tension bars using modal sensitivities [J]. International Journal of Solids and Structures, 2006, 43(10): 3185-3196.

[7] Kim B H,Park T. Estimation of cable tension force using the frequency-based system identification method [J]. Journal of Sound and Vibration, 2007,304(3-5): 660-676.

[8] Kim B H, Park T, Shin H,et al. A comparative study of the tension estimation methods for cable supported bridges [J]. International Journal of Steel Structures, 2007,7(1): 77-84.

[9] Ceballos M A , Prato C A. Determination of the axial force on stay cables accounting for their bending stiffness and rotational end restraints by free vibration tests [J]. Journal of Sound and Vibration, 2008; 317(1-2): 127-141.

[10] Ren W X, Liu H L , Chen G. Determination of cable tensions based on frequency differences [J]. Engineering Computations, 2008, 25(2): 172-189.

[11] Fang Z ,Wang J. Practical formula for cable tension estimation by vibration method [J]. Journal of Bridge Engineering, 2012,17: 161-164.

[12] Irvine H M. Cable Structures [M]. Cambridge : The MIT Press, 1981.

[13] Li R H. Galerkin Finite Element Method for Boundary Value Problems (in Chinese) [M].Beijing: Science Press, 2005.

[14] Franklin G F, Workman M L , Powell D. Digital Control of Dynamic Systems[M]. Boston: Addison-Wesley Longman Publishing Co Inc, 1997.