# Quantifying Population Movement Using a Novel Implementation of Digital
Image Correlation in the ICvectorfields package

Devin W. Goodsman

## Abstract

Movements in imagery captivate the human eye and imagination. They are also of
interest in variety of scientific disciplines that study spatiotemporal
dynamics. Popular methods for quantifying movement in imagery include particle
image velocimetry and digital image correlation. Both methods are widely
applied in engineering and materials science, but less applied in other
disciplines. This paper describes an implementation of a basic digital image
correlation algorithm in R [R Core Team, 2021] as well as an extension
designed to quantify persistent movement velocities in sequences of three or
more images. Algorithms are applied in the novel arena of landscape ecology to
quantify population movement and to produce vector fields for easy
visualization of complex movement patterns across space. Functions to
facilitate analyses are available in ICvectorfields [Goodsman, 2021]. These
methods and functions are likely to produce novel insights in theoretical and
landscape ecology because they facilitate visualization and comparison of
theoretical and observed data in complex and heterogeneous environments.

## Introduction

Living organisms move through space in complex ways that have inspired many
branches of spatial pattern analysis from Turing instabilities [Alonso et al.,
2002, Ruan, 1998], to complex systems analysis of the emergent properties of
individual-level behaviours when organisms live in groups [Parrish and
Edelstein-Keshet, 1999, Johnson, 2009]. Moreover, in mathematical ecology
there is a long history of deriving analytic expressions for traveling wave
speeds from mathematical models of biological systems based on partial
differential equations [Kolmogorov, 1937, Skellam, 1951] and integrodifference
equations [Kot et al., 1996]. In addition to standard traveling waves and
wave-trains, simulation studies have revealed more unusual patterns of
population level movement can arise from the way organism interact and move on
the landscape [Hassell et al., 1994]. Spiral waves are one example of
surprising spatiotemporal dynamics that can arise in biological systems
[Hassell et al., 1991].

Travelling waves and spiral waves emerge from mathematical models of
population expansion, which are often based on partial differential equations
and integrodifference equations. These types of models, which represent the
movement patterns of populations of organisms, are sometimes classified as
Eulerian approaches to distinguish them from Lagrangian approaches that focus
on the trajectories of individuals. The majority of R packages that quantify
organismal movement, however, are Lagrangian as they pertain to the analysis
of the tracks or trajectories of individual animals with tracking collars or
tracking devices. Integrated step selection models [Avgar et al., 2016], such
as those in the amt package [Signer et al., 2019], which incorporate the
impact of spatially variable habitat or environmental variables on movement of
individuals modeled using a discrete time and discrete space framework, are an
example of a Lagrangian approach when fitted to movement data from
individuals. Because my focus in this work is on population-level movements
that are evident in imagery, I will forego further discussion of Lagrangian
models and instead refer the interested reader to a review of R packages for
modeling animal movement [Joo et al., 2020]. At the time of writing, R
packages that focus on the Eulerian approach include IDE [Zammit-Mangion,
2019], deSolve [Soetaert et al., 2010], and ReacTran [Soetaert and Meysman,
2012]. These packages are designed primarily to obtain numerical solutions to
Eulerian models, analyze their dynamics, and fit them to data. Recently,
movement modeling based on stochastic differential equations, stochastic
partial differential equations [Krainski et al., 2018], and other stochastic
process models [Buderman et al., 2016] has proliferated. Computationally
efficient Bayesian statistical approaches are often required to fit these
stochastic models to data due to the ubiquity of noise in spatiotemporal time
series in combination with nonlinear dynamical processes [Krainski et al.,
2018].

In contrast to the R packages and approaches I have cited above, this work is
focused on the description of an empirical method for quantifying spatially
heterogeneous rates of movement in sequences of images without fitting a
model–although I do approach the quantification problem from an Eulerian
perspective. Empirically quantifying spread rates without imposing a specific
mathematical model allows the user to abandon many of the assumptions implicit
in mathematical modeling of population expansion. For example, tractable
mathematical models of consumer-resource systems that generate traveling
waves, wave-trains, and spiral waves, often rely on assumptions of a
homogeneous spatial environment with respect to resources or one in which
there are no discontinuities in resources. In contrast, many organisms spread
in environments that are spatially heterogeneous [Urban et al., 2008], and in
environments subject to persistent directional flows that impact organism
movement [Hoffman et al., 2006]. The ramifications of this claim are more
easily understood using a meteorological analogy. In meteorology, vector
fields are frequently used to illustrate the impacts of high and low pressure
systems on wind speed and direction, and thus on the movement of weather
systems. In such meteorological systems, wind speed and direction are
complicated functions of topography and complex atmospheric dynamics. As a
result, vector fields representing movement in such systems are often variable
at the regional scale, with winds flowing in one direction on one side of a
map and possibly in an opposing direction on the other side. In ecology,
populations of organisms are like the weather systems in that their movement
on the landscape is what is of primary interest to researchers; variable wind
causes spatially variable movement of weather systems similarly to how
persistent directional fluid flows, including wind, in an organism’s
environment impact dispersal, and therefore population movement.

At the time of writing, the only tool in R [R Core Team, 2021] designed to
empirically estimate spreading speed or the speed of wave-trains in
populations without fitting a mathematical model is implemented in the ncf R
package [Bjornstad, 2020]. The ncf package relies on lagged non-parametric
cross-correlation functions to estimate spreading speed of traveling waves
[Bjørnstad and Falck, 2001]. To do so, it takes two spatiotemporal data sets
that differ from one another only in that one is a time-lagged version of the
other, and projects their planar coordinates onto lines of varying angles that
can be specified using function arguments. After projection onto a line,
cross-correlation is estimated using a spline-correlogram approach [Bjørnstad
and Falck, 2001] and the location of maximum cross-correlation gives an
estimate of displacement along the direction of the projection line. This
approach was used to estimate the velocity of traveling wave-trains in the
larch budmoth system in the European Alps [Bjørnstad et al., 2002].

Projecting population data from a domain with two spatial coordinates onto a
domain that has only one spatial coordinate and then using a correlogram
approach precludes quantification of more complicated patterns of movement on
the landscape. For example if two spatially separated populations are moving
towards one another at the same speed, such an approach will yield enigmatic
correlograms. Similarly, if several populations move radially around a central
fulcrum, the correlogram will be difficult to translate to an estimate of
directional movement.

In this paper I present an approach for estimating vector-fields in systems
with spatially variable movement that is inspired by a technique from
engineering and materials science called Digital Image Correlation or DIC
[Anuta, 1970, Sutton et al., 2009]. Among other things, Digital Image
Correlation is used to estimate displacement based on photographs of a planar
material before, during, and after a force has been applied to warp its
surface [Sutton et al., 2009]. A typical DIC approach as well as the
extensions described in this paper are implemented in the ICvectorfields
package [Goodsman, 2021], in which the IC is the abbreviation for Image
Correlation. I demonstrate these approaches using the ICvectorfieds package to
analyze a simulated data set as well as the larch budmoth data set provided
with the ncf R package [Bjornstad, 2020, Bjørnstad et al., 2002].

## Mathematical and Computational Details

Here I provide mathematical and computational details of the algorithms used
in the ICvectorfields R package starting with a standard digital image
correlation approach, and following with extensions to estimate persistent
movement and to account for spatial variability in persistent movement. The
ICvectorfields package capitalizes on the algorithms written in C under the
title FFTW which stands for Fastest Fourier Transform in the West [Frigo and
Johnson, 2005], and a convenient wrapper package in R called fftwtools [Rahim,
2021]. Input raster images and raster stacks are read and manipulated using
the terra package [Hijmans, 2021].

### Digital image correlation

One of the earliest applications of cross-correlation in image analysis was to
align images taken from different sensors or at different times using
satellites or aircraft [Anuta, 1970]. The theoretical and computational
details I present here loosely follow those in this pioneering application. I
will provide the mathematical underpinning of two-dimensional cross-
correlation, and then elaborate on its computational implementation, which
involves some additional complexity due to the circular nature of discrete
fast Fourier transforms. In all descriptions below, I do not normalize the
cross-correlation function to obtain Pearson correlation coefficients and
therefore, I follow the convention of using the terms cross-correlation and
cross-covariance interchangeably.

Given two images that have been converted to square matrices $f$ and $g$ of
dimension $m\times m$, two-dimensional cross-correlation can be defined in
terms of a convolution:

$\left(f\star
g\right)\left(x_{j},y_{i}\right)=\left(\overline{f\left(-x_{j},-y_{i}\right)}*g\left(x_{j},y_{i}\right)\right)\left(x_{j},y_{i}\right),$
(1)

in which $\left(f\star g\right)$ is the two-dimensional cross-correlation
matrix, the $*$ operator denotes convolution,
$\overline{f\left(-x_{j},-y_{i}\right)}$ is the complex conjugate of the
$f\left(x_{j},y_{i}\right)$ matrix, $i$ is the matrix row index, and $j$ is
the matrix column index $i,j\in\mathbb{N}=\\{1,2,...\\}$. Note that I use
array indices that start at one rather than zero. The coordinates of the
centroids of each pixel are given by $x_{j}$ and $y_{i}$.

Based on the convolution theorem, equation (1) can be rewritten as

$\left(f\star
g\right)\left(x_{j},y_{i}\right)=\mathbb{F}^{-1}\left(\overline{\mathbb{F}\left(f\left(x_{j},y_{i}\right)\right)}\mathbb{F}\left(g\left(x_{j},y_{i}\right)\right)\right)\left(x_{j},y_{i}\right),$
(2)

wherein $\mathbb{F}$ denotes the two-dimensional Fourier transform,
$\mathbb{F}^{-1}$ denotes its inverse, and
$\overline{\mathbb{F}\left(f\left(x_{j},y_{i}\right)\right)}$ is the complex
conjugate of $\mathbb{F}\left(f\left(x_{j},y_{i}\right)\right)$. Because
$\overline{\mathbb{F}\left(f\left(x_{j},y_{i}\right)\right)}=\mathbb{F}\left(\overline{f\left(-x_{j},-y_{i}\right)}\right)$,
and because $f\left(x_{j},y_{i}\right)$ contains only real numbers, I forego
calculating the complex conjugate of the Fourier transform of the $f$ matrix
and choose instead capitalize on matrix multiplication to compute the complex
conjugate:

$\overline{f\left(-x_{j},-y_{i}\right)}=r\times f\times r,$ (3)

in which the r matrix is a $m\times m$ matrix that has zeros everywhere except
for along the diagonal that runs from its lower left to upper right corners,
which contains ones:

$r=\begin{pmatrix}0&0&\cdots&0&1\\\ 0&0&\cdots&1&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&1&\cdots&0&0\\\
1&0&\cdots&0&0\end{pmatrix}.$ (4)

The matrix calculations in equations (3) and (4) substitute for calculation of
the complex conjugate only in the case where all values of the $f$ matrix are
real, which is the case in most natural science applications.

