在测量、信号处理和控制等学科中，为了求得某些未知参数，常常需要进行一系列的观测。由于人为观测存在许多不确定的因素，观测得到的数据通常只是未知量的某些函数，且得到的观测值夹杂着随机误差，如何利用含有误差的观测值求得未知参数估值是测量平差需要解决的问题。经典的Gauss-Markov模型中仅考虑了观测向量y的随机误差，忽略或假定系数矩阵A不受随机误差的影响，采用最小二乘方法(least squares，LS)便可求得模型参数解。变量误差模型(errors-in-variables，EIV)中既考虑了观测向量y的随机误差，同时考虑到系数矩阵A也可能受到随机误差的影响，采用总体最小二乘方法(total least squares，TLS)，便可求得模型参数解^[1-2]。但在大地测量和工程测量的实际应用中^[3]，许多情况下系数矩阵只有部分含有随机误差，这类问题宜采用Partial EIV模型^[4]进行未知参数的求解。

Partial EIV模型由文献[4]提出，该模型相对于EIV模型更适合于解决实际问题。目前，针对Partial EIV模型的研究成果相当有限，文献[4]导出了Partial EIV模型的加权总体最小二乘方法(weighted total least squares，WTLS)，但该方法收敛速度较慢，所需的迭代次数过大；文献[5]基于Partial EIV模型导出了适用于海量数据处理问题的WTLS方法；文献[6]利用Partial EIV模型求解附有不等式约束的EIV模型问题。文献[7]在Partial EIV模型的基础上进行了方差分量估计问题的研究。但在以上研究中，均对随机模型进行了特定的假设即假定观测向量和系数矩阵中的元素不存在相关性^[8]。很明显，以上假设在实际问题中无法得到保证，因此如何能将Partial EIV模型的数据处理方法更加有效地应用于实际问题中，而不需要受到特殊假设的限制，也是本文主要研究的问题。本文将特殊假定条件扩展到不限定观测数据相关性的一般情况，并以Partial EIV模型为基础，详细推导了观测向量和系数矩阵元素相关且不等精度情况下的加权总体最小二乘方法；通过算例试验验证本文算法的可行性，并与目前已有的解决EIV模型相关观测情况下的方法^[9-13]进行了比较。

1 相关观测情形下的Partial EIV模型

不同的平差问题，函数模型中系数矩阵的结构也不同，EIV模型中系数矩阵的结构主要分为3类：① 由常数元素列和随机元素列组成，如直线拟合模型等；② 由随机元素零散分布组成，且随机元素重复出现，如自回归模型等；③ 由常数元素和随机元素零散分布组成，且随机元素重复出现，如坐标转换模型、地壳应变参数反演模型^[14]等。对于这3种情况，平差过程中必需保证常数元素不进行改正，重复出现的随机元素改正数相同。

文献[4]提出了Partial EIV模型，即式(1)，将上述情况都考虑在内，是一种概括性模型。但文献[4]中假定观测向量与系数矩阵中的观测数据不相关，下面探讨观测数据相关性情况下的Partial EIV模型。在实际测量数据处理中，系数矩阵的元素并非全都是随机的，如线性回归模型和坐标转换模型，部分是由0或者1等常数组成，此时的平差函数模型即为Partial EIV模型^[4]，即

式中，y∈R^n×1表示观测向量；e_y∈R^n×1表示观测向量y的随机误差；表示系数矩阵真值，h∈R^nm×1主要由系数矩阵A中非随机元素组成；B∈R^nm×t是固定矩阵，其阶数取决于系数矩阵A中随机元素的数目(此处假设为t)；a∈R^t×1是系数矩阵A中随机元素组成的列向量(若随机元素重复出现，则只提取一次)，其真值用表示；e_a∈R^t×1是a中包含的随机误差；β∈R^m×1表示未知参数向量。符号⊗表示Kronecker积，·β=(β^T⊗I_n)vec()；vec(·)表示矩阵拉直运算。e_a和e_y的协因数阵和权阵分别表示为Q_a、P_a和Q_y、P_y，即

假设观测值与系数矩阵中随机元素是相关的，即

式中，cov表示协方差。

2 相关观测的Partial EIV模型求解方法

通过观测值y和系数矩阵随机元素a估计Partial EIV模型中和β的加权总体最小二乘平差准则为^[15]

