载荷识别主要用于计算恶劣服役环境下结构主承力构件所受载荷或者无法通过传感器获取的外部激励^[1]。动载荷识别一般分为直接测量法和间接识别法2种^[2]，直接测量法是通过在飞行器表面设置大量测压孔感知整个表面的压力分布。此方式特点是测量精度高、计算量小，但缺点是打孔给飞行器结构引入损伤，一定程度上影响了飞行器结构力学性能。

间接识别算法是在结构表面布置传感器网络直接采集载荷响应信息来反求载荷，此种方法不需破坏结构、操作较简单。国内外许多学者对此进行了大量深入的研究，Hong等提出了一种倒谱域的机械动态载荷谱识别方法，基于高阶谱的谐波恢复，建立了ARMA(Auto Regvessive Moving Avevage Model)模型，确定格林函数的系数，求出动载荷，但只是进行了仿真验证，缺乏工程验证^[3]。Maes等基于杆件的模态特性，假设了Timoshenko梁理论，在考虑杆的转动惯量和传感器质量的基础上，提出了一种估算组合结构杆件轴向力的方法，但不适用于变截面等复杂结构的动载荷识^[4]。Ma等提出了一种基于卡尔曼滤波和递推最小二乘算法的在线递推动载荷方法。此算法对高斯白噪声具有很好的抑制作用，通过正弦、三角脉冲、矩形脉冲、一系列脉冲和随机脉冲仿真，验证了该方法对由噪声测量值估计梁结构外部激励具有良好的性能，但未实现载荷作用位置辨识^[5]。间接识别法属于数学上的逆问题，因而存在矩阵病态导致的识别精度不高和计算量偏大的问题^[6]。目前关于间接识别法的研究中，除了矩阵病态问题，还有传感器的安装和布局等难点^[7-8]。

由于光纤传感技术能够实现针对应变、温度、振动等多种物理量的直接测量，具有便于分布式组网、灵敏度高、芯径细、重量轻、柔韧性好、适于强电磁干扰以及腐蚀性危险环境使用等独特优点^[9]。基于此，本文研究了一种基于光纤光栅传感器与卡尔曼滤波器的变截面悬臂梁结构动载荷识别算法。该方法基于自适应卡尔曼滤波器反馈修正的思想进行载荷识别，一定程度上解决了间接识别法中矩阵求逆所导致的病态问题，给出了载荷所用位置判断依据，并通过仿真与实验验证了此算法具有较好的抗噪能力、较高的载荷辨识精度和工程适用性。

1 光纤光栅应变传感原理

光纤光栅传感器(Fiber Bragg Grating，FBG)基于光纤周期或纤芯折射率的变化，当光纤Bragg光栅感受到外界温度、应变变化时，会引起光栅栅区折射率的变化，从而使得其对应中心波长发生改变^[10]。光纤光栅传感原理，如图 1所示。

由光纤光栅传感器原理可知，假定环境温度恒定，此时FBG仅受轴向应变因素影响，当中心波长变化量为Δλ_B时，其应变表达式为

$ \varepsilon = \frac{1}{{1 - {P_{\rm{e}}}}} \cdot \frac{{\Delta {\lambda _{\rm{B}}}}}{{{\lambda _{\rm{B}}}}} $	（1）

式中：P_e为光纤光栅有效弹光系数，取为常数0.22; λ_B为FBG反射光谱中心波长; Δλ_B为中心波长偏移量^[11]。

2 动载荷识别原理 2.1 变截面梁有限元分析

针对如图 2所示长度为l、厚度为h、上底边宽度为b₀、下底边宽度为b₁的锥形变截面梁单元形式，其宽度b(x)、截面积A(x)和惯性矩I(x)可以表示为^[12]

