基于IMPSCO和改进Newmark-β算法的结构系统及激励辨识研究

麻胜兰， 姜绍飞， 陈志刚

(福州大学 土木工程学院，福州 350108)

基于IMPSCO和改进Newmark-β算法的结构系统及激励辨识研究

麻胜兰， 姜绍飞， 陈志刚

(福州大学 土木工程学院，福州 350108)

基于振动的土木工程结构系统识别中，常用的识别方法如基于Newmark-β算法的方法，常常需要已知激励，而事实上激励和完整的输出响应是很难测得的。基于此，在现有的Newmark-β算法中引入多项式拟合来解决结构系统识别误差带来的漂移现象，进而提出了改进的Newmark-β算法，并结合改进协同粒子群算法(IMPSCO)，给出了仅用部分结构响应的系统识别和激励辨识方法。最后，通过一七层钢框架数值算例和实验室模型试验验证了方法的有效性和可行性，并探讨了噪声、输出响应完整性的影响。研究表明，所提方法能准确地实现未知激励和部分实测响应状态下的结构系统和未知激励的识别，而且具有较强的容噪性和鲁棒性。

未知激励；系统识别；改进的Newmark-β算法

结构系统参数识别是结构健康监测的核心部分，因此国内外专家提出了很多系统识别方法[1-4]。其中，时域下基于结构振动响应的识别方法受到了广泛关注，通过直接识别结构系统的物理参数，如刚度、质量或阻尼等，可以有效确定损伤位置、损伤程度和评估结构状态。这些方法如最小二乘法(Least Square,LS)[5]、卡尔曼滤波(Kalman Filter,KF)[6]和遗传算法(Genetic Algorithm,GA)[7]等。但分析发现，这些方法大多仅适用于输入与输出都已知的情形，而实际工程中，由于环境的复杂性，一般难以测得结构实际受到的激励大小，甚至也难以测得完整的结构输出响应。因此，研究未知激励及部分实测输出响应状态下的结构系统辨识方法，意义重大。

国内外专家学者开展了一系列研究，取得了一定成果。Chen等[8]利用改进的迭代最小二乘法同时识别结构系统参数和未知输入力；此后，Chen等[9]又提出利用蚁群算法同时识别结构系统参数和未知激励，并通过试验验证了方法的有效性；此外，Sun等[10]将改进Newmark-β算法与人工蜂群算法结合，来同时识别未知激励及结构系统参数。但研究发现Chen的方法并不适用于仅有部分实测响应的情形；而Sun的方法通过将加速度积分成速度和位移，进而根据动力方程求得未知激励，但实际工程中，加速度积分会由于“平移项”和“趋势项”的影响[11]使速度、位移和识别激励发生漂移，导致识别错误。

因此，本文的目的在于研究一种未知激励及部分实测响应下，能同时识别激励与结构系统参数的时域下的辨识策略。时域下的结构系统识别通常需要与各种数值积分方法相结合，其中最常用的数值方法为Newmark-β算法[12-13]，但是常规Newmark-β算法需要已知激励，为此，通常做法是通过易测得的加速度积分得到速度和位移响应，再根据结构动力方程求得激励，然后进行系统参数识别。但如上所述，加速度积分存在趋势项误差，会使速度、位移和识别激励发生漂移，并导致系统参数识别错误。鉴于此，本文借鉴多项式去趋势的思想，对常规Newmark-β算法进行改进，利用多项式拟合来修正加速度积分导致的识别激励的漂移现象；在此基础上，结合改进协同粒子群算法(IMPSCO)，提出了一种适用于未知激励及部分实测响应下的结构系统参数识别与激励辨识方法，并通过一7层框架数值仿真和试验研究验证了方法的有效性和可行性。

1 改进的Newmark-β算法

本文对Newmark-β算法进行了改进，具体步骤如下：

步骤1 初始计算加速度、速度、位移和力：

①确定结构的质量、阻尼和刚度矩阵M,C,K,及时间步长Δt；

步骤2 由每个时间步p计算第p+1个时间步的结构响应和激励：

③按照传统Newmark-β算法计算结构响应：

步骤3 对下一个时间步进行循环令p=p+1，重复步骤2中的①～④；

步骤6 识别激励修正：fiden=f-ffit。

步骤7 将fiden作为已知激励，并重复执行步骤2中的③步骤直到所有时间步循环结束。

