粒子群算法在基于頻率的兩端固結吊桿索力識別中的應用

李 睿， 李曉章， 鄭祥隆， 周亦唐

(1. 昆明理工大學 建筑工程學院， 昆明 650504； 2. 浙江大學 建筑工程學院， 杭州 310058)

吊桿是中下承式拱橋的重要傳力構件，吊桿的索力測試是在拱橋施工以及運營過程中都必不可少的一項工作。索力測試方法包括油壓表法、傳感器法、磁通量法、頻率法等。其中，頻率法以簡便、經濟、實用等優點成為普遍使用的測試方法。

吊桿的頻率和索力的關系可用頻率方程表示，頻率法求解索力實際就是求解吊桿頻率方程的逆問題。傳統的頻率法忽略吊桿的彎曲剛度和端部約束作用，假定吊桿是一根張緊的弦，這時根據吊桿頻率方程可以建立吊桿基頻與索力間的顯式關系

T=4ml2f21

(1)

式中：T,m,l,f1分別為吊桿索力、線密度、長度、基頻。對于較細長拉索而言，由于其彎曲剛度和端部約束對頻率的影響較小，采用弦模型式(1)計算索力能滿足工程要求。但對于長度相對較短并且兩端通常錨固于拱肋和橋面混凝內的拱橋吊桿，其彎曲剛度以及兩端固結約束對頻率的影響已不可忽略。此時，對兩端固結吊桿的頻率方程是一個超越方程，無法建立索力與頻率的顯式關系。

不少學者對頻率法進行研究與探索，提出若干實用方法。Zui等[1]針對斜拉索提出索力求解的分段擬合實用公式；Ren等[2]考慮了索垂度和彎曲剛度對斜拉索頻率的影響，采用能量法和曲線擬合方法，建立了利用基頻計算索力的實用公式；方志等[3]基于兩端固結梁在軸向拉力作用的橫向振動方程，通過擬合得出適用于拉索和吊桿張力測試計算公式。李新生等[4-5]基于能量原理，將兩階以上近似振型函數代入Rayleigh-Ritz公式，得到聯立方程組然后消去彎曲剛度求解索力。他們的不同在于吊桿邊界條件的假定以及近似振型函數的選擇。上述學者都采用不同的策略避免對復雜超越方程的直接求解以實現索力計算簡便性，但是計算簡便性往往以犧牲計算精度為代價。

此外，還有不少學者則通過采用數值優化方法求解超越方程實現索力識別。Sagüés等[6]通過試算法進行體外預應力索索力求解； Lagomarsino等[7]利用牛頓梯度法對誤差函數進行最小化優化，用三階頻率進行金屬拉桿張力、彎曲剛度、端部轉動剛度識別；Kim等[8]利用基于參數靈敏度修正的迭代算法識別斜拉索的索力、彎曲剛度和軸向剛度；Ceballos等[9]也提出用多階頻率進行斜拉索索力的迭代識別的方法；Xie等[10]利用有限元模型對吊桿進行模擬，并通過特征值分析建立吊桿張力、彎曲剛度和邊界條件與頻率的超越方程，然后采用自適應選擇的遺傳算法作進行吊桿張力、彎曲剛度以及邊界條件的識別。采用數值方法求解超越的頻率方程，能夠在計算機輔助下實現吊桿索力精確識別。但直接數值求解法常見于國外研究；國內方面，除Xie等利用遺傳算法進行吊桿多參數識別外，尚未見到其他相關研究。

另一方面，粒子群算法是由Kennedy等[11]于1995年提出的一種數值優化算法。它與遺傳算法一樣具有優秀的全局搜索能力，并且粒子群算法還具有思路簡潔、容易編程的特點。目前已有不少學者應用粒子群算法來解決復雜工程問題。Elegbede[12]最早將粒子群算法用于求解結構可靠度指標；賈義鵬等[13]利用粒子群算法進行巖爆預測；陳志軍等[14]基于粒子群算法對斜拉橋成橋索力進行優化。

