Simulation of Pedestrian Evacuation Considering Cohesion of Social Groups on University Campus
-
摘要: 通过大学校园内通道及楼梯上的观测实验，获取行人群组运动的相关数据，基于Pathfinder软件构建群组疏散模型，研究在人员疏散过程中出现的瓶颈以及群组行为对疏散效率的影响。结果表明：瓶颈处密度与速度呈负相关；较大的群组比例和较强的群组凝聚力均会延长疏散时间，行人成群组运动也会增加疏散过程中的堵塞次数和拥堵持续时间，所以后期人员疏散研究中需要考虑群组行为对疏散结果的影响, 在没有其他行人需要帮助的情况下，鼓励行人尽量独自行走，减少对同伴的等待，同时在建筑物出口设计时需要注意群组对疏散效率的负面影响，提高设计规格。
-
关键词:
- 交通工程 /
- 疏散模拟 /
- Pathfinder /
- 群组
Abstract:Objectives Pedestrian and evacuation dynamics have attracted much attention in recent decades, specifically with respect to crowd risk analysis and pedestrian safety. In our daily life, individuals do not walk alone and they often walk with their friends or families because of the social relationship. However, the effect of group behavior observed in reality is rarely investigated in current studies.Methods In this paper, movement data of pedestrian social groups on corridors and stairs are obtained through observation experiments on the university campus. Evacuation model for social groups is built based on Pathfinder software to study the bottlenecks in the process of egress.Results The impact of group behavior on evacuation efficiency is investigated, and the results show that the density and speed in the bottleneck area are negatively correlated. More groups and stronger cohesion among group members will significantly extend the evacuation time. Besides, jam time increases due to the existence of social groups.Conclusions The group behavior must be considered during the evacuation process in the future study. In the absence of other pedestrians needing help, pedestrians are encouraged to walk independently as much as possible to reduce the waiting times for others. At the same time, the negative effect of social groups on evacuation efficiency should be considered in the design of pedestrian facilities or public buildings.-
Keywords:
- traffic engineering /
- evacuation simulation /
- Pathfinder /
- pedestrian group
-
整体最小二乘(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算法虽然总体表现良好，但随着粗差个数的增加，对于个别仿真数据抗差失效，而本文算法仍表现出了良好的抗差性能，总体抗差效果更稳定。
-
![]()
图 7 有群组和无群组情况下不同时刻疏散模拟对比
Figure 7. Comparison of Evacuation Simulation with and Without Groups at Different Time
![]()
图 13 有群组和无群组时瓶颈处速度及密度对比
Figure 13. Comparison of Speed and Density at the Bottleneck with and Without Groups
等级 人均面积/(m2·人-1） 概述 A >3.3 自由流，通行方便，可以避免碰撞 B (2.3, 3.3] 小范围碰撞，通行和速度受到限制 C (1.4, 2.3] 较拥挤但移动顺畅，通行受到限制，交叉或反向交会困难 D (0.9, 1.4] 碰撞明显，速度受到影响 E (0.5, 0.9] 交叉困难，断断续续的停滞 F ≤0.5 临界密度，低流率，频繁停滞，更多的身体接触 
下载: 导出CSV
表 2 有群组与无群组时人员疏散参数对比/s
Table 2 Comparison of Evacuation Parameters with and Without Groups/s
有无群组 最大连续拥堵时间 累计拥堵时间 运动时间 无群组 2.48 11.63 78.21 有群组 3.15 15.61 91.44
下载: 导出CSV
-
[1] Zhang J, Cao S C, Salden D, et al. Homogeneity and Activeness of Crowd on Aged Pedestrian Dynamics[J]. Procedia Computer Science, 2016, 83: 361-368 doi: 10.1016/j.procs.2016.04.137
[2] 杨雪, 陈立, 田欢欢, 等. 单向和双向行人流经过通道瓶颈的实验[J]. 上海大学学报(自然科学版), 2015, 21(3): 356-363 https://www.cnki.com.cn/Article/CJFDTOTAL-SDXZ201503010.htm Yang Xue, Chen Li, Tian Huanhuan, et al. Experiments of Unidirectional and Bidirectional Pedestrian Flows Through a Bottleneck in a Channel[J]. Journal of Shanghai University (Natural Science Edition), 2015, 21(3): 356-363 https://www.cnki.com.cn/Article/CJFDTOTAL-SDXZ201503010.htm
[3] Ma Y P, Li L H, Zhang H, et al. Experimental Study on Small Group Behavior and Crowd Dynamics in a Tall Office Building Evacuation[J]. Physica A: Statistical Mechanics and Its Applications, 2017, 473: 488-500 doi: 10.1016/j.physa.2017.01.032
[4] Cao S C, Seyfried A, Zhang J, et al. Fundamental Diagrams for Multidirectional Pedestrian Flows[J]. Journal of Statistical Mechanics: Theory and Experiment, 2017, 2 017(3): 033404 doi: 10.1088/1742-5468/aa620d
[5] 张鑫龙, 陈秀万, 李怀瑜, 等. 一种改进元胞自动机的人员疏散模型[J]. 武汉大学学报·信息科学版, 2017, 42(9): 1 330-1 336 doi: 10.13203/j.whugis20150763 Zhang Xinlong, Chen Xiuwan, Li Huaiyu, et al. An Improved Cellular Automata Model for Simulating Pedestrian Evacuation[J]. Geomatics and Information Science of Wuhan University, 2017, 42(9): 1 330-1 336 doi: 10.13203/j.whugis20150763
[6] Kallianiotis A, Papakonstantinou D, Arvelaki V, et al. Evaluation of Evacuation Methods in Underground Metro Stations[J]. International Journal of Disaster Risk Reduction, 2018, 31: 526-534 doi: 10.1016/j.ijdrr.2018.06.009
[7] Cao S C, Fu L B, Song W G. Exit Selection and Pedestrian Movement in a Room with Two Exits Under Fire Emergency[J]. Applied Mathematics and Computation, 2018, 332: 136-147 doi: 10.1016/j.amc.2018.03.048
[8] Cao S C, Wang P, Yao M, et al. Dynamic Analysis of Pedestrian Movement in Single-File Experiment Under Limited Visibility[J]. Communications in Nonlinear Science and Numerical Simulation, 2019, 69: 329-342 doi: 10.1016/j.cnsns.2018.10.007
[9] Wei X G, Lü W, Song W G, et al. Survey Study and Experimental Investigation on the Local Behavior of Pedestrian Groups[J]. Complexity, 2015, 20(6): 87-97 doi: 10.1002/cplx.21633
[10] 张泽天, 黄升, 屈楷博, 等. 基于Pathfinder的人员疏散中小群体行为研究[C]. 第30届全国高校安全科学与工程学术年会, 合肥, 中国, 2018 Zhang Zetian, Huang Sheng, Qu Kaibo, et al. Study on Small Group Behavior in Human Evacuation Based on Pathfinder[C]. The 30th National Confe-rence on Safety Science and Engineering in Universities, Heifei, China, 2018
[11] Thornton C, O'Konski R, Hardeman B, et al. Pathfinder: An Agent-Based Egress Simulator[C]//Pedestrian and Evacuation Dynamics, Boston: Springer, 2011
[12] Lam T N. Pedestrian Planning and Design by John J. Fruin[J]. Transportation Science, 1972, 6(2): 214-215 doi: 10.1287/trsc.6.2.214
-
期刊类型引用(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)
计量
- 文章访问数: 842
- HTML全文浏览量: 150
- PDF下载量: 71
- 被引次数: 11
