在实际的测量工作中，观测数据中除含有随机误差外，在一定程度上还含有粗差与系统误差等。粗差的存在使得以随机误差作为假设前提的平差模型无法获得最优的参数估值。为降低粗差对平差结果的影响，常采用选权迭代^[1-2]的方法降低粗差在平差中的权重，进而获得稳健的参数估计结果。

在测绘领域，三维基准转换模型^[3-4]、边长变化反演应变参数模型^[5-6]以及遥感影像混合像元解混模型^[7]等都可以看作是变量含误差(errors-in-variables, EIV)模型，因为其观测向量与系数矩阵均由含有误差的变量组成。观测值的增加提高了模型存在粗差的可能性，因此将总体最小二乘(total least squares, TLS)方法和稳健估计方法引入到EIV模型的解算中。当前，TLS理论已有了长足的发展，包括顾及系数矩阵结构性的混合TLS方法^[8]、加权TLS方法^[9-10]、顾及先验约束信息的约束TLS方法^[11]和顾及观测值粗差的稳健总体最小二乘方法(robust total least squares, RTLS)^[12-17]等。RTLS方法主要包括基于常用权函数的RTLS^[12-16]和基于中位数法的RTLS^[17]。基于常用权函数的RTLS方法的讨论热点是稳健估计策略的选择，即是否对系数矩阵和观测向量采用同一个单位权中误差进行降权处理。文献[12-14]采用同一个单位权中误差分别进行Huber与IGGⅢ权函数的降权处理；文献[15-16]对系数矩阵和观测向量分别采用各自的单位权中误差进行Huber与IGGⅢ权函数的降权处理。

采用RTLS^[12-16]方法对含有多列独立变量(以下简称为多变量)的EIV模型进行平差计算时，多变量EIV模型的各列变量的改正数受其对应的参数估值与观测向量先验精度的联合影响。两者的乘积越大，则该列改正数越大。因此，使用现有RTLS^[12-16]方法处理多变量EIV模型时，会优先对多变量EIV模型的某一列向量进行降权处理，造成平差结果不合理甚至错误，本文称该现象为虚假稳健估计现象。鉴于此，提出多变量稳健总体最小二乘(RTLS for multivariable EIV model, MRTLS)平差方法。该方法将多变量稳健估计策略与RTLS方法相结合，解决含有粗差的多变量EIV模型的虚假稳健估计问题。仿真算例验证了本文方法的正确性和有效性。

1 多变量稳健总体最小二乘估计 1.1 MRTLS方法

已有文献报道的RTLS方法存在多种不同的函数模型^[12-15]。本文结合特征矩阵^[18]的思路，将其函数模型统一为概括总体最小二乘函数模型：

${\mathit{\boldsymbol{L}}_{\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}} = {\mathit{\boldsymbol{A}}_{\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}}\mathit{\boldsymbol{X}}$

(1)

式中，${\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} \in {\mathit{\boldsymbol{R}}^{k \times 1}}$为EIV模型的独立观测值；$\mathit{\boldsymbol{e}} \in {{\bf{R}}^{k \times 1}}$为${\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}}$的随机误差；${\mathit{\boldsymbol{A}}_{\left({{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}} \in {{\bf{R}}^{m \times n}}$表示系数矩阵是$\left({{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)$的函数；${\mathit{\boldsymbol{L}}_{\left({{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}} \in {{\bf{R}}^{m \times 1}}$表示观测向量为$\left({{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)$的函数；$\mathit{\boldsymbol{X}} \in {{\bf{R}}^{n \times 1}}$为待估参数。

结合特征矩阵^[18]的思路，式(1)经线性化后得：

$\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \left( {{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{e}}$

(2)

式中，$\mathit{\boldsymbol{A}} \in {{\bf{R}}^{m \times n}}$为系数矩阵；$\mathit{\boldsymbol{L}} \in {{\bf{R}}^{m \times 1}}$为观测向量；$\mathit{\boldsymbol{B}} = {\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_m}$，$ \otimes $表示Kronecker直积；${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}$是根据系数矩阵观测量与EIV模型独立观测量之间的线性关系构造的特征函数矩阵；${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}}$是根据观测向量与EIV模型独立观测量之间的线性关系构造的特征函数矩阵。

特征矩阵C_A、C_L的计算公式为：

${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}} = {\left. {\frac{{\partial {\rm{vec}}\left( {{\mathit{\boldsymbol{A}}_{\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}}} \right)}}{{\partial \left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}}} \right|_{{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}}}},{\rm{\;\;\;\;}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} = {\left. {\frac{{\partial {\mathit{\boldsymbol{L}}_{\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}}}}{{\partial \left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}} - \mathit{\boldsymbol{e}}} \right)}}} \right|_{{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{e}}}}}$

(3)

式中，vec(·)表示矩阵的拉直算子。

根据式(3)可以建立独立观测量误差与观测向量误差和系数矩阵误差的关系式：

${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{L}}} = {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}}\mathit{\boldsymbol{e}},{\rm{\;}}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{A}}} = {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{e}}$

