宏观人群恐慌传播动力学模型建立方法及动力学模型应用与流程

本发明涉及人群疏散仿真技术领域，尤其是涉及一种宏观人群恐慌传播动力学模型建立方法及动力学模型应用。

背景技术：

疏散人群的群体行为中最具灾难性的是由恐慌所致的人群踩踏，其踩踏和过度拥挤常常导致大量人员伤亡。为此研究者们进行大量的疏散模拟演练，希望得出恐慌传播特点，然而模拟演练都是预先组织事件，不能完全反映灾难现场的疏散人群的恐慌情绪，加上恐慌传播具有复杂性、时变性、无法复制性。因此，有必要建立描述恐慌传播的动力学模型，研究恐慌传播的特点，从而为防止过度拥挤和踩踏提供方法支持。此外，经典的社会力模型和helbingd提出的心理和行为恐慌波动模型都为恐慌测量提供了基础性方法。但是迄今为止，大多数研究人员在微观模型中将恐慌程度映射为疏散个体速度的变化；在宏观模型中仅通过对动物和昆虫的模拟实验来研究恐慌，未涉及恐慌程度的定量分析。然而，宏观模型对于大规模人群疏散问题更具有说服力。因此，需要对宏观模型开展进一步研究。

技术实现要素：

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种宏观人群恐慌传播动力学模型建立方法及动力学模型应用。

本发明的目的可以通过以下技术方案来实现：

一种宏观人群恐慌传播动力学模型建立方法，包括以下步骤：

1)将疏散空间划分为多个离散化的单位网格，基于信息熵理论构建每一单位网格的恐慌熵；

2)构建各相邻单位网格的恐慌熵值梯度；

3)基于所述恐慌熵和恐慌熵值梯度构建宏观人群恐慌传播动力学模型，具体表示为：

其中，ex，y(k+1)是速度梯度比，ξ是0和1之间的随机系数，是时间步长k时网格(x，y)与相邻网格(a，b)的恐慌熵值的梯度，∑grad(k)为整体恐慌熵梯度值，vx，y(k)是时间步长为k时网格(x，y)中的速度，ex，y(k)是时间步长为k时网格(x，y)中的恐慌熵值，σ1和σ2是可调节常数，e1和e2是恐慌熵特征点，dir(ex，y(k+1))、||ex，y(k+1)||是时间步长为k+1时整个疏散空间中的恐慌熵方向和恐慌熵大小。

进一步地，所述每一单位网格中，将疏散人群的速度方向离散化为八个等间距的主方向。

进一步地，构建所述恐慌熵时，考虑在网格(x，y)中疏散人群的恐慌熵方向和恐慌熵大小

疏散人群的恐慌熵方向表示为：

疏散人群的恐慌熵大小表示为：

其中，n是每个单位网格中疏散个体的总数，nc是每个方向间隔中疏散个体的总数，md是每个离散速度区间中疏散个体的总数，vmax为疏散个体最大速度。

进一步地，所述恐慌熵值梯度表示为：

其中，ea，b(k)是时间步长为k时与网格(x，y)相邻的网格(a，b)中的恐慌熵值，a＝x-1，x，x+1，b＝y-1，y，y+1。

进一步地，所述整体恐慌熵梯度值表示为：

本发明还提供一种基于所述的宏观人群恐慌传播动力学模型的人群疏散仿真方法，包括以下步骤：

a)设置疏散空间以及疏散个体的初始位置和初始速度；

b)计算每个网格的恐慌熵和网格中人群的移动方向，实时更新疏散人员的位置信息；

c)循环利用所述宏观人群恐慌传播动力学模型计算每一时间步长对应的疏散人群速度和恐慌熵值，直至仿真结束。

进一步地，所述步骤b)中，根据人群的速度梯度比确定所述移动方向。

进一步地，该方法还包括：

d)将对恐慌传播过程中人群恐慌态势演化过程以3d空间展示。

与现有技术相比，本发明基于宏观人群疏散模型，重点研究疏散过程中人群恐慌因素，提出恐慌熵概念，实时计算人群恐慌程度，进而提出恐慌传播动力学模型，经过人群疏散状态参数初始化后，实时计算宏观人群恐慌程度，利用恐慌传播特性自适应获取下一时刻的人群恐慌程度分布，动态地分析恐慌对人群疏散的影响，本发明具有以如下有益效果：