Together, equations (1) through (4) constitute an elegant way to compute two
dimensional cross-correlation. Computer implementation of these, however,
requires some additional complexity due to the use of discrete fast Fourier
transforms to efficiently compute convolutions. Discrete fast Fourier
transforms are inherently circular, which means that what happens on the outer
edges of matrices will impact their discrete fast Fourier transform on the
opposite side. In order to mitigate this problem, zeros are added to the outer
edge on all sides of both the $f$ and $g$ matrices [Anuta, 1970]. In
ICvectorfields, the $f$ and $g$ matrices are padded with as many zeros as
there are rows and columns in the original matrix and then additional zeros
are added to ensure that both matrices are square with an even number of rows
and columns.

When matrices are padded as described above and discrete fast Fourier
transforms are used as in equation (2), the non-cyclic components of the
convolution end up in the outer corners of the cross-correlation matrix
$\left(f\star g\right)\left(x_{j},y_{i}\right)$ [Anuta, 1970]. Thus, to obtain
a correct estimate of cross-correlation, $\left(f\star
g\right)\left(x_{j},y_{i}\right)$ must be divided into four quadrants and each
quadrant must be flipped along its horizontal and vertical axes using matrix
multiplication. For example, if the zero-padded matrices have dimension
$k\times k$, each quadrant of $\left(f\star g\right)\left(x_{j},y_{i}\right)$
will have dimension $k/2\times k/2$ and the following operation flips each
quadrant matrix:

$q_{f}=r\times q\times r,$ (5)

where $r$ is a $k/2\times k/2$ matrix as in equation (4). Then the four
quadrants can be reassembled into the $k\times k$ cross-correlation matrix.
The mathematical operations in equations (1) through (5) are implemented in
the _Xcov2D_ function in the ICvectorfields package.

Once a cross-correlation matrix has been estimated, it can be used to compute
displacement in the horizontal and vertical directions in terms of the
horizontal and vertical shifts in pixel numbers that maximize cross-
correlation. In ICvectorfields, shifts to the right and up are designated as
positive, whereas shifts to the left and down are considered negative.

A typical implementation of DIC will define a region of interest within the
input images or their corresponding matrices wherein displacement vectors are
sought using a bounding box as in the _DispFieldbb_ function in ICvectorfields
or using a sub-grid of equal sized regions of interest as in the _DispField_
function in ICvectorfields. Note that all of the functions in ICvectorfields
that use DIC or variations of it, translate displacement or velocities in
terms of pixel shifts to the spatial units defined in the projection
information of the original input rasters. The coordinate information required
for translation of pixel shifts to the correct spatial units is obtained using
functions in the terra R package [Hijmans, 2021].

### Extending DIC to quantify persistent movement

In applications of DIC in earth systems with persistent directional flows that
influence movement, it is valuable to determine directional movement of
populations of interest that persist for more than one time step. In such
situations, a spatiotemporal array of images with two space dimensions and one
time dimension is required. Often these can be formulated as stacks of raster
images, with each layer in the stack representing spatially referenced
observations for one time step (step one in Figure 1). A variation of DIC
which I call Spatiotemporal Image Correlation (STIC) permits estimation of
persistent directional movement in terms of orthogonal velocity vectors.

In STIC, the three dimensional array is first lagged by duplicating it and
then removing an integer number of layers from the top of one duplicate and
the bottom of the other (steps two and three in Figure 1). The integer lag is
user defined and serves to minimize estimates of zero movement which always
would occur in the absence of a lag. To differentiate the duplicate lagged
arrays, I will refer to the first as the reference array, and the second as
the lagged array. Regions of interest in the reference array are selected and
locations outside the region of interest in the reference array are assigned
values of zero (grey shaded region in steps two and three in Figure 1
represent regions of interest). The reference array and the lagged array are
then dimension reduced by averaging along rows to obtain one pair of two-
dimensional matrices and by averaging along columns to obtain a second pair of
two dimensional matrices (step 4 in Figure 1). The first pair of matrices
comprises row averaged reference and lagged matrices. The second pair of
matrices comprises column averaged reference and lagged matrices. Each matrix
in the two pairs has one space dimension and one time dimension.

Figure 1: The STIC algorithm: The input array in step one is a raster stack of
images in which each image layer represents the phenomenon of interest in
planar space at a different time instance. In step two, the input array is
duplicated and based on a user specified time lag, layers are removed from the
top of one array and the bottom of the other. In addition, a region of
interest is defined in one of the duplicate arrays represented by the grey
shaded region at the top of the prism on the left. In step four the rows are
dimension reduced by averaging along one of the axes (either rows or columns).
This produces a pair or row-averaged matrices and a pair of column-averaged
matrices that are analyzed using cross-correlation to estimate orthogonal
velocity vectors.

Cross-correlation between the pairs of reference and lagged matrices is then
computed as described for DIC. Recall that one dimension of each of the row or
column-averaged matrices is spatial while the other is temporal, which enables
calculation of two orthogonal velocity vectors based on space shifts and time
shifts obtained by application of DIC:

$\displaystyle v_{x}$ $\displaystyle=s_{x}/s_{tx},s_{tx}\neq 0$ (6)
$\displaystyle v_{y}$ $\displaystyle=s_{y}/s_{ty},s_{ty}\neq 0$ (7)

in which $v_{x}$ and $v_{y}$ are velocity in the horizontal and vertical
directions, $s_{x}$ and $s_{y}$ are shifts in the horizontal and vertical
direction, $s_{tx}$ is the time shift that corresponds to spatial shifts in
the horizontal direction, and $s_{ty}$ is the time shift that corresponds to
spatial shifts in the vertical direction. Note that due to the time shift, the
user-defined time lag does not necessarily pre-determine the denominator of
the orthogonal velocity vectors.

### Spatially variable velocities

When the magnitudes of movement velocities are highly spatially variable, a
single time lag is not optimal for quantifying orthogonal velocity vectors.
For these scenarios a variation on the STIC algorithm called STIC+ allows the
user to specify a maximum time lag. The algorithm then repeats the steps
described for STIC for each integer time lag from one to the maximum time lag.
For each repetition and each location of interest, the total velocity
magnitude (speed) is calculated as

$|v|=\sqrt{v_{x}^{2}+v_{y}^{2}}.$ (8)

For each region of interest, the horizontal and vertical velocity vectors are
determined by the time lag STIC calculation that maximizes equation (8).

A summary table describing which functions in ICvectorfields use each of the
algorithms described above is provided (Table 1). Two functions in
ICvectorfields use a standard implementation of DIC similar to that described
by Anuta [1970], two functions use the STIC extension and two functions use
the STIC+ extension (Table 1).

function | algorithm | Images | ROI | Velocities
---|---|---|---|---
_DispField_ | DIC | 2 | grid | variable or not
_DispFieldbb_ | DIC | 2 | bounding box | variable or not
_DispFieldST_ | STIC | 3+ | grid | less variable
_DispFieldSTbb_ | STIC | 3+ | bounding box | less variable
_DispFieldSTall_ | STIC+ | 3+ | grid | more variable
_DispFieldSTbball_ | STIC+ | 3+ | bounding box | more variable
Table 1: ICvectorfields functions, algorithms, and use contexts to facilitate
decisions on which function is most applicable. ROI stands for region of
interest, which is defined either using a grid or a bounding box, Velocities
refers to whether the magnitudes of velocities in the vector field are
presumed to be spatially variable or not.

## Application

In this section I demonstrate the use of DIC and extensions implemented in
ICvectorfields [Goodsman, 2021] using an example in which data were simulated
based on a partial differential equation and using the classic larch budmoth
defoliation data that are embedded in the ncf package [Bjornstad, 2020]. The
data from the simulated example are embedded in ICvectorfields. For
visualization of results, the demonstrations call functions in several R
packages: These include ggplot2 [Wickham, 2016], ggnewscale [Campitelli,
2021a], metR [Campitelli, 2021b], and terra [Hijmans, 2021].

### Demonstration using simulated data

The model used to simulate data to test and demonstrate the functionality of
ICvectorfields is a convection-diffusion equation, which is a partial
differential equation with terms for diffusion, advection, and reaction:

$\frac{\partial u}{\partial t}=\nabla\cdot\left(D\nabla
u\right)-\nabla\cdot\left(\mathbf{v}u\right)+ru,$ (9)

in which $r$ is the growth rate with units of per unit time, $D$ is the
diffusion coefficient with squared spatial units per unit time, $\mathbf{v}$
is the bivariate orthogonal velocity vector in units of space per unit time,
$\nabla$ represents the gradient, and $\nabla\cdot$ represents divergence. The
orthogonal velocity vector is spatially dependent in the simulations that
follow:

$\mathbf{v}=\begin{cases}(0,0.2),&x<0,y\geq 0\\\ (0.2,0),&x\geq 0,y\geq 0\\\
(0,-0.2),&x\geq 0,y<0\\\ (-0.2,0),&x<0,y<0\end{cases}$ (10)

Note that by convention in equations (9) and (10), movement to the right and
up has a negative sign, and movement to the left and down has a positive sign.
This is the opposite convention used in ICvectorfields. Note that the
discontinuities at $x=0$ and $y=0$ in the advection term in equation (10)
create strange model behaviours once concentration reaches $x=0$ or $y=0$, and
so the model was simulated for only 6 time steps to minimize encounters with
these axes. Parameter values for the diffusion coefficient and the growth rate
were $D=(0.01,0.01)$ squared spatial units per unit time and $r=0.5$ per unit
time.

The model in equations (9) and (10) was simulated using the ReacTran R package
[Soetaert and Meysman, 2012], using a finite differencing scheme with backward
differencing on a square domain of 202 cells in each direction, each with a
width of 0.049 spatial units. The initial condition was a concentration of one
units per arbitrary unit of volume in the centre of each quadrant of the
spatial domain. Boundary conditions were zero flux (reflecting) on all four
sides of the spatial domain. The simulation data are saved in table format
within ICvectorfields.

The data are imported and then converted from table format to a raster stack
using the _RastStackData_ helper function. They can then be visualized as
rasters as shown below.

    # import simulated data
    data(SimData, package = "ICvectorfields")

    # convert to raster stack
    SimStack <- ICvectorfields::RastStackData(SimData)

    # confirming dimension
    dim(SimStack)

    # plotting
    layout(matrix(1:6, 2, 3, byrow = TRUE))
    layout.show(6)
    terra::plot(SimStack[[1]], legend = FALSE, main = "t1")
    terra::plot(SimStack[[2]], legend = FALSE, main = "t2")
    terra::plot(SimStack[[3]], legend = FALSE, main = "t3")
    terra::plot(SimStack[[4]], legend = FALSE, main = "t4")
    terra::plot(SimStack[[5]], legend = FALSE, main = "t5")
    terra::plot(SimStack[[6]], legend = FALSE, main = "t6")