(4)

式中，e_L为观测向量的误差向量；e_A为系数矩阵的误差向量。

EIV模型的观测量个数相较于经典平差模型有所增加，提高了模型存在粗差的可能性。因此在稳健估计时，需采用相应的稳健估计权函数。根据拉格朗日条件极值法原理，建立多变量稳健总体最小二乘平差准则函数：

$\psi = {\mathit{\boldsymbol{e}}^{\rm{T}}}\bar Q_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{e}} - 2{\mathit{\boldsymbol{\eta }}^{\rm{T}}}\left\{ {\mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}} + \left( {{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{e}}} \right\} $

(5)

式中，${\mathit{\boldsymbol{\eta }}_{m{\rm{\; }} \times 1}}$为联系数向量；$\mathit{\boldsymbol{\bar Q}}_\mathit{\boldsymbol{e}}^{}$为等价协因数阵；$\mathit{\boldsymbol{\bar P}}_\mathit{\boldsymbol{e}}^{}$为等价权矩阵，且$\mathit{\boldsymbol{\bar Q}}_\mathit{\boldsymbol{e}}^{ - 1} = \mathit{\boldsymbol{\bar P}}_\mathit{\boldsymbol{e}}^{}$。通过多变量稳健估计权函数计算降权矩阵$\mathit{\boldsymbol{W}}_\mathit{\boldsymbol{e}}^{}$，进而根据$\mathit{\boldsymbol{\bar P}}_\mathit{\boldsymbol{e}}^{} = \mathit{\boldsymbol{W}}_\mathit{\boldsymbol{e}}^{}\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{e}}^{}$计算等价权矩阵($\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{e}}^{}$为初始的观测值权阵)。

对准则函数分别求关于${\mathit{\boldsymbol{e}}^{\rm{T}}}$、${\mathit{\boldsymbol{X}}^{\rm{T}}}$、$\mathit{\boldsymbol{\eta }}$的偏导，可得：

$\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial \psi }}{{\partial {\mathit{\boldsymbol{e}}^{\rm{T}}}}} = \mathit{\boldsymbol{e}}_{}^{\rm{T}}\mathit{\boldsymbol{\bar Q}}_e^{ - 1} - {\mathit{\boldsymbol{\eta }}^{\rm{T}}}\left( {{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}} \right)}\\ {\frac{{\partial \psi }}{{\partial {\mathit{\boldsymbol{X}}^{\rm{T}}}}} = - {\mathit{\boldsymbol{\eta }}^{\rm{T}}}\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}^{}} \right)}\\ {\frac{{\partial \psi }}{{\partial \mathit{\boldsymbol{\eta }}}} = \mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}} + \left( {{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{e}}} \end{array}} \right.$

(6)

约定上标“^”表示变量估值，则整理可得参数估计式为：

$ \mathit{\boldsymbol{\hat X}} = {\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}{{\mathit{\boldsymbol{\bar M}}}^{ - 1}}\mathit{\boldsymbol{\tilde A}}} \right)^{ - 1}}{\mathit{\boldsymbol{\tilde A}}^{\rm{T}}}{\mathit{\boldsymbol{\bar M}}^{ - 1}}\mathit{\boldsymbol{\tilde L}} $

(7)

