The Stability of GPS/BDS Inter-system Biases at the Receiver End and its Effect on Ambiguity Resolution
-
摘要: GNSS接收机端的UPD与接收到的信号频率有关，这导致GPS和BDS系统间的双差模糊度不具有整数特性，为了恢复其整数特性，两系统间的系统偏差需要进行估计或改正。在顾及GPS和BDS之间的时间系统、坐标系统和频率间偏差的基础上，推导出GPS/BDS系统偏差计算模型，并利用不同实验对系统偏差的稳定性进行验证。实验结果表明，不同品牌接收机在GPS/BDS系统偏差方面在一定条件下均具有稳定性；天线类型和天线连接线长度没有对GPS/BDS系统偏差产生显著影响。加入系统偏差改正的GPS/BDS紧组合定位在恶劣环境下表现良好，可将模糊度固定平均所需时间缩短33%，模糊度固定成功率提高31%。Abstract: Uncalibrated phase delay at the receiver end may depend on the frequency of the signal received, which leads to mixed double difference ambiguity between GPS and BDS and, lost integer character. In order to solve this problem, inter-system biases between GPS and BDS need to be estimated or corrected. Taking the biases of time system, coordinate system and frequency between GPS and BDS into consideration, a computing model of GPS-BDS inter-system bias is deduced, based on which, different experiments verify the stability of inter-system biases. Experimental results show that GPS-BDS inter-system biases of different brand receivers have stability under certain conditions. GPS-BDS inter-system biases are nott influenced by antenna type and antenna cable length significantly. GPS-BDS tight combining with inter-system biases corrected performs well under poor observing condition, and improves the mean time-to-ambiguity-fixing by 33%, and the success rate of ambiguity-fixing by 31%.
-
Keywords:
- GPS /
- BDS /
- uncalibrated phase delay /
- inter-system bias /
- ambiguity resolution
-
整体最小二乘(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算法虽然总体表现良好，但随着粗差个数的增加，对于个别仿真数据抗差失效，而本文算法仍表现出了良好的抗差性能，总体抗差效果更稳定。
-
![]()
图 2 良好观测环境下卫星个数和PDOP值
Figure 2. Number of Satellites and PDOP under Good Observing Condition
![]()
图 4 恶劣观测环境下卫星个数和PDOP值
Figure 4. Number of Satellites and PDOP under Poor Observing Condition
表 1 2014年6月19-21日利用A接收机短基线数据计算ISB/周
Table 1 The Result of ISB Using Short Baseline Data of Receiver A on June 19, 2014/Cycle
时段编号 WHU1-WHU2 WHU1-WHU3 WHU2-WHU3 6月19日 6月20日 6月21日 6月19日 6月20日 6月21日 6月19日 6月20日 6月21日 1 0.00 0.05 0.05 0.01 0.05 0.05 0.00 0.06 0.00 2 0.01 0.04 0.05 0.01 0.03 0.04 0.01 0.05 0.01 3 0.01 0.01 0.06 0.01 0.04 0.04 0.00 0.04 0.01 平均值 0.01 0.03 0.05 0.01 0.04 0.04 0.00 0.05 0.01 
下载: 导出CSV
表 2 2014年12月15日利用A接收机短基线数据计算ISB/周
Table 2 The Result of ISB Using Short Baseline Data of Receiver A on December 15, 2014/Cycle
时段编号 WHU1-WHU2 WHU1-WHU3 WHU2-WHU3 1 0.01 0.00 0.01 2 0.02 0.01 0.01 3 0.01 0.00 0.01 平均值 0.01 0.00 0.01
下载: 导出CSV
表 3 2014年7月1-3日利用B接收机短基线数据计算ISB/周
Table 3 The Result of ISB Using Short Baseline Data of Receiver B on July 1, 2014/Cycle
时段编号 WHU1-WHU2 WHU1-WHU3 WHU2-WHU3 7月1日 7月2日 7月3日 7月1日 7月2日 7月3日 7月1日 7月2日 7月3日 1 0.16 0.00 0.18 0.28 0.22 0.40 0.14 0.23 0.28 2 0.12 0.00 0.15 0.27 0.23 0.41 0.16 0.21 0.27 3 0.11 0.00 0.16 0.24 0.23 0.43 0.13 0.21 0.30 平均值 0.13 0.00 0.16 0.26 0.23 0.41 0.14 0.22 0.28
下载: 导出CSV
表 4 不同品牌接收机短基线数据计算ISB/周
Table 4 The Result of ISB Using Short Baseline Data of Different Brand Receivers/Cycle
接收机品牌 基线长度/m 时段1 时段2 时段3 C 4.27 -0.01 -0.06 -0.01 D 54.21 -0.12 0.48 0.48 E 8.42 -0.01 -0.01 -0.01 F 25.47 0.05 -0.29 -0.05 G 25.17 0.02 -0.12
下载: 导出CSV
表 5 2014年7月3日、4日和7日利用A和B接收机短基线数据计算ISB/周
Table 5 The Result of ISB Using Short Baseline Data of A and B Receivers on July 3, 2014/Cycle
时段编号 WHU1-WHU2 WHU1-WHU3 WHU2-WHU3 7月3日 7月4日 7月7日 7月3日 7月4日 7月7日 7月3日 7月4日 7月7日 1 0.03 0.49 0.23 0.18 0.23 0.19 0.11 -0.28 -0.06 2 0.04 0.46 0.23 0.13 0.23 0.14 0.08 -0.26 -0.06 3 0.02 0.46 0.22 0.15 0.24 0.14 0.09 -0.27 -0.06 平均值 0.03 0.47 0.23 0.15 0.23 0.16 0.09 -0.27 -0.06
下载: 导出CSV
表 6 混合接收机短基线数据计算ISB/周
Table 6 The Result of ISB Using Short Baseline Data of Mixed Receivers/Cycle
混合接收机 基线长度/m 时段1 时段2 时段3 C-E 5.27 -0.09 -0.07 -0.09 D-F 44.59 0.33 0.38 -0.44
下载: 导出CSV
表 7 不同天线类型条件下ISB计算结果/周
Table 7 The Result of ISB Under the Condition of Different Antenna Types/Cycle
时段编号 WHU1-WHU2 WHU1-WHU3 WHU2-WHU3 1 0.03 0.02 0.02 2 0.03 0.04 0.01 3 0.08 0.07 0.02 平均值 0.05 0.04 0.02
下载: 导出CSV
表 8 不同天线连接线长度条件下
Table 8 The Result of ISB Under the Condition of Different Antenna Cable Length/Cycle
时段编号 WHU1-WHU2 WHU1-WHU3 WHU2-WHU3 1 0.07 0.03 0.02 2 0.01 0.03 0.01 3 0.03 0.06 0.06 平均值 0.04 0.04 0.03
下载: 导出CSV
表 9 良好观测环境下不同定位模式对模糊度固定影响
Table 9 Impact of Ambiguity Resolution on Different Positioning Mode under Good Observing Condition
组合类型 模糊度固定平均所需时间/历元 计算时长为30 s的固定成功率/% 松组合 2.38 88 紧组合 2.29 91
下载: 导出CSV
表 10 恶劣观测环境下不同定位模式对模糊度固定影响
Table 10 Impact of Ambiguity Resolution on Different Positioning Mode under Poor Observing Condition
组合类型 模糊度固定平均所需时间/历元 计算时长为5 s的固定成功率/% 松组合 6 65 紧组合 4 85
下载: 导出CSV
-
[1] 王世进, 秘金钟, 李得海, 等. GPS/BDS的RTK定位算法研究[J].武汉大学学报·信息科学版, 2014, 39(5):621-625 http://ch.whu.edu.cn/CN/abstract/abstract2980.shtml Wang Shijin, Bei Jinzhong, Li Dehai, et al. Real-time Kinematic Positioning Algorithm of GPS/BDS[J]. Geomatics and Information Science of Wuhan University, 2014, 39(5), 621-625 http://ch.whu.edu.cn/CN/abstract/abstract2980.shtml
[2] Zhang W, Cannon M E, Julien O, et al. Investigation of Combined GPS/GALILEO Cascading Ambiguity Resolution Schemes[C]. Proceedings of the 16th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS/GNSS '03), Portland, Ore, USA, 2003
[3] Julien O, Cannon M E, Alves P, et al. Triple Frequency Ambiguity Resolution Using GPS/Galileo[J]. European Journal of Navigation, 2004, 2(2):51-57 doi: 10.1007/s10291-013-0333-9
[4] Deng C, Tang W, Liu J, et al. Reliable Single-epoch Ambiguity Resolution for Short Baselines Using Combined GPS/BeiDou System[J]. GPS Solutions, 2014, 18(3):375-386 doi: 10.1007/s10291-013-0337-5
[5] Teunissen P J G, Odolinski R, Odijk D. Instantaneous BeiDou/GPS RTK Positioning with High cut-off Elevation Angles[J]. Journal of Geodesy, 2014, 88(4):335-350 doi: 10.1007/s00190-013-0686-4
[6] Odijk D, Teunissen P J G. Characterization of Between-receiver GPS-Galileo Inter-system Biases and Their Effect on Mixed Ambiguity Resolution[J]. GPS Solutions, 2013, 17(4):521-533 doi: 10.1007/s10291-012-0298-0
[7] Paziewski J, Wielgosz P. Accounting for Galileo-GPS Inter-system Biases in Precise Satellite Positioning[J]. Journal of Geodesy, 2015, 89(1):81-93 doi: 10.1007/s00190-014-0763-3
[8] Montenbruck O, Hauschild A, Hessels U. Characterization of GPS/GIOVE Sensor Stations in the CONGO Network[J]. GPS Solutions, 2011, 15(3):193-205 doi: 10.1007/s10291-010-0182-8
[9] 楼益栋, 龚晓鹏, 辜声峰, 等. GPS/BDS混合双差分RTK定位方法及结果分析[J].大地测量与地球动力学, 2016, 36(1):1-5 Lou Yidong, Gong Xiaopeng, Gu Shengfeng, et al. An Algorithm and Results Analysis for GPS+BDS Inter-System Mix Double-Difference RTK[J]. Journal of Geodesy and Geodynamics, 2016, 36(1), 1-5
[10] Wanninger L. Carrier-phase Inter-frequency Biases of GLONASS Receivers[J]. Journal of Geodesy, 2012, 86(2):139-148 doi: 10.1007/s00190-011-0502-y
[11] 程鹏飞, 文汉江, 成英燕, 等. 2000国家大地坐标系椭球参数与GRS80和WGS84的比较[J].测绘学报, 2009, 38(3):189-194 http://jz.docin.com/p-503223503.html Cheng Pengfei, Wen Hanjiang, Cheng Yingyan, et al. Parameters of the CGCS2000 Ellipsoid and Comparisons with GRS80 and WGS84[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(3):189-194 http://jz.docin.com/p-503223503.html
[12] Hofmann-wellenhof B, Lichtenegger H, Wasle E. GNSS-Global Navigation Satellite Systems:GPS, GLONASS, Galileo, and More[M]. Berlin:Springer, 2008
[13] Blewitt G. Carrier Phase Ambiguity Resolution for the Global Positioning System Applied to Geodetic Baselines up to 2000 km[J]. Journal of Geophysical Research, 1989, 94(B8):10187-10203 doi: 10.1029/JB094iB08p10187
[14] Zhang H, Li X, Zhu L, et al. Research on GNSS System Time Offset Monitoring and Prediction[C]. China Satellite Navigation Conference (CSNC) 2014, Nanjing, China, 2014
[15] Dach R, Schaer S, Lutz S, et al. Combining the Observations from Different GNSS[C]. Proceedings of EUREF 2010 Symposium, Gävle, Sweden, 2010
[16] Wang J. An Approach to GLONASS Ambiguity Resolution[J]. Journal of Geodesy, 2000, 74(5):421-430 doi: 10.1007/s001900000096
[17] Teunissen P J G. The Least-squares Ambiguity Decorrelation Adjustment:A Method for Fast GPS Integer Ambiguity Estimation[J]. Journal of Geodesy, 1995, 70(1-2):65-82 doi: 10.1007/BF00863419
[18] Nadarajah N, Teunissen P J G, Raziq N. BeiDou Inter-satellite-type Bias Evaluation and Calibration for Mixed Receiver Attitude Determination[J]. Sensors, 2013, 13(7):9435-9463 doi: 10.3390/s130709435
-
期刊类型引用(7)
1. 倪育德，闫苗玉，刘瑞华. 基于DOA-BP神经网络的电离层TEC短期预测. 航空学报. 2024(04): 192-205 . 
百度学术
2. 吴文坛，王新广，石少坚，史进志，郧晓光. 基于河北省CORS的区域电离层模型精度分析. 地理空间信息. 2023(01): 148-151 . 百度学术
3. 蒋磊，孙蕊，刘正午，徐成，梁的达，胡德振. 基于GA-BP的中欧GNSS电离层误差建模与精度分析. 北京航空航天大学学报. 2023(06): 1533-1542 . 百度学术
4. 赵奕源，吴文坛，赵春梅，秘金钟，莫雁寒，田时雨，李得海. 不同活跃状态下香港区域电离层建模分析. 测绘科学. 2023(08): 72-80 . 百度学术
5. 孙帆，宁一鹏，邢建平，王森，代培培. 单站区域实时电离层异常改正拟合模型研究. 电子器件. 2023(05): 1294-1299 . 百度学术
6. 章淑君，邱蕾，陆浩楠，夏朋飞，叶世榕. 利用北斗GEO卫星观测监测深圳市周围地区电离层空间环境. 测绘地理信息. 2022(05): 12-16 . 百度学术
7. 张兴汉. Klobuchar与GIM模型对GNSS单点定位精度影响分析. 测绘与空间地理信息. 2022(12): 184-187 . 百度学术
其他类型引用(4)
计量
- 文章访问数: 1982
- HTML全文浏览量: 290
- PDF下载量: 473
- 被引次数: 11
