Multi-source venue flow anomaly detection method and system based on spatial-temporal feature fusion

Abstract

The invention provides a multi-source venue flow anomaly detection method and system based on spatio-temporal feature fusion, and the method comprises the following steps: constructing a tetrahedral mesh space of a venue in combination with a CIM and an indoor model, and mapping audience individuals in multi-source monitoring data into a dynamic trajectory particle swarm in a mesh; a real-time crowd rheological viscosity field is constructed based on a shear thickening fluid model by calculating the velocity gradient and the strain rate of a particle swarm. The viscosity field and the particle velocity field are combined, a dynamic abnormal area is recognized by solving a fluid vorticity evolution equation, a path failure risk area is recognized by calculating virtual shear stress of a channel boundary, and finally the two risk areas are fused to generate a final flow anomaly detection result. According to the invention, the disordered congestion phenomenon in the high-density crowd evacuation process can be accurately identified.

Description

Technical Field

[0001] This invention belongs to the field of anomaly detection technology, specifically relating to a multi-source venue traffic anomaly detection method and system based on spatiotemporal feature fusion. Background Technology

[0002] Crowd monitoring in large venues during major events is of paramount importance for public safety management. Modern security systems can use multi-source surveillance cameras to achieve real-time counting of crowds and visualization of spatial distribution of crowd density. By setting warning thresholds, such as when the number of people per unit area in a certain area exceeds a critical value, an alarm will be automatically triggered to alert management. However, this crowd density-based detection method cannot effectively distinguish between two distinct crowd states when facing large-scale crowd evacuation scenarios, as the single density scalar field is insufficient.

[0003] Specifically, an extremely high crowd density could correspond to either a high-density but uniformly oriented, stable-speed, and low-internal-friction orderly evacuation flow—a state that, while crowded, carries relatively manageable risks—or a high-density but disorderly crowding state characterized by head-on collisions, stagnation, and pushing and shoving among individuals. The latter is highly prone to causing serious accidents such as crushing and stampedes. However, these two states may exhibit identical or extremely similar values ​​on the density map, making it impossible for traditional density threshold alarm systems to effectively distinguish between them. Summary of the Invention

[0004] This invention provides a method and system for detecting abnormal traffic flow in multi-source venues based on spatiotemporal feature fusion, in order to solve the above-mentioned technical problems.

[0005] In a first aspect, the present invention provides a method for detecting anomalies in multi-source venue traffic based on spatiotemporal feature fusion, the method comprising the following steps: Obtain the CIM data and indoor model of the target venue, combine the CIM data and indoor model of the venue and use the tetrahedral subdivision algorithm to construct a grid space model of the target venue containing multiple tetrahedral grid units; By accessing multi-source monitoring data within the target venue and aligning spatiotemporal coordinates, the individual audience members identified in the multi-source monitoring data are mapped onto a grid space model, generating a dynamic trajectory particle swarm containing location and velocity information. The velocity gradient and strain rate of a dynamic trajectory particle swarm within a tetrahedral mesh cell are calculated, and a nonlinear mapping relationship between strain rate and flow resistance is constructed based on a shear-thickening fluid model to generate a real-time crowd rheological viscosity field. By combining the rheological viscosity field of the crowd and the particle velocity field constructed based on the dynamic trajectory particle swarm, a fluid vortex evolution equation is constructed. By solving the fluid vortex evolution equation, dynamic anomaly regions in the grid space model are identified. The virtual shear stress of the dynamic trajectory particle swarm on the channel boundary of the target venue is calculated by combining the rheological viscosity field of the crowd and the structural conductivity properties of the tetrahedral mesh unit, and the path failure risk area in the mesh space model is identified based on the virtual shear stress. The flow anomaly detection results for the target venue are generated by integrating the dynamic anomaly region and the path failure risk region.

[0006] Optionally, the step of combining CIM data and the venue's interior model and using a tetrahedral partitioning algorithm to construct a grid spatial model of the target venue containing multiple tetrahedral grid cells includes the following steps: Analyze the geometric topology information in CIM data and venue interior models to identify the boundary surface equations of walls, columns, fences and stairs in the target venue; The effective access range of the target venue in three-dimensional space is defined based on the boundary surface equation, and the effective access range is divided into a set of unstructured tetrahedral mesh elements using the tetrahedral subdivision algorithm. Based on the centroid coordinates and volume parameters of each tetrahedral mesh element in the tetrahedral mesh element set, the spatial positional relationship between the tetrahedral mesh element and the boundary surface equation is calculated, and the impassable normal vectors existing in the tetrahedral mesh element are determined according to the spatial positional relationship. Based on the impassable normal vector, a structural conduction tensor in the local coordinate system is constructed as the structural conduction property of the tetrahedral mesh element. The conduction coefficient along the direction of the impassable normal vector in the structural conduction tensor is set to a preset stagnation value, and the conduction coefficient perpendicular to the plane of the impassable normal vector in the structural conduction tensor is set to a unit passability value. Laplace smoothing is applied to the structural conduction tensor of all tetrahedral mesh elements to form a continuously distributed mesh space model of the target venue.

[0007] Optionally, the process of accessing multi-source monitoring data within the target venue and mapping individual audience members identified in the multi-source monitoring data to a grid space model through spatiotemporal coordinate alignment to generate a dynamic trajectory particle swarm containing location and velocity information includes the following steps: Acquire video surveillance data containing audience pixel coordinates collected by surveillance cameras in the target venue, and crowd density heat map data containing the number of audience members in the area collected by people counting cameras in the target venue; Based on a preset time synchronization protocol, video surveillance data and crowd density thermal data are aligned to the same sampling time. The viewer's pixel coordinates are back-projected onto the grid space model using a perspective transformation matrix to generate a three-dimensional observation point representing the position of an individual viewer in the real world; The three-dimensional observation points are input into a preset Kalman filter. The state transition matrix is ​​used to predict the state of the audience individual at the next moment. The optimal estimated position and optimal estimated velocity vector of the audience individual are updated in combination with the current three-dimensional observation points. The entity particles representing the audience individual are generated in the grid space model. By combining video surveillance data and crowd density thermal data, target areas in the target venue with only crowd density thermal data are identified. Based on the macroscopic density distribution function within the target area, supplementary virtual particles are generated in the grid space model using the Monte Carlo sampling method. The virtual particles are given an initial velocity that conforms to the characteristics of the local flow field. By merging physical and virtual particles, a dynamic trajectory particle swarm containing audience flow information within the target venue is generated.

[0008] Optionally, the step of generating supplementary virtual particles in the grid space model based on the macroscopic density distribution function within the target region and using the Monte Carlo sampling method, and assigning initial velocities to the virtual particles that conform to the characteristics of the local flow field, includes the following steps: The target region is mapped to the mesh space model, and the total number of target particles in each tetrahedral mesh cell within the target region is calculated by integration. Generate a candidate set of random coordinate points that follow a uniform distribution within the tetrahedral geometry of each tetrahedral mesh cell; The structural conductivity property of the tetrahedral mesh element is used to validate the random coordinate point candidate set and invalid coordinate points that fall inside physical obstacles are removed. Based on the total number of target particles, a corresponding number of coordinate points are extracted from the verified candidate set of random coordinate points as the generation positions of virtual particles. Search for the velocity vectors of existing physical particles within the neighborhood of the tetrahedral mesh cell where the virtual particle is located, calculate the velocity vectors of the physical particles using inverse distance weighted interpolation, and assign the velocity vectors to the virtual particle as its initial velocity.

[0009] Optionally, the calculation of the velocity gradient and strain rate of the dynamic trajectory particle swarm within the tetrahedral mesh element, and the construction of a nonlinear mapping relationship between strain rate and flow resistance based on the shear-thickening fluid model to generate a real-time crowd rheological viscosity field includes the following steps: The velocity vector of the dynamic trajectory particle swarm within the tetrahedral mesh element is mapped to the vertex nodes of the tetrahedral mesh element using shape function interpolation to obtain the nodal velocity field. A velocity gradient tensor describing local deformation characteristics is constructed by calculating the spatial partial derivatives of the nodal velocity field, and the symmetric part of the velocity gradient tensor is extracted as the strain rate tensor of the dynamic trajectory particle swarm. The shear rate amplitude of crowd flow within a tetrahedral mesh cell is quantified by calculating the second invariant of the strain rate tensor. The apparent viscosity of the tetrahedral mesh element is calculated by substituting the shear rate amplitude into a preset shear-thickening fluid model that includes the rheological index and consistency coefficient. After traversing all tetrahedral mesh elements and calculating the corresponding apparent viscosity values, and performing spatial smoothing, the rheological viscosity field of the target venue is generated.

[0010] Optionally, the step of substituting the shear rate amplitude into a preset shear-thickening fluid model that includes the rheological index and consistency coefficient to calculate the apparent viscosity value of the tetrahedral mesh element includes the following steps: Collect environmental decibel data within the target venue and calculate the acceleration variance of the dynamic trajectory particle swarm; A panic index model based on environmental decibel data and acceleration variance was constructed, and a positive correlation mapping function between the panic index and the rheological index output by the panic index model was established. A modified shear-thickening fluid equation is constructed using the updated rheological index and the preset consistency coefficient; The shear rate amplitude is input into the shear thickening fluid equation to calculate the apparent viscosity of the tetrahedral mesh element.

