动载荷识别是结构动力学研究中的关键问题之一，其原理是根据已知的系统特性和动响应得到结构所受的动载荷。动载荷识别技术主要发展了频域和时域两类识别方法，其中时域方法可以更直观地表现载荷随时间历程的变化规律，更适用于工程实践。近数十年来，国内外专家学者针对时域动载荷识别的精度问题提出了许多有效的算法和研究，其中，比较典型的方法包括：逆卷积方法^[1]、逆系统方法^[2]、函数拟合法^[3-4]、 模态滤波器方法^[5]、SWAT方法^[6]及一些交叉学科方法，如神经网络^[7]和遗传算法^[8]等。对于所识别的载荷形式而言，也从简单的单点输入载荷扩展到了多源载荷^[9]、随机载荷^[10]、移动载荷^[11]，分布载荷 ^[12]以及高频载荷^[13]等。

本文为了解决由于累积误差造成识别结果不稳定甚至发散的情况，提出了基于数值迭代的改进算法，该算法能够有效地提高识别精度，避免识别发散问题，同时通过典型力学模型的仿真算例证实了识别效果。文中较为系统深入地研究了该动载荷识别算法的若干影响因素，包括区间放大系数、区间分割系数、精确度指标。通过大量算例和函数拟合，定性定量地描述了上述参数对与识别效率和精度的影响，为使用数值修正识别算法提供参考。

1 载荷识别的数值修正算法

对于已知载荷作用位置的情形，以往一步到位得到结果的载荷识别算法难免会遇到误差累积的问题。而数值修正算法在每一个时间步长内部迭代修正，可以有效地减小由于累积误差导致的发散，得到收敛的识别结果，并且有着理想的抗噪性能，具备一定的工程实用价值。

1.1 载荷初值的获得

使用数值迭代算法的前提之一是要获得可靠的初值范围，这个范围必须包含真实值且大小合适。计算初值可由一些已有的载荷识别算法获得，如正交多项式方法^[14]、Wilson-θ反分析法^[15]等。本文介绍了一种拟静态载荷法，其计算过程简单，所得结果不发散。

加速度信号是载荷识别计算中常见的输入形式，Wilson-θ法在计算由加载力得到加速度响应时具有良好的数值稳定性，取θ>1.37即可保证结果无条件收敛，本文取θ=1.4。

对于离散多自由度的系统，其动力学平衡方程为

$M\ddot x\left( t \right) + C\dot x\left( t \right) + Kx\left( t \right) = f\left( t \right)$

(1)

在Wilson-θ法^[16]的计算中，假设外力在t至t+θΔt时间内是线性变化的，于是可以得到t+θΔt时刻的系统等效静力方程

$\hat Kx\left( {t + \theta \Delta t} \right) = \hat f\left( {t + \theta \Delta t} \right)$

(2)

式中：$\hat K$为等效刚度；$\hat f$为等效载荷。在已知作用有外力的自由度上施加单位结点载荷${F^i}\left( {t + \theta \Delta t} \right)$,则有

${f^i}\left( {t + \theta \Delta t} \right) = \lambda _{t + \theta \Delta t}^i{F^i}\left( {t + \theta \Delta t} \right)i = 1,2,\cdots ,n$

(3)

式中：n为结构上作用有载荷的自由度数；${f^i}\left( {t + \theta \Delta t} \right)$时刻第i个自由度处施加的动载荷；λ_t+θΔtⁱ为载荷系数。由于${F^i}\left( {t + \theta \Delta t} \right)$是已知的单位载荷，所以求动载荷${f^i}\left( {t + \theta \Delta t} \right)$的问题就转变为求解载荷系数λ_t+θΔtⁱ。

引入位移向量dⁱ(t+θΔt)，其表示在(t+θΔt)时刻，当第i自由度处单独施加的单位节点力ⁱ(t+θΔt)时，各自由度所对应的位移，于是有

$\hat K{d^i}\left( {t + \theta \Delta t} \right) = {F^i}\left( {t + \theta \Delta t} \right)$

(4)

解dⁱ(t+θΔt)，可得线性关系式

$x\left( {t + \theta \Delta t} \right) = \sum\limits_{i = 1}^n {\left( {\lambda _{t + \theta \Delta t}^i{d^i}\left( {t + \theta \Delta t} \right)} \right)} $

