GB2588556A - A method for detecting traffic anomaly based on travel time distribution 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 travel time 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 time probability distributions. Jensen-Shannon divergence is then used to measure a difference between said historical and real-time travel time probability distributions and an anomaly index of traffic conditions is determined from said difference. The accuracy of said anomaly index is then evaluated.

Description

A Method for Detecting Traffic Anomaly Based on Travel Time Distribution with Equidistant Spatial-temporal Division This is a divisional application of the following 'parent' application: Application Number: GB1905907.0 Application Source: NJ'] Filing Date: 30 December 2017 POT Application No.. PC1V1B2017/058536 Priority Claimed: 30 December 2016 in W1PO -PC7T1132016/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 or traffic anomaly. By GNSS positioning devices of the floating cars, the spatial-temporal position infommtion can be obtained. After data preprocessing, map matching and data fusion, the probability distribution of travel time for specific spatial-temporal region will be derived. According to the changes of time distribution, the city road traffic incidents can be effectively identified.

Technological background

Traffic incident 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 are the following representative techniques: 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 arid 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: I 0, 0: 10-0:20 arid so on. Each time segment is called a time subzone. Divide the implemented area into several spatial segments, e.g., the area of 121.58°E-121.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 incident 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 incident 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) = 17,P,(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 fommla is: BIC = -2111_1,-Hk-inn where L is the maximum of the likelihood function. K is the number of unknown parameters, and Nis 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(81x) is equal to the probability of a given parameter 0 after the variable X. That is, L(01x-)=P(X=x10).

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 Q. 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: I) 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 incident 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, 1RNSS 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 abnomml 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 die 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 is based on travel time of floating vehicles. The urban road traffic incident detection is based on the travel time. Therefore, data sampling refers to the sampling of travel time.

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 time. In the past, many models used a single indicator to represent the traffic characteristics of a certain road scction or region, such as the historical average speed. Although these methods are easy to implement, thc 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 time.

7) Real-time trajectory data analysis and feature extraction The so-called real-time trajectory data refers to the trajectory data 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 time 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 figure are based on it.

In recent years, in the development of anomaly detection teclmology, 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 detennine which features are 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 sonic 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 abnommlity.

The location of the anomaly is projected onto the electronic map and is publicly released by the smart mobile device APR 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 200mx200m 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 200mx200m 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 flow rate --when the peak hour traffic flow rate is greater than or equal to 1000 vehicles per hour, take 30min as one time segment and 200m x200m as one space segment; when peak hour traffic flow rate is less than 1000 vehicles per hour, take 30min as one time segment and 400m x400m as one 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 As = {(xs, JO xs E, Xr+,), E [y7, yr +1)1, each grid area contains several road sections, the set of these road sections can be denoted as /2, each road section of the abovementioned set of road section Rs can be denoted as and assign an index to each road section; 32) Determining the grid area where the positioning point is located, and searching for the section where the positioning point A is located in the set Rs 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 is less than the threshold, i.e. 161, -81 <8, , then the matching process is completed, the abovementioned threshold can be 2.5' 10° etc if9 A!, < 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 as /6=1//i); denote the point adjacent to A in time as P(ti-to), then P(ti-v/i)) can be defined as the 1-adjacent point of A, and P(tA-2t0),P(:4-12t0) can be defined as the 2-adjacent point of A, and so on, then P(14-kt0), P(171-Fkin) can be defined as the k-adjacent point of A. Whenfi<lHz, k=1 or 2. Take the road segment ü with the smallest distance from the k-adjacent point and A, and calculate I9A., the mean value of the direction angle of the k-adjacent point of A and A. If 15A, -Od <8 is met, the matching process is completed. If not, search in other road segments until the condition - <go is met; 33) 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): y -y, = k(x- ) Where the slope is: xi -The projection line equation is: Y 1, Yt = -T(x -x4) ky -ky, + +x Solve the projected coordinate P as: lc' +1 yr, - 'A + Y, IccA k 2 + 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 spatial-temporal sub-zone.