式中，$\mathit{\boldsymbol{\tilde A}} = \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{\hat E}}_\mathit{\boldsymbol{A}}^{}{\rm{\; }}$；$\mathit{\boldsymbol{\tilde L}} = \mathit{\boldsymbol{L}} - \mathit{\boldsymbol{\hat E}}_\mathit{\boldsymbol{A}}^{}\mathit{\boldsymbol{\hat X}}$；$\mathit{\boldsymbol{\bar M}} = \left({{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} - \mathit{\boldsymbol{\hat B}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{\bar Q}}_\mathit{\boldsymbol{e}}^{}{\left({{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} - \mathit{\boldsymbol{\hat B}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}} \right)^{\rm{T}}}$；$\mathit{\boldsymbol{\hat E}}_\mathit{\boldsymbol{A}}^{} = {\rm{ivec}}\left({{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{A}}}} \right)$表示重构后的系数矩阵改正数矩阵，ivec(·)表示拉直算子的逆变换。

改正数估计式为：

$\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\hat e}} = \mathit{\boldsymbol{\bar Q}}_\mathit{\boldsymbol{e}}^{}{{\left( {{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}} - \mathit{\boldsymbol{\hat B}}{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}} \right)}^{\rm{T}}}{{\mathit{\boldsymbol{\bar M}}}^{ - 1}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A\hat X}}} \right)}\\ {{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{L}}} = {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}}\mathit{\boldsymbol{\hat e}}}\\ {{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{A}}} = {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{\hat e}}} \end{array}} \right. $

(8)

单位权中误差估计式为：

$ \hat \sigma _0^2 = \frac{{{{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{\bar Q}}_\mathit{\boldsymbol{e}}^{ - 1}\mathit{\boldsymbol{\hat e}}}}{{m - n}} $

(9)

参数方差-协方差阵估计式为：

$ {\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{\hat X}}}} = \hat \sigma _0^2 \cdot {\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}{{\mathit{\boldsymbol{\bar M}}}^{ - 1}}\mathit{\boldsymbol{\tilde A}}} \right)^{ - 1}} $

(10)

从式(7)和式(10)可看出，本文导出的MRTLS方法的参数估计式与方差-协方差及TLS方法的参数估计式有着相似的形式。需要指出的是，当稳健估计策略采用本文的多变量稳健估计策略时，由此导出的方法称之为MRTLS方法；当采用文献[12-16]的稳健估计策略时，导出的方法称为RTLS方法；当不采用稳健估计方法时，导出的方法称为TLS方法。

多变量稳健总体最小二乘迭代计算步骤如下：

1) 输入迭代初值${\mathit{\boldsymbol{\hat X}}^{\left(0 \right)}}$、A、L、$\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^{}$；

2) 根据式(7)、式(8)进行总体最小二乘平差计算直至收敛，获得${\mathit{\boldsymbol{\hat X}}^{\left(i \right)}}$、${\mathit{\boldsymbol{\hat e}}^{\left(i \right)}}$；

3) 根据${\mathit{\boldsymbol{\hat e}}^{\left(i \right)}}$进行稳健估计，获得$\mathit{\boldsymbol{\bar P}}_\mathit{\boldsymbol{e}}^{\left(i \right)}$；

4) 重复步骤2)、3)，直至迭代收敛($\left| {{{\mathit{\boldsymbol{\hat X}}}^{\left({i + 1} \right)}} - {{\mathit{\boldsymbol{\hat X}}}^{\left(i \right)}}} \right| < 1 \times {10^{ - 10}}$)，输出${\mathit{\boldsymbol{\hat X}}^{\left({i + 1} \right)}}$、${\mathit{\boldsymbol{\hat e}}^{\left({i + 1} \right)}}$、$\hat \sigma _0^2$、${\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{\hat X}}}}$。

MRTLS方法的迭代计算步骤嵌套两个循环：内层循环是步骤2)的总体最小二乘循环；外层循环是步骤4)的稳健总体最小二乘循环。

1.2 现有RTLS稳健估计策略问题分析

本节依据§1.1导出的式(7)~式(10)分析指出文献[12-16]中的稳健估计策略在应用到多变量EIV模型时存在的虚假稳健估计问题。以Huber为例，RTLS方法的权函数^[16]为：