总体而言，改进的目的在于解决以下两个问题：

(1) 激励未知下，程序的执行问题及未知激励的识别；

(2) 修正可能发生的识别激励的漂移现象。对于前者，本文借鉴Sun的方法，利用已知的实测加速度积分成速度和位移，进而根据动力方程求得未知激励。然而，Sun在计算下一时间步的速度响应时采用的上一时间步得到的速度和位移经常规Newmark-β算法公式(类似本文所提算法步骤步骤2步骤中的①～③)计算得到的数值加速度，这无疑会增加计算误差，因而本文将其改用为实测加速度(步骤2步骤中的①～②)；此外，本文在计算中加入了识别激励的更新(步骤2步骤中的④)，以进一步减少可能存在的由于积分引起的识别误差。对于后者，本文利用多项式拟合去除可能发生的激励漂移现象，得到修正的识别激励，并将其作为已知值重新计算结构响应(步骤5～8)。

2 系统识别与激励辨识策略

如上所述，本文将IMPSCO与改进Newmark-β算法结合起来，提出了一种适用于未知激励及部分实测响应状态下的结构系统识别与激励辨识策略，其流程如图1所示；IMPSCO是对MPSCO的改进，通过对局部最优粒子的替换处理，提高了算法的稳定性和鲁棒性，具体过程可详见文献[14]。此外，考虑到加速度相较容易获得，本文基于加速度响应进行研究。因此，在激励位置、激励部位及其相邻部位加速度已知时，即可利用本策略识别结构系统参数和未知激励。具体步骤如下：

(1) 对结构系统参数进行编码并确定参数搜索范围，本文考虑质量已知情形，对于阻尼，则采用Rayleigh阻尼模型，因此需要编码的结构系统参数θ为刚度ki(i=1,2,…,n)和阻尼比ξj(j=1,2)；

(2) 设置IMPSCO的基本算法参数，本文设置如下：子种群数量为Ns=3；子种群大小为ms=30；进化代数me=500；学习因子c1=c2=2；最差次数限值Ie=6；w=[0.4,0.9]，线性递减，且当迭代次数大于45时，取0.4；

(3) 在各参数搜索范围内随机产生若干随机粒子(即不同参数值的组合)，并按第一节，利用改进Newmark-β算法计算每个随机粒子下的结构响应及未知激励(fiden)；

(4) 以实测加速度响应和计算所得加速度响应为基础建立适应度函数F，如式(1)所示，并计算每个随机粒子的适应度值；

(1)

式中，amea和acom分别为结构实测加速度响应和数值计算的加速度响应;k是结构测点数;l为加速度响应数据的个数。适应度计算值越大，识别结果越接近实际值。

(5) 替换掉适应度值最差的粒子，并对每个粒子按IMPSCO的设定进行粒子进化；

(6) 同样按第一节，利用改进Newmark-β算法计算每个进化粒子组合下的结构响应及未知激励(fiden)，并重复步骤(4)和(5)；

(7) 重复步骤(6)，进行粒子循环迭代进化，直到达到设定的最大迭代次数，迭代结束，输出适应度值最高的粒子，即最优系统参数和未知激励(fiden)。

图1 结构系统识别与激励辨识策略流程图

3 数值算例

3.1 结构模型

图2为一7层剪切型框架结构，刚度为k1=…=k7=375 kN/m;质量为m1=…=m6=3.78 kg和m7=3.31 kg;采用Rayleigh阻尼模型，前两阶阻尼比为ξ1=ξ2=2%。顶层受随机激励，采用常规Newmark-β算法计算结构动力响应，并取加速度响应用于分析，采样频率1 000 Hz，采样时间5 s。

3.2 识别过程与结果

本例中，质量已知，以结构各层层间刚度及阻尼比为未知，进行编码识别，即θ=[k1,k2,…,k7,ξ1,ξ2]。参数搜索范围取理论值的0.5倍～2倍，即ki∈[187.5,700]kN/m,i=1,2,…,7;ξj∈[0.01,0.04],j=1,2。按图1流程，识别结构参数θ及未知激励fiden，改进Newmark-β算法中多项式拟合阶数取n=4。识别时，仅考虑第2,3及5～7层加速度响应以模拟响应不完整情形。无噪声下，阻尼比的识别值分别为2.05%和2.09%，结果较为精确；刚度识别结果如表2所示；激励识别结果如图3所示，为简便，图中仅给出0～0.5 s识别结果。