The step 5) may specifically adopt one of the following methods: 51) A Full sample scheme of travel time: All travel time data of the same path between the OD pairs constitutes the whole. The implementation method is to calculate the travel time of each vehicle in the spatial-temporal sub-zone, where t1.2...tn-1,n is die distance of the 1st and 2nd GNSS positioning point in sub-zone "E" the distance between the n-lth and the nth GNSS positioning point, tl...tn is the 1st,..., nth in the spatial-temporal sub-zone, ..., the time stamp of the nth GNSS positioning point; the data in each spatial-temporal sub-zone is not filtered and constitutes a set T4 for subsequent processing.

52) Time-smooth sampling method of travel time information: Specify the length of the time segment, set the upper limit of the number of segments of the same time; search for the travel time data in each time segment of a time-space sub-zone, if the number of travel time data in the time segment exceeds the upper limit, The data of the upper limit number is randomly used for subsequent processing; the implementation method is to calculate the travel time of each vehicle in the spatial-temporal sub-zone: /, =I. +19,3 +...+ç 1. where ti,2...62_,J2is the travel time between the first and second GNSS anchor points in the spatial-temporal sub-zone, ..., the travel time between the n-lth and nth GNSS anchor points; specify the time segment length //22, the upper limit of the number of segment data at the same time Ant"; search for the travel time data in the time segment of the i-1 th time segment in a time-space sub-zone, if the number of travel time data in the time segment exceeds the upper limit pmax, randomly take The Ana strip data constitutes a set 7'2, and is 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 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. 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 travel time. The number of components is manually specified according to the distribution pattern of travel time in a typical case. In order to ensure the reliability of the probability distribution, the component number K cannot be too small. Generally, K=42-26 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 follows: 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: fj.() (v) n E r-1 totally S mixed model. At the same time, calculate the BIC of the five models which is defined as: BIC = -2 hi L inn, where L is the maximum likelihood function value, kis 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 ij, t, a are recorded, where ii 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: Select A,/ kinds of probability distribution models as the distribution types of sub-components, including but not limited to: nomml distribution, gamma distribution. Weibull distribution. When using a normal distribution, the sub-distribution function takes: f(v) = cxp 2c When using a gamma distribution, the sub-distribution function takes: v -p)2 Mr' [(v,) = fc: e When using a Weibull distribution, the sub-distribution function takes: Assuming that the distribution types of all sub-components of the hybrid model are the same, determine the number K of possible maximum components. For the selection of M seed component distribution type and K component quantity, a total of Mt. 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-zones, and the quantized differences arc used for clustering. Temperature, precipitation, visibility and traffic control measures are used as characteristic factors to perfonu 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), a traffic feature model is established and parameter estimation is performed.

632) According to the previous parameter estimation results of the finite-mix model, the probability density function p;(4 of the travel time 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: = yq; a).

where K represents the number of subcomponents of the travel time distribution, q represents the proportion of a subcomponent in the travel time distribution, IA represents the mean value of a subcomponent in the travel time distribution, and a represents the standard division of a subcomponent in the travel time distribution.

633) Calculate the Jensen-Shannon divergence between the two distributions du: = ISD(1' -D(1' M)+ -111(0 M) , 2 2 where P and Q are two different probability distributions, Al =-(1' + Q), D is the Kullback-Le bier divergence: D(1' H o)=Ip(x)log) 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: Dmt(i II g)-1 -Dog f(x) >D(f 11g) n g(x) where Dme represents the Kullback-Leibler divergence approximated by Monte Carlo sampling, and/. and g represent any two distribution functions.