式中，W_y、W_ya、W_ay、W_a为逆矩阵Q^-1=中的子块矩阵；Q_y为观测向量的协因数阵；Q_ya为观测向量关于系数矩阵中随机元素的协因数阵。

当观测向量y和系数矩阵随机元素a中包含有相同的观测值时(例如自回归模型)，平差准则式(5) 中的矩阵Q奇异，此时无法得到它的逆矩阵，进而无法利用平差准则式(5) 进行未知参数的求解。此时需要对原有的Partial EIV模型进行一些调整，并根据扩展的模型重新构造目标函数进行参数求解。调整方式如下

将式(1) 进行改写，得到函数模型

式中，；e是由观测向量y和系数矩阵中随机元素a组成的列向量中非重复元素的随机误差组成的s×1列向量；C₁∈R^(n+t)×s为给定的转换矩阵。

以观测数据随机误差的平方和最小为平差准则

根据式(7) 构造拉格朗日条件极值的目标函数

式中，C₂=[I_n-(β^T⊗I_n)B]C₁；W为向量e的权阵，且满足W=Q_e^-1；Q_e为向量e的协因数阵；λ为拉格朗日乘子。

对目标函数式(8) 求偏导并令其等于零得

式(10) 中，Ivec表示矩阵拉直的逆变换，并记。

由式(9) 可得

将式(12) 代入式(11) 得

将式(13) 代入式(10)，得

将式(13) 代入式(10) 后进行变化，在y-(β^T⊗I_n)(h+Ba)中同时加上和减去，记y₁=y+，则式(14) 又可以写成

迭代过程如下：

(1) 给定已知或者观测数据a、y、P_y、Q_e；

(2) 根据函数模型中系数矩阵结构确定矩阵B和向量h, 以及构造转换矩阵C₁；

(3) 计算最小二乘结果；

(4) 初始化C₂=[I_n-(β^T⊗I_n)B]C₁；

(5) 通过式(12)—式(14) 计算得到，更新矩阵C₂；

(6) 当与前一次迭代值的差值2范数小于设定的阈值，迭代终止，得到最终的参数解，否则跳到步骤(5)，重复迭代计算。

3 相关观测的Partial EIV模型算法精度评定 3.1 单位权方差估值公式

根据文献[8]，相关观测的Partial EIV模型算法的单位权方差估值公式为

式中, r为自由度，其数值大小与观测值个数和待求参数个数有关，在相关观测的Partial EIV模型中，观测值个数为n+t；待求参数个数为t+m。

3.2 参数协因数阵

由式(14) 可知，相关观测的Partial EIV模型下的参数解具有最小二乘解的形式，根据协因数传播律得到，Q_EA表示对系数矩阵拉直后的协因数阵，Q_EA=BQ_aB^T，且经过化简后Q_y1=C₂(W)^-1C₂^T。因此采用经典最小二乘精度评定理论中参数协因数阵求解公式直接求得相关观测的Partial EIV模型下的参数协因数阵公式^[16]

式中，Q₁=C₂(W)^-1C₂^T。

其他相关量的协因数阵可以通过协因数传播律求得。

4 算例和分析 4.1 算例1直线拟合

引用文献[17]中的试验数据，坐标观测值(x_i, y_i)及其相应的权W_xi、W_yi见表 1。为了模拟系数矩阵A和观测向量y相关的情况，利用表 2中的相关系数^[13]给观测值加入相关性。分别采用最小二乘方法(LS)、不考虑相关性的总体最小二乘(TLS)^[16]、本文方法和文献[12]的方法(记Fang方法)、文献[13]的方法(记Snow方法)进行参数求解，图 1给出了Fang方法和Snow方法的流程图，图 2给出了本文方法的流程图。参数计算结果见表 3。表 1给出了利用本文方法计算得到的观测向量和系数矩阵随机元素的残差值。本文将迭代的收敛阈值均设为10^-10。

假设直线方程为

直线拟合的Partial EIV模型为式(1)，式中，固定矩阵。

通过公式σ_xiyi=ρ_xiyiσ_xiσ_yi计算得到观测数据的协方差，相关情况下的协因数阵为(其中σ₀²=1)

从表 3可以看出，本文方法与Fang方法、Snow方法可以得到相同的参数估值、单位权方差及参数协方差，不同的是本文方法的迭代次数要少于其他两种方法；参数估值表达式的不同是影响迭代次数的主要原因。

