基于形函数方法快速识别结构动态荷载的试验验证

张青霞，段忠东，Lukasz Jankowski，王 丰

(1.大连民族学院 土木建筑工程学院，大连 116600;2.哈尔滨工业大学 深圳研究生院，深圳 518055;3.哈尔滨工业大学 土木工程学院，哈尔滨 150090;4.Smart－Tech Centre，Institute of Fundamental Technological Research，Polish Academy of Sciences，02 －106，Warsaw，Poland)

准确的结构动态荷载信息是确定结构动态行为的关键因素之一，对于结构损伤识别、安全评定等具有重要的作用。此外，它为灾难性事故的分析评估提供了科学依据，有益于预防和避免类似事故的重发生。荷载识别属于结构动力学中的反问题分析，由结构的动力响应和其动态特性确定。

目前的荷载识别方法［1－4]中，大多数属于频域［5，6]、时域［7，8]或小波域［9]内分析，基本思想是通过计算实测的结构响应和提前获取的系统脉冲响应函数的反卷积计算荷载，也就是所谓的反卷积方法［2]。此外，还有时间有限元法［10]，基于鲁棒观测器(robust observer)法［11]，逆结构滤波法(ISF)［12]等。每种方法具有不同的优缺点，但均面临反问题分析中常存在的病态性问题，如解的唯一性、稳定性等。反卷积法以其形式简单、可灵活采用时域和频域内分析的特点应用较广泛;但若测量时间长或采样频率高时，则导致脉冲响应函数(传递函数)组成的系数矩阵维数很大，从而使该方法数值求解困难。

文献［8] 提出荷载形函数方法来改善荷载识别反卷积方法中的不足，利用有限元中的形函数概念逼近动态荷载曲线，从而把识别未知荷载的时间历程转变为识别数目有限的荷载形函数的权重或系数，大大降低了系数矩阵的维数，提高计算效率。另外，荷载形函数是连续的，在一定程度上能够过滤噪声的影响，提高逆问题对噪声的鲁棒性。本文对荷载形函数的相关参数设置进行了理论补充，并通过试验验证了该方法的有效性。

1 反卷积方法

假定所考虑的结构为线弹性而且在荷载作用期内为小变形即忽略几何非线性，那么结构的响应线性依赖于外部荷载作用。设结构初始状态为零，根据杜哈姆积分，结构的响应可由外部荷载通过卷积表示为:

式中:α为测点编号;fi为作用在结构上的第i个外部荷载，hαi(t)为测点α和荷载fi之间的脉冲响应函数。

实际应用中，结构的响应一般通过测量或数值模拟得到，从而是离散化的，因此考虑测量时间内所有的离散的结构测点响应，并按一定顺序排列，式(1)可以表示为如下的矩阵形式:

式中yM汇集了所有测点在测量时间内的离散响应;f汇集了所有外部荷载的离散时间历程;矩阵h是分块矩阵，子矩阵是由测点和外部荷载之间的相关脉冲响应函数构成的Toeplize矩阵。

由式(1)、(2)知，通过解实测结构响应和系统脉冲响应函数计算反卷积可以得到未知荷载，即求解方程(2)。这是一个线性方程组，为保障所求解的唯一性，需要独立测点(传感器)的数目至少等于可能的未知荷载数目。在实际应用中，需要预先知道或假定荷载可能存在的位置信息以确定传感器的数目和构造对应的脉冲响应函数。

解式(2)就是解线性方程组Ax=b或求解min Ax=b2，往往面临病态问题。实测响应中一个小的扰动，比如无法避免的测量误差，将会引起很大的荷载识别误差，即方程的解对扰动非常敏感［3]。因此，有必要对解进行数值正则化，常用的正则化技术有截断奇异值分解方法(TSVD)或 Tikhonov方法［3]。

另外，式(2)中方程组和未知数的数目与时间步数成正比，当采样频率高或测量时间较长时，比如离线荷载识别中，线性方程组可能会很大，导致计算操作难以实现。利用提出的荷载形函数方法可以改善荷载识别方法中存在的病态问题，减小计算工作，提高识别效率。

2 荷载形函数方法

2.1 基本原理

荷载形函数的基本思想源于有限元中的形函数。不同的荷载形式决定了不同类型的形函数类型，如荷载只作用在一点时，其时间历程是一个曲线，可以采用梁单元的形函数，本质上是采用曲线插值的方法对荷载曲线进行逼近。如果结构上作用一个线荷载，其时间历程就是一个曲面，需要选用板的形函数进行曲面插值。对于面荷载，则需要借用实体单元的形函数进行插值。因此选择形函数类型时，需要根据经验判断荷载的类型。

本文假定荷载的时间历程是一个曲线，如图 1，可以比拟为一个“时间梁”的变形，其中荷载的时间轴比拟为

图1 荷载的时间历程曲线Fig.1 The load time history

“梁”长度方向的位置坐标。只考虑梁的弯曲变形，即每个结点有竖向和转动两个位移。假如把荷载历程按时间划分为四段，那么该“时间梁”单元有五个结点，十个形函数，定义为荷载形函数，如图2，其意义为“时间梁”某个结点在竖向或转动方向发生单位位移引起的荷载历程。记荷载形函数矩阵N=［N1N2…Nnl] ，nl是形函数的个数。

图2 荷载形函数Fig.2 Load shape function

根据形函数的意义可知，动态荷载时间历程f(t)可由若干个形函数近似表示:

式中:α=［α1α2… αnl]T，αi是对应第 i个形函数的“时间梁”结点的位移，称为该形函数的系数。

2.2 形函数的选取

图2给出了形函数的形式，具体表达式参考有限元中梁单元形函数。在实际应用中，荷载历程一般被均匀划分为若干段。设相邻每段两个结点的时间长度定义为形函数的半个周期Tf/2，则形函数的频率为ff=1/Tf，两个结点之间为lnt=Tf/(2Δ)=fs/(2ff)个时间步，其中fs为采样频率，采样时间间隔为Δ=1/fs。若分析的离散荷载时间为t，共含nt个时间步，则形函数的个数nl:

由上可知，形函数内时间步数取决于采样频率fs和形函数的频率ff，而形函数个数取决于形函数的频率ff和测量时间t。测量时间一定时，形函数个数由形函数频率决定，形函数的频率选取会影响识别精度。结构动态荷载往往含有多个频率成分，为能够准确逼近动态荷载曲线，要求形函数的频率大于荷载所分析的最高截断频率。形函数的频率小于荷载的截断频率时，它不能准确地逼近荷载;而形函数的频率较高时，取极限类似单位脉冲，识别结果对噪声敏感。另外，较高的形函数频率意味着形函数个数多，这时系数矩阵的维数大，计算效率低。因此，选取合适的形函数频率是基于形函数法识别荷载的关键之一。对未知荷载而言，鉴于线性结构的响应线性依赖于外部荷载，确定其所需的形函数频率可以通过对结构实测响应的频谱分析决定。

设F(ω)代表响应的傅里叶变换，定义下式:

对于一般荷载，可取r(ωc)=95%对应的频率ωc为截断频率，视ωc为荷载可能的最高频率，即把ωc作为逼近未知荷载所需的形函数的参考频率。实际操作中根据荷载形式，r(ωc)的取值可适当增减。那么两个相邻结点含的时间步lnt=fs/2ωc，代入式(4)可确定所需的形函数个数nl。

2.3 荷载识别

实际应用中，根据式(3)，式(2)中离散的荷载历程可以表示为:

式中:Nf=Inf×nf⊗N 是 nfnt× nfnl的矩阵，N 是 nt× nl维的离散的荷载形函数矩阵。是nfnl维的列向量，由各个荷载的形函数的系数组成，αi(j)是荷载i的第j个形函数的系数。

把式(7)代入式(2)得到:

式(8)中的脉冲响应函数矩阵、荷载形函数矩阵均可以提前构造。形函数的个数nl一般远小于测量时间步数nt，因此把式(2)直接计算荷载转化为先求解式(8)中荷载形函数的系数，然后计算荷载，使未知数个数大大减少，显著提高计算效率。另外，荷载形函数是连续光滑的，根据式(8)计算荷载在一定程度上能削减噪声的影响。

利用式(8)计算荷载，其系数矩阵的维数得以减小，采用截断奇异值分解方法(TSVD)对解进行数值正则化更简单快捷，即:

式中:ui是酉矩阵 U=［u1，…，unfnt] 的第 i个列向量，vi是酉矩阵V=［v1，…，vnfnl] 的第i个列向量，酉矩阵U，V由式(8)中系数矩阵B的奇异值分解得到。σi是矩阵B的第i个正奇异值。q指选用前q个正奇异值求解。奇异值截断数目可以根据L－曲线计算。通过L－曲线的拐点来权衡响应残差和荷载光滑度。

根据公式(9)，系数矩阵很小的奇异值会放大结构测量中的响应误差的影响，从而使识别的结果严重偏离真实值。利用反卷积方法识别荷载时，公式(2)中系数矩阵由脉冲响应函数组成，而脉冲响应的振动一般比较剧烈，频率成分含量比较丰富，特别是高频成分，系数矩阵中小的奇异值可认为是由脉冲响应中的高频成分导致的，因此在计算反卷积时需要通过截断奇异值分解等正则化方法去掉较小的奇异值，以改善解的识别精度。当采用荷载形函数识别荷载时，系数矩阵由脉冲响应函数组成的矩阵和形函数矩阵相乘得到(见式(8))，该系数矩阵中的每列元素相当于对应每个荷载形函数引起的结构响应。因为荷载形函数是连续光滑的，它引起的结构响应也相对光滑，相对于单位脉冲引起的结构响应而言，它含有高频成分很少或几乎没有，使系数矩阵中只存在少数几个或不存在特别小的奇异值，因此采用荷载形函数识别荷载时，系数矩阵只截断少数几个奇异值或不需截断奇异值。

