CN107991060B - 基于自适应和迭代算法的载荷分布式光纤辨识方法 - Google Patents

Abstract

本发明公开一种基于自适应和迭代算法的载荷分布式光纤辨识方法，属于结构健康监测领域。包括：步骤一：采用有限元方法，得到梁结构的离散化状态方程；步骤二：基于分布式光纤传感器的气动载荷‑应变响应信号采集；步骤三：基于卡尔曼滤波器和载荷估计器的气动载荷分布状态反演；步骤四：采用Sage‑Husa自适应和迭代实现对步骤三气动载荷分布反演过程中噪声特性参数Q，R和收敛特性的调节，分别得到Q，R参数优化值；步骤五：将步骤四优化所得参数Q，R代入步骤三所述算法，作为下一次采样时刻的基准参量；步骤六：按照步骤二至步骤五顺序，依次反复循环上述过程。本发明提高了收敛速度以及对动态载荷的实时估计精度，具有简单方便，实时性强等特点。

Description

基于自适应和迭代算法的载荷分布式光纤辨识方法

技术领域

本发明属于结构健康监测领域，提出了一种基于自适应和迭代算法的载荷分布式光纤辨识方法。

背景技术

在航空航天结构设计和健康状态监测中，载荷识别可以确保结构的设计安全，实时监测结构的健康状态，为航空航天器安全服役提供可靠保障。

动载荷识别的方法分为直接测量法和间接识别法两种，前者即通过传感器直接测量载荷本身或者与载荷相关的参数进而确定载荷的大小，但在大多数实际工程应用中，动态载荷无法通过直接测量获得，如火箭飞行状态中所受的推力，爆炸冲击的载荷，飞机飞行状态中机翼所受的气动载荷等。此时必须采用间接方法对动载荷进行识别，即载荷识别技术。

目前载荷识别方法中应用最广泛的是频域法，即频响函数矩阵求逆，最小二乘法和模态坐标变换法。频响函数求逆法只要获得频率响应函数矩阵及响应谱，即可识别动态载荷。最小二乘法是频响函数矩阵求逆法的补充和深入，当响应点数目多于激励点数目时，频响函数矩阵不能直接求逆，求最小二乘解获得载荷。这两种方法的缺点在于低频的识别精度较差。目前可以通过奇异值分解技术改善，还可以通过相干函数作为阈值控制矩阵的病态，但是阈值的选取尚无成熟方法，识别精度有待进一步提高。模态坐标变化法必须知道系统的固有频率及模态振型才能识别载荷特性，参与计算的模态的选取及截断对载荷识别精度影响较大，选取不同阶数的模态，测试噪声等因素对识别精度的影响较大，在工程应用中存在较大的局限性。时域法克服了频域法只能对平稳载荷识别的缺点，但传统的杜哈梅积分方法存在着积分计算量大，对边值条件和初值条件敏感，数值不稳定等问题。时域法探讨的模型多数建立在仿真模型的基础上，识别模型与实际模型存在差异，与实际应用还有一段距离，实际工程中常见的分布载荷的识别仍然比较罕见。

目前载荷识别算法在实时估计方面尚未得到充分研究，为了实时估计气动载荷，在此提出了一种新的载荷识别算法，此算法拥有良好的噪声处理能力和实时显示的特点，并可以与系统控制完美结合，具有良好的工程应用价值。在此基础上，提出利用Sage-Husa自适应和迭代对噪声特性和收敛性能进行实时调节，提高了系统的实时性和辨识精度。

由于分布式光纤光栅传感器可以实现对温度、应变等物理量的直接测量，且具有不受电磁干扰、信号带宽大、灵敏度高、易于复用、重量轻、适于在高温、腐蚀性的危险环境中使用等优点。因此提出了基于分布式光纤布拉格光栅传感器应变信息实时采集的气动载荷反演辨识方法。

发明内容

技术问题：本发明目的在于提供一种识别精度高，收敛速度快的基于自适应和迭代算法的载荷分布式光纤辨识方法。

该方法采用分布式光纤布拉格光栅传感器网络感知梁结构中不同位置的载荷响应信号，通过研究传感器测得的应变响应信号，计算气动载荷大小，利用迭代算法加快收敛性能，利用自适应算法识别噪声特性。