$ {w_{i,L}} = \left\{ {\begin{array}{*{20}{c}} {1,}&{\left| {{\mathit{\boldsymbol{e}}_{i,L}}} \right| \le 2{\sigma _L}}\\ {2{\sigma _L}/\left| {{\mathit{\boldsymbol{e}}_{i,L}}} \right|,}&{\left| {{\mathit{\boldsymbol{e}}_{i,L}}} \right| > 2{\sigma _L}} \end{array}} \right. $

(11)

$ {w_{j,A}} = \left\{ {\begin{array}{*{20}{c}} {1,}&{\left| {{\mathit{\boldsymbol{e}}_{j,A}}} \right| \le 2{\sigma _A}}\\ {2{\sigma _A}/\left| {{\mathit{\boldsymbol{e}}_{j,A}}} \right|,}&{\left| {{\mathit{\boldsymbol{e}}_{j,A}}} \right| > 2{\sigma _A}} \end{array}} \right. $

(12)

式中，w_i，L、w_j，A分别为观测向量和系数矩阵的权因子，下标i表示观测向量的第i个观测值，下标j表示系数阵的第j个观测值；e_i，L表示其改正数；σ_L表示其单位权中误差；e_j，A表示系数阵的改正数向量；σ_A表示系数阵的单位权中误差。

分析式(11)、式(12)可以发现，现有稳健估计策略对系数矩阵和观测向量分别进行降权处理，但对系数矩阵的观测量进行降权时采用同一个σ_A，从而忽视了系数矩阵中不同列观测量的改正数之间的比例关系，导致不合理的稳健估计结果。这里以多元线性回归模型为例进行具体的阐述。有线性回归模型：

${\mathit{\boldsymbol{L}}_0} = {X_1} \cdot {\mathit{\boldsymbol{L}}_1} + {X_2} \cdot {\mathit{\boldsymbol{L}}_2} + \cdots + {X_n} \cdot {\mathit{\boldsymbol{L}}_n}$

(13)

式中，L₀为观测向量；${\mathit{\boldsymbol{L}}_i}\left({i \in \left\{ {1 \cdots n} \right\}} \right)$为组成系数矩阵的变量; ${X_i}\left({i \in \left\{ {1 \cdots n} \right\}} \right)$为模型待估参数。本文假设${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_n}$为独立观测获得的不等精度数据，即非等精度的独立观测值，${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_n}$对应的先验协因数阵分别为$\mathit{\boldsymbol{Q}}_i^{\rm{*}}{\rm{\; }}, {\rm{\; }}i \in \left\{ {0, 1 \cdots n} \right\}$。

显然，式(13)所示的多元线性回归模型是EIV模型，因此可使用TLS方法进行平差计算。根据上述信息按照式(8)计算模型的改正数：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat e}}}^ * } = \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^ * {{\left( {C_L^ * - {{\hat B}^ * }C_A^ * } \right)}^{\rm{T}}}{M^{ * - 1}}\left( {{L_0} - {A^ * }{{\hat X}^ * }} \right) = }\\ { - \left[ {\begin{array}{*{20}{c}} {{S_0}}\\ {{S_1}}\\ \vdots \\ {{S_n}} \end{array}} \right]\left( {{L_0} - {A^ * }{{\hat X}^ * }} \right)}\\ {式中,{{\hat B}^ * } = {{\hat X}^{ * {\rm{T}}}} \otimes I;C_A^ * = \left[ {\begin{array}{*{20}{c}} 0&I \end{array}} \right];}\\ {{M^ * } = \left( {C_L^ * - {{\hat B}^ * }C_A^ * } \right)\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^ * {{\left( {C_L^ * - {{\hat B}^ * }C_A^ * } \right)}^{\rm{T}}};}\\ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}^ * = I \otimes \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}_0^ * }& \cdots &{\mathit{\boldsymbol{Q}}_n^ * } \end{array}} \right];C_L^ * = \left[ {\begin{array}{*{20}{c}} I& \cdots &0 \end{array}} \right];}\\ {{\mathit{\boldsymbol{A}}^ * } = \left[ {\begin{array}{*{20}{c}} {{L_1}}& \cdots &{{L_n}} \end{array}} \right];{{\mathit{\boldsymbol{\hat X}}}^ * } = {{\left[ {\begin{array}{*{20}{c}} {{{\hat X}_1}}& \cdots &{{{\hat X}_n}} \end{array}} \right]}^{\rm{T}}};}\\ {{\mathit{\boldsymbol{S}}_i} = {{\mathit{\boldsymbol{\hat X}}}_i}\mathit{\boldsymbol{Q}}_i^ * {{\left( {\mathit{\boldsymbol{Q}}_0^ * + \sum\limits_{i = 1}^n {\mathit{\boldsymbol{\hat X}}_i^2\mathit{\boldsymbol{Q}}_i^ * } } \right)}^{ - 1}};}\\ {{\mathit{\boldsymbol{S}}_0} = - \mathit{\boldsymbol{Q}}_0^ * {{\left( {\mathit{\boldsymbol{Q}}_0^ * + \sum\limits_{i = 1}^n {\mathit{\boldsymbol{\hat X}}_i^2\mathit{\boldsymbol{Q}}_i^ * } } \right)}^{ - 1}}。} \end{array} $

