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

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

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

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

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

1 理論推導(dǎo)

1.1 斜拉索的運(yùn)動(dòng)方程

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

圖1 斜拉索模型

圖2 微元體受力分析

(1)

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

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

(2)

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

(3)

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

(4)

(5)

邊界條件：

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

(6)

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

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

(7)

1.2 索運(yùn)動(dòng)方程的近似級(jí)數(shù)解

索的動(dòng)撓度可以用有限級(jí)數(shù)表示為

(8)

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

(9)

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

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

(10)

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

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

(11)

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

(12)

(13)

(14)

(15)

1.3 索力的識(shí)別

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

i=1,2,…,r

(16)

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

(17)

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

AP=Q

(18)

其中

An×(l+1)=

(19)

(20)

(21)

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

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

(22)

1.4 由有限觀(guān)測(cè)進(jìn)行索力識(shí)別

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

當(dāng)只有有限觀(guān)測(cè)時(shí)，如{v(ξ1,t)，…，v(ξi,t)}T，而{v(ξi+1,t)，…，v(ξr,t)}T沒(méi)有被觀(guān)測(cè)到。我們采用擴(kuò)展卡爾曼濾波來(lái)根據(jù)已觀(guān)測(cè)的部分對(duì)未觀(guān)測(cè)的部分進(jìn)行估計(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)下的運(yùn)動(dòng)方程：

(24)

(25)

(26)

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

(27)

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

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

(28)

(29)

其中

(30)

(31)

(32)

式中，Δt是采樣周期。

(2)根據(jù)過(guò)程方程，有對(duì)tn+1時(shí)的先驗(yàn)估計(jì)

(33)

(34)

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

(35)

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

(36)

(37)

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

2 數(shù)值試驗(yàn)

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

2.1 背景資料

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

表1 N26索的基本信息

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

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

2.2 通過(guò)對(duì)一個(gè)模態(tài)所有狀態(tài)變量觀(guān)測(cè)而得到的索力識(shí)別結(jié)果

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

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

2.3 通過(guò)有限觀(guān)測(cè)得到的索力識(shí)別結(jié)果

沿索長(zhǎng)觀(guān)測(cè)所有的位移、速度和加速度是不現(xiàn)實(shí)的，因此，這里采用擴(kuò)展卡爾曼濾波來(lái)進(jìn)行有限觀(guān)測(cè)條件下的索力識(shí)別。在下面的識(shí)別中，只觀(guān)測(cè)到索上兩點(diǎn)的位移。

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

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

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

3 結(jié) 論

本文提出了一種基于小波級(jí)數(shù)結(jié)合卡爾曼濾波的時(shí)變索力識(shí)別方法。以濱州黃河公路大橋?yàn)楸尘?，進(jìn)行了數(shù)值模擬研究，結(jié)果表明：(1)在觀(guān)測(cè)全局信息時(shí)，本文所提出的方法對(duì)于較小的索力變化亦有較高的精度；(2)在有限觀(guān)測(cè)條件下本文方法依然可行，但識(shí)別結(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.