技术方案：为解决上述技术问题，本发明的所采用的技术方案包括以下步骤：

步骤一：Bernoulli-Euler梁结构状态方程的建立及离散化

利用有限元方法将梁结构进行离散化，得到n个有限元单元，在每个单元内粘贴两个光纤布拉格光栅传感器，这两个光纤布拉格光栅传感器依次沿梁结构表面中心线等间距布置，得到光纤布拉格光栅传感器所采集的应变值与梁结构位移转角的关系如下所示：

式中，为不同位置光纤布拉格光栅传感器采集所得应变值，w₁,w₂,…w_2n为不同单元节点的位移，θ₁,θ₂,…θ_2n为不同单元节点的转角，l为单元长度，h为梁厚度，ξ_i由光纤布拉格光栅传感器粘贴位置决定。

对于2n个自由度的梁结构振动系统，矩阵形式的微分方程如下所示：

其中M是2n×2n的质量矩阵，C是2n×2n的阻尼矩阵，K是2n×2n的刚度矩阵，F(t)是2n×1的梁结构所受气动载荷向量，Y(t),分别是梁结构振动对应的位移，速度和加速度。将上述矩阵形式的微分方程转换为状态空间方程，其中根据公式(2)，状态方程可以改写为如下形式：

以光纤布拉格光栅传感器实时采集的应变值作为观测值，根据应变值与节点位移和转角的关系，观测方程改写成如下形式：

Z(t)＝HX(t) (4)

其中：

F(t)＝[F₁,F₂,F₃,…F_2n]表示气动载荷向量。A是状态矩阵，由梁结构材料的质量矩阵，刚度矩阵和阻尼矩阵决定。H是状态方程的观测矩阵，由观测值和状态变量的关系决定。Z(t)是观测序列，由光纤布拉格光栅传感器采集的应变值确定。将采样间隔设为ΔT，将方程(3)、(4)离散化得到：

X(k+1)＝ΦX(k)+Γ(F(k)+w(k)) (5)

Z(k)＝HX(k)+v(k) (6)

Φ＝exp(A*ΔT) (7)

其中X(k)代表k时刻的状态向量，表示梁结构的位移转角的状态。Φ代表状态转移矩阵，τ为噪声积分变量，Γ代表驱动矩阵，H为状态方程的观测矩阵，F(k)为气动载荷序列。w(k)和v(k)为互不相关的高斯白噪声序列，将其噪声特性描叙为：

E[w(k)]＝0,E[w(k)w^T(m)]＝Q(k)δ_kl,Q＝Q_wI_2n*2n (9)

E[v(k)]＝0,E[v(k)v^T(m)]＝R(k)δ_kl,R＝R_vI_2n*2n,R_v＝σ² (10)

式中E为代表均值，Q和R是协方差矩阵，I_2n*2n为单位矩阵，k和m为不同的采样时刻。

步骤二：基于分布式光纤布拉格光栅应变传感器的气动载荷-应变响应信号采集；

光纤布拉格光栅应变传感器数目等于梁结构单元数目的二倍。将光纤布拉格光栅应变传感器接入光纤光栅解调仪，传输信号至上位机。通过位于不同位置的光纤布拉格光栅应变传感器采样得到梁结构所受气动载荷与传感器所在位置应变的对应关系。

步骤三：基于卡尔曼滤波器和载荷估计器的气动载荷分布状态反演

载荷识别算法利用系统响应来估计梁结构的载荷。本发明所述估计算法包含两部分：卡尔曼滤波器和载荷估计器。载荷估计器利用卡尔曼滤波器产生的增益矩阵、新息序列和协方差矩阵实时估计梁结构所受载荷，其步骤如下所示：

卡尔曼滤波器：

P(k/k-1)＝ΦP(k-1/k-1)Φ^T+ΓQΓ^T (12)

S(k)＝HP(k/k-1)H^T+R (13)

K_a(k)＝P(k/k-1)H^TS^-1(k) (14)

P(k/k)＝[I-K_a(k)H]P(k/k-1) (15)

