近年来，加权总体最小二乘平差是测绘数据处理领域研究的热点问题。目前,处理加权总体最小二乘问题的算法有Schaffrin-Wieser算法^[1]、Mahboub算法^[2]和Amiri-Simkooei算法^[3]等。通过对加权总体最小二乘已有算法的进一步分析发现，在平差前随机模型要已知验前方差才能精确地给观测数据定权。在过去很长时间，平差都在单一的同类观测量中进行，例如测角网平差、水准网平差；定权都是采用测量平差中常用方法定权，例如水准高差按路线长度倒数定权等。随着平差对象从单一同类观测扩展为不同类的多种观测量，一般它们的验前方差都不能精确已知，而且实践证明，原先采用的估计验前方差确定各类观测量权的方法，在许多情况下不够准确。为了提高方差估计的精度，20世纪70年代开始出现了用验后信息估计各类观测量方差，然后定权的方法。方差分量估计(VCE)主要有赫尔默特估计(Helmert)^[4]、最小范数二次无偏估计(MINQUE)、最优不变二次无偏估计(BIQUE)及VCE的扩展研究^{[5, 6]}，这些都是在最小二乘准则的基础上进行的，没有顾及系数矩阵的误差。对于以上存在的问题，文献[7, 8]将总体最小二乘问题和VCE进行了综合研究，分别采用LS-VCE和MINQUE-VCE对总体最小二乘问题中观测向量和系数矩阵元素的方差分量进行求解，从而得到更准确的参数估计值。针对EIV模型中存在随机模型不明确，参数估计不准确的问题，本文提出了一种加权总体最小二乘平差随机模型的验后估计方法，既考虑到观测向量 y和系数矩阵A 中元素可能具有不同精度的情况，又顾及到在给观测向量 y和系数矩阵A 中元素定权时所选取的初始单位权方差可能不同导致定权不准确的问题。在算例部分，第1个算例采用文献[9]中的实验数据，将本文算法与文献[3]、[7]中的算法进行了比较，验证本文算法的可行性；第2个直线拟合和第3个线性参数估计算例，主要对本文的验后估计算法和本文提出的加权总体最小二乘快速算法，以及最小二乘算法进行了比较，实验表明本文提出的算法在解决随机模型不准确的总体最小二乘问题时更合理、有效。

当同时考虑系数矩阵和观测向量误差时，有EIV模型为：

式中,y ∈ R ^n×1为观测向量； e _y 为观测向量的随机误差； A ∈ R ^n×m为列满秩的系数矩阵； E _A 为系数矩阵的随机误差； X ∈ R ^m×1为待估参数向量；其随机模型分别为：

式中,e _A =vec( E _A )，即 e _A 是将 E _A 按列向量化后得到的nm×1列向量； P _y 和 P _A 分别为 e _y 和 E _A 的列向量化，即vec( E _A )的权阵；σ₀₁²为观测向量的单位权方差；σ₀₂²为系数矩阵元素的单位权方差。

假设 e _y 、 E _A 的估值分别为 ，则加权总体最小二乘平差准则为：

一般的加权总体最小二乘中观测值和系数阵中的元素是相互独立的，且系数矩阵中的元素通常是由其他观测值经过计算得到的，也就是观测值和系数阵中元素的随机模型是不相同的；在总体最小二乘平差之前给定的观测值的权阵 P _y 和系数矩阵元素的权阵 P _A 是不恰当的，或者说它们所对应的单位权方差σ₀₁²和σ₀₂²不相等，因此对这个问题必须利用平差后的结果进行处理，这就是加权总体最小二乘平差随机模型的验后估计。

设参数 X 的估值为， 其中 X ⁰表示参数的近似值，为未知参数近似值的改正数，观测值的估值，总体最小二乘平差除了考虑观测值 y 的误差外，同时还考虑系数矩阵 A 的误差 。根据式(1)有：

令 ，则式(5)变为：

所以,加权总体最小二乘准则式(3)变为：

结合式(6)、式(7)得加权总体最小二乘的估值：

式中， 。

下面给出 Q_VV 的求解方法。由 ，得^{[10, 11]}：

式中， ，其中 I _n表示n阶单位阵，是Kronecker积算子;

。其中 表示列向量的第i个元素，表示矩阵 的第(i,j)个元素。

若已知列向量中的所有元素 、之间的互协因数及协因数组成的矩阵，则由式(9)根据协因数传播律得：

单位权方差的估计为：

参数估值的协因数阵为：

式中，Q _y 表示观测向量 y 的协因数阵。

将式(6)和式(8)代入式(9)得：

由最小二乘条件平差得：

由此可以得到观测值及系数阵随机误差的估值和 ,对应的单位权方差估值分别为：

式中，r _y 为观测向量 y 对应的自由度；r _A 为系数矩阵 A 对应的自由度。

加权总体最小二乘快速算法(简称fast-WTLS)迭代过程如下。

1) 给定已知或者观测数据 A 、 y 、 P _A 、 P _y ；

2) 设置估计值的初值为：

3) 利用式(8)、(10)计算 Q_VV 、，从而得到；

