CN109716414B - Multi-mode road traffic abnormity detection method - Google Patents

Abstract

A multi-mode road traffic abnormity detection method is characterized in that traffic states are represented by probability distribution of travel speeds, traffic state differences are reflected by probability distribution difference measurement indexes, different modes are classified and processed by introducing environmental states and traffic management variables, spatial position information of floating vehicles at different moments can be obtained by utilizing vehicle-mounted GNSS positioning devices of the floating vehicles, and intelligent detection of urban road traffic abnormity events can be achieved by analyzing and mining massive floating vehicle tracks. The method has the characteristics of clear principle, simple and convenient implementation and high detection rate.

Description

Multi-mode road traffic abnormity detection method

Technical Field

The invention belongs to the technical field of traffic detection. In particular, the invention relates to a real-time detection method for urban road traffic abnormity. The method comprises the steps that spatial position information of a floating vehicle at different moments can be obtained through a vehicle-mounted GNSS positioning device of the floating vehicle, and travel vehicle speed probability distribution in a specific space-time range is obtained through data preprocessing, map matching and data fusion; according to the change condition of the speed distribution, the urban road traffic abnormal event can be effectively identified.

Background

The detection of traffic abnormal events is an important component of urban traffic management and is also one of the core functions of an intelligent traffic system. Traffic anomalies primarily include traffic accidents, vehicle breakdown, truck droppings, damage or malfunction to road traffic facilities, and other special events that cause traffic flow disturbances. The events easily cause traffic jam and reduce road section traffic capacity, and influence the normal operation of the whole road traffic system in serious cases. Through the detection of the traffic abnormal events, traffic managers can know the traffic abnormal information in time and take proper induction and control measures to reduce the adverse effect of the traffic abnormal events.

The traffic abnormal event detection can be divided into a manual mode and an automatic mode. The manual mode comprises patrol cars, emergency call reporting, video monitoring and the like, and the requirements of traffic management cannot be met due to the fact that manpower and material resources are consumed and the real-time performance is poor. The automatic mode is realized by means of an automatic event Detection (AID) algorithm, and the basic principle is to identify traffic abnormal events by detecting the change of road traffic flow at different positions. The AID algorithms commonly used at present include pattern recognition algorithms (such as California algorithm and monica algorithm), statistical prediction algorithms (such as exponential smoothing method and kalman filtering algorithm), traffic flow model algorithms (such as McMaster algorithm) and intelligent recognition algorithms (such as artificial neural network and fuzzy logic algorithm).

However, the current detection method has the defects of high facility requirement, high calculation complexity, incapability of further judging the situation of the abnormal condition and the like. The invention utilizes the track data returned by the GNSS positioning devices on the taxi and the bus to establish a historical traffic state database and a real-time traffic state database, and identifies the traffic abnormal event by analyzing the traffic flow characteristic difference reflected by the historical traffic state database and the real-time traffic state database. The method has the characteristics of good real-time performance, capability of parallel processing, high recognition rate, low requirement on detection facilities and the like, and is suitable for detecting the urban road traffic abnormal events in the data environment with real-time floating car positioning data.

Currently, for monitoring abnormal traffic events, there are the following representative technologies:

one united states patent application, US 20160148512, discloses the principle of construction and implementation of a traffic anomaly detection and reporting system. The system consists of a sensor, a communication module, a mobile processing module and a user interaction module. The sensor is used for collecting relevant data of the periphery of the vehicle; the communication module is used for sending the data of the vehicle and receiving the data of the surrounding vehicles; the mobile processing module is used for processing and analyzing data of related vehicles in a certain area and generating a traffic incident report; the user interaction module can provide traffic event reports like a user. The scheme is a traffic abnormal event detection technology based on a vehicle-vehicle and vehicle-road communication network, and abnormal events can be judged by utilizing various information collected by a sensor. However, the sensor and the communication unit need to be separately installed and debugged, so that the implementation difficulty is high; mobile processing unit processing capacity is limited; meanwhile, a mobile and fixed message receiving end is needed, the system has failure probability, and the reliability is poor.

A chinese patent application, CN 104809878A, discloses a method for detecting abnormal traffic conditions on urban roads by using GPS data of buses. According to the scheme, a road section delay time index is obtained according to GPS historical data, an instantaneous speed, a periodic average speed, a weighted moving average speed and a multi-vehicle average speed are obtained according to GPS current data, and an abnormal condition is detected by using a standard variable analysis algorithm. The scheme does not need to add a detection facility, and is convenient to implement. However, the representation of the traffic situation is too simplified, and the characteristics and causes of the abnormal traffic condition cannot be analyzed; the method has the advantages that the method lacks basis for dividing the traffic scene, and cannot consider the influence of factors such as weather on the change of the traffic situation.

Disclosure of Invention

In order to more clearly illustrate the contents of the present invention, first, the terms of art involved are explained as follows:

floating car: also known as probe cars. The vehicle-mounted positioning device is mounted on the bus and the taxi which run on the urban road.

GNSS: global Navigation Satellite System (Global Navigation Satellite System). Including GPS, GLONASS, GALILEO, and Beidou satellite navigation systems.

Space-time subregion: the fragment area is divided according to two dimensions of time and space, and reflects the situation in a certain space range within a period of time. Divide a day into several time segments, e.g. 0: 00-0: 10,0: 10-0: 20 … …, each time slice being referred to as a time sub-section; dividing an implementation area for detecting the abnormal urban road traffic into a plurality of space segments, such as areas with the longitude of 121.58 degrees E-121.59 degrees E and the latitude of 31.16 degrees N-31.17 degrees N, wherein each space segment is called a space subregion; the spatio-temporal segments formed by the intersection of any one time sub-region and any one space sub-region, called spatio-temporal sub-regions, have a longitude of 121.58 ° E-121.59 ° E, and a latitude of 31.16 ° N-31.17 ° N in the region 0: 00-0: 10, in space-time segments.

Historical track data: the historical trajectory data is trajectory data that is accumulated over a long period of time and stored in a database. The historical track data is dynamically changed data and needs to be updated in time and reprocessed and analyzed periodically to ensure the accuracy of extracting the historical traffic characteristics. The data for each spatio-temporal subregion can be processed in parallel to improve efficiency. Which may be referred to herein simply as historical data.

Real-time trajectory data: the real-time trajectory data is a set of trajectory data in a time segment closest to the current time. Which may be referred to herein simply as real-time data.

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

And (3) traffic abnormity: traffic flow disorder caused by traffic accidents, vehicle breakdown, truck dropping, damage or failure of road traffic facilities and other events.

Severity of traffic abnormalities: that is, the severity of the traffic flow disorder is the difference between the normal traffic flow and the traffic flow characteristics after the occurrence of the traffic abnormality.

Traffic anomaly index: a measure of the severity of a traffic abnormality. The range is 0-10, and the larger the numerical value is, the more serious the traffic abnormality is.

The traffic environment is as follows: the sum of all external influences and forces acting on road traffic participants. Including road conditions, traffic facilities, terrain and topography, weather conditions, and other traffic activities of traffic participants.

Map matching: a process of associating geographic coordinates with a city road network.

Peak hour flow: the maximum value of hour traffic flow in a day of a certain urban road section.

Finite mixture model: a mathematical method for modeling complex densities with simple densities. The finite mixture model with the variable set as y and the number of components as K can be expressed as:

response variable: variables that change in accordance with an independent variable are also referred to as dependent variables.

Bayesian information criterion: under incomplete information, the method estimates the partially unknown state by subjective probability and then corrects the occurrence probability by a Bayes formula to obtain an evaluation index of the reliability of the result. The calculation method comprises the following steps:

BIC＝-2ln L+k·ln n

in the formula, L is the maximum value of the likelihood function, k is the number of unknown parameters, and n is the sample size.

Likelihood function: the likelihood function is a function of the parameters of the statistical model. The probability that the likelihood function L (θ | X) for a parameter θ is (numerically) equal to the variable X after the given parameter θ, given an output X: l (θ | X) ═ P (X ═ X | θ).

Parameter estimation: a method of estimating an unknown parameter contained in a population distribution from samples extracted from the population.

The EM algorithm: the Expectation Maximization Algorithm (Expectation Maximization Algorithm) is an iterative Algorithm used for maximum likelihood estimation or maximum a posteriori probability estimation of probability parameter models with hidden variables.