(1)恐慌传播定量分析与恐慌态势演化展示

本发明基于二维空间的人群疏散aw-rascle模型，在恐慌传播定量分析中，提出相邻网格熵值梯度计算公式，同时计算出整体的熵值梯度值。利用离散网格熵值梯度值与整体熵值梯度值之比，提出人群速度梯度比计算公式，获得速度方向上的最大梯度比，利用速度梯度比定量反应疏散人群k+1时刻的人群疏散状态，基于恐慌熵值梯度以及k时刻人群速度，提出k+1时刻的恐慌熵大小与方向的计算公式，实现了恐慌传播过程的定量分析，提高仿真可靠性。

由于恐慌熵值直接反应人群恐慌程度，因此本发明实现了对恐慌传播过程中人群恐慌态势演化过程的3d展示，使得仿真结果更加切实直观。

(2)恐慌传播动力学仿真迭代方法

本发明与现有技术相比，最大的优势在于利用动力学传播特性实现人群恐慌传播特性分析。恐慌熵值的计算离不开人群速度，所以本发明利用k时刻的人群速度，结合恐慌大小对速度的影响特点，提出k+1时刻速度的公式。另外在恐慌熵方向的迭代方法中，首先提出k时刻的恐慌熵值梯度，由于恐慌熵值直接反应人群恐慌程度，恐慌程度又与人群密度与速度有关，所以恐慌熵值梯度最大时，其相邻两网格的人群密度与速度差最大，可以确定人群移动方向，并通过迭代计算实现实时更新疏散空间中的恐慌熵值，。

附图说明

图1为t型路口行人流速度方向重叠情况；

图2为疏散个体速度方向划分和离散化网格的恐慌熵，其中，(a)为速度的方向示意图，(b)为离散化网格的恐慌熵图；

图3为时步为70时恐慌熵3d空间分布仿真图；

图4为时步为70时恐慌熵2d空间分布仿真图；

图5为t型街道交叉口疏散过程中密度和恐慌熵值随时间步长变化的趋势图；

图6为疏散人群运动速度与恐慌熵值的关系。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本发明提供一种宏观人群恐慌传播动力学模型建立方法，包括以下步骤：1)将疏散空间划分为多个离散化的单位网格，基于信息熵理论构建每一单位网格的恐慌熵；2)构建各相邻单位网格的恐慌熵值梯度；3)基于所述恐慌熵和恐慌熵值梯度构建宏观人群恐慌传播动力学模型。同时，本发明还基于上述宏观人群恐慌传播动力学模型实现人群恐慌传播，进而实现人群疏散模拟。

一、人群疏散宏观动力学模型

基于一维流动模型，a.aw和m.rascle提出了二维空间的人群疏散模型(以下简称为aw-rascle模型)。经典质量守恒方程由偏微分方程的守恒确定，如公式(1)所示。

其中ρ和v分别是行人密度和速度。是时间t的偏导数，是距离x的偏导数。假设ρ和v是相互独立的，压力项p由公式计算，c0是密度的预期响应系数。用于反映ρ和v关系的偏微分方程，如公式(2)所示。

其中p(ρ)是压力项，v是水平速度。人群一维动态aw-rascle模型由公式(1)和公式(2)的两个部分双曲偏微方程(pde)组成。利用公式(3)、公式(4)和公式(5)将一维aw-rascle模型转换为二维aw-rascle模型。

其中v是水平速度，u是垂直速度。基于人群疏散aw-rascle模型的t型街道人群疏散动态模型增加了矢量双向叠加的影响，垂直方向为u，矢量水平方向为v。s1和s2是松弛项目，它们使速度和密度保持平衡。ph和pv分别是考虑交叉区域的影响矩阵的水平压力项和垂直压力项，它使交叉口的行人密度分布更加合理。为了解决双向流动叠加问题，引入交叉区域的影响矩阵mimp，如公式(6)所示。c是最大影响系数，i和j是交叉区域的疏散个体的网格坐标。基于人群疏散的一维模型，构建了如图1所示的t型路口的二维流动模型。交叉区域的影响矩阵mimp如公式(6)所示。

其中i和j是交叉区域的水平和垂直坐标，是交叉区域(i，j)坐标处的影响因子。

