基于监测加速度的拉索时变索力历程识别的数据驱动方法

摘要

本发明提出了一种基于监测加速度的拉索时变索力历程识别的数据驱动方法，运用复杂度寻踪这种高效的无监督学习算法，利用拉索上布置的多通道加速度传感器的监测信息，实现了对索力时程进行实时辨识。复杂度寻踪算法可以自主地将拉索的加速度响应分解为拉索的单模态响应，进而通过极短时间的加速度信息辨识拉索的实时频率，通过张紧弦理论计算索力时程。通过结合实测索力和实测风速的斜拉桥的模拟分析以及拉索的模型试验，证明了所提出的复杂度寻踪算法能够对时变索力历程进行准确的实时识别。本发明是一种直接有效的时变索力历程辨识方法，简单易用，索力辨识精度高，时效性强且能够实现在线实时辨识，尤其适用于拉索的在线评估。

说明书

技术领域

本发明涉及一种土木工程结构响应辨识的方法，具体涉及一种基 于监测加速度的拉索时变索力历程识别的数据驱动方法。

背景技术

斜拉桥由于其跨度大、造型美观、施工方便等特点，是目前世界 上应用最广泛的桥型。斜拉桥一般会成为一个地区的标志性建筑，而 且往往是一个地区的交通枢纽，对该地区的政治经济有着至关重要的 影响。

斜拉索作为斜拉桥主要的承重构件，由高强度钢丝束和PE护套 组成。斜拉索在长达几十年使用期内，环境侵蚀、材料老化和荷载的 长期效应、疲劳效应与突变效应等灾害因素的耦合作用下，将不可避 免地导致结构和系统的损伤积累和抗力衰减，从而抵抗自然灾害、甚 至正常环境作用的能力下降，极端情况下将引发灾难性的突发事故。 为保证桥梁结构的正常运营，最好的解决方法就是对斜拉索进行定期 的检测维护，且需要依据检测结果对其进行安全评估。但是由于无损 检测方法(漏磁检测、X射线检测、超声检测、基于振动法的索力检 测等)的局限性，在斜拉桥实际运营过程中，只有极小部分斜拉索都 得到检测。鉴于斜拉索病害的普遍性，以及斜拉索修复和更换的高昂 费用，为应对斜拉索结构运营状况下进行腐蚀疲劳评估的需要，目前 亟待认识斜拉索结构在服役全过程中的结构索力历程响应。

目前存在的索力测试装置有：压力环、磁通量传感器、光纤光栅 智能拉索等，其中压力环和光纤光栅智能拉索能够直接测试索力历 程，广泛的安装和应用于新建桥梁上。但是，这些索力监测传感器一 般价格昂贵，安装复杂(只能用于新建桥梁)，更主要的是传感器的耐 久性比较差，这些固有缺点限制了上述传感器的大规模应用。由于更 换索力监测传感器耗时耗力，价格昂贵，因而亟需发展一种省时、省 力、经济的实时索力历程监测方法。

发明内容

基于以上不足之处，本发明提供一种基于多通道监测加速度辨识 时变拉索索力历程的数据驱动方法，解决拉索时变索力历程识别的问 题。

本发明采用如下技术方案实现：一种基于监测加速度的拉索时变 索力历程识别的数据驱动方法，步骤如下：

步骤1：对待测试拉索同一平面内布设2-3个的加速度传感器， 测试斜拉索在环境激励下的加速度响应，并通过一个通道的10秒内 的加速度信号，辨识拉索基频及该时间段内能够激励起来的最高频 率，拉索的基频用f₁表示，监测加速度信号中能够辨识的拉 索的最高频率用f_i表示，f_i≈i×f₁，其中，表示拉 索在索力T作用下，第i阶模态的圆频率，μ为单位索长的密度，L 为斜拉索长度为L；

步骤2：设计高通滤波器，截止频率(i-1.5)f₁，将拉索的多通道 加速度信号中出现的最高两阶频率f_i和f_i-1提取出来；

步骤3：选择3秒的窗函数，对窗函数内的多通道加速度信号进 行滤波预处理，预处理后加速度信号通过复杂度寻踪算法分离得到单 模态响应信号，通过快速傅里叶变换或功率谱计算单模态响应信号的 频率，利用张紧弦理论，通过辨识的频率计算拉索的时变索力；

步骤4：滑动窗函数，对窗函数内的加速度信号重复步骤3，辨 识拉索的时变索力历程。

本发明还具有如下技术特征：

1、应通过多通道加速度信号辨识拉索的时变频率，进而辨识拉索 的时变索力历程。

2、对多通道加速度信号进行高通滤波预处理，提取拉索能够被激 励起的最高的2阶频率。

3、利用复杂度寻踪处理算法，将滑动窗内多模态响应分解为单模 态响应，频率分辨率Δf≤0.025HZ的频率分辨率准确辨识单模态响应的 频率，利用张紧弦理论计算时变索力历程。

本发明基于成熟的加速度传感器测试技术，相对于目前发展的其 他类型的索力监测传感器，本发明采用的加速度传感器成熟可靠，测 试精度高，价格便宜，传感器安装以及更换都很方便，从而使得本发 明提出的拉索时变索力历程辨识系统具有极高的可靠性和耐久性。本 发明监测加速度信息，通过辨识时变频率计算索力，方法简单易用， 索力辨识精度高，时效性强且能够实现在线实时辨识，本发明方法的 鲁棒性及可靠性强，尤其适用于拉索的在线评估。

附图说明

图1为Benchmark桥梁的智能拉索制作构造和C8拉索位置简图；

图2为Benchmark桥梁实测风速和索力时程曲线图(测试时间为 2008年1月17日1:00-2:00AM，红框内为选择的待辨识索力)；

图3为Benchmark桥梁拉索加速度时程曲线图(开始时间为2008 年1月17日1:40AM)；

图4为C8拉索八分点通道10秒加速度功率谱图；

图5为C8拉索倍频特性图；

图6为C8拉索八分点通道和四分点通道3秒加速度滤波后的时 程曲线及其功率谱图；

图7为复杂度寻踪得到的单模态响应时程曲线及其功率谱图；

图8为辨识的C8拉索时变频率图；

图9为本发明算法辨识C8拉索索力与实测索力对比图；

图10为拉索时变索力历程辨识试验装置简图；

图11为S#1测试加速度10秒功率谱图；

图12为试验拉索倍频关系图

图13为试验拉索S#1和S#3传感器3秒加速度滤波后的时程曲 线及其功率谱；

图14为复杂度寻踪得到的单模态响应时程曲线及其功率谱图；

图15为辨识的试验拉索时变频率图；

图16为本发明算法辨识试验拉索索力与实测索力对比工况1图；

图17为本发明算法辨识试验拉索索力与实测索力对比工况2图。

具体实施方式

本发明的具体实施方案，通过结合实测索力和实测风速的斜拉桥 的模拟分析以及拉索的模型试验进行说明。

对于小垂度的斜拉索(长度为L，截面面积为A，杨氏模量为E)， 其平面外的第i阶模态的横向振动运动方程可以表示为

$m_i{\stackrel{··}{q}}_i\left(t\right)+2{\zeta }_i{m}_i{\omega }_i{\stackrel{·}{q}}_i\left(t\right)+m_i{\omega }_i^2{q}_i\left(t\right)+{\alpha }_i{q}_i\left(t\right)+\underset{k=1}{\overset{n}{\Sigma }}{\beta }_{\mathrm{ik}}{q}_i\left(t\right)q_k^2\left(t\right)=F_i\left(t\right)---\left(1\right)$

式中，m_i＝μL／2表示单位索长的重量，μ表示单位索长的密度；ζ_i为 第i阶模态阻尼比；α_i=i²π²T／2L，β_ik=EAπ⁴i²k²／8L³；F_i(t)为外界激 励；表示拉索在索力T作用下，第i阶模态的圆频率，ω_i的 单位是rad／s。的计算公式表明，拉索各阶圆频率与基频之 间存在倍频关系，即ω_i≈i×ω₁；同时拉索的各阶模态圆频率ω_i都和拉 索索力T之间存在直接的映射关系。依据张紧弦理论，拉索索力T可 以表示为

$T=\frac{{\mathrm{\mu L}}^2}{{\pi }^2}{\left(\frac{{\omega }_i}{i}\right)}^2\approx \frac{{\mathrm{\mu L}}^2{\omega }_1^2}{{\pi }^2}---\left(2\right)$

本发明所述的基于监测加速度的拉索时变索力历程识别的数据 驱动方法的核心即是时变模态圆频率ω_i的辨识算法。依据模态叠加理 论，环境激励下拉索的位移向量x(t)=[x₁(t)，…，x_n(t)]^T可以表示为

