Robust Weight Total Least Squares Algorithm of Correlated Observation Based on Median Parameter Method
-
摘要: 在抗差加权整体最小二乘算法中，抗差模型的抗差性与初值的好坏关系极大，若以最小二乘或整体最小二乘估值作为初值，必定会受到粗差污染而影响其抗差性。考虑到观测向量和系数矩阵存在相关性，首先推导了部分变量误差（partial errors-in-variables，Partial EIV）模型的加权整体最小二乘算法，在此基础上提出了一种利用中位参数法求解抗差迭代初值的相关观测抗差加权整体最小二乘算法。然后采用中位参数法确定抗差初值，考虑到可能出现的粗差对观测空间与结构空间的综合影响，基于标准化残差构造权因子函数，实现其抗差解法。仿真实验结果表明，此算法具有良好的抗差性能，其参数估计结果比传统算法精度更高，且随着粗差个数的增加，其抗差稳定性较好。
-
关键词:
- 抗差加权整体最小二乘 /
- 中位参数法 /
- Partial EIV模型 /
- 相关观测 /
- 标准化残差
Abstract: In the robust weighted total least squares(RWTLS) algorithm, its robustness of the robust model is highly related to the initial values. If the least squares or total least squares estimates is used as the initial value, it will be affected by gross error, and certainly impacted the robust characteristics of RWTLS estimates. Considering the correlation between the observed vector and the coefficient matrix, we first deduce the weighted least-squares solution of Partial-EIV model, and a new RWTLS algorithm of correlated observation is proposed to solve the initial values of robust iterations by using the median parameter method. Then the median parameter method is used to determine the initial value, and on this basis we propose a new robust estimated method, which is based on the standardized residual error and considered the influence of gross error both on observation and structure spaces. The experiment results show that the proposed estimated method has a good performance to resist gross error, and the presented solution is more accurate than the traditional method for line fitting, and with the increase of the number of gross errors, the stability of the algorithm is superior to the traditional method. -
整体最小二乘(total least squares，TLS)估计可以有效解决变量误差(errors in variables，EIV)模型中系数矩阵含有随机误差的问题，国内外学者对此展开了广泛的研究，产生了大量研究成果。文献[1]通过扩展最小二乘(least squares，LS)估计，最早提出了将观测向量和系数矩阵的残差平方和最小化的TLS估计准则；文献[2]提出了著名的奇异值分解法，并正式命名了TLS估计方法的通用术语，也使得TLS估计引起了各领域的广泛关注；文献[3]提出了附加等式约束的TLS估计，自此，TLS估计方法才真正应用到大地测量领域。近年来，国内外学者针对加权整体最小二乘(weighted total least squares，WTLS)算法进行了深入的研究。文献[4]提出了基于拉格朗日函数的WTLS迭代算法；文献[5]基于Gauss-Helmert模型研究了平面坐标转换的WTLS估计；文献[6]提出了非线性条件下WTLS的高斯-牛顿迭代法；文献[7-9]研究了一般性权矩阵、非线性、系数矩阵部分元素为固定量等情形下的WTLS算法；文献[10]从非线性LS入手，将WTLS算法扩展到部分变量误差(Partial EIV)模型中，推导了有限样本下WTLS算法的精度评定公式，可以将各种形式下的系数矩阵纳入统一的模型下求解，且有效减少了未知数的数目。
上述讨论仅考虑了观测向量和系数矩阵中观测值包含的随机误差，在实际测量中观测值往往会受到粗差污染，而LS和WTLS估计并不具备抵御粗差的能力，因此，有必要进一步研究抗差加权整体最小二乘算法(robust weighted total least squares，RWTLS)。在LS估计中，粗差的处理方法一般可归纳为两类：(1)基于均值漂移模型，把粗差归入函数模型，研究粗差探测方法[11]，(2)基于方差-协方差膨胀模型，将粗差归入随机模型，研究抗差估计方法[12-14]。在TLS估计中，国内外学者进行了一些关于RWTLS的探讨，文献[15]基于Gauss-Helmert模型和Huber权函数提出了三维坐标转换的稳健WTLS方法；文献[16]提出了WTLS的选权迭代算法，并用于直线拟合和坐标转换中；文献[17]将文献[15]中的Huber权函数替换成了IGG函数，并验证了该方法的优势；文献[18-19]基于Gauss-Helmert模型研究了三维坐标转换和GPS高程拟合中的RWTLS算法；文献[20]认识到在上述关于RWTLS的研究中均采用残差构造权因子函数，无法顾及结构空间的抗差性，对于加权情况下的抗差效果不理想，因此基于标准化残差提出了一种RWTLS算法，并推导了残差协因数阵的表达式。然而，上述算法均存在共同的缺陷：(1)相关算法未考虑观测向量和系数矩阵的相关性；(2)由于等价权抗差估计的抗差性与初值好坏的关系极大，而抗差迭代前的初值往往采用LS和TLS估计的结果，受到粗差的影响而容易失真。文献[21]提出采用基于中位数法的抗差总体最小二乘算法，但最终只用到少量观测数据计算参数估值，丢失大量的有效信息，按文中分组选取子样并不十分合理，且抗差权因子函数由残差构造，并不能顾及结构空间抗差。
鉴于上述分析，本文首先结合Partial EIV模型和标准化残差各自的优势，基于高斯-牛顿迭代法推导了基于Partial EIV模型相关观测的RWTLS算法，然后利用中位数法求解抗差迭代的初值，代入RWTLS算法中重新定权，进行参数求解。
1 相关观测的Partial EIV模型的高斯-牛顿迭代解法
观测向量和系数矩阵完全相关的Partial EIV模型的函数模型和随机模型可表示为：
$$ \left. \begin{array}{l} \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} = \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\left( {\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B}}\left( {\mathit{\boldsymbol{a}} - {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{a}}}} \right)} \right)\\ \mathit{\boldsymbol{v}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}}\\ {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{a}}}} \end{array}} \right] \sim \left( {\left[ {\begin{array}{*{20}{c}} 0\\ 0 \end{array}} \right]\sigma _0^2\mathit{\boldsymbol{Q}}} \right) \end{array} \right\} $$ (1) 式中，y和vy分别表示m×1维观测向量及其对应的随机误差向量；X为t维参数向量；h为非随机元素构成的mt维向量；B为由常数构成的mt×n维系数矩阵；a为随机元素对应的n维列向量；va为a的随机误差向量；随机模型中σ0为单位权中误差；Q为v对应的协因数阵，
$$ \mathit{\boldsymbol{Q}} = \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) 式中，Qy和Qa分别表示vy和va对应的协因数阵；Qya和Qay分别表示互协因数阵。
将式(1)中的函数模型在(X(i), v a(i))处进行一阶泰勒级数展开可得：
$$ \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} + {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{X}} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{a}}} $$ (3) 将式(3)整理可得：
$$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} + {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{X}} + \left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]\mathit{\boldsymbol{v}} $$ (4) 式中，
$$ \mathit{\boldsymbol{A}} = {\rm{ivec}}\left( {\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}}} \right) $$ (5) $$ {\mathit{\boldsymbol{A}}_{\left( i \right)}} = {\rm{ivec}}\left( {\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B}}\left( {\mathit{\boldsymbol{a}} - {\mathit{\boldsymbol{v}}_{\mathit{\boldsymbol{a}}\left( i \right)}}} \right)} \right) = \mathit{\boldsymbol{A}} - {\rm{ivec}}\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{v}}_{\mathit{\boldsymbol{a}}\left( i \right)}}} \right) $$ (6) 式中，ivec表示将列向量恢复为指定维数矩阵的算子。
构造WTLS的拉格朗日目标函数：
$$ \begin{array}{*{20}{c}} {\mathit{\Phi }\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{K}}} \right) = {\mathit{\boldsymbol{v}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}^{ - 1}}\mathit{\boldsymbol{v}} + 2{\mathit{\boldsymbol{K}}^{\rm{T}}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - } \right.}\\ {\left. {{\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{X}} - \left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]\mathit{\boldsymbol{v}}} \right)} \end{array} $$ (7) 对各变量求偏导并令偏导数为零：
$$ \frac{1}{2}\frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{v}}}} = {\mathit{\boldsymbol{Q}}^{ - 1}}\mathit{\boldsymbol{\hat v}} - {\left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]^{\rm{T}}}\mathit{\boldsymbol{\hat K}} = 0 $$ (8) $$ \frac{1}{2}\frac{{\partial \mathit{\Phi }}}{{\partial {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{X}}}} = - \mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{\hat K}} = 0 $$ (9) $$ \begin{array}{*{20}{c}} {\frac{1}{2}\frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{K}}}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat X}} - }\\ {\left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]\mathit{\boldsymbol{\hat v}} = 0} \end{array} $$ (10) 由(8)可得：
$$ \mathit{\boldsymbol{\hat v}} = \mathit{\boldsymbol{Q}}{\left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]^{\rm{T}}}\mathit{\boldsymbol{\hat K}} $$ (11) 将式(11)代入式(10)可得：
$$ \begin{array}{l} \mathit{\boldsymbol{\hat K}} = {\left( {\mathit{\boldsymbol{MQ}}{\mathit{\boldsymbol{M}}^{\rm{T}}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat X}}} \right) = \\ \;\;\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat X}}} \right) \end{array} $$ (12) 式中，
$$ \begin{array}{l} \mathit{\boldsymbol{M}} = \left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]\\ \;\;\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{c\left( i \right)}} = \left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]\mathit{\boldsymbol{Q}}\left[ {{\mathit{\boldsymbol{I}}_m} - } \right.}\\ {{{\left. {\left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]}^{\rm{T}}}} \end{array} \end{array} $$ (13) 将式(12)代入式(9)可以推导出：
$$ {\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}} = {\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)^{ - 1}} \cdot \mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}}} \right) $$ (14) 则在(i+1)次迭代后更新的参数估值为：
$$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}} = {\mathit{\boldsymbol{X}}_{\left( i \right)}} + {\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}} = {{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1} \cdot }\\ {\left( {\mathit{\boldsymbol{y}} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{v}}_{\mathit{\boldsymbol{a}}\left( i \right)}}} \right)} \end{array} $$ (15) $$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat v}}}_{\left( {i + 1} \right)}} = {\mathit{\boldsymbol{Q}}^{ - 1}}{{\left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]}^{\rm{T}}} \cdot }\\ {\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}}} \right)} \end{array} $$ (16) $$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat v}}}_{\mathit{\boldsymbol{y}}\left( {i + 1} \right)}} = \left( {{\mathit{\boldsymbol{Q}}_y} + {\mathit{\boldsymbol{Q}}_{ya}}{{\left( { - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right)}^{\rm{T}}}} \right) \cdot }\\ {\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}}} \right)} \end{array} $$ (17) $$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat v}}}_{\mathit{\boldsymbol{a}}\left( {i + 1} \right)}} = \left( {{\mathit{\boldsymbol{Q}}_{ay}} + {\mathit{\boldsymbol{Q}}_a}{{\left( { - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right)}^{\rm{T}}}} \right) \cdot }\\ {\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}}} \right)} \end{array} $$ (18) 将最小二乘的参数估值作为初值，按照式(14)至式(18)进行迭代计算，当满足$ \parallel \delta \mathit{\boldsymbol{\hat X}}{_{(i + 1)}}\parallel < \theta $ (θ为给定小量)时，停止迭代，此时可以获得Partial EIV的WTLS解。
当观测向量和系数矩阵不相关时，Qya = 0，Qay = 0，式(17)和式(18)可变为：
$$ {{\mathit{\boldsymbol{\hat v}}}_{\mathit{\boldsymbol{y}}\left( {i + 1} \right)}} = {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}}} \right) $$ (19) $$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat v}}}_{\mathit{\boldsymbol{a}}\left( {i + 1} \right)}} = - {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_{\left( i \right)}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{y}} - } \right.}\\ {\left. {\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}}} \right)} \end{array} $$ (20) 2 基于中位参数法的RWTLS算法
2.1 基于标准化残差的降权策略
为了同时顾及观测向量和结构空间抗差，在利用等价权抗差原理进行参数的抗差估计时，需要构建检验量，其关键在于残差向量的协因数阵推导。根据式(17)和式(18)得到第(i+1)次迭代的预测残差向量：
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{v}}_{i + 1}} = \left[ \begin{array}{l} {{\mathit{\boldsymbol{\hat v}}}_{\mathit{\boldsymbol{y}}\left( {i + 1} \right)}}\\ {{\mathit{\boldsymbol{\hat v}}}_{\mathit{\boldsymbol{a}}\left( {i + 1} \right)}} \end{array} \right] = }\\ {\left[ \begin{array}{l} \left( {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}} + {\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}{{\left( { - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right)}^{\rm{T}}}} \right)\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\\ \left( {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}} + {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}{{\left( { - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right)}^{\rm{T}}}} \right)\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1} \end{array} \right] \cdot }\\ {\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}}} \right) = \left[ \begin{array}{l} {\mathit{\boldsymbol{\omega }}_i}\\ {\mathit{\boldsymbol{\psi }}_i} \end{array} \right]{\mathit{\boldsymbol{\varphi }}_{i + 1}}} \end{array} $$ (21) 式中，vi+1为第(i+1)次迭代的残差向量；$ {{\mathit{\boldsymbol{\hat v}}}_{\mathit{\boldsymbol{y}}(i + 1)}}$和$ \mathit{\boldsymbol{\hat v}}{_{\mathit{\boldsymbol{a}}{\rm{ }}(i + 1)}}$分别表示观测向量和随机元素向量对应的第(i+1)次迭代的残差向量；$\mathit{\boldsymbol{\omega }}{_i} = \left( {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}} + {\rm{ }}{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}\left( { - \left( {\mathit{\boldsymbol{X}}_{(i)}^{\rm{T}} \otimes {\rm{ }}\mathit{\boldsymbol{I}}{_m}} \right)\mathit{\boldsymbol{B}}} \right){^{{\rm{T}}}}} \right){\rm{ }}\mathit{\boldsymbol{Q}}_{c(i)}^{ - 1};{\rm{ }} $ $ \mathit{\boldsymbol{\psi }}{_i} = {\rm{ }}\left( {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}\mathit{\boldsymbol{}} + {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}\mathit{\boldsymbol{}}\left( { - \left( {\mathit{\boldsymbol{X}}_{(i)}^{\rm{T}} \otimes {\rm{ }}\mathit{\boldsymbol{I}}{_m}} \right)\mathit{\boldsymbol{B}}} \right){^{{\rm{T}}}}} \right)\mathit{\boldsymbol{Q}}_{c(i)}^{ - 1};$ $ \mathit{\boldsymbol{\varphi }}{_{i + 1}} = {\rm{ }}\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{{\rm A}X}}{_{(i)}} - \mathit{\boldsymbol{A}}{_{(i)}}\mathit{\boldsymbol{\delta \hat X}}{_{(i + 1)}}$。
则第(i+1)次迭代的预测残差向量的协因数阵为：
$$ {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{v}}_{i + 1}}}} = \left[ \begin{array}{l} {\mathit{\boldsymbol{\omega }}_i}\\ {\mathit{\boldsymbol{\psi }}_i} \end{array} \right] \cdot {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{\varphi }}_{i + 1}}}} \cdot \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\omega }}_i^{\rm{T}}}&{\mathit{\boldsymbol{\psi }}_i^{\rm{T}}} \end{array}} \right] $$ (22) 将式(14)代入φi+1得：
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\varphi }}_{i + 1}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat X}}}_{\left( {i + 1} \right)}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}} - }\\ {{\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}} \cdot \mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}}} \right) = }\\ {\left( {{\mathit{\boldsymbol{I}}_m} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}} \right) \cdot \left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}}} \right)} \end{array} $$ (23) 根据协因数传播率可以求得[20]：
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{\varphi }}_{i + 1}}}} = \left( {{\mathit{\boldsymbol{I}}_m} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}} \right) \cdot }\\ {{\mathit{\boldsymbol{Q}}_{\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}}} \right)}}{{\left( {{\mathit{\boldsymbol{I}}_m} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}} \right)}^{\rm{T}}} = }\\ {\left( {{\mathit{\boldsymbol{I}}_m} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}} \right) \cdot }\\ {{\mathit{\boldsymbol{Q}}_{c\left( i \right)}}{{\left( {{\mathit{\boldsymbol{I}}_m} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}} \right)}^{\rm{T}}} = }\\ {{\mathit{\boldsymbol{Q}}_{c\left( i \right)}} - {\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}} \end{array} $$ (24) 将式(24)代入式(22)可得：
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{v}}_{i + 1}}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\omega }}_i}}\\ {{\mathit{\boldsymbol{\psi }}_i}} \end{array}} \right] \cdot {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{\varphi }}_{i + 1}}}} \cdot \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\omega }}_i^{\rm{T}}}&{\mathit{\boldsymbol{\psi }}_i^{\rm{T}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\omega }}_i}}\\ {{\mathit{\boldsymbol{\psi }}_i}} \end{array}} \right]\left( {{\mathit{\boldsymbol{Q}}_{c\left( i \right)}} - } \right.}\\ {\left. {{\mathit{\boldsymbol{A}}_{\left( i \right)}}{{\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}} \right) \cdot \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\omega }}_i^{\rm{T}}}&{\mathit{\boldsymbol{\psi }}_i^{\rm{T}}} \end{array}} \right]} \end{array} $$ (25) 由于上述推导获得的残差协因数阵与迭代过程中的参数和系数阵残差有关，故残差协因数阵需要在迭代中不断更新。
在此基础上，对式(21)的残差进行标准化：
$$ {{\hat v}_i} = \frac{{{v_i}}}{{{\sigma _0}\left( {\sqrt {{q_{{v_i}}}} } \right)}}\left( {{q_{{v_i}}} \ne 0} \right) $$ (26) 式中，${{\hat v}_i} $为标准化残差向量${\mathit{\boldsymbol{\hat v}}} $中的第i个元素；vi表示残差向量v中的第i个元素；σ0为单位权中误差；qvi为残差协因数阵Qv中对角线上第i个元素。σ0由中位函数计算，可表示为：
$$ {\sigma _0} = \mathop {{\rm{med}}}\limits_{i = 1}^t \left( {\left| {{v_i}/\sqrt {{q_{{v_i}}}} } \right|} \right) \cdot 1.482\;6\left( {{q_{{v_i}}} \ne 0} \right) $$ (27) 基于标准化残差进行抗差估计，可选的权因子函数众多，本文采用IGG 3权因子函数，为了计算方便，采用其对应的协因数因子函数：
$$ {R_{ii}} = \left\{ \begin{array}{l} 1.0,\left| {{{\hat v}_i}} \right| \le {k_0}\\ \frac{{\left| {{{\hat v}_i}} \right|}}{{{k_0}}}{\left( {\frac{{{k_1} - {k_0}}}{{{k_1} - \left| {{{\hat v}_i}} \right|}}} \right)^2},{k_0} < \left| {{{\hat v}_i}} \right| \le {k_1}\\ {10^{10}},{k_1} < \left| {{{\hat v}_i}} \right| \end{array} \right. $$ (28) 式中，第3段理论上应该为无穷大，为了方便计算，可采用极大值替代(如1×1010)，满足实际要求；k0取为2.0~3.0，k1取为4.0~8.0[20]。
采用双因子模型，由式(28)可得观测向量和系数矩阵中观测值的等价协因数：
$$ {{\bar q}_{ij}} = {q_{ij}}\sqrt {{R_{ii}}} \sqrt {{R_{jj}}} $$ (29) 式中，${{\bar q}_{ij}} $和${q_{ij}} $分别为观测值的等价协因数阵${\mathit{\boldsymbol{\bar Q}}} $和协因数阵Q中的对应位置元素。采用${\mathit{\boldsymbol{\bar Q}}} $替代Q，构建类似于式(7)的加权整体最小二乘问题的拉格朗日目标函数，对各变量求导并令导数为零，类似于式(8)至式(10)，可得第(i+1)次迭代的待求参数抗差估计的改正数为：
$$ {\rm{ \mathsf{ δ} }}{{\mathit{\boldsymbol{\hat {\bar X}}}}_{\left( {i + 1} \right)}} = {\left( {\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{\bar Q}}_{c\left( i \right)}^{ - 1}{\mathit{\boldsymbol{A}}_{\left( i \right)}}} \right)^{ - 1}}\mathit{\boldsymbol{A}}_{\left( i \right)}^{\rm{T}}\mathit{\boldsymbol{\bar Q}}_{c\left( i \right)}^{ - 1} \cdot \left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\left( i \right)}}} \right) $$ (30) 式中，
$$ {{\mathit{\boldsymbol{\bar Q}}}_{c\left( i \right)}} = \left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]\mathit{\boldsymbol{\bar Q}} \cdot {\left[ {{\mathit{\boldsymbol{I}}_m} - \left( {\mathit{\boldsymbol{X}}_{\left( i \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{B}}} \right]^{\rm{T}}} $$ 采用类似于式(14)至式(18)的迭代过程，在给定的阈值条件下，可获得参数的最终估值及其精度评定指标。
2.2 基于中位参数法确定抗差初值
文献[21]中基于中位数法的思路为：(1)采用中位数法求解分组后的中位数参数估值，并用其计算观测向量与系数矩阵中含有的误差向量的估值；(2)依据绝对误差中位数与中误差的关系，分类确定观测向量与系数矩阵中的观测元素对应的中误差，并应用抗差权函数对观测元素进行重新分类定权；(3)根据重新定义的权值分组求解参数估值，并求取其参数的中位数值，计算各组模型参数的估值与中位参数的差值向量，以各组中差值向量的最小范数对应的参数估值作为参数的最终估值。
然而，文献[21]方法在解算最终参数估值的过程中使用的是分组的观测数据，即最终只用到少量观测数据计算参数估值，损失了大量的有效信息。而本文只采用中位数法确定抗差迭代前的初始值，在抗差迭代过程中使用了全部的观测数据，理论上更加合理。
式(1)函数模型中有m×1维观测向量，待估参数$ {\mathit{\boldsymbol{\hat X}}}$为t维向量(t < m)，因此，采用相关观测的Partial EIV模型可以获得Cmt组参数估值，可表示为：
$$ \mathit{\boldsymbol{\hat X}} = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_1}}&{{{\mathit{\boldsymbol{\hat X}}}_2}}& \cdots &{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{C}}_m^t}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {x_1^1}&{x_2^1}& \cdots &{x_{{\rm{C}}_m^t}^1}\\ {x_1^2}&{x_2^2}& \cdots &{x_{{\rm{C}}_m^t}^2}\\ \vdots &{}&{}& \vdots \\ {x_1^t}&{x_2^t}& \cdots &{x_{{\rm{C}}_m^t}^t} \end{array}} \right] $$ (31) 式中，x表示待估参数元素。
对式(31)中获得的Cmt组参数估值(行向量)求中位数，可得到待估参数${\mathit{\boldsymbol{\hat X}}} $的t个参数的中位数值，表示为：
$$ {{\mathit{\boldsymbol{\hat X}}}_{{\rm{med}}}} = {\left[ {\begin{array}{*{20}{c}} {x_{{\rm{med}}}^1}&{x_{{\rm{med}}}^2}& \cdots &{x_{{\rm{med}}}^t} \end{array}} \right]^{\rm{T}}} $$ (32) 式中，下标med表示中位数值。
计算Cmt组估值与中位参数向量的差值，并求二次范数$ \parallel {{\rm{d}}_{{{\mathit{\boldsymbol{\hat X}}}_i}}}\parallel $，可表示为：
$$ \left\| {{{\rm{d}}_{{{\mathit{\boldsymbol{\hat X}}}_i}}}} \right\| = \left\| {{{\mathit{\boldsymbol{\hat X}}}_i} - {{\mathit{\boldsymbol{\hat X}}}_{{\rm{med}}}}} \right\|,i \in \left[ {1,2 \cdots {\rm{C}}_m^t} \right] $$ (33) 式中，$ \parallel {{\rm{d}}_{{{\mathit{\boldsymbol{\hat X}}}_i}}}\parallel $的最小值对应的参数估值${{\mathit{\boldsymbol{\hat X}}}_i} $为所求的中位参数解，即抗差迭代前的初值，根据式(17)和式(18)可以分别求出观测向量和系数矩阵对应的残差向量，再根据式(21)至式(29)进行重新定权，重新定权后根据标准化残差的降权策略(见§2.1)进行抗差估计即可。
应用中位参数法进行参数估值的求取，可以有效降低受观测值含有的粗差污染的影响。众所周知，单参数中位数法的崩溃污染率接近50%，当一半以上的观测值正常时，其中位数也是正常的。因此，当观测值中含有k个粗差时，则参数的Cmt组参数估值中，不含粗差的解向量为Cmt个。为保证中位参数估值的抗差性，需要满足Cm-kt≥Cmt/2，此时，可确保中位参数解由不含粗差的观测值求得。
3 算例与分析
假设一直线方程为y=5x+9，采用均匀分布函数在(0，18)之间选择18个x值，并计算对应y值，将对应的18个点位坐标(x，y)视为真值，并加入非等精度的随机误差，其数学期望为0，标准差在0~0.05间随机产生。考虑坐标点存在的相关性，设同一点位内x和y的相关系数$ {\rho _{{x_i}{y_i}}}$均为0.6，点间的相关系数$ {\rho _{{x_i}{x_j}}}$和${\rho _{{y_i}{y_j}}} $均为0.3，在实际计算中，通过${\sigma _{{x_i}{y_i}}} = {\rho _{{x_i}{y_i}}}{\sigma _{{x_i}}}{\sigma _{{y_i}}}, {\sigma _{{x_i}{x_j}}} = {\rho _{{x_i}{x_j}}}{\sigma _{{x_i}}}{\sigma _{{x_j}}}, {\sigma _{{y_i}{y_j}}} = {\rho _{{y_i}{y_j}}}{\sigma _{{y_i}}}{\sigma _{{y_j}}} $计算观测数据的协方差，完全相关观测情况下的协方差阵为(其中σ02=1)：
$$ \mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\sigma _0^2}}{{\sigma _{{x_1}}^2}}}& \cdots &{\frac{{\sigma _0^2}}{{{\sigma _{{x_1}{x_{18}}}}}}}&{\frac{{\sigma _0^2}}{{{\sigma _{{x_1}{y_1}}}}}}& \cdots &{\frac{{\sigma _0^2}}{{{\sigma _{{x_1}{y_{18}}}}}}}\\ \vdots &{}&{}&{}&{}& \vdots \\ {\frac{{\sigma _0^2}}{{{\sigma _{{x_{18}}{x_1}}}}}}& \cdots &{\frac{{\sigma _0^2}}{{\sigma _{{x_{18}}}^2}}}&{\frac{{\sigma _0^2}}{{\sigma _{{x_{18}}{y_1}}^2}}}& \cdots &{\frac{{\sigma _0^2}}{{\sigma _{{x_{18}}{y_{18}}}^2}}}\\ {\frac{{\sigma _0^2}}{{\sigma _{{y_1}{x_1}}^2}}}& \cdots &{\frac{{\sigma _0^2}}{{\sigma _{{y_1}{x_{18}}}^2}}}&{\frac{{\sigma _0^2}}{{\sigma _{{y_1}}^2}}}& \cdots &{\frac{{\sigma _0^2}}{{\sigma _{{y_1}{y_{18}}}^2}}}\\ \vdots &{}&{}&{}&{}& \vdots \\ {\frac{{\sigma _0^2}}{{\sigma _{{y_{18}}{x_1}}^2}}}& \cdots &{\frac{{\sigma _0^2}}{{\sigma _{{y_{18}}{x_{18}}}^2}}}&{\frac{{\sigma _0^2}}{{\sigma _{{y_{18}}{y_1}}^2}}}& \cdots &{\frac{{\sigma _0^2}}{{\sigma _{{y_{18}}}^2}}} \end{array}} \right] $$ (34) 加入的粗差个数分别为1、2、3，粗差大小介于标准差绝对值的5~20倍之间，位置随机产生。反复进行500次模拟，方案设计如下。
1) 方案1：加入粗差前，使用WTLS估计。
2) 方案2：加入粗差后，使用WTLS估计。
3) 方案3：加入粗差后，使用基于标准化残差的RWTLS方法(文献[20])。
4) 方案4：加入粗差后，使用基于中位参数法的RWTLS方法。
上述方案中，方案1和方案2分别是在未加入粗差前和加入粗差后使用相关观测的Partial EIV模型进行求解；方案3是基于文献[20]的算法；方案4的算法是将文献[20]的算法推广到相关观测的情形，且在此基础上利用了中位参数法相关观测的RWTLS算法。
现有的RWTLS方法在进行抗差迭代前使用的初值一般采用WTLS的结果，而受到粗差污染的WTLS方法所求得的抗差迭代初值会严重失真。当粗差个数为3时，方案3和方案4抗差迭代初值的偏差序列分别如图 1和图 2所示，可以看出方案4中使用中位参数法求出的抗差迭代初值明显优于方案3中使用WTLS方法的结果。
![]()
图 1 方案3和方案4抗差迭代初值(斜率a)偏差序列对比Figure 1. Deviation Sequence Comparison of Resistance Values (Slope a) Between Schemes 3 and 4![]()
图 2 方案3和方案4抗差迭代初值(截距b)偏差序列对比Figure 2. Deviation Sequence Comparison of Resistance Values (Intercept b) Between Schemes 3 and 4当粗差个数分别为1、2、3时，采用4种方案进行参数估计，分析不同方案的效果如下。
表 1至表 3分别给出了在不同粗差个数时各方案直线拟合斜率a与截距b的均方根误差(root mean square error, RMSE)和最大绝对偏差(max)的对比情况。
表 1 1个粗差时各方案的统计结果Table 1. Statistical Results of Each Scheme When There Exists One Gross Error方案 RMSE(a) RMSE(b) max (|a|) max (|b|) 1 0.009 6 0.104 0 0.008 1 0.512 1 2 0.022 8 0.213 9 0.155 8 1.444 7 3 0.010 1 0.111 1 0.073 3 0.635 2 4 0.009 8 0.107 2 0.072 4 0.629 9 表 2 2个粗差时各方案的统计结果Table 2. Statistical Results of Each Scheme When There Exists Two Gross Errors方案 RMSE(a) RMSE(b) max (|a|) max (|b|) 1 0.008 8 0.123 2 0.031 2 0.480 5 2 0.029 9 0.426 0 0.183 7 2.819 0 3 0.013 4 0.292 7 0.105 1 2.819 0 4 0.010 5 0.128 6 0.054 8 1.342 5 表 3 3个粗差时各方案的统计结果Table 3. Statistical Results of Each Scheme When There Exists Three Gross Errors方案 RMSE(a) RMSE(b) max (|a|) max (|b|) 1 0.008 9 0.118 3 0.027 8 0.349 7 2 0.037 9 0.428 4 0.210 4 2.749 9 3 0.021 1 0.309 8 0.197 6 2.729 5 4 0.014 2 0.161 4 0.064 9 0.796 2 分析表 1至表 3得出以下结论：(1)由于方案1中没有混入粗差，所求得的均方根误差和最大绝对偏差均最小；(2)由于方案2中混入了粗差，均方根误差和最大绝对偏差均明显大于方案1，说明混入粗差的观测值使得WTLS方法求得的参数估值严重失真，WTLS方法缺乏抵御粗差的能力；(3)当粗差个数为1时，方案3采用了基于标准化残差的RWTLS方法，相较方案2，其均方根误差和最大绝对偏差均显著减小，说明方案3在粗差个数较少时，由于同时顾及了观测空间和结构空间抗差，其模型具有较好的抵御粗差的能力；但随着粗差个数增加到2和3时，方案3中最大绝对偏差相较方案2并未发生显著变化，说明在考虑相关观测的情况下，方案3对于个别仿真数据，抗差会失效，整体抗差能力在减弱；(4)当粗差个数分别为1、2、3时，方案4的均方根误差相较于方案3均进一步减小，且随着粗差个数的增加，效果更加明显。当粗差个数为3时，方案4相比于方案3，斜率a和截距b分别减小了约32.70%和47.90%，抗差效果改善明显，尤其是最大绝对偏差显著减小，分别减小了约67.16%和70.83%，说明方案4的抗差方案稳定性较好，相较其他方案，效果总体最佳。
图 3和图 4给出了当粗差个数为3时，方案3和方案4直线拟合斜率a与截距b的解算结果与真值的偏差序列。从图 3、图 4中可以看出，方案3中个别偏差值较大，部分仿真数据抗差效果不佳，而方案4的拟合结果与真值偏差序列总体比方案3更稳定，即对于方案3抗差失效的数据，方案4仍能实现有效抗差估计，抗差性能稳定。
![]()
图 3 方案3和方案4拟合结果(斜率a)对比Figure 3. Fitting Results (Slope a) Comparison Between Schemes 3 and 4![]()
图 4 方案3和方案4拟合结果(截距b)对比Figure 4. Fitting Results (Intercept b) Comparison Between Schemes 3 and 44 结语
本文在推导相关观测Partial-EIV模型的WTLS算法的基础上，结合标准化残差提出了一种基于中位参数法求解抗差迭代初值的相关观测RWTLS迭代算法。本文算法解决了抗差迭代初值受粗差污染会严重失真的问题，并顾及了观测空间和结构空间的抗差性，同时考虑了坐标点位的相关性。算例结果表明：(1)与基于标准化残差的RWTLS算法相比，本文算法得到的斜率与截距的均方根误差显著减小，算法的精度较现有的算法有一定的提高；(2)考虑到点位之间的相关性，通过比较基于标准化残差方案的RWTLS算法和本文算法，可以看出基于标准化残差的RWTLS算法虽然总体表现良好，但随着粗差个数的增加，对于个别仿真数据抗差失效，而本文算法仍表现出了良好的抗差性能，总体抗差效果更稳定。
-
![]()
图 1 方案3和方案4抗差迭代初值(斜率a)偏差序列对比
Figure 1. Deviation Sequence Comparison of Resistance Values (Slope a) Between Schemes 3 and 4
![]()
图 2 方案3和方案4抗差迭代初值(截距b)偏差序列对比
Figure 2. Deviation Sequence Comparison of Resistance Values (Intercept b) Between Schemes 3 and 4
![]()
图 3 方案3和方案4拟合结果(斜率a)对比
Figure 3. Fitting Results (Slope a) Comparison Between Schemes 3 and 4
![]()
图 4 方案3和方案4拟合结果(截距b)对比
Figure 4. Fitting Results (Intercept b) Comparison Between Schemes 3 and 4
表 1 1个粗差时各方案的统计结果
Table 1 Statistical Results of Each Scheme When There Exists One Gross Error
方案 RMSE(a) RMSE(b) max (|a|) max (|b|) 1 0.009 6 0.104 0 0.008 1 0.512 1 2 0.022 8 0.213 9 0.155 8 1.444 7 3 0.010 1 0.111 1 0.073 3 0.635 2 4 0.009 8 0.107 2 0.072 4 0.629 9 
下载: 导出CSV
表 2 2个粗差时各方案的统计结果
Table 2 Statistical Results of Each Scheme When There Exists Two Gross Errors
方案 RMSE(a) RMSE(b) max (|a|) max (|b|) 1 0.008 8 0.123 2 0.031 2 0.480 5 2 0.029 9 0.426 0 0.183 7 2.819 0 3 0.013 4 0.292 7 0.105 1 2.819 0 4 0.010 5 0.128 6 0.054 8 1.342 5
下载: 导出CSV
表 3 3个粗差时各方案的统计结果
Table 3 Statistical Results of Each Scheme When There Exists Three Gross Errors
方案 RMSE(a) RMSE(b) max (|a|) max (|b|) 1 0.008 9 0.118 3 0.027 8 0.349 7 2 0.037 9 0.428 4 0.210 4 2.749 9 3 0.021 1 0.309 8 0.197 6 2.729 5 4 0.014 2 0.161 4 0.064 9 0.796 2
下载: 导出CSV
-
[1] Adcock R J. Note on the Method of Least Squares[J]. The Analyst, 1877, 4(6):183-184 doi: 10.2307/2635777
[2] Golub H G, Vanloan F C. An Analysis of the Total Least Squares Problem[J]. SIAM Journal on Numerical Analysis, 1980, 17(6):883-893 doi: 10.1137/0717073
[3] Schaffrin B, Felus Y A. On Total Least-Squares Adjustment with Constraints[J]. IAG-Symp, 2005, 128:417-421 doi: 10.1007/3-540-27432-4_71
[4] Schaffrin B, Felus Y A. On the Multivariate Total Least-Squares Approach to Empirical Coordinate Transformations[J]. Journal of Geodesy, 2008, 82(6):373-383 doi: 10.1007/s00190-007-0186-5
[5] Neitzel F. Generalization of Total Least-Squares on Example of Unweighted and Weighted 2D Similarity Transformation[J]. Journal of Geodesy, 2010, 84(12):751-762 doi: 10.1007/s00190-010-0408-0
[6] Shen Y Z, Li B F, Chen Y. An Iterative Solution of Weighted Total Least-Squares Adjustment[J]. Journal of Geodesy, 2011, 85(4):229-238 doi: 10.1007/s00190-010-0431-1
[7] Fang X. Weighted Total Least Squares Solutions for Applications in Geodesy[D]. Germany: Leibniz University Hannover, 2011
[8] Fang X. Weighted Total Least Squares:Necessary and Sufficient Conditions, Fixed and Random Parameters[J]. Journal of Geodesy, 2013, 87(8):733-749 doi: 10.1007/s00190-013-0643-2
[9] 邓才华, 周拥军, 朱建军, 等.一类新函数模型及通用加权总体最小二乘平差方法[J].中国矿业大学学报, 2015, 44(5):952-958 http://d.old.wanfangdata.com.cn/Periodical/zgkydxxb201505025 Deng Caihua, Zhou Yongjun, Zhu Jianjun, et al. General Weighted Total Least Squares Method for a Type of New Functional Model[J]. Journal of China University of Mining andTechnology, 2015, 44(5):952-958 http://d.old.wanfangdata.com.cn/Periodical/zgkydxxb201505025
[10] Xu P L, Liu J N, Shi C. Total Least Squares Adjustment in Partial Errors-in-Variables Models:Algorithm and StatisticalAnalysis[J]. Journal of Geo-desy, 2012, 86(8):661-675 doi: 10.1007/s00190-012-0552-9
[11] 欧吉坤.粗差的拟准检定法(QUAD法)[J].测绘学报, 1999, 28(1):15-20 doi: 10.3321/j.issn:1001-1595.1999.01.004 Ou Jikun. Quasi-Accurate Detection of Gross Errors (QUAD)[J]. Acta Geodaetica et Cartographica Sinica, 1999, 28(1):15-20 doi: 10.3321/j.issn:1001-1595.1999.01.004
[12] Yang Y X. Robust Estimation of Geodetic Datum Transformation[J]. Journal of Geodesy, 1999, 73(5):268-274 doi: 10.1007/s001900050243
[13] 周江文, 黄幼才, 杨元喜, 等.抗差估计最小二乘法[M].武汉:华中理工大学出版社, 1997 Zhou Jiangwen, Huang Youcai, Yang Yuanxi, et al. Robust Least Squares Method[M]. Wuhan:Huazhong University of Science and Technology Press, 1997
[14] 高井祥, 张华海, 余学祥.矿区GPS网坐标转换的抗差模型[J].中国矿业大学学报, 1999, 28(2):99-103 doi: 10.3321/j.issn:1000-1964.1999.02.001 Gao Jingxiang, Zhang Huahai, Yu Xuexiang. Study on Robust Estimation Model for Coordinate Transformation of GPS Network in Ming Areas[J]. Journal of China University of Mining and Technology, 1999, 28(2):99-103 doi: 10.3321/j.issn:1000-1964.1999.02.001
[15] 陈义, 陆珏.以三维坐标转换为例解算稳健总体最小二乘方法[J].测绘学报, 2012, 41(5):715-722 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201205237128 Chen Yi, Lu Jue. Peforming 3D Similarity Transformation by Robust Total Least Squares[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(5):715-722 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201205237128
[16] MahboubV, Amiri-Simkooei A R, Sharifi M A.Iteratively Reweighted Total Least Squares:A Robust Estimation in Errors-in-Variables Models[J]. Surv Rev, 2013, 45(329):92-99 doi: 10.1179/1752270612Y.0000000017
[17] 龚循强, 李志林.稳健加权总体最小二乘法[J].测绘学报, 2014, 43(9):888-894 http://d.old.wanfangdata.com.cn/Periodical/kckxjs201301005 Gong Xunqiang, Li Zhilin. A Robust Weighted Total Least Squares Method[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(9):888-894 http://d.old.wanfangdata.com.cn/Periodical/kckxjs201301005
[18] Lu J, Chen Y, Li B, et al. Robust Total Least Squares with Reweighting Iteration for Three-Dimensional Similarity Transformation[J]. Surv Rev, 2014, 46(334):28-36 doi: 10.1179/1752270613Y.0000000050
[19] Tao Y Q, Gao J X, Yao Y F. TLS Algorithm for GPS Height Fitting Based on Robust Estimation[J]. Surv Rev, 2014, 46(336):184-188 doi: 10.1179/1752270613Y.0000000083
[20] Wang B, Li J C, Liu C. A Robust Weighted Total Least Squares Algorithm and Its Geodetic Applications[J]. Studia Geophysica et Geodaetica, 2016, 60(2):177-194 doi: 10.1007/s11200-015-0916-8
[21] 陶叶青, 高井祥, 姚一飞.基于中位数法的抗差总体最小二乘估计[J].测绘学报, 2016, 45(3):297-301 http://d.old.wanfangdata.com.cn/Periodical/chxb201603008 Tao Yeqing, Gao Jingxiang, Yao Yifei. Solution for Robust Total Least Squares Estimation Based on Median Method[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(3):297-301 http://d.old.wanfangdata.com.cn/Periodical/chxb201603008
-
期刊类型引用(8)
1. 齐志军，方兴，吕志鹏. 两种适用于线性回归EIV模型的高崩溃污染率算法. 武汉大学学报(信息科学版). 2025(01): 74-82 . 
百度学术
2. 王中，彭飞，韩玉超，孟庆旭，邓为耀. 非均匀采样加权最小二乘圆拟合的潜艇承压检测数据处理方法. 武汉大学学报(信息科学版). 2024(06): 970-976 . 百度学术
3. Jianjun ZHU，Leyang WANG，Jun HU，Bofeng LI，Haiqiang FU，Yibin YAO. Recent Advances in the Geodesy Data Processing. Journal of Geodesy and Geoinformation Science. 2023(03): 33-45 . 必应学术
4. 郭绍禹. 三原则整体最小二乘用于平面自由网平差计算. 测绘工程. 2022(01): 30-34 . 百度学术
5. 刘江，谭思伦，蔡伯根，王剑. 基于车路信息交互的车辆卫星定位协同定权方法. 交通运输系统工程与信息. 2022(05): 85-96 . 百度学术
6. 张志友，许长辉，杨映新，李冰峰. 不动产测量数据最小二乘融合技术. 工程勘察. 2022(11): 63-67 . 百度学术
7. 宋迎春，宋采薇，左廷英. 带有界不确定性的加权混合估计方法. 武汉大学学报(信息科学版). 2020(07): 949-955 . 百度学术
8. 傅英坤. 应力-渗流-化学多场耦合作用下混凝土蠕变特性试验研究——以龙口港工程为例. 长江科学院院报. 2020(09): 135-141 . 百度学术
其他类型引用(13)
计量
- 文章访问数: 1525
- HTML全文浏览量: 183
- PDF下载量: 257
- 被引次数: 21
