GB2582531A - Method for detecting traffic anomally of urban road with equidistant spatial-temporal division - Google Patents

Abstract

Equidistant spatial-temporal division is used in the detection of a traffic anomaly. A spatial-temporal subzone is created by dividing one day into several time segments of 30 minutes, each defined as a time subzone; dividing one implemented area into several spatial segments of 200m x 200m, each defined as a spatial subzone; and defining an intersection of one time subzone and one spatial subzone as a spatial-temporal subzone. Historical and real-time GNSS positioning data of floating vehicles is respectively pre-processed into sampled historical and real-time vehicle speed data relating to historical and real-time trajectories. These are used in a finite mixed modelling method to establishing respective historical and real-time travel speed probability distributions. Jensen-Shannon divergence is then used to measure a difference between said historical and real-time travel speed probability distributions and an anomaly index of traffic conditions is determined from said difference. The accuracy of said anomaly index is then evaluated.

Description

Method for Detecting Traffic Anomaly of Urban Road with Equidistant Spatial-temporal Division This is a divisional application of the following 'parent' application: Application Number: GB1909405.1 Application Source: NPI Filing Date: 30 December 2017 PCT Application No.: PC1782017/058531 Priority Claimed: 30 December 2016 in UFO -PC7/182016/058105

Technical Field

The invention belongs to the technical field of traffic detection. In particular, the present invention relates to a real-time detection method for urban road traffic incidents. Through GNSS positioning devices of the floating cars, the spatial-temporal position information can be obtained. After data preprocessing, map matching and data fusion, the probability distribution of travel speed for specific spatial-temporal region will be derived. According to the changes of velocity distribution, the city road traffic incidents can be effectively identified.

Technological background

Traffic incident detection or traffic anomaly detection is an important part of urban traffic management, and it is also one of the core functions of intelligent transportation system. Traffic accidents mainly include traffic accidents, vehicle breakdown, goods falling, road traffic facilities damage or failure, and other special events that cause traffic flow disorder. This kind of events can easily lead to traffic congestion and reduce the traffic capacity of the road, which seriously affects the normal operation of the whole road traffic system. Traffic incident detection enables traffic managers to understand traffic information in a timely manner and take appropriate guidance and control measures to reduce adverse effects of traffic incidents.

Traffic incident detection can be divided into manual and automatic methods. Manual methods include patrol cars, emergency telephone reporting and video surveillance, which cannot meet the needs of traffic management because of the consumption of human resources and poor real-time ability Automatic methods depend on Automated Incidence Detection (AID) algorithm. The basic principle is to identify traffic abnormal events by detecting traffic flow changes at different locations. The commonly used AID algorithm including pattern recognition algorithms (e.g., California algorithm, Monica algorithm)., statistical prediction algorithm (e.g., exponential smoothing method, Kalman filtering algorithm), algorithm based on traffic flow theory (e.g., McMaster algorithm) and intelligent recognition algorithms (e.g., artificial neural network, fuzzy logic algorithm).

However, the current detection methods have the disadvantages of high requirements for facilities and high computational complexity. Furthermore, these methods cannot make a further judgment for incidents. The invention uses the track data of the GNSS positioning device transmitted by taxi and bus to establish historical traffic status database and real-time traffic state database and identifies traffic incident by analyzing the difference of traffic flow characteristics reflected by them. This method has the characteristics of real-time and parallel processing, with high recognition rate and low computation requirements for detection facilities, and thus, is suitable for detecting incidents in urban road traffic data under real-time floating car positioning data environment.

At present, in view of the traffic incident detection, there arc the following representative teclmiques: A U.S. patent application, US 20160148512, disclosed the principle and implementation of a traffic incident detection and reporting system. The system consists of sensors, the communication module, the mobile processing module and the user interaction module. The sensors are used for data collection around the vehicle; the communication module is used for sending and receiving data of the vehicle; the processing module is used for processing and analysis of the relevant vehicle within a certain range of data and generating traffic incident reports; the user interaction module can provide the user traffic incident report. The scheme is a traffic incident detection technology based on vehicle to vehicle and vehicle to infrastructure communication network. It can detect various kinds of traffic events by using the information collected by sensors. However, due to the installation and experiment of sensors and communication units, it is difficult to implement, and the processing capacity of mobile processing units is limited. Moreover, the mobile and fixed message receiver ends are required, so the system itself has relative higher failure probability and poor reliability.

A China patent application, CN 104809878 A, disclosed a method of using bus GPS data to detect the abnormal state of urban road traffic. The scheme obtains the link delay time index based on the GPS historical records, then obtains the instantaneous speed, the periodic average speed, the weighted moving average speed and the multi vehicle average speed through the GPS real-time data. The canonical variable analysis algorithm is used to detect the traffic incidents. This scheme does not require additional detection facilities and is convenient to implement. However, the representation of traffic situation is too simplistic to analyze the characteristics and causes of traffic incidents. Furtherly, the division of traffic scene is lack of evidence and fails to consider the influence of weather and other factors on the change of traffic status.

Summary of the invention

To clarify the content of the invention more clearly, the technical terms involved are explained as follows.

Floating vehicle: also known as detection vehicle. The vehicle including the city bus and taxi on which the positioning device is installed.

GNSS: Global Navigation Satellite System, e.g., GPS. GLONASS, GALILEO and Compass.

Spatial-temporal subzone: the area divided by the two dimensions of time and space which reflects the traffic status in a certain period and a certain range of space. Divide the day into several time fragments, such as 0:00-0:10, 0:10-0:20 and so on. Each time segment is called a time subzone. Divide the implemented area into several spatial segments, e.g., the area of I21.58°E-I 21.59°E and 31.16°N-31.17°N. Each spatial fragment is called a spatial subzone. The intersection of an arbitrary time subzone and spatial subzone is called a spatial-temporal subzone, such as the spatial-temporal segment in the area of 121.58°E-121.59°E and 31.16°N-31.17°N together with the time range of 0:00-0:10.

Historical trajectory data: historical trajectory data is the trajectory data accumulated a long time and stored in the database. Historical trajectory data are dynamically changing data, which need updating in a timely manner then reprocessing and analyzing regularly, so as to ensure the accuracy of the extraction of historical traffic characteristics. The data within each spatial-temporal subzone can be processed in parallel to improve efficiency. In the present invention, historical trajectory data can be referred to as historical data for short.

Real-time trajectory data: real-time trajectory data is a set of trajectory data within the latest time segment. In the present invention, real-time trajectory can be referred to as real-time data for short.

Traffic situation: the general term of the comprehensive situation of traffic operation in a certain time and in a certain space.

Traffic anomalies: traffic flow disorder caused by traffic accidents, vehicle breakdown, goods falling, damage or failure of road traffic facilities.

Traffic anomaly severity: the severity of traffic flow disorder. In other words, the difference between the traffic flow characteristics in normal status and in anomaly status.

Traffic anomaly index: measure of traffic abnormality severity. The range is 0-10. The greater value means the more serious of anomaly severity.

Traffic environment: all of the external effects and forces acting on the road traffic participants. It includes road conditions, traffic facilities, landforms, weather conditions, and other traffic activities.

Map matching: the process of association of geographic coordinates with urban road networks. Peak hour flow: the maximum amount of traffic flow in one day of a city road section.

Finite mixed model: a mathematical method of simulating complex density function with simple density function. The finite mixed model of the set of variables y and the number of components K can be expressed as: P(Y)= (Y) Response variables: variables that change according to the independent variable, also known as the dependent variable.

Bayesian information criterion: a criterion for model selection among a finite set of models, the model with the lowest BIC is preferred. The computing formula is: BIC = -21n L + k * In n wherein L is the maximum of the likelihood function, K is the number of unknown parameters, and N is the sample size.

