# Graph Neural Network-Based Anomaly Detection for River Network Systems

Katie Buchhorn  Queensland University of Technology, School of Mathematical
Sciences  Queensland University of Technology, Centre for Data Science  Edgar
Santos-Fernandez  Queensland University of Technology, School of Mathematical
Sciences  Queensland University of Technology, Centre for Data Science  Kerrie
Mengersen  Queensland University of Technology, School of Mathematical
Sciences  Queensland University of Technology, Centre for Data Science  Robert
Salomone  Queensland University of Technology, School of Mathematical Sciences
Queensland University of Technology, School of Computer Science

###### Abstract

Water is the lifeblood of river networks, and its quality plays a crucial role
in sustaining both aquatic ecosystems and human societies. Real-time
monitoring of water quality is increasingly reliant on in-situ sensor
technology. Anomaly detection is crucial for identifying erroneous patterns in
sensor data, but can be a challenging task due to the complexity and
variability of the data, even under normal conditions. This paper presents a
solution to the challenging task of anomaly detection for river network sensor
data, which is essential for accurate and continuous monitoring. We use a
graph neural network model, the recently proposed Graph Deviation Network
(GDN), which employs graph attention-based forecasting to capture the complex
spatio-temporal relationships between sensors. We propose an alternate anomaly
threshold criteria for the model, GDN+, based on the learned graph. To
evaluate the model’s efficacy, we introduce new benchmarking simulation
experiments with highly-sophisticated dependency structures and subsequence
anomalies of various types. We further examine the strengths and weaknesses of
this baseline approach, GDN, in comparison to other benchmarking methods on
complex real-world river network data. Findings suggest that GDN+ outperforms
the baseline approach in high-dimensional data, while also providing improved
interpretability. We also introduce software called gnnad.

## Keywords

Anomaly Detection, Graph Deviation Network, Graph Neural Network, Multivariate
Time Series, Graph Attention Forecasting, Spatio-temporal Data, Complex
Systems

## Introduction

River network systems play a vital role as freshwater habitats for aquatic
life, and as support for terrestrial ecosystems in riparian zones, but are
particularly sensitive to the anthropogenic impacts of climate change, water
pollution and over-exploitation, among other factors. As a United Nations
Sustainable Development Goal [1], water quality is a major environmental
concern worldwide. The use of in-situ111in-situ refers to an instrument in
direct contact with the medium of observation. sensors for data collection on
river networks is increasingly prevalent [2, 3], generating large amounts of
data that allow for the identification of fine-scale spatial and temporal
patterns, trends, and extremes, as well as potential sources of pollutants and
their downstream impacts. However, such sensors are susceptible to technical
errors relating to the equipment, herein defined as anomalies, for example due
to miscalibration, biofouling, electrical interference and battery failure. In
contrast, extreme events in rivers occur as result of heavy rain and floods.
Technical anomalies must be identified before the data are considered for
further analysis, as they can introduce bias in model parameters and affect
the validity of statistical inferences, which then confounds the
identification of true changes in water variables. Trustworthy data is needed
to produce reliable and accurate assessments of water quality, for enhanced
environmental monitoring, and for guiding management decisions in the
prioritisation of ecosystem health.

Anomaly detection in river networks is challenging due to the highly dynamic
nature of river water even under typical conditions [4], as well as the
complex spatial relationships between sensors. The unique spatial
relationships between neighbouring sensors on a river network are
characterised by a branching network topology with flow direction and
connectivity, embedded within the 3-D terrestrial landscape. Common anomalies
from data obtained from in-situ sensors are generally characterised by
multiple consecutive observations (subsequence or persistent, [5]), including
sensor drift and periods of unusually high or low variability, which may
indicate the necessity for sensor maintenance or calibration [6]. Such
anomalies are difficult to detect and often associated with high false
negative rates [7].

Earlier statistical studies have focused on developing autocovariance models
based on within-river relationships to capture the unique spatial
characteristics of rivers [8]. Although these methods are adaptable to
different climate zones, and have recently been extended to take temporal
dependencies into account [9], data sets generated by in-situ sensors still
pose significant computational challenges with such prediction methods due to
the sheer volume of data [10]. Autocovariance matrices must be inverted when
fitting spatio-temporal models and making predictions, and the distances
between sites must be known. Previous work [11], aimed to detect drift and
high-variability anomalies in water quality variables, by studying a range of
neural networks calibrated using a Bayesian multi-objective optimization
procedure. However, the study was limited to analyzing univariate time series
data, and the supervised methods required a significant amount of labeled data
for training, which are not always available.

There are limited unsupervised anomaly detection methods for subsequence
anomalies (of variable length), and even less so for multivariate time series
anomaly detection [5]. One such method uses dynamic clustering on learned
segmented windows to identify global and local anomalies [12]. However, this
algorithm requires the time series to be well-aligned, and is not suitable for
the lagged temporal relationships observed with river flow. Another method to
detect variable-length subsequence anomalies in multivariate time series data
uses dimensionality reduction to construct a one-dimensional feature, to
represent the density of a local region in the recurrence representation,
indicating the recurrence of patterns obtained by a sliding window [13]. A
similarity measure is used to classify subsequences as either non-anomalous or
anomalous. Only two summary statistics were used in this similarity measure
and the results were limited to a low-dimensional simulation study. In
contrast, the technique introduced by [14], DeepAnT, uses a deep convolutional
neural network (CNN) to predict one step ahead. This approach uses Euclidean
distance of the forecast errors as the anomaly score. However, an anomaly
threshold parameter must be provided.

Despite the above initial advances, challenges still remain in detecting
persistent variable-length anomalies within high-dimensional data exhibiting
noisy and complex spatial and temporal dependencies. With the aim of
addressing such challenges, we explore the application of the recently-
proposed Graph Deviation Network (GDN) [15], and explore refinements with
respect to anomaly scoring that specifically address the needs of
environmental monitoring. The GDN approach [15] is a state-of-the-art model
that uses sensor embeddings to capture inter-sensor relationships as a learned
graph, and employs graph attention-based forecasting to predict future sensor
behaviour. Anomalies are flagged when the error scores are above a calculated
threshold value. By learning the interdependencies among variables and
predicting based on the typical patterns of the system in a semi-supervised
manner, this approach is able to detect deviations when the expected spatial
dependencies are disrupted. As such, GDN offers the ability to detect even the
small-deviation anomalies generally overlooked by other distance based and
density based anomaly detection methods for time series [16, 17], while
offering robustness to lagged variable relationships. Unlike the commonly-used
statistical methods that explicitly model covariance as a function of
distance, GDN is flexible in capturing complex variable relationships
independent of distance. GDN is also a semi-supervised approach, eliminating
the need to label large amounts of data, and offers a computationally
efficient solution to handle the ever increasing supply of sensor data.
Despite the existing suite of methods developed for anomaly detection, only a
limited number of corresponding software packages are available to
practitioners. In summary, we have identified the following gaps in the
current literature and research on this topic:

  1. 1.