[0011] Optionally, the step of constructing a fluid vortex evolution equation by combining the rheological viscosity field of the crowd and the particle velocity field based on dynamic trajectory particle swarm, and identifying the dynamic anomaly region in the grid space model by solving the fluid vortex evolution equation includes the following steps: A particle velocity field is constructed based on a dynamic trajectory particle swarm. The curl operator is applied to the particle velocity field to calculate the vorticity vector of each tetrahedral mesh element. The fluid vorticity evolution equation is constructed based on the time change rate of the vorticity vector. Decompose the fluid vorticity evolution equation and calculate the vortex stretching term caused by the velocity gradient and the viscous diffusion term caused by the rheological viscosity field of the crowd, respectively. When the magnitude of the vortex stretching term is greater than that of the viscous diffusion term, the flow field is determined to be in an unstable growth state of vorticity accumulation. By integrating the square of the modulus of the vorticity vector within the tetrahedral grid cell by volume, the pseudo-vortex energy characterizing the rotational intensity of the local flow field is calculated. The tetrahedral mesh element set whose pseudo-vortex energy exceeds the preset growth threshold and is in an unstable growth state is marked as a dynamic anomaly region.

[0012] Optionally, the step of calculating the virtual shear stress of the dynamic trajectory particle swarm on the channel boundary of the target venue by combining the rheological viscosity field of the crowd and the structural conductivity properties of the tetrahedral mesh elements, and identifying the path failure risk area in the mesh space model based on the virtual shear stress, includes the following steps: In a mesh space model, a mesh whose structural conductivity properties of tetrahedral mesh elements change abruptly is defined as a channel boundary layer. The virtual shear stress of the passenger flow on the tangential plane of the channel boundary layer is calculated by combining the apparent viscosity and strain rate tensor at the channel boundary layer. The erosion resistance coefficient is set according to the physical width and geometry of the channel boundary layer; Based on dynamic trajectory particle swarm optimization and combined with virtual shear stress and erosion resistance coefficient, the attenuation rate of the channel's conduction function corresponding to the channel boundary layer is calculated. The grid area where the channel boundary layer is located, where the rate of attenuation of the dredging function exceeds the preset failure threshold, is identified as a path failure risk area.

[0013] In a second aspect, the present invention also provides a multi-source venue traffic anomaly detection system based on spatiotemporal feature fusion, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion as described in any one of the first aspects.

[0014] Thirdly, the present invention also provides a computer-readable storage medium storing instructions, characterized in that, when executed by a processor, the instructions cause the processor to be configured to perform the multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to any one of the first aspects.

[0015] The beneficial effects of this invention are: This invention constructs a tetrahedral mesh spatial model with physical constraints and maps individuals from multi-source monitoring data as dynamic trajectory particle swarms. It utilizes a nonlinear crowd rheological viscosity field based on a shear-thickening fluid model to characterize the dynamic changes in flow resistance at different strain rates, effectively distinguishing between benign high-density ordered flow and dangerous stagnation, collision, and other high internal friction states. Next, it captures dynamic anomalies caused by crowd collision and swirling by solving the fluid vortex evolution equation. Simultaneously, it assesses the path failure risk areas caused by overcrowding on boundary structures such as passageways and fences by calculating virtual shear stress. Finally, by fusing the two anomaly detection results, it accurately identifies disorderly crowding in large venues and issues timely warnings, thereby preventing stampedes caused by disorderly crowding. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating a multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion in one embodiment of this application.

[0017] Figure 2 This is a schematic diagram illustrating the process of constructing a grid space model in one embodiment of this application. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0019] The terms "first," "second," etc., used in the specification and claims of this application are used to distinguish similar objects and not to describe a specific order or sequence. It should be understood that such use of data can be interchanged where appropriate so that embodiments of this application can be implemented in orders other than those illustrated or described herein, and the objects distinguished by "first," "second," etc., are generally of the same class and the number of objects is not limited; for example, a first object can be one or more. Furthermore, in the specification and claims, "and / or" indicates at least one of the connected objects, and the character " / " generally indicates that the preceding and following objects are in an "or" relationship.

[0020] Figure 1 This is a flowchart illustrating a multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion in one embodiment. It should be understood that, although... Figure 1 The steps in the flowchart are shown sequentially as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified herein, there is no strict order in which these steps are executed, and they can be performed in other orders. Figure 1 At least some steps in the process may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps. For example Figure 1 As shown, the multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion disclosed in this invention specifically includes the following steps: S101. Obtain the CIM data and indoor model of the target venue, combine the CIM data and indoor model of the venue and use the tetrahedral subdivision algorithm to construct a grid space model of the target venue containing multiple tetrahedral grid units.

[0021] The process begins by analyzing the geometric parameters in the City Information Model (CIM) data and the refined interior model of the venue. This allows for the extraction of boundary surface equations defining key elements of the venue's physical structure, such as walls, columns, fences, and stairs. These equations, composed of coordinate variables in three-dimensional space and corresponding normal vector coefficients, precisely delineate the outlines of physical obstacles. Based on these boundary surface equations, the effective passageway for pedestrians is defined in a three-dimensional Cartesian coordinate system. An unstructured tetrahedral partitioning algorithm is then introduced to discretize this irregular passageway into a massive and tightly connected set of tetrahedral mesh cells. This mesh form perfectly conforms to the complex and varied architectural structure within the venue. For each tetrahedral mesh cell in the set, its centroid coordinates and volume parameters are calculated, thereby quantifying the spatial relationship between the mesh cell and the boundary surface equations and identifying the directions of inaccessible normal vectors within the mesh due to physical obstructions. Based on this direction, a structural conduction tensor is constructed in a local coordinate system. This tensor serves as a mathematical matrix describing the anisotropic conduction capability of space. The conduction coefficient along the impassable normal vector direction is forcibly set to a near-zero blocking value, while the conduction coefficient perpendicular to the plane of this normal vector is assigned a unit passability value. This mathematically characterizes the obstructive effect of walls and other obstacles on pedestrian flow. Laplace smoothing is performed on the structural conduction tensor of all tetrahedral mesh elements across the entire venue to eliminate parameter abrupt changes between adjacent meshes, forming a continuous, smooth target venue mesh space model with physical constraints.

[0022] S102. Access multi-source monitoring data within the target venue, and map the individual audience members identified in the multi-source monitoring data to the grid space model through spatiotemporal coordinate alignment, generating a dynamic trajectory particle swarm containing position and velocity information.

[0023] The process involves acquiring video surveillance data covering key areas of the venue and crowd density thermal data collected by people counting cameras. Network time protocols are used to strictly align the timestamps of both types of data to the same sampling time, ensuring data consistency across time. For the video data, a pre-calibrated perspective transformation matrix is ​​applied to project the two-dimensional viewer pixel coordinates on the image plane into a three-dimensional mesh space model, calculating the three-dimensional observation points of the viewers in the real physical space. To eliminate observation noise and smooth the trajectory, the sequence of three-dimensional observation points is input into a Kalman filter. A state transition matrix is ​​used to predict the current state based on the optimal estimate from the previous moment, and combined with the current observation values ​​to update the optimal estimated position and velocity vector of each viewer in three-dimensional space, thereby generating physically-possessing entity particles in the mesh model. For blind spots covered only by thermal data of population density, a Monte Carlo sampling method is used for probabilistic completion based on the macroscopic density distribution function within the area. Specifically, uniformly distributed coordinate points are randomly generated within the corresponding tetrahedral grid cells. Invalid points falling inside obstacles are removed after structural continuity property verification, and valid points are retained as virtual particles. An inverse distance weighted interpolation method is then used to assign initial velocities to the virtual particles, conforming to the characteristics of the local flow field, referencing the velocity vectors of neighboring physical particles. Finally, the physical and virtual particles are merged to form a dynamic trajectory particle swarm containing precise position and velocity information.

[0024] S103. Calculate the velocity gradient and strain rate of the dynamic trajectory particle swarm within the tetrahedral mesh element, and construct the nonlinear mapping relationship between strain rate and flow resistance based on the shear thickening fluid model to generate a real-time crowd rheological viscosity field.