$ \left\{ {\begin{array}{*{20}{l}} {b(x) = {b_0}\left( {1 - \frac{x}{l}} \right) + {b_1}\frac{x}{l}}\\ {A(x) = h\left[ {{b_0}\left( {1 - \frac{x}{l}} \right) + {b_1}\frac{x}{l}} \right]}\\ {I(x) = \frac{{{h^3}}}{{12}}\left[ {{b_0}\left( {1 - \frac{x}{l}} \right) + {b_1}\frac{x}{l}} \right]} \end{array}} \right. $	（2）

根据结构力学的有限元分析可知，等截面梁单元刚度矩阵K^e可以表示为

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{K}}^e} = \int\limits_0^l {{\mathit{\boldsymbol{B}}^{\rm{T}}}} EI\mathit{\boldsymbol{B}}{\rm{d}}x = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{EI}}{{{l^3}}}\left[ {\begin{array}{*{20}{c}} {12}&{6l}&{ - 12}&{6l}\\ {6l}&{4{l^2}}&{ - 6l}&{2{l^2}}\\ { - 12}&{ - 6l}&{12}&{ - 6l}\\ {6l}&{2{l^2}}&{ - 6l}&{4{l^2}} \end{array}} \right]} \end{array} $

（3）

对于等截面梁单元的一致质量矩阵M^e可以表示为

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{M}}^e} = \int_0^l {{\mathit{\boldsymbol{N}}^{\rm{T}}}} \rho A\mathit{\boldsymbol{N}}{\rm{d}}x = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\rho Al}}{{420}}\left[ {\begin{array}{*{20}{c}} {156}&{22l}&{54}&{ - 13l}\\ {22l}&{4{l^2}}&{13l}&{ - 3{l^2}}\\ {54}&{13l}&{156}&{ - 22l}\\ { - 13l}&{ - 3{l^2}}&{ - 22l}&{4{l^2}} \end{array}} \right]} \end{array} $

（4）

式中：EI为梁单元抗弯刚度；A为梁的截面面积；b为梁截面宽度^[13]；N为梁单元的形函数矩阵；ρ为材料密度；B为几何矩阵。将式(2)代入式(3)、式(4)，可以得到变截面梁单元刚度矩阵与质量矩阵。

2.2 单元节点划分

通过有限元法将梁离散成n个单元，在每个单元内选取2个应变测量点，共计2n个应变测量点。由此，可通过形函数分别建立起每个单元内所测应变与单元节点的位移、转角之间的映射关系:

$ \left[ {\begin{array}{*{20}{c}} {{\varepsilon _1}}\\ {{\varepsilon _2}}\\ \vdots \\ {{\varepsilon _{2n - 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{B}}(1,{\xi _1})}&{}&{}&{}&{}&{}\\ {\mathit{\boldsymbol{B}}(2,{\xi _2})}&{}&{}&{}&{}&{}\\ {}&{\mathit{\boldsymbol{B}}(3,{\xi _3})}&{}&{}&{}&{}\\ {}&{\mathit{\boldsymbol{B}}(4,{\xi _4})}&{}&{}&{}&{}\\ {}&{}&\mathit{\boldsymbol{O}}&{}&{}&{}\\ {}&{}&{}&{\mathit{\boldsymbol{B}}(i,{\xi _i})}&{}&{}\\ {}&{}&{}&{}&\mathit{\boldsymbol{O}}&{}\\ {}&{}&{}&{}&{}&{\mathit{\boldsymbol{B}}(2n - 1,{\xi _{2n - 1}})}\\ {}&{}&{}&{}&{}&{\mathit{\boldsymbol{B}}(2n,{\xi _{2n}})} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _1}}\\ {{\theta _1}}\\ {{\omega _2}}\\ {{\theta _2}}\\ \vdots \\ {{\omega _{n - 1}}}\\ {{\theta _{n - 1}}}\\ {{\omega _n}}\\ {{\theta _n}} \end{array}} \right] $

（5）

式中：ε₁~ε_2n为n个单元对应的2n个应变值，实验中可由光纤光栅传感器测得；ω₁~ω_n为n个单元节点的位移；θ₁~θ_n为n个单元节点的转角；B(i, ξ_i)=[-6+12ξ_i, l(-4+6ξ_i), 6-12ξ_i, l(-2+6ξ_i)]×y/l², ξ为光纤光栅传感器粘贴位置，y=h/2为梁中性轴距离其上表面的距离。