本文針對中下承式拱橋中比較常見的兩端固結吊桿，利用粒子群算法對超越的頻率方程進行多維空間搜索求解，實現了利用兩階實測頻率進行吊桿索力和彎曲剛度的精確識別。和Xie等的方法相比，本文基于兩端固結吊桿的頻率解析方程進行求解，避免對吊桿動力特性的有限元建模分析，使索力求解方法更易于編程實現。

1 兩端固結吊桿的頻率方程

圖1為中、下承式拱橋中常見的吊桿形式，吊桿上下兩端分別固結于拱肋和橋面混凝土內。圖中E為吊桿的彈性模量，I為截面慣性矩，l為長度，m為單位長度質量，T為吊桿索力。

不考慮吊桿大變形和二次力影響，將吊桿當作一根受常軸力作用的均勻截面梁進行分析，其自由振動方程為

圖1 兩端固結吊桿分析模型

(2)

式中：v(x,t)為吊桿各點在t時刻的橫向位置；m、EI和T分別為吊桿的線密度，彎曲剛度和張力。式(2)為線性微分方程，其解具有

v(x,t)=φ(x)·Y(t)

(3)

將式(3)代入式(2)，分離變量之后，可得到兩個獨立方程

EIφIV(x)-Tφ″(x)-ω2φ(x)=0

(4a)

(4b)

式中:ω為吊桿的自振圓頻率。式(4a)的通解為

φ(x)=A1cos(δx)+A2sin(δx)+

A3sinh(εx)+A4cosh(εx)

(5)

式中：A1、A2、A3、A4為任意實常數。對兩端固結吊桿，已知的邊界條件為

φ(x)|x=0=0;φ′(x)|x=0=0

φ(x)|x=l=0;φ′(x)|x=l=0

(6)

將邊界條件式(6)代入式(5)得到關于常數A1、A2、A3、A4的線性方程組，然后根據方程組具有非零解的充要條件，可獲得兩端固結吊桿的頻率方程

1-cos(δl)cosh(εl)+

(7)

式中：參數δ、ε的計算見式(5)。

頻率方程式(7)建立了吊桿索力和頻率的關系。同時，應注意到吊桿是由多根平行或半平行鋼絲并聯組成的構件，其彎曲剛度EI受鋼絲間的黏結程度、吊桿的索力等的影響，不能準確預知。因此，在利用頻率求解索力時，吊桿的彎曲剛度EI也應該和索力T一樣當做未知量來處理。此時，頻率方程式(7)為二元超越方程，無法通過代入兩階頻率實現索力和彎曲剛度的直接求解。本文采用粒子群算法實現求解過程。

2 利用粒子群算法進行吊桿索力數值求解

粒子群算法(Particle Swarm Optimization, PSO)通過模擬鳥群的捕食行為來達到優化問題的求解。首先在解空間內隨機初始化鳥群粒子，這些“粒子”在解空間內以某種規律移動，經過若干次迭代后找到最優解。在迭代過程中，粒子通過跟蹤兩個“優粒子”來更新自己：一個是粒子個體歷史最優解Pid，另一個是整個粒子群歷史最優解Pgd。迭代終止判據為最新迭代步的粒子群最優解滿足精度要求或者迭代數超過限值。

利用粒子群算法求解吊桿索力的迭代求解過程通過自編程序實現，圖2為程序流程框圖。對整個求解流程，本文進行幾點詳細說明。

2.1 種群粒子{T,EI}T的解空間

吊桿張力T和彎曲剛度EI可以預先設定一個取值范圍，在這個取值范圍內進行搜索求解能有效提高優化效率。本文對吊桿張力和彎曲剛度的取值范圍參考Xie等的研究，設定如下

(8)

