g-RIPS Sendai 2023
IHI Group

Final Report

Authors:

Mentors:

Melissa Adrian1

Natsuo Miyatake+

Namito Asano2

Masao Ono*

Ann Mwangi3
Tomoki Okamoto2
Aaron Scheiner4

1

The University of Chicago

2

Musashino University

3

African Institute for Mathematical Sciences

4

Northwestern University+

Academic Mentor, Tohoku University*

Industrial Mentor, IHI Corporation

Date of presentation August 8, 2023.
Finalized by industrial May 16, 2024.
mentorTitle: Mathematics for Trajectory Extrapolation Using Vehicle and Human Traffic Data
Towards Zero Traffic Fatalities
IHI document number: HAPD117219

Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 1 of 19

G-RIPS Sendai 2023 IHI Group

Contents
1

Introduction2

2

Objectives

2

3

Data Description

3

4

Methods

45

4.1

Car-Following Models ................................................................................................................4

4.2

Gipps Model ...............................................................................................................................5

4.3

Ensemble Kalman Smoother ......................................................................................................7
4.3.1

Assumptions ...................................................................................................................7

4.3.2

Algorithms .....................................................................................................................8

4.3.3

Specifying Known Quantities .........................................................................................9Results10

5.1

Objective 1: Assigning Labels to Unknown Trajectory Directions .......................................... 10

5.2

Objective 2: Reducing Noise in the LiDAR Sensors ................................................................ 11

5.3

Objective 3: Extrapolating Trajectories ................................................................................... 14

6

Conclusion

14

7

Further Directions

14

7.1

7.2

Improvements to the Forward Model ...................................................................................... 14
7.1.1

2D Car-Following Models ............................................................................................ 15

7.1.2

Joint Forward Model Learning and State Space Estimation ....................................... 15

Improvements to the Data Assimilation Framework ............................................................... 15
7.2.1

Covariance Tapering .................................................................................................... 15

7.2.2

Assimilating Temporally Intermediate Observations ................................................... 15

References
Appendix A.

16
Data Preprocessing

17

A.1

Encoding Reference Trajectories ............................................................................................... 17

A.2

Spatial Region of Interest .......................................................................................................... 18

A.3

2D to 1D Coordinate Conversion .............................................................................................. 18

1
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 2 of 19

G-RIPS Sendai 2023 IHI Group

1

Introduction

Traffic sensors are increasingly being applied to intelligent transportation systems (ITS). In particular,Light Detection And Ranging (LiDAR) can precisely capture shapes of vehicles and pedestrians as 3D
point cloud data. Owing to this capability, the LiDAR technology holds significant potential for improving
both traffic safety and sustainability to a great extent.
Traditional traffic measurements rely on video cameras for object detection. However, poor visibility ininclement weather or halation at night often deteriorates the detection accuracy. To solve this problem,
IHI has been developing more reliable sensing systems using 3D Laser RadarTM [1], which is a proprietary
LiDAR sensor of IHI. One of the strengths of these systems is high performance object tracking by detection
algorithms refined over 20 years of the development and operation of 3D Laser Radar [2][3].
The LiDAR sensors, including 3D Laser Radar, use lasers instead of radio waves to scan objects on theroad. The shorter wavelengths of lasers provide a finer resolution for observing the shapes of objects
compared to conventional radar sensors. The LiDAR sensing typically gathers point cloud data byleveraging the Time-of-Flight (ToF) principle, which measures the time taken for laser light to return. The
ToF sensing facilitates object detection even in environments imperceptible to the human eye.
In addition, the detection algorithms implemented in 3D Laser Radar systems enable to accurately track
and classify each object along its trajectory. To date, over 2,400 units of 3D Laser Radar sensors have beenmonitoring vehicles and pedestrians across the world, especially at railroad crossings [3]. Recently, thedevelopment has been extended to various ITS applications such as intersections and highways [4]. Thisadvancement helps both drivers and autonomous vehicles to drive safely as eyes on the roadside.
While the LiDAR technology is in high demand in accident-prone areas, object tracking using the ToF
sensors can be a challenge, particularly in high-traffic intersections and merging sections. For example, if
one vehicle is obscuring another from the view of a roadside sensor, the sensor cannot detect the hidden
one. This is called occlusion and often makes tracking performance more vulnerable to detection errors.
Therefore, to fully observe all trajectories of vehicles and pedestrians, we need to cover the measurement
areas with a large number of sensors and eliminate any blind spots.
However, a significant disadvantage of the LiDAR, its high implementation cost, restricts the number oflocations where the sensors can be installed. Consequently, our primary objective in this study is tominimize the requisite number of LiDAR sensors to capture accurate trajectories by harnessing
mathematics. Our team proposes a novel approach to smooth partial LiDAR data and estimate complete
trajectories by using mathematical car-following models and offline smoothing methods. This methodology
can pave the way for broader applications of LiDAR sensors, contributing to the realization of safe andsustainable transportation systems.

2

Objectives

The primary goals in this project are as follows:1.

Assign labels that represent traveling directions to vehicles labeled as ‘unknown’ directions by using
geographic and geometric information on a specified intersection as well as given LiDAR data.

2.

Reduce noises and generate probable trajectories by eliminating measurement errors caused by data
fusion of multiple LiDAR sensors.

2
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 3 of 19
3.

G-RIPS Sendai 2023 IHI Group

Design mathematical models to extrapolate the trajectories of vehicles.
If time permits, we can explore the following additional goals related to this project:

4.

Develop an alternative approach for assigning direction labels to vehicles labeled as ‘unknown’directions.

5.

Evaluate safety levels of the intersection by analyzing trajectories and labels generated by our
proposed methods.

3

Data Description

In this project, we studied mathematical models by using traffic flow datasets that were observed by four
LiDAR systems temporarily installed on the roadside. These measurements were taken in real traffic
situations (and not on specific test courses), such as multilane traffic intersections with traffic lights in
downtown Tokyo (Figure 1).

Figure 1: Traffic situations of datasets.
The measurements were conducted under the following conditions:⚫ Date: 11/24/2021 8:00 am to 5:00 pm (JST)⚫ Number of LiDAR sensors: four units➢ The data was generated by fusing all LiDAR sensors.
➢ The data contains measurement noises caused by synchronization errors at multiple time periods.
The dataset consists of time series data of each object. The columns are shown in Figure 2.
⚫ The following data are assigned to each vehicle:➢ Type: normal car, large car, bike, and pedestrian➢ Size: height [m], length[m], and width [m]⚫ The following data are recorded in each index (frame):➢ x-coordinate of vehicle’s head [m] (on local coordinate system)➢ y-coordinate of vehicle’s head [m] (on local coordinate system)

3
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 4 of 19

G-RIPS Sendai 2023 IHI Group★ Although the dataset includes pedestrian trajectories, the mathematical modeling is excluded fromour research objects mentioned in Section 2.
Figure 2: Measurement items of datasets.

4

Methods

In achieving our research objectives, we explore a variety of car-following models, especially onedimensional models mentioned in [5]. Furthermore, we use a framework of Ensemble Kalman Smoother(EnKS) [6] to integrate the car-following models and actual LiDAR data.

4.1

Car-Following Models

Car-following models are commonly classified into three categories: General Motors (GM)-type models,safety-distance models, and psychophysical car-following models. In the GM-type models, follower vehicle’s
acceleration is proportional to three factors: (1) follower’s speed, (2) difference in speed between the follower
vehicle and the leading one, and (3) spacing between the follower vehicle and the leading one. In the safetydistance models, the followers maintain their safe distance to the leading vehicles. Lastly, in thepsychophysical car-following models, the follower’s driving behavior is determined by “perceptual
thresholds.” All drivers can react to changes in spacing or relative velocity only when certain variables
exceed their own thresholds. For instance, one of the thresholds is the minimum value of relative velocity
that the followers can perceive [5][7].
The GM-type car-following models are based on the simple GM model, in which the acceleration of thefollowing vehicle is calculated by subtracting the current speed of the leader from that of the follower. The
Gazis-Herman-Rothery model incorporates the concept of distance into the simple GM model and can
describe more detailed driving behavior [5][7].
The Gipps car-following model, which is the most famous safety-distance model, assumes the following

4
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 5 of 19

G-RIPS Sendai 2023 IHI Grouptwo constraints: (1) drivers should not exceed their desired speed; (2) they should always maintain theirown safe distance [5][7][8].
The Wiedemann model, a commonly referenced psychophysical car-following model, postulates thatcharacteristics of each driver depend on the following four factors: (1) driving capabilities for perception,
reaction, and estimation of surrounding traffic environment, (2) safety requirements, (3) desired speed, and(4) aggressiveness towards the maximum acceleration and deceleration. Based on these premises, eachdriver determines subsequent driving behavior, as determined by perceptual thresholds, and transitions
between the following four regimes: (1) free driving, (2) approaching/closing, (3) following/closing, and (4)emergency [5][7].
Additionally, we also address other models and several applications.
Optimal velocity model (OVM) describes the dynamics of traffic congestion using the equation of motionof each vehicle. In the OVM, all drivers possess constant sensitivities. Each driver's velocity alters toward
the optimal velocity, which is determined by a function of the following headway, representing the distanceto the preceding vehicle. The optimal velocity, 𝑉(𝛥𝑥) , is empirically known to follow 𝑉(𝛥𝑥) =tanh(∆𝑥 − 2) + tanh⁡(2), where 𝛥𝑥 represents the following headway [9].
The car-following models described above have been utilized for driving control systems of vehicles such
as Adaptive Cruise Control (ACC) and Cooperative Adaptive Cruise Control (CACC). The ACC serves
as driver assistance systems, while the CACC is an advancement of the ACC system, incorporating vehicleto-vehicle (V2V) communication systems [5].
In addition to the car-following models, a lane-change model was proposed by Peter G. Gipps, theoriginator of the Gipps car-following model. The lane-change model entails several rules governing lanechanging behavior. Lane changes are further categorized into two categories: mandatory lane changes,
which are necessary, and discretionary lane changes, which are optional [5][7].
In this study, we opted for the Gipps model for three reasons. The foremost reason pertains to theexpressive power of the Gipps model. This model can closely emulate real-world driving behavior withrelatively simple mathematical model. Furthermore, this concept of the safe distance inherent in this model
appears to align with car-following behavior within the intersections shown in Section 3. Finally, theparameters of the Gipps model can be identified from our intersection data. From these reasons, we believethat the Gipps model, or the safety-distance model, is appropriate for our research.

4.2

Gipps Model