目前已有的加权总体最小二乘算法中，文献[8, 16, 18-19]在迭代中收敛阈值都设置为10^-10，文献[20]收敛阈值设置为10^-12，文献[21]收敛阈值设置为10^-18，收敛阈值设置为10^-10较为普遍；为了进一步说明本文方法的优势，在算例1中通过设置不同的收敛阈值比较不同方法的迭代次数，得到的结果见表 4。

从表 4可以看出，收敛阈值的不同对迭代次数存在影响，通过比较可以发现本文方法的迭代次数在不同阈值下都是最少的。

在该直线拟合算例中，数据的相关性存在于观测向量和系数矩阵的随机元素中，从相关系数中可以看出，各系数矩阵中随机元素x_i对观测向量中元素y_i所起作用并非相同。数据之间的相关性主要从协因数阵中体现出来，若不考虑相关性，则式(18) 为对角阵，很明显若缺少了非对角线元素的约束，则无法体现出数据之间的相互作用关系，最终导致求得的参数结果存在偏差。

4.2 算例2自回归模型

时间序列法是一种利用按时间顺序排列的数据预测未来的方法，它通过对各种动态数据来建立对应的数学模型，并对所建立的模型进行分析研究，从而了解这些观测数据的特性和相互联系，对未来数据的变化趋势作出合理的分析和预测。自回归模型的参数估计问题是时间序列求解中的基本问题^[22-23]。

p阶自回归过程{Y_n}满足^[22]

将式(19) 改写成矩阵的形式

式中

采用文献[24]中的观测数据，该数据为某一建筑物在某个时间段定期进行的36次沉降观测。建立一个阶数为3的自回归模型，分别采用最小二乘方法(LS)、不考虑相关性的总体最小二乘方法(TLS)^[16]、本文方法和Fang方法^[12]、Snow方法^[13]进行参数求解，参数计算结果见表 5。本文将迭代的收敛阈值均设为10^-10。在自回归模型中，数据之间的相关性是由于存在重复元素而产生的，系数矩阵中存在元素重复出现的情况，且观测向量与系数矩阵中存在相同元素，导致元素之间存在自相关，此时在进行模型解算时需要保证重复元素具有相同的改正数。若不考虑相关性，则会出现同一元素具有不同的改正数，与实际情况不相符，最终导致参数解算结果存在偏差。

从表 5可以看出，本文方法与Fang方法、Snow方法计算结果的等价性，不同的是迭代次数，本文方法在一定程度上可以提高计算效率。

同算例1，通过设置不同的收敛阈值来比较不同方法的迭代次数，得到的结果见表 6。

从表 6可以看出，迭代次数随着收敛阈值的变小而增多，当收敛阈值取为10^-13时，Fang方法与Snow方法迭代次数分别为7462、11 233，而本文方法只需要286次，在不同的收敛阈值下本文方法的迭代次数都是最少的，体现了本文方法可以提高计算效率这一优势。

4.3 算例3模拟直线拟合

在算例1直线拟合中，由式(18) 可以看出未考虑观测量自身与系数矩阵观测量自身的相关性，其协因数阵为对角阵，互协因数阵也为对角阵；因此为了说明本文方法的一般性，考虑观测量自身的相关性，则式(18) 形式变为

模拟一条直线y_i=ax_i+b, a与b分别表示斜率与截距，其真值为2、1.5，x_i从2至9.2每隔0.8取一个值，并计算相应的y_i值；其中y_i与x_i的协因数阵及y_i与x_i相关的互协因数阵由Matlab函数randn生成，其形式为Q={[randn(100, 2*n)]^T×[randn(100, 2*n)]})/1000，令单位权方差σ₀²=0.1，再根据Matlab函数mvnrnd函数生成随机误差进行平差解算，保持协因数阵Q及单位权方差不变模拟计算1000次，迭代中收敛阈值设置为10^-10，用Fang方法、Snow方法及本文方法计算得到参数估值均值与迭代次数均值结果见表 7。

从表 7可以看出，当考虑观测量自身的相关性及观测量与系数矩阵的相关性时，协因数阵式(22) 与式(19) 形式相同，使用不同的方法可以得到相同的计算结果，但本文方法的迭代次数均值要少于Fang方法与Snow方法；在构造协因数阵方面，本文方法只需考虑随机元素的协因数阵及它们之间的互协因数阵，简单方便，而Fang方法与Snow方法需要获得增广矩阵[A y]的协因数阵及互协因数阵，维度较大，比本文更复杂，尤其是考虑了观测量自身的相关性及观测量与系数矩阵元素的相关性时，本文方法一定程度上简化了计算程序。