载荷估计器：

B_s(k)＝H[ΦM_s(k-1)+I]Γ (18)

M_s(k)＝[I-K_a(k)H][ΦM_s(K-1)+I] (19)

K_b(k)＝γ^-1P_b(k-1)B^T _s(k)[B_s(k)γ^-1P_b(k-1)B^T _s(k)+S(k)]^-1

(20)

P_b(k)＝[I-K_b(k)B_s(k)]γ^-1P_b(k-1) (21)

其中I代表单位矩阵，P代表协方差矩阵，S(k)代表新息协方差，B_s(k)和M_s(k)代表灵敏度矩阵，代表新息，K_b(k)是用来估计气动载荷的增益矩阵，代表k时刻气动载荷的估计值，γ是衰减因子，用于调节自适应能力。

步骤四：采用Sage-Husa自适应和迭代实现对步骤三气动载荷分布反演过程噪声特性和收敛性能的调节，分别得到Q，R的参数优化值

由于基于分布式光纤应变信号采集的气动载荷估计算法存在延迟特性和收敛性受限制的缺陷，为了在基于动态应变采集的气动载荷突变信号等特殊情况下得到更准确的识别结果，需要对识别结果进行循环迭代修正。修正算法流程图如图2所示：

其中ΔJ用来衡量新息误差的大小，η表示测量误差允许范围，i是迭代次数，n代表最大迭代次数。当ΔJ>η且没有达到最大迭代次数时，基于分布式光纤应变数据采集的气动载荷反演模型将循环迭代，加快收敛速度。

利用Sage-Husa时变噪声统计估计器对Q(k)，R(k)进行估计，得到：

其中d_k＝(1-b)/(1-b^k+1),0＜b＜1,b是遗忘因子，是公式(9)(10)中k时刻的Q和R的估计值，K(k)为公式(14)中的K_a(k),ε(k/k-1)为公式(16)中的

步骤五：气动载荷监测系统观测噪声和量测噪声特性参数的实时更新

将步骤四优化所得前一次采样时刻对应的Q，R直接代入步骤三所述的梁结构气动载荷分布反演算法，作为下一次采样时刻基于分布式光纤布拉格光栅传感器应变分布信息采集的气动载荷反演方法的基准参量。

步骤六：基于循环计算的气动载荷实时估计

按照步骤二至步骤五顺序，依次反复循环上述过程，实现梁结构所受气动载荷的实时估计，直至测量过程结束。

本发明的优点是：

(1)传统电阻应片传感方式具有布线困难，易受电磁干扰的缺陷，分布式光纤应变传感器具有线路布置简单、质量轻、便于分布式组网、抗电磁干扰能力强等优点。基于FBG的载荷识别技术解决了载荷识别领域传感器安装困难的缺陷。

(2)常规算法实时性较差，往往具有迟后性，而本发明提出的载荷估计器能够实现每个时刻估计一次载荷，具有实时估计的优点。

(3)为了解决工程中噪声特性未知，从而影响估计效果的问题，在此提出Sage-Husa自适应算法对噪声特性进行实时估计。为了解决采样频率较低情况下的收敛速度较慢的缺陷，还提出迭代修正算法加快系统估计载荷的收敛速度。

附图说明

图1是系统装配图。其中：1 光源，2 耦合器，3 FBG，4 计算机，5 解调仪。

图2是迭代算法流程图

具体实施方式

1、一种基于Sage-Husa自适应迭代算法的梁结构气动载荷分布式光纤辨识方法。包括以下步骤：

步骤一：Bernoulli-Euler梁结构状态方程的建立及离散化

利用有限元方法将梁结构进行离散化，得到n个有限元单元，在每个单元内粘贴两个光纤布拉格光栅传感器，这两个光纤布拉格光栅传感器依次沿梁结构表面中心线等间距布置，得到光纤布拉格光栅传感器所采集的应变值与梁结构位移转角的关系如下所示：

式中，为不同位置光纤布拉格光栅传感器采集所得应变值，w₁,w₂,…w_2n为不同单元节点的位移，θ₁,θ₂,…θ_2n为不同单元节点的转角，l为单元长度，h为梁厚度，ξ_i由光纤布拉格光栅传感器粘贴位置决定；

