《测绘学报》

构建与学术的桥梁 拉近与权威的距离

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

王乐洋1,2,3, 许光煜1,4, 温贵森1,2

1. 东华理工大学测绘工程学院, 江西 南昌 330013;

2. 流域生态与地理环境监测国家测绘地理信息局重点实验室, 江西 南昌 330013;

3. 江西省数字国土重点实验室, 江西 南昌 330013;

4. 武汉大学测绘学院, 湖北 武汉 430079

收稿日期：2016-09-06；修回日期：2017-06-02

基金项目：国家自然科学基金（41664001；41204003）；江西省杰出青年人才资助计划项目（20162BCB23050）；国家重点研发计划（2016YFB0501405）；江西省教育厅科技项目（GJJ150595）

第一作者简介：王乐洋(1983-), 男, 博士, 副教授, 主要研究方向为大地测量反演及大地测量数据处理.E-mail:wleyang@163.com

摘要：Partial errors-in-variables（Partial-EIV）模型作为EIV模型的扩展形式，其构造方式更有规律，解算方法更为简便，能有效应用于实际情况。针对已有Partial EIV模型方法未考虑观测向量和系数矩阵存在相关性这一情况，通过提取观测向量和系数矩阵组成的增广矩阵中非重复出现的随机元素，构建更具一般适用性的Partial EIV模型，在该模型的基础上，将特殊假定条件扩展到不限定观测数据相关性的一般情况，详细推导了观测向量和系数矩阵元素相关且不等精度情况下的加权总体最小二乘方法，通过算例试验，并与目前已有的解决EIV模型相关观测情况下的方法进行了比较分析，研究表明本文方法可以提高计算效率，更具一般性，特别是对于观测向量和系数矩阵中存在常数元素和重复元素的情况。

A Method for Partial EIV Model with Correlated Observations

WANG Leyang1,2,3, XU Guangyu1,4, WEN Guisen1,2

Abstract: As an extended form of the errors-in-variables(EIV) model, partial errors-in-variables(Partial EIV) model has more advantages than the previous one, such as regular structure, simple solving method, which make it has a wide range of applications. Considering the situation that the correlation between the observations and elements in coefficient matrix is not taken into account in the existed algorithms derived from Partial EIV model, the non-repetitive random elements in the augmented matrix consisting of observation vector and coefficient matrix are extracted to build a more suitable partial EIV model. Based on this model, the special assumptions are extended to the general case where the observations are correlated, a new weighted total least squares(WTLS)algorithms is derived when the observations and elements in coefficient matrix are heteroscedastic and correlated. Through two examples, the algorithm proposed in this paper and the existed algorithms which consider the correlation of the observation in EIV model are compared and analyzed. Research shows that these algorithms can improve the calculation efficiency and more general, especially for the situation that coefficient matrix consists of constant elements and repeated elements.

Key words: total least squares correlated observations partial errors-in-variables autoregression model

在测量、信号处理和控制等学科中，为了求得某些未知参数，常常需要进行一系列的观测。由于人为观测存在许多不确定的因素，观测得到的数据通常只是未知量的某些函数，且得到的观测值夹杂着随机误差，如何利用含有误差的观测值求得未知参数估值是测量平差需要解决的问题。经典的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]，即

(1)

式中，y∈Rn×1表示观测向量；ey∈Rn×1表示观测向量y的随机误差；表示系数矩阵真值，h∈Rnm×1主要由系数矩阵A中非随机元素组成；B∈Rnm×t是固定矩阵，其阶数取决于系数矩阵A中随机元素的数目(此处假设为t)；a∈Rt×1是系数矩阵A中随机元素组成的列向量(若随机元素重复出现，则只提取一次)，其真值用

表示；ea∈Rt×1是a中包含的随机误差；β∈Rm×1表示未知参数向量。符号⊗表示Kronecker积，·β=(βT⊗In)vec；vec(·)表示矩阵拉直运算。ea和ey的协因数阵和权阵分别表示为Qa、Pa和Qy、Py，即

(2)

(3)

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

(4)

式中，cov表示协方差。

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

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

(5)

式中，Wy、Wya、Way、Wa为逆矩阵Q-1=

中的子块矩阵；Qy为观测向量的协因数阵；Qya为观测向量关于系数矩阵中随机元素的协因数阵。

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

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

(6)

式中，

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

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

(7)

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

(8)

式中，C2=[In-(βT⊗In)B]C1；W为向量e的权阵，且满足W=Qe-1；Qe为向量e的协因数阵；λ为拉格朗日乘子。

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

(9)

(10)

(11)

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

。

由式(9) 可得

(12)

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

(13)

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

(14)

将式(13) 代入式(10) 后进行变化，在y-(βT⊗In)(h+Ba)中同时加上和减去

，记y1=y+，则式(14) 又可以写成

(15)

迭代过程如下：

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

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

(3) 计算最小二乘结果

；

(4) 初始化C2=[In-(βT⊗In)B]C1；

(5) 通过式(12)—式(14) 计算得到

，更新矩阵C2；(6) 当

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

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

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

(16)

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

3.2 参数协因数阵

由式(14) 可知，相关观测的Partial EIV模型下的参数解

