部分变量误差^[1](partial errors-in-variables，PEIV)模型是针对变量误差^{[2, 3]}(errors-in-vari-ables，EIV)模型中系数矩阵只有部分元素存在误差而提出的模型，将系数矩阵中的随机元素提取，减少了模型中待改正量的个数；顾及了系数矩阵误差的EIV模型近几年被广泛应用^[4-8]，相比于只考虑观测量误差的Gauss-Markov(GM)模型，EIV模型在某些程度上提高了参数的解算精度。自文献[1]提出PEIV模型以来，文献[9]推导了适用于海量数据处理的加权总体最小二乘(total least squares, TLS)方法，文献[10]使用PEIV模型求解附有不等式约束的EIV模型问题，文献[11]对PEIV模型进行了移项，将参数表达式写成标准化的最小二乘(least squares, LS)形式，简化了文献[1]的参数表达式，且在运算速率上有所提高，文献[12]给出了改进目标函数的PEIV模型解法，文献[13]研究了PEIV的加权总体最小二乘算法及其应用，文献[14, 15]研究了PEIV模型的方差分量估计(variance-covariance component estimation, VCE)，文献[16, 17]分别针对多项式拟合及水准测量尺度比参数的附加系统参数模型的特点，使用PEIV模型进行解算，扩展了模型的应用范围，文献[18]研究了附有相对权比的PEIV模型平差方法。上述文献推导都是在观测量和系数矩阵不相关且只考虑相同的方差分量基础上得到的，文献[19]提出了一种相关观测情况下的PEIV模型解法，处理了相关观测情况的解算。观测量与系数矩阵元素之间相关的模型有自回归模型^{[20, 21]}、高程拟合模型^[22]等。

平差模型包括函数模型与随机模型，一个合理的随机模型对平差结果有很大的影响，然而初始平差时给定的随机模型往往是不准确的，方差分量估计VCE方法^[23-27]可以对随机模型进行修正。针对最小二乘，文献[28]分析了方差-协方差分量估计的统一理论，文献[29]针对GPS双差观测值的方差-协方差阵进行了方差-协方差分量估计，文献[30]分析了Helmert型方差-协方差的通用公式，文献[31]基于概括函数平差模型分析了方差分量估计的通用公式；而针对总体最小二乘，文献[14]分析了PEIV模型的Helmert方差分量估计(Helmert-VCE)方法，文献[32]分析了EIV模型的最小二乘方差分量估计(LS-VCE)^[33-35]方法，文献[15]分析了PEIV模型的非负最小二乘方差分量估计，处理了方差分量估计中出现的负方差情况。Helmert方差分量估计方法与LS-VCE方法可以得到相同的参数估值，然而以上研究并未考虑观测向量与系数矩阵相关的情况，即随机模型的方差-协方差分量估计。

基于上述分析，本文针对相关观测且随机模型不准确情形下的PEIV模型，根据最小二乘方差分量估计原理，对PEIV模型的随机模型进行分析，转换得到3个方差分量的表达式，并给出迭代算法；对公式进行进一步转换，且与已有的方差分量估计公式进行对比，得到更具一般性的表达式，并通过算例说明了本文方法的可行性。

1 相关PEIV模型解算

PEIV的数学模型表示如下：函数模型为：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{\tilde A\beta }} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}} = ({\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n})(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\tilde a}}) - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}}\\ {\mathit{\boldsymbol{a}} = \mathit{\boldsymbol{\tilde a}} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}} \end{array}} \right. $

(1)

随机模型为：

$ \left[\begin{array}{l} {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}\\ {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}} \end{array} \right] \sim (\mathit{\boldsymbol{0}}, \sigma _0^2\left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}} \end{array}} \right]) $

(2)

式(1)、式(2)中，y∈R^n×1表示观测向量；e_y∈R^n×1表示观测向量的随机误差；$ \mathit{\boldsymbol{\tilde A}} \in {\mathit{\boldsymbol{R}}^{n \times m}} $表示系数矩阵真值；h∈R^nm×1由系数矩阵中的非随机元素组成；B∈R^nm×t表示固定矩阵，t表示系数矩阵随机元素的个数；a∈R^t×1表示系数矩阵中随机元素组成的列向量，$ \mathit{\boldsymbol{\tilde a}} $表示其真值，e_a∈R^t×1表示其随机误差；β∈R^m×1表示未知参数向量；符号⊗表示Kronecker积；$ \mathit{\boldsymbol{\tilde A\beta }} = ({\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{n}}}){\rm{vec}}(\mathit{\boldsymbol{\tilde A}}) $，且$ {\rm{vec}}(\mathit{\boldsymbol{\tilde A}}) = \mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\tilde a}} $，vec(·)表示矩阵拉直运算；Q_y、Q_a、Q_ay分别表示观测量y、a的协因数阵及互协因数阵；σ₀²表示单位权方差。

