CN104713673A - 基于扩展卡尔曼滤波器的拉索时变索力历程识别方法 - Google Patents

Abstract

一种基于扩展卡尔曼滤波器的拉索时变索力历程识别方法，如下：利用建立的考虑拉索端部位移和抗弯刚度的拉索运动方程，选择若干阶拉索振动的控制模态，将拉索的横向振动离散为振型函数和广义坐标的形式，利用包含拉索索力的扩展状态变量，将拉索振动微分方程转化为状态空间方程，基于拉索加速度的观测方程和实际桥梁拉索的监测加速度时程曲线，采用扩展卡尔曼滤波器的预测过程和更新过程，实时识别拉索时变索力历程，同时本发明可以实现在已知或未知外界激励(风荷载和其他环境荷载)监测数据的情况下时变索力历程的辨识。本方法能够对时变索力历程进行准确的实时识别，尤其适用于拉索的在线评估。

Description

基于扩展卡尔曼滤波器的拉索时变索力历程识别方法

技术领域

本发明涉及一种基于扩展卡尔曼滤波器的拉索时变索力历程识别方法。

背景技术

由于拉索具有承载能力强、施工建造方便、造型美观等优点，被广泛应用于大跨度桥梁中。据不完全统计，截至目前为止，我国已建成斜拉桥100余座、跨径50m以上的钢管混凝土拱桥230余座，特别是近十年来，我国大跨度桥梁建设技术已进入世界先进行列，建桥速度之快、数量之多、跨越能力之大为世人所瞩目，尤其是近年来建成及在建的一批跨海、跨江大桥，如昂船洲大桥、杭州湾大桥、东海大桥、苏通大桥、巫山长江大桥、重庆菜园坝长江大桥等的主通航孔均为斜拉桥或钢管混凝土拱桥。这些结构的安全问题已经成为一个国家(地区)重大政治、经济、社会问题，一旦出现破坏，不仅会对桥梁本身造成严重的经济损失，交通的中断更将导致灾难性的损失和影响。

拉索作为大跨度桥梁(斜拉桥、吊杆拱桥等)的主要承重构件，在长达几十年的使用期内，在环境侵蚀、材料老化和荷载的长期效应、疲劳效应与突变效应等灾害因素的耦合作用下，将不可避免地导致结构、系统的损伤累积和抗力衰减，从而使其抵抗自然灾害、甚至正常荷载作用的能力下降，极端情况下更会引发灾难性的突发事故。据不完全统计，20世纪70年代至90年代初，我国修建的30余座斜拉桥中，已经加固修复的桥占65％，有4座斜拉桥已拆除或改用其他桥型，有35％的斜拉桥已全部或部分更换了斜拉索，最近几年内尚有10余座90年代后修建的斜拉桥需要换索，而对于钢管混凝土拱桥，由于拉索(吊杆)断裂造成的事故不下10余起。因此，为了保障结构的安全性、完整性、适用性与耐久性，已建成使用的许多重大工程结构和基础设施亟需采用有效的手段监测和评定其安全状况、修复和控制损伤。结构健康监测是当前土木工程领域的研究热点课题之一，它对土木工程结构、特别是大型和超大型结构的设计、建造、维护和运营安全具有重要的意义，为揭示结构真实服役环境、荷载、响应和性能演化规律提供了现场试验手段。近年来，桥梁结构健康监测技术在世界范围内得到快速的发展和应用，我国在许多大跨度桥梁上安装了包括多种和较大规模传感器的健康监测系统。这些健康监测系统运行积累的大量数据，为开展结构健康监测基础科学问题的研究奠定了基础。

目前健康监测系统中存在的索力监测方法主要有两类：一类是利用索力监测装置直接监测实时时变索力，如压力传感器、磁通量传感器、光纤光栅智能拉索等监测装置，上述监测装置已经集成到某些新建桥梁的结构健康监测系统中；另一类是分析加速度信号的频率成分，通过解析求解、近似求解或经验拟合等方法建立频率与索力之间的关系，然后进行索力识别。

各类索力监测方法的原理和优缺点列于下表：

表1拉索索力监测方法及其优缺点

对于第一类索力监测方法，索力监测传感器(压力环、磁通量传感器、光纤光栅智能拉索等)及其数据采集系统一般价格昂贵、安装复杂(只能用于新建桥梁)，更主要的是目前该类传感器的耐久性不能满足健康监测系统长期监测的要求，索力监测传感器的更换耗时耗力，这些固有缺点限制了这类监测方式的大规模应用；而对于第二类索力监测方法，目前基于振动频率量测法的索力识别均通过建立频率与索力之间的关系来实现，而对于拉索振动频率的辨识需要一段时间的加速度监测数据，因此辨识得到的拉索振动频率和相应索力是某种意义上的平均值，无法识别拉索的时变索力时程。因而亟需发展一种省时、省力、经济的实时索力历程监测方法。

发明内容

针对以上问题，本发明公开一种基于扩展卡尔曼滤波器的拉索时变索力历程识别方法，采用监测加速度信息和拉索振动的状态空间系统方程，通过预测与更新两个过程辨识时变索力历程，索力辨识精度高、时效性强、方法简单易用，能够实现在线实时辨识，具有较高的鲁棒性和可靠性。

本发明所采用的技术如下：一种基于扩展卡尔曼滤波器的拉索时变索力历程识别方法，包括以下步骤：

步骤1：传感器安装及控制模态辨识；在拉索的平面内或平面外布设一个或多个加速度传感器，针对加速度传感器一个通道的一段加速度时程曲线，计算其功率谱，辨识拉索基频以及其他频率成分，利用功率谱的幅值辨识拉索振动的控制模态；

步骤2：建立拉索状态空间系统方程和观测方程；根据拉索的物理特性：拉索长度、单位长度质量、恒荷载索力，拉索振动控制模态和加速度传感器安装信息，建立拉索状态空间系统方程和观测方法；