(14)

式(14)即为改正数的计算公式，${\mathit{\boldsymbol{S}}_0}\left({{\mathit{\boldsymbol{L}}_0} - {\mathit{\boldsymbol{A}}^{\rm{*}}}{{\mathit{\boldsymbol{\hat X}}}^{\rm{*}}}} \right)$即为观测向量${\mathit{\boldsymbol{L}}_0}$的改正数，${\mathit{\boldsymbol{S}}_i}{\rm{\; }}\left({{\mathit{\boldsymbol{L}}_0} - {\mathit{\boldsymbol{A}}^{\rm{*}}}{{\mathit{\boldsymbol{\hat X}}}^{\rm{*}}}} \right)$为第i列观测向量${\mathit{\boldsymbol{L}}_i}$的改正数。由式(14)可以看出各列观测向量的改正数受到${\mathit{\boldsymbol{S}}_i}{\rm{\; }}$的影响，${\mathit{\boldsymbol{S}}_i}{\rm{\; }}$越大，则改正数越大。

进一步分析改正数的大小关系，设同一列的观测值为等精度观测值，不同列的观测值精度不同，即$\mathit{\boldsymbol{Q}}_i^{\rm{*}} = {q_i} \cdot \mathit{\boldsymbol{I}}{\rm{\; \; }}\left({{q_i} > 0} \right)$，其中${q_i}$为各列变量的协因数。将$\mathit{\boldsymbol{Q}}_i^{\rm{*}}$代入式(14)可得：

$ {\mathit{\boldsymbol{\hat e}}^ * } = - \frac{1}{s}\left[ {\begin{array}{*{20}{c}} { - {q_0}\left( {{\mathit{\boldsymbol{L}}_0} - {\mathit{\boldsymbol{A}}^ * }{{\mathit{\boldsymbol{\hat X}}}^ * }} \right)}\\ {{q_1}{{\hat X}_1}\left( {{\mathit{\boldsymbol{L}}_0} - {\mathit{\boldsymbol{A}}^ * }{{\mathit{\boldsymbol{\hat X}}}^ * }} \right)}\\ \vdots \\ {{q_n}{{\hat X}_n}\left( {{\mathit{\boldsymbol{L}}_0} - {\mathit{\boldsymbol{A}}^ * }{{\mathit{\boldsymbol{\hat X}}}^ * }} \right)} \end{array}} \right] $

(15)

式中，$s = {q_0} + \mathop {\mathop \sum \limits_{i = 1} }\limits^n \hat X_i^2{q_i}$。