文献[36]针对系数矩阵中的元素不再是自变量本身而是关于变量的函数的情况提出了非线性总体最小二乘平差算法，对模型进行线性化求解。根据文献[36]的求解思想，PEIV模型可以视为非线性函数，将其表示为非线性GM模型：

$ \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{y}}\\ \mathit{\boldsymbol{a}} \end{array}} \right] + \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] = \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}) $

(3)

式中，$ f\left( \mathit{\boldsymbol{x}} \right) = \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde A\beta }}}\\ {\mathit{\boldsymbol{\tilde a}}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\left( {\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\tilde a}}} \right)}\\ {\mathit{\boldsymbol{\tilde a}}} \end{array}} \right];\mathit{\boldsymbol{x}} = \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{\beta }}\\ {\mathit{\boldsymbol{\tilde a}}} \end{array}} \right] $。

考虑到模型的非线性，将模型在参数近似值x₀处展开，$ \mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{x}}_0} + \Delta \mathit{\boldsymbol{x}} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\beta }}_0}}\\ {{\mathit{\boldsymbol{a}}_0}} \end{array}} \right] + \left[{\begin{array}{*{20}{c}} {\Delta \mathit{\boldsymbol{\beta }}}\\ {\Delta \mathit{\boldsymbol{a}}} \end{array}} \right] $，则式(3)可写为：

$ \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{y}}\\ \mathit{\boldsymbol{a}} \end{array}} \right] + \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] = f({\mathit{\boldsymbol{x}}_0}) + \mathit{\boldsymbol{\hat J}}\Delta \mathit{\boldsymbol{x}} $

(4)

式中，$ f({\mathit{\boldsymbol{x}}_0}) = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_0}{\mathit{\boldsymbol{\beta }}_0}}\\ {{\mathit{\boldsymbol{a}}_0}} \end{array}} \right] $；$ \mathit{\boldsymbol{\hat J}} $为f(x)在x₀^T处的Jacobian矩阵，其形式为：

$ \mathit{\boldsymbol{\hat J}} = \frac{{\partial f(\mathit{\boldsymbol{x}})}}{{\mathit{\boldsymbol{x}}_{\rm{0}}^{\rm{T}}}} = \left[{\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat A}}}_i}}&{{\rm{(}}\mathit{\boldsymbol{\hat \beta }}_i^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_n}{\rm{)}}\mathit{\boldsymbol{B}}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{I}}_t}} \end{array}} \right] $

(5)

式中，$ {{{\mathit{\boldsymbol{\hat A}}}_i}} $与$ \mathit{\boldsymbol{\hat \beta }}_i $表示第i次迭代计算时的系数矩阵估值与参数估值，当i=0时，表示在x₀处展开得到的初值。

将式(4)改写成最小二乘间接平差形式有：

$ \mathit{\boldsymbol{\hat L}} + \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{\hat J}}\Delta \mathit{\boldsymbol{\hat x}} $

(6)

式中，$ \mathit{\boldsymbol{\hat L}} = \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{y}}-{\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{\beta }}_i}}\\ {\mathit{\boldsymbol{a}}-{\mathit{\boldsymbol{a}}_i}} \end{array}} \right], \mathit{\boldsymbol{e}} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] $，其协因数阵$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] $。

对式(6)进行最小二乘解算得：

$ \Delta \mathit{\boldsymbol{\hat x}} = {({\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})^{ - 1}}{\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat L}} $

(7)

通过迭代计算得到参数估值为：

$ {\mathit{\boldsymbol{\hat x}}^{i + 1}} = {\mathit{\boldsymbol{\hat x}}^i} + \Delta \mathit{\boldsymbol{\hat x}} $

(8)

观测值改正数的表达式为^[33]：