[0025] In this process, shape function interpolation is used to map the velocity vectors of discrete particles within a tetrahedral mesh cell to the four vertex nodes of the mesh cell, forming a continuously distributed nodal velocity field. Spatial partial derivatives of this velocity field are calculated to construct a velocity gradient tensor describing the local fluid deformation characteristics. The symmetric portion of this tensor is extracted as a strain rate tensor, which physically describes the rate of expansion, contraction, and shear deformation of the crowd fluid micro-clusters. The second invariant of the strain rate tensor is then calculated to quantify the amplitude of the shear rate of the crowd flow within the tetrahedral mesh cell, reflecting the intensity of the relative motion between the crowd. Based on a shear-thickening fluid model in non-Newtonian fluid mechanics, a nonlinear mapping relationship between strain rate and flow resistance is established. Considering that a crowd in a panic state may exhibit a congestion-solidification effect similar to that of a non-Newtonian fluid, environmental decibel data within the venue is collected to calculate the particle acceleration variance to construct a panic index. A positive correlation function between the panic index and the rheological index is then established. The apparent viscosity is calculated using a constitutive equation constructed using the modified rheological index and consistency coefficient. The apparent viscosity values ​​of all grid cells are calculated and spatially smoothed. The resulting viscosity field can dynamically characterize the flow resistance of crowds under different densities and psychological states. High viscosity regions indicate extremely poor crowd mobility and a high risk of crushing.

[0026] S104. By combining the rheological viscosity field of the crowd and the particle velocity field constructed based on the dynamic trajectory particle swarm, the fluid vorticity evolution equation is constructed. By solving the fluid vorticity evolution equation, the dynamic anomaly region in the grid space model is identified.

[0027] In this study, the Euler velocity field constructed based on dynamic trajectory particle swarm optimization is used to calculate the vorticity vector within each tetrahedral grid cell using the curl operator. This vector physically describes the rotational intensity and direction of the local fluid micro-particles. Based on the vorticity transport principle in fluid mechanics, a fluid vorticity evolution equation is constructed, decomposing the time-varying rate of vorticity into a vortex stretching term caused by the velocity gradient and a viscous diffusion term caused by the rheological viscosity field of the crowd. The vortex stretching term reflects the nonlinear process in the flow field where uneven velocity distribution causes vortices to be elongated, leading to a sharp increase in rotational angular velocity, and is a key factor causing flow field instability. The viscous diffusion term reflects the dissipative effect of viscosity on vortices. By comparing the magnitudes of the two terms through numerical calculations, when the vortex stretching term is significantly larger than the viscous diffusion term, the local flow field is determined to be in an unstable growth state of vorticity accumulation, i.e., self-excited and enhanced disordered turbulence has appeared in the crowd. Furthermore, by volume integration of the square of the modulus of the vorticity vector within the tetrahedral grid cell, the pseudo-vortex energy, representing the total rotational kinetic energy of the local flow field, is calculated. By monitoring the temporal variation characteristics of pseudo-vortex energy, grid areas that are in an unstable growth state and whose pseudo-vortex energy values ​​exceed a preset safety threshold are marked as dynamically abnormal areas. Such areas usually indicate the precursors to high-risk dynamic events such as people falling or trampling due to collision or crowding.

[0028] S105. Calculate the virtual shear stress of the dynamic trajectory particle swarm on the channel boundary of the target venue by combining the rheological viscosity field of the crowd and the structural conductivity properties of the tetrahedral mesh unit, and identify the path failure risk area in the mesh space model based on the virtual shear stress.

[0029] In the scanning mesh space model, the set of mesh cells whose structural conductivity properties undergo drastic abrupt changes—that is, whose conductivity coefficient drops sharply from a unit value to a stagnation value—is defined as the channel boundary layer, representing the edges of walls or fences on both sides of the corridor or exit. Based on fluid dynamics boundary layer theory, and combined with the high apparent viscosity and strain rate tensor at the channel boundary layer, the virtual shear stress exerted by the flow of people on the tangential plane of the channel boundary layer is calculated. This stress physically simulates the tangential friction and compression caused by a high-density crowd moving against lateral walls or obstacles. An erosion resistance coefficient is set according to the physical width, geometry, and material properties of the channel boundary layer. This coefficient characterizes the upper limit of the channel's ability to maintain normal traffic flow. Utilizing the flow characteristics of dynamic trajectory particle swarms, combined with the calculated virtual shear stress and erosion resistance coefficient, a channel traffic flow attenuation model is constructed to calculate the rate of attenuation of the channel's traffic flow function. When the rate of decay of the channel boundary layer exceeds the preset failure threshold, it means that although the physical channel exists geometrically, it has actually lost its effective passage capacity due to excessive lateral pressure and friction from the crowd, and is very likely to form deadlock or arch effect blockage. Therefore, it is necessary to accurately identify the grid area where the channel boundary layer is located as the path failure risk area, indicating that the exit will fail due to physical compression.

[0030] S106. Generate the traffic anomaly detection results for the target venue by fusing dynamic anomaly areas and path failure risk areas.

[0031] The system considers dynamic anomaly areas as sources of instability within the population, representing the risk of conflict between people, such as pushing, falling, or turbulence. Path failure risk areas are viewed as structural bottlenecks under environmental constraints, representing the risk of conflict between people and the environment, such as exit blockage or passage congestion. Under the unified coordinate system of the grid spatial model, these two types of risk areas are spatially superimposed for analysis. If an area is simultaneously marked as both types of anomalies, it is identified as an extremely high-risk composite congestion point, requiring immediate triggering of the highest-level alarm. If only a single type of anomaly exists, it is graded and marked according to the anomaly's intensity index. The system output includes not only the three-dimensional spatial coordinate range of the anomaly area but also the anomaly's type attributes, formation mechanism, and suggested mitigation strategies. The final detection results are overlaid on the CIM visual interface in the form of a three-dimensional heat map or vector field, which intuitively displays potential dangerous hotspots in the venue. This allows security management personnel to clearly distinguish which areas are chaotic due to crowd panic and which areas are congested due to insufficient facility capacity, so as to take targeted diversion or expansion measures and realize the perception and early warning of complex traffic conditions in large venues.

[0032] In one embodiment, reference is made to Figure 2 The process of constructing a grid spatial model of the target venue containing multiple tetrahedral mesh elements by combining CIM data and the venue's interior model and using a tetrahedral subdivision algorithm includes the following steps: S201. Analyze the geometric topology information in the CIM data and the indoor model of the venue, and identify the boundary surface equations of the walls, columns, fences and stairs in the target venue; S202. The effective passage range of the target venue in three-dimensional space is defined based on the boundary surface equation, and the effective passage range is divided into an unstructured tetrahedral mesh unit set using the tetrahedral subdivision algorithm; S023. Based on the centroid coordinates and volume parameters of each tetrahedral mesh element in the tetrahedral mesh element set, calculate the spatial positional relationship between the tetrahedral mesh element and the boundary surface equation, and determine the impassable normal vectors existing in the tetrahedral mesh element according to the spatial positional relationship. S204. Based on the impassable normal vector, construct the structural conduction tensor in the local coordinate system as the structural conduction property of the tetrahedral mesh element. Set the conduction coefficient along the direction of the impassable normal vector in the structural conduction tensor to a preset stagnation value, and set the conduction coefficient perpendicular to the plane of the impassable normal vector in the structural conduction tensor to a unit passability value. S205. Perform Laplace smoothing on the structural conduction tensor of all tetrahedral mesh elements to form a continuously distributed mesh space model of the target venue.

[0033] In this embodiment, an IFC or CityGML format file containing the building structure is read, and the geometric data stream within the file is traversed. A feature recognition algorithm filters out decorative components, identifying physical entities within the venue that substantially obstruct crowd flow. These primarily include load-bearing walls, supporting columns, barriers, and stair treads. For each identified physical entity, the vertex coordinate sequence of its constituent surfaces is extracted, and the least squares method is used for geometric fitting of the spatial plane to calculate the boundary surface equations describing these physical interfaces. For any j-th obstruction plane, its boundary surface equation is established as follows: The form is given by x, y, z, where x, y, z are three-dimensional spatial coordinate variables. These constitute the normal vector components of the plane. This is the intercept parameter. Based on the boundary surface equation, the effective passage range for personnel is defined in a three-dimensional Cartesian coordinate system. Specifically, it is the overall bounding box space of the venue minus the space occupied by all obstacles defined by the boundary surface equation.

[0034] For effective passageways with irregular shapes and complex topologies, an unstructured tetrahedral partitioning algorithm is employed for spatial discretization. This algorithm, based on the Delaunay triangulation principle, fills the continuous passageway space with a massive and tightly connected set of tetrahedral mesh cells. Compared to structured meshes, tetrahedral mesh cells offer significant advantages in geometric adaptability, perfectly fitting the winding corridor boundaries, circular column surfaces, and irregular stair corners within the venue, avoiding jagged approximation errors. During the partitioning process, the algorithm dynamically adjusts the mesh density, automatically densifying the mesh in geometrically complex areas while maintaining moderate sparseness in open areas to balance computational accuracy and efficiency.

[0035] After mesh generation, all tetrahedrons in the set are traversed, and their geometric centroid coordinates and volume parameters are calculated based on the coordinates of their four vertices. These two parameters represent the spatial centroid location and physical scale of the mesh, respectively. Then, using the centroid coordinates as a reference point, the minimum Euclidean distance from the centroid to the plane defined by the equations of each boundary surface is calculated to quantify the proximity of the mesh cell to physical obstacles. When the calculated minimum distance is less than a critical threshold set based on the mesh volume parameters, the tetrahedral mesh cell is determined to be adjacent to a physical obstacle, and the normal vector of the nearest boundary surface is immediately extracted and defined as the impassable normal vector within that mesh cell. The impassable normal vector indicates which direction is absolutely forbidden at that grid location (i.e., the direction pointing inwards towards a wall or pillar). Next, based on the identified impassable normal vectors, a structural connectivity tensor is constructed in the local coordinate system for each tetrahedral grid cell. The structural connectivity tensor is a second-order symmetric matrix used to describe the ease of passage for people in different directions within space.

