当前位置:网站首页>[GNSS data processing] Helmert variance component estimation analysis and code implementation
[GNSS data processing] Helmert variance component estimation analysis and code implementation
2022-07-06 12:53:00 【Proletarians】
One 、 background
helmert Variance component estimation is variance - A special case of covariance component estimation , In modern surveying adjustment , There are obvious changes in the data we process , The data sequence includes gross errors , The data type is no longer single . When we encounter two or more types of data , Consider using Helmert variance component estimation , The premise of application is that different types of data need to be independent , So our weight matrix is diagonal matrix .
Two 、 introduce
Let's first introduce the game of King glory , The king is a sharp weapon for flirting with his sister , ha-ha , It's like you have two mobile phones ( Android and apple ), When I take my sister, I probably use apples , Because Apple is better than Android in terms of game experience , Don't drop the frame when taking a sister , No, carton , Don't back off , To operate, there are operations , There is skin, there is skin , This is to give apple a high evaluation , In other words , Give Apple more weight , This is a priori .
that , In the field of Surveying and Mapping , When we measure ( Although I am surveying and mapping , But I am a thief in traditional surveying and mapping ), When measuring the corner of triangulation , There are two types of data ; stay gnss in , When processing multimode data , for example gps and bds There are also two types of data .
3、 ... and 、matlab Code implementation
Data use 《 Generalized survey adjustment 》( The second edition )P106 Example , I use matlab The code implements a wave .
I believe it is the people who have read the theory , Briefly summarize the data processing process :
(1) A priori weighted adjustment ,
(2) Variance component estimation ,
(3) Then reset the weight adjustment ,
(4) When the processed unit weight variance is equal to or less than the threshold .
The formula is not posted , Just post the code
close all
clear all
clc
% function :helmert Variance component estimation
% step1: According to the error equation (v=Bx-l) Determine the constant matrix
% Angle measurement coefficient matrix 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];
% Edge coefficient matrix 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];
% coefficient matrix B 18*4
B=[B1;B2];
% Angle measurement constant term 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;
% Boundary constant term l2 6*1
l2=[ -0.84 1.54 -3.93 2.15 -12.58 -8.21]';
l2=-l2;
% Constant term l 18*1
l=[l1;l2];
% Measure the number of angles and sides
n1=12;
n2=6;
% step2: Classify according to the source of observations , Fixed right
cnt=0;
% Initial weighting , Mean square error of angle measurement =1.5 Mean square error of side measurement =2.0, Save the angle measurement and side measurement variance
sita=[1.5*1.5 2.0*2.0];
while(1)
% Fixed right , Get the power matrix
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)]);
% seek N
N=B'*P*B;
N1=B1'*P1*B1;
N2=B2'*P2*B2;
% seek W
W=B'*P*l;
W1=B1'*P1*l1;
W2=B2'*P2*l2;
% Adjustment
x=inv(N)*W;
v1=B1*x-l1;
v2=B2*x-l2;
% according to hermet The estimation formula determines the matrix
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;% Redefining
cnt=cnt+1;% Used for statistics helmert Number of executions
disp(['【0】helmet Number of times ',num2str(cnt),' 【1】 Unit weight variance of angle measurement and side measurement ',num2str(sita(1)),' ',num2str(sita(2)),'',' 【2】 Ratio of unit weight variance 1:',num2str(sita(2)/sita(1))]);
if(abs(1-min(sita)/max(sita))<0.01)
break;
else
% Angle measurement variance and side measurement variance
sita(1)=sita(1)*inv(sita1);
sita(2)=sita(2)*inv(sita2);
end
end
disp('finish');
Code comments all have , The process is also written , Then put the results I printed
The example in the book is iteration 3 Time out , My is 4 Time . My right does not take two decimal places , Even if you take the last two , You will find in the book N Erroneous , Ha ha ha !
边栏推荐
- Teach you to release a DeNO module hand in hand
- Unity scene jump and exit
- FairyGUI增益BUFF数值改变的显示
- Fabrication d'un sac à dos simple fairygui
- Affichage du changement de valeur du Buff de gain de l'interface graphique de défaillance
- 【rtklib】在rtk下使用抗差自适应卡尔曼滤波初步实践
- Matlab读取GNSS 观测值o文件代码示例
- Unity3d makes the registration login interface and realizes the scene jump
- Special palindromes of daily practice of Blue Bridge Cup
- Particle system for introduction to unity3d Foundation (attribute introduction + case production of flame particle system)
猜你喜欢
[algorithm] sword finger offer2 golang interview question 1: integer division
KF UD分解之UD分解基础篇【1】
[算法] 剑指offer2 golang 面试题2:二进制加法
The master of double non planning left the real estate company and became a programmer with an annual salary of 25W. There are too many life choices at the age of 25
Excel导入,导出功能实现
[algorithm] sword finger offer2 golang interview question 5: maximum product of word length
RTKLIB: demo5 b34f.1 vs b33
Combination of fairygui check box and progress bar
闇の連鎖(LCA+树上差分)
SVN更新后不出现红色感叹号
随机推荐
Prove the time complexity of heap sorting
【干货】提升RTK模糊度固定率的建议之周跳探测
Latex learning
[899] ordered queue
It has been solved by personal practice: MySQL row size too large (> 8126) Changing some columns to TEXT or BLOB or using ROW_ FORMAT
MySQL shutdown is slow
[algorithm] sword finger offer2 golang interview question 2: binary addition
Fairygui joystick
MySQL error warning: a long semaphore wait
堆排序【手写小根堆】
Force buckle 1189 Maximum number of "balloons"
Derivation of logistic regression theory
[算法] 剑指offer2 golang 面试题13:二维子矩阵的数字之和
GPS高程拟合抗差中误差的求取代码实现
Knowledge system of digital IT practitioners | software development methods -- agile
Meanings and differences of PV, UV, IP, VV, CV
[算法] 剑指offer2 golang 面试题5:单词长度的最大乘积
Easy to use shortcut keys in idea
FairyGUI簡單背包的制作
Matlab读取GNSS 观测值o文件代码示例