在许多情况下，系数矩阵同样由观测量组成，而这些观测量是通过测量手段获取的，就不可避免地存在误差，即系数矩阵也含有误差。当系数矩阵和观测向量同时含有误差时，应建立EIV模型，采用总体最小二乘求解。近年来，已有大量关于总体最小二乘的研究，其理论日趋成熟，成为测绘数据处理领域研究的热点问题。文献[1-2]对加权总体最小二乘算法进行研究，文献[3-4]推导病态总体最小二乘算法，文献[5-6]研究附加约束条件的总体最小二乘，以及文献[7]研究多元总体最小二乘问题。以上研究都认为各观测量的精度是相同的，然而在测绘数据处理领域中，如GPS网、三角网和导线网，经常含有不同类型或不同精度的观测量。平差计算时，一般根据先验方差进行定权，但是先验方差并不能精确已知，此时确定的各类观测值的权不够准确，需要采用验后方差分量估计的方法进行二次定权。验后方差分量估计(VCE)方法主要包括赫尔默特估计(Helmert)、最小范数二次无偏估计(MINQUE)、最优不变二次无偏估计(BIQUE)及VCE的扩展研究等。针对总体最小二乘中含有不同精度的观测量，文献[8-9]将总体最小二乘问题和VCE进行综合研究，分别采用LS-VCE和MINQUE-VCE对总体最小二乘问题中观测向量和系数矩阵中元素的方差分量进行求解，从而得到更准确的参数估计值；文献[10-11]将Helmert方差分量估计用于总体最小二乘中，同样对观测向量和系数矩阵中元素的方差分量进行求解。

另外，观测量中通常不可避免地含有粗差，如何对粗差进行探测和定位，消除其对平差结果的影响，是提高平差结果精度的一个关键问题。基于最小二乘理论，周江文^[12]提出IGG权函数，杨元喜^[13]进行了扩充，提出IGGⅡ权函数。已有许多学者将稳健估计理论与总体最小二乘结合，陈义、陆珏^[14]以三维坐标转换为例，阐述一种稳健总体最小二乘方法，实际是将Huber权函数用于Gauss-Helmert加权总体最小二乘；龚循强、李志林^[15]进行改进，提出基于IGG权函数的稳健总体最小二乘，并与文献[14]提出的方法作详细的比较分析，得出IGG权函数的稳健总体最小二乘计算结果更为可靠；文献[16-17]引入标准化残差和中位数来提高稳健总体最小二乘的抗差性能。

总体最小二乘方差分量估计并不具有稳健性，当观测数据中含有粗差时，其计算结果往往面目全非，得不到准确的参数估值。因此，本文在进行总体最小二乘方差分量估计前，先根据稳健估计的原理选取合适的权函数，分别构建系数矩阵和观测向量的等价权，再以构建的等价权作为初始权，按照总体最小二乘Helmert方差分量估计进行计算。为了方便直接引用Helmert方差分量估计，本文从最小二乘配置的角度，根据广义最小二乘原理，将总体最小二乘转换为最小二乘的形式求解。

1 加权总体最小二乘算法

总体最小二乘函数模型如下^[1]：

(1)

式中，是观测值向量，是观测值向量的误差改正值，是系数矩阵，E_A是系数矩阵误差改正矩阵，是待估计参数。

采用高斯-牛顿法将非线性模型(1)转化为线性模型，迭代求解。设第j次迭代后的估值为ξ_(j)，则参数ξ表示为：

(2)

将式(2)代入式(1)中，为了避免迭代发散或收敛错误，采用E_A的第j次迭代值E_A(j)代替E_A计算系数矩阵，记A_(j)=A+E_A(j)，得：

(3)

式中，e_A=vec(E_A)，vec(·)表示将矩阵按列拉直变换。

考虑到系数矩阵中可能同时含有无误差的常数项和随机元素，并且不同的位置有相同的元素存在，为了保证常数项的改正数为0，不同位置的相同元素改正数相同，将系数矩阵中的随机元素提取出。假设系数矩阵中含有u个不同的随机元素，则可以组成u×1的列向量a，用e_a表示随机元素a的误差改正数，则有e_A=Me_a，M是由常数构成的a的系数矩阵，为mn×u矩阵，可通过求矩阵vec(A)对随机元素a的雅可比矩阵获得。

代入式(3)得：

(4)

可以看出，式(4)为最小二乘配置的标准形式，将δξ看作不考虑随机信息的倾向参数，e_a看作考虑随机信息的信号参数。根据广义最小二乘原理，将e_a的期望作为虚拟观测值，而e_a的期望为0，则可以得到：

(5)

将总体最小二乘转换为最小二乘求解，整理得：

(6)

此时总体最小二乘目标函数为V^TP_LV+V_ea^TP_aV_ea=min，P_L为观测向量L的权阵，P_a为随机元素a的权阵，令，根据最小二乘原理得：

(7)

通过迭代计算，直到满足收敛条件为止。

2 稳健总体最小二乘Helmert方差分量估计 2.1 总体最小二乘的Helmert方差分量估计