对于2n个自由度的梁结构振动系统，矩阵形式的微分方程如下所示：

其中M是2n×2n的质量矩阵，C是2n×2n的阻尼矩阵，K是2n×2n的刚度矩阵，F(t)是2n×1的梁结构所受气动载荷向量，Y(t),分别是梁结构振动对应的位移，速度和加速度；将上述矩阵形式的微分方程转换为状态空间方程，其中根据公式(2)，状态方程可以改写为如下形式：

以光纤布拉格光栅传感器实时采集的应变值作为观测值，根据应变值与节点位移和转角的关系，观测方程改写成如下形式：

Z(t)＝HX(t) (4)

其中：

F(t)＝[F₁,F₂,F₃,…F_2n]表示气动载荷向量；A是状态矩阵，由梁结构材料的质量矩阵，刚度矩阵和阻尼矩阵决定；H是状态方程的观测矩阵，由观测值和状态变量的关系决定；Z(t)是观测序列，由光纤布拉格光栅传感器采集的应变值确定；将采样间隔设为ΔT，将方程(3)、(4)离散化得到：

X(k+1)＝ΦX(k)+Γ(F(k)+w(k)) (5)

Z(k)＝HX(k)+v(k) (6)

Φ＝exp(A*ΔT) (7)

其中X(k)代表k时刻的状态向量，表示梁结构的位移转角的状态；Φ代表状态转移矩阵，τ为噪声积分变量，Γ代表驱动矩阵，H为状态方程的观测矩阵，F(k)为k时刻载荷序列；w(k)和v(k)为互不相关的高斯白噪声序列，将其噪声特性描叙为：

E[w(k)]＝0,E[w(k)w^T(l)]＝Q(k)δ_kl,Q＝Q_wI_2n*2n (9)

E[v(k)]＝0,E[v(k)v^T(l)]＝R(k)δ_kl,R＝R_vI_2n*2n,R_v＝σ² (10)

式中Q和R是协方差矩阵，I_2n*2n为单位矩阵，k和m为不同的采样时刻；

步骤二：基于分布式光纤布拉格光栅应变传感器的气动载荷-应变响应信号采集；

光纤布拉格光栅应变传感器数目等于梁结构单元数目的二倍；将光纤布拉格光栅应变传感器接入光纤光栅解调仪，传输信号至上位机；通过位于不同位置的光纤布拉格光栅应变传感器采样得到梁结构所受气动载荷与传感器所在位置应变的对应关系。

步骤三：基于卡尔曼滤波器和载荷估计器的气动载荷分布状态反演

载荷识别算法利用系统响应来估计梁结构的载荷；本发明所述估计算法包含两部分：卡尔曼滤波器和载荷估计器；载荷估计器利用卡尔曼滤波器产生的增益矩阵、新息序列和协方差矩阵实时估计梁结构所受载荷，其步骤如下所示：

卡尔曼滤波器：

P(k/k-1)＝ΦP(k-1/k-1)Φ^T+ΓQΓ^T

(12)

S(k)＝HP(k/k-1)H^T+R (13)

K_a(k)＝P(k/k-1)H^TS^-1(k) (14)

P(k/k)＝[I-K_a(k)H]P(k/k-1) (15)

载荷估计器：

B_s(k)＝H[ΦM_s(k-1)+I]Γ (18)

M_s(k)＝[I-K_a(k)H][ΦM_s(K-1)+I] (19)

K_b(k)＝γ^-1P_b(k-1)B^T _s(k)[B_s(k)γ^-1P_b(k-1)B^T _s(k)+S(k)]^-1 (20)

P_b(k)＝[I-K_b(k)B_s(k)]γ^-1P_b(k-1) (21)

其中I代表单位矩阵，P代表协方差矩阵，S(k)代表新息协方差，B_s(k)和M_s(k)代表灵敏度矩阵，代表新息，K_b(k)是用来估计气动载荷的增益矩阵，代表k时刻气动载荷的估计值，γ是衰减因子，用于调节自适应能力。