2.3 建立状态方程

离散化后的梁结构为二阶振动系统，其2n个自由度系统的振动方程为

$ \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}^{\prime \prime }}(t) + \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{X}}^\prime }(t) + \mathit{\boldsymbol{KX}}(t) = \mathit{\boldsymbol{F}}(t) $	（6）

式中：M、C和K分别为2n×2n的质量矩阵、阻尼矩阵和刚度矩阵；F(t)为2n×1的载荷向量；X(t)为系统的位移；X'(t)为系统的速度；X″(t)为系统的加速度。

对于阻尼矩阵在实模态空间内不能解耦的情况，引入状态坐标变换^[14]:

$ \mathit{\boldsymbol{Y}}(t) = [\mathit{\boldsymbol{X}}(t)\quad {\mathit{\boldsymbol{X}}^\prime }(t)] $	（7）

将式(7)代入式(6)，可把上述矩阵形式的微分方程转换为状态空间方程，表示为

$ {\mathit{\boldsymbol{Y}}^\prime }(t) = \left[ {\begin{array}{*{20}{c}} {{{\bf{0}}_{2n \times 2n}}}&{{\mathit{\boldsymbol{I}}_{2n \times 2n}}}\\ { - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{K}}}&{ - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{C}}} \end{array}} \right]\mathit{\boldsymbol{Y}}(t) + \left[ {\begin{array}{*{20}{c}} {{{\bf{0}}_{2n \times 2n}}}\\ {{\mathit{\boldsymbol{M}}^{ - 1}}} \end{array}} \right]\mathit{\boldsymbol{F}}(t) $

（8）

令$ \mathit{\boldsymbol{D}} = \left[ \begin{gathered} {{\bf{0}}_{2n \times 2n}}\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{I}}_{2n \times 2n}} \hfill \\ - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{K}}\;\;\;\;\; - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{C}} \hfill \\ \end{gathered} \right], \mathit{\boldsymbol{G}} = \left[ \begin{gathered} {{\bf{0}}_{2n \times 2n}} \hfill \\ {\mathit{\boldsymbol{M}}^{ - 1}} \hfill \\ \end{gathered} \right]$，则状态方程可以表示为

$ {\mathit{\boldsymbol{Y}}^\prime }(t) = \mathit{\boldsymbol{DY}}(t) + \mathit{\boldsymbol{GF}}(t) $	（9）

式中：F(t)=[F₁  F₂  …  F_2n]^T; D为状态矩阵，由梁的质量矩阵、刚度矩阵和阻尼矩阵决定。系统观测方程可表示为

$ \mathit{\boldsymbol{Z}}(t) = \mathit{\boldsymbol{HY}}(t) $	（10）

式中：H为观测矩阵；Z(t)为观测序列，即结构的应变响应^[15]。

$ \mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}_1}}&{}&{}&{}&{}\\ {{\mathit{\boldsymbol{B}}_2}}&{}&{}&{}&{}\\ {}&{{\mathit{\boldsymbol{B}}_3}}&{}&{}&{}\\ {}&{{\mathit{\boldsymbol{B}}_4}}&{}&{}&{}\\ {}&{}&\mathit{\boldsymbol{O}}&{}&{}\\ {}&{}&{}&\mathit{\boldsymbol{O}}&{}\\ {}&{}&{}&{}&{{\mathit{\boldsymbol{B}}_{2n - 1}}}\\ {}&{}&{}&{}&{{\mathit{\boldsymbol{B}}_{2n}}} \end{array}} \right] $

（11）

采样时间间隔设置为ΔT，将式(9)和式(10)离散化，得到线性离散系统的状态方程和观测方程：

$ {\mathit{\boldsymbol{Y}}(k + 1) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} Y}}(k) + \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}(\mathit{\boldsymbol{F}}(k) + \mathit{\boldsymbol{w}}(k))} $	（12）