步骤3：在已知监测风荷载输入和加速度响应输出的情况下，根据k-1步的最优估计值和系统方程计算第k步状态变量初步估计值，初步估计状态变量的协方差估计值通过系统方程的线性化得到，采用观测的拉索加速度响应，计算观测变量的初步预测误差，根据卡尔曼增益更新第k步状态变量及其协方差的最优估计值；

若拉索外界激励信息未知，将风荷载等外界激励视为系统噪声的一部分，改写拉索振动状态空间系统方程，重复步骤1-3，同样采用扩展卡尔曼滤波器的预测过程和更新过程识别拉索时变索力。

本发明还具有如下技术特征：通过监测加速度信号和拉索振动方程构建状态空间方程，通过离散扩展卡尔曼滤波方法识别拉索时变索力，具体过程如下：

忽略索的垂度和轴向惯性影响，拉索的振动方程表示为：

$\frac{d^2{\tilde{q}}_n}{d{\mathrm{\τ}}^2}+2{\mathrm{\ξ}}_{n}n\frac{d{\tilde{q}}_n}{\mathrm{d\τ}}+n^2[1+\tilde{u}\left(\mathrm{\τ}\right)]{\tilde{q}}_n+\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}{\widetilde{\mathrm{\β}}}_{\mathrm{nk}}{\tilde{q}}_n{\tilde{q}}_k^2={\tilde{U}}_n\left(\mathrm{\τ}\right)---\left(1\right)$

式中：q_n为第n阶模态的广义位移；ξ_n为拉索第n阶模态阻尼比；u₀为拉索在初始索力T₀作用下的伸长量，E为拉索的弹性模量，A为拉索横截面积；L为拉索无应力长度；τ＝ω₁t，t为时间，τ为无量纲时间，ω₁为拉索的第一阶自振频率，μ为拉索单位长度质量；u为拉索支座沿轴向运动位移；为非线性振动系数；为无量纲外荷载，其中U_n为拉索上的外荷载，f_w(x，t)为作用在索上的横向分布时变荷载；m为拉索模态质量，即

上述模态坐标下的无量纲运动方程描述拉索沿轴向存在变形时的振动情况，公式(1)变为：

$m\frac{d^2{q}_n}{d{t}^2}+c_n\frac{d{q}_n}{\mathrm{dt}}+\mathrm{\α}{K}_n{q}_n+\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nk}}{q}_n{q}_k^2=U_n---\left(2\right)$

式中：c_n＝2mξ_nω_n，ω_n表示第n阶自振圆频率；K_n为第n阶模态的刚度初始值，即是由于初始索力T₀引起的第n阶模态刚度，α是索力随时间变化的参数，其中T_u为由于支座轴向移动引起的索力变化， $T_u=\frac{\mathrm{EAu}}{L},{\mathrm{\β}}_{\mathrm{nk}}=\frac{\mathrm{EA}{k}^2{n}^2{\mathrm{\π}}^4}{8L^3};$

设索的振动只包含r₀阶模态，将前r₀阶的振动方程写成如下矩阵形式：

$M_{r_0}\stackrel{\·\·}{q}+C_{r_0}\stackrel{\·}{q}+\mathrm{\α}{K}_{r_0}q+{\mathrm{\Λ}}_{r_0}q=U_{r_0}---\left(1\right)$

式中：分别为r₀阶模态相应的质量矩阵、阻尼矩阵、刚度矩阵和非线性刚度矩阵；

$U_{r_0}={\left[\begin{array}{c}{U}_1\\ U_2\\ .\\ .\\ .\\ U_{r_0}\end{array}\right]}_{r_0\×N_t}---\left(6\right)$

式中：N_t为采样点数；

在进行拉索索力识别时，设采用的模态阶数为r，定义如下状态空间向量

$Z\left(t\right)={\left[\begin{array}{ccccccc}{q}_1& ...& q_r& {\stackrel{\·}{q}}_1& ...& {\stackrel{\·}{q}}_r& \mathrm{\α}\end{array}\right]}^T---\left(7\right)$

则公式(3)改写为如下状态空间方程：

$\stackrel{\·}{Z}\left(t\right)=g\left(Z\left(t\right)\right)=\left[\begin{array}{c}\stackrel{\·}{q}\\ M_r^{-1}(U_r-C_r\stackrel{\·}{q}-\mathrm{\α}\left(t\right)K_{r}q-{\mathrm{\Λ}}_{r}q)\\ 0_{1\×1}\end{array}\right]+w_0\left(t\right)---\left(10\right)$

式中：w₀(t)为系统噪声，式中：M_r、C_r、K_r、Λ_r分别为r阶模态相应的质量矩阵、阻尼矩阵、刚度矩阵和非线性刚度矩阵；

非线性项Λ_rq视为高斯白噪声过程，将其归入过程噪声项，公式(10)改写为：

$\stackrel{\·}{Z}\left(t\right)=A\left(t\right)Z\left(t\right)+B{U}_r\left(t\right)+w\left(t\right)---\left(8\right)$

式中： $A\left(t\right)=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& 0_{r\×1}\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right],B=\left[\begin{array}{c}{0}_{r\×r}\\ M_r^{-1}\\ 0_{1\×r}\end{array}\right];$ w(t)为系统的过程噪声，设w(t)为均值为0、方差为Q(t)的白噪声；

将公式(11)离散，得到

$\begin{array}{c}{Z}_k=f(Z_{k-1},U_{r,k})\\ =\mathrm{\Φ}(t_k,t_k-\mathrm{\Δt})Z_{k-1}+{\∫}_{t_k-\mathrm{\Δt}}^{t_k}\mathrm{\Φ}(t_k,t)B\left(t\right)U_r\left(t\right)\mathrm{dt}+w_{k-1}\end{array}---\left(9\right)$

