[GNSS] robust estimation (robust estimation) principle and program implementation

Robust estimation

The principle of robust estimation

Robust estimation is a category of modern measurement adjustment , Also known as robust estimation （robust estimate）, According to academician Yang, the CAS system likes to call it robust estimation , Wu Da likes to call it robust estimation . Our measurements are random variables , According to the normal distribution , If there is gross error （gross error） Words , We will make the result deviate from the true value when we apply the least squares or Kalman filter （ Filter divergence ） The phenomenon of . When we solve gross errors or systematic errors , It can be understood from two aspects , Mean shift or variance inflation , Robust estimation belongs to variance inflation model , That is, the mean value remains unchanged , The phenomenon of variance change . We can reduce the weight of observations with gross errors .
I compare lazy , I don't want to write a formula , You should have seen a similar derivation , In fact, the formula is the same as the minimum second difference , The place of change is the power matrix , Replaced by equivalent weight , Then what affects the equivalent weight is the equivalent weight function , Common are huber、IGG III etc. , I recommend IGG III, The code is easy to implement .
When we do it , This power matrix is very important , The weight matrix is the inverse of what we call the variance matrix , It is divided into independent and non independent , That is, whether the observed values are related , Whether it is relevant also affects the form of the equivalent weight function . Take the simplest example of independent observation distance

Examples of robust estimation

Let's assume that for a distance 10 Sub independent observation ,10 The observed values are ：5.09、5.10、5.13、5.09、5.12、5.08、5.46、7.81、5.10、5.11（ Unit is m）, I made up the data casually .
From the above paragraph , We can know the number of measurements n=10, Necessary observation m=1, The redundant observation is n-m=9. Independent observation, then the weight matrix is diagonal , Prior weight matrix is unit diagonal matrix . And what is visible to the naked eye is 7.81 This observation is a gross error ,5.46 Well, I'm not sure if .

Robust program implementation

matlab Realized

%%  function 3- Robust estimation 
% %  Suppose a length is observed ,5.09 5.10 5.13 5.09 5.12 5.08 5.46 7.81 5.10 5.11
close all
clear all
clc
k0=1.0;k1=2.5;k=1.0;
B=ones(10,1);%  Design matrix 
l=[5.09 5.10 5.13 5.09 5.12 5.08 6.46 7.81 5.10 5.11]';
% P=diag(ones(10,1));%  A priori weight matrix 1, Equivalent weight 
P=diag(1./l);%  A priori weight matrix 2, Weight according to the reciprocal of length 
x_prev=0;
x0=5.10;%  Assign initial value to 
l=l-x0;

while(1)

    W=B'*P*l;
    N=B'*P*B;
    x=inv(N)*W;%  Parameter vector to be evaluated 
    v=B*x-l;%  Residual vector 
    sigma0=sqrt(v'*P*v/(10-1));%  Mean square error of unit weight 
    disp([' Parameter vector to be evaluated ：',num2str(x),'  Mean square error of unit weight ：',num2str(sigma0)]);
    Q=inv(P);
    Qvv=Q-B*inv(N)*B';

    for i=1:10
        v_=v(i)/(sigma0*sqrt(Qvv(i,i)));
        if abs(v_)<=k0
            k=1.0;
        elseif  abs(v_)>k1
            k=1e-8;
        else
            k=(k0/abs(v_))*((k1-abs(v_))/(k1-k0))^2;
        end
        P(i,i)=P(i,i)*k;
    end

    if x_prev==0
        x_prev=x;
        continue;
    elseif abs(x_prev-x)<0.01
        break;
    end
    x_prev=x;
end
x=x0+x_prev;% Find the final length   initial value + Correct number 
disp(['length is ',num2str(x)])

The results are as follows ：
The result of using equivalent weight is ：

 Parameter vector to be evaluated ：0.309  Mean square error of unit weight ：0.8512
 Parameter vector to be evaluated ：0.042222  Mean square error of unit weight ：0.11331
 Parameter vector to be evaluated ：0.0025  Mean square error of unit weight ：0.01472
 Parameter vector to be evaluated ：-0.00031553  Mean square error of unit weight ：0.0084525
length is 5.1025

Using the reciprocal of length to determine the weight, the result is ：

 Parameter vector to be evaluated ：0.2218  Mean square error of unit weight ：0.31128
 Parameter vector to be evaluated ：0.03985  Mean square error of unit weight ：0.048702
 Parameter vector to be evaluated ：0.0024523  Mean square error of unit weight ：0.0065128
 Parameter vector to be evaluated ：-0.0003396  Mean square error of unit weight ：0.0037381
length is 5.1025

The first line in the figure is the result of our least squares ,10 Only one of the results is close to the value obtained by the least square , Obviously, it is inconsistent with our measured results , Therefore, the existence of a gross error will have a devastating effect on the result of the least squares . So if we apply robust estimation , after 3 The final result and unit weight mean square error are obtained in the second iteration .
If I'm on the seventh value , Making a change , by 5.46 become 6.46. The results are shown below ：
The result of using the unit right is ：

 Parameter vector to be evaluated ：0.409  Mean square error of unit weight ：0.91426
 Parameter vector to be evaluated ：0.1257  Mean square error of unit weight ：0.38585
 Parameter vector to be evaluated ：0.0025  Mean square error of unit weight ：0.01472
 Parameter vector to be evaluated ：-0.00031551  Mean square error of unit weight ：0.0084526
length is 5.1025

Using the reciprocal of length to determine the weight, the result is ：

 Parameter vector to be evaluated ：0.2218  Mean square error of unit weight ：0.31128
 Parameter vector to be evaluated ：0.03985  Mean square error of unit weight ：0.048702
 Parameter vector to be evaluated ：0.0024523  Mean square error of unit weight ：0.0065128
 Parameter vector to be evaluated ：-0.0003396  Mean square error of unit weight ：0.0037381
length is 5.1025

Conclusion
The accuracy of parameter estimation is improved by using robust estimation , In addition, the results of different prior weight matrices are inconsistent in the iterative process , For this example, the reciprocal of length is used to determine the weight , Whether or not the minimum two difference is applied , The results are close to the real value , The mean square error of unit weight is also small .

Robust estimation of correlation between observations

The main thing is the three-step iteration ：
1、 Standardized residuals , Mean square error of unit weight
2、 Variance amplification factor
3、 Equivalent variance

Welcome to exchange ～

Updated on 2020/11/10 14:20:00
In the above example , Use standardized residuals , Thank you for reminding me by Mr. Xiao of China Academy of testing , I haven't seen classic adjustment for a long time , The covariance matrix propagation formula is really forgotten , Ha ha ha .
Expand ： In fact, we are gnss spp in , There is also no correlation between the observations , So you see , Find the corresponding parameters and apply them directly .