基于振動監(jiān)測的斜拉索時變索力識別方法

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

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

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

Casas[1]、Hiroshi[2]各自采用不同的索力與基頻之間的關系，建立了基于基頻的索力識別方法；Ressel和Lardner在他們的索力識別中重點強調了模態(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]采用曲線擬合的方法建立了兩端固支索的索力與基頻之間的關系，并基于該關系提出了索力識別的方法。

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

1 理論推導

1.1 斜拉索的運動方程

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

圖1 斜拉索模型

圖2 微元體受力分析

(1)

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

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

(2)

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

(3)

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

(4)

(5)

邊界條件：

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

(6)

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

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

(7)

1.2 索運動方程的近似級數解

索的動撓度可以用有限級數表示為

(8)

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

(9)

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

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

(10)

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

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

(11)

式中，質量矩陣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子空間的小波尺度函數表示為

(17)

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

AP=Q

(18)

其中

An×(l+1)=

(19)

(20)

(21)

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

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

(22)

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

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

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

q0=Φ-1v0

(23)

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

(24)

(25)

(26)

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

(27)

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

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

(28)

(29)

其中

(30)

(31)

(32)

式中，Δt是采樣周期。

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

(33)

(34)

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

(35)

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

(36)

(37)

圖3 小波-擴展卡爾曼濾波識別數據流程

2 數值試驗

研究方法如下：首先，由ANSYS有限元程序對橋梁在車輛荷載下的索力響應進行模擬；然后，根據所得到的時變索力采用Newmark-β法求解方程(11)，得到索在環(huán)境激勵和軸向激勵共同作用下的響應；最后，選取部分索上的響應作為觀測，進行索力的識別。

2.1 背景資料

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

表1 N26索的基本信息

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

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

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

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

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

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

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

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

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

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

3 結 論

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

[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.