式中： $\mathrm{\Φ}(t_2,t_1)=\exp \left({\∫}_{t_1}^{t_2}A\left(t\right)\mathrm{dt}\right);$

设拉索上安装了加速度传感器，则拉索的观测方程为：

式中：p为加速度传感器的数目，h(·)为观测函数，为第j个加速度传感器的位置；为振型函数，且v_k是第k个时间步的观测噪声，设它不随时间变化，是满足均值为0、方差为R_k的高斯白噪声；

采用公式(12)及公式(13)通过离散扩展卡尔曼(EKF)滤波方法识别拉索的时变索力，识别方法分为预测与更新两个过程：

(1)预测

首先，根据系统方程和前一步的状态变量最优估计值有

${\hat{Z}}_{k|k-1}=f({\hat{Z}}_{k-1|k-1},U_{r,k})---\left(11\right)$

式中：是k-1步的状态变量估计值，U_r，k是k步的系统输入，为根据k-1步的最优估计值和系统方程计算得到的第k步状态变量初步估计值；

通过系统方程的线性化得到初步估计状态变量的协方差估计值，结果如下：

${\hat{P}}_{k|k-1}=F_{k-1}{\hat{P}}_{k-1|k-1}{F}_{k-1}^T+Q_{k-1}---\left(12\right)$

式中：是k-1步的协方差矩阵估计值，U_r，k是k步的系统输入，P_k|k-1为根据k-1步的最优估计值和观测数据计算得到的第k步协方差矩阵估计值；

$G\left(t\right)=\frac{\∂g}{\∂Z}=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& -M_r^{-1}{K}_{r}q\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right];$

(2)更新

采用观测的拉索加速度响应，得到观测变量的初步预测误差

${\tilde{y}}_k=y_k-h\left({\hat{Z}}_{k|k-1}\right)---\left(13\right)$

式中：y_k是第k步的观测量；

得到卡尔曼增益如下：

${\mathrm{\κ}}_k=P_{k|k-1}{H}_k^T{S}_k^{-1}---\left(14\right)$

式中：

$S_k=H_k{P}_{k|k-1}{H}_k^T+R_k$

为观测噪声的协方差矩阵；

从而根据卡尔曼增益得到更新的第k步状态变量及其协方差的最优估计值

${\hat{Z}}_{k|k}={\hat{Z}}_{k|k-1}+{\mathrm{\κ}}_k{\tilde{y}}_k---\left(15\right)$

P_k|k＝(I-κ_kH_k)P_k|k-1      (16)

将上述过程在整个观测时间段内重复计算，则得到整个时间段内拉索的索力与时间的关系，即拉索索力时程。

对未安装风速仪的情况，将风荷载作为过程噪声，从而系统方程应改写为

$\stackrel{\·}{Z}\left(t\right)=g_2\left(Z\left(t\right)\right)=A\left(t\right)Z\left(t\right)+w_2\left(t\right)---\left(17\right)$

其离散形式如下：

Z_k＝f₂(Z_k-1)＝Φ(t_k，t_k-Δt)Z_k-1+w_k-1        (18)

式中： $\mathrm{\Φ}(t_2,t_1)=\exp \left({\∫}_{t_1}^{t_2}A\left(t\right)\mathrm{dt}\right)$

前节采用的加速度观测方程中包含有外荷载项，因此也需改为下式：

综上，无风速观测时的索力识别采用如下过程：

首先，根据系统方程和前一步的状态变量最优估计值有

${\hat{Z}}_{k|k-1}=\mathrm{\Φ}(t_k,t_k-\mathrm{\Δt}){\hat{Z}}_{k-1|k-1}---\left(20\right)$

式中：是来自上一步的状态变量估计值，为根据系统方程计算得到的第k步状态变量初步估计值；

通过系统方程的线性化得到协方差初步估计值

${\hat{P}}_{k|k-1}=F_{k-1}{\hat{P}}_{k-1|k-1}{F}_{k-1}^T+Q_{k-1}---\left(21\right)$

式中： $F_{k-1}=\frac{\∂f_2}{\∂Z}{|}_{{\hat{Z}}_{k-1|k-1}}=\frac{\∂f}{\∂Z}{|}_{{\hat{Z}}_{k-1|k-1}}=\exp \left({\∫}_{(k-1)\mathrm{\Δt}}^{\mathrm{k\Δt}}G\left(t\right)\mathrm{dt}\right)$

$G\left(t\right)=\frac{\∂g}{\∂Z}=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& -M_r^{-1}{K}_{r}q\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right]$

$Q_{k-1}=E\left[w_{k-1}{w}_{k-1}^T\right]$

采用观测到的加速度响应，得到观测值的估计误差

${\tilde{y}}_k=y_k-h_2({\hat{Z}}_{k|k-1},{\hat{Z}}_{k-1|k-1})---\left(22\right)$

继而得到卡尔曼增益如下：

${\mathrm{\κ}}_k={\hat{P}}_{k|k-1}{H}_k^T{S}_k^{-1}---\left(23\right)$

式中： $S_k=H_k{\hat{P}}_{k|k-1}{H}_k^T+R_k;R_k=E\left[v_k{v}_k^T\right]$ 为观测噪声的协方差矩阵；

从而根据卡尔曼增益得到更新的状态变量及其协方差的最优估计值

${\hat{Z}}_{k|k}={\hat{Z}}_{k|k-1}+{\mathrm{\κ}}_k{\tilde{y}}_k---\left(25\right)$

${\hat{P}}_{k|k}=(I-{\mathrm{\κ}}_k{H}_k){\hat{P}}_{k|k-1}---\left(26\right)$

。

本发明的有益效果和优点：

本发明利用监测加速度信息和拉索振动的状态空间系统方程，通过预测与更新两个过程实时识别拉索时变索力历程，同时可实现在已知或未知外界激励(风载荷和其他环境载荷)监测数据的情况下进行时变索力历程的辨识，方法直接有效、辨识精度高、尤其适用于拉索在线评估。

附图说明

图1南京某斜拉桥索编号及梁段划分图；

图2南京某斜拉桥模拟风速曲线；

图3南京某斜拉桥模拟单辆车作用索力时程曲线；

图4南京某斜拉桥模拟多辆车作用索力时程曲线；

图5南京某斜拉桥J03拉索单辆车作用下L/6通道加速度时程曲线；

图6南京某斜拉桥J03拉索多辆车作用下L/6通道加速度时程曲线；

图7南京某斜拉桥J03拉索单辆车作用下L/6通道30s加速度功率谱；

图8南京某斜拉桥J03拉索多辆车作用下L/6通道30s加速度功率谱；

图9南京某斜拉桥J03拉索单辆车作用下算法辨识与真实时变索力对比图；

图10南京某斜拉桥J03拉索单辆车作用下观测加速度与识别加速度对比图；

图11南京某斜拉桥J03拉索多辆车作用下算法辨识与真实时变索力对比图；

图12南京某斜拉桥J03拉索多辆车作用下观测加速度与识别加速度对比图；

图13拉索时变索力历程辨识试验装置；

图14试验拉索加速度传感器30s平面内加速时程图；

图15试验拉索加速度传感器30s平面外加速时程图；

图16试验拉索加速度传感器30s平面内加速度频谱图；

图17试验拉索加速度传感器30s平面外加速度频谱图；

图18基于平面内拉索加速度的算法辨识试验拉索索力与实测索力对比图；

图19基于平面外拉索加速度的算法辨识试验拉索索力与实测索力对比图；

具体实施方式

下面根据附图举例进一步说明：

实施例1

通过监测加速度信号和拉索振动方程构建状态空间方程，通过离散扩展卡尔曼滤波方法识别拉索时变索力，具体过程如下：

忽略索的垂度和轴向惯性影响，拉索的振动方程表示为：

$\frac{d^2{\tilde{q}}_n}{d{\mathrm{\τ}}^2}+2{\mathrm{\ξ}}_{n}n\frac{d{\tilde{q}}_n}{\mathrm{d\τ}}+n^2[1+\tilde{u}\left(\mathrm{\τ}\right)]{\tilde{q}}_n+\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}{\widetilde{\mathrm{\β}}}_{\mathrm{nk}}{\tilde{q}}_n{\tilde{q}}_k^2={\tilde{U}}_n\left(\mathrm{\τ}\right)---\left(1\right)$

式中：q_n为第n阶模态的广义位移；ξ_n为拉索第n阶模态阻尼比；u₀为拉索在初始索力T₀作用下的伸长量，E为拉索的弹性模量，A为拉索横截面积；L为拉索无应力长度；τ＝ω₁t，t为时间，τ为无量纲时间，ω₁为拉索的第一阶自振频率，μ为拉索单位长度质量；u为拉索支座沿轴向运动位移；为非线性振动系数；为无量纲外荷载，其中U_n为拉索上的外荷载，f_w(x，t)为作用在索上的横向分布时变荷载；m为拉索模态质量，即

上述模态坐标下的无量纲运动方程描述拉索沿轴向存在变形时的振动情况，公式(1)变为：

$m\frac{d^2{q}_n}{d{t}^2}+c_n\frac{d{q}_n}{\mathrm{dt}}+\mathrm{\α}{K}_n{q}_n+\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nk}}{q}_n{q}_k^2=U_n---\left(2\right)$

式中：c_n＝2mξ_nω_n，ω_n表示第n阶自振圆频率；K_n为第n阶模态的刚度初始值，即是由于初始索力T₀引起的第n阶模态刚度，α是索力随时间变化的参数，其中T_u为由于支座轴向移动引起的索力变化， $T_u=\frac{\mathrm{EAu}}{L},{\mathrm{\β}}_{\mathrm{nk}}=\frac{\mathrm{EA}{k}^2{n}^2{\mathrm{\π}}^4}{8L^3};$

设索的振动只包含r₀阶模态，将前r₀阶的振动方程写成如下矩阵形式：

$M_{r_0}\stackrel{\·\·}{q}+C_{r_0}\stackrel{\·}{q}+\mathrm{\α}{K}_{r_0}q+{\mathrm{\Λ}}_{r_0}q=U_{r_0}---\left(27\right)$

式中：分别为r₀阶模态相应的质量矩阵、阻尼矩阵、刚度矩阵和非线性刚度矩阵；

$U_{r_0}={\left[\begin{array}{c}{U}_1\\ U_2\\ .\\ .\\ .\\ U_{r_0}\end{array}\right]}_{r_0\×N_t}---\left(32\right)$

式中：N_t为采样点数；

在进行拉索索力识别时，设采用的模态阶数为r，定义如下状态空间向量

$Z\left(t\right)={\left[\begin{array}{ccccccc}{q}_1& ...& q_r& {\stackrel{\·}{q}}_1& ...& {\stackrel{\·}{q}}_r& \mathrm{\α}\end{array}\right]}^T---\left(33\right)$

则公式(3)改写为如下状态空间方程：

$\stackrel{\·}{Z}\left(t\right)=g\left(Z\left(t\right)\right)=\left[\begin{array}{c}\stackrel{\·}{q}\\ M_r^{-1}(U_r-C_r\stackrel{\·}{q}-\mathrm{\α}\left(t\right)K_{r}q-{\mathrm{\Λ}}_{r}q)\\ 0_{1\×1}\end{array}\right]+w_0\left(t\right)---\left(10\right)$

式中：w₀(t)为系统噪声，式中：M_r、C_r、K_r、Λ_r分别为r阶模态相应的质量矩阵、阻尼矩阵、刚度矩阵和非线性刚度矩阵；