The Gipps model focuses on the concept of the safety distance from a following vehicle to its preceding
one. As mentioned above, this model has two constraints:(1) Drivers’ velocity should not exceed their desired velocity.
(2) Drivers should maintain their own safety distance.
These constraints appear to align with typical behavior of actual drivers following the preceding vehicles.
In this model, the velocity of the following vehicle after reacting to any behavior of the preceding one, isdescribed as follows:𝑣𝑛 (𝑡 + 𝑇) = min⁡{𝑣𝑛𝑎 (𝑡 + 𝑇), 𝑣𝑛𝑏 (𝑡 + 𝑇)},
where 𝑣𝑛𝑎 (𝑡 + 𝑇)⁡and 𝑣𝑛𝑏 (𝑡 + 𝑇) are given by5
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

(4.1)

Page 6 of 19

G-RIPS Sendai 2023 IHI Group

𝑣 (𝑡)

𝑣 (𝑡)

𝑛
𝑛
𝑣𝑛𝑎 (𝑡 + 𝑇) = 𝑣𝑛 (𝑡) + 2.5⁡𝑎𝑛max ⁡𝑇⁡ {1 − desired√0.025 + desired },

𝑣𝑛

(4.2)

𝑣𝑛

𝑣2

(𝑡)

𝑣𝑛𝑏 (𝑡 + 𝑇) = 𝑑𝑛max ⁡𝑇 + √(𝑑𝑛max ⁡𝑇)2 − 𝑑𝑛max [2{𝑣𝑛−1 (𝑡) − 𝑠𝑛−1 − 𝑥𝑛 (𝑡)} − 𝑣𝑛 (𝑡)⁡𝑇⁡ 𝑛−1
],
̂
𝑑𝑛−1

(4.3)

where
⚫

𝑇 is a reaction time of the following vehicle.

⚫

𝑛 is the follower’s index and 𝑛 − 1 is the leader’s index as shown in Figure 3.

⚫

𝑣𝑛 (𝑡)⁡ and 𝑣𝑛−1 (𝑡) are vehicles’ velocities of the follower and the leader at the time 𝑡, respectively.

⚫

𝑣𝑛desired ⁡is the follower’s maximum desired velocity.

⚫

𝑎𝑛max is the follower’s maximum acceleration.

⚫

𝑑𝑛max is the follower’s deceleration at the most severe braking.

⚫

𝑑̂𝑛−1 ⁡is the leader’s 𝑑𝑛max estimated by the follower 𝑛.

⚫

𝑥𝑛 (𝑡)⁡and 𝑥𝑛−1 (𝑡) are locations of the follower and the leader at the time 𝑡, respectively.

⚫

𝑠𝑛−1 ⁡ is inter-vehicle spacing at a stop, which is the follower’s desired inter-vehicle spacing at the stop
including the length of the leader’s vehicle.

Figure 3: Indexes of leader and follower.

Each vehicle’s velocity 𝑣𝑛 (𝑡 + 𝑇) at the time 𝑡 + 𝑇 depends on the previous value at 𝑡 + 𝑇 − 1. The
constraint (1) is considered by 𝑣𝑛𝑎 (𝑡 + 𝑇) and drivers’ velocities 𝑣𝑛 (𝑡 + 𝑇)⁡ accelerate toward their desired
velocities. In other words, 𝑣𝑛𝑎 (𝑡 + 𝑇) adjusts the accelerating behavior of the follower 𝑛. On the other hand,𝑣𝑛𝑏 (𝑡 + 𝑇) takes into account the constraint (2) and decelerates the follower 𝑛 when there is a possibility
that the vehicle will break its safe distance and approach the leading vehicle 𝑛 − 1.
The transition between 𝑣𝑛𝑎 (𝑡 + 𝑇) and 𝑣𝑛𝑏 (𝑡 + 𝑇) occurs smoothly because 𝑣𝑛𝑏 (𝑡 + 𝑇) is put into effect6
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 7 of 19

G-RIPS Sendai 2023 IHI Group

before the maximum braking is required. The following vehicles can continue to accelerate again even after
a transition from accelerating behavior to decelerating one [8].

4.3

Ensemble Kalman Smoother

Our team considers the ensemble Kalman smoother (EnKS) presented in [6] to be a useful approach foraddressing objective 2. Kalman filters can be utilized in an online setting, whereas Kalman smoothers are
better suited to be used in an offline setting. Since the setting for our denoising task is offline, Kalmansmoothing methods are more appropriate.
We propose to utilize EnKS for two particular reasons: (1) this approach can handle cases where ourforward model is a nonlinear operator, therefore avoiding the need for a linear adjoint as in the case of ageneric Kalman smoother, and (2) we would like to utilize information from future time points to inform
our estimate at each time point. In the case of the ensemble Kalman filter (EnKF), this method estimates
quantities of interest across time only using information from [0, 𝑡] for each 𝑡, does not utilize observational
data across the entire timespan and therefore leads to a less smooth fit [11].
In our setting, the forward model that we plan to pair with Kalman methodology is a nonlinear operatorof the latent state space, namely the car-following models mentioned in Section 4.1. Therefore, we need toutilize the EnKS to smooth the data using information across the entire timespan without constructing an
operator mapping backward in time.
To achieve the same goal of denoising the vehicle data, our group will investigate results obtained byusing a cubic spline and compare these results to the EnKS.

4.3.1 Assumptions
We can formalize the assumptions underlying this Kalman smoother with the following Gaussian statespace model [12]:𝑋𝑡 = 𝐹(𝑋𝑡−1 ) + 𝜉𝑡 , 𝜉𝑡 ~𝑁(0, 𝑄𝑡 ), 1 ≤ 𝑡 ≤ 𝑇,

(4.4)

𝑌𝑡 = ⁡ 𝐻𝑡 𝑋𝑡 + 𝜂𝑡 ,⁡⁡⁡𝜂𝑡 ~𝑁(0, 𝑅𝑡 ), 1 ≤ 𝑡 ≤ 𝑇,

(4.5)

𝑋0 = 𝑝0 (𝑋0 ),(4.6)

where 𝑋𝑡 ∈ ℝ𝑑𝑥𝑡 is a vector of all quantities of interest to reconstruct at time 𝑡, including vehicle locations,velocities, and accelerations at time 𝑡. 𝑌𝑡 ∈ ℝ𝑑𝑦𝑡 is a vector of all observed quantities of interest at time𝑡, which in our setting, includes positional readings from the LiDAR sensors for vehicles in the intersection
at time 𝑡. In our application setting, the dimension of the space we would like to estimate is larger thanthe dimension of our observational space (𝑑𝑥𝑡 ≫ 𝑑𝑦𝑡 ). We emphasize that the dimension of the problem isalso changing with time, denoted by the subscripts 𝑡 in the dimension notations 𝑑𝑥𝑡 ⁡ and⁡ 𝑑𝑦𝑡 .
the car-following models that estimate vehicle variables at timerefers tousing traffic and pedestrian information

at time⁡ 𝑡 − 1. 𝐻𝑡 ∈ ℝ𝑑𝑥𝑡 ×𝑑𝑦𝑡 is a linear observation operator mapping the latent quantities of interest to
the observation space. The initial distribution 𝑝0 is assumed to be known, meaning we assume that we

7
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 8 of 19

G-RIPS Sendai 2023 IHI Group

have access to a complete characterization of the latent space at time zero [12]. Lastly, we assume
independent Gaussian noise statistics for the forward model errors and the observation errors withrespective known covariance matrices 𝑄𝑡 and 𝑅𝑡 [12].

4.3.2 Algorithms
Algorithm 1 and Algorithm 2 provide pseudo-code for EnKF and EnKS, respectively. The EnKF algorithm
alternates between two steps at each time point: forecasting into the future based on the previous time
point, and correcting this forecast based on observational data. EnKS extends on this idea utilizing future
observations; first, the EnKF solution is computed across time points 1 ≤ 𝑡 ≤ 𝑇, then this solution is
improved utilizing observations up to 𝐿 time points in the future.
Ensemble Kalman Filter (EnKF) [12]Algorithm 1.
1:

Input: 𝑌1:𝑇 , 𝑋01:𝑁 . (If 𝑋01:𝑁 is not specified, draw 𝑋0𝑛 ~i.i.d. 𝑝0 (𝑋0 ))2:for 𝑡 = 1, … , 𝑇⁡do

3:𝑛 )Set 𝑋̂𝑡𝑛 = 𝐹(𝑋𝑡−1+ 𝜉𝑡𝑛 , where⁡⁡𝜉𝑡𝑛 ~i.i.d. 𝑁(0, 𝑄𝑡 ).

▷ Forecast−14:̂𝑡 = 𝐶̂𝑡 𝐻𝑡T (𝐻𝑡 𝐶̂𝑡 𝐻𝑡T + 𝑅𝑡 ) .
Compute 𝑚
̂ 𝑡 and 𝐶̂𝑡 using (4.7) and (4.8), and set 𝐾

5:

̂𝑡 (𝑌𝑡 + 𝜂𝑡𝑛 − 𝐻𝑡 𝑋̂𝑡𝑛 ), where 𝜂𝑡𝑛 ~i.i.d. 𝑁(0, 𝑅𝑡 ).
Set 𝑋𝑡𝑛 = 𝑋̂𝑡𝑛 + 𝐾

6:

end for7:

1:𝑁
Output: 𝑋1:𝑇

▷ Analysis

𝑁

𝑚
̂𝑡 =

1
∑ 𝑋̂𝑡𝑛
𝑁

(4.7)

𝑛=1

𝑁

1
𝐶̂𝑡 = ⁡∑(𝑋̂𝑡𝑛 − 𝑚̂ 𝑡 )(𝑋̂𝑡𝑛 − 𝑚̂ 𝑡 )T
𝑁−1

(4.8)

𝑛=1

Ensemble Kalman Smoother (EnKS) [6]Algorithm 2.

𝑌1:𝑇 ,⁡ 𝑋01:𝑁 . (If

𝑋01:𝑁 is not specified, draw 𝑋0𝑛 ~i.i.d. 𝑝0 (𝑋0 ))1:

Input:

2:

1:𝑁
Compute 𝑋̂1:𝑇
using Algorithm 1.

3:

for 𝑡 = 1, … , 𝑇⁡do

▷ EnKF

4:

Compute 𝑚
̂ 𝑡 and 𝐶̂𝑡 using (4.7) and (4.8).

5:TT
̂𝑡,𝑡+𝑙 = 𝐶̂𝑡,𝑡+𝑙 𝐻𝑡+𝑙
Compute 𝐾(𝐻𝑡+𝑙 𝐶̂𝑡,𝑡+𝑙 𝐻𝑡+𝑙
+ 𝑅𝑡+𝑙 ) .

6:

𝑛
𝑛
𝑛
̂𝑡,𝑡+𝑙 (𝑌𝑡+𝑙 + 𝜂𝑡+𝑙
Set 𝑋𝑡𝑛 = 𝑋̂𝑡𝑛 + ∑𝐿𝑙=1 𝐾− 𝐻𝑡+𝑙 𝑋̂𝑡+𝑙), where⁡ 𝜂𝑡+𝑙
~i.i.d. 𝑁(0, 𝑅𝑡+𝑙 ).

−1

7:

end for8:

1:𝑁
Output: 𝑋1:𝑇

▷Look ahead
time steps

8
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 9 of 19

G-RIPS Sendai 2023 IHI Group

4.3.3 Specifying Known Quantities
One main challenge in applications of EnKF and EnKS is that some quantities are assumed to be knowna priori, but in practice, these quantities are potentially not known by the researchers. Particularly, theerror covariances 𝑄𝑡 and 𝑅𝑡 are assumed to be known, as well as the parameterization of the forward
model 𝐹 and the complete initial condition⁡𝑋0 . In practice, these quantities are estimated prior to the
EnKF/EnKS analysis.
In this work, we do not have a prior sense of the magnitudes of 𝑄𝑡 and 𝑅𝑡 , and the forward model
parameters. We can fully specify the initial condition using our observation data.
Dimensionality of the latent and observation states.
The dimensionality of the system in terms of latent states and observations is dynamic in this trafficsystem since vehicles are entering and exiting our region of interest. The latent state at time 𝑡 is comprisedof all vehicles’ positions, speeds, and accelerations in the intersection. This vector contains all quantities
that we would like to estimate at the given time 𝑡. The number of quantities we would like to estimate attime 𝑡 is denoted by 𝑑𝑥𝑡 . The observations at time

are the vehicles’ positions detected by the LiDAR

sensor. The number of vehicles’ positions we detect at time 𝑡 is denoted by 𝑑𝑦𝑡 . We note in this setting,𝑑𝑥𝑡 >>𝑑𝑦𝑡 .
Choosing the ensemble size.
Since the ensemble size is directly tied to how well 𝐶̂𝑡 in (4.8) is estimated, we would like to choose thisuser-defined hyper-parameter as large as possible. However, setting this hyper-parameter to a considerable
number increases the storage and computational space needed for the problem. This trade-off is problem
specific and depends on the size of 𝑑𝑥𝑡 . For this project, we set 𝑁 = 100 as a starting point to assess theperformance.
Creating the initial ensemble.
To create the initial ensemble, we choose to initialize the system using observational data that corresponds
to the first time point in the dataset that we have both velocity and acceleration information, whichcorresponds to the first timepoint in the dataset where we have both velocity and acceleration information
for all observed vehicles in the dataset. We denote this initial condition as 𝑋0 .
To create the initial ensemble, we add Gaussian noise to the initial condition in a way that describes our
uncertainty of this initial condition. More formally,
𝑋0𝑛 = ⁡ 𝑋0 + 𝜀 𝑛 ,⁡⁡⁡𝜀 𝑛 ~𝑁(0, 𝐵),⁡

(4.9)

where 𝐵⁡ is specified a priori.
Forward model.
The forward model is a user-chosen model that maps the latent states 𝑋𝑡 to a prediction at time⁡ 𝑡 + 𝑇.
There are a variety of different choices for forward models for our system, and our choice of forward modelis the Gipps model described in Section 4.2. Since this model describes traffic behavior for a 1D systemalong a straight line, some preprocessing steps were taken to create a dataset compatible with this forward
model. Appendix A will describe the preprocessing steps that were applied to the dataset. We will mention
alternative car-following models in Section 7.1.1 that describe 2D traffic models, which may be more

9
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 10 of 19

G-RIPS Sendai 2023 IHI Group

suitable for our dataset.
Forward model error covariance.
We assume that the forward model error covariance in equation (4.4), 𝑄𝑡 , is a matrix with nonzero
diagonal variances and nonzero covariances between vehicles that are spatially close to each other.

Observation noise error covariance.
We assume that the observation error covariance in equation (4.5), 𝑅𝑡 , is a diagonal matrix withdimension 𝑑𝑦𝑡 × 𝑑𝑦𝑡 . More specifically, in this work, we assume this matrix has the form 𝑅𝑡 = 𝑟𝐼𝑑𝑦 , where𝑡𝑟 ∈ ℝ . This choice of matrix assumes that the LiDAR sensor’s measurement errors for different objects
are independent (based on the diagonal specification), and the errors across vehicles are distributed withthe same magnitude (based on the scalar 𝑟 specification).
Observation operator.
In our setting, this matrix is in {0,1}𝑑𝑦𝑡 ×𝑑𝑥𝑡 . This matrix has the following form:
1⁡⁡if⁡𝑌𝑖,𝑡 ⁡is⁡the⁡observation⁡of⁡quantity⁡𝑋𝑗,𝑡 ,𝐻𝑖𝑗,𝑡 = {⁡0⁡⁡otherwise,(4.10)

where 𝑖 = 1, … , 𝑑𝑦𝑡 ⁡and 𝑗 = 1, … , 𝑑𝑥𝑡 . For example, if 𝑌𝑖,𝑡 is the observation for a certain vehicle’s position
and 𝑋𝑗,𝑡 is the true position for that same vehicle, then 𝐻𝑖𝑗,𝑡 = 1. On the other hand, if 𝑋𝑗,𝑡 were thevelocity for that same vehicle, 𝐻𝑖𝑗,𝑡 would be zero since 𝑌𝑖,𝑡 would not be a measurement of the quantity
𝑋𝑗,𝑡 .

5
5.1

Results

Objective 1: Assigning Labels to Unknown Trajectory Directions

The procedure for assigning labels is as follows:(1) We extracted data with unknown trajectory directions from the given dataset. The extracted datacontains the following information of each vehicle: ID, labels of traveling directions (from, to), the firstmeasured time, the last measured time, initial 𝑥, initial 𝑦, last 𝑥, last 𝑦, initial velocity, and initialacceleration.
(2) We quantified the similarities between the observed trajectory and all reference trajectories by
calculating mean squared error (MSE) described below:
𝑁

MSE𝑗

1
= ∑{𝑦𝑖 − 𝑓𝑗 (𝑥𝑖 )}2 ,⁡
𝑁
𝑖

where
⚫

𝑦𝑖 ⁡and 𝑥𝑖 are the 𝑖th observed coordinate of 𝑦⁡and 𝑥, respectively.

⚫

𝑓𝑗 (𝑥𝑖 )⁡is the y-coordinate of the 𝑗th reference trajectory on the 𝑖th observed x-coordinate⁡𝑥𝑖 .

10
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

(5.1)

Page 11 of 19

G-RIPS Sendai 2023 IHI Group

Please refer to Appendix A for details on 𝑓.
⚫

is the number of observed points (𝑥𝑖 , 𝑦𝑖 ) of any objects. (𝑥1 , 𝑦1 ) and (𝑥𝑁 , 𝑦𝑁 ) are the initial point
and the last point, respectively.

(3) Figure 4 shows an example of label assignment. Using a vehicle trajectory (orange line) starting from
"East" and ending at "Unknown", we replaced “Unknown” with the direction label of a reference trajectory.
The orange line represents a partial trajectory of the vehicle, which was lost from the LiDAR sensors’ view
before exiting the intersection. After calculating MSE𝑗 for all reference trajectories, we assigned the vehicle
as the label “North” of the reference trajectory with the minimum MSE𝑗 (blue line).

Figure 4: An example of label assignments.

5.2

Objective 2: Reducing Noise in the LiDAR Sensors

To reduce noise of the LiDAR sensor data, we can utilize the EnKS framework. Results shown in thissection contain estimated positions, velocities, and accelerations over time for one vehicle. The estimates
for this vehicle were computed independently from any other vehicles in the dataset, which does not accountfor information about surrounding cars that would aid in estimating this vehicle’s information. We note
this limitation in mind when interpreting our results.
The results show EnKF, EnKS, and Gipps forecasting compared to the observational LiDAR data. In
these results, for EnKS, we set 𝐿 = 7, and in both EnKF and EnKS, we set 𝑅 = 𝑟𝐼, where 𝑟 = 0.05, and𝑄 = 𝑞𝐼, where 𝑞 = 0.5, and 𝑁 = 100. The parameters of the Gipps model in this example are shown in
Table 1.

11
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 12 of 19

G-RIPS Sendai 2023 IHI Group

Table 1: The parameter values in this study [10]Parameter

Value

Unit

𝑣𝑛desired

15.3

m/s

𝑎𝑛max

1.77

m/s2

|𝑑𝑛max |

3.51

m/s2

|𝑑̂𝑛−1 |

6.3

m/s2

𝑠𝑛−1

5.9

m

𝑇

0.6

s

The estimated time series of positions, velocities, and accelerations are shown in Figure 5, Figure 6, andFigure 7, respectively. Figure 5 and Figure 6 demonstrate that Gipps forecast and EnKF/EnKS analysis
can be assimilated into the trend of positions and velocities. As shown in (4.1)-(4.3), the forward model has
no information on second or higher derivative. Thus, the acceleration cannot be estimated, and the noiseis amplified compared to the complete observation.
Meanwhile, Figure 6 indicates that the forward model highly contaminated the smoothed velocities due
to its large uncertainty. Based on the current specification of the Gipps model, specifying 𝑞 to be large
and 𝑟 to be small would be most optimal since we see that the Gipps forecasts are unreliable. If we wereto change 𝑞 and⁡𝑟 in this way, we would receive EnKS and EnKF results that closely match theobservations. To obtain results that are less noisy compared to the observational data, we need a forward
model whose errors are on a smaller scale compared to the observations. With this in mind, we suggest
future directions for consideration for improving the forward model in Section 7.1.

Figure 5: Comparison of time series of positions.

12
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 13 of 19

G-RIPS Sendai 2023 IHI Group

Figure 6: Comparison of time series of velocities.

Figure 7: Comparison of time series of accelerations.

13
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 14 of 19

G-RIPS Sendai 2023 IHI Group

5.3

Objective 3: Extrapolating Trajectories

We present pseudo code below in Algorithm 3 for how EnKS can extrapolate vehicle trajectories that are
hidden in the dataset, and how to account for the dynamic nature of the latent space and observation
space dimensions. We provided code to IHI related to the steps of Algorithm 3, however, we note that wehave not integrated all steps of this algorithm in our code.
Extrapolating Trajectories

Algorithm 3.
𝑋01:𝑁 .

1:

Initialize

2:

for 𝑡 = 1, … , 𝑇⁡do3:

1:𝑁
Apply the Gipps forward model from equation (4.1)-(4.3) to 𝑋𝑡−1.

4:

Check if the latent dimension is correct.
(remove vehicles no longer in our region of interest, unless their observed position is still in theregion, and add vehicles that enter the region)

5:

Construct a new 𝐻𝑡 based on which latent states are observed versus unobserved including
hidden vehicles.

6:

Construct 𝑅𝑡 and 𝑄𝑡 based on 𝑑𝑦𝑡 and 𝑑𝑦𝑡 , and apply the EnKS from Algorithm 2.

7:

end for8:

1:𝑁
Output: 𝑋1:𝑇
, which includes estimated speeds, positions, and accelerations for both hidden

and observed vehicles.

6

Conclusion

We have implemented a method for assigning labels to trajectories with “unknown” travel directions. Thisapproach applies to our data preprocessing steps, which includes mapping the data from a 2D format to aset of 1D trajectories. Regarding our EnKS results, we acknowledge that work still needs to be done in
optimizing the user-chosen parameters, as well as specifying the best setup for EnKS. We provide moredetails on further directions in Section 7.

7

Further Directions

This preliminary attempt to denoise observed traffic trajectories and extrapolate missing trajectories canbe extended to account for more realistic assumptions for an intersection. For example, this approach canbe extended to include a 2D car-following model as the forward model, which would capture information
about vehicle lane changes in the intersection. Additionally, we suggest methods for jointly learning theforward model parameters and estimating the state space quantities of interest to further improve theresults.

7.1

Improvements to the Forward Model

14
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 15 of 19

G-RIPS Sendai 2023 IHI Group

Our team sees two avenues to improve upon the EnKS results regarding the specification of the forward
model, the first of which would be to consider vehicle prediction models that can produce predictions in a
2D space, and the second of which would be to consider implementing an approach to jointly learn or tune
the forward model dynamics while estimating the latent state space. We consider these two potential
avenues to be the most fruitful in improving the results.

7.1.1 2D Car-Following Models
There is less research literature on 2D car-following models compared to 1D models. However, newmodeling methods have been discussed in recent years, and further expansion is expected in the future [13].

7.1.2 Joint Forward Model Learning and State Space Estimation
This framework can be further improved to utilize the observational data to improve the forward model
parameterization. One such method is the auto-differentiable ensemble Kalman filter [14], which alternates
between estimating the latent states 𝑋𝑡 given a set of parameters and updating this parameter set by
autodifferentiating the loss function and taking a gradient step to minimize the loss function. This method
illustrated superior performance compared to an expectation-maximization approach to learning theforward model parameters.

7.2

Improvements to the Data Assimilation Framework

We see a few areas of improvement that could be made regarding our application of EnKS, specificallyrelated to including covariance tapering and assimilating temporally intermediate observations. We
consider these two directions to be important to consider improving upon the current results.

7.2.1 Covariance Tapering
Covariance tapering is a standard technique in applications of data assimilation where the latent states
represent some sort of spatial information, such as in our case of tracking information about individual carsin an intersection. In our setting, we can define a tapering matrix 𝜌𝑡 ∈ {0,1}𝑑𝑥𝑡 ×𝑑𝑥𝑡 that describes which
values of 𝐶̂𝑡 that we would like to estimate [15]. Covariance tapering is applied in the EnKF and EnKS
formula by replacing 𝐶̂𝑡 with 𝜌𝑡 ∘ 𝐶̂𝑡 , where denotes the Hadamard product [15]. The goal of this matrix
is to zero out correlations between states that we would expect to be zero, such as two vehicles that are
spatially far apart or on different roads, which increases our effective ensemble size and leads to better
approximations [15].

7.2.2 Assimilating Temporally Intermediate Observations
We note in this work that we discard observations that are temporally in between the forward model time
jump. For example, the forward model maps the latent states from time

to time 𝑡 + 𝑇 . Observations

from time points⁡ 𝑡 and⁡ 𝑡 + 𝑇 are assimilated, but the intermediate observations at 𝑡 +

𝑇
6

, 𝑡+

𝑇
3

, etc. are

not assimilated. A further extension of this work could be to extend the current implementation to be able
to assimilate these temporally intermediate observations.

15
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 16 of 19

G-RIPS Sendai 2023 IHI Group

References[1]

Realizing a sustainable mobility society with traffic flow data. IHI Engineering Review, 56(1), 2023.
https://www.ihi.co.jp/en/technology/techinfo/contents_no/1199408_13586.html[2]

3D Laser Radar. IHI Corporation.

https://www.ihi.co.jp/3DLaserRadar/en/index.html.
[3]

Level crossing obstacle detection system. IHI Corporation.
https://www.ihi.co.jp/3DLaserRadar/en/products/01.html.

[4]Intelligent Transport Systems(ITS) Equipment. IHI Corporation.
https://www.ihi.co.jp/3DLaserRadar/en/products/02.html.

[5]

Hafiz Usman Ahmed, Ying Huang, and Pan Lu. A review of car-following models and modeling tools
for human and autonomous-ready driving behaviors in micro-simulation. Smart Cities, 4(1):314–335,
2021.

[6]

Geir Evensen, Peter Jan, and Van Leeuwen. An ensemble kalman smoother for nonlinear dynamics.
Monthly Weather Review, 128:1852–1864, 2000.

[7]

Johan Janson Olstam and Andreas Tapani. Comparison of Car-following models. Swedish National
Road and Transport Research Institute: Linkoping, Sweden, 2004.

[8]

P.G. Gipps. A behavioural car-following model for computer simulation. Transportation ResearchPart B: Methodological, 15:105–111, 1981.

[9]The Mathematical Society of Traffic Flow. Optimal velocity model.
http://traffic.phys.cs.is. nagoya-u.ac.jp/~mstf/sample/ov_e.html.

[10] Luis Vasconcelos, Ana B. Silva, Alvaro M. Seco, Paulo Fernandes, and Margarida C. Coelho.
Turboroundabouts: Multicriterion assessment of intersection capacity, safety, and emissions.
Transportation Research Record, 2402(1):28–37, 2014.
[11] Roger R Labbe Jr. Kalman and Bayesian Filters in Python. 2020.
[12] Geir Evensen. The Ensemble Kalman Filter: Theoretical formulation and practical implementation.
Ocean Dynamics, 53(4):343–367, 2003.
[13] Jing Zhao, Victor L Knoop, and Meng Wang. Two-dimensional vehicular movement modelling at
intersections based on optimal control. Transportation Research Part B: Methodological, 138:1–22,
2020.
[14] Yuming Chen, Daniel Sanz-Alonso, and Rebecca Willett. Auto-differentiable ensemble kalman filters,

16
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 17 of 19

G-RIPS Sendai 2023 IHI Group2021.
[15] Thomas M. Hamill, Jeffrey S. Whitaker, and Chris Snyder. Distance-dependent filtering of backgrounderror covariance estimates in an ensemble kalman filter. Monthly Weather Review, 129(11):2776 –
2790, 2001.
[16] T.D. Barfoot and C.M. Clark. Motion planning for formations of mobile robots. Robotics andAutonomous Systems, 46(2):65–78, 2004.

Appendix A.

Data Preprocessing

Because the car-following model of interest operates only in scenarios where the road is a straight pathwith only one lane, we needed to preprocess our 2D dataset so that it is compatible with the methodology
we would like to use. This preprocessing does throw away useful traffic information, such as lane changes
between vehicles, but we noted this limitation in this document and suggest remedies for this limitation in
Section 7.

A.1 Encoding Reference Trajectories
We experimented with linear interpolation, cubic spline, and quadratic spline techniques to produce
suitable smooth curves representing reference trajectories. Figure 8 shows an example of interpolated resultsof the reference trajectory shown in Figure 4.

y [m]

x [m]

17
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 18 of 19

G-RIPS Sendai 2023 IHI Group
Figure 8: An example of label assignments.

A.2

Spatial Region of Interest

While we have traffic data for each of the four directions prior to the intersection, we are primarily
concerned with the behavior of the cars in the intersection. Thus, we restrict our data to the points
contained in the box representing the intersection, in which each edge is a crosswalk. The below image asshown in Figure 9 demonstrates the lanes on which the cars can travel in the intersection (in blue) as wellas the restricted box (in orange).

y [m]

x [m]
Figure 9: Spatial region of interest in this project.

A.3 2D to 1D Coordinate Conversion
Our smoothing and project goals necessitate one-dimensional coordinates as inputs to the forward model.
Thus, we convert from Cartesian coordinates to curvilinear coordinates. However, we specifically choose
the length along the reference trajectory, 𝑠, as our one-dimensional coordinate of interest. We initially usedthe curvilinear conversion seen in [16]. This entailed the conversion from Cartesian to curvilinear
coordinates in which:
𝑠=

1
𝜅𝑏 𝑥
arctan()
𝜅𝑏
1 − 𝜅𝑏 𝑦

,

(A.1)

1
𝑥
−
𝜅𝑏 sin(𝜅𝑏 𝑠)

.

(A.2)

and
𝑛=

18
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.

Page 19 of 19

G-RIPS Sendai 2023 IHI Group

This resulted in⁡(𝑠,⁡ 𝑛) coordinates that closely resembled (𝑥, 𝑦) coordinates. We expected 𝑠 to be thedistance along the reference trajectory, and we expected 𝑛 to be the closest distance from a point to thereference trajectory. This would imply only positive 𝑠, for example. However, our 𝑠 values and 𝑥 values
were equivalent, and the same was true for our 𝑛 and 𝑦 values. Even before this issue, we intended toonly consider the 𝑠 coordinate. In the intersection box, there are rarely lane changes, so 𝑛 can be
discounted without significant impact.
However, our s results were not appearing as we had expected them to, so we found an alternative way
to calculate 𝑠. For a given Cartesian car trajectory coordinate, we find the closest point 𝛼 along thereference trajectory and calculate the distance along the reference trajectory from its start to 𝛼. As wepreviously discussed encoding the reference trajectories, we measure distance along the reference trajectory
via summing linear distances between discretely sampled points along the reference trajectory. We use thismethod to generate 𝑠 coordinates.

19
Copyright © 2024 IHI Corporation and Tohoku University All Rights Reserved.