Traffic traceability and flow statistics method based on hidden Markov

Abstract

The invention discloses a traffic traceability and flow statistics method based on hidden Markov, and belongs to the field of intelligent traffic. The method comprises the following steps: step 1, cleaning and standardizing floating car data; step 2, establishing a path matching model based on hidden Markov; step 3, carrying out parallel optimization on road network partitioning and path matching through multiple threads; and step 4, path traceability and flow analysis. The problems of low precision, poor efficiency and insufficient traceability in vehicle track traceability in the prior art are solved, spatial-temporal feature modeling, multi-thread parallel optimization and dynamic path backtracking technologies are fused, the traffic traceability and flow statistics method based on the hidden Markov model is provided, traffic traceability and flow statistics can be accurately achieved, and the traffic traceability and flow statistics efficiency are improved. And powerful support is provided for optimization of the intelligent traffic system.

Description

Technical Field

[0001] The present invention relates to the field of intelligent transportation, and in particular to a floating vehicle trajectory path matching method based on a hidden Markov model, which is used for traffic flow analysis and congestion cause analysis. Background Art

[0002] Vehicle trajectory tracing technology is key to analyzing traffic conditions in the intelligent transportation sector. It can restore actual vehicle paths based on floating vehicle latitude and longitude data, providing support for traffic flow analysis, congestion cause investigation, and intelligent transportation system optimization. With the development of intelligent transportation, floating vehicle trajectory data has become a crucial source for addressing traffic monitoring and management challenges. Its high timeliness and wide coverage facilitate tasks such as traffic flow monitoring. However, accurately tracing vehicle paths and analyzing the causes of traffic problems from massive amounts of data remains a challenge for intelligent transportation. Existing vehicle trajectory tracing methods often rely on traditional path matching algorithms, which suffer from limited accuracy and efficiency in complex traffic environments and struggle to accurately restore actual vehicle paths. Traditional geometric projection and topological matching methods fail to consider the temporal characteristics of vehicle travel, are susceptible to GPS noise, and are time-consuming to process large amounts of data, making them inadequate for city-level real-time traffic analysis. Currently, research on vehicle trajectory tracing technology focuses primarily on improving accuracy or efficiency. Key approaches include describing the tracing process based on probabilistic models and establishing tracing models based on vehicle travel characteristics. Strategies such as selecting the shortest path are often employed in these applications. However, these studies fail to consider the uncertainty and real-time nature of attribution in complex traffic scenarios, and are unable to accurately describe the attribution process. Existing vehicle trajectory attribution research has not proposed a reasonable attribution method for complex scenarios, and cannot achieve accurate simulation and dynamic tracking of vehicle trajectory attribution in complex traffic scenarios. Summary of the Invention

[0003] To address the above issues, the present invention proposes a traffic tracing and flow statistics method and system based on a hidden Markov model (HMM). By integrating spatiotemporal feature modeling, multi-threaded parallel optimization, and dynamic path backtracking technology, the present invention addresses the problems of low accuracy, poor efficiency, and insufficient tracing capabilities in existing technologies. The specific technical solution is as follows:

[0004] Step 1: Data preprocessing

[0005] This article uses floating vehicle data mainly from floating vehicles equipped with positioning and communication devices.

[0006] Table 1 Requirement fields

[0007] Field Name meaning DeviceId Equipment number of the floating vehicle Longitude Longitude recorded by the floating vehicle equipment Latitude Latitude recorded by the floating vehicle equipment GPStime Time standard used by the Global Positioning System

[0008] "DeviceID" represents the device ID of the floating vehicle and is used to uniquely identify the device.

[0009] "Longitude" refers to the longitude recorded by the floating vehicle's equipment, indicating the east-west coordinate of the vehicle's location. "Latitude" refers to the latitude recorded by the floating vehicle's equipment, indicating the north-south coordinate of the vehicle's location. "GPStime" is the time standard used by the Global Positioning System, providing an accurate time reference for the location information recorded by the device, starting at 00:00:00 on January 6, 1980.

[0010] 1. Data cleaning and standardization

[0011] In raw GPS datasets, vehicle numbers are often presented in a mixed format of letters and numbers. To enhance data security and storage efficiency, a hash function is used to process these numbers. This process converts the complex mixed format into a unique digital identifier.

[0012] The GPS timestamps of raw data are often presented in a concise format. This format is limited when conducting in-depth time series analysis and makes it difficult to effectively integrate and compare with other time series data. It is necessary to uniformly convert these GPS timestamps to the more universal "YYYY-MM-DDHH:MM:SS" format.

[0013] At the same time, when pre-processing GPS data, the data records need to be sorted according to vehicle ID and time sequence. In this process, focus on and eliminate duplicate data reported by the same vehicle at almost the same time point. For the same vehicle, if the latitude and longitude in the reported data are exactly the same and the time difference is less than 1 second, it is considered a duplicate record and will be eliminated. The specific operation is as follows: for two adjacent data with the same number in the data set R1 = (ID1, t1, lat1, lon1) and R2 = (ID2, t2, lat2, lon2)

[0014] ID1=ID2,

[0015] lat1=lat2

[0016] lon1=lon2

[0017] |t1-t2|<1,

[0018] Then remove the data point, where ID1 and ID2 represent the unique identifiers of the vehicle; t1 and t2 are the corresponding timestamps; lat1 and lon1 and lat2 and lon2 represent the latitude and longitude coordinates of the two records respectively.

[0019] The sliding window algorithm is used to filter the vehicle travel path and the vehicle stop status. The vehicle trajectory dataset is S = {s1, s2, ..., s n} Each trajectory point s i =(xi ,y i ,t i ) contains the longitude and latitude coordinates (x i ,y i )

[0020] and time t i For any time t k Construct a sliding time window:

[0021] W k ={s i ||t i -t k ∣≤Δt},Δt=5min

[0022] The maximum spatial displacement within the window is:

[0023]

[0024] D k <θ,θ=50m

[0025] Among them, S a and S b That is, the time window W k Any two trajectory points within. (x a ,y a ) and (x b ,y b ) is the trajectory point S a and S b The longitude and latitude coordinates of the vehicle are shown in Figure 1. The displacement threshold θ is set to 50 m according to the Urban Traffic Operation Evaluation Standard. When the window meets the above conditions, that is, the vehicle displacement within the window is less than 50 meters, it is marked as a stop point.

[0026] At the same time, the drift points in the GPS data are accurately eliminated. Specifically, the drift points are identified by calculating the instantaneous speed between adjacent points. Assume that a vehicle is at a continuous time point t i and t i+1 Corresponding to the spatial position point S i (x i ,y i ) and S i+1 (x i+1 ,y i+1 )

[0027] Calculate the difference in latitude and longitude:

[0028] Δφ=x2-x1

[0029] Δλ=y2-y1

[0030] Haversine formula calculates the distance between two data points

[0031]

[0032] Here, R is the radius of the Earth, which is 6371 km based on actual conditions. d is the average speed between two points, and x1 and x2 are the latitudes of the two points:

[0033]

[0034] According to the Road Traffic Safety Law of the People's Republic of China, we identify data with v < 120km / h, which meets the standard of urban road speed limit, as GPS drift points and remove them; the time interval is

[0035] Δt=t i+1 -t i ;

[0036] 2 Coordinate system conversion

[0037] When processing raw GPS data, given that it uses the WGS-84 coordinate system (EPSG:4326), directly projecting the longitude and latitude in this coordinate system onto a flat map will inevitably cause distortion errors due to the inherent properties of the Earth's curved surface and flat surface transformations. To effectively address this issue, the present invention converts the data to the Beijing local coordinate system (EPSG:2416) and utilizes the Gauss-Krüger projection method. This projection method significantly reduces the degree of distortion during the projection process.