4) 如果与上轮迭代的差值范数小于设定的阈值，则终止迭代得到最终解；否则更新参数近似值，转到步骤3)，重复迭代计算。

将式(4)进行调整，结合公式

得：

式中，表示fast-WTLS算法迭代过程中第i次迭代计算得到的系数矩阵和未知参数的值；Δ A _i、Δ X _i为迭代过程中相应的改正数。 舍去二次项Δ A _iΔ X _i，式(17)、式(18)可以改写为^[8]：

令 ，，。

式(19)可以表示为：

且有下列关系式：

式中，权阵。

此时便可以借鉴Welsch推证的赫尔默特方差估计严密公式^[4]得到EIV模型的方差分量估计公式(Helmert-VCE)。

在EIV模型中，赫尔默特方差分量估计公式可以表示为：

式中，

表示迹算子。

通过对式(21)求解，可以得到观测向量和系数矩阵的方差分量：

本文提出的方差分量求解方法是将函数模型式(19)表示成间接平差的形式，通过间接平差模型下的Helmert-VCE公式推导得到EIV模型中方差分量估计公式，因此EIV模型的精度评定和多余观测数的求解工作完全可以在此基础上按间接平差的原理继续进行。EIV模型的多余观测数估计公式为：

EIV模型的观测数据分为两个部分，分别为n个观测向量元素和t个系数矩阵中观测元素(随机变量)，共n+t个。求解的参数也分为两个部分，分别是m个模型参数和t个系数矩阵的待估计量(随机变量)，共m+t个。因此，EIV模型和高斯-马尔科夫模型相比，自由度没有变化^{[8, 12]}。

同样可以得到观测向量和系数矩阵的多余观测分量：

加权总体最小二乘平差随机模型验后估计的迭代计算步骤为:

1) 给定观测值和系数矩阵权阵的初值 P _y1和 P A₁，求得 Q _y1和 Q _A1，然后组成；

2) 进行第一次加权总体最小二乘平差，由式(14)求得，从而得到观测值及系数阵第一次平差后的和；

3) 根据式(15)、式(16)、式(24)、式(25)分别得到观测值及系数阵第一次平差后的单位权方差估值和；

4) 依式(26)重新定权：

式中，c为任一常数，一般是选或中的某个值，i≥1为迭代次数。求得( Q _y )_i和( Q _A )_i，然后组成 ；

5) 进行步骤2)~4)，即平差—单位权方差估计—定权后再平差，直到 为止。

为了检验本文提出的加权总体最小二乘平差随机模型验后估计方法的效果，采用文献[9]中的直线拟合数据，利用本文提出的算法进行参数求解，并与文献[3]中的加权总体最小二乘算法(WTLS)，文献[7]中提出的最小二乘方差分量估计(LS-VCE)结合加权总体最小二乘算法，以及文献[3]中的WTLS算法结合本文提出的Helmert-VCE方法所求得的参数结果进行比较，以分析这四种算法的差异。结果见表 1,上述方法所求得的单位权方差估值计算结果见表 2。图 1表示迭代计算过程中两个方差分量的收敛图。

图 1 算例2.1中线性回归模型迭代过程中 Fig. 1 Estimates of Two Variance Factors and of Linear Regression Model over Different Iteration Steps

图选项

假设直线方程为：

直线拟合的模型为式(27)，系数矩阵所对应的协因数阵 Q _A 为式(28)：

1) 在加权总体最小二乘算法^[3]中，把观测向量和系数阵所对应的单位权方差当作相同，得到的参数结果与文献[9]中提供的“参数真值”基本一致。在采用方差分量估计结合WTLS中，考虑了观测向量和系数阵的验前单位权方差一般不能准确已知或二者是不同的，例如在直线拟合模型中，就考虑了x和y具有不同的单位权方差分量 和；因此，参数的估计结果就得到了进一步改正。

2) 从表 1可以看出，采用方差分量估计结合加权总体最小二乘方法(包括本文方法)求解的参数结果与文献[9]提供的“参数真值”有差别，这是合理的，因为文献[9]中的数据和权值是由不同人给出的，无法确定权值的准确性，采用WTLS算法进行参数求解，显然不合理。文献[7]中也对文献[9]提供的“参数真值”的准确性提出了质疑，认为在确定观测向量和系数阵权阵时采用了不同的单位权方差，因此权阵并不合理。

3) 本文提出的加权总体最小二乘平差随机模型验后估计方法实质是fast-WTLS和赫尔默特方差分量估计的结合。其中fast-WTLS算法是一种近似的WTLS算法^[13]，从算法结构可以看出，该方法的求解更接近最小二乘算法，在进行数据处理时更快速、便捷。从表 1中可以看出，本文算法和WTLS结合LS-VCE算法求解得到的参数结果有差异，而WTLS结合Helmert-VCE算法与WTLS结合LS-VCE算法求解得到的参数结果基本相同，从而可以判断本文算法与WTLS结合LS-VCE算法的差异存在于求解参数的加权总体最小二乘算法中，而方差分量估计部分Helmert-VCE与LS-VCE等价。