[0036] During the construction process, a local orthogonal coordinate system is first established with the impassable normal vector as one axis. Under this coordinate system, the structural conduction tensor exhibits significant anisotropic characteristics. Specifically, the main diagonal component of the tensor along the direction of the impassable normal vector, i.e., the conduction coefficient, is forcibly set to a preset hindrance value. This value is a positive number approaching zero, physically simulating the infinite resistance in the direction of the obstacle; simultaneously, the other two main diagonal components of the tensor located perpendicular to the impassable normal vector plane are set to a unit passability value of 1, representing that passage along the tangential direction of the wall is completely free and unimpeded. The mathematical construction of the structural passability tensor T can be expressed as follows: Where I is the unit tensor, This represents the dyadic product operation. Through tensor quantization, any velocity component attempting to penetrate the wall is naturally suppressed in numerical calculations, while the velocity component sliding along the wall is preserved, thus simulating the physical phenomenon of a crowd walking close to a wall.

[0037] The original structure conduction tensor field often exhibits drastic numerical jumps at obstacle boundaries. Therefore, a Laplace smoothing algorithm is used to iteratively update the structure conduction tensor of each grid cell. Specifically, for each central grid cell, all its topologically connected neighboring grid cells are searched, the weighted average of the neighboring cell tensors is calculated, and this average is mixed with the current tensor value of the central cell. The smoothing update formula can be described as follows: ,in and These are the tensors before and after the update, respectively. As a smoothing factor, Let M be the tensor of the neighborhood grid k, and M be the number of neighborhoods. The smoothed model not only ensures the stability of the flow field calculation numerically, but also physically simulates the expected behavior of a crowd gradually adjusting its direction of movement when approaching an obstacle, making the final generated flow anomaly detection results more consistent with the real movement patterns of crowds.

[0038] In one implementation, accessing multi-source monitoring data within the target venue and mapping individual audience members identified from the multi-source monitoring data to a grid space model through spatiotemporal coordinate alignment to generate a dynamic trajectory particle swarm containing location and velocity information includes the following steps: Acquire video surveillance data containing audience pixel coordinates collected by surveillance cameras in the target venue, and crowd density heat map data containing the number of audience members in the area collected by people counting cameras in the target venue; Based on a preset time synchronization protocol, video surveillance data and crowd density thermal data are aligned to the same sampling time. The viewer's pixel coordinates are back-projected onto the grid space model using a perspective transformation matrix to generate a three-dimensional observation point representing the position of an individual viewer in the real world; The three-dimensional observation points are input into a preset Kalman filter. The state transition matrix is ​​used to predict the state of the audience individual at the next moment. The optimal estimated position and optimal estimated velocity vector of the audience individual are updated in combination with the current three-dimensional observation points. The entity particles representing the audience individual are generated in the grid space model. By combining video surveillance data and crowd density thermal data, target areas in the target venue with only crowd density thermal data are identified. Based on the macroscopic density distribution function within the target area, supplementary virtual particles are generated in the grid space model using the Monte Carlo sampling method. The virtual particles are given an initial velocity that conforms to the characteristics of the local flow field. By merging physical and virtual particles, a dynamic trajectory particle swarm containing audience flow information within the target venue is generated.

[0039] In this embodiment, high-definition video surveillance cameras distributed at key nodes of the target venue are respectively connected to crowd density thermal cameras or infrared counters specifically used for counting the number of people in an area. For video surveillance data, target detection algorithms from computer vision are used to identify the bounding boxes of individual audience members in real time on the two-dimensional pixel plane of each frame image, and the center point at the bottom of the bounding box is extracted as the pixel coordinates (u,v) representing the audience member's position, thereby converting the continuous video stream into a discrete sequence of pixel positions. At the same time, the structured data uploaded by the crowd counting cameras is parsed. The structured data exists in the form of a gridded matrix or total area, reflecting the degree of crowd gathering in a specific physical area, i.e., crowd density thermal data. Since different types of monitoring devices are often connected in different sub-networks, and their data acquisition frequencies, clock references, and network transmission delays are significantly different, direct fusion will lead to serious spatiotemporal misalignment. Therefore, a unified global time reference axis must be established based on the network time protocol. Specifically, a time alignment module with buffering function can monitor the timestamp tags of various data packets in real time, and select a standard sampling time as a reference frame according to a preset time synchronization protocol. For video surveillance data, find the image frame with the closest timestamp; for crowd density thermal data, if the update frequency is low, use linear interpolation or zero-order hold to estimate the density value at the standard sampling time.