式中:EI0為按吊桿鋼絲截面為剛性計算得到的彎曲剛度，上限系數取1.2則是考慮鋼絲外層PE套管影響；T0為按式(1)計算得到的張力值。在隨機生成初始種群時，種群粒子按式(8)給出的取值范圍進行取值；在迭代過程中，某粒子的位置超出取值空間時，將超越邊界點位置作為當前粒子的更新位置進行處理。

圖2 粒子群算法求解索力程序框圖

2.2 粒子適應度

粒子群算法作為一種優化算法，需要設定優化目標函數，目標函數值的大小體現粒子適應程度的優劣，所以也稱為適應度。本文選取識別頻率與實測頻率差異最小的粒子作為最優解。因此，設定的目標函數為頻率相對誤差平方和開方，即

(9)

2.3 粒子速度和位置的更新

在每一個迭代步中，粒子位置通過跟蹤兩個“優粒子”來更新自己。為保證識別方法的全局收斂，粒子速度和位置采用如下公式進行迭代更新[16]

Vid=K×[Vid+c1R1(Pid-xid)+c2R2(Pgd-xid)]

xid=xid+Vid

(10)

式中:下標i為粒子編號；下標d為解空間的維數，對于本文求解問題，d=1代表張力T，d=2代表彎曲剛度EI；Vid為粒子速度，xid為當前粒子位置；pid和pgd分別表示粒子個體歷史最優解和粒子群歷史最優解；R1和R2為兩個在[0,1]范圍內取值的隨機數；c1和c2稱為學習因子，取正常數；K為收斂因子，其計算為

(11)

式中：φ=c1+c2，且φ>4。根據Eberhart等的研究結果[16]，取學習因子c1=c2=2.05，所以φ=4.1，K=0.729。

通常粒子速度需要限定最大值Vmax。本文識別過程中，當粒子更新位置超出預設{T,EI}T的取值空間時，以超越邊界點位置作為當前粒子更新位置，因而未設置Vmax。

2.4 種群粒子多樣性

為了提高識別算法的收斂速度，本文借鑒了遺傳算法中的“變異”思想，以保持種群粒子多樣性，確保粒子迅速向全局最優解靠攏。具體操作是在種群迭代過程中，根據各個粒子適應度排名按比例區分優等粒子和劣等粒子，優等粒子保留到下一迭代步，劣等粒子則直接淘汰，并由隨機生成的新粒子替代，即完成劣等粒子變異。

3 方法驗證和誤差分析

為驗證方法的可行性，本文選取具有代表性的10根虛擬吊桿進行索力識別。虛擬吊桿的理論計算參數如表1。表中的吊桿截面參數參照實際吊桿取值，吊桿截面由139根Ф7鋼絲組成，按鋼絲完全黏結考慮時其剛性截面彎曲剛度為EI0=656.0 kN·m。表中2列給出其線密度m=45.7 kg/m；表中第3列為吊桿的長度，3～15 m為常見的拱橋吊桿長度區間。表中第4列、第5列為吊桿索力和彎曲剛度的理論設定值，上標“*”表示該值為真值，與之后的識別值加以區分；第6列、第7列為吊桿前二階頻率值，為驗證方法本身的識別精度，表1中近似以精確至小數點后五位的頻率作為精確頻率，減少頻率舍入誤差帶來的影響。吊桿前二階頻率值通過將吊桿物理參數代入式(7)，然后分區間利用弦割法進行求解獲得。

現假設表1中吊桿的索力和彎曲剛度未知，以第6列、第7列中的精確頻率及其他已知參數作為識別條件進行識別。設定種群內個體規模數量為100，種群優等粒子占比為0.2。

表1 虛擬吊桿參數

圖3給出粒子群歷史最優解pgd隨迭代數的變化過程。橫坐標表示識別計算迭代步，圖3(a)中豎坐標T/T*為張力的識別值與理論設定值的比值，圖3(b)中豎坐標EI/EI*為彎曲剛度的識別值與理論設定值的比值?？梢钥闯?，所有吊桿的識別索力和彎曲剛度均能在100個迭代步以內迅速收斂。