非线性项Λ_rq视为高斯白噪声过程，将其归入过程噪声项，公式(10)改写为：

$\stackrel{\·}{Z}\left(t\right)=A\left(t\right)Z\left(t\right)+B{U}_r\left(t\right)+w\left(t\right)---\left(34\right)$

式中： $A\left(t\right)=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& 0_{r\×1}\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right],B=\left[\begin{array}{c}{0}_{r\×r}\\ M_r^{-1}\\ 0_{1\×r}\end{array}\right];$ w(t)为系统的过程噪声，设w(t)为均值为0、方差为Q(t)的白噪声；

将公式(11)离散，得到

$\begin{array}{c}{Z}_k=f(Z_{k-1},U_{r,k})\\ =\mathrm{\Φ}(t_k,t_k-\mathrm{\Δt})Z_{k-1}+{\∫}_{t_k-\mathrm{\Δt}}^{t_k}\mathrm{\Φ}(t_k,t)B\left(t\right)U_r\left(t\right)\mathrm{dt}+w_{k-1}\end{array}---\left(35\right)$

式中： $\mathrm{\Φ}(t_2,t_1)=\exp \left({\∫}_{t_1}^{t_2}A\left(t\right)\mathrm{dt}\right);$

设拉索上安装了加速度传感器，则拉索的观测方程为：

式中：p为加速度传感器的数目，h(·)为观测函数，为第j个加速度传感器的位置；为振型函数，且v_k是第k个时间步的观测噪声，设它不随时间变化，是满足均值为0、方差为R_k的高斯白噪声；

采用公式(12)及公式(13)通过离散扩展卡尔曼(EKF)滤波方法识别拉索的时变索力，识别方法分为预测与更新两个过程：

(1)预测

首先，根据系统方程和前一步的状态变量最优估计值有

${\hat{Z}}_{k|k-1}=f({\hat{Z}}_{k-1|k-1},U_{r,k})---\left(37\right)$

式中：是k-1步的状态变量估计值，U_r，k是k步的系统输入，为根据k-1步的最优估计值和系统方程计算得到的第k步状态变量初步估计值；

通过系统方程的线性化得到初步估计状态变量的协方差估计值，结果如下：

${\hat{P}}_{k|k-1}=F_{k-1}{\hat{P}}_{k-1|k-1}{F}_{k-1}^T+Q_{k-1}---\left(38\right)$

式中：是k-1步的协方差矩阵估计值，U_r，k是k步的系统输入，P_k|k-1为根据k-1步的最优估计值和观测数据计算得到的第k步协方差矩阵估计值；

$G\left(t\right)=\frac{\∂g}{\∂Z}=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& -M_r^{-1}{K}_{r}q\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right];$

(2)更新

采用观测的拉索加速度响应，得到观测变量的初步预测误差

${\tilde{y}}_k=y_k-h\left({\hat{Z}}_{k|k-1}\right)---\left(39\right)$

式中：y_k是第k步的观测量；

得到卡尔曼增益如下：

${\mathrm{\κ}}_k=P_{k|k-1}{H}_k^T{S}_k^{-1}---\left(40\right)$

式中：

$S_k=H_k{P}_{k|k-1}{H}_k^T+R_k$

为观测噪声的协方差矩阵；

从而根据卡尔曼增益得到更新的第k步状态变量及其协方差的最优估计值

${\hat{Z}}_{k|k}={\hat{Z}}_{k|k-1}+{\mathrm{\κ}}_k{\tilde{y}}_k---\left(41\right)$

P_k|k＝(I-κ_kH_k)P_k|k-1        (42)

将上述过程在整个观测时间段内重复计算，则得到整个时间段内拉索的索力与时间的关系，即拉索索力时程。

对未安装风速仪的情况，将风荷载作为过程噪声，从而系统方程应改写为

$\stackrel{\·}{Z}\left(t\right)=g_2\left(Z\left(t\right)\right)=A\left(t\right)Z\left(t\right)+w_2\left(t\right)---\left(43\right)$

其离散形式如下：

Z_k＝f₂(Z_k-1)＝Φ(t_k，t_k-Δt)Z_k-1+w_k-1        (44)

式中： $\mathrm{\Φ}(t_2,t_1)=\exp \left({\∫}_{t_1}^{t_2}A\left(t\right)\mathrm{dt}\right)$

前节采用的加速度观测方程中包含有外荷载项，因此也需改为下式：

综上，无风速观测时的索力识别采用如下过程：

首先，根据系统方程和前一步的状态变量最优估计值有

${\hat{Z}}_{k|k-1}=\mathrm{\Φ}(t_k,t_k-\mathrm{\Δt}){\hat{Z}}_{k-1|k-1}---\left(46\right)$

式中：是来自上一步的状态变量估计值，为根据系统方程计算得到的第k步状态变量初步估计值；

通过系统方程的线性化得到协方差初步估计值

${\hat{P}}_{k|k-1}=F_{k-1}{\hat{P}}_{k-1|k-1}{F}_{k-1}^T+Q_{k-1}---\left(47\right)$

式中： $F_{k-1}=\frac{\∂f_2}{\∂Z}{|}_{{\hat{Z}}_{k-1|k-1}}=\frac{\∂f}{\∂Z}{|}_{{\hat{Z}}_{k-1|k-1}}=\exp \left({\∫}_{(k-1)\mathrm{\Δt}}^{\mathrm{k\Δt}}G\left(t\right)\mathrm{dt}\right)$

$G\left(t\right)=\frac{\∂g}{\∂Z}=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& -M_r^{-1}{K}_{r}q\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right]$

$Q_{k-1}=E\left[w_{k-1}{w}_{k-1}^T\right]$

采用观测到的加速度响应，得到观测值的估计误差

${\tilde{y}}_k=y_k-h_2({\hat{Z}}_{k|k-1},{\hat{Z}}_{k-1|k-1})---\left(48\right)$

继而得到卡尔曼增益如下：

${\mathrm{\κ}}_k={\hat{P}}_{k|k-1}{H}_k^T{S}_k^{-1}---\left(49\right)$

式中： $S_k=H_k{\hat{P}}_{k|k-1}{H}_k^T+R_k;R_k=E\left[v_k{v}_k^T\right]$ 为观测噪声的协方差矩阵；

从而根据卡尔曼增益得到更新的状态变量及其协方差的最优估计值

${\hat{Z}}_{k|k}={\hat{Z}}_{k|k-1}+{\mathrm{\κ}}_k{\tilde{y}}_k---\left(51\right)$

${\hat{P}}_{k|k}=(I-{\mathrm{\κ}}_k{H}_k){\hat{P}}_{k|k-1}---\left(52\right)$

。

步骤1：传感器安装及控制模态辨识。在拉索的平面内或平面外布设1个或多个加速度传感器，针对拉索加速度传感器一个通道的一段加速度时程曲线(取30秒)，计算其功率谱，辨识拉索基频以及其他频率成分，拉索的基频用f₁表示利用功率谱的幅值辨识拉索振动的控制模态；

步骤2：建立拉索状态空间系统方程和观测方程。根据拉索的物理特性(拉索长度、单位长度质量、恒荷载索力等)、拉索振动控制模态和拉索加速度传感器安装信息，建立拉索状态空间系统方程和观测方法；

步骤3：在已知监测风荷载输入和加速度响应输出的情况下，根据k-1步的最优估计值和系统方程计算第k步状态变量初步估计值，初步估计状态变量的协方差估计值通过系统方程的线性化得到，采用观测的拉索加速度响应，计算观测变量的初步预测误差，根据卡尔曼增益更新第k步状态变量及其协方差的最优估计值。

实施例2

如图1所示，南京某斜拉桥上的J03拉索，由109根7mm钢丝组成，索长L＝112.029m，截面面积A＝4.195×10^-3m2，单位长度的密度为μ＝32.93kg/m。

步骤一：本发明算法需要利用拉索监测加速度传感器信息辨识拉索索力，通过数值模拟计算J03拉索在图2所示风速和图3、图4所示索力作用下L/6位置的加速度响应。图5、图6分别为J03拉索在单辆车作用下和多辆车作用下L/6通道的30秒加速度时程，对其进行功率谱分析结果如图7、图8所示，从图中可以看出主要是前11阶模态参与了振动，其他模态的作用非常小，拉索的振动控制模态取为前20阶是可行的。

步骤二：利用步骤一选择的20阶控制模态和拉索的物理参数，建立拉索状态空间系统方程和观测方程，其中过程噪声和观测噪声的协方差矩阵分别取为Q＝diag([0_2r×1；1e-6])，R＝1e-3，式中r表示识别中所采用的模态阶数。

步骤三：采用扩展卡尔曼滤波器的预测过程和更新过程识别拉索时变索力。单辆车和多辆车作用下索力时程曲线和加速度时程曲线的对比结果如图9-12所示。模拟索力时程曲线与辨识索力时程曲线很好的吻合在一起，验证了本发明所提算法的准确性。

实施例3

试验拉索长14.02m，直径为1.5cm，单位长度质量为1.33kg/m。，试验拉索一端固定，另一端采用螺纹杆调节索的拉伸长度，进而调节索力的大小，同时在试验拉索端部布置测力计，监测索力时程的变化，试验拉索采用两个风机作为外界激励产生振动，模拟拉索外界激励信息未知条件下拉索索力的识别。

步骤一：分别在2.43m、3.60m处布设加速度传感器测试平面内及平面外振动响应，采样频率为200Hz，试验装置、传感器布置以及几何尺寸如图13所示。对如图14、图15所示的试验拉索加速度传感器的30s平面内和平面外加速度时程进行功率谱分析，结果如图16、图17所示，试验拉索激励起来的最高阶频率为第8阶；

步骤二：利用步骤一选择的20阶控制模态和拉索的物理参数，建立拉索状态空间系统方程和观测方程。时变索力识别程序中，过程噪声的均值为0，过程噪声的协方差矩阵由加速度响应的频谱幅值确定，取为Q＝diag([298 849 436 42 599 894 201 31 396 542 648])×10^-5；观测噪声的均值和方差分别为0和0.1；

步骤三：采用扩展卡尔曼滤波器的预测过程和更新过程识别拉索时变索力。基于平面内和平面外加速度的索力时程曲线对比如图18、图19所示。

外界激励位置的情况下，试验拉索的监测索力时程曲线与本发明所提算法辨识索力时程曲线很好的吻合在一起，从试验角度验证了本发明所提算法的准确性。

Claims (3)

1.一种基于扩展卡尔曼滤波器的拉索时变索力历程识别方法，其特征在于，包括以下步骤：

步骤1：传感器安装及控制模态辨识；在拉索的平面内或平面外布设一个或多个加速度传感器，针对加速度传感器一个通道的一段加速度时程曲线，计算其功率谱，辨识拉索基频以及其他频率成分，利用功率谱的幅值辨识拉索振动的控制模态；

步骤2：建立拉索状态空间系统方程和观测方程；根据拉索的物理特性：拉索长度、单位长度质量、恒荷载索力，拉索振动控制模态和加速度传感器安装信息，建立拉索状态空间系统方程和观测方法；

步骤3：在已知监测风荷载输入和加速度响应输出的情况下，根据k-1步的最优估计值和系统方程计算第k步状态变量初步估计值，初步估计状态变量的协方差估计值通过系统方程的线性化得到，采用观测的拉索加速度响应，计算观测变量的初步预测误差，根据卡尔曼增益更新第k步状态变量及其协方差的最优估计值；