[0040] After time alignment, the pre-calibrated camera intrinsic matrix (containing focal length and principal point coordinates) and extrinsic matrix (containing the camera's rotation and translation vectors in the world coordinate system) are invoked to construct the inverse process of the projection model from 3D world space to the 2D image plane. Due to the lack of image depth information, it is usually assumed that the viewer's feet are on the ground, i.e., the ground plane equation Z=0 is used as an auxiliary constraint. A perspective transformation matrix is ​​constructed, and the inverse projection formula is used to map the pixel coordinates back to 3D space, generating 3D observation points representing the viewer's position in the real world. This point data falls directly within the pre-constructed tetrahedral mesh space model coordinate system, giving the pixels that originally only existed on the screen a concrete physical location in the digital twin space. The original 3D observation points often contain measurement errors caused by sensor noise, occlusion, or lighting changes, and only contain position information without velocity information. Therefore, a Kalman filter is also needed for state estimation and smoothing.

[0041] Specifically, each identified audience member is treated as an independent state system with its own state vector. The prior state at the next moment is predicted using the state transition matrix based on a uniform motion model. Then, the 3D observation point obtained through perspective transformation at the current moment is used as the observation vector. The Kalman gain is calculated in conjunction with the observation matrix, and the prior state is corrected using the observation residuals. Recursive iteration effectively filters out high-frequency jitter noise at the observation position, obtaining the optimal estimated position for each audience member. Furthermore, the smooth optimal estimated velocity vector is implicitly calculated through temporal difference and filtering of the position. Finally, the entity particles generated in the grid space model not only have accurate positions but also possess stable kinematic properties, realistically reflecting the audience's direction and speed of movement.

[0042] In large venues, blind spots inevitably exist where video surveillance cannot cover the area, but density statistical sensors are deployed. To ensure the continuity of the overall flow model, it is necessary to combine two types of data to identify target areas with only crowd density thermal data. For these target areas, a local macroscopic density distribution function is constructed based on the total number of people and thermal distribution map fed back by the sensors. Using Monte Carlo sampling, random sampling is performed within the grid space model volume corresponding to the target area based on the probabilistic characteristics of this density distribution function. Several candidate points are generated within each tetrahedral grid cell. Points falling inside obstacles are removed using a rejection sampling strategy, retaining valid points as supplementary virtual particles to ensure that the total number of generated particles is consistent with the number of people in the area. Since these areas lack individual trajectory observations and cannot directly measure velocity, it is necessary to refer to the flow field information of physical particles in adjacent areas and, using the continuity assumption in fluid mechanics, assign initial velocities to virtual particles that conform to the local flow field characteristics. Finally, the set of physical particles generated by Kalman filtering in the video coverage area is merged with the set of virtual particles generated by Monte Carlo sampling in the blind spots. During the merging process, spatial deduplication logic is executed to ensure that there are no double counts of particles in areas where the video and heatmap overlap, prioritizing the more accurate video data. The resulting dynamic trajectory particle swarm is a unified dataset containing information on all audience flow within the venue, where each particle carries its three-dimensional position coordinates, three-dimensional velocity vector, and attribute label identifying its origin. The particle swarm constitutes a Lagrange-style flow field description, dynamically displaying the gathering, dispersal, flow, and stagnation of crowds within the venue.

[0043] In one implementation, supplementary virtual particles are generated in the grid space model based on the macroscopic density distribution function within the target region and using the Monte Carlo sampling method. Assigning initial velocities to these virtual particles that conform to the characteristics of the local flow field includes the following steps: The target region is mapped to the mesh space model, and the total number of target particles in each tetrahedral mesh cell within the target region is calculated by integration. Generate a candidate set of random coordinate points that follow a uniform distribution within the tetrahedral geometry of each tetrahedral mesh cell; The structural conductivity property of the tetrahedral mesh element is used to validate the random coordinate point candidate set and invalid coordinate points that fall inside physical obstacles are removed. Based on the total number of target particles, a corresponding number of coordinate points are extracted from the verified candidate set of random coordinate points as the generation positions of virtual particles. Search for the velocity vectors of existing physical particles within the neighborhood of the tetrahedral mesh cell where the virtual particle is located, calculate the velocity vectors of the physical particles using inverse distance weighted interpolation, and assign the velocity vectors to the virtual particle as its initial velocity.

[0044] In this embodiment, all tetrahedral mesh cells in the mesh space model are traversed, and a geometric inclusion test algorithm is used to determine which mesh cells' centroids or the entire set of cells are located within the target region, thereby constructing a subset of mesh cells corresponding to the target region. Then, based on the macroscopic density field provided by the population density thermal data, the volume integral of each tetrahedral mesh cell in the subset is calculated to determine the theoretically required number of particles within that cell.

[0045] After determining the number of particles required to be generated for each tetrahedral mesh cell, candidate particle positions need to be generated within the cell. This involves random point generation techniques under tetrahedral geometric constraints. For each tetrahedral mesh cell, there are four vertex coordinates... Define a candidate set of random coordinate points that follows a uniform distribution, using the principle of the barycentric coordinate system. Introduce three random variables that are uniformly distributed in the interval [0,1]. First, the first two variables are transformed to ensure that the points are uniform on the triangle surface. Then, the third variable is used to extend the points to the interior of the tetrahedron. The specific generation formula is: If Then let ;like Then, the weight coefficients are adjusted according to specific rules, and the final weight vector is calculated. This makes the candidate point coordinates This process is repeated to generate far more candidate points than the total number of target particles, forming a high-density set of random coordinate candidate points. This centroid-based generation method ensures that the probability density of candidate points within the tetrahedron is completely uniform, avoiding artificial clustering at vertices or edges.

[0046] Although the generated set of candidate random coordinate points is geometrically located inside the tetrahedron, it does not consider the physical obstacle constraints within the venue. Therefore, the structural connectivity properties of the mesh cells must be used for rigorous validity verification. The structural connectivity tensor of the current tetrahedral mesh cell is read, which encodes information about impassable normal vectors. For each coordinate point in the candidate set, its position is checked to see if it lies within the obstruction region defined by the structural connectivity tensor or within a buffer zone extremely close to the physical boundary. In practice, this can be combined with the pre-calculated impassable normal vectors and obstacle distance fields of the mesh cells. If the distance of a candidate point from the obstacle surface is less than a preset human radius threshold, or if the point lies in invalid space in the direction where the connectivity coefficient is the obstruction value, then the point is determined to be an invalid coordinate point. This eliminates invalid coordinate points that are physically located within wall cavities, inside columns, or in fenced restricted areas.

[0047] After validity verification, the remaining candidate points constitute the set of valid coordinate points. Next, the final selection and instantiation are performed based on the previously calculated total number of target particles. For the i-th tetrahedral mesh element, the total number of target particles is... The number of valid candidate points after verification is A random sampling algorithm without replacement is used to sample from... Randomly selected from 10 valid candidate points These points will serve as the final virtual particle generation locations. If... This indicates that the effective space of the grid cell is too small to accommodate the expected number of particles, in which case all of them should be used. The remaining ungenerated particle quota is transferred to neighboring grid cells for generation. These selected coordinate points are formally instantiated as virtual particle objects, and their spatial distribution follows both macroscopic density statistics and strictly adheres to microscopic physical boundary constraints. For each newly generated virtual particle, a search radius is set centered on its position to search for existing physical particles within the neighborhood of the grid space model. If physical particles exist within the neighborhood, the set of velocity vectors of these physical particles is extracted. The velocity vector of the virtual particle is calculated using the inverse distance weighted interpolation method. The calculation formula is as follows: The weight Distance from the virtual particle to the j-th physical particle Inversely proportional, take (p is a power exponent). If there are no physical particles in the neighborhood, the value is assigned by referring to the historical average velocity vector of the grid cell or the preset guiding velocity along the structural conduction direction.

[0048] In one implementation, the following steps are taken to calculate the velocity gradient and strain rate of the dynamic trajectory particle swarm within the tetrahedral mesh cells, and to construct a nonlinear mapping relationship between strain rate and flow resistance based on a shear-thickening fluid model to generate a real-time crowd rheological viscosity field: The velocity vector of the dynamic trajectory particle swarm within the tetrahedral mesh element is mapped to the vertex nodes of the tetrahedral mesh element using shape function interpolation to obtain the nodal velocity field. A velocity gradient tensor describing local deformation characteristics is constructed by calculating the spatial partial derivatives of the nodal velocity field, and the symmetric part of the velocity gradient tensor is extracted as the strain rate tensor of the dynamic trajectory particle swarm. The shear rate amplitude of crowd flow within a tetrahedral mesh cell is quantified by calculating the second invariant of the strain rate tensor. The apparent viscosity of the tetrahedral mesh element is calculated by substituting the shear rate amplitude into a preset shear-thickening fluid model that includes the rheological index and consistency coefficient. After traversing all tetrahedral mesh elements and calculating the corresponding apparent viscosity values, and performing spatial smoothing, the rheological viscosity field of the target venue is generated.

[0049] In this embodiment, for each tetrahedral mesh element in the mesh space model, all dynamically tracing particles falling into it are identified, and the instantaneous velocity vectors of these particles are extracted. This can be achieved using shape function interpolation techniques from the finite element method. Specifically, the velocity contribution of the internal particles is distributed to the four vertex nodes of the mesh element using the linear shape function of the tetrahedral element. If a node is shared by multiple tetrahedral elements, the contributions from all adjacent elements are accumulated and averaged. This effectively smooths the random motion noise of individual particles, constructing a continuous and smooth nodal velocity field defined on the mesh vertices.

[0050] Within a tetrahedral mesh element, the velocity field is assumed to be linear, therefore its spatial partial derivatives are constant. The velocity gradient tensor L at the center of the mesh element is directly calculated using the nodal velocities and the spatial derivatives of the shape functions. The velocity gradient tensor is a... A matrix whose elements Describes velocity components Along spatial direction The rate of change. To separate the rotational and deformational motion of the fluid, the symmetric part of the velocity gradient tensor is extracted, and the strain rate tensor D is constructed. The calculation formula is as follows: ,in This is the transpose matrix. The strain rate tensor physically describes the expansion and contraction rates and shear deformation rates of a micro-fluid mass within a crowd per unit time, eliminating the pure rigid body rotation component. This quantitative description captures the compression (negative diagonal elements) and displacement (off-diagonal elements) characteristics of the crowd within local regions.

[0051] As a symmetric tensor, the strain rate tensor's second invariant remains unchanged under coordinate transformations. Under the assumption of incompressible fluids, the concept of generalized shear rate or effective shear rate is adopted, and its calculation is based on the deformation form of the second invariant. The shear rate amplitude within a tetrahedral mesh element is calculated. The formula is ,in The trace operation represents the sum of the diagonal elements of a matrix. The intensity of crowd flow within the grid cells was quantified; a higher value indicates a greater difference in relative speed between crowds, suggesting strong friction, pushing, or rapid queue changes. Next, based on the physical properties of shear-thickening fluids, a nonlinear mapping relationship between strain rate and flow resistance needs to be established to simulate the frictional lock-up effect in high-density crowds: when crowds attempt to move quickly or experience violent conflicts, the physical resistance between people increases sharply, manifested as an instantaneous increase in fluid viscosity, thus restricting the degrees of freedom of flow.

[0052] Finally, all tetrahedral mesh cells in the grid space model are traversed, and their corresponding apparent viscosity values ​​are calculated one by one. Global spatial smoothing is then performed to generate the final crowd rheological viscosity field. Specifically, Laplace smoothing or Gaussian filtering algorithms are used to iteratively smooth the viscosity data across the entire field. For each mesh cell, its viscosity value is updated to a weighted average of its own and its neighboring mesh viscosity values. The processed viscosity field is a continuously distributed scalar field in three-dimensional space, clearly showing the distribution of flow resistance in different areas of the venue. In the visualization, high-viscosity areas correspond to extremely crowded, difficult-to-move solidified zones, while low-viscosity areas correspond to freely flowing, unobstructed zones.

[0053] In one embodiment, substituting the shear rate amplitude into a preset shear-thickening fluid model that includes the rheological index and consistency coefficient, and calculating the apparent viscosity value of the tetrahedral mesh element includes the following steps: Collect environmental decibel data within the target venue and calculate the acceleration variance of the dynamic trajectory particle swarm; A panic index model based on environmental decibel data and acceleration variance was constructed, and a positive correlation mapping function between the panic index and the rheological index output by the panic index model was established. A modified shear-thickening fluid equation is constructed using the updated rheological index and the preset consistency coefficient; The shear rate amplitude is input into the shear thickening fluid equation to calculate the apparent viscosity of the tetrahedral mesh element.

[0054] In this embodiment, ambient decibel (dB) data is collected in real time by acoustic sensors deployed throughout the venue, and short-time energy analysis or root mean square (RMS) values ​​are calculated to quantify the intensity level of ambient noise, as unusual screaming or shouting often foreshadows sudden events. Simultaneously, for the particle set within each grid cell, the statistical dispersion of its acceleration vector, i.e., the acceleration variance, is calculated. First, the instantaneous acceleration is calculated using the time difference of the particle's velocity. Then, the deviation of the acceleration of all particles within the unit from the average acceleration is statistically analyzed. Acceleration variance reflects the disorder and abrupt changes in crowd movement. When a crowd is in a state of panic or disturbance, individuals often frequently stop abruptly, turn sharply, or collide, leading to a significant increase in acceleration variance. Then, based on the extracted environmental decibel data and acceleration variance, a panic index model is constructed, defining the panic index. Let be a normalized scalar, with values ​​ranging from [0,1]. A weighted fusion method is used to construct the model, expressed by the formula: ,in and These are the weighting coefficients. and These are normalized mapping functions for decibels and variance, respectively, used to map the original physical quantities to dimensionless probability intervals. Furthermore, a positive correlation mapping function is established between the panic index and the rheological index n in the rheological model. Under normal conditions, the value of n for crowd movement is close to 1 or slightly greater than 1; however, under panic conditions, the shear thickening effect is sharply enhanced. The mapping function is defined as follows. ,in Based on the basic rheological index, This is the sensitivity coefficient.

[0055] Using the updated rheological index and the preset consistency coefficient, a modified shear-thickening fluid equation is constructed. In the constitutive equation of the shear-thickening fluid equation, the consistency coefficient K is kept as a constant or slowly varying variable characterizing the basic environmental drag, while the rheological index n is replaced with a dynamic value. The corrected constitutive equation is updated to $. This equation implies that under the same physical compression, if the crowd is in a state of high panic, the flow resistance calculated by the model will be much greater than in a calm state. Finally, the previously calculated shear rate amplitude is input into this modified shear-thickening fluid equation to complete the final solution for the apparent viscosity value. Specifically, for each tetrahedral mesh element, its current shear rate and panic index are read simultaneously, and then a power operation is performed. The apparent viscosity value at that time and location is obtained. A high viscosity value directly indicates that the area is at an extremely dangerous rheological solidification critical point, meaning that the crowd has lost its ability to move due to physical compression and psychological panic. A dynamically changing viscosity distribution map will be generated through grid-by-grid calculations.

[0056] In one implementation, a fluid vortex evolution equation is constructed by combining the rheological viscosity field of the crowd and the particle velocity field based on dynamic trajectory particle swarm. The identification of dynamically anomalous regions in the grid space model by solving the fluid vortex evolution equation includes the following steps: A particle velocity field is constructed based on a dynamic trajectory particle swarm. The curl operator is applied to the particle velocity field to calculate the vorticity vector of each tetrahedral mesh element. The fluid vorticity evolution equation is constructed based on the time change rate of the vorticity vector. Decompose the fluid vorticity evolution equation and calculate the vortex stretching term caused by the velocity gradient and the viscous diffusion term caused by the rheological viscosity field of the crowd, respectively. When the magnitude of the vortex stretching term is greater than that of the viscous diffusion term, the flow field is determined to be in an unstable growth state of vorticity accumulation. By integrating the square of the modulus of the vorticity vector within the tetrahedral grid cell by volume, the pseudo-vortex energy characterizing the rotational intensity of the local flow field is calculated. The tetrahedral mesh element set whose pseudo-vortex energy exceeds the preset growth threshold and is in an unstable growth state is marked as a dynamic anomaly region.

[0057] In this embodiment, the irregularly distributed particle velocities are projected onto the mesh nodes using the shape function interpolation technique of tetrahedral mesh elements, forming a continuously defined particle velocity field throughout the entire three-dimensional space. Based on this, the curl operator from vector analysis is introduced. Spatial differentiation is performed on the velocity field to calculate the vorticity vector within each tetrahedral mesh element. The vorticity vector physically describes precisely the instantaneous rotational intensity and direction of a fluid particle around its own axis, and its calculation formula is as follows: Furthermore, to explore the dynamic evolution trend of crowd flow, a partial differential equation describing the evolution of vorticity over time, namely the fluid vorticity evolution equation, was constructed based on the vorticity transport principle in fluid mechanics. This equation closely links the time-varying rate of vorticity with the spatial structure of the flow field, no longer limiting itself to static density or velocity observations, but examining the stability of crowd flow from the perspective of dynamic evolution.

[0058] Next, the fluid vorticity evolution equation is decomposed into independent physical mechanism terms driving the flow field changes. The focus is on calculating the vortex stretching term caused by the velocity gradient and the viscous diffusion term caused by the crowd rheological viscosity field. (Vortex stretching term) The process of fluid particles being elongated and twisted under the influence of a velocity gradient is described. According to the law of conservation of angular momentum, when a vortex tube is elongated, its cross-section decreases, and its rotational angular velocity inevitably increases dramatically. This is the core nonlinear mechanism leading to local concentration of kinetic energy and enhanced turbulence in the flow field. Simultaneously, the viscous diffusion term is calculated. ,in This is the apparent viscosity value calculated based on the shear thickening model. This value represents the dissipation effect of viscous friction on vortices within the fluid. Physically, it manifests as the uniform diffusion of concentrated vorticity to the surrounding area, thus providing a damping effect to smooth disturbances and restore flow field order.

[0059] Based on the comparison of the magnitudes of the physical terms, a criterion for determining the dynamic stability of the flow field is established. Within each computational step, tetrahedral mesh elements are scanned one by one, and the magnitudes of the vortex stretching term and the viscous diffusion term are compared. When the calculation results show that the magnitude of the vortex stretching term is significantly greater than that of the viscous diffusion term, i.e. Physically, this means that the rate at which new vortices are generated due to velocity shear in the flow field within this region exceeds the rate of viscous dissipative vortices. In this scenario, minute disturbances within the flow field cannot be absorbed by its own damping mechanism; instead, they are amplified through positive feedback, leading to an exponential increase in local rotational motion and disordered counter-movement. Based on this, the local flow field is determined to be in an unstable growth state of vorticity accumulation. To quantitatively assess the destructive energy that the flow field in this unstable state may release, the vorticity vector within each tetrahedral mesh element is... The square of the modulus is used for full-volume integration. Pseudo-vortex energy, as a scalarized energy measure, is calculated using the following formula: ,in dV represents the geometric volume of a tetrahedral mesh element. This physical quantity intuitively reflects the total rotational kinetic energy contained within the mesh element; the higher the value, the more intense the rotation, collision, and turbulence of the population within that region.

[0060] Finally, the entire venue's grid space model is traversed. First, it checks whether the pseudo-vortex energy of the grid cells exceeds a preset empirical growth threshold, which represents the upper limit of rotational energy that can be tolerated for normal crowd flow. Second, it confirms whether the grid cell was determined by previous steps to be in an unstable growth state. Only tetrahedral grid cell sets that simultaneously meet the conditions of excessive energy and unstable trend are marked as dynamically abnormal regions. These regions are highlighted in the 3D model, indicating high-risk locations in the venue where crowds are experiencing violent collisions, potentially forming crowd vortices, or where trampling is imminent.

[0061] In one embodiment, applying a curl operator to the particle velocity field to calculate the vorticity vector of each tetrahedral mesh element, and constructing a fluid vorticity evolution equation based on the time-varying rate of change of the vorticity vector, includes the following steps: Collect environmental decibel data within the target venue and calculate the acceleration variance of the dynamic trajectory particle swarm; The spatial distribution of a dynamic trajectory particle swarm is transformed into a continuous macroscopic density scalar field using the kernel density estimation method, and the spatial density gradient vector of the macroscopic density scalar field is calculated. A pressure scalar field is constructed by combining acceleration variance and ambient decibel data, and the pressure gradient vector of the pressure scalar field is calculated. By using the cross product of the spatial density gradient vector and the pressure gradient vector, the baroclinic torque term is constructed, and a fluid vorticity evolution equation containing the baroclinic torque term is constructed.

[0062] In this implementation, acoustic sensors deployed throughout the venue collect ambient decibel (dB) data in real time, and perform short-time energy analysis or calculate the root mean square (RMS) value to quantify the intensity level of ambient noise, as unusual screaming or shouting often foreshadows sudden events. Simultaneously, for each set of particles within a grid cell, the statistical dispersion of its acceleration vector, i.e., the acceleration variance, is calculated. Next, kernel density estimation is used to infer the overall probability density function from a limited data sample. Specifically, a kernel function with a specific bandwidth is placed centered on the position of each dynamically tracing particle, and all kernel functions are superimposed in space to generate a continuously distributed macroscopic density scalar field in three-dimensional space. Each point in this field corresponds to a specific density value, smoothly representing the degree of crowd aggregation within the venue and eliminating numerical noise caused by the random positions of individual particles. Based on this, the gradient operator in vector calculus is used to perform spatial partial derivative operations on the macroscopic density scalar field to calculate the spatial density gradient vector. This vector points in the direction of the fastest density increase, and its magnitude reflects the spatial drasticness of density change, revealing the transition characteristics between the core and peripheral areas of crowd aggregation.

[0063] Combining the acceleration variance and environmental decibel data obtained in previous steps, a pressure scalar field is constructed to characterize the psychological tension and physical compression intensity of a crowd, and its gradient is then calculated. Here, pressure is a generalized concept that integrates physical congestion and psychological panic. The functional relationship of the pressure scalar field is typically constructed using a weighted summation, with the normalized values ​​of the environmental decibel data and the normalized values ​​of the acceleration variance used as input variables to quantify the potential energy of the crowd's desire to escape or be "compressed" within the area. Subsequently, the gradient operator is applied to the pressure scalar field to calculate the pressure gradient vector. This vector physically represents the direction of the driving force of pressure, pointing towards the direction of the fastest pressure drop (e.g., evacuation exits or open areas), and also reveals the distribution location of the high-pressure core area. Finally, using the calculated spatial density gradient vector and pressure gradient vector, a crucial baroclinic torque term in fluid mechanics is constructed through a cross product operation. In an ideal barotropic fluid, density isosurfaces and pressure isosurfaces are parallel. However, in complex crowd flow fields, the non-uniformity of panic propagation (pressure gradient) and the non-uniformity of physical aggregation (density gradient) often do not coincide, resulting in an angle between them. Therefore, a vector cross product operation is needed to generate a new vector field, namely the baroclinic torque term. The physical meaning of this term represents the spontaneously generated rotational torque due to the mismatch between the density and pressure distribution within the fluid. When the density gradient direction is perpendicular to the pressure gradient direction, this torque reaches its maximum value, forcibly driving fluid particles to rotate, thereby creating vortices in a flow field that might otherwise be merely translational.

[0064] Finally, by synthesizing all the above fluid dynamic components, a fluid vorticity evolution equation including a baroclinic torque term is constructed. This equation is a partial differential equation, and its specific form is as follows:

[0065] in, The vorticity vector. For particle velocity field, This is the apparent viscosity value. For macroscopic density scalar fields, For pressure scalar field, The first term on the right-hand side represents the time rate of change of the vorticity vector. The term "vortex stretching" describes the enhanced vorticity stretching caused by the velocity gradient; the second term on the right is... The term represents viscous diffusion, describing the eddy current dissipation caused by apparent viscosity; the third term on the right-hand side... This is the newly introduced baroclinic torque term. By solving this equation, the distribution and intensity of the vorticity field at future moments can be predicted.

[0066] In one implementation, the virtual shear stress of the dynamic trajectory particle swarm on the channel boundary of the target venue is calculated by combining the rheological viscosity field of the crowd and the structural conductivity properties of the tetrahedral mesh elements, and the path failure risk areas in the mesh space model are identified based on the virtual shear stress, including the following steps: In a mesh space model, a mesh whose structural conductivity properties of tetrahedral mesh elements change abruptly is defined as a channel boundary layer. The virtual shear stress of the passenger flow on the tangential plane of the channel boundary layer is calculated by combining the apparent viscosity and strain rate tensor at the channel boundary layer. The erosion resistance coefficient is set according to the physical width and geometry of the channel boundary layer; Based on dynamic trajectory particle swarm optimization and combined with virtual shear stress and erosion resistance coefficient, the attenuation rate of the channel's conduction function corresponding to the channel boundary layer is calculated. The grid area where the channel boundary layer is located, where the rate of attenuation of the dredging function exceeds the preset failure threshold, is identified as a path failure risk area.

[0067] In this embodiment, all tetrahedral mesh cells are traversed, and the structural conductivity tensor stored within each cell is examined. This tensor has previously been used to quantify the passageability of the space. In open areas far from obstacles, the conductivity properties exhibit uniform isotropic characteristics. However, near walls, columns, or exit edges, the conductivity properties undergo a drastic jump from full conductivity to complete blockage within a very short spatial distance. The rate of change of conductivity coefficients between adjacent mesh cells is calculated using a spatial gradient algorithm, and the set of mesh cells with rates of change exceeding a preset gradient is defined as the channel boundary layer. Within the defined channel boundary layer, the real-time apparent viscosity values ​​within the channel boundary layer mesh are retrieved, and the strain rate tensor describing the fluid deformation intensity is extracted. Based on the constitutive equations of non-Newtonian fluids, the shear stress vector acting on the boundary tangent plane is calculated. Since the channel boundary layer is adjacent to impassable physical surfaces, the boundary normal vector is used. A projection operator is constructed to project the stress tensor onto a direction parallel to the wall surface, thereby obtaining the virtual shear stress. The calculation formula can be expressed as follows: Virtual shear stress simulates the tangential resistance generated when a high-density crowd passes through a narrow passage due to friction and compression between the shoulders and torso of individuals and walls or railings. Higher values ​​indicate more intense interaction between the crowd and the facility; this friction not only slows down the movement of people at the edges but also reduces the effective width of the entire passage through force chains.

[0068] To assess the inherent properties of different passage structures in resisting congestion and maintaining traffic flow, it is necessary to set erosion resistance coefficients based on their physical characteristics. This depends on the physical width of the passage boundary layer, the complexity of its geometry, and the properties of its surface material. Analyzing the geometric descriptions of passages in CIM data reveals that wider main roads have greater tolerance and stronger resistance to interference; however, narrow passages such as funnel-shaped, right-angle turns, or stepped passages are prone to deadlock. The erosion resistance coefficients are then constructed. The calculation model is as follows: ,in Based on material factors, The effective physical width of the channel. The geometry factor (higher values ​​for straight channels, lower values ​​for irregular channels) and the erosion resistance coefficient characterize the upper limit of the channel's theoretically maximum boundary shear stress without functional failure. Next, the particle flux flowing through the channel boundary layer is statistically analyzed. A functional decay model is constructed by combining virtual shear stress and erosion resistance coefficient. This model indicates that the channel's throughput capacity is not constant but gradually decreases with the accumulation of boundary friction. The rate of conduction function decay is calculated, and its functional form is: ,in This is the time-scale adjustment constant. The formula indicates that when the pedestrian flow is very high and boundary friction is extremely strong, while the passage itself is narrow or irregular, the congestion reduction function will decrease extremely rapidly.

[0069] Finally, based on the calculated attenuation index, the final path failure risk area detection results are generated. Specifically, a failure threshold is set, which represents the bottom line for maintaining the channel's function. The system iterates through the grid areas containing all channel boundary layers, comparing the calculated attenuation rate of the guiding function with the failure threshold in real time. When the attenuation rate of the guiding function in a certain area exceeds the failure threshold, it is determined that although the physical structure of the channel is intact, its guiding function has substantially failed due to severe boundary friction and arching effects. The system then marks this area as a path failure risk area. These areas are typically presented as dark red warning bands on the visualization interface, covering specific exits or corridor sections.

[0070] The present invention also discloses a multi-source venue traffic anomaly detection system based on spatiotemporal feature fusion, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion as described above.

[0071] The processor can be a central processing unit (CPU). Of course, depending on the actual use, it can also be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), off-the-shelf programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor can be a microprocessor or any conventional processor, etc., and this application does not limit it.