图2 结构模型

从表2可以发现，结构刚度的识别结果与理论值非常接近，最大误差仅为0.20%；而从图3(a)可以发现，激励的识别值也与理论值非常吻合。此外，图4为未经多项式修正时的激励识别结果(即识别过程没有加入第一节中步骤5～8)，可以发现即使未利用多项式进行修正时，识别激励也没有明显发生漂移，这是由于本例为数值算例，在正向计算结构数值响应与反向识别结构物理参数时，实际上都是基于Newmark-β算法的基本假设进行的，因此识别中不会出现明显的激励漂移现象，因此，识别时多项式的修正值为0，但实际工程中，结构响应是实际测得的，并不完全满足这一假设，因而容易产生趋势误差，导致识别激励发生漂移，这将在下一节试验中得到验证。以上分析表明，无噪声下，本文所提策略可以精确识别结构系统参数和未知激励。

表2 刚度识别结果

(a) 无噪声

(b) SNR=20 dB

图4 未修正下激励识别结果

3.3 比较与讨论

为验证本文方法的适用性和鲁棒性，从噪声水平和输出响应的完整性两个方面进行分析讨论。

(1) 噪声影响

按式(2)向原始加速度响应中添加噪声，信噪比(SNR)分别为25 dB和20 dB。

SNR=20lg(1/ε)=20lg(Asignal/Anoise)

(2)

式中:yia和yi分别代表原始位移响应信号和噪声污染后的信号;R是均值为0，偏差为1的正态分布随机数;ε为噪声水平;A为信号幅值。

同样，按图1流程进行识别, 刚度和激励的识别结果如表2和图3所示。可以看出，SNR=25 dB和20 dB时，刚度的最大识别误差也仅为0.81%和2.20%，而激励的识别结果也仍然非常精确。此外，SNR=25 dB时，阻尼比的识别值分别为2.12%和2.17%；SNR=20 dB时，则分别为2.24%和2.25%；实质上Rayleigh阻尼比为2%。可见阻尼比的识别误差与刚度识别一致，随噪声增加，其中阻尼比误差更大一些，这是由于本例为小阻尼结构，阻尼对系统响应影响较小所致，但总体上其误差尚在可接受范围内。这些表明，本文所提策略具有较好的容噪性和鲁棒性。

(2) 输出响应的完整性

为验证输出完整性对本文方法的影响，无噪声下，分别利用第5～7层、第2,3,5～7层和第1～7层的加速度响应进行识别。刚度和激励识别结果如表3和图5所示。

1～7层加速度都参与识别时，阻尼比分别为2.08%和2.02%，最大误差4.00%，而仅有5～7层的加速度响应参与识别时，其值分别为2.19%和2.25%，最大误差为12.50%。可见参与识别的加速度响应个数对于阻尼比的识别影响较大，但从表3可以发现，其对刚度的识别影响较小。当只有5～7层的加速度响应参与识别时，刚度的最大识别误差也仅增大到0.74%，激励的识别结果也仍然非常精确(图6(b))。这些表明，虽然加速度响应的完整性会影响参数及激励识别结果，但影响较小，本文策略仍具有较高识别精度。

表3 不同测点组合下的辨识结果

(a) 第1～7层

(b) 第5～7层

4 试验验证

4.1 试验模型

如图6所示，为2跨×1跨的7层钢框架缩尺模型，相应的数学简化模型如图2所示。梁、柱间距均为200 mm。梁为空钢管，柱为薄钢板，且均为300 W级钢材。构件特性如表4所示,据此可求得结构各层质量为m1=m2=…=m6=3.78 kg，m7=3.31 kg，完整状态下结构各层理论刚度均为375 kN/m。利用激振器在结构顶层施加随机激励，并利用力传感器记录激振力大小，采样频率5 000 Hz。在结构各层布置传感器以测量各层加速度响应。

图6 实验室钢框架模型

参数梁柱截面尺寸/mm25×25×325×4.6截面面积A/m2286×10-6115×10-6惯性矩I/m42.41×10-87.78×10-10杨氏模量E/Pa206×109206×109体密度ρ/(kg·m-3)78507850

4.2 结构系统修正

利用结构刚度理论值计算出结构的理论频率并与试验中测得实测频率进行对比，如表5所示。

表5 结构固有频率