$x\left(t\right)=\mathrm{\Phi q}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\phi }_i{q}_i\left(t\right)\iff \tilde{q}\left(t\right)={\tilde{\Phi }}^{-1}x\left(t\right)---\left(3\right)$

式中，表示模态振型矩阵，q(t)=[q₁(t)，…，q_n(t)]^T为模态响应向量， 第i阶振型(Φ矩阵的第i列)和模态响应分别表示为φ_i和q_i(t)。上述 位移向量的分解计算通过盲元分离算法实现，令表示模态 响应向量估计值，表示模态振型向量估计值，y_i(t)信号的复杂 度定义为

$F\left(y_i\right)=\log \frac{V(w_i,x)}{U(w_i,x)}=\log \frac{w_i\bar{R}{w}_i^T}{w_i\hat{R}{w}_i^T}---\left(4\right)$

式中，和表示n×n的长期预测指标和短期预测指标的方差矩阵， 和的矩阵元素分别表示为

$\left(\begin{array}{c}\overline{r_{\mathrm{ij}}}=\underset{t=1}{\overset{N}{\Sigma }}[x_i\left(t\right)-\overline{x_i}\left(t\right)][x_j\left(t\right)-\overline{x_j}\left(t\right)]\\ \widehat{r_{\mathrm{ij}}}=\underset{t=1}{\overset{N}{\Sigma }}[x_i\left(t\right)-\widehat{x_i}\left(t\right)][x_j\left(t\right)-\widehat{x_j}\left(t\right)]\end{array}\right)---\left(5\right)$

长期预测值和短期预测值计算公式为

$\left(\begin{array}{cc}\overline{x_i}\left(t\right)={\lambda }_L\overline{x_i}(t-1)+(1-{\lambda }_L)x_i(t-1)& 0\le {\lambda }_L\le 1\\ \widehat{x_i}\left(t\right)={\lambda }_S\widehat{x_i}(t-1)+(1-{\lambda }_S)x_i(t-1)& 0\le {\lambda }_S\le 1\end{array}\right)---\left(6\right)$

本发明中，${\lambda }_S=2^{-1/h_s}=1/2,{\lambda }_S=2^{-1/h_L}=2^{-1/900000,}$h_L＞＞h_S，和矩阵 可以通过快速卷积算法计算获得。如果信号y_i(t)为理想的单模态响应 信号，信号复杂度函数值F(y_i)最大，因此只需要对信号复杂度函数 F(y_i)求极值，最优响应信号y_i(t)就能准确描述模态响应向量的重构 值。

复杂度函数F(y_i)求极值的过程可以通过经典的梯度优化算法实 现，复杂度函数F(y_i)对w_i的微分可以表示为

${▿}_{w_i}F=\frac{{2w}_i}{V_i}\bar{R}-\frac{{2w}_i}{U_i}\hat{R}---\left(7\right)$

当F=0时，复杂度函数F(y_i)即可达到最优点，上式可以改写 为

${▿}_{w_i}F=\frac{{2w}_i}{V_i}\bar{R}-\frac{{2w}_i}{U_i}\hat{R}=0\Rightarrow w_i\bar{R}=\frac{V_i}{U_i}{w}_i\hat{R}---\left(8\right)$

上式变形为一个广义特征值问题，所以w_i即为的特征向量，特征 值为γ_i=V_i／U_i。获得w_i后即可有公式(3)得到重构的单模态响应向量 上述通过优化信号复杂度函数F(y_i)，通过监测位移响应 向量x(t)获得为单模态响应q(t)=[q₁(t)，…，q_n(t)]^T的过程称为复杂度寻踪 算法。对于加速度响应和速度响应可以通过同样的方式进行计算。

对于单模态响应信号q(t)=[q₁(t)，…，q_n(t)]^T，由于每个模态响应中只 有控制模态的频率成分占优，所以只需要对分离出的 q(t)=[q₁(t)，…，q_n(t)]^T进行快速傅里叶(FFT)变换或功率谱分析，即可求得 信号的模态频率，通过公式(2)所示拉索索力与拉索频率之间的关系， 求得拉索索力。

对于单模态响应信号q_i(t)而言，频率分辨率Δf为选择的窗函数时 间τ的倒数Δf=1／τ，存在误差项的索力计算公式为