(a) 索力

(b) 彎曲剛度

表2給出了第100迭代步時十根虛擬吊桿的索力識別結果。表中第5、第6列為利用精確頻率進行識別的結果，可以看到索力識別精度可達0.01%，幾乎沒有誤差。這驗證了本方法的準確性。

此外，表中第2、第3列還給出了利用弦模型式(1)計算得到的索力值以及誤差，可以看到索力計算誤差最小也達到15%以上。并且吊桿越短，采用弦模型公式計算的索力誤差越大。這是因為吊桿越短，其彎曲剛度對頻率的影響越大，弦模型公式因忽略吊桿彎曲剛度產生的誤差也就越大。這表明弦模型公式不適用于兩端固結吊桿尤其是兩端固結短吊桿的索力計算。

考慮到現場實測時頻率不可能精確至小數點后很多位，現將表1中頻率四舍五入保留兩位小數，然后再進行索力識別，粒子群算法參數設定與前文一致。表2中第7、第8列分別給出利用舍入頻率進行索力識別的結果，可以看到最大誤差達到0.612%。與采用精確頻率進行識別的結果進行比較，采用舍入頻率的進行識別的精度有所降低，這表明索力識別誤差主要由頻率誤差導致。

表2 虛擬吊桿索力識別結果

為分析頻率誤差對索力識別誤差的影響，表2中第2列給出了頻率誤差指標Ef，Ef計算方法為

(12)

式中：Ef為頻率誤差指標；ft,j、fe,j分別為精準頻率和舍入頻率；下標j為頻率階數。

圖4給出10根虛擬吊桿的頻率誤差指標Ef和索力識別誤差的散點關系圖，這里索力識別誤差用絕對值表示。從圖中可以看到，索力識別誤差隨頻率誤差增加大致呈遞增關系。這進一步表明當采用舍入頻率進行識別時，索力識別誤差主要是由頻率誤差引起的。

圖4 頻率誤差與索力識別誤差的關系

4 實例應用

為進一步驗證本文方法的工程實用性，以文獻[10]中給出實測吊桿數據進行索力識別，文獻中給出的索力準確值為1 635.5kN。實測吊桿線密度m=45.7 kg/m，長度l=9.00 m，吊桿實測前2階頻率為f1=12.286 Hz、f1=26.004 Hz。粒子群算法參數設定與前文一致。圖5給出實測吊桿的索力識別過程。圖中，實線表示索力準確值；虛線表示弦模型公式計算的索力值，散點實線表示索力識別值?？梢钥吹?，弦模型公式計算結果與索力準確值相差很大，弦模型公式索力計算值為2 235.0 kN，誤差可達36.66%；索力識別值可在100個迭代步內收斂至索力準確值附近，第100迭代步的索力識別值為1 615.0kN，誤差為-1.25%。實測吊桿的索力識別誤差要比虛擬吊桿的誤差要大。這需要考慮到在對實測吊桿進行索力識別時，吊桿參數取值還存在誤差，包括吊桿的線密度取值誤差，長度取值誤差以及上述的頻率誤差。這些誤差的累積導致了索力識別誤差的增加。盡管如此，實測吊桿索力識別誤差仍然小于2%，這表明本文方法能夠滿足工程精度要求。

圖5 實測吊桿索力識別收斂過程

5 結 論

(1) 采用工程中常用的弦模型公式計算拱橋兩端固結吊桿的索力誤差較大，難以滿足工程要求。并且吊桿越短，索力計算誤差越大，應盡量避免使用。

(2) 利用粒子群算法進行頻率超越方程逆向求解思路簡潔，易于編程實現，本文基于此建立了兩端固結吊桿的索力識別方法。通過對虛擬吊桿的索力識別以及誤差分析，驗證了本文方法能夠實現對吊桿索力的迅速準確識別；通過對實測吊桿的索力識別應用，表明本文方法滿足工程精度要求，能夠用于工程實測。

