Adaptive sampling for learning gaussian processes using mobile sensor networks.

This paper presents a novel class of self-organizing sensing agents that adaptively learn an anisotropic, spatio-temporal gaussian process using noisy measurements and move in order to improve the quality of the estimated covariance function. This approach is based on a class of anisotropic covariance functions of gaussian processes introduced to model a broad range of spatio-temporal physical phenomena. The covariance function is assumed to be unknown a priori. Hence, it is estimated by the maximum a posteriori probability (MAP) estimator. The prediction of the field of interest is then obtained based on the MAP estimate of the covariance function. An optimal sampling strategy is proposed to minimize the information-theoretic cost function of the Fisher Information Matrix. Simulation results demonstrate the effectiveness and the adaptability of the proposed scheme.

Introduction

In recent years, due to global climate changes, more environmental scientists are interested in the change of ecosystems over vast regions in lands, oceans, and lakes. For instance, for certain environmental conditions, rapidly reproducing harmful algal blooms in the Great Lakes can produce cyanotoxins. Besides such natural disasters, there exist growing ubiquitous possibilities of the release of toxic chemicals and contaminants in the air, lakes, and public water systems. This resulted in the rising demands to utilize autonomous robotic systems that can perform a series of tasks such as estimation, prediction, monitoring, tracking and removal of a scalar field of interest undergoing often complex transport phenomena (Common examples are diffusion, convection, and advection).

Significant enhancements have been made in the areas of mobile sensor networks and mobile sensing vehicles such as unmanned ground vehicles, autonomous underwater vehicles, and unmanned aerial vehicles. Emerging technologies have been reported on the coordination of mobile sensing agents [1-6]. Mobile sensing agents form an ad-hoc wireless communication network in which each agent usually operates under a short communication range, with limited memory and computational power. Mobile sensing agents are often spatially distributed in an uncertain surveillance environment.

The mobility of mobile agents can be designed in order to perform the optimal sampling of the field of interest. Recently in [5], Leonard et al. developed mobile sensor networks that optimize ocean sampling performance defined in terms of uncertainty in a model estimate of a sampled field. In [6], distributed learning and cooperative control were developed for multi-agent systems to discover peaks of the unknown field based on the recursive estimation of an unknown field. In general, we design the mobility of sensing agents to find the most informative locations to make observations for a spatio-temporal phenomenon. To find these locations that predict the phenomena best, one needs a model of the spatio-temporal phenomenon itself. In our approach, we focus on Gaussian processes to model fields undergoing transport phenomena. A Gaussian process (or kriging in geostatistics) has been widely used as a nonlinear regression technique to estimate and predict geostatistical data [7-11]. A Gaussian process is a natural generalization of the Gaussian probability distribution. It generalizes the Gaussian distribution with a finite number of random variables to a Gaussian process with an infinite number of random variables in the surveillance region. Gaussian process modeling enables us to predict physical values, such as temperature and plume concentration, at any of spatial points with a predicted uncertainty level efficiently. For instance, near-optimal static sensor placements with a mutual information criterion in Gaussian processes were proposed by [12,13]. Distributed kriged Kalman filter for spatial estimation based on mobile sensor networks are developed by [14]. A distributed adaptive sampling approach was proposed by [15] for sensor networks to find locations that maximize the information contents by assuming that the covariance function is known up to a scaling parameter. Multi-agent systems that are versatile for various tasks by exploiting predictive posterior statistics of Gaussian processes were developed by [16,17].

The motivation of our work is as follows. Even though there have been efforts to utilize Gaussian processes to model and predict the spatio-temporal field of interest, most of recent papers assume that Gaussian processes are isotropic, implying that the covariance function only depends on the distance between locations. Many studies also assume that the corresponding covariance functions are known a priori for simplicity. However, this is not the case in general as pointed out in literature [12,13,18], in which they treat the non-stationary process by fusing a collection of isotropic spatial Gaussian processes associated with a set of local regions. Hence our motivation is to develop theoretically-sound algorithms for mobile sensor networks to learn the anisotropic covariance function of a spatio-temporal Gaussian process. Mobile sensing agents can then predict the Gaussian process based on the estimated covariance function in a nonparametric manner.

The contribution of this paper is to develop covariance function learning algorithms for the sensing agents to perform nonparametric prediction based on a properly adapted Gaussian process for a given spatio-temporal phenomenon. By introducing a generalized covariance function, we expand the class of Gaussian processes to include the anisotropic spatio-temporal phenomena. The maximum a posteriori probability (MAP) estimator is used to find hyperparameters for the associated covariance function. The proposed optimal navigation strategy for autonomous vehicles minimizes the information-theoretic cost function such as D-or A-optimality criterion using the Fisher Information Matrix (or Cramér-Rao lower bound (CRLB)[19], improving the quality of the estimated covariance function. A Gaussian process with a time-varying covariance function has been proposed to demonstrate the adaptability of the proposed scheme.

This paper is organized as follows. In Section 2, we briefly review the mobile sensing network model and the notation related to a graph. A nonparametric approach to predict a field of interest based on measurements is presented in Section 3. Section 4 introduces a covariance function learning algorithm for an anisotropic, spatio-temporal Gaussian process. An optimal navigation strategy is described in Section 5. In Section 6, simulation results illustrate the usefulness of our approach and its adaptability for unknown and/or time-varying covariance functions.

The standard notation will be used in the paper. Let ℝ, ℝ≥0, ℤ denote, respectively, the set of real, non-negative real, and integer numbers. The positive semi-definiteness of a matrix A is denoted by A ≽ 0. Let |B| denotes the determinant of a matrix B. 𝔼 denote the expectation operator.

Mobile Sensor Networks

First, we explain the mobile sensing network and the measurement model used in this paper. Let N be the number of sensing agents distributed over the surveillance region 𝒬 ∈ ℝ2. Assume that 𝒬 is a compact set. The identity of each agent is indexed by ℐ := {1, 2,⋯, N}. Let q(t) ∈ 𝒬 be the location of the i-th sensing agent at time t ∈ ℝ≥0. We assume that the measurement y(q(t), t) of agent i is the sum of the scalar value of the Gaussian process z(q(t), t) and sensor noise w(t), at its position q(t) and some measurement time t,

The communication network of mobile agents can be represented by a graph with edges. Let G(t) = (ℐ, ℰ(t)) be an undirected communication graph such that an edge (i, j) ∈ ℰ(t) if and only if agent i can communicate with agent j ≠ i. We define the neighborhood of agent i at time t by N(t) := {j | (i, j) ∈ ℰ(t), i ∈ ℐ}. We also define the closed neighborhood of agent i at time t by the union of its index and its neighbors, i.e., N̄ (t) := {i} ∪ N(t).

The Nonparametric Approach

With the spatially distributed sampling capability, agents need to estimate and predict the field of interest by fusing the collective samples from different space and time indices. We show a nonparametric approach to predict a field of interest based on measurements. We assume that a field undergoing a physical transport phenomenon can be modeled by a spatio-temporal Gaussian process, which can be used for nonparametric prediction.

Consider a spatio-temporal Gaussian process: where s, s′ ∈ 𝒬, t, t′ ∈ ℝ≥0 and μ(s, t) denotes the mean value at location s and time t. We then propose the following generalized covariance function 𝒦(s, t, s′, t′;Ψ) with a hyperparameter vector : where s is the l-th entry of s. {σ, σ} and σ are kernel bandwidths for space and time, respectively. Equation (2) shows that points close in the measurement space and time indices are strongly correlated and produce similar values. In reality, the larger temporal distance two measurements are taken with, the less correlated they become, which strongly supports our generalized covariance function in Equation (2). This may also justify the truncation (or windowing) of the observed time series data to limit the size of the covariance matrix for reducing the computational cost.

In the case that the global coordinates are different from the local model coordinates, a similarity transformation can be used to address this issue. For instance, a rotational relationship between the model basis {e⃗, e⃗} and the global basis {E⃗, E⃗} is: where θ represents the angle of rotation. We can then use the following relationship to change the coordinates: where x and y indicate coordinates in the local basis and X and Y indicate their counterparts in the global basis. Equation (2) can then be rewritten in terms of global coordinates as where s and s are the coordinates in the global basis. In this case, the parameter vector Ψ is redefined as .

Up to time t, agent i has noisy collective data {y(q(t), t) | m ∈ ℤ, j ∈ N̄(t), 1 ≤ m ≤ k}, where N̄(t) denotes the closed neighborhood of agent i at time t. The measurements y(q(t), t) = z(q(t), t) + w(t) are taken at different positions q(t) ∈ 𝒬 and different times t ∈ ℝ≥0. The measurements are corrupted by the sensor and communication noises represented by Gaussian white noise . For the case in which the noise level σ is not known and needs to be estimated, the hyperparameter vector can be expanded to include σ, i.e., . The column-vectorized measurements collected by agent i is denoted by with a joint distribution where n is the total number of observations up to time t, μ := 𝔼(Y≤) is the mean vector of Y≤, Σ:= 𝔼((Y≤ − μ) (Y≤ − μ)) is the covariance matrix of Y≤ obtained by in which δ denotes the Kronecker delta function.

If the covariance function is known a priori, the prediction of the random field z(s, t) at location s and time t is then obtained by where ẑ(s, t|t) := 𝔼 (z(s, t|t)) is and the prediction error variance is where Σ is the covariance of z, obtained by 𝒦(s, t, s, t;Ψ), is the covariance matrix between z and Y, obtained by [Σ] = 𝒦(s, t, s, t; Ψ). Each agent can then predict the field of interest at any location and time with the associated uncertainty in a nonparametric way. In the next section, we present a learning approach for unknown covariance functions.

The MAP Estimate of the Hyperparameter Vector

Without loss of generality, we use a zero mean Gaussian process z(s, t) ∼ 𝒢𝒫(0, 𝒦(s, t, s′, t′)) for modeling the field undergoing a physical transport phenomenon. This is not a strong limitation since the mean of the posterior process is not confined to zero [11].

If the covariance function of a Gaussian process is not known a priori, mobile agents need to estimate parameters of the covariance function (Ψ) based on the observed samples. Using Bayes’ rule, the posterior p(Ψ|Y≤) is proportional to the likelihood p(Y≤|Ψ) times the prior p(Ψ), i.e., At time t, the maximum a posteriori (MAP) estimate Ψ̂ of the hyperparameter vector can be obtained by This is equivalent to maximize the logarithm of the posterior p(Ψ|Y≤), i.e., The log likelihood function is given by where n is the size of Y≤. Notice that if no prior information is given, the MAP estimate in Equation (4) is equal to the maximum likelihood (ML) estimate.

A gradient ascent algorithm is used to find a MAP estimate of Ψ: where is a small positive number which can be obtained by using a backtracking line search and ∇(x) is the partial derivative of f(x) with respect to x. The partial derivative of the log likelihood function with respect to a hyperparameter ψ is given by where . Alternatively, a simplex search method [20] can be used to find a MAP estimate of Ψ. This is a direct search method that does not use numerical or analytic gradients.

After finding a MAP estimate of Ψ, agents can proceed the prediction of the field of interest using Equation (3).

An Adaptive Sampling Strategy

Agents should find new sampling positions to improve the quality of the estimated covariance function in the next iteration at time t. For instance, to precisely estimate the anisotropic phenomenon, i.e., processes with different correlations along x-axis and y-axis directions, sensing agents need to explore and sample measurements along different directions.

To this end, we consider a centralized scheme. Suppose that a leader agent (or a central station) knows the communication graph at the next iteration time t and also has access to all measurements collected by agents. Let Y and Y be the measurements at time t and the collective measurements up to time t, respectively, i.e.,

To derive the optimal navigation strategy, we compute the log likelihood function of observations of Y≤: where n≤ is the size of Y≤.

Since the locations of observations in Y≤ were already fixed, we represent the log likelihood function in terms of a vector of future sampling points q̃ at time t only and the hyperparameter vector Ψ:

Now consider the Fisher Information Matrix (FIM) that measures the information produced by measurements Y≤ for estimating the hyperparameter vector at time t. The Cramér-Rao lower bound (CRLB) theorem states that the inverse of the FIM is a lower bound of the estimation error covariance matrix [19,21]: where Ψ̂ represents the estimation of Ψ at time t. The FIM [19] is given by where the expectation is taken with respect to p(Y≤ | Ψ). The analytical closed-form of FIM is given by Since the true value of Ψ is not available, we will evaluate the FIM in Equation (7) at the currently available best estimate Ψ̂. This has been an effective practical solution when we evaluate the FIM and estimate Ψ simultaneously [22,23]. The term due to the MAP estimation error in evaluating the FIM in Equation (7) will decrease as the number of samples increases.

We can expect that minimizing the CRLB results in a decrease of uncertainty in estimating Ψ [22]. Using the D-optimality criterion [24,25], the objective function J is given by Minimizing J in Equation (8) corresponds to minimizing the volume of the ellipsoid which represents the maximum confidence region for the estimated hyperparameters. However, if one hyperparameter has a much larger variance compared to the others, minimizing the volume may not be very useful [25]. As an alternative, the A-optimality which minimizes the sum of the variances may be used. The objective function J based on the A-optimality criterion is A control law for the mobile sensor network can be formulated as follows: A gradient descent strategy can be used to find the next optimal sampling positions: where α( is a small positive number which can be obtained by using a backtracking line search. Alternatively, a control law for the mobile sensor network can be formulated as follows: where where n = 2 denotes the dimension of the surveillance region 𝒬 and δ is the maximum distance for each agent to move in x and y directions.

However, optimization on ln p(Y≤|Ψ) in Equation (6) and J(q̃, Ψ̂) in Equation (8) or (9) can be numerically costly due to the increasing size of Σ. One way to deal with this problem is to use a truncated date set instead of using Y≤. In addition, this approach based on the truncated observations can be viewed as a strategy to deal with a slowly time-varying parameter vector Ψ, which will further investigated in Section 6.2.

The overall protocol for the sensor network is summarized as in Table 1.

Table 1.

An adaptive sampling strategy for mobile sensor networks.

Learning: At time tk, the sensor network updates Ψ̂k using a MAP estimate Equation (4) for a data set Y≤k. Start this MAP optimization with the initial point Ψ̂k−1. Prediction: For given Y≤k and Ψ̂k, agents can compute prediction at any point and time using Equation (3), i.e., p(z(s, t)|Y≤k; Ψ̂k). Sampling: Based on {Ψ̂k, Y≤k}, the sensor network computes the control Equation (10) in order to maximize J(q̃, Ψ̂k). Update the positions of agents accordingly and collect measurements at time tk+1. Repeat the steps 1–3 until Ψ converges.

Simulation Results

In this section, we evaluate the proposed approach for a spatio-temporal Gaussian process (Section 6.1) and an advection-diffusion process (Section 6.3). For both cases, we compare the simulation results using the proposed optimal sampling strategy with results using a benchmark random sampling strategy. In this random sampling strategy, each agent was initially randomly deployed in the surveillance region. At each time step, the next sampling position for agent i is generated randomly with the same mobility constraint, viz. a random position within a square region with length 2 centered at the current position q. For fair comparison, the same values were used for all other conditions. In Section 6.2, our approach based on truncated observations has been applied to a Gaussian process with a time-varying covariance function to demonstrate the adaptability of the proposed scheme.

A Spatio-Temporal Gaussian Process

We apply our approach to a spatio-temporal Gaussian process. The Gaussian process was numerically generated for the simulation [11]. The hyperparameters used in the simulation were chosen such that . Two snap shots of the realized Gaussian random field at time t = 1 and t = 20 are shown in Figure 1. In this case, N = 5 mobile sensing agents were initialized at random positions in a surveillance region . The initial values for the algorithm were given to be . A prior of the hyperparameter vector has been selected as where p(σ) = p(σ) = p(σ) = p(σ) = Γ(5, 2), and p(σ) = Γ(5, 0.2). Γ(a, b) is a Gamma distribution with mean ab and variance ab2 in which all possible values are positive. The gradient method was used to find the MAP estimate of the hyperparameter vector.

Figure 1.

Snap shots of the realized Gaussian process at (a) t = 1 and (b) t = 20.

For simplicity, we assumed that the global basis is the same as the model basis. We considered a situation where at each time, measurements of agents are transmitted to a leader (or a central station) that uses our Gaussian learning algorithm and sends optimal control back to individual agents for next iteration to improve the quality of the estimated covariance function. The maximum distance for agents to move in one time step was chosen to be 1 for both x and y directions. The A-optimality criterion was used for optimal sampling.

For both proposed and random strategies, Monte Carlo simulations were run for 100 times and the statistical results are shown in Figure 2. The estimates of the hyperparameters (shown in circles and error bars) tend to converge to the true values (shown in dotted lines) for both strategies. As can be seen, the proposed scheme (Figure 2(a)) outperforms the random strategy (Figure 2(b)) in terms of the A-optimality criterion.

Figure 2.

Monte Carlo simulation results (100 runs) for a spatio-temporal Gaussian process using (a) the random sampling strategy, and (b) the adaptive sampling strategy. The estimated hyperparameters are shown in blue circles with error-bars. The true hyperparameters that used for generating the process are shown in red dashed lines.

Figure 3 shows the predicted field along with agents’ trajectories at time t = 1 and t = 20 for one trial. As shown in Figure 1(a) and Figure 3(a), at time t = 1, the predicted field is far from the true field due to the inaccurate hyperparameters estimation and small number of observations. As time increases, the predicted field will be closer to the true field due to the improved quality of the estimated the covariance function and the cumulative observations. As expected, at time t = 20, the quality of the predicted field is very well near the sampled positions as shown in Figure 3(b). With 100 observations, the running time is around 30s using Matlab, R2008a (MathWorks) in a PC (2.4 GHz Dual-Core Processor). No attempt has been made to optimize the code. After converging to a good estimate of Ψ, agents can switch to a decentralized configuration and collect samples for other goals such as peak tracking and prediction of the process [6,16,17].

Figure 3.

The predicted fields along with agents’ trajectories at (a) t = 1 and (b) t = 20.

Time-Varying Covariance Functions

To illustrate the adaptability of the proposed strategy to time-varying covariance functions, we introduce a Gaussian process defined by the following covariance function. The time-varying covariance function is modeled by a time-varying weighted sum of two known covariance functions 𝒦1(·) and 𝒦2(·) such as where λ(t) ∈ [0, 1] is a time-varying weight factor that needs to be estimated. In the simulation study, 𝒦1(·) is constructed with σ = 1, σ = 0.2, σ = 0.1, σ = 8, and σ = 0.1; and 𝒦2(·) is with σ = 1, σ = 0.1, σ = 0.2, σ = 8, and σ = 0.1. This Gaussian process defined in (12) with theses particular 𝒦1 and 𝒦2 effectively models hyperparameter changes in x and y directions.

To improve the adaptability, the mobile sensor network uses only observations sampled during the last 20 iterations for estimating λ(t) online. The true λ(t) and the estimated λ(t) are shown in Figure 4(a,b), respectively. From Figure 4, it is clear that the weighting factor λ(t) can be estimated accurately after some delay about 5–8 iterations. The delay is due to using the truncated observations that contain past observations since the time-varying covariance function changes continuously in time.

Figure 4.

(a) The weighting factor λ(t) and (b) the estimated λ(t).

Fitting a Gaussian Process to an Advection-Diffusion Process

We apply our approach to a spatio-temporal process generated by physical phenomena (advection and diffusion). This work can be viewed as a statistical modeling of a physical process, i.e., as an effort to fit a Gaussian process to a physical advection-diffusion process in practice. The advection-diffusion model developed in [26] was used to generate the experimental data numerically. An instantaneous release of Qkg of gas occurs at a location (x0, y0, z0). This is then spread by the wind with mean velocity Assuming that all measurements are recorded at a level z = 0, and the release occurs at a ground level (i.e., z0 = 0), the concentration C at an arbitrary location (x, y, 0) and time t is described by the following analytical solution [27]: where Δx = x − x0, Δy = y − y0, and Δt = 5(t − 1) + t0. The parameters used in the simulation study are shown in Table 2. Notice that this process generates an anisotropic concentration field with parameters K = 20m2/min and K = 10m2/min as in Table 2. The fields at time t = 1 and t = 20 are shown in Figure 5. Notice the center of the concentration moved. In this case, N = 5 mobile sensing agents were initialized at random positions in a surveillance region .

Table 2.

Parameters used in simulation.

Parameter	Notation	Unit	Value
Number of agents	Ns	-	5
Sampling time	ts	min	5
Initial time	t0	min	100
Gas release mass	Q	kg	106
Wind velocity in x axis	ux	m/min	0.5
Eddy diffusivity in x axis	Kx	m2/min	20
Eddy diffusivity in y axis	Ky	m2/min	10
Eddy diffusivity in z axis	Kz	m2/min	0.2
Location of explosion	x0	m	2
Location of explosion	y0	m	5
Location of explosion	z0	m	0
Sensor noise level	σw	kg/m3	0.1

Figure 5.

Snap shots of the advection-diffusion process at (a) t = 1 and (b) t = 20.

The initial values for the algorithm was chosen to be where we assumed σ = 1 and σ = 0.1. For this application, we did not assume any prior knowledge about the covariance function. Hence, the MAP estimator was the same as the ML estimator. The gradient method was used to find the ML estimate.

We again assumed that the global basis is the same as the model basis and assumed all agents have the same level of measurement noises for simplicity. In our simulation study, agents start sampling at t0 = 100min and take measurements at time t with a sampling time of t = 5min as in Table 2.

Monte Carlo simulations were run for 100 times, and Figure 6 shows the estimated σ, σ, and σ with (a) the random sampling strategy and (b) the optimal sampling strategy, respectively. With 100 observations, the running time at each time step is around 20s using Matlab, R2008a (MathWorks) in a PC (2.4 GHz Dual-Core Processor). No attempt has been made to optimize the code. As can be seen in Figure 6, the estimates of the hyperparameters tend to converge to similar values for both strategies. Clearly, the proposed strategy outperforms the random sampling strategy in terms of the estimation error variance.

Figure 6.

Simulation results (100 runs) for a advection-diffusion process. The estimated hyperparameters with (a) random sampling and (b) optimal sampling.

Summary

In this paper, we presented a novel class of self-organizing sensing agents that learn an anisotropic, spatio-temporal Gaussian process using noisy measurements and move in order to improve the quality of the estimated covariance function. The MAP estimator was used to estimate the hyperparameters in the unknown covariance function and the prediction of the field of interest was obtained based on the MAP estimates. An optimal navigation strategy was proposed to minimize the information-theoretic cost function of the Fisher Information Matrix for the estimated hyperparameters. The proposed scheme was applied to both a spatio-temporal Gaussian process and a true advection-diffusion field. Simulation study indicated the effectiveness of the proposed scheme and the adaptability to time-varying covariance functions. The trade-off between a precise estimation and computational efficiency using truncated observations will be studied in the future work.