An urgent need exists for a flexible approach that can effectively capture
complex spatial relationships in river networks without the specification of
an autocovariance model, and the ability to learn from limited labeled data,
in a computationally efficient manner.

  2. 2.

Lack of data and anomaly generation schemes on which to benchmark methods,
that exhibit complex spatial and temporal dependencies, as observed across
river networks.

  3. 3.

Lack of open-source software for anomaly detection, which hinders the
accessibility and reproducibility of research in this field, and limits the
ability for individuals and organizations with limited resources to implement
effective anomaly detection strategies.

Our work makes four primary contributions to the field:

  1. 1.

An improvement of the GDN approach via the threshold calculation based on the
learned graph is presented, and shown to detect anomalies more accurately than
GDN while improving the ability to locate anomalies across a network.

  2. 2.

New simulation data generation methods for highly-sophisticated spatio-
temporal structures, contaminated with various types of persistent anomalies
are provided as benchmarking data.

  3. 3.

Experimentation on benchmark data, as well as a real-world river network
dataset, is conducted to evaluate the strengths and limitations of GDN (and
its variants) in increasingly challenging settings.

  4. 4.

User-friendly, free open-source software for the GDN/GDN+ approach is
available on the pip repository as gnnad, with data and anomaly generation
modules.

The structure of the remainder of the paper is as follows: The next section
details the methods of GDN and the model extension GDN+, and describes the
methodology of the simulated data and anomaly generation. In the Results
section we present an extensive simulation study on the benchmarking data, as
well as a real-world case study. In each case, the performance of GDN/GDN+ is
assessed against other state-of-the-art anomaly detection models. Further
details and example code for the newly-released software are also provided.
The paper concludes with a discussion of the findings, and the strengths and
weaknesses of this model.

## Methods

Consider multivariate time series data
$\mathrm{Y}=\left[\boldsymbol{y}^{(1)},\ldots,\boldsymbol{y}^{(T)}\right]$,
obtained from $n$ sensors over $T$ time ticks. The (univariate) data collected
from sensor $i=1,\ldots,n$ at time $t=1,\ldots,T$ are denoted as
$y_{i}^{(t)}$. Following the standard semi-supervised anomaly detection
approach [18, 19], non-anomalous data are used for training, while the test
set may contain anomalous data. That is, we aim to learn the sensor behaviour
using data obtained under standard operational conditions throughout the
training phase and identify anomalous sensor readings during testing, as those
which deviate substantially from the learned behaviour. As the algorithm
output, each test point $\boldsymbol{y}^{(t)}$ is assigned a binary label,
$a(t)\in\\{0,1\\}$, where $a(t)=1$ indicates an anomaly at time $t$, anywhere
across the full sensor network.

The GDN approach [15] for anomaly detection is composed of two aspects:

  1. 1.

Forecasting-based time series model: a non-linear autoregressive multivariate
time series model that involves graph neural networks is trained, and

  2. 2.

Threshold-based anomaly detection: transformations of the individual
forecasting errors are used to determine if an anomaly has occurred, if such
errors exceed a calculated threshold.

The above components are described in ore detail below.

Forecasting-based Time Series Model. To predict $\boldsymbol{y}^{(t)}$, the
model takes as input $p\in\mdmathbb{N}$ lags of the multivariate series,

$\mathrm{X}^{(t)}:=\left[\boldsymbol{y}^{(t-1)},\ldots,\boldsymbol{y}^{(t-p)}\right].$

The $i$-th row (containing sensor $i$’s measurements for the previous $p$
lags) of the above input matrix is represented by the column vector,
$\boldsymbol{x}_{i}^{(t)}=(y_{i}^{(t-1)},\ldots,y_{i}^{(t-p)})$. Prior to
training, the practitioner specifies acceptable candidate relationships via
the sets $\mathcal{C}_{1},\ldots,\mathcal{C}_{n}$, where each
$\mathcal{C}_{i}\subseteq\\{1,2,\ldots,n\\}$ and does not contain $i$. These
sets specify which nodes are allowed to be considered to be connected from
node $i$ (noting that the adjacency graph connections are not necessarily
symmetric).

The model implicitly learns a graph structure via training sensor embedding
parameters $\boldsymbol{v}_{i}\in\mdmathbb{R}^{d}$ for $i=1,\ldots,n$ which
are used to construct a graph. The intuition is that the embedding vectors
capture the inherent characteristics of each sensor, and that sensors which
are “similar” in terms of the angle between their vector embeddings are
considered connected. Formally, the quantity $e_{ji}$ is defined as the cosine
similarity between the vector embeddings of sensors $i$ and $j$:

$e_{ji}=\frac{\boldsymbol{v}_{i}^{\top}\boldsymbol{v}_{j}}{||\boldsymbol{v}_{i}||\,||\boldsymbol{v}_{j}||}\mdmathbb{I}\\{j\in{\mathcal{C}}_{i}\\},\quad
i,j\in\\{1,\ldots,n\\},$

with $||\cdot||$ denoting the Euclidean norm, and indicator function
$\mdmathbb{I}$ which equals 1 when node $j$ belongs to set
${\mathcal{C}}_{i}$, and 0 otherwise. Note that the similarity is forced to be
zero if a connecting node is not in the permissible candidate set. Next, let
$e_{j,(i)}$ be the $i$-th largest value in $(e_{j1},\ldots,e_{jn})$. A graph-
adjacency matrix (and in turn a graph itself) is then constructed from the
sensor similarities via:

$A_{ji}=\mdmathbb{I}\\{\\{e_{ji}\geq e_{{j,(K)}}\\}\cup\\{i=j\\}\\},$

for user-specified $K\in\\{1,\ldots,n\\}$ which determines the maximum number
of edges from a node, referred to as the “Top-K" hyperparameter.

The above describes how the trainable parameters
$\\{\boldsymbol{v}_{k}\\}_{k=1}^{n}$ yield a graph. Next, the lagged series
are fed individually through a shallow Graph Attention Network [20] that uses
the previously constructed graph. Here, each node corresponds to a sensor, and
the node features for node $i$ are the lagged (univariate) time-series values,
$\boldsymbol{x}_{i}^{(t)}\in\mdmathbb{R}^{p}$. Allow a parameter weight matrix
$\mathrm{W}\in\mdmathbb{R}^{d\times p}$ to apply a shared linear transform to
each node. Then, the output of the network is given by