(5)

设有m个测点，若已经测出其加速度响应，根据Wilson-θ法的假设，有

$\begin{gathered} x\left( {t + \theta \Delta t} \right) = x\left( t \right) + \theta \Delta t\dot x\left( t \right) + \hfill \\ \frac{{{{\left( {\theta \Delta t} \right)}^2}}}{6}\left( {\ddot x\left( {t + \theta \Delta t} \right) + 2\ddot x\left( t \right)} \right) \hfill \\ \end{gathered} $

(6)

式中$\ddot x\left( t \right),\dot x\left( t \right),x\left( t \right)$为得到的已知值。$\ddot x\left( {t + \Delta t} \right)$中m个测点加速度已知，故由式(6)可求出这m个已知测点在t+θΔt时刻的位移x^j(t+θΔt)，j＝1，2，…，m，由式(5)可得

$\begin{gathered} {x^j}\left( {t + \theta \Delta t} \right) = \sum\limits_{i = 1}^n {\left( {\lambda _{t + \theta \Delta t}^i{d^i}{{\left( {t + \theta \Delta t} \right)}^i}} \right)} \hfill \\ j = 1,2,\cdots ,m \hfill \\ \end{gathered} $

(7)

式(7)实际是含有n个未知数λ_t+θΔtⁱ(i＝1，2，…，n)，m个方程组成的线性方程组，该方程组的解可分为以下3种情况讨论：

(1) 当m≤n时，即测点数小于未知载荷数，不能得到唯一解，无法识别；

(2) 当m＝n时，可以求得唯一解；

(3) 当m≥n时，对位移矩阵求广义逆，可得到一组最小二乘解。

解出｛λ｝_t+θΔt后根据式(3)，可以求出t+θΔt时刻fⁱ(t+θΔt)的大小，进而计算下一时间点的动载荷值

${f^i}\left( {t + \Delta t} \right) = \left( {{f^i}\left( {t + \theta \Delta t} \right) - {f^i}\left( t \right)} \right)/\theta + {f^i}\left( t \right)$

(8)

1.2 利用数值迭代算法进行修正

拟静态法计算简便，解得载荷的计算结果不会出现发散。但该方法将每个时间间隔点独立计算，没有考虑前后时间点之间载荷与响应的相互影响，所以这种方法对于动载荷识别问题并不适用。为进一步得到比较精确的结果，可以用数值迭代修正拟静态法求得的载荷初值。

由Wilson-θ法和上述推导，可以将$\ddot x$(t+Δt)的表达化简为

$\ddot x\left( {t + \Delta t} \right) = \alpha {\hat K^{ - 1}}f\left( {t + \Delta t} \right) + c\left( t \right)$

(9)

式中：α为常数；c(t)为前一时间点已知量。若在第b₁，…，b_i，…，b_n自由度上施加载荷，在第a₁，…，a_i，…，a_m自由度上测得响应加速度，则有

$\begin{gathered} f\left( {t + \Delta t} \right) = {\left[\begin{gathered} 0,\cdots ,0,f{\left( {t + \Delta t} \right)_{b1}},0,\cdots ,0,\hfill \\ f{\left( {t + \Delta t} \right)_{bj}},0,\cdots ,0,f{\left( {t + \Delta t} \right)_{bn}},0,\cdots ,0 \hfill \\ \end{gathered} \right]^T} = \hfill \\ \sum\limits_{i = 1}^n {f{{\left( {t + \Delta t} \right)}_{bj}} \cdot {e_{bj}}} \hfill \\ \end{gathered} $

(10)

于是

$\ddot x\left( {t + \Delta t} \right) = \alpha {\hat K^{ - 1}} \cdot \left( {\sum\limits_{i = 1}^n {f{{\left( {t + \Delta t} \right)}_{bj}} \cdot {e_{bj}}} } \right) + c\left( t \right)$

(11)

在该方程组中取出n个方程，设取前n个方程，可得到如下方程式

