Analysis of Sea and Land Water Storage Changes in the Western Pacific Under the Influence of ENSO
-
摘要: 利用2005-01-2016-12的GRACE卫星数据，结合测高、ARGO（array for real-time geostrophic oceanography）、TRMM（tropical rainfall measurement mission）降水数据和MEI（multivariate ENSO index），从海平面、海水质量、海水比容和陆地水储量的非季节性变化分析了厄尔尼诺-南方涛动（El Niño/Southern Oscillation，ENSO）对西太平洋及其沿岸地区的影响。在GRACE数据处理过程中考虑了地震构造信号的影响，并采用正演模型法来恢复时变信号。结果表明，西太平洋地区平均海平面（western Pacific mean sea level，WPMSL）变化对ENSO的响应表现为：厄尔尼诺期间海平面下降，拉尼娜期间海平面上升；尤其在2014-2016年厄尔尼诺期间，WPMSL下降了近26.2 mm。比容变化为影响WPMSL非季节性变化的主导因素。澳大利亚和中南半岛水储量异常在年际变化上也较好响应了ENSO，与MEI的相关系数分别为-0.61（滞后3个月）和-0.65（滞后8个月）。长江流域对ENSO的响应主要在中下游区域，2014-2016年厄尔尼诺造成了该时期长江流域水储量达到近几年的极大值。Abstract: In this paper, the gravity recovery and climate experiment (GRACE), altimetry, ARGO, TRMM (tropical rainfall measurement mission) data and multivariate ENSO index (MEI) from January 2005 to December 2016 are adopted to analyze the effects of ENSO on the Western Pacific and its coastal areas, from the aspects of sea level, oceanic mass, steric and the terrestrial water storage change. During the GRACE data processing, we remove the effect of tectonic signals, and use an iterative forward modeling method to recover time-varying signal. The result shows, the changes in the Western Pacific mean sea level (WPMSL) in response to ENSO are as follows:The sea level is lower during El Niño and higher during La Niña; in particular, during 2014-2016 El Nino, the WPMSL dropped by nearly 26.2 mm. The steric change is the dominant factor affecting the nonseasonal anomalies of the WPMSL. The nonseasonal anomalies of water storage in Australia and Indochina Peninsula also respond well to the ENSO, and the correlation coefficients with the MEI are -0.61 (lag three months) and -0.65 (lag eight months). The response of the Yangtze River basin to ENSO is mainly in the middle and lower basin, 2014-2016 El Niño caused the water storage in the Yangtze River basin to reach the maximum value in recent years.
-
Keywords:
- GRACE /
- forward modeling /
- sea level change /
- land water storage change /
- the Western Pacific /
- ENSO
-
整体最小二乘(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算法虽然总体表现良好，但随着粗差个数的增加，对于个别仿真数据抗差失效，而本文算法仍表现出了良好的抗差性能，总体抗差效果更稳定。
致谢: 感谢CSR提供的GRACE数据、AVISO提供的卫星测高数据、JAMSTC、SIO和IPRC提供的ARGO数据、GES DISC提供的TRMM数据以及NOAA提供的MEI数据；感谢武汉大学高性能计算中心提供的计算支持。 -
![]()
图 1 西太平洋及沿海地区地表质量变化速率分布
Figure 1. Distributions of Mass Change Rates in the Western Pacific and Coastal Areas
![]()
图 2 西太平洋地区海平面年际变化时间序列及MEI
Figure 2. Time Series of Interannual Variations of Sea Level in the Western Pacific and MEI
![]()
图 3 西太平洋地区海平面异常非趋势性变化序列与MEI比较
Figure 3. Comparison of Non-trending Interannual Variations of the Sea Level Change in the Western Pacific with MEI
![]()
图 4 西太平洋沿岸陆地水储量和降水量年际变化
Figure 4. Time Series of Interannual Variations of Terrestrial Water Storage and Land Precipitation in the Western Pacific Coastal Areas
![]()
图 5 ENSO影响下的西太平洋沿岸陆地水储量和降水量异常
Figure 5. Anomalies of Terrestrial Water Storage and Precipitation in the Western Pacific Coastal Areas Under the Influence of ENSO
![]()
图 6 不同滞后月份下陆地水储量变化与MEI的相关系数
Figure 6. Correlation Coefficients of Terrestrial Water Storage and MEI Under Different Lagging Months
表 1 不同时间段西太平洋地区平均海平面、海水质量、海水比容变化趋势/(mm∙a-1)
Table 1 Trends of Sea Level Change, Ocean Mass Change and Steric Change with Different Periods in the Western Pacific /(mm∙a-1)
海平面变化 时间段 与MEI的相关
系数2005-01—2016-12 2005-01—2010-12 2011-01—2016-12 GMSL 3.88±0.17 2.88±0.25 5.59±0.41 0.62 WPMSL 3.01±0.56 4.43±1.26 -1.65±1.49 -0.75 质量（GRACE）+比容（ARGO） 2.38±0.50 3.65±1.17 -1.61±1.28 -0.73 质量（GRACE） 2.72±0.19 1.91±0.45 3.83±0.54 0.44 比容（ARGO） -0.33±0.54 1.74±1.15 -5.44±1.35 -0.78 
下载: 导出CSV
-
[1] Nerem R S, Chambers D P, Choe C, et al. Estimating Mean Sea Level Change from the TOPEX and Jason Altimeter Missions[J]. Marine Geodesy, 2010, 33(s1):435-446 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.1080/01490419.2010.491031
[2] Boening C, Willis J K, Landerer F W, et al. The 2011 La Niña:So Strong, the Oceans Fell[J]. Geophysical Research Letters, 2012, 39(19):1-5 doi: 10.1029/2012GL053055/abstract
[3] 袁媛, 周宁芳, 李崇银.中国华北雾霾天气与超强El Nino事件的相关性研究[J].地球物理学报, 2017, 60(1):11-21 http://www.cnki.com.cn/Article/CJFDTotal-DQWX201701003.htm Yuan Yuan, Zhou Ningfang, Li Chongyin. Correlation Between Haze in North China and Super El Niño Events[J]. Chinese J Geophys, 2017, 60(1):11-21 http://www.cnki.com.cn/Article/CJFDTotal-DQWX201701003.htm
[4] 冯伟, 钟敏, 许厚泽.联合卫星测高、卫星重力和海洋浮标资料研究2005-2013年的全球海平面变化[J].地球物理学进展, 2014, 29(2):471-477 http://www.cnki.com.cn/Article/CJFDTotal-DQWJ201402001.htm Feng Wei, Zhong Min, Xu Houze. Global Sea Level Changes Estimated from Satellite Altimetry, Satellite Gravimetry and Argo Data During 2005-2013[J]. Progress in Geophysics, 2014, 29(2):471-477 http://www.cnki.com.cn/Article/CJFDTotal-DQWJ201402001.htm
[5] 张保军, 王泽民.联合卫星重力、卫星测高和海洋资料研究全球海平面变化.武汉大学学报·信息科学版, 2015, 40(11):1453-1459 http://ch.whu.edu.cn/CN/abstract/abstract3363.shtml Zhang Baojun, Wang Zemin. Global Sea Level Variations Estimated from Satellite Altimetry, GRACE and Oceanographic Data[J]. Geomatics and Information Science of Wuhan University, 2015, 40(11):1453-1459 http://ch.whu.edu.cn/CN/abstract/abstract3363.shtml
[6] Chambers D P, Cazenave A, Champollion N, et al. Evaluation of the Global Mean Sea Level Budget Between 1993 and 2014[J]. Surveys in Geophysics, 2017, 38(1):1-19 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=29764a83e008fcfc0187da7080dc107d
[7] Willis J K, Chambers D P, Nerem R S. Assessing the Globally Averaged Sea Level Budget on Seasonal to Interannual Timescales[J]. Journal of Geophysical Research:Oceans, 2008, 113(C6):205-210 doi: 10.1029/2007JC004517/pdf
[8] Chambers D P. Evaluation of New GRACE Time-Variable Gravity Data over the Ocean[J]. Geophysical Research Letters, 2006, 331(17):1-5 doi: 10.1029-2006GL027296/
[9] Chambers D P, Bonin J A. Evaluation of Release-05 GRACE Time-Variable Gravity Coefficients over the Ocean[J]. Ocean Science, 2012, 8(5):859-868 doi: 10.5194/os-8-859-2012
[10] Landerer F W, Swenson S C. Accuracy of Scaled GRACE Terrestrial Water Storage Estimates[J]. Water Resources Research, 2012, 48(4):W04531 http://d.old.wanfangdata.com.cn/NSTLQK/NSTL_QKJJ0227162423/
[11] Yi S, Wang Q, Sun W. Basin Mass Dynamic Changes in China from GRACE Based on a Multibasin Inversion Method[J]. Journal of Geophysical Research:Solid Earth, 2016, 121(5):3782-3803 doi: 10.1002/2015JB012608
[12] Han S, Riva R, Sauber J, et al. Source Parameter Inversion for Recent Great Earthquakes from a Decade-Long Observation of Global Gravity Fields[J]. Journal of Geophysical Research:Solid Earth, 2013, 118(3):1240-1267 doi: 10.1002/jgrb.50116
[13] Swenson S, Chambers D, Wahr J. Estimating Geocenter Variations from a Combination of GRACE and Ocean Model Output[J]. Journal of Geophysical Research:Atmospheres, 2008, 113(B8):194-205 doi: 10.1029-2007JB005338/
[14] Wahr J, Molenaar M, Bryan F. Time Variability of the Earth's Gravity Field:Hydrological and Oceanic Effects and Their Possible Detection Using GRACE[J]. Journal of Geophysical Research, 1998, 103(B12):30205-30230 doi: 10.1029/98JB02844
[15] Peltier W R, Argus D F, Drummond R. Space Geodesy Constrains Ice Age Terminal Deglaciation:The Global ICE-6G_C (VM5a) Model[J]. Journal of Geophysical Research:Solid Earth, 2015, 120(1):450-487 doi: 10.1002/2014JB011176
[16] Broerse T, Riva R, Simons W, et al. Postseismic GRACE and GPS Observations Indicate a Rheology Contrast Above and Below the Sumatra Slab[J]. Journal of Geophysical Research:Solid Earth, 2015, 120(7):5343-5361 doi: 10.1002/2015JB011951
[17] Chen J L, Wilson C R, Tapley B D. Contribution of Ice Sheet and Mountain Glacier Melt to Recent Sea Level Rise[J]. Nature Geoscience, 2013, 6(7):549-552 doi: 10.1038/ngeo1829
[18] Yi S, Sun W, Heki K, et al. An Increase in the Rate of Global Mean Sea Level Rise Since 2010[J]. Geophysical Research Letters, 2015, 42(10):3998-4006 doi: 10.1002/2015GL063902
[19] Yi S, Sun W, Feng W, et al. Anthropogenic and Climate-Driven Water Depletion in Asia[J]. Geophysical Research Letters, 2016, 43:9061-9069 doi: 10.1002/2016GL069985
[20] Riser S C, Freeland H J, Roemmich D, et al. Fifteen Years of Ocean Observations with the Global Argo Array[J]. Nature Climate Change, 2016, 6(2):145-153 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=958878ab7fb8df8bdfc6482ad3863c32
[21] Masters D, Nerem R S, Choe C, et al. Comparison of Global Mean Sea Level Time Series from TOPEX/Poseidon, Jason-1, and Jason-2[J]. Marine Geodesy, 2012, 35(s1):20-41 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.1080/01490419.2012.717862
[22] Piecuch C G, Quinn K J. El Niño, La Niña, and the Global Sea Level Budget[J]. Ocean Science, 2016, 12(6):1-18
[23] Zhang Z, Chao B F, Chen J, et al. Terrestrial Water Storage Anomalies of Yangtze River Basin Droughts Observed by GRACE and Connections with ENSO[J]. Global & Planetary Change, 2015, 126:35-45 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=54f378c4e6eb6a2aded09085714fb370
[24] Fasullo J T, Boening C, Landerer F W, et al. Australia's Unique Influence on Global Sea Level in 2010-2011[J]. Geophysical Research Letters, 2013, 40(16):4368-4373 doi: 10.1002/grl.50834
-
期刊类型引用(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)
计量
- 文章访问数:
- HTML全文浏览量:
- PDF下载量:
- 被引次数: 11