另外，已知结构的有限元模型，荷载形函数方法结合“移动时间窗”［13]可以推广到在线识别，也就是式(8)中的实测响应yM被替代为，其中，yM(n)是第n个时间窗的采样数据;假定前面时间窗内的荷载已被识别，该时间窗内的结构初始状态和非零初始状态引起的桥体的自由振动根据结构模型易求。

3 数值算例

首先借助一个两跨连续梁模型(图3)来验证所提出的荷载识别方法的有效性。

梁每跨5 m，工字型组合截面55 mm×100 mm×7.2 mm ×4.5 mm。弹性模量 E=210 GPa，密度 ρ=7 800 kg/m3，截面面积为1 200 mm2，受弯惯性距为1 943 800 mm4。有限元梁模型被均匀划分为20个单元，在距左端2.26 m、5.26 m的截面底部分别布一个应变传感器，记为 S1，S2。取前两阶瑞利阻尼比为0.01。

假定结构承受两个竖向荷载f1，f2(图4)，分别作用在距梁左端4.5 m、6.5 m截面上。结构动力响应计算采用Newmark逐步积分方法，Newmark参数α=0.25，β=0.5。积分步长0.001 s，即采样频率1 000 Hz。测量时间为0.9 s。考虑响应含有5%的高斯随机噪声。

假定外部荷载未知，采用荷载形函数法进行识别，并与荷载识别的基本方法—反卷积法(解式(2))进行比较。鉴于荷载识别等反问题往往是病态的，采用TSVD正则化方法求解。奇异值截断的个数k即正则化程度通过L－曲线确定。

考虑5%的高斯噪声，反卷积法识别荷载的L－曲线(图5)的拐点对应k=470，也就是小于最大奇异值0.14%的奇异值都被截断，识别的相应荷载(见图6)与真实值相比有明显振荡，说明直接利用反卷积方法识别荷载对噪声敏感。

利用形函数进行荷载识别时，荷载形函数通过测点响应进行傅里叶变换确定。两个测点响应的各频率含量相似，以测点S1的傅里叶变换为例进行分析，如图7。选取ωmax=100 Hz，根据式(6)得到的形函数的频率ωc=14.43 Hz，见图8，根据此确定每个荷载需要62个形函数。考虑5%的高斯噪声，得到的L－曲线如图9。该曲线没有典型的L－曲线拐点，说明解式(8)条件数好，k=1即可得到不错的识别结果，如图10，与真实值非常接近。

从上面的分析可知荷载识别的主要计算工作来自系数矩阵的奇异值分解，对一个m×n维的矩阵(m≥n)的计算量为O(mn2)阶［14]。基本反卷积法的系数维数矩阵为1 802×1 802，利用形函数方法后变为1 802×124，因此这里形函数方法能使计算量降低几乎两个数量级。基本反卷积方法识别荷载对噪声敏感，需要通过合适的奇异值截断来进行，而且仍存在明显的识别误差。形函数方法不但降低计算量，而且能够在一定程度上过滤噪声的影响，识别结果抗噪能力强。

4 试验验证

图11 试验装置Fig.11 Experimental setup

试验装置如图11所示，试件为一个悬臂铝梁，长136.15 cm，矩形截面 2.7 cm ×0.31 cm。固定端夹在一个稳定框架上。梁的弹性模量为70 GPa，密度是2 700 kg/m3。采用压电激励器APA(Amplified Piezo Actuator)对试件施加一个纯弯矩激励，结构动态响应利用贴在梁上的三个压电应变片(PVDF)测量，记为传感器S1～S3。试验中设计了光滑连续荷载、三角形荷载，见图12。为保证测得的样本数据含有必需的结构响应信息，采样频率为2 500 Hz。两种工况对应的测点传感器的响应见图12。

图12 两种工况下的实测结构响应和激励Fig.12 Measured responses and excitations

利用移动时间窗进行荷载在线识别，每个窗内采用形函数方法识别。首先以测点S1为代表，对各工况的前2 000个点的实测响应进行傅里叶变化，根据式(6)计算形函数的参考频率，从而确定每个窗内形函数的个数。两种工况取lnt分别为30、40，定每个窗含的时间步数分别为420、400。移动时间窗前后重叠一半。两种工况的总识别时间步分别为8 821、8 801。