步骤四：采用Sage-Husa自适应和迭代实现对步骤三气动载荷分布反演过程噪声特性和收敛性能的调节，分别得到Q，R的参数优化值

由于基于分布式光纤应变信号采集的气动载荷估计算法存在延迟特性和收敛性受限制的缺陷，为了在基于动态应变采集的气动载荷突变信号等特殊情况下得到更准确的识别结果，需要对识别结果进行循环迭代修正；修正算法如图2所示：

其中ΔJ用来衡量新息误差的大小，η表示测量误差允许范围，i是迭代次数，n代表最大迭代次数。当ΔJ>η且没有达到最大迭代次数时，基于分布式光纤应变数据采集的气动载荷反演模型将循环迭代，加快收敛速度。

利用Sage-Husa时变噪声统计估计器对Q(k)，R(k)进行估计，得到：

其中d_k＝(1-b)/(1-b^k+1),0＜b＜1,b是遗忘因子，是公式(9)(10)中k时刻的Q和R的估计值，K(k)为公式(14)中的K_a(k),ε(k/k-1)为公式(16)中的

步骤五：气动载荷监测系统观测噪声和量测噪声特性参数的实时更新

将步骤四优化所得前一次采样时刻对应的Q，R直接代入步骤三所述的梁结构气动载荷分布反演算法，作为下一次采样时刻基于分布式光纤布拉格光栅传感器应变分布信息采集的气动载荷反演方法的基准参量。

步骤六：基于循环计算的气动载荷实时估计

按照步骤二至步骤五顺序，依次反复循环上述过程，实现梁结构所受气动载荷的实时估计，直至测量过程结束。

Claims (1)

1.一种基于自适应和迭代算法的载荷分布式光纤辨识方法，包括以下步骤：

步骤一：Bernoulli-Euler梁结构状态方程的建立及离散化

利用有限元方法将梁结构进行离散化，得到n个有限元单元，在每个单元内粘贴两个光纤布拉格光栅传感器，这两个光纤布拉格光栅传感器依次沿梁结构表面中心线等间距布置，得到光纤布拉格光栅传感器所采集的应变值与梁结构位移转角的关系如下所示：

式中，i＝1,2…2n,ε₁…ε_2n为不同位置光纤布拉格光栅传感器采集所得应变值，w₁,w₂,…w_2n为不同单元节点的位移，θ₁,θ₂,…θ_2n为不同单元节点的转角，l为单元长度，h为梁厚度，ξ_i由光纤布拉格光栅传感器粘贴位置决定；

对于2n个自由度的梁结构振动系统，矩阵形式的微分方程如下所示：

其中M是2n×2n的质量矩阵，C是2n×2n的阻尼矩阵，K是2n×2n的刚度矩阵，F(t)是2n×1的梁结构所受气动载荷向量，Y(t),分别是梁结构振动对应的位移，速度和加速度；将上述矩阵形式的微分方程转换为状态空间方程，其中根据公式(2)，状态方程可以改写为如下形式：

以光纤布拉格光栅传感器实时采集的应变值作为观测值，根据应变值与节点位移和转角的关系，观测方程改写成如下形式：

Z(t)＝HX(t) (4)

其中：

F(t)＝[F₁,F₂,F₃,…F_2n]表示气动载荷向量；A是状态矩阵，由梁结构材料的质量矩阵，刚度矩阵和阻尼矩阵决定；H是状态方程的观测矩阵，由观测值和状态变量的关系决定；Z(t)是观测序列，由光纤布拉格光栅传感器采集的应变值确定；将采样间隔设为ΔT，将方程(3)、(4)离散化得到：

X(k+1)＝ΦX(k)+Γ(F(k)+w(k)) (5)

Z(k)＝HX(k)+v(k) (6)

Φ＝exp(A*ΔT) (7)

其中X(k)代表k时刻的状态向量，表示梁结构的位移转角的状态,Φ代表状态转移矩阵，τ为噪声积分变量，Γ代表驱动矩阵，H为状态方程的观测矩阵，F(k)为气动载荷序列,w(k)和v(k)为互不相关的高斯白噪声序列，将其噪声特性描叙为：

