1
Scalable Pythagorean Mean based Incident Detection in
Smart Transportation Systems
MD. JAMINUR ISLAM∗ , Western Michigan University, USA
JOSE PAOLO TALUSAN∗ , Vanderbilt University, USA
SHAMEEK BHATTACHARJEE‡ , Western Michigan University, USA
FRANCIS TIAUSAS, Nara Institute of Sc. and Tech., JapanABHISHEK DUBEY, Vanderbilt University, USA
KEIICHI YASUMOTO, Nara Institute of Sc. and Tech., JapanSAJAL K. DAS, Missouri University of Science and Technology, USA
Modern smart cities need smart transportation solutions to quickly detect various traffic emergencies andincidents in the city to avoid cascading traffic disruptions. To materialize this, roadside units and ambienttransportation sensors are being deployed to collect speed data that enables the monitoring of traffic conditions
on each road segment. In this paper, we first propose a scalable data-driven anomaly-based traffic incident
detection framework for a city-scale smart transportation system. Specifically, we propose an incremental
region growing approximation algorithm for optimal Spatio-temporal clustering of road segments and their
data; such that road segments are strategically divided into highly correlated clusters. The highly correlated
clusters enable identifying a Pythagorean Mean-based invariant as an anomaly detection metric that is highly
stable under no incidents but shows a deviation in the presence of incidents. We learn the bounds of theinvariants in a robust manner such that anomaly detection can generalize to unseen events, even whenlearning from real noisy data. Second, using cluster-level detection, we propose a folded Gaussian classifier
to pinpoint the particular segment in a cluster where the incident happened in an automated manner. We
perform extensive experimental validation using mobility data collected from four cities in Tennessee, comparewith the state-of-the-art ML methods, to prove that our method can detect incidents within each cluster inreal-time and outperforms known ML methods.
CCS Concepts: • Computing methodologies → Learning latent representations; • Applied computing→ Transportation.
Additional Key Words and Phrases: Weakly Unsupervised Learning, Anomaly Detection, Smart Transportation,
Graph Algorithms, Cluster Analysis, Regression, Incident Detection, Approximation Algorithm.
∗ Both authors contributed equally to this research.
‡ Corresponding Author

Authors’ addresses: Md. Jaminur Islam, Western Michigan University, Kalamazoo, MI, USA, mohammadjaminur.islam@wmich.edu; Jose Paolo Talusan, Vanderbilt University, Nashville, TN, USA, jose.paolo.talusan@vanderbilt.edu; Shameek
Bhattacharjee‡ , Western Michigan University, Kalamazoo, MI, USA, shameek.bhattacharjee@wmich.edu; Francis Tiausas,
Nara Institute of Sc. and Tech., Nara, Japan, tiausas.francis_jerome.ta5@is.naist.jp; Abhishek Dubey, Vanderbilt University,
Nashville, TN, USA, abhishek.dubey@vanderbilt.edu; Keiichi Yasumoto, Nara Institute of Sc. and Tech., Nara, Japan,
yasumoto@is.naist.jp; Sajal K. Das, Missouri University of Science and Technology, Rolla, MO, USA, sdas@mst.edu.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee
provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and thefull citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored.
Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires
prior specific permission and/or a fee. Request permissions from permissions@acm.org.
© 2023 Copyright held by the owner/author(s). Publication rights licensed to ACM.
2378-962X/2023/1-ART1 $15.00
https://doi.org/10.1145/3603381
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:2

Md. J. Islam and J.P. Talusan, et al.

ACM Reference Format:
Md. Jaminur Islam, Jose Paolo Talusan, Shameek Bhattacharjee‡ , Francis Tiausas, Abhishek Dubey, Keiichi
Yasumoto, and Sajal K. Das. 2023. Scalable Pythagorean Mean based Incident Detection in Smart Transportation
Systems. ACM Trans. Cyber-Phys. Syst. 1, 1, Article 1 (January 2023), 25 pages. https://doi.org/10.1145/3603381

1

INTRODUCTION

Rapid urbanization has proliferated the number of vehicles in cities leading to increasing congestion
and a higher number of traffic accidents. For any traffic accident, delayed detection and response
from first responders or emergency management agencies can worsen into heavy city-wide congestion and even in the loss of life. This delay is one of the most important challenges faced bycommunities across the globe [17].
To monitor the transportation infrastructure, three approaches have emerged to increase thevisibility of real-time road conditions: (i) vehicular crowdsourcing, (ii) video-based anomaly detection; and (iii) sensor-based data collection. Vehicular crowdsourcing involves cities leveraging
commercial crowdsourcing platforms (such as Waze), to gather content reported by citizen users onthese platforms to get real-time observations on traffic events. However, traffic incident detection
is often unreliable and strong verification of the human-reported data cannot be guaranteed inreal-time. Video anomaly detection[5] leverages cameras and sensors deployed by the city to detect
traffic emergencies. This approach requires expensive edge devices, longer model training times,
and continuous maintenance. Many environmental and connectivity constraints also negatively
influence video quality and real-time availability. The computational resources needed to monitor
and identify traffic incidents are high and not community scalable.
To avoid the above problems in these two paradigms, smart cities are deploying traffic sensors
and Road Side Units (RSU) along roads and highways that collect traffic data from speed sensors or
smart cars [11]. The RSU infrastructure is a typical IoT network that is decentralized, low-powered,and resource-constrained in nature. However, given the ubiquity and number of devices, the RSU
infrastructure can be utilized to work together in a distributed capacity, to design intelligent
lightweight anomaly-based traffic incident detection in real-time that would otherwise be too
computationally intensive, geographically impossible, or costly.
Challenges: We view traffic incidents as anomalies that occur between otherwise normal traffic
patterns. However, characterizing a normal traffic pattern that works at a large city scale is notstraightforward due to (i) day-to-day variability of traffic, (ii) local neighborhood dependencies,(iii) a large number of speed sensors and road segments. Hence, the nature of the problem falls
under smart living CPS, which, unlike industrial CPS, is not just bound by tightly defined laws ofphysics. Therefore, the anomaly detection problem is much more challenging and requires novel
advances compared to existing theories of anomaly detection in CPS.
Furthermore, previous works on smart metering [1] have attempted to solve the anomaly detection challenge in smart living CPS. However, such efforts used data collected from small experimental
testbeds. Thus, the scale of the problem was smaller, and training data was free from noise. In
contrast, our transportation CPS setting includes data collected from the wild, across a whole city.
This needs to be accounted for in the design. Specifically, geospatial factors need to be blended
with causal factors of the underlying structure of the data that characterizes benign situations.
While many prior works exist in this area, the effort in this paper takes the challenge for thewhole city with a dataset analyzed over one year to account for all seasonal and human behavioral
effects. The validation and the performance metrics reported are very robust compared to existing
works in [10, 20, 21].
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:3

Paper Contributions: We propose an unsupervised time series-based anomaly detection framework for large-scale smart transportation networks that detects traffic incidents in real time while
maintaining a low false alarm rate. The framework automatically pinpoints the area of the incident.
Specifically, we first show theoretical parallelism between the transportation problem and an
existing anomaly detection metric (Harmonic to Arithmetic Mean ratios) previously developed foranomaly detection in smart energy systems. Second, we propose a region-growing approximation
algorithm that allows the strategic partition of smart transportation CPS into clusters where thedata is highly positively correlated. The strategic partitioning guarantees 1) high invariance of
the anomaly detection metric and 2) decentralized cluster-wise implementation of our detection
framework which enables the framework to pinpoint the area of the incident. Third, we propose adata cleaning and augmentation technique to enable learning the underlying structure of benign
conditions from the data collected from the wild to reduce false alarms. Fourth, we give a techniqueto learn the bounds of the anomaly detection metric in each of the strategic partitions undernormal traffic conditions to establish the anomaly detection criterion. Finally, we validate ourapproach through extensive large-scale experiments on real mobility datasets from four cities inTennessee. Results show that our model is able to detect traffic incidents in a cluster, in real-time.
We extend our work in [12] by identifying the segment within the cluster, where the incident
originated. Pinpointing the incident at this level lets responsible agencies make decisions faster.
The performance is measured by comparing our framework’s decisions with a separate groundtruth dataset containing actual incidents recorded by the Nashville Fire and Safety Department.
The rest of the paper is organized as follows: Section 2 discusses the related work. Section 3
introduces the transportation system model. Section 4 discusses the proposed framework. Section 5extends the framework for segment-level detection. Experimental results are discussed in Section 6
followed by conclusions.
2

PREVIOUS WORK

Existing research on automatic incident detection for cyber-physical transportation systems broadly
falls into two classes. They can be classified into model-based and data-driven approaches. Modelbased approaches [8, 14] include probabilistic models [24], fuzzy C-means clustering [23], andstate-based methods which used Kalman Filtering [13] to describe the state of monitored traffic sothat usual traffic behaviors can be learned and unusual incidents can be detected. However, thesemethods require realistic assumptions for the target area and assume that their forecasting models
are representative of true uncertainty in the data. Thus, requiring extensive time-series validation.
Data-driven approaches, on the other hand, include classification methods which typically
include nearest neighbors [18], neural networks, and support vector machines. These methods
require labeled data to train and introduce new challenges regarding user data privacy. Techniquessuch as convolutional-LSTM models and neural networks [25], consume a lot of time comparing
real-time data with historical data and have high computational costs limiting their effectiveness indecentralized deployments.
3 SYSTEM MODEL DESCRIPTION
The smart transportation CPS monitors the physical world of road conditions via TMC sensors thatare deployed in each road segment. In our setting, there is one TMC sensor per road segment, so
the number of TMC sensors equals the number of road segments. The data collected from TMC is
used for various operational decisions that can control the appropriate volume of traffic to reducethe disturbance in mobility and travel times.
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:4

Md. J. Islam and J.P. Talusan, et al.

The TMC sensors are small computational units with minimal memory. Hence, each captured
information is sent to a Road Side unit [22] (RSU). Each RSU receives data from multiple TMCs andhas a larger computational power and memory. The RSUs usually have a wired back-haul link to anedge or cloud server, where data from all TMCs of an area of interest is accumulated. Dependingon implementation variations, the RSU itself could also serve as a decentralized edge server foredge analytics. However, the system-level implementation is out of the scope of our paper. We
provide a framework that can run on the edge or fog, based on the computational and networking
capabilities available to the smart city.
For this paper, a traffic incident is an anomalous event such as vehicular accidents, crime, or
man-made disaster affecting traffic flow, fire, non-recurring high duration congestion to which thepolice, emergency, and fire safety required a response. The ground truth information on incidents
was collected from Nashville Fire and Safety Department. This ground truth information contains
the location, timestamp, and date of each incident responded by the City of Nashville in 2019.
Our goal in this paper is to develop a framework and learn the parameters that automaticallydetect congestion in real-time in the test/deployment stage. The ground truth information during the testing set is used to measure the incident detection accuracy of our anomaly detection
framework. The ground truth information during the training phase is used to cross-reference
for data augmentation and cleaning that enables efficient learning of the underlying structure ofdata corresponding to benign conditions in the transportation CPS. Each TMC at the end of atime window 𝑡 sends the following information to the RSU: timestamp, road segment ID, mean
speed over the 𝑡-th time window). The TMC sensor is located at the center of each road segment.
Therefore, the distance between two road segments is the distance between the midpoint of any
two road segments. The TMCs capture ambient speeds as vehicles pass by over a particular road.
4

PROPOSED FRAMEWORK

Fig. 1. Information architecture of the framework and the control flow of interaction between components.
Aside from the region growing clustering, all components run in parallel. The cluster detection framework is
generated in the detection training phase and is used on real-world data upon deployment.