Likelihood function: likelihood function is a function of the parameters of a statistical model. When the output x is given, the likelihood function L(01x) is equal to the probability of a given parameter 0 after the variable X. That is, L(0x)=P(X=x0).

Parameter estimation: a method of estimating the unknown parameters in the overall distribution based on the sample extracted from the population.

EM Algorithm: the abbreviation of Expectation Maximization Algorithm. It is an iterative algorithm for the maximum likelihood estimation or the maximum posteriori probability estimation of the probabilistic parameter model with implicit variables.

Kullback-Leibler divergence: a measure of the difference between two probability distributions of P and Jensen-Shannon divergence: a symmetric form of Kullback-Leibler divergence.

K-Medoids algorithm: a clustering algorithm. In each iteration the algorithm selects a point from the current category with the minimum distance to all other points in the current category, as a new central point.

The purpose of the invention is to set up a scheme based on floating car track recording system, which uses historical GNSS positioning data and real-time GNSS positioning data combined with traffic environment information to identify road traffic anomalies. In order to achieve the above purpose, the present invention provides the following technical scheme: The premise of the invention is the floating vehicle (taxi, bus, etc.) equipped with GNSS track recorder, and data center with large-scale storage, computation and real-time task processing ability.

The applicable scope of the invention is as follows: the urban road (including the ground road and the elevated road) with the floating vehicles.

The implementation steps of the present invention include: 1) Determining the spatial-temporal range of detection and setting up spatial-temporal subzones Based on the practical application requirements, the time range and space range of traffic anomaly detection should be determined. The time range can be set for all day long, namely 0:00-24:00; it can also be set for a specific period of time, such as 17:00-20:00. The spatial range can be set up according to administrative divisions to a certain city area, such as Beijing, Shanghai, Huangpu District, etc.; and it can also be set to a certain urban functional area according to the urban spatial structure, for example, a city's central business district and industrial area The establishment of spatial-temporal subzones refers to dividing the detection time span into several smaller time segments, dividing the detection space, namely the implementation area of traffic incident detection, into several smaller space segments. The establishment of spatial-temporal subzones can be implemented by a variety of empirical methods, including equidistant spatial-temporal division and non-equidistant spatial-temporal division.

2) Data pre-processing GNSS positioning data should be processed by data cleaning, data integration, data conversion and data reduction, in order to improve structural level of data. GNSS, or Global Navigation Satellite System, is aerial radio navigation and positioning system that can provide full-time three-dimensional coordinates, velocity and time information at any place on the earth's surface or near-Earth space. It mainly includes GPS (Global Positioning System) of the United States, GLONASS (Global Navigation Satellite System) of Russia, GALILEO of the European Union and Beidou Satellite Navigation System of China. It also includes QZSS of Japan, TRNSS of India and United States WASS, Japan's MSAS and other satellite positioning augmentation systems. In order to establish a unified data distribution standard in different navigation and positioning system equipment, the National Marine Electronics Association (NEMA) formulated the unified communication protocol to regulate GNSS data broadcasting. Therefore, although each member system of GNSS, such as GPS, GLONASS, is established and maintained by different countries and institutions, it does not need to transform the data format thanks to the consistent data distribution format.

Within the selected spatial range, there are many vehicles equipped with GNSS positioning equipment, such as taxis, buses, freight cars, private cars and so on. Based on the current situation of urban traffic data application, urban taxi is usually used as floating vehicle as the data source of traffic incident detection system.

The collected GNSS location information contains some unreasonable information. In order to ensure the accuracy of traffic incident detection and discrimination results, it is necessary to screen the abnormal data to ensure the reliability of the data. These abnormal data include: the data falling outside the detection spatial-temporal range, and the leap of location obviously beyond the reasonable range. The so-called "leap of location obviously beyond the reasonable range" is illustrated below. If the positioning point uploaded by the positioning equipment of a floating car is marked as A at 10:30:00 on a certain day-, the positioning point uploaded by the positioning equipment of the floating car is marked as B at 10:30:30 on the same day, and the distance between position A and position B is 1500 meters, then the driving speed of the floating car is calculated to be at least 180 km/h, which is beyond the common sense, so it is a kind of abnormal spatial position leap, which should be eliminated in data processing.

3) Fast map matching After pre-processing, the GNSS positioning data need to be combined with the urban road network data by projecting the GNSS positioning points to the urban map through map matching algorithm to establish the matching relationship between the positioning points and the road segments, and correct the errors caused by positioning drift.

At present, the electronic map around the world has been detailed, this electronic map can be derived from the city's geographical information system, of course, can also be derived from other ways and means. These electronic maps depict the urban road information in detail, and several road sections can be obtained by dividing them. By means of distance, angle and other information, the positioning points are matched to the road sections, then the positioning information is matched to the actual geographical environment.

4) Representing the path of floating car and the matching of different vehicle paths Given a set of starting points, the vehicle path may not be unique. Complex urban traffic network contains a number of road sections, these different sections are numbered, for example, as LI, L2, and so on. A road section may have two different directions. in this case, the two different directions should be represented as two different road sections with different section number.

The the starting point and end point can usually adopt the intersection point of the urban road network. Given the path of a floating car, it is necessary to select the same path from the path information sent by other floating cars, so as to obtain the same path group between the starting point and the end point.

5) Data sampling The location data of floating cars include location coordinates, instantaneous speed and recording time. In the urban road traffic incident detection method based on floating car data proposed in this patent, data sampling is to select part of the data from all floating car data for subsequent analysis and processing, which depends on the data center's computing power and the pre-proposed accuracy requirements. Data sampling methods can be chosen based on computing power and accuracy requirements. For example, when the data center has strong computing power and requires high detection accuracy, all floating vehicle positioning data can be processed and analyzed comprehensively, while when the data center has limited computing power, it is assumed that the current data center can process 5 for each spatial subzone within I minute, in fact, however, 2000 float car positioning record can be generated in each spatial subzone in one minute. Then 500 float car positioning records can be randomly extracted from 2000 float car positioning records for analysis, so as to obtain the processing results with limited precision within the computing capacity of the data center.

According to the different ways of using floating car data, different attributes of floating car data can be sampled, such as travel speed and travel time. In this patent, the urban road traffic incident detection method based on floating vehicle data is proposed. The urban road traffic incident detection is based on the travel speed. Therefore, data sampling refers to the sampling of vehicle speed.

6) Historical trajectory data analysis and feature extraction The so-called historical trajectory data refer to the floating car trajectory data accumulated in the long-term operation. Using historical floating car trajectory data, a traffic characteristic model of urban road can be established to reflect the general characteristics of urban traffic operation. In this patent, the urban road traffic characteristic model can refer to some specific indicators, such as average speed, weighted average speed, etc. It can also refer to a variety of statistical models, such as the probability distribution of travel speed. In the past, many models used a single indicator to represent the traffic characteristics of a certain road section or region, such as the historical average speed. Although these methods are easy to implement, the accuracies are not high enough and the sensitivities are relatively poor, thus they cannot play a good role in traffic incident detection. Therefore, for each spatial-temporal subzone, this patent uses the probability distribution of traffic characteristic variables to describe the traffic characteristics. establishes the traffic characteristic model and estimates the parameters.

Traffic characteristic variables, including travel speed and travel time, can be collected. The probability distribution of traffic characteristic variables described in this patent refers to the probability distribution of travel speed.

7) Real-time trajectory data analysis and feature extraction The so-called real-time trajectory data refers to the trajectorydata of floating cars in a period of time not far from the current time. Using real-time floating vehicle trajectory data, we can grasp the dynamic changes of traffic characteristics, which can reflect the real-time characteristics of the current traffic operation. This patent describes the current traffic characteristics using the travel speed of the current time and spatial subzone.