$\left| \begin{gathered} {{\ddot x}_1} - {c_1} \hfill \\ {{\ddot x}_2} - {c_2} \hfill \\ {\text{ }} \vdots \hfill \\ {{\ddot x}_n} - {c_n} \hfill \\ \end{gathered} \right| = \left| \begin{gathered} {a_{11}}{a_{12}} \cdots {a_{1n}} \hfill \\ {a_{21}}{a_{22}} \cdots {a_{2n}} \hfill \\ {\text{ }} \vdots {\text{ }} \vdots {\text{ }} \ddots {\text{ }} \vdots \hfill \\ {a_{n1}}{a_{n2}} \cdots {a_{nn}} \hfill \\ \end{gathered} \right|\left| \begin{gathered} {f_1} \hfill \\ {f_2} \hfill \\ \vdots \hfill \\ {f_n} \hfill \\ \end{gathered} \right|$

(12)

式中：${\ddot x_i} = \ddot x{\left( {t + \Delta t} \right)_{ai}};{c_i} = c{\left( t \right)_{ai}};{f_i} = f{\left( {t + \Delta t} \right)_{bj}}$。求解方程(12)，消去未知量，只保留f₁得

$\sum\limits_{i = 1}^m {{A_{i1}} \cdot {{\ddot x}_i} = \left| A \right| \cdot {f_1} + } \sum\limits_{i = 1}^m {{A_{i1}} \cdot {c_i}} $

(13)

式中A_i1为系数矩阵去掉第i行第1列的代数余子式。而其余未知量用f₁表示，则有

${f_j} = \frac{{ - \sum\limits_{i = 1}^n {{A_{in,1j}} \cdot {x_i} + {A_{nj}} \cdot {f_1}} }}{{{A_{n1}}}}j = 2,3,\cdots ,n$

(14)

式中A_in，1j为系数矩阵去掉第i，n行第1,j列的代数余子式。令

$g\left( {{f_1}} \right) = \sum\limits_{i = 1}^m {{A_{i1}} \cdot {{\ddot x}_i} - } \sum\limits_{i = 1}^m {{A_{i1}} \cdot \ddot x_i^0} $

(15)

式中$\ddot x_i^0$表示真实加速度值，定义ｇ(f₁)用以衡量计算载荷所得加速度的误差，对其求导得

$g'\left( {{f_1}} \right) = \left| A \right|$

(16)

A的行列式值为零的结点，无论载荷如何变化其响应不变，本文中不考虑这种情况。那么式(16)中导数的值是常数，所以g(f₁)是一个单调函数。