另外，行人流进入交叉区域，道路环境对人流的影响是不同的。当人们第一次进入交叉口区域时，靠近墙壁的行人流量影响也不同。是主干路和支路上速度的偏分后的影响矩阵因子。t型街道恐慌人群疏散的宏观模型如公式(7)、公式(8)和公式(9)所示。

ph(ρ，v，u)和pv(ρ，v，u)分别是考虑交叉区域的影响矩阵后的水平压力项和垂直压力项。

是在水平方向的速度偏分后的影响因子，是在垂直方向的速度偏分后的影响因子。β和γ是用于将密度转换为压力的功率因数。

二、恐慌传播动力学模型

(1)初始化恐慌熵的定义

本发明利用恐慌熵来定量分析个体行为和整个群体混乱程度之间的关系。恐慌熵是用来度量疏散人群恐慌程度的量，即疏散人群恐慌程度越大恐慌熵值越大，恐慌熵定义如下：

假设u是一个任意的疏散个体，其属于集合g＝{u1，u2，…}，该集合代表一个单位网格(x，y)，(x，y)为该网格中心点坐标，其中的元素有限可数。设概率值相应局部区域中单位网格空间(x，y)的恐慌熵定义为：

其中ui表示疏散个体，表示速度分布概率。如图2所示，x为水平坐标值，y为垂直坐标值，单位网格中的人群看作一个整体，疏散个体的行为特征由单位网格人群速度的大小和方向来表示。图2中，将每个网格中的每个疏散个体的速度方向划分为8个方向，

在网格(x，y)中疏散人群的恐慌熵方向dir(ex，y)和恐慌熵大小||ex，y||的计算公式如下：

其中n是每个离散网格中疏散个体的总数，nc是每个方向间隔中疏散个体的总数，md是每个离散速度区间中疏散个体的总数。因此，可以通过疏散空间中人群速度的大小和方向的分布来计算恐慌熵的初始值。

在初始情况的零时刻，离散网格(x，y)中某一e方向(总共8个方向)的疏散群体的速度为：

γ是以e方向为基准的-22.5°和+22.5°之间的扇形区域，是γ方向的速度值。根据疏散人群的初始速度分布，同理初始化每个网格的疏散速度。

(2)恐慌传播模型的建立

根据恐慌熵的公式，可以获得时间步长k时离散网格(x，y)与相邻网格(a，b)的恐慌熵值的梯度，如公式(16)所示。

因此整体的恐慌熵梯度值如公式(17)所示：

∑grad(k)＝∑a，b(ex，y(k)-ea，b))(17)

相邻网格(a，b)为位于离散网格(x，y)周围的8个方向对应的8个网格，其中a＝x-1，x，x+1和b＝y-1，y，y+1，离散空间的中心点的位置是(x，y)。在时间步长为k时，离散网格(x，y)的恐慌熵是ex，y(k)；离散网格(a，b)的恐慌熵是ea，b(k)。

在恐慌熵值梯度公式的基础上，提出疏散人群速度的梯度比，即网格(x，y)与网格(a，b)的恐慌熵值梯度与整体恐慌熵值梯度的比值。在时间步长为k+1时，人群的速度梯度比如下：

其中ξ是0和1之间的随机系数。在梯度比中使用随机系数ξ来代表速度方向的随机性，获得速度方向上的最大梯度比，以及与最大比率对应的方向，则此方向的行人最有可能进入下一个网格，从而确定人群速度方向。

在时间步长k+1处离散化网格(x，y)中的人群的速度大小如下：

其中vx，y(k)是时间步长为k时网格(x，y)中的速度，ex，y(k)是时间步长为k时网格(x，y)中的恐慌熵值。σ1和σ2是可调节常数，e1和e2是恐慌熵的特征点。当离散网格中人群的恐慌熵小于e1时，恐慌对疏散人群的速度几乎没有影响。当离散网格中的人群恐慌熵在e1和e2之间时，恐慌因素促进人群疏散，人群疏散速度增加并减少了疏散时间。当离散网格中的人群恐慌熵大于e2时，恐慌因素会减缓人群的疏散，不利于人群尽快撤离。在时间步长为k+1时整个疏散空间中的恐慌熵方向和恐慌熵大小由公式(20)和公式(21)给出。

其中a＝x-1，x，x+1和b＝y-1，y，y+1。离散空间的中心点的位置为(x，y)，并且在时间步长为k时离散网格(a，b)的速度为va，b(k)。因此，恐慌传播动力学模型由公式(18)、公式(19)、公式(20)和公式(21)组成。

