经典最小二乘理论是以Gauss-Markov模型为基础，以观测向量的加权残差平方和最小为准则，获得统计意义上具有无偏性、有效性、一致性的最优参数估值^[1]，参数估计的形式简单，易于进行精度评定。在现代大地测量及相关数据处理领域(如坐标转换^[2]、直线和平面拟合^[3]、地震参数反演^[4])，设计矩阵元素也含有误差，由此建立了EIV模型及相应的整体最小二乘平差理论^[5]。Gauss-Markov模型中，解可以表示成观测量的显式表达，EIV模型只有当系数矩阵元素和观测向量等权不相关时，才能得到解的显式表达式。自Golub等^[6]利用奇异值分解技术推导EIV模型的整体最小二乘解以来，众多学者针对整体最小二乘模型算法、解的统计性质和模型可靠性展开大量研究。Mahboub^[7]采用Lagrange乘子法分别研究了系数矩阵元素和观测量不等权且两者互不相关情况下的整体最小二乘解。Amiri-Simkooei等^[8]在经典最小二乘理论框架下设计了形式简单的迭代算法，并给出解的近似方差。Fang^[9]将EIV模型扩展到系数矩阵元素和观测量不等权，在两者相关的情况下分别利用Lagrange乘子法、GHM方法和数值分析方法得到EIV模型的解。本文在此基础上，采用GHM方法推导一般权阵下EIV模型的近似方差矩阵，并将EIV模型重新用类似于Gauss-Markov模型的形式表达，在经典最小二乘理论下推导EIV模型的迭代算法和近似方差矩阵，其形式简便，易于算法的执行，且与已有算法等价。

1 经典最小二乘原理概述

在测量数据处理及相关领域，为获得待求参数的最优估值，通常建立观测值和待求参数之间的线性或线性化函数，即Gauss-Markov模型:

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{A\xi }} + \mathit{\boldsymbol{e}} $

(1)

式中，A为n×u设计矩阵，所有元素精确已知; y和e分别为n维观测向量和误差向量，观测值误差e ~N(0, D(y))，其中$D\left( {{\rm{ }}\mathit{\boldsymbol{y}}{\rm{ }}} \right) = \sigma _0^2\mathit{\boldsymbol{Q}}{_y} = \sigma _0^2\mathit{\boldsymbol{P}}_y^{ - 1} $。令观测值误差的加权平方和最小，可以获得参数的最优无偏估计为：

$ \mathit{\boldsymbol{\hat \xi }} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_y^{ - 1}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_y^{ - 1}\mathit{\boldsymbol{y}} $

(2)

观测值和误差向量的估计值(改正数)表示为：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\widehat y}} = {\mathit{\boldsymbol{P}}_A}\mathit{\boldsymbol{y}}}\\ {\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{\widehat e}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\widehat y}} = {\mathit{\boldsymbol{P}}_{\frac{1}{A}}}\mathit{\boldsymbol{y}}} \end{array}} \right. $

(3)

式中，$ \mathit{\boldsymbol{P}}{_A} = {\rm{ }}\mathit{\boldsymbol{A}}{\rm{ }}{(\mathit{\boldsymbol{A}}{^{\rm{T}}}\mathit{\boldsymbol{Q}}_y^{ - 1}\mathit{\boldsymbol{A}}{\rm{ }})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_y^{ - 1}$和${\rm{ }}{\mathit{\boldsymbol{P}}_{\frac{1}{A}}} = {\rm{ }}\mathit{\boldsymbol{I}}{\rm{ }} - {\rm{ }}\mathit{\boldsymbol{A}}{\rm{ }}{({\rm{ }}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_y^{ - 1}\mathit{\boldsymbol{A}}{\rm{ }})^{1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_y^{ - 1} $是两个n×n正交投影矩阵。式(3)将观测值的估值和误差改正数通过投影矩阵与观测向量建立起联系，参数估值、观测值估值和误差估值的协因数阵分别为$\mathit{\boldsymbol{ }}{\mathit{\boldsymbol{Q}}_{\hat x}} = {\left( {{\rm{ }}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_y^{ - 1}\mathit{\boldsymbol{A}}{\rm{ }}} \right)^{ - 1}} $、${\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat y}}}} = {\rm{ }}{\mathit{\boldsymbol{P}}_A}\mathit{\boldsymbol{Q}}{_y} $和${\rm{ }}{\mathit{\boldsymbol{Q}}_{\hat e}} = {\rm{ }}{\mathit{\boldsymbol{P}}_{\frac{1}{A}}}\mathit{\boldsymbol{Q}}{_y} $。单位权方差估值为：

