总体最小二乘^[1]是顾及了系数矩阵误差的平差方法，是Errors-in-Variables(EIV)^[2-3]模型的严密估计方法。在EIV模型的基础上，文献[4]将系数矩阵中不同的随机元素提取并作为参数进行求解，提出了Partial Errors-in-Variables(Partial EIV)模型。文献[5]将Partial EIV模型进行转换，进行两步间接平差法得到参数估值，在形式上比文献[4]简单。文献[6]分析了Partial EIV模型算法与一般算法的优势分别为：① Partial EIV模型是EIV模型的更一般的表达式；② Partial EIV模型是将系数矩阵中不同的随机元素提取并作为参数进行求解，因此Partial EIV模型系数矩阵中待改正量的个数要小于或等于对应EIV模型中待改正量的个数，提高了计算效率；③ 便利了后续的精度评定。总体最小二乘法在近年来发展迅速^[7-13]，在平差时随机模型的不准确对参数估值有很大的影响，方差分量估计可以对随机模型进行修正从而得到更加合理的参数估值。方差分量估计(VCE)方法主要有赫尔默特估计^[14-15](Helmert)、最小范数二次无偏估计^[16](MINQUE)、最优不变二次无偏估计^[17-18](BIQUE)及最小二乘方差分量估计^[19-22](LS-VCE)；最小二乘方差分量估计(LS-VCE)方法是Teunissen提出，该方法使用的是最小二乘准则，经过转换将方差分量作为参数进行解算，得到的方差分量估值具有无偏性，且便利了后续工作。文献[20]给出了不同情形下的最小二乘方差分量估计方法并探索其特性。文献[21]将最小二乘方差分量估计应用于EIV模型中，将加权总体最小二乘参数解的表达式写成标准化的最小二乘形式，结合最小二乘方差分量估计方法确定方差分量估值，然而算例部分函数模型的系数矩阵既有随机元素又有非随机元素，使用更具有一般性的Partial EIV模型解算更具有代表性。与其他方差分量估计方法相似，在计算中可能出现负方差，函数模型多余观测量的不足和随机模型结构的不正当是出现负方差的主要原因，在最小二乘方差分量估计方法的基础上，文献[22]结合非负最小二乘理论，将非负最小二乘方差分量估计应用到GPS时间序列中，方差分量估计方法在本质上是相同的。文献[15]指出由MINQUE可以导出严密的Helmert估计公式。文献[23]针对概括函数平差模型进行公式推导，总结了不同方差分量估计方法并指出最小二乘方差分量估计方法的一般性，通过公式推导得到了MINQUE、Helmert和BIQUE方法均是最小二乘方差分量估计的特例。相比于其他方法的方差分量计算公式推导，LS-VCE方法过程更加简单、易于理解、具有较强的应用价值。

本文在文献[5]参数估计方法的基础上，结合最小二乘方差分量估计方法计算系数矩阵与观测向量的方差分量估值，在迭代过程中更新协因数阵，从而达到修正参数估值的效果。针对估计中出现的负方差，增加非负约束条件，使用非负最小二乘方差分量估计方法计算方差分量估值。通过算例试验对本文方法进行验证，并与已有方差分量估计方法进行比较。与文献[14]相比，本文方法不需要引入权比因子，在计算参数估值时迭代次数要更少，且方便了后续方差分量估值的精度评定，当计算中出现负方差时，文献[14]的方法则会出现不可估的现象。与文献[21]相比，本文方法继承了原有Partial EIV模型的优势。

1 Partial EIV模型算法

在某些实际应用中，EIV^[24-25]模型中系数矩阵由一些非随机元素与随机元素组成，文献[4]对系数矩阵进行处理，提取了系数矩阵中的随机元素，将EIV模型改写成Partial EIV模型^[4]

函数模型

随机模型

式中，表示n×m系数矩阵真值；h是一个确定的nm×1常数向量；由中非随机元素和零组成；B是一个给定的nm×t矩阵；t表示系数矩阵中随机元素的个数；y表示n×1的观测向量；e_y表示观测向量误差；a是系数矩阵中随机元素组成的t×1列向量；表示其真值；e_a表示a的随机误差；Q_y表示观测量误差的协因数阵；Q_a表示系数矩阵误差的协因数阵；σ₀²表示单位权方差；I_n表示n×n单位矩阵。

文献[5]对式(1) 进行变形得到

式中，令是n+t列向量；；。

将观测值a作为参数进行平差，计算得到^[5]

用参数估值重新构造系数矩阵估值，并计算系数矩阵改正数，A表示由观测向量组成的矩阵，计算协因数阵Q_C^{[5, 26]}

经计算得到^{[5, 26]}

为将法矩阵变成对称正定矩阵，对式(6) 加上和减去并进行解算得到