从表5可以看出，理论频率与实测频率偏差较大，可见理论刚度值不符合结构实际条件，这是由于试验材料自身的缺陷、焊接的影响、试验设备的测量误差、试验人员人为的因素及测量时噪声的影响等多种因素引起的，因此需要对其进行修正。

为此，在激励已知下，利用IMPSCO进行刚度识别与修正。首先对结构各层层间刚度及阻尼比进行编码，即θ=[k1,k2,…,k7,ξ1,ξ2]。参数搜索范围为ki∈[187.5,700]kN/m,i=1,2,…,7;ξj∈[0,1],j=1,2。刚度的识别结果如表6所示。而后利用识别刚度计算结构的固有频率，如表5所示，可以看出，识别频率与实测频率更加吻合，从而表明识别的刚度值更符合实际结构状况。因此，后文分析时以此刚度识别值作为激励未知下的刚度基准值，以便对比。

表6 不同测点组合下的刚度识别结果

4.3 识别过程与结果

以第2，3及5～7层的加速度响应作为已知输出响应，以模拟响应不完整情形，对各层刚度及阻尼比进行编码，搜索范围与4.2节一致，并按图1流程进行识别，多项式阶数同样取n=4。阻尼比的识别结果分别为0.68%和2.02%；刚度和激励的结果分别如表6和图7(b)所示。图7(a)为未经多项式修正时的激励识别结果(第一节中步骤5～8步)。

(a) 未修正(不完整响应,第2,3,5～7层加速度)

(b) 修正(不完整响应,第2,3,5～7层加速度)

(c) 未修正(完整响应,第1～7层加速度)

(d) 修正(完整响应,第1～7层加速度)

图7 激励识别结果

Fig.7 The identified results of excitation

从表6可以看出，激励未知下，刚度的识别结果与激励已知时的识别结果较为接近,最大误差仅为9.84%；此外，从图7(a)可以发现，识别激励由于加速度积分的趋势误差发生了明显漂移；但经本文所提策略修正后，激励的识别结果明显没有了漂移现象，且与实测值较为吻合，如图7(b)所示。以上分析表明，本文策略可以有效识别结构系统参数和未知激励，并消除激励的漂移现象。

4.4 比较与讨论

为了验证本文方法的鲁棒性和有效性，对同样的结构模型进行了激励未知及完整响应状态(7层加速度都已知)下的系统辨识。首先对结构各层层间刚度及阻尼比进行编码，参数搜索范围与4.2节一致，同样按图2流程进行识别，多项式拟合阶数同样取n=4。阻尼比的识别值分别为0.52%和1.86%，与4.3节的识别结果较为接近，且都比较小；刚度和激励的结果如表6和图7(d)所示，图7(c)同样为未经多项式修正时的激励识别结果。

从表6可以发现，相较完整输出响应情形，仅有部分实测响应时，刚度的最大误差从5.85%增加到9.84%，可见识别误差虽有所增大，但仍有较高的识别精度；此外，从图7(b)和图7(d)对比发现，相较完整输出响应情形，部分实测响应下激励识别值与激励实测值的吻合度虽有所不如，但总体上仍然较为吻合。以上对比表明，虽然加速度响应的不完整性对参数和激励的识别有所影响，但本文所提策略仍表现出较高的识别精度和鲁棒性，能满足工程需要。

5 结 论

本文通过研究，得出了以下结论：

(1) 通过改进常规Newmark-β算法并结合IMPSCO，提出了系统辨识策略，有效解决了仅有部分实测输出响应下的结构系统识别和激励辨识问题。

(2) 响应完整性确会影响识别结果，但本文策略对于仅有部分实测输出响应下的系统参数识别仍具有较高的识别精度和较强的鲁棒性。

(3) 研究表明，所提改进Newmark-β算法确实能有效消除识别激励的漂移现象，保证了系统识别结果的准确性。

数值模拟和试验研究表明，所提策略能有效识别结构系统参数和激励，初步验证了该方法的可行性和有效性，且具有较好的容噪性，但仍需要更多的试验和工程实践来进一步检验。此外，本文策略仅适用于激励位置已知的情形，对于激励位置未知情形，还有待进一步研究和扩展。

[1] 张延年, 范鹤, 董锦坤. 基础隔震结构的系统参数优化设计[J]. 沈阳建筑大学学报, 2006,22(4):529-533.