8) Incident detection The idea of system state anomaly detection was first put forward by Denning, that is, by monitoring the anomaly of system usage on system audit records, we can detect the events that violate security and may cause system anomalies. Denning's model is independent of any specific system, application environment, system weakness, fault type, and therefore is a general anomaly detection model. The model includes 5 parts: subject, object, audit record, outline, abnormal record and activity rule. Contour is the normal behavior of the subject relative to the object expressed in terms of measurement and statistical models. Denning's model defines three metrics, namely event counter, interval timer and resource measurer. Five statistical models are proposed, namely operation model, mean and standard deviation model, multivariable model, Markov process model and time series model. Denning's model establishes a statistical-based normal behavior profile of the system subject by analyzing the system audit data. When detecting, the audit data in the system is compared with the normal behavior profile of the established subject. If the difference part exceeds a certain threshold, it is considered as an abnormal event. The model lays the foundation of anomaly detection, and many anomaly detection methods and systems developed in the future are based on it.

In recent years, in the development of anomaly detection technology, more artificial intelligence methods have been introduced to improve the performance of anomaly detection. These artificial intelligence methods mainly include data mining, artificial neural network, fuzzy evidence theory and so on. Data mining is used to determine which features arc the most important in a large number of data sets. This technique is mainly used in anomaly detection to find a more concise definition of normal mode, rather than simply enumerate all normal modes as the traditional anomaly detection method. The introduction of data mining makes it possible for the detection system to generalize the normal patterns not included in the training data only by recognizing the main features of the normal patterns. Artificial neural network anomaly detection problem can be regarded as a general data classification problem. in the statistical anomaly detection mentioned above, user behavior data are divided into two categories according to some statistical criteria: anomaly behavior and normal behavior. Because the statistical method has some difficulties in extracting and abstracting audit instances, which may cause large errors, it must rely on some probability distribution assumptions. Generally, it needs to describe the measurement of user behavior based on experience and feeling, so the clustering method of artificial neural network is introduced. Artificial neural network has the ability of self-learning and self-adapting. The neural network is trained by the sample points representing the normal user behavior. Through repeated learning, the neural network can extract the normal user or system activity patterns from the data and encode them into the network structure. The network can determine whether the system is normal. Because of the fuzziness of anomaly evaluation criteria, fuzzy evidence theory is introduced into anomalies. For example, an intrusion detection framework model based on fuzzy expert system can reduce the false alarm rate and the false alarm rate.

This patent proposes an anomaly detection scheme based on statistical features. The basic idea is to measure the difference between historical traffic characteristics and real-time traffic characteristics by Jensen-Shannon divergence, so as to detect abnormal traffic conditions. The scheme has the advantages of good EXPLANABILITY and low computational burden. It not only overcomes the weakness of inaccurate and untimely detection by single statistics, but also avoids the shortcomings of heavy computational burden and high hardware requirements by artificial neural network.

9) Abnormal severity quantitative characterization and abnormal information release The severity of traffic abnormalities should be publicized to the public in a concise and clear manner to avoid possible congestion areas and improve the operational efficiency of urban transport. The severity of the abnormal condition was characterized by traffic abnormality index, which ranged from 0 to 10, with 0 indicating no abnormality and 10 indicating high abnormality.

The location of the anomaly is projected onto the electronic map and is publicly released by the smart mobile device APP.

10) System performance evaluation The evaluation of system performance refers to the evaluation of the accuracy of traffic abnormal state detection, and its evaluation indexes include the false alarm rate and the false alarm rate. The lower the false positive rate and the false negative rate, the better the performance of the system.

In the step 1), the division of the spatial-temporal subzones may specifically adopt the following methods: 11) Equidistant spatial-temporal division method. Determining the segment size of the time dimension, the time segment span is a fixed value, usually taking 30min as a time segment; determining the segment size of the spatial dimension, the spatial segment span is a fixed value, and usually takes a spatial grid of 200m x 200m as a spatial segment; 12) Non-equidistant spatial-temporal division method based on road network density: the segmentation is based on road network density --when the road network density is greater than or equal to 2km/km2, take 30min time segment and 200m x200m space segment; when the road network density is less than 2km/km2, take 30min time segment and 400m x400m space segment; 13) Non-equidistant spatial-temporal division method based on peak hour traffic volume: the segmentation is based on peak hour traffic volume --when the peak hour traffic flow rate is greater than or equal to 1000 vehicles per hour, take 30min time segment and 200mx200m space segment; when peak hour traffic flow rate is less than 1000 vehicles per hour, take 30min time segment and 400mx400m space segment.

The step 3) specifically includes the following steps: 31) Divide the space area to be processed into a grid of a certain size, then the range of each grid area can be expressed as A = {(x,,Y, )1)c, G r)31, E Dit,y,, )1, each grid area contains several road sections, the set of these road sections can be denoted as R. each road section of the abovementioned set of road section Rican be denoted as ij, and assign an index to each road section; 32) Determining the grid area where the positioning point is located, and searching for the section if where the positioning point A is located in the set R" of the road segment by using the distance and the azimuth angle, the matching schemes include: 321) Single point matching scheme: Searching for the road segment closest to point A, if the difference between the direction of travel of the point A and the direction of the section i j is less than the threshold, i.e. 1614 -0,)1< , then the matching process is completed, the abovementioned threshold can be 2.5° , 5° , 10° etc.; if 104 -Girl< go is not satisfied, delete the road segment ij in the search space and continue to search for other road segments until the conditions are met; see Fig. 3 for the matching method; 322) Point sequence matching scheme: This scheme is applicable for high frequency floating vehicle data. Denote the floating vehicle GNSS data acquisition frequency asA=1/in; denote the point adjacent to A in time as I(i"-6), then P(l,4+6) can be defined as the 1-adjacent point of A, and 1)(44-210), 1-)(44+210) can be defined as the 2-adjacent point of A, and so on, then P(t-1-ktp), P(t l+kto) can be defined as the k-adjacent point of A. Whenti,<IHz, k=1 or 2. Take the road segment ij with the smallest distance from the k-adjacent point and A, and calculate 6,r, the mean value of the direction angle of the k-adjacent point of A and A. If 1/9,1. - c5(, is met, the matching process is completed. if not, search in other road segments until the condition PA' 0,3,1 <8, is met; 33) Use the straight line equation of the road segment (if it is a cuned road segment, it is roughly split into straight lines), calculate the projection coordinates of the GNSS positioning point on the road segment, and reduce the error caused by the GNSS positioning drift. The specific method is descript as follows.

Determine the straight line equation of the road segment ij (if the road segment is a curve, it is divided into several straight line segments): 12-Yi= k(x-x,) Where the slope is: k = Y Y ^, The projection line equation is: y Y A = Solve the projected coordinate P as: ky, +k x, +x, k2 +1 + y, + Ia, -a, k2 +1 After the map matching process, combined with the timestamp data of the coordinates of the positioning point, the positioning point is matched to the space-time sub-region.

The step 5) may specifically adopt one of the following methods:

Y

51) A full sample scheme for speed information. The overall travel speed data of each floating vehicle in a time and space sub-zone constitutes the whole. The implementation method is to calculate the travel speed d d --of each vehicle in the spatial-temporal sub-zone 2f,: = " , where di,,...da-ta is the distance of the 1 and 2"6 GNSS positioning point in sub-zone cf., ..., the distance between the n-1h and the nth GNSS positioning point, is the Pt" nth in the space-time sub-region, ..., the time stamp of the nth GNSS positioning point; the data in each spatial-temporal sub-zone is not filtered and constitutes a set 1-7: for subsequent processing.