$ \hat \sigma _0^2 = \frac{{{{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_y}\mathit{\boldsymbol{\hat e}}}}{{n - u}} $

(4)

2 WTLS方法及精度分析

当设计矩阵元素含误差时，EIV模型的函数表达式为:

$ \mathit{\boldsymbol{y}} = \left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_A}} \right)\mathit{\boldsymbol{\xi }} + {\mathit{\boldsymbol{e}}_y} $

(5)

相应的随机模型为：

$ \begin{array}{*{20}{c}} {D\left( \mathit{\boldsymbol{l}} \right) = D\left( {\begin{array}{*{20}{c}} {{\rm{vec}}\left( \mathit{\boldsymbol{A}} \right)}\\ \mathit{\boldsymbol{y}} \end{array}} \right) = D\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_A}}\\ {{\mathit{\boldsymbol{e}}_y}} \end{array}} \right) = D\left( \mathit{\boldsymbol{e}} \right) = }\\ {\sigma _0^2{\mathit{\boldsymbol{Q}}_{ll}} = \sigma _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{AA}}}&{{\mathit{\boldsymbol{Q}}_{Ay}}}\\ {{\mathit{\boldsymbol{Q}}_{yA}}}&{{\mathit{\boldsymbol{Q}}_{yy}}} \end{array}} \right] = \sigma _0^2{\mathit{\boldsymbol{P}}^{ - 1}}} \end{array} $

(6)

式中，vec(·)为矩阵按列进行向量化运算算子，l =(vec^T(A)y^T)^T，e =(e_A^T e_y^T)^T，Q_AA为nu×nu系数阵元素的协因数阵，Q_yy为n×n观测值协因数阵，Q_Ay= Q_yA^T为nu×n系数阵元素和观测值向量的互协因数阵，Q_ll为对称正定矩阵，σ₀²为单位权方差。根据整体最小二乘准则，组成拉格朗日目标函数：

$ \mathit{\Phi }(\mathit{\boldsymbol{\xi }}) = {\mathit{\boldsymbol{e}}^{\mathit{d }}}\mathit{\boldsymbol{Pe}} + 2{\mathit{\boldsymbol{\lambda }}^{\rm{T}}}\left( {\mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{e}}_y} - \mathit{\boldsymbol{A\xi }} + {\mathit{\boldsymbol{E}}_A}\mathit{\boldsymbol{\xi }}} \right) $

(7)

令${\mathit{\boldsymbol{\hat \xi }}} $、${\mathit{\boldsymbol{\hat \lambda }}} $分别为参数和Lagrange乘子的估值，$ {\mathit{\boldsymbol{\hat e}}}$为观测误差的预测值，且$\mathit{\boldsymbol{\hat B}}{\rm{ }} = [{{\mathit{\boldsymbol{\hat \xi }}}^{\rm{T}}} \otimes {\rm{ }}\mathit{\boldsymbol{I}}{_n}, {\rm{ }} - {\rm{ }}\mathit{\boldsymbol{I}}{_n}] $。根据文献[12]，将式(7)分别对ξ、e、λ求偏导，并令其值为0，可以得到如下关系式：

$ \boldsymbol{A}^{\mathrm{T}} \hat{\boldsymbol{\lambda}}-\boldsymbol{E}_{A}^{\mathrm{T}} \hat{\boldsymbol{\lambda}}=0 $

(8)

$ \hat{\boldsymbol{e}}=-\boldsymbol{Q}_{l l} \hat{\boldsymbol{B}}^{\mathrm{T}} \hat{\lambda} $

(9)