634) The divergence between the two distributions is expressed as a distance matrix: D = The matrix satisfies du-di, and do=0(i=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 (F-1) between the traffic environment E 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 sct.

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 time in the spatial-temporal sub-zone, which constitutes the overall real-time travel time 722) Establish a travel time probability distribution model prt(t,.")= Eq, * f ( o-),) and carry out f-i 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: 81) 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 travel time distributions according to the description parameters iv, pri, an of the current traffic characteristics and the description parameters ii, g. a of the historical traffic characteristics: cliff 1(ii",tt " , ") ,(11, Et, 05)1= JSD(P, HP) Among them, IF, is the proportional vector of each sub-component in the real-time traffic feature model, ftrt is the mean vector of each sub-component in the real-time traffic feature model, an is the standard deviation vector of each subcomponent in the real-time traffic feature model; n is the proportional vector of each sub-component in the history traffic feature model, tt 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 time distribution and real-time travel time 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, see Figure 7 for details.

The step 9) specifically includes the following steps: 91) Normalize the difference in travel time distribution of each spatial-temporal subzone to a normalized value aeof 0-1: dij4, -min (diff) a max (thio_min(diff) 92) Calculate the traffic incident index of each spatial-temporal subzone A e-ael0; 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 forin 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: n, a, =-x100% n" 102) Calculate the false positive rate of traffic anomaly: a2 x100% no In the above two formulas, fiir is the total number of missed events per unit time, Pie is the total number of false positive events per unit time, and n, 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: (1) The method can make full 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 incident 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 are 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 incident 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, QZSS and any similar navigation satellite positioning system. GNSS track recorders equipped with floating cars, buses, etc., record the position information of vehicles at various points in time with a certain sampling frequency.1. (general requirementsP0.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 time, use a large number of travel time data in the same spatial-temporal subzonc, establish a probability distribution model of travel time, 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 travel time data in the current time period, and also establish the current travel time 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 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 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 4 = {(x" ) Ix E [X" -1), ,yr_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 69=2.5°. If the difference between the azimuth angle of the point A and the direction of the road section 11 is less than the threshold bo, i.e. 19.4 - <ö9, then the matching process is completed; if 19A 8,/ < is not satisfied, delete the road segment U 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 &script as follows.

Determine the straight-line equation of the road segment ij (if the road segment is a curve, it s divided into several straight-line segments): y -y, = k(x -

V P

where the slope is: The projection line equation is: Solve the projected coordinate P as: xi -y -y, = -1(x -x") ki -ky, + k1x, + x-k2 +1 h k2 After the map matching process, combined with the timestamp data of the coordinates of the positioning point, the positioning point is matched to the spatial-temporal sub-zone.

Step 14. All travel time data of the same path between the OD pairs constitutes the whole. The implementation method is to calculate the travel time of each vehicle in the spatial-temporal sub-zone E,: where t1,2...tn-1,n is the distance of the 1st and 2nd GNSS positioning point in sub-zone..., the distance between the n-I th and the nth GNSS positioning point, t I...tn is the 1st..., nth in the spatial-temporal sub-zone, ..., the time stamp of the nth GNSS positioning point; the data in each spatial-temporal sub-zone is not filtered and constitutes a set re, 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 + +(cc,, -kx, models, the best model is determined by Bayesian infontion Criterion (BIC). Calculate: f ('3 =Lit) n e (1,2,.... 51 totally 5 mixed model. At the same time calculate the BIC of the five models which is defined as: MC = -2 In L + k In n 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 ii. g, 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 fiunction of the current traffic condition. The method is the same as Step 15. Record the parameter vectors II, Step 17. Calculate the difference between the two travel time distributions according to the description parameters iv_ an of the current traffic characteristics and the description parameters ii, ft. a of the historical traffic characteristics: cliff [(Tin,r+." arD It, = JSD(P, HP).

Step 18. Normalize the difference in travel time distribution of each spatial-temporal subzonc to a normalized value ay of 0-1: a mqx. (diff)-min (diff) Calculate the traffic incident index of each spatial-temporal subzone Ar-ce,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 x200m 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 4 ={(x" yc) I e [x" e [ye, ; 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 mid the azimuth angle: diff; rntin (cliff) Search for the road segment closest to point A, and take (50=2.5°, if the difference between the azimuth angle of the point A and the direction of the road section is less than the threshold 60, i.e. l94 -eul< 59, then the matching process is completed; if < is not satisfied, delete the road segment 1/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 /f (if the road segment is a curve, it is divided into several straight-line segments): y-j = k(x -where the slope is: k = -3', The projection line equation is: X1 -=--(x-xj

P

Solve the projected coordinate P as: P = k2 After the map matching process, combined with the timestamp data of the coordinates of the positioning point, the positioning point is matched to the spatial-temporal sub-zone.

Step 24. Calculate calculate the travel time of each vehicle in the spatial-temporal sub-zone: 14, = tfl*t2,3 where 112... is the travel time between the first and second GNSS anchor points in the spatial-temporal sub-zone, ..., the travel time between the n-lth and nth GNSS anchor points; specify the time segment length fp, the upper limit of the number of segment data at the same time pmax; search for the travel time data in the time segment of the i-1 th time segment in a time-space sub-zone, if the number of travel time data in the time segment exceeds the upper limit pmat, randomly take The strip data constitutes a set T4: and is used for subsequent processing.

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: 1)(i 5) =Diu u(t n E 0,2 5} =1 totally 5 mixed model. At the same time, calculate the BIC of the five models, which is defined as: BIC = -21n 1. + Jr * In n 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 ii. 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 Mx) of the travel time 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: pi( ..,)= I * . f, (t,, Calculate the Jensen-Shannon divergence between the two distributions di,: Jd = SD(P 0) = 1 -D(P An+ 1 -2D(Qi Al) 2 where P and Q are two different probability distributions. M = + Q), I) is the Kullback-Le bier divergence: P(x) D H = IP(;)log:) k-1 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: ow( /II g)= itiog f(r1) >D(.111g) n g(xf) The divergence between the two distributions is expressed as a distance matrix: D = The matrix satisfies ciy=c4; and cL-0(1=1).

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 (F-fl between the traffic environment Hand 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 t el time in the spatial-temporal sub-zone, which constitutes the overall real-time travel time Establish a travel time probability ity distribution model p" (1") = IT * f,(1. ,uf,CT i) and carry 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-:-1).

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 travel time distributions according to the description parameters ikt, pre, an of the current traffic characteristics and the description parameters ii. g, a of the historical traffic characteristics: rity[(fft, ff,,, cc),(niff, a)] = J-521)(P, 11 P) Step 28. Normalize the difference in travel time distribution of each spatio-temporal subzone to a normalized value aj of 0-1: cuff -min (cliff) 5: max -min (cliff) Calculate the traffic incident index of each spatial-temporal subzone 2t=CfrX 10.

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 = v3)1 x, E [XT Xt E Denote the floating vehicle GNSS data acquisition frequency as_fi)=1/to, denote the point adjacent to A in time as P(tz-to), then PQA-lto) can be defined as the 1-adjacent point of A. and P(/.4-2t0), P(t4+2t0) can be defined as the 2-adjacent point of A, and so on, then P(14-1(10), P(t.,+kto) can be defined as the k-adjacent point of A. Whenfo< IHz, k I or 2. Take the road segment //with the smallest distance from the k-adjacent point and A, and calculate 19.45 the mean value of the direction angle of the k-adjacent point of A and A. if -old < c5;, is met, the matching process is completed. If not, search in other road segments until the condition 9A' A' - So is 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 t].. (if the road segment is a curve, it is divided into several straight-line segments): y-= k(x -where the slope is: The projection line equation is: Y A = -1( iYc -x.4) Solve the projected coordinate P as: -Icy, +k2 x1 x k1 +1 1(2Y YE k1+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 spatial-temporal sub-zone.

Step 34. Calculate the travel time of each vehicle in the spatial-temporal sub-zone: t,2 = ...+ where 62... is the travel time between the first and second GNSS anchor points in the spatial-temporal sub-zone. ..., the travel time between the n-lth and nth GNSS anchor points; specify,' the time segment length t, the upper limit of the number of segment data at the same time p","; search for the travel time data in the time segment of the i-1 th time segment in a time-space sub-zone, if the number of travel time data in the time segment exceeds the upper limit p" randomly take The p"", strip data constitutes a set Tc, and is used for subsequent processing.

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, ..., K components are parameterized separately; for the K models, the best model is determined by Bayesian Information Criterion (BIC). Calculate: 1,(t)=E17,/,,( ) 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 = -21n 1, + k * Inn where L is the maximum likelihood function value, k is the number of parameters in the model, and P7 is the total amount of data.

After that, the minimum hybrid model of BIC is selected and its parameter vectors th g, 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 functionMx) of the travel time 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: = Eq.; * (lj "PP 'Cr 7_7 Calculate the Jensen-Shannon divergence between the two distributions = JSD(P 0) = -1 D(P It/1)+ -2,0(Q Al) -where P and Q are two different probability distributions, Al = -1(P + 0) , D is the Kullback-Leibler divergence: P(;)

K

D(P H = P(; )log Qock) A =1 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: DKr (f g) = II M. 11g) The divergence between the two distributions is expressed as a distance matrix: D = The matrix satisfies clu-d,, and do=0(i=1).

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. -7) 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 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 time in the spatial-temporal sub-zone, which constitutes the overall real-time travel time Establish a trip travel time probability distribution model p(1/4) = (ti; ,c) j) and carry out J=1 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--1).

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 travel time distributions according to the description parameters rin, ELI, an of the current traffic characteristics and the description parameters'', Et, a of the historical traffic characteristics: rieff r(rt" p.", an), (4 = J,STD(P, P)* Step 38. Normalize the difference in travel time distribution of each spatial-temporal subzone to a normalized value air of 0-1: -min (dill) a. -max (thin -min (dW.) Calculate the traffic incident index of each spatial-temporal subzone A,r-ct,rx 10.

Claims (11)

1. CLAIMS1. A method for detecting traffic anomaly based on travel time distribution with equidistant epatial-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) Pre-processing historical trajectory data: processing GNSS positioning historical data of a floating vehicle into sampled vehicle travel time data of historical trajectory; Preprocessing real-time trajectory data: processing GNSS positioning real-time data of a floating vehicle into sampled vehicle travel time data of real-time trajectory; 3) Historical trajectory data analysis: using said sampled vehicle travel time data of historical trajectory, establishing a historical travel time probability distribution, and then obtaining a historical traffic feature model Ph; Real-time trajectory data analysis: using said sampled vehicle travel time data of real-time trajectory, establishing a real-time travel time probability distribution, and then obtaining a real-time traffic feature model Prt; 4) Difference calculation: Jensen-Shannon divergence is used to measure difference between said historical travel time probability distribution and said real-time travel time 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 step I) adopts the following division method: Equidistant spatial-temporal division method based on fixed segment size of time segment and spatial segment: take 30min as one time segment; take 200mx200m as one spatial segment.
2. 2. A method according to claim I, wherein the preprocessing of the historical trajectory data in step 2) comprises: 2a) Data structumlizing: GNSS positioning data is processed by data cleaning, data integration, data conversion and data reduction, in order to get stnicturalized GNSS positioning historical data; 2b) Map matching: combined with urban road network data, project GNSS positioning points to 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) Travel time calculation and sampling of historical trajectory: select a plurality of OD pairs in each spatial-temporal subzone, calculate a travel time between the OD pairs according to the structuralized GNSS positioning historical data, obtain travel time data of historical trajectory, and perform data sampling on said travel time data of historical trajectory to obtain sampled travel time data of historical trajectory.
3. 3 A method according to claim 2, wherein said map matching in step 2b) comprises the following sub-steps: 2b1) Dividing 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.= {(x, , ) I E [-r" Y, E [yr, IA, 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 Rs can be denoted as ij, and assign an index to each road section; 262) Searching for road section ü where a positioning point A is located in said set of road sections Rs according to its distance and azimuth angle; 2b3) GNSS positioning point linear projection method is used to calculate projection coordinates of the GNSS positioning point on the road segment.
4. 4. A method according to claim 3, wherein the 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 094 and direction angle of said road section if bf) is less than a threshold (69), 164 <80' then said map matching is completed; if condition 1634 -9"1< is not met, search in other road sections until said condition is met; 2b22) Point sequence matching scheme, which is applicable for high frequency floating vehicle data" denote floating vehicle GNSS data acquisition frequency as fu I Ito; one point adjacent to A in time as P(t4-to) or P(ti+to) is defined as 1-point-adjacent to A, and P(ti-2t0) or P(14-12t0) is defined as 2point-adjacent to A, and so on, then P(ti-kto) or P(ii-lkto) is defined as k-point-adjacent to A; when fi<lHz, k=1 or 2, take said road section ij with the smallest distance from said k-point-adjacent to A to said point A, and calculate 9, which is mean value of direction angles of said k-point-adjacent to A and said point A; if condition 1194 -<89 is met, then said map matching is completed; if not, search in other road sections until said condition 1M, - < t' is met.
5. 5. A method according to claim 2, wherein the travel time calculation and sampling of the historical trajectory in step 2c) adopts one of the following methods: 2c1) Full sample scheme of travel time: calculate travel time of each vehicle in a spatial-temporal sub-zone t; = t2.3 where 112...6,_," is the distance between the lst and 2" GNSS positioning point in subzone..., the distance between the n-Ith and the nth GNSS positioning point, is the l' 11th time stamp of GNSS positioning point in said space-time subzone, respectively: travel time data in each spatial-temporal subzone is not filtered out and constitutes a set T; for subsequent processing; 2c2) Time-smooth sampling scheme: specify' time segment size, set up an upper limit: maximum number of data pieces in said time segment; search for travel data in each time segment in a spatial-temporal subzone; if the number of travel time data in said time segment exceeds said upper limit, maximum number of travel time data is randomly selected for subsequent processing; then calculate travel time of each vehicle in said special-temporal: 15 =11,2 +12,7, +... ± , where ti.. is the travel time between lr and 2111 GNSS positioning point in subzone..., the travel time between the n_lthl and the nth GNSS positioning point, 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.; then search for travel time data in the ith time segment of said special-temporal subzone, if the number of travel time data in the time segment exceeds said upper limit p a, then p ma r pieces of travel time data are randomly selected and constitutes a set T4: for subsequent processing.
6. 6. A method according to any of claims 1 to 5, wherein one of the following methods is used to establish said historical travel time probability distribution in step 3): 3a) Historical trajectory data classification method: According to temperature, precipitation, visibility and traffic control measures, historical data under no traffic anomalies are divided into different categories; 3b) Historical data clustering method: For historical data, differences between any two spatial-temporal subzoncs are used to obtain quantitative characterization of different spatial-temporal subzoncs, and the quantitative characterization are used for clustering; using temperature, precipitation, visibility and traffic control measures as traffic environment characteristic factors, carrying out multiple Logit regressions to establish a mapping relationship between traffic environment characteristic factors and different categories.
7. 7. A method according to claim 6, wherein the historical data clustering method according to step 3b) comprises: 3b1) Using EM algorithm to estimate the parameters of the historical traffic feature model Ph, and to obtain parameters ii. g, a, wherein i is the proportional vector of each sub-component in the historical traffic feature model, g is the mean vector of each sub-component in the historical traffic feature model, a is a standard deviation vector of each sub-component in the historical traffic feature model: 3b2) According to the parameters estimated, establishing a probability density function pi(x) of the historical travel time probability distributions corresponding to spatial-temporal subzones of different dates; 3b3) Calculating Jensen-Shannon divergence du between any two of said historical travel time probability distributions: = CO= -2D(P 11 Al) + -2D(Q11 Ai) Al 1 = -(P + 0.) Wherein P and Q represent two different probability distributions; 2; .1S7) is Jensen-Shannon divergence, and D is Kullback-Leibler divergence: DO' H = 1'(x)log PGY0A-k=1In case of finite mixing modeling, the value of D cannot be explicitly expressed, but Monte Carlo sampling method can be used to approximate the calculation: I kre ( Hg) 1 = -Dog >I)(f 11g) n,=, g(x) wherein, Rig: represents Kullback-Leibler divergence approximated by the Monte Carlo sampling method, while f and g represent any two historical travel time probability distributions; 3b4) Combining all the divergence dij to form a distance matrix: d" D= the distance matrix satisfies d"-d", d"-0 for i=j; 3b5) using the distance matrix as an input of K-Medoids algorithm, obtaining clustering results and indexing the different categories; 3b6) taking index of said different categories category as dependent variable, and taking traffic environment characteristic factors as independent variable, performing multiple logit regression to obtain mapping relationship between traffic environment E and traffic category T indicated as 3b7) Aggregating historical data from the same category, re-establishing a new historical travel time probability distribution, and obtaining said historical traffic feature model Ph.
8. 8. A method according to any of claims 1 to 5, the real-time trajectory data analysis in step 3) includes the following steps: 3b3) Calculating travel time Lyi between any two OD pairs in the spatial-temporal subzone, forming a real-time travel time set T;-,n; 3b4) Establishing a trip time probability distribution; obtaining a real-time traffic feature model Eq, * ft (t; wherein K represents the number of sub-components of real-time A-1 traffic feature model, i represents the proportion of a certain sub-component of the real-time traffic feature model, ,tt represents the mean of a certain sub-component of the real-time traffic feature model, and a represents the standard deviation of a sub-component of the real-time traffic feature model; implement parameter estimation to obtain description parameters Tin, pin, art, wherein, in is the proportional vector of each sub-component in the real-time traffic feature model, ix, is the mean vector of each sub-component in the real-time traffic feature model, an is the standard deviation (2(xk) vector of each sub-component in the real-time traffic feature model; 3b5) taking index of said different categories category as dependent variable, and taking traffic environment characteristic factors as independent variable, performing multiple logit regression to obtain mapping relationship between traffic environment E and traffic category T indicated as 9.
9. A method according to any of claims I to 5, wherein traffic anomaly detection in step 4) comprises: 4a) if a traffic category T is identified, then locate historical trajectory data under category T; 4b) calculate the difference between historical and real-time traffic features: dijr[01,11",t3"),(n,11.03)1= ./.57)(p, 11 Po wherein irt is the proportional vector of each sub-component in the real-time traffic feature model, Etn is the mean vector of each sub-component in the real-time traffic feature model, ari is the standard deviation vector of each sub-component in the real-time traffic feature model; ri is the proportional vector of each sub-component in the history traffic feature model, 1.1 is the mean vector of each sub-component in the historical traffic feature model, and o is the standard deviation vector of each sub-component in the historical traffic feature model; JSD is Jensen-Shannon divergence.ID.
10. A method according to claims 1 to 5, wherein step 5) comprises: 5a) Normalizing differences among travel time distributions of all spatial-temporal subzoncs to a normalized value 05 of 0-1: diff" - (diall[max (cliff)-min (dig' )1 wherein, dillrefers to difference between travel time distributions of any two spatial-temporal subzones; 5b) Calculating a traffic incident index for each spatial-temporal subzone A ci-acfix10.
11. 11. A method according to claims 1 to 5, wherein step 6) comprises: 6a) Calculating a false negative rate of traffic anomaly events, a1=(1") x 100% ,where in nr represents the total number of false negative events, and a, represents the total number of traffic anomaly events that should be reported; 6b) Calculating a false positive rate of the traffic incident state, a2 = "26) 1 00%, wherein n, is the total number of false positive events and fl" is the total number of traffic anomaly events that should be reported.