赫尔默特方差分量估计

1 赫尔默特方差分量估计
我们知道，平差前观测值向量的方差阵一般是未知的，因此平差时随机模型都是使用观测值向量的权阵。

而权的确定往往都是采用经验定权，也称为随机模型的验前估计，对于同类观测值可按第一章介绍的常用定权方法定权；对于不同类的观测值，就很难合理地确定各类观测值的权。

为了合理地确定不同类观测值的权，可以根据验前估计权进行预平差，用平差后得到的观测值改正数来估计观测值的方差，根据方差的估计值重新进行定权，以改善第一次平差时权的初始值，再依据重新确定的观测值的权再次进行平差，如此重复，直到不同类观测值的权趋于合理，这种平差方法称为验后方差分量估计。

此概念最早由赫尔默特（F.R.Helmert ）在1924年提出，所以又称为赫尔默特方差分量估计。

一、赫尔默特方差分量估计公式
为推导公式简便起见，设观测值由两类不同的观测量组成，不同类观测值之间认为互不相关，按间接平差时的数学模型为
2221
11~
~
∆-=∆-=X B L X B L （函数模型） （8-4-1）
0),(()()()()(212112
2022112011=∆∆==∆==∆=--D L L D P D L D P D L D ）
，σσ （随机模型） （8-4-2）
其误差方程为
111ˆl x
B V -= 权阵1P （8-4-3）
222ˆl x
B V -= 权阵2P （8-4-4） 作整体平差时，法方程为
0ˆ=-W x
N （8-4-5）
式中
2222111121B P B N B P
B N N N N T
T
==+=，， 2222111121l P B W l P
B W W W W T
T ==+=，，
一般情况下，由于第一次给定的权1P 、2P 是不恰当的，或者说它们对应的单位权方差
是不相等的，设为201σ和202σ，则有
1
220221
12
011)()(--==P L D P L D σσ
（8-4-6）
但只有2
0202201σσσ==才认为定权合理。

方差分量估计的目的就是根据事先初定的权1P 、2P 进行预平差，然后利用平差后两类观测值的111V P V T 、222V P V T
来求估计量202201ˆˆσσ、，再根据（8-4-6）式求出)(ˆ)(ˆ21L D L D 、，由这个方差估值再重新定权，再平差，直到202201ˆˆσσ=为止。

为此需要建立111V P V T 、222V P V T
与估计量202201
ˆˆσσ、之间的关系式。

由数理统计知识可知，若有服从任一分布的q 维随机变量1⨯q Y
，已知其数学期望为1⨯q η
，方差阵为q
q ⨯∑
，则1⨯q Y
向量的任一二次型的数学期望可以表达为：
ηηB B tr BY Y E T T +∑=)()(
（8-4-7）
式中B 为任意q 阶的对称可逆阵。

现用V 向量代替上式中的Y 向量，则其中η的应换为)(V E ，∑应换为)(V D ，B 阵可以换成权阵P ，于是有
)()())(()(V PE V E V PD tr PV V E T
T += （8-4-8） 前面已经证明0)(=V E ，于是有：
))(()(111V PD tr PV V E T
= （8-4-9）
而
11
11l W N B V -=- 12111)(l W W N B -+=-
12221111111l l P B N B l P B N B T
T -+=--
2221
111111)(l P B N B l I P B N B T
T --+-= 对上式应用协因数传播律得
+--=--T
T T I P B N B L D I P B N B V D ))(()()(1111111111 T T
B N B P L D P B N B 112222211)(--
将1
22022112011)()(--==P L D P L D σσ、代入上式，整理后得
T
T T B N N N B P B N B B N N N B V D 11
21120211111111112011)2()(------++-=σσ 将上式代入（8-4-9）式，得 ))(()(11111V D P tr V P V E T
=
)
()2(11211120211111111111112
01T T
T
B N N N B P tr P P B N B P B N N N B P tr ------++-=σσ
顾及矩阵迹的性质，上式可写为
)()]()(2[)(1
2112021111111201111-----++-=N N N N tr N N N N tr N N tr n V P V E T
σσ
同理可得
)()]()(2[)(12112011212122202222-----++-=N N N N tr N N N N tr N N tr n V P V E T
σσ 去掉上面两式的期望符号，相应的单位权方差202201σσ、也改用估值符号2
02201ˆˆσσ、表示，整
理顺序后得
1112
0212112011111111ˆ)(ˆ)]()(2[V P V N N N N tr N N N N tr N N tr n T =++------σσ
（8-4-10）
2222
0212121222011211ˆ)]()(2[ˆ)(V P V N N N N tr N N tr n N N N N tr T =+-+-----σσ
（8-4-11）
其矩阵形式可写为
θ
θW S =⨯⨯1
222ˆ （8-4-12） θ
θW S 11
2ˆ-⨯=
（8-4-13）
式中
⎥⎦⎤⎢⎣
⎡+-+-=----------)()(2)()
()()(21212122121112111111111N N N N tr N N tr n N N N N tr N N N N tr N N N N tr N N tr n S
[
]T
202201ˆˆˆσσθ=
[]
T
T
T V P V V P V W 222111=θ （8-4-12）、（8-4-13）两式即为赫尔默特方差分量估计的严密公式。

由此式可以求得两类观测值的单位权方差估值，从而可以根据（8-4-6）式求得观测值方差的估值，以此方差估值再次定权，再次平差，直至满足要求为止。

现将以上推导扩展至m 组观测值。

误差方程为
i i i l x B V -=ˆ )21
(m i ，= 令
1
20)(-=i i i P L D σ
∑==m
i
i i T
i
N N B P B N 1
，
∑==m
i
i i T i
W W l P B W 1
，
则得参数的估值为
W N x
1ˆ-= 按照上述类似的推导，则有
)
()]()(2[)(1120,11
112
0--≠=---∑+
+-=N N N N tr N N N N tr N N tr n V P V E j i j m
i
j j i i i i i
i i T
i σ
σ
去掉期望符号，相应的单位权方差20i σ也改为用估值符号2
0ˆi σ，则有
θ
θW S m m m =⨯⨯1
ˆ （8-4-14）
式中
⎥⎥
⎥
⎥⎥⎦⎤
⎢⎢⎢⎢⎢⎣⎡+-+-+-=---------------------)()(2)()()()()(2)()()()()(21
1112111112
121212122111211112111111111N N N N tr N N tr n N N N N tr N N N N tr N N N N tr N N N N tr N N tr n N N N N tr N N N N tr N N N N tr N N N N tr N N tr n S m m m m m m m ，，，
[]T
m 20202201ˆˆˆˆσσσθ =
[]
T
m m T
m T T V P V V P V V P V W 222111=θ
二、计算步骤
1.将观测值分类，并进行验前权估计，即确定各类观测值的权的初值m P P P 21、；
2.进行第一次平差，求得
i i T i V P V ； 3.按（8-4-14）式求各类观测值单位权方差估值2
0ˆi σ；
4.按（8-4-6）式计算各类观测值方差的估值；
5.依据定权公式再次定权，再次平差，如此反复，直到各类单位权方差的估值相等或接近相等为止。

2 秩亏自由网平差
在前面介绍的经典平差中，都是以已知的起算数据为基础，将控制网固定在已知数据上。

如水准网必须至少已知网中某一点的高程，平面网至少要已知一点的坐标、一条边的边长和一条边的方位角。

当网中没有必要的起算数据时，我们称其为自由网，本节将介绍网中没有起算数据时的平差方法，即自由网平差。

在经典间接平差中，网中具备必要的起算数据，误差方程为
11
1
ˆ⨯⨯⨯⨯-=n t t n n l x
B V （8-2-1）
式中系数阵B 为列满秩矩阵，其秩为t B R =)(。

在最小二乘准则下得到的法方程为
0ˆ1
1
=-⨯⨯⨯t t t
t bb W x
N （8-2-2）
由于其系数阵的秩为
t B R PB B R N R T
bb ===)()()(，所以bb N 为满秩矩阵，即为非奇异阵，具有凯利逆bb N
1
-，因此具有唯一解，即
W N x
bb 1ˆ-= （8-2-3）
当网中无起算数据时，网中所有点均为待定点，设未知参数的个数为u ，误差方程为
11
1
ˆ⨯⨯⨯⨯-=n u u n n l x
B V （8-2-4）
式中
d t u +=
d 为必要的起算数据个数。

尽管增加了d 个参数，但B 的秩仍为必要观测个数，即
u t B R <=)(
其中B 为不满秩矩阵，称为秩亏阵，其秩亏数为d 。

组成法方程
0ˆ1
1
=-⨯⨯⨯u u u u W x
N （8-2-5）
式中
Pl
B W PB B N T u T u
u ==⨯⨯1
，，且
u t B R PB B R N R T <===)()()(，所以N 也为秩亏阵，秩亏数为：
t u d -=
（8-2-6）
由上式知，不同类型控制网的秩亏数就是经典平差时必要的起算数据的个数。

即有：
⎪⎩⎪
⎨⎧=测角网网
测边网、边角网、导线水准网、测站平差,4,3,
1d
在控制网秩亏的情况下，法方程有解但不唯一。

也就是说仅满足最小二乘准则，仍无法
求得x
ˆ的唯一解，这就是秩亏网平差与经典平差的根本区别。

为求得唯一解，还必须增加新的约束条件，来达到求唯一解的目的。

秩亏自由网平差就是在满足最小二乘min =PV V T
和
最小范数min ˆˆ=x x
T 的条件下，求参数一组最佳估值的平差方法。

下面将推导自由网平差常用两种解法的有关计算公式。

一、直接解法
根据广义逆理论，相容方程组0ˆ1
1=-⨯⨯⨯u u u u W x
N 虽然具有无穷多组解，但它有唯一的最小
范数解，即：
W N x m r 1
ˆ-= （8-2-7）
式中--=)(T T m NN N N ，称为矩阵N 的最小范数g 逆。

-)(T NN 称为矩阵T
NN 的g 逆。

代入（8-2-7）式得
W NN N x T T r -
=)(ˆ （8-2-8）
上式就是根据广义逆理论直接求解参数的唯一最小范数解的公式。

由于广义逆计算较为复杂，下面将公式做进一步改化： 令
⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭
⎫
⎝⎛=⨯⨯⨯⨯⨯2122211211
N N N N N N N d d t d d t t
t u
u （8-2-9）
⎪⎪
⎭⎫
⎝⎛=⨯⨯⨯12111
d t u W W W （8-2-10） 式中1N 行满秩，即t N R =)(1，于是有
()⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛=T T
T T
T T
T
N N N N N N N N N N N N NN 221
221112
1
21 （8-2-11）
而t N R N N R T ==)()(111，所以
)(11T
N N 为满秩方阵，按照降阶法求矩阵广义逆的方法，即：如果有矩阵
⎪⎪⎪⎭⎫ ⎝⎛=-⨯-⨯--⨯⨯⨯)()(22)(21)
(1211
r n r m r
r m r n r r r n
m A A A
A A
其中1111)(A r A R ，=存在凯利逆，则有n m A
⨯的g 逆
⎪⎪⎭⎫ ⎝
⎛=⨯-⨯-
00011
1r r n
m A A
（8-2-12）
根据上式可得
⎪⎪⎭⎫
⎝⎛=⎪⎪⎭⎫ ⎝⎛=--
00
0000)()(111
111
1Q N N N N T T
（8-2-13）
代入（8-2-8）式，得
()111121112
1
000ˆW Q N W W Q N N x
T
T T
=⎪
⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛=
（8-2-14）
或写成
11111)(ˆW N N N x
T T -= （8-2-15）
未知参数的协因数阵为：
11111111111111ˆˆ)(11N Q N Q N Q N Q Q N Q T
T T W W T X
X == （8-2-16）
二、附加条件法（伪观测值法）
前面已提及，秩亏自由网平差就是在满足最小二乘min =PV V T
和最小范数
min ˆˆ=x x
T 的条件下，求参数一组最佳估值的平差方法，实际上就是求相容方程组0ˆ1
1
=-⨯⨯⨯u u u u W x
N 的最小范数解。

附加条件法的基本思想：由于网中没有起算数据，平差时多
选了d 个未知参数，因此在u 个参数之间必定满足d 个附加条件式，即在原平差函数模型中需要加入d 个未知参数间的限制条件方程，从而可以按附有条件的间接平差法求解。

问题的
关键是如何导出等价于min ˆˆ=x x
T 的限制条件方程的具体形式。

为叙述方便，我们先给出该限制条件方程，然后再推导平差计算公式，最后证明，在给定的限制条件方程下所求得的解，就是相容方程组0ˆ1
1
=-⨯⨯⨯u u u u W x
N 的最小范数解。

设等价于约束条件min ˆˆ=x x
T 的限制条件方程为
0ˆ1
=⨯⨯u T
u
d x
S
（8-2-17）
式中，d S R =)(且满足，
0=BS S 称为附加阵。

故秩亏自由网平差的函数模型为
11
1
ˆ⨯⨯⨯⨯-=n u u n n l x
B V
权阵为P
0ˆ1
=⨯⨯u T
u
d x
S
按照附有条件的间接平差可得法方程
00ˆ0=⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝
⎛W K x S S N
s T
（8-2-18）
式中Pl
B W PB B N T u T u
u ==⨯⨯1
，，且u t B R PB B R N R T
<===)()()(，唯一不同的是这里
N 为秩亏阵。

为解决秩亏问题，将（8-2-18）中的第二式左乘S 矩阵后，再加到第一组中得：
00ˆ0=⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝
⎛W K x S S N
s T
（8-2-19）
式中T
SS N N +=，且u N R =)(
根据附有条件的间接平差原理，上式的解为
W N S S N S K T s 1
)(-= （8-2-20）
)(ˆ1s SK W N x -=- （8-2-21）
由于上述解是通过增加未知参数间满足的d 个附加条件，按照附有条件的间接平差法而实现的，因此人们把此法称为附加条件法。

但它又不同于经典的附有条件的间接平差法，其主要表现为：当S 阵满足0=BS 时，必定有下式成立（证明从略）
0=s K （8-2-22）
将（8-2-22）式代入（8-2-21）式，可得参数的解为
Pl B N W N x
T 11ˆ--== （8-2-23）
现在只需证明，按（8-2-23）式求得的解Pl B N W N x
T 11ˆ--==就是法方程0ˆ11
=-⨯⨯⨯u u u u W x
N 的最小范数解。

为此只需证明1
-N
是N 的最小范数g 逆中的一个即可，即只
需证明1
-N
满足以下两式：
N N N N N N N
N T 111
)(---==和
（8-2-24）
现证明如下：因为 T
SS N N += ，所以有
I SS N N N N T =+=--)(11 右乘S 阵并展开，则有
S S SS N NS N S N N T =+=---111
而0==PBS B NS T
，所以有
S S S S N T
=-)(1
（8-2-25）
由于
d S S R T =)(，存在逆阵，则有 1
1)(--=S S S S N T
（8-2-26）
所以有 T T T T S S S S I SS N I SS N N N N
1111
)()()(-----=-=-=
（8-2-27） N S S S NS N N N N T =-=--11)()(
（8-2-28）
因此（8-2-24）第一式得到验证。

由（8-2-27）式得
T T T T T T S S S S I S S S S I N N
111
)())(()(----=-=
考虑到（8-2-26）式，则上式为
N N N N N I SS N I N N
T T 1111
)()(----=--=-=
（8-2-29）
（8-2-28）、（8-2-29）两式说明1
-N 是N 的最小范数g 逆中的一个，因此按（8-2-23）式求
得的x
ˆ一定是相容方程组0ˆ1
1=-⨯⨯⨯u u u u W x
N 的最小范数解。

三、精度评定
单位权中误差估值的计算
r PV
V T ±=0ˆσ
（8-2-30）
式中PV V T
可以直接计算，也可以按下式求得
x
W Pl l PV V T
T T ˆ-= （8-2-31）
未知参数的协因数阵为
T
T T X X P B N Q P B N Q )(11ˆˆ--⋅⋅= 1
111----=⋅⋅=N N N N PB Q P B N T
)()(1111
-----=-=N SS I N N SS N N
T T
1
11----=N SS N N T
（8-2-32）
实际计算时，通常要对S 进行标准化，设标准化后的S 阵用G 表示，即不仅要求满足
0=BG ，还要求满足I G G T =，此时（8-2-26）式变成G G G G G N T ==--1
1)(，转置后有T
T G N G =-1，因此（8-2-32）式将变成如下形式
T
X X GG N Q -=-1ˆˆ
（8-2-33）
四、两点说明 ①若将
0=s K 代入法方程，则法方程变为
0ˆ)0ˆ=-+=-Pl B x SS PB B W x N T
T T 或（ 上式相当于下列误差方程联合组成的法方程
⎪⎭⎪
⎬⎫=-=x G V l x B V T g ˆˆ
上式的第一式为观测值L 的误差方程，第二式可以看作是为求最小范数解而人为增设的d 个虚拟误差方程，因此附加条件法又叫伪观测值法。

②该方法的特点就是用求凯利逆替代了求广义逆，因此便于计算和计算机编程，但首要条件是必须知道附加阵S ，关于附加阵的确定问题，本教材不准备作详细讨论，下面直接给出常见控制网的附加阵S 及其标准化后的矩阵G 的具体形式：
水准网（设有u 个点）
()1111
=⨯T
S
μ ；
⎪⎪⎭⎫
⎝⎛=⨯μμ
μμ11
11
T
G
（8-2-34）
测边网（设有m 个点）
⎪⎪⎪⎭⎫
⎝⎛---=⨯00020
2
1012310
1
10010101
m m
m
T
x y x y x y S
（8-2-35）
式中0
0i i y x 、为第I 点的近似坐标
⎪
⎪⎪⎪⎪⎪⎪⎭⎫
⎝⎛---∑=⨯00
020
2
10
12
23(1)10
1
10(1)010101(1
m m
m T x y x y x y R m m G
（8-2-36）
式中00i i y x 、是以中心坐标为原点的第I 点的近似坐标，它们的计算如下：
∑-=m i i
i
x m x x 1000
1，∑-=m i
i i y m y y 100
01
测角网（设有m 个点）
∑+=
∑m
i i y x R 1
2
0202
)
(
只需在（8-2-35）式中增加一行()
00
02020101m
m
y x y x y x
元素、在（8-2-36）
式中增加一行()
0020201012
1
m
m y x y x y x
R
∑元素即可得到相应的S 阵和G
阵。

例[8-3] 如图8-2水准网，C B A 、、点全为待定点，同精度独立高差观测值为
m h 345.121=，m h 478.32=，m h 817.153-=，平差时选取C B A 、、三个待定点的高
程平差值为未知参数321ˆ
ˆˆX X X 、、，并取近似值
)
823.25345.22100302010
m X X X X （⎪
⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛=
试分别用直接法和附加条件法求解参数的平差值及其协因数阵。

解：1．直接解法 误差方程为
⎪⎪⎪
⎭⎫ ⎝⎛-⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛---=600ˆˆˆ101110011321x x
x V 法方程为
0606ˆˆˆ211121112321=⎪⎪⎪
⎭⎫ ⎝⎛--⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛------x x
x
由法方程易知
⎪⎪⎭⎫ ⎝⎛--=211211N ， ⎪⎪⎭⎫ ⎝⎛----=1211121N ， ⎪
⎪⎭⎫ ⎝⎛=061W
所以有
⎪⎪⎭⎫
⎝⎛=
=-6336271)(11111T N N Q
未知参数的改正数为
图8-2
A
C
3
)
(202)(ˆ111111111mm W Q N W N N N x T
T T ⎪⎪⎪⎭⎫ ⎝⎛-===-
未知参数的平差值为
)
821.25345.22002.10ˆˆˆˆˆˆ321030201321m x x x X X X X X X （⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎭⎫ ⎝⎛
未知参数的协因数阵为
⎪⎪
⎪⎭⎫
⎝⎛------==2111211129111111111ˆ
ˆN Q N Q N Q T X X
2．附加条件法
解法一中已求得法方程为0ˆ=-W x
N 的具体形式为： 0606ˆˆˆ211121112321=⎪⎪⎪
⎭⎫ ⎝⎛--⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛------x x
x
该水准网有3个待定点，所以附加阵为
()
1111
3=⨯T
S
⎪⎪⎭⎫ ⎝⎛=⨯313
13
11
3T
G
则有
T
GG N N +=
⎪
⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝
⎛------=11111111131211121112
⎪⎪⎪⎭⎫
⎝⎛------=72227222731
⎪⎪⎪⎭⎫
⎝⎛=-522252225911
N
所以有
)(20260652225222591ˆ1
mm W N x
⎪⎪⎪⎭⎫
⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛-⎪⎪⎪⎭⎫ ⎝⎛==-
未知参数的的协因数阵为
⎪⎪
⎪⎭⎫
⎝⎛------=-=-211121112911ˆ
ˆT X X GG N Q
结果与直接解法完全相同。