$ \mathit{\boldsymbol{\hat \lambda }} = {\left( {\mathit{\boldsymbol{\widehat B}}{\mathit{\boldsymbol{Q}}_{ll}}{{\mathit{\boldsymbol{\widehat B}}}^{\rm{T}}}} \right)^{ - 1}}(\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A\hat \xi }}) $

(10)

将式(10)代入式(9)得：

$ \mathit{\boldsymbol{\widehat e}} = - {\mathit{\boldsymbol{Q}}_{ll}}{\mathit{\boldsymbol{\widehat B}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{\widehat B}}{\mathit{\boldsymbol{Q}}_{ll}}{{\mathit{\boldsymbol{\widehat B}}}^{\rm{T}}}} \right)^{ - 1}}(\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}\mathit{\boldsymbol{\widehat \xi }}) $

(11)

根据式(6)中Q_ll的定义，将$ {{\mathit{\boldsymbol{\hat e}}}_A}$从${\mathit{\boldsymbol{\hat e}}} $中提取出来：

$ {\mathit{\boldsymbol{\hat e}}_A} = - \left[ {{\mathit{\boldsymbol{Q}}_{AA}}{\mathit{\boldsymbol{Q}}_{Ay}}} \right]\left( {{{\mathit{\boldsymbol{\hat B}}}^{\rm{T}}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{ll}}{{\mathit{\boldsymbol{\hat B}}}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A\hat \xi }}} \right)} \right) $

(12)

对式中$ {{\mathit{\boldsymbol{\hat e}}}_A}$进行逆向量化计算，得到矩阵形式的$ {{\mathit{\boldsymbol{\hat e}}}_A}$=ivec($ {{\mathit{\boldsymbol{\hat e}}}_A}$)，其中算子ivec(·)是vec(·)的逆运算。将式(12)逆向量化计算的结果代入式(8)，并结合式(10)，可以得到：

$ \begin{array}{*{20}{c}} {\left( {{{\left( {\mathit{\boldsymbol{A}} - {{\mathit{\boldsymbol{\widehat e}}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{\widehat B}}{\mathit{\boldsymbol{Q}}_{ll}}{{\mathit{\boldsymbol{\widehat B}}}^{\rm{T}}}} \right)}^{ - 1}}} \right)\mathit{\boldsymbol{A}}\mathit{\boldsymbol{\widehat \xi }} = }\\ {{{\left( {\mathit{\boldsymbol{A}} - {{\mathit{\boldsymbol{\widehat e}}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{\widehat B}}{\mathit{\boldsymbol{Q}}_{ll}}{{\mathit{\boldsymbol{\widehat B}}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{y}}} \end{array} $

(13)

经过简单变换，可得参数${\mathit{\boldsymbol{\hat \xi }}} $的表达式：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat \xi }} = {{\left( {{{\left( {\mathit{\boldsymbol{A}} - {{\mathit{\boldsymbol{\hat e}}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{\hat B}}{\mathit{\boldsymbol{Q}}_{ll}}{{\mathit{\boldsymbol{\hat B}}}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{A}} - {{\mathit{\boldsymbol{\hat e}}}_A}} \right)} \right)}^{ - 1}} \cdot }\\ {{{\left( {\mathit{\boldsymbol{A}} - {{\mathit{\boldsymbol{\hat e}}}_A}} \right)}^{\rm{T}}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{ll}}{{\mathit{\boldsymbol{\hat B}}}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{y}} - {{\mathit{\boldsymbol{\hat e}}}_A}\mathit{\boldsymbol{\hat \xi }}} \right)} \end{array} $

(14)

由于式(14)的右端$ {{\mathit{\boldsymbol{\hat e}}}_A}$和$ {\mathit{\boldsymbol{\hat B}}}$都是$ {\mathit{\boldsymbol{\hat \xi }}}$的函数，所以需要选定恰当的初始值${{\mathit{\boldsymbol{\hat \xi }}}^0} $进行迭代计算，具体迭代步骤见文献[12]中算法2, 它属于Gauss-Newton算法，具有一阶线性收敛性。为得到解的精度，将EIV模型视为GHM模型，线性化后推导解的一阶近似方差矩阵。设

$ f\left( {\mathit{\boldsymbol{l}} + \mathit{\boldsymbol{e}},\mathit{\boldsymbol{\xi }}} \right) = \left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_A}} \right)\mathit{\boldsymbol{\xi }} - \mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{e}}_y} $