式中，根据协方差传播律可得到。进而可得到标准化最小二乘形式的参数估值表达式^{[5, 26]}

进而得到观测量估值及残差^[21]

式中，表示投影阵。

2 Partial EIV模型的非负最小二乘方差分量估计 2.1 最小二乘方差分量估计

最小二乘方差分量估计^[19-23]由Teunissen提出，使用的是最小二乘准则，解算出的方差分量估值具有无偏性；文献[14]在Partial EIV模型中引入权比因子，使用赫尔默特方差分量估计方法确定权比因子，并通过算例验证了与其他方差分量估计方法的等价性。文献[14]在参数计算过程中使用的是文献[4]的算法，在表达形式上相对较复杂，且文献[5]指出该算法收敛较慢，影响计算效率；方差分量估计的本质是相同的，都是对随机模型进行修正，即修正观测值的权，从而达到修正参数估值的效果，最小二乘方差分量估计是更为一般的方差分量估计方法，而且最小二乘方差分量估计有利于后续对方差分量估值的精度评定及其他方面的应用，如出现负方差时可以增加约束条件使其成为附有约束条件的最小二乘平差，处理出现负方差的情况，而出现负方差时，其他方差分量估计方法会出现不可估的情况。针对总体最小二乘中系数矩阵与观测量有不同方差分量的情况，使用LS-VCE对随机模型进行修正更加合理。

由式(1)、式(2) 及式(8) 可知，期望可以表示为，引入矩阵R和向量t，矩阵R满足，t是一个均值为0的n-m维向量，参数和观测向量y₁与t、R的关系式^{[20, 22]}为

式中，表示参数估值的协因数阵。

由式(10) 可知，根据协因数传播律可得到向量t的协因数阵，因此参数及向量t的期望与协因数阵又可以表示为

式中，协因数阵Q_C如式(5) 所示，经过相应的变换可以得到方差分量的关系^[19-22]为

式中，A_vh和y_vh均表示对矩阵进行vh^[21]操作后的矩阵和向量，σ表示方差分量组成的向量。根据文献[20]可得

式中，σ为方差分量组成的向量；Q_vh^-1为σ的协因数阵；为p×p的矩阵；L=A_vh^TQ_vh^-1y_vh为p×1的列向量；p表示方差分量的个数，本文中p取值为2。根据协方差传播律可得到方差分量估值的方差阵为Q_σ=N^-1，经过进一步的化简，将其用Partial EIV模型中的协因数阵、残差等表示，其形式^[19-23]为

式中，n_kj、l_k分别表示矩阵N, L中的某一元素；Q_k、Q_l都表示为系数矩阵或观测量的协因数阵；Q₀为已知的矩阵，一般为零矩阵。

在Partial EIV模型中，观测向量和系数矩阵在某些情况具有不同的方差分量，对于式(2) 的随机模型存在

式中，σ₁、σ₂分别为观测向量和系数矩阵误差的方差分量；Q_y1、Q_a1为给定的观测量和系数矩阵的协因数阵。

将式(5) 改写成线性求和形式，即

Partial EIV模型是顾及了系数矩阵误差的总体最小二乘模型，文献[20]从经典的Gauss-Markov模型出发分析了最小二乘方差分量估计方法，相比于Partial EIV模型的数学模型，Gauss-Markov的数学模型可以表示为

对式(17) 进行平差解算可得到参数的表达式

误差向量e_y可以表示成多个分量求和的形式。在Partial EIV模型的解算中，式(8) 和式(18) 与Gauss-Markov模型解算得到的参数表达式及协因数阵的形式相同，这与文献[21]中EIV模型的最小二乘方差分量估计是一致的。

迭代计算分为参数估计和方差分量估计，具体的迭代步骤为：数据准备：A、y、Q_y、Q_a、h、a、B、收敛条件ε=10^-10。

(1) 用最小二乘法求解模型参数初值，给定初值a⁽⁰⁾，构造矩阵l、D；

(2) 给定方差分量初值σ⁽⁰⁾=[σ₁²⁽⁰⁾ σ₂²⁽⁰⁾]^T；设计迭代次数i=0；

外循环：

(3) 用式(15) 更新Q_y和Q_a；设计迭代次数j=0；

内循环：

(1) 更新Q_e，由式(4) 计算；

(2) 由式(5) 更新Q_C；

(3) 由重新构造矩阵，计算，更新；

(4) 由式(8) 计算参数，更新矩阵l，D，j=j+1；

(5) 重复步骤① 至步骤④ 直至满足时迭代终止。

(4) 更新；

(5) 由式(16) 更新Q_C；

(6) 计算和；

(7) 由式(14) 更新N和L，计算σ⁽ⁱ⁺¹⁾=N^-1L，i=i+1；