$\boldsymbol{z}_{i}^{(t)}=\max\left\\{\boldsymbol{0},\left(\sum_{j:A_{ji}>0}{}_{ij}\mathrm{W}\boldsymbol{x}_{i}^{(t)}\right)\right\\},$

where $\boldsymbol{z}_{i}^{(t)}$ is called the node representation, and
coefficients ij are the attention paid to node $j$ when computing the
representation for node $i$, with:

ij $\displaystyle={\rm
LeakyReLU}\left(\boldsymbol{a}^{\top}\left(\boldsymbol{v}_{i}\oplus\mathrm{W}\boldsymbol{x}_{i}^{(t)}+\boldsymbol{v}_{j}\oplus\mathrm{W}\boldsymbol{x}_{j}^{(t)}\right)\right),$
(1) $\displaystyle\text{where
}{}_{ij}=\frac{\exp({}_{ij})}{\sum_{k:A_{ki}>0}\exp({}_{ik})}.$

with learnable parameters $\boldsymbol{a}\in\mdmathbb{R}^{2d}$, where $\oplus$
denotes concatenation, and ${\rm
LeakyReLU}(\boldsymbol{x}):=\max\\{\delta\boldsymbol{x},\boldsymbol{x}\\}$ for
$\delta>0$, with the maximum operation applied elementwise. Note the
addition222It seems that the authors of [15] mistakenly denote this addition
as concatenation in the paper, however, their corresponding reference code
computes addition. in Equation (1) and that $\sum_{j=1}^{n}{}_{ij}=1$.
Intuitively, the above is an automated mechanism to aggregate information from
a node itself and neighbouring nodes (whilst simultaneously assigning a weight
of how much information to take from each neighbour) to compute a vector
representing extracted information about node $i$ itself and its neighbours’
interaction with it. The final model output (prediction) is given by,

$\widehat{\boldsymbol{y}}^{(t)}=f_{\boldsymbol{\eta}}\left(\left[\boldsymbol{v}_{1}\odot\boldsymbol{z}_{1}^{(t)},\ldots,\boldsymbol{v}_{n}\odot\boldsymbol{z}_{n}^{(t)}\right]\right),$

where $f_{\boldsymbol{\eta}}:\mdmathbb{R}^{d\times n}\to\mdmathbb{R}^{n}$ is a
feedforward neural network with parameters , and $\odot$ denotes element-wise
multiplication. The model is trained by optimizing the parameters
$\\{\boldsymbol{v}_{i}\\}_{i=1}^{n}$, $\mathrm{W}$, $\boldsymbol{a}$, and to
minimize the mean squared error loss function

${\mathcal{L}}=\frac{1}{T-p}\sum_{t=p+1}^{T}\bigg{|}\bigg{|}\widehat{\boldsymbol{y}}^{(t)}-\boldsymbol{y}^{(t)}\bigg{|}\bigg{|}^{2}.$

Threshold-based Anomaly Detection. Given the learned inter-sensor and temporal
relationships, we are able to detect anomalies as those which deviate from
these interdependencies. An anomalousness score is computed for each time
point in the test data. For each sensor $i$, we denote the prediction error at
time, $t$, as,

${}_{i,t}=|y_{i}^{(t)}-\widehat{y}_{i}^{(t)}|,$

with $|\cdot|$ denoting the absolute value, and the vector of prediction error
for each sensor is, $\boldsymbol{\epsilon}_{i}\in\mdmathbb{R}^{T-p}$. Since
the error values of different sensors may vary substantially, we perform a
robust normalisation of each sensor’s errors to prevent any one sensor from
overly dominating the others, that is,

$\tilde{\boldsymbol{\epsilon}}_{i}=\left(\frac{\boldsymbol{\epsilon}_{i}-{\rm
Median}(\boldsymbol{\epsilon}_{i})}{{\rm
IQR}(\boldsymbol{\epsilon}_{i})}\right),$

where IQR denotes inter-quartile range. In the original work by [15], a time
point $t$ is flagged as anomalous if,

$A(t)=\mdmathbb{I}\\{\max_{i}(\tilde{\epsilon}_{i,t})>\kappa\\},$

using the notation $\tilde{\epsilon}_{i,t}$ for the error value at the $t$-th
index. Alternatively, the authors recommend using a simple moving average of
$\tilde{\boldsymbol{\epsilon}}_{i}$ and flagging time $t$ as anomalous if the
maximum of that moving average exceeds . The authors specify as the maximum of
the normalised errors observed on some (non-anomalous) validation data,
denoted by variant epsilon, $\tilde{\boldsymbol{\varepsilon}}$. However, this
is applied to all sensors as a fixed threshold.

### Sensor-based Anomaly Threshold, GDN+

The behavior of water quality variables may differ across space, for example,
water level at high-altitude river network branches are generally
characterised by rainfall patterns, whereas water level downstream near a
river outlet can also be influenced by tidal patterns. A fixed threshold
across the network does not allow for any local adaptations in error
sensitivity. For this reason, this work also considers the novel sensor-
specific threshold calculation,

$A_{i}(t)=\mdmathbb{I}\\{\tilde{\epsilon}_{i,t}>{}_{i}\\},$

where i is chosen such that,

$\frac{1}{|\\{\tilde{\varepsilon}_{j,t}\\}_{j:A_{ji}>0}|}\sum_{j:A_{ji}>0}\mdmathbb{I}\\{\tilde{\varepsilon}_{j,t}<{}_{i}\\}=\tau,$

for some user-specified quantile, $\tau\in(0,1)$, and where $|\cdot|$ is the
cardinality. In other words, the threshold for each sensor, i, is set as the
-th quantile of the normalised error scores across the neighbourhood of $i$,
on the validation data set. Unless otherwise stated, we set $\tau=0.99$. In
this way, the sensor threshold is based only on its direct neighbourhood, as
opposed to the original method which uses the global maximum, and is thus more
robust as a percentile and in tune with the local behaviour of the system.

We refer to the GDN model using this variant of the threshold-based anomaly
detection as GDN+.

### New Class of Benchmarking Data