First, we provide a high-level overview of the framework, its architecture can be seen in Fig. 1,
followed by a summary of the notations used in this paper in Table 1. There are five logical modules
in which the contribution is divided:
Theoretical Intuition: We discuss the choice of harmonic mean to arithmetic mean ratio
metric [1] as an anomaly detection metric, its relevance to the problem, and its advantages andmodifications necessary to fit the transportation application.
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:5

Region Growing Approximation: For the metric to achieve invariance, we need spatial andtemporal partitions of the high dimensional data at which the positive correlation within eachpartition is maximized, which is achieved through a region growing approximation algorithm.
Invariant Design: Involves metric derivation after the region growing approximation.
Pre-processing and Augmentation: Due to the characteristics of real-world traffic data, theinvariant contains the effects of accidents. This poses a practical problem for unsupervised learning
problems such as anomaly detection. Therefore, our framework invokes a data cleaning andsanitization technique to augment synthetic benign samples of the invariant.
Learning normal operating range of invariant: Once the cleaning has been done, we obtaina low dimensional invariant that is a suitable candidate for pattern recognition of this invariant
that remains stable when there are no incidents.
Anomaly Detection Criterion: We identify the best hyperparameter inputs to the training
algorithm that gives the best output.
Table 1. List of symbols.
SymbolTotal clusters in the target area𝑐𝑘

𝑘 𝑡ℎ cluster within the set of clusters 𝐶𝑆

Set of segments in the target area𝑛Number of segments

𝑆𝑐 𝑘

Set of segments located in cluster 𝑐𝑘

𝑠𝑐𝑙 𝑘𝑙 𝑡ℎ segment in the 𝑘 𝑡ℎ cluster𝑝

Speed correlation𝑝 (𝑚𝑖𝑛)Correlation threshold

𝑝𝑐𝑢𝑡

Cut off correlation value𝑡Time slot, based on temporal granularity𝑚

Number of time slots

𝑑𝑠𝑙 (𝑡 )Mean speed 𝑑 at segment 𝑙 within 𝑆𝑐𝑘 at time 𝑡

𝐻 𝑀𝑐𝑘 (𝑡 )Harmonic Mean of cluster 𝑐𝑘 at time index 𝑡

𝐴𝑀𝑐𝑘 (𝑡 )Arithmetic Mean of cluster 𝑐𝑘 at time index 𝑡

𝑄𝑐𝑘 (𝑡 )Q-ratio metric of cluster 𝑐𝑘 at time index 𝑡(𝑡 ), Γ𝑐𝑙𝑜𝑤(𝑡 )𝑘Upper and lower Safe margins for the ratio of cluster 𝑐𝑘 at time 𝑡

𝑐𝑘

Γ𝑐𝑘

ℎ𝑖𝑔ℎ

𝜏𝑐𝑚𝑎𝑥
(ℎ), 𝜏𝑐𝑚𝑖𝑛(ℎ)
𝑘𝑘∇𝑐𝑘 (𝑡 )4.1

Description

𝐶

Upper and lower standard limits of cluster 𝑐𝑘 over historical data ℎResiduals for the ratios of cluster 𝑐𝑘 , a non-zero residual indicates a possible anomaly𝜅

Scalar Factor Hyperparameter𝑆𝐹

Sliding Frame Size Hyperparameter

Theoretical Intuition

For a large-scale CPS application such as smart transportation, the anomaly detection metric shouldhave the following properties:(1) Invariance Under Benign Conditions: Under no incidents, the metric should show minimalchange across time and across history. This is important to reduce false alarms given the low baserate of incident occurrence.
(2) Deviation Under Incidents: Under incidents, the metric should have properties that cause quickand discernible deviation in the metric. This is important to increase detection accuracy.
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:6

Md. J. Islam and J.P. Talusan, et al.

As a starting point, we leverage a recent result from [1] that showed that a collection of positively
correlated random variables sampled repeatedly over time can be represented as a time series ofratio between the harmonic to the arithmetic mean of the aggregate data; and can be used as ananomaly detection metric. This is because this the metric is stationary in its time series as longas a positive co-variance structure can be preserved. Any unforeseen data falsification attack that
disturbs the space-time covariance structure will cause deviations in the otherwise stationary
time series of Harmonic Means to Arithmetic Means. In the following, we explain the theoreticalexplanation of why the HM to AM ratio is a good starting point for our problem and examine whatnovel theoretical and applied contributions are necessary to make it work for incident detection
for a transportation CPS.

55

0.9

50

0.8

45
40

Qratio per cluster

1.0

