Discharge Prediction in Tidal Reach Using Harmonic Analysis and VMD-BP Neural Network
-
摘要: 针对感潮河段受径流和潮汐共同影响导致流量预报精度较低的问题，结合变分模态分解（variational mode decomposition，VMD）处理非平稳信号的能力与反向(back propagation，BP)神经网络处理非线性问题的优势，提出了基于潮流调和分析及VMD-BP神经网络组合模型的感潮河段流量预报方法。首先，采用潮流调和分析方法对原始流量进行潮流、余流分离；然后，根据误差逆向传播算法构建BP神经网络，并对潮流数据和经VMD处理后的余流数据进行仿真训练；最后，将仿真训练输出的潮流和余流分量叠加重构进而得到最终的流量预报结果。在长江口徐六泾断面开展流量预报实验，结果表明，单独采用BP神经网络方法相对于传统潮流调和分析方法的流量预报精度提高了约3 400 m3/s，相对精度提高了约6%；所提组合模型方法的流量预报精度相对于传统潮流调和分析方法提高了约5 000 m3/s，相对精度提高了约9%。基于调和分析及VMD-BP神经网络的组合模型可以有效提高感潮河段流量预报精度，同时也为流况多变水域的流量预报提供新思路。Abstract:Objectives In view of the low accuracy of flow forecast in tidal reach due to the combined influence of runoff and tide, combining the ability of variational mode decomposition (VMD) to process non-stationary signals and the advantages of back propagation (BP) neural network to deal with nonlinear problems, this paper proposes a new method for discharge prediction in tidal reach based on the combination model of tidal current harmonic analysis and VMD-BP neural network.Methods First, the harmonic analysis method is used to separate the tidal discharge and residual discharge from the original discharge. Second, a BP neural network is constructed based on error backpropagation algorithm. Then, the BP neural network is used to simulate and train the power flow data and residual current data processed by VMD. Finally, the output tidal discharge and residual discharge components are superimposed and reconstructed to obtain the final discharge prediction results.Results Experiments performed at the Xuliujing section of the Yangtze Estuary show that the discharge prediction accuracy of BP neural network method is better than that of traditional harmonic analysis method, with the root mean square error (RMSE) is decreased by about 3 400 m3/s, and the relative standard accuracy (RSD) is increased by about 6%. In addition, compared with the traditional method, RMSE of the proposed method in this paper is decreased by about 5 000 m3/s, and the RSD is increased by about 9%.Conclusions The combined model based on harmonic analysis and VMD-BP neural network can effectively improve the prediction accuracy of discharge in tidal reach. At the same time, it also provides a new idea for the discharge prediction in the waters with variable flow conditions.
-
整体最小二乘(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算法虽然总体表现良好，但随着粗差个数的增加，对于个别仿真数据抗差失效，而本文算法仍表现出了良好的抗差性能，总体抗差效果更稳定。
-
![]()
图 4 短期、中期及长期的原始流量调和分析结果
Figure 4. Harmonic Analysis of Raw Flow Results of Short-Term, Medium-Term and Long-Term
![]()
图 6 2020⁃01⁃03全天原始、潮流、余流数据集的BP神经网络预报结果
Figure 6. Prediction Results of BP Neural Network of Origin, Current and Residual Datasets in Jan.3, 2020
![]()
图 7 趋势、高频、低频数据集的BP神经网络预报
Figure 7. Prediction Results of BP Neural Network of Trend, High Frequency and Residual Low Frequency Datasets
![]()
图 8 VMD重构后的余流BP神经网络预报结果
Figure 8. BP Neural Network Prediction Results of Residual Discharge After VMD Reconstruction
![]()
图 9 基于潮流调和分析及VMD⁃BP神经网络组合模型的原始流量预报结果
Figure 9. Prediction Results of Original Discharge Based on Current Harmonic Analysis and VMD-BP Neural Network Combination Model
![]()
图 10 小潮、中潮及大潮时段3种方法的流量预报结果
Figure 10. Forecast Results Using Three Methods in Three Periods
表 1 BP神经网络的参数设置
Table 1 Parameter Setting of BP Neural Network
网络结构（n×l×m） 参数 参数设置 96×15×48 迭代次数 1 000 学习速率 0.01 期望误差 0.000 1 传递函数 tansig（） 
下载: 导出CSV
表 2 3种方法的流量预报结果
Table 2 Discharge Prediction of Three Methods
方法 RMSE/（m3·s-1） RSD/% 潮流调和分析 11 819 19.90 BP神经网络* 8 427 14.18 本文组合方法* 6 786 11.42 注：*表示模型运行100次后取得的均值。
下载: 导出CSV
-
[1] 吴玲莉. 感潮河段设计水位方法确定与水位预报研究[D]. 南京: 河海大学, 2006. Wu Lingli. Study on Determination of Design Water Level Method and Water Level Forecast in Tidal Reach[D]. Nanjing: Hohai University, 2006.
[2] Xiong L H, O'Connor K M. Comparison of Four Updating Models for Real-Time River Flow Forecasting[J]. Hydrological Sciences Journal, 2002, 47(4): 621-639. doi: 10.1080/02626660209492964
[3] 陈志高. 长江口局域时空流场构建及实时流量精确估计[D]. 武汉: 武汉大学, 2015. Chen Zhigao. Construction of Local Spatiotemporal Flow Field and Accurate Estimation of Real-Time Flow in the Yangtze River Estuary[D]. Wuhan: Wuhan University, 2015.
[4] 方国洪. 潮汐和潮流的分析和预报[M]. 北京: 海洋出版社, 1986. Fang Guohong. Analysis and Forecast of Tides and Tidal Currents[M]. Beijing: Ocean Press, 1986.
[5] 唐岩, 暴景阳, 刘雁春, 等. 短期潮汐潮流数据的正交潮响应分析研究[J]. 武汉大学学报(信息科学版), 2010, 35(10): 1151-1156. https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201010004.htm Tang Yan, Bao Jingyang, Liu Yanchun, et al. Research on Orthogonal Response Analysis of Tide and Current Using Short Term Data[J]. Geomatics and Information Science of Wuhan University, 2010, 35(10): 1151-1156. https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201010004.htm
[6] Zhao J H, Chen Z G, Zhang H M, et al. Multiprofile Discharge Estimation in the Tidal Reach of Yangtze Estuary[J]. Journal of Hydraulic Engineering, 2016, 142(12): 04016056. doi: 10.1061/(ASCE)HY.1943-7900.0001201
[7] 戚德虎, 康继昌. BP神经网络的设计[J]. 计算机工程与设计, 1998, 19(2): 48-50. https://www.cnki.com.cn/Article/CJFDTOTAL-SJSJ202207010.htm Qi Dehu, Kang Jichang. On Design of the BP Neural Network[J]. Computer Engineering and Design, 1998, 19(2): 48-50. https://www.cnki.com.cn/Article/CJFDTOTAL-SJSJ202207010.htm
[8] Criswell J A., Lin E B. River Flow Forecasting Using An Inverse Wavelet Transform Neural Network Approach [J]. International Journal of Applied Mathematics, Computational Science and Systems Engineering, 2021, 3: 67-70.
[9] 要震, 许继平, 孔建磊, 等. 基于GA-Elman的河流水位预测方法研究[J]. 长江科学院院报, 2018, 35(9): 34-37. https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB201809009.htm Yao Zhen, Xu Jiping, Kong Jianlei, et al. Prediction of River Water Level by GA-Elman Model[J]. Journal of Yangtze River Scientific Research Institute, 2018, 35(9): 34-37. https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB201809009.htm
[10] Hecht-Nielsen R. Theory of the Backpropagation Neural Network [M]. Amsterdam : Elsevier, 1992.
[11] 王炎, 王船海, 王妮, 等. BP神经网络在感潮河段水动力洪水预报中的应用[J]. 水力发电, 2016, 42(2): 21-25. https://www.cnki.com.cn/Article/CJFDTOTAL-SLFD201602013.htm Wang Yan, Wang Chuanhai, Wang Ni, et al. Application of BP Neural Network in Tidal River Hydrodynamic Flood Forecasting[J]. Water Power, 2016, 42(2): 21-25. https://www.cnki.com.cn/Article/CJFDTOTAL-SLFD201602013.htm
[12] 秦思远, 李进军, 龙冰心, 等. 基于GPOS-BP神经网络模型的潮汐预报[J]. 海洋信息, 2020, 35(2): 1-5. https://www.cnki.com.cn/Article/CJFDTOTAL-HTXX202002001.htm Qin Siyuan, Li Jinjun, Long Bingxin, et al. Tide Forecast Model Based on GPOS-BP Neural Network[J]. Marine Information, 2020, 35(2): 1-5. https://www.cnki.com.cn/Article/CJFDTOTAL-HTXX202002001.htm
[13] 张峰, 王琪, 卢美, 等. 基于人工神经网络的海流预报研究[J]. 海洋预报, 2018, 35(4): 41-46. https://www.cnki.com.cn/Article/CJFDTOTAL-HYYB201804006.htm Zhang Feng, Wang Qi, Lu Mei, et al. Study of Tidal Current Prediction Based on Artificial Neural Network[J]. Marine Forecasts, 2018, 35(4): 41-46. https://www.cnki.com.cn/Article/CJFDTOTAL-HYYB201804006.htm
[14] Elganiny M A, Eldwer A E. Enhancing the Forecasting of Monthly Streamflow in the Main Key Stations of the River Nile Basin[J]. Water Resources, 2018, 45(5): 660-671.
[15] Hidayat H, Hoitink A J F, Sassi M G, et al. Prediction of Discharge in a Tidal River Using Artificial Neural Networks[J]. Journal of Hydrologic Engineering, 2014, 19(8): 04014006.
[16] Dragomiretskiy K, Zosso D. Variational Mode Decomposition[J]. IEEE Transactions on Signal Processing, 2014, 62(3): 531-544.
[17] Huang N E, Shen Z, Long S R, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis[J]. Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903-995.
[18] 万猛. 象山港余流特征及其影响因素研究[D]. 杭州: 浙江大学, 2015. Wan Meng. Study on Residual Current Characteristics and Its Influencing Factors in Xiangshan Port[D]. Hangzhou: Zhejiang University, 2015.
[19] 赵力学. 基于混合BP神经网络的河流水位流量预测方法研究[D]. 武汉: 武汉理工大学, 2019. Zhao Lixue. Study on Forecasting Method of River Water Level and Discharge Based on Hybrid BP Neural Network[D]. Wuhan: Wuhan University of Technology, 2019.
[20] 张红梅, 黄家勇, 赵建虎, 等. 基于梯度训练法的径向基函数潮流分离方法[J]. 武汉大学学报(信息科学版), 2017, 42(12): 1811-1817. https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201712020.htm Zhang Hongmei, Huang Jiayong, Zhao Jianhu, et al. An Improved Tidal Current Separation Method of Radial Basis Function Using Gradient Training[J]. Geomatics and Information Science of Wuhan University, 2017, 42(12): 1811-1817. https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201712020.htm
[21] 罗亦泳, 黄城, 张静影. 基于变分模态分解的变形监测数据去噪方法[J]. 武汉大学学报(信息科学版), 2020, 45(5): 784-790. https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH202005019.htm Luo Yiyong, Huang Cheng, Zhang Jingying. Denoising Method of Deformation Monitoring Data Based on Variational Mode Decomposition[J]. Geomatics and Information Science of Wuhan University, 2020, 45(5): 784-790. https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH202005019.htm
[22] 汪亚平, 潘少明, Wang H V, 等. 长江口水沙入海通量的观测与分析[J]. 地理学报, 2006, 61(1): 35-46. https://www.cnki.com.cn/Article/CJFDTOTAL-DLXB200601003.htm Wang Yaping, Pan Shaoming, Wang H V, et al. Measurements and Analysis of Water Discharges and Suspended Sediment Fluxes in Changjiang Estuary[J]. Acta Geographica Sinica, 2006, 61(1): 35-46. https://www.cnki.com.cn/Article/CJFDTOTAL-DLXB200601003.htm
[23] 沈花玉, 王兆霞, 高成耀, 等. BP神经网络隐含层单元数的确定[J]. 天津理工大学学报, 2008, 24(5): 13-15. https://www.cnki.com.cn/Article/CJFDTOTAL-TEAR200805006.htm Shen Huayu, Wang Zhaoxia, Gao Chengyao, et al. Determining the number of BP neural network hidden layer units[J]. Journal of TianJin University of Technology, 2008, 24(5): 13-15. https://www.cnki.com.cn/Article/CJFDTOTAL-TEAR200805006.htm
-
期刊类型引用(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