具有最小二乘解的形式，根据协因数传播律得到

，QEA表示对系数矩阵拉直后的协因数阵，QEA=BQaBT，且经过化简后Qy1=C2(W)-1C2T。因此采用经典最小二乘精度评定理论中参数协因数阵求解公式直接求得相关观测的Partial EIV模型下的参数协因数阵公式[16]

(17)

式中，Q1=C2(W)-1C2T。

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

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

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

图 1 Fang方法和Snow方法流程[12-13]Fig. 1 Fang's method and Snow's method[12-13]注：*QAA表示系数矩阵A的协因数阵；QyA表示观测向量y关于系数矩阵A的协因数阵；Qll表示观测向量和系数矩阵组成的增广矩阵[A,y]的协因数阵。

图选项

表 1 坐标观测数据及相应的权阵和残差[17]Tab. 1 Coordinate observations and corresponding weights as well as their residuals[17]

点号xiWxiyiWyieaey10.010005.91.00.002 611-0.543 92620.910005.41.8-0.008 749-0.452 01531.85004.44.0-0.012 180.136 25642.68004.68.00.015 800-0.443 97853.32003.520-0.074 3580.375 96564.4803.7200.156 868-0.435 37275.2602.870-0.038 5480.186 98686.1202.870-0.208 059-0.148 47696.51.82.4100-0.058 783-0.000 719107.41.01.55000.982 8340.007 635

表选项

表 2 坐标观测值之间的相关系数[13]Tab. 2 Correlation coefficient between coordinate observations[13]

相关系数ρx1y1ρx2y2ρx3y3ρx4y4ρx5y5ρx6y6ρx7y7ρx8y8ρx9y9ρx10y10数值-0.165 9560.440 649-0.999 771-0.395 335-0.706 488-0.815 323-0.627 480-0.308 879-0.206 4650.077 633

表选项

图 2 本文方法的流程图Fig. 2 The method in this paper

图选项

表 3 算例1不同方法所得结果及迭代次数Tab. 3 The results from different methods and the number of iterations of the first example

LSTLSFang方法Snow方法本文方法-0.610 812 97-0.480 533 41-0.459 228 67-0.459 228 67-0.459 228 676.100 109 325.479 910 225.357 272 565.357 272 565.357 272 564.293 150 941.483 294 152.090 685 972.090 685 972.090 685 970.003 886 390.004 987 220.005 969 980.005 969 980.005 969 980.179 826 420.129 058 060.156 490 230.156 490 230.156 490 23迭代次数—71097

表选项

假设直线方程为

(18)

直线拟合的Partial EIV模型为式(1)，式中

，固定矩阵

。

通过公式σxiyi=ρxiyiσxiσyi计算得到观测数据的协方差，相关情况下的协因数阵为(其中σ02=1)

(19)

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

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

表 4 算例1不同收敛阈值下不同算法的迭代次数Tab. 4 The number of iterations of different algorithms under different thresholds of the first example

收敛阈值Fang方法Snow方法本文方法10-787510-887610-998710-10109710-111110810-121210910-13131110

表选项

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

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

4.2 算例2自回归模型

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

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

(20)

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

(21)

式中

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

表 5 算例2不同方法的解算结果Tab. 5 The results from different methods of the second example

LSTLSFang方法Snow方法本文方法0.635 059 321.077 117 921.179 081 341.179 081 341.179 081 340.327 809 390.518 098 670.041 899 550.041 899 550.041 899 550.041 086 62-0.589 297 81-0.214 448 09-0.214 448 09-0.214 448 090.647 620 490.329 331 340.413 991 220.413 991 220.413 991 220.022 699 910.048 607 440.012 658 060.012 658 060.012 658 060.017 724 900.020 957 500.009 481 770.009 481 770.009 481 770.021 176 810.047 205 560.009 074 550.009 074 550.009 074 55迭代次数—334310933

表选项

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

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

表 6 算例2不同收敛阈值下不同方法的迭代次数Tab. 6 The number of iterations of different algorithms under different thresholds of the second example

收敛阈值Fang方法Snow方法本文方法10-732772410-836872710-939963010-10431093310-11461223610-12591445810-13746211 233286

表选项

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

4.3 算例3模拟直线拟合

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

(22)

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

表 7 算例3不同方法得到的参数估值均值与迭代次数均值Tab. 7 The average results of estimated parameters and the number of iterations from different methods of the third example

Fang方法Snow方法本文方法真值2.000 290 772.000 290 772.000 290 7721.496 678 241.496 678 241.496 678 241.5迭代次数均值4.8084.9523.808—

表选项

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

4.4 算例分析

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

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

(3) 在自回归模型算例中，文献[4, 15]中基于Partial EIV模型提出的方法无法使用，这是由于这两种方法的结构导致的，在自回归模型中，观测向量和系数矩阵中存在重复出现的元素，此时得到的观测向量和系数矩阵的增广矩阵[Ay]的协因数阵奇异[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模型方法的实际应用这方面还需要进一步研究。

【引文格式】王乐洋，许光煜，温贵森。一种相关观测的Partial EIV模型求解方法[J]. 测绘学报，2017，46(8)：978-987. DOI: 10.11947/j.AGCS.2017.20160430

更多精彩内容：