[0072] The memory can be an internal storage unit of a computer device, such as a hard disk or RAM, or an external storage device, such as a plug-in hard disk, smart memory card (SMC), secure digital card (SD), or flash memory card (FC) provided on the computer device. Furthermore, the memory can be a combination of internal storage units and external storage devices of a computer device. The memory is used to store computer programs and other programs and data required by the computer device. The memory can also be used to temporarily store data that has been output or will be output. This application does not limit this.

[0073] The present invention also discloses a computer-readable storage medium storing instructions that, when executed by a processor, configure the processor to perform the multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion described in any of the above embodiments.

[0074] The computer program can be stored in a machine-readable medium. The computer program includes computer program code, which can be in the form of source code, object code, executable file, or certain middleware. The machine-readable medium includes any entity or device capable of carrying computer program code, recording media, USB flash drive, portable hard drive, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signals, telecommunication signals, and software distribution media, etc. It should be noted that the machine-readable medium includes, but is not limited to, the above-mentioned components.

[0075] The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion described in the above embodiments is stored in the computer-readable storage medium and loaded and executed on the processor to facilitate the storage and application of the above method.

[0076] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of protection of this application is limited to these examples; within the framework of this application, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of one or more embodiments of this application as described above, which are not provided in detail for the sake of brevity.