(15)

将式(15)分别在观测值和近似解处线性化得：

$ \begin{array}{*{20}{c}} {f\left( {\mathit{\boldsymbol{l}} + \mathit{\boldsymbol{e}},\mathit{\boldsymbol{\xi }}} \right) = \frac{{\partial f}}{{\partial {\mathit{\boldsymbol{\xi }}^{\rm{T}}}}}\left| {_{l + {\mathit{\boldsymbol{e}}^{\left( i \right)}},{\xi ^{\left( i \right)}}}} \right.{\rm{d}}{\xi ^{i + 1}} + \frac{{\partial f}}{{\partial {\mathit{\boldsymbol{l}}^{\rm{T}}}}}\left| {_{l + {e^{\left( i \right)}},{\xi ^{\left( i \right)}}}} \right. \cdot }\\ {{{\left( {\mathit{\boldsymbol{e}} - \mathit{\boldsymbol{e}}} \right)}^{\left( i \right)}} + f\left( {\mathit{\boldsymbol{l}} + {\mathit{\boldsymbol{e}}^{\left( i \right)}},{\mathit{\boldsymbol{\xi }}^{\left( i \right)}}} \right) = }\\ {\frac{{\partial f}}{{\partial {\mathit{\boldsymbol{\xi }}^{\rm{T}}}}}\left| {_{l + {e^{\left( i \right)}},{\xi ^{\left( i \right)}}}} \right.{\rm{d}}{\xi ^{i + 1}} + \frac{{\partial f}}{{\partial {\mathit{\boldsymbol{l}}^{\rm{T}}}}}\left| {_{l + {e^{\left( i \right)}},{\xi ^{\left( i \right)}}}} \right. \cdot }\\ {\mathit{\boldsymbol{e}} + f\left( {\mathit{\boldsymbol{l}},{\mathit{\boldsymbol{\xi }}^{\left( i \right)}}} \right) = \left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}}_A^{\left( i \right)}} \right){\rm{d}}{\mathit{\boldsymbol{\xi }}^{\left( {i + 1} \right)}} + }\\ {\left[ {\begin{array}{*{20}{c}} {{{\left( {{\mathit{\boldsymbol{\xi }}^{\left( i \right)}}} \right)}^{\rm{T}}} \otimes {\mathit{\boldsymbol{l}}_n}}&{ - {\mathit{\boldsymbol{l}}_n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{e}}_A^{\left( {i + 1} \right)}}\\ {\mathit{\boldsymbol{e}}_y^{\left( {i + 1} \right)}} \end{array}} \right] + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\xi }}^{\left( i \right)}} - \mathit{\boldsymbol{y}}} \end{array} $

(16)

式(16)在经典平差中是附有参数的条件平差模型，其解为^[1]:

$ \begin{array}{l} {\rm{d}}{{\mathit{\boldsymbol{\hat \xi }}}^{\left( {i + 1} \right)}} = \\ {\left( {{{\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)}^{\rm{T}}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}{\mathit{\boldsymbol{Q}}_{ll}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}} \right)}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)} \right)^{ - 1}} \cdot \\ {\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)^{\rm{T}}}{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}{\mathit{\boldsymbol{Q}}_{ll}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}} \right)}^{\rm{T}}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat \xi }}}^{\left( i \right)}}} \right) \end{array} $

(17)

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat \xi }}}^{\left( {i + 1} \right)}} = {{\mathit{\boldsymbol{\hat \xi }}}^{\left( i \right)}} + {\rm{d}}{{\mathit{\boldsymbol{\hat \xi }}}^{\left( {i + 1} \right)}} = }\\ {{{\left( {{{\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)}^{\rm{T}}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}{\mathit{\boldsymbol{Q}}_{ll}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}} \right)}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)} \right)}^{ - 1}} \cdot }\\ {{{\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)}^{\rm{T}}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}{\mathit{\boldsymbol{Q}}_{ll}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}} \right)}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}{{\mathit{\boldsymbol{\hat \xi }}}^{\left( i \right)}}} \right)} \end{array} $