$ {\mathit{\boldsymbol{Z}}(k) = \mathit{\boldsymbol{HY}}(k) + \mathit{\boldsymbol{v}}(k)} $	（13）

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = {\rm{exp}}(\mathit{\boldsymbol{D}}*\Delta T)} $	（14）

$ {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = \int_{(k - 1)\Delta T}^{k\Delta T} {{\rm{exp}}} [\mathit{\boldsymbol{D}}(k\Delta T - \tau )]\mathit{\boldsymbol{G}}{\rm{d}}\tau } $	（15）

式中：Y(k)为状态向量；Φ为状态转移矩阵(描述动态系统从时间k到时间k+1状态的转移)；Ψ为驱动矩阵；w(k)和v(k)分别为互不相干的系统白噪声和观测白噪声，其系统特性为^[16]：数学期望(系统白噪声)为0，即E[w(k)]=0；系统白噪声自适应函数为E[w(k)w^T(k)]=Q_kd_kl；数学期望(观测白噪声)为0，即E[v(k)]=0；观测白噪声自适应函数为E[v(k)v^T(k)]=R_kd_kl, d_kl为自相关函数；Q、R为协方差矩阵, Q=Q_wI_2n×2n，R=R_vI_2n×2n, Q_w和R_v为w(k)白噪声给定的协方差初始值。

2.4 卡尔曼滤波器原理

卡尔曼滤波器可分为2个部分：时间更新预测和测量更新修正^[17-18]。

时间更新：将光纤光栅传感器所测应变响应作为输入变量进行时间方程更新。

$ \mathit{\boldsymbol{\bar Y}}(k/(k - 1)) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} \bar Y}}(k - 1/(k - 1)) $	（16）

当前时刻的状态预计：利用上一时刻的最优解Y(k-1/(k-1))与状态转移矩阵Φ，得到当前时刻的状态估计值Y(k/(k-1))。

$ \mathit{\boldsymbol{P}}(k/(k - 1)) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} P}}(k - 1/(k - 1)){\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}} + \mathit{\boldsymbol{ \boldsymbol{\varPsi} Q}}{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}} $	（17）

根据上一状态的协方差P(k-1/(k-1))和系统过程中的协方差Q，求出当前状态的协方差预测值P(k/(k-1))。

测量更新：

$ \mathit{\boldsymbol{S}}(k) = \mathit{\boldsymbol{HP}}(k/(k - 1)){\mathit{\boldsymbol{H}}^{\rm{T}}} + \mathit{\boldsymbol{R}} $	（18）

利用式(17)得到当前时刻的协方差预测值P(k/(k-1))，再结合观测矩阵H，求得新息协方差S(k)。

$ {\mathit{\boldsymbol{K}}_a}(k) = \mathit{\boldsymbol{P}}(k/(k - 1)){\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{S}}^{ - 1}}(k) $	（19）

利用新息协方差S(k)以及观测矩阵H，得到卡尔曼增益矩阵K_a(k)。

$ {\mathit{\boldsymbol{P}}(k/k) = [\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_a}(k)\mathit{\boldsymbol{H}}]\mathit{\boldsymbol{P}}(k/(k - 1))} $	（20）

$ {\mathit{\boldsymbol{\bar Z}}(k) = \mathit{\boldsymbol{Z}}(k) - \mathit{\boldsymbol{H\bar Y}}(k/(k - 1))} $	（21）

利用观测方程Z(k)、观测矩阵H以及当前时间的状态估计值Y(k/(k-1))，求得当前时间的新息序列Z(k)。

$ \mathit{\boldsymbol{\bar Y}}(k/k) = \mathit{\boldsymbol{\bar Y}}(k/(k - 1)) + {\mathit{\boldsymbol{K}}_a}(k)\mathit{\boldsymbol{\bar Z}}(k) $	（22）

利用当前时间的状态估计值Y(k/(k-1))、卡尔曼增益矩阵K_a(k)以及当前时间的新息序列Z(k)，求得当前时间的最优解Y(k/k)^[19]。