$ \mathit{\boldsymbol{\hat e}} = \mathit{\boldsymbol{\hat J}}\Delta \mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\hat L}} = \mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot \mathit{\boldsymbol{\hat L}} $

(9)

式中，$ \mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot = {\mathit{\boldsymbol{I}}_{n + t}} - \mathit{\boldsymbol{\hat J}}({\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}){\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1} $表示为投影阵，其满足$ \mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot = {(\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot )^{\rm{T}}} $。

相应的迭代算法为：

1) 给定数据A、Q_a、Q_ya、Q_y、h、B、a；

2) 计算最小二乘初值$ {\mathit{\boldsymbol{\beta }}_{{\rm{LS}}}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}^{ - 1}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}^{ - 1}\mathit{\boldsymbol{y}} $，令$ {\mathit{\boldsymbol{a}}_0} = \mathit{\boldsymbol{a}}, {\mathit{\boldsymbol{\beta }}_0} = {\mathit{\boldsymbol{\beta }}_{{\rm{LS}}}}, {\mathit{\boldsymbol{x}}^i} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\beta }}_0}}\\ {{\mathit{\boldsymbol{a}}_0}} \end{array}} \right] $；

3) 构造Jacobian矩阵$ \mathit{\boldsymbol{\hat J}} $及观测量$ \mathit{\boldsymbol{\hat L}} $，用式(7)计算$ \Delta \mathit{\boldsymbol{\hat x}} $；

4) 利用式(8)计算$ {\mathit{\boldsymbol{\hat x}}^{i + 1}} $，并令$ {\mathit{\boldsymbol{\hat x}}^i} = {\mathit{\boldsymbol{\hat x}}^{i + 1}} $，重复步骤3)，直至$ \left\| {{{\mathit{\boldsymbol{\hat x}}}^{i + 1}}\mathit{\boldsymbol{ - }}{{\mathit{\boldsymbol{\hat x}}}^i}} \right\| $小于给定阈值。

从该算法可以看出，与已有的迭代算法相比，该迭代步骤比较简单，易于实现。

2 最小二乘方差分量估计

最小二乘方差分量估计(LS-VCE)方法是将方差分量作为最小二乘平差参数进行解算，得到的结果具有无偏性，且更易于精度评定，将PEIV模型进行泰勒展开后，可以得到最小二乘间接平差形式，但此时需要进行迭代计算。由式(6)可得：

$ E\left\{ {(\mathit{\boldsymbol{\hat L}} - \mathit{\boldsymbol{\hat J}}\Delta \mathit{\boldsymbol{\hat x}}){{(\mathit{\boldsymbol{\hat L}} - \mathit{\boldsymbol{\hat J}}\Delta \mathit{\boldsymbol{\hat x}})}^{\rm{T}}}} \right\} = \sum\limits_{i = 1}^p {\sigma _i^2{\mathit{\boldsymbol{Q}}_i}} $

(10)

考虑平差随机模型为：

$ {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{e}}} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{y}}}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{ya}}}}}\\ {{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] = {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}} = \left[{\begin{array}{*{20}{c}} {\sigma _1^2{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}}&{\sigma _3^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}\\ {\sigma _3^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] $

(11)

式中，Q_e此时表示的是协方差阵；σ₁²、σ₂²、σ₃²表示为未知方差分量，因此将式(11)改写成线性求和形式：

$ {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{e}}} = {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}} = \sum\limits_{i = 1}^3 {\sigma _i^2{\mathit{\boldsymbol{Q}}_i}} $

(12)

式中，Q_i的形式为：

$ {\mathit{\boldsymbol{Q}}_i} = \left\{ \begin{array}{l} \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}} \end{array}} \right], i = 1\\ \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}} \end{array}} \right], i = 2\\ \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&\mathit{\boldsymbol{0}} \end{array}} \right], i = 3 \end{array} \right. $

(13)

在式(7)的迭代计算中，$ {\mathit{\boldsymbol{\hat Q}}_e} $为初始的协因数阵，通过方差分量估计确定未知方差分量估值来修正初始平差时所给定的权，从而达到修正参数估值的效果；存在矩阵H满足$ {\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{\hat J}} = 0 $或$ \mathit{\boldsymbol{H}}{\mathit{\boldsymbol{\hat J}}^{\rm{T}}} = 0 $，在式(10)等号两边左乘H^T和右乘H得到：