(18)

显然，由GHM方法推导的式(18)和式(14)等价。由于式(16)为线性化模型，在近似值处可以求得参数近似协因数矩阵及单位权方差:

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{\hat \xi \hat \xi }} = \left( {{{\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)}^{\rm{T}}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}{\mathit{\boldsymbol{Q}}_{ll}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}} \right)}^{\rm{T}}}} \right)}^{ - 1}}} \right. \cdot }\\ {{{\left. {\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_A^{\left( {i + 1} \right)}} \right)} \right)}^{ - 1}}} \end{array} $

(19)

$ \begin{array}{*{20}{c}} {\sigma _0^2 = \frac{{{{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{P\hat e}}}}{{n - u}} = }\\ {\frac{{{{\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat \xi }}}^{\left( {i + 1} \right)}}} \right)}^{\rm{T}}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}{\mathit{\boldsymbol{Q}}_{ll}}{{\left( {{{\mathit{\boldsymbol{\hat B}}}^{\left( i \right)}}} \right)}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat \xi }}}^{\left( {i + 1} \right)}}} \right)}}{{n - u}}} \end{array} $

(20)

由于EIV模型固有的非线性特征，式(19)、式(20)为有偏估计。以上参数估计基于Lagrange乘子法，而精度评定是将EIV模型当作非线性GHM模型进行一阶线性化后，根据方差协方差传播律求取。下面根据经典最小二乘理论研究EIV模型的参数估计和精度评定问题。

3 基于经典LS原理的WTLS方法

对于式(5)表示的EIV模型，两边分别减去E_Aξ，得到模型(5)的变式：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{E}}_A}\mathit{\boldsymbol{\xi }} = \left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_A}} \right)\mathit{\boldsymbol{\xi }} + \left( {{\mathit{\boldsymbol{e}}_y} - {\mathit{\boldsymbol{E}}_A}\mathit{\boldsymbol{\xi }}} \right) = }\\ {\left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_A}} \right)\mathit{\boldsymbol{\xi }} - \mathit{\boldsymbol{Be}}} \end{array} $

(21)

令$\mathit{\boldsymbol{\tilde y}} = {\rm{ }}\mathit{\boldsymbol{y}}{\rm{ }} - \mathit{\boldsymbol{E}}{_A}\mathit{\boldsymbol{\xi }}, \mathit{\boldsymbol{\tilde A}} = {\rm{ }}\mathit{\boldsymbol{A}}{\rm{ }} - {\rm{ }}\mathit{\boldsymbol{E}}{_A}, {\rm{ }}\mathit{\boldsymbol{\tilde e}} = - {\rm{ }}\mathit{\boldsymbol{Be}} $，则式(21)可写成类似式(1)的形式：

$ \mathit{\boldsymbol{\tilde y}} = \mathit{\boldsymbol{\tilde A\xi }} + \mathit{\boldsymbol{\tilde e}} $

(22)

由方差协方差传播律可得：

$ {\mathit{\boldsymbol{Q}}_{\tilde y}} = {\mathit{\boldsymbol{Q}}_{\tilde e}} = \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{ll}}{\mathit{\boldsymbol{B}}^{\rm{T}}} $

(23)

对式(22)、式(23)直接采用最小二乘准则可得：

$ \mathit{\boldsymbol{\hat \xi }} = {\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{\tilde A}}} \right)^{ - 1}}{\mathit{\boldsymbol{\tilde A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{\tilde y}} $

(24)

式(24)和式(14)等价，且与式(2)在形式上一致。参数相应的协因数矩阵为：

$ {\mathit{\boldsymbol{Q}}_{\hat \xi \hat \xi }} = {\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{\tilde A}}} \right)^{ - 1}} $

(25)