由式(15)可得系数矩阵各列变量的改正数的表达式为$\mathit{\boldsymbol{\hat e}}_i^{\rm{*}} = - {s^{ - 1}}{q_i}{\hat X_i}\left({{\mathit{\boldsymbol{L}}_0} - {\mathit{\boldsymbol{A}}^{\rm{*}}}{{\mathit{\boldsymbol{\hat X}}}^{\rm{*}}}} \right)$。当$\left({{\mathit{\boldsymbol{L}}_0} - {\mathit{\boldsymbol{A}}^{\rm{*}}}{{\mathit{\boldsymbol{\hat X}}}^{\rm{*}}}} \right)$固定时，各列变量的改正数仅受$ - {s^{ - 1}}{q_i}{\hat X_i}$项的影响。这揭示了一个重要现象：在TLS平差中，改正数的大小与改正数所对应的待估参数值和观测值先验精度存在一定的比例关系。换言之，待估参数值与对应列的先验协因数的乘积越大，则该列的改正数越大；反之亦然。当q_i中至少有两个协因数相等时，不妨假设i=1列与i=2列的先验协因数相等(q₁=q₂)，则根据式(15)可得i=1列与i=2列观测向量的改正数比值约为${\hat X_1}/{\hat X_2}$。若${\hat X_1}$与${\hat X_2}$的差异非常大，则采用常规的稳健估计策略(即对设计阵的观测量采用同一个方差因子)会产生虚假稳健估计的现象，即${\mathit{\boldsymbol{L}}_1}$与${\mathit{\boldsymbol{L}}_2}$虽然是等精度，但${\mathit{\boldsymbol{L}}_1}$与${\mathit{\boldsymbol{L}}_2}$的改正数却因为${\hat X_1}$与${\hat X_2}$的不同呈现出不同的数量关系。因此在传统的RTLS方法^[12-16]中，会优先将X_max($\left| {{{\hat X}_1}} \right|$与$\left| {{{\hat X}_2}} \right|$之间的最大值)所对应列的观测值进行降权处理，从而导致不合理的稳健估计结果，本文称之为虚假稳健估计现象。

1.3 本文多变量稳健估计策略

针对具有多列独立变量的平差模型，已有的RTLS方法或是对全部的变量统一采用一个单位权中误差，或是对观测向量与系数矩阵分开采用两个单位权中误差，但均无法解决多变量情况下的虚假估计问题。因此，本文提出了适用于多变量EIV模型的多变量权函数计算策略(以Huber权函数为例)，称之为多变量稳健估计策略：

$ {w_{i,k}} = \left\{ {\begin{array}{*{20}{c}} {1,}&{\left| {{e_{i,k}}} \right| \le {c_k}}\\ {{c_k}/\left| {{e_{i,k}}} \right|,}&{\left| {{e_{i,k}}} \right| > {c_k}} \end{array}} \right. $

(16)

式中，${w_{i, k}}$为权因子，其中k表示第k列观测量，i表示第i个观测值；${e_{i, k}}$表示其改正数；${c_k}$为常数，一般取${c_k} = 2{\sigma _k}$，${\sigma _k}$表示第k列观测量的或然中误差。

分析式(16)可知，多变量权函数计算策略实际上是对不同列的独立观测量分别进行降权处理。这样处理是考虑到参数估值对改正数的影响，从而将降权处理仅限定在每个参数估值所对应的列，则在稳健估计时可避免形如§1.2所分析的虚假稳健估计的发生。

本文导出的RTLS方法因采用式(16)所示的多变量权函数计算策略，称之为MRTLS平差方法。结合式(7)~(10)和式(16)可知，MRTLS方法适用于EIV模型含有多列独立变量的情形。因采用了特征矩阵${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{A}}}$和${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{L}}}$的表达方式，因此MRTLS方法对于具有重复观测值、观测值间有函数关系等复杂结构的EIV模型同样适用。在实际应用中，具有复杂结构EIV模型的有高程异常拟合、遥感影像叶面积指数反演以及地面激光扫描标靶球定位等。对于独立观测量处于同一列，即系数矩阵不存在多类观测值的问题，如三维基准转换模型^[4]，其独立观测量集中于平差模型系数矩阵的第4列，使用本文方法时只需对三维基准转换模型的观测向量与系数矩阵的第4列向量分别进行降权处理即可。

2 虚假稳健估计现象及MRTLS方法的算例及分析

