粒子群算法在基于频率的两端固结吊杆索力识别中的应用

李 睿， 李晓章， 郑祥隆， 周亦唐

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