2.5 载荷识别器原理

基于卡尔曼滤波器的载荷估计器是以光纤光栅传感器测得的应变响应作为观测信号，通过卡尔曼滤波器生成的增益矩阵、新息序列以及协方差矩阵，得到灵敏度矩阵和估计力的增益矩阵。在此基础上，利用广义回归模型及其最小二乘算法估计载荷特征。

基于卡尔曼滤波器的动载荷识别器方程如下所示^[20]：

$ {{\mathit{\boldsymbol{B}}_s}(k) = \mathit{\boldsymbol{H}}[\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{M}}_{\rm{s}}}(k - 1) + \mathit{\boldsymbol{I}}]\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} $	（23）

$ {{\mathit{\boldsymbol{M}}_{\rm{s}}}(k) = [\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_a}(k)\mathit{\boldsymbol{H}}][\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{M}}_{\rm{s}}}(k - 1) + \mathit{\boldsymbol{I}}]} $	（24）

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_{\rm{b}}}(k) = {\gamma ^{ - 1}}{\mathit{\boldsymbol{P}}_{\rm{b}}}(k - 1)\mathit{\boldsymbol{B}}_{\rm{s}}^{\rm{T}}(k)[{\mathit{\boldsymbol{B}}_{\rm{s}}}(k){\gamma ^{ - 1}} \cdot }\\ {{\mathit{\boldsymbol{P}}_{\rm{b}}}(k - 1)\mathit{\boldsymbol{B}}_{\rm{s}}^{\rm{T}}(k) + \mathit{\boldsymbol{S}}(k){]^{ - 1}}} \end{array} $	（25）

$ {{\mathit{\boldsymbol{P}}_{\rm{b}}}(k) = [\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_{\rm{b}}}(k){\mathit{\boldsymbol{B}}_{\rm{s}}}(k)]{\gamma ^{ - 1}}{\mathit{\boldsymbol{P}}_{\rm{b}}}(k - 1)} $	（26）

$ {\mathit{\boldsymbol{\hat F}}(k) = \mathit{\boldsymbol{\hat F}}(k - 1) + {\mathit{\boldsymbol{K}}_{\rm{b}}}(k)[\mathit{\boldsymbol{\bar Z}}(k) - {\mathit{\boldsymbol{B}}_{\rm{s}}}(k)\mathit{\boldsymbol{\hat F}}(k - 1)]} $	（27）

式中：B_s(k)和M_s(k)代表灵敏度矩阵；K_b(k)是用来估计力的增益矩阵；$ \mathit{\boldsymbol{\hat F}}\left( k \right)$为此时刻的激振力的估计值；γ是衰减因子用来调节自适应能力。

基于光纤光栅传感器和卡尔曼滤波器的动载荷识别流程，如图 3所示。

3 有限元数值仿真 3.1 变截面悬臂梁结构建模

利用Ansys Workbench对锥形梁单元形式的变截面梁结构进行有限元数值仿真，材料密度为2 770 kg/m³，弹性模量为71 GPa，泊松比为0.3，上底边长为100 mm，下底边长250 mm，模型长度为1 200 mm，厚度为6 mm。

将变截面悬臂梁结构模型平均划分为3个单元，从左到右依次定义相应3个单元节点：A~ C。在每个单元内的1/3处和2/3处设置应变监测提取点，则共有6个应变提取点，从固支端开始依次将其定义为1~6。单元节点划分与应变提取点位置，如图 4所示，其中矩形标记为单元内应变提取点，圆形标记为单元节点。

为验证此载荷辨识算法的抗噪能力，在提取应变响应后加上σ=1×10^-9，Q_w=1×10^-15的高斯白噪声作为系统的观测序列代入相关算法。设置衰减因子γ=0.65，采样频率为500 Hz。

假设变截面悬臂梁结构系统阻尼为比例阻尼，阻尼矩阵是质量矩阵和刚度矩阵的线性组合，可以表示为