若拉索外界激励信息未知，将风荷载等外界激励视为系统噪声的一部分，改写拉索振动状态空间系统方程，重复步骤1-3，同样采用扩展卡尔曼滤波器的预测过程和更新过程识别拉索时变索力。

2.根据权利要求1所述的一种基于扩展卡尔曼滤波器的拉索时变索力历程识别方法，其特征在于，通过监测加速度信号和拉索振动方程构建状态空间方程，通过离散扩展卡尔曼滤波方法识别拉索时变索力，具体过程如下：

忽略索的垂度和轴向惯性影响，拉索的振动方程表示为：

$\frac{d^2{\tilde{q}}_n}{d{\mathrm{\τ}}^2}+2{\mathrm{\ξ}}_n\stackrel{\·}{n}\frac{d{\tilde{q}}_n}{\mathrm{d\τ}}+n^2[1+\tilde{u}\left(\mathrm{\τ}\right)]{\tilde{q}}_n+\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}{\widetilde{\mathrm{\β}}}_{\mathrm{nk}}{\tilde{q}}_n{\tilde{q}}_k^2={\tilde{U}}_n\left(\mathrm{\τ}\right)---\left(1\right)$

式中：q_n为第n阶模态的广义位移；ξ_n为拉索第n阶模态阻尼比；u₀为拉索在初始索力T₀作用下的伸长量，E为拉索的弹性模量，A为拉索横截面积；L为拉索无应力长度；τ＝ω₁t，t为时间，τ为无量纲时间，ω₁为拉索的第一阶自振频率，μ为拉索单位长度质量；u为拉索支座沿轴向运动位移；为非线性振动系数；为无量纲外荷载，其中U_n为拉索上的外荷载，f_w(x，t)为作用在索上的横向分布时变荷载；m为拉索模态质量，即

上述模态坐标下的无量纲运动方程描述拉索沿轴向存在变形时的振动情况，公式(1)变为：

$m\frac{d^2{q}_n}{d{t}^2}+c_n\frac{d{q}_n}{\mathrm{dt}}+\mathrm{\α}{K}_n{q}_n+\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nk}}{q}_n{q}_k^2=U_n---\left(2\right)$

式中：c_n＝2mξ_nω_n，ω_n表示第n阶自振圆频率；K_n为第n阶模态的刚度初始值，即是由于初始索力T₀引起的第n阶模态刚度，α是索力随时间变化的参数，其中T_u为由于支座轴向移动引起的索力变化， ${\mathrm{\β}}_{\mathrm{nk}}=\frac{\mathrm{EA}{k}^2{n}^2{\mathrm{\π}}^4}{8L^3};$

设索的振动只包含r₀阶模态，将前r₀阶的振动方程写成如下矩阵形式：

$M_{r_0}\stackrel{\·\·}{q}+C_{r_0}\stackrel{\·}{q}+\mathrm{\α}{K}_{r_0}q+{\mathrm{\Λ}}_{r_0}q=U_{r_0}---\left(1\right)$

式中：分别为r₀阶模态相应的质量矩阵、阻尼矩阵、刚度矩阵和非线性刚度矩阵；

$U_{r_0}={\left[\begin{array}{c}{U}_1\\ U_2\\ .\\ .\\ .\\ U_{r_0}\end{array}\right]}_{r_0\×N_t}---\left(6\right)$

式中：N_t为采样点数；

在进行拉索索力识别时，设采用的模态阶数为r，定义如下状态空间向量

$Z\left(t\right)={\left[\begin{array}{ccccccc}{q}_1& ...& q_r& {\stackrel{\·}{q}}_1& ...& {\stackrel{\·}{q}}_r& \mathrm{\α}\end{array}\right]}^T---\left(7\right)$

则公式(3)改写为如下状态空间方程：

$\stackrel{\·}{Z}\left(t\right)=g\left(Z\left(t\right)\right)=\left[\begin{array}{c}\stackrel{\·}{q}\\ M_r^{-1}(U_r-C_r\stackrel{\·}{q}-\mathrm{\α}\left(t\right)K_{r}q-{\mathrm{\Λ}}_{r}q)\\ 0_{1\×1}\end{array}\right]+w_0\left(t\right)---\left(10\right)$

式中：w₀(t)为系统噪声，式中：M_r、C_r、K_r、Λ_r分别为r阶模态相应的质量矩阵、阻尼矩阵、刚度矩阵和非线性刚度矩阵；

非线性项Λ_rq视为高斯白噪声过程，将其归入过程噪声项，公式(10)改写为：

$\stackrel{\·}{Z}\left(t\right)=A\left(t\right)Z\left(t\right)+B{U}_r\left(t\right)+w\left(t\right)---\left(8\right)$

式中： $A\left(t\right)=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& 0_{r\×1}\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right],B=\left[\begin{array}{c}{0}_{r\×r}\\ M_r^{-1}\\ 0_{1\×r}\end{array}\right];$ w(t)为系统的过程噪声，设w(t)为均值为0、方差为Q(t)的白噪声；

将公式(11)离散，得到

$\begin{array}{c}{Z}_k=f(Z_{k-1},U_{r,k})\\ =\mathrm{\Φ}(t_k,t_k-\mathrm{\Δt})Z_{k-1}+{\∫}_{t_k-\mathrm{\Δt}}^{t_k}\mathrm{\Φ}(t_k,t)B\left(t\right)U_r\left(t\right)\mathrm{dt}+w_{k-1}\end{array}---\left(9\right)$

式中： $\mathrm{\Φ}(t_2,t_1)=\exp \left({\∫}_{t_1}^{t_2}A\left(t\right)\mathrm{dt}\right);$

设拉索上安装了加速度传感器，则拉索的观测方程为：