Kullback-Leibler divergence: a measure of the difference between the two probability distributions P and Q.

Jensen-Shannon divergence: is a symmetric form of the Kullback-Leibler divergence.

K-Medoids algorithm: a clustering algorithm selects, for each iteration, a point from the current class whose sum of distances to all other points (in the current class) is the smallest as a new center point.

The invention aims to establish a scheme for identifying road traffic abnormal events by combining historical GNSS positioning data and real-time GNSS positioning data and traffic environment information based on a floating car track recording system. In order to achieve the purpose, the invention provides the following technical scheme:

the implementation premise of the invention is as follows: floating cars (taxis, buses, etc.) carrying GNSS trajectory recorders; the data center has large-scale storage, calculation and real-time task processing capacity.

The application range of the invention is as follows: urban roads (including ground roads and elevated roads) through which the floating vehicles pass.

The implementation steps of the invention comprise:

1) determining the space-time range of detection and establishing space-time subareas.

And determining the time range and the space range in which the traffic abnormal event detection is required based on the actual application requirement. The time range may be set to all days, i.e. 0: 00-24: 00; it may also be set to a certain period of time, for example to detect 17: 00-20: 00, the detection time range is set to 17: 00-20: 00, to name just one specific example, there are many other scenarios that will not be described here. The spatial range can be set to a certain market domain according to administrative divisions, such as Beijing, Shanghai, Huangpu district, etc.; or a certain city functional area, such as a certain city central business area, an industrial area, etc., may be set according to the city space structure.

The establishment of the space-time subarea refers to dividing the detection time range into a plurality of smaller time segments and dividing the detection space range, namely the implementation area of the urban road traffic abnormality detection, into a plurality of smaller space segments. The space-time sub-regions can be established by adopting various empirical division methods, including an equidistant space-time division method and a non-equidistant space-time division method.

2) And (4) preprocessing data.

The GNSS positioning data is subjected to data cleaning, data integration, data conversion and data reduction, and the structuralization degree of the data is improved. GNSS, a global navigation satellite system positioning system, is a space-based radio navigation positioning system that can provide all-weather three-dimensional coordinates and velocity and time information at any location on the earth's surface or in near-earth space. The Satellite Positioning system mainly comprises four global Navigation Positioning systems, namely a GPS (global Positioning system) in the United states, a GLONASS (global Navigation Satellite system) in Russia, a GALILEO in European Union and a Beidou Satellite Navigation system in China, and also comprises regional Navigation Positioning systems such as a QZSS in Japan, an IRNSS in India and the like, and Satellite Positioning enhancement systems such as a WASS in the United states, an MSAS in Japan and the like. In order to establish a unified data distribution standard among different navigation positioning system devices, the american national Marine Electronics association has established a unified nema (national Marine Electronics association) communication protocol to standardize the data broadcasting of GNSS. Therefore, the member systems in GNSS, such as GPS, GLONASS, etc., have consistent data distribution formats although established and maintained by different countries and organizations, respectively, and thus do not need to be transformed.

In the selected space range, there are many vehicles equipped with GNSS positioning devices, and there are common taxis, buses, freight cars, private cars, and the like. Based on the current situation of urban traffic data application, in practical application, an urban taxi is usually selected as a floating car as a data source of a traffic abnormality detection system.

The acquired GNSS positioning information contains some unreasonable information, and in order to ensure the accuracy of the detection and judgment result of the abnormal traffic state, the GNSS positioning information needs to be screened firstly to provide abnormal data, so that the reliability of the data is ensured. These anomaly data include: data falling outside the detection spatio-temporal range, spatial position jumps significantly beyond a reasonable range. The "spatial position jump significantly out of the reasonable range" is exemplified below. If 10: 30: the locating point uploaded by the locating equipment of a certain floating car at the time of 00 is recorded as A, and the current day is 10: 30: and at 30 moment, the positioning point uploaded by the positioning equipment of the floating car is marked as B, the distance between the position A and the position B is 1500 meters, and the running speed of the floating car is calculated to be at least 180km/h according to the calculation, which exceeds the common knowledge, so that the floating car is an abnormal spatial position jump and should be removed in data processing.

3) And (4) fast map matching.

The preprocessed GNSS positioning data need to be combined with urban road network data, a GNSS positioning point is projected to an urban map through a map matching algorithm, a matching relation between the positioning point and a road section is established, and errors caused by positioning drift are corrected.

At present, electronic maps of various geographic areas are detailed, and the electronic maps can be sourced from geographic information systems of cities, and can also be sourced from other modes and paths. The electronic maps describe the urban road information in detail, and a plurality of road sections can be obtained through division. By means of information such as distance and angle, positioning points are matched to the road sections, and therefore the positioning information is matched to the actual geographic environment.

4) A representation of a floating vehicle path and matching of different vehicle paths.

The path of the vehicle may not be unique given a set of starting and ending points. A complex urban traffic network comprises several links, the different links being numbered, for example, as L1, L2, etc. A road may have two different directions of travel, in which case the two different directions of travel should be represented as two different road segments, given different road segment numbers.

Given a starting point and an end point, intersections of road segments in the urban road network may generally be used. When a path on which a certain floating vehicle travels is known, it is now necessary to select the same path as that of the floating vehicle from the path information already transmitted from other floating vehicles, thereby obtaining the same path group between the starting point and the ending point.

5) And (6) sampling data.

The positioning data of the floating car includes information such as position coordinates, instantaneous speed, recording time, and the like. In the method for detecting the urban road traffic abnormality based on the floating car data, the data sampling refers to screening partial data from all floating car data for subsequent analysis and processing, and the screening is performed based on the computing capacity of a data center and the accuracy requirement proposed in advance. Different data sampling methods may be employed based on different computational power and accuracy requirements. For example, when the computing power of the data center is strong and the requirement on the detection precision is high, all floating car positioning data can be used as processing objects to perform comprehensive processing analysis; when the computing power of the data center is limited, assuming that the current data center can process 500 pieces of data in 1 minute for each space subregion, and the actual situation is that 2000 pieces of floating car positioning data can be generated in 1 minute for each space subregion, 500 pieces of data can be randomly extracted from 2000 pieces of data to be analyzed, so that the processing result with limited accuracy can be obtained within the computing power range of the data center.

According to different utilization modes of the floating car data, different attributes of the floating car data can be sampled, such as travel speed, travel time and the like. The urban road traffic abnormity detection method based on floating car data provided by the patent adopts the travel speed as the basis to detect the urban road traffic abnormity. Thus, data sampling refers to sampling the trip vehicle speed.

6) Analyzing historical track data and extracting characteristics.

The historical track data refers to track data of the floating cars accumulated in long-term urban road traffic operation. By utilizing historical floating vehicle track data, an urban road traffic characteristic model can be established to reflect the general characteristics of urban traffic operation. The urban road traffic characteristic model can refer to some specific indexes, such as average speed, weighted average speed and the like; but also to some kind of statistical model, such as the probability distribution of the speed of the journey. Many conventional models use a single index to represent traffic characteristics (such as historical average vehicle speed) of a certain road section or area, but this method is simple and convenient to apply, but has low accuracy and low sensitivity, and often cannot achieve good effect in detecting abnormal traffic conditions. Therefore, the patent proposes that for each space-time subregion, traffic characteristics are described by using probability distribution of traffic characteristic variables, a traffic characteristic model is established, and parameter estimation is carried out.

The traffic characteristic variables which can be collected comprise travel speed, travel time and the like, and the probability distribution of the traffic characteristic variables in the patent refers to the probability distribution of the travel speed.

7) And analyzing real-time track data and extracting characteristics.

The real-time track data refers to track data of a floating car in traffic operation within a period of time short of the current moment. By utilizing the real-time floating car track data, the change dynamics of the traffic characteristics can be mastered so as to reflect the instant characteristics of the current traffic operation. The patent describes current traffic characteristics by using the travel speed of a current time-space subregion.

8) And (4) detecting the abnormality.

