桥梁施工监控卡尔曼滤波数据处理源代码

2022-12-17 09:53 刘泉

主跨195米连续刚构施工监控项目,想如何进行监测数据的科学处理,为合理的桥梁设计线形保驾护航,比较了用“最小二乘法”、“卡尔曼滤波法”来处理监测数据的优劣。决定用”自适应卡尔曼滤波法”,使用Matlab软件编了代码,(另“最小二乘法简单,很多期刊上都有,就不在这分享了”)现分享下:
卡尔曼滤波法的原理如图:

微信图片_20230529100157.png

程序代码如下:-----------------------------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------------------------------------------------------
clear;clc;clearvars;%由1#、2#块的数据开始
%预拱度:范指 成桥后,无车辆荷载状况下,桥面标高高于设计标高的值
n=12;
%理论预拱度计算值(由1#块的数据开始)
h_z_y=[23,16,21,26,34,44,50,54,55,50,38,29;15,4,5,6,10,18,22,27,29,27,19,16];
%预拱度实测值Y=(实测立模标高 - 实测张拉后标高 - 挂篮变形值 - 节段挠度计算理论值) + 节段预拱度计算理论值;        理论计算立模标高
%预拱度实测值(由1#块的数据开始)
Sh_z_y=[21,18,27,31,39,47,52,59,63,55,39,31;14,7,10,13,19,30,35,42,48,41,32,24];
%取 1#块两悬臂端理论预拱度与实测预拱度差值的平方为协方差的初始值P1)
P1=[1,0;0,1];
%Q为i#块(左i#块理论立模标高(或理论预拱度) - 实测立模标高(或实测预拱度);右i#块理论立模标高 (或理论预拱度)- 实测立模标高(或实测预拱度))即Q=[。。。;。。。];,Q为系统误差自协方差,它反映的是计算模型与实际施工状态之间的偏差(由1#块的数据开始)
%R为测量值均方差,取从2块开始的(i#预拱度实测值的均方差),即R=[。。。。],根据误差传播定律计算的梁段悬臂长为s时预拱度实测值的均方差σ = k= √2s/1000,从2#块开始~5#号结束(由2#块的数据开始)
R=[0.53,0.77,1.02,1.26,1.50,1.31,1.56,1.36,1.61,1.41,1.66];
Q(1,1)=2;
Q(2,1)=1;
q=[Q(1,1);Q(2,1)];
Q1=[Q(1,1)^2,0;0,Q(2,1)^2];
%遗忘因子b=0.98
b=0.98;
for i =1:n-2
X11=[h_z_y(1,i);h_z_y(2,i)];
Z2=[Sh_z_y(1,i+1);Sh_z_y(2,i+1)];
A2=[h_z_y(1,i+1)/h_z_y(1,i) , 0;   0 , h_z_y(2,i+1)/h_z_y(2,i)] ;    %A2是本段与前一段函数模型状态转移阵
%r=[sqrt(R(i));sqrt(R(i))]; %r=0是本段测量噪声
R2=[R(i),0;0,R(i)];         %R2是本段测量噪声方差阵
P2=A2*P1*A2'+Q1; %P1是前一段方差阵;P2是本段方差阵
K2=P2*inv(P2+R2) ; %K2是本段增益阵
X21=A2*X11+q;   %X11是前一段函数模型值;X21是本段函数模型值
EE=Z2-X21;           % EE是本段残差值;        Z2是本段实测值
X22=X21+K2*EE ;   %X22是本段卡尔曼预测值
P3=(eye(2)-K2)*P2;   %P3是本段方差阵卡尔曼修正更新
%以下为噪声估计器
dk=(1-b)/(1-b^(i+1));
q=(1-dk)*q+dk*(X22-A2*X11);
Q1=(1-dk)*Q1+dk*(K2*EE*EE'*K2'+P3-A2*P1*A2');
%以下为预测未浇筑段卡尔曼滤波预测结果值
A3=[h_z_y(1,i+2)/h_z_y(1,i+1) , 0;   0 , h_z_y(2,i+2)/h_z_y(2,i+1)]; %A3是本段与未浇筑段需卡尔曼预测段模型状态转移阵
X32=A3*X22 %X32是未浇筑段卡尔曼滤波预测结果值
%以下为递推
P1=P2;
P2=P3;
%打印数据
%Yc1(i)=X32(1,1)
%Yc2(i)=X32(2,1)
plot(i,X32(2,1),'-.b*');
plot(i,Sh_z_y(2,i+2),'-.g*');
plot(i,h_z_y(2,i+2),'-.k*');
%plot(i,X32(2,1),'--mo');
hold on
end