$ E\left\{ {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{\hat L}}{{({\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{\hat L}})}^{\rm{T}}}} \right\} = \sum\limits_{i = 1}^p {\sigma _i^2{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_i}\mathit{\boldsymbol{H}}} $

(14)

令$ \mathit{\boldsymbol{t}} = {\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{\hat L}} $，并对式(14)作vh变换(表示取出矩阵对角线的下三角元素并拉直)得到：

$ E({\mathit{\boldsymbol{L}}_{{\rm{vh}}}}) = {\mathit{\boldsymbol{J}}_{{\rm{vh}}}}\mathit{\boldsymbol{\sigma }} $

(15)

式中，L_vh=vh(tt^T)；σ为方差分量组成的向量；$ {\mathit{\boldsymbol{J}}_{{\rm{vh}}}} = [{\rm{vh}}({\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{H}}), {\rm{vh}}({\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{H}}), \cdots, {\rm{vh}}({\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_n}\mathit{\boldsymbol{H}})] $。

对式(15)进行最小二乘解算得到：

$ \mathit{\boldsymbol{\hat \sigma }} = {\rm{(}}\mathit{\boldsymbol{J}}_{{\rm{vh}}}^{\rm{T}}{\mathit{\boldsymbol{W}}_{{\rm{vh}}}}{\mathit{\boldsymbol{J}}_{{\rm{vh}}}}{\rm{)}}\mathit{\boldsymbol{J}}_{{\rm{vh}}}^{\rm{T}}{\mathit{\boldsymbol{W}}_{{\rm{vh}}}}{\mathit{\boldsymbol{L}}_{{\rm{vh}}}} = {\mathit{\boldsymbol{N}}^{ - 1}}\mathit{\boldsymbol{l}} $

(16)

考虑方差分量的权阵W_vh与向量t权阵的关系$ {\mathit{\boldsymbol{W}}_{{\rm{vh}}}} = {\mathit{\boldsymbol{D}}^{\rm{T}}}{\rm{(}}{\mathit{\boldsymbol{W}}_\mathit{\boldsymbol{t}}} \otimes {\mathit{\boldsymbol{W}}_\mathit{\boldsymbol{t}}}{\rm{)}}\mathit{\boldsymbol{D}} $，其中矩阵D满足vec(W_t)= Dvh(W_t)，将矩阵J_vh、L_vh、W_vh代入式(16)并进行转换得到^[34]：

$ \left\{ \begin{array}{l} {n_{kj}} = {\rm{tr}}(\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{k}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{l}}})\begin{array}{*{20}{c}} {, k, l = 1, ..., p}&{}&{\begin{array}{*{20}{c}} {}&{}&{} \end{array}}&{} \end{array}\\ {l_k} = {{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{k}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat e}}\begin{array}{*{20}{c}} {, k = 1, ..., p}&{} \end{array} \end{array} \right. $

(17)

考虑到相关观测情形下的PEIV模型，结合式(12)得到法矩阵N及向量l的形式为：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{N}} = \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{-1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{-1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_1}}&{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{-1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_2}}&{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_3}}\\ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_1}}&{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_2}}&{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_3}}\\ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_1}}&{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_2}}&{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot {\mathit{\boldsymbol{Q}}_3}} \end{array}} \right]\\ \mathit{\boldsymbol{l}} = {\left[{\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{-1}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{-1}\mathit{\boldsymbol{\hat e}}}&{{{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{-1}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat e}}}&{{{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}{\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat e}}} \end{array}} \right]^{\rm{T}}} \end{array} \right. $

(18)

文献[37]给出了最小二乘平差的最小范数二次无偏估计(the minimum norm quadratic unbiased estimator，MINQUE)方法，针对相关观测PEIV模型，将其化简为式(6)的形式可认为线性GM模型，因此MINQUE的公式为：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{S}} = {\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_k}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_l}) = \\ \left[{\begin{array}{*{20}{c}} {{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1})}&{{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_2})}&{{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_3})}\\ {{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1})}&{{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_2})}&{{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_3})}\\ {{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1})}&{{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_2})}&{{\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_3})} \end{array}} \right]\\ \mathit{\boldsymbol{W}} = {\left[{\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{C\hat e}}}&{{{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{e}}^{-1}\mathit{\boldsymbol{\hat e}}}&{{{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_3}\mathit{\boldsymbol{C\hat e}}} \end{array}} \right]^{\rm{T}}} \end{array} \right. $