观测值和残差的估值为：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\hat y}} = {\mathit{\boldsymbol{P}}_{\tilde A}}\mathit{\boldsymbol{\tilde y}} = \mathit{\boldsymbol{\tilde A\hat \xi }} = \left( {\mathit{\boldsymbol{A}} - {{\mathit{\boldsymbol{\tilde e}}}_A}} \right)\mathit{\boldsymbol{\hat \xi }}\\ \mathit{\boldsymbol{\hat e}} = P_{\tilde A}^ \bot \mathit{\boldsymbol{\tilde y}} = \mathit{\boldsymbol{\tilde y}} - \mathit{\boldsymbol{\hat y}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A\hat \xi }} = {\mathit{\boldsymbol{Q}}_{\tilde y}}\mathit{\boldsymbol{\hat \lambda }} \end{array} \right. $

(26)

式中，$ {\mathit{\boldsymbol{P}}_{\tilde A}} = \mathit{\boldsymbol{\tilde A}}{\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{\tilde A}}} \right)^{ - 1}}{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}, \mathit{\boldsymbol{P}}_{\tilde A}^ \bot = \mathit{\boldsymbol{I}} - \mathit{\boldsymbol{\tilde A}} \cdot {\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{\tilde A}}} \right)^{ - 1}}{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}$。需要指出的是，$\mathit{\boldsymbol{\hat e}} \ne {{\mathit{\boldsymbol{\hat e}}}_y} $。对式(26)运用方差传播律，可以得到观测值和残差估值的协因数阵分别为：

$ {\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat y}}}} = {\mathit{\boldsymbol{P}}_{\tilde A}}{\mathit{\boldsymbol{Q}}_{\tilde y}} = \mathit{\boldsymbol{\tilde A}}{\mathit{\boldsymbol{Q}}_{\hat \xi \hat \xi }}{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}} $

(27)

$ {\mathit{\boldsymbol{Q}}_v} = {\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat e}}}} = \mathit{\boldsymbol{P}}_{\tilde A}^ \bot {\mathit{\boldsymbol{Q}}_{\tilde y}} = {\mathit{\boldsymbol{Q}}_{\tilde y}} - {\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat y}}}} $

(28)

单位权方差为：

$ \sigma _0^2 = \frac{{{{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{\hat e}}}}{{n - u}} = \frac{{{{\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A\hat \xi }}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A\hat \xi }}} \right)}}{{n - u}} $

(29)

忽略上标可以看出，式(25)、式(29)和式(19)、式(20)等价。基于经典最小二乘原理加权EIV模型的迭代步骤为：

1) 给定系数阵A、观测向量y及协因数矩阵Q_ll，计算初始值${{\mathit{\boldsymbol{\hat \xi }}}^0} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1}\mathit{\boldsymbol{y}} $;

2) 令i=0, 1, …k，k为迭代次数，依次计算${{\mathit{\boldsymbol{\hat B}}}^{(i + 1)}} = [{({{\mathit{\boldsymbol{\hat \xi }}}^{\left( i \right)}})^{\rm{T}}} \otimes {\rm{ }}\mathit{\boldsymbol{I}}{_n} - {\rm{ }}\mathit{\boldsymbol{I}}{_n}], \mathit{\boldsymbol{\hat \lambda }}{^{(i + 1)}} = {(\mathit{\boldsymbol{\hat B}}{^{(i + 1)}}\mathit{\boldsymbol{Q}}{_{ll}}\cdot{(\mathit{\boldsymbol{\hat B}}{^{(i + 1)}})^{\rm{T}}})^{ - 1}}\left( {{\rm{ }}\mathit{\boldsymbol{y}} - {\rm{ }}\mathit{\boldsymbol{A\hat \xi }}{^{(i)}}} \right), $ $ \mathit{\boldsymbol{\hat e}}_A^{\left( {i + 1} \right)} = - \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}{_{AA}}}&{\mathit{\boldsymbol{Q}}{_{Ay}}} \end{array}} \right]\cdot{({{\mathit{\boldsymbol{\hat B}}}^{(i + 1)}})^{\rm{T}}}{{\mathit{\boldsymbol{\hat \lambda }}}^{(i + 1)}}, \mathit{\boldsymbol{\hat e}}_A^{\left( {i + 1} \right)} = {\rm{ivec}}\left( {\mathit{\boldsymbol{\hat e}}_A^{\left( {i + 1} \right)}} \right), {{\mathit{\boldsymbol{\tilde A}}}^{(i + 1)}} = {\rm{ }}\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat e}}_A^{\left( {i + 1} \right)}, $ ${({\rm{ }}\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1})^{(i + 1)}} = {({\rm{ }}\mathit{\boldsymbol{B}}{^{(i + 1)}}\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{}}_{ll}}{({\rm{ }}\mathit{\boldsymbol{B}}{^{(i + 1)}})^{\rm{T}}})^{ - 1}}, {{\mathit{\boldsymbol{\tilde y}}}^{(i + 1)}} = {\rm{ }}\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\hat e}}_A^{\left( {i + 1} \right)}\mathit{\boldsymbol{\xi }}{^{(i)}}, \mathit{\boldsymbol{Q}}_{\hat \xi \hat \xi }^{(i + 1)} = {({({{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}})^{(i + 1)}}{(\mathit{\boldsymbol{Q}}_{\tilde y}^{ - 1})^{(i + 1)}}{{\mathit{\boldsymbol{\tilde A}}}^{(i + 1)}})^{1}}, $ $ \mathit{\boldsymbol{\hat \xi }}{^{(i + 1)}} = {\rm{ }}\mathit{\boldsymbol{Q}}_{\hat \xi \hat \xi }^{(i + 1)}{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}^{\left( {i + 1} \right)}{\rm{ }}{(Q_{\tilde y}^{ - 1})^{(i + 1)}};$