The figure showing the simulated spatiotemporal dynamics is identical to
Figure 2.

Figure 2: Visualization of simulation data for six time steps. The initial
condition at t0 is not shown. Green colours represent the highest
concentrations.

To analyze displacement based on a pair of images, I chose to use the standard
implementation of DIC in the _DispField_ function of ICvectorfields. The first
two arguments of the _DispField_ function are the input rasters. The first
input raster is treated as the reference image and the second is treated as
the shifted image. In this case, the first image is the raster layer
corresponding to the first time step (t1) and the second image is the raster
layer corresponding to the sixth time step (t6). The function selects regions
of interests based on a grid of dimensions given in the _factv1_ and _facth1_
arguments, which represent to number of rows and columns in each sub-grid.
Sub-grids start in the upper left corner and as many sub-grids as fit within
the original domain are constructed. In the code below, sub-grids are
$101\times 101$, which is approximately the size of one quarter of the
original spatial domain of the simulation. The restricted argument is by
default set to _FALSE_. In that case, the DIC algorithm cross-correlates each
region of interest in the first image with the entirety of the second image.
When _restricted = TRUE_ The algorithm cross correlates both images only
within the region of interest. If the user has reason to believe that movement
is predominantly occurring within sub-grids the _restricted = TRUE_ option has
the added benefit of speeding up computation.

    VFdf1 <- DispField(SimStack[[1]], SimStack[[6]], factv1 = 101, facth1 = 101,
                        restricted = TRUE)
    VFdf1

The output of _DispField_ is in data table format. Because the data table is
small, I’ve printed it to the console in the code example above. For ease of
viewing in this paper, the output is duplicated in Table 2.

row | col | frowmin | frowmax | fcolmin | fcolmax | centx | centy | dispx | dispy
---|---|---|---|---|---|---|---|---|---
51 | 51 | 1 | 101 | 1 | 101 | -2.5 | 2.5 | 0.00 | -0.98
152 | 51 | 102 | 202 | 1 | 101 | -2.5 | -2.5 | 0.98 | 0.00
51 | 152 | 1 | 101 | 102 | 202 | 2.5 | 2.5 | -0.98 | 0.00
152 | 152 | 102 | 202 | 102 | 202 | 2.5 | -2.5 | 0.00 | 0.98
Table 2: ICvectorfields output from a call of the _DispField_ function using
simulated data. The table is a duplicate of the data table returned after
calling the function except that the headings for the first two columns have
been shortened from rowcent to row and from colcent to col so that the table
fits the page. In addition centx and centy have been rounded to a single
decimal place and dispx and dispy have been rounded to two decimal places.

The directions of movement coincide with the directions of advection in the
simulation with movement downwards in the upper left quadrant (first row of
Table 2), movement to the right in the lower left quadrant (second row of
Table 2), movement to the left in the upper right quadrant (third row of Table
2), and upwards movement in the lower right quadrant (fourth row of Table 2).
Speed of movement can be computed by dividing displacement by the number of
time steps that passed $0.98/5=0.196$, which is slightly slower than the
simulated advection speed of 0.2 spatial units per time step. The discrepancy
is likely due to the blurring effect of diffusion in the partial differential
equation (equation (9)).

In situations where the speed is constant, velocity can be estimated from
pairs of images as I have demonstrated above. However, the _DispFieldST_
algorithm is designed to return orthogonal velocity vectors and so for
confirmation purposes, I demonstrate it below:

    VFdf2 <- DispFieldST(SimStack, lag1 = 1, factv1 = 101, facth1 = 101, restricted = TRUE)
    VFdf2

The data table that is printed to the console after running the code above
looks similar to the data table duplicated in Table 2 except that under the
heading dispx and dispy the algorithm returns horizontal and vertical
velocities rather than displacement vectors. The directions of movement are
the same as those shown in Table 2, but the speed is 0.196 spatial units per
unit time as previously estimated.

The vector field and the raw data can be visualized simultaneously using
plotting functionality in ggplot2 with extensions in the ggnewscale and metR
packages.

    SimVF = ggplot() +
      xlim(c(-5, 5)) +
      ylim(c(-5, 5)) +
      geom_raster(data = SimData,
                  aes(x = xcoord, y = ycoord, fill = t1)) +
      scale_fill_gradient(low = "white", high = "blue", na.value = NA) +
      new_scale("fill") +
      geom_raster(data = SimData,
                  aes(x = xcoord, y = ycoord, fill = t6), alpha = 0.5) +
      scale_fill_gradient(low = "white", high = "red", na.value = NA) +
      geom_vector(data = VFdf2,
                  aes(x = centx, y = centy,
                      mag = Mag(dispx, dispy),
                      angle = Angle(dispx, dispy))) +
      theme_bw()
    SimVF

Figure 3: Vector field for radial movement simulated using a convection-
diffusion equation. the orthogonal velocity vectors are estimated using the
_DispFieldST_ function in ICvectorfields.

The resulting figure is duplicated in Figure 3. The velocity vectors in the
vector field (Figure 3) are consistent with the simulated advection vectors,
although they slightly underestimate movement speed.

Before proceeding to the next demonstration I will illustrate one of the
potential pitfalls of estimating movement based on cross-correlation. If the
argument of the _DispFieldST_ function is left its default _restricted =
FALSE_ configuration, the algorithm will search the entire domain for shifts
that maximize cross-correlation. Because the simulations in each quadrant of
the spatial domain are quite similar, cross-correlation is in fact maximized
by shifts that cross quadrants, even though simulated movement was not that
large. Therefore, calling _DispFieldST_ with _restricted = FALSE_ produces
incorrect output (Table 3): The simulated advection speed is not at all close
to the estimated maximum orthogonal advection speed of 3.87 spatial units per
unit time.

row | col | frowmin | frowmax | fcolmin | fcolmax | centx | centy | dispx | dispy
---|---|---|---|---|---|---|---|---|---
51 | 51 | 1 | 101 | 1 | 101 | -2.5 | 2.5 | 0.10 | -3.87
152 | 51 | 102 | 202 | 1 | 101 | -2.5 | -2.5 | 3.87 | 0.10
51 | 152 | 1 | 101 | 102 | 202 | 2.5 | 2.5 | -3.87 | -0.10
152 | 152 | 102 | 202 | 102 | 202 | 2.5 | -2.5 | -0.10 | 3.87
Table 3: ICvectorfields output from a call of the _DispFieldST_ function using
simulated data. This call is meant to demonstrate a potential pitfall in using
the cross-correlation approach because when _restricted = FALSE_ , the
algorithm finds positive cross-correlations that are not caused by movement.
The table is a duplicate of the data table returned after calling the function
except that the headings for the first two columns have been shortened from
rowcent to row and from colcent to col so that the table fits the page. In
addition centx and centy have been rounded to a single decimal place and dispx
and dispy have been rounded to two decimal places.

### Demonstration using larch budmoth data

Larch budmoths are lepidopteran defoliators that exhibit periodic outbreaks
every 8 to 9 years in the European Alps [Bjørnstad et al., 2002]. The larch
budmoth data originated from survey information collected by the forest
administrative agencies of France, Italy, Switzerland, and Austria from 1961
to 1998. The data record the presence of defoliation by larch budmoth
caterpillars within 1 $\times$ 1 km grid cells (a binary variable). These data
were aggregated up to 20 $\times$ 20 km grid cells so that records at this
spatial scale were population proxies for larch budmoth caterpillar abundance
[Bjørnstad et al., 2002] based on the assumption that defoliation damage is
proportional to the abundance of the causal agents. Grid cells were excluded
from the data set if they exhibited less than one percent defoliation or if
more than ninety percent of years in which data were collected at that
location exhibited no defoliation by larch budmoth [Bjørnstad et al., 2002].
The larch budmoth defoliation data exhibit directional traveling wave-trains
that travel from the southwest to the northeast along the European Alps
[Bjørnstad et al., 2002]. These data are embedded in the ncf R package.

After loading the ncf library as well as ICvectorfields, the data can be
loaded, converted to a raster stack and visualized as follows:

    # import larch budmoth data
    data(lbm, package = "ncf")

    # convert to raster stack
    LBMStack <- ICvectorfields::RastStackData(lbm)

    # confirming dimension
    dim(LBMStack)

    # visualizing
    layout(matrix(1:6, 2, 3, byrow = TRUE))
    layout.show(6)
    terra::plot(LBMStack[[1]], legend = FALSE, main = "1961")
    terra::plot(LBMStack[[2]], legend = FALSE, main = "1962")
    terra::plot(LBMStack[[3]], legend = FALSE, main = "1963")
    terra::plot(LBMStack[[4]], legend = FALSE, main = "1964")
    terra::plot(LBMStack[[5]], legend = FALSE, main = "1965")

Figure 4: The first five years of the larch budmoth defoliation data included
in ncf. Green colours represent the highest level of defoliation.

This code plots the first five years of the data set (Figure 4), which show a
standard progression of outbreaks from the southwest corner of the Alps to the
northeast. This pattern repeats relatively regularly every 8 to 10 years in
the data set.

The study region covers a large geographic area and so it is likely the
population movement speeds vary geographically. For this reason, I elected to
use the STIC+ algorithm to analyze the data using _DispFieldSTall_. In the
code below I analyze the first 23 years of the time series (1961 to 1983) as
defoliation patterns from 1984 to 1998 are less regular.

    VFdf3 <- DispFieldSTall(LBMStack[[1:23]], lagmax = 3, factv1 = 3, facth1 = 3, restricted = FALSE)

    LBMVF1 = ggplot() +
      geom_tile(data = lbm,
                  aes(x = x, y = y,
                      fill = X1962)) +
      scale_fill_gradient(low = "white", high = "blue", na.value = NA) +
      new_scale("fill") +
      geom_tile(data = lbm,
                  aes(x = x, y = y, fill = X1964), alpha = 0.5) +
      scale_fill_gradient(low = "white", high = "red", na.value = NA) +
      geom_vector(data = VFdf3,
                  aes(x = centx, y = centy,
                      mag = Mag(dispx, dispy),
                      angle = Angle(dispx, dispy))) +
      theme_bw()
    LBMVF1

Calling _DispFieldSTall_ returns a data frame object that is convenient for
plotting the vector field. The vector field reveals that moth movement is to
the north on the southwestern side of the Alps and then to the east on the
northern side of the Alps (Figure 5). It also reveals deceleration as
outbreaks turn the corner and then acceleration as outbreaks move eastward
(Figure 5).

Figure 5: Vector field for Larch Budmoth persistent movement. The orthogonal
velocity vectors are estimated using the _DispFieldSTall_ function in
ICvectorfields. Blue colours show the locations and intensities of defoliation
in 1962 and red colours show the locations and intensities of defoliation in
1964. Vectors have their own scale that is distinct from the scale of the map.