(19)

式中，$ \mathit{\boldsymbol{C}} = \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1} - \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}{({\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})^{ - 1}}{\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1} $，等价于LS-VCE中的$ \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot = \mathit{\boldsymbol{C}} $，因此矩阵S与N是等价的。

顾及式(9)中观测值改正数$ \mathit{\boldsymbol{\hat e}} = \mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot \mathit{\boldsymbol{\hat L}} $，且$ \mathit{\boldsymbol{\hat L}} = \mathit{\boldsymbol{\hat J}}\Delta \mathit{\boldsymbol{\hat x}} - \mathit{\boldsymbol{\hat e}} $，又存在$ \mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot \mathit{\boldsymbol{\hat J}} = 0 $，即$ \mathit{\boldsymbol{\hat e}} = \mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot \mathit{\boldsymbol{\hat e}} $，将其代入式(9)中得：

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat e}} = {(\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot \mathit{\boldsymbol{\hat e}})^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}(\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{\hat J}}}^ \bot \mathit{\boldsymbol{\hat e}}) = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{\mathit{\boldsymbol{\hat e}}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{C\hat e}} \end{array} $

(20)

因此可得到N = S，l = W，LS-VCE与MINQUE是等价的。

考虑到式(12)中Q_i的形式，即存在$ \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{i}}} = {\mathit{\boldsymbol{I}}_i} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_n}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}} \end{array}} \right] $，令$ \mathit{\boldsymbol{N}} = ({\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}), {\mathit{\boldsymbol{N}}_i} = {({\mathit{\boldsymbol{\hat J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})_i} $，根据文献[37]的推导可得：

$ \left\{ \begin{array}{l} {S_{ii}} = {\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_i}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_i}) = {\rm{tr}}\left\{ {({\mathit{\boldsymbol{I}}_i} - \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}^{ - 1}}\mathit{\boldsymbol{J}}_i^{\rm{T}})({\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{n}}} - \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})\mathit{\boldsymbol{J}}_i^{\rm{T}})} \right\}\\ \;\;\;\; = {n_i} - 2{\rm{tr}}\left\{ {{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}^{ - 1}}{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}_i}} \right\} + {\rm{tr}}{\left\{ {{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}^{ - 1}}{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}_i}} \right\}^2}\\ \;\;\;\; = {n_i} - 2{\rm{tr(}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}_i}{\rm{)}} + {\rm{tr(}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}_i}{{\rm{)}}^{\rm{2}}}\\ {S_{ij}} = {\rm{tr}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_i}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_i}) = {\rm{tr}}\left\{ {({\mathit{\boldsymbol{I}}_i} - \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}^{ - 1}}\mathit{\boldsymbol{J}}_i^{\rm{T}})({\mathit{\boldsymbol{I}}_j} - \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})\mathit{\boldsymbol{J}}_j^{\rm{T}})} \right\}\\ \;\;\;\; = {\rm{tr}}\left\{ {{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}^{ - 1}}{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}_i}{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}^{ - 1}}{{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})}_j}} \right\} \end{array} \right. $

(21)