The idea of system state anomaly detection is firstly proposed by Denning, that is, by monitoring the abnormal conditions of system use on the audit record of the system, the events which violate safety and possibly cause system anomaly can be detected. The model established by Denning is independent of any specific system, application environment, system weakness and fault type, and is an anomaly detection model in a general sense. The model comprises 5 parts of a subject, an object, an audit record, a contour, an exception record and an activity rule. A contour is the normal behavior of a subject relative to an object, represented by a metric and statistical model. The Denning model defines 3 metrics, namely event counters, interval timers, resource measurers, and proposes 5 statistical models, namely, operational models, mean and standard deviation models, multivariate models, markov process models, and time series models. The model proposed by Denning establishes a normal behavior feature profile of a system main body based on statistics by analyzing audit data of the system, and during detection, the audit data in the system is compared with the established normal behavior feature profile of the main body, and if the different parts exceed a certain threshold, an abnormal event is considered. The model lays a foundation for anomaly detection, and a plurality of anomaly detection methods and systems developed in the future are developed on the basis of the anomaly detection method and the system.

In the development process of the anomaly detection technology in recent years, more artificial intelligence methods are introduced to improve the performance of anomaly detection. The artificial intelligence methods mainly comprise data mining, artificial neural networks, fuzzy evidence theory and the like. The method of data mining is used to determine what features are most important in a large collection of data. The technology is mainly used for anomaly detection, and a more concise definition of a normal mode is sought, rather than simply listing all normal modes as in the conventional anomaly detection method. The introduction of the data mining method enables the detection system to generally include normal patterns not included in the training data by merely identifying key features in the normal patterns. In the statistical anomaly detection referred to above, user behavior data is divided into two categories according to some statistical criteria: namely abnormal behavior and normal behavior. Because the statistical-based method has certain difficulty in extracting and abstracting the audit instance, large errors can be caused, the method must depend on some probability distribution assumptions, and the measurement of the user behavior is generally described by experience and feeling, so that the clustering method of the artificial neural network is introduced. The artificial neural network has self-learning self-adaptive capacity, the neural network is trained by using sample points representing normal user behaviors, the neural network can extract normal user or system activity patterns from data through repeated learning for many times and codes the patterns into a network structure, and when the artificial neural network is detected, the audited data passes through the learned neural network to judge whether the system is normal. Because the judgment standard of the abnormity has certain ambiguity, the fuzzy evidence theory is introduced into the abnormity, for example, an intrusion detection framework model based on a fuzzy expert system is established, and the false alarm rate can be reduced better.

The patent provides an anomaly detection scheme based on statistical characteristics, and the basic idea is to measure the difference between historical traffic characteristics and real-time traffic characteristics through Jensen-Shannon divergence to realize the detection of abnormal traffic conditions. The scheme has the advantages of good interpretability and low calculation burden, overcomes the defects of inaccurate and untimely detection by adopting single statistic, and avoids the defects of heavy calculation burden and high hardware requirement of methods such as an artificial neural network and the like.

9) And (4) quantifying and representing the severity of the abnormality and issuing abnormal information.

The severity of the abnormal traffic condition should be released to the public in a concise and clear manner to avoid possible areas of congestion and improve the operating efficiency of urban traffic. The severity of the abnormal condition is characterized by a traffic anomaly index ranging from 0 to 10, where 0 indicates no anomaly and 10 indicates a high degree of anomaly.

And projecting the abnormal occurrence position on an electronic map, and publishing the abnormal occurrence position in forms of intelligent mobile equipment APP and the like.

10) And (5) evaluating the system performance.

The evaluation of the system performance refers to the evaluation of the accuracy of the detection of the abnormal traffic state, and the evaluation indexes comprise a false alarm rate and a missing report rate. The lower the false alarm rate and the false negative rate, the better the performance of the system.

In the step 1), the space-time sub-regions may be divided specifically by the following method:

11) equidistant space-time division method. Determining the segment scale of a time dimension, wherein the span of the time segment is a fixed value, and usually 30min is taken as a time segment; determining the segment scale of the space dimension, wherein the span of the space segment is a fixed value, and a space grid of 200m × 200m is usually taken as a space segment;

12) non-equidistant space-time division method based on road network density: based on the density of road network as the judgment index, when the density of road network is more than or equal to 2km/km²Taking a time segment of 30min and a space segment of 200m multiplied by 200 m; when the density of the road network is less than 2km/km²Taking a time segment of 30min and a space segment of 400m multiplied by 400 m;

13) non-equidistant space-time division method based on peak hourly flow: taking a time segment of 30min and a space segment of 200m multiplied by 200m when the peak hour flow is more than or equal to 1000 vehicles/hour based on the peak hour flow as a judgment index; when the peak hour flow is less than 1000 vehicles/hour, a time segment of 30min and a space segment of 400m × 400m are taken.

The step 3) specifically comprises the following steps:

31) the spatial region to be processed is divided into grids of a certain size, and the range of each grid region can be represented as A_s＝{(x_s，y_s)|x_s∈[x_r，x_r+1)，y_s∈[y_r，y_r+1) -each grid area contains a number of road segments, the set of these road segments being denoted R_sSet R of said road sections_sEach road section in (1) is represented as ij, and each road section is assigned with a number;

32) judging the grid area where the positioning point is positioned, and utilizing the distance and the azimuth angle to assemble R in the road section_sSearching a road section ij where a certain positioning point A is located, wherein the matching scheme comprises the following steps:

321) the single-point matching scheme comprises the following steps:

searching a road section closest to the point A, and when the difference value between the driving direction angle meeting the point A and the direction angle of the road section ij is smaller than a threshold value, namely, meeting the condition of theta_A-θ_ij|＜δ_θCompleting the matching, said threshold value delta_θ2.5 degrees, 5 degrees, 10 degrees and the like can be selected; if not satisfy | theta_A-θ_ij|＜δ_θAnd deleting the road sections ij in the search space, and continuing to search other road sections until the condition is met. The matching method is shown in fig. 3.

322) Point sequence matching scheme:

this scheme is applicable to high frequency floating car data. Representation of floating car GNSS data acquisition frequency as f₀＝1/t₀A point P (t) adjacent to A in time_A-t₀)，P(t_A+t₀) 1-neighbor, P (t), defined as A_A-2t₀)，P(t_A+2t₀) Defined as 2-neighbors of A, and so on, then P (t)_A-kt₀)，P(t_A+kt₀) Defined as k-neighbors of A. At f₀If the frequency is less than 1Hz, k is 1 or 2. Taking the road section ij with the minimum distance between the distance A and the k-adjacent point of the distance A, and calculating the average value of the driving direction angles of the k-adjacent points of the distance A and the distance A

If it satisfies

Completing matching; otherwise, searching other road sections until meeting

33) And calculating the projection coordinate of the GNSS positioning point on the road section by using a linear equation of the road section (approximately splitting into a straight line if the road section is a curve road section), thereby reducing the error caused by GNSS positioning drift. The specific method adopts a GNSS positioning point linear projection method as follows:

determining a linear equation of the road section ij (if the road section is a curve, dividing the road section into a plurality of linear sections):

y-y_i＝k(x-x_i)

wherein the slope is:

the projection line equation is:

solving the projection coordinate P as:

after the map matching process, matching the positioning points to the space-time sub-area by combining the timestamp data of the coordinates of the positioning points.

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

51) the overall scheme of speed information. The whole travel speed data of each sub-floating car in a space-time sub-area form a whole. The implementation method is that the travel speed of each vehicle in a space-time subregion xi is calculated:

wherein d is_1，2...d_n-1，nIs the distance between the 1 st and the 2 nd GNSS positioning points in the space-time sub-area xi, the distance between the n-1 st and the n th GNSS positioning point, and t₁...t_nThe timestamp of the nth GNSS positioning point is the 1 st within the space-time sub-region xi. Data in each space-time subregion are not screened to form a set V_ξFor subsequent processing.

52) A time-smoothed sampling scheme of velocity information. Appointing the length of the time segment, and setting the upper limit of the number of data of the same time segment; searching the speed data in each time segment in one time-space subregion, and randomly taking the data of the upper limit number for subsequent places if the number of the speed data in the time segment exceeds the upper limitAnd (6) processing. The implementation method is that the travel speed of each vehicle in a space-time subregion xi is calculated:

wherein d is_1，2...d_n-1，nIs the distance between the 1 st and the 2 nd GNSS positioning points in the space-time sub-area xi, the distance between the n-1 st and the n th GNSS positioning point, and t₁...t_nThe timestamp of the nth GNSS positioning point is the 1 st within the space-time sub-region xi. Specifying a time slice length t_pUpper limit of number of data pieces p of same time segment_max(ii) a Searching the speed data in the ith time segment of a time in a space-time subregion, if the number of the speed data in the time segment exceeds the upper limit p_maxRandom access of p_maxStrip data Add V_ξAnd used for subsequent processing.

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