算例1  为验证§1.2中关于TLS方法的改正数比例关系分析的正确性，设计仿真算例：多元线性回归模型为${\mathit{\boldsymbol{L}}_0} = \mathop {\mathop \sum \limits_{i = 1} }\limits^6 {X_i}{\mathit{\boldsymbol{L}}_i}$，${X_1} \sim {X_6}$为待估参数，其值见表 1。${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_6}$为7组观测向量，其值由以下步骤生成：首先生成6组观测向量${\mathit{\boldsymbol{L}}_1} \sim {\mathit{\boldsymbol{L}}_6}$，每组观测向量由30个服从[0, 10]区间上均匀分布的观测值组成，然后代入${\mathit{\boldsymbol{L}}_0} = \mathop {\mathop \sum \limits_{i = 1} }\limits^6 {X_i}{\mathit{\boldsymbol{L}}_i}$求出${\mathit{\boldsymbol{L}}_{0{\rm{\; \; }}30 \times 1}}$，最后对观测向量${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_6}$分别加入不等精度的相互独立的高斯白噪声(同一组观测向量下的观测值等精度，不同组观测向量下的观测值不等精度)。加入的噪声与参数真值有关，观测向量${\mathit{\boldsymbol{L}}_1} \sim {\mathit{\boldsymbol{L}}_6}$加入的噪声信息见表 1，为方便分析，将加入的噪声转化为平差时的协因数${q_i}$，$i \in \left\{ {1, 2 \cdots 6} \right\}$。采用TLS方法计算${\mathit{\boldsymbol{L}}_1} \sim {\mathit{\boldsymbol{L}}_6}$的改正数组(平差时的随机模型根据加入的先验噪声而定)，将前5组误差方程的改正数结果列于表 2。

从表 1可以看出，1~6列的参数与协因数有不同的组合形式，包括参数与协因数完全相同(1、2列)，参数相同但协因数不同(1、4列)，协因数阵相同但参数不同(1、3列)，参数与协因数均不同、但两者乘积相同(1、5列)，1~6列的乘积比为1:1:2:2:1:10。对比分析表 2结果可以发现，${\mathit{\boldsymbol{L}}_1} \sim {\mathit{\boldsymbol{L}}_6}$观测向量的改正数比值近似为1:1:2:2:1:10，与表 1给出的比值完全相同，表明无论在哪种组合情形下，最终影响各列变量的改正数比例关系的是参数与协因数的乘积。表 2所示的平差结果验证了§1.2关于改正数比例关系分析的正确性。需要指出的是，这里的改正数比值近似而非完全为1:1:2:2:1:10，是因为参数估计结果较真值有偏差。

此外，该回归模型中，${\mathit{\boldsymbol{L}}_6}$的改正数是${\mathit{\boldsymbol{L}}_1} \sim {\mathit{\boldsymbol{L}}_5}$改正数的5~10倍，此时从理论上来讲，采用RTLS方法将会优先对第6列进行降权处理，从而造成虚假稳健估计现象。

算例2  为验证虚假估计现象的存在并比较分析本文多变量稳健估计方法的正确性和有效性，设计仿真算例如下。采用算例1的函数模型${\mathit{\boldsymbol{L}}_0} = \mathop {\mathop \sum \limits_{i = 1} }\limits^6 {X_i}{\mathit{\boldsymbol{L}}_i}$，随机模型采用以下4种不同方案：①对${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_6}$分别加入σ=0.001的相互独立高斯白噪声；②对${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_6}$分别加入σ=0.01的相互独立高斯白噪声；③对${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_6}$分别加入σ=0.1的相互独立高斯白噪声；④对${\mathit{\boldsymbol{L}}_0} \sim {\mathit{\boldsymbol{L}}_6}$分别加入σ=0.5的相互独立高斯白噪声。同时，对${\mathit{\boldsymbol{L}}_0}$的第1个观测值、${\mathit{\boldsymbol{L}}_1}$的第10个观测值以及${\mathit{\boldsymbol{L}}_0}$、${\mathit{\boldsymbol{L}}_2}$的第23个观测值分别加入50σ的粗差。方案①、②、③均重复仿真步骤1 000次，并分别采用基于Huber权函数的MRTLS(记为MRTLS-Huber)与基于Huber权函数的RTLS(记为RTLS-Huber)两种方法分别进行平差计算，统计两种平差方法获得合理平差结果的比例，结果见表 3。

从表 3可以看出，在不同的噪声情况下，RTLS-Huber方法的平差结果合理比例均为0，而MRTLS-Huber方法的平差结果合理比例均大于90%，高于RTLS-Huber方法。这表明当各列变量的改正数之间具有较大的倍数差距时，RTLS-Huber方法因为对各列变量的改正数统一进行降权处理，确实存在虚假稳健估计的现象，从而导致平差结果不合理。反观MRTLS-Huber方法因对各列变量的改正数分别进行降权处理而避免了虚假稳健估计现象，表明了本文方法的正确性和有效性。