式中，i, j=1, 2, …, p。同理，将矩阵C代入到$ {\mathit{\boldsymbol{\hat e}}^{\rm{T}}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_i}\mathit{\boldsymbol{C\hat e}} $中，且$ {\mathit{\boldsymbol{J}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat e}} = 0 $，可以得到：

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_i}\mathit{\boldsymbol{C\hat e}} = {{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}(\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1} - \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})^{ - 1}}{{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}) \times \\ \;\;\;\;\;\;\;{\mathit{\boldsymbol{Q}}_i}(\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1} - \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}}{({{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat J}})^{ - 1}}{{\mathit{\boldsymbol{\hat J}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1})\mathit{\boldsymbol{\hat e}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}{\mathit{\boldsymbol{Q}}_i}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat e}} = {{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_i^{ - 1}\mathit{\boldsymbol{\hat e}} \end{array} $

(22)

式中，Q_i的形式见式(13)。

由式(18)~式(22)可以看出，LS-VCE、MINQUE和Helmert方法是等价的，若不考虑两类观测量的相关性，即式(12)中Q₃为零矩阵，此时进行最小二乘迭代计算方法与已有文献是等价的，且该方差分量估计公式与文献[32]方法是等价的，计算出方差分量后，可以对观测量权进行修正，从而修正参数估值。相关观测PEIV模型的最小二乘方差分量估计的迭代流程图见图 1。

3 算例分析 3.1 算例1

数据采用文献[32]的直线拟合数据，在该算例数据中模拟相关系数，其中坐标观测值与相关系数见表 1。

直线拟合模型为：

$ \left[{\begin{array}{*{20}{c}} {{y_1}}\\ {{y_2}}\\ \vdots \\ {{y_n}} \end{array}} \right] - \left[{\begin{array}{*{20}{c}} {{e_{{y_1}}}}\\ {{e_{{y_2}}}}\\ \vdots \\ {{e_{{y_n}}}} \end{array}} \right] = \left( {\left[{\begin{array}{*{20}{c}} {{x_1}}&1\\ {{x_2}}&1\\ \vdots & \vdots \\ {{x_n}}&1 \end{array}} \right] - \left[{\begin{array}{*{20}{c}} {{e_{{x_1}}}}&0\\ {{e_{{x_2}}}}&0\\ \vdots & \vdots \\ {{e_{{x_n}}}}&0 \end{array}} \right]} \right) \left[{\begin{array}{*{20}{c}} {{\beta _1}}\\ {{\beta _2}} \end{array}} \right] $

(23)

式中，[β₁, β₂]^T表示待估参数。

相关系数与坐标元素的协方差之间的关系为：

$ {\rho _{{x_i}{y_i}}} = \frac{{{\sigma _{{x_i}{y_i}}}}}{{\sqrt {\sigma _{{x_i}}^2\sigma _{{y_i}}^2} }}(i = 1, 2, \cdots, 10) $

(24)

由式(24)可知，当观测向量与系数矩阵的相关协因数阵的方差分量为σ₃²时，通过已知的相关系数和协方差可以得到互协因数Q_xiyi：

$ \sigma _3^2{\mathit{\boldsymbol{Q}}_{{x_i}{y_i}}} = {\rho _{{x_i}{y_i}}}{\sigma _1}\sqrt {{\mathit{\boldsymbol{Q}}_{{x_i}}}} {\sigma _2}\sqrt {{\mathit{\boldsymbol{Q}}_{{y_i}}}} $

(25)

继而得到：

$ {\mathit{\boldsymbol{Q}}_{{x_i}{y_i}}} = \frac{{{\rho _{{x_i}{y_i}}}{\sigma _1}\sqrt {{\mathit{\boldsymbol{Q}}_{{x_i}}}} {\sigma _2}\sqrt {{\mathit{\boldsymbol{Q}}_{{y_i}}}} }}{{\sigma _3^2}} $

(26)

因此，当相关系数ρ_xiyi和协因数阵Q_x、Q_y不变时，可通过改变3个方差分量确定互协因数阵Q_xy。

由上述分析可以得到相关观测时协因数阵Q_e的形式：

$ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] = \\\left[{\begin{array}{*{20}{c}} {\sigma _{{y_1}}^2/\sigma _1^2}& \cdots &0&{{\sigma _{{y_1}{x_1}}}/\sigma _3^2}& \cdots &0\\ \vdots & & \vdots & \vdots & & \vdots \\ 0& \cdots &{\sigma _{{y_n}}^2/\sigma _1^2}&0& \cdots &{{\sigma _{{y_n}{x_n}}}/\sigma _3^2}\\ {{\sigma _{{x_1}{y_1}}}/\sigma _3^2}& \cdots &0&{\sigma _{{x_1}}^2/\sigma _2^2}& \cdots &0\\ \vdots & & \vdots & \vdots & & \vdots \\ 0& \cdots &{{\sigma _{{x_n}{y_n}}}/\sigma _3^2}&0& \cdots &{\sigma _{{x_n}}^2/\sigma _2^2} \end{array}} \right] $

(27)

现分别使用最小二乘(LS)、不考虑相关性的总体最小二乘(TLS)、不考虑相关性的总体最小二乘方差分量估计(TLS-VCE)、相关观测PEIV模型总体最小二乘(PEIV-TLS)及相关观测PEIV模型的方差分量估计(PEIV-VCE)进行解算，$ \hat \sigma _1^2、\hat \sigma _2^2、\hat \sigma _3^2 $分别表示协因数阵Q_y、Q_a及其互协因数阵对应的方差分量估值，得到的结果见表 2。算例1方差分量估值变化图见图 2。

从表 2可以看出，考虑了方差分量估计的平差解算结果与不考虑方差分量估计的结果有所差异，当不考虑数据相关性时，本文给出的方法与已有方法结果相同；而考虑了相关性时，协因数阵Q_y与Q_a对应的方差分量估值与不考虑相关观测时接近，但此时可以得到互协因数阵对应的方差分量估值，从而修正参数估值。

3.2 算例2

假设一组数据从4~11.2每隔0.8模拟一个数据，并组成真值$ \mathit{\boldsymbol{\tilde x}} $，第二组数据从9.5~23.9每隔0.6模拟一个数据，并组成真值$ \mathit{\boldsymbol{\tilde y}} $，两组数据通过不同手段获取，并且考虑其相关性，其中模拟的权值为：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{x}}} = \\ {\rm{diag}}({[1000,300,600,100,100,800,400,300,600,100]^{\rm{T}}})\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{y}}} = \\ {\rm{diag}}({[100,180,200,100,400,400,200,300,200,200]^{\rm{T}}}) \end{array} $