61) and (3) a simple historical track data fusion method. And taking the historical data under the condition of no traffic abnormality as a whole to establish a traffic characteristic model and estimate parameters. The method utilizes a finite mixture model to establish a traffic characteristic model and carries out parameter estimation. Specifically, one of the following three schemes can be adopted:

611) gaussian mixture model with fixed composition

The scheme adopts a mixed Gaussian model with fixed component quantity K to describe the probability distribution of the vehicle speed. The component amounts are 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 number of components K cannot be too small. Generally, K is 4 to 6.

612) Gaussian mixture model with variable component number

The scheme adopts a model evaluation-based method to select proper component quantity, and the method comprises the following steps:

determining the possible maximum component quantity K, and respectively carrying out parameter estimation on a mixed Gaussian model with n being 1, 2, … and K components; for the K models, the best model is determined by Bayesian Information Criterion (BIC). The number of largest components is generally chosen according to the accuracy requirement, but it must be noted that the larger the number of components, the slower the convergence of the maximization algorithm is expected to be. The maximum ingredient amount selected here is K5, i.e. it needs to be calculated:

there are 5 mixed models. Meanwhile, BIC was calculated for 5 models, which is defined as:

BIC＝-2ln L+k·ln n

in the formula, L is the maximum likelihood function value, k is the number of parameters in the model, and n is the total amount of data.

And then, selecting a mixed model with the minimum BIC, and recording parameter vectors eta, mu and sigma of the mixed model, wherein eta is a proportion vector occupied by each subcomponent in the historical traffic characteristic model, mu is a mean vector of each subcomponent in the historical traffic characteristic model, and sigma is a standard deviation vector of each subcomponent in the historical traffic characteristic model, and the mixture model is used as a characteristic record of the space-time subregion. The density profile of the mixed model is shown in fig. 6.

613) Finite mixed model with variable component quantity and distribution type

The scheme adopts the same model-based evaluation method as 612), but the distribution form of the subcomponents and the quantity of the subcomponents are variable, and the method is as follows:

selecting M probability distribution models as the distribution types of the subcomponents, including but not limited to: normal distribution, gamma distribution, weibull distribution. When a normal distribution is used, the sub-distribution function takes:

when using a gamma distribution, the sub-distribution function takes:

wherein

When using a weibull distribution, the sub-distribution function takes:

assuming that the distribution types of all the subcomponents of the mixture model are the same, the possible maximum number of components K is determined. And selecting the distribution type of the M seed components and the quantity of the K components to jointly form MK combinations, respectively calculating the BIC value, and taking the model with the minimum BIC as the optimal model.

62) And (4) classifying historical track data by context. According to the air temperature, the precipitation, the visibility and the traffic control measures, the historical data under the condition of no traffic abnormality is divided into different categories, and a model is respectively established and parameter estimation is carried out. The implementation method comprises the following steps:

according to the difference of air temperature, precipitation, visibility and traffic control measures, traffic environments are divided into 5-8 categories, and historical data are classified into the categories according to the difference of the traffic environments corresponding to the historical data. For each category, the processing as described in 5) is performed, so that a mapping relation R (E → T) is established, where E is the traffic environment and T is the traffic situation.

63) Historical data clustering. And for historical data, obtaining difference quantitative representations of different space regions through comparison between every two space-time sub-regions, and clustering by using the quantized differences. And performing multiple Logit regressions by taking the air temperature, the precipitation, the visibility and the traffic control measures as characteristic factors, and establishing a mapping relation between the traffic environment and the category. The implementation scheme is shown in figure 4. The implementation steps are as follows:

631) according to the method of 5), a traffic characteristic model is established, and parameter estimation is carried out.

632) According to the previous finite mixture model parameter estimation result, a probability density function p of the travel vehicle speed distribution corresponding to different dates of the space-time subarea is represented_i(x) The parameters are given by a Gaussian mixture model as an example:

in the formula, K represents the number of subcomponents of the stroke vehicle speed distribution, η represents the proportion of a certain subcomponent in the stroke vehicle speed distribution, μ represents the average value of a certain subcomponent in the stroke vehicle speed distribution, and σ represents the standard deviation of a certain subcomponent in the stroke vehicle speed distribution.

633) Calculating Jensen-Shannon divergence d between every two distributions_ij：

Where P, Q are two different probability distributions,

d is Kullback-Leibler divergence:

in the case of the finite mixture model, the value cannot be expressed explicitly, but can be approximated by a monte carlo sampling method, which is:

in the formula, D_MCThe Kullback-Leibler divergence obtained by approximate calculation by adopting a monte carlo sampling method is shown, and f and g represent any two distribution functions.

634) The divergence between each two distributions is represented as a distance matrix:

the matrix satisfies d_ij＝d_ji，d_ij＝0(i＝j)。

635) And taking the distance matrix as the input of a K-Medoids algorithm to obtain a clustering result, and establishing indexes for categories.

636) And (4) performing multiple Logit regression by taking the category index as a response variable and taking the traffic environment data (including air temperature, precipitation, visibility and the like) as independent variables to obtain a mapping relation R (E → T) between the traffic environment E and the traffic situation category T.

637) And aggregating the data of the same category, reestablishing a mixed model by utilizing the new aggregated data set, and performing parameter estimation to obtain a final historical traffic characteristic data set.

The step 7) may specifically adopt the following method:

71) simple real-time data processing. The method is carried out simultaneously with 61). And carrying out model establishment and parameter estimation on the real-time traffic data to obtain a characteristic function of the current traffic condition. The implementation steps of the method are the same as 61), except that the adopted data are real-time traffic data.

72) And (4) a classification treatment method. The method is carried out simultaneously with 62) or 63). And acquiring a characteristic function of the traffic condition, acquiring information such as current air temperature, precipitation, visibility, traffic control measures and the like, and judging the category of the current traffic condition. The implementation scheme is shown in figure 5. The implementation steps are as follows:

721) calculating travel speed in the space-time subarea to form a real-time travel speed overall V_ξ，rt；

722) Establishing travel vehicle speed probability distribution model

And estimating parameters;

723) and taking the current traffic environment data (including air temperature, precipitation, visibility and the like) as input parameters, and obtaining the category T of the current traffic situation by utilizing a mapping relation R (E → T).

The step 8) specifically comprises the following steps:

81) when step 72) is adopted, positioning historical traffic characteristic data under the category according to the category T to which the current traffic situation belongs, otherwise, not processing;

82) description parameter eta according to current traffic characteristics_rt、μ_rt、σ_rtAnd the description parameters eta, mu and sigma of the historical traffic characteristics calculate the difference between the two speed distributions: diff [ (eta) ]_rt，μ_rt，σ_rt)，(η，μ，σ)]＝JSD(P_rtP). Wherein eta is_rtIs the proportion vector, mu, of each sub-component in the real-time traffic characteristic model_rtIs the mean vector, sigma, of each subcomponent in the real-time traffic feature model_rtThe standard deviation vector of each sub-component in the real-time traffic characteristic model is obtained; eta is a proportion vector occupied by each subcomponent in the historical traffic characteristic model, mu is a mean vector of each subcomponent in the historical traffic characteristic model, and sigma is a standard deviation vector of each subcomponent in the historical traffic characteristic model. When the historical traffic characteristics are similar to the real-time traffic characteristics (namely the historical travel vehicle speed distribution and the real-time travel vehicle speed distribution), a smaller Jensen-Shannon divergence value is obtained, namely the difference between the Jensen-Shannon divergence value and the Shannon divergence value is smaller; when the historical traffic characteristics and the real-time traffic characteristics are greatly different, a larger Jensen-Shannon divergence value is obtained, that is, the difference between the two is larger, that is, the probability of abnormality is larger, see fig. 7.

The step 9) specifically comprises the following steps:

91) standardizing the speed distribution difference of each space-time subregion to a normalized value a of 0-1_ξ：

92) Calculating the traffic abnormality index A of each space-time subregion_ξ＝a_ξ×10；

93) The regional position with the abnormal index higher than 5 is projected on an electronic map, and the intelligent mobile device APP and other forms are published to the society in a public mode, so that drivers can avoid potential congestion points, and the traffic efficiency of urban road traffic is improved.

The step 10) specifically comprises the following steps:

101) calculating the rate of missing report of the abnormal traffic state:

102) calculating the false alarm rate of the abnormal traffic state:

in the above two formulae, n_lIs the total number of missed events per unit time, n_eIs the total number of false alarm events in unit time, n_aIs the total number of times that an anomaly actually occurs within a unit time.

Compared with similar technologies in the same field, the invention has the following advantages:

(1) the method fully utilizes the existing floating car operation data (GNSS track data), detects the change of the traffic state through historical traffic feature extraction and real-time traffic situation analysis, and can realize real-time, low-cost and intelligent detection of the urban road traffic abnormal event;

(2) the probability distribution of the traffic characteristic parameters is used as the description of the traffic characteristics, the reflected characteristics are more comprehensive, the one-sidedness and instability of the traffic characteristics represented by a single index are avoided, and the detection reliability is higher;

(3) aiming at the characteristic that the traffic characteristics are influenced by the traffic environment (such as weather conditions), a clustering-multinomial Logit regression joint algorithm is introduced, and the mapping relation between the traffic environment characteristics and the traffic situation categories is established;

(4) through the inspection of actual data, the urban road traffic anomaly detection technology based on floating car data provided by the invention can realize the detection of an abnormal event with higher accuracy, the detection rate exceeds 90%, the missing report rate is lower than 15%, the false report rate is lower than 20%, a good detection effect is obtained, and the technology can be applied to the intelligent management and service of urban traffic.

Drawings

The details and advantages of the present invention will become apparent and readily appreciated when considered in conjunction with the following drawings, in which:

FIG. 1 shows a schematic diagram of the components and basic principles of the present invention;

FIG. 2 shows a general flow diagram of the present invention in its implementation;

FIG. 3 is a diagram illustrating a fast map matching algorithm embodiment of the present invention;

FIG. 4 is a flow diagram illustrating a historical traffic feature extraction scheme implemented in accordance with the present invention;

FIG. 5 is a schematic flow diagram illustrating a real-time traffic feature extraction scheme implemented in accordance with the present invention;

FIG. 6 shows a morphological schematic of a Gaussian mixture model probability distribution;

FIG. 7 illustrates a graphical representation of the difference in comparison of historical traffic characteristics to real-time traffic characteristics.

Detailed description of the preferred embodiments

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

As shown in fig. 1, the overall system architecture of the present invention comprises: the vehicle-mounted GNSS track recorder is carried by the floating vehicle, and comprises a data center, a GNSS satellite and a communication system. The GNSS here includes any similar navigation satellite positioning system such as GPS, GLONASS, GALILEO, Beidou, IRNSS, QZSS, etc. The GNSS track recorder carried by the floating vehicles such as taxies, buses and the like records the position information of the vehicles at each time point in the driving process at a certain sampling frequency f (generally, f is required to be more than 0.1Hz), and transmits the position information to the data center in real time through a GPRS mobile communication network (wireless network communication technologies such as WCDMA, TD-LTE and the like can also be adopted, but the cost is correspondingly improved). The data center establishes a historical road traffic characteristic database through data preprocessing and data fusion and a specific algorithm; establishing a real-time traffic characteristic database for recently received real-time data; and judging whether the current traffic characteristics are abnormal or not through the mapping relation between the historical database and the real-time database, and performing visual display through the processing terminal to generate a traffic abnormal event report.

The general flow of the scheme is shown in fig. 2, and comprises the steps of collecting and storing GNSS track data, establishing a space-time subregion, extracting historical traffic characteristics, extracting real-time traffic characteristics, identifying anomalies and the like. The acquisition and storage of GNSS trajectory data is a data basis of the whole scheme, and due to the huge data magnitude, a distributed storage scheme is adopted, and the content of the invention is not the content of the existing mature technology for distributed storage. The method is characterized in that space-time sub-areas are established, the basic assumption is that the same traffic characteristics exist in a specific area and a specific time period, and the assumption is universally applicable after long-term observation. The historical traffic characteristic extraction method is characterized by utilizing GNSS track data to calculate travel speed, utilizing a large amount of travel speed data of the same time-space subregion to establish a probability distribution model of the speed, carrying out parameter estimation and representing traffic characteristics by a small amount of parameters. The principle of the real-time traffic feature extraction is that speed data in the current time period is processed and analyzed, and a current vehicle speed probability distribution model is built. The difference measure index is adopted to judge the change degree of the real-time characteristic compared with the historical characteristic, and whether the traffic abnormal event occurs is determined according to whether the change degree reaches the threshold value.

According to a combination of the implementation methods described in the summary of the invention, the examples are given below.

Example one

Step 11, determining the segment scale of the time dimension by adopting an equidistant space-time division method, wherein the time segment span is a fixed value, and generally taking 30min as a time segment; the segment size of the spatial dimension is determined, the span of the spatial segment is a fixed value, and a spatial grid of 200m × 200m is taken as a spatial segment.

And step 12, preprocessing data, and performing data cleaning, data integration, data conversion and data reduction on the GNSS positioning data to improve the structuralization degree of the data.

Step 13, dividing the spatial region to be processed into grids of a certain size, where the range of each grid region can be represented as a_s＝{(x_s，y_s)|x_s∈[x_r，x_r+1)，y_s∈[y_r，y_r+1) }; judging a grid area where the positioning point is located, and searching a road section where the positioning point is located by using the distance and the azimuth angle; searching a road section nearest to the point A, and taking a threshold value delta_θWhen the difference between the driving direction angle at the point a and the direction angle at the link ij is smaller than the threshold value δ 2.5 °_θWhen, i.e. | θ is satisfied_A-θ_ij|＜δ_θCompleting matching; if not satisfy | theta_A-θ_ij|＜δ_θDeleting the road sections ij in the search space, and continuing to search other road sections until the condition is met; the method comprises the following steps of calculating the projection coordinates of the GNSS positioning points on the road section by using a linear equation of the road section (approximately splitting into a straight line if the road section is a curve road section), and reducing errors caused by GNSS positioning drift, wherein the specific method comprises the following steps:

determining a linear equation of the road section ij (if the road section is a curve, dividing the road section into a plurality of linear sections):

y-y_i＝k(x-x_i)

wherein the slope is:

the projection line equation is:

solving the projection coordinate P as:

after the map matching process, matching the positioning points to the space-time sub-area by combining the timestamp data of the coordinates of the positioning points.

And step 14, forming a whole by all travel speed data of each sub-floating vehicle in a space-time sub-area. Calculating the travel speed of each vehicle in a space-time sub-area xi:

wherein d is_1，2...d_n-1，nIs the distance between the 1 st and the 2 nd GNSS positioning points in the space-time sub-area xi, the distance between the n-1 st and the n th GNSS positioning point, and t₁...t_nThe timestamp of the nth GNSS positioning point is the 1 st within the space-time sub-region xi. Data in each space-time subregion are not screened to form a set V_ξFor subsequent processing.

And step 15, taking the historical data under the condition of no traffic abnormality as a whole to establish a traffic characteristic model and estimate parameters. The method utilizes a finite mixture model to establish a traffic characteristic model and carries out parameter estimation. Taking the maximum component number K as 5, and respectively carrying out parameter estimation on a Gaussian mixture model with n as 1, 2, … and K components; for the K models, the best model is determined by Bayesian Information Criterion (BIC). And (3) calculating:

there are 5 mixed models. Meanwhile, the BIC of 5 models was calculated:

BIC＝-2ln L+k·ln n

in the formula, L is the maximum likelihood function value, k is the number of parameters in the model, and n is the total amount of data.

And then, selecting a mixed model with the minimum BIC, and recording parameter vectors eta, mu and sigma of the mixed model as the characteristic record of the space-time sub-area.

Step 16, carrying out model establishment and parameter estimation on the real-time traffic data to obtain a characteristic function of the current traffic condition, wherein the method has the same steps as the step five, and recording a parameter vector eta_rt、μ_rt、σ_rt。

Step 17, according to the description parameter eta of the current traffic characteristics_rt、μ_rt、σ_rtAnd the description parameters eta, mu and sigma of the historical traffic characteristics calculate the difference between the two speed distributions: diff [ (eta) ]_rt，μ_rt，σ_rt)，(η，μ，σ)]＝JSD(P_rt||P)。

Step 18, standardizing the speed distribution difference of each space-time subregion to a normalized value a of 0-1_ξi：