这里激励器施加的荷载是一个纯弯矩，因此，至少需要一个传感器才能准确识别荷载。分别利用单个测点S1，两个测点S1－S2，三个测点S1－S3进行识别。为避免赘述，本文仅针对工况1给出了由不同测点数目识别的结果，见图13。其中当没有截断奇异值时，由S1识别的荷载尾部发散，使荷载计算收敛的奇异值截断数目为4;利用多个测点识别时不需奇异值截断。

由图13可以看出，利用单个传感器可以识别荷载，但识别精度稍差。利用两个传感器和三个传感器识别的结果基本上是一致的，识别精度更高。采用不同传感器识别的荷载误差曲线见图14，考虑所有测量时间内的荷载识别相对误差 F计算－F实测2/F实测2分别为51.83%，32.46%，26.6%。针对工况二只给出了利用两个测点S1－S2识别的结果，见图15，未采用奇异值截断，荷载识别相对误差为22.93%。

可以看出，两种工况下利用提出的荷载形函数识别方法均能较好地识别未知荷载。识别的荷载历程曲线与实际曲线吻合较好，能够反映实际荷载信息。与多个传感器相比，利用单个传感器识别荷载，结构的动态响应信息相对较少，由于噪声和测量误差的影响，识别精度相对较低。移动时间窗在线识别荷载时，利用反卷积方法需要每个窗内重新计算L－曲线，确定奇异值截断数目，这很耗时，难于实现。而荷载形函数在一定程度上能够过滤噪声的影响，不需要或只需少数几个奇异值截断，避免计算L－曲线。这进一步验证在线荷载识别中形函数方法的重要角色。

5 结论

本文介绍并理论补充了基于形函数的荷载识别方法，通过数值算例和一个悬臂梁的试验验证了该方法的有效性，并得到以下主要结论:

(1)针对反卷积方法识别荷载中存在的测量时间长或采样频率高时数值求解困难的问题，利用有限的形函数逼近未知荷载，将识别离散的荷载历程转化为识别有限的荷载形函数系数，提高计算效率。

(2)荷载形函数是连续光滑的，在一定程度上能够过滤噪声的影响，提高逆问题对噪声的鲁棒性。

(3)通过实测结构响应的频谱分析可以确定逼近未知荷载所需的形函数个数。

(4)形函数与移动时间窗相结合，每个窗内不需或只需截断少数几个奇异值就可以有效地识别未知荷载，从而实现荷载的在线识别。

［1] 瞿伟廉，王锦文.振动结构动态荷载识别综述［J] .华中科技大学学报(城市科学版)，2004，21(4):94－99.

［2] Inoue H，Harrigan JJ，Reid SR.Review of inverse analysis for indirect measurement of impact force［J] .Applied Mechanics Reviews，2001，54(6):503 －524.

［3] Jacquelin E，Bennani A，Hamelin P.Force reconstruction:analysis and regularization of a deconvolution problem［J] .Journal of Sound and Vibration，2003，265(1):81－107.

［4] Uhl T.The inverse identification problem and its technical application［J] .Archive of Applied Mechanics， 2007，77(5):325－337.

［5] 智 浩，文祥荣，缪龙秀，等.动态载荷的频域识别方法［J] .北方交通大学学报，2004(4):5－10.

［6] Inoue H，Ishida H，Kishimoto K.et al.Measurement of impact load by using analysis technique［C] .JSME Int.J.Series 1，1991，34(4):453 －458.

［7] JankowskiŁ.Off－line identification of dynamic loads［J] .Structural and Multidisciplinary Optimization，2009，37(6):609－623.

［8] Zhang Q，JankowskiŁ，Duan Z.Identification of coexistent load and damage based on virtual distortion method［C] .In:Proceedings of the 4th European Workshop on Structural Health Monitoring，Cracow，Poland，DEStech，2008:1121－1128.

［9] Doyle J F.A wavelet deconvolution method for impact force identification［J] .Experimental mechanics ，1997，37(4):403－408.

［10] 张 方，朱德懋，张福祥.动荷载识别的时间有限元模型理论及其应用［J] .振动与冲击，1998，17(2):1－4.

［11] Ha QP，Trinh H.State and input simultaneous estimation for a class of nonlinear systems［J] .Automatica，2004，40:1779－1785.

［12] Allen M S，Carne T G.Delayed，multi－step inverse structural filter for robust force identification［J] .Mechanical Systems and Signal Processing，2008，22:1036－1054.

［13] Hansen P C.Numerical tools for analysis and solution of Fredholm integral equations of the first kind［J] .Inverse Problems 8，1992:849 －872.

［14] Zhang Q，JankowskiŁ，Duan Z.Identification of coexistent load and damage［J] .Structural and Multidisciplinary Optimization，2010，41(2):243－253.

［15] Dahlquist G，Björck Ä.Numerical methods in scientific computing，volⅡ［M] .2008，to be published by SIAM，URL http://www.mai.liu.se/～ akbjo/NMbook.html.