0 引言

频率法是目前索力测试最常用方法。对于单跨拉索，已有大量学者提出了基于频率的拉索索力实用计算公式^[1-6]，精度满足工程需要。出于减振目的，拉索一般带有一个或多个中间弹性支承，如拉索两端设置的减振阻尼器、斜拉索或吊杆间设置的抗风用辅助索、拱桥系杆的支承架或转向块等。这类索结构的振动频率与索力对应关系不甚明确，难以考虑中间支承的影响。王修勇等^[7]研究了单个拉索减振器对拉索振动频率的影响，提出了视减振器为刚性情形下采用高阶频率进行索力测试的策略。柯红军等^[8]根据拉索减振器安装前后频率的变化关系修正拉索的计算长度，并进行了实例验证。王荣辉等^[9]基于能量变分原理，采用动态精确梁单元，提出了识别具有多个弹性支承的拱桥系杆索力的有限元法。李平杰等^[10]将弹性支承反力看作是作用于梁上的强迫力，构造了带中间弹性支承及轴向力的梁模型，得到其横向振动频率解析算法。

动力刚度法是一种精确算法，本文采用该法求解带多个中间支承拉索的频率问题。将拉索视为无限自由度体系，导出精确的单元动力刚度矩阵(Dynamic Stiffness Matrix)，通过集成及求解得到整体动力刚度矩阵和频率方程，根据Wittrick-Williams算法^[11]采用二分法得到拉索自振频率。已有研究表明，基于拉索测试频率拟合的参数识别方法可得到较理想的索力测试精度^[12-13]。粒子群算法^[14]是求解参数识别问题的有效工具，本文提出了基于动力刚度法与粒子群算法的拉索参数识别，并通过仿真分析和一座斜拉桥拉索索力测试实例进行验证。

1 动力刚度矩阵的建立 1.1 单元动力刚度矩阵

单元动力刚度即为引起单元单位振幅所需要的动态力。带有多个中间支承拉索模型如图 1所示，不计拉索垂度，拉索抗弯刚度为EI，单位长度质量为m，两端张力为T，拉索总长为L；拉索两端边界简化为固结和铰接两种情况；拉索中间有n-1个刚度为K_i(i=1, 2, …, n-1) 的线弹性支承，将拉索分为n跨，每跨长度为L_i。拉索第i跨横向自由振动微分方程为：

(1)

其中v(x_i, t)为横向动位移，根据分离变量法，有

(2)

式中ω为拉索的固有振动圆频率。

将式(2) 代入式(1)，得

(3)

式(3) 的一般解为

(4)

式中，A=[A₁ A₂ A₃ A₄]^T为待定系数；

如图 2 (b)所示，单元杆端位移向量v=[v_A θ_AL v_B θ_BL]^T=[v(0) v′(0) v(L_i) v′(L_i)]^T，引入边界条件，式(4) 可写为矩阵形式：

(5)

其中

(6)

式中，s=sin(αu_i)；c=cos(αu_i)；S=sinh(βu_i)；C=cosh(βu_i)；u_i=L_i/L。

若单元两端固定，则由式(5) 可得

(7)

式(7) 有非零解则

(8)

展开并化简后得频率方程

(9)

满足式(9) 的频率称为固端频率，记为ω_F。当ω≠ω_F时，由式(5) 可得到系数矩阵，再进一步推导单元动力刚度矩阵。

式(5) 给出了用式(4) 中系数表示节点位移的表达式。由图 2(a)可知，单元杆端力向量可表示为

(10)

展开并化简后可得

(11)

式中

(12)

式中。

由式(5) 和(12) 可得

(13)

由式(13) 可得单元i的动力刚度矩阵：

(14)

式中,

1.2 总体动力刚度矩阵

根据式(14)，当弹性支承位置作为单元划分的节点时，就可以按照传统的结构力学方法将单元动力刚度矩阵集成为总体动力刚度矩阵，并对号入座将弹性支撑刚度叠加在总体刚度矩阵中相应节点竖向位移所对应的自由度主对角线位置上，得到拉索总体动力刚度矩阵如式(15) 所示。

(15)

式中。

对单元i的边界条件，引入取值区间为[0, 1]的无量纲参数^[15]：

式中，当k_i-1, k_i=0时，表示单元两端无竖向支承；当k_i-1, k_i=1时，表示单元两端竖向固结。

2 频率计算 2.1 Wittrick-Williams算法

Wittrick-Williams算法^[11](以下简称W-W算法)是结构自由振动问题精确分析计算的核心算法。这一算法并不直接计算频率，它是一种计数方法，它给出了结构低于某个给定值ω^*的频率数目的计算公式。

在结构的所有频率中，低于某个给定值ω^*的频率个数J由式(16) 给出：

(16)