$ \mathit{\boldsymbol{C}} = \alpha \mathit{\boldsymbol{M}} + \beta \mathit{\boldsymbol{K}} $	（28）

式中:α为质量比例常数，β为刚度比例常数，设置α=0.01，β=0.02。

为进一步描述动载荷识别效果，量化动载荷识别误差，采用平均绝对误差MAE和均方根误差RMSE来描述动载荷识别误差。

$ {{\rm{MAE}} = \frac{1}{N} \cdot \sum\limits_{i = 1}^N | {\mathit{\boldsymbol{F}}_i} - {{\mathit{\boldsymbol{\hat F}}}_i}|} $	（29）

$ {{\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^N {{{({\mathit{\boldsymbol{F}}_i} - {{\mathit{\boldsymbol{\hat F}}}_i})}^2}} }}{{N - 1}}} } $	（30）

式中：F_i为在梁结构上施加的模拟真实激励；${\mathit{\boldsymbol{\hat F}}_i} $为算法识别的载荷；N为样本总数。

3.2 有限元仿真结果验证 3.2.1 单点加载下动载荷辨识效果

1) 动载荷大小识别

依次在单元节点A、B和C施加正弦激励：F=5sin(3πt)，识别效果如图 5所示。

根据图 5可知，单点载荷作用下，单元节点识别载荷效果良好，载荷识别曲线与仿真施加载荷曲线基本吻合，各单元节点载荷识别误差相差不大。单元节点A处的动载荷识别精度较好，单元节点B和单元节点C处的载荷识别曲线在峰值处略微偏小，这主要是因为悬臂梁模型在外加载荷激励过程中，距离固支端越远的测点，提取的应变响应受到外部激励扰动越明显，使得载荷辨识精度相对于固支端有所降低。

单点激励下各单元节点载荷识别误差，如表 1所示。

2) 动载荷激励位置识别

根据式(5)函数关系，本文提出通过提取变截面模型每个单元内2个位置的应变响应来计算相应单元节点处的位移和转角作为载荷识别算法的输入变量，而输出值是每个单元节点对应的载荷与弯矩。因此，算法识别出来的载荷与其相应单元节点呈一一对应关系，也即在载荷激励位置，将会存在显著的载荷响应，据此可以判断出动载荷施加位置。

以在单元节点B处施加激励载荷为例，单元节点B识别的有效载荷和单元节点A、单元节点C识别到的扰动载荷响应曲线，如图 6所示。可以发现，在单元节点B(距离固支端0.8 m)处识别的载荷显著大于单元节点A(距离固支端0.4 m)、单元节点C(距离固支端1.2 m)，且单元节点A、单元节点C处识别的载荷波动很小，无明显变化规律，因此可以判断载荷施加位置为单元节点B处，这与模拟仿真所加载荷位置相符。

3.2.2 三单元节点同时加载下动载荷辨识效果

为进一步验证此算法在复合工况下动载荷识别效果，选择在单元节点A、单元节点B和单元节点C处同时施加3个正弦激励，具体加载情况如表 2所示。

3个单元节点同时加载条件下动载荷识别效果，如图 7所示。从图 7可以看出，3个单元节点同时加载工况下，单元节点A、B、C处对应的载荷识别曲线与有限元仿真结果基本吻合，仅在峰值处均与存在一定偏差。单元节点C处识别的载荷受其他激励载荷影响，识别曲线出现一定的波动。相对于单点激励工况，由于加载幅值的增加，3个单元节点的载荷识别误差略有增大。

3个单元节点同时加载条件下动载荷识别误差，如表 3所示。

4 实验验证 4.1 变截面悬臂梁结构动载荷识别实验系统

变截面悬臂梁结构单元节点划分与传感器布局示意图，如图 8所示。变截面悬臂梁结构分布式光纤动态载荷识别实验系统，如图 9所示。

实验平台主要由光纤光栅传感器、光纤光栅解调仪、计算机、33500B型信号发生器、HEAS-50型功率放大器、力显示器、CZL-1012型力传感器、NI USB-6341数据采集卡、HEV-50型激振器等组成。