$\left(\begin{array}{c}\tilde{T}=\frac{{\mathrm{\mu L}}^2}{{\pi }^2}{\left(\frac{{\tilde{\omega }}_i}{i}\right)}^2=\frac{{\mathrm{\mu L}}^2}{{\pi }^2}{\left(\frac{{\tilde{\omega }}_i\pm 2\mathrm{\pi \Delta f}}{i}\right)}^2\approx \frac{{\mathrm{\mu L}}^2}{{\pi }^2}[{\left(\frac{{\tilde{\omega }}_i}{i}\right)}^2\pm \frac{4\mathrm{\pi \Delta f}{\tilde{\omega }}_i}{i^2}+{\left(\frac{2\mathrm{\pi \Delta f}}{i}\right)}^2]\\ \approx T\pm \frac{{4\mathrm{\mu L}}^2}{\pi }\left(\frac{{\tilde{\omega }}_i}{i}\right)\frac{\mathrm{\Delta f}}{i}\end{array}\right)---\left(9\right)$

式中，约等于拉索的基频，误差项的大小取决于Δf／i。从公式 (9)可以看出，对于时间长度比较短的窗函数，只要拉索的辨识频率 阶数(一般取i>10)比较高，也能达到很高的频率分辨率。

实施例1：

Benchmark斜拉桥(图1上的C8斜拉索)，由139根5mm钢丝组 成，索长L=100.95m，截面面积A=2.73×10^-3m²，单位长度的密 度为μ=21.43kg／m。2006年，C8斜拉索在换索工程中，更换为智 能拉索，能够实时监测索力历程。

步骤1：本发明算法需要利用2个或最多3个通道的加速度传感 器信息辨识拉索索力，通过数值模拟计算C8斜拉索在监测风速(图2 (a))和监测索力(图2(b))作用下L／8和L／4位置的加速度响应。 对C8拉索L／8通道的10秒加速度时程(图3)进行功率谱分析(图4)， C8拉索激励起来的最高阶频率为第30阶，同时可以看出拉索频率之 间存在明显的倍频关系(图5)。

步骤2：设计高通滤波器，截止频率约等于(30-1.5)f₁≈38Hz，将第 29阶和30阶频率成分过滤出来；

步骤3：选择3秒的窗函数，对窗函数内的加速度曲线应用步骤2 设计的滤波器进行滤波，滤波后的L／8和L／4通道的加速度时程曲线 及其功率谱见图6。对滤波后的L／8和L／4通道的加速度时程曲线进 行复杂度寻踪算法处理，分离出单模态响应信号，并计算其功率谱(图 7)，辨识单模态信号的占主导模态的频率，并利用张紧弦理论计算 索力。

步骤4：逐点滑动3秒的时间窗，重复步骤3，对滑动窗内加速度 信号滤波，分离单模态响应信号，通过傅里叶变换或者功率谱分析辨 识频率，从而得到频率随时间变化的曲线(图8)，进而得到索力时程 曲线(图9)。

监测索力时程曲线与辨识索力时程曲线很好的吻合在一起，验证 了本发明所提算法的准确性。

实施例2：

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

步骤1：分别在L／4、L／2和3／4L布设加速度传感器S#1-S#3测 试平面外响应，采样频率为200Hz，试验装置、传感器布置以几何尺 寸如图8所示。对试验拉索S#1传感器的10秒加速度时程进行功率 谱分析(图11)，试验拉索激励起来的最高阶频率为第14阶，同时 可以看出试验拉索频率之间存在明显的倍频关系(图12)；

步骤2：设计高通滤波器，截止频率约等于(14-1.5)f₁≈28Hz，能够 将第13阶和14阶频率成分过滤出来；

步骤3：选择3秒的窗函数，对窗函数内的加速度曲线应用步骤2 设计的滤波器进行滤波，滤波后的S#1和S#3传感器的加速度时程曲 线及其功率谱见图13。对滤波后的S#1和S#3传感器的加速度时程 曲线进行复杂度寻踪算法处理，分离出单模态响应信号，并计算其功 率谱(图14)，辨识单模态信号的占主导模态的频率，并利用张紧 弦理论计算索力。

步骤4：逐点滑动3秒的时间窗，重复步骤3，得到频率随时间变 化的曲线(图15)，进而得到索力时程曲线(图16-17)。试验拉索2个 工况的监测索力时程曲线与辨识索力时程曲线很好的吻合在一起，从 试验角度验证了本发明所提算法的准确性。