从实施上来看，W-W算法中J_K的计算(常规Gauss消元)并无任何困难，而J₀的计算对于索结构来说则要困难得多。但是，可采用迂回的办法来避免J₀的计算。假设欲求的某一阶频率在区间(ω_l, ω_u)里，该方法就是将杆件划分成若干个单元，使得每个单元的第1阶固端频率都高于所求频率的上界ω_u，这样在利用W-W算法时就有J₀=0，避免了J₀的计算。

对于等截面拉索而言，索力T和抗弯刚度EI均为常数，单元长度L决定了每个单元固端频率的大小，也就是说每个频率上界ω_u都对应一个最大的单元长度L_max。

Hiroshi Zui^[2]提出了拉索固端基频的近似公式：

(17)

式中。

根据式(17) 可以反算出已知频率时拉索的长度，从而确定最大单元长度L_max，然后与所求频率的阶次(单元划分数目需大于所求频率阶次)相结合对拉索进行单元划分，便可直接取J₀=0，回避了J₀的复杂计算。

2.2 基于W-W算法的二分法计算拉索频率

用W-W算法求结构的频率时，最简单的是采用二分法，对索结构也不例外，只是需要根据前述理论对单元划分加以控制。拉索第k阶频率的计算流程图如图 3所示。

3 拉索参数识别 3.1 参数识别原理

通过测试可以方便地获取拉索的多阶自振频率，根据拉索的设计参数，采用本文推导的动力刚度矩阵结合W-W算法则可确定各阶频率的计算值。因此，拉索的参数识别问题实则为一优化问题，其目标函数为：

(18)

式中, f_DSM^j和f_MEAS^j分别为动力刚度法计算和实测所得第j阶固有频率；N为固有频率阶次。

拉索频率的实测值和计算值之间不可避免地存在偏差^[12]，因此，可认为当式(18) 足够小时所对应的识别结果接近各待识别参数真实值。式(18) 可采用多种方法进行求解，本文采用带变异算子的粒子群优化算法(PSO)求解得到拉索参数。

3.2 带变异算子的PSO算法

PSO算法^[14]模拟鸟群飞行觅食的行为，通过鸟之间的集体协作使群体达到最优的目的。PSO系统中每个备选解被称为一个"粒子"，多个粒子共存、合作寻优，每个粒子根据它自身的"经验"和相邻粒子群的最佳"经验"在问题空间中向更好的位置(最优解)"飞行"，每个粒子的速度v_i=(v_i1, v_i2, …, v_iD)和位置x_i=(x_i1, x_i2, …, x_iD)根据如下公式进行更新：

(19)

(20)

式中，N为粒子个数；D为粒子的维数；c₁, c₂为跟踪个体及群体最优值的权重系数，称为加速因子；rand()为[0, 1]之间的随机数；w为惯性因子；p_id为第i个粒子历史最优值；p_gd为第全部粒子搜索到的最优值。

文献[16]针对PSO易陷入局部最优的缺点，提出了带变异算子的PSO算法。该算法后期引入变异算子，使其摆脱后期易陷入局部极优点的束缚，又保持前期搜索速度快的特点。定义带变异算子PSO算法如下：

如果LogjamStep > =MaxStep

  如果SwarmDist < BorderDist

   对邻居子群内所有粒子的所有维，按几率ρ进行重新随机初始化位置和速度

   LogjamStep=0

  否则

  按式(19) 更新粒子位置和速度

  结束

否则

  按式(19) 更新粒子位置和速度

结束

其中，LogjamStep为子群历史最优粒子位置P_g连续不变化或变化极小的迭代次数；MaxStep为连续不变化次数的阀值；SwarmDist为群内所有粒子到历史最佳位置P_g的欧几里得空间距离，可采用最大聚集距离MaxDist，定义如下：

(20)

BorderDist为判断群内粒子聚集程度的距离阈值。

根据文献[16]的研究成果，本文在计算过程中取MaxStep=10，ρ=1/(2D)，BorderDist=0.5^k+1R，k为粒子群变异次数，R为全空间半径。

3.3 识别算法验证

为验证算法的有效性，取三根长度分别为L=10, 40 m和100 m，索力T=2.4 MN，单位长度质量m=34.8 kg/m，抗弯刚度EI=0.442 MN·m²，带有两个中间弹性支承的两端固结拉索(对应的ξ分别为21.4，85.7，214)，中间弹性支承距锚固端的距离为L₁=L₃=3 m，弹性支承刚度K₁=K₂=0.1 MN/m，具体参数真实值及参数识别范围如表 1所示。拉索计算模型如图 4所示。采用ANSYS建立拉索模型，分别用beam3单元和combin14单元模拟拉索和弹性支承，并将计算得到的拉索频率作为初始条件输入本文提出的算法对应的计算程序中，对拉索参数(T，EI，K₁，K₂，L₁，L₃)进行了识别与验证。