[0038] Because basic 4D maps are extremely slow to process, we considered deeply analyzing and converting their SHP layers. This layer was parsed and constructed into a directed graph structure, which was stored and presented in GraphML format. In this directed graph model, nodes are assigned specific geographic meanings, representing road intersections or locations of key geographic significance. Edges represent road segments, and each edge is assigned detailed attribute information, including key parameters such as segment length and number of lanes.

[0039] Step 2. HMM-based path matching model

[0040] Path matching based on a hidden Markov model (HMM) first identifies a set of candidate road segments based on the location of the trajectory point and the surrounding road network information. Each candidate segment is abstracted as a hidden state in a Markov chain and assigned an observation state probability. This probability is used to quantify the likelihood of a match between the GPS point and the candidate segment. If the GPS point is very close to a road segment, it is assigned a higher output observation probability. Subsequently, state transition probabilities are calculated for edges between adjacent hidden states in the Markov chain, taking into account factors such as road connectivity, vehicle direction, and speed. Finally, a method for solving the HMM, such as the Viterbi algorithm, is used to obtain the maximum likelihood path.

[0041] 1. Generation of candidate road segments

[0042] By dynamically adjusting the search radius, with each GPS point as the center, the search radius r of the candidate road segment is generated and adaptively adjusted according to the vehicle's real-time speed v and the sampling interval Δt. The formula is:

[0043] r=3σ+v·Δt

[0044] Where σ is the GPS positioning standard deviation (5 meters), and v·Δt is the theoretical maximum displacement of the vehicle within the sampling interval.

[0045] 2. Candidate road section screening:

[0046] To further improve the accuracy and efficiency of path matching within the area defined by the calculated search radius r, candidate road segments need to be screened. Priority is given to road segments with a heading angle deviation of less than 30 degrees from the vehicle. The vehicle heading angle is accurately calculated using consecutive GPS points and reflects the vehicle's direction of travel. This screening criterion is based on the assumption that vehicle direction is consistent in real traffic scenarios, meaning that vehicles do not make drastic turns in a short period of time.

[0047]

[0048] Δθ=min(|α-θ|,2π-|α-θ|)

[0049]

[0050] 3 Observation probability calculation

[0051] For the data points to be matched, the matching degree between a certain measurement point and multiple candidate matching points is judged based on the observation probability. Specifically, given a certain road section k i , the measurement point s observed at time t t , and its corresponding observation probability p(z t ∣k i ) can be calculated using the following formula:

[0052]

[0053] Where: p(z t ∣k i ) indicates that the vehicle is actually located at section k i In the case of t The probability of measurement,

[0054] z t It is the two-dimensional latitude and longitude position measurement representative at time t;

[0055] ki Represents the discrete state in HMM, i.e., the i-th road segment;

[0056] s t,i For a given measurement value z t and section k i , which is the road segment k i Upper distance measurement point z t The nearest point.

[0057] |z t -s t,i | great_circle : For a given measurement value z t and section k i Upper closest point x t,i The great circle distance on the Earth's surface between the two points, which is used to measure the proximity of the measurement point to the road segment.

[0058] σ z : The standard deviation of GPS measurements, used to quantify the noise in GPS measurements. This value is estimated based on analysis of actual data.

[0059] 4 State transition probability

[0060] In order to obtain the probability of the vehicle reaching the next state from the current state at time point t, consider constructing the state transition probability. The specific formula is as follows:

[0061]

[0062] Among them, e is a natural constant, approximately 2.71828, which is the base of the exponential function. t is the distance between the current state and the next state:

[0063] d t =|z t -z t+1 | great_circle -|s t,j -s t+1,j | route

[0064] Among them, |z t -z t+1 | great_circle Represents the measurement point z at time t t and the measurement point z at time t+1 t+1 The great circle distance on the surface of the Earth.

[0065] |s t,j -s t+1,j | route represents the road segment k at time t i Upper closest point s t,i and the road section s at time t+1i Upper closest point s t+1 ,j * The driving distance between * and j * Represents the actual road section in the actual driving route. It is the parameter of the exponential distribution, which is used to adjust the transition probability with the distance difference d t The rate of change.