52) Time-smooth sampling scheme for speed information. Specify the length of the time segment, set the upper limit of the number of data pieces in the same time segment search for the velocity data in each time segment in a spatial-temporal sub-zone: if the number of speed data in the time segment exceeds the upper limit the data of the upper limit is randomly used for subsequent processing. The implementation method is d + 2,; + ...+ d i." to calculate the travel speed of each vehicle in the special-temporal sub-zone v, -- * where 6/12...cifia, is the distance of the I' and rd GNSS positioning point in sub-zone..., the distance between the n-1h and the nth GN SS positioning point, ti...ta is the ls,..., nth in the space-time sub-region, ..., the time stamp of the nth GNSS positioning point; specify the length of time segment as i", the maximum number of fragment data at the same time asp,",,.; searching for speed data in the in' time segment of a special-temporal sub-zone, if the number of speed data in the time segment exceeds the upper limit pmaa, the p max pieces of data are randomly added to V. and used for subsequent processing.

The step 6) may specifically adopt one of the following methods: 61) Simple historical trajectory data fusion method. Regard the historical data under the condition of nonnal condition as a whole, the traffic feature model establishment and parameter estimation. The method uses a finite hybrid model to establish a traffic feature model and perform parameter estimation. Specifically, one of the following three options can be used: 611) Mixed Gaussian model of fixed composition This scheme uses a mixed Gaussian model with a fixed component quantity K to describe the probability distribution of vehicle speed. The number of components is manually specified according to the distribution pattern of the vehicle speed in a typical case. in order to ensure the reliability of the probability distribution, the component number K cannot be too small. Generally, K-4-6 is acceptable.

612) Mixed Gaussian model with variable composition This program uses a model-based evaluation method to select the appropriate number of components, as fol lows: The possible maximum component quantity K is determined, and the mixed Gaussian model of n=1, 2, K components is parameterized separately; for K models, the best model is determined by Bayesian Information Criterion (BIC). The maximum number of components is generally selected according to the accuracy requirements, but it must be noted that the more components; the slower the expectation maximization algorithm converges. The maximum number of components selected here is K=5, which requires calculation: J, (ye) = /id (ye) n E 11,2. ,51, totally 5 mixed model. At the same time, calculate the BIC of the five models, which is defined as: BIC = -21riL + k Inn, where L is the maximum likelihood function value, k is the number of parameters in the model, and n is the total amount of data.

After that, the minimum hybrid model of BIC is selected, and its parameter vectors rl, ft, a are recorded, where ri is the proportional vector occupied by each sub-component in the historical traffic feature model, It is the mean vector of each sub-component in the historical traffic feature model, and a is the standard deviation vector of each sub-component in the historical traffic feature model, as the feature of the present temporal sub-zone. The morphology of the density curve of the hybrid model is shown in Figure 6.

613) Finite mixed model with variable composition and distribution type This scheme uses the same model-based evaluation method as 612), but the distribution morphology and the number of sub-components are variable. The detailed steps are shown as follows: SelectMkinds of probability distribution models as the distribution types of sub-components, including but not limited to: normal distribution, gamma distribution, Weibull distribution. When using a normal distribution, the sub-distribution function takes: 1 (-1'. -Pu)2 f, (v, ) _ exp 2 27o- 2cr ^ 0 i When using a gamma distribution, the sub-distribution function takes: f (12)= e. ,1-(»;) PI'Rec,) When using a Weibull distribution, the sub-distribution function takes: ( \ -1 k v.

Assuming that the distribution types of all sub-components of the hybrid model are the same, detemfine the number K of possible maximum components. For the selection of NI seed component distribution type and K component quantity, a total of NI.K combinations were formed, and the BIC values were calculated respectively, and the model with the smallest BIC was taken as the best model.

62) Situational historical trajectory data classification method. According to temperature, precipitation, visibility and traffic control measures, the historical data without traffic incidents are divided into different categories, and models and parameter estimates are established. The implementation method is as follows: According to the temperature, precipitation, visibility and traffic control measures, the traffic environment is divided into 5-8 categories, and the historical data is classified into the above categories according to the different traffic environments corresponding to historical data. For each category, the processing as described in 5) is performed separately, thereby establishing a mapping relationship R(E--T), where E is the traffic environment and T is the traffic situation.

63) Historical data clustering method. For the historical data, the difference quantitative representation of different spatial-temporal regions is obtained by comparison between the time and space sub-regions, and the quantized differences are used for clustering. Temperature, precipitation, visibility and traffic control measures are used as characteristic factors to perform multiple Logit regressions to establish a mapping relationship between traffic environment and categories. See Figure 4 for the implementation process. The implementation steps are as follows: 631) According to the method described in 5), atraffic feature model is established and parameter estimation is performed.

632) According to the previous parameter estimation results of the (mite-mix model, the probability density function/Mx) of the travel speed distribution corresponding to the spatial-temporal subzone on different dates is written. The parameters of mixed Gaussian model, as an example, is shown as follows: P (K)=IO *_f(ve:11;,a1), j-1 where K represents the number of subcomponents of the travel speed distribution, t7 represents the proportion of a subcomponent in the travel speed distribution, p represents the mean value of a subcomponent in the travel speed distribution, and a represents the standard division of a subcomponent in the travel speed distribution.

633) Calculate the Jensen-Shannon divergence between the two distributions dl: d =IND(P 0)=-1D(P111)+ -1D(0, 2 2 -where P and Q are two different probability distributions, M I =1(P + 0) D is the Kullback-Le bler 2 -divergence: P(KK) D(PH Q.) = XPO'Olog 0(5) In the case of a finite mixed model, the value cannot be explicitly expressed, but Monte Carlo sampling method can be used to approximate the value. The calculation method is shown as follows: DMC(_f H g) = 1 (x,) log n(f II g) n g(1;) where DN1C represents the Kullback-Leibler divergence approximated by Monte Carlo sampling, and f and g represent any two distribution functions.

634) The divergence between the two distributions is expressed as a distance matrix: ii d," I) = The matrix satisfies 6/=4 and du=0(1=j).

635) The distance matrix is used as the input of the K-Medoids algorithm to obtain clustering results and index the categories.

636) Taking the category index as the response variable, the traffic environment data (including temperature, precipitation, visibility, etc.) is used as an independent variable to perform multiple logit regression to obtain the mapping relationship R (E-..-7) between the traffic environment F and the traffic situation category T. 637) The same type of data is aggregated, and the hybrid model is re-established with the new data set after aggregation, and the parameter estimation is performed to obtain the final historical traffic characteristic data set.

The step 7) may specifically adopt one of the following methods: 71) Simple real-time data processing. This method is implemented simultaneously with 61). The real-time traffic data is modeled and parameter estimation is performed to obtain the characteristic function of the current traffic condition. The implementation steps of the method are exactly the same as 61), except that the data used is real-time traffic data.

72) Classification processing. This method is carried out simultaneously with 62) or 63). Obtain the characteristic function of the traffic condition, and obtain the current information such as temperature, precipitation, visibility, traffic control measures, etc., and judge the current traffic condition category. See Figure 5 for the implementation process. The implementation steps are as follows: 721) Calculate the travel speed in the spatial-temporal sub-zone, which constitutes the overall real-time travel speed Eccr; 722) Establish a trip speed probability distribution model pr, )= ) and carry out c=1 parameter estimation; 723) The current traffic environment data (including temperature, precipitation, visibility, etc.) is taken as an input parameter, and the category T of the current traffic situation is obtained by using the mapping relationship R(E-7).