4.2 单点加载工况下动载荷辨识效果

在单元节点C和单元节点B进行单点激励，通过算法识别出来的动载荷与力传感器直接测得的真实激励进行数据分析与对比，激励工况及动载荷识别效果如下。

4.2.1 正弦激励

单元节点B正弦激励(7 Hz)动载荷识别效果，如图 10所示。由图 10(b)可以判别载荷激励位置为单元节点B。

4.2.2 方波激励

单元节点C处方波激励(1 Hz)动载荷识别效果，如图 11所示。由图 11(b)可以判别载荷激励位置为单元节点C。

4.2.3 锯齿波激励

单元节点C处锯齿波激励(1 Hz)动载荷识别效果，如图 12所示。由图 12(b)可以判别载荷激励位置为单元节点C。

根据图 10~图 12可知，在单点加载工况下，载荷识别曲线与真实激励曲线吻合度较好，只是在最大幅值处出现一定偏差。由于在施加方波、锯齿波时激振器会产生一定的顿挫，在顿挫瞬间，激励本身出现振荡，使得载荷识别曲线相较于实际施加载荷曲线存在一定波动。

通过各单元节点载荷识别曲线，根据每个单元节点对应的识别载荷可以判断出真实激励单元节点。

单点激励下，不同频率与波形激励载荷对应的动载荷辨识误差，如表 4所示。

4.3 两单元节点同时加载工况下动载荷辨识效果

在单元节点B(激励频率2 Hz)和单元节点C(激励频率4 Hz)同时施加外部激励荷载，载荷识别效果如图 13所示。

根据图 13所示，在两点同时激励的复合动载荷加载工况中，载荷辨识曲线与实际加载曲线轮廓基本吻合，算法能够较好地识别出单元节点B和单元节点C处对应的外部激励载荷特征，具有良好的识别效果。根据图 14可以判断出单元节点载荷激励位置，分别为单元节点B和单元节点C。

动载荷识别结果中呈现2个激励频率波形相互叠加的情况，这是因为在单元节点B和单元节点C处同时激励的过程中，2个激振器通过变截面悬臂梁结构相互作用，激励本身出现耦合所致。这表明该算法针对复合载荷同时加载工况，具有较好的动载荷识别能力。两单元节点同时激励实验误差，如表 5所示。

综合上述分析，动载荷识别结果也存在一定偏差，分析其原因可能包括：一是系统结构离散的单元数较少；二是光纤光栅传感器与变截面悬臂梁结构上表面未均匀胶结或粘贴方向有偏差，引起应变测量误差；三是光纤光栅传感器中心波长偏移量受结构振动扰动影响较大，尤其是在最大振幅时识别精度有所下降。

5 结论

本文针对飞行器结构载荷监测需求，推导了变截面梁模型质量矩阵、刚度矩阵等参数信息，提出了基于分布式光纤传感器与卡尔曼滤波器的变截面悬臂梁动载荷辨识方法。

1) 基于自适应卡尔曼滤波器反馈修正的思想进行载荷识别，一定程度上解决了常规间接识别法中矩阵求逆所导致的病态问题。

2) 借助ANSYS有限元分析软件，开展针对结构在不同载荷激励状态下的数值仿真，验证了该算法对于多种不同载荷形式辨识的可行性以及良好的噪声抑制能力。

3) 基于形函数建立的变截面悬臂梁单元所测应变与单元节点的位移、转角之间映射关系，提出了动载荷激励位置判别方法。

4) 构建了基于分布式光纤传感器的变截面悬臂梁应变监测与动载荷识别系统，试验结果表明，该算法对于复合载荷工况具有较好的识别能力和工程适用性。

5) 本文所述方法采用全分布式光纤传感网络，相较于电学式测量方法，具有抗电磁干扰、质量轻不改变被测结构动力学特征、适用于狭小空间机载环境的非视觉测量等特点，能够为风洞和飞行试验中飞行器载荷辨识提供技术支撑。