E[w(k)]＝0,E[w(k)w^T(m)]＝Q(k)δ_kl,Q＝Q_wI_2n*2n (9)

E[v(k)]＝0,E[v(k)v^T(m)]＝R(k)δ_kl,R＝R_vI_2n*2n,R_v＝σ² (10)

式中E为代表均值，Q和R是协方差矩阵，I_2n*2n为单位矩阵，k和m为不同的采样时刻；

步骤二：基于分布式光纤布拉格光栅应变传感器的气动载荷-应变响应信号采集

光纤布拉格光栅应变传感器数目等于梁结构单元数目的二倍,将光纤布拉格光栅应变传感器接入光纤光栅解调仪，传输信号至上位机，通过位于不同位置的光纤布拉格光栅应变传感器采样得到梁结构所受气动载荷与传感器所在位置应变的对应关系；

步骤三：基于卡尔曼滤波器和载荷估计器的气动载荷分布状态反演

载荷识别算法利用光纤布拉格光栅应变传感器测量值来估计梁结构的载荷；本方法所述估计算法包含两部分：卡尔曼滤波器和载荷估计器；载荷估计器利用卡尔曼滤波器产生的增益矩阵、新息序列和协方差矩阵实时估计梁结构所受载荷，其步骤如下所示：

卡尔曼滤波器：

P(k/k-1)＝ΦP(k-1/k-1)Φ^T+ΓQΓ^T (12)

S(k)＝HP(k/k-1)H^T+R (13)

K_a(k)＝P(k/k-1)H^TS^-1(k) (14)

P(k/k)＝[I-K_a(k)H]P(k/k-1) (15)

载荷估计器：

B_s(k)＝H[ΦM_s(k-1)+I]Γ (18)

M_s(k)＝[I-K_a(k)H][ΦM_s(K-1)+I] (19)

K_b(k)＝γ^-1P_b(k-1)B^T _s(k)[B_s(k)γ^-1P_b(k-1)B^T _s(k)+S(k)]^-1 (20)

P_b(k)＝[I-K_b(k)B_s(k)]γ^-1P_b(k-1) (21)

其中I代表单位矩阵，P代表协方差矩阵，S(k)代表新息协方差，B_s(k)和M_s(k)代表灵敏度矩阵，代表新息，K_b(k)是用来估计气动载荷的增益矩阵，代表k时刻气动载荷的估计值，γ是衰减因子，用于调节自适应能力；

步骤四：采用Sage-Husa自适应和迭代实现对步骤三气动载荷分布反演过程噪声特性和收敛性能的调节，分别得到Q，R的参数优化值

由于基于分布式光纤应变信号采集的气动载荷估计算法存在延迟特性和收敛性受限制的缺陷，为了在基于动态应变采集的气动载荷突变信号等特殊情况下得到更准确的识别结果，需要对识别结果进行循环迭代修正，修正算法根据的值进行判断；如果ΔJ在误差允许范围之内，则停止迭代，如果ΔJ在误差允许范围之外，则进行迭代修正，其中ΔJ用来衡量新息误差的大小，η表示测量误差允许范围，i是迭代次数，n代表最大迭代次数；当ΔJ>η且没有达到最大迭代次数时，基于分布式光纤应变数据采集的气动载荷反演模型将循环迭代，加快收敛速度；

利用Sage-Husa时变噪声统计估计器对Q(k)，R(k)进行估计，得到：

其中d_k＝(1-b)/(1-b^k+1),0＜b＜1,b是遗忘因子，是公式(9)(10)中k时刻的Q和R的估计值，K(k)为公式(14)中的K_a(k),ε(k/k-1)为公式(16)中的

步骤五：气动载荷监测系统观测噪声和量测噪声特性参数的实时更新

将步骤四优化所得前一次采样时刻对应的Q，R直接代入步骤三所述的梁结构气动载荷分布反演算法，作为下一次采样时刻基于分布式光纤布拉格光栅传感器应变分布信息采集的气动载荷反演方法的基准参量；

步骤六：基于循环计算的气动载荷实时估计

按照步骤二至步骤五顺序，依次反复循环上述过程，实现梁结构所受气动载荷的实时估计，直至测量过程结束。