本文从最小二乘配置的角度，根据广义最小二乘原理将总体最小二乘转换为最小二乘形式，进而可以直接引用Helmert方差分量估计。

令，则总体最小二乘Helmert方差分量估计计算公式为^[11]:

(8)

(9)

(10)

式中，表示观测向量的方差分量估值，表示系数矩阵中随机元素的方差分量估值，，tr(·)表示求矩阵的迹。

通过对式(9)求解，可以得到观测向量和系数矩阵中随机元素的方差分量估值：

(11)

总体最小二乘多余观测数r与最小二乘多余观测数是一样的：r=(m+u)－(n+u)=m－n，但总体最小二乘多余观测数在观测向量和系数矩阵的随机元素之间进行分配，即r=r_L+r_a。这里可以直接得出观测向量和系数矩阵中随机元素的多余观测数为：

(12)

则观测向量和系数矩阵中随机元素的方差分量估值可简化为：

(13)

总结Helmert方差分量估计的计算步骤如下:

1) 按最小二乘求参数初值ξ₍₀₎=(A^TP_LA)^－1A^TP_LL，E_A0=0；

2) 进行第一次总体最小二乘平差计算，求得V和V_ea，并计算V^TP_LV和V_ea^TP_aV_ea；

3) 按式(13)进行方差分量估计，求得观测向量和系数矩阵中随机元素的方差分量估值、，再依下式重新定权：

(14)

式中，c为任意常数，一般选中的一个值；

4) 重复步骤2)、3)，直至为止。

2.2 稳健估计

为了使总体最小二乘方差分量估计具有稳健性，本文继续引入稳健估计理论，在总体最小二乘方差分量估计前先进行稳健估计。已有大量文献证明，IGG权函数更适于测量计算。IGG权函数为^[13]：

(15)

(16)

式中，P_{L, i}为观测向量的等价权P_L的第i个对角线元素，P_{a, i}为随机元素a的等价权P_a的第i个对角线元素；P_{L, i}为P_L的第i个对角线元素，P_{a, i}为P_a的第i个对角线元素；V_i为V的第i个元素，V_{ea, i}为V_ea的第i个元素；k₀、k₁为常数(一般取k₀=1.5, k₁=2.5)。

本文提出的稳健总体最小二乘Helmert方差分量估计分为两步来计算：第一步是按照总体最小二乘计算残差V和V_ea，按式(13)计算观测向量和系数矩阵中随机元素的单位权中误差估值，然后由式(14)构建等价权，通过迭代计算，获得最终的等价权P_L和P_a；第二步以P_L和P_a作为观测向量和系数矩阵中随机元素的初始权，代替P_L和P_a，按总体最小二乘Helmert方差分量估计的计算步骤求解参数估值。

3 算例分析

以六参数平面坐标转换为例，采用文献[14]中的数据，并在u_s1上加入大小为1的粗差，在v_t5上加入大小为1.5的粗差。数据如表 1所示，(u_s, v_s)表示原坐标系的坐标，(u_t, v_t)表示目标坐标系的坐标，参数真值。为了进行比较分析，采用不同的方法进行计算，计算方法列于表 2，并将不同方法计算得到的坐标转换参数估值和估值与真值的差值范数列于表 3。图 1为方案4的迭代计算过程中两个方差分量的收敛图，图 2为方案6的迭代计算过程中两个方差分量的收敛图。

1) 从表 3的计算结果可以看出，方案2与方案1的计算结果基本相同，方案4与方案3的计算结果基本相同，说明本文的加权总体最小二乘算法和总体最小二乘Helmert方差分量估计均是切实可行的。加权总体最小二乘认为观测向量和系数矩阵所对应的单位权方差是相同的，然而实际上表 1中原坐标系坐标观测值与目标坐标系坐标观测值的精度并不相同，导致观测向量和系数矩阵分别对应不同的单位权方差，此时采用加权总体最小二乘算法显然是不合理的，需要通过验后估计的方法进行二次定权。本文的总体最小二乘Helmert方差分量估计将原坐标系坐标观测值与目标坐标系坐标观测值看作具有不同的单位权方差分量，通过方差分量改变权值，最终达到方差分量一致，因此参数的计算结果得到了进一步的改正，与真值最为接近。

2) 从方案5的计算结果可以看出，其参数估值与真值差异较大，这是因为总体最小二乘Helmert方差分量估计不具有抵抗粗差的能力。方案6的计算结果与真值接近，因为方案6通过先构建等价权，再进行方差分量估计，使得算法具有了稳健性，说明本文的稳健总体最小二乘Helmert方差分量估计是可行的，在观测数据中含有粗差的情况下，可以很好地消除粗差的影响，提高方差分量估计计算结果的可靠性。

4 结语

本文从最小二乘配置的角度，根据广义最小二乘原理将总体最小二乘转换为最小二乘求解，从而可以引用Helmert方差分量估计的原理，处理不同精度或者不同类型的观测量数据，得到更可靠的计算结果。同时还考虑到观测数据中可能含有粗差的情况，提出稳健总体最小二乘Helmert方差分量估计，使得总体最小二乘方差分量估计具有稳健性。