The average speed of larch budmoth movement can be computed from the data
frame output of _DispFieldSTall_ as follows:

    VFdf3$speed <- sqrt((VFdf3$dispx^2) + VFdf3$dispy^2)

    # subsetting to remove locations where speed is zero
    VFdf4 <- subset(VFdf3, speed > 0)

    # computing mean, standard deviation and dimension of data frame
    # to obtain sample size
    mean(VFdf4$speed)
    sd(VFdf4$speed)
    dim(VFdf4)

    # upper and lower Wald-type 95 percent confidence interval on average speed
    mean(VFdf4$speed)/1000 + qt(0.975, dim(VFdf4)[1] - 1)*sd(VFdf4$speed)/1000/sqrt(dim(VFdf4)[1] - 1)
    mean(VFdf4$speed)/1000 + qt(0.025, dim(VFdf4)[1] - 1)*sd(VFdf4$speed)/1000/sqrt(dim(VFdf4)[1] - 1)

Using the approach above, the average movement speed is estimated as $176\pm
43$ $\text{km}(\text{Yr})^{-1}$, an estimate that is less than the previous
speed estimates for northeastern spread of 220 $\text{km}(\text{Yr})^{-1}$
[Bjørnstad et al., 2002] and 254 $\text{km}(\text{Yr})^{-1}$ [Johnson et al.,
2004]. The difference between estimates in the literature and estimates
produced here are likely due to the direction of movement. The vectors in the
larch budmoth vector field point predominantly north and east (Figure 5). In
other words they are orthogonal. In contrast the movement speeds estimated by
Bjørnstad et al. [2002] and Johnson et al. [2004] are projected along lines
that run to the northeast. A simple application of geometry reveals that an
average speed of 176 $\text{km}(\text{Yr})^{-1}$ in the north and east
directions corresponds to an estimated speed of 249
$\text{km}(\text{Yr})^{-1}$ in the northeastern direction (Pythagorean theorem
$\sqrt{176^{2}+176^{2}}$). This estimate is consistent with prior speed
estimates for larch budmoth population movement [Bjørnstad et al., 2002,
Johnson et al., 2004].

## Summary

The ICvectorfields R package implements standard Digital Image Correlation
algorithms in addition to a novel extension that permits estimation of
orthogonal velocities of persistent movement in series of three or more
images. Here I demonstrate the usefulness of DIC and the extension implemented
in ICvectorfields in a new arena: Whereas DIC is often applied in engineering
and materials science to quantify the effects of force application on
materials [Sutton et al., 2009], it has not been used in landscape ecology. In
this field, the approach has potential to provide new insights into how
populations move across landscapes and to demonstrate the untenable nature of
assumptions of homogeneity inherent in most analyses based on the traveling
wave paradigm. Even when models of sufficient complexity to capture
environmental heterogeneity can be used, I expect that the methods in
ICvectorfields will be useful because they facilitate comparison between
modeled and empirical population movement data as demonstrated in the partial
differential equation example in this study. Approaches such as this one that
estimate movement based on cross-correlation, however, have a weakness: Under
certain circumstances, they are prone to finding cross-correlations that are
unrelated to movement as was demonstrated in this paper. For this reason,
users must exercise vigilance in interpreting the results of vector field
analyses like those demonstrated herein. If possible, results should be
checked against a standard or against prior published results regarding
movement propensity. Nevertheless, the methods described here hold promise for
exploratory analyses, hypothesis generation, and synoptic pattern analyses of
population movements.

## References

  * Alonso et al. [2002] David Alonso, Frederic Bartumeus, and Jordi Catalan.  Mutual interference between predators can give rise to turing spatial patterns.  _Ecology_ , 83(1):28–34, 2002.
  * Anuta [1970] Paul E Anuta.  Spatial registration of multispectral and multitemporal digital imagery using fast fourier transform techniques.  _IEEE transactions on Geoscience Electronics_ , 8(4):353–368, 1970.
  * Avgar et al. [2016] Tal Avgar, Jonathan R Potts, Mark A Lewis, and Mark S Boyce.  Integrated step selection analysis: bridging the gap between resource selection and animal movement.  _Methods in Ecology and Evolution_ , 7(5):619–630, 2016.
  * Bjornstad [2020] Ottar N. Bjornstad.  _ncf: Spatial Covariance Functions_ , 2020.  URL https://CRAN.R-project.org/package=ncf.  R package version 1.2-9.
  * Bjørnstad and Falck [2001] Ottar N Bjørnstad and Wilhelm Falck.  Nonparametric spatial covariance functions: estimation and testing.  _Environmental and Ecological Statistics_ , 8(1):53–70, 2001.
  * Bjørnstad et al. [2002] Ottar N Bjørnstad, Mikko Peltonen, Andrew M Liebhold, and Werner Baltensweiler.  Waves of larch budmoth outbreaks in the european alps.  _Science_ , 298(5595):1020–1023, 2002.
  * Buderman et al. [2016] Frances E Buderman, Mevin B Hooten, Jacob S Ivan, and Tanya M Shenk.  A functional model for characterizing long-distance movement behaviour.  _Methods in Ecology and Evolution_ , 7(3):264–273, 2016.
  * Campitelli [2021a] Elio Campitelli.  _ggnewscale: Multiple Fill and Colour Scales in ’ggplot2’_ , 2021a.  URL https://CRAN.R-project.org/package=ggnewscale.  R package version 0.4.5.
  * Campitelli [2021b] Elio Campitelli.  _metR: Tools for Easier Analysis of Meteorological Fields_ , 2021b.  URL https://github.com/eliocamp/metR.  R package version 0.9.2.
  * Frigo and Johnson [2005] Matteo Frigo and Steven G. Johnson.  The design and implementation of FFTW3.  _Proceedings of the IEEE_ , 93(2):216–231, 2005\.  Special issue on “Program Generation, Optimization, and Platform Adaptation”.
  * Goodsman [2021] Devin Goodsman.  _ICvectorfields: Vector Fields from Spatial Time Series of Population Abundance_ , 2021.  URL https://CRAN.R-project.org/package=ICvectorfields.  R package version 0.0.2.
  * Hassell et al. [1991] Michael P Hassell, Hugh N Comins, and Robert M Mayt.  Spatial structure and chaos in insect population dynamics.  _Nature_ , 353(6341):255–258, 1991.
  * Hassell et al. [1994] Michael P Hassell, Hugh N Comins, and Robert M May.  Species coexistence and self-organizing spatial dynamics.  _Nature_ , 370(6487):290–292, 1994.
  * Hijmans [2021] Robert J. Hijmans.  _terra: Spatial Data Analysis_ , 2021.  URL https://CRAN.R-project.org/package=terra.  R package version 1.2-10.
  * Hoffman et al. [2006] Aaron L Hoffman, Julian D Olden, Jeremy B Monroe, N LeRoy Poff, Todd Wellnitz, and John A Wiens.  Current velocity and habitat patchiness shape stream herbivore movement.  _Oikos_ , 115(2):358–368, 2006.
  * Johnson et al. [2004] Derek M Johnson, Ottar N Bjørnstad, and Andrew M Liebhold.  Landscape geometry and travelling waves in the larch budmoth.  _Ecology Letters_ , 7(10):967–974, 2004.
  * Johnson [2009] Neil Johnson.  _Simply complexity: A clear guide to complexity theory_.  Simon and Schuster, 2009.
  * Joo et al. [2020] Rocio Joo, Matthew E Boone, Thomas A Clay, Samantha C Patrick, Susana Clusella-Trullas, and Mathieu Basille.  Navigating through the r packages for movement.  _Journal of Animal Ecology_ , 89(1):248–267, 2020\.
  * Kolmogorov [1937] Andrei N Kolmogorov.  Étude de l’équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique.  _Bull. Univ. Moskow, Ser. Internat., Sec. A_ , 1:1–25, 1937\.
  * Kot et al. [1996] Mark Kot, Mark A Lewis, and Pauline van den Driessche.  Dispersal data and the spread of invading organisms.  _Ecology_ , 77(7):2027–2042, 1996.
  * Krainski et al. [2018] Elias Krainski, Virgilio Gómez-Rubio, Haakon Bakka, Amanda Lenzi, Daniela Castro-Camilo, Daniel Simpson, Finn Lindgren, and Håvard Rue.  _Advanced spatial modeling with stochastic partial differential equations using R and INLA_.  Chapman and Hall/CRC, 2018.
  * Parrish and Edelstein-Keshet [1999] Julia K Parrish and Leah Edelstein-Keshet.  Complexity, pattern, and evolutionary trade-offs in animal aggregation.  _Science_ , 284(5411):99–101, 1999.
  * R Core Team [2021] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2021.  URL https://www.R-project.org/.
  * Rahim [2021] Karim Rahim.  _fftwtools: Wrapper for ’FFTW3’ Includes: One-Dimensional, Two-Dimensional, Three-Dimensional, and Multivariate Transforms_ , 2021.  URL https://CRAN.R-project.org/package=fftwtools.  R package version 0.9-11.
  * Ruan [1998] SHIGUI Ruan.  Turing instability and travelling waves in diffusive plankton models with delayed nutrient recycling.  _IMA journal of applied mathematics_ , 61(1):15–32, 1998.
  * Signer et al. [2019] Johannes Signer, John Fieberg, and Tal Avgar.  Animal movement tools (amt): R package for managing tracking data and conducting habitat selection analyses.  _Ecology and Evolution_ , 9:880–890, 2019.  URL https://doi.org/10.1002/ece3.4823.
  * Skellam [1951] John Gordon Skellam.  Random dispersal in theoretical populations.  _Biometrika_ , 38(1/2):196–218, 1951.
  * Soetaert and Meysman [2012] Karline Soetaert and Filip Meysman.  Reactive transport in aquatic ecosystems: Rapid model prototyping in the open source software r.  _Environmental Modelling & Software_, 32:49–60, 2012\.
  * Soetaert et al. [2010] Karline Soetaert, Thomas Petzoldt, and R. Woodrow Setzer.  Solving differential equations in R: Package deSolve.  _Journal of Statistical Software_ , 33(9):1–25, 2010.  doi: 10.18637/jss.v033.i09.
  * Sutton et al. [2009] Michael A Sutton, Jean Jose Orteu, and Hubert Schreier.  _Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications_.  Springer Science & Business Media, 2009.
  * Urban et al. [2008] Mark C Urban, Ben L Phillips, David K Skelly, and Richard Shine.  A toad more traveled: the heterogeneous invasion dynamics of cane toads in australia.  _The American Naturalist_ , 171(3):E134–E148, 2008.
  * Wickham [2016] Hadley Wickham.  _ggplot2: Elegant Graphics for Data Analysis_.  Springer-Verlag New York, 2016.  ISBN 978-3-319-24277-4.  URL https://ggplot2.tidyverse.org.
  * Zammit-Mangion [2019] Andrew Zammit-Mangion.  _IDE: Integro-Difference Equation Spatio-Temporal Models_ , 2019.  URL https://CRAN.R-project.org/package=IDE.  R package version 0.3.0.