The step 8) specifically includes the following steps: 8I) When step 72) is adopted, historical traffic characteristic data under the category is located according to the category T of the current traffic situation, otherwise no processing is performed; 82) Calculate the difference between the two speed distributions according to the description parameters nrr, nn, a,1 of the current traffic characteristics and the description parameters lb it, a of the historical traffic characteristics: dill " , art), (q, It, a)] = H P). Among them, s, is the proportional vector of each sub-component in the real-time traffic feature model, Int is the mean vector of each sub-component in the real-time traffic feature model, a,, is the standard deviation vector of each sub-component in the real-time traffic feature model; 11 is the proportional vector of each sub-component in the history traffic feature model, ft is the mean vector of each sub-component in the historical traffic feature model, and a is the standard deviation vector of each sub-component in the historical traffic feature model. When historical traffic characteristics and real-time traffic characteristics (i.e., historical travel speed distribution and real-time travel speed distribution) are similar, a smaller Jensen-Shannon divergence value will be obtained, that is, the difference between the two value is smaller; when the difference between the historical traffic characteristics and the real-time traffic characteristics is large, a larger Jensen-Shannon divergence value will be obtained, that is, the difference between the two value is large, that is, the probability of abnormality is large, sec Figure 7 for details.

The step 9) specifically includes the following steps: 91) Normalize the difference in speed distribution of each spatio-temporal subzone to a normalized value cte of 0-1: duff -min(diff) max (duff) -min (Off) 92) Calculate the traffic anomaly index of each spatial-temporal subzone*-ax10; 93) The location of the region with an anomaly index higher than 5 is projected onto the electronic map, and is released to the public in the form of a smart mobile device APP, so that the driver can avoid potential congestion points and improve the efficiency of urban road traffic.

The step 10) specifically includes the following steps: 101) Calculate the false negative rate of traffic anomaly: a, = x100% 102) Calculate the false positive rate of traffic anomaly: a2 = X I 00% In the above two formulas, n1 is the total number of missed events per unit time, n, is the total number of false positive events per unit time, and is the total number of abnormal times actually occurred in the unit time.

The present invention has the following advantages over similar technologies in the same field: ( I) The method can make -HI use of existing floating vehicle operation data (GNSS trajectory data) detect historical traffic state changes through historical traffic feature extraction and real-time traffic situation analysis, and realize real-time, low-cost and intelligent urban road traffic anomaly events detection; (2) The method takes the probability distribution of traffic characteristic parameters as the description of traffic characteristics, so that the characteristics reflected are more comprehensive, avoiding the use of a single index to characterize the one-sidedness and instability of traffic features, and the reliability of detection is higher; (3) Aiming at the characteristics that the traffic characteristics arc affected by the traffic environment (such as weather conditions), a cluster-multiple Logit regression joint algorithm is introduced to establish the mapping relationship between traffic environment characteristics and traffic situation categories.

(4) According to the test by actual data, the urban road traffic anomaly detection technology based on floating car data proposed by the present invention can realize the detection of abnormal events with high accuracy, the detection rate exceeds 90%, the false negative rate is less than 15%, and the false positive rate is less than 20%. It has achieved good detection results and can be applied to intelligent management and service of urban traffic.

Drawings The specific content and advantages of the present invention will become apparent and readily understood in conjunction with the following drawings: Figure 1 shows a schematic diagram of the components and basic principles of the present invention; Figure 2 is a schematic view showing the overall flow of the present invention in the implementation process; Figure 3 is a schematic diagram showing an implementation manner of a fast map matching algorithm of the present invention; Figure 4 is a schematic flow chart showing a historical traffic feature extraction scheme implemented by the present invention; Figure 5 is a schematic flowchart diagram of a real-time traffic feature extraction scheme implemented by the present invention; Figure 6 is a schematic diagram showing the morphology of a Gaussian mixture model probability distribution; Figure 7 shows a measurement of the difference in the comparison between historical traffic characteristics and real-time traffic characteristics.

Specific implementation In order to clarify the objects, technical solutions, and advantages of the present invention more clearly, the specific embodiments of the present invention are described in detail below.

As shown in Figure 1, the overall system architecture of the present invention includes an onboard GNSS track recorder, data center, GNSS satellites, and communication system carried by a floating vehicle. The GNSS here includes GPS, GLONASS, GALILEO, Beidou, IRNSS, QZ SS and any similar navigation satellite positioning system. GNSS track recorders equipped with floating cars, buses, etc., record the position information of vehicles al various points in time with a certain sampling frequency f (general requirements 7("A-1.1Hz), and send the location information to the data center in real time through the GPRS mobile communication network (other wireless network communication technology such as WCDMA and TD-LTE can also be adopted, but the cost will be correspondingly improved). The data center establishes a historical road traffic feature database through data preprocessing and data fusion, and establishes a real-time traffic feature database for the recently received real-time data; and determines whether the current traffic feature is abnormal through the mapping relationship between the historical database and the real-time database. Finally, the incident information is visualized through the processing terminal and a traffic incident report is generated.

The overall process of the scheme is shown in Figure 2. it includes the steps of collecting and storing GNSS trajectory data, establishing spatial-temporal subzones, historical traffic feature extraction, real-time traffic feature extraction, and incident identification. Collecting and storing GNSS trajectory data is the data foundation of the whole scheme. Due to the huge amount of data, a distributed storage scheme should be adopted. For distributed storage, there are mature technologies, which are not the content of the present invention. The basic assumption of establishing a spatial-temporal subzone is that it has the same traffic characteristics in a certain area and a specific time period. This assumption is universally applicable after long-term observation. Historical traffic feature extraction, the principle is to use GNSS trajectory data, calculate the travel speed, use a large number of travel speed data in the same spatial-temporal subzone, establish a probability distribution model of vehicle speed, and estimate the parameters, with a small number of parameters to characterize the traffic characteristics. For real-time traffic feature extraction, the principle is to process and analyze the speed data in the current time period, and also establish the current vehicle speed probability distribution model. The abnormality identification is to use the difference measurement index to judge the degree of change of the real-time feature compared with the historical feature and determine whether a traffic incident occurs according to whether it reaches the threshold.

According to the combination of the embodiments of the invention, the implementation is given below.

Embodiment 1 Step 11. Implement equidistant spatial-temporal division method and detemnne the segment size of the time dimension, the time segment span is a fixed value, usually taking 30min as a time segment; determine the segment size of the spatial dimension, the spatial segment span is a fixed value, and usually takes a spatial grid of 200m x200m as a spatial segment; Step 12. Implement data pre-processing. GNSS positioning data should be processed by data cleaning, data integration, data conversion and data reduction, in order to improve structural level of data.

Step 13. Divide the space area to be processed into a grid of a certain size, then the range of each grid area can be expressed as = ) [Xr X711),Y5 [Yr, Yr determine the grid area where the positioning point is located, and searching for the road section where the positioning point is located by using the distance and the azimuth angle; Search for the road segment closest to point A, and take 69=2.5°. if the difference between the azimuth angle of the point A and the direction of the road section ij is less than the threshold (59, i.e. <8 then the matching process is completed; if PA -6101 <150 is not satisfied, delete the road segment if in the search space and continue to search for other road segments until the conditions are met; use the straight-line equation of the road segment (if it is a curved road segment, it is roughly split into straight lines), calculate the projection coordinates of the GNSS positioning point on the road segment, and reduce the error caused by the GNSS positioning drift. The specific method is descript as follows.

Determine the straight-line equation of the road segment ij (if the road segment is a curve, it is divided into several straight-line segments): = k(x -k= y; -y, x-y -y, k, kV kY, ±k2 ± X4

-

P k2 +1

where the slope is: The projection line equation is: Solve the projected coordinate P as: y4+ j + la,

P k2 +1

After the map matching process, combined with the timestamp data of the coordinates of the positioning point, the positioning point is matched to the space-time sub-region.

Step 14. The overall travel speed data of each floating vehicle in a time and space sub-zone constitutes the whole. The implementation method is to calculate the travel speed of each vehicle in the spatial-temporal d, , + + ... + d," sub-zone dr: = d,2 -3 " , where dm. etfri," is the distance of the 1" and 2111 GNSS positioning t;2-ti point in sub-zone dr, ..., the distance between the n-lth and the nth GNSS positioning point, is the 1", in the space-time sub-region, ..., the time stamp of the GNSS positioning point; the data in each spatial-temporal sub-zone is not filtered and constitutes a set 127, for subsequent processing.

Step 15. Regard the historical data under the condition of normal condition as a whole, the traffic feature model establishment and parameter estimation. The method uses a finite hybrid model to establish a traffic feature model and perform parameter estimation. Take the component quantity K=5 and the mixed Gaussian models with the component quantity n=1, 2, K components are parameterized separately; for the K models, the best model is determined by Bayesian information Criterion (Bin. Calculate: f(v,)=EriuMv;) n E{1,2,...,5} totally 5 mixed model. At the same time, calculate the BIC of the five models, which is defined as: BIC=-21nL+k*lnn where L is the maximum likelihood function value, k is the number of parameters in the model, and n is the total amount of data After that, the minimum hybrid model of BIC is selected, and its parameter vectors q, p, a are recorded, as the feature of the present temporal sub-zone.

Step 16. The real-time traffic data is modeled, and parameter estimation is performed to obtain the characteristic function of the current traffic condition. The method is the same as Step 15. Record the parameter vectors in. prh Step 17. Calculate the difference between the two speed distributions according to the description parameters gn, an of the current traffic characteristics and the description parameters it p, a of the historical traffic characteristics: cliff [(11," it, " ),(q,p, a) 1= .15D(P" H P).

Step 18. Normalize the difference in speed distribution of each spatial-temporal subzone to a normalized value tk of 0-1: difJ-min (diff) a. -max (dig) -mM @of) Calculate the traffic anomaly index of each spatial-temporal subzone 2.4.--a,rx 10.

Embodiment 2 Step 21. Implement equidistant spatial-temporal division method and determine the segment size of the time dimension, the time segment span is a fixed value, usually taking 30min as a time segment; determine the segment size of the spatial dimension, the spatial segment span is a fixed value, and usually takes a spatial grid of 200m x 200m as a spatial segment; Step 22. implement data pre-processing. GNSS positioning data should be processed by data cleaning, data integration, data conversion and data reduction, in order to improve structural level of data.

Step 23. Divide the space area to be processed into a grid of a certain size, then the range of each grid area can be expressed as A, = {(x"ys) xs. G[)C, c [y "Y,1)}; determine the grid area where the positioning point is located, and searching for the road section where the positioning point is located by using the distance and the azimuth angle; Search for the road segment closest to point A, and take 60=2.5°. If the difference between the azimuth angle of the point A and the direction of the road section if is less than the threshold 80, i.e. 16' -0 -l<8 then the matching process is completed; if PA -001<80 is not satisfied, delete the road segment ij in the search space and continue to search for other road segments until the conditions are met; use the straight-line equation of the road segment (if it is a curved road segment, it is roughly split into straight lines), calculate the projection coordinates of the GNSS positioning point on the road segment, and reduce the error caused by the GNSS positioning drift. The specific method is descript as follows.

Determine the straight-line equation of the road segment ij (if the road segment s a curve, it is divided into several straight-line segments): where the slope is: y -v, = k(x - k Y, Y, The projection line equation is: ' I Solve the projected coordinate P as: - -x,21 k, k3;, +ex, x; x +x4 kz +1 y + y, + acv d-fix, YR r k-+1 After the map matching process, combined with the timestamp data of the coordinates of the positioning point, the positioning point is matched to the space-time sub-region.

Step 24. Calculate the travel speed of each vehicle in the special-temporal sub-zone j: + d + ...+ d v, = where is the distance of the PI and 2' GNSS positioning point in sub-zone s, ..., the distance between the n-11'and the nth GNSS positioning point, is the 1s nth in the space-time sub-region, ..., the time stamp of the GNSS positioning point; specify the length of time segment as 1k,, the maximum number of fragment data at the same time as p",,,r; searching for speed data in the it time segment of a special-temporal sub-zone, if the number of speed data in the time segment exceeds the upper limit p.x, the p.m pieces of data are randomly added to V. Step 25. Regard the historical data under the condition of normal condition as a whole, the traffic feature model establishment and parameter estimation. The method uses a finite hybrid model to establish a traffic feature model and perform parameter estimation. Take the component quantity K=5 and the mixed Gaussian models with the component quantity n=1, 2, K components are parameterized separately; for the K models, the best model is determined by Bayesian Information Criterion (BIC). Calculate: /(v.) = ./0. (v; ) n E 2,.. 51 =I totally 5 mixed model. At the same time, calculate the BIC of the five models, which is defined as: BIC = -21nL+k*Inn where L is the maximum likelihood function value. It is the number of parameters in the model, and n is the total amount of data After that, the minimum hybrid model of BIC is selected, and its parameter vectors q, it, a are recorded, as the feature of the present temporal sub-zone.

According to the previous parameter estimation results of the finite-mix model, the probability density function p;(x) of the travel speed distribution corresponding to the spatial-temporal subzone on different dates is written. The parameters of mixed Gaussian model, as an example, is shown as follows: P, =1"0, f, a 1=1 Calculate the Jensen-Shannon divergence between the two distributions di": = D(P Q)= -1D(P M) + -1D(C.),11,1), 2 2 where P and Q are two different probability distributions, M = -(P +O) , D is the Kullback-Leibler divergence:

POO

D(P e) = Po:01°g 0,0 k-I In the case of a finite mixed model, the value cannot be explicitly expressed, but Monte Carlo sampling method can be used to approximate the value. The calculation method is shown as follows: I f(x) Ric(/' H log. >N./ II g) n g(x,) The divergence between the two distributions is expressed as a distance matrix: d" D= * : The matrix satisfies do=do and dif-0(i-j).

The distance matrix is used as the input of the K-Medoids algorithm to obtain clustering results and index the categories.

Taking the category index as the response variable, the traffic environment data (including temperature, precipitation, visibility, etc.) is used as an independent variable to perform multiple logit regression to obtain the mapping relationship R 1) between the traffic environment E and the traffic situation category T. The same type of data is aggregated, and the hybrid model is re-established with the new data set after aggregation, and the parameter estimation is performed to obtain the final historical traffic characteristic data set.

Step 26. Obtain the characteristic function of the traffic condition, and obtain the current information such as temperature, precipitation, visibility, traffic control measures, etc., and judge the current traffic condition category.

Calculate the travel speed in the spatial-temporal sub-zone, which constitutes the overall real-time travel speed -17,,,,,q; Establish a trip speed probability distribution model pn(1,,)=Iri, * ,6,) and cam( out parameter estimation; The current traffic environment data (including temperature, precipitation, visibility, etc.) is taken as an input parameter, and the category T of the current traffic situation is obtained by using the mapping relationship R(E--7).

Step 27. Historical traffic characteristic data under the category is located according to the category T of the current traffic situation; calculate the difference between the two speed distributions according to the description parameters t,,, an of the current traffic characteristics and the description parameters q, II, a of the historical traffic characteristics: duff,g," a, ) = ISTI(Pr, H P).

Step 28. Normalize the difference in speed distribution of each spatio-temporal subzone to a normalized value tk of 0-1: IEEE -min(diff) a, -max Off) -min (cliff Calculate the traffic anomaly index of each spatial-temporal subzone A;2=ae10.

Embodiment 3 Step 31. Implement non-equidistant spatial-temporal division method. When the peak hour traffic flow rate is greater than or equal to 1000 vehicles per hour, take 30min time segment and 200m x200m space segment; when peak hour traffic flow rate is less than 1000 vehicles per hour, take 30min time segment and 400mx400m space segment.

Step 32. Implement data pre-processing. GNSS positioning data should be processed by data cleaning, data integration, data conversion and data reduction, in order to improve structural level of data.

Step 33. Divide the space area to be processed into a grid of a certain size, then the range of each grid area can be expressed as 4 = {(xs, ys.) x, E 1), ys [Y, 'Yr i)}-Denote the floating vehicle GNSS data acquisition frequency as fu=1/to; denote the point adjacent to A in time as P(tA-to), then P(tA+to) can be defined as the 1-adjacent point of A, and P(2.4-2t0), P(t4+2to) can be defined as the 2-adjacent point of A, and so on, then P(ki-k-to), P(214+kto) can be defined as the k-adjacent point of A. When fo. <1Hz, k=1 or 2. Take the road segment with the smallest distance from the k-adjacent point and A, and calculate ot, the mean value of the direction angle of the k-adjacent point of A and A. If 1974. -ed< 59 is met, the matching process is completed. If not, search in other road segments until the condition 1/9;,* -Of/I< 3, is met.

Usc the straight-line equation of the road segment (if it is a curved road segment, it is roughly split into straight lines), calculate the projection coordinates of the GNSS positioning point on the road segment, and reduce the error caused by the GNSS positioning drift. The specific method is descript as follows.

Determine the straight-line equation of the road segment U (if the road segment is a curve, it is divided into several straight-line segments): where the slope is: y-Y,=4)c- 1 k= Y Y, The projection line equation is: Y YA A124) Solve the projected coordinate P as: + X -x +xA, k2 +1 e y4+ y, + la,

Y P +1

After the map matching process, combined with the timestamp data of the coordinates of the positioning point, the positioning point is matched to the space-time sub-region.

Step 34. Calculate the travel speed of each vehicle in the special-temporal sub-zone + +...+ d v, - " , where 42...(4,-0, is the distance of the 1st and rd GNSS positioning point in sub-

-

zone 5, the distance between the n-1u and the nth GNSS positioning point, ti...t" is the 1s nth in the space-time sub-region, the time stamp of the nth GNSS positioning point; specify the length of time segment as 1", the maximum number of fragment data at the same time as p","-; searching for speed data in the ith time segment of a special-temporal sub-zone, if the number of speed data in the time segment exceeds the upper limit pnth, the p"", pieces of data are randomly added to V. Step 35. Regard the historical data under the condition of normal condition as a whole, the traffic feature model establishment and parameter estimation. The method uses a finite hybrid model to establish a traffic feature model and perform parameter estimation. Take the component quantity K=5 and the mixed Gaussian models with the component quantity n=1, 2, Kcomponents are parameterized separately; for the K models, the best model is determined by Bayesian Information Criterion (BIC). Calculate: = n E totally 5 mixed model. At the same time, calculate the BIC of the five models, which is defined as: BIC'=-21nL+k*hin where L is the maximum likelihood function value, k is the number of parameters in the model, and n is the total amount of data After that, the minimum hybrid model of BIC is selected, and its parameter vectors 11, it, a are recorded, as the feature of the present temporal sub-zone.

According to the previous parameter estimation results of the finite-mix model, the probability density function p;(x) of the travel speed distribution corresponding to the spatial-temporal subzone on different dates is written. The parameters of mixed Gaussian model, as an example, is shown as follows: f(ve-*11,-(Tj)* j-1 Calculate the Jensen-Shannon divergence between the two distributions cl": = ISD(P 0)= -I D(P A/1)+ -I 1.)(() d 2 2 where P and Q are two different probability distributions, M = 1(P +0) , I) is the Kullback-Leibler divergence: P(X) D(P110= IP(X0113g /-Q(x,c) In the case of a finite mixed model, the value cannot be explicitly expressed, but Monte Carlo sampling method can be used to approximate the value. The calculation method is shown as follows: Dme(f H g) = 1 -I log(. >D(/' II g) n =I g(x,) The divergence between the two distributions is expressed as a distance matrix: D= The matrix satisfies did, and dii=0(i=j).

The distance matrix is used as the input of the K-Medoids algorithm to obtain clustering results and index the categories.

Taking the category index as the response variable, the traffic environment data (including temperature, precipitation, visibility, etc.) is used as an independent variable to perform multiple logit regression to obtain the mapping relationship R (E-T) between the traffic environment E and the traffic situation category T. The same type of data is aggregated, and the hybrid model is re-established with the new data set after aggregation, and the parameter estimation is performed to obtain the final historical traffic characteristic data sct.

Step 36. Obtain the characteristic function of the traffic condition, and obtain the current information such as temperature, precipitation, visibility, traffic control measures, etc., and judge the current traffic condition category.

Calculate the travel speed in the spatial-temporal sub-zone, which constitutes the overall real-time travel speed C.;,,,t; Establish a trip speed probability distribution model = ri * ff (Yr Jr; p f,o-i) and cam/ out parameter estimation; The current traffic environment data (including temperature, precipitation, visibility, etc.) is taken as an input parameter, and the category T of the current traffic situation is obtained by using the mapping relationship R(E-29) . Step 37. Historical traffic characteristic data under the category is located according to the category T of the current traffic situation; calculate the difference between the two speed distributions according to the description parameters lin, run, an of the current traffic characteristics and the description parameters q, it, a of the historical traffic characteristics: Off [(qgr, , a), (Th it, a)] = JSD(P,., 11P).

Step 38. Normalize the difference in speed distribution of each spatio-temporal subzone to a normalized value tk of 0-1: dill; - (cl cle max (deli) - (di,iy) Calculate the traffic anomaly index of each spatial-temporal subzonc.A.-ti,rx 10.

Claims (9)

1. Claims 1. A method for detecting traffic anomaly of urban road with equidistant spatial-temporal division, comprising the following steps: 1) Creating spatial-temporal subzone: divide one day into several time segments; each time segment is defined as a time subzone; divide one implemented area into several spatial segments; each spatial segment is defined as a spatial subzone; intersection of one time subzone and one spatial subzone is defined as a spatial-temporal subzone; 2) Preprocessing of historical trajectory data: processing GNSS positioning historical data of floating vehicles into sampled vehicle speed data of historical trajectory; Preprocessing of real-time trajectory data: processing GNSS positioning real-time data of floating vehicles into sampled vehicle speed data of real-time trajectory; 3) Historical trajectory data analysis and feature extraction: Using said sampled vehicle speed data of historical trajectory, establishing a historical travel speed probability distribution, to obtain a historical traffic feature model Ph: by means of taking said sampled vehicle speed data of historical trajectory without traffic anomaly as a whole, and using a finite mixed modelling method to establish said historical traffic feature model Ph; Real-time trajectory data analysis and feature extraction: Using said sampled vehicle speed data of real-time trajectory, establishing a real-time travel speed probability distribution, to obtain a real-time traffic feature model Pri; 4) Difference calculation: Jensen-Shannon divergence is used to measure difference between said historical travel speed probability distribution and said real-time travel speed probability distribution; 5) Quantitative characterization of anomaly severity: Using said difference between historical and real-time traffic features, calculate an anomaly index of traffic conditions; 6) Accuracy evaluation: evaluate accuracy of said anomaly index; wherein said time segment size takes 30min as one time segment; said spatial segment size takes a spatial grid of 200m x200m as one spatial segment.
2. 2. A method for detecting traffic anomaly of urban road according to claim 1, wherein said preprocessing of the historical trajectory data in step 2) comprises: 2a) Data structuralizing: GNSS positioning data is processed by data cleaning, data integration, data conversion and data reduction, in order to get structuralized GNSS positioning historical data; 2b) Map matching: combined with urban road network data, project GNSS positioning points to an urban map through map matching algorithm to establish a matching relationship between said GNSS positioning points and road segments, and correct errors caused by positioning drift; 2c) Vehicle speed calculation and sampling of historical trajectory: calculate traffic operation characteristic parameters according to said structural i zed GNSS positioning history data, obtain vehicle speed data of historical trajectory, and perform data sampling on said vehicle speed data of historical trajectory to obtain sampled vehicle speed data of historical trajectory.
3. 3. A method for detecting traffic anomaly of urban road according to claim 2, wherein said map matching in step 2b) comprises the following sub-steps: 2b1) Divide space area to be processed into grid areas of a certain size, the range of each grid area is expressed as = {(x"y5) x,. E [x" x" ), ys [y"y, ,,)} , each grid area contains several road sections, a set of said road sections is denoted as R,, each road section of said set of road sections Kris denoted as ij, then assign an index to each road section; 2b2) Searching for road section if where a positioning point A is located in said set of road segments Rs according to its distance and azimuth angle; 2b3) GNSS positioning point linear projection method is used to calculate projection coordinates of GNSS positioning point on said road section.
4. 4 A method for detecting traffic anomaly of urban road according to claim 3, wherein sub-step 2b2) adopts one of the following methods: 2b21) Single point matching scheme: searching for a road section that is closest to said point A, if difference between direction angle of travel of said point A (OA) and direction angle of said road section 0 (Of) is less than a threshold (Se),16j A 61,11<6;i' then said map matching is completed; if condition 6jfil< is not met, search in other road sections until said condition is met; 2b22) Point sequence matching scheme: denote floating vehicle GNSS data acquisition frequency asin=1/to; one point adjacent to A in time as 1'(t4-to) or P(tt+to) is defined as 1-point-adjacent to A, and P(t1-26) or P(ti+2to) is defined as 2-point-adjacent to A, and so on, then P(tt-kto) or P(ti-Ekto) is defined as k-point-adjacent to A; when fo<lHz, Ic=1 or 2, take said road section ij with the smallest distance from said k-point-adj acent to A to said point A, and calculate which is mean value of direction angles of said k-point-adjacent to A and said point A; if condition 1:9A, - 8, is met, then said map matching is completed; if not, search in other road sections until said condition 15A' -Bid <S. is met.
5. 5. A method for detecting traffic anomaly of urban road according to claim 2, wherein the vehicle speed calculation and sampling of the historical trajectory in step 2c) adopts one of the following methods: 2c1) Full sample scheme: calculate travel speed of each vehicle in a spatial-temporal sub-zone d + ...+ dy, zone: , , = , where is distance between the 1st and 2nd -tt GNSS positioning point in subzone..., distance between the n-lth and the nth GNS S positioning point, ti...!" is the 1",..., nth time stamp of GNSS positioning point in said spatial-temporal subzone, respectively; speed data in each spatial-temporal subzone is not filtered out, and constitutes a speed set Vi; for subsequent processing; 2c2) Time-smooth sampling scheme: calculate travel speed of each vehicle in said special- 2+ dn, + ...+ temporal subzone v. - , where d,1_,," is distance between the 1" and 2"d GNSS positioning point in subzone... distance between the n-lth and the nth GNSS positioning point, is the nth time stamp of GNSS positioning point in said spatial-temporal subzone, respectively; then define time segment size as 1,, define an upper limit, which is maximum number of data pieces in said time segment, as p,,,,,x; then search for speed data in the time segment of said special-temporal subzone, if the number of speed data in said time segment exceeds said upper thenp,,,,,N pieces of speed data are randomly selected and added to a speed set V,: for subsequent processing.
6. 6 A method for detecting traffic anomaly of urban road according to any of previous claims, wherein said finite mixed modelling method in step 3) adopts one of the following methods: 3a) Mixed Gaussian modelling with fixed composition: using a mixed Gaussian model with a fixed component number K to describe probability distribution of vehicle speed; said component number K is from 4 to 6; then real-time traffic data is modeled and parameter estimation is performed to obtain a characteristic function of current traffic condition; 3b) Variable Mixed Gaussian modelling: using variable component numbers, or variable component numbers and variable distribution types; using classification processing to obtain a characteristic function of current traffic condition, and obtaining information of temperature, precipitation, visibility, and traffic control measures, then determining category of said current traffic condition.
7. 7. A method for detecting traffic anomaly of urban road according to claim 6, wherein 3b) comprises one of the following sub-methods: 3b 1) Mixed Gaussian modelling with variable composition numbers: using a model-based evaluation method to select component number, determine maximum component number K and mixed Gaussian model of n=1, 2, ..., K components is parameterized separately; for K models, the best model is determined by Bayesian Information Criterion; 3b2) Finite mixed modelling with variable composition numbers and distribution types: both distribution of subcomponents and component number are variables; Wherein said classification processing comprises: 3b3) Calculate travel speed in sai spatial-temporal subzone, which constitutes real-time travel speed set V (,,t; Ii 3b4) Establish a travel speed probability distribution model P,T(; ,r) I, * I p cri), where K represents the number of subcomponents of real-time traffic feature, II represents the proportion of a certain sub-component of real-time traffic feature, ft represents the mean of a certain sub-component of real-time traffic feature, and a represents the standard deviation of a sub-component of real-time traffic feature; implement parameter estimation to obtain description parameters mn, tun. an of current real-time traffic features; wherein, tin is the proportional vector of each subcomponent in said real-time traffic feature model, grt is the mean vector of each subcomponent in said real-time traffic feature model, an is the standard deviation vector of each sub-component in said real-time traffic feature model; rl is the proportional vector of each sub-component in historical traffic feature model, ti is the mean vector of each sub-component in said historical traffic feature model, and a is the standard deviation vector of each sub-component in said historical traffic feature model; 3b5) Current traffic environment data, including temperature, precipitation, and visibility, is taken as an input parameter, and category T of current traffic situation is obtained by using mapping relationship R(/ 7), where E represents said current traffic environment data
8. 8. A method for detecting traffic anomaly of urban road according to any of claims 1 to 5, wherein the step 4) comprises: 4a) According to category T of current traffic situation, locate historical traffic feature data under said category T; if there is no category division, it is not necessary to distinguish; 4b) Based on description parameters of current description parameters lin, prt, an and description parameters in historical traffic feature model sl, it, a to calculate difference between two speed distributions: dill,f) (it it. (7)1 = JSD(p, ph), where lin is the proportional vector of each sub-component in said real-time traffic feature model, grt is the mean vector of each sub-component in said real-time traffic feature model, an is the standard deviation vector of each sub-component in said real-time traffic feature model; sl is the proportional vector of each sub-component in said historical traffic feature model, g is the mean vector of each sub-component in said historical traffic feature model, and a is the standard deviation vector of each sub-component in said historical traffic feature model.
9. 9 A method for detecting traffic anomaly of urban road according to any of claims 1 to 5, wherein the step 5) comprises: 5a) Normalize differences among speed distributions of all spatial-temporal subzones to a normalized value gi of 0-1: a =cliff -min(diff)1/ max (diff)- @tiff)] ; wherein, cliff refers to difference between speed distributions of any two spatial-temporal subzones; 5b) Calculate traffic anomaly index of each spatial-temporal subzone: AL;=a0 x 10.