ZHANG Yannian, FAN He, DONG Jinkun. Parameter optimization of base-isolated structure[J]. Journal of Shenyang Construction University, 2006,22(4):529-533.

[2] TANG H, ZHANG W, XIE L, et al. Multi-stage approach for structural damage identification using particle swarm optimization[J]. Smart Structures and Systems, 2013, 11(1): 69-86.

[3] WANG F G, LING X Z, XU X, et al. Structural stiffness identification based on the extended kalman filter research[J].Abstract and Applied Analysis, 2014(2):267-290.

[4] YANG J N, HUANG H, LIN S. Sequential non-linear least-square estimation for damage identification of structures[J]. International Journal of Non-Linear Mechanics, 2006, 41(1):124-140.

[5] LIN J W, BETTI R, SMYTH A W, et al. On-line identification of non-linear hysteretic structural systems using a variable trace approach[J]. Earthquake Engineering & Structural Dynamics, 2001, 30(9):1279-1303.

[6] HOSHIYA M, SAITO E. Structural identification by extended Kalman filter[J]. Journal of Engineering Mechanics, 1984,

110 (12):1757-1770.

[7] PERRY M J, KOH C G, CHOO Y S. Modified genetic algorithm strategy for structural identification[J]. Computers & Structures, 2006, 84(8/9):529-540.

[8] CHEN J, LI J. Simultaneous identification of structural parameters and input time history from output-only measurements[J]. Computational Mechanics, 2004, 33(5): 365-374.

[9] CHEN J, CHEN X, LIU W.Complete inverse method using ant colony optimization algorithm for structural parameters and excitation identification from output only measurements[J]. Mathematical Problems in Engineering, 2014, 18(1):249-269.

[10] SUN H, BETTI R. Simultaneous identification of structural parameters and dynamic input with incomplete output-only measurements[J]. Structural Control and Health Monitoring, 2014, 21: 868-889.

[11] JUNHEE K, KIYOUNG K, HOON S. Autonomous dynamic displacement estimation from data fusion of acceleration and intermittent displacement measurements[J]. Mechanical Systems and Signal Processing, 2014,42(1/2):194-205.

[12] MAUNG T S, CHEN H P, ALANI A M. Real time structural damage assessment from vibration measurements[J]. Journal of Physics: Conference Series, 2011, 305(1):12096-12103.

[13] SHIN S B, OH S H. Damage assessment of shear buildings by synchronous estimation of stiffness and damping using measured acceleration[J]. Smart Structures and Systems, 2007, 3(3):245-261.

[14] JIANG S F, WU S Y, DONG L Q.A time-domain structural damage detection method based on improved multi-particle swarm coevolution optimization algorithm[J].Mathematical Problems in Engineering, 2014, 44(1):77-85．

Identification of a structural system and its excitation based on IMPSCO and modified Newmark-βalgorithm

MA Shenglan, JIANG Shaofei, CHEN Zhigang

(College of Civil Engineering, Fuzhou University, Fuzhou 350108, China)

In vibration-based structural system identification (SI) of a civil structure, the commonly used SI method, such as, the one based on Newmark-βalgorithm usually needs the known excitation. In fact, it is difficult to measure the actual excitation and full output responses of a structure system. Here, a modified Newmark-βalgorithm was proposed through introducing polynomial fitting into Newmark-βalgorithm to effectively solve the drift phenomenon brought by the error of SI. Furthermore, combining the modified Newmark-βalgorithm with the improved multi-particle swarm co-evolution optimization (IMPSCO), an identification approach of SI and excitation adopting responses of parts of a structure system was presented. Finally, numerical simulation for a seven-story steel frame and its model tests were conducted to verify the effectiveness and feasibility of the proposed method. The effects of noise and completeness of structural response were also discussed. The results showed that the proposed method can be used to effectively identify a structural system and its excitation only adopting responses of parts of the structure system; it has stronger noise tolerance and robustness．

unknown excitation; system identification; modified Newmark-βalgorithm

国家十二五科技支撑计划(2015BAK14B02)国家自然科学基金(51278127)

2015-09-22 修改稿收到日期：2016-06-28

麻胜兰 女，博士,1986年生

姜绍飞 男,闽江学者特聘教授，博士，博士生导师,1969年生 E-mail：cejsf@fzu.edu.cn

TU317

A

10.13465/j.cnki.jvs.2017.15.004