Dynamic stochastic modeling for adaptive sampling of environmental variables using an AUV

Gunhild Elisabeth Berget; Jo Eidsvik; Morten Omholt Alver; Tor Arne Johansen

doi:10.1007/s10514-023-10095-8

Dynamic stochastic modeling for adaptive sampling of environmental variables using an AUV

Berget, Gunhild Elisabeth; Eidsvik, Jo; Alver, Morten Omholt; Johansen, Tor Arne 2023-04-01 00:00:00 Discharge of mine tailings signiﬁcantly impacts the ecological status of the sea. Methods to efﬁciently monitor the extent of dispersion is essential to protect sensitive areas. By combining underwater robotic sampling with ocean models, we can choose informative sampling sites and adaptively change the robot’s path based on in situ measurements to optimally map the tailings distribution near a seaﬁll. This paper creates a stochastic spatio-temporal proxy model of dispersal dynamics using training data from complex numerical models. The proxy model consists of a spatio-temporal Gaussian process model based on an advection–diffusion stochastic partial differential equation. Informative sampling sites are chosen based on predictions from the proxy model using an objective function favoring areas with high uncertainty and high expected tailings concentrations. A simulation study and data from real-life experiments are presented. Keywords Adaptive sampling · Gaussian process · Stochastic modeling · AUV · Mine tailings 1 Introduction The world’s oceans are heavily impacted by human activities, such as pollution and overﬁshing (Halpern et al., 2008). One example is the harmful effects of mine tailings, drill cuttings, We acknowledge support from the Norwegian Research Council and mud discharge. Along the Norwegian continental shelf, (RCN) through the Centers of Excellence funding scheme, Project Number 223254 - Centre for Autonomous Marine Operations and there are sensitive species like cold-water coral reefs and Systems (NTNU-AMOS), and the INDORSE Project 267793. The sponges which are negatively impacted by such deposition authors would like to thank Tore Mo-Bjørkelund from NTNU and (Trannum et al., 2010), but still, there exist several active Frédéric Py from SINTEF Digital for help related to the ﬁeld mining seaﬁlls and new seaﬁlls are in the planning phase. experiments, and Tor Nordam and Finn Are Michelsen from SINTEF Ocean for running the numerical models SINMOD and DREAM and Improved methods for monitoring the extent of the contam- supplying the data sets used in this paper. ination from seaﬁlls are important to ensure that sensitive areas are protected. In this paper, we use a seaﬁll located in B Gunhild Elisabeth Berget Frænfjorden (see Fig. 1) as a case study where we validate gunhild.berget@ntnu.no our proposed method for monitoring the tailings distribution Jo Eidsvik close to the outlet. jo.eidsvik@ntnu.no When mapping the tailings distribution, models of ocean Morten Omholt Alver dynamics can assist us with helpful information. How- morten.alver@ntnu.no ever, building such models is challenging because of the Tor Arne Johansen oceans’ large-scale nonlinear processes and high spatio- tor.arne.johansen@ntnu.no temporal variability. Existing models use numerical methods Department of Engineering Cybernetics and Centre for to improve their accuracy (Grifﬁes et al., 2000) continuously, Autonomous Marine Operations and Systems, NTNU, but they still have errors due to simpliﬁcations, limitations Trondheim, Norway in numerical implementation, uncertain parameters, and ini- Sintef Community, Trondheim, Norway tial and boundary conditions. Thus, data assimilation using Department of Mathematical Sciences, NTNU, Trondheim, in situ measurements can help improve the model accuracy. Norway 123 Autonomous Robots This function ensures that the area is explored by choosing locations assumed to be information-rich and also consid- ers the travel length limitations of the AUV, leading to an adaptive sampling strategy. In our case study, the method is tested in simulation using the Frænfjorden area (Norway), a fjord containing a seaﬁll for mine tailings. The fjord can be seen in Fig. 1, where the location is marked with a red circle. We compare the proposed method to four simpler strategies in the simulation study. As a proof of concept, real-life experiments are carried out in the same area. Fig. 1 Map of the Frænfjorden area. The pink contours in the ﬁg- The paper is organized as follows: Sect. 2 presents back- ure show the location of the factory Omya Hustadmarmor, and the ground information on the Frænfjorden area, the ocean black lines indicate the pipes of the seaﬁll. The numbers (2.1–2.4) and models, and adaptive sampling. In Sect. 3, the development (1.1–1.2) mark possible outlets for the tailings. At the time of the ﬁeld of the proxy model is presented. Section 4 presents a com- experiment, outlet 2.3 was the operational discharge. The blue rectan- gle represents the operational area for the simulation study and ﬁeld bined state space model explaining how the data assimilation experiments (Color ﬁgure online) is done before the sampling strategies are discussed in Sect. 5. Section 6 contains simulation results comparing our pro- posed method to four simpler strategies, and results from experimental ﬁeld trials are presented in Sect. 7. Section 8 Measurements are commonly obtained using remote sensing, ships, or buoys, but this data is expensive to collect, making discusses the results before Sect. 9 gives a conclusion and it hard to observe the environment in detail. Recently, robotic future work. vehicles like autonomous underwater vehicles (AUVs) have become more robust and less expensive, and by equipping 1.1 Related work and contribution them with sensors, one can collect measurements (sample) more efﬁciently. Still, the possible observations are sparsely With the improved technology and robustness of AUVs, distributed, leading to signiﬁcant gaps in our understanding adaptive sampling using these robotic platforms has had of the ocean (Stewart, 2008). This measurement limitation increased focus lately, and a review of the existing adap- prompts the need for intelligent and targeted observation tive sampling methods for AUV missions can be found in designs, and an important question is when and where to Hwang et al. (2019). Most adaptive sampling methods uti- sample. lize GPs (Cressie & Wikle, 2011), which are widely used to This paper addresses this question and suggests an adap- create non-parametric and computationally efﬁcient models. tive sampling method utilizing prior data from ocean models They are popular not only for adaptive sampling but used in and real-time in situ measurements from an AUV. Because various applications like neuroimaging (Hyun et al., 2016), the already existing ocean dynamics models have a high com- machine learning (Rasmussen & Williams, 2005) and sen- putational load, making them unﬁt for running on embedded sor placement (Krause et al., 2008). Combining GPs with robotic systems, an onboard low-complexity proxy model robotic vehicles is amongst others explored in Zhang et al. is created so that real-time adaption of the mission can be (2012) where an AUV is used to track an upwelling front, obtained. The proxy model represents the state of the ocean and in Luo and Sycara (2018), which explore a method of at the time and can be updated when new information is multiple vehicles using a mixture of GP models. Foss et al. available. We use a spatio-temporal Gaussian process (GP) (2022) use a spatio-temporal model for the same case we are with a kernel describing spatial dependencies and a temporal studying here, but with another objective function relying on model built using an advection–diffusion stochastic partial specifying a threshold for the particle concentration. differential equation (SPDE). The linear SPDE model and This work presents an implementation of an adaptive sam- the Gaussian measurement model are written as a state-space pling method for environmental sensing of mine tailings model, and the Kalman ﬁlter is used to update the state recur- close to a seaﬁll. We show how to utilize prior information sively over time. Two existing numerical models, SINMOD from ocean models to create GP models. Similar approaches, (ocean model) and DREAM (particle dynamics model) are also utilizing prior knowledge to learn GP models, can be used to train the compact GP model creating a proxy model found in Das et al. (2015) who use an AUV to collect samples of the tailing distribution. Sensor readings on the AUV can for ex-situ analysis, selecting the sampling locations based on then be used to update the proxy model onboard in real-time. previous mission data and maximizing a utility function, and An objective function for waypoint selection is presented to in Ma et al. (2018), who create a model of a spatio-temporal maximize the value of information from the measurements. environment by learning and reﬁning the hyperparameters of 123 Autonomous Robots a GP kernel using prior ocean model data. In our approach, discharge for the data in the simulation study and the ﬁeld the ocean model data are used to develop a non-stationary experiments is on the western pipe. This outlet is marked as correlation kernel as in Berget et al. (2018) and Fossum et al. 2.3 in the ﬁgure. (2018) and set the parameters in a temporal model based on the stochastic advection–diffusion partial differential equa- 2.2 Ocean models tion (Sigrist et al., 2015). This allows the model to capture critical dynamic phenomena such as tides. Estimated solu- Numerical oceanographic models are used in this paper to tions of the SPDE can be written in a linear form so that the build a prior belief of the state of the ocean and to predict Kalman ﬁlter can be used to update the model state. An objec- the movement of particles with the currents. Ocean mod- tive function based on variance reduction, hotspot detection, els provide information on salinity, currents, temperature, and AUV operation time constraints is used; see also Berget density, and pressure and describe the state of the ocean at et al. (2018). The approach sequentially selects waypoints a given time. The models use a set of hydrodynamic and based on the current model state, and the method is tested in thermodynamic equations called the primitive equations, and both simulations and ﬁeld experiments near a mining seaﬁll these are commonly solved numerically using ﬁnite differ- in a Norwegian fjord. ence methods. This makes the models computer intensive; The main contributions of this paper are: hence, high-resolution predictions can only be computed for limited areas, depending on the available computer resources. – Propose an adaptive sampling method for environmental Input to ocean models typically includes tides, sea-level sensing of mine tailings. pressure, wind, heat exchange, bathymetry, and freshwater – Use ocean model data to set the parameters in a lin- runoff. Ocean model outputs are uncertain due to limits in ear temporal model based on a stochastic advection– model resolution, limits to our knowledge of the processes diffusion partial differential equation. resolved by the model, errors in initial values, boundary con- – Present results from real-life experiments using the pro- ditions and parameter values, and inaccuracies in numerical posed adaptive sampling method. implementations (Lermusiaux et al., 2006). It is, therefore, necessary to validate the output against observations, and data assimilation can be used to correct the model state esti- 2 Background mates. For monitoring of releases of mine tailings or other particles or dissolved material, there is a need to develop enabling technology for efﬁcient and targeted ocean sam- This section presents deﬁnitions and background informa- tion on ocean models, the spatio-temporal model, and data pling. Observations from different platforms such as buoys, assimilation. ship-based sampling, or robotic sampling from, for example, AUVs can be used to evaluate or, through data assimilation, 2.1 The Frænfjorden area improve the model estimates. In this paper, we use data from two speciﬁc models to In Norway, several active and planned seaﬁlls have impli- build the stochastic onboard model, SINMOD (ocean dynam- cations for the marine environment (Ramirez-Llodra et al., ics model) (Slagstad & McClimans, 2005; Wassmann et al., 2015; Morello et al., 2016). The tailings typically consist of 2006) and DREAM (particle transportation model) (Rye et particulates and residual chemicals from processing, and this al., 1998, 2008). Nepstad et al. (2020) described the SIN- waste will smother the sea bed close to the sea ﬁll. In addi- MOD and DREAM setups for Frænfjorden and evaluated the models against ﬁeld measurements. The relation between tion, some of the waste may be transported by ocean currents away from the outlet, thus affecting a larger area of the sea SINMOD, DREAM, and the onboard model can be seen in bed and the water column. Fig. 2. More details on the two ocean models are given below. As the case study for this paper, we focus on a seaﬁll located in Frænfjorden, Norway (see Fig. 1). A calcium car- 2.2.1 SINMOD: ocean dynamics model bonate processing plant in the factory has been operating for several decades. This processing of ore from mining causes SINMOD is a numerical ocean model system that connects ﬁne-ground waste (tailings) (Kvassnes & Iversen, 2013), dis- and simulates physical and biological processes in the ocean posed of on the ocean ﬂoor, in a seaﬁll close to the factory. (Slagstad & McClimans, 2005). The model is based on The outlets in Frænfjorden are located at approximately the Navier–Stokes equations and uses a nesting technique 20 m depth, and tailings are continuously discharged with a where high resolution models obtain their boundary con- time-varying discharge rate. A typical discharge rate from the ditions from larger model domains with lower resolution. outlets is around 400 m /h. The discharge alternates between SINMOD is established in conﬁgurations with horizontal two pipes, illustrated by black lines in Fig. 1. The operational resolutions ranging from 20 km to 32 m. Input to the model 123 Autonomous Robots the onboard proxy model. This data is also used to estimate parameters in the SPDE model, including the GP kernel. 2.3 Adaptive sampling Observing the ocean by taking measurements or samples is crucial to increasing our understanding of ocean processes. Ocean sampling can be seen as the action of sensing or mea- suring one or more physical properties in the ocean at a speciﬁc point in space and time. Traditional methods include mounting sensors to buoys which limits the sampling to only one location, or immersing instruments in the water from boats, commonly using a predetermined plan of sampling locations. Exploiting robotic platforms and especially AUVs can intensify ocean sampling, giving more data and better insight. We need more sophisticated sampling strategies to exploit the beneﬁts that robotic platforms have regarding autonomy and increased mobility. To obtain the most scientiﬁcally rel- evant measurements effectively, adaptive approaches need Fig. 2 The relation between SINMOD, DREAM, and the adaptive sam- to be developed. In contrast to static/pre-scripted schemes, pling framework. SINMOD provides ocean currents data as input to the adaptive/data-driven strategies can react to measured condi- particle transportation model DREAM. Both SINMOD and DREAM tions, having access to both prior and in situ information; provides data used to build the prior belief of the onboard model. Ocean selecting sampling locations depending on past observations current predictions are used in the temporal part of the model and parti- cle concentration predictions on the boundary are used in the updating taken during exploration. of the model The beneﬁts of adaptive sampling will depend on the situ- ation. There are numerous studies on the value of additional (sequential) sampling efforts in different contexts, say clin- includes atmospheric forcing, freshwater runoff, and bound- ical trials in medical research (Jennison & Turnbull, 1999) ary conditions (temperature, salinity, tidal forcing, currents). or geoscience data (Eidsvik et al., 2015). In our case, there Ocean dynamics data from SINMOD is used in this case would be limited to gain by adaptation if a priori information study as input to DREAM and the onboard proxy model (see makes it possible to pre-script a high-value AUV plan before Fig. 2). Atmospheric forcing for daily forecast simulations deployment. However, if there is high prior variability, the was obtained from NOAA’s Global Forecast System with ability to adapt the initial plan one way or the other, based on 0.5 resolution (for large-scale model domains) and from onboard measurements, would bring added value by guiding the Norwegian Meteorological Institute’s MetCoOp Ensem- the AUV to more relevant places. ble Prediction System with 2.5 km resolution (for local model domains). 3 Modeling 2.2.2 DREAM: particle transportation model This section describes the advection–diffusion SPDE model DREAM is a Lagrangian particle transport model which and its associated spatio-temporal GP model representation. can be used to simulate the behavior and fate of marine We also outline parameters in the model. pollutants, including particulate discharges from drilling operations (Rye et al., 1998, 2008). It provides time series 3.1 Spatio-temporal onboard model of concentrations of released materials in the water column and deposition of these materials onto the sea ﬂoor. Input to To select optimal sampling locations in a dynamic domain, the DREAM model includes hydrodynamic data delivered by information about the spatial conditions at all times is essen- SINMOD and information about the release (amount, rate, tial. Since the computational demands using ocean models densities, grain size distribution). DREAM is often used as a such as SINMOD and DREAM are too high for real-time decision support tool for managing operational discharges to onboard modeling, we suggest a stochastic spatio-temporal the marine environment. This paper uses data from DREAM proxy model based on an advection–diffusion SPDE. When to set the initial belief for the mine tailings distribution on this SPDE is discretized, it represents a spatio-temporal GP 123 Autonomous Robots model. Apart from having a smaller computational footprint, B(t , s) describes external information, e.g. boundary infor- GPs are conventional tools for statistical modeling of spatial mation. data and have been widely adopted in oceanographic appli- This paper assumes horizontal, isotropic diffusion lead- cations (Binney et al., 2013). In this paper, the proxy model ing to a constant diffusion parameter D. We denote the t t T is used to approximate the distribution of the particle con- 2-dimensional drift vector as c(t , s) = (c , c ) at time step e n centration near the seaﬁll in Frænfjorden. t. Then we get There exist various types of spatio-temporal GP mod- ∂ X (t , s) ∂ X (t , s) ∂ X (t , s) els, see e.g. Chapter 6 of Cressie and Wikle (2011)for a t t = λX (t , s) − c + c e n ∂t ∂s ∂s review. For statisticians, it has been important to understand e n 2 2 the spatio-temporal covariance function of such processes, ∂ X (t , s) ∂ X (t , s) + D + + B(t , s) + (t , s). which is crucial in spatio-temporal Kriging. In the simplest 2 2 ∂s ∂s e n case, the covariance function is separable in space and time, (2) but there has also been much work on non-separable covari- ance functions (Gneiting, 2002). Stochastic differential and The SPDE is discretized using the upwind differencing difference equations in the spatial domain ﬁt naturally into scheme to ﬁnd an approximate solution. We deﬁne the length state-space models, which enable spatio-temporal Kalman of a time step as dt. Forward differences are used in time, ﬁltering using in-situ observations. The spatial version of the giving autoregressive process is one such example of a linear dif- ∂ X (t + dt , s) − X (t , s) ference equation. So are integro-difference equations which X (t , s) ≈ . (3) weight spatial variables from the previous time step to pre- ∂t dt dict the spatial variables at the current time step (Storvik et We deﬁne a regular grid D with a total of N grid points al., 2002). Recently, extensions of the advection–diffusion with resolution ds in the east direction and ds in the north e n equation have been developed (Sigrist et al., 2015; Richard- direction. The set of spatial locations (east, north) is then son, 2017). These models provide more interpretability in s = (s , s ) ∈ D ⊂ R , i = 1,..., N.For theﬁrst i e,i n,i s the spatial weighting as well as opportunities for embedding derivative in space, we alternate between forward and back- time-dependent drift. We use one of these models in the cur- ward differences based on the direction of the drift vector c . rent paper. The advection–diffusion equation is also endorsed In the eastward direction, we have in the physical modeling community, and in our setting the advection can be speciﬁed from the tidal currents and the ∂ X (t , s ± (ds , 0) ) − X (t , s) X (t , s) ≈ , (4) diffusion coefﬁcients can be set from the ocean models SIN- ∂s ds e e MOD and DREAM. We consider a real-valued stochastic process {X (t , s), t ≥ t where the upper sign is used for c < 0 and the lower is used 0, s ∈ D}, representing the uncertain particle concentration t when c ≥ 0. Similarly, in the northward direction, we have on log scale in (east,north) location s = (s , s ) ∈ D ⊂ R e n at time t. T ∂ X (t , s ± (0, ds ) ) − X (t , s) X (t , s) ≈ , (5) The spatio-temporal process of transport of a medium in ∂s ds n n a ﬂuid can be expressed by the advection–diffusion SPDE (Sigrist et al., 2015) where the upper version is used for c < 0 and the lower is used when c ≥ 0. The second derivatives are discretized as ∂ X (t , s) follows, = λX (t , s) − c(t , s) ∇ X (t , s) ∂t +∇ D∇ X (t , s) + B(t , s) + (t , s), (1) X (t , s) ≈ ∂s T T X (t , s + (ds , 0) ) − 2X (t , s) + X (t , s − (ds , 0) ) T e e with c(t , s) = (c (t , s), c (t , s)) is the drift vector ﬁeld for e n (6) ds advection, D the diffusion matrix, (t , s) a noise term, ∇ e ∂ ∂ the gradient operator , and λ ∈[−1, 0] a damping ∂s ∂s e n X (t , s) ≈ ∂s constant controlling the autoregressive relationship between T T state vectors (Richardson, 2017). The noise term (t , s) is X (t , s + (0, ds ) ) − 2X (t , s) + X (t , s − (0, ds ) ) n n . (7) assumed to be a zero-mean Gaussian noise term that is 2 ds uncorrelated in the temporal domain but spatially correlated, Corr((t , s), (t , s )) = R(|s − s |), where R() represents By scaling the time discretization, a change in time can be a correlation function for distance input. The source term denoted t +1 instead of t +dt.Wehave {t = t + j ×dt ; j = 123 Autonomous Robots 0, 1, 2 ..., T }. When inserting all the discretized terms into The prior covariance matrix Σ is then given by eq. (1), the scheme can be written as a linear process on the ⎡ ⎤ form Σ Σ ... Σ 11 12 1N ⎢ ⎥ Σ Σ ... Σ 21 22 2N ⎢ ⎥ Σ = , (10) ⎢ ⎥ 0 . . . . . . . ⎣ ⎦ . . . X = A X + B + , ∼ N (0, W ), t = 1,..., T , t t t −1 t t t t Σ Σ ... Σ N 1 N 2 NN (8) where Σ = σ σ R and R represents the spatial cor- ij i j ij ij relation. Hence, the diagonal of the covariance matrix in where X is the size N ×1 vector of particle concentrations at equation (10) contains the location-wise variances σ , and time step t for all spatial locations, A is the N × N propaga- the off-diagonal elements describe the covariance between tor matrix including information about advection, diffusion the locations. The covariance function is chosen by com- and damping, B is a matrix adding external information, in paring the empirical covariance of the ocean model training our case boundary conditions, and W is a positive deﬁnite data with known valid parametric covariance functions. The variance-covariance matrix for the noise term . Matern (3/2) kernel (Matérn, 2013) is ﬁtted to such covari- The particle concentrations at time t = 0 areassumedtobe ance function here; Gaussian distributed with mean vector μ =[μ ,...,μ ] 1 N and N × N variance-covariance matrix Σ . For the boundary conditions in the SPDE, we use Dirichlet R = (1 + φh ) exp(−φh ), (11) ij ij ij BCs on the east side of our survey area. This can be consid- ered as a constant forcing along the east side and is chosen where h =|s − s | is the Euclidean distance between two ij i j because of the constant ﬂow of mine tailings that are let locations s and s and φ> 0 is a constant meta-parameter i j out on this side of the rectangle. We assume no drift across regulating the correlation decay with the distance. The best the boundary on the remaining sides of the survey area and value for φ is speciﬁed using training data. therefore use Neumann BCs. The boundary conditions are set To specify covariance model parameters and visualize using predicted particle concentration data from DREAM. how the spatial correlation decays as a function of the distance between points in the domain, it is common in geostatistics to use the variogram of the data, see e.g. 3.2 Parameter specification in the SPDE model Cressie and Wikle (2011). The variogram works on differ- ence X (s ) − X (s ) between spatial data. Doing so, it tends t j t i The advection and diffusion parameters in the SPDE (1)must to be more robust to background trends than the covariance be speciﬁed. The velocity ﬁeld is set from physical assump- function. The theoretical relation between the variogram γ() t t T tions, where c (s) = (c (s), c (s)) as the forcing of the e n and the covariance is convection in our case is treated as the currents in the ocean. The values are set using predicted currents from the ocean γ(s − s ) = Var(X (s ) − X (s )) = i j t i t j model SINMOD. The isotropic diffusion D is settobethe same as for the computations in the DREAM model D = 0.1 (Var(X (s )) + Var(X (s )) − 2Cov(X (s ), X (s ))). t i t j t i t j m /s. The initial distribution is mainly estimated from the (12) DREAM data. Since ocean processes are affected by cur- rents, wind patterns, bathymetry, and freshwater runoff, we Figure 3 shows results from the variogram analysis. The will have elevated variability in some locations. Thus, we empirical variogram (blue) is computed by averaging the choose a non-stationary prior model (Jun & Stein, 2008) square differences of variable pairs having spatial distance where the mean and variance in each location are speciﬁed within bins. We use regular bins of distances (ﬁrst axis). empirically from the training data available by DREAM: Three different theoretical variograms were ﬁtted to the train- ing data. All of them ﬁt the empirical variogram adequately, but the Matern variogram (green) was chosen here because 2 2 it is closer to the empirical variogram for very small spa- σ = (y − μ ) . (9) i ,m i M − 1 tial distances. This indicates that a Matern spatial ﬁeld has a m=1 smoothness that ﬁts the DREAM data better. We see that the variogram increases with distance and then ﬂattens out. Here, Here, μ = y is the empirical mean concentra- the effective spatial correlation length is about 500 m. For the i i ,m m=1 tion in location s over all M realizations of training data. Matern correlation function in equation (11) this means that 123 Autonomous Robots represented by the initial onboard model, and observations are taken from sensors on the AUV. We deﬁne our observation model as Y = G X + ω , ω ∼ N (0,τ I ), (13) t t t t t where Y is a length m vector of observations at time t, G is t t a time-varying observation matrix of size m × N that maps the state variables to the observations, and ω represent the observational error which is modeled as a zero-mean Gaus- sian with standard deviation τ , and I representing the identity matrix. Each row of G has a one-entry on the vector position associated with an observed grid cell, and the other entries are zero. This method constraints the sampling locations to Fig. 3 The empirical variogram is shown in blue, and a Matern(3/2), an exponential, and a Gaussian variogram is ﬁtted to the data. We ﬁnd be at the predeﬁned grid cells. Continuous path sampling, as an average variance of 0.6 and an effective spatial correlation length of in Choi and How (2010), could improve the sampling method about 500 m (Color ﬁgure online) and should be considered in future work. Since the SPDE model deﬁnes a GP with a linear process representation, the spatial Kalman ﬁlter can be used to update φ ≈ 0.01. The variogram converges to the average variance our state at every time step. Letting Z = (Y ,..., Y ) rep- t 1 t which is about 0.6. resent all observations and prior information up to time t, The advection–diffusion model parameters are further set the ﬁltering distribution is X |Z ∼ N (μ , Σ ) with mean t t t |t t |t from the physical assumptions in the ocean model. From μ and variance-covariance Σ deﬁned by t |t t |t visual inspection, this GP model appears to give reasonable simulations of the ﬁeld but with more patchy randomness Σ = A Σ A + W than what is seen in the ocean model. As commonly seen t |t −1 t t −1|t −1 t T 2 in physical models data, the variogram computed from the Q = G Σ G + τ I t t |t −1 t t ocean models data starts to ﬂuctuate after a while. In many T −1 K = Σ G Q t t |t −1 t t cases this can be related to physical principles. Notably, the μ = A μ + B t t GP model is a statistical proxy model that can be run onboard t |t −1 t −1|t −1 but is potentially missing some large-scale oceanographic μ = μ + K (Y − G μ ) t t t t |t t |t −1 t |t −1 features. The ocean model is a numerical description with Σ = Σ − K Q Σ . (14) t |t t |t −1 t t |t −1 much skill but potentially getting systematic bias due to mis- speciﬁed inputs growing with time. They both have their This recursion starts by the initial distribution with mean μ merits. and Σ , where X ∼ N (μ , Σ ). The boundary conditions 0 0 0 The covariance matrix W in the noise term of the equation and the model parameters going into the matrices and vectors is set to have the same correlation structure as the covariance in equation (14) are all speciﬁed as described in Sect. 3.2. matrix for the initial state Σ ,giving W = V Σ , where V 0 t 0 is a constant scaling parameter. Considering this in relation to the original continuous-time SPDE, the discretization in 5 Sampling strategy time would also impact this variance, but in our work, we use a ﬁxed time discretization identical to the one used onboard. An objective function is suggested to obtain an informative There does not seem to be any damping over time in the ocean path for the AUV. The function is created based on three model data, but for stability reasons (Richardson, 2017)we criteria, as in Berget et al. (2019): specify a minimal level λ =−0.0010. 1. Locations with high variance are preferred 2. Locations with high predicted concentration are pre- 4 State space model and data assimilation ferred 3. Locations leading to a suitable travel length within the The state of concentrations is assumed to be modeled by next time interval for the AUV are preferred. equation (8). The aim of data assimilation is to combine observations with this prior process knowledge to obtain an optimal sequential posterior or ﬁltering estimate of the cur- The ﬁrst criterion is chosen because sampling in areas rent state (Sar, 2013). In this study, the prior knowledge is with a high model variance will reduce the total variance, 123 Autonomous Robots hence creating a more accurate model. This criterion also egy. Making the AUV ﬂy with the currents may increase ensures that the AUV travels to unexplored areas. The sec- endurance and reduce travel time between waypoints. This ond criterion makes the method adaptive. When studying the topic is amongst others explored in Pereira et al. (2013) and simulation results of the particle transport from the complex should be considered in future work. model, it is clear that the variability is highest where there is a high concentration of particles. Hence, the criterion is inspired by this observation and assumes that locations with 6 Summary of the suggested approach high predicted concentration will be rich with information. The last criterion comes from the travel length limitations of We now provide an overview of the suggested approach. At the AUV. When choosing the next sample location, it is also the beginning of the deployment, the GP model is initialized. essential that the AUV does not travel too far, assuring that The following steps are executed: the risk of the vehicle drifting outside the survey area is kept low. – Kernel hyperparameters are estimated using prediction The suggested objective function is then created by having and hindcast data from the particle transportation model a term for the ﬁrst two criteria. At time step t for location s , DREAM. (Described in Sect. 3.2.) the objective function is given by – The initial GP state deﬁned by μ and Σ are found. (Described in Sect. 3.2.) f (s ) = θ σ + θ μ , (15) t i 1 2 i ,t |t – Ocean current prediction data from SINMOD is used to i ,t |t deﬁne the advection ﬁelds c , which is used together with where the constant parameters θ =[θ ,θ ] deﬁne the weight- 1 2 the diffusion constant D to deﬁne transition matrices A . ing for each criterion. Deciding the weighting of criteria 1 (Described in Sect. 3.1.) and 2 depends on design choices. When the goal is to obtain – Particle transportation data from DREAM deﬁnes the a good map of the overall region, then criterion 1 should be boundaries that are incorporated in the vector B . dominating. However, if hotspots are the main focus, crite- (Described in Sect. 3.1.) rion 2 should dominate. In our case, the parameters θ are tuned in simulation to obtain a good balance between the The sampling method starts out either with the initial two criteria. GP model or with the latest updated GP model. The AUV The chosen sampling location s at time step t is deter- t performs adaptive sampling choosing the next waypoint as mined as the location that maximizes the objective function described in Sect. 5. The onboard model of the AUV is then f (s ) for s ∈[s ,..., s ]. Given the previous sampling t i i 1 N updated; the ﬁeld is moved with the currents, and new infor- location s , we choose the next sampling location as t −1 mation from observations is added. Updating is executed whenever a waypoint is reached, using data from observa- s = arg max f (s ) t i tions taken along the AUV transect in the assimilation. An overview of the adaptive sampling method is given in Algo- s.t. |s − s |≥ d i min t −1 rithm 1. |s − s |≤ d . (16) i max t −1 Algorithm 1 Adaptive Sampling method The constraints on the objective function ensure that the third Initialize GP model criterion in the objective function is adhered to by choosing while running do d and d such that a suitable travel length is obtained. min max Select new waypoint using (16) |s − s | is the Euclidean distance between the previous i Travel to new waypoint and collect observations t −1 sampling location and s . Equation (16) is solved using brute- if waypoint is reached then Update GP model with observations using (14) force enumeration; the value of the objective function (15) is found for all the grid cells within the region deﬁned by the constraints in (16). The next sampling waypoint s is chosen as the grid cell with the highest objective value. In this paper, ocean currents are included in the spatio-temporal proxy model, but they are not considered in the AUV nav- 7 Simulation study igation plans for selecting sampling locations. Rather, the AUV is set to surface often to avoid risk and reduce navi- A simulation study compares the proposed model with four gation uncertainty, as described in the experimental section. simpler strategies. We vary both the process model and In a larger-scale mission, it would certainly be beneﬁcial to the sampling strategy. The temporal model based on the include more of the ocean currents in the sampling strat- advection–diffusion SPDE is compared to a stationary pro- 123 Autonomous Robots Table 1 Overview of the ﬁve methods used for comparison in the sim- 7.1 Results and discussion ulation study In this section, we present results from the simulation study Method Process model Sampling strategy and discuss the performance of the methods. SPDEObj SPDE Objective SPDELawn SPDE Lawnmower 7.1.1 Performance measures SPDE SPDE No sampling Obj Stationary Objective To evaluate the different methods, three performance mea- Lawn Stationary Lawnmower sures are studied. The mean absolute error (MAE) is com- puted as the mean of the absolute error between the predicted concentration μ and the true ﬁeld from DREAM (April t |t 2nd) for all the grid cells and time steps. The mean standard Table 2 Performance measures for the ﬁve methods deviation (MSD) is computed similarly but then focuses on Method MAE MSD MO the standard deviation in the proxy model in all grid cells at SPDEObj 0.76 2.38 209.3 all time steps. Lastly, we compute the mean objective value SPDELawn 0.87 2.67 272.8 (see Eq. (15)) (MO) in all grid cells at every time step. Table SPDE 1.11 4.70 848.1 2 shows these performance measures computed as the mean Obj 0.85 2.52 230.2 over the total time steps. Figure 4 shows how the mean objec- tive value (MO) (see Eq. (15)) evolves through time, plotting Lawn 0.99 2.67 280.5 the mean value for every time step. The table shows the mean absolute error (MAE), the mean standard Considering the MAE values in Table 2, we see that deviation (MSD), and the mean objective function value (MO) SPDEObj performs best in this case. Including a non- stationary process model increases the method’s perfor- mance, as this is the case using both the Obj sampling strategy cess model where the next time step is equal to the previous and the Lawn sampling strategy. The objective function also time step, as in Berget et al. (2018). Different sampling increases the method’s performance compared to the pre- strategies are tested with sampling strategies based on the planned lawnmower strategy. Considering the MSD and MO suggested objective function, a predeﬁned sampling strat- measures, the SPDE method stands out with signiﬁcantly egy following a lawnmower pattern, and a method that uses higher values than the other methods. This is due to the only the SPDE model with no sampling (no observations are taken). The SPDE and Stationary process model are the same for all sampling strategies using the same parameter values given in the text. For the strategy using the objective function we choose the travel distance restrictions d = 180 meters min and d = 220 meters and set the time between each sample max to be 5 min. This gives a suitable travel speed for the AUV, and a constant time step simpliﬁes the updating of the tempo- ral model. When using the lawnmower pattern, the distance between each sampling location is set to 200 ms, and the same time step of 5 min is used. The parameters in the objec- tive function are set to θ = 15 and θ = 23. The different 1 2 strategies are shown in Table 1. Ocean model data from two consecutive days, 1st and 2nd of April 2013, is used in the simulation study. The data from April 1st is used to train the SPDE parameters as described in sections 3 and 4. Data from April 2nd is used as the test data, meaning that the predicted ﬁeld on April 2nd is our ground truth and the ﬁeld from which observations (samples) are drawn. We use realizations of the GP prior as the initial belief and run 100 replicate runs with the presented sampling methods in Table 1. We set the measurement error standard Fig. 4 The mean objective value (MO) (see Eq. (15)) computed over all deviation τ = 0.5 log(ppb). We simulate for a total of 6 h, the grid cells for each of the methods. Results are shown for the whole simulation time period spanning from 15:00 to 21:00 aiming to map the ﬁeld on April 2nd from 15:00 to 21:00. 123 Autonomous Robots Fig. 5 Hourly plot of the model standard deviation at the given time sampling location at the current time. The arrow indicates the direction together with the path for the different sampling strategies up until that of the AUV (Color ﬁgure online) time. The blue dot shows the path’s starting point, and the yellow is the method doing no sampling at all, and hence there is no way for bilize at a mean objective value between 200 and 300. For this method to decrease the proxy models’ variance, affecting both lawnmower methods the standard deviation and the objective values. (SPDELawn and Lawn), periodical ﬂuctuations in the objec- The effect of increasing variance with no sampling can tive function can be observed. This occurs when the sampling also be seen in Fig. 4. Again, the SPDE method stands out is done close to the east end of the survey area, where the with steadily increasing values while the other methods sta- assimilation with new data affects a minor part of the modeled ﬁeld compared to adding information in the center. Gener- 123 Autonomous Robots ally, the Obj and SPDEObj methods have a lower overall mean objective function compared with the other methods, indicating an intelligent sampling design. In Fig. 4, we see that the Lawn method has a lower MO value at some time steps compared to the SPDELawn method. This occurs at times when both Lawn approaches get high MO. Using the lawnmower pattern for sampling (see Fig. 5a) is the reason for this unexpected behavior. Sam- pling is performed in boxes, and some locations in the middle of the boxes get little inﬂuence from measurements. The SPDELawn method solves this by having the temporal model move neighboring sample values into this box. The Lawn method does not have a solution for this, and the values inside this box will be similar to the initial values. Since the initial mine tailings concentration in our example (see Fig. 7a) has low values all over the area, the tailing concentration inside this box will stay low throughout the whole simulation for the Lawn method. This will affect the MO value, leading to Fig. 6 The distribution in one single grid cell using the SPDEObj artiﬁcially low values, even though the performance (here, method. Showing mean (blue line) and 95% conﬁdence interval (red better illustrated by the MAE) can be poor. dotted lines) together with the true value (green line). The selected grid cell is located 340 ms east and 200 ms north in the survey area (Color ﬁgure online) 7.1.2 Sampling design and ﬁeld predictions Figure 5 shows the AUV path and waypoints for the different sampling strategies plotted on top of the standard deviation at every hour of the simulation study. The same lawnmower pattern is used for SPDELawn and Lawn and can be seen in out the whole mission, relying on observations to change its the ﬁrst column of the ﬁgure. The path chosen by the Obj predicted ﬁeld. This is not the case with SPDEObj, which has and the SPDEObj method can be seen in respectively the a non-stationary process model that moves the ﬁeld with the second and third columns. All the methods cover the area currents. Although data from two different days are used as quite well, but both Obj and SPDEObj focus on the eastern training and test data, both datasets are from the same ocean area. This region is closer to the outlet, and the true mine models. This may cause the comparison with ground-truth tailings distribution (seen in the ﬁrst column of Fig. 7)shows to be optimistic. Still, based on the simulation results, it is elevated concentrations in this area. Hence, more frequent reasonable to assume that SPDEObj can map the ﬁeld quite sampling here is expected because of the proposed objective well. The MAE scores of the two methods are indicated for function (see Eq. (15)). different times (Fig. 7). With the relatively high dynamical Figure 6 gives insight into how the GP methods behave, uncertainty, the MAE increases with time. The MAE is gen- showing the distribution in one single grid cell over time erally lower for the SPDEObj approach. using SPDEObj. The blue line shows the method’s mean Simulation studies with different AUV velocities have value, the green line is the true value, and the two red dotted shown that the SPDEObj method’s advantages are higher lines are the 95% conﬁdence interval. The point we are con- when the AUV takes fewer samples in the survey area over sidering is 340 ms east and 200 ms north in the survey area. the same time period. Naturally, lower AUV velocities led When sampling and assimilation are done close to this point, to fewer samples (comparable to having a bigger survey we see a reduction in the conﬁdence interval, and the mean area), while higher velocities led to an increase in the num- value tends to be corrected towards the actual value. ber of samples. When doubling the number of samples, the In Fig. 7 we focus on the two methods Obj and SPDEObj. difference in the performance between the approaches was The updated model mean for every hour is plotted in columns minimal, making it hard to decide which method was best. 2 and 3, and the true ﬁeld (the ﬁeld where the observations However, when operating with half the number of samples are drawn from) can be seen in the ﬁrst column. Comparing (compared to the simulation study above), the differences the methods with the true ﬁeld, they capture the trend quite became more evident, and the SPDEObj method showed well. The display illustrates an important difference between improved results. Hence, it is reasonable to think that the the two methods; because of its stationary process model, difference between the methods will increase with a larger the Obj method inherits the structure from the GP through- survey area. 123 Autonomous Robots Fig. 7 Hourly plot of the “true” value and the model mean value for the Obj and SPDEObj method. The MAE scores are included in the captions 7.1.3 Adaptive versus pre-planned designs pling design for one hour before the end time is shown in the variance plots. We now give an example comparing an the adaptive approach In the example, the adaptive approach gives a better esti- (SPDEObj) with a pre-planned method (SPDELawn) (see mate of the true ﬁeld than the pre-planned method. From the Fig. 8). The example shows the predicted results at the end true ﬁeld, we observe increased concentrations on the eastern of a mission run of total 5 h. The true ﬁeld is seen in Fig. 8a. side. The adaptive approach captures this information, and Figure 8b, c show the model mean and variance for the adap- because of elevated values in the area, continues to investi- tive approach, and Fig. 8d, e present the model state for the gate the eastern side. The pre-planned lawnmower approach, pre-planned approach. The AUV path describing the sam- on the other hand, is set to sample in the western area at the 123 Autonomous Robots Fig. 8 Example run comparing a pre-planned approach with an adap- approach, and Fig. 8d, e present the model state for the pre-planned tive approach. The true ﬁeld after ﬁve hours of sampling can be seen in approach. The AUV path showing the sampling design for one hour Fig. 8a. Figure 8b, c show the model mean and variance for the adaptive before the current time is shown in the variance plots time, and the samples taken are not as informative. Adap- MOD models, will ﬁnd high concentrations of mine tailings. tive approaches excel when conditions are not as expected. This is veriﬁed by our sensor data as well. Figure 11 shows In such cases, the ﬂexibility of an adaptive approach is ben- one transect (waypoint-to-waypoint) for the AUV including eﬁcial since the sampling design is automatically adjusted. sensor data values. It is clear from the ﬁgure that a layer with A pre-planned approach will always stick to the pre-scripted increased turbidity is reached. sampling design. 8.1 Robotic platform and implementation 8 Experimental results The robotic platform in our experiments consisted of a Light AUV from OceanScan (see Fig. 9) equipped with a Wetlabs As a proof of concept, ﬁeld experiments were conducted in EcoPuck sensor measuring total suspended matter (TSM), Frænfjorden. Two experiments were carried out on Decem- which is a measure of optical back-scattering in the form ber 17th and 18th, 2020. The operational discharge outlet at of an attenuation coefﬁcient. The measured sensor values the time of the experiments is marked as 2.1 in Fig. 1, and of TSM indicate turbidity, and this observation value was a rectangular operational area of size 1150 m × 250 m (also used to assimilate our model. It was converted to a mea- marked in the ﬁgure) was used. We focus on the depth layer sure of (log) concentration to adjust the value to our onboard at around 22–25 ms since this is a layer with a lot of variation model. This conversion is not straightforward, and we need and a layer where we, according to the DREAM and SIN- to determine the uncertainty τ . The Wetlabs EcoPuck sen- 123 Autonomous Robots a series of actions (e.g., Goto and Arrive_at). The set of all those actions forms a plan built while ensuring operational constraints of the mission (e.g., the vehicle should never dive deeper than a certain depth or leave a deﬁned area.) The plan is then sent to DUNE, which handles low-level control and execution. Before deploying the AUV, our method was tested in a simulated environment similar to our embedded system. In simulations, sensor readings from the AUV were replaced with data from SINMOD and DREAM, and an AUV simula- tor in DUNE was used to simulate the behavior of the AUV. Fig. 9 The light AUV (Harald) produced by OceanScan was used in The layout of the simulator is shown in Fig. 10. our experiments. The length of the AUV is 240 cm, and it is equipped with multiple sensors like CTD, Dissolved Oxygen, and Fluorometer 8.2 AUV path Inside the operational area, a 2-dimensional grid of size 58 × 12 was deﬁned, giving N = 684. The distance between each node was approximately 20 m in both grid directions, and each grid point was a possible waypoint for the AUV. The AUV was set to surface whenever reaching a new way- point to reduce the uncertainty caused by navigation errors and vehicle drift and for safety reasons. The AUV did a yoyo maneuver between 20 and 26 ms between each way- Fig. 10 The architecture of the embedded system (on top) and the point. As in the simulation study, we set d = 180 meters min simulator (bottom). The Uniﬁed Navigation Environment (DUNE) runs and d = 220 meters. The AUV speed is set to 1.2m/s, max onboard the AUV, and the agent architecture T-REX sits on top of this, and as opposed to the simulation study, where a constant allowing the embedding of multiple decision processes such as planning time of 5 min between two waypoints was used, we now let the time vary. The time used between two waypoints is found using the start and end times of the AUV transect. The temporal model and the updating are done when a way- sor user manual suggests an attenuation error of about 4 % −1 at attenuation coefﬁcient 1 m . With concentrations varying point is reached, updating the model with the amount of time between approximately 3 and 8 log(ppb), assuming measure- used for that particular transect. Figure 11 shows one tran- ment standard deviation of 4% gives τ between 0.12 and 0.32 sect (waypoint-to-waypoint) from our survey, starting with a log(ppb). However, both navigation errors and vehicle drift dive until reaching the selected depth layer, doing the yoyo can lead to uncertainty in the sampling location. Hence, we maneuver, and then surfacing in a corkscrew maneuver. The set τ = 0.5 log(ppb) for the experiments. corkscrew maneuver was chosen to obtain a sampling value Figure 10 shows the layout of the agent architecture used as close as possible to the waypoint. The yoyo maneuver was in the experimental survey. The Uniﬁed Navigation Envi- chosen to increase the possibility of ﬁnding the layer of ele- ronment (DUNE) (Pinto et al., 2012) runs onboard the AUV vated mine tailings concentration. The sensor reading used and is used for control, navigation, vehicle supervision, com- for updating the model was chosen as the maximum value in the layer 22–25 ms (the layer marked with red surfaces) for munication, and interaction with actuators. On top of this sits the agent architecture T-REX (Teleo-Reactive EXecu- each of the separate yoyos. For the transect in Fig. 11,this tive) (Rajan & Py, 2012), which enables an adaptive mission. would result in 6 updating values. Because the transect length T-REX allows the embedding of multiple complex decision could vary for different waypoints, the number of yoyos was processes (including planning). The communication between not ﬁxed but mostly alternated between 5 or 6 yoyos per DUNE and T-REX was handled by the LSTS toolchain (Pinto transect, giving 5 or 6 updating values per transect. et al., 2013), which provides the back-seat driver API to DUNE. This allows external controllers, such as T-REX, 8.3 Results to provide desired tasks for our platform while receiving progress updates on the current state. Our sampling algo- The survey area is located near the coast in a fjord arm, so rithm was written as a python script and was implemented the currents are highly inﬂuenced by the tidal cycle. To get as a reactor in the T-REX framework. A reactor is a com- a better insight into how the ocean currents behave at the ponent of T-REX acting as an internal control loop and is time of the missions, we consider the tidal cycle at the time capable of producing goals that the planner integrates into in Fig. 12. Our temporal model considers the tides’ varia- 123 Autonomous Robots before a tidal high, and the second is when the tides turn at the top. The period of the two missions is marked in red in the ﬁgure. Hourly predicted ocean currents between 09:00 and 13:00 from the SINMOD model can be seen in Fig. 13. As expected, the currents are highly inﬂuenced by the tides showing an inward ﬂow (eastward) before the tidal high (09:00–11:00) and weaker ocean currents during the tidal high (11:00–13:00). These predictions are used as input to our onboard model as the forcing in the SPDE model. Fig. 11 Showing the AUV transect from waypoint to waypoint, where the distance between the waypoints was 218 ms. The color of the transect 8.3.1 Mission 1: 08:53–10:44 indicates normalized mine tailings concentration values, and the two red surfaces show the depth layer (22–25 ms) where samples are collected (Color ﬁgure online) The ﬁrst mission was conducted between 08:53 and 10:44 on December 17th, before a tidal high. The estimated cur- rent from the SINMOD model can be seen in Figs. 13a, 13b, 13c. Figure 14 presents the onboard model results. The left column shows the updated mean of the onboard model, rep- resenting the mine tailings log(concentration) for selected time steps. The right column shows the onboard model stan- dard deviation at selected time steps and the AUV path up Fig. 12 The blue line shows the tidal cycle on December 17th, showing until the given time. The red line represents the AUV path, low tides at around 07:00 and high tides around 13:00. The red boxes and the red dots are the waypoints selected by our sampling indicate the time intervals for the two missions strategy. The results show increased tailing concentration in the eastern area, especially at the beginning of the mis- tions, and we choose to focus on two missions occurring sion (see Fig. 14a,14c and 14e). We observe a decay in the at different intervals in the tidal cycle. The ﬁrst mission is concentration in this area after 1.5 h (see Fig. 14g). This is Fig. 13 Hourly ocean current predictions from SINMOD between 09:00 and 13.00. The currents are highly inﬂuenced by the tides showing an inward ﬂow (eastward) before the tidal high (09:00–11:00) and weaker ocean currents during the tidal high (11:00–13:00) 123 Autonomous Robots Fig. 14 Results from mission 1. The left column shows the path up until the given time. The red line represents the AUV path, and updated mean of the onboard model, representing the mine tailings the red dots are the waypoints selected by our sampling strategy (Color log(concentration) for selected time steps. The right column shows the ﬁgure online) onboard model standard deviation at selected time steps and the AUV when the ocean currents get weaker due to the tidal cycle. In (see Fig. 13d and 13e). Figure 15 presents the onboard model the beginning, the chosen sampling strategy seems to focus results. The left column shows the updated mean of the on exploring the area, covering large parts of the whole onboard model, and the right column shows the standard rectangle. After the exploration, the focus switches to the deviation together with the selected AUV path. This mission eastern area, where high mine tailing concentrations have does not contain areas of clearly increased tailings concen- been observed. tration, and we see that the AUV chooses a different sampling strategy. Since no area shows increased concentration, there is no extra focus on any special area. Instead, the waypoints 8.3.2 Mission 2: 11:52–13:24 are chosen to cover as much as possible of the whole survey area during the whole mission. The second mission was conducted during a tidal high. Figure 15 shows results from this survey. The predicted SIN- MOD ocean currents are much weaker than for mission 1 123 Autonomous Robots Fig. 15 Results from mission 2. The left column shows the path up until the given time. The red line represents the AUV path, and updated mean of the onboard model, representing the mine tailings the red dots are the waypoints selected by our sampling strategy (Color log(concentration) for selected time steps. The right column shows the ﬁgure online) onboard model standard deviation at selected time steps and the AUV 9 Discussion with our assumptions on how the tailings concentration behave. Since the outlet is located on the east side of our The chosen sampling paths (see second column of Fig. 14 and area, we assume elevated concentration in the east. It is also 15) follow two different strategies. In the beginning, they are reasonable to assume that the concentration of mine tailings quite similar, focusing on exploring the area, but after the in the survey area will be higher when the currents ﬂow west- exploration, we see differences in waypoint selection. The ward, as is the case for mission 1. We have also compared sampling path in mission 1 (Fig. 14) focuses on the eastern the results with the ocean model DREAM data and observed area where increased concentration is observed, while the the same tendencies. However, this is not an accurate com- sampling path in mission 2, where there are no areas with parison since the proxy model is trained with data from this increased tailings concentration, continues to cover the whole model. region also after the initial exploration. For both surveys, the standard deviation is reduced at the end of the mission, indicating that the method manages to explore the area. The 10 Closing remarks selection of different strategies shows that the method is adaptive and will change depending on prior data and obser- In this work, we propose a method for adaptive sampling vations and that the sampling is focused on areas with high of ocean processes with an AUV. Prior information from tailings concentration. The AUV data are undoubtedly useful ocean models is combined with information from in situ for reﬁning the DREAM results. measurements to obtain the best possible sampling strategy. It is hard to rate the ﬁeld experiments’ performance since A spatio-temporal GP proxy model with a non-stationary there is no ground truth for this case. Still, the model mean covariance function and a process model based on ﬁnite dif- concentration (ﬁrst column of Fig. 14 and 15 coincide well ference solutions of an advection–diffusion SPDE are built. 123 Autonomous Robots This is a low-complexity model which can run onboard the References AUV. We can choose informative sampling locations by using Berget, G. E., Eidsvik, J., Alver, M., Py, F., Grøtli, E. I., & Johansen, an objective function favoring locations with high uncertainty T. A. (2019). Adaptive underwater robotic sampling of dispersal and high predicted particle concentration. The main contribu- dynamics in the coastal ocean. In Proceedings of the international tion of the paper is to have a dynamic (spatio-temporal) model symposium on robotics research Accepted. onboard a dynamic agent. This enables reliable uncertainty Berget, G. E., Fossum, T., Johansen, T. A., Eidsvik, J., & Rajan, K. (2018). Adaptive sampling of ocean processes using an AUV with quantiﬁcation in spatial predictions while actively search- a gaussian proxy model. IFAC-PapersOnLine, 51, 238–243. ing for ’hot-spots’ in space and time. The method is tested Binney, J., Krause, A., & Sukhatme, G. S. (2013). Optimizing waypoints in simulations and compared with simpler models and sam- for monitoring spatiotemporal phenomena. The International pling strategies. A simulation study shows that the SPDE Journal of Robotics Research, 32(8), 873–888. Choi, H.-L., & How, J. P. (2010). Continuous trajectory planning of process model gives improved prediction results in a dynamic mobile sensors for informative forecasting. Automatica, 46(8), environment, and the sampling strategy using the suggested 1266–1275. objective function performs better than the sampling strategy Cressie, N., & Wikle, C. (2011). Statistics for Spatio-Temporal Data. using a simpler lawnmower design. Since there is no easy CourseSmart Series. Hoboken: Wiley. Das, J., Py, F., Harvey, J. B., Ryan, J. P., Gellene, A., Graham, R., way to obtain a ground truth for ﬁeld experiments, it is difﬁ- Caron, D. A., Rajan, K., & Sukhatme, G. S. (2015). Data-driven cult to determine how well the method performs in the ﬁeld. robotic sampling for marine ecosystem monitoring. The Interna- Still, our results agree well with our assumptions of the mine tional Journal of Robotics Research, 34(12), 1435–1452. tailings distribution at the time, and the selection of differ- Eidsvik, J., Mukerji, T., & Bhattacharjya, D. (2015). Value of infor- mation in the earth sciences: Integrating spatial modeling and ent sampling strategies for two separate missions illustrates decision analysis. Cambridge: Cambridge University Press. the method’s adaptivity. Future work includes expanding the Foss, K. H., Berget, G. E., & Eidsvik, J. (2022). Using an autonomous model to 3D (with depth coordinates), adding current estima- underwater vehicle with onboard stochastic advection-diffusion tion from measurements onboard the AUV, including ocean models to map excursion sets of environmental variables. Envi- ronmetrics, e2702. currents in the waypoint selection method, and exploration Fossum, T., Eidsvik, J., Ellingsen, I., Alver, M., Fragoso, G., Johnsen, using multiple AUVs. It is also interesting to test such meth- G., Mendes, R., Ludvigsen, M., & Rajan, K. (2018). Information- ods for extended time periods, potentially combining the driven robotic sampling in the coastal ocean. Journal of Field AUV measurements with satellite data or data from drones. Robotics, 35, 1101–1121. Gneiting, T. (2002). Nonseparable, stationary covariance functions for Funding Open access funding provided by SINTEF. This work is sup- space-time data. Journal of the American Statistical Association, ported by the norges forskningsråd (Grand Nos. 223254, 267793) 97(458), 590–600. Grifﬁes, S. M., Böning, C., Bryan, F. O., Chassignet, E. P., Gerdes, Data availability The datasets generated during and/or analysed dur- R., Hasumi, H., Hirst, A., Treguier, A.-M., & Webb, D. (2000). ing the current study are available from the corresponding author on Developments in ocean climate modelling. Ocean Modelling, 2, reasonable request. 123–192. Halpern, B. S., Walbridge, S., Selkoe, K. A., Kappel, C. V., Micheli, F., D’Agrosa, C., Bruno, J. F., Casey, K. S., Ebert, C., Fox, H. Declarations E., Fujita, R., Heinemann, D., Lenihan, H. S., Madin, E. M. P., Perry, M. T., Selig, E. R., Spalding, M., Steneck, R., & Watson, R. (2008). A global map of human impact on marine ecosystems. Conﬂict of interest The authors declare that they have no conﬂict of Science, 319(5865), 948–952. interest. Hwang, J., Bose, N., & Fan, S. (2019). AUV adaptive sampling methods: Areview. Applied Sciences, 9(15). Open Access This article is licensed under a Creative Commons Hyun, J. W., Li, Y., Huang, C., Styner, M., Lin, W., & Zhu, H. (2016). Attribution 4.0 International License, which permits use, sharing, adap- STGP: Spatio-temporal gaussian process models for longitudinal tation, distribution and reproduction in any medium or format, as neuroimaging data. NeuroImage, 134, 550–562. long as you give appropriate credit to the original author(s) and the Jennison, C., & Turnbull, B. W. (1999). Group sequential methods with source, provide a link to the Creative Commons licence, and indi- applications to clinical trials. Boca Raton: CRC Press. cate if changes were made. The images or other third party material Jun, M., & Stein, M. L. (2008). Nonstationary covariance models for in this article are included in the article’s Creative Commons licence, global data. The Annals of Applied Statistics, 2(4), 1271–1289. unless indicated otherwise in a credit line to the material. If material Krause, A., Singh, A., & Guestrin, C. (2008). Near-optimal sensor is not included in the article’s Creative Commons licence and your placements in gaussian processes: Theory, efﬁcient algorithms and intended use is not permitted by statutory regulation or exceeds the empirical studies. Journal of Machine Learning and Research, 9, permitted use, you will need to obtain permission directly from the copy- 235–284. right holder. To view a copy of this licence, visit http://creativecomm Kvassnes, A., & Iversen, E. (2013). Waste sites from mines in Norwe- ons.org/licenses/by/4.0/. gian fjords. Mineralproduksjon, 3, A27–A38. Lermusiaux, P., Chiu, C.-S., Gawarkiewicz, G., Abbot, P., Robinson, A., Miller, R., Haley, P., Leslie, W., Majumdar, S., Pang, A., & Lekien, F. (2006). Quantifying uncertainties in ocean predictions. Oceanography, 19, 92–105. 123 Autonomous Robots Luo, W. & Sycara, K. (2018). Adaptive sampling and online learning in Storvik, G., Frigessi, A., & Hirst, D. (2002). Stationary space-time multi-robot sensor coverage with mixture of gaussian processes. In Gaussian ﬁelds and their time autoregressive representation. Sta- 2018 IEEE international conference on robotics and automation tistical Modelling, 2(2), 139–161. (ICRA) (pp. 6359–6364). Trannum, H. C., Nilsson, H. C., Schaanning, M. T., & Øxnevad, S. Ma, K.-C., Liu, L., Heidarsson, H., & Sukhatme, G. (2018). Data-driven (2010). Effects of sedimentation from water-based drill cuttings learning and planning for environmental sampling. Journal of Field and natural sediment on benthic macrofaunal community structure Robotics, 35(5), 643–661. and ecosystem processes. Journal of Experimental Marine Biology Matérn, B. (2013). Spatial variation. Meddelanden från Statens Skogs- and Ecology, 383(2), 111–121. forskningsinstitut, 36(5), 1–144. Wassmann, P., Slagstad, D., Riser, C. W., & Reigstad, M. (2006). Mod- Morello, E., Haywood, M., Brewer, D., Apte, S., Asmund, G., Kwong, elling the ecosystem dynamics of the barents sea including the Y., & Dennis, D. (2016). The ecological impacts of submarine marginal ice zone: II. Carbon ﬂux and interannual variability. Jour- tailings placement. In Hughes, R., Hughes, D., Smith, I. P., & nal of Marine Systems, 59(1): 1–24. Dale, A. (eds.), Oceanography and Marine Biology: An Annual Zhang, Y., Godin, M. A., Bellingham, J. G., & Ryan, J. P. (2012). Using Review, volume 54 (pp. 315–366). CRC Press, 1 edition. an autonomous underwater vehicle to track a coastal upwelling Nepstad, R., Liste, M., Alver, M. O., Nordam, T., Davies, E., & Glette, front. IEEE Journal of Oceanic Engineering, 37(3), 338–347. T. (2020). High-resolution numerical modelling of a marine mine tailings discharge in Western Norway. Regional Studies in Marine Publisher’s Note Springer Nature remains neutral with regard to juris- Science, 39, 101404. dictional claims in published maps and institutional afﬁliations. Pereira, A., Binney, J., Hollinger, G., & Sukhatme, G. (2013). Risk- aware path planning for autonomous underwater vehicles using predictive ocean models. Journal of Field Robotics, 30(5), 741– Gunhild Elisabeth Berget is Pinto, J., Calado, P., Braga, J., Dias, P., Martins, R., Marques, E., & a Ph.D. Student in Engineering Sousa, J. (2012). Implementation of a control architecture for Cybernetics at The Norwegian networked vehicle systems. IFAC Proceedings Volumes, 45(5), University of Science and Tech- 100–105. nology (NTNU). She received her Pinto, J., Dias, P., Martins, R., Fortuna, J., Marques, E., & Sousa, J. M.Sc. degree in Applied Mathe- (2013). The LSTS toolchain for networked vehicle systems. In matics in 2017 from NTNU. Her OCEANS 2013 MTS/IEEE Bergen: The Challenges of the Northern Ph.D. project is related to adap- Dimension. tive sampling using AUVs focus- Rajan, K. & Py, F. (2012). T-REX: Partitioned inference for AUV mis- ing on modeling, data assimila- sion control. Further advances in unmanned marine vehicles (pp. tion and path planning. 171–199). Ramirez-Llodra, E., Trannum, H., Evenset, A., Levin, L., Andersson, M., Finne, T., Hilário, A., Flem, B., Christensen, G., Schaanning, M., & Vanreusel, A. (2015). Submarine and deep-sea mine tailing placements: A review of current practices, environmental issues, natural analogs and knowledge gaps in norway and internationally. Jo Eidsvik is Professor of Statistics Marine Pollution Bulletin, 97(1), 13–35. in the Department of Mathemat- Rasmussen, C. E., & Williams, C. K. I. (2005). Gaussian processes for ical Sciences at the Norwegian machine learning (adaptive computation and machine learning). University of Science and Tech- Cambridge: The MIT Press. nology (NTNU). He got his M.Sc. Richardson, R. A. (2017). Sparsity in nonlinear dynamic spatiotemporal (1997) in Applied Mathematics models using implied advection. Environmetrics, 28(6), e2456. from the University of Oslo and Rye, H., Reed, M., & Ekrol, N. (1998). The PARTRACK model for his Ph.D. (2003) in Statistics from calculation of the spreading and deposition of drilling mud, chem- NTNU. His research interest are icals and drill cuttings. Environmental Modelling Software, 13(5), spatiotemporal statistics and com- 431–441. putational statistics, often applied Rye, H., Reed, M., Frost, T., Smit, M., Durgut, I., Johansen, O., & to the Earth sciences. Along with Ditlevsen, M. (2008). Development of a numerical model for cal- methods for parameter estimation, culating exposure to toxic and nontoxic stressors in the water prediction and data assimilation in column and sediment from drilling discharges. Integrated Envi- these settings, he also works with ronmental Assessment and Management, 4, 194–203. questions related to the design of experiments, sequential uncertainty Särkkä, S. (2013). Bayesian ﬁltering and smoothing. Cambridge Uni- reduction and value of information analysis in this context of spatio- versity Press: Cambridge. temporal applications. He co-authored a book on ‘Value of informa- Sigrist, F., Künsch, H., & Stahel, W. (2015). Stochastic partial differen- tion in the Earth Sciences’, Cambridge University press, 2015. He is tial equation based modeling of large space-time data sets. Journal Associate Editor of Statistics and Computing and Mathematical Geo- of the Royal Statistical Society: Series B (Statistical Methodology), sciences. 77, 3–33. Slagstad, D., & McClimans, T. (2005). Modeling the ecosystem dynam- ics of the barents sea including the marginal ice zone: I. Physical and chemical oceanography. Journal of Marine Systems, 58, 1–18. Stewart, R. H. (2008). Introduction To physical oceanography.Texas A&MUniversity. 123 Autonomous Robots Morten Omholt Alver received the M.Sc. and Ph.D .in Engineer- ing Cybernetics from The Nor- wegian University of Science and Technology (NTNU) in 2002 and 2007, respectively. He is currently an Associate Professor at the Department of Engineering Cyber- netics at NTNU. His research inter- ests include ocean modelling and data assimilation, modelling of bio- logical systems, and instrumen- tation and automation in marine aquaculture. Tor Arne Johansen received the M.Sc. degree in 1989 and the Ph.D. degree in 1994, both in electrical and computer engineer- ing, from the Norwegian Univer- sity of Science and Technology, Trondheim, Norway. From 1995 to 1997, he worked at SINTEF as a researcher before he was appointed Associated Professor at the Norwegian University of Sci- ence and Technology in Trond- heim in 1997 and Professor in 2001. He has published several hundred articles in the areas of control, estimation and optimization with applications in the marine, aerospace, automotive, biomedical and process industries. In 2002 Johansen co-founded the company Marine Cybernetics AS where he was Vice President until 2008. Prof. Johansen received the 2006 Arch T. Colwell Merit Award of the SAE, and is currently a prin- cipal researcher within the Center of Excellence on Autonomous Marine Operations and Systems (NTNU-AMOS) and director of the Unmanned Aerial Vehicle Laboratory at NTNU and the SmallSat Lab- oratory at NTNU. He recently co-founded the spin-off companies Scout Drone Inspection, UBIQ Aerospace, Zeabuz and SentiSystems. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Autonomous Robots Springer Journals http://www.deepdyve.com/lp/springer-journals/dynamic-stochastic-modeling-for-adaptive-sampling-of-environmental-kDw1OPYnAC

Loading next page...

References (44)

J. Pinto, P. Calado, José Braga, P. Dias, R. Martins, Eduardo Marques, Joel Sousa (2012)
Implementation of a Control Architecture for Networked Vehicle Systems
IFAC Proceedings Volumes, 45
J. Hyun, Yimei Li, Chao Huang, M. Styner, Weili Lin, Hongtu Zhu (2016)
STGP: Spatio-temporal Gaussian process models for longitudinal neuroimaging data
NeuroImage, 134
G. Berget, T. Fossum, T. Johansen, J. Eidsvik, K. Rajan (2018)
Adaptive Sampling of Ocean Processes Using an AUV with a Gaussian Proxy Model
IFAC-PapersOnLine, 51
E. Morello, M. Haywood, David Brewer, S. Apte, G. Asmund, Y. Kwong, D. Dennis (2016)
The Ecological Impacts of Submarine Tailings Placement
D. Slagstad, T. McClimans (2005)
Modeling the ecosystem dynamics of the Barents sea including the marginal ice zone: I. Physical and chemical oceanography
Journal of Marine Systems, 58
W. Arx (1962)
Introduction to Physical Oceanography
R. Nepstad, M. Liste, M. Alver, T. Nordam, E. Davies, Tormod Glette (2020)
High-resolution numerical modelling of a marine mine tailings discharge in Western Norway
Regional Studies in Marine Science, 39
C. Robert (2014)
Statistics for Spatio-Temporal Data
CHANCE, 27
Yanwu Zhang, M. Godin, J. Bellingham, J. Ryan (2012)
Using an Autonomous Underwater Vehicle to Track a Coastal Upwelling Front
IEEE Journal of Oceanic Engineering, 37
M. Jun, M. Stein (2008)
Nonstationary covariance models for global data
The Annals of Applied Statistics, 2
S. Ounpraseuth (2008)
Gaussian Processes for Machine Learning
Journal of the American Statistical Association, 103
G. Storvik, A. Frigessi, D. Hirst (2002)
Stationary space-time Gaussian fields and their time autoregressive representation
Statistical Modeling, 2
Kai-Chieh Ma, Lantao Liu, H. Heidarsson, G. Sukhatme (2017)
Data‐driven learning and planning for environmental sampling
Journal of Field Robotics, 35
Andreas Krause, A. Singh, Carlos Guestrin (2008)
Near-Optimal Sensor Placements in Gaussian Processes: Theory, Efficient Algorithms and Empirical Studies
J. Mach. Learn. Res., 9
K. Rajan, F. Py (2012)
T-REX: partitioned inference for AUV mission control
E. Ramirez-Llodra, H. Trannum, A. Evenset, L. Levin, M. Andersson, T. Finne, A. Hilário, B. Flem, G. Christensen, M. Schaanning, A. Vanreusel (2015)
Submarine and deep-sea mine tailing placements: A review of current practices, environmental issues, natural analogs and knowledge gaps in Norway and internationally.
Marine pollution bulletin, 97 1-2
J. Eidsvik, T. Mukerji, D. Bhattacharjya (2016)
Value of Information in the Earth Sciences: Integrating Spatial Modeling and Decision Analysis
Fabio Sigrist, H. Künsch, W. Stahel (2012)
Stochastic partial differential equation based modelling of large space–time data sets
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77
H. Trannum, H. Nilsson, M. Schaanning, S. Øxnevad (2010)
Effects of sedimentation from water-based drill cuttings and natural sediment on benthic macrofaunal community structure and ecosystem processes
Journal of Experimental Marine Biology and Ecology, 383
(1994)
Johansen received the M.Sc. degree in 1989 and the Ph.D. degree in 1994, both in electrical and computer
S. Särkkä (2013)
Bayesian Filtering and Smoothing
, 3
H. Rye, M. Reed, N. Ekrol (1998)
The PARTRACK model for calculation of the spreading and deposition of drilling mud, chemicals and drill cuttings
Environmental Modelling and Software, 13
P. Wassmann, D. Slagstad, C. Riser, M. Reigstad (2006)
Modelling the ecosystem dynamics of the Barents Sea including the marginal ice zone: II. Carbon flux and interannual variability
Journal of Marine Systems, 59
B. Halpern, Shaun Walbridge, K. Selkoe, Carrie Kappel, F. Micheli, C. D'Agrosa, J. Bruno, K. Casey, C. Ebert, H. Fox, R. Fujita, D. Heinemann, H. Lenihan, E. Madin, Matthew Perry, E. Selig, M. Spalding, R. Steneck, R. Watson (2008)
A Global Map of Human Impact on Marine Ecosystems
Science, 319
A. Pereira, J. Binney, Geoffrey Hollinger, G. Sukhatme (2013)
Risk‐aware Path Planning for Autonomous Underwater Vehicles using Predictive Ocean Models
Journal of Field Robotics, 30
J. Das, F. Py, J. Harvey, J. Ryan, A. Gellene, Rishi Graham, D. Caron, K. Rajan, G. Sukhatme (2015)
Data-driven robotic sampling for marine ecosystem monitoring
The International Journal of Robotics Research, 34
S. Griffies, C. Böning, F. Bryan, E. Chassignet, R. Gerdes, H. Hasumi, A. Hirst, A. Treguier, D. Webb (2000)
Developments in ocean climate modelling
Ocean Modelling, 2
C. Jennison, B. Turnbull (1999)
Group Sequential Methods with Applications to Clinical Trials
Jimin Hwang, N. Bose, Shuangshuang Fan (2019)
AUV Adaptive Sampling Methods: A Review
Applied Sciences
J. Binney, Andreas Krause, G. Sukhatme (2013)
Optimizing waypoints for monitoring spatiotemporal phenomena
The International Journal of Robotics Research, 32
T. Fossum, J. Eidsvik, I. Ellingsen, M. Alver, G. Fragoso, G. Johnsen, R. Mendes, M. Ludvigsen, K. Rajan (2018)
Information‐driven robotic sampling in the coastal ocean
Journal of Field Robotics, 35
(2008)
Development of a Numerical Model for Calculating Exposure to Toxic and Nontoxic Stressors in the Water Column and Sediment from Drilling Discharges
Wenhao Luo, K. Sycara (2018)
Adaptive Sampling and Online Learning in Multi-Robot Sensor Coverage with Mixture of Gaussian Processes
2018 IEEE International Conference on Robotics and Automation (ICRA)
J. Pinto, P. Dias, R. Martins, João Fortuna, Eduardo Marques, J. Sousa (2013)
The LSTS toolchain for networked vehicle systems
2013 MTS/IEEE OCEANS - Bergen
T. Gneiting (2002)
Nonseparable, Stationary Covariance Functions for Space–Time Data
Journal of the American Statistical Association, 97
F. Lutscher (2019)
Spatial Variation
Interdisciplinary Applied Mathematics
G. Berget, J. Eidsvik, M. Alver, F. Py, E. Grøtli, T. Johansen (2019)
Adaptive Underwater Robotic Sampling of Dispersal Dynamics in the Coastal Ocean
Robert Richardson (2017)
Sparsity in nonlinear dynamic spatiotemporal models using implied advection
Environmetrics, 28
Han-Lim Choi, J. How (2010)
Continuous trajectory planning of mobile sensors for informative forecasting
Autom., 46
Carl Rasmussen, Christopher Williams (2005)
Gaussian processes for machine learning
Pierre Lermusiaux, C. Chiu, G. Gawarkiewicz, P. Abbot, A. Robinson, Robert Miller, P. Haley, W. Leslie, S. Majumdar, A. Pang, F. Lekien (2006)
Quantifying Uncertainties in Ocean Predictions
Oceanography, 19
Karine Foss, G. Berget, J. Eidsvik (2021)
Using an autonomous underwater vehicle with onboard stochastic advection‐diffusion models to map excursion sets of environmental variables
Environmetrics, 33
KyungMann Kim (2001)
Group Sequential Methods with Applications to Clinical Trials
Controlled Clinical Trials, 22
(2013)
Waste sites from mines in Norwegian fjords

Publisher: Springer Journals
Copyright: Copyright © The Author(s) 2023
ISSN: 0929-5593
eISSN: 1573-7527
DOI: 10.1007/s10514-023-10095-8
Publisher site: See Article on Publisher Site

Abstract

Discharge of mine tailings signiﬁcantly impacts the ecological status of the sea. Methods to efﬁciently monitor the extent of dispersion is essential to protect sensitive areas. By combining underwater robotic sampling with ocean models, we can choose informative sampling sites and adaptively change the robot’s path based on in situ measurements to optimally map the tailings distribution near a seaﬁll. This paper creates a stochastic spatio-temporal proxy model of dispersal dynamics using training data from complex numerical models. The proxy model consists of a spatio-temporal Gaussian process model based on an advection–diffusion stochastic partial differential equation. Informative sampling sites are chosen based on predictions from the proxy model using an objective function favoring areas with high uncertainty and high expected tailings concentrations. A simulation study and data from real-life experiments are presented. Keywords Adaptive sampling · Gaussian process · Stochastic modeling · AUV · Mine tailings 1 Introduction The world’s oceans are heavily impacted by human activities, such as pollution and overﬁshing (Halpern et al., 2008). One example is the harmful effects of mine tailings, drill cuttings, We acknowledge support from the Norwegian Research Council and mud discharge. Along the Norwegian continental shelf, (RCN) through the Centers of Excellence funding scheme, Project Number 223254 - Centre for Autonomous Marine Operations and there are sensitive species like cold-water coral reefs and Systems (NTNU-AMOS), and the INDORSE Project 267793. The sponges which are negatively impacted by such deposition authors would like to thank Tore Mo-Bjørkelund from NTNU and (Trannum et al., 2010), but still, there exist several active Frédéric Py from SINTEF Digital for help related to the ﬁeld mining seaﬁlls and new seaﬁlls are in the planning phase. experiments, and Tor Nordam and Finn Are Michelsen from SINTEF Ocean for running the numerical models SINMOD and DREAM and Improved methods for monitoring the extent of the contam- supplying the data sets used in this paper. ination from seaﬁlls are important to ensure that sensitive areas are protected. In this paper, we use a seaﬁll located in B Gunhild Elisabeth Berget Frænfjorden (see Fig. 1) as a case study where we validate gunhild.berget@ntnu.no our proposed method for monitoring the tailings distribution Jo Eidsvik close to the outlet. jo.eidsvik@ntnu.no When mapping the tailings distribution, models of ocean Morten Omholt Alver dynamics can assist us with helpful information. How- morten.alver@ntnu.no ever, building such models is challenging because of the Tor Arne Johansen oceans’ large-scale nonlinear processes and high spatio- tor.arne.johansen@ntnu.no temporal variability. Existing models use numerical methods Department of Engineering Cybernetics and Centre for to improve their accuracy (Grifﬁes et al., 2000) continuously, Autonomous Marine Operations and Systems, NTNU, but they still have errors due to simpliﬁcations, limitations Trondheim, Norway in numerical implementation, uncertain parameters, and ini- Sintef Community, Trondheim, Norway tial and boundary conditions. Thus, data assimilation using Department of Mathematical Sciences, NTNU, Trondheim, in situ measurements can help improve the model accuracy. Norway 123 Autonomous Robots This function ensures that the area is explored by choosing locations assumed to be information-rich and also consid- ers the travel length limitations of the AUV, leading to an adaptive sampling strategy. In our case study, the method is tested in simulation using the Frænfjorden area (Norway), a fjord containing a seaﬁll for mine tailings. The fjord can be seen in Fig. 1, where the location is marked with a red circle. We compare the proposed method to four simpler strategies in the simulation study. As a proof of concept, real-life experiments are carried out in the same area. Fig. 1 Map of the Frænfjorden area. The pink contours in the ﬁg- The paper is organized as follows: Sect. 2 presents back- ure show the location of the factory Omya Hustadmarmor, and the ground information on the Frænfjorden area, the ocean black lines indicate the pipes of the seaﬁll. The numbers (2.1–2.4) and models, and adaptive sampling. In Sect. 3, the development (1.1–1.2) mark possible outlets for the tailings. At the time of the ﬁeld of the proxy model is presented. Section 4 presents a com- experiment, outlet 2.3 was the operational discharge. The blue rectan- gle represents the operational area for the simulation study and ﬁeld bined state space model explaining how the data assimilation experiments (Color ﬁgure online) is done before the sampling strategies are discussed in Sect. 5. Section 6 contains simulation results comparing our pro- posed method to four simpler strategies, and results from experimental ﬁeld trials are presented in Sect. 7. Section 8 Measurements are commonly obtained using remote sensing, ships, or buoys, but this data is expensive to collect, making discusses the results before Sect. 9 gives a conclusion and it hard to observe the environment in detail. Recently, robotic future work. vehicles like autonomous underwater vehicles (AUVs) have become more robust and less expensive, and by equipping 1.1 Related work and contribution them with sensors, one can collect measurements (sample) more efﬁciently. Still, the possible observations are sparsely With the improved technology and robustness of AUVs, distributed, leading to signiﬁcant gaps in our understanding adaptive sampling using these robotic platforms has had of the ocean (Stewart, 2008). This measurement limitation increased focus lately, and a review of the existing adap- prompts the need for intelligent and targeted observation tive sampling methods for AUV missions can be found in designs, and an important question is when and where to Hwang et al. (2019). Most adaptive sampling methods uti- sample. lize GPs (Cressie & Wikle, 2011), which are widely used to This paper addresses this question and suggests an adap- create non-parametric and computationally efﬁcient models. tive sampling method utilizing prior data from ocean models They are popular not only for adaptive sampling but used in and real-time in situ measurements from an AUV. Because various applications like neuroimaging (Hyun et al., 2016), the already existing ocean dynamics models have a high com- machine learning (Rasmussen & Williams, 2005) and sen- putational load, making them unﬁt for running on embedded sor placement (Krause et al., 2008). Combining GPs with robotic systems, an onboard low-complexity proxy model robotic vehicles is amongst others explored in Zhang et al. is created so that real-time adaption of the mission can be (2012) where an AUV is used to track an upwelling front, obtained. The proxy model represents the state of the ocean and in Luo and Sycara (2018), which explore a method of at the time and can be updated when new information is multiple vehicles using a mixture of GP models. Foss et al. available. We use a spatio-temporal Gaussian process (GP) (2022) use a spatio-temporal model for the same case we are with a kernel describing spatial dependencies and a temporal studying here, but with another objective function relying on model built using an advection–diffusion stochastic partial specifying a threshold for the particle concentration. differential equation (SPDE). The linear SPDE model and This work presents an implementation of an adaptive sam- the Gaussian measurement model are written as a state-space pling method for environmental sensing of mine tailings model, and the Kalman ﬁlter is used to update the state recur- close to a seaﬁll. We show how to utilize prior information sively over time. Two existing numerical models, SINMOD from ocean models to create GP models. Similar approaches, (ocean model) and DREAM (particle dynamics model) are also utilizing prior knowledge to learn GP models, can be used to train the compact GP model creating a proxy model found in Das et al. (2015) who use an AUV to collect samples of the tailing distribution. Sensor readings on the AUV can for ex-situ analysis, selecting the sampling locations based on then be used to update the proxy model onboard in real-time. previous mission data and maximizing a utility function, and An objective function for waypoint selection is presented to in Ma et al. (2018), who create a model of a spatio-temporal maximize the value of information from the measurements. environment by learning and reﬁning the hyperparameters of 123 Autonomous Robots a GP kernel using prior ocean model data. In our approach, discharge for the data in the simulation study and the ﬁeld the ocean model data are used to develop a non-stationary experiments is on the western pipe. This outlet is marked as correlation kernel as in Berget et al. (2018) and Fossum et al. 2.3 in the ﬁgure. (2018) and set the parameters in a temporal model based on the stochastic advection–diffusion partial differential equa- 2.2 Ocean models tion (Sigrist et al., 2015). This allows the model to capture critical dynamic phenomena such as tides. Estimated solu- Numerical oceanographic models are used in this paper to tions of the SPDE can be written in a linear form so that the build a prior belief of the state of the ocean and to predict Kalman ﬁlter can be used to update the model state. An objec- the movement of particles with the currents. Ocean mod- tive function based on variance reduction, hotspot detection, els provide information on salinity, currents, temperature, and AUV operation time constraints is used; see also Berget density, and pressure and describe the state of the ocean at et al. (2018). The approach sequentially selects waypoints a given time. The models use a set of hydrodynamic and based on the current model state, and the method is tested in thermodynamic equations called the primitive equations, and both simulations and ﬁeld experiments near a mining seaﬁll these are commonly solved numerically using ﬁnite differ- in a Norwegian fjord. ence methods. This makes the models computer intensive; The main contributions of this paper are: hence, high-resolution predictions can only be computed for limited areas, depending on the available computer resources. – Propose an adaptive sampling method for environmental Input to ocean models typically includes tides, sea-level sensing of mine tailings. pressure, wind, heat exchange, bathymetry, and freshwater – Use ocean model data to set the parameters in a lin- runoff. Ocean model outputs are uncertain due to limits in ear temporal model based on a stochastic advection– model resolution, limits to our knowledge of the processes diffusion partial differential equation. resolved by the model, errors in initial values, boundary con- – Present results from real-life experiments using the pro- ditions and parameter values, and inaccuracies in numerical posed adaptive sampling method. implementations (Lermusiaux et al., 2006). It is, therefore, necessary to validate the output against observations, and data assimilation can be used to correct the model state esti- 2 Background mates. For monitoring of releases of mine tailings or other particles or dissolved material, there is a need to develop enabling technology for efﬁcient and targeted ocean sam- This section presents deﬁnitions and background informa- tion on ocean models, the spatio-temporal model, and data pling. Observations from different platforms such as buoys, assimilation. ship-based sampling, or robotic sampling from, for example, AUVs can be used to evaluate or, through data assimilation, 2.1 The Frænfjorden area improve the model estimates. In this paper, we use data from two speciﬁc models to In Norway, several active and planned seaﬁlls have impli- build the stochastic onboard model, SINMOD (ocean dynam- cations for the marine environment (Ramirez-Llodra et al., ics model) (Slagstad & McClimans, 2005; Wassmann et al., 2015; Morello et al., 2016). The tailings typically consist of 2006) and DREAM (particle transportation model) (Rye et particulates and residual chemicals from processing, and this al., 1998, 2008). Nepstad et al. (2020) described the SIN- waste will smother the sea bed close to the sea ﬁll. In addi- MOD and DREAM setups for Frænfjorden and evaluated the models against ﬁeld measurements. The relation between tion, some of the waste may be transported by ocean currents away from the outlet, thus affecting a larger area of the sea SINMOD, DREAM, and the onboard model can be seen in bed and the water column. Fig. 2. More details on the two ocean models are given below. As the case study for this paper, we focus on a seaﬁll located in Frænfjorden, Norway (see Fig. 1). A calcium car- 2.2.1 SINMOD: ocean dynamics model bonate processing plant in the factory has been operating for several decades. This processing of ore from mining causes SINMOD is a numerical ocean model system that connects ﬁne-ground waste (tailings) (Kvassnes & Iversen, 2013), dis- and simulates physical and biological processes in the ocean posed of on the ocean ﬂoor, in a seaﬁll close to the factory. (Slagstad & McClimans, 2005). The model is based on The outlets in Frænfjorden are located at approximately the Navier–Stokes equations and uses a nesting technique 20 m depth, and tailings are continuously discharged with a where high resolution models obtain their boundary con- time-varying discharge rate. A typical discharge rate from the ditions from larger model domains with lower resolution. outlets is around 400 m /h. The discharge alternates between SINMOD is established in conﬁgurations with horizontal two pipes, illustrated by black lines in Fig. 1. The operational resolutions ranging from 20 km to 32 m. Input to the model 123 Autonomous Robots the onboard proxy model. This data is also used to estimate parameters in the SPDE model, including the GP kernel. 2.3 Adaptive sampling Observing the ocean by taking measurements or samples is crucial to increasing our understanding of ocean processes. Ocean sampling can be seen as the action of sensing or mea- suring one or more physical properties in the ocean at a speciﬁc point in space and time. Traditional methods include mounting sensors to buoys which limits the sampling to only one location, or immersing instruments in the water from boats, commonly using a predetermined plan of sampling locations. Exploiting robotic platforms and especially AUVs can intensify ocean sampling, giving more data and better insight. We need more sophisticated sampling strategies to exploit the beneﬁts that robotic platforms have regarding autonomy and increased mobility. To obtain the most scientiﬁcally rel- evant measurements effectively, adaptive approaches need Fig. 2 The relation between SINMOD, DREAM, and the adaptive sam- to be developed. In contrast to static/pre-scripted schemes, pling framework. SINMOD provides ocean currents data as input to the adaptive/data-driven strategies can react to measured condi- particle transportation model DREAM. Both SINMOD and DREAM tions, having access to both prior and in situ information; provides data used to build the prior belief of the onboard model. Ocean selecting sampling locations depending on past observations current predictions are used in the temporal part of the model and parti- cle concentration predictions on the boundary are used in the updating taken during exploration. of the model The beneﬁts of adaptive sampling will depend on the situ- ation. There are numerous studies on the value of additional (sequential) sampling efforts in different contexts, say clin- includes atmospheric forcing, freshwater runoff, and bound- ical trials in medical research (Jennison & Turnbull, 1999) ary conditions (temperature, salinity, tidal forcing, currents). or geoscience data (Eidsvik et al., 2015). In our case, there Ocean dynamics data from SINMOD is used in this case would be limited to gain by adaptation if a priori information study as input to DREAM and the onboard proxy model (see makes it possible to pre-script a high-value AUV plan before Fig. 2). Atmospheric forcing for daily forecast simulations deployment. However, if there is high prior variability, the was obtained from NOAA’s Global Forecast System with ability to adapt the initial plan one way or the other, based on 0.5 resolution (for large-scale model domains) and from onboard measurements, would bring added value by guiding the Norwegian Meteorological Institute’s MetCoOp Ensem- the AUV to more relevant places. ble Prediction System with 2.5 km resolution (for local model domains). 3 Modeling 2.2.2 DREAM: particle transportation model This section describes the advection–diffusion SPDE model DREAM is a Lagrangian particle transport model which and its associated spatio-temporal GP model representation. can be used to simulate the behavior and fate of marine We also outline parameters in the model. pollutants, including particulate discharges from drilling operations (Rye et al., 1998, 2008). It provides time series 3.1 Spatio-temporal onboard model of concentrations of released materials in the water column and deposition of these materials onto the sea ﬂoor. Input to To select optimal sampling locations in a dynamic domain, the DREAM model includes hydrodynamic data delivered by information about the spatial conditions at all times is essen- SINMOD and information about the release (amount, rate, tial. Since the computational demands using ocean models densities, grain size distribution). DREAM is often used as a such as SINMOD and DREAM are too high for real-time decision support tool for managing operational discharges to onboard modeling, we suggest a stochastic spatio-temporal the marine environment. This paper uses data from DREAM proxy model based on an advection–diffusion SPDE. When to set the initial belief for the mine tailings distribution on this SPDE is discretized, it represents a spatio-temporal GP 123 Autonomous Robots model. Apart from having a smaller computational footprint, B(t , s) describes external information, e.g. boundary infor- GPs are conventional tools for statistical modeling of spatial mation. data and have been widely adopted in oceanographic appli- This paper assumes horizontal, isotropic diffusion lead- cations (Binney et al., 2013). In this paper, the proxy model ing to a constant diffusion parameter D. We denote the t t T is used to approximate the distribution of the particle con- 2-dimensional drift vector as c(t , s) = (c , c ) at time step e n centration near the seaﬁll in Frænfjorden. t. Then we get There exist various types of spatio-temporal GP mod- ∂ X (t , s) ∂ X (t , s) ∂ X (t , s) els, see e.g. Chapter 6 of Cressie and Wikle (2011)for a t t = λX (t , s) − c + c e n ∂t ∂s ∂s review. For statisticians, it has been important to understand e n 2 2 the spatio-temporal covariance function of such processes, ∂ X (t , s) ∂ X (t , s) + D + + B(t , s) + (t , s). which is crucial in spatio-temporal Kriging. In the simplest 2 2 ∂s ∂s e n case, the covariance function is separable in space and time, (2) but there has also been much work on non-separable covari- ance functions (Gneiting, 2002). Stochastic differential and The SPDE is discretized using the upwind differencing difference equations in the spatial domain ﬁt naturally into scheme to ﬁnd an approximate solution. We deﬁne the length state-space models, which enable spatio-temporal Kalman of a time step as dt. Forward differences are used in time, ﬁltering using in-situ observations. The spatial version of the giving autoregressive process is one such example of a linear dif- ∂ X (t + dt , s) − X (t , s) ference equation. So are integro-difference equations which X (t , s) ≈ . (3) weight spatial variables from the previous time step to pre- ∂t dt dict the spatial variables at the current time step (Storvik et We deﬁne a regular grid D with a total of N grid points al., 2002). Recently, extensions of the advection–diffusion with resolution ds in the east direction and ds in the north e n equation have been developed (Sigrist et al., 2015; Richard- direction. The set of spatial locations (east, north) is then son, 2017). These models provide more interpretability in s = (s , s ) ∈ D ⊂ R , i = 1,..., N.For theﬁrst i e,i n,i s the spatial weighting as well as opportunities for embedding derivative in space, we alternate between forward and back- time-dependent drift. We use one of these models in the cur- ward differences based on the direction of the drift vector c . rent paper. The advection–diffusion equation is also endorsed In the eastward direction, we have in the physical modeling community, and in our setting the advection can be speciﬁed from the tidal currents and the ∂ X (t , s ± (ds , 0) ) − X (t , s) X (t , s) ≈ , (4) diffusion coefﬁcients can be set from the ocean models SIN- ∂s ds e e MOD and DREAM. We consider a real-valued stochastic process {X (t , s), t ≥ t where the upper sign is used for c < 0 and the lower is used 0, s ∈ D}, representing the uncertain particle concentration t when c ≥ 0. Similarly, in the northward direction, we have on log scale in (east,north) location s = (s , s ) ∈ D ⊂ R e n at time t. T ∂ X (t , s ± (0, ds ) ) − X (t , s) X (t , s) ≈ , (5) The spatio-temporal process of transport of a medium in ∂s ds n n a ﬂuid can be expressed by the advection–diffusion SPDE (Sigrist et al., 2015) where the upper version is used for c < 0 and the lower is used when c ≥ 0. The second derivatives are discretized as ∂ X (t , s) follows, = λX (t , s) − c(t , s) ∇ X (t , s) ∂t +∇ D∇ X (t , s) + B(t , s) + (t , s), (1) X (t , s) ≈ ∂s T T X (t , s + (ds , 0) ) − 2X (t , s) + X (t , s − (ds , 0) ) T e e with c(t , s) = (c (t , s), c (t , s)) is the drift vector ﬁeld for e n (6) ds advection, D the diffusion matrix, (t , s) a noise term, ∇ e ∂ ∂ the gradient operator , and λ ∈[−1, 0] a damping ∂s ∂s e n X (t , s) ≈ ∂s constant controlling the autoregressive relationship between T T state vectors (Richardson, 2017). The noise term (t , s) is X (t , s + (0, ds ) ) − 2X (t , s) + X (t , s − (0, ds ) ) n n . (7) assumed to be a zero-mean Gaussian noise term that is 2 ds uncorrelated in the temporal domain but spatially correlated, Corr((t , s), (t , s )) = R(|s − s |), where R() represents By scaling the time discretization, a change in time can be a correlation function for distance input. The source term denoted t +1 instead of t +dt.Wehave {t = t + j ×dt ; j = 123 Autonomous Robots 0, 1, 2 ..., T }. When inserting all the discretized terms into The prior covariance matrix Σ is then given by eq. (1), the scheme can be written as a linear process on the ⎡ ⎤ form Σ Σ ... Σ 11 12 1N ⎢ ⎥ Σ Σ ... Σ 21 22 2N ⎢ ⎥ Σ = , (10) ⎢ ⎥ 0 . . . . . . . ⎣ ⎦ . . . X = A X + B + , ∼ N (0, W ), t = 1,..., T , t t t −1 t t t t Σ Σ ... Σ N 1 N 2 NN (8) where Σ = σ σ R and R represents the spatial cor- ij i j ij ij relation. Hence, the diagonal of the covariance matrix in where X is the size N ×1 vector of particle concentrations at equation (10) contains the location-wise variances σ , and time step t for all spatial locations, A is the N × N propaga- the off-diagonal elements describe the covariance between tor matrix including information about advection, diffusion the locations. The covariance function is chosen by com- and damping, B is a matrix adding external information, in paring the empirical covariance of the ocean model training our case boundary conditions, and W is a positive deﬁnite data with known valid parametric covariance functions. The variance-covariance matrix for the noise term . Matern (3/2) kernel (Matérn, 2013) is ﬁtted to such covari- The particle concentrations at time t = 0 areassumedtobe ance function here; Gaussian distributed with mean vector μ =[μ ,...,μ ] 1 N and N × N variance-covariance matrix Σ . For the boundary conditions in the SPDE, we use Dirichlet R = (1 + φh ) exp(−φh ), (11) ij ij ij BCs on the east side of our survey area. This can be consid- ered as a constant forcing along the east side and is chosen where h =|s − s | is the Euclidean distance between two ij i j because of the constant ﬂow of mine tailings that are let locations s and s and φ> 0 is a constant meta-parameter i j out on this side of the rectangle. We assume no drift across regulating the correlation decay with the distance. The best the boundary on the remaining sides of the survey area and value for φ is speciﬁed using training data. therefore use Neumann BCs. The boundary conditions are set To specify covariance model parameters and visualize using predicted particle concentration data from DREAM. how the spatial correlation decays as a function of the distance between points in the domain, it is common in geostatistics to use the variogram of the data, see e.g. 3.2 Parameter specification in the SPDE model Cressie and Wikle (2011). The variogram works on differ- ence X (s ) − X (s ) between spatial data. Doing so, it tends t j t i The advection and diffusion parameters in the SPDE (1)must to be more robust to background trends than the covariance be speciﬁed. The velocity ﬁeld is set from physical assump- function. The theoretical relation between the variogram γ() t t T tions, where c (s) = (c (s), c (s)) as the forcing of the e n and the covariance is convection in our case is treated as the currents in the ocean. The values are set using predicted currents from the ocean γ(s − s ) = Var(X (s ) − X (s )) = i j t i t j model SINMOD. The isotropic diffusion D is settobethe same as for the computations in the DREAM model D = 0.1 (Var(X (s )) + Var(X (s )) − 2Cov(X (s ), X (s ))). t i t j t i t j m /s. The initial distribution is mainly estimated from the (12) DREAM data. Since ocean processes are affected by cur- rents, wind patterns, bathymetry, and freshwater runoff, we Figure 3 shows results from the variogram analysis. The will have elevated variability in some locations. Thus, we empirical variogram (blue) is computed by averaging the choose a non-stationary prior model (Jun & Stein, 2008) square differences of variable pairs having spatial distance where the mean and variance in each location are speciﬁed within bins. We use regular bins of distances (ﬁrst axis). empirically from the training data available by DREAM: Three different theoretical variograms were ﬁtted to the train- ing data. All of them ﬁt the empirical variogram adequately, but the Matern variogram (green) was chosen here because 2 2 it is closer to the empirical variogram for very small spa- σ = (y − μ ) . (9) i ,m i M − 1 tial distances. This indicates that a Matern spatial ﬁeld has a m=1 smoothness that ﬁts the DREAM data better. We see that the variogram increases with distance and then ﬂattens out. Here, Here, μ = y is the empirical mean concentra- the effective spatial correlation length is about 500 m. For the i i ,m m=1 tion in location s over all M realizations of training data. Matern correlation function in equation (11) this means that 123 Autonomous Robots represented by the initial onboard model, and observations are taken from sensors on the AUV. We deﬁne our observation model as Y = G X + ω , ω ∼ N (0,τ I ), (13) t t t t t where Y is a length m vector of observations at time t, G is t t a time-varying observation matrix of size m × N that maps the state variables to the observations, and ω represent the observational error which is modeled as a zero-mean Gaus- sian with standard deviation τ , and I representing the identity matrix. Each row of G has a one-entry on the vector position associated with an observed grid cell, and the other entries are zero. This method constraints the sampling locations to Fig. 3 The empirical variogram is shown in blue, and a Matern(3/2), an exponential, and a Gaussian variogram is ﬁtted to the data. We ﬁnd be at the predeﬁned grid cells. Continuous path sampling, as an average variance of 0.6 and an effective spatial correlation length of in Choi and How (2010), could improve the sampling method about 500 m (Color ﬁgure online) and should be considered in future work. Since the SPDE model deﬁnes a GP with a linear process representation, the spatial Kalman ﬁlter can be used to update φ ≈ 0.01. The variogram converges to the average variance our state at every time step. Letting Z = (Y ,..., Y ) rep- t 1 t which is about 0.6. resent all observations and prior information up to time t, The advection–diffusion model parameters are further set the ﬁltering distribution is X |Z ∼ N (μ , Σ ) with mean t t t |t t |t from the physical assumptions in the ocean model. From μ and variance-covariance Σ deﬁned by t |t t |t visual inspection, this GP model appears to give reasonable simulations of the ﬁeld but with more patchy randomness Σ = A Σ A + W than what is seen in the ocean model. As commonly seen t |t −1 t t −1|t −1 t T 2 in physical models data, the variogram computed from the Q = G Σ G + τ I t t |t −1 t t ocean models data starts to ﬂuctuate after a while. In many T −1 K = Σ G Q t t |t −1 t t cases this can be related to physical principles. Notably, the μ = A μ + B t t GP model is a statistical proxy model that can be run onboard t |t −1 t −1|t −1 but is potentially missing some large-scale oceanographic μ = μ + K (Y − G μ ) t t t t |t t |t −1 t |t −1 features. The ocean model is a numerical description with Σ = Σ − K Q Σ . (14) t |t t |t −1 t t |t −1 much skill but potentially getting systematic bias due to mis- speciﬁed inputs growing with time. They both have their This recursion starts by the initial distribution with mean μ merits. and Σ , where X ∼ N (μ , Σ ). The boundary conditions 0 0 0 The covariance matrix W in the noise term of the equation and the model parameters going into the matrices and vectors is set to have the same correlation structure as the covariance in equation (14) are all speciﬁed as described in Sect. 3.2. matrix for the initial state Σ ,giving W = V Σ , where V 0 t 0 is a constant scaling parameter. Considering this in relation to the original continuous-time SPDE, the discretization in 5 Sampling strategy time would also impact this variance, but in our work, we use a ﬁxed time discretization identical to the one used onboard. An objective function is suggested to obtain an informative There does not seem to be any damping over time in the ocean path for the AUV. The function is created based on three model data, but for stability reasons (Richardson, 2017)we criteria, as in Berget et al. (2019): specify a minimal level λ =−0.0010. 1. Locations with high variance are preferred 2. Locations with high predicted concentration are pre- 4 State space model and data assimilation ferred 3. Locations leading to a suitable travel length within the The state of concentrations is assumed to be modeled by next time interval for the AUV are preferred. equation (8). The aim of data assimilation is to combine observations with this prior process knowledge to obtain an optimal sequential posterior or ﬁltering estimate of the cur- The ﬁrst criterion is chosen because sampling in areas rent state (Sar, 2013). In this study, the prior knowledge is with a high model variance will reduce the total variance, 123 Autonomous Robots hence creating a more accurate model. This criterion also egy. Making the AUV ﬂy with the currents may increase ensures that the AUV travels to unexplored areas. The sec- endurance and reduce travel time between waypoints. This ond criterion makes the method adaptive. When studying the topic is amongst others explored in Pereira et al. (2013) and simulation results of the particle transport from the complex should be considered in future work. model, it is clear that the variability is highest where there is a high concentration of particles. Hence, the criterion is inspired by this observation and assumes that locations with 6 Summary of the suggested approach high predicted concentration will be rich with information. The last criterion comes from the travel length limitations of We now provide an overview of the suggested approach. At the AUV. When choosing the next sample location, it is also the beginning of the deployment, the GP model is initialized. essential that the AUV does not travel too far, assuring that The following steps are executed: the risk of the vehicle drifting outside the survey area is kept low. – Kernel hyperparameters are estimated using prediction The suggested objective function is then created by having and hindcast data from the particle transportation model a term for the ﬁrst two criteria. At time step t for location s , DREAM. (Described in Sect. 3.2.) the objective function is given by – The initial GP state deﬁned by μ and Σ are found. (Described in Sect. 3.2.) f (s ) = θ σ + θ μ , (15) t i 1 2 i ,t |t – Ocean current prediction data from SINMOD is used to i ,t |t deﬁne the advection ﬁelds c , which is used together with where the constant parameters θ =[θ ,θ ] deﬁne the weight- 1 2 the diffusion constant D to deﬁne transition matrices A . ing for each criterion. Deciding the weighting of criteria 1 (Described in Sect. 3.1.) and 2 depends on design choices. When the goal is to obtain – Particle transportation data from DREAM deﬁnes the a good map of the overall region, then criterion 1 should be boundaries that are incorporated in the vector B . dominating. However, if hotspots are the main focus, crite- (Described in Sect. 3.1.) rion 2 should dominate. In our case, the parameters θ are tuned in simulation to obtain a good balance between the The sampling method starts out either with the initial two criteria. GP model or with the latest updated GP model. The AUV The chosen sampling location s at time step t is deter- t performs adaptive sampling choosing the next waypoint as mined as the location that maximizes the objective function described in Sect. 5. The onboard model of the AUV is then f (s ) for s ∈[s ,..., s ]. Given the previous sampling t i i 1 N updated; the ﬁeld is moved with the currents, and new infor- location s , we choose the next sampling location as t −1 mation from observations is added. Updating is executed whenever a waypoint is reached, using data from observa- s = arg max f (s ) t i tions taken along the AUV transect in the assimilation. An overview of the adaptive sampling method is given in Algo- s.t. |s − s |≥ d i min t −1 rithm 1. |s − s |≤ d . (16) i max t −1 Algorithm 1 Adaptive Sampling method The constraints on the objective function ensure that the third Initialize GP model criterion in the objective function is adhered to by choosing while running do d and d such that a suitable travel length is obtained. min max Select new waypoint using (16) |s − s | is the Euclidean distance between the previous i Travel to new waypoint and collect observations t −1 sampling location and s . Equation (16) is solved using brute- if waypoint is reached then Update GP model with observations using (14) force enumeration; the value of the objective function (15) is found for all the grid cells within the region deﬁned by the constraints in (16). The next sampling waypoint s is chosen as the grid cell with the highest objective value. In this paper, ocean currents are included in the spatio-temporal proxy model, but they are not considered in the AUV nav- 7 Simulation study igation plans for selecting sampling locations. Rather, the AUV is set to surface often to avoid risk and reduce navi- A simulation study compares the proposed model with four gation uncertainty, as described in the experimental section. simpler strategies. We vary both the process model and In a larger-scale mission, it would certainly be beneﬁcial to the sampling strategy. The temporal model based on the include more of the ocean currents in the sampling strat- advection–diffusion SPDE is compared to a stationary pro- 123 Autonomous Robots Table 1 Overview of the ﬁve methods used for comparison in the sim- 7.1 Results and discussion ulation study In this section, we present results from the simulation study Method Process model Sampling strategy and discuss the performance of the methods. SPDEObj SPDE Objective SPDELawn SPDE Lawnmower 7.1.1 Performance measures SPDE SPDE No sampling Obj Stationary Objective To evaluate the different methods, three performance mea- Lawn Stationary Lawnmower sures are studied. The mean absolute error (MAE) is com- puted as the mean of the absolute error between the predicted concentration μ and the true ﬁeld from DREAM (April t |t 2nd) for all the grid cells and time steps. The mean standard Table 2 Performance measures for the ﬁve methods deviation (MSD) is computed similarly but then focuses on Method MAE MSD MO the standard deviation in the proxy model in all grid cells at SPDEObj 0.76 2.38 209.3 all time steps. Lastly, we compute the mean objective value SPDELawn 0.87 2.67 272.8 (see Eq. (15)) (MO) in all grid cells at every time step. Table SPDE 1.11 4.70 848.1 2 shows these performance measures computed as the mean Obj 0.85 2.52 230.2 over the total time steps. Figure 4 shows how the mean objec- tive value (MO) (see Eq. (15)) evolves through time, plotting Lawn 0.99 2.67 280.5 the mean value for every time step. The table shows the mean absolute error (MAE), the mean standard Considering the MAE values in Table 2, we see that deviation (MSD), and the mean objective function value (MO) SPDEObj performs best in this case. Including a non- stationary process model increases the method’s perfor- mance, as this is the case using both the Obj sampling strategy cess model where the next time step is equal to the previous and the Lawn sampling strategy. The objective function also time step, as in Berget et al. (2018). Different sampling increases the method’s performance compared to the pre- strategies are tested with sampling strategies based on the planned lawnmower strategy. Considering the MSD and MO suggested objective function, a predeﬁned sampling strat- measures, the SPDE method stands out with signiﬁcantly egy following a lawnmower pattern, and a method that uses higher values than the other methods. This is due to the only the SPDE model with no sampling (no observations are taken). The SPDE and Stationary process model are the same for all sampling strategies using the same parameter values given in the text. For the strategy using the objective function we choose the travel distance restrictions d = 180 meters min and d = 220 meters and set the time between each sample max to be 5 min. This gives a suitable travel speed for the AUV, and a constant time step simpliﬁes the updating of the tempo- ral model. When using the lawnmower pattern, the distance between each sampling location is set to 200 ms, and the same time step of 5 min is used. The parameters in the objec- tive function are set to θ = 15 and θ = 23. The different 1 2 strategies are shown in Table 1. Ocean model data from two consecutive days, 1st and 2nd of April 2013, is used in the simulation study. The data from April 1st is used to train the SPDE parameters as described in sections 3 and 4. Data from April 2nd is used as the test data, meaning that the predicted ﬁeld on April 2nd is our ground truth and the ﬁeld from which observations (samples) are drawn. We use realizations of the GP prior as the initial belief and run 100 replicate runs with the presented sampling methods in Table 1. We set the measurement error standard Fig. 4 The mean objective value (MO) (see Eq. (15)) computed over all deviation τ = 0.5 log(ppb). We simulate for a total of 6 h, the grid cells for each of the methods. Results are shown for the whole simulation time period spanning from 15:00 to 21:00 aiming to map the ﬁeld on April 2nd from 15:00 to 21:00. 123 Autonomous Robots Fig. 5 Hourly plot of the model standard deviation at the given time sampling location at the current time. The arrow indicates the direction together with the path for the different sampling strategies up until that of the AUV (Color ﬁgure online) time. The blue dot shows the path’s starting point, and the yellow is the method doing no sampling at all, and hence there is no way for bilize at a mean objective value between 200 and 300. For this method to decrease the proxy models’ variance, affecting both lawnmower methods the standard deviation and the objective values. (SPDELawn and Lawn), periodical ﬂuctuations in the objec- The effect of increasing variance with no sampling can tive function can be observed. This occurs when the sampling also be seen in Fig. 4. Again, the SPDE method stands out is done close to the east end of the survey area, where the with steadily increasing values while the other methods sta- assimilation with new data affects a minor part of the modeled ﬁeld compared to adding information in the center. Gener- 123 Autonomous Robots ally, the Obj and SPDEObj methods have a lower overall mean objective function compared with the other methods, indicating an intelligent sampling design. In Fig. 4, we see that the Lawn method has a lower MO value at some time steps compared to the SPDELawn method. This occurs at times when both Lawn approaches get high MO. Using the lawnmower pattern for sampling (see Fig. 5a) is the reason for this unexpected behavior. Sam- pling is performed in boxes, and some locations in the middle of the boxes get little inﬂuence from measurements. The SPDELawn method solves this by having the temporal model move neighboring sample values into this box. The Lawn method does not have a solution for this, and the values inside this box will be similar to the initial values. Since the initial mine tailings concentration in our example (see Fig. 7a) has low values all over the area, the tailing concentration inside this box will stay low throughout the whole simulation for the Lawn method. This will affect the MO value, leading to Fig. 6 The distribution in one single grid cell using the SPDEObj artiﬁcially low values, even though the performance (here, method. Showing mean (blue line) and 95% conﬁdence interval (red better illustrated by the MAE) can be poor. dotted lines) together with the true value (green line). The selected grid cell is located 340 ms east and 200 ms north in the survey area (Color ﬁgure online) 7.1.2 Sampling design and ﬁeld predictions Figure 5 shows the AUV path and waypoints for the different sampling strategies plotted on top of the standard deviation at every hour of the simulation study. The same lawnmower pattern is used for SPDELawn and Lawn and can be seen in out the whole mission, relying on observations to change its the ﬁrst column of the ﬁgure. The path chosen by the Obj predicted ﬁeld. This is not the case with SPDEObj, which has and the SPDEObj method can be seen in respectively the a non-stationary process model that moves the ﬁeld with the second and third columns. All the methods cover the area currents. Although data from two different days are used as quite well, but both Obj and SPDEObj focus on the eastern training and test data, both datasets are from the same ocean area. This region is closer to the outlet, and the true mine models. This may cause the comparison with ground-truth tailings distribution (seen in the ﬁrst column of Fig. 7)shows to be optimistic. Still, based on the simulation results, it is elevated concentrations in this area. Hence, more frequent reasonable to assume that SPDEObj can map the ﬁeld quite sampling here is expected because of the proposed objective well. The MAE scores of the two methods are indicated for function (see Eq. (15)). different times (Fig. 7). With the relatively high dynamical Figure 6 gives insight into how the GP methods behave, uncertainty, the MAE increases with time. The MAE is gen- showing the distribution in one single grid cell over time erally lower for the SPDEObj approach. using SPDEObj. The blue line shows the method’s mean Simulation studies with different AUV velocities have value, the green line is the true value, and the two red dotted shown that the SPDEObj method’s advantages are higher lines are the 95% conﬁdence interval. The point we are con- when the AUV takes fewer samples in the survey area over sidering is 340 ms east and 200 ms north in the survey area. the same time period. Naturally, lower AUV velocities led When sampling and assimilation are done close to this point, to fewer samples (comparable to having a bigger survey we see a reduction in the conﬁdence interval, and the mean area), while higher velocities led to an increase in the num- value tends to be corrected towards the actual value. ber of samples. When doubling the number of samples, the In Fig. 7 we focus on the two methods Obj and SPDEObj. difference in the performance between the approaches was The updated model mean for every hour is plotted in columns minimal, making it hard to decide which method was best. 2 and 3, and the true ﬁeld (the ﬁeld where the observations However, when operating with half the number of samples are drawn from) can be seen in the ﬁrst column. Comparing (compared to the simulation study above), the differences the methods with the true ﬁeld, they capture the trend quite became more evident, and the SPDEObj method showed well. The display illustrates an important difference between improved results. Hence, it is reasonable to think that the the two methods; because of its stationary process model, difference between the methods will increase with a larger the Obj method inherits the structure from the GP through- survey area. 123 Autonomous Robots Fig. 7 Hourly plot of the “true” value and the model mean value for the Obj and SPDEObj method. The MAE scores are included in the captions 7.1.3 Adaptive versus pre-planned designs pling design for one hour before the end time is shown in the variance plots. We now give an example comparing an the adaptive approach In the example, the adaptive approach gives a better esti- (SPDEObj) with a pre-planned method (SPDELawn) (see mate of the true ﬁeld than the pre-planned method. From the Fig. 8). The example shows the predicted results at the end true ﬁeld, we observe increased concentrations on the eastern of a mission run of total 5 h. The true ﬁeld is seen in Fig. 8a. side. The adaptive approach captures this information, and Figure 8b, c show the model mean and variance for the adap- because of elevated values in the area, continues to investi- tive approach, and Fig. 8d, e present the model state for the gate the eastern side. The pre-planned lawnmower approach, pre-planned approach. The AUV path describing the sam- on the other hand, is set to sample in the western area at the 123 Autonomous Robots Fig. 8 Example run comparing a pre-planned approach with an adap- approach, and Fig. 8d, e present the model state for the pre-planned tive approach. The true ﬁeld after ﬁve hours of sampling can be seen in approach. The AUV path showing the sampling design for one hour Fig. 8a. Figure 8b, c show the model mean and variance for the adaptive before the current time is shown in the variance plots time, and the samples taken are not as informative. Adap- MOD models, will ﬁnd high concentrations of mine tailings. tive approaches excel when conditions are not as expected. This is veriﬁed by our sensor data as well. Figure 11 shows In such cases, the ﬂexibility of an adaptive approach is ben- one transect (waypoint-to-waypoint) for the AUV including eﬁcial since the sampling design is automatically adjusted. sensor data values. It is clear from the ﬁgure that a layer with A pre-planned approach will always stick to the pre-scripted increased turbidity is reached. sampling design. 8.1 Robotic platform and implementation 8 Experimental results The robotic platform in our experiments consisted of a Light AUV from OceanScan (see Fig. 9) equipped with a Wetlabs As a proof of concept, ﬁeld experiments were conducted in EcoPuck sensor measuring total suspended matter (TSM), Frænfjorden. Two experiments were carried out on Decem- which is a measure of optical back-scattering in the form ber 17th and 18th, 2020. The operational discharge outlet at of an attenuation coefﬁcient. The measured sensor values the time of the experiments is marked as 2.1 in Fig. 1, and of TSM indicate turbidity, and this observation value was a rectangular operational area of size 1150 m × 250 m (also used to assimilate our model. It was converted to a mea- marked in the ﬁgure) was used. We focus on the depth layer sure of (log) concentration to adjust the value to our onboard at around 22–25 ms since this is a layer with a lot of variation model. This conversion is not straightforward, and we need and a layer where we, according to the DREAM and SIN- to determine the uncertainty τ . The Wetlabs EcoPuck sen- 123 Autonomous Robots a series of actions (e.g., Goto and Arrive_at). The set of all those actions forms a plan built while ensuring operational constraints of the mission (e.g., the vehicle should never dive deeper than a certain depth or leave a deﬁned area.) The plan is then sent to DUNE, which handles low-level control and execution. Before deploying the AUV, our method was tested in a simulated environment similar to our embedded system. In simulations, sensor readings from the AUV were replaced with data from SINMOD and DREAM, and an AUV simula- tor in DUNE was used to simulate the behavior of the AUV. Fig. 9 The light AUV (Harald) produced by OceanScan was used in The layout of the simulator is shown in Fig. 10. our experiments. The length of the AUV is 240 cm, and it is equipped with multiple sensors like CTD, Dissolved Oxygen, and Fluorometer 8.2 AUV path Inside the operational area, a 2-dimensional grid of size 58 × 12 was deﬁned, giving N = 684. The distance between each node was approximately 20 m in both grid directions, and each grid point was a possible waypoint for the AUV. The AUV was set to surface whenever reaching a new way- point to reduce the uncertainty caused by navigation errors and vehicle drift and for safety reasons. The AUV did a yoyo maneuver between 20 and 26 ms between each way- Fig. 10 The architecture of the embedded system (on top) and the point. As in the simulation study, we set d = 180 meters min simulator (bottom). The Uniﬁed Navigation Environment (DUNE) runs and d = 220 meters. The AUV speed is set to 1.2m/s, max onboard the AUV, and the agent architecture T-REX sits on top of this, and as opposed to the simulation study, where a constant allowing the embedding of multiple decision processes such as planning time of 5 min between two waypoints was used, we now let the time vary. The time used between two waypoints is found using the start and end times of the AUV transect. The temporal model and the updating are done when a way- sor user manual suggests an attenuation error of about 4 % −1 at attenuation coefﬁcient 1 m . With concentrations varying point is reached, updating the model with the amount of time between approximately 3 and 8 log(ppb), assuming measure- used for that particular transect. Figure 11 shows one tran- ment standard deviation of 4% gives τ between 0.12 and 0.32 sect (waypoint-to-waypoint) from our survey, starting with a log(ppb). However, both navigation errors and vehicle drift dive until reaching the selected depth layer, doing the yoyo can lead to uncertainty in the sampling location. Hence, we maneuver, and then surfacing in a corkscrew maneuver. The set τ = 0.5 log(ppb) for the experiments. corkscrew maneuver was chosen to obtain a sampling value Figure 10 shows the layout of the agent architecture used as close as possible to the waypoint. The yoyo maneuver was in the experimental survey. The Uniﬁed Navigation Envi- chosen to increase the possibility of ﬁnding the layer of ele- ronment (DUNE) (Pinto et al., 2012) runs onboard the AUV vated mine tailings concentration. The sensor reading used and is used for control, navigation, vehicle supervision, com- for updating the model was chosen as the maximum value in the layer 22–25 ms (the layer marked with red surfaces) for munication, and interaction with actuators. On top of this sits the agent architecture T-REX (Teleo-Reactive EXecu- each of the separate yoyos. For the transect in Fig. 11,this tive) (Rajan & Py, 2012), which enables an adaptive mission. would result in 6 updating values. Because the transect length T-REX allows the embedding of multiple complex decision could vary for different waypoints, the number of yoyos was processes (including planning). The communication between not ﬁxed but mostly alternated between 5 or 6 yoyos per DUNE and T-REX was handled by the LSTS toolchain (Pinto transect, giving 5 or 6 updating values per transect. et al., 2013), which provides the back-seat driver API to DUNE. This allows external controllers, such as T-REX, 8.3 Results to provide desired tasks for our platform while receiving progress updates on the current state. Our sampling algo- The survey area is located near the coast in a fjord arm, so rithm was written as a python script and was implemented the currents are highly inﬂuenced by the tidal cycle. To get as a reactor in the T-REX framework. A reactor is a com- a better insight into how the ocean currents behave at the ponent of T-REX acting as an internal control loop and is time of the missions, we consider the tidal cycle at the time capable of producing goals that the planner integrates into in Fig. 12. Our temporal model considers the tides’ varia- 123 Autonomous Robots before a tidal high, and the second is when the tides turn at the top. The period of the two missions is marked in red in the ﬁgure. Hourly predicted ocean currents between 09:00 and 13:00 from the SINMOD model can be seen in Fig. 13. As expected, the currents are highly inﬂuenced by the tides showing an inward ﬂow (eastward) before the tidal high (09:00–11:00) and weaker ocean currents during the tidal high (11:00–13:00). These predictions are used as input to our onboard model as the forcing in the SPDE model. Fig. 11 Showing the AUV transect from waypoint to waypoint, where the distance between the waypoints was 218 ms. The color of the transect 8.3.1 Mission 1: 08:53–10:44 indicates normalized mine tailings concentration values, and the two red surfaces show the depth layer (22–25 ms) where samples are collected (Color ﬁgure online) The ﬁrst mission was conducted between 08:53 and 10:44 on December 17th, before a tidal high. The estimated cur- rent from the SINMOD model can be seen in Figs. 13a, 13b, 13c. Figure 14 presents the onboard model results. The left column shows the updated mean of the onboard model, rep- resenting the mine tailings log(concentration) for selected time steps. The right column shows the onboard model stan- dard deviation at selected time steps and the AUV path up Fig. 12 The blue line shows the tidal cycle on December 17th, showing until the given time. The red line represents the AUV path, low tides at around 07:00 and high tides around 13:00. The red boxes and the red dots are the waypoints selected by our sampling indicate the time intervals for the two missions strategy. The results show increased tailing concentration in the eastern area, especially at the beginning of the mis- tions, and we choose to focus on two missions occurring sion (see Fig. 14a,14c and 14e). We observe a decay in the at different intervals in the tidal cycle. The ﬁrst mission is concentration in this area after 1.5 h (see Fig. 14g). This is Fig. 13 Hourly ocean current predictions from SINMOD between 09:00 and 13.00. The currents are highly inﬂuenced by the tides showing an inward ﬂow (eastward) before the tidal high (09:00–11:00) and weaker ocean currents during the tidal high (11:00–13:00) 123 Autonomous Robots Fig. 14 Results from mission 1. The left column shows the path up until the given time. The red line represents the AUV path, and updated mean of the onboard model, representing the mine tailings the red dots are the waypoints selected by our sampling strategy (Color log(concentration) for selected time steps. The right column shows the ﬁgure online) onboard model standard deviation at selected time steps and the AUV when the ocean currents get weaker due to the tidal cycle. In (see Fig. 13d and 13e). Figure 15 presents the onboard model the beginning, the chosen sampling strategy seems to focus results. The left column shows the updated mean of the on exploring the area, covering large parts of the whole onboard model, and the right column shows the standard rectangle. After the exploration, the focus switches to the deviation together with the selected AUV path. This mission eastern area, where high mine tailing concentrations have does not contain areas of clearly increased tailings concen- been observed. tration, and we see that the AUV chooses a different sampling strategy. Since no area shows increased concentration, there is no extra focus on any special area. Instead, the waypoints 8.3.2 Mission 2: 11:52–13:24 are chosen to cover as much as possible of the whole survey area during the whole mission. The second mission was conducted during a tidal high. Figure 15 shows results from this survey. The predicted SIN- MOD ocean currents are much weaker than for mission 1 123 Autonomous Robots Fig. 15 Results from mission 2. The left column shows the path up until the given time. The red line represents the AUV path, and updated mean of the onboard model, representing the mine tailings the red dots are the waypoints selected by our sampling strategy (Color log(concentration) for selected time steps. The right column shows the ﬁgure online) onboard model standard deviation at selected time steps and the AUV 9 Discussion with our assumptions on how the tailings concentration behave. Since the outlet is located on the east side of our The chosen sampling paths (see second column of Fig. 14 and area, we assume elevated concentration in the east. It is also 15) follow two different strategies. In the beginning, they are reasonable to assume that the concentration of mine tailings quite similar, focusing on exploring the area, but after the in the survey area will be higher when the currents ﬂow west- exploration, we see differences in waypoint selection. The ward, as is the case for mission 1. We have also compared sampling path in mission 1 (Fig. 14) focuses on the eastern the results with the ocean model DREAM data and observed area where increased concentration is observed, while the the same tendencies. However, this is not an accurate com- sampling path in mission 2, where there are no areas with parison since the proxy model is trained with data from this increased tailings concentration, continues to cover the whole model. region also after the initial exploration. For both surveys, the standard deviation is reduced at the end of the mission, indicating that the method manages to explore the area. The 10 Closing remarks selection of different strategies shows that the method is adaptive and will change depending on prior data and obser- In this work, we propose a method for adaptive sampling vations and that the sampling is focused on areas with high of ocean processes with an AUV. Prior information from tailings concentration. The AUV data are undoubtedly useful ocean models is combined with information from in situ for reﬁning the DREAM results. measurements to obtain the best possible sampling strategy. It is hard to rate the ﬁeld experiments’ performance since A spatio-temporal GP proxy model with a non-stationary there is no ground truth for this case. Still, the model mean covariance function and a process model based on ﬁnite dif- concentration (ﬁrst column of Fig. 14 and 15 coincide well ference solutions of an advection–diffusion SPDE are built. 123 Autonomous Robots This is a low-complexity model which can run onboard the References AUV. We can choose informative sampling locations by using Berget, G. E., Eidsvik, J., Alver, M., Py, F., Grøtli, E. I., & Johansen, an objective function favoring locations with high uncertainty T. A. (2019). Adaptive underwater robotic sampling of dispersal and high predicted particle concentration. The main contribu- dynamics in the coastal ocean. In Proceedings of the international tion of the paper is to have a dynamic (spatio-temporal) model symposium on robotics research Accepted. onboard a dynamic agent. This enables reliable uncertainty Berget, G. E., Fossum, T., Johansen, T. A., Eidsvik, J., & Rajan, K. (2018). Adaptive sampling of ocean processes using an AUV with quantiﬁcation in spatial predictions while actively search- a gaussian proxy model. IFAC-PapersOnLine, 51, 238–243. ing for ’hot-spots’ in space and time. The method is tested Binney, J., Krause, A., & Sukhatme, G. S. (2013). Optimizing waypoints in simulations and compared with simpler models and sam- for monitoring spatiotemporal phenomena. The International pling strategies. A simulation study shows that the SPDE Journal of Robotics Research, 32(8), 873–888. Choi, H.-L., & How, J. P. (2010). Continuous trajectory planning of process model gives improved prediction results in a dynamic mobile sensors for informative forecasting. Automatica, 46(8), environment, and the sampling strategy using the suggested 1266–1275. objective function performs better than the sampling strategy Cressie, N., & Wikle, C. (2011). Statistics for Spatio-Temporal Data. using a simpler lawnmower design. Since there is no easy CourseSmart Series. Hoboken: Wiley. Das, J., Py, F., Harvey, J. B., Ryan, J. P., Gellene, A., Graham, R., way to obtain a ground truth for ﬁeld experiments, it is difﬁ- Caron, D. A., Rajan, K., & Sukhatme, G. S. (2015). Data-driven cult to determine how well the method performs in the ﬁeld. robotic sampling for marine ecosystem monitoring. The Interna- Still, our results agree well with our assumptions of the mine tional Journal of Robotics Research, 34(12), 1435–1452. tailings distribution at the time, and the selection of differ- Eidsvik, J., Mukerji, T., & Bhattacharjya, D. (2015). Value of infor- mation in the earth sciences: Integrating spatial modeling and ent sampling strategies for two separate missions illustrates decision analysis. Cambridge: Cambridge University Press. the method’s adaptivity. Future work includes expanding the Foss, K. H., Berget, G. E., & Eidsvik, J. (2022). Using an autonomous model to 3D (with depth coordinates), adding current estima- underwater vehicle with onboard stochastic advection-diffusion tion from measurements onboard the AUV, including ocean models to map excursion sets of environmental variables. Envi- ronmetrics, e2702. currents in the waypoint selection method, and exploration Fossum, T., Eidsvik, J., Ellingsen, I., Alver, M., Fragoso, G., Johnsen, using multiple AUVs. It is also interesting to test such meth- G., Mendes, R., Ludvigsen, M., & Rajan, K. (2018). Information- ods for extended time periods, potentially combining the driven robotic sampling in the coastal ocean. Journal of Field AUV measurements with satellite data or data from drones. Robotics, 35, 1101–1121. Gneiting, T. (2002). Nonseparable, stationary covariance functions for Funding Open access funding provided by SINTEF. This work is sup- space-time data. Journal of the American Statistical Association, ported by the norges forskningsråd (Grand Nos. 223254, 267793) 97(458), 590–600. Grifﬁes, S. M., Böning, C., Bryan, F. O., Chassignet, E. P., Gerdes, Data availability The datasets generated during and/or analysed dur- R., Hasumi, H., Hirst, A., Treguier, A.-M., & Webb, D. (2000). ing the current study are available from the corresponding author on Developments in ocean climate modelling. Ocean Modelling, 2, reasonable request. 123–192. Halpern, B. S., Walbridge, S., Selkoe, K. A., Kappel, C. V., Micheli, F., D’Agrosa, C., Bruno, J. F., Casey, K. S., Ebert, C., Fox, H. Declarations E., Fujita, R., Heinemann, D., Lenihan, H. S., Madin, E. M. P., Perry, M. T., Selig, E. R., Spalding, M., Steneck, R., & Watson, R. (2008). A global map of human impact on marine ecosystems. Conﬂict of interest The authors declare that they have no conﬂict of Science, 319(5865), 948–952. interest. Hwang, J., Bose, N., & Fan, S. (2019). AUV adaptive sampling methods: Areview. Applied Sciences, 9(15). Open Access This article is licensed under a Creative Commons Hyun, J. W., Li, Y., Huang, C., Styner, M., Lin, W., & Zhu, H. (2016). Attribution 4.0 International License, which permits use, sharing, adap- STGP: Spatio-temporal gaussian process models for longitudinal tation, distribution and reproduction in any medium or format, as neuroimaging data. NeuroImage, 134, 550–562. long as you give appropriate credit to the original author(s) and the Jennison, C., & Turnbull, B. W. (1999). Group sequential methods with source, provide a link to the Creative Commons licence, and indi- applications to clinical trials. Boca Raton: CRC Press. cate if changes were made. The images or other third party material Jun, M., & Stein, M. L. (2008). Nonstationary covariance models for in this article are included in the article’s Creative Commons licence, global data. The Annals of Applied Statistics, 2(4), 1271–1289. unless indicated otherwise in a credit line to the material. If material Krause, A., Singh, A., & Guestrin, C. (2008). Near-optimal sensor is not included in the article’s Creative Commons licence and your placements in gaussian processes: Theory, efﬁcient algorithms and intended use is not permitted by statutory regulation or exceeds the empirical studies. Journal of Machine Learning and Research, 9, permitted use, you will need to obtain permission directly from the copy- 235–284. right holder. To view a copy of this licence, visit http://creativecomm Kvassnes, A., & Iversen, E. (2013). Waste sites from mines in Norwe- ons.org/licenses/by/4.0/. gian fjords. Mineralproduksjon, 3, A27–A38. Lermusiaux, P., Chiu, C.-S., Gawarkiewicz, G., Abbot, P., Robinson, A., Miller, R., Haley, P., Leslie, W., Majumdar, S., Pang, A., & Lekien, F. (2006). Quantifying uncertainties in ocean predictions. Oceanography, 19, 92–105. 123 Autonomous Robots Luo, W. & Sycara, K. (2018). Adaptive sampling and online learning in Storvik, G., Frigessi, A., & Hirst, D. (2002). Stationary space-time multi-robot sensor coverage with mixture of gaussian processes. In Gaussian ﬁelds and their time autoregressive representation. Sta- 2018 IEEE international conference on robotics and automation tistical Modelling, 2(2), 139–161. (ICRA) (pp. 6359–6364). Trannum, H. C., Nilsson, H. C., Schaanning, M. T., & Øxnevad, S. Ma, K.-C., Liu, L., Heidarsson, H., & Sukhatme, G. (2018). Data-driven (2010). Effects of sedimentation from water-based drill cuttings learning and planning for environmental sampling. Journal of Field and natural sediment on benthic macrofaunal community structure Robotics, 35(5), 643–661. and ecosystem processes. Journal of Experimental Marine Biology Matérn, B. (2013). Spatial variation. Meddelanden från Statens Skogs- and Ecology, 383(2), 111–121. forskningsinstitut, 36(5), 1–144. Wassmann, P., Slagstad, D., Riser, C. W., & Reigstad, M. (2006). Mod- Morello, E., Haywood, M., Brewer, D., Apte, S., Asmund, G., Kwong, elling the ecosystem dynamics of the barents sea including the Y., & Dennis, D. (2016). The ecological impacts of submarine marginal ice zone: II. Carbon ﬂux and interannual variability. Jour- tailings placement. In Hughes, R., Hughes, D., Smith, I. P., & nal of Marine Systems, 59(1): 1–24. Dale, A. (eds.), Oceanography and Marine Biology: An Annual Zhang, Y., Godin, M. A., Bellingham, J. G., & Ryan, J. P. (2012). Using Review, volume 54 (pp. 315–366). CRC Press, 1 edition. an autonomous underwater vehicle to track a coastal upwelling Nepstad, R., Liste, M., Alver, M. O., Nordam, T., Davies, E., & Glette, front. IEEE Journal of Oceanic Engineering, 37(3), 338–347. T. (2020). High-resolution numerical modelling of a marine mine tailings discharge in Western Norway. Regional Studies in Marine Publisher’s Note Springer Nature remains neutral with regard to juris- Science, 39, 101404. dictional claims in published maps and institutional afﬁliations. Pereira, A., Binney, J., Hollinger, G., & Sukhatme, G. (2013). Risk- aware path planning for autonomous underwater vehicles using predictive ocean models. Journal of Field Robotics, 30(5), 741– Gunhild Elisabeth Berget is Pinto, J., Calado, P., Braga, J., Dias, P., Martins, R., Marques, E., & a Ph.D. Student in Engineering Sousa, J. (2012). Implementation of a control architecture for Cybernetics at The Norwegian networked vehicle systems. IFAC Proceedings Volumes, 45(5), University of Science and Tech- 100–105. nology (NTNU). She received her Pinto, J., Dias, P., Martins, R., Fortuna, J., Marques, E., & Sousa, J. M.Sc. degree in Applied Mathe- (2013). The LSTS toolchain for networked vehicle systems. In matics in 2017 from NTNU. Her OCEANS 2013 MTS/IEEE Bergen: The Challenges of the Northern Ph.D. project is related to adap- Dimension. tive sampling using AUVs focus- Rajan, K. & Py, F. (2012). T-REX: Partitioned inference for AUV mis- ing on modeling, data assimila- sion control. Further advances in unmanned marine vehicles (pp. tion and path planning. 171–199). Ramirez-Llodra, E., Trannum, H., Evenset, A., Levin, L., Andersson, M., Finne, T., Hilário, A., Flem, B., Christensen, G., Schaanning, M., & Vanreusel, A. (2015). Submarine and deep-sea mine tailing placements: A review of current practices, environmental issues, natural analogs and knowledge gaps in norway and internationally. Jo Eidsvik is Professor of Statistics Marine Pollution Bulletin, 97(1), 13–35. in the Department of Mathemat- Rasmussen, C. E., & Williams, C. K. I. (2005). Gaussian processes for ical Sciences at the Norwegian machine learning (adaptive computation and machine learning). University of Science and Tech- Cambridge: The MIT Press. nology (NTNU). He got his M.Sc. Richardson, R. A. (2017). Sparsity in nonlinear dynamic spatiotemporal (1997) in Applied Mathematics models using implied advection. Environmetrics, 28(6), e2456. from the University of Oslo and Rye, H., Reed, M., & Ekrol, N. (1998). The PARTRACK model for his Ph.D. (2003) in Statistics from calculation of the spreading and deposition of drilling mud, chem- NTNU. His research interest are icals and drill cuttings. Environmental Modelling Software, 13(5), spatiotemporal statistics and com- 431–441. putational statistics, often applied Rye, H., Reed, M., Frost, T., Smit, M., Durgut, I., Johansen, O., & to the Earth sciences. Along with Ditlevsen, M. (2008). Development of a numerical model for cal- methods for parameter estimation, culating exposure to toxic and nontoxic stressors in the water prediction and data assimilation in column and sediment from drilling discharges. Integrated Envi- these settings, he also works with ronmental Assessment and Management, 4, 194–203. questions related to the design of experiments, sequential uncertainty Särkkä, S. (2013). Bayesian ﬁltering and smoothing. Cambridge Uni- reduction and value of information analysis in this context of spatio- versity Press: Cambridge. temporal applications. He co-authored a book on ‘Value of informa- Sigrist, F., Künsch, H., & Stahel, W. (2015). Stochastic partial differen- tion in the Earth Sciences’, Cambridge University press, 2015. He is tial equation based modeling of large space-time data sets. Journal Associate Editor of Statistics and Computing and Mathematical Geo- of the Royal Statistical Society: Series B (Statistical Methodology), sciences. 77, 3–33. Slagstad, D., & McClimans, T. (2005). Modeling the ecosystem dynam- ics of the barents sea including the marginal ice zone: I. Physical and chemical oceanography. Journal of Marine Systems, 58, 1–18. Stewart, R. H. (2008). Introduction To physical oceanography.Texas A&MUniversity. 123 Autonomous Robots Morten Omholt Alver received the M.Sc. and Ph.D .in Engineer- ing Cybernetics from The Nor- wegian University of Science and Technology (NTNU) in 2002 and 2007, respectively. He is currently an Associate Professor at the Department of Engineering Cyber- netics at NTNU. His research inter- ests include ocean modelling and data assimilation, modelling of bio- logical systems, and instrumen- tation and automation in marine aquaculture. Tor Arne Johansen received the M.Sc. degree in 1989 and the Ph.D. degree in 1994, both in electrical and computer engineer- ing, from the Norwegian Univer- sity of Science and Technology, Trondheim, Norway. From 1995 to 1997, he worked at SINTEF as a researcher before he was appointed Associated Professor at the Norwegian University of Sci- ence and Technology in Trond- heim in 1997 and Professor in 2001. He has published several hundred articles in the areas of control, estimation and optimization with applications in the marine, aerospace, automotive, biomedical and process industries. In 2002 Johansen co-founded the company Marine Cybernetics AS where he was Vice President until 2008. Prof. Johansen received the 2006 Arch T. Colwell Merit Award of the SAE, and is currently a prin- cipal researcher within the Center of Excellence on Autonomous Marine Operations and Systems (NTNU-AMOS) and director of the Unmanned Aerial Vehicle Laboratory at NTNU and the SmallSat Lab- oratory at NTNU. He recently co-founded the spin-off companies Scout Drone Inspection, UBIQ Aerospace, Zeabuz and SentiSystems.

Journal

Autonomous Robots – Springer Journals

Published: Apr 1, 2023

Keywords: Adaptive sampling; Gaussian process; Stochastic modeling; AUV; Mine tailings

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Dynamic stochastic modeling for adaptive sampling of environmental variables using an AUV

Dynamic stochastic modeling for adaptive sampling of environmental variables using an AUV

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Dynamic stochastic modeling for adaptive sampling of environmental variables using an AUV

Dynamic stochastic modeling for adaptive sampling of environmental variables using an AUV

References (44)

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies