当前位置:网站首页>【GNSS数据处理】赫尔默特(helmert)方差分量估计解析及代码实现
【GNSS数据处理】赫尔默特(helmert)方差分量估计解析及代码实现
2022-07-06 09:18:00 【Proletarians】
一、背景
helmert方差分量估计是方差-协方差分量估计的一个特例,在近现代测量平差中,我们处理的数据有明显的变化,数据序列中包括粗差,数据类型不再单一等问题。当我们遇到数据类型两类及以上时,可以考虑使用赫尔默特方差分量估计,应用的前提条件是不同类别的数据需要相互独立,这样我们的权阵就是对角阵。
二、引入
我们先用王者荣耀这款游戏引入一下,王者可是撩妹利器,哈哈,就像你有两部手机(安卓和苹果),带妹的时候大概率会用苹果,因为游戏体验方面苹果要优于安卓,带妹的时候不掉帧,不卡顿,不闪退,要操作有操作,有皮肤有皮肤,这就是给苹果的评价要高,换言之,给苹果的权重要大,这是先验的。
那么,在测绘领域,我们在进行测量时(虽然我是测绘的,但是我传统测绘方面贼菜),进行三角网边角测量时,就是两类数据;在gnss中,进行多模数据处理时,例如gps和bds也是两类数据。
三、matlab代码实现
数据采用《广义测量平差》(第二版)P106的例子,我用matlab代码实现了一波。
相信都是看了理论的人,简单总结下数据处理流程:
(1)先验定权平差,
(2)方差分量估计,
(3)再重新定权平差,
(4)当处理后的单位权方差相等或小于阈值即可。
公式就不贴了,直接贴代码了
close all
clear all
clc
% 功能:helmert方差分量估计
% step1:根据误差方程(v=Bx-l)确定常量矩阵
% 测角系数矩阵B1 12*4
B1=[0.5532 -0.8100 0 0;0.2434 0.5528 0 0;-0.7966 0.2572 0 0;-0.2434 -0.5528 0 0;
0.6298 0.6368 0 0;-0.3864 -0.0840 0 0;0.7966 -0.2572 -0.2244 -0.3379;-0.8350 -0.1523 0.0384 0.4095;
0.0384 0.4095 0.1860 -0.0716;-0.0384 -0.4095 0.2998 0.1901;-0.3480 0.3255 -0.0384 -0.4095;0.3864 0.0840 -0.2614 0.2194];
% 测边系数矩阵B2 6*4
B2=[0.3072 0.9516 0 0;-0.9152 0.4030 0 0;0.2124 -0.9722 0 0;
0 0 -0.6429 -0.7660;0 0 -0.8330 0.5532;0.9956 -0.0934 -0.9956 0.0934];
% 系数矩阵B 18*4
B=[B1;B2];
% 测角常数项l1 12*1
l1=[0.18 -0.53 3.15 0.23 -2.44 1.01 2.68 -4.58 2.80 -3.10 8.04 -1.14]';
l1=-l1;
% 测边常数项l2 6*1
l2=[ -0.84 1.54 -3.93 2.15 -12.58 -8.21]';
l2=-l2;
% 常数项l 18*1
l=[l1;l2];
% 测角和测边个数
n1=12;
n2=6;
% step2:按观测值来源分类,定权
cnt=0;
% 初始定权,测角中误差=1.5 测边中误差=2.0,保存测角和测边方差
sita=[1.5*1.5 2.0*2.0];
while(1)
% 定权,得到权阵
sita1=sita(1)/sita(1);
sita2=sita(1)/sita(2);
% sita1=1;
% sita2=0.56;
P1=diag(sita1*ones(1,12));
P2=diag(sita2*ones(1,6));
P=diag([sita1*ones(1,12) sita2*ones(1,6)]);
% 求N
N=B'*P*B;
N1=B1'*P1*B1;
N2=B2'*P2*B2;
% 求W
W=B'*P*l;
W1=B1'*P1*l1;
W2=B2'*P2*l2;
% 平差
x=inv(N)*W;
v1=B1*x-l1;
v2=B2*x-l2;
% 按照hermet估算公式确定矩阵
w1=v1'*P1*v1;
w2=v2'*P2*v2;
w=[w1 w2]';
S=[n1-2*trace(inv(N)*N1)+trace(inv(N)*N1*inv(N)*N1) trace(inv(N)*N1*inv(N)*N2);
trace(inv(N)*N1*inv(N)*N2) n2-2*trace(inv(N)*N2)+trace(inv(N)*N2*inv(N)*N2)];
sita=inv(S)*w;% 重新定权
cnt=cnt+1;%用于统计helmert执行次数
disp(['【0】helmet进行次数 ',num2str(cnt),' 【1】测角和测边单位权方差 ',num2str(sita(1)),' ',num2str(sita(2)),'',' 【2】单位权方差之比 1:',num2str(sita(2)/sita(1))]);
if(abs(1-min(sita)/max(sita))<0.01)
break;
else
% 测角方差和测边方差
sita(1)=sita(1)*inv(sita1);
sita(2)=sita(2)*inv(sita2);
end
end
disp('finish');
代码注释都有,流程也写了,再放上我打印的结果
书上的例子是迭代了3次退出,我的是4次。我的权没有取小数后两位,即使取了后两位,会发现书上N有错误,哈哈哈!
边栏推荐
- Office prompts that your license is not genuine pop-up box solution
- Mixed use of fairygui button dynamics
- CUDA C programming authoritative guide Grossman Chapter 4 global memory
- FGUI工程打包发布&导入Unity&将UI显示出来的方式
- [leetcode15] sum of three numbers
- JS Title: input array, exchange the largest with the first element, exchange the smallest with the last element, and output array.
- HCIP Day 12
- 1041 be unique (20 points (s)) (hash: find the first number that occurs once)
- 编译原理:源程序的预处理及词法分析程序的设计与实现(含代码)
- Force buckle 1189 Maximum number of "balloons"
猜你喜欢
Database course design: college educational administration management system (including code)
Unity3D制作注册登录界面,并实现场景跳转
Classification, understanding and application of common methods of JS array
基於Redis的分布式ID生成器
Office prompts that your license is not genuine pop-up box solution
Fairygui gain buff value change display
FairyGUI复选框与进度条的组合使用
Problèmes avec MySQL time, fuseau horaire, remplissage automatique 0
MySQL takes up too much memory solution
There is no red exclamation mark after SVN update
随机推荐
Combination of fairygui check box and progress bar
FGUI工程打包发布&导入Unity&将UI显示出来的方式
Derivation of logistic regression theory
Learning notes of JS variable scope and function
Problèmes avec MySQL time, fuseau horaire, remplissage automatique 0
燕山大学校园网自动登录问题解决方案
How to improve the deletion speed of sequential class containers?
Theoretical derivation of support vector machine
Fabrication d'un sac à dos simple fairygui
InnoDB dirty page refresh mechanism checkpoint in MySQL
Agile development helps me
Vulnhub target: hacknos_ PLAYER V1.1
如何给Arduino项目添加音乐播放功能
341. Flatten nested list iterator
Database course design: college educational administration management system (including code)
记一次云服务器被密码爆破的经历——关小黑屋、改密码、改端口
KF UD分解之UD分解基础篇【1】
Fairygui joystick
Knowledge system of digital IT practitioners | software development methods -- agile
Special palindromes of daily practice of Blue Bridge Cup