The following is a method for simulating synthetic datasets with persistent
anomalies inspired by the statistical models recently used to model river
network data [9]. Let $\mathcal{S}$ denote an arbitrary set of individual
spatial locations, with locations
$\boldsymbol{s}_{1},\boldsymbol{s}_{2},\ldots,\boldsymbol{s}_{n}\in\mathcal{S}$
chosen by experimental design [21] or otherwise. Consider a linear mixed model
with $n\times 1$ response $\boldsymbol{Y}$, and $n\times m$ design matrix
$\boldsymbol{X}$ of explanatory variables spatially indexed at locations
$\boldsymbol{s}_{1},\boldsymbol{s}_{2},\ldots,\boldsymbol{s}_{n}$,

$\boldsymbol{Y}_{t}=\boldsymbol{\beta}_{0}+\boldsymbol{X}_{t}\boldsymbol{\beta}+\boldsymbol{Z}+\boldsymbol{{}_{0}},\quad
t=1,\ldots,T,$ (2)

with time-homogeneous spatially-correlated random effects
$\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{\boldsymbol{Z}})$
and vector of independent noise terms
$\boldsymbol{{}_{0}}{\sim}\mathcal{N}(\boldsymbol{0},{}_{0}^{2}\mathrm{I})$,
yielding
$\mbox{Cov\,}[\boldsymbol{Y}_{t}|\boldsymbol{X}_{t}=\boldsymbol{x}_{t}]=\boldsymbol{\Sigma}_{\boldsymbol{Z}}+{}_{0}^{2}\mathrm{I}$,
where $\mathrm{I}$ denotes the $n\times n$ identity matrix and
$\boldsymbol{\epsilon}_{0}$ is an error term. The covariates
$\boldsymbol{X}_{t}$ for $t=1,\ldots,T$ are simulated according to an
autoregressive process based on an underlying sequence of independent random
fields,

$\boldsymbol{X}_{t}={}^{p}_{i=0}{}_{i}\widetilde{\boldsymbol{X}}_{t-i},$

where $p$ is the order of the autoregressive process, and

$\widetilde{\boldsymbol{X}}_{t}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}{\mathcal{N}}(\boldsymbol{0},\boldsymbol{\Sigma}_{\boldsymbol{X}}),\quad
t=1,\ldots,T.$

Note that other distributions may be used. Above,

$(\boldsymbol{\Sigma}_{\boldsymbol{X}})_{ij}=k(\boldsymbol{s}_{i},\boldsymbol{s}_{j}^{\prime})$

for some covariance kernel $k$. For example,

$k(\mathbf{s},\mathbf{s}^{\prime};\sigma,\alpha)={}^{2}\exp\left(-\frac{||\mathbf{s}-\mathbf{s}^{\prime}||^{2}}{\alpha}\right),$
(3)

where $||\cdot||$ denotes the Euclidean norm, ${}^{2}>0$ is the covariance-
scaling parameter, and $\alpha\in\mdmathbb{R}$ is the range parameter that
controls the rate of decay of correlation between points over distance. Figure
1 illustrates an example of a generated Gaussian random field evolving over
time.

Figure 1:  Smooth random Gaussian field used to generate covariate values,
$\boldsymbol{X}$, in the simulation studies. Examples of the field are shown
for $t=1,2,3$. Sensor locations are shown as white dots.

We consider two scenarios in generating simulated data: 1) spatial
relationships characterised by Euclidean distance only, where the covariance
matrix, $\boldsymbol{\Sigma}_{Z}$, is constructed via the kernel function
given in Equation (3), and 2) data simulated on a river network, where
$\boldsymbol{\Sigma}_{Z}$ is constructed via the kernel function given in
Equation (4). Together, this approach provides a variety of simulated data
with highly-sophisticated dependency structure, see Figure 2.

Two points $\boldsymbol{s}_{i}$, $\boldsymbol{s}_{j}$ on a river network are
said to be _flow-connected_ if they share water flow, and _flow-unconnected_
otherwise. We define stream distance,
$h_{Riv}(\boldsymbol{s}_{i},\boldsymbol{s}_{j})$, as the shortest distance
separating $\boldsymbol{s}_{i}$ and $\boldsymbol{s}_{j}$ when travelling along
a given river network. Tail-up covariance models for river networks,
introduced in [8], effectively represent spatial relationships when variables
are dominated by flow (e.g. pollutants enter a stream and only impact
downstream locations). By construction, the tail-up covariance function only
allows for correlation between flow-connected sites:

$k_{\rm
TU}(\boldsymbol{s}_{i},\boldsymbol{s}_{j};\sigma,\alpha)={}_{ij}{}^{2}\exp\left(-\frac{h_{\rm
Riv}(\boldsymbol{s}_{i},\boldsymbol{s}_{j})}{\alpha}\right)\mathcal{F}_{ij},$
(4)

where ${\mathcal{F}}_{ij}$ is equal to one if $s_{i}$ and $s_{j}$ are flow-
connected, and zero otherwise, and ij is a weighting attributed to each stream
segment to account for the upstream branching network structure and ensure
stationarity in variance (for full details, see [22, 23]). The weightings
corresponding to each segment may incorporate flow volume, or the area of the
catchment, or a proxy such as stream order [24]. Note that there are various
choices of covariance models. Tail-down models allow correlation between both
flow-connected and flow-unconnected locations, and may be more suitable for
water variables such as temperature, or organisms that can move both upstream
and downstream [25]. See further covariance model examples [8], here we use
the exponential function for decay.

Once the base multivariate time series is constructed as above, it is modified
to include river-network type anomalies as follows. Hyperparameters are,
$n_{\rm anomaly}\geq 0$, the number of subsequence anomalies and, ${}_{\rm
anomaly}>0$, the average length of an anomalous subsequence, for each anomaly
type. In this example, we consider two types of anomalies, high-variability
and drift, see Algorithm 1.

input : Time series data
$\mathrm{Y}=\left[\boldsymbol{y}^{(1)},\ldots,\boldsymbol{y}^{(T)}\right]$,
locations
$\boldsymbol{s}_{1},\boldsymbol{s}_{2},\ldots,\boldsymbol{s}_{n}\in\mathcal{S}$,
expected length of each anomaly type, anomaly, and the number of anomalies,
$n_{\rm anomaly}$.

for  _${\rm anomaly}\in\\{{\rm drift},{\rm variability}\\}$_ do

for  _$i=1,\ldots,n_{\rm anomaly}$_ do

Draw location $S\sim{\rm Uniform}(\\{1,2,\ldots,n\\})$

Draw time $t\sim{\rm Uniform}(\\{1,\ldots,T\\})$

Draw length $L\sim{\rm Poisson}({}_{\rm anomaly})$

if  _${\rm anomaly}={\rm drift}$_ then