(8) 重复步骤(3) 到步骤(7) 直至时迭代终止。

2.2 非负最小二乘方差分量估计

最小二乘方差分量估计方法是根据最小二乘准则使得其误差平方和取最小，与其他方差分量估计方法类似，在进行方差分量解算时也会出现负方差，这与方差分量的自身含义是相互冲突的，出现负方差的原因有两种^[22]：① 函数模型观测量的不足，即多余观测量的不足。多余观测量也叫自由度，增加模型的多余观测量可以提高数据解算精度。② 结构不正当的随机模型，随机模型又表现在观测值的权中，然而观测数据初始的权往往是不精确的，严重影响了平差结果的精度。非负最小二乘方法^[27-28]是在最小二乘准则下增加非负约束条件而得，在最小二乘方差分量估计的基础上增加约束条件使得方差分量非负，因此可以将非负最小二乘方差分量估计作为附有约束条件的最小二乘问题，其准则^{[22, 28]}为

令，F(σ)是一个凸函数。式(19) 就等价于求解凸函数最小值问题，在Karush-Kuhn-Tucker(KKT)^[28]条件下构造拉格朗日函数得到

式中，u是拉格朗日乘子，对式(20) 求偏导计算得到

采用文献[28]的非负最小二乘方法，令σ_k是向量中的第k个元素，去除k之外剩余的元素用K={1, 2, …, p}\{k}表示；用式(13) 计算出方差分量，并对里面的元素进行约束，即

式中，表示中的第k个元素; 表示中的第k个元素; u_kⁱ表示中的第k个元素。

对中的某一元素进行约束，剩下的元素保持不变，因此式(21) 改写成^[27]

式中，n(:, k)表示矩阵N的第k列。

由式(23) 可以看出，迭代计算出的方差分量估值，若某一分量出现负值，通过式(23) 进行约束。非负最小二乘方差分量估计可以当成附有约束条件的最小二乘方差分量估计，其约束条件可以表达成

式中，C₁^T是z×p的矩阵；z表示σ中0元素的个数。

在计算参数协因数阵时，根据附有约束条件的平差方法，由协因数传播律可得到方差分量估值的协因数阵为

式中，。

具体迭代步骤为：

数据准备同最小二乘方差分量估计，在计算出参数估值后，在最小二乘方差分量估计的基础上使用非负最小二乘原理。

(1) 给定方差分量初值，设计迭代次数i=0；

外循环：

(2) 计算和；

(3) 用式(14) 更新N和L；令σ⁽⁰⁾=0，u⁽⁰⁾=-L；设计迭代次数j=0；

内循环：

(1) 用式(22) 计算σ^(j+1)；式(23) 计算u^(j+1)；

(2) 更新u^(j)=u^(j+1)，j=j+1；

(3) 重复步骤① 到② 直至时迭代终止。

(4) 更新；

(5) 重复步骤(2) 到(4) 直至ε时迭代终止；

(6) 计算矩阵C₁^T和P_1/C1；

(7) 由式(25) 计算方差分量估值的协因数阵。

3 算例与分析 3.1 算例1

数据采用文献[5, 21]直线拟合数据，已知坐标观测值(x_i, y_i)和相应的权值(p_xi, p_yi)，见表 1。

直线拟合的模型^{[5, 21]}为

为待估参数，直线拟合模型是常见的含有非随机元素与随机元素的模型，在使用Partial EIV模型进行解算时还需要给出向量h和固定矩阵B，由式(26) 可知，向量h和矩阵B的形式为

权阵p_a=，p_y=，收敛条件为ε=10^-10，分别用文献[5]的方法(WTLS)、文献[21]方法(方法1)、文献[14]方法(方法2) 以及本文方法计算参数估值及方差分量估值，得到的结果见表 2，其中表示单位权方差分量估值，表示系数矩阵与观测量的方差分量估值，、表示方差分量估值的方差，可以发现得到的方差分量估值相等，且相对应方差分量估值的方差也相等；方法1与本文方法的方差分量迭代过程变化相同(如图 1所示)，两种方法迭代次数都为21次，且文献[14]方差分量估计方面迭代次数也为21(迭代中收敛条件都取为ε=10^-10)，而参数估计方面文献[14]需要267次，本文只需要94次。

文献[14]与文献[15]推导了方差分量估值的方差公式为

式中, ，表示方差分量估值组成的向量，文献[14]再根据确定的权比因子计算方差；而在最小二乘方差分量估计中，由式(13) 根据协方差传播律可得，更简单且更易获取，其形式为

