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, ¹ water quality is a major environmental concern worldwide. The use of in-situ [ 1] 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, confounding 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, ⁴ 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 ⁵ ), including sensor drift and periods of unusually high or low variability, which may indicate the necessity for sensor maintenance or calibration. ⁶ Such anomalies are difficult to detect and often associated with high false negative rates. ⁷

Earlier statistical studies have focused on developing autocovariance models based on within-river relationships to capture the unique spatial characteristics of rivers. ⁸ Although these methods are adaptable to different climate zones, and have recently been extended to take temporal dependencies into account, ⁹ data sets generated by in-situ sensors still pose significant computational challenges with such prediction methods due to the sheer volume of data. ¹⁰ Autocovariance matrices must be inverted when fitting spatio-temporal models and making predictions, and the distances between sites must be known. Previous work, ¹¹ 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 optimisation 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. ⁵ One such method uses dynamic clustering on learned segmented windows to identify global and local anomalies. ¹² 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. ¹³ 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 Ref. 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 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), ¹⁵ and explore refinements with respect to anomaly scoring that address the needs of environmental monitoring. The GDN approach ¹⁵ 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.

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.

Our work makes four primary contributions to the field: 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.

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. 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 the considered models.

Consider multivariate time series data Y = y 1 … y T , obtained from n sensors over T time ticks. The (univariate) data collected from sensor i = 1 , … , n at time t = 1 , … , 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 y t is assigned a binary label, a t ∈ 0 1 , where a t = 1 indicates an anomaly at time t , anywhere across the full sensor network.

The GDN approach ¹⁵ for anomaly detection is composed of two aspects: 1.

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

The above components are described in more detail below.

Forecasting-based Time Series Model. To predict y t , the model takes as input w ∈ ℕ lags of the multivariate series, X t ≔ y t − 1 … y t − w .

The i -th row (containing sensor i ’s measurements for the previous w lags) of the above input matrix is represented by the column vector, x i t = y i t − 1 … y i t − w . Prior to training, the practitioner specifies acceptable candidate relationships via the sets C 1 , … , C n , where each C i ⊆ 1 2 … 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 v i ∈ ℝ d for i = 1 , … , 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 = v i Τ v j ‖ v i ‖ ‖ v j ‖ I j ∈ C i , i , j ∈ 1 … n , with ‖ ⋅ ‖ denoting the Euclidean norm, and indicator function I which equals 1 when node j belongs to set 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 j 1 … e jn . A graph-adjacency matrix (and in turn a graph itself) is then constructed from the sensor similarities via: A ji = I e ji ≥ e j , K ∪ i = j , for user-specified K ∈ 1 … 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 v k k = 1 n yield a graph. Next, the lagged series are fed individually through a shallow Graph Attention Network ²⁰ 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, x i t ∈ ℝ w . Allow a parameter weight matrix W ∈ ℝ d × w to apply a shared linear transform to each node. Then, the output of the network is given by z i t = max 0 ∑ j : A ji > 0 α ij W x i t , where 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 = LeakyReLU a Τ v i ⊕ W x i t + v j ⊕ W x j t , where α ij = exp π ij ∑ k : A ki > 0 exp π ik , (1) with learnable parameters a ∈ ℝ 2 d , where ⊕ denotes concatenation, and LeakyReLU x ≔ max δ x x for δ > 0 , with the maximum operation applied elementwise. Note the addition [ 2] in Equation 1 and that ∑ 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, y ̂ t = f η v 1 ⊙ z 1 t … v n ⊙ z n t , where f η : ℝ d × n → ℝ n is a feedforward neural network with parameters η , and ⊙ denotes element-wise multiplication. The model is trained by optimizing the parameters v i i = 1 n , W , a , and η to minimize the mean squared error loss function L = 1 T − w ∑ t = w + 1 T y ̂ t − y t 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 − y ̂ i t | , with | ⋅ | denoting the absolute value, and the vector of prediction error for each sensor is, ε i ∈ ℝ T − w . 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, ε ˜ i = ε i − Median ε i IQR ε i , where IQR denotes inter-quartile range. In the original work by Ref. 15, a time point t is flagged as anomalous if, A t = I max i ε ˜ i , t > κ , using the notation ε ˜ i , t for the error value at the t -th index. Alternatively, the authors recommend using a simple moving average of ε ˜ 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, ε ˜ . However, this is applied to all sensors as a fixed threshold.

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 TP TP + 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 TP TP + FP , accuracy TP + TN TP + TN + FP + FN and specificity TN TN + FP .

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 Refs. 7, 26 ^{; 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.}

Multivariate anomaly detection and prediction models for spatio-temporal sensor data have the potential to transform water quality observation, modelling, and management. ^{7 , 8} 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 Ref. 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. ^{28 – 30} 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, ¹⁵ 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. ⁹ 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. A replication study with multiple benchmarking data sets demonstrated consistency in these results. 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 appeared as flat-lining, due to the river drying out. 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.

There are two recent methodological extensions to the GDN approach. First, the Fused Sparse Autoencoder and Graph Net ³¹ which extends the GDN approach by augmenting the prediction-based loss with a reconstruction-based term arising from the output of a sparse autoencoder, and a further extension that allows the individual sensors to have multivariate data. Second, a probabilistic (latent variable) extension is trained using variational inference. ³² Since these were published contemporaneously to the present research, and only the latter provided accompanying research code, these approaches were not considered in the paper. Future extensions of this work could consider incorporating the above methodological extensions, as well as developing on the existing software package.

Other applications 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. The extensive body of literature on signal processing for source identification has the potential to inspire future solutions in addressing this issue. ^{33 , 34}

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.

This work studied the application of Graph Deviation Network (GDN) based approaches for anomaly detection on the challenging setting on river network data, which often feature sensors that generate high-dimensional data with complex spatio-temporal relationships. We introduced alternative defection criteria for the model (GDN+/GDN++), and their practicality was explored on both real and simulated benchmark data. The findings indicated that GDN and its variants were effective in correctly (and conservatively) identifying anomalies. Benchmark data were generated via an approach that was also introduced in this paper, along with open-source software, and may serve useful in the development and testing of other anomaly-detection methods. In short, we found that graph neural network based approaches to anomaly detection offer a flexible framework, able to capture and model non-standard, highly dynamic, complex relationships over space and time, with the ability to flag a variety of anomaly types. However, the task of anomaly detection on river network sensor data remains a considerable challenge.

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