式中：p为加速度传感器的数目，h(·)为观测函数，为第j个加速度传感器的位置；为振型函数，且v_k是第k个时间步的观测噪声，设它不随时间变化，是满足均值为0、方差为R_k的高斯白噪声；

采用公式(12)及公式(13)通过离散扩展卡尔曼(EKF)滤波方法识别拉索的时变索力，识别方法分为预测与更新两个过程：

(1)预测

首先，根据系统方程和前一步的状态变量最优估计值有

${\hat{Z}}_{k|k-1}=f({\hat{Z}}_{k-1|k-1},U_{r,k})---\left(11\right)$

式中：是k-1步的状态变量估计值，U_r，k是k步的系统输入，为根据k-1步的最优估计值和系统方程计算得到的第k步状态变量初步估计值；

通过系统方程的线性化得到初步估计状态变量的协方差估计值，结果如下：

${\hat{P}}_{k|k-1}=F_{k-1}{\hat{P}}_{k-1|k-1}{F}_{k-1}^T+Q_{k-1}---\left(12\right)$

式中：是k-1步的协方差矩阵估计值，U_r，k是k步的系统输入，P_k|k-1为根据k-1步的最优估计值和观测数据计算得到的第k步协方差矩阵估计值；

$G\left(t\right)=\frac{\∂g}{\∂Z}=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& -M_r^{-1}{K}_{r}q\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right];$

(2)更新

采用观测的拉索加速度响应，得到观测变量的初步预测误差

${\tilde{y}}_k=y_k-h\left({\hat{Z}}_{k|k-1}\right)---\left(13\right)$

式中：y_k是第k步的观测量；

得到卡尔曼增益如下：

${\mathrm{\κ}}_k=P_{k|k-1}{H}_k^T{S}_k^{-1}---\left(14\right)$

式中：

$S_k=H_k{P}_{k|k-1}{H}_k^T+R_k$

为观测噪声的协方差矩阵；

从而根据卡尔曼增益得到更新的第k步状态变量及其协方差的最优估计值

${\hat{Z}}_{k|k}={\hat{Z}}_{k|k-1}+{\mathrm{\κ}}_k{\tilde{y}}_k---\left(15\right)$

P_k|k＝(I-κ_kH_k)P_k|k-1          (16)

将上述过程在整个观测时间段内重复计算，则得到整个时间段内拉索的索力与时间的关系，即拉索索力时程。

3.根据权利要求1或2所述的一种基于扩展卡尔曼滤波器的拉索时变索力历程识别方法，其特征在于，

对未安装风速仪的情况，将风荷载作为过程噪声，从而系统方程应改写为

$\stackrel{\·}{Z}\left(t\right)=g_2\left(Z\left(t\right)\right)=A\left(t\right)Z\left(t\right)+w_2\left(t\right)---\left(17\right)$

其离散形式如下：

Z_k＝f₂(Z_k-1)＝Φ(t_k，t_k-Δt)Z_k-1+w_k-1         (18)

式中： $\mathrm{\Φ}(t_2,t_1)=\exp \left({\∫}_{t_1}^{t_2}A\left(t\right)\mathrm{dt}\right)$

前节采用的加速度观测方程中包含有外荷载项，因此也需改为下式：

综上，无风速观测时的索力识别采用如下过程：

首先，根据系统方程和前一步的状态变量最优估计值有

${\hat{Z}}_{k|k-1}=\mathrm{\Φ}(t_k,t_k-\mathrm{\Δt}){\hat{Z}}_{k-1|k-1}---\left(20\right)$

式中：是来自上一步的状态变量估计值，为根据系统方程计算得到的第k步状态变量初步估计值；

通过系统方程的线性化得到协方差初步估计值

${\hat{P}}_{k|k-1}=F_{k-1}{\hat{P}}_{k-1|k-1}{F}_{k-1}^T+Q_{k-1}---\left(21\right)$

式中： $F_{k-1}=\frac{\∂f_2}{\∂Z}{|}_{{\hat{Z}}_{k-1|k-1}}=\frac{\∂f}{\∂Z}{|}_{{\hat{Z}}_{k-1|k-1}}=\exp \left({\∫}_{(k-1)\mathrm{\Δt}}^{\mathrm{k\Δt}}G\left(t\right)\mathrm{dt}\right)$

$G\left(t\right)=\frac{\∂g}{\∂Z}=\left[\begin{array}{ccc}{0}_{r\×r}& I_{r\×r}& 0_{r\×1}\\ -M_r^{-1}\left(\mathrm{\α}\left(t\right)K_r\right)& -M_r^{-1}{C}_r& -M_r^{-1}{K}_{r}q\\ 0_{1\×r}& 0_{1\×r}& 0_{1\×1}\end{array}\right]$

$Q_{k-1}=E\left[w_{k-1}{w}_{k-1}^T\right]$

采用观测到的加速度响应，得到观测值的估计误差

${\tilde{y}}_k=y_k-h_2({\hat{Z}}_{k|k-1},{\hat{Z}}_{k-1|k-1})---\left(22\right)$

继而得到卡尔曼增益如下：

${\mathrm{\κ}}_k={\hat{P}}_{k|k-1}{H}_k^T{S}_k^{-1}---\left(23\right)$

式中： $S_k=H_k{\hat{P}}_{k|k-1}{H}_k^T+R_k;R_k=E\left[v_k{v}_k^T\right]$ 为观测噪声的协方差矩阵；

从而根据卡尔曼增益得到更新的状态变量及其协方差的最优估计值

${\hat{Z}}_{k|k}={\hat{Z}}_{k|k-1}+{\mathrm{\κ}}_k{\tilde{y}}_k---\left(25\right)$

${\hat{P}}_{k|k}=(I-{\mathrm{\κ}}_k{H}_k){\hat{P}}_{k|k-1}---\left(26\right)$