4.4 算例分析

(1) 从表 3、表 5和表 7中可以看出，本文提出的解决相关观测的Partial EIV模型方法与已有的解决EIV模型相关观测情况下的方法在参数解算方面具有相同的效果，无论是在直线拟合或是自回归模型问题中，求解得到的参数结果一致，计算得到的单位权方差也一致，而且迭代次数要少于已有的算法，一定程度上提高了计算效率，从表 4和表 6可以看出，设置不同的收敛阈值，本文方法的迭代次数都是最少的，进一步说明了本文方法的优势。相比较其他相关观测的WTLS方法，本文提供了一种可供选择的其他方法。主要体现Partial EIV模型的优势，将其用于解决观测值相关情况下的总体最小二乘问题。

(2) 由于直线拟合模型的系数矩阵是由常数元素列和随机元素列组成的，系数矩阵中并不存在随机元素重复出现的情况，因此构造EIV模型或是Partial EIV模型进行参数求解并无太大差异，但算例2自回归模型中含有重复的随机元素明显地体现出了本文方法的优势。

(3) 在自回归模型算例中，文献[4, 15]中基于Partial EIV模型提出的方法无法使用，这是由于这两种方法的结构导致的，在自回归模型中，观测向量和系数矩阵中存在重复出现的元素，此时得到的观测向量和系数矩阵的增广矩阵[A y]的协因数阵奇异^[25]，因此无法得到该增广矩阵的权阵，而文献[4, 15]中均需要提供增广矩阵的权阵。这也是本文针对目前已有Partial EIV模型方法存在的不足之处所进行的改进。

(4) 本文是在Partial EIV模型的基础上对该模型进行了扩展，在标准的Partial EIV模型中，仅对系数矩阵中的常数元素和随机元素进行了分离，在本文方法中则是将观测向量和系数矩阵组成的增广矩阵中的常数元素和随机元素进行了分离，依照Partial EIV模型的构建思想，对该模型进行了扩展。从表 5中可以看出，本文方法在求解自回归模型参数时只需迭代33次，而Fang方法需要43次，Snow方法需要109次，这主要是由不同的参数求解表达式引起的。

(5) 在本文方法中，遵循了Partial EIV模型的构造思想，将重复出现的元素进行了提取，因此在方法求解过程中，无论一个元素重复出现多少次，都不影响求解参数的个数，而Fang方法和Snow方法中对增广矩阵中出现的元素进行了求解，因此在求解参数的个数方面远多于本文方法所需求解的参数个数，因此这两种方法的计算量更大。同时，必需注意在平差模型构建过程中元素的重复出现，只是同一个观测量的重复使用，对其观测精度并无影响。

(6) 通过本文方法与图 1中Fang方法和Snow方法的迭代表达式的对比，可以发现本文方法表达形式更为简便，组成结构更有规律。由于本文方法是基于扩展的Partial EIV模型推导得到的，因此继承了Partial EIV模型的优势；从式(19)、式(22) 可以看出，本文方法解决了已有Partial EIV方法未考虑相关性的问题。该模型考虑到方法的集成性发展，根据不同平差问题预先设计固定矩阵B的形式，可以形成一个统一的求解模式，使其具有更为广阔的实际应用前景。

5 结论

本文以Partial EIV模型为基础，推导了适用于相关观测情况下的方法，进一步完善了Partial EIV模型的理论研究。该方法在结构和理论推导方面更加简便、易懂，且更适合解决系数矩阵中存在元素重复出现的问题，对总体最小二乘理论的发展具有一定的理论与实际价值。通过算例表明了本文方法与目前已有的解决EIV模型相关观测情况下的方法在参数求解结果方面的一致性，相较于EIV模型下的方法，本文方法更具普适性，可以形成统一的求解模式，实际应用前景更为广阔。本文还比较了不同方法的迭代次数，并对同一问题中不同方法具有不同的迭代次数进行了分析，但未能给出更为全面的解释，对方法机理的进一步探究能够加深对方法迭代过程的深入理解，有利于提高方法的效果，对于这一问题以及Partial EIV模型方法的实际应用这方面还需要进一步研究。