$\boldsymbol{v}\leftarrow(\alpha,2\alpha,\ldots,L\alpha)$

// drift

else

$\boldsymbol{v}\sim\mathcal{N}(\boldsymbol{0},{}^{2}\mathrm{I}_{L\times L})$

// variability

$\boldsymbol{y}_{S,t:(t+L)}\leftarrow\boldsymbol{y}_{S,t:(t+L)}+\boldsymbol{v}$

Algorithm 1 Two-type Anomaly Generation

Figure 2:  SimEuc site locations across space (left) and SimRiver site
locations along a river network (right). Both simulations use the same (x, y)
coordinates for sensor locations. The direction of flow is from top to bottom
for SimRiver. Red dots indicate sites for which time series have been shown in
Figure 3 and Figure 4.

### Python Package

The Python package gnnad introduced in this paper extends and generalises the
research code originally implemented by [15], which is incompatible with newer
package dependencies and offers only a command line interface. The code is
refactored to be modular and user-friendly, with a scikit-inspired interface,
and extended to include visualisation, data and anomaly generation modules, as
well as the GDN+ model extension. A continuous integration/continuous
deployment (CI/CP) pipeline is established with unit testing to ensure that
changes to the code are tested and deployed efficiently. Comprehensive
documentation now accompanying the codebase enhances readability for future
developers, facilitating maintenance, reuse, and modification. Furthermore,
rigorous error handling is implemented to improve the software experience. The
software developments have resulted in a more robust, user-friendly and easily
distributable package that is available via https://github.com/KatieBuc/gnnad
and the pip repository, gnnad. See below for example code that shows the
process of fitting the GDN+ model.

⬇

from gnnad.graphanomaly import GNNAD

from gnnad.generate import GenerateGaussian,

… GenerateAnomaly

# generate sample data

gengauss = GenerateGaussian()

X = gengauss.generate()

# split train test

X_train, X_test = X[:P,:], X[P:, :]

# generate anomalies on test set

anoms = GenerateAnomaly(X_test)

X_test = anoms.generate(anoms.variability,

… lam = 3, prop_anom = 0.07, seed=45)

X_test = anoms.generate(anoms.drift,

… lam = 11, prop_anom = 0.07, seed=234)

y_test = anoms.get_labels()

# instantiate and fit GDN model object

model = GNNAD(threshold_type="max_validation",

… topk=6, slide_win=200)

fitted_model = model.fit(X_train, X_test,

… y_test)

# sensor based threshold based on GDN+

pred_label = fitted_model.sensor_threshold_

… preds(tau = 99)

# print evaulation metrics

fitted_model.print_eval_metrics(pred_label)

## Results

This section presents a summary of the main findings for anomaly detection
using GDN/GDN+ on both simulated and real-world data. To ensure quality of
data, the aim for practitioners is to maximise the ability to identify
anomalies of different types, while minimising false detection rates. We
define the following metrics in terms of true positive (TP), true negative
(TN), false positive (FP) and false negative (FN) classifications. In other
words, the main priority is to minimise FN, while maintaining a reasonable
number of FP, such that it is not an operational burden to check the total
number of positive flags. Accordingly, we use recall, defined by $\frac{{\rm
TP}}{{\rm TP}+{\rm FN}}$, to evaluate the performance on the test set and to
select hyperparameters. That is, the proportion of actual positive cases that
were correctly identified by the model(s).We also report the performance using
precision $\big{(}\frac{{\rm TP}}{{\rm TP}+{\rm FP}}\big{)}$, accuracy
$\big{(}\frac{{\rm TP}+{\rm TN}}{{\rm TP}+{\rm TN}+{\rm FP}+{\rm FN}}\big{)}$
and specificity $\big{(}\frac{{\rm TN}}{{\rm TN}+{\rm FP}}\big{)}$.

To evaluate model performance, three existing anomaly detection models are
used as benchmarks: 1. The naive (random walk) Autoregressive Integrated
Moving Average Model (ARIMA) prediction model from [26, 7]; 2. HDoutliers
[17], an unsupervised algorithm designed to identify anomalies in high-
dimensional data, based on a distributional model that allows for probability
assignment to an anomaly; and 3. DeepAnT [14], an unsupervised, deep learning-
based approach to detecting anomalies in time series data.

Dataset #Train #Test %Anomalies SimEuc $3,000$ $1,000$ $13.2$ SimRiver $3,000$
$1,000$ $13.2$ Herbert $12,745$ $3,499$ 58.0

Table 1: Details of the three data sets used in the case studies.

### Simulation Study: Benchmark data

Data are generated using the linear-mixed model described in Equation (2),
with differing spatial dynamics: SimEuc where the random effect,
$\boldsymbol{Z}$, is characterised by Euclidean distance only, and SimRiver
where $\boldsymbol{Z}$ simulates complex river network dynamics, using the
same site locations and covariate values, $\boldsymbol{X}$ (see accompanying
code for details). Detecting anomalies that involve multiple consecutive
observations is a difficult task that often requires user intervention, and is
the focus of this study. We consider scenarios with drift and high-variability
anomaly types, which together contaminate $13.2\%$ of the test data, given
$n_{\rm drift}=5,{}_{\rm drift}=11,n_{\rm variability}=24,{}_{\rm
variability}=3$, see Table 1.

Figure 3:  A selection of time series from the SimEuc data set; sensor 8 and
sensor 37 separated by a short Euclidean distance share high correlation
(Pearson‘s coefficient of 0.65). A Savitzky-Golay filter smoothens the time
series (purple line). On the bottom, sensor 23 illustrates high-variability
and drift anomalies, consecutively (red dots). Figure 4:  A selection of
series from the SimRiver data set; sensor 8 and sensor 37 are not flow-
connected sites, and share low correlation (Pearson‘s coefficient of 0.07).
Drift anomalies (red dots) are shown in data from sensor 4. A Savitzky-Golay
filter is used to smoothen the time series for visual representation (purple
line).

Figure 3 visualises aspects of the SimEuc dataset, where sensor 37 and sensor
8 are in close (Euclidean) proximity, resulting in a high correlation between
the locations (0.65), as anticipated. Note that sensor 23 exhibits anomalous
behavior, high-variability and drift, consecutively, over time. Compared to
the SimRiver dataset, shown in Figure 4, we note how the time series from
sensor 37 and sensor 8 are no longer strongly correlated (0.07), despite their
close proximity, as they are not flow-connected in the simulated river
network.

