[GNSS data processing] Helmert variance component estimation analysis and code implementation

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 ！