【GNSS数据处理】赫尔默特(helmert)方差分量估计解析及代码实现

一、背景
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有错误，哈哈哈！