图 5示出了L=10 m时拉索参数基于带变异算子的PSO的识别过程。可以看出，随机生成的20个"粒子"在400次的更新"飞行"中迅速接近目标值，并进行了多次变异，摆脱了解陷入局部最优点的束缚。由表 1所示识别结果可知，尽管各参数的识别范围达到±20%甚至更多，基于PSO算法的频率拟合法仍能够较高精度地识别拉索各参数值，其中EI, lgK₁, lgK₂, L₁, L₃等参数识别误差最大为4.24%，但最重要的拉索索力T的识别误差均小于1%，完全满足工程要求。

4 实例分析

广州海印桥是一座双塔单索面预应力混凝土斜拉桥，主桥跨径布置为35+(85.5+175+85.5)+35 m=416 m。全桥拉索共186根，每根拉索两端均设有高阻尼橡胶阻尼器，以减少拉索振动。采用东京测器的动态应变仪(DC-204R)对全桥拉索进行了频率测试，采样频率为100 Hz。选取北塔边跨3根典型拉索进行研究，拉索参数如表 2所示，典型拉索振动时程频谱如图 6所示。

用动力刚度法及带变异算子的PSO算法对拉索的索力T、减振器刚度K_d1, K_d2以及拉索抗弯刚度EI进行识别。参数识别过程如图 7所示，可以看出，随机生成的20个"粒子"经过近200次的"飞行"，索力T的识别值已经趋于收敛，3根拉索在参数识别过程中发生了1~4次不等的变异，其中NB15-1索变异较为明显。

拉索的参数识别结果如表 3所示，由于该桥为旧桥，拉索两端没有设置穿心式压力传感器，缺乏拉索的真实索力值，表 5中给出了拉索的成桥索力设计值，注意该值并不是目前索力真实值，此处仅做参考。同时本文还用另外两种方法对表 2中3根拉索的索力值进行了计算以与本文算法得到的索力值进行对比。方法一：根据文献[18]提出的采用高阶振动模态来评估拉索索力的方法，计算得到考虑拉索抗弯刚度时的索力值列于表 5中；方法二：采用ANSYS建立各拉索有限元模型，通过不断修改拉索张力来计算拉索的前10阶频率，取计算得到频率与实测频率最接近时对应的拉索张力为其索力值，计算结果见表 5，模型中阻尼器刚度取值参考文献[19]。

同时，表 5中还给出了上述两种计算方法得到的索力值与本文算法得到的索力值之间的差值百分比(括号中数据)，由表 5中的计算结果可以看出，3种方法计算结果非常接近，最大差值百分比仅为6.81%，说明了本文算法的可靠性；同时，3种索力计算结果均与成桥索力设计值相差较小，也能部分反映本文计算结果的有效性。

对于两端带减振器或阻尼器的斜拉索，在索力测试时计算索长的确定一直是个棘手的问题，特别是对于短索及两端减振器或阻尼器安装位置到锚固点距离占总索长比例较大的情形。文献[18](方法1) 计算时假定阻尼器刚度无穷大，取计算索长为两端减振器中心之间的距离，为近似计算，方法2则需要反复试算，且需要先根据经验事先假定阻尼器刚度，而本文的算法则仅需在测得拉索各阶频率的情况下即可较高精度地识别出拉索的索力。

上述待识别参数中，一般情况下索力在实际工程最为重要，而其计算精度又受其他参数的影响。文献[20]对此作了相应的参数分析，得出拉索减振器安装距相应端部较远时，索力计算不能忽略减振器刚度的影响，而对于不同ξ值的拉索，拉索边界条件对索力计算精度的影响也将随之变化，因此一般的近似计算公式将出现局限性，而本文提出的算法则可以在任意拉索边界条件以及减振器安装位置情况下对拉索进行参数识别，因此较其他方法具有更好的适应性。

5 结论

(1) 推导了拉索单元的单元动力刚度矩阵，基于动力刚度法，提出采用控制单元划分方式的W-W算法求解拉索的振动频率，拉索参数识别实例结果表明，基于动力刚度矩阵的W-W算法能够较高精度地得到拉索的振动频率；

(2) 根据拉索振动频率拟合的参数识别方法，采用带变异算子的粒子群优化算法进行拉索参数识别，仿真分析及实桥试验表明，基于动力刚度法与PSO算法的参数识别方法能够获得较好索力测试精度，适应性较为广泛。

(3) 对于两端带减振器或阻尼器的斜拉索，在索力测试时计算索长的确定一直是个棘手的问题，应用本方法可以识别出拉索的抗弯刚度、减震器刚度、拉索索力等参数，比常规的偏经验选取拉索计算索长进行索力识别更为合理。