其中，Q_x= P_x^－1；Q_y= P_y^－1。

构造直线拟合模型(式(23))，其中参数真值为2和1.5，考虑两组数据的相关性，相关系数由MATLAB函数rand随机生成，互协因数阵由相关系数生成，计算方式同式(27)。现给两组数据模拟均值为0、协方差阵为$ \left[{\begin{array}{*{20}{c}} {3{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}}&{2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}\\ {2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] $的随机误差，观测数据真值、观测值及相关系数见表 3。

分别使用LS、PEIV-TLS及PEIV-VCE模型解算得到的结果如表 4所示，其中，$ \left\| {\mathit{\boldsymbol{\hat \beta }} - \mathit{\boldsymbol{\tilde \beta }}} \right\| $表示参数估值与真值的差值范数。算例2方差分量估值变化图见图 3。

从表 4与图 3可以看出，方差分量估计方法可以得到更加合理的参数结果，且估计得到的方差分量估值与真值比较接近。

3.3 算例分析

1) 从算例1可以看出，将PEIV视为非线性最小二乘，并在近似值处展开化为间接平差形式进行迭代求解，可以得到与已有方法相同的参数估值结果，并在不考虑观测值间相关性时，使用方差分量估计方法可以得到与文献[14, 27, 32]相同的参数估值与方差分量估值，说明了方差分量估计方法可以对平差时的随机模型进行修正，从而得到在修正后随机模型下更加合理的参数结果。

2) 考虑观测数据的相关性，并针对随机模型的不准确进行相关观测PEIV模型的方差分量估计，在算例1中随机模拟相关系数生成观测值之间的互协因数，使用方差分量估计确定3个方差分量估值，从而修正参数估值。从表 2可以看出，协因数阵Q_y与Q_a对应的方差分量估值与未考虑相关性得到的估值比较接近，考虑了观测数据的相关性后，可以得到互协因数阵对应的方差分量估值，修正协因数阵，从而修正参数估值。

3) 从算例2模拟数据可以看出，考虑了方差分量估计得到的参数估值结果更接近于参数真值，得到的估值与真值间的差值范数更小，且得到的方差分量估值也接近于真值，进一步说明了方差分量估计方法的有效性。

4 结束语

随机模型的验后估计即方差分量估计方法是平差中重要的一部分，本文针对相关观测的PEIV模型并考虑平差时随机模型的不准确，将PEIV模型视为非线性最小二乘进行线性化得到间接平差形式进行迭代计算。对平差模型进行转换，根据最小二乘方差分量估计方法，推导了相关观测PEIV模型的最小二乘方差分量估计，将方差分量估值作为平差参数进行解算，通过方差分量估值修正协因数阵，从而得到更加合理的参数估值。对最小二乘方差分量估计公式转换可以得到与MINQUE及Helmert方差分量估计相同的表达式，进一步说明了方差分量估计方法的等价性，且本文是针对相关观测情形下，得到的公式更具有一般性。在算例实验中，通过改变协因数阵进行迭代计算，此时可能会出现不收敛或负方差的情况，且总体最小二乘实质上是有偏的，在迭代计算中，偏差是否会对解算结果有影响是今后要研究的工作。