[0077] One or more embodiments in this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of this application. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of one or more embodiments in this application should be included within the protection scope of this application.

Claims

1. A method for detecting anomalies in venue traffic from multiple sources based on spatiotemporal feature fusion, characterized in that, Includes the following steps: Obtain the CIM data and indoor model of the target venue, combine the CIM data and indoor model of the venue and use the tetrahedral subdivision algorithm to construct a grid space model of the target venue containing multiple tetrahedral grid units; By accessing multi-source monitoring data within the target venue and aligning spatiotemporal coordinates, the individual audience members identified in the multi-source monitoring data are mapped onto a grid space model, generating a dynamic trajectory particle swarm containing location and velocity information. The velocity gradient and strain rate of a dynamic trajectory particle swarm within a tetrahedral mesh cell are calculated, and a nonlinear mapping relationship between strain rate and flow resistance is constructed based on a shear-thickening fluid model to generate a real-time crowd rheological viscosity field. By combining the rheological viscosity field of the crowd and the particle velocity field constructed based on the dynamic trajectory particle swarm, a fluid vortex evolution equation is constructed. By solving the fluid vortex evolution equation, dynamic anomaly regions in the grid space model are identified. The virtual shear stress of the dynamic trajectory particle swarm on the channel boundary of the target venue is calculated by combining the rheological viscosity field of the crowd and the structural conductivity properties of the tetrahedral mesh unit, and the path failure risk area in the mesh space model is identified based on the virtual shear stress. The flow anomaly detection results for the target venue are generated by integrating the dynamic anomaly region and the path failure risk region.