4) 在EIV模型中，LS-VCE首先设定方差分量初值，一般设为1，经过迭代得到不同的方差分量估值，而迭代过程中相应的协因数阵不变；而Helmert-VCE估计是通过方差分量修改权阵，达到最终方差分量一致，其中权阵变了。两种方法最终都能达到统一单位权方差的目的。从表 2中单位权方差估值的计算结果也可以看出这一点。

5) 对于非线性模型，fast-WTLS算法在推导过程中出现了非线性化的情况，通过对二次项的舍弃达到了模型的线性化，因此对于非线性问题也可以较好的处理。

利用Matlab软件编程进行数值模拟实验，模拟普通直线 y=β₁+β₂x,其不含误差的x、y值见表 3。待估的未知参数分别为直线的斜率β₂和截距β₁，真值为β₂=5、β₁=8。在观测值x和y上分 别加入均值为0，方差为σ_x²=σ₀₁²· P^－1_x和σ² _y =σ₀₂²· P_y ^－1 的随机误差，其中σ₀₁²=1，σ₀₂²=0.25；采用本文提出的加权总体最小二乘平差随机模型的验后估计进行求解，并与最小二 乘(LS) 、fast-WTLS所求解的参数结果进行比较。参数估值结果见表 4。图 2表示普通直线拟合算例中两个方差分量的收敛图。

表 4 模型参数计算结果 Tab. 4 The Calculation Results of Model Parameters

	真值	LS	fast-WTLS	fast-WTLS结合 Helmert-VCE
	5	4.826 204 402	4.841 768 064	4.907 255 311
	8	8.734 877 616	8.725 097 232	8.546 508 439
	—	0.755 149 005	0.742 161 264	0.554 322 155

表选项

图 2 算例2.2中线性回归模型迭代计算 过程中两个方差分量估计值 Fig. 2 Estimates of Two Variance Factorsand of Linear Regression Model over Different Iteration Steps

图选项

模拟一个线性参数估计 问题，有 Ax = b ，其不含误差的设计矩阵 A ∈ R ^n×m和观测值 b ∈ R ^n×1见表 5。待估的未知参数有5个，x =[x₁ x₂ x₃ x₄ x₅]^T ，真值=[1 1 1 1 1]^T。观测值 b的噪声e ~N( 0 ,σ₀₁²P_b ^－1 )，其中σ₀₁²=0.1，P _b =diag(1,2,…,n)。设计矩阵元素的误差为 e _A ~N(0,σ₀₂² P_A ^－1 )，其中σ₀₂²=0.5，P _A =diag(n×m,n×m－1,…,2,1)；采用本文算法求解参数，并与 LS、fast-WTLS所求解的参数进行比较，参数估值结果见表 6。图 3表示线性参数估计算例中两个方差分量的收敛图。

表 6 模型参数计算结果 Tab. 6 The Calculation Results of Model Parameters.

	1	2	3	4	5	‖Δ‖
真值	1	1	1	1	1	/
LS	0.946 093 655	1.093 801 465	0.970 298 037	0.986 810 384	1.066 541 909	0.131 105 328
fast-WTLS	0.960 154 916	1.071 157 298	0.994 505 507	0.997 673 115	1.050 579 762	0.096 150 444
fast-WTLS结合 Helmert-VCE	0.966 018 915	1.063 779 361	0.995 815 997	1.000 280 826	1.047 417 512	0.086 536 271

表选项

图 3 线性参数估计迭代计算过程中两 个方差分量估计值 Fig. 3 Estimates of Two Variance Factors and of Linear Parameter Estimation Model over Different Iteration Steps

图选项

从表 4和表 6中可以看出，在普通直线拟合和线性参数求解中使用本文算法进行求解，拟合的结果比fast-WTLS、LS方法拟合的结果都要好。在本文提出的方法中，不仅考虑到系数矩阵中某些元素是常数，不存在误差，而且考虑到各元素之间是不等精度的，以及各元素之间存在相互联系，同时又考虑到不同类型观测值的验前方差一般都很难精确估计，这样就无法给观测向量和系数阵进行精确定权，采用随机模型验后估计方法，能够解决这些问题。所以采用本文方法解算的结果更加接近真值，也更为合理。

在线性模型参数估计中随机模型的验后估计应用很广泛，本文提出的加权总体最小二乘随机模型的验后估计分为两部分内容，第一部分是针对函数模型参数的求解；第二部分针对观测向量 y 和系数矩阵 A 的方差分量的求解，因此引入了随机模型的验后估计方法，对观测向量和系数矩阵进行权的重新分配，直至它们所对应的方差分量相等。本文算法推导过程依托了间接平差中赫尔墨特方差分量估计的原理，其优点为可以利用其理论体系中已有的公式推导适用于EIV模型的精度评定和求解自由度的公式，但同时也会继承其缺点。由于算法结构的特点，求解方差分量的迭代过程会出现不收敛，或虽然收敛但是出现负方差分量的情况；本文所模拟的实验仅考虑观测向量 y和系数矩阵A 相互独立的情况，并且假设协因数阵均为对角阵，相关情况未进行数值模拟实验，这两点无疑对本文方法的适用范围进行了限制，今后会针对这两方面做进一步工作。