Table 2 shows the performance of the different anomaly detection methods. The
best performance in terms of recall is highlighted in bold, while the second-
best performance is underlined. We observe that GDN/GDN+ outperformed most
other models in terms of recall, which is the fraction of true positives among
all actual positive instances. Specifically, GDN has a recall score of 83.3%
on SimEuc and 72.7% on SimRiv, while GDN+ has the second highest recall score
of 85.6% on SimEuc and 78.0% on SimRiv. Although ARIMA performed best in terms
of recall, the high percentage of detected anomalies, 81.6% and 85.6%, is
impractical to use (discussed below) and results in low accuracy scores of
28.6% and 25.3%, respectively. HDoutliers did not flag any time point as
anomalous in any test case. DeepAnT tends to classify most samples as
negative, resulting in a low recall score. These results suggest that GDN and
its variants are best able to detect a high proportion of actual anomalies in
the datasets.

Figure 5:  Anomaly detection performance of ARIMA, GDN, and GDN+ for the
simulation study, in terms of true positive (TP), true negative (TN), false
positive (FP) and false negative (FN).

Figure 5 shows the classifications of anomalies in the simulation study in
terms of true positive (TP), true negative (TN), false positive (FP) and false
negative (FN). For reasons mentioned, recall is the performance metric of
interest, but we also consider the trade-off between recall and precision.
Lower precision means that the model may also identify some normal instances
as anomalies, leading to false positives. In the context of river network
anomaly detection, FP may be manually filtered, but it is critical to minimise
FN. Note that GDN+ improves in comparison to GDN in minimising the FN count,
but at the cost of increasing FP, in both data sets. Such a trade-off is
acceptable and considered an improvement. However, while ARIMA has the highest
recall score, the number of FP is impractical for practitioners. We also note
that drift anomalies are harder to detect than high-variability anomalies,
with drift as the majority of FN counts, in all cases.

The authors of the GDN model demonstrated its efficacy in detecting anywhere-
within-system failures at time $t$ by applying a threshold to all sensors
within a system. However, the use of sensor-based thresholds in GDN+ has the
advantage of indicating anomalies at the individual sensor level. In the
context of monitoring river networks, it is crucial to identify the anomalous
sensor, $i$, at a given time $t$. The percentage of true positives detected at
the correct sensor, $i$, using the sensor-based anomaly threshold, $A_{i}(t)$,
in GDN+, was 92% and 89% for SimEuc and SimRiver, respectively. Similarly, the
rate of true positives detected in the neighbourhood of the correct sensor $i$
were 96% and 91%, respectively. This granularity of information is essential
for large networks consisting of independent sensors that are separated by
significant spatial distances, where the cost of time and travel for sensor
replacement or maintenance is substantial.

Data Model Rec Prec Acc Spec SimEuc HDoutliers 0.0 0.0 86.8 100.0 ARIMA 88.6
14.4 28.6 19.4 DeepAnT 3.8 11.9 83.5 95.7 GDN 83.3 55.3 88.9 89.7 GDN+ 85.6
48.1 85.9 85.9 SimRiv HDoutliers 0.0 0.0 86.8 100.0 ARIMA 91.7 14.2 25.3 15.1
DeepAnT 0.8 9.1 85.9 98.8 GDN 72.7 54.2 88.3 90.6 GDN+ 78.0 43.1 83.5 84.3

Table 2: Anomaly detection performance in terms of recall (%), precision (%),
accuracy (%), and specificity (%) of GDN and its variants and baseline methods
for the simulation study.

### Case Study: Herbert River

This case study examines water-level data collected from eight sites located
across the Herbert river, a major river system located in the tropical region
of Australia, as shown in Figure 6. The time series data is highly non-
stationary, characterised by river events caused by abnormal rainfall
patterns, with some coastal sites exhibiting shorter periodicity trends which
can be attributed to tidal patterns, see Figure 7. The spatial relationships
are complex, and depend on the surrounding water catchment areas, spatial
rainfall patterns, dams, and other impediments. In-situ sensors are prone to
various anomalies such as battery failure, biofouling (accumulation of
microorganisms, plants, algae, or small animals), and damage. In some cases,
anomalies can manifest as the absence of a water event (i.e., flatlining)
rather than the presence of abnormal time series patterns (i.e., spikes,
variability, drift). In real-world scenarios, anomalies can persist for
extended periods, and resolving them may require traveling to remote locations
to inspect and repair sensors. As seen in Figure 7, anomalies at time $t$ are
largely attributed to Sensor 4, which was out of water for long periods of
time.

Figure 6:  The Herbert river system and sensor locations.

The Herbert river is a challenging data set for all of the anomaly detection
models, due to the sparse placement of sensors across the network, i.e., fewer
sensors at greater distances apart resulting in weaker spatial relationships,
and the test set contains a high proportion of anomalies (58%; see Table 1).
GDN applied to the real-world dataset yields a recall of $29.2\%$, with GDN+
improving recall to $34.8\%$, see Table 3. Model performance suffers primarily
due to the failure to detect the large anomaly spanning across 22-01-10 in
Figure 7. This may be attributed to the learned graph relationships being
characterised by river events in the training data, and without such events,
it is difficult to identify when a sensor is flat-lining. However, the model
successfully identified anomalies spanning across 2021-12-27 and 22-01-03,
which coincided with river events.

Data Model Rec Prec Acc Spec Herbert HDoutliers 0.0 0.0 42.0 100.0 ARIMA 30.5
62.9 51.3 77.5 DeepAnT 1.6 39.7 44.0 97.0 GDN 29.2 60.6 50.1 76.2 GDN+ 34.8
59.2 50.4 70.0 GDN++ 100.0 80.7 86.7 70.0

Table 3: Anomaly detection performance in terms of recall (%), precision (%),
accuracy (%), and specificity (%) of GDN and its variants and baseline methods
for the Herbert river case study.

Figure 8 shows the learned graph adjacency matrix, $A$. Sensor 1, separated
geographically by a large distance from the other sensors, has weaker
relationships with the other nodes. Interestingly, the attention weights
indicate that sensor 3 is strongly influenced by sensor 6, with tidal patterns
from sensor 6 being evident in the predictions of sensor 3. Large error scores
indicating an anomaly spanning across 2021-12-27 are primarily observed in
sensor 2, impacted by anomalous sensors 3 and 4, where the predicted values
are low. However, due to the small network of sensors in this case study, it
is difficult to determine which sensor had the anomaly.

Adjusting the threshold in the GDN+ model can only enhance performance to a
limited extent, since some anomalies may have insignificant error scores due
to the underlying time series model. For instance, the anomaly spanning across
22-01-10 has negligible error scores because the prediction model was
performing well and no river events had occurred, making it challenging to
detect the anomaly in this particular scenario. Therefore, relying solely on
the model may not be sufficient, and we recommend implementing some basic
expert rules. We introduce another model variant GDN++, by applying a simple
filter to the time series ensuring that all values are positive, used in
conjunction with adjusting the threshold calculation (GDN+). That is, an
anomaly at time, $t$, on sensor, $i$, is flagged if,
$\max\\{\mdmathbb{I}\\{y_{i}^{(t)}<0\\},A_{i}(t)\\}=1.$ GDN++ successfully
detects all anomalies.

Figure 7:  Test data of water level collected from sensors across the Herbert
river (light blue), and the corresponding predictions from the GDN model (dark
blue). Actual anomalies (red dots) are shown, along with the predicted
anomalies (orange line) on the bottom. Figure 8:  Graph with learned adjacency
matrix, $A$. The edge weightings are determined by ij, which indicates
attention and depend on model input, ${\rm X}^{(t)}$.

## Discussion

Multivariate anomaly detection and prediction models for spatio-temporal
sensor data have the potential to transform water quality observation,
modelling, and management [8, 7]. The provision of trustworthy sensor data has
four major benefits: 1. It enables the production of finer-scale, reliable and
more accurate estimates of sediment and nutrient loads, 2. It provides real-
time feedback to landholders and managers, 3. It guides compliance with water
quality guidelines, and 4. It allows for the assessment of ecosystem health
and the prioritisation of management actions for sustaining aquatic
ecosystems. However, technical anomalies in the data provided by the sensors
can occur due to factors such as low battery power, biofouling of the probes,
and sensor miscalibration. As noted by [7], there are a wide variety of
anomaly types present within in-situ sensor data (e.g., high-variability,
drift, spikes, shifts). Most anomaly detection methods used for water quality
applications tend to target specific anomalies, such as sudden spikes or
shifts [27, 28, 29]. Detecting persistent anomalies, such as sensor drift and
periods of abnormally high variability, remains very challenging for
statistical and machine learning research. Such anomalies are often overlooked
by distance and kernel based methods, yet must be detected before the data can
be used, because they confound the assessment of status and trends in water
quality. Understanding the relationships among, and typical behaviours of,
water quality variables and how these differ among climate zones is thus an
essential step in distinguishing anomalies from real water quality events.

We investigated the graph-based neural network model, GDN [15], for its
ability to capture complex interdependencies between different variables in a
semi-supervised manner. As such, GDN offered the ability to capture deviations
from expected behaviour within a high-dimensional setting. We developed novel
bench-marking data sets for subsequence anomaly detection (of variable length
and type), with a range of spatio-temporal complexities, inspired by the
statistical models recently used for river network data [9]. Results showed
that GDN tended to outperform the benchmarks in anomaly detection. We
developed a model extension, GDN+, by adjusting the threshold calculation.
GDN+ was shown to further improve performance. Sensor-based thresholds also
proved useful in terms of identifying which neighbourhood the anomaly
originated from in the simulation study.

We used a real-world case study of water level in the Herbert river, with non-
stationary time series characterised by multiple river events caused by
abnormal rainfall patterns. In this case, most of the anomalies arose due to
the sensors being out-of-water. Considering an individual time series, such
anomalies may not appear obvious, as it is the failure to detect river events
(in a multivariate setting) that is indicative of an anomalous sensor. Despite
the challenges in the data, GDN+ was shown to successfully detect technical
anomalies when river events occurred, and combined with a simple expert-based
rule, GDN++, all anomalies were successfully detected.

Future applications of this work could consider the separation of river events
from technical anomalies. In terms of interpretability, since the prediction
model aggregates feature data from neighbouring sensors, anomalous data can
affect the prediction of any other sensor within the neighbourhood. Therefore,
large error scores originating from one sensor anomaly can infiltrate through
to the entire neighbourhood, and impairs the ability to attribute an anomaly
to a sensor.

In summary, this work extends and examines the practicality of the GDN
approach when applied to an environmental monitoring application of major
international concern, with complex spatial and temporal interdependencies.
Successfully addressing the challenge of anomaly detection in such settings
can facilitate the wider adoption of in-situ sensors and could revolutionise
the monitoring and management of air, soil, and water.

## Conclusions

In this work, we used the Graph Deviation Network (GDN) approach for anomaly
detection on river networks, based on sensors that generate high-dimensional
data with complex spatio-temporal relationships. We introduce an extension of
the model (GDN+/GDN++) and assess the practicality of this approach on various
data sets. Simulation studies showed that GDN was effective in accurately
identifying anomalies, further improved by GDN+/GDN++, with the improved
ability to locate anomalies across a network. Experiments on a real-world data
set indicated that GDN, and its variants, are most effective in capturing
river network dynamics when compared to other methods, but challenges still
remain for anomaly detection on river network sensor data.

### Author contributions

KB; Investigation, Formal Analysis, Methodology, Software, Validation,
Visualization, Writing – Original Draft Preparation. RS; Methodology,
Supervision, Writing – Review & Editing. KM; Funding Acquisition,
Conceptualization, Methodology, Supervision, Writing – Review & Editing. ES;
Data Curation, Visualization, Writing – Review & Editing.

### Grant information

This work was supported by the Australian Research Council (ARC) Linkage
Project (LP180101151) titled “Revolutionising water-quality monitoring in the
information age". Case study data were provided by the Department of
Environment and Science, Queensland, and is available as part of the Python
package gnnad.

### Acknowledgements

Łukasz Mentel; Software. James Mcgree; Writing – Review & Editing.