2. The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to claim 1, characterized in that, The process of combining CIM data and the venue's indoor model and using a tetrahedral subdivision algorithm to construct a grid spatial model of the target venue containing multiple tetrahedral grid cells includes the following steps: Analyze the geometric topology information in CIM data and venue interior models to identify the boundary surface equations of walls, columns, fences and stairs in the target venue; The effective access range of the target venue in three-dimensional space is defined based on the boundary surface equation, and the effective access range is divided into a set of unstructured tetrahedral mesh elements using the tetrahedral subdivision algorithm. Based on the centroid coordinates and volume parameters of each tetrahedral mesh element in the tetrahedral mesh element set, the spatial positional relationship between the tetrahedral mesh element and the boundary surface equation is calculated, and the impassable normal vectors existing in the tetrahedral mesh element are determined according to the spatial positional relationship. Based on the impassable normal vector, a structural conduction tensor in the local coordinate system is constructed as the structural conduction property of the tetrahedral mesh element. The conduction coefficient along the direction of the impassable normal vector in the structural conduction tensor is set to a preset stagnation value, and the conduction coefficient perpendicular to the plane of the impassable normal vector in the structural conduction tensor is set to a unit passability value. Laplace smoothing is applied to the structural conduction tensor of all tetrahedral mesh elements to form a continuously distributed mesh space model of the target venue.

3. The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to claim 1, characterized in that, The process of accessing multi-source monitoring data within the target venue and mapping individual audience members identified from the multi-source monitoring data to a grid space model through spatiotemporal coordinate alignment to generate a dynamic trajectory particle swarm containing location and velocity information includes the following steps: Acquire video surveillance data containing audience pixel coordinates collected by surveillance cameras in the target venue, and crowd density heat map data containing the number of audience members in the area collected by people counting cameras in the target venue; Based on a preset time synchronization protocol, video surveillance data and crowd density thermal data are aligned to the same sampling time. The viewer's pixel coordinates are back-projected onto the grid space model using a perspective transformation matrix to generate a three-dimensional observation point representing the position of an individual viewer in the real world; The three-dimensional observation points are input into a preset Kalman filter. The state transition matrix is ​​used to predict the state of the audience individual at the next moment. The optimal estimated position and optimal estimated velocity vector of the audience individual are updated in combination with the current three-dimensional observation points. The entity particles representing the audience individual are generated in the grid space model. By combining video surveillance data and crowd density thermal data, target areas in the target venue with only crowd density thermal data are identified. Based on the macroscopic density distribution function within the target area, supplementary virtual particles are generated in the grid space model using the Monte Carlo sampling method. The virtual particles are given an initial velocity that conforms to the characteristics of the local flow field. By merging physical and virtual particles, a dynamic trajectory particle swarm containing audience flow information within the target venue is generated.

4. The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to claim 3, characterized in that, The process of generating supplementary virtual particles in the grid space model based on the macroscopic density distribution function within the target region and using the Monte Carlo sampling method, and assigning initial velocities to the virtual particles that conform to the characteristics of the local flow field, includes the following steps: The target region is mapped to the mesh space model, and the total number of target particles in each tetrahedral mesh cell within the target region is calculated by integration. Generate a candidate set of random coordinate points that follow a uniform distribution within the tetrahedral geometry of each tetrahedral mesh cell; The structural conductivity property of the tetrahedral mesh element is used to validate the random coordinate point candidate set and invalid coordinate points that fall inside physical obstacles are removed. Based on the total number of target particles, a corresponding number of coordinate points are extracted from the verified candidate set of random coordinate points as the generation positions of virtual particles. Search for the velocity vectors of existing physical particles within the neighborhood of the tetrahedral mesh cell where the virtual particle is located, calculate the velocity vectors of the physical particles using inverse distance weighted interpolation, and assign the velocity vectors to the virtual particle as its initial velocity.

5. The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to claim 1, characterized in that, The process of calculating the velocity gradient and strain rate of the dynamic trajectory particle swarm within the tetrahedral mesh element, and constructing a nonlinear mapping relationship between strain rate and flow resistance based on a shear-thickening fluid model to generate a real-time crowd rheological viscosity field includes the following steps: The velocity vector of the dynamic trajectory particle swarm within the tetrahedral mesh element is mapped to the vertex nodes of the tetrahedral mesh element using shape function interpolation to obtain the nodal velocity field. A velocity gradient tensor describing local deformation characteristics is constructed by calculating the spatial partial derivatives of the nodal velocity field, and the symmetric part of the velocity gradient tensor is extracted as the strain rate tensor of the dynamic trajectory particle swarm. The shear rate amplitude of crowd flow within a tetrahedral mesh cell is quantified by calculating the second invariant of the strain rate tensor. The apparent viscosity of the tetrahedral mesh element is calculated by substituting the shear rate amplitude into a preset shear-thickening fluid model that includes the rheological index and consistency coefficient. After traversing all tetrahedral mesh elements and calculating the corresponding apparent viscosity values, and performing spatial smoothing, the rheological viscosity field of the target venue is generated.

6. The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to claim 5, characterized in that, The step of substituting the shear rate amplitude into a preset shear-thickening fluid model that includes the rheological index and consistency coefficient to calculate the apparent viscosity value of the tetrahedral mesh element includes the following steps: Collect environmental decibel data within the target venue and calculate the acceleration variance of the dynamic trajectory particle swarm; A panic index model based on environmental decibel data and acceleration variance was constructed, and a positive correlation mapping function between the panic index and the rheological index output by the panic index model was established. A modified shear-thickening fluid equation is constructed using the updated rheological index and the preset consistency coefficient; The shear rate amplitude is input into the shear thickening fluid equation to calculate the apparent viscosity of the tetrahedral mesh element.

7. The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to claim 1, characterized in that, The process of constructing a fluid vortex evolution equation by combining the rheological viscosity field of the population and the particle velocity field based on dynamic trajectory particle swarm, and identifying dynamic anomaly regions in the grid space model by solving the fluid vortex evolution equation includes the following steps: A particle velocity field is constructed based on a dynamic trajectory particle swarm. The curl operator is applied to the particle velocity field to calculate the vorticity vector of each tetrahedral mesh element. The fluid vorticity evolution equation is constructed based on the time change rate of the vorticity vector. Decompose the fluid vorticity evolution equation and calculate the vortex stretching term caused by the velocity gradient and the viscous diffusion term caused by the rheological viscosity field of the crowd, respectively. When the magnitude of the vortex stretching term is greater than that of the viscous diffusion term, the flow field is determined to be in an unstable growth state of vorticity accumulation. By integrating the square of the modulus of the vorticity vector within the tetrahedral grid cell by volume, the pseudo-vortex energy characterizing the rotational intensity of the local flow field is calculated. The tetrahedral mesh element set whose pseudo-vortex energy exceeds the preset growth threshold and is in an unstable growth state is marked as a dynamic anomaly region.

8. The multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to claim 5, characterized in that, The process of calculating the virtual shear stress of the dynamic trajectory particle swarm on the channel boundary of the target venue by combining the rheological viscosity field of the crowd and the structural conductivity properties of the tetrahedral mesh elements, and identifying the path failure risk area in the mesh space model based on the virtual shear stress, includes the following steps: In a mesh space model, a mesh whose structural conductivity properties of tetrahedral mesh elements change abruptly is defined as a channel boundary layer. The virtual shear stress of the passenger flow on the tangential plane of the channel boundary layer is calculated by combining the apparent viscosity and strain rate tensor at the channel boundary layer. The erosion resistance coefficient is set according to the physical width and geometry of the channel boundary layer; Based on dynamic trajectory particle swarm optimization and combined with virtual shear stress and erosion resistance coefficient, the attenuation rate of the channel's conduction function corresponding to the channel boundary layer is calculated. The grid area where the channel boundary layer is located, where the rate of attenuation of the dredging function exceeds the preset failure threshold, is identified as a path failure risk area.

9. A multi-source venue traffic anomaly detection system based on spatiotemporal feature fusion, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion as described in any one of claims 1 to 8.

10. A computer-readable storage medium storing instructions thereon, characterized in that, When executed by a processor, the instruction causes the processor to be configured to perform the multi-source venue traffic anomaly detection method based on spatiotemporal feature fusion according to any one of claims 1 to 8.

Images