Calculating the traffic abnormality index A of each space-time subregion_ξ＝a_ξ×10。

Example two

Step 21, determining the segment scale of the time dimension by adopting an equidistant space-time division method, wherein the time segment span is a fixed value, and generally taking 30min as a time segment; the segment size of the spatial dimension is determined, the span of the spatial segment is a fixed value, and a spatial grid of 200m × 200m is taken as a spatial segment.

And step 22, preprocessing data, and performing data cleaning, data integration, data conversion and data reduction on the GNSS positioning data to improve the structuralization degree of the data.

Step 23, dividing the spatial region to be processed into grids of a certain size, where the range of each grid region can be represented as a_s＝{(x_s，y_s)|x_s∈[x_r，x_r+1)，y_s∈[y_r，y_r+1) }; judging a grid area where the positioning point is located, and searching a road section where the positioning point is located by using the distance and the azimuth angle; searching a road section nearest to the point A, and taking a threshold value delta_θWhen the difference between the driving direction angle at the point a and the direction angle at the link ij is smaller than the threshold value δ 2.5 °_θWhen, i.e. | θ is satisfied_A-θ_ij|＜δ_θCompleting matching; if not satisfy | theta_A-θ_ij|＜δ_θDeleting the road sections ij in the search space, and continuing to search other road sections until the condition is met; calculating the projection coordinate of the GNSS positioning point on the road section by using the linear equation of the road section (approximately splitting into a straight line if the road section is a curve road section), and reducing the GNSS positioning driftThe error is specifically determined as follows:

determining a linear equation of the road section ij (if the road section is a curve, dividing the road section into a plurality of linear sections):

y-y_i＝k(x-x_i)

wherein the slope is:

the projection line equation is:

solving the projection coordinate P as:

after the map matching process, matching the positioning points to the space-time sub-area by combining the timestamp data of the coordinates of the positioning points.

Step 24, calculating the travel speed of each vehicle in the space-time sub-area xi:

wherein d is_{1， 2}...d_n-1，nIs the distance between the 1 st and the 2 nd GNSS positioning points in the space-time sub-area xi, the distance between the n-1 st and the n th GNSS positioning point, and t₁...t_nThe timestamp of the nth GNSS positioning point is the 1 st within the space-time sub-region xi. Specifying a time slice length t_pUpper limit of number of data pieces p of same time segment_max(ii) a Searching the speed data in the ith time segment of a time in a space-time subregion, if the number of the speed data in the time segment exceeds the upper limit p_maxRandom access of p_maxStrip data Add V_ξ。。

And 25, taking the historical data under the condition of no traffic abnormality as a whole to establish a traffic characteristic model and estimate parameters. The method utilizes a finite mixture model to establish a traffic characteristic model and carries out parameter estimation. Taking the maximum component number K as 5, and respectively carrying out parameter estimation on a Gaussian mixture model with n as 1, 2, … and K components; for the K models, the best model is determined by Bayesian Information Criterion (BIC). And (3) calculating:

there are 5 mixed models. Meanwhile, the BIC of 5 models was calculated:

BIC＝-2ln L+k·ln n

in the formula, L is the maximum likelihood function value, k is the number of parameters in the model, and n is the total amount of data.

And then, selecting a mixed model with the minimum BIC, and recording parameter vectors eta, mu and sigma of the mixed model as the characteristic record of the space-time sub-area.

According to the parameter estimation result, a probability density function p of the travel vehicle speed distribution of the space-time subarea corresponding to different dates is shown_i(x)：

Calculating Jensen-Shannon divergence d between every two distributions_ij：

Where P, Q are two different probability distributions,

d is Kullback-Leibler divergence:

under the condition of adopting a finite mixture model, adopting a Monte Carlo sampling method to approximate calculation, wherein the calculation method comprises the following steps:

the divergence between each two distributions is represented as a distance matrix:

the matrix satisfies d_ij＝d_ji，d_ij＝0(i＝j)。

And taking the distance matrix as the input of a K-Medoids algorithm to obtain a clustering result, and establishing indexes for categories.

And (4) performing multiple Logit regression by taking the category index as a response variable and taking the traffic environment data (including air temperature, precipitation, visibility and the like) as independent variables to obtain a mapping relation R (E → T) between the traffic environment E and the traffic situation category T.

And aggregating the data of the same category, reestablishing a mixed model by utilizing the new aggregated data set, and performing parameter estimation to obtain a final historical traffic characteristic data set.

And step 26, acquiring a characteristic function of the traffic condition, acquiring information such as current air temperature, precipitation, visibility and traffic control measures, and judging the type of the current traffic condition.

Calculating travel speed in the space-time subarea to form a real-time travel speed overall V_ξ，rt；

Establishing travel vehicle speed probability distribution model

And estimating parameters;

and taking the current traffic environment data (including air temperature, precipitation, visibility and the like) as input parameters, and obtaining the category T of the current traffic situation by utilizing a mapping relation R (E → T).

Step 27, positioning historical traffic characteristic data under the category according to the category T to which the current traffic situation belongs; description parameter eta according to current traffic characteristics_rt、μ_rt、σ_rtAnd the description parameters eta, mu and sigma of the historical traffic characteristics calculate the difference between the two speed distributions: diff [ (eta) ]_rt，μ_rt，σ_rt)，(η，μ，σ)]＝JSD(P_rt||P)。

Step 28, standardizing the speed distribution difference of each space-time subregion to a normalized value a of 0-1_ξ：

Calculating the traffic abnormality index A of each space-time subregion_ξ＝a_ξ×10。

EXAMPLE III

Step 31, adopting a non-equidistant space-time division method to obtain a road network with the density of more than 2km/km²Or in the urban central area with the peak hour flow rate more than 1000 vehicles/hour, taking a time segment of 30min and a space segment of 200m multiplied by 200m, and regarding the road network density less than 2km/km²Or suburban areas with peak hourly traffic less than 1000 vehicles/hour, and a time segment of 30min and a space segment of 400m × 400m are taken.

And step 32, preprocessing data, performing data cleaning, data integration, data conversion and data reduction on the GNSS positioning data, and improving the structuralization degree of the data.

Step 33, dividing the spatial region to be processed into meshes of a certain size, where the range of each mesh region can be represented as a_s＝{(x_s，y_s)|x_s∈[x_r，x_r+1)，y_s∈[y_r，y_r+1)}；

Representation of floating car GNSS data acquisition frequency as f₀＝1/t₀A point P (t) adjacent to A in time_A-t₀)，P(t_A+t₀) 1-neighbor, P (t), defined as A_A-2t₀)，P(t_A+2t₀) Defined as 2-neighbors of A, and so on, then P (t)_A-kt₀)，P(t_A+kt₀) Defined as k-neighbors of A. At f₀If the frequency is less than 1Hz, k is 1 or 2. Taking the road section ij with the minimum distance between the distance A and the k-adjacent point of the distance A, and calculating the average value of the driving direction angles of the k-adjacent points of the distance A and the distance A

Taking the threshold value delta_θIf it is 5 °, then

Completing matching; otherwise, searching other road sections until the condition is met.

And calculating the projection coordinate of the GNSS positioning point on the road section by using a linear equation of the road section (approximately splitting into a straight line if the road section is a curve road section), thereby reducing the error caused by GNSS positioning drift. The specific method comprises the following steps:

determining a linear equation of the road section ij (if the road section is a curve, dividing the road section into a plurality of linear sections): y-y_i＝k(x-x_i)

Wherein the slope is:

the projection line equation is:

solving the projection coordinate P as:

after the map matching process, matching the positioning points to the space-time sub-area by combining the timestamp data of the coordinates of the positioning points.

Step 34, calculating each vehicle in space-time sub-area xiTravel vehicle speed of (1):

wherein d is_{1， 2}...d_n-1，nIs the distance between the 1 st and the 2 nd GNSS positioning points in the space-time sub-area xi, the distance between the n-1 st and the n th GNSS positioning point, and t₁...t_nThe timestamp of the nth GNSS positioning point is the 1 st within the space-time sub-region xi. Specifying a time slice length t_pUpper limit of number of data pieces p of same time segment_max(ii) a Searching the speed data in the ith time segment of a time in a space-time subregion, if the number of the speed data in the time segment exceeds the upper limit p_maxRandom access of p_maxStrip data Add V_ξ。。