算例3  为验证本文抗差方法的粗差定位能力以及参数估计结果的正确可靠性，设计仿真算例如下。采用算例1给出的函数模型${\mathit{\boldsymbol{L}}_0} = \mathop {\mathop \sum \limits_{i = 1} }\limits^6 {X_i}{\mathit{\boldsymbol{L}}_i}$，随机模型采用算例2的方案③。在不同的位置对观测量加入粗差，方案如下：⑤对观测向量${\mathit{\boldsymbol{L}}_0}$的第1个观测值加入50σ的粗差; ⑥在方案⑤的基础上再对系数矩阵${\mathit{\boldsymbol{L}}_1}$的第10个观测值加入50σ的粗差; ⑦在方案⑥的基础上再对${\mathit{\boldsymbol{L}}_0}$、${\mathit{\boldsymbol{L}}_2}$的第23个观测值加入50σ的粗差。

仿真1 000次，分别采用MRTLS-Huber、基于Huber权函数的稳健最小二乘方法(记为LS-Huber)和TLS平差方法进行计算。图 1(a)、图 1(b)和图 1(c)为MRTLS-Huber的粗差位置结果，图 1(d)、图 1(e)和图 1(f)为平差方法的单位权中误差结果。3种平差方法的参数估值与真值的均方根误差结果见表 4。

从表 4方案⑤结果可以发现，3种平差方法的参数均方根误差大小关系为：E_MRTLS-Huber < E_LS-Huber < E_TLS，并且随着粗差个数的增多，方案⑤、⑦呈现同样的大小关系。E_MRTLS-Huber < E_TLS以及E_LS-Huber < E_TLS体现了稳健估计的优越性，而E_MRTLS-Huber < E_LS-Huber则体现了TLS方法的优越性。本文MRTLS-Huber方法结合了TLS与稳健估计的特点，较LS-Huber和TLS方法能够获得更接近真值的参数估计结果，进一步表明本文方法的正确性和有效性。

图 1(a)、图 1(b)、图 1(c)只列出了方案⑤~⑦前100组的平差结果。可以看出，随着粗差类型与个数的增加，图 1(a)中在组号为1处存在一条由点组成的线，在图 1(b)的组号为1与10处、图 1(c)的组号为1、10和23处均出现了此现象。这表明在该次仿真中，当粗差处于观测向量、系数矩阵或是两者都存在的情况下，本文MRTLS-Huber方法均能够定位出粗差所对应的误差方程，说明本文方法正确、有效。

需要说明的是，本文方法能够探测出的粗差所对应的误差方程组号而非定位出粗差的具体位置，是因为TLS方法所获得的观测值改正数受到参数估值与协因数乘积的影响，若误差方程的某个观测值存在粗差^[19]，则经平差后的该行观测值均出现较大的改正数，从而对该行误差方程的每个观测值都会进行降权处理，造成只能探测出粗差所对应的误差方程的现象。

图 1(d)、图 1(e)、图 1(f)分别对应方案⑤、⑥、⑦下3种平差方法的内符合精度。分析结果可知，在不同的粗差位置下，MRTLS-Huber方法的内符合精度均高于LS-Huber与TLS方法的内符合精度。随着加入粗差个数的增多，MRTLS-Huber方法的精度明显高于TLS方法，体现本文方法在稳健估计方面的优越性，进一步表明本文方法的正确性和有效性。

3 结语

本文分析了现有RTLS方法的稳健估计策略在处理多变量EIV模型时产生的虚假稳健估计问题。为避免该问题，提出了MRTLS平差方法，并通过仿真算例验证了MRTLS方法的正确性和有效性。

1) 采用RTLS方法对多变量EIV模型进行平差计算时，模型的各列变量的改正数受其对应的参数估值与观测量先验精度的联合影响。两者的乘积越大，则该列改正数越大。因此，使用现有RTLS方法处理多独立变量EIV模型时，会优先对多变量EIV模型的改正数大的列向量进行降权处理，造成平差结果不合理甚至错误。

2) 提出了多变量稳健总体最小二乘方法。该方法对多变量EIV模型的各列变量分别进行降权处理，从而避免了虚假稳健估计现象的发生。仿真计算结果表明，MRTLS方法能够有效地减少虚假稳健估计现象的发生，定位出粗差所对应的误差方程以及给出可靠的参数估计结果。