[0066] In addition, β is the parameter of the exponential distribution, which determines the sensitivity of the state transition probability to the distance difference in the map matching algorithm.

[0067]

[0068] Among them, |z t -z t+1 | great_circle ,|s t,j -s t+1,j | route Consistent with the above formula, median() represents the median operation. is a constant coefficient that plays a specific scaling role.

[0069] 5. Path backtracking

[0070] In the framework of Hidden Markov Model (HMM), Viterbi algorithm solves the maximum a posteriori probability path through dynamic programming. Let the trajectory point sequence be S = {s1, s2, ..., s T}, the candidate road segment set is k={k1,k2,...,k j}, assign initial probabilities to all candidate road segments at the first observation point s1 (t=1):

[0071]

[0072] in:

[0073] π(k j ) is the initial state probability, and the prior information is uniformly distributed (π(k j )=1 / O), O is the total number of road sections in the candidate road section set k, indicating that when there is no observation information, the vehicle is on the candidate road section k j The probability of P(s1|k j ) is the above observation probability, which refers to the vehicle actually being in the candidate road section k j The probability of observing the first observation point s1 under the condition .

[0074] For each time step t (t = 2, ..., T) and candidate road segment k i , calculate the maximum cumulative probability and the optimal predecessor node:

[0075]

[0076] δ t (k j ) refers to the candidate road segment k at time step t j It is calculated by considering the cumulative probability of each candidate road segment in the previous time step, the state transition probability and the current observation probability.

[0077] δ t-1 (k i ):At time step t-1, candidate road segment k i The maximum cumulative probability.

[0078] P(k j ∣k i ):State transition probability, which refers to the vehicle’s transition from road section k i Transfer to section k j It describes the transfer relationship between road segments.

[0079] P(s t ∣k j ): At time step t, the vehicle is on candidate road segment k j Observation point s t The probability of observation.

[0080]

[0081] Where: P(k j ∣k i ) is the state transition probability, which refers to the vehicle leaving the road section k at time t. i Transfer to section k j The probability of ψ t (k j ) record reaches k j The optimal front-wheel drive section.

[0082] At the last time step T, the candidate segment with the largest cumulative probability is selected as the path end point

[0083]

[0084] Among them, δ T (k j ) refers to the candidate road segment k at time step T j The maximum cumulative probability.

[0085] Recursively obtain the optimal path from t=T-1,T-2,...,1:

[0086]

[0087] in, The road segment on the optimal path at time step t.

[0088] At time step t+1, the road segment is reached The optimal predecessor section is obtained by continuously backtracking this value and recursively deducing from the end point of the path to obtain the entire optimal path.

[0089] The optimal path sequence is Step 3. Multi-threaded parallel optimization

[0090] 1. Map data segmentation

[0091] The entire road network for the matching map is a graph G = (V, E), where V is the node set and E is the edge set. The road network G is divided into n subgraphs G according to administrative divisions or traffic grids. i =(V i ,E i ),i=1,2,,n satisfies and Each subgraph G i Assign independent threads to process to achieve parallel computing.

[0092] 2. Parallel Operation

[0093] Matching each floating car trajectory requires complex computations such as calculating state transition probabilities and observing sequence probabilities, which takes a long time to process. Consider using Python's Concurrent.Futures library to dynamically manage thread pools to avoid resource contention and deadlock. Specifically, the thread pool size M can be dynamically adjusted based on the number of system cores C and the number of tasks N:

[0094] M=min(C,N)

[0095] Where M is the size of the thread pool, C is the number of CPU cores, and N is the number of sub-regions. By adjusting the thread pool size to balance the computing load,

[0096] In order to ensure load balancing in parallel operations and ensure that the workload of each thread is balanced, we use the standard deviation σ to measure the load balance. Assume there are T threads, and the workload of each thread is m1,m2,...,m T , the average load is:

[0097]

[0098] Load balancing meets the following requirements:

[0099]