And step 35, taking the historical data under the condition of no traffic abnormality as a whole to establish a traffic characteristic model and estimate parameters. The method utilizes a finite mixture model to establish a traffic characteristic model and carries out parameter estimation. Taking the maximum component number K as 5, and respectively carrying out parameter estimation on a Gaussian mixture model with n as 1, 2, … and K components; for the K models, the best model is determined by Bayesian Information Criterion (BIC). And (3) calculating:

there are 5 mixed models. Meanwhile, the BIC of 5 models was calculated:

BIC＝-2ln L+k·ln n

in the formula, L is the maximum likelihood function value, k is the number of parameters in the model, and n is the total amount of data.

And then, selecting a mixed model with the minimum BIC, and recording parameter vectors eta, mu and sigma of the mixed model as the characteristic record of the space-time sub-area.

According to the parameter estimation result, a probability density function p of the travel vehicle speed distribution of the space-time subarea corresponding to different dates is shown_i(x)：

Calculating Jensen-Shannon divergence d between every two distributions_ij：

Where P, Q are two different probability distributions,

d is Kullback-Leibler divergence:

under the condition of adopting a finite mixture model, adopting a Monte Carlo sampling method to approximate calculation, wherein the calculation method comprises the following steps:

the divergence between each two distributions is represented as a distance matrix:

the matrix satisfies d_ij＝d_ji，d_ij＝0(i＝j)。

And taking the distance matrix as the input of a K-Medoids algorithm to obtain a clustering result, and establishing indexes for categories.

And (4) performing multiple Logit regression by taking the category index as a response variable and taking the traffic environment data (including air temperature, precipitation, visibility and the like) as independent variables to obtain a mapping relation R (E → T) between the traffic environment E and the traffic situation category T.

And aggregating the data of the same category, reestablishing a mixed model by utilizing the new aggregated data set, and performing parameter estimation to obtain a final historical traffic characteristic data set.

And step 36, acquiring a characteristic function of the traffic condition, acquiring information such as current air temperature, precipitation, visibility and traffic control measures, and judging the category of the current traffic condition.

Calculating travel speed in the space-time subarea to form a real-time travel speed overall V_ξ，rt；

Establishing travel vehicle speed probability distribution model

And estimating parameters;

and taking the current traffic environment data (including air temperature, precipitation, visibility and the like) as input parameters, and obtaining the category T of the current traffic situation by utilizing a mapping relation R (E → T).

Step 37, positioning historical traffic characteristic data under the category according to the category T to which the current traffic situation belongs; description parameter eta according to current traffic characteristics_rt、μ_rt、σ_rtAnd the description parameters eta, mu and sigma of the historical traffic characteristics calculate the difference between the two speed distributions: diff [ (eta) ]_rt，μ_rt，σ_rt)，(η，μ，σ)]＝JSD(P_rt||P)。

Step 38, standardizing the speed distribution difference of each space-time subregion to a normalized value a of 0-1_ξ：

Calculating the traffic abnormality index A of each space-time subregion_ξ＝a_ξ×10。

Claims (11)

1. A multi-modal road traffic anomaly detection method comprises the following steps:

1) establishing a space-time subregion: dividing one day into a plurality of time segments, wherein each time segment is called a time subarea; dividing an implementation area for urban road traffic anomaly detection into a plurality of space segments, wherein each space segment is called a space sub-area; the intersection of any one time subregion and any one space subregion is called a space-time subregion; the space-time sub-area is established by adopting an equidistant space-time division method and a non-equidistant space-time division method.

2) Preprocessing historical track data: processing the GNSS positioning historical data of the floating car into sampling speed data of a historical track;

preprocessing real-time track data: processing the GNSS positioning real-time data of the floating car into sampling speed data of a real-time track;

3) analyzing historical track data and extracting characteristics: describing traffic characteristics by using the probability distribution of traffic characteristic variables for each space-time subregion by using the sampling vehicle speed data of the historical track, establishing a traffic characteristic model and performing parameter estimation, establishing historical travel speed probability distribution, and obtaining a historical traffic characteristic model P_h；

Real-time trajectory data analysis and feature extraction: using the sampling speed data of the real-time track to perform model establishment and parameter estimation on the real-time traffic data, acquiring a characteristic function of the current traffic condition, establishing real-time travel speed probability distribution, and obtaining a real-time traffic characteristic model P_rt；

4) Abnormality detection: measuring the difference between the historical traffic characteristic model and the real-time traffic characteristic model through the Jensen-Shannon divergence, and calculating to obtain a historical and real-time traffic characteristic difference value through the Jensen-Shannon divergence calculated by the historical traffic characteristic model and the Jensen-Shannon divergence calculated by the real-time traffic characteristic model; when the historical traffic characteristics are similar to the real-time traffic characteristics, a smaller Jensen-Shannon divergence value is obtained; when the difference between the historical traffic characteristics and the real-time traffic characteristics is larger, a larger Jensen-Shannon divergence value is obtained, namely, the probability of abnormality is larger;

5) quantitative characterization of the severity of the abnormality: calculating an abnormal index of the traffic condition by using the difference value of the historical traffic characteristics and the real-time traffic characteristics; the severity of the abnormal condition is characterized by a traffic anomaly index ranging from 0 to 10, where 0 indicates no anomaly and 10 indicates a high degree of anomaly;

6) evaluating the system performance: and calculating the missing report rate of the abnormal traffic state, calculating the false report rate of the abnormal traffic state, evaluating the detection accuracy of the abnormal traffic state and measuring the running stability of the system.

2. The road traffic abnormality detection method according to claim 1, characterized in that step 1) employs one of the following methods:

1a) equidistant space-time division method: determining the segment scale of the time dimension, wherein the time segment span is a fixed value, and taking 30min as a time segment; determining the segment scale of the space dimension, wherein the span of the space segment is a fixed value, and taking a space grid of 200m multiplied by 200m as a space segment;

1b) non-equidistant space-time division method based on road network density: based on the density of road network as the judgment index, when the density of road network is more than or equal to 2km/km²Taking a time segment of 30min and a space segment of 200m multiplied by 200 m; when the density of the road network is less than 2km/km²Taking a time segment of 30min and a space segment of 400m multiplied by 400 m;

1c) non-equidistant space-time division method based on peak hourly flow: taking a time segment of 30min and a space segment of 200m multiplied by 200m when the peak hour flow is more than or equal to 1000 vehicles/hour based on the peak hour flow as a judgment index; when the peak hour flow is less than 1000 vehicles/hour, a time segment of 30min and a space segment of 400m × 400m are taken.

3. The method for detecting road traffic abnormality according to claim 1, wherein the preprocessing of the historical trajectory data of step 2) includes:

2a) data structuring: performing data cleaning, data integration, data conversion and data reduction on the floating vehicle GNSS positioning historical data to obtain structured GNSS positioning historical data;

2b) fast map matching: projecting the structured GNSS positioning historical data to an urban road network by combining urban road network data through a map matching algorithm, establishing a matching relation between a positioning point in the structured GNSS positioning historical data and a road section, obtaining a matching relation table between the positioning point in the structured GNSS positioning historical data and the road section, and correcting errors caused by positioning drift;

2c) vehicle speed calculation and sampling of historical tracks: and calculating traffic running characteristic parameters according to the structured GNSS positioning historical data to obtain vehicle speed data of the historical track, and performing data sampling on the vehicle speed data of the historical track to obtain sampled vehicle speed data of the historical track.

4. The method of detecting road traffic abnormality according to claim 3, characterized in that the rapid map matching of step 2b) includes:

2b1) the spatial region to be processed is divided into grids of a certain size, and the range of each grid region can be represented as A_s＝{(x_s,y_s)|x_s∈[x_r,x_r+1),y_s∈[y_r,y_r+1) -each grid area contains a number of road segments, the set of these road segments being denoted R_sSet R of said road sections_sEach road section in (1) is represented as ij, and each road section is assigned with a number;

2b2) judging the grid area where the positioning point is positioned, and utilizing the distance and the azimuth angle to assemble R in the road section_sSearching a road section ij where a certain positioning point A is located;

2b3) and calculating the projection coordinates of the GNSS positioning points on the road section by utilizing a GNSS positioning point linear projection method.

5. The method for detecting road traffic abnormality according to claim 4, characterized in that step 2b2) employs one of the following methods:

2b21) the single point matching method comprises the following steps: searching a road section nearest to a certain positioning point A, wherein the implementation method comprises the following steps: set R for road sections_sWhen the traveling direction angle theta of the point A is satisfied, a certain segment ij_ADirection angle theta with road section ij_ijIs less than a threshold value delta_θWhen, i.e. | θ is satisfied_A-θ_ij|<δ_θWhen the matching is finished, the matching is finished; if not satisfy | theta_A-θ_ij|<δ_θContinuing to search the set of road segments R_sUntil | θ is satisfied_A-θ_ij|<δ_θ；

2b22) The point sequence matching method comprises the following steps: the scheme is suitable for high-frequency floating car data; the GNSS data time interval of the floating car at every two adjacent times is represented as t₀The floating car GNSS data acquisition frequency is expressed as f₀＝1/t₀The time record of a certain anchor point A is denoted t_AA point P (t) temporally adjacent to the anchor point A_A-t₀)，P(t_A+t₀) 1-neighbor, P (t), defined as A_A-2t₀)，P(t_A+2t₀) Defined as the 2-neighbor of a certain anchor point A, and so on, P (t)_A-kt₀)，P(t_A+kt₀) Defining k-adjacent points of a certain positioning point A; at f₀<Taking k as 1 or 2 when 1 Hz; taking a road section ij with the minimum distance from a certain positioning point A and a k-adjacent point of the positioning point, and calculating the mean value of the driving direction angles of the positioning point and the k-adjacent point thereof

If it satisfies

Completing matching; otherwise, searching other road sections until meeting

6. The method of detecting road traffic abnormality according to claim 3, wherein the vehicle speed calculation and sampling of the history track in step 2c) is performed by one of:

2c1) the method comprises the following steps: the total travel speed data of each secondary floating vehicle in a space-time sub-area xi form a whole, and the implementation method is that the travel speed of each vehicle in the space-time sub-area xi is calculated:

wherein d is_1,2…d_n-1,nIs the distance between the 1 st and the 2 nd GNSS positioning points in the space-time sub-area xi, … …, the distance between the n-1 st and the n th GNSS positioning point, t₁…t_nTime stamp of the 1 st, … … nth GNSS positioning point in space-time sub-area xi; the travel speed data in each space-time subregion are not screened to form a set V_ξFor subsequent processing;

2c2) time-smoothed sampling method: appointing the length of the time segment, and setting the upper limit of the number of data of the same time segment; searching speed data in each time segment in a space-time sub-area, and if the number of the speed data in the time segment exceeds an upper limit, randomly taking the data of the upper limit number for subsequent processing, wherein the implementation method is to calculate the travel speed of each vehicle in a space-time sub-area xi:

wherein d is_1,2…d_n-1,nIs the distance between the 1 st and the 2 nd GNSS positioning points in the space-time sub-area xi, … …, the distance between the n-1 st and the n th GNSS positioning point, t₁…t_nTime stamp of the 1 st, … … nth GNSS positioning point in space-time sub-area xi; specifying a time slice length t_pUpper limit of number of data pieces p of same time segment_max(ii) a Searching the speed data in the ith time segment of a time in a space-time subregion, if the number of the speed data in the time segment exceeds the upper limit p_maxRandom access of p_maxStrip data Add V_ξAnd used for subsequent processing.

7. The method for detecting road traffic abnormality according to one of claims 1 to 6, characterized in that the historical trajectory data analysis and feature extraction of step 3) employs one of the following methods:

3a) situation-based historical track data classification: according to air temperature, precipitation, visibility and traffic control measures, dividing historical data under the condition of no traffic abnormality into different categories, respectively establishing a model and carrying out parameter estimation;

3b) historical data clustering method: for historical data, difference quantitative representations of different space-time regions are obtained through comparison between every two space-time sub-regions, and clustering is carried out by utilizing quantized differences; and performing multiple Logit regressions by taking the air temperature, the precipitation, the visibility and the traffic control measures as characteristic factors, and establishing a mapping relation between the traffic environment and the category.

8. The method for detecting road traffic abnormality according to claim 7, wherein the historical data clustering method of step 3b) includes:

3b1) performing parameter estimation on a historical traffic characteristic model by adopting an EM (effective noise) algorithm to obtain description parameters eta, mu and sigma of historical traffic characteristics, wherein eta is a proportion vector occupied by each subcomponent in the historical traffic characteristic model, mu is a mean vector of each subcomponent in the historical traffic characteristic model, and sigma is a standard deviation vector of each subcomponent in the historical traffic characteristic model;

3b2) according to the parameter estimation result, a probability density function p of the travel vehicle speed distribution of the space-time subarea corresponding to different dates is shown_i(x)；

3b3) Calculating Jensen-Shannon divergence d between every two distributions_ij：

Where P, Q are two different probability distributions,

JSD is Jensen-Shannon divergence, D is Kullback-Leibler divergence:

in the case of the finite mixture model, the value cannot be expressed explicitly, but can be approximated by a monte carlo sampling method, which is:

in the above formula, D_MCRepresenting the Kullback-Leibler divergence obtained by approximate calculation by adopting a Monte Carlo sampling method, wherein f and g represent any two distribution functions;

3b4) distributing the travel time by divergence d between each two_ijAnd combining to form a distance matrix:

the matrix satisfies d_ij＝d_ji，d_ij＝0(i＝j)；

3b5) Taking the distance matrix as the input of a K-Medoids algorithm to obtain a clustering result, and establishing indexes for categories;

3b6) performing multiple Logit regression by taking the category index as a response variable and taking the traffic environment data as an independent variable to obtain a mapping relation R (E → T) between the traffic environment E and the traffic situation category T;

3b7) and aggregating the data of the same category, reestablishing a mixed model by utilizing the new aggregated data set, and performing parameter estimation to obtain a final historical traffic characteristic data set.

9. The method for detecting road traffic abnormality according to claim 8, wherein the real-time trajectory data analysis and feature extraction of step 3) includes the steps of:

3c) calculating travel speed in the space-time subarea to form a real-time travel speed overall V_ξ,rt；

3d) Establishing travel vehicle speed probability distribution model

In the formula, KRepresenting the quantity of the subcomponents of the real-time traffic characteristics, eta representing the proportion of a certain subcomponent of the real-time traffic characteristics, mu representing the mean value of a certain subcomponent of the real-time traffic characteristics model, sigma representing the standard deviation of a certain subcomponent of the real-time traffic characteristics model, and performing parameter estimation to obtain the description parameter eta of the current real-time traffic characteristics_rt、μ_rt、σ_rtWherein η_rtIs the proportion vector, mu, of each sub-component in the real-time traffic characteristic model_rtIs the mean vector, σ, of each subcomponent in the real-time traffic feature model_rtThe standard deviation vectors of all the subcomponents in the real-time traffic characteristic model are obtained;

3e) and taking the current traffic environment data (including air temperature, precipitation, visibility and the like) as input parameters, and obtaining the category T of the current traffic situation by utilizing a mapping relation R (E → T), wherein E represents the traffic environment data.

10. The road traffic abnormality detection method according to claim 8, characterized in that the abnormality detection of step 4) includes:

4a) according to the category T of the current traffic situation, positioning historical traffic characteristic data under the category, and if no category division exists, the category does not need to be distinguished;

4b) description parameter eta according to current traffic characteristics_rt、μ_rt、σ_rtAnd the description parameters eta, mu and sigma of the historical traffic characteristics calculate the difference between the two speed distributions: diff [ (eta) ]_rt,μ_rt,σ_rt),(η,μ,σ)]＝JSD(P_rt||P_h) Wherein η_rtIs the proportion vector, mu, of each sub-component in the real-time traffic characteristic model_rtIs the mean vector, σ, of each subcomponent in the real-time traffic feature model_rtThe standard deviation vectors of all the subcomponents in the real-time traffic characteristic model are obtained; eta is a proportion vector occupied by each subcomponent in the historical traffic characteristic model, mu is a mean vector of each subcomponent in the historical traffic characteristic model, and sigma is a standard deviation of each subcomponent in the historical traffic characteristic modelAnd (5) vector quantity.

11. The road traffic abnormality detection method according to one of claims 1 to 6, characterized in that the step 5) quantitative characterization of abnormality severity comprises:

5a) dividing a space-time subregion xi_iIs normalized to a normalized value a of 0 to 1_ξi：

5b) Calculating the traffic abnormality index A of each space-time subregion_ξi＝a_ξi×10。

Images