方差分量估计是针对随机模型不准确进而修正随机模型并再次平差的方法，本文以Partial EIV模型为基础的的最小二乘方差分量估计方法重新确定了两类观测值的权阵，得到的权值见表 3，可以发现与文献[14]和文献[21]重新确定的权值是相等的。

3.2 算例2

模拟一个二维仿射变换，对文献[29]的数据进行改化，已知6个转换参数真值a₁=0.9，b₁=-0.8，c₁=1，a₂=0.6，b₂=0.7，c₂=0.5和原始坐标、目标坐标的真值，相应的协因数阵为Q_a=0.005·diag([1 1 2 2 3 3 1 1 5 5 4 4 2 2 7 7 1 1 8 8 3 3 6 6]^T)，Q_y=0.005·diag([1 1 3 3 6 6 1 1 1 1 8 8 4 4 3 3 6 6 5 5 4 4 5 5 2 2]^T)，用Matlab软件mvnrnd函数给真值加上均值为0，协方差分别为Q_y和3Q_a的随机误差，得到随机一组坐标值见表 4。

二维仿射变换模型^[29]为

表 4中提供了26组坐标，向量h和确定矩阵B、y、a的形式分别为

式中，。

数据准备：h、B、y、a、p_a、p_y，权阵为协因数阵的逆，收敛条件为ε=10^-10；计算得到的方差分量估值见表 5，因为该算例出现负方差，导致文献[14]出现不可估现象，且本文与文献[21]方法等价，因此表 5中只列出了最小二乘方差分量估计(LSVCE)及非负最小二乘方差分量估计(NNLSVCE)进行计算，相应的方差分量估值变化图见图 2。其中表示目标坐标系(观测量)的方差分量，表示原始坐标系(系数矩阵)的方差分量，表示的方差，表示的方差。

图中，LS-VCE的迭代次数都为37次，NNLSVCE的迭代次数为14次，而在使用文献[14]方法进行求解时，因为方差分量出现负值的现象，计算时则出现不收敛的情况。

3.3 算例分析

(1) 由表 2可以看出本文方法与文献[14、21]得到的待估参数结果及方差分量估值相等，对应的方差分量估值方差也相等，图 1显示了方差分量估值的变化图，两种方法的迭代次数都为21次，表 3显示了修正方差分量所对应的权值。在不考虑方差分量估计情况下，文献[5]的参数估计需要迭代7次，文献[30]需要迭代8次；文献[14]计算中进行参数估计时需要迭代267次，而本文只需要94次，在一定程度提高了计算效率，且在计算方差时本文只需要根据协方差传播律获得，如式(25)，相对更简单。

(2) 用最小二乘方法计算算例2方差分量估值时出现负值，而文献[14]的方法则出现不可估的情况，使用非负最小二乘理论得到的方差分量估值变化图见图 2。在非负约束下，通过解算其余方差分量估值从而达到准则下的整体最优解。非负最小二乘方差分量估计可以看成为附有限制条件的间接平差，在结合LS-VCE后，非负方差分量估计实质上就为非负最小二乘估计，因此本文方法与文献[22]方法是等价的。

(3) Partial EIV模型是将系数矩阵中的重复的随机元素提取，使得模型需要改正的待改正量个数相对更少。对于Partial EIV模型的随机模型，其协因数阵构造简单，更具有一般性的Partial EIV模型在解算海量数据时更能体现其优势。在式(5) 计算协因数阵时，BQ_aB^T与EIV模型的系数矩阵的协因数阵Q_A相等，与文献[30]求解公式一致，相应的在进行方差分量估计时是等价的。

4 结论

本文以Partial EIV模型为基础，当观测数据的随机模型不准确时，结合最小二乘方差分量估计方法对随机模型进行修正，推导了以Parial EIV模型为基础的最小二乘方差分量估计公式，将观测向量误差和系数矩阵误差分别作为一类数据，从而计算出两个方差分量估值；并且当方差分量估值出现负的时候，使用非负最小二乘理论，增加非负约束条件，可以处理方差分量出现负值的现象。计算过程分为参数估计与方差分量估计，Partial EIV模型的优势是减少了待估参数的个数，如算例1直线拟合在不考虑方差分量估计时Partial EIV模型解算的迭代次数比EIV模型少，虽然对Partial EIV模型进行解算之后参数的表达式与文献[29]等价，但是使用更具一般性的Partial EIV模型更能体现其代表性。将Partial EIV模型与最小二乘方差分量估计结合可以得到与EIV模型一样的结果，但在参数估计上使用Partial EIV模型一定程度上提高运算效率，该方法继承了Partial EIV模型原有的优点，对总体最小二乘理论进行了必要的完善。本文只讨论了Partial EIV模型系数矩阵与观测向量不相关时出现负方差的处理，而针对相关观测情况及负方差产生的具体原因及分析是今后需要研究的工作。