(3) 利用本文方法進行吊桿索力識別時，方法本身不引起誤差。索力識別誤差主要是由于吊桿參數取值誤差的累積造成的。因此在進行實測吊桿索力識別時，首先要確保吊桿線密度、吊桿長度以及吊桿頻率等參數的實測精度，以確保索力的準確識別。

參 考 文 獻

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

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

[3] 方志, 汪建群, 顏江平. 基于頻率法的拉索及吊桿張力測試[J]. 振動與沖擊, 2007, 26(9): 78-82.

FANG Zhi, WANG Jianqun, YAN Jiangping. Tension test of cable and hanger based on vibration method[J]. Journal of Vibration and Shock, 2007, 26(9): 78-82.

[4] 李新生, 項貽強. 基于撓度曲線振型函數的系桿拱橋柔性吊桿索力測量公式[J]. 工程力學, 2010, 27(8): 174-178.

LI Xinsheng, XIANG Yiqiang. Tension measurement formula of flexible hanger rods in tied-rods arch bridges based on vibration shape function of deflection curve[J]. Engineering Mechanics, 2010, 27(8): 174-178.

[5] 劉文峰，應懷樵，柳春圖. 考慮剛度及邊界條件的索力精確求解[J]. 振動與沖擊, 2003, 22(4): 12-14.

LIU Wenfeng, YING Huaiqiao, LIU Chuntu. Precise solution of cable tensile force in consideration of cable stiffness and its boundary conditions[J]. Journal of Vibration and Shock, 2003, 22(4): 12-14.

[6] SAGüéS A A, KRANC S C, EASON T G. Vibrational tension measurement of external tendons in segmental posttensioned bridges[J]. Journal of Bridge Engineering,

2006, 11(5): 582-589.

[7] LAGOMARSINO S, CALDERINI C. The dynamical identi?cation of the tensile force in ancient tie-rods[J]. Engineering Structures, 2005, 27(7): 846-856.

[8] 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/4/5): 660-676.

[9] CEBALLOS M A, PRATO C A. Determination of the axial force on stay cables accounting for their bending Stiffness and rotational end Restraints[J]. Journal of Sound and Vibration, 2008, 317(1): 127-141.

[10] XIE X, LI X Z. Genetic algorithm-based tension identi-fication of hanger by solving inverse eigenvalue problem[J]. Inverse Problems in Science & Engineering, 2014, 22(6):966-987.

[11] KENNEDY J, EBERHART R. Particle swarm optimization[C]// Proceedings of International Conference on Neural Networks. Perth: IEEE, 1995: 1942-1948.

[12] ELEGBEDE C. Structural reliability assessment based on particles swarm optimization[J]. Structural Safety, 2005, 27(2): 171-186.

[13] 賈義鵬, 呂慶, 尚岳全. 基于粒子群算法和廣義回歸神經網絡的巖爆預測[J]. 巖石力學與工程學報, 2013, 32(2): 343-348.

JIA Yipeng, Lü Qing, SHANG Yuequan. Rockburst prediction using particle swarm optimization algorithm and general regression neural network[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(2): 343-348.

[14] 陳志軍, 劉洋, 楊立飛,等. 基于粒子群優化算法的獨塔斜拉橋成橋索力優化[J]. 橋梁建設, 2016, 46(3): 40-44.

CHEN Zhijun, LIU Yang, YANG Lifei, et al. Optimization of stay cable tension of completed bridge of single-pylon cable-stayed bridge based on particle swarm optimization algorithm[J]. Bridge Construction, 2016, 46(3): 40-44.

[15] 丁麗娟, 程杞元. 數值計算方法[M]. 北京: 北京理工大學出版社, 2005.

[16] EBERHART R C, SHI Y. Comparing inertia weights and constriction factors in particle swarm optimization[C]// Proceedings of Congress on Evolutionary Computation. San Diego: IEEE, 2000: 84-88.