假设由拟静态法得到识别载荷的初始时间序列f(t)。对于时间点t_m，设其对应的力为f(t_m)，令a₀=-r·f(t_m)，b₀=r·f(t_m)，其中，r为区间放大倍数，且能满足g(a₀)·g(b₀) < 0，则[a₀,b₀]为初始含根区间(假设f(t_m>0)。根据单调函数性质可知，在此区间内，g(f₁)的零点即为真实的载荷。一般地，如果已计算得到含根区间，则令

$f_1^{k + 1} = \left( {1 - q} \right) \cdot {a_k} + q \cdot {b_k}{\text{ }}0{\text{ < }}q{\text{ < }}1$

(17)

式中q为区间分割系数，对于不同的数值迭代算法，取各自对应的q值。若

$ g\left( {{a_k}} \right)g\left( {f_1^{k + 1}} \right)\left\{ \begin{gathered} {\text{ > }}0{\text{ }}取\left[{{a_{k + 1}},{b_{k + 1}}} \right] = \left[{f_1^{k + 1}{b_k}} \right] \hfill \\ {\text{ < }}0{\text{ }}取\left[{{a_{k + 1}},{b_{k + 1}}} \right] = \left[{{a_k}f_1^{k + 1}} \right] \hfill \\ = 0{\text{ }}则f_1^{ * 0} = f_1^{k + 1},终止计算 \hfill \\ \end{gathered} \right. $

实际计算中，所用的终止原则为b_k-a_k≤ε(一般取ε=10^-p)，其中p为正整数。

2 计算参数对算法的影响

由以上分析可以看到，这种数值修正算法增加的计算量主要取决于迭代次数的多少，所以计算耗时和结果的精确性与区间放大倍数r、区间分割系数q以及精确度指标ε有很大关系。

以无噪声情况下的两端简支矩形梁为例，该模型长度为1m，截面尺寸为0.08m×0.05 m，各阶阻尼比为0.02，弹性模量为7.2×10¹⁰，材料密度为2700kg/m³，划分为12个有限元单元，则该平面结构具有24个自由度，在第16和18自由度上分别加载f₁=10sin(4πt)N，f₂=15sin(8πt)N的正弦力，已知数据为第8，12和20自由度上的加速度响应。

取时间步长Δt=0.001s，计算时间t=1s，令r=100,q=0.5,ε=10^-7，对于上述计算模型，其拟静态初值和修正后的识别力结果如图 1(a，b)所示。

2.1 计算参数1：区间放大倍数r

区间放大倍数是使用数值迭代算法的一个重要选取参数，若选取值太小，有可能出现初始计算区间不包含真实值的情况，从而导致结果不准确甚至发散，如图 2所示。但若是r选取过大，则会带来计算量的增加，表 1分析了不同r值选取时对于迭代步数以及误差的影响，其中第314时间点为第18自由度识别力曲线峰值点(取q=0.5,ε=10^-7)。

表 1说明在其余参数选取合适时，区间放大倍数对计算精度影响并不大。误差小于0.5%在实际运用中可作为真实值的近似解，r值的影响主要体现在对于计算步数的增加上。为进一步探讨二者关系，以第314时间点为例，设其计算步数为K^*，取r=60,70,…,190,200,300,…,900,1000，分别得出对应的K^*值，经过线性插值近似去除取整和函数拟合处理后得K^*与r大致成如下对数关系

${K^ * } = \left[ {a + b \cdot \lg r} \right]$

(18)

式中[]为取整符号。经拟合计算，a=30.28,b=3.32，见图 3。

2.2 计算参数2：区间分割系数q

区间分割系数的选取也是使用数值迭代算法的重要步骤之一，一般而言该系数的选取范围为q∈0,1。不同的数值迭代算法其q值的选取也不一样，若使用二分法，则取q=0.5，若使用黄金分割法，则取q=0.618，表 2分析了不同分割系数对计算效率与精度的影响(取r=100,ε=10^-7)。

从表 2数据来看，在其余参数合适时，区间分割系数对计算精度影响也不大，误差依然小于0.5%，但q值对于计算效率的影响很大。同上以第314时间点为例，设其计算步数为K^*，取q=0.1,0.2,…,0.9，分别得出对应的K^*值，二者大致关系如图 4所示。由计算数据可知K^*与q成复杂的非线性对数关系，其具体函数关系有待进一步研究论证。

2.3 计算参数3：精确度指标ε

精确度指标的选取会直接决定计算结果的可靠度，不合适的ε值可能会导致结果严重发散(如ε=10^-2)。表 3和图 5列举了一些不同的精确度指标对于计算的影响(取r=100,q=0.5)。

从图 5和表 3中可以看出，精度指标对于计算效率和准度的影响都很大。为方便表述，使用p来表示精确度指标(ε=10^-p，即p=-igε)，以第314时间点为例，设其计算步数为K^*，误差为E，取p=3,4,…,10，分别得出对应的K^*值和E值，经过去除取整处理和函数拟合得出K^*与p大致成式(19)线性关系，如图 6所示。

${K^ * } = \left[ {a + b \cdot p} \right]a = 13.29,b = 3.32$

(19)

而E与p大致成如式(20)的负指数关系，如图 7所示，在p不断变大时，E逐渐趋于零。

$E = a \cdot {e^{bp}}{\text{ }}a = 207699.66,b = - 2.14$

(20)

3 结论

本文阐述了一种结合数值修正思想的载荷识别算法，以多输入简支梁模型为算例，分析了区间放大倍数、区间分割系数以及精确度指标3个计算参数对数值修正算法效率和精度的影响，总结出了如下结论：

(1) 数值修正算法在每一个时间步长内利用数值算法进行迭代修正，可以有效减少载荷识别过程中误差累积导致的发散。在参数选取合适时，该算法可以得到稳定的识别结果，基本收敛于真实加载。

(2) 区间放大系数选取不当有可能引起结果发散，在初始区间包含真实值的情况下，该系数对于计算结果的精确性影响不大，对于计算步数有一定影响，但不如其余参数作用大。

(3) 改变区间分割系数对于结果精度的影响也不大，但是对于计算效率的影响很大，两者具体数值关系有待进一步探索研究。

(4) 精确度指标对于计算效率和精度的影响在这3个参数中最大，不合适的精确度指标也可能导致结果发散。该参数与计算步数成线性关系，与误差成负指数关系，在实际计算中，该参数的选择十分重要。