Helmert方差分量估计
在融合多源异构信息时，不同观测值应具有不同的权阵。但是往往在数据处理之前，我们不知道不同观测值之间的关系，进而无法对不同的观测值进行定权。举个栗子，在GNSS数据处理时，不同卫星系统因其星座设计、观测频率设计导致观测值噪声不同，对于不同的卫星系统的观测值应设计不同的权。即使我们有一些先验信息可以对权阵进行设计，但是随着观测环境的变化会导致先验权不能准确地反映观测值之间的关系，Helmert后验方差分量估计方法更新权阵弥补了这点的不足。
接下来，本篇博文将以两类不同观测值的Helmert方差分量估计作为粒子，推导Helmert背后的原理及其应用流程。
间接平差模型
间接平差函数模型和随机模型为
L
=
B
X
+
Δ
(1)
L=BX+\Delta \tag{1}
L=BX+Δ(1)
D
(
L
)
=
D
(
Δ
)
=
σ
2
P
−
1
(2)
D(L)=D(\Delta)=\sigma^2P^{-1} \tag{2}
D(L)=D(Δ)=σ2P−1(2)
其中，
L
L
L为观测值，
B
B
B为系数矩阵。
误差方程为
V
=
B
X
^
−
L
(3)
V=B\hat X-L \tag{3}
V=BX^−L(3)
其中，
X
^
\hat X
X^为未知状态，则
V
V
V为观测残差，(3)的法方程为
B
T
P
B
X
^
−
B
T
P
L
=
N
X
^
−
W
=
0
(4)
B^TPB\hat X-B^TPL=N\hat X-W=0 \tag{4}
BTPBX^−BTPL=NX^−W=0(4)
X
^
\hat X
X^解为
X
^
=
N
−
1
W
(5)
\hat X=N^{-1}W \tag{5}
X^=N−1W(5)
Helmert方差分量估计原理
方差分量估计的原理：不同观测值的单位权方差相等。
假设有两类相互独立的观测值，则其函数模型为
L
1
=
B
1
X
+
Δ
1
(5)
L_1=B_1X+\Delta_1 \tag{5}
L1=B1X+Δ1(5)
L
2
=
B
2
X
+
Δ
2
(6)
L_2=B_2X+\Delta_2 \tag{6}
L2=B2X+Δ2(6)
其随机模型为
D
(
L
1
)
=
D
(
Δ
1
)
=
σ
1
2
P
1
−
1
(7)
D(L_1)=D(\Delta_1)=\sigma_1^2P_1^{-1} \tag{7}
D(L1)=D(Δ1)=σ12P1−1(7)
D
(
L
2
)
=
D
(
Δ
2
)
=
σ
2
2
P
2
−
1
(8)
D(L_2)=D(\Delta_2)=\sigma_2^2P_2^{-1} \tag{8}
D(L2)=D(Δ2)=σ22P2−1(8)
两类观测值的误差方程可以写为
V
1
=
B
1
X
^
−
L
1
(9)
V_1=B_1\hat X-L_1 \tag{9}
V1=B1X^−L1(9)
V
2
=
B
2
X
^
−
L
2
(10)
V_2=B_2\hat X-L_2 \tag{10}
V2=B2X^−L2(10)
对应权阵分别为
P
1
P_1
P1和
P
2
P_2
P2，合并成矩阵形式
[
V
1
V
2
]
=
[
B
1
B
2
]
X
^
−
[
L
1
L
2
]
(11)
\begin {bmatrix} V_1 \\ V_2 \end {bmatrix} = \begin {bmatrix} B_1 \\ B_2 \end {bmatrix}\hat X - \begin {bmatrix} L_1 \\ L_2 \end {bmatrix} \tag{11}
[V1V2]=[B1B2]X^−[L1L2](11)
由于两类观测值互相独立，计算法方程
(
B
1
T
P
1
B
1
+
B
2
T
P
2
B
2
)
X
^
−
(
B
1
T
P
1
L
1
+
B
2
T
P
2
L
2
)
=
0
(12)
(B_1^TP_1B_1+B_2^TP_2B_2)\hat X-(B_1^TP_1L_1+B_2^TP_2L_2) =0\tag{12}
(B1TP1B1+B2TP2B2)X^−(B1TP1L1+B2TP2L2)=0(12)
其中，
N
=
N
1
+
N
2
=
B
1
T
P
1
B
1
+
B
2
T
P
2
B
2
N=N1+N2=B_1^TP_1B_1+B_2^TP_2B_2
N=N1+N2=B1TP1B1+B2TP2B2，
W
=
B
1
T
P
1
L
1
+
B
2
T
P
2
L
2
W=B_1^TP_1L_1+B_2^TP_2L_2
W=B1TP1L1+B2TP2L2，解得
X
^
=
N
−
1
W
\hat X=N^{-1}W
X^=N−1W。
按照不同观测值单位权方差相等的原则
σ
1
2
=
σ
2
2
\sigma_1^2=\sigma_2^2
σ12=σ22，需要建立观测残差
V
i
T
P
i
V
i
V_i^TP_iV_i
ViTPiVi和单位权方差之间的关系。根据二次型期望公式
E
(
V
1
T
P
1
V
1
)
=
t
r
(
P
1
D
(
V
1
)
)
(13)
E(V_1^TP_1V_1)=tr(P_1D(V_1))\tag{13}
E(V1TP1V1)=tr(P1D(V1))(13)
又由于
V
1
=
B
1
X
^
−
L
1
=
B
1
N
−
1
W
−
L
1
=
B
1
N
−
1
(
W
1
+
W
2
)
−
L
1
=
B
1
N
−
1
(
B
1
T
P
1
L
1
+
B
2
T
P
2
L
2
)
−
L
1
=
(
B
1
N
−
1
B
1
T
P
1
−
E
)
L
1
+
B
1
N
−
1
B
2
T
P
2
L
2
(14)
\begin{aligned} V_1&=B_1\hat X-L_1 \\&=B_1N^{-1}W-L_1 \\&=B_1N^{-1}(W_1+W_2)-L_1 \\&=B_1N^{-1}(B_1^TP_1L_1+B_2^TP_2L_2)-L_1 \\&=(B_1N^{-1}B_1^TP_1-E)L_1+B_1N^{-1}B_2^TP_2L_2 \end{aligned} \tag{14}
V1=B1X^−L1=B1N−1W−L1=B1N−1(W1+W2)−L1=B1N−1(B1TP1L1+B2TP2L2)−L1=(B1N−1B1TP1−E)L1+B1N−1B2TP2L2(14)
根据协方差传播律，得到
V
1
V_1
V1的方差
D
(
V
1
)
=
(
B
1
N
−
1
B
1
T
P
1
−
E
)
D
L
1
(
B
1
N
−
1
B
1
T
P
1
−
E
)
T
+
B
1
N
−
1
B
2
T
P
2
D
L
2
P
2
B
2
N
−
1
B
1
T
=
σ
1
2
(
B
1
N
−
1
N
1
N
−
1
B
1
T
−
2
B
1
N
−
1
B
1
T
+
P
1
−
1
)
+
σ
2
2
B
1
N
−
1
N
2
N
−
1
B
1
T
(15)
\begin{aligned} D(V_1)&=(B_1N^{-1}B_1^TP_1-E)D_{L_1}(B_1N^{-1}B_1^TP_1-E)^T+B_1N^{-1}B_2^TP_2D_{L_2}P_2B_2N^{-1}B_1^T \\&=\sigma_1^2(B_1N^{-1}N_1N^{-1}B_1^T-2B_1N^{-1}B_1^T+P_1^{-1})+\sigma_2^2B_1N^{-1}N_2N^{-1}B_1^T \end{aligned} \tag{15}
D(V1)=(B1N−1B1TP1−E)DL1(B1N−1B1TP1−E)T+B1N−1B2TP2DL2P2B2N−1B1T=σ12(B1N−1N1N−1B1T−2B1N−1B1T+P1−1)+σ22B1N−1N2N−1B1T(15)
将(15)代入(13)得
E
(
V
1
T
P
1
V
1
)
=
t
r
(
P
1
D
(
V
1
)
)
=
σ
1
2
t
r
(
P
1
B
1
N
−
1
N
1
N
−
1
B
1
T
−
2
P
1
B
1
N
−
1
B
1
T
+
P
1
P
1
−
1
)
+
σ
2
2
t
r
(
P
1
B
1
N
−
1
N
2
N
−
1
B
1
T
)
(16)
\begin{aligned} E(V_1^TP_1V_1)&=tr(P_1D(V_1)) \\&=\sigma_1^2tr(P_1B_1N^{-1}N_1N^{-1}B_1^T-2P_1B_1N^{-1}B_1^T+P_1P_1^{-1})+\sigma_2^2tr(P_1B_1N^{-1}N_2N^{-1}B_1^T) \end{aligned}\tag{16}
E(V1TP1V1)=tr(P1D(V1))=σ12tr(P1B1N−1N1N−1B1T−2P1B1N−1B1T+P1P1−1)+σ22tr(P1B1N−1N2N−1B1T)(16)
考虑矩阵迹的性质，(16)可以写为
E
(
V
1
T
P
1
V
1
)
=
σ
1
2
[
n
1
−
2
t
r
(
N
1
N
−
1
)
+
t
r
(
N
1
N
−
1
)
2
]
+
σ
2
2
t
r
(
N
1
N
−
1
N
2
N
−
1
)
(17)
E(V_1^TP_1V_1)=\sigma_1^2[n_1-2tr(N_1N^{-1})+tr(N_1N^{-1})^2]+\sigma_2^2tr(N_1N^{-1}N_2N^{-1}) \tag{17}
E(V1TP1V1)=σ12[n1−2tr(N1N−1)+tr(N1N−1)2]+σ22tr(N1N−1N2N−1)(17)
同理可得
E
(
V
2
T
P
2
V
2
)
=
σ
1
2
t
r
(
N
1
N
−
1
N
2
N
−
1
)
+
σ
2
2
[
n
2
−
2
t
r
(
N
2
N
−
1
)
+
t
r
(
N
2
N
−
1
)
2
]
(18)
E(V_2^TP_2V_2)=\sigma_1^2tr(N_1N^{-1}N_2N^{-1})+ \sigma_2^2[n_2-2tr(N_2N^{-1})+tr(N_2N^{-1})^2]\tag{18}
E(V2TP2V2)=σ12tr(N1N−1N2N−1)+σ22[n2−2tr(N2N−1)+tr(N2N−1)2](18)
去除(17)和(18)中的期望，并将单位权方差
σ
i
2
\sigma_i^2
σi2改写为估计值
σ
^
i
2
\hat \sigma_i^2
σ^i2
[
n
1
−
2
t
r
(
N
1
N
−
1
)
+
t
r
(
N
1
N
−
1
)
2
]
σ
^
1
2
+
t
r
(
N
1
N
−
1
N
2
N
−
1
)
σ
^
2
2
=
V
1
T
P
1
V
1
t
r
(
N
1
N
−
1
N
2
N
−
1
)
σ
^
1
2
+
[
n
2
−
2
t
r
(
N
2
N
−
1
)
+
t
r
(
N
2
N
−
1
)
2
]
σ
^
2
2
=
V
2
T
P
2
V
2
(19)
\begin{aligned} &[n_1-2tr(N_1N^{-1})+tr(N_1N^{-1})^2]\hat \sigma_1^2+tr(N_1N^{-1}N_2N^{-1})\hat \sigma_2^2=V_1^TP_1V_1\\ &tr(N_1N^{-1}N_2N^{-1})\hat \sigma_1^2+ [n_2-2tr(N_2N^{-1})+tr(N_2N^{-1})^2]\hat \sigma_2^2=V_2^TP_2V_2 \end{aligned}\tag{19}
[n1−2tr(N1N−1)+tr(N1N−1)2]σ^12+tr(N1N−1N2N−1)σ^22=V1TP1V1tr(N1N−1N2N−1)σ^12+[n2−2tr(N2N−1)+tr(N2N−1)2]σ^22=V2TP2V2(19)
矩阵形式写为
S
θ
^
=
W
(20)
S\hat \theta=W\tag{20}
Sθ^=W(20)
θ
^
=
S
−
1
W
(21)
\hat \theta=S^{-1}W\tag{21}
θ^=S−1W(21)
式中
θ
^
=
[
σ
1
2
σ
2
2
]
\hat\theta=\begin {bmatrix} \sigma_1^2 \\\sigma_2^2 \end {bmatrix}
θ^=[σ12σ22]，
W
=
[
V
1
T
P
1
V
1
V
2
T
P
2
V
2
]
W=\begin {bmatrix} V_1^TP_1V_1 \\ V_2^TP_2V_2 \end {bmatrix}
W=[V1TP1V1V2TP2V2]，
S
=
[
n
1
−
2
t
r
(
N
1
N
−
1
)
+
t
r
(
N
1
N
−
1
)
2
t
r
(
N
1
N
−
1
N
2
N
−
1
)
t
r
(
N
1
N
−
1
N
2
N
−
1
)
n
2
−
2
t
r
(
N
2
N
−
1
)
+
t
r
(
N
2
N
−
1
)
2
]
S=\left[\begin{matrix} n_1-2tr(N_1N^{-1})+tr(N_1N^{-1})^2 & tr(N_1N^{-1}N_2N^{-1})\\tr(N_1N^{-1}N_2N^{-1})&n_2-2tr(N_2N^{-1})+tr(N_2N^{-1})^2 \end{matrix}\right]
S=[n1−2tr(N1N−1)+tr(N1N−1)2tr(N1N−1N2N−1)tr(N1N−1N2N−1)n2−2tr(N2N−1)+tr(N2N−1)2]
方差分量估计应用流程
1）根据先验信息构建不同观测值之间的关系，即单位权方差，可设置两类观测值单位权方差相等，在同一类观测值内部采用相同的权阵构建方法。如GNSS定位过程中对于GPS和BDS权的计算采用高度角方法。
2）参数估计（最小二乘、卡尔曼滤波）
X
^
\hat X
X^，并计算残差
V
V
V；
3）通过式（21）计算单位权方差；
4）更新权阵：
P
1
=
c
σ
1
2
P
1
P_1=\frac{c}{\sigma_1^2}P_1
P1=σ12cP1，
P
2
=
c
σ
2
2
P
2
P_2=\frac{c}{\sigma_2^2}P_2
P2=σ22cP2，其中c为常数，往往取
c
=
σ
1
2
c=\sigma_1^2
c=σ12。（举个例子，在GPS和BDS融合定位时，将GPS卫星观测值作为基准，调节BDS的权重）
5）若
σ
1
2
=
σ
2
2
\sigma_1^2=\sigma_2^2
σ12=σ22(理论是严格等于，实际运算中约等即可)，则中止迭代，并将更新后的权阵纳入到参数估计中重新结算，否则重复步骤2）-4）至迭代超出次数或单位权方差相等。
参考文献
[1] 崔希璋, 於宗俦, 陶本藻, 等. 广义测量平差[Ｍ]. 武汉: 武汉大学出版社, 2009.
[2] 刘金海, 涂锐, 张睿, 张鹏飞, 卢晓春. Helmert方差分量估计在GPS/GLONASS/BDS组合定位权比确定中的应用[J]. 大地测量与地球动力学, 2018, 38(06): 568-570+576.
扫一扫
2255
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