QratioSpeed (kph1.1

speeds per segment

60

0.7
0.6
0.5
0.4

35
03:00 06:00 09:00 12:00 15:00 18:00 21:00
Time

0.3

03:00 06:00 09:00 12:00 15:00 18:00 21:00
Time

(a)

(b)

Fig. 2. Comparison of (a) Raw speed data and (b) Harmonic to the arithmetic mean ratio within a cluster.

4.1.1 Invariance Under Benign Conditions. We explain why the harmonic to the arithmetic
mean ratio is a candidate for an anomaly detection metric that is invariant under benign conditions.
The basic premise is that humans react with some shared driving behavior based on the time of theday, traffic level on the road, road type (highway or city lanes), and road width, etc. Such shared
driving behavior in the absence of incidents causes driving speeds to increase or decrease together,or remain similar that in turn manifests itself as having high positive correlation among data points.
One of the achievements of [1] is that it proved that the upper bound on the absolute difference
between the arithmetic mean and harmonic mean of the data collected from a positively correlated
system depends on two things: (1) minimum possible value of the data (denoted by 𝑑𝑚𝑖𝑛 ) and (2)the average difference in data observed between any two arbitrary sensing end-points averaged
over an appropriate time granularity (say, 𝑇 ), (denoted by 𝜉 (𝑇 )).
As long as it can be guaranteed that 𝑑𝑚𝑖𝑛 and 𝜉 (𝑇 ) do not change with time, the invariance inthe harmonic to the arithmetic mean ratio is guaranteed. We identified, however, that 𝜉 (𝑇 ) does
not change only under strategic spatial and temporal partitions which is nontrivial to achieve for atransportation CPS. We show in Fig. 2a that relying only on raw data from the sensors (harmonicmean) can lead to noisy data. However, the corresponding harmonic to arithmetic mean ratio is still
stable under such conditions, Fig.2b. Note that, the studies [1, 2] worked with a small experimental
micro-grid. Furthermore, weather affects all areas equally which implicitly preserves similar citywide power consumption patterns. For the above reasons, a positive correlation was implicitly
guaranteed in the advanced metering infrastructure (AMI) application. However, this is not the casewith transportation applications. In a transportation CPS, data is collected from the wild, and cities
are a complex mix of narrower lanes and highways. The data is also affected by the uniqueness
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:7of the neighborhoods (e.g. downtown vs uptown) and thus a positive co-variance structure is not
implicitly guaranteed. Therefore, a computationally tractable clustering method is required toachieve invariance.
4.1.2 Deviation Under Incidents. Another key achievement of [1] and [2] is that it proved thatany short-lived disturbance on the covariance structure will lead to deviation in any metric thatcombines Harmonic and Arithmetic Mean calculated from a highly positively correlated set of
random variables. This is attributed to an asymmetry in Schur Concavity properties. The Harmonic
Mean is strictly Schur Concave while Arithmetic Mean is Schur Convex. This imbalance causes
deviations in the HM to AM ratio metric whenever any event triggers a decrease in the correlation.
4.1.3 Domain-Specific Challenges: We need the following domain-specific adaptations: First,
in [1], the strength of the positive correlation was implicit in smart metering CPS. However, intransportation CPS, several localized factors affect traffic data patterns in sub-areas of the city. Thisrequires intelligent clustering that preserves a high space-time covariance structure strategically.
Second, [1] was designed for power consumption data from smart meters for a small experimental
micro-grid. In such applications, geospatial factors play little role, which is not the case withcity-wide smart transportation CPS. This requires bounding the clustering region size. Third, AMI
application had only one observation per hour and the framework proposed was suitable for attack
detection and not incident detection. The time for detection of attacks was in the order of hours.
In our CPS use case, incident detection needs to happen within minutes. This requires too manydetection rounds, increases the false alarm reduction challenge. Fourth, in [1] the data was free fromanomalies due to a controlled environment of an experimental micro-grid. Instead, this application
contains data from the wild from a real city and therefore framework adaptions are necessary tolearn the underlying structure of the benign pattern of the CPS.
4.2

Region Growing Clustering Algorithm

This is the main theoretical core of the contribution which mainly addresses the first two challenges.
We need to strategically group the road segments into spatial clusters such that the speed data
has maximum positive correlation which leads to the highest invariance. At the same time, theclustering needs to be geographically proximate for disturbances in the co-variance structure tohave a causal link to the traffic incidents.
All the road segments exhibiting correlations above a threshold may be grouped together to forma cluster. Thereafter, if 𝐶 = {𝑐 1, ...𝑐𝑘 , ...𝑐 𝐾 } is a candidate cluster set and 𝑠𝑖 and 𝑠 𝑗 are any two roadsegments where 𝑖 ≠ 𝑗, 1 ≤ 𝑖, 𝑗 ≤ 𝑛 such that 𝑠𝑖 and 𝑠 𝑗 are in the same cluster 𝑐𝑘 , we can formalize
the problem as the following:
∑︁
∑︁
max
𝐶𝑜𝑟 (𝑠𝑖 , 𝑠 𝑗 )𝑐 ∈𝐶{𝑠𝑖 ,𝑠 𝑗 } ∈𝑐(1)s.t. 𝐶𝑜𝑟 (𝑠𝑖 , 𝑠 𝑗 ) > 𝑝 (𝑚𝑖𝑛)In the above optimization 𝐶𝑜𝑟 (𝑠𝑖 , 𝑠 𝑗 ) represents the correlation between two road segments and𝑝 (𝑚𝑖𝑛) is a threshold. The above optimization problem is 𝑁 𝑃 hard since with |𝑆 | number of roadsegments, there is an exponential number of possible solutions which is computationally intractable.
We need an approximation to the exact solution. This is done by first converting the clustering
problem into a graph problem.
Reformulation into a Graph Problem We convert our optimal clustering problem into agraph problem, where we visualize each road segment as a vertex on the graph 𝐺 ′ and the roadsegment connections as an edge. The weight of an edge is equal to the correlation between the
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:8

Md. J. Islam and J.P. Talusan, et al.

Fig. 3. Problem reformulation to graph problem.

road segments (vertices) it connects. Fig. 3 shows the remapped graph abstraction from the originalroad network to our reformulated graph mapping.
Theoretically, a correlation may exist between any pair of road segments. Therefore the initialgraph 𝐺 ′ is a complete graph. However, since all road segments are not necessarily positively
correlated (e.g. geographically distant, city roads to highways in the same geographical area), there
will be edges with negative or zero weights and relatively low weights. Let there be a bound on theminimum correlation value 𝑝𝑐𝑢𝑡 > 0 necessary to be considered a feasible edge of the graph. Alledges whose weights are less than 𝑝𝑐𝑢𝑡 are pruned from the complete graph. A low 𝑝𝑐𝑢𝑡 affects thelevel of invariance in the ratio invariant which is key to low false alarms and improved detection.
Formally, this reduced graph is denoted as 𝐺 = (𝑉 , 𝐸), where 𝑉 is the set of vertices and 𝐸 is theset of edges. The set of vertices and edges are indexed by 𝑣 ∈ 𝑉 and 𝑒 ∈ 𝐸. Each edge 𝑒 is assigned
a weight 𝑝𝑒 that is equal to the speed correlation between the two road segments which are further
can be denoted in the transformed graph G as the two vertices, such as 𝑣𝑖 and 𝑣 𝑗≠𝑖 , where 𝑣𝑖 , 𝑣 𝑗 ∈ 𝑉 .
Edge Preference Hyper-parameter (𝑝𝑐𝑢𝑡 ): From the feasible set, we introduce a notion of desirability
to form the strongest grouping of clusters. Note, only using Euclidean distance is not appropriate forcausal link because a narrow lane may have a highway road segment running over it. Geographicallythey are close, but even if one incident is affecting a ramp connecting the two, the correlation will
not be as strong due to inherent differences in their physical characteristics. Also, some roads arelong and can see an incident’s effect quickly propagate, and segments not geographically very close
still become affected by that same incident when not geographically close. Hence, we bring in anotion of edge preference hyper-parameter 𝑝𝑐𝑢𝑡 < 𝑝 (𝑚𝑖𝑛) < 1. Using 𝑝 (𝑚𝑖𝑛) separates all edges 𝐸′
into two subsets. We let the set 𝐸𝑠 include all edges whose 𝑝𝑒 < 𝑝 (𝑚𝑖𝑛) while the set 𝐸𝑠 include
those edges whose 𝑝𝑒 ≥ 𝑝 (𝑚𝑖𝑛) . This separation improves causal linkage.
Distance Weight Variable (𝑥𝑒 ): As explained earlier, geographically closer road segments will beaffected by the same incident. Hence, the distance should be factored in the clustering too. Eachedge 𝑒 ∈ 𝐸 can be visualized as associated with a weight variable, 𝑥𝑒 ∈ (0, 1]. The weight 𝑥𝑒 equals𝑒the normalized distance between two vertices (road segments) such that 𝑥𝑒 = 𝑑𝑑𝑑𝑑𝑚𝑎𝑥, where 𝑑𝑑𝑒
is the distances between two vertices of edge 𝑒 and 𝑑𝑑𝑚𝑎𝑥 is the maximum distance among all
distances between any pair of vertices.
Many optimization problems are formulated as error minimization problems where error is anunfavorable outcome that needs to be minimized. In our setting, two kinds of errors happen for
any candidate solution (cluster). First, the two end vertices {𝑣𝑖 , 𝑣 𝑗 } of an edge 𝑒 has correlation
value 𝑝𝑒 > 𝑝 (𝑚𝑖𝑛) but they are in two different candidate clusters (positive error). Second, the two
end vertices {𝑣𝑖 , 𝑣 𝑗 } of an edge 𝑒 has correlation value 𝑝𝑒 < 𝑝𝑚𝑖𝑛 but they are in the same cluster
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:9

negative error. By minimizing these two errors, the optimal clustering can be achieved, maximizing
the correlations in a cluster.
Transformed Optimization Problem: In the graph-theoretic mapping of the original network,
the original optimization can then be re-written as the following:

(a)

(b)

(c)

Fig. 4. (a) A region around the initial node. (b) Increased radius for the region based on the nearest segment
whose one connecting end is inside the previous region. (c) A region where the volume is greater than cut.

arg minh ∑︁𝑝𝑒 𝑥𝑒 +

∑︁𝑝𝑒 (1 − 𝑥𝑒 )i(2)s.t. 𝑥𝑒 ∈ (0, 1]
In the above optimization, the first term includes all positive errors for a given candidate solution
𝐶, when 𝑒 ∈ 𝐸𝑠 (they have 𝑝𝑒 ≥ 𝑝 (𝑚𝑖𝑛) ) and {𝑣𝑖 , 𝑣 𝑗 } are not in the same candidate cluster 𝐶. Similarly,′
the second term includes all negative errors for the same candidate solution 𝐶, when 𝑒 ∈ 𝐸𝑠 (theyhave 𝑝𝑒 < 𝑝 (𝑚𝑖𝑛) ) and {𝑣𝑖 , 𝑣 𝑗 } are in the same candidate solution 𝐶. We need to find the solution 𝐶
which jointly minimizes both errors.
In [4], a clustering problem for a weighted graph which has the same form as Eqn. 2, was solved.
Their study revealed that the relaxed form of the integer problem has an Ω(log 𝑛) integrality gap.
Hence, the best-known factor of the approximation can be 𝑂 (log 𝑛). Hence, it ensures theoreticalguarantees to our approach.
Approximation Algorithm: To understand the core idea of the approximation algorithm, which
is based on growing a region with radius 𝑟 from some random starting point, we first need to definesome key elements.
Region: A 𝑅𝑒𝑔𝑖𝑜𝑛(𝑣𝑖𝑛𝑖𝑡 , 𝑟 ) is the set of road segments 𝑣 ∈ 𝑉 that are within the area with radius 𝑟from an initial vertex 𝑣𝑖𝑛𝑖𝑡 .
Cut: A 𝑐𝑢𝑡 (𝑐), where 𝑐 ∈ 𝐶, is the sum of the weights of the edges 𝑒 ∈ 𝐸𝑠 with 𝑝𝑒 > 𝑝 (𝑚𝑖𝑛) where
each edge 𝑒 has one vertex 𝑣𝑖 , within 𝑐 and the other 𝑣 𝑗 , is outside.
∑︁𝑐𝑢𝑡 (𝑐) =
𝑝𝑒
(3)
𝐶

𝑒 ∈𝐸𝑠

𝑒 ∈𝐸𝑠

′

{𝑣𝑖 ,𝑣 𝑗 }∩𝑐=1
𝑒 ∈𝐸𝑠

Volume: 𝑣𝑜𝑙 (𝑐) is also the total sum of the weights of the edges 𝑒 ∈ 𝐸𝑠 , with 𝑝𝑒 > 𝑝 (𝑚𝑖𝑛) , where atleast one vertex 𝑣𝑖 or 𝑣 𝑗 is inside the cluster.
∑︁𝑣𝑜𝑙 (𝑐) =𝑝𝑒 𝑥𝑒(4){𝑣𝑖 ,𝑣 𝑗 } ∈𝑐𝑒 ∈𝐸𝑠

ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:10

Md. J. Islam and J.P. Talusan, et al.

Algorithm 1: Approximation Algorithm
Input: 𝐺 = (𝑉 , 𝐸 )Output: 𝐶 = 𝑐 1 , ..., 𝑐 𝐾
Initialize: 𝑘 = 1, 𝐶 = { }1 begin2
while 𝐺 𝑖𝑠 𝑛𝑜𝑡 ∅ do3
𝑠𝑒𝑙𝑒𝑐𝑡 𝑟𝑎𝑛𝑑𝑜𝑚 𝑣𝑖𝑛𝑖𝑡 ∈ 𝑉
4𝑟 → 𝑟𝑖𝑛𝑖𝑡
5
while 𝑐𝑢𝑡 (𝑅𝑒𝑔𝑖𝑜𝑛 (𝑣𝑖𝑛𝑖𝑡 , 𝑟 ) ) ≥ 𝑣𝑜𝑙 (𝑅𝑒𝑔𝑖𝑜𝑛 (𝑣𝑖𝑛𝑖𝑡 , 𝑟 ) ) do
6𝑟 → 𝑟 + min(𝑑 𝑣𝑖𝑛𝑖𝑡 ,𝑣 − 𝑟 )𝑣 ∈𝑉𝑣∉𝑅𝑒𝑔𝑖𝑜𝑛 (𝑣𝑖𝑛𝑖𝑡 ,𝑟 )

7
8
9
10
11

𝑐𝑘 = 𝑅𝑒𝑔𝑖𝑜𝑛 (𝑣𝑖𝑛𝑖𝑡 , 𝑟 )𝑖𝑛𝑠𝑒𝑟𝑡 (𝑐𝑘 , 𝐶 )𝑟𝑒𝑚𝑜𝑣𝑒 (𝑐𝑘 , 𝐺 )𝑘 →𝑘 +1
return 𝐶 = 𝑐 1 , 𝑐 2 , ...𝑐 𝐾

The algorithm returns the set of different clusters. The formation of one cluster happens via theregion growing process. Algorithm 1 starts by checking if the graph 𝐺 has more nodes to cluster inline 1. Nodes that are chosen as part of a cluster are then removed from 𝐺. The core of the region
growing approximation are lines 3-7 in Algorithm 1 which decides what is included in one cluster,while the rest of the algorithm repeats the process of finding clusters for the whole graph. It starts
at a random vertex 𝑣𝑖𝑛𝑖𝑡 at line 3 with an initial radius 𝑟 1 = 𝑟𝑖𝑛𝑖𝑡 > 0 in line 4, that forms an initialregion 𝑅𝑒𝑔𝑖𝑜𝑛(𝑣𝑖𝑛𝑖𝑡 , 𝑟 1 ).
In line 5 the algorithm checks if a cluster has met the terminating conditions based on bothcut and volume. In Line 6 the algorithm grows the radius to include the nearest node outside of
the cluster and repeats the process until reaches the terminating criteria. In Line 7 it forms thecluster 𝑐𝑘 , and the formed cluster is added to the final cluster set 𝐶. The nodes in the cluster arethen removed from the graph 𝐺. Then, the algorithm repeats the entire process until all nodes are
clustered. To illustrate the algorithm 1, Fig. 4a represents an initial region centered on a random
vertex (shaded in red) with radius 𝑟 1 . Next, it finds the nearest vertex 𝑣, in the neighborhood of
the initial region. To find the nearest vertex, it lists all other vertices 𝑣𝑛𝑒𝑎𝑟 such that each vertex
in set 𝑣𝑛𝑒𝑎𝑟 includes only vertices that are directly connected to the region 𝑅𝑒𝑔𝑖𝑜𝑛(𝑣𝑖𝑛𝑖𝑡 , 𝑟 1 ). The
term “directly connected to the region” implies that an edge exists between a vertex in 𝑣𝑛𝑒𝑎𝑟 andany vertex inside 𝑅𝑒𝑔𝑖𝑜𝑛(𝑣𝑖𝑛𝑖𝑡 , 𝑟 1 ) (any vertex shaded red). For illustration, in Fig. 4a, 𝑣 1 and 𝑣 2 are
the only vertices directly connected to the 𝑅𝑒𝑔𝑖𝑜𝑛(𝑣𝑖𝑛𝑖𝑡 , 𝑟 1 ).
It then calculates the euclidean distance from 𝑣𝑖𝑛𝑖𝑡 to each vertex in 𝑣𝑛𝑒𝑎𝑟 and selects the vertex
𝑣 that is nearest 𝑣𝑖𝑛𝑖𝑡 . Let this smallest distance (from 𝑣𝑖𝑛𝑖𝑡 ) be denoted as 𝑑 (𝑣𝑖𝑛𝑖𝑡 ,𝑣) . This distanceis used as the radius of the new region such that 𝑟 2 = 𝑑 (𝑣𝑖𝑛𝑖𝑡 ,𝑣) . Since 𝑟 2 > 𝑟 1 , the region grows
from the initial region, hence the term region growing approximation. This is illustrated in Fig. 4b,
where 𝑟 2 is formed by the distance between the nearest vertex 𝑣 1 and 𝑣𝑖𝑛𝑖𝑡 .
Note, that with the new region, the set 𝑣𝑛𝑒𝑎𝑟 now includes 𝑣 2 and 𝑣 3 . For simplicity, we drop thesuffix of the radius parameter such that the radius at any iteration of region growth is simply 𝑟 .
Then, the above process of region growing happens continuously The region then continuously
grows until the stopping condition is met and we have our first cluster 𝑐𝑘 = 𝑐 1 where 𝑘 = 1. We
insert this first cluster into our final cluster set denoted by 𝐶 (See line 8 in Algorithm 1).
All the vertices in cluster 𝑐 1 are then removed from the graph 𝐺 to avoid duplication while
generating other clusters. Finally, 𝑘 is increased by 1 for the next iteration. The process starts again
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:11

with a new initial region centered around a new random vertex 𝑣𝑖𝑛𝑖𝑡 . It generates the cluster 𝑐𝑘 = 𝑐 2
by executing the lines 3 − 7 which is added to the cluster set 𝐶 before removing the vertices incluster 𝑐 2 from graph 𝐺. Algorithm 1 continues until there are no vertices left to cluster. Once thegraph 𝐺 is empty, the set of clusters 𝐶 is returned. The approximation algorithm takes polynomial
time to cluster all the segments. We proved the complexity in our previous work in [12] which is
𝑂 (𝑛) where 𝑛 is the number of road segments.
4.3

Ratio Invariant per Cluster

The clustering process ensures clusters that maximize the correlation strategically. Let any cluster
𝑐𝑘 have |𝑆𝑐𝑘 | number of road segments. Then, we calculate a ratio metric 𝑄𝑐𝑘 (𝑡) for every cluster 𝑐𝑘
at each time index 𝑡, which is the invariant. The ratio metric is defined as the ratio of the harmonicmean 𝐻 𝑀𝑐𝑘 (𝑡) and arithmetic mean 𝐴𝑀𝑐𝑘 (𝑡) of data collected from all road segments within acluster such that:Í |𝑆𝑐𝑘 |
𝑑𝑆𝑐𝑙
𝑆𝑐 𝑘
𝑙=0
𝑘
𝐻 𝑀𝑐𝑘 (𝑡) = Í |𝑆 |𝐴𝑀𝑐𝑘 (𝑡) =
(5)
𝑐𝑘
𝑆𝑐
1
𝑘

𝑙=0

𝑑𝑆 𝑙

𝑐𝑘

where 𝑑𝑆𝑐𝑙 is the aggregate speed reported by the 𝑙-th TMC for a time index 𝑡. Consequently, the𝑘ratio sample of the cluster 𝑐𝑘 at any time index 𝑡 is calculate by the following:
𝑄𝑐𝑘 (𝑡) =𝐻 𝑀𝑐𝑘 (𝑡)𝐴𝑀𝑐𝑘 (𝑡)(6)

To illustrate the importance of the approximation algorithm for clustering to maximize positive
correlation strategically, we compare the plots of the time series of the ratio samples for the sametime and area in Fig. 5. Fig. 5a is a cluster with a high data correlation (0.87) and Fig. 5b is a clusterwith a low data correlation (0.37). Observe that the time series of ratio samples in Fig. 5a is highly
stable under benign traffic conditions (stationarity and low variance) and shows a sharp deviationon the incident that happened at 13:00 hrs. In contrast, Fig. 5b that did not maximize correlation
has poor stability and does not show clear deviation in its time series when the incident happens.

1.0

0.9

0.9

0.8

0.8

0.7

0.4

incident

0.7

0.6
0.5

QratioQ-ratio

1.1

1.0

Q-ratio

1.1

0.6
0.5

Qratio0.4

incident

0.3
00:00 03:00 06:00 09:00 12:00 15:00 18:00 21:00
01-Aug0.3
00:00 03:00 06:00 09:00 12:00 15:00 18:00 21:00
01-Aug

(a)

(b)

Time

Time

Fig. 5. Effect of cluster level correlation on invariance: (a) High correlation (0.87). (b) Low correlation (0.35).

An important thing to note is that every incident is unique in its manifestation and the method
has to generalize for various clusters. Hence, we need to learn the underlying general structure of
the ratio time series. However, the training data collected from the wild have incidents, and the datacollected from connected transportation is very noisy due to human behavioral randomness. ThisACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:12

Md. J. Islam and J.P. Talusan, et al.

is unlike traditional industrial CPS where the data patterns are only governed by tightly modeledlaws of physics. Hence, we cannot simply learn the ratio samples themselves.
4.4

Data Pre-processing and Augmentation

Real-world mobility data pose a practical problem for unsupervised learning problems such as
anomaly detection, due to the presence of various incidents in the training phase. This preventsthe learning of the underlying structure of benign data patterns. We need a mechanism to bypass
this problem which we discuss here.

1.0

0.008

0.8

Probability

0.010

Q-ratio

0.006

0.6

0.004

0.4

0.002
0.000

incident training
replacement

Ground Truth Incident

Region A

Region B

19:00

20:00

0.2

−100

−50

0

50

Duration from GT(t)
(a)

100

0.0

Time

21:00

22:00

(b)

Fig. 6. (a) Distribution of temporal neighborhood of disturbances across all incidents (b) Data augmentation.

The intuition is to use the time and location stamp of the ground truth incidents and superimpose
them on the ratio time series of the cluster which falls under the location of a particular incident.
Then we identify the neighborhood of the time series of 𝑄𝑐𝑘 (𝑡) around all incidents to learn the
portions of the time series that were disturbed. Unless these disturbances are cleaned out, it will
prevent learning the structure of the benign behavior.
Note, ground truth incident recording itself is noisy due to human-in-the-loop issues. We observedin many cases, they are recorded much after the physical world has been affected by the incident. In
other cases, the incident is reported and recorded instantly but it takes some time for the physicalworld to get really affected (e.g. in sparse traffic scenarios).
Temporal Disturbance Period Selection: We know that prominent incidents in the city cause
large congestion that gets captured in the congestion factor metric available with the dataset.
Additionally, the moving average of the invariant decreases near incidents. We utilize the decrease
to differentiate between benign and noisy ratios (invariant) which are then used to select a neighborhood around the incident ground truth timestamp 𝐺𝑇 (𝑡). This can be visualized through Fig. 6bwhere region 𝐵 is such a neighborhood.
From the time stamp where the incident was recorded (minutes=0), we check how many incidents
showed a low moving average of ratio samples for a window before and after (minutes=0). One
can find 𝑎𝑟 % of the incidents create a decrease in the ratio time series for less than Minutes =±𝑦minutes. Hence, the temporal neighborhood of ratio time series sample around the 𝐺 (𝑡) timestampthat needs to be cleaned and discarded is on average y minutes before and after the 𝐺𝑇 (𝑡) shownin Fig. 6a. In Fig. 6b, this region is marked as region B. In Fig. 6b, we showed region markings
assuming, 𝑎𝑟 = 60%, the 𝑦 is around 30 minutes before and after any corresponding ground truth
time stamps. Similarly, any confidence interval can be used for cleaning.
Ratio Sample Cleaning and Augmentation: To clean the incident neighborhood, we discard
the ratios of region B from the training examples of ratio samples and replace them with the
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:13

cumulative moving average (CMA) of an equal length of time just before the start of region B(temporal disturbance window of a cluster).
As an illustration, the cumulative sliding moving average of the ratio samples from region A are
copied into the discarded ratio samples from region B, as demonstrated by Fig. 6b. The CMA for
any cluster 𝑐𝑘 at time 𝑡, is calculated by the following:
𝑄𝑐𝑀𝐴
(𝑡) =
𝑘

(𝑡 − 1)𝑄𝑐𝑀𝐴(𝑡 − 1) + 𝑄𝑐𝑘 (𝑡)
𝑘

(7)

𝑡

This process is executed for all ratio sample neighborhoods of ground truth incidents found in
all clusters during the training phase. The CMA of region A is then used to replace the signature in
region B. Figure 6b shows the incident signature being replaced by the cleaned data. This allows
the model to learn the underlying structure of the data without incidents.
4.5

Detection Framework Design

After cleaning effects of ground truth recorded incidents, there are other behavioral randomness
and noise that make lowering false alarms challenging without sacrificing the detection accuracy.
Therefore, a two-tier approach (NIST recommended [7]) to learning the thresholds and an appropriate anomaly detection criterion is essential. The two-tier principal mandates short-term andlong-term errors of any underlying detection metric. We adapt this idea in our context in thefollowing manner:
4.5.1 First Tier Stateless Residuals. The first tier uses the time series distribution of the ratios
𝑄𝑐𝑘 to set up a varying threshold that follows the ratio distribution for each cluster 𝑐𝑘 where𝑘 ∈ {1, · · · , 𝐾 }. A particular ratio 𝑄𝑐𝑘 (𝑡) can be greater than or less than the mean ratio 𝑄𝑐𝑀𝑒𝑎𝑛(𝑡).
𝑘
The acceptable margin creates the upper and lower side boundary using the mean ratio of a clusterℎ𝑖𝑔ℎ
𝑄𝑐𝑀𝑒𝑎𝑛(𝑡) and the standard deviation 𝜎 𝑐𝑘 . The upper boundary is denoted as Γ𝑐𝑘 (𝑡) and the lower𝑘boundary is denoted as Γ𝑐𝑙𝑜𝑤(𝑡). The boundaries are termed as safe margins which can be calculated𝑘using the following equations:
Γ𝑐𝑘

ℎ𝑖𝑔ℎ

(𝑡) = 𝑄𝑐𝑀𝑒𝑎𝑛(𝑡) + 𝜅𝜎 𝑐𝑘
𝑘

(8)

Γ𝑐𝑙𝑜𝑤
(𝑡) = 𝑄𝑐𝑀𝑒𝑎𝑛(𝑡) − 𝜅𝜎 𝑐𝑘
𝑘
𝑘

(9)

4.5.2 Second Tier Stateful Residuals. The second tier consists of two thresholds. These thresholds
are termed as standard limits in [1]. The upper side standard limit is 𝜏𝑐𝑚𝑖𝑛
(ℎ) and the lower side𝑘standard limit is 𝜏𝑐𝑚𝑎𝑥
(ℎ).
Setting
up
these
two
thresholds
is
not
as
straightforward
as tier 1. To𝑘calculate the thresholds, the first step is to get the residuals ∇𝑐𝑘 (𝑡) of each time index using thesafe margin. This residual will be used to again calculate the residual under curve 𝑅𝑈 𝐶𝑐𝑘 (𝑡) over asliding frame of size 𝑆𝐹 for each time index and the sub-region. Finally the 𝑅𝑈 𝐶𝑐𝑘 (𝑡)’s are used tolearn the standard limits by the given algorithm 2.
Residual is defined as the difference between the safe margin and the ratio. If a ratio is higherℎ𝑖𝑔ℎthan Γ𝑐𝑘 (𝑡), the residual will be positive and if a ratio is less than Γ𝑐𝑙𝑜𝑤(𝑡), the residual will be𝑘negative. The following equation calculates the residual:ℎ𝑖𝑔ℎ = 𝑄𝑐𝑘 (𝑡) − Γ𝑐𝑘 (𝑡),𝑙𝑜𝑤∇𝑐𝑘 (𝑡) : = 𝑄𝑐𝑘 (𝑡) − Γ𝑐 (𝑡),𝑘
 = 0,
ℎ𝑖𝑔ℎ

if 𝑄𝑐𝑘 (𝑡) > Γ𝑐𝑘 (𝑡); 𝑙𝑜𝑤
if 𝑄𝑐𝑘 (𝑡) < Γ𝑐𝑘 (𝑡); otherwise;

(10)

Residual Under Curve A non-zero residual indicates the possible presence of an anomaly. However,
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Md. J. Islam and J.P. Talusan, et al.
Average difference between speeds (mph)

1:14
10

0.10

RUC

0.05
0.00

−0.05
−0.10

Residual
τmax
τmin

−0.15

Ground Truth Incidents

06:00

08:00

10:00

12:00

14:00

Time

16:00

18:00

20:00

Fig. 7. Detection illustration: RUC of 𝑖 𝑡ℎ cluster.

i th cluster

8
6
4
2
0
06:00

08:00

10:00

12:00

14:00

Time

16:00

18:00

20:00

Fig. 8. Average difference in speeds.

to confirm, a sum of residuals is calculated over a fixed optimal time window size which can becalled as sliding frame size. The summation is termed as residuals under curve. It is calculated by:
𝑅𝑈𝐶𝑐𝑘 (𝑡) =𝑡∑︁∇𝑐𝑘 (𝑘)(11)

𝑗=𝑡 −𝐹𝑆

Algorithm 2: Calculate 𝜏𝑐𝑚𝑎𝑥
(ℎ)
𝑘
1
2

for 𝑐𝑘 , 𝑡, 𝜏 do
if (𝑅𝑈 𝐶𝑐𝑘 (𝑡 ) < 𝜏 then|𝜏 −𝑅𝑈 𝐶𝑐 (𝑡 ) |𝑘𝑐𝑐𝑚𝑎𝑥 :
2
CC ← 𝑐𝑐𝑚𝑎𝑥

3
4
5
6
7
8

else
𝑝𝑝𝑚𝑎𝑥 = |𝑅𝑈 𝐶𝑐𝑘 (𝑡 ) − 𝜏 |2
PP ← 𝑝𝑝𝑚𝑎𝑥
𝜏𝑐𝑚𝑎𝑥(ℎ) = 𝜂1+ argmin𝜏𝑘ÍCC 𝑐𝑐𝑚𝑎𝑥 −ÍPP 𝑝𝑝𝑚𝑎𝑥

Learning Standard Limit The computed 𝑅𝑈 𝐶𝑐𝑘 is later used to learn the standard limit using
Algo. 2. The algorithm treats the interior and exterior RUC differently by multiplying two different
weights. An interior-point contributes less to the overall loss and an exterior point contributes
more. The algorithm minimizes the difference between the loss of interior and exterior points tolearn the optimal standard limit both for the higher and lower sides. Eventually, both of the learned
thresholds use the same algorithm, here we have shown only for the 𝜏𝑐𝑚𝑎𝑥
(ℎ). For 𝜏𝑐𝑚𝑖𝑛(ℎ) only the𝑘𝑘𝑚𝑎𝑥
negative 𝑅𝑈 𝐶𝑐𝑘 (−) are used whereas for 𝜏𝑐𝑘 only the positive 𝑅𝑈 𝐶𝑐𝑘 (+) are used.
4.5.3 Anomaly Detection Criterion. Similarly, the RUC can be calculated at every time index
in the test set. Let 𝑅𝑈 𝐶𝑐𝑘 (𝑡 𝑐 ) is the RUC value for the cluster 𝑐𝑘 in the test set at the current timeindex 𝑡 𝑐𝑢𝑟𝑟𝑒𝑛𝑡 . Then, the incident detection criterion is given∈ [𝜏𝑐𝑚𝑖𝑛(ℎ), 𝜏𝑐𝑚𝑎𝑥
(ℎ)] No Incident;
𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝑘𝑘
𝑅𝑈 𝐶𝑐𝑘 (𝑡
):
(12)
𝑚𝑖𝑛
∉ [𝜏𝑐𝑘 (ℎ), 𝜏𝑐𝑚𝑎𝑥
(ℎ)],
Incident Inferred;𝑘Fig. 7 illustrates the incident detection where the vertical lines are the ground truth incidents
and the horizontal lines represent the standard limits. we can see that RUC(T) metric goes beyond
the learned standard limit near the growth truth time stamps.
Why do ratios and RUC deviate? As and when an accident occurs within a subarea of a city,
the immediate neighborhood of the location where the accident happened, experiences a decrease
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems
0.104

Missed Detection Rate
False Positive Rate

0.8

0.102
0.100

MAD Value

0.6

Percent

0.098

0.4

0.096
0.094

0.2
0.0

1:15

0.092
0.090

0.25

0.5

0.75

1.0
1.25
Kappa

1.5

1.75

0.1

2.0

Fig. 9. Effect of varying 𝜅 on the detection of a given cluster.

0.2

0.3

0.4

0.5

pcut

0.6

0.7

0.8

0.9

Fig. 10. MAD over 𝑝𝑐𝑢𝑡 for 𝑘 𝑡ℎ cluster.

in vehicle speeds instantly. However, this reduction in speed takes time to propagate beyond this
immediate neighborhood until it affects the whole cluster. This delay in propagation causes thedeviation in the signature and can be detected by the metric as an anomaly. Fig. 8, is an illustrationthat shows that the average difference between any two pairs of TMC values within the identified
clusters over time 𝜉 (𝑇 ) do not vary much, which is required for ratio stability.
4.6

Hyperparameter Tuning

There are four different hyperparameters. The first set of hyperparameters is 𝑝𝑐𝑢𝑡 and 𝑝 (𝑚𝑖𝑛) which
affect the clustering process and the distribution of ratios. The second parameter set includes 𝜅 and𝑆𝐹 values which affect the standard limits. The distribution will produce different standard limits
for every combination of 𝜅 and 𝑆𝐹 the same ratio. To learn the 𝑝𝑐𝑢𝑡 , we measure the deviation ininvariance (ratios) 𝑄 𝑀𝐴𝐷 for different margins of 𝑝𝑐𝑢𝑡 . The 𝑄 𝑀𝐴𝐷 is used to select the 𝑝𝑐𝑢𝑡 value
since it directly affects the level of invariance in the ratiometric. Since the lowest median absolutedeviation in the series imply the most stability, it implies that the smallest 𝑝𝑐𝑢𝑡 for which theminimum value of 𝑄 𝑀𝐴𝐷 stops decreasing across consecutive values of candidate 𝑝𝑐𝑢𝑡 is desirable.
The smallest value of 𝑝𝑐𝑢𝑡 , shown in Fig. 10, is recommended since too much positive correlation
reduces the sensitivity to smaller incidents.
As we increase 𝑝𝑐𝑢𝑡 from 0.0 to 0.99, the mean absolute deviation of the ratios in a cluster
decreases. This trend continues until a certain point where the deviation stabilizes. Accounting this,we settle on 𝑝𝑐𝑢𝑡 = 0.7 as a lower bound. The hyper-parameter 𝑝 (𝑚𝑖𝑛) controls the area coverageand the performance of the cluster-wise incident detection. We learn it by the following:arg max (𝐶𝑜𝑣𝑅 + 𝑇 𝑃𝑅 − 𝐹 𝑃𝑅)
𝑝 (𝑚𝑖𝑛)(13)

where 𝐶𝑜𝑣𝑅 is the coverage rate of road segments while clustering, 𝑇 𝑃𝑅 true positive rate ofdetection, 𝐹 𝑃𝑅 are the false positive rate, for incident detection. The above equation ensures the
maximization of performance by reducing the false positives at the same time and increasing thecoverage percentage. The approximate solution cover a majority of the area which had a total of
6,928 road segments. As we increase the correlation threshold, 𝑝 (𝑚𝑖𝑛) , the radius becomes smallerwith fewer segments being included in each cluster. However, as a result, there are now moreclusters generated, resulting a larger coverage of the target area. The performance for different
values for 𝑝 (𝑚𝑖𝑛) are given in Table 2.
For the incident detection model we learn the optimal value for the hyper parameters 𝜅 and 𝑆𝐹 .
The parameter 𝜅 is a value in (0, 3) and 𝑆𝐹 is a sliding frame size which varies among the integer
values in the set {3, 5, 7, 9} . The optimal values of 𝜅 and 𝑆𝐹 are learnt byarg min
𝜅,𝑆𝐹(𝑀𝐷𝑅 + 𝐹 𝑃𝑅)(14)

ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:16

Md. J. Islam and J.P. Talusan, et al.

where 𝑀𝐷𝑅 is the missed detection rate and 𝐹 𝑃𝑅 is the false positive rate. The 𝜅 and 𝑆𝐹 from theabove equation are selected as the optimal hyperparameters for the considered cluster. Each cluster
𝑐𝑘 has its own optimal selection of 𝜅 and 𝑆𝐹 . Figure 9 shows the effect of the hyperparameter 𝜅on the detection performance of the 𝑘 𝑡ℎ cluster. A detection model with large 𝜅 reduces the total
number of false alarms however, it increases the number of missed incidents detected. The equation
above ensures an acceptable performance.
Table 2. Cluster Information.

Cluster Info

5

0.75

𝑝 (𝑚𝑖𝑛)0.85

Table 3. 𝜇𝑐𝑀𝑅
based rating levels.
𝑘
Scenario

0.95

𝑆𝑐𝑙

Discrete Rating Level (rl)

count
mean
min
max
ave. data correlation
ave. radius (m)area coverage

317
16.06
4
161
0.863
611.11
73%

354
14.54
4
132
0.862
610.07
74%

472
11.81
4
112
0.867
490.0380%

2Δ𝑎𝑏𝑠 < 𝜃𝑑𝑖 𝑘𝑓 𝑓 ≤ 3Δ𝑎𝑏𝑠

2

True Positive Rate
False Positive Rate

0.916092
0.010727

0.924653
0.032051

0.926340
0.030957

otherwise1

𝜃𝑑𝑖 𝑘𝑓 𝑓 ≤ Δ𝑎𝑏𝑠
𝑆𝑐𝑙𝑘

Δ𝑎𝑏𝑠 < 𝜃𝑑𝑖 𝑓 𝑓 ≤ 2Δ𝑎𝑏𝑠
𝑆𝑐𝑙

4
3

SEGMENT LEVEL DETECTION

The primary motivation of segment-level incident detection is to provide optimal resource dispatch
for agents that manage traffic incidents. Though cluster-level detection [12] helps to identify
the incidents in a city block, it can not provide the exact location. Without this information, theemergency response will be more time-consuming since the department would need to manually
locate the road affected or wait for citizens to call in. For example, one cluster can have 15 roadsegments within 600 meter radius. The agents would have to allot resources to locate the incident
area before providing the service. In this context, our segment-level detection actually brings
locality which allows them to respond and dispatch optimally. Another motivation was to address
the scalability issue. A city or a region can have a vast number of road segments. For example,
a region with 500 roads would require 500 detection models introducing plenty of overhead. In
contrast, the cluster level first detects whether there is a problem with only one detection model.
Then, only if there is an anomaly at the segment level would detection be initiated. Hence, thesegment-level detection runs only when needed, it saves computation and networking resourcesthat support such computations.
The incident detection framework detects incidents at a cluster level where each cluster 𝑐𝑘 is
made up of a set of 𝑆𝑐𝑘 segments. Any segment 𝑆𝑐𝑙 𝑘 (𝑙 𝑡ℎ in cluster 𝑐𝑘 ) can be the origin of an actual
incident that cannot be identified in the current framework. This motivated us to propose an
extended architecture that would locate the origin of the incident after detecting the incident at acluster level. The extended framework is a significant amendment to the prior framework that canhelp commuters save time by identifying segments and generating alternative routes.
For segment-level detection, first a cluster 𝑐𝑘 under a detected incident at time 𝑡 is considered
to narrow down the origin (segment 𝑆𝑐𝑙 𝑘 ) of the incident. Intuitively, in the cluster 𝑐𝑘 , only theaffected segment will show a change(decrease) in speed values. Therefore,the average speed 𝜇𝑐𝑘 (𝑡)at time 𝑡 also changes(decreases) for the cluster 𝑐𝑘 . However, an approximation of the unaffected
mean 𝜇𝑐𝑘 (𝑡) of the cluster 𝑐𝑘 can identify the affected road segment(where (segment)the incident
occurred) by comparing the road segment speed 𝑑𝑆𝑐𝑙 (𝑡) with the mean. Therefore, we approximate𝑘the value of mean 𝜇𝑐𝑘 (𝑡) and denote it as corrected mean by 𝜇𝑐𝑀𝑅(𝑡). Further, depending on the𝑘
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:17

distance of the speed of each segment 𝑑𝑆𝑐𝑙 (𝑡) with 𝜇𝑐𝑀𝑅(𝑡) at time 𝑡, each road segment is assigned𝑘𝑘a rating level between 1 to 4. We calculate the rating for a time window of length 𝑇 for a segment.
A key point is that throughout the window 𝑇 the approximation of mean 𝜇𝑐𝑀𝑅(𝑡) is calculated only𝑘once at 𝑡. Using the rating values, a trust score per segment is calculated which lets us apply binaryclassification to identify impacted and non-impacted road segments in cluster 𝑐𝑘 under the incident.
Overall, segment-level detection is divided into two major parts, the first is approximating thecluster mean and the second is the trust scoring model. The entire architecture is redefined todetect the road segments under incident for smart transportation problems from [3], where a trustscore-based framework was first used to detect individual smart meters under attack.
5.1

Mean Approximation at Incident Detection Time

It is proven in [3] that if the individual values are decreased due to a data poisoning attack orother means, the means (arithmetic and harmonic) are also decreased. In such cases, the generalmean 𝜇 can be corrected by using 𝜇 = 𝐴𝑀 + (𝐴𝑀 − 𝐻 𝑀). Similarly, due to incidents, we assume
that the speed values are decreased in a smart transportation system. Thereafter, the same mean
correction can be utilized to get the approximation of the mean. Therefore, the correction is doneby 𝜇𝑐𝑀𝑅(𝑡) = 𝐴𝑀𝑐𝑘 (𝑡) + (𝐴𝑀𝑐𝑘 (𝑡) − 𝐻 𝑀𝑐𝑘 (𝑡)), and is calculated each time 𝑡 when an incident is𝑘detected. Further, to measure the distance of the speed values 𝑑𝑆𝑐𝑙 (𝑡) of each segments 𝑆𝑐𝑙 𝑘 (𝑡) the𝑘unbiased standard deviation is approximated. For corrected standard deviation, we simply takethe standard deviation in the window just before the occurrence of the incident. That means forstandard deviation correction we will take 𝜎𝑐𝑀𝑅(𝑡) = 𝜎𝑐𝑘 (𝑡 − 1). Again, both the approximated𝑘𝑀𝑅𝑀𝑅mean 𝜇𝑐𝑘 (𝑡) and standard deviation 𝜎𝑐𝑘 (𝑡) are calculated only once for the window length 𝑇 .
5.2

Trust Scoring Model

Trust scores for each segment are computed in three steps. First, the corrected or approximatedmean 𝜇𝑐𝑀𝑅(𝑡) and standard deviation 𝜎𝑐𝑀𝑅(𝑡) are used to calculate ratings per segment for a fixed𝑘𝑘window of length 𝑇 = 30 minutes. Second, each segment is assigned a weight based on the ratingsreceived in the first step. Lastly, the weights are converted to a trust score using an inverse powerlaw (IPL) based kernel function. In the following, each of these three steps is presented in detail.
5.2.1 Segment Rating. A rating of a segment 𝑆𝑐𝑙 𝑘 at time 𝑡 is an integer value that points out thecloseness of the speed value 𝑑𝑆𝑐𝑙 (𝑡) of the segment with the corrected mean value 𝜇𝑐𝑀𝑅(𝑡) of a𝑘𝑘cluster 𝑐𝑘 . A segment can have a different rating value based on time 𝑡 within a window 𝑇 . As
mean value 𝜇𝑐𝑀𝑅(𝑡) stays same for such window 𝑇 , we remove the 𝑡 from 𝜇𝑐𝑀𝑅(𝑡) thus it becomes𝑘𝑘𝑀𝑅
𝜇𝑐𝑘 . Now, we calculate the absolute difference between the 𝑑𝑆𝑐𝑙 (𝑡) and 𝜇𝑐𝑀𝑅
which
is denoted as𝑘𝑘Δ𝑎𝑏𝑠 = 𝜎𝑀𝑅 for distance under one standard deviation. Thereafter, if the distance for a segment
𝑆𝑐𝑙𝑘
𝑆𝑐𝑙 𝑘 in cluster 𝑐𝑘 has absolute difference denoted by 𝜃𝑑𝑖 𝑘𝑓 𝑓 = 𝑑𝑆𝑐𝑙 (𝑡) − 𝜇𝑐𝑀𝑅
, the segment get the𝑘rating according to the policy presented by Table 3. The Table 3 shows 4 different rating values
where the highest rating is 4 which is closest to 𝜇𝑐𝑀𝑅
, and 1 is the lowest rating that denotes the𝑘furthest distance from the 𝜇𝑐𝑀𝑅.
The
ratings
forallsegments
in a cluster 𝑐𝑘 are represented by a
𝑘

𝑆𝑐𝑙

𝑙

𝑘
rating vector 𝑟 𝑆𝑐𝑘 . This rating vector is sorted as 𝑟𝑠𝑜𝑟𝑡
in ascending order to assign weights to eachsegment depending on the rating value.

𝑆𝑐𝑙𝑘
5.2.2 Segment Weight. Here each rating 𝑟𝑠𝑜𝑟𝑡, for each segment 𝑆𝑐𝑙 𝑘 in cluster 𝑐𝑘 at time 𝑡 gets a𝑙weight based on the Gaussian rating distribution . The weight for each rating are denoted as 𝑊 𝑆𝑐𝑘
The lower rating gets lower weights and vice versa. After getting the rating, they are sorted which
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:18

Md. J. Islam and J.P. Talusan, et al.

makes it easier to give lower weights to smaller ratings through Eqn 15 by dividing the rating spaceover the considered time window. In the following, a weight parameter 𝑥 (𝑡) distributed between 1
to 4 is calculated as:(𝑇 𝐾 − 1)𝑡
𝑥 (𝑡) = 1 +∀ 𝑡 = 0, 1, ...,𝑇 − 1
(15)
𝑇 −1
In Eqn 15, 𝑇 𝐾 = 4 denotes the total number of discrete rating levels in the system and 𝑇 is thelength of the window. Now, 𝑥 (𝑡) from Eqn 15 is used to calculate the distance from highest ratinglevel 𝜇𝐵𝑅 = 4. The larger the distance, the smaller the density value on the shape of the Gaussian
distribution. The following equation is used to calculate the density as weights
𝑐𝑤 𝑆𝑐𝑙 (𝑡) =
𝑘

1
√
𝑑𝑟

𝜎 𝑙

𝑆𝑐 𝑘

2𝜋

exp −

(𝑥 (𝑡) − 𝜇𝐵𝑅 ) 2
2
2 𝜎 𝑑𝑟𝑙

(16)

𝑆𝑐 𝑘

In Eqn 16, 𝜎 𝑑𝑟𝑙 represents the standard deviation of discrete ratings of each segment from 𝜇𝐵𝑅 = 4
𝑆𝑐 𝑘

in a window length 𝑇 . The 𝜎 𝑑𝑟𝑙 for each segment will be different based on different observations
𝑆𝑐 𝑘compared to common mixture data, which captures certain individual differences in speed. The raw
𝑐𝑤𝑆 𝑙 (𝑡 )𝑐weights from Eqn. 16, are normalized using 𝑤 𝑆𝑐𝑙 (𝑡) = Í𝑇 −1 𝑐𝑤𝑘
𝑡 =0

𝑘

𝑆𝑐𝑙
𝑘

(𝑡 )

Thereafter through an indicatorfunction which is denoted as 𝐼 (𝑟𝑙, 𝑡), all the weights 𝑤 𝑆𝑐𝑙 (𝑡) corresponding to each unique rating𝑘level 𝑟𝑙 (between 1 to 4) withing 𝑇 are added up to get the cumulative weight associated with eachrating level. The cumulative weight is by 𝑊 𝐷 (𝑟𝑙) which is calculated by the following equation
𝑊 𝐷 (𝑟𝑙) =

𝑇∑︁
−1

𝑤𝑡

𝐼 (𝑟𝑙, 𝑡)

(17)

𝑡 =0

In Eqn 17, the indicator function 𝐼 (𝑟𝑙, 𝑡) indicates whether any particular rating level 𝑟𝑙 occurs in
the time slot 𝑡 and It is written as below(
1 if rl occurred in time slot t
𝐼 (𝑟𝑙, 𝑡) =(18)
0 otherwise
The cumulative weights for each rating level are considered together to get the final weight
for each segment 𝑆𝑐𝑙 𝑘 . The weight is presented by 𝑅𝑆𝑐𝑙 of the 𝑙 𝑡ℎ segment of 𝑘 𝑡ℎ cluster which is a𝑘continuous value between 1 and 4 and is given by:
5.2.3 Trust Score. The weight ( 𝑅𝑆𝑐𝑙 ) is used to get the trust score by injecting into an inverse𝑘power law (IPL) based kernel function. IPL transforms the weight value based on magnitude such astransforming larger weight to increase it, and vice versa. If a weight 𝑅𝑆𝑐𝑙 is 4, it means the highest𝑘trustworthiness which is followed by an exponential ‘discounting’ of trust, as 𝑅𝑆𝑐𝑙 decreases. The𝑘final transformed value which is the trust score is between 0 and 1. The following Eqn 19 is showing
IPL based trust score calculation.
(𝑅𝑆𝑐𝑙 )𝜂𝑘
, 𝑇 𝑅𝑆𝑐𝑙 ∈ {0, 1}(19)𝑘(𝑇 𝐾)𝜂
In Eqn 19, 𝜂 is a scaling factor controlling the rate of discounting. The Eqn 19, gives exponentiallyless trust to 𝑅𝑆𝑐𝑙 as it decreases from the maximum value of 4, in adherence to the Folded Gaussian𝑘shape of the rating distribution of road segments. The scaling factor 𝜂 depends on the skewness offolded Gaussian in the benign data set. The Eqn. 19 produces trust values such that impacted and

𝑇 𝑅𝑆𝑐𝑙 =𝑘

ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:19

non-impacted road segments by any incident are linearly separable, which enables the calculationof an unsupervised threshold for classification.
5.2.4 Hyper Parameter Learning. The segment-level detection performance depends on the parameter 𝜂 which is a scaling factor learned by maximizing the following objective𝜂 ∗ = arg max (𝑇 𝑃𝑅)

(20)

∀𝜂>0

In Eq.20 we select 𝜂 that maximizes TPR. Further, the optimal value of 𝜂 is data specific, hence
we tune the parameter for each separate data. Parameter tuning was done by searching througha range from 0.25 to 10.0 in steps of 0.25 with 5 cross-fold validation. We selected the 𝜂 with thehighest overall accuracy across the cross-validation set for four cities.
5.3

Classification and Evaluation

Here, the trust scores are considered as input to an unsupervised binary classifier. The binaryclassifier learns a threshold based on which it divides the input set into two classes. The threshold islearned through a binary classification algorithm, K-medoids [19], which is quick and sufficient for
the use case. The binary classification has two outputs either 0 or 1 where 0 refers to non-impacted
and 1 denotes an impacted road segment. After getting the classification results of the road segments,
the performance is evaluated utilizing ground truth information. In the following neighborhood
selection, validation, and performance evaluation are presented as part of the evaluation.
Neighborhood Selection: The neighborhood area is identified by locating the real road segment
where the incident occurred utilizing ground truth data. The neighborhood consists of several roadsegments depending on max hop distance ℎ𝑑 which is an integer denoting the distance from theorigin. Here, the max ℎ𝑑 is set to 4 to get a neighborhood around the incident. Each segment in
the neighborhood has a level of ℎ𝑑 between 1 to 4 based on the hop distance from the incident
segment. The lowest hop distance is 0 which denotes the incident segment itself. The highest hop
distance is ℎ𝑑 = 4 denoting the furthest segments. In the following, the neighborhood segments
are, compared against the segments classified as impacted by the binary classifier.
Validation: To validate, the segments with hop distances 0 and 1 are considered to keep it more
realistic. Intuitively, it requires significant time to propagate the incident impact up to ℎ𝑑 = 4.
Therefore, considering level>1 will unnecessarily create a bias in validation. The classification ofa road segment is considered as success if it is classified as impacted and the road segment in
the neighborhood has a hop distance of either 0 or 1. On the contrary, the classification of a roadsegment is considered as failure if the road segment is classified as impacted but the road segment
is either not found or has a hop distance greater than 1.
Performance Evaluation: For performance evaluation, we simply count the successful classification while validating.For example if a cluster 𝑐𝑘 has |𝑆𝑐𝑘 | number of segment, the binaryclassification will have 𝑆𝑐𝑘 predictions. Now, out of |𝑆𝑐𝑘 |, if the number of successful classification
𝑆𝑢𝑐𝑐
is 𝑐𝑜𝑢𝑛𝑡𝑆𝑢𝑐𝑐 , the success rate is defined as SR = 𝑐𝑜𝑢𝑛𝑡
|𝑆𝑐 |
𝑘

6

EXPERIMENTAL EVALUATIONIn this section, we introduce the dataset and the performance of our framework in terms of incident
detection rate, false alarm rate, time to detection of incidents, and the impact of undetected incidents
on the CPS application.
6.1

Details of Dataset

To evaluate our framework, we use one year-long traffic data collected from the city of Nashville,
Tennessee by roadside sensors at five-minute intervals. This dataset is bigger in duration and
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:20

Md. J. Islam and J.P. Talusan, et al.

in terms of coverage area compared to validation used in existing works [10, 21]. For groundtruth incidents, we use another dataset collected from Nashville’s Fire and Emergency Response
Department during the same year. In all phases of the experiment, we only consider weekdays of
2019, focusing on the period between 6:00 AM to 9:00 PM since during the weekend and late night,
the traffic has discrete patterns. Details are shown in Table 4.
Table 4. Detail of datasets.
Data Sources

Properties

Nashville

MemphisKnoxville

Chattanooga

Road network# intersections# streets

6,928
19,493

12,07527,144

5,513
11,786

5,121
11,532

Traffic incidents
Sensors# instances# count

8,116
6,928

12,11312,081

4,745
6,073

3,470
5,560

2019

2021

Collection period

2020

Experimental Setup: The twelve months of data is divided across training, cross-validation,and testing sets. The training phase learns the model for a combination of hyperparameters The
cross-validation set is used to find the best hyperparameters that give the best outcome. The best
hyperparameters are fitted to the model to find the final learned model that is used for testing. The
first eight months (Jan. to Aug.) are used for training. The next two months, (Sept. and Oct.), used
for cross-validation. The final two months (Nov. and Dec.) are used to test the model itself.
6.1.1 Training Dataset Details. To train, we focus on Southwest(−87.0506, 35.9895) and Northeast(−86.5275, 36.4168). Then, the segments inside the area are clustered following the clustering
process discussed in Sect. 4.2. To cluster, the road segments from the transformed graph problem
where correlation 𝑝𝑒 is assigned as the weights for each edge 𝑒 in the graph. We consider a cutoff
correlation value 𝑝𝑐𝑢𝑡 and a minimum level of correlation 𝑝𝑚𝑖𝑛 which leads to more invariance.
The clustering generated 354 clusters. However, we found that most incidents in the ground truth
were limited to fewer clusters that correspond to the busiest parts of the city. Hence, we selected
25 clusters with the most reported incidents from the ground truth and used it to evaluate theperformance of incident detection. Limiting the dataset gives us the opportunity to validate ourframework quickly since testing on clusters with little or no incidents would not prove whether
our method can actually detect various kinds of incidents in accident-prone areas. For the temporal
disturbance due to the ground truth incidents, the 𝑎𝑟 = 60% was used from the distribution ofduration from 𝐺𝑇 (𝑡) variable. This corresponds to 𝑦 = ±30 minutes around the neighborhood fromall 𝐺𝑇 (𝑡). Since 𝑡 is slotted every 5 minutes, there are 12 ratio samples augmented per incident.
6.1.2 Testing Dataset Details. This section evaluates our decentralized implementation of a lightweight anomaly detection framework. We tested two months (Nov. and Dec.) from the dataset. In
these two months, there were a total of 851 incidents recorded in 580 active segments. We present
how our technique enables us to detect these incidents which can lead to actionable information.
6.2

Performance Results

This section shows a sensitivity analysis of our performance to changing hyperparameters insteadof learned hyperparameters. Then, we report the performance of the optimal learned model withhyperparameters using the fitness function described in Section 4.6.
6.2.1 Sensitivity Analysis of Performance. Here we give the sensitivity analysis of performance
where the 𝑘𝑎𝑝𝑝𝑎 is not learned but varied as a free parameter to check its effect on the changing
performance. The performance metrics include time to detection, true positive rate, false-positive
rate, expected time between false alarms, and impact of undetected incidents.
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

0.9

True Positive Rate

0.8
0.7
0.6
0.5
0.4

ROC

0.3
0.005 0.010 0.015 0.020 0.025 0.030

Time between false alarms (minutes)

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:21

500
400
300
200
100

False Positive Rate

0

0.25

0.5

0.75

(a)

1.0

1.25

1.5

Kappa

1.75

2.0

(b)

Fig. 11. (a) Average ROC curve for 25 clusters (b) Mean time between false positives with different 𝜅 for 𝑐𝑘 .
1.0

Minutes saved per year

5000

0.8

4000

0.6

CDF

3000

0.4

2000

0.2
0.0

1000
0

10
20
30
40
50
60
Time to detect incidents (minutes)

0.25

(a)

0.5

0.75

1.0

1.25

Kappa (κ)

1.5

1.75

2.0

(b)

Fig. 12. (a) Time to detect incidents (b) Impact on route travel times under different 𝜅.

Detection Rate and False Alarms Fig. 11a shows the average ROC curve across the 25 clusters,
underscoring the performance of our framework. One can see that at 90% true positive detection
rate, the false alarm/false positive rate (FPR) is only 0.030. The low FPR is a significant achievement because: (1) anomaly detection methods are prone to false alarms (2) due to lower rates of
emergency/incidents, the cost of FPR is usually high for any CPS. Each cluster has 16,560 detection
attempts and the false positives are few. Fig. 11b shows the rarity of these false positives even whenusing 𝜅 = 0.25 which has the best overall detection rate. Specifically, Fig. 11b gives an idea on the
expected time between two false alarms for various 𝜅.
Mean Time to Detection A key performance indicator of usability in a CPS application is time todetection. Fig. 12a shows that 78% incidents were detected in the first 5 minutes and 90% incidents
were detected within 30 minutes. Quick time to detection is essential to warn commuters earlier
and control the flow of traffic to prevent congestion spread.
Impacts of Undetected Incidents By successfully detecting an incident, a traffic congestion will
be prevented or mitigated and the commuters will be diverted to different routes, avoiding travel
delays. We generated 200 routes that pass through segments with incidents and computed the travel
time of each route with incidents. Using the true positive rates, we calculated what percentage ofincidents will be detected. Detected incidents allow the system to notify commuters early and thus
will experience less delay in their travel times. The time saved per vehicle is given by the following:△𝑡 = 𝑅𝐷 (

1
1
−
)
𝐼𝑆 𝐹𝑆

(21)

where 𝑅𝐷 is the total distance in a given route, 𝐼𝑆 is the speed due to incidents and 𝐹𝑆 is thefree-flow speed which is experienced when affected areas are avoided. Assuming on average that
10,000 vehicles pass by any segment per year, we can identify the impact of our anomaly detector
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:22

Md. J. Islam and J.P. Talusan, et al.

on the travel time saved over a year. Figure 12b show the amount of travel time saved on a macrolevel depending on the hyperparameter and granularity used respectively.
6.2.2 Overall Performance with Learnt Hyperparameter. We applied the learnt values of 𝜅, 𝑆𝐹 , 𝑝𝑐𝑢𝑡
and 𝑝𝑚𝑖𝑛 for each of the 25 clusters. The final result is an average from all 25 clusters for all traffic
incidents over the 2 months. The average true positive detection rate 𝑇 𝑃𝑅 = 0.90 and 𝐹 𝑃𝑅 = 0.03.
6.3

Comparative Analysis

We evaluate our framework further by comparing it with different baseline models and by using
the additional dataset from three of the largest cities in the state of Tennessee.
We compare our cluster-level incident detection against prediction-based (e.g., AR, DeepLogLSTM) and reconstruction-based (e.g., AE, VAE, EncDec-AD) anomaly detectors. The models aretrained on harmonic speed data instead of the ratio metric we previously derived. We identify
incidents on the time series data per segment per cluster and augment them following the method
described in Section 4.4. Figures 13a and 13b compare the true positive rates and false positive rates
of cluster-level detection of the framework and baseline models. All results are obtained using the
optimal hyperparameters 𝑆𝐹 and 𝜅 for each cluster. We used the same clustering for all methods.
Hence, clustering-wise there is no performance gap. Later, we show that our detection part of
the framework outperforms the other methods. And it is worth noting that the performance of
our framework depends entirely on the quality of the in-variance which is the end product of theclustering strategy. Our detection model blends with the clustering approach resulting in improved
performance compared to other methods. A summary for each model is as follows:(1) AR [15] models use linear regression to calculate a sample’s deviance from its predictedvalue which is then used to identify outliers. This model can handle multivariate time series
data by training independent linear regression models for each dimension and computes
the anomaly score for each sample. The score can be based on the mean, maximization, or
median weighted deviance of each dimension.
(2) DeepLog-LSTM [6] is a deep neural network that uses LSTM to model a system log as a natural
language sequence. It learns the log patterns from normal executions to detect anomalies.
(3) AE [9] is an autoencoder model base on a multi-layer perceptron (MLP) that encodes data intoa lower dimensional latent space which is then reconstructed back to the original structure
by a decoder. Anomalies are identified based on the reconstruction errors.
(4) VAE is similar to an AE. However, it outputs parameters of a pre-defined distribution in thelatent space for every input.
(5) EncDec-AD [16] replaces the MLP layers in the AE with LSTM layers.
We extend our framework to Memphis, Knoxville, and Chattanooga. These three cities, alongwith Nashville are the four largest cities in Tennessee by population. We use one year of long traffic
data for each city in our training and experimental evaluation. The data provider is the same as
the Nashville dataset. Each contains "real-time" harmonic mean flow speeds. We perform the samedivision of data for training, cross-validation, and testing. Each city generated varying numbers ofclusters, however, we only select 25 clusters with the most number of incidents for each city.
In contrast to the Nashville dataset, incidents for these three cities were gathered from crowdsourced incident reports through Waze. Incidents were limited to only include accident reports.
Data collection periods were also varied across three years, to show the ability of our framework tohandle spatially and temporally different data. Fig. 14a and 14b show the ability of the frameworkto adapt and scale to other areas with minimal to no customization.
Overall, the results show the advantage of our framework against baseline models and whenused in other cities and during different periods. Our framework offers higher TPRs and lower
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:23

0.6

1.0

0.5

0.8

0.4
FPR

TPR

0.6
0.4

0.3
0.2

0.2

0.1

0.0
Ratio metricLSTM AE ANN AE ANN VAE

Cluster level detection

AR

0.0

DeepLog

Ratio metricANN AE ANN VAE

AR

LSTM AE DeepLog

Cluster level detection models

(a)

(b)

Fig. 13. Cluster level detection with various methods: (a) True positive rates and (b) False positive rates.
0.175

1.0

0.150
0.125

0.6

FPR

TPR

0.8

0.4

0.075
0.050

0.2
0.0

0.100

0.025
0.000
Nashville

Chattanooga

CitiesMemphisKnoxville

Nashville

Chattanooga(a)CitiesMemphisKnoxville(b)

Segment Level Detection Accuracy

Fig. 14. Cluster level detection in various cities, TN: (a) True positive rates and (b) False positive rates.
1.0
0.8
0.6
0.4
0.2
0.0

(a)

Chattanooga

Knoxville

City

MemphisNashville(b)Fig. 15. (a) Cluster and segment level mapping. (b) Segment level detection accuracy cities in TN.

FPRs across the board. Proving that it can be effective in the first step to not only detect incidents
but to identify the extent of its effects both temporally and spatially.
6.4

Segment Level Incident Detection

After an incident is detected in any cluster, the framework must then be able to detect incidents atthe segment level. However, incident reports are inherently noisy, not always aligning perfectly withsegment information (harmonic speed data) which in contrast, is gathered in real-time. Reported
incidents may be incomplete or incorrectly logged, they may also be reported late. Thus, to beable to correctly evaluate the framework, we have to specify a region around the reported incident
origin to search. Detections are considered successful if the framework can detect segments locatedACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

1:24

Md. J. Islam and J.P. Talusan, et al.

at or one hop away from the ground truth as discussed in Sec. 5.3. The left side of Fig. 15a shows anincident detected in one cluster within Nashville and on the right a group of segments in that cluster,
we highlight the reported point of origin and its adjacent segments in red. We use the results fromthe cluster-level detection and select 100 successful incident detections from 25 clusters uniformly
at random. We set the time window as the duration of time from the start of detection to the end,
defined as the duration when the ratio decreased and then returned to normal. We do the samefor each of the four cities and then calculate the rating levels of each segment for each incident
window in each offending cluster. Figure 15b shows the accuracy of the proposed framework at
detecting the affected segments given a successful detection at the cluster level. Detection at thesegment level is only considered successful when it can detect both the origin and the adjacentsegments. Given these requirements and the scale of the problem, where each dataset has 100
incidents across 1000s segments, the framework can identify the affected segments >75% of thetime. The framework is lightweight and quick, and detection time is less than a second per incident,
which enables it to be placed even on resource-constrained devices.
7

CONCLUSION

We proposed an unsupervised time series-based anomaly detection framework for city-scale smart
transportation CPS. We discuss how an existing anomaly detection metric (Harmonic to Arithmetic
Mean ratios) can be applied to a transportation problem, by using a strategic partitioning of cityarea into positively correlated clusters that guarantee high invariance in detection metric. We utilize
a data augmentation technique to enable unsupervised learning of the anomaly detection techniqueand learn the bounds of the technique under sanitized, normal traffic conditions to establish anomaly
detection criteria. Results show that our proposed unsupervised anomaly detection framework
allows strategic partitions to independently generate, sanitize, learn and detect anomalies withhigh accuracy and low false-positive rates. Further, we extend our incident detection framework to
enable individual road segment detection under incident by trust score assignment. These enable
our approach to be deployed in a decentralized manner while maintaining high-performanceincident detection in a real-time manner.
ACKNOWLEDGMENTS
This research was supported by National Science Foundation, USA grants, SATC-2030611, CNS1818942, CNS-1818901, SATC-2030624, OAC-2017289. We thank the Tennessee department oftransportation and Nashville fire department for help in obtaining the relevant data for thisresearch.
REFERENCES[1] Shameek Bhattacharjee and Sajal K. Das. 2021. Detection and Forensics against Stealthy Data Falsification in Smart
Metering Infrastructure. IEEE Transactions on Dependable and Secure Computing 18, 1 (2021), 356–371. https://doi.org/10.1109/TDSC.2018.2889729[2] Shameek Bhattacharjee, Venkata Praveen Kumar Madhavarapu, Simone Silvestri, and Sajal K. Das. 2021. Attack
Context Embedded Data Driven Trust Diagnostics in Smart Metering Infrastructure. ACM Trans. Priv. Secur. 24, 2,
Article 9 (Jan. 2021), 36 pages. https://doi.org/10.1145/3426739[3] Shameek Bhattacharjee, Aditya Thakur, and Sajal K. Das. 2018. Towards Fast and Semi-Supervised Identification of
Smart Meters Launching Data Falsification Attacks. In Proceedings of the 2018 on Asia Conference on Computer andCommunications Security (Incheon, Republic of Korea) (ASIACCS ’18). Association for Computing Machinery, New
York, NY, USA, 173–185. https://doi.org/10.1145/3196494.3196551[4] Erik D. Demaine, Dotan Emanuel, Amos Fiat, and Nicole Immorlica. 2006. Correlation clustering in general weighted
graphs. Theoretical Computer Science 361, 2 (2006), 172–187. https://doi.org/10.1016/j.tcs.2006.05.008 Approximation
and Online Algorithms.
ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.

Scalable Pythagorean Mean based Incident Detection in Smart Transportation Systems

1:25[5] Keval Doshi and Yasin Yilmaz. 2021. An Efficient Approach for Anomaly Detection in Traffic Videos. arXiv preprint
arXiv:2104.09758 (2021).
[6] Min Du, Feifei Li, Guineng Zheng, and Vivek Srikumar. 2017. DeepLog: Anomaly Detection and Diagnosis from System
Logs through Deep Learning. In Proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications
Security (Dallas, Texas, USA) (CCS ’17). Association for Computing Machinery, New York, NY, USA, 1285–1298.
https://doi.org/10.1145/3133956.3134015[7] David I. Urbina Jairo Giraldo Alvaro A. Cardenas Junia Valente Mustafa Faisal. 2016. Survey and New Directions for
Physics-Based Attack Detection in Control Systems. NIST NIST GCR 16-010 (2016).
[8] Yaser E. Hawas. 2007. A fuzzy-based system for incident detection in urban street networks. Transportation ResearchPart C: Emerging Technologies 15, 2 (2007), 69–95. https://doi.org/10.1016/j.trc.2007.02.001[9] Simon Hawkins, Hongxing He, Graham Williams, and Rohan Baxter. 2002. Outlier Detection Using Replicator Neural
Networks. In Data Warehousing and Knowledge Discovery, Yahiko Kambayashi, Werner Winiwarter, and Masatoshi
Arikawa (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 170–180.
[10] Yue Hu and Daniel B. Work. 2021. Robust Tensor Recovery with Fiber Outliers for Traffic Events. ACM Trans. Knowl.
Discov. Data 15, 1, Article 6 (Dec. 2021), 27 pages. https://doi.org/10.1145/3417337[11] Zafar Iqbal and Majid Iqbal Khan. 2018. Automatic incident detection in smart city using multiple traffic flow
parameters via V2X communication. International Journal of Distributed Sensor Networks 14, 11 (2018), 1550147718815845.
https://doi.org/10.1177/1550147718815845[12] Jaminur Islam, Jose Paolo Talusan, Shameek Bhattacharjee, Francis Tiausas, Sayyed Mohsen Vazirizade, Abhishek
Dubey, Keiichi Yasumoto, and Sajal K. Das. 2022. Anomaly based Incident Detection in Large Scale Smart Transportation
Systems. In 2022 ACM/IEEE 13th International Conference on Cyber-Physical Systems (ICCPS). 215–224. https://doi.org/
10.1109/ICCPS54341.2022.00026[13] S.J. Julier and J.K. Uhlmann. 2004. Unscented filtering and nonlinear estimation. Proc. IEEE 92, 3 (2004), 401–422.
https://doi.org/10.1109/JPROC.2003.823141[14] Akira Kinoshita, Atsuhiro Takasu, and Jun Adachi. 2015. Real-time traffic incident detection using a probabilistic topic
model. Information Systems 54 (2015), 169–188. https://doi.org/10.1016/j.is.2015.07.002[15] Kwei-Herng Lai, Daochen Zha, Guanchu Wang, Junjie Xu, Yue Zhao, Devesh Kumar, Yile Chen, Purav Zumkhawaka,
Minyang Wan, Diego Martinez, and Xia Hu. 2020. TODS: An Automated Time Series Outlier Detection System. CoRR
abs/2009.09822 (2020). arXiv:2009.09822 https://arxiv.org/abs/2009.09822[16] Pankaj Malhotra, Anusha Ramakrishnan, Gaurangi Anand, Lovekesh Vig, Puneet Agarwal, and Gautam Shroff. 2016.
LSTM-based Encoder-Decoder for Multi-sensor Anomaly Detection. CoRR abs/1607.00148 (2016). arXiv:1607.00148
http://arxiv.org/abs/1607.00148[17] Ayan Mukhopadhyay, Geoffrey Pettet, Sayyed Mohsen Vazirizade, Di Lu, Alejandro Jaimes, Said El Said, Hiba Baroud,
Yevgeniy Vorobeychik, Mykel Kochenderfer, and Abhishek Dubey. 2022. A Review of Incident Prediction, Resource
Allocation, and Dispatch Models for Emergency Management. Accident Analysis & Prevention 165 (2022), 106501.
https://doi.org/10.1016/j.aap.2021.106501[18] Murat Ozbayoglu, Gokhan Kucukayan, and Erdogan Dogdu. 2016. A real-time autonomous highway accident detection
model based on big data processing and computational intelligence. In 2016 IEEE International Conference on Big Data
(Big Data). 1807–1813. https://doi.org/10.1109/BigData.2016.7840798[19] Hae-Sang Park and Chi-Hyuck Jun. 2009. A simple and fast algorithm for K-medoids clustering. Expert Systems with
Applications 36, 2, Part 2 (2009), 3336–3341. https://doi.org/10.1016/j.eswa.2008.01.039[20] Yasas Senarath, Ayan Mukhopadhyay, Sayyed Vazirizade, hemant Purohit, Saideep Nannapaneni, and Abhishek Dubey.
2021. Practitioner-Centric Approach for Early Incident Detection Using Crowdsourced Data for Emergency Services.
In 21st IEEE International Conference on Data Mining (ICDM 2021).
[21] Yasas Senarath, Saideep Nannapaneni, Hemant Purohit, and Abhishek Dubey. 2020. Emergency Incident Detection from
Crowdsourced Waze Data using Bayesian Information Fusion. In 2020 IEEE/WIC/ACM International Joint Conference on
Web Intelligence and Intelligent Agent Technology (WI-IAT). 187–194. https://doi.org/10.1109/WIIAT50758.2020.00029[22] Siemens. 2018. Connected Vehicle Roadside Unit (RSU). Brochure. https://www.mobotrex.com/download/downloaddatasheet-for-siemens-connected-vehicle-roadside-unit-rsu/[23] Hong Yang, Zhenyu Wang, Kun Xie, and Dong Dai. 2017. Use of ubiquitous probe vehicle data for identifying secondary
crashes. Transportation Research Part C: Emerging Technologies 82 (2017), 138–160. https://doi.org/10.1016/j.trc.2017.
06.016[24] Kun Zhang and Michael A.P. Taylor. 2006. Effective arterial road incident detection: A Bayesian network based algorithm.
Transportation Research Part C: Emerging Technologies 14, 6 (2006), 403–417. https://doi.org/10.1016/j.trc.2006.11.001[25] Lin Zhu, Rajesh Krishnan, Aruna Sivakumar, Fangce Guo, and John W. Polak. 2019. Traffic Monitoring and Anomaly
Detection based on Simulation of Luxembourg Road Network. In 2019 IEEE Intelligent Transportation Systems Conference(ITSC). 382–387. https://doi.org/10.1109/ITSC.2019.8917015

ACM Trans. Cyber-Phys. Syst., Vol. 1, No. 1, Article 1. Publication date: January 2023.