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

王乐洋^1,2,3, 温贵森^1,2     

1. 东华理工大学测绘工程学院, 江西 南昌 330013; 
2. 流域生态与地理环境监测国家测绘地理信息局重点实验室, 江西 南昌 330013; 
3. 江西省数字国土重点实验室, 江西 南昌 330013

收稿日期：2016-10-11；修回日期：2017-05-27

基金项目：国家自然科学基金（41664001；41204003）；江西省杰出青年人才资助计划项目（20162BCB23050）；国家重点研发计划（2016YFB0501405）；江西省教育厅科技项目（GJJ150595）；江西省数字国土重点实验室开放研究基金资助项目（DLLJ201705）；东华理工大学研究生创新专项资金资助项目（DHYC-2016005）

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

摘要：Partial Errors-in-Variables（Partial EIV）模型是EIV模型的扩展形式，权阵构造简单，当系数矩阵中存在非随机元素和随机元素时，Partial EIV模型的适用性更强。针对Partial EIV模型中随机模型不准确的情况，将系数矩阵和观测向量分别作为一类数据，本文在该模型的基础上，使用最小二乘方差分量估计方法，推导相关计算公式及迭代算法，分别估计出相应的方差分量估值。并对出现的负方差使用非负最小二乘理论，增加约束条件，对随机模型进行修正，得到更加合理的参数估值。试实验结果表明，本文的方法与其他方差分量估计方法等价。

关键词：Partial EIV模型    EIV模型    最小二乘方差分量估计    非负最小二乘    

Non-negative Least Squares Variance Component Estimation of Partial EIV Model

WANG Leyang^1,2,3, WEN Guisen^1,2     

Abstract: As an extended form of the errors-in-variables (EIV) model, and the weight matrix is easy to structure, the applicability is stronger when both non-random elements and random elements exist in the coefficient matrix. According to the inaccurate stochastic model in the Partial EIV model, the coefficient matrix and observation vector are used as a kind of data respectively. The least squares variance component estimation method based on Partial EIV model is used and the relevant calculated formulas and iterative algorithm are derived. Then the corresponded variance components are estimated. The non-negative least squares is used when the negative variance appears, then add constraint condition to correct the rand model, so the estimated parameters are more reasonable. The experiments show that the results obtained by the method of this paper and other methods are equivalent.

Key words: Partial EIV model     EIV model     least squares variance component estimation     non-negative least squares    

总体最小二乘^[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]

函数模型

 (1)

随机模型

 (2)

式中，表示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) 进行变形得到

 (3)

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

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

 (4)

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

 (5)

经计算得到^[5, 26]

 (6)

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

 (7)

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

 (8)

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

 (9)

式中，表示投影阵。

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)

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

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

 (11)

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

 (12)

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

 (13)

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

 (14)

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

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

 (15)

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

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

 (16)

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

 (17)

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

 (18)

误差向量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]为

 (19)

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

 (20)

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

 (21)

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

 (22)

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

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

 (23)

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

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

 (24)

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

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

 (25)

式中，。

具体迭代步骤为：

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

(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。

表 1 坐标观测值及相应的权值^[5, 21]Tab. 1 Coordinate observations and corresponding weights^[5, 21]

表选项

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

 (26)

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

 (27)

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

表 2 算例1中不同方法的解算结果Tab. 2 Results from different methods of the first example

表选项

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

 (28)

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

 (29)

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

表 3 修正后的观测值权值Tab. 3 The observations weights of corrected

表选项

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。

表 4 原始坐标和目标坐标模拟坐标值Tab. 4 Simulated coordinates in start system and target system

表选项

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

 (30)

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

 (31)

式中，。

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

表 5 算例2中不同方法得到的结果Tab. 5 Results from different methods of the second example

表选项

图中，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模型系数矩阵与观测向量不相关时出现负方差的处理，而针对相关观测情况及负方差产生的具体原因及分析是今后需要研究的工作。

【引文格式】王乐洋，温贵森。Partial EIV模型的非负最小二乘方差分量估计[J]. 测绘学报，2017，46(7)：857-865. DOI: 10.11947/j.AGCS.2017.20160501



更多精彩内容：

科技部发布“重大自然灾害监测预警与防范”重点专项2017年度项目申报指南

《战狼2》撤侨到底发生在哪？破译电影里一幅近代地图的秘密

《测绘学报》2017年第7期网刊发布

2017无人机测绘技术及应用论坛在线报名已开通

惊艳的骗局：原来你每天坐的地铁线路图这样欺骗你的双眼！

论文推荐| 张力：基于多视影像匹配模型的倾斜航空影像自动连接点提取及区域网平差方法

论文推荐| 张春森：顾及曝光延迟的无人机GPS辅助光束法平差方法

2017中国地理信息产业优秀工程公示

78岁院士刘先林穿旧鞋坐高铁 二等座上仍不忘备课

如何评价 Google推出的室内定位技术VPS？

北斗专辑|张小红：BeiDou B2/Galileo E5b短基线紧组合相对定位模型及性能评估

想着为国家和人民做事，就充满活力——首个全国科技工作者日专访78岁李德仁双院士

武大测绘学院李星星教授获欧洲地学联合会杰出青年科学家奖

再不去把握这九大空间信息趋势, 等着被颠覆吧！