3) 给定一个足够小的正常数ε，若$\parallel {{\mathit{\boldsymbol{\hat \xi }}}^{\left( {i + 1} \right)}} - {{\mathit{\boldsymbol{\hat \xi }}}^{\left( i \right)}}\parallel < \varepsilon $，停止迭代，$ \mathit{\boldsymbol{\hat \xi }} = {{\mathit{\boldsymbol{\hat \xi }}}^{\left( {i + 1} \right)}}$。

上述迭代步骤与基于Lagrange乘子法的迭代方法实质上是相同的。

4 算例分析

采用文献[12]中的直线拟合算例，10个点的观测值及相应的权列于表 1，相应的EIV模型见文献[12]，根据表 1中数据估计斜率和截距。

按照本文算法，以加权最小二乘解为初值，令ε=10^－10，经7次迭代得到拟合的2个参数估计值ξ =(5.479 910 224 －0.480 533 407)^T。根据式(25)得到解的协因数矩阵为：

$ {\mathit{\boldsymbol{Q}}_{\hat \xi \hat \xi }} = {10^{ - 2}} \times \left( {\begin{array}{*{20}{c}} {8.700\;8}&{ - 1.6473}\\ { - 1.6473}&{0.3362} \end{array}} \right) $

根据式(29)得到单位权方差估值σ₀²=1.483 3。则参数估值的一阶近似方差矩阵为：

$ {\mathit{\boldsymbol{D}}_{\hat \xi \hat \xi }} = \sigma _0^2{\mathit{\boldsymbol{Q}}_{\hat \xi \hat \xi }} = \left( {\begin{array}{*{20}{c}} {0.129\;1}&{ - 0.024\;4}\\ { - 0.024\;4}&{0.005\;0} \end{array}} \right) $

另外，观测量和残差的估值及相应的协因数矩阵可以由式(26)~(28)求得。可见，加权EIV模型可以用经典最小二乘理论进行迭代求解，相应的单位权方差、参数估值、观测值和残差的估值及协因数矩阵也可以用经典最小二乘法导出，且表达形式与Gauss-Markov模型一致。

5 结语

本文推导了权矩阵是一般对称正定矩阵EIV模型的近似单位权方差及一阶协因数矩阵。将EIV模型当作包含观测值向量和模型参数的GHM模型并线性化，以最优解作为线性化的近似值，根据方差协方差传播定律获得解的精度信息。将EIV模型用类似于Gauss-Markov模型的形式表达，以经典最小二乘理论为基础，推导出EIV模型的迭代算法和近似方差矩阵，与已有算法得到的参数估计和GHM模型推导的方差是等价的。同时，建立了EIV模型的精度评价体系。