# A Large-Scale Annotated Multivariate Time Series Aviation Maintenance
Dataset from the NGAFID

Hong Yang
Rochester Institute of Technology
Rochester, NY 14623
<EMAIL_ADDRESS>
Travis Desell
Rochester Institute of Technology
Rochester, NY 14623
<EMAIL_ADDRESS>

###### Abstract

This paper presents the largest publicly available, non-simulated, fleet-wide
aircraft flight recording and maintenance log data for use in predicting part
failure and maintenance need. We present 31,177 hours of flight data across
28,935 flights, which occur relative to 2,111 unplanned maintenance events
clustered into 36 types of maintenance issues. Flights are annotated as before
or after maintenance, with some flights occurring on the day of maintenance.
Collecting data to evaluate predictive maintenance systems is challenging
because it is difficult, dangerous, and unethical to generate data from
compromised aircraft. To overcome this, we use the National General Aviation
Flight Information Database (NGAFID), which contains flights recorded during
regular operation of aircraft, and maintenance logs to construct a part
failure dataset. We use a novel framing of Remaining Useful Life (RUL)
prediction and consider the probability that the RUL of a part is greater than
2 days. Unlike previous datasets generated with simulations or in laboratory
settings, the NGAFID Aviation Maintenance Dataset contains real flight records
and maintenance logs from different seasons, weather conditions, pilots, and
flight patterns. Additionally, we provide Python code to easily download the
dataset and a Colab environment to reproduce our benchmarks on three different
models. Our dataset presents a difficult challenge for machine learning
researchers and a valuable opportunity to test and develop prognostic health
management methods.

##  1 Introduction

Maintenance related issues are a notable risk in the operation single-engine
fixed-wing aircraft Goldman et al., (2002) for private use, known as General
Aviation (GA). With improvements in data collection and processing
capabilities, there is now an opportunity to indirectly monitor the condition
of a variety of components in GA aircraft for the purposes of Predictive
Maintenance (PM).

The domain of prognostics and health management (PHM) has shifted towards the
use of more data driven approaches Tsui et al., (2015), sometimes also
involving deep learning methods. These methods aim to predict the RUL, detect
faults, or monitor the condition of a system and its parts. However, it is
particularly difficult to collect data for these purposes. Rezaeianjouybari
and Shang, (2020) states that the data collection process is often too time
consuming and too perfect (the data collected in laboratory conditions may not
translate well to the real world). It is clear that there is a lack of high
quality, real world data, yet there is a great need for such data to evaluate
new PHM methodologies.

Existing datasets for PHM often use simulated data Liu et al., (2012) or model
systems in controlled environments, such as a power plant data in the 2015
Prognostics and Health Management Society Data Challenge. However, PHM
methodologies can be extended to real world systems in largely uncontrolled
environments, such as cars or aircraft. As vehicles collect more data, see
Arena et al., (2021), there is an opportunity to combine said data with
unplanned maintenance records to create a novel RUL estimation dataset for
predictive maintenance.

GA flight data is particularly challenging for PHM due to the nature of the
aircraft sensors. Unlike commercial aircraft, GA aircraft only possess basic
sensors for monitoring critical systems, such as the engine and battery, and
flight instruments, for measuring air speed and altitude. A predictive model
relying on such basic sensors for predicting RUL of specific parts must be
both highly sensitive to the most subtle of changes and highly robust to
noise. This is because a significant portion of the variance in flight data is
explained by pilot action and not the condition of the airplane’s parts. This
poses a difficult problem for conventional approaches, such as auto encoders.

In this paper we present the NGAFID Aviation Maintenance Dataset for use in
binary RUL estimation. This dataset is constructed using flight records and
unplanned maintenance logs of a single flight school over multiple years. For
each unplanned breakdown of a part, which lead to unplanned maintenance, there
are multiple flights days before and after the maintenance event. From here,
we create two groups of flight data, one with parts that are going to break
within $x$ days and ones with new parts, which are not expected to break
within $x$ days. This challenging PHM problem uses one of the largest real-
world PHM datasets and it also poses a particularly difficult machine learning
problem for time series classification and time series anomaly detection.

###  1.1 Our Contributions

Data: https://doi.org/10.5281/zenodo.6624956 and
https://www.kaggle.com/datasets/hooong/aviation-maintenance-dataset-from-the-
ngafid A dask dataframe containing 31,177 hours of flight data across 28,935
flights, with a header csv that describes each flight and the associated
maintenance file. The data was collected automatically by a flight school
after each flight and uploaded to the NGAFID database. Flights occurred in a
variety of seasons and weather conditions. All flights were flown with Cessna
172 aircraft. 23 sensors record data every second, resulting in more than 100
million rows of flight data with a total size of 4.3 GB.

Benchmarks: https://github.com/hyang0129/NGAFIDDATASET A repository containing
helper code to automatically download and process the dataset and 2 Colab
notebooks for replicating the benchmark experiments. Anyone can run the
benchmark experiments with one click using the Colab notebooks in a web
browser, with a free Linux environment including GPUs and TPUs provided by
Google. This demonstrates how to use the dataset and replicate experiments,
regardless of the replicator’s hardware and software limitations.

The above files are licensed under GNU General Public License V3.0.

##  2 Related Work

###  2.1 Aircraft and Aircraft Maintenance Datasets

Publicly available aircraft maintenance datasets often use simulated data
rather than data collected from real life events. Commercial Modular Aero-
Propulsion System Simulation (C-MAPSS) Liu et al., (2012) is a commonly used
simulation program for aircraft engines. These simulations are sometimes
supplemented by real flight conditions, e.g. Arias Chao et al., (2021), where
flight data is used as the input for the simulation. However, these approaches
rely on synthetic data, which may not reflect the noisy nature of data
collection. There are papers that reference real data, but they do not
publicly release their datasets, e.g. Celikmih et al., (2020) and Dangut et
al., (2022). The task of detecting anomalies from flight data was investigated
by Chu et al., (2010), however the dataset was generated from a flight
simulator. Yildirim and Kurt, (2018) collects real world data from aircraft,
but only during normal operation of the aircraft and not when there is a
fault. A common difficulty in collecting real world failure data is the
associated cost and safety risk; it would be unethical to ask pilots to fly
faulty aircraft. Yang et al., (2021) presented a smaller version of this
dataset for aviation maintenance.

Aircraft flight data has also been used for other tasks, aside from
maintenance prediction. Klein, (1989) used real data to estimate aerodynamic
properties of aircraft. Khadilkar and Balakrishnan, (2012) used aircraft
flight recordings to estimate fuel usage during the taxing stage. Aircraft
flight data can also be used to improve estimation of the remaining range for
the aircraft, Randle et al., (2011). In most PHM dataests, Liu et al., (2012),
the goal is to predict an estimate or probability density function of the RUL,
prior to a fault occuring. Prior to this work, there were no publicly
available fleet-wide maintenance records of this size and scope coupled with
flight data recordings in both fault free and faulty states.

###  2.2 Time Series Datasets

There are many time series classification datasets, such as Bagnall et al.,
(2018). These datasets are used to evaluate both deep learning methods, such
as Inception Time, Fawaz et al., (2020), and non deep methods, such as Mini
Rocket, Dempster et al., (2021) and HIVE-COTE, Lines et al., (2018). Time
series classification could also include audio datasets, such as Vincent et
al., (2013), but such datasets belong in their own domain of audio processing.
As shown by Bagnall et al., (2018), a collection of time series datasets will
contain very diverse types of data, such as motion sensor data,
electroencephalogram data, traffic data and more. Few time series datasets
exist that contain as many training examples with so many data points per
training example.

###  2.3 Classification Methods

Several methods have been developed for MTS classification, for a review see
Fawaz et al., (2019). Notable non-deep learning methods include distance based
k-nearest neighbors by Orsenigo and Vercellis, (2010) and Dynamic Time Warping
KNN by Seto et al., (2015). For deep learning methods, well performing MTS
classifiers tend to utilize some combination of Recurrent Neural Network (RNN)
and Convolutional Neural Network (CNN) methods, Karim et al., (2017), or CNN
only methods Wang et al., (2017). However, RNN methods struggle with long
sequences due to the vanishing gradient problem, as mentioned by Le and
Zuidema, (2016). CNN only methods can provide strong results for MTS
classification, as shown by Assaf et al., (2019), but may struggle when
relevant features are temporally sparse and related. Yang et al., (2021) also
used multi-headed self-attention methods.

##  3 Dataset Description

###  3.1 Data Acquisition

The full dataset is constructed using flight recordings from the NGAFID
database and maintenance records from MaintNet. Each maintenance record
describes the type of maintenance performed on an aircraft and the time period
when maintenance occurred. From this information, we extract that aircraft’s
flights occurring before and after the maintenance period. Only the first 5
flights before and after and any flights during maintenance are extracted. All
maintenance activities were unplanned and done on request. However, when
issues are detected, FAA policy forbids the plane from flying. This means that
issues are always detected sometime during or after the first flight before
maintenance. Each flight record contains readings from 23 different sensors
every second. A description of the sensors can be found in Table 1. All flight
data recordings come from the same aircraft model, the Cessna 172.

Column Name | Description
---|---
volt1 | Main electrical system bus voltage (alternators and main battery)
volt2 | Essential bus (standby battery) bus voltage
amp1 | Ammeter on the main battery (+ charging, - discharging)
amp2 | Ammeter on the standby battery (+ charging, - discharging)
FQtyL | Fuel quantity left
FQtyR | Fuel quantity right
E1 FFlow | Engine fuel flow rate
E1 OilT | Engine oil temperature
E1 OilP | Engine oil pressure
E1 RPM | Engine rotations per minute
E1 CHT1 | 1st cylinder head temperature
E1 CHT2 | 2nd cylinder head temperature
E1 CHT3 | 3rd cylinder head temperature
E1 CHT4 | 4th cylinder head temperature
E1 EGT1 | 1st Exhaust gas temperature
E1 EGT2 | 2nd Exhaust gas temperature
E1 EGT3 | 3rd Exhaust gas temperature
E1 EGT4 | 4th Exhaust gas temperature
OAT | Outside air temperature
IAS | Indicated air speed
VSpd | Vertical speed
NormAc | Normal acceleration
AltMSL | Altitude miles above sea level
Table 1: Description of the data collected by aircraft sensors

The NGAFID serves as a repository for general aviation flight data, with a web
portal for viewing and tracking flight safety events for individual pilots as
well as for fleets of aircraft Karboviak et al., (2018). The NGAFID currently
contains over 900,000 hours of flight data generated by over 780,000 flights
by 12 different types of aircraft, provided by 65 fleets and individual users,
resulting in over 3.15 billion per second flight data records across 103
potential flight data recorder parameters. Five years of textual maintenance
records from a fleet which provided data to the NGAFID have been clustered by
maintenance issue type and then validated by domain experts for the MaintNet
project, see Akhbardeh et al., (2020).

MaintNet’s maintenance record logbook data was clustered into 36 different
maintenance issue types. The count of flights per issue type is located in
table 3 in Appendix A. Although some maintenance issues occur very rarely, all
maintenance issues with flight data are included in the full dataset. It is
important to note that the NGAFID collects real flight data from aircraft
flying with potentially faulty parts (as is the case for any real world fleet
of aircraft). This is because individual components may fail without causing
catastrophic failure of the aircraft during regular operation. The collection
of this data poses no additional safety risk to the pilots because data
collection occurs for all flights performed.

###  3.2 Dataset Subset for Benchmarking

Attempting to train a model on the full dataset is particularly difficult.
This is because of two main issues. The first is the under representation of
certain classes. For example, the spark plug related issue contains only 15
flights. These classes are included in the dataset to allow for future
researchers to train models that address the class imbalance issue, as class
imbalance is an important area of machine learning research Johnson and
Khoshgoftaar, (2019). The second is that each flight contains a significant
amount of data, which makes it hard to regularize a model trained on the data.
In the spark plug example, each of the 15 flights contain more than 40,000
data points, which would be used to predict a single class out of 36.
Regularization is another important research area of machine learning, Kukačka
et al., (2017), but it is beyond the scope of this paper. These two challenges
are inevitable when collecting real world data, as it is impossible to
guarantee that each part in a real world aircraft fails at the same rate.

A subset of the full dataset is used in this paper to benchmark time series
classification methods and provide a baseline. Since addressing class
imbalance and regularization are beyond the scope of this paper, the subset is
designed to minimize the impact of those two issues. By limiting the subset to
only classes that have at least 50 flights, it removes flights where
generalization is a problem and reduces the impact of class imbalance. The
subset also defines the binary RUL problem as $P(RUL>x)$ with $x$ as 2 days.
We label flights within 2 days before maintenance as $P(RUL>x)=0$ and flights
after 2 days after maintenance as $P(RUL>x)=1$, as such flights contain brand
new parts.

The data subset is defined as all flights within two days of the maintenance
period, excluding any flights during the maintenance period. Furthermore, any
flights belonging to classes that would have fewer than 50 flights before or
after the maintenance period are also excluded. This leaves us with a total of
5844 flights after maintenance and 5602 flights before maintenance in 19
different classes of maintenance issues. The count of flights per issue type
is located in table 3 in Appendix A.

###  3.3 Data Preprocessing

The full dataset is provided without any preprocessing and contains the full
flight data for each flight. Some sensors recorded NaN values at certain time
steps. This can be caused by many factors, including but not limited to
initialization of the sensor at the beginning of the flight, failure of the
sensor during flight, and failure of the recording system during flight.
Approximately 1% of all datapoints contain NaN values. We decided to present
future researchers with NaN values as replacing such values with another, such
as 0, would change the meaning.

For training a model in our benchmark experiments, the data subset is scaled.
All values are scaled via MinMax, with a scaled minimum value of 0 and a
scaled maximum value of 1. Minimums and maximums for each channel of data were
calculated using all of the data. Please note that the full dataset is not
distributed with its values scaled, which allows for examination of different
normalization techniques by users of the dataset.

###  3.4 Visualization

Readings from 4 sensors from two flights are provided in figure 1. Only the
output from 4 sensors are included to keep the graph readable. Even with only
4 sensors, it should be clear that determining the probability of part failure
for a human would be extremely difficult. Even with domain knowledge, it would
only be possible to detect that an issue has occurred, but not that an issue
is going to occur.

Figure 1: Exhaust Gas Temperature before and after Maintenance for Intake
Gasket Issue. This shows that the sensor readings before an issue occurs and
after it is fixed is largely the same. Notice that the EGT values are mostly
in line with each other, but diverge at certain times. It is not obvious where
an anomaly occurs, nor is it obvious whether or not there is an issue. The
difference in temperature values indicate a difference in flight plan as well.

###  3.5 Data Privacy

To ensure privacy, information regarding the serial number of the aircraft of
the flight and the date of the flight have been removed from the publicly
released dataset. Latitude and longitude are also removed from the flight. The
specific text of the maintenance logs are also withheld.

###  3.6 Labels

Labels are assigned to flights based on the date of the flight and the date of
the maintenance. As per the Federal Aviation Administration regulations,
aircraft are only permitted to operate in a safe manner. This means that
maintenance issues and part failures that occurred prior to the maintenance
date will be fixed after the maintenance date.

Part failures in the aircraft can be described as acute or chronic. Acute part
failure describes a complete failure of a part due to an unexpected event. For
example, if the aircraft collides with a small object, such as a bird, this
would create data representing an acute failure and an instantaneous change in
the remaining useful life. Chronic part failure describes gradual wear and
tear that renders a part unsafe for flight, before scheduled maintenance can
replace the part. For example, an intake gasket may wear down more quickly
than expected, leading to a leakage that negatively impacts the safety of the
aircraft.

It should be noted that dramatic decreases in remaining useful life are rare.
Based on analysis of the intake gasket leak/damage class, there were 9085
cases were described as leaking and only 15 cases were described as torn. This
may negatively impact anomaly detection methods and will be discussed in a
later section.

Because all maintenance records are related to unplanned maintenance, it is
clear that the related parts reached the end of their remaining useful life
prematurely. These unplanned maintenance events occurred outside of scheduled
maintenance and can be considered anomalies in the expect wear and tear of the
associated parts.

###  3.7 Problem Structure

####  3.7.1 Detection of Maintenance Issues

The detection of maintenance issues and their associated part failure is quite
challenging. We define this as $\min_{\forall i\in I}P(RUL_{i}>x)$, where $I$
represents the set of all parts and $RUL_{i}$ is the remaining useful life of
part $i$.

While this problem may resemble anomaly detection, it is important to note
that a significant portion of variance in the data is caused by outside
factors. Pilot actions can dramatically alter readings for every sensor on the
aircraft, causing no two flights to be the same. Weather conditions, which may
vary with altitude, are also inconsistent. For these reasons, an anomaly
detection approach may find a large number of anomalous sections on any given
flight. These anomalies may be caused by a maintenance issue or irregular
flight activities.

It should be clear that the data collected by aircraft sensors pose a greater
challenge than sensors on laboratory or industrial assets. In those
situations, the environment is much more controlled, as opposed to the
environment of an aircraft in flight. They may also utilize dedicated sensors
designed to collect data for the purpose of detecting specific issue. This is
not the case for our dataset, as the sensors collect general flight data,
which were not originally designed to collect data for maintenance issue
prediction and binary RUL estimation.

####  3.7.2 Classification of Maintenance Issues

An extension of the previous problem is to also classify the type of
maintenance issue based on the maintenance record left by the mechanic. This
problem is significantly more difficult for two reasons. First is the class
imbalance, with some classes having thousands of flights and others having
less than one hundred. Second is the similarity in how maintenance issues
affect flight characteristics. This is because two different issues may cause
very similar changes in flight characteristics.

Formally, this is defined as $\forall i\in I[P(RUL_{i}>x)]$. This is
significantly more difficult than detection of the maintenance issue.

####  3.7.3 Hierarchical Classification

In addition to the class labels for specific maintenance issues, this dataset
is provided with a simple hierarchy for maintenance issues. There are 5
groupings of maintenance issues: engine, baffle, oil, cylinder, and other.
Argyriou et al., (2006) has shown that hierarchies can provide benefits in
terms of regularization, as they improve the information provided by the
labels. These hierarchies are included for potential future research.

##  4 Benchmark Experiments

###  4.1 Dataset and Task Definition

Three different tasks are defined and all three tasks use the same subset of
data. The first task is maintenance issue detection. Here, we assign after
maintenance flights with a negative label and before maintenance flights with
a positive label. The second task is maintenance issue classification, where
we assign after maintenance flights with class 0 and before maintenance
flights with their maintenance issue class (from 1 to 19). The third is
combined detection and classification, where a network is trained on both
tasks simultaneously.

###  4.2 Models

####  4.2.1 Mini Rocket

Mini Rocket Dempster et al., (2021) is an improvement over the original Rocket
classifier, Dempster et al., (2020). While the original Rocket transforms the
input time series using a large number of random convolutional kernels and
trains a linear classifier on the transformed features, Mini Rocket improves
upon this by using a smaller, fixed set of kernels. This method has performed
quite well on many time series datasets and trains much more quickly than deep
learning methods. We evaluate a GPU implementation of Mini Rocket by Oguiza,
(2022).

####  4.2.2 HIVECOTEv2

We attempted to test HIVECOTEv2 Middlehurst et al., (2021) on the NGAFID
maintenance dataset using the sktime implementation by Löning et al., (2022).
However, the training time for the model exceeded the maximum time allocated
for Google Colab instances (24hrs). Shifaz et al., (2020) notes that HIVECOTE
scales polynomially with data, which could explain why it could not train in
time on the NGAFID maintenance dataset. Not only is the NGAFID maintenance
dataset large, it also has many time steps and many channels. No results for
HIVECOTEv2 are included.

####  4.2.3 Inception Time

Convolutional networks are very popular in the domain of computer vision, but
can also perform well in time series classification. Fawaz et al., (2020)
proposed the InceptionTime model as an ensemble of five Inception models. Each
Inception model is composed of two residual blocks, each containing three
Inception modules, followed at the very end by a global average pooling layer
and a dense classification head layer. Inception modules contain convolutions
of various kernel sizes and a bottleneck layer. The residual blocks help
mitigate the vanishing gradient issue by allowing for direct gradient flow He
et al., (2016). Fawaz et al., (2020) noted that ensembling was necessary due
to the high standard deviation in accuracy of single Inception models and the
small size of time series datasets.

For this study, we evaluate the Inception model without ensembling, so
references to InceptionTime refer to just the Inception model, without
ensembling. This is because the compute cost is reduced five-fold allowing for
a more efficient use of limited resources.

####  4.2.4 ConvMHSA

Multi-Headed Self Attention (MHSA) modules were popularized by Devlin et al.,
(2018) for usage in Natural Language Processing and by Dosovitskiy et al.,
(2020) for Computer Vision. The ConvMHSA model used for benchmarking is the
same as the one in Yang et al., (2021). The model implements attention layers
that mimic the functionality of the encoder layers present in BERT Devlin et
al., (2018). Instead of token embeddings, the model generates sequence
embeddings with the use of 1D convolutions along the temporal dimension. These
learnable sequence embeddings capture local relationships and to compress the
MTS to a shorter length. It uses a series of 1D convolutions to reduce the
temporal resolution from 4096 to 512 and then employs 4 stacked MHSA encoder
layers with 8 heads each and 64 dense units per head. The output is globally
average pooled and fed to a dense layer for classification.

####  4.2.5 Model Configuration

