基本思想

先对各类观测值定初权，进行预平差，利用预平差得到的相应各类观测值的残差平方和，依据一定的原则对各类观测值的验前方差作出估计，重新定权。通过这样不断计算迭代，使得不同类观测值的单位权中误差趋向一致，从而达到最佳平差效果。

赫尔默特方差分量估计的计算步骤（二元为例）

实际应用

有一组模拟的GPS和GLONASS观测数据，采用伪距单历元单点定位，根据经验，ＧＰＳ的伪距中误差为０.３ｍ，GLONASS的伪距中误差为０.８ｍ，已给出误差观测方程的Ａ阵和Ｌ阵，试按间接平差方差进行验后Helmert方差分量估计，并求待定点的坐标值。

import numpy as np

# 设初始权
# 令GPS伪距中误差（0.3m）为单位权中误差
p1_0 = np.eye(10)  # GPS权阵(10*10)
p2_0i = (0.3 / 0.8) ** 2
array_1 = np.array([p2_0i] * 9)
p2_0 = np.diag(array_1)  # glonass权阵(9*9)

# 读取数据
data = np.genfromtxt(r'C:\data.csv', delimiter=',')
B_GPS = data[:10, :-1]  # GPS系数阵
L_GPS = data[:10, -1]  # GPS L阵
B_GLONASS = data[10:, :-1]  # GLONASS系数阵
L_GLONASS = data[10:, -1]  # GLONASS L阵


def cal_canshu(B, L, P):
    # 计算Ni,Wi
    Ni = np.dot(B.T, P).dot(B)  # Ni=BT*P*B
    Wi = np.dot(B.T, P).dot(L)  # Wi=BT*P*L
    return Ni, Wi


def cal_(B_GPS, L_GPS, p1_0, B_GLONASS, L_GLONASS, p2_0):
    # 法方程：
    # N1*x=W1,N2*x=W2,即(N1+N2)*x=W1+W2
    # 即：Nx=W
    # 计算X,V,VTPV
    N1, W1, = cal_canshu(B_GPS, L_GPS, p1_0)
    N2, W2, = cal_canshu(B_GLONASS, L_GLONASS, p2_0)
    N = N1 + N2
    W = W1 + W2
    x = np.dot(np.linalg.inv(N), W)  # 参数改正数x=inv(N)*W
    v_GPS = np.dot(B_GPS, x) - L_GPS  # gps改正数v=b_gps*x-l
    v_GLONASS = np.dot(B_GLONASS, x) - L_GLONASS  # glonass改正数v=b_glonass*x-l
    V1TPV1 = np.dot(v_GPS.T, p1_0).dot(v_GPS)  # vt*p*v
    V2TPV2 = np.dot(v_GLONASS, p2_0).dot(v_GLONASS)
    return N, N1, N2, V1TPV1, V2TPV2, x


def Helmert(N, N1, N2, V1TPV1, V2TPV2, n1, n2):
    #     helmert方差分量计算
    inv_N = np.linalg.inv(N)
    s2 = s3 = np.trace(inv_N.dot(N1).dot(inv_N).dot(N2))  # s2=tr(inv(N)*N1*inv(N)*N2)
    s1 = n1 - 2 * np.trace((np.dot(inv_N, N1))) + np.trace(
        np.dot(inv_N, N1).dot(inv_N).dot(N1))  # s1=n1-2*tr(inv(N)*N1)+tr[(inv(N)*n1)**2]
    s4 = n2 - 2 * np.trace((np.dot(inv_N, N2))) + np.trace(np.dot(inv_N, N2).dot(inv_N).dot(N2))
    S = np.array([s1, s2, s3, s4]).reshape(2, 2)  # S矩阵
    M = np.array([V1TPV1, V2TPV2]).T  # M矩阵
    sigma = np.dot(np.linalg.inv(S), M)  # 方差估计值=inv(S)*M
    sigma[0] = round(sigma[0], 4)  # 保留五位小数
    sigma[1] = round(sigma[1], 4)
    return sigma, S, M

if __name__ == '__main__':
    N, N1, N2, V1TPV1, V2TPV2, x = cal_(B_GPS, L_GPS, p1_0, B_GLONASS, L_GLONASS, p2_0)
    sigma, S, M = Helmert(N, N1, N2, V1TPV1, V2TPV2, 10, 9)
    k = 1  # 迭代次数
    row = []
    while sigma[0] != sigma[1]:
        m = sigma[0] / sigma[1]  # 权阵转化因子m：设m=sigma[0]
        p2_0 = p2_0 * m
        # 或者令m=1
        # p1_0*=1/sigma[0]
        # p2_0*=1/sigma[1]
        N, N1, N2, V1TPV1, V2TPV2, x = cal_(B_GPS, L_GPS, p1_0, B_GLONASS, L_GLONASS, p2_0)
        sigma, S, M = Helmert(N, N1, N2, V1TPV1, V2TPV2, 10, 9)
        k += 1
        print('S:', S, 'M:', M, 'sigma:', sigma, 'sigma比值：', sigma[0] / sigma[1], '迭代次数k:%d' % (k),
              '参数改正数(坐标改正数，钟差改正数)：', x,
              sep='')

结果

S:[[6.07339446 0.45040852]
 [0.45040852 8.02578849]]
M:[6.54536641 8.50379278]
 sigma:[1.0033 1.0033]
sigma比值：1.0
 迭代次数k:10
 参数改正数(坐标改正数，钟差改正数)：
[0.03935843  1.04245476  0.28245052  0.52865248]