format long;
clear all;
load XT.txt;
t = XT(:,1);
L = XT(:,2)*1000;
ab = polyfit(t,L,1); %线性拟合求a和b的初始值
X0=[ab(2);ab(1);3.0;-pi/3;3.0;-pi/3]; %初始近似值a和b为拟合得到 其余为经验值
P=eye(length(L));
[QXk_1,Xk_1]=LS(P,X0,t,L); %第一次最小二乘平差计算
[Q,X]=XuG_LS(P,QXk_1,Xk_1,t,L); %根据第一次平差结果做序贯平差运算
for i = 1: length(t)
L_new(i,1)=L(i)-X(3)*sin(2*pi*t(i)+X(4))+X(5)*sin(4*pi*t(i)+X(6)); %计算扣除季节性非线性变化后的时间序列
L_m(i,1)=X(1)+X(2)*t(i)+X(3)*sin(2*pi*t(i)+X(4))+X(5)*sin(4*pi*t(i)+X(6)); %建模得到序列
end
%绘图原始时间序列
figure;
plot(t,L,'o');
grid on
xlabel('时间/y');
ylabel('X坐标');
%绘图 扣除季节性变化的时间序列
figure;
plot(t,L_new,'*');
grid on
xlabel('时间/y');
ylabel('X坐标');
%绘图 建模序列与实际序列差值
L_L=L-L_m;
figure;
plot(t,L_L,'.');
grid on
xlabel('时间/y');
ylabel('坐标差值');
fprintf('---X------Sigma---(mm)\n');
for i=1:length(X)
fprintf('%8.4f %8.4f\n',X(i),Q(i,i)); %输出X的平差值以及对应的中误差
end
fprintf('线性速度: v=%8.4f mm/a\n中误差: σ=%8.4f mm\n',X(2),Q(2,2));
function [QX,X]=LS(P,X0,t,L) %LS 最小二乘平差函数
A = zeros(length(t),6);
%X0 =[a0;b0;A10;fi10;A20;fi20]
for i = 1: length(t)
A(i,1) = 1;
A(i,2) = t(i);
A(i,3) = sin(2*pi*t(i)+X0(4));
A(i,4) = X0(3)*cos(2*pi*t(i)+X0(4));
A(i,5) = sin(4*pi*t(i)+X0(6));
A(i,6) = X0(5)*cos(4*pi*t(i)+X0(6));
l(i,1)= L(i)-(X0(1)+X0(2)*t(i)+X0(3)*sin(2*pi*t(i)+X0(4))+X0(5)*sin(4*pi*t(i)+X0(6)));
end
N=inv(A'*P*A);
X=N*(A'*P*l)+X0;
QX=N;
end
function [Q,X]=XuG_LS(P,QX,X0,t,L)
A = zeros(length(t),6);
for i = 1: length(t)
A(i,1) = 1;
A(i,2) = t(i);
A(i,3) = sin(2*pi*t(i)+X0(4));
A(i,4) = X0(3)*cos(2*pi*t(i)+X0(4));
A(i,5) = sin(4*pi*t(i)+X0(6));
A(i,6) = X0(5)*cos(4*pi*t(i)+X0(6));
l_add(i)= L(i)-(X0(1)+X0(2)*t(i)+X0(3)*sin(2*pi*t(i)+X0(4))+X0(5)*sin(4*pi*t(i)+X0(6)));
end
l=l_add';
Qvk=inv(P)+A*QX*A';
M=A*QX;
J=M'*inv(Qvk);
X=X0+J*l;
Q=QX-J*M;
U=l'*P*l;
u=sqrt(U/length(t));
Q=u*Q;
end %Su_LS 序贯最小二乘平差函数 输出为协方差矩阵
没有合适的资源?快使用搜索试试~ 我知道了~
资源推荐
资源详情
资源评论
收起资源包目录
vel.rar (1个子文件)
vel.m 2KB
共 1 条
- 1
资源评论
小贝德罗
- 粉丝: 70
- 资源: 1万+
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功