All results reported were generated using a Google Colab instance with a v2-8
TPU. All models were trained for 200 epochs with 200 steps per epoch using a
batch size of 128, with 5-fold cross validation. Flights are truncated to the
last 4096 time steps and padded to be of the same size. Models used an Adam
optimizer with a learning rate of 3e-5 for CONVMHSA and 1e-4 for
InceptionTime.

Mini-rocket was trained using a Google Colab instance with a Nvidia V100 GPU.
Models were trained for 200 epochs with 143 steps and a batch size of 64, with
5 fold cross validation. Flights are truncated to the last 4096 time steps and
padded to be of the same size. Models used an Adam optimizer with a learning
rate of 2.5e-5.

Notebooks used to train the models are available on Github and can be run in
Colab through a web browser. Colab notebooks can be exported as regular
Jupyter notebooks to run the experiments locally on a GPU. See
https://github.com/hyang0129/NGAFIDDATASET for benchmark instructions.

###  4.3 Results

A summary of results can be found in Table 2. Overall, models tend to overfit
the data, but especially so on classifying the specific maintenance issue.
This is because certain classes contain only 75 examples. MiniRocket performed
relatively poorly compared to the deep learning models. If we compare the
training loss and the validation loss, it appears that MiniRocket may have
difficulty in generalization in the multi class case. It also seems to suffer
from under fitting in the binary case. These results suggest that deep
learning methods may have an advantage for this type of problem.

Model | Task | Binary Acc. | Multiclass Acc. | Loss | Train Loss
---|---|---|---|---|---
ConvMHSA | Binary | 76.0% | N/A | 0.526 | 0.003
Multi | N/A | 52.8% | 2.168 | 1.097
Both |  76.1% |  56.3% | 1.756 | 1.377
InceptionTime | Binary | 75.5% | N/A | 0.569 | 0.214
Multi | N/A | 54.1% | 2.251 | 1.365
Both | 74.0% | 55.4% | 1.667 | 1.038
MiniRocket | Binary | 59.8% | N/A | 0.667 | 0.395
Multi | N/A | 50.4% | 1.800 | 0.424
Table 2: Validation metrics for model and task combinations, averaged across 5
folds. Train loss is included to measure overfitting

###  4.4 Early Detection Testing

We can test out of fold early detection performance by creating an early
detection dataset, consisting of only validation flights before maintenance in
the intake gasket leak/damage class. Note that this does not account for false
positives, since they would not be taken in for maintenance. This test is
repeated for each fold using the best InceptionTime model, trained on both
tasks simultaneously. The results are summarized in Figure 2. Based on the 5
fold validation, one can conclude that there is no statistically significant
difference in recall across the number of flights before maintenance. This
suggests that the model is capable of predicting a part failure or the need
for maintenance before the problem arises in the flight immediately before
maintenance.

Figure 2: Min, mean, and max early detection accuracy of maintenance issues on
validation flights. Only out of fold pre maintenance flights are included in
this calculation, done using each fold’s best trained model based on
validation accuracy.

##  5 Discussion of Dataset Research Potential

##### Flight Safety and Predictive Maintenance Systems

The Cessna 172 is the most produced aircraft at over 43,000 units Goyer et
al., (2022) and is widely used in general aviation. This includes flight
schools, recreation, and many other applications. However, general aviation is
one of the most dangerous forms of civil aviation according to Board, (2011).
The NGAFID Aviation Maintenance dataset can be used for training and testing
systems that detect maintenance issues early to improve aircraft safety.

##### Class Imbalance Research

One of the major challenges in this dataset is the under-representation of
many classes. Half of all flights are recorded after maintenance and two
damages classes, intake gasket leak/damage and rocker cover leak/loose/damage,
represent the majority of all damage classes. This is so problematic that
models trained to predict the specific damage class for a flight will predict
that the flight is post maintenance or that the flight belongs to one of the
two majority classes more than 95% of the time. Future researchers can use
this dataset to test methods that mitigate problems caused by imbalanced
training data. This is particularly relevant for prognostic health management,
where it is often easy to gather data of a system in normal operation, but
extremely difficult to gather abnormal operation data.

##### Contrastive Representation Learning

Contrastive representation learning is a technique that attempts to learn a
representation of the input data in order to provide a benefit to a downstream
task, as described by Wang and Isola, (2020). While some flights were excluded
from the subset of data used for benchmarking, those flights could be used to
learn representations of flight data, which may be useful for improving
performance on the downstream detection and classification task.

##### Time Series Classification Benchmarking

Benchmark results suggest that the NGAFID Aviation Maintenance dataset is
particularly challenging for non deep learning methods. Of the two main tasks,
deep learning models seem to perform reasonable well for maintenance issue
detection, but all models perform quite poorly on maintenance issue
classification. One possible explanation is that the data is more similar to
audio data, where deep methods perform well, and less similar to the time
series data benchmarks by Bagnall et al., (2018), where methods such as Mini
Rocket and HIVECOTE perform well. While we were unable to test HIVECOTE, we
invite other researchers to evaluate their methods on the NGAFID Aviation
Maintenance dataset.

##### Anomaly Detection

While this dataset does not contain localized annotations of failures, one
could separate flights into regular operation and compromised operation and
attempt anomaly detection. As noted in earlier sections, the authors of this
paper believe that anomaly detection in this dataset would be extremely
difficult. This is because the data collected contains a significant degree of
variance caused by other factors, such as the pilot’s experience, the weather,
the flight plan, the payload carried, and so on.

Furthermore, anomaly detection works best in cases of acute part failure,
where the failure of a part is dramatic or has a dramatic impact on the
overall system. As noted in earlier sections, chronic part failure is more
common and does not create a dramatic impact on the overall system. While such
failures make the aircraft unsafe to fly, it does not immediately cause the
aircraft to crash or cease operation.

An auto encoder system, such as An and Cho, (2015), would struggle with this
dataset. First, the encoder must accommodate significant degrees of variance
when there are no anomalies. Second, the encoder must detect the most subtle
of changes in flight characteristics, which may be easily overshadowed by
other sources of variance. Imagine a pilot performing a barrel roll with a
leaky gasket; the variance in data caused by the barrel roll would eclipse any
variance caused by a leaky gasket.

##### Transfer Learning

This dataset also presents opportunities for transfer learning, Torrey and
Shavlik, (2010). While this dataset uses sensor data from the Cessna 172
collected from a flight school, the same sensors can be placed on other Cessna
172 aircraft for other uses. For example, the Cessna 172 can be used for
passenger, cargo, or military purposes. It is also possible that these sensors
can be mounted in similar single engine aircraft. This would allow researchers
to develop PHM systems without needing to collect as much data from the other
aircraft; the data from the other aircraft can be supplemented with the NGAFID
data.

The authors believe that any transfer learning for a military purpose will not
directly endanger more lives. This is because this dataset only contains
flight and maintenance information, which can be used to improve the safety
and reliability of systems.

##  6 Future Research

##### Flight Event Detection Datasets

Currently, the NGAFID web application provides detection services for certain
flight events, such as a aircraft stall event. Using the existing flight data
and expert labeling of said flight data, it is possible to create a time
series localization dataset, with the goal of detecting both the presence and
timing of such events.

##### Unsupervised and Self Supervised Datasets

The NGAFID database contains more than 900,000 hours of flight data from
various fleets and aircraft models. While most of the data is unlabeled, it
can still be used for self supervised and unsupervised learning, such as
contrastive representation learning or masked data modeling, similar to Devlin
et al., (2018).

##### Dataset Expansion

The authors are actively working with the Federal Aviation Administration and
existing flights schools to obtain more maintenance data to expand this
dataset to multiple airframes and fleets of aircraft. This process is both
slow and costly due to data governance and legal issues. Given the value of
this dataset even from a single fleet, the authors will provide the current
available data for the Cessna 172 and work on a future data release to include
additional air frames and maintenance issues.

##  7 Conclusion

In this paper we present 31,177 hours of flight data across 28,935 flights,
which occur relative to 2,111 unplanned maintenance events clustered into 36
types of maintenance issues. Each flight records information from 23 different
sensors every second on the Cessna 172 aircraft, during normal operation of a
flight school. Our paper makes the significant contribution of providing non-
simulated, compromised aircraft flight data, collected ethically at no
additional danger to the pilots involved. This dataset is made easily
accessible at the links in Section 1.

The large amount of flight data involving a compromised aircraft is
particularly valuable to prognostic health management and predictive
maintenance. Because aircraft in question, the Cessna 172, is often used in
flight schools, recreation, agriculture, and more, this dataset can help
create systems that can greatly improve flight safety.

Finally, the NGAFID Aviation Maintenance dataset will be of particular
interest to machine learning researchers working with time series data. It is
our aim, by releasing this dataset and identifying areas for future research,
to encourage further work in the detection and classification of maintenance
issues. We hope this in turn leads to improved future detection and
classification algorithms.