[0100] By minimizing the standard deviation σ, we can optimize the distribution of tasks so that the workload of each thread is as close to m as possible. avg , thus achieving load balancing.

[0101] Ensure optimal utilization of system resources Multi-threading technology reduces the matching time of 8408 tracks from 4 hours to 1877 seconds, increasing efficiency by 3.4 times

[0102] Step 4. Tracing and Traffic Analysis

[0103] The trajectory sequence of each vehicle obtained in step 3

[0104] 1. Path tracing

[0105] When a user selects a specific route, the trajectory of all vehicles along that route is retrieved based on the route's geometry ID. By comparing the vehicle's position with the selected route's geometry, the database is searched for all routes within that geometry and the set of geometry IDs for the vehicle trajectories that meet the criteria is returned. The path geometry range is [x1, x2] on the x-axis and [y1, y2] on the y-axis.

[0106] G j ={K i ∣x j ∈[x1,x2]∧y j ∈[y1,y2]}

[0107] On the selected path, the set K of all vehicle trajectories that meet the following conditions i Included in the geometry ID set G j middle:

[0108] B j ={k i ∣G j ∩t i ∈Δt}

[0109] The system is based on the geometry ID G j The time limit Δt is used to filter out the trajectories of all corresponding vehicles related to the path within the restricted time period, understand the source of the vehicle trajectory, and realize trajectory tracing.

[0110] 2. Traffic aggregation

[0111] The number of trajectories passing through each path within the map range is B j , generate dynamic heat maps, identify high traffic source areas (such as attached Figure 2 ).

[0112] 3. Visualization and Interaction

[0113] Through the PyQt front-end interface (see attached Figure 3 ), users can flexibly select specific road sections, and the system will display in real time the vehicle trajectories and traffic source distribution matching the selected route. This graphical display allows traffic management departments to intuitively identify traffic flow distribution on the road. Furthermore, through different color-coding (such as yellow for traffic source), managers can quickly locate potential sources of congestion, optimize traffic scheduling and management strategies, and improve the efficiency and accuracy of emergency response. This real-time, dynamic path tracing function contributes to the development of intelligent transportation systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0115] Figure 1 A schematic diagram of a path matching process for constructing a matching model provided by an embodiment of the present invention.

[0116] Figure 2 This is a diagram of traffic statistics results based on the matching model provided by an embodiment of the present invention.

[0117] Figure 3 This is a diagram of the path tracing results based on the matching model provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0118] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0119] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions of the present invention are described in detail based on the floating vehicle data within 1 hour on October 10, 2023 in Tongzhou District, Beijing (a total of 8408 trajectories).

[0120] Step 1: Data Preprocessing

[0121] Input data: Floating vehicle data from 00:00 to 01:00 on October 10, 2023 in Tongzhou District, Beijing, including the following fields:

[0122] Table 1. Floating car data format

[0123]

[0124] Duplicate records were removed: 213 records were removed for the same vehicle (DeviceID) with a time difference of less than 1 second and identical latitude and longitude. A sliding window (10-minute window duration, 50-meter displacement threshold) was used to mark 124 stop points. The instantaneous speeds of adjacent points were calculated, and outliers with speeds > 120 km / h (37 in total) were removed. The WGS-84 coordinate system was converted to the Beijing local coordinate system (EPSG:2416). The Tongzhou District SHP map data was parsed, and a directed graph structure was constructed (nodes were intersections, edges were road sections), which was stored in GraphML format.

[0125] Step 2: HMM-based path matching

[0126] Dynamic formula for search radius: r = 3σ + v·Δt, where σ = 5m, Δt = 10s. Set heading angle deviation threshold: 30°. Observation probability calculation: Based on Gaussian distribution,

[0127]

[0128] Step 3: Multi-threaded parallel optimization

[0129] The Tongzhou District road network is divided into 8 sub-regions (based on administrative divisions), and each sub-region is processed by an independent thread. The thread pool size is dynamically adjusted and set to 16 threads. The standard deviation of the number of trajectories processed by each thread is <