## References

  * ass [2015] UN General Assembly Transforming our World: The 2030 Agenda for Sustainable Development.  _United Nations_ , (A/RES/70/1), 2015.
  * Silva et al. [2022] Gabriel Marinho e Silva, Daiane Ferreira Campos, José Artur Teixeira Brasil, Marcel Tremblay, Eduardo Mario Mendiondo, and Filippo Ghiglieno.  Advances in technological research for online and in situ water quality monitoring—A review.  _Sustainability_ , 14(9):5059, 2022.
  * Ritchie et al. [2003] Jerry C Ritchie, Paul V Zimba, and James H Everitt.  Remote sensing techniques to assess water quality.  _Photogrammetric engineering & remote sensing_, 69(6):695–704, 2003.
  * Kang et al. [2009] James M Kang, Shashi Shekhar, Michael Henjum, Paige J Novak, and William A Arnold.  Discovering teleconnected flow anomalies: A relationship analysis of dynamic neighborhoods (RAD) approach.  In _International Symposium on Advances in Spatial and Temporal Databases_ , pages 44–61. Springer, 2009.
  * Blázquez-García et al. [2021] Ane Blázquez-García, Angel Conde, Usue Mori, and Jose A Lozano.  A review on outlier/anomaly detection in time series data.  _ACM Computing Surveys (CSUR)_ , 54(3):1–33, 2021\.
  * Bourgeois et al. [2003] Wilfrid Bourgeois, Anne-Claude Romain, Jacques Nicolas, and Richard M Stuetz.  The use of sensor arrays for environmental monitoring: interests and limitations.  _Journal of Environmental Monitoring_ , 5(6):852–860, 2003.
  * Leigh et al. [2019] Catherine Leigh, Omar Alsibai, Rob J Hyndman, Sevvandi Kandanaarachchi, Olivia C King, James M McGree, Catherine Neelamraju, Jennifer Strauss, Priyanga Dilini Talagala, Ryan DR Turner, et al.  A framework for automated anomaly detection in high frequency water-quality data from in situ sensors.  _Science of the Total Environment_ , 664:885–898, 2019\.
  * Ver Hoef and Peterson [2010] Jay M Ver Hoef and Erin E Peterson.  A moving average approach for spatial statistical models of stream networks.  _Journal of the American Statistical Association_ , 105(489):6–18, 2010.
  * Santos-Fernandez et al. [2022] Edgar Santos-Fernandez, Jay M Ver Hoef, Erin E Peterson, James McGree, Daniel J Isaak, and Kerrie Mengersen.  Bayesian spatio-temporal models for stream networks.  _Computational Statistics & Data Analysis_, 170:107446, 2022.
  * Porter et al. [2012] John H Porter, Paul C Hanson, and Chau-Chin Lin.  Staying afloat in the sensor data deluge.  _Trends in Ecology & Evolution_, 27(2):121–129, 2012.
  * Rodriguez-Perez et al. [2020] Javier Rodriguez-Perez, Catherine Leigh, Benoit Liquet, Claire Kermorvant, Erin Peterson, Damien Sous, and Kerrie Mengersen.  Detecting technical anomalies in high-frequency water-quality data using artificial neural networks.  _Environmental Science & Technology_, 54(21):13719–13730, 2020.
  * Wang et al. [2018] Xing Wang, Jessica Lin, Nital Patel, and Martin Braun.  Exact variable-length anomaly detection algorithm for univariate and multivariate time series.  _Data Mining and Knowledge Discovery_ , 32:1806–1844, 2018\.
  * Hu et al. [2019] Min Hu, Xiaowei Feng, Zhiwei Ji, Ke Yan, and Shengchen Zhou.  A novel computational approach for discord search with local recurrence rates in multivariate time series.  _Information Sciences_ , 477:220–233, 2019.
  * Munir et al. [2018] Mohsin Munir, Shoaib Ahmed Siddiqui, Andreas Dengel, and Sheraz Ahmed.  DeepAnT: A deep learning approach for unsupervised anomaly detection in time series.  _Ieee Access_ , 7:1991–2005, 2018.
  * Deng and Hooi [2021] Ailin Deng and Bryan Hooi.  Graph neural network-based anomaly detection in multivariate time series.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 4027–4035, 2021.
  * Schmidl et al. [2022] Sebastian Schmidl, Phillip Wenig, and Thorsten Papenbrock.  Anomaly detection in time series: a comprehensive evaluation.  _Proceedings of the VLDB Endowment_ , 15(9):1779–1797, 2022.
  * Wilkinson [2017] Leland Wilkinson.  Visualizing big data outliers through distributed aggregation.  _IEEE transactions on visualization and computer graphics_ , 24(1):256–266, 2017.
  * Goldstein and Uchida [2016] Markus Goldstein and Seiichi Uchida.  A comparative evaluation of unsupervised anomaly detection algorithms for multivariate data.  _PloS one_ , 11(4):e0152173, 2016.
  * Nassif et al. [2021] Ali Bou Nassif, Manar Abu Talib, Qassim Nasir, and Fatima Mohamad Dakalbab.  Machine learning for anomaly detection: A systematic review.  _Ieee Access_ , 9:78658–78700, 2021.
  * Velivcković et al. [2018] Petar Velivcković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  _International Conference on Learning Representations (ICLR)_ , 2018\.
  * Buchhorn et al. [2022] Katie Buchhorn, Kerrie Mengersen, Edgar Santos-Fernandez, Erin E Peterson, and James M McGree.  Bayesian design with sampling windows for complex spatial processes.  _arXiv preprint arXiv:2206.05369_ , 2022.
  * Cressie et al. [2006] Noel Cressie, Jesse Frey, Bronwyn Harch, and Mick Smith.  Spatial prediction on a river network.  _Journal of Agricultural, Biological, and Environmental Statistics_ , 11(2):127, 2006.
  * Ver Hoef et al. [2006] Jay M Ver Hoef, Erin Peterson, and David Theobald.  Spatial statistical models that use flow and stream distance.  _Environmental and Ecological Statistics_ , 13(4):449–464, 2006.
  * Shreve [1966] Ronald L Shreve.  Statistical law of stream numbers.  _The Journal of Geology_ , 74(1):17–37, 1966\.
  * Peterson and Fausch [2003] Douglas P Peterson and Kurt D Fausch.  Upstream movement by nonnative brook trout (salvelinus fontinalis) promotes invasion of native cutthroat trout (oncorhynchus clarki) habitat.  _Canadian Journal of Fisheries and Aquatic Sciences_ , 60(12):1502–1516, 2003.
  * Hyndman and Athanasopoulos [2021] Rob J Hyndman and George Athanasopoulos.  _Forecasting: principles and practice_.  OTexts, third edition, 2021.
  * Talagala et al. [2019] Priyanga Dilini Talagala, Rob J Hyndman, Catherine Leigh, Kerrie Mengersen, and Kate Smith-Miles.  A feature-based procedure for detecting technical outliers in water-quality data from in situ sensors.  _Water Resources Research_ , 55(11):8547–8568, 2019.
  * Hill and Minsker [2010] David J Hill and Barbara S Minsker.  Anomaly detection in streaming environmental sensor data: A data-driven modeling approach.  _Environmental Modelling & Software_, 25(9):1014–1022, 2010.
  * Ba and McKenna [2015] Amadou Ba and Sean A McKenna.  Water quality monitoring with online change-point detection methods.  _Journal of Hydroinformatics_ , 17(1):7–19, 2015.