假设两网格分别为(x，y)、(a，b)，则k+1时刻，由公式(18)可获得人群从网格(x，y)涌向网格(a，b)的人数最多，所以此方向则为该时刻恐慌熵值的方向，此推导过程符合动力学传播特性。利用公式(19)，(20)和(21)进行算法迭代，可以实现实时更新疏散空间中的恐慌熵值。

三、人群疏散仿真

在每个时间步骤中，恐慌传播模型必须计算恐慌熵并更新人群流速。本实施例以t型街道为例，仿真执行过程如下：

步骤1：在t型街道处设置疏散空间，并且设置两条街道的宽度和两条街道的交叉口。

第2步：设置撤离人员的初始位置和速度。假设所有撤离者随机分布在空间中，并且假设所有撤离者都具有初始速度。

步骤3：根据所有撤离人员的初始速度，计算每个网格的恐慌熵和网格中人群的移动方向。

步骤4：使用公式(18)计算时间步长为k时每个网格中疏散人群的速度和密度(速度梯度比)，并实时更新疏散人员的位置信息。

步骤5：根据公式(19)、(20)和(21)进行公式迭代，实时更新疏散空间中恐慌熵值。

步骤6：更新时间步长为k＝k+1，然后返回步骤4，循环步骤直到人群密度达到对应踩踏事件的最大值。

实施例

运用本发明的宏观人群恐慌传播动力学建模方法，以2015年麦加踩踏事件为背景，通过仿真重现发生在204号街道和223号街道交汇处的丁字路口踩踏事件。该踩踏事件发生在06时(格林尼治标准时间09时09分)，在204号街道和223号街道的丁字路口，朝圣者走向五层楼，围绕着被称为jamarat大桥的柱子。在两条相反方向的街道，大量人群聚集在丁字路口的交汇处，后面的疏散人群，不了解前方的挤塞情况，不断前进，导致了朝圣者之间的严重踩踏事件，该事件导致至少2177人死亡。

根据204号街道与223号街道交汇处的踩踏事件的实际参数，在街道上加载了2000名疏散人员作为初始条件，右侧壁上的出口宽度设为d，地图中红色标志的区域为实际人群蜂拥的位置，如图4所示。其中204号街道是主干路，223号街道是支路。204号街道的宽度设置为wmain＝10m，223号街道的宽度设置为wbranch＝9m，204街道的模拟长度设置为lmain＝500m，并设置223号街道的模拟长度为lbranch＝300m。街道204和街道223的宽度比为1.1∶1。人群的高恐慌值主要分布在初始位置和丁字路口的交汇处。初始化位置中人群的集中会导致恐慌熵值的增加，之后随着204号街道和233号街道的人群移动到交汇区域，每个网格的人群密度迅速增加，导致人群恐慌程度也迅速增加，人群密度与人群恐慌熵的关系如图5所示。

另外，在恐慌传播模型中，人群运动速度和人群恐慌程度是相互影响的。不同的运动速度可以影响下一个时间步长下人群的恐慌熵值，并且恐慌熵值的大小影响之后的人群运动速度。利用公式(13)和公式(14)已知人群疏散速度矢量计算恐慌熵值，接着再利用恐慌传播动力学模型的公式(19)计算出下一时间步长下的疏散人群速度，再通过公式(20)和公式(21)计算此时的恐慌熵值。循环此过程得出人群疏散速度与恐慌熵值关系图，如图6所示。综合图5和图6，恐慌传播的动力学特征以及疏散人群密度，速度和恐慌熵之间的关系是显而易见的。

通过matlabr2013b的程序仿真，得到t型街道口疏散人群恐慌熵分布情况。根据模拟结果，它可视地反映了人群的高恐慌熵值主要分布在初始位置和t型街道的交叉口。最初，初始化位置中由于人群的突然集中导致较大的恐慌熵，此时街道的其他位置的恐慌熵值较低。随着人群的恐慌传播，整个街道的恐慌熵值越来越高，直到模拟结束。

利用恐慌传播定量分析公式，加载t型街道口疏散人群数量及速度初值，采用matlabr2013b进行程序仿真，得到t型街道口疏散人群恐慌熵值变化情况，并且通过3d空间进行展示，如图3所示，3d空间中的恐慌熵分布在仿真时间步长为70时达到更高的恐慌程度。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。