## References

  * Akhbardeh et al., (2020) Akhbardeh, F., Desell, T., and Zampieri, M. (2020).  MaintNet: A collaborative open-source library for predictive maintenance language resources.  In Proceedings of the 28th International Conference on Computational Linguistics: System Demonstrations, pages 7–11, Barcelona, Spain (Online). International Committee on Computational Linguistics (ICCL).
  * An and Cho, (2015) An, J. and Cho, S. (2015).  Variational autoencoder based anomaly detection using reconstruction probability.  Special Lecture on IE, 2(1):1–18.
  * Arena et al., (2021) Arena, F., Collotta, M., Luca, L., Ruggieri, M., and Termine, F. G. (2021).  Predictive maintenance in the automotive sector: A literature review.  Mathematical and Computational Applications, 27(1):2.
  * Argyriou et al., (2006) Argyriou, A., Evgeniou, T., and Pontil, M. (2006).  Multi-task feature learning.  Advances in neural information processing systems, 19.
  * Arias Chao et al., (2021) Arias Chao, M., Kulkarni, C., Goebel, K., and Fink, O. (2021).  Aircraft engine run-to-failure dataset under real flight conditions for prognostics and diagnostics.  Data, 6(1):5.
  * Assaf et al., (2019) Assaf, R., Giurgiu, I., Bagehorn, F., and Schumann, A. (2019).  Mtex-cnn: Multivariate time series explanations for predictions with convolutional neural networks.  In 2019 IEEE International Conference on Data Mining (ICDM), pages 952–957. IEEE.
  * Bagnall et al., (2018) Bagnall, A., Dau, H. A., Lines, J., Flynn, M., Large, J., Bostrom, A., Southam, P., and Keogh, E. (2018).  The uea multivariate time series classification archive, 2018.  arXiv preprint arXiv:1811.00075.
  * Board, (2011) Board (2011).  Review of U.S. civil aviation accidents: review of aircraft accident data, 2007-2009.  U.S. National Transportation Safety Board.
  * Celikmih et al., (2020) Celikmih, K., Inan, O., and Uguz, H. (2020).  Failure prediction of aircraft equipment using machine learning with a hybrid data preparation method.  Scientific Programming, 2020.
  * Chu et al., (2010) Chu, E., Gorinevsky, D., and Boyd, S. (2010).  Detecting aircraft performance anomalies from cruise flight data.  In AIAA Infotech@ Aerospace 2010, page 3307.
  * Dangut et al., (2022) Dangut, M. D., Jennions, I. K., King, S., and Skaf, Z. (2022).  A rare failure detection model for aircraft predictive maintenance using a deep hybrid learning approach.  Neural Computing and Applications, pages 1–19.
  * Dempster et al., (2020) Dempster, A., Petitjean, F., and Webb, G. I. (2020).  Rocket: exceptionally fast and accurate time series classification using random convolutional kernels.  Data Mining and Knowledge Discovery, 34(5):1454–1495.
  * Dempster et al., (2021) Dempster, A., Schmidt, D. F., and Webb, G. I. (2021).  Minirocket: A very fast (almost) deterministic transform for time series classification.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pages 248–257.
  * Devlin et al., (2018) Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K. (2018).  Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805.
  * Dosovitskiy et al., (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al. (2020).  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv preprint arXiv:2010.11929.
  * Fawaz et al., (2019) Fawaz, H. I., Forestier, G., Weber, J., Idoumghar, L., and Muller, P.-A. (2019).  Deep learning for time series classification: a review.  Data mining and knowledge discovery, 33(4):917–963.
  * Fawaz et al., (2020) Fawaz, H. I., Lucas, B., Forestier, G., Pelletier, C., Schmidt, D. F., Weber, J., Webb, G. I., Idoumghar, L., Muller, P.-A., and Petitjean, F. (2020).  Inceptiontime: Finding alexnet for time series classification.  Data Mining and Knowledge Discovery, 34(6):1936–1962.
  * Goldman et al., (2002) Goldman, S. M., Fiedler, E. R., and King, R. E. (2002).  General aviation maintenance-related accidents: A review of ten years of ntsb data.
  * Goyer et al., (2022) Goyer, B. I., Staff, F., Staff, I. G., and Flying (2022).  Cessna 172: Still relevant today?
  * He et al., (2016) He, K., Zhang, X., Ren, S., and Sun, J. (2016).  Identity mappings in deep residual networks.  In European conference on computer vision, pages 630–645. Springer.
  * Johnson and Khoshgoftaar, (2019) Johnson, J. M. and Khoshgoftaar, T. M. (2019).  Survey on deep learning with class imbalance.  Journal of Big Data, 6(1):1–54.
  * Karboviak et al., (2018) Karboviak, K., Clachar, S., Desell, T., Dusenbury, M., Hedrick, W., Higgins, J., Walberg, J., and Wild, B. (2018).  Classifying aircraft approach type in the national general aviation flight information database.  In International Conference on Computational Science, pages 456–469. Springer.
  * Karim et al., (2017) Karim, F., Majumdar, S., Darabi, H., and Chen, S. (2017).  Lstm fully convolutional networks for time series classification.  IEEE access, 6:1662–1669.
  * Khadilkar and Balakrishnan, (2012) Khadilkar, H. and Balakrishnan, H. (2012).  Estimation of aircraft taxi fuel burn using flight data recorder archives.  Transportation Research Part D: Transport and Environment, 17(7):532–537.
  * Klein, (1989) Klein, V. (1989).  Estimation of aircraft aerodynamic parameters from flight data.  Progress in Aerospace Sciences, 26(1):1–77.
  * Kukačka et al., (2017) Kukačka, J., Golkov, V., and Cremers, D. (2017).  Regularization for deep learning: A taxonomy.  arXiv preprint arXiv:1710.10686.
  * Le and Zuidema, (2016) Le, P. and Zuidema, W. (2016).  Quantifying the vanishing gradient and long distance dependency problem in recursive neural networks and recursive lstms.  arXiv preprint arXiv:1603.00423.
  * Lines et al., (2018) Lines, J., Taylor, S., and Bagnall, A. (2018).  Time series classification with hive-cote: The hierarchical vote collective of transformation-based ensembles.  ACM Transactions on Knowledge Discovery from Data, 12(5).
  * Liu et al., (2012) Liu, Y., Frederick, D. K., DeCastro, J. A., Litt, J. S., and Chan, W. W. (2012).  User’s guide for the commercial modular aero-propulsion system simulation (c-mapss): Version 2.  Technical report.
  * Löning et al., (2022) Löning, M., Bagnall, T., Király, F., Middlehurst, M., Ganesh, S., Oastler, G., Lines, J., ViktorKaz, Walter, M., Mentel, L., RNKuhns, chrisholder, Tsaprounis, L., Owoseni, T., Rockenschaub, P., Khrapov, S., jesellier, danbartl, Bulatova, G., eenticott shell, Lovkush, Take, K., Meyer, S. M., AidenRushbrooke, Gilbert, C., Schäfer, P., oleskiewicz, Xu, Y.-X., Ansari, A., and Sakshi, A. (2022).  alan-turing-institute/sktime: v0.11.4.
  * Middlehurst et al., (2021) Middlehurst, M., Large, J., Flynn, M., Lines, J., Bostrom, A., and Bagnall, A. (2021).  Hive-cote 2.0: a new meta ensemble for time series classification.  Machine Learning, 110(11):3211–3243.
  * Oguiza, (2022) Oguiza, I. (2022).  tsai - a state-of-the-art deep learning library for time series and sequential data.  Github.
  * Orsenigo and Vercellis, (2010) Orsenigo, C. and Vercellis, C. (2010).  Combining discrete svm and fixed cardinality warping distances for multivariate time series classification.  Pattern Recognition, 43(11):3787–3794.
  * Randle et al., (2011) Randle, W. E., Hall, C. A., and Vera-Morales, M. (2011).  Improved range equation based on aircraft flight data.  Journal of Aircraft, 48(4):1291–1298.
  * Rezaeianjouybari and Shang, (2020) Rezaeianjouybari, B. and Shang, Y. (2020).  Deep learning for prognostics and health management: State of the art, challenges, and opportunities.  Measurement, 163:107929.
  * Seto et al., (2015) Seto, S., Zhang, W., and Zhou, Y. (2015).  Multivariate time series classification using dynamic time warping template selection for human activity recognition.  In 2015 IEEE Symposium Series on Computational Intelligence, pages 1399–1406. IEEE.
  * Shifaz et al., (2020) Shifaz, A., Pelletier, C., Petitjean, F., and Webb, G. I. (2020).  Ts-chief: a scalable and accurate forest algorithm for time series classification.  Data Mining and Knowledge Discovery, 34(3):742–775.
  * Torrey and Shavlik, (2010) Torrey, L. and Shavlik, J. (2010).  Transfer learning.  In Handbook of research on machine learning applications and trends: algorithms, methods, and techniques, pages 242–264. IGI global.
  * Tsui et al., (2015) Tsui, K. L., Chen, N., Zhou, Q., Hai, Y., and Wang, W. (2015).  Prognostics and health management: A review on data driven approaches.  Mathematical Problems in Engineering, 2015.
  * Vincent et al., (2013) Vincent, E., Barker, J., Watanabe, S., Le Roux, J., Nesta, F., and Matassoni, M. (2013).  The second ‘chime’ speech separation and recognition challenge: Datasets, tasks and baselines.  In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 126–130.
  * Wang and Isola, (2020) Wang, T. and Isola, P. (2020).  Understanding contrastive representation learning through alignment and uniformity on the hypersphere.  In III, H. D. and Singh, A., editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 9929–9939. PMLR.
  * Wang et al., (2017) Wang, Z., Yan, W., and Oates, T. (2017).  Time series classification from scratch with deep neural networks: A strong baseline.  In 2017 International joint conference on neural networks (IJCNN), pages 1578–1585. IEEE.
  * Yang et al., (2021) Yang, H., LaBella, A., and Desell, T. (2021).  Predictive maintenance for general aviation using convolutional transformers.  arXiv preprint arXiv:2110.03757.
  * Yildirim and Kurt, (2018) Yildirim, M. T. and Kurt, B. (2018).  Aircraft gas turbine engine health monitoring system by real flight data.  International Journal of Aerospace Engineering, 2018.

##  Appendix A Appendix

Class Name | Overall Count of Flights | Subset Count of Flights
---|---|---
aircraft start/external issue | 90 | 0
baffle bracket loose/damage | 35 | 0
baffle crack/damage/loose/miss | 545 | 304
baffle mount loose/damage | 42 | 0
baffle plug need repair/replace | 465 | 254
baffle rivet loose/miss/damage | 92 | 0
baffle screw miss/loose | 348 | 211
baffle seal loose/damage | 336 | 197
baffle spring damage | 79 | 0
baffle tie/tie rod loose or damage | 303 | 0
cowling miss/loose/damage | 89 | 0
cylinder compression issue | 302 | 143
cylinder crack/fail/need part repair | 196 | 108
cylinder exhaust valve/stuck valve issue | 48 | 0
cylinder head/exhaust gas temperature issue | 71 | 0
cylinder/exhaust push rod/tube damage | 106 | 0
drain line/tube damage | 127 | 0
engine crankcase/crankshaft/firewall near repair | 99 | 0
engine failure/fire/time out | 236 | 161
engine idle/rpm issue | 150 | 93
engine need repair/reinstall/clean | 148 | 80
engine run rough | 311 | 141
engine seal/tube/bolt loose or damage | 144 | 76
engine/propeller overspeed or damage | 137 | 94
induction damage/hardware fail | 51 | 0
intake gasket leak/damage | 4244 | 2098
intake tube/bolt/seal/boot loose or damage | 556 | 269
magneto failure | 51 | 0
mixture fail/need adjust | 17 | 0
oil cooler need maintenance | 123 | 75
oil dipstick/tube need repair | 20 | 0
oil leak/pressure issue | 53 | 0
oil return line issue | 37 | 0
pilot/in-flight noticed issue | 317 | 75
rocker cover leak/loose/damage | 2157 | 1024
spark plug need repair/replace | 15 | 0
flights recorded after maintenance | 10291 | 5844
flights recorded during maintenance | 6504 | 0
Table 3: Count of flights by class in the full dataset and the subset used for
the benchmark experiments.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes] Seciton 3

    3. (c)

Did you discuss any potential negative societal impacts of your work? [No] The
data is collected from purpose built sensors on aircraft. We believe that this
data is so specialized as to be only applicable in the aviation industry and
only to single engine aircraft, which very few people interact with regularly.

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [N/A]

    2. (b)

Did you include complete proofs of all theoretical results? [N/A]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]
As a repository with links to Colab Notebooks

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes] See the code in the repo

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes] But only for the early detection in section
4.4. The overall results in 4.3 only report mean values. However, you can
reproduce results as the splits are availalbe in the dataset.

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See section
4.2.5

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes]

    2. (b)

Did you mention the license of the assets? [Yes]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [Yes]

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [No] See data sheet supplementary material

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [Yes] Section 3.3

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]