[0130] 5% to ensure that tasks are evenly distributed.

[0131] Step 4: Source tracing and traffic analysis

[0132] In the matching trajectory, a single vehicle trajectory generates an independent tag ID, and the time information corresponding to the longitude and latitude points is intercepted and associated with the trajectory, generating the matching trajectory structure as follows:

[0133] Table 2. Path matching results

[0134]

[0135] The vehicle trajectories passing through Xinhua Avenue between 0:30 and 0:45 were screened, and a total of 123 were matched.

[0136] Typical trajectory examples:

[0137] Table 3. Path tracing results

[0138]

[0139] Traffic metering:

[0140] The Tongzhou District road network traffic situation and heat map are generated through path aggregation, and the expression is as follows:

[0141] Table 4. Flow metering results

[0142]

[0143] The real-time heat map is displayed through the PyQt interface, and users can select an area to query the vehicle trajectory distribution and traffic source area of a certain road section.

Claims

1. A traffic tracing and flow statistics method based on Hidden Markov Model. "DeviceID" represents the device ID of a floating vehicle, which is used to uniquely identify the device; "Longitude" refers to the longitude recorded by the floating vehicle device, indicating the east-west coordinate of the vehicle's location; "Latitude" is the latitude recorded by the floating vehicle device, indicating the north-south coordinate of the vehicle's location; "Gpstime" is the time standard used by the Global Positioning System, which provides an accurate time reference for the location information recorded by the device. It is characterized by: The following steps are involved: 1) Data cleaning and standardization In the original GPS dataset, vehicle numbers are often presented in a mixed format of letters and numbers. A hash function is used to process them and convert the complex mixed format numbers into unique digital identifiers. GPS timestamps are uniformly converted to the more common "YYYY-MM-DDHH:MM:SS" format; At the same time, when preprocessing GPS data, the data records need to be sorted according to vehicle ID and time sequence. In this process, focus on and eliminate duplicate data reported by the same vehicle at almost the same time point. For the same vehicle, if the latitude and longitude in the reported data are exactly the same and the time difference is less than 1 second, it is considered a duplicate record and is eliminated. The specific operation is as follows: for two adjacent data with the same number in the data set R1 = (ID1, t1, lat1, lon1) and R2 = (ID2, t2, lat2, lon2) ID1=ID2, lat1=lat2 lon1=lon2 |t1-t2|<1, Then the data point is eliminated, where ID1 and ID2 represent the unique identifier of the vehicle; t1, t2 are the corresponding timestamps; lat1, lon1 and lat2, lon2 represent the longitude and latitude coordinates of the two records respectively; The sliding window algorithm is used to filter the vehicle travel path and the vehicle stop status. The vehicle trajectory dataset is S = {s1, s2, ..., s n } Each trajectory point s i =(x i ,y i ,t i ) contains the longitude and latitude coordinates (x i ,y i ) and time t i For any time t k Construct a sliding time window: W k ={s i ||t i -t k ∣≤Δt},Δt=5min The maximum spatial displacement within the window is: D k <θ,θ=50m Among them, S a and S b That is, the time window W k Any two trajectory points within (x a ,y a ) and (x b ,y b ) is the trajectory point S a and S b The longitude and latitude coordinates of the vehicle are as follows: the displacement threshold θ is set to 50m. When the vehicle displacement in the window is less than 50 meters, it is marked as a stop point. At the same time, the drift points in the GPS data are eliminated; specifically, the drift points are identified by calculating the instantaneous speed between adjacent points; a vehicle is set at a continuous time point t i and t i+1 Corresponding to the spatial position point S i (x i ,y i ) and S i+1 (x i+1 ,y i+1 ) Calculate the latitude and longitude differences: Δφ=x2-x1 Δλ=y2-y1 Haversine formula calculates the distance between two data points Where R is the radius of the Earth, d is the average velocity between two points, and x1 and x2 are the latitudes of the two points: According to the data with v < 120km / h, which meets the standard of urban road speed limit, it is determined as GPS drift point and removed; among them, the time interval is Δt=t i+1 -t i ; 2) Coordinate system conversion When processing the original GPS data, since the data uses the WGS-84 coordinate system, the data is converted to the Beijing local coordinate system (EPSG:2416) and the Gauss-Krüger projection method is used; Because basic 4D maps are extremely slow to process, we considered deeply analyzing and converting their SHP layers. This layer was parsed and constructed into a directed graph structure, which was stored and presented in GraphML format. In this directed graph model, nodes were assigned specific geographic meanings, representing road intersections or locations of key geographic significance. Edges represented road segments, and each edge was assigned detailed attribute information, including segment length and number of lanes. Step 2. HMM-based path matching model 2.1 Generation of candidate road segments By dynamically adjusting the search radius, with each GPS point as the center, the search radius r of the candidate road segment is generated and adaptively adjusted according to the vehicle's real-time speed v and the sampling interval Δt. The formula is: r=3σ+v·Δt Where σ is the GPS positioning standard deviation, which is 5 meters, and v·Δt is the theoretical maximum displacement of the vehicle within the sampling interval; 2.2 Candidate road section screening: Candidate road sections are screened within the area defined by the calculated search radius r; sections with a heading angle deviation of less than 30 degrees from the vehicle are prioritized; the vehicle heading angle is accurately calculated using consecutive GPS points, reflecting the vehicle's direction of travel; Δθ=min(|α-θ|,2π-|α-θ|) (Convert to radians) 3) Observation probability calculation For the data points to be matched, the matching degree between a certain measurement point and multiple candidate matching points is judged based on the observation probability; specifically, given a certain road section k i , the measurement point s observed at time t t , and its corresponding observation probability p(z t ∣k i ) is calculated using the following formula: Where: p (z t ∣k i ) indicates that the vehicle is actually located at section k i In the case of t The probability of measurement, z t It is the two-dimensional latitude and longitude position measurement representative at time t; k i Represents the discrete state in HMM, i.e., the i-th road segment; s t,i For a given measurement value z t and section k i , which is the road segment k i Upper distance measurement point z t nearest point; |z t -s t,i | great_circle : For a given measurement value z t and section k i Upper closest point x t,i The great circle distance on the surface of the earth between the two points is used to measure the proximity of the measuring point to the road segment; σ z : standard deviation of GPS measurement, used to quantify the size of GPS measurement noise; This value is estimated based on analysis of actual data; 4) State transition probability In order to obtain the probability of the vehicle reaching the next state from the current state at time point t, consider constructing the state transition probability. The specific formula is as follows: Among them, e is a natural constant and the base of the exponential function; d t is the distance between the current state and the next state: d t =||z t -z t+1 | great_circle -|s t,j -s t+1,j | route | Among them, |z t -z t+1 | great_circle Represents the measurement point z at time t t and the measurement point z at time t+1 t+1 The great circle distance on the surface of the Earth; |s t,j -s t+1,j | route represents the road segment k at time t i Upper closest point s t,i and the road section s at time t+1 i Upper closest point s t+1 ,j * The driving distance between * and j * Represents the actual road section in the actual driving route; is the parameter of the exponential distribution, used to adjust the transition probability with the distance difference d t rate of change; In addition, β is the parameter of the exponential distribution, which determines the sensitivity of the state transition probability to the distance difference in the map matching algorithm; Among them, |z t -z t+1 | great_circle ,|s t,j -s t+1,j | route Consistent with the above formula, median() represents the median operation. is a constant coefficient; 5) Path backtracking In the framework of Hidden Markov Model (HMM), the Viterbi algorithm solves the maximum a posteriori probability path through dynamic programming; let the trajectory point sequence be S = {s1, s2, ..., s T }, the candidate road segment set is k={k1,k2,...,k j }, Assign initial probabilities to all candidate road segments at the first observation point s1 (t=1): Where: π(k j ) is the initial state probability, and the prior information is uniformly distributed (π(k j )=1 / O), O is the total number of road sections in the candidate road section set k, which represents the probability that the vehicle is on the candidate road section kj when there is no observation information; P(s1|k j ) is the above observation probability, which refers to the vehicle actually being in the candidate road section k j , the probability of observing the first observation point s1; For each time step t (t = 2, ..., T) and candidate road segment k i , calculate the maximum cumulative probability and the optimal predecessor node: δ t (k j ) refers to the maximum cumulative probability of the candidate link kj at time step t; it is calculated by considering the cumulative probability of each candidate link in the previous time step, the state transition probability, and the current observation probability; δ t-1 (k i ):At time step t-1, candidate road segment k i The maximum cumulative probability of P(k j ∣k i ):State transition probability, which refers to the vehicle’s transition from road section k i Transfer to section k j The probability of ; it describes the transfer relationship between road segments; P(s t ∣k j ): At time step t, the vehicle is on the candidate road segment kj and observes the observation point s t The probability of , that is, the observation probability; Where: P(k j ∣k i ) is the state transition probability, which refers to the vehicle leaving the road section k at time t. i The probability of transferring to the road segment kj; ψ t (k j ) record reaches k j The optimal front-drive section; At the last time step T, the candidate segment with the largest cumulative probability is selected as the path end point Among them, δ T (k j ) refers to the candidate road segment k at time step T j The maximum cumulative probability of Recursively obtain the optimal path from t=T-1,T-2,...,1: in, The road segment on the optimal path at time step t; At time step t+1, the road segment is reached The optimal predecessor section is obtained by continuously backtracking this value and recursively deducing from the end point of the path to obtain the entire optimal path. The optimal path sequence is Step 3. Multi-threaded parallel optimization 3.1 Map Data Blocking The entire road network of the matching map is a graph G = (V, E), where V is the node set and E is the edge set; the road network G is divided into n subgraphs G according to administrative divisions or traffic grids. i =(V i ,E i ), i=1,2,…,n satisfies and Each subgraph G i Assign independent threads to process to achieve parallel computing; 3.

2. Parallel Operation Matching each floating car trajectory requires complex operations such as calculating state transition probabilities and observation sequence probabilities, which takes a long time to process. Consider using Python's Concurrent.Futures library to dynamically manage thread pools to avoid resource contention and deadlock. Specifically, the thread pool size M can be dynamically adjusted based on the number of system cores C and the number of tasks N: M=min(C,N) in, M is the size of the thread pool, C is the number of CPU cores, and N is the number of sub-regions; by adjusting the thread pool size to balance the computing load, In order to ensure load balance in parallel operations and ensure that the workload of each thread is balanced, we use the standard deviation σ to measure the balance of the load; suppose there are T threads, and the workload of each thread is m1,m2,...,m T , the average load is: Load balancing meets the following requirements: By minimizing the standard deviation σ, we optimize the distribution of tasks and achieve load balancing; Step 4. Tracing and Traffic Analysis The trajectory sequence of each vehicle obtained in step 3 4.

1. Path tracing When a user selects a specific route, the trajectory of all vehicles on the route is obtained based on the geometric ID of the route. By comparing the vehicle position with the geometric range of the selected route, the database is searched for all routes within the geometric range and the geometric ID set of the vehicle trajectories that meet the conditions is returned. The geometric range of the route is [x1, x2] and [y1, y2] on the x-axis and y-axis respectively. G j ={K i ∣x j ∈[x1,x2]∧y j ∈[y1,y2]} On the selected path, the set K of all vehicle trajectories that meet the following conditions i Included in the geometry ID set G j middle: B j ={k i ∣G j ∩t i ∈Δt} The system is based on the geometric IDG j The time limit Δt is used to filter out all corresponding vehicle trajectories related to the path within the restricted time period, understand the source of the vehicle trajectory, and realize trajectory tracing; 4.

2. Traffic aggregation The number of trajectories passing through each path within the map range is B j ,generate dynamic heat maps to identify high traffic source areas; 4.

3. Visualization and interaction.

Images