A Low-Complexity Approach for Improving the Accuracy of Sensor Networks

Abstract

The paper addresses the problem of improving the accuracy of the measurements collected by a sensor network, where simplicity and cost-effectiveness are of utmost importance. An adaptive Bayesian approach is proposed to this aim, which allows improving the accuracy of the delivered estimates with no significant increase in computational complexity. Remarkably, the resulting cooperative algorithm does not require prior knowledge of the (hyper)parameters and is able to provide a “denoised” version of the monitored field without losing accuracy in detecting extreme (less frequent) values, which can be very important for a number of applications. A novel performance metric is also introduced to suitably quantify the capability to both reduce the measurement error and retain highly-informative characteristics at the same time. The performance assessment shows that the proposed approach is superior to a low-complexity competitor that implements a conventional filtering approach.

1. Introduction and Motivations

In the last years, sensor networks have started to be deployed for an increasing number of different applications [1]. The availability of low-cost commercial off-the-shelf nodes, fostered by significant advances in wireless communication technologies and size scaling of integrated circuits, has enabled the deployment of small low-cost sensor nodes with increased lifetime [2]. Typical applications are sensing/estimation of some parameters [3, 4] such as temperature, pollution level [5], electromagnetic exposure [6, 7], or field reconstruction [8, 9]. Such problems are particularly important in environmental monitoring [10], ecology [11], meteorology, agriculture, and related fields as reported in a number of case studies [12, 13]; see also [14] and references therein. More in general, sensing capabilities are currently regarded as a key enabler for smart applications in contexts as diverse as transportation systems [15–17], cyber-physical systems [18–20], and ad hoc networks [21] and in (opportunistic) applications like position estimation for location awareness [22–25]. Finally, interconnection of standalone systems can lead to advanced sensing capabilities, for example, in radar applications [26].

Both centralized and distributed approaches can be adopted to process the information collected through a sensor network [27]. In the centralized approach, data are sent to a fusion center (FC) performing the whole computation [28], while in the distributed one neighboring nodes cooperate in a peer-to-peer fashion until convergence [29, 30]. Regardless of implementation aspects, the goal can be formalized as an inference problem, based on sensor observations, about an underlying physical phenomenon. Clearly, observations are affected by errors introduced in the sensing/measurement process, for instance, due to thermal noise, atmospheric effects, and residual sensor calibration errors. This is especially true when data are collected through general-purpose devices, even smartphones, [31] instead of dedicated (expensive) sophisticated sensors. Therefore, techniques aimed at improving the accuracy of sensor measurements are highly desirable [32, 33].

Unfortunately, a peculiarity of sensor networks is that each sensor has quite limited power supply and computation capabilities; since advanced processing techniques cannot be often implemented with reduced effort, novel low-complexity approaches are needed to actually make sensor network applications feasible in most real-world contexts. As a matter of fact, systems designed to be effectively deployable have to face a number of practical difficulties; thus things are kept as simple and cheap as possible; this, however, may negatively have an impact on the final accuracy [12, 34]. In particular, a coarse granularity may be reported on reconstructed sensing maps, with resulting sharp edges between contiguous levels (as, e.g., in [16]). In other cases, values are averaged [5] or, if their spatial distribution is of concern, smoothed via suitable low-pass filters [35]. A simple moving average is often used in practice, which can be easily implemented as a correlation with a weight mask through a sliding window approach.

In this work, we propose a low-complexity Bayesian approach for improving the accuracy of the measurements, so that a more reliable value of the monitored field is obtained without the need for a sophisticated processing. We will show that, with other things being equal, few lines of code can be effective to improve the final accuracy if the information coming from other sensors is exploited in a cooperative way. The advantage of this approach is that the “filtering” procedure takes into account the statistical relationships in the data at hand. We consider a quite general observation model, where the measurement error is modeled through a Gaussian law. The latter is a versatile model for measurement errors and other random effects [36], supported by the central limit theorem (CLT). A Bayesian approach is used to estimate the value of the field by minimizing the mean squared error at each monitoring point. Different from conventional Bayesian techniques, which require full prior knowledge of the data distribution, we follow an Empirical Bayes approach [37], where the parameters of the prior (hyperparameters) are unknown. We derive their Maximum Likelihood (ML) estimator and show that it has a simple closed form amenable to practical implementation in low-cost devices. An application to spatial field monitoring and reconstruction is reported to highlight the performance improvement compared to a conventional “denoising” technique based on low-pass two-dimensional filtering (moving average).

The rest of the paper is organized as follows. In Section 2, we formulate the problem within the reference scenario. Then, in Section 3, we introduce the proposed filtering approach, which includes the derivation of the Minimum Mean Square Error (MMSE) estimator of the field and the Maximum Likelihood (ML) estimator of the hyperparemters. Besides the mathematical derivation, we provide also a scheme (and the corresponding algorithm's pseudocode) for practical application. In Section 4, we evaluate the proposed approach, showing that it can improve the estimation accuracy without significantly increasing the complexity. To better spotlight the ability to represent correctly specific characteristics of the monitored field, a novel metric is preliminarily introduced. Finally, Section 5 contains the conclusions of the work.

2. Problem Formulation

A general scenario is considered, where N sensor nodes observe a given phenomenon. We denote by $x_i$ , $i=\mathrm{1},\dots ,N$ , the measurement at sensor (equivalently, location) i relative to an unknown local parameter $m_i$ ; that is,

$$\begin{array}{cc}& x_i=m_i+{ϵ}_i,\end{array}$$

(1)

where ${ϵ}_i$ is a zero-mean error term with variance ${\sigma }^{\mathrm{2}}$ . Lacking specific models, it is customary to consider ${ϵ}_i$ as normally distributed result of the CLT. The reference scenario is depicted in Figure 1. The monitored field is represented with shades proportional to the intensity of the parameter under investigation. Circles indicate the locations of the deployed sensors, with (some of the) links reported as edges of a (time-varying) graph, which depends on the communication policies of the deployed network and on environmental conditions. Measurements may be sent to a fusion center if the processing is centralized; in a distributed setup, conversely, nodes will exchange status updates by combining their measurement with the status updates of neighboring nodes until convergence.

Figure 1

Reference scenario of a monitoring sensor network.

One can exploit the fact that sensors are deployed over a “continuous” field for monitoring purposes; hence, they measure a variable which represents a “sampling” of the underlying whole process. As a consequence, some correlation can be expected according to the spatial proximity (distance) between sensors. However, correlations are difficult to model exactly, since they require deep knowledge of the process at hand. Moreover, they are very site-specific and change with time. Complicated models, if available, require in turn computational-intensive techniques; conversely, as mentioned, simple approaches are needed in low-cost sensor networks.

To this aim, we propose modelling the relationships between the sensed points of the field by means of a prior distribution on the measurements with unknown hyperparameters. In particular, given the Gaussian model for $x_i$ , we consider a conjugate prior distribution for the mean $m_i$ , which is again a Gaussian distribution. Thus, the model can be written as

$$\begin{array}{cc}& x_i\hspace{0.17em}|\hspace{0.17em}{m}_i~\mathcal{N}\left(m_i,{\sigma }^{\mathrm{2}}\right),\\ m_i~\mathcal{N}\left(\mu ,{\nu }^{\mathrm{2}}\right).\end{array}$$

(2)

The goal is to recover each $m_i$ , based on the statistic $x_i$ , through the lens of the Bayesian hierarchy, that is, to enhance the estimation quality by considering the probability of observing specific field values and computing the estimator that minimizes the mean squared error (MMSE). However, the fact that the hyperparameters $(\mu ,{\nu }^{\mathrm{2}})$ are unknown makes this approach different from the classical Bayesian framework, where the prior distribution is completely known. In particular, in Section 3.2, we derive in closed form the Maximum Likelihood (ML) estimators for $(\mu ,{\nu }^{\mathrm{2}})$ , following the so-called Empirical Bayes (EB) rationale [37]. To this aim, we consider the whole statistic $\mathbf{x}=[x_{\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\cdots \hspace{0.17em}\hspace{0.17em}{x}_N{]}^T$ and exploit the fact that observations share the same probability distribution. This allows fitting the Bayesian model by the most suitable (hyper)parameters, which are “learned” adaptively from the data, thus allowing for generality. Beforehand, we derive below the MMSE estimator for the value of the field at the sensor locations.

3. Empirical Bayes-Based Measurement Filtering

3.1. Minimum Mean Square Error Field Estimation

The conditional distribution of $x_i$ given $m_i$ is

$$\begin{array}{cc}& p\left(x_i\hspace{0.17em}|\hspace{0.17em}{m}_i\right)=\frac{\mathrm{1}}{\sqrt{\mathrm{2}\pi {\sigma }^{\mathrm{2}}}}\text{exp}\left\{-\frac{{\left(x_i-m_i\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}}\right\}\end{array}$$

(3)

and, as mentioned, we consider a normal prior for $m_i$ so that

$$\begin{array}{cc}& p\left(m_i\right)=\frac{\mathrm{1}}{\sqrt{\mathrm{2}\pi {\nu }^{\mathrm{2}}}}\text{exp}\left\{-\frac{{\left(m_i-\mu \right)}^{\mathrm{2}}}{\mathrm{2}{\nu }^{\mathrm{2}}}\right\}.\end{array}$$

(4)

We can derive the joint probability distribution of $x_i$ and $m_i$ :

$$\begin{array}{l}p(x_i,m_i)=p(x_i|m_i)=p(x_i|m_i)p(m_i)=\frac{1}{2\pi \sigma v}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\cdot exp\{-\frac{1}{2}[\left(\frac{1}{v^2}-\frac{1}{{\mu }^2}\right)m_i^2-2\left(\frac{x_i}{{\sigma }^2}+\frac{{\mu }^2}{v^2}\right)m_i+\frac{{\mu }^2}{v^2}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{x_i^2}{{\sigma }^2}]\}.\hfill \end{array}$$

(5)

The Minimum Mean Square Error (MMSE) estimator of $m_i$ is obtained by computing the conditional mean [37]; therefore, it is first necessary to calculate the posterior distribution of $x_i$ , according to Bayes’ theorem:

$$\begin{array}{cc}& p\left(m_i\hspace{0.17em}|\hspace{0.17em}{x}_i\right)=\frac{p\left(x_i,m_i\right)}{p\left(x_i\right)}.\end{array}$$

(6)

The unconditional distribution of $x_i$ , $p(x_i)$ , can be obtained from the joint distribution $p(x_i,m_i)$ by integrating $m_i$ out:

$$\begin{array}{l}p\left(x_i\right)={{\displaystyle \int }}_{-\infty }^{+\infty }p\left(x_i,m_i\right)\mathrm{d}{m}_i\\ =\frac{\mathrm{1}}{\mathrm{2}\pi \sigma \nu }\text{exp}\left\{-\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{\mu }^{\mathrm{2}}}{{\nu }^{\mathrm{2}}}+\frac{x_i^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}}\right)\right\}{I}_{x_i},\end{array}$$

(7)

where

$$\begin{array}{l}{I}_{x_i}\\ ={{\displaystyle \int }}_{-\infty }^{+\infty }\text{exp}\left\{-\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathrm{1}}{{\nu }^{\mathrm{2}}}+\frac{\mathrm{1}}{{\sigma }^{\mathrm{2}}}\right)m_i^{\mathrm{2}}+\left(\frac{\mu }{{\nu }^{\mathrm{2}}}+\frac{x_i}{{\sigma }^{\mathrm{2}}}\right)m_i\right\}\mathrm{d}{m}_i.\end{array}$$

(8)

The integral above can be computed by completing the square to a Gaussian kernel. After some calculation, the result is

$$\begin{array}{cc}& I_{x_i}=\text{exp}\left\{\frac{{\left(\mu {\sigma }^{\mathrm{2}}+x_i{\nu }^{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}{\nu }^{\mathrm{2}}\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)}\right\}\frac{\sqrt{\mathrm{2}\pi }}{\sqrt{\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)/{\nu }^{\mathrm{2}}{\sigma }^{\mathrm{2}}}}.\end{array}$$

(9)

Substituting (9) into (7) yields

$$\begin{array}{l}p\left(x_i\right)=\frac{\mathrm{1}}{\sqrt{\mathrm{2}\pi \left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)}}\\ ·\text{exp}\left\{-\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{\mu }^{\mathrm{2}}}{{\nu }^{\mathrm{2}}}+\frac{x_i^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}}\right)+\frac{{\left(\mu {\sigma }^{\mathrm{2}}+x_i{\nu }^{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}{\nu }^{\mathrm{2}}\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)}\right\}.\end{array}$$

(10)

Using (6) (together with (5) and (10)) we obtain the final expression of the posterior distribution:

$$\begin{array}{l}p(m_i|x_i)=\frac{1}{\sqrt{2\pi ({\sigma }^2{v}^2/({\sigma }^2+v^2))}}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\cdot exp\{-\frac{{(\mu {\sigma }^2+x_i{v}^2)}^2}{2{\sigma }^2{v}^2({\sigma }^2+v^2)}-\frac{1}{2}\left(\frac{1}{v^2}+\frac{1}{{\sigma }^2}\right)m_i^2\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\left(\frac{\mu }{v^2}+\frac{x_i}{{\sigma }^2}\right)m_i\}=\frac{1}{\sqrt{2\pi ({\sigma }^2{v}^2/({\sigma }^2+v^2))}}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\cdot exp\left\{-\frac{{(m_i-(\mu {\sigma }^2+x_i{v}^2)/({\sigma }^2/v^2))}^2}{2({\sigma }^2{v}^2/({\sigma }^2+v^2))}\right\}\hfill \end{array}$$

(11)

which can be recognized as a Gaussian distribution with mean $\left({\sigma }^{\mathrm{2}}\mu +{\nu }^{\mathrm{2}}{x}_i\right)/\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)$ and variance ${\sigma }^{\mathrm{2}}{\nu }^{\mathrm{2}}/\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)$ . Such a result is important because it yields a very tractable posterior distribution.

Given the result above, the MMSE is obtained in a simple closed form as

$$\begin{array}{cc}& {\widehat{m}}_i=\mathrm{E}\left[m_i\hspace{0.17em}|\hspace{0.17em}{x}_i\right]=\frac{{\sigma }^{\mathrm{2}}\mu +{\nu }^{\mathrm{2}}{x}_i}{{\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}}=\gamma \mu +\left(\mathrm{1}-\gamma \right)x_i,\end{array}$$

(12)

where

$$\begin{array}{cc}& \gamma =\frac{{\sigma }^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}}\in \left(\mathrm{0,1}\right).\end{array}$$

(13)

Equation (12) has a revealing interpretation as convex combination of the measurement $x_i$ and the prior mean μ, according to the parameter γ. The latter depends on the relative importance of the measurement uncertainty (quantified by the variance ${\sigma }^{\mathrm{2}}$ ) with respect to the overall uncertainty (i.e., including also the prior variance ${\nu }^{\mathrm{2}}$ ). It is obvious, therefore, that the hyperparameters $(\mu ,{\nu }^{\mathrm{2}})$ cannot be assumed known a priori as in the conventional Bayesian approach; we propose instead fitting them to the observed data by leveraging cooperation between the nodes, that is, using the whole statistic x. In particular, in the next section, we derive the Maximum Likelihood (ML) estimator of the hyperparameters, finally obtaining the so-called Empirical Bayes (MMSE) estimator for $m_i$ by replacing the estimates $(\widehat{\mu },{\widehat{\nu }}^{\mathrm{2}})$ into (12).

3.2. Maximum Likelihood Hyperparameter Estimation

The joint probability distribution of the sensor measurements is obtained from (10) as

$$\begin{array}{l}p\left(\mathbf{x}\hspace{0.17em}|\hspace{0.17em}\mu ,{\nu }^{\mathrm{2}}\right)={\displaystyle \prod }_{i=\mathrm{1}}^{N}p\left(x_i\hspace{0.17em}|\hspace{0.17em}\mu ,{\nu }^{\mathrm{2}}\right)={\left[\mathrm{2}\pi \left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)\right]}^{-N/\mathrm{2}}\\ ·\text{exp}\left\{{\displaystyle \sum }_{i=\mathrm{1}}^N\left[-\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{\mu }^{\mathrm{2}}}{{\nu }^{\mathrm{2}}}+\frac{x_i^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}}\right)+\frac{{\left(\mu {\sigma }^{\mathrm{2}}+x_i{\nu }^{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}{\nu }^{\mathrm{2}}\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)}\right]\right\},\end{array}$$

(14)

where the conditioning stresses the fact that the hyperparameters are the unknown quantities to be estimated.

It is easy to verify that the expression between square brackets above is simplified as $-\left(\mathrm{1}/\mathrm{2}({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}})\right)(x_i-\mu {)}^{\mathrm{2}}$ , so that (14) can be rewritten as

$$\begin{array}{l}p\left(\mathbf{x}\hspace{0.17em}|\hspace{0.17em}\mu ,{\nu }^{\mathrm{2}}\right)\\ ={\left[\mathrm{2}\pi \left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)\right]}^{-N/\mathrm{2}}\text{exp}\left\{-\frac{{\sum }_{i=\mathrm{1}}^N{\left(x_i-\mu \right)}^{\mathrm{2}}}{\mathrm{2}\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)}\right\}\end{array}$$

(15)

to be maximized with respect to $(\mu ,{\nu }^{\mathrm{2}})$ . It is convenient to rewrite this optimization problem after a logarithmic transformation; negating the result and omitting an irrelevant term, we obtain the equivalent problem:

$$\begin{array}{l}\left(\widehat{\mu },{\widehat{\nu }}^{\mathrm{2}}\right)\\ =\underset{\mu \in \mathbb{R},{\nu }^{\mathrm{2}}\ge \mathrm{0}}{\text{argmin}}\left\{N\text{log}\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)+\frac{{\sum }_{i=\mathrm{1}}^N{\left(x_i-\mu \right)}^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}}\right\}.\end{array}$$

(16)

We can write the Karush-Kuhn-Tucker conditions [38] for this constrained problem as follows:

$$\begin{array}{l}\mathrm{2}{\displaystyle \sum }_{i=\mathrm{1}}^N\left(x_i-\mu \right)=\mathrm{0}\\ \frac{N}{{\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}}-\frac{{\sum }_{i=\mathrm{1}}^N{\left(x_i-\mu \right)}^{\mathrm{2}}}{{\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)}^{\mathrm{2}}}-\lambda =\mathrm{0},\\ -{\nu }^{\mathrm{2}}\le \mathrm{0},\lambda {\nu }^{\mathrm{2}}=\mathrm{0},\hspace{0.17em}\lambda \ge \mathrm{0},\end{array}$$

(17)

where the first two equations are the derivative of the Lagrangian function

$$\begin{array}{cc}& \mathcal{L}\left(\mu ,{\nu }^{\mathrm{2}},\lambda \right)=N\text{log}\left({\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}\right)+\frac{{\sum }_{i=\mathrm{1}}^N{\left(x_i-\mu \right)}^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}+{\nu }^{\mathrm{2}}}-\lambda {\nu }^{\mathrm{2}}\end{array}$$

(18)

with respect to μ and ${\nu }^{\mathrm{2}}$ , the third line is the positivity constraint on the variance (primal feasibility), the fourth equation is the complementary slackness condition, and the last inequality is the dual feasibility condition. From the first equation, we obtain

$$\begin{array}{cc}& \widehat{\mu }=\frac{\mathrm{1}}{N}{\displaystyle \sum }_{i=\mathrm{1}}^N{x}_i,\end{array}$$

(19)

that is, the sample mean, denoted by $\bar{x}$ . From the feasibility conditions, it is easy to realize that the Lagrange multiplier λ is null at the stationary point if ${\nu }^{\mathrm{2}}>\mathrm{0}$ ; in that case, from the second equation, we obtain ${\widehat{\nu }}^{\mathrm{2}}=s^{\mathrm{2}}-{\sigma }^{\mathrm{2}}$ , where

$$\begin{array}{cc}& s^{\mathrm{2}}=\frac{\mathrm{1}}{N}{\displaystyle \sum }_{i=\mathrm{1}}^N{\left(x_i-\widehat{\mu }\right)}^{\mathrm{2}}\end{array}$$

(20)

denotes the sample variance. Any $\lambda >\mathrm{0}$ necessarily implies ${\widehat{\nu }}^{\mathrm{2}}=\mathrm{0}$ in force of the complementary slackness. The ML estimator of ${\nu }^{\mathrm{2}}$ is therefore given by the following expression:

$$\begin{array}{cc}& {\widehat{\nu }}^{\mathrm{2}}=\text{max}\left(\mathrm{0},s^{\mathrm{2}}-{\sigma }^{\mathrm{2}}\right).\end{array}$$

(21)

This automatically accounts for the case $s^{\mathrm{2}}-{\sigma }^{\mathrm{2}}<\mathrm{0}$ , which would lead to a (unfeasible) negative value for ${\widehat{\nu }}^{\mathrm{2}}$ .

3.3. Practical Application of the Algorithm

The approach developed above can be applied in a straightforward way to a monitoring network deployed in a given area. A schematic representation is depicted in Figure 2, which details the scheme of Figure 1. Some of the nodes are exploded to indicate the local processing, which can be summarized as a function f of the local measurement and of the “state” $(\bar{x},s^{\mathrm{2}})$ . The latter is the result of the cooperation in the network and as mentioned can be computed by a fusion center by means of a distributed processing. Clearly, we have that

$$\begin{array}{cc}& f\left(x_i;\bar{x},s^{\mathrm{2}}\right)=\gamma \bar{x}+\left(\mathrm{1}-\gamma \right)x_i,\gamma =\text{min}\left(\mathrm{1},\frac{{\sigma }^{\mathrm{2}}}{s^{\mathrm{2}}}\right),\end{array}$$

(22)

where the expression for γ has been obtained by simple manipulations of (13).

Figure 2

Schematic representation of the proposed estimation approach.

It is reasonable to expect that, except for small-area networks, nodes may not be able to communicate in a completely meshed way, but rather they have limited connectivity dictated by their communication range. For each node i, we can define the set of neighbors $N_i(t)$ according to the coverage area allowed by the communication technology available in the network; then, the proposed algorithm can be applied locally, with each node exchanging a limited number of communication packets to neighboring nodes in order to calculate the value of the “state” $(\bar{x},s^{\mathrm{2}})$ for its local area, that is, by interacting with nodes $k\in N_i(t)$ . Clearly, each node may contribute (as neighbor) to the calculation of different nodes, according to the extent of the coverage areas. At limit, one could even have a single area enclosing all N nodes; otherwise, a certain number of “clusters” will arise, where the (possibly distributed) computation is performed independently, though it may share some information with gateway nodes. This means that the approach is consistent and very scalable.

Summing up, despite the lenghty calculation in Sections 3.1 and 3.2 (necessary for a rigorous derivation), the result is very handy and can be easily implemented in a few lines of code for a generic cluster as follows.

Pseudocode of the Proposed Algorithm

•  

compute $\bar{x}=\left(\mathrm{1}/N\right){\sum }_{i=\mathrm{1}}^N{x}_i$ and $s^{\mathrm{2}}=\left(\mathrm{1}/N\right){\sum }_{i=\mathrm{1}}^N(x_i-\bar{x}{)}^{\mathrm{2}}$

 

compute $\gamma =\mathrm{m}\mathrm{i}\mathrm{n}(\mathrm{1},{\sigma }^{\mathrm{2}}/s^{\mathrm{2}})$

 

compute ${\widehat{m}}_i=\gamma \bar{x}+(\mathrm{1}-\gamma )x_i$ at each node i

4. Performance Assessment

In this section, we show how the proposed approach can be used to improve the accuracy of a sensor network without significantly increasing the computational cost. We resort to simulations to control the ground truth; that is, we simulate the true value of the field (of some physical quantity, namely, temperature) plus additive noise that models measurement errors. Performance are assessed as function of the power of the noise, that is, the variance ${\sigma }^{\mathrm{2}}$ of the Gaussian disturbance.

4.1. A Novel Metric for Accuracy Evaluation

In order to reduce the error introduced in the measurement process, a smoothing filter is typically used. However, simple approaches actually used in real networks just rely on techniques that ignore the statistical properties of the data; in particular, data are often processed through a sliding window where measurements are low-pass filtered to reduce the disturbance, similarly to a basic image denoising algorithm. Although this approach provides reasonable results on average, it may negatively have an impact on the values that deviate from the mean. The latter are conversely the most interesting data in a number of monitoring applications, for instance, to detect extreme events. As a result, it is important to measure the ability of a smoothing approach to retain the low-probability characteristics of the monitored field.

On the other hand, an algorithm that focuses too much on extreme events tends to lose its ability to “denoise” areas where the field is almost stationary, that is, plateaux with very similar values that fluctuate just because of the intrinsic uncertainty introduced by the measurement process. To correctly evaluate the overall performance, thus, it is necessary that the performance metric takes into account both these conflicting objectives, that is, ability to retain low-probability values and ability to simultaneously ensure a satisfactory smoothing.

To this aim, in the following, we propose a novel compound metric based on a weighted version of the Frobenius norm. More precisely, denoting by M the matrix of the true field values $m_{i,j}$ , where $(i,j)$ indicates a possible location of the monitored area (in this section we use matrix notation to link more clearly the measured value to the corresponding location where it has been measured; therefore we need a pair of indices, instead of the single index i formerly used to denote node-i-related quantities.), the overall error in the reconstructed field $\widehat{\mathbf{M}}$ (which is function of the measurement matrix X) can be evaluated as

$$\begin{array}{cc}& D_{\mathrm{1}}\left(\widehat{\mathbf{M}}\right)={‖\widehat{\mathbf{M}}-\mathbf{M}‖}_F,\end{array}$$

(23)

where ${‖·‖}_F$ is the Frobenius norm, that is, the sum of the squares of all elements of the matrix argument. In order to spotlight deviations in particular points, however, a weighting matrix must be introduced. We therefore define the following metric:

$$\begin{array}{cc}& D_W\left(\widehat{\mathbf{M}}\right)={‖\mathbf{W}\circ \left(\widehat{\mathbf{M}}-\mathbf{M}\right)‖}_F,\end{array}$$

(24)

where ∘ is the Hadamard (entry-wise) matrix product and W is a matrix of weights. Clearly, the metric $D_{\mathrm{1}}(·)$ is obtained as particular case of $D_W(·)$ for all-ones matrix W.

Based on the two metrics above, we propose the following compound metric:

$$\begin{array}{cc}& D\left(\widehat{\mathbf{M}}\right)={‖\widehat{\mathbf{M}}-\mathbf{M}‖}_F{‖\mathbf{I}\left(\mathbf{M}\right)\circ \left(\widehat{\mathbf{M}}-\mathbf{M}\right)‖}_F,\end{array}$$

(25)

where $\mathbf{I}(\mathbf{M})$ denotes the matrix of self-information for the field, that is, $I_{i,j}=\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{1}/P_{i,j}\right)=-\mathrm{l}\mathrm{o}\mathrm{g}{P}_{i,j}$ , where $\mathbf{P}=[P_{i,j}]$ is (an estimate of) the probability distribution of the field. This means that P is such that ${\sum }_{i=\mathrm{1}}^N{\sum }_{j=\mathrm{1}}^N{P}_{i,j}=\mathrm{1}$ . More precisely, P can be taken as the normalized histogram of the true values of the field M resulting from some binning rule. (Since this information is not available to the algorithm, it is not possible to design an estimation scheme that minimizes the metric D. The use of D here is only for the purpose of performance assessment.)

The definition in (25) can be rewritten as the product of the distance between the real field and the reconstructed one in both the original and the transformed (weighted) spaces:

$$\begin{array}{cc}& D\left(\widehat{\mathbf{M}}\right)={‖\widehat{\mathbf{M}}-\mathbf{M}‖}_F{‖{\widehat{\mathbf{M}}}^{\prime }-{\mathbf{M}}^{\prime }‖}_F,\end{array}$$

(26)

where ${\widehat{\mathbf{M}}}^{\prime }=\mathbf{I}(\mathbf{M})\circ \widehat{\mathbf{M}}$ and ${\mathbf{M}}^{\prime }=\mathbf{I}(\mathbf{M})\circ \mathbf{M}$ can be interpreted as the transformed versions (through $\mathbf{I}(\mathbf{M})$ ) of $\widehat{\mathbf{M}}$ and M, respectively; this highlights the ability of the metric to account for both aspects discussed above.

In Figure 3(a), we show a clarifying example: the true values of the field M change very slowly but for a region near the left-bottom corner, where an event has occurred that perturbs the field (e.g., accidental overheating). The histogram of the field values is reported in Figure 3(b): it shows that more than 90% of the values fall in the first bin (ranging from 0 to about 2.5), while only few points are distributed over the rest of the span (the inset reports a zoomed version of the distribution tail for a better visualization). In this example, an algorithm that completely misses to accurately represent higher values but retains the lower ones would have a very good score for the metric $D_{\mathrm{1}}$ ; however, it would perform badly through the lens of the metric $D_I$ (i.e., $D_W$ with $\mathbf{W}=\mathbf{I}(\mathbf{M})$ ), which weighs more heavily the points where the information is high. Accurate estimation of values deviating from the expected level is in fact of utmost importance in many cases, namely, for early detection of extreme events, with alarms typically triggered when exceeding a threshold. On the other hand, what we can measure is not M but rather a noisy version X as the one reported in Figure 3(c): the algorithm, thus, must be able to also smooth out the variations due to noise so as to not increase the false alarm rate due to values accidentally exceeding the threshold. To this aim, the metric $D_{\mathrm{1}}$ is more appropriate as indicator of the goodness of the smoothing process. As mentioned, these aspects are both taken into account in D, resulting in a compound metric that reflects both the estimation root mean square error and the accuracy in representing less likely (hence more informative) values.

Figure 3

Example of monitoring application scenario: (a) ground truth (matrix M), (b) histogram of the values (leading to matrix P), and (c) noisy version (measurement matrix X) for ${\sigma }^{\mathrm{2}}=\mathrm{1}$ .

We have applied the proposed Empirical Bayes approach to a sensor network deployed on the field shown in Figure 3(a). Following the cluster-based approach described above, we have run the algorithm on the noisy version resulting from the measurement process as function of the noise level ${\sigma }^{\mathrm{2}}$ . We use normalized units in order to have a direct mapping between locations and index pairs in the matrix notation. Nodes embed a technology with communication range limited to a few units: in particular, we considered measurements from nodes closer than 3 units in range, which resulted in a minimum number of neighbors equal to 11; that is, $N_i(t)\ge \mathrm{11}$ .

As mentioned, as competitor, we consider the conventional moving-average filter, which is a typical smoothing (“denoising”) technique, to improve the quality of sensed images [39]. Such an algorithm uses a spatial mask centered in each location $(x,y)$ , with constant weights, and performs the local average by moving across the whole field; this corresponds to the two-dimensional convolution between the matrix X representing the measured field and the mask and produces a low-pass filtered output where noise has been reduced, while sharp edges have been better preserved compared to nonconstant weight masks.

The result of the filtering for the case of Figure 3(c), which is related to ${\sigma }^{\mathrm{2}}=\mathrm{1}$ , is reported in Figure 4 together with the ground truth and the competitor moving average (on the same data and neighbor sets). The figure clearly shows that the proposed approach is able to smooth the measured data, reconstructing (an estimate of) the field in a way that also preserves the information about low-probability events, that is, the ones with higher values but occurring much less frequently. On the contrary, the moving average algorithm treats all locations in the same way, smoothing out too much the left-bottom corner; as a result, it completely misses the higher values of the field. This behavior is reflected by the metric D, which takes on a much smaller value for the Empirical Bayes than for the competitor. The proposed approach is able to adaptively “learn” from the data how the field varies in the local region from a statistical point of view, an information that is ultimately condensed in the parameter γ. The values of the latter at the different locations, reported in Figure 5, reveal that the algorithm can adaptively identify the points where the measurement must be prevalent (i.e., $\gamma \approx \mathrm{0}$ , dark locations), treating them differently from the ones where a smoothing can be safely performed due to the limited variability (i.e., $\gamma \approx \mathrm{1}$ , light locations). Intermediate conditions are automatically accommodated (shades of gray).

Figure 4

Result comparison for the proposed algorithm versus the moving average, ${\sigma }^{\mathrm{2}}=\mathrm{1}$ .

Figure 5

Values of γ for the different points of the field, ${\sigma }^{\mathrm{2}}=\mathrm{1}$ .

4.2. Numerical Results

It is worth noticing that the proposed Empirical Bayes approach leads to a better value of D, but it is also superior on $D_{\mathrm{1}}$ and $D_I$ individually, as it can be noticed by comparing the individual scores in Figure 4.

By increasing the noise level, it turns out that the proposed algorithm still outperforms the competitor and, furthermore, remains able to retain the highly informative (low-probability) values while ensuring a satisfactory smoothing of the noise. Results for ${\sigma }^{\mathrm{2}}=\mathrm{4}$ are reported in Figures 6 and 7. The value of D has increased, of course, but is still significantly lower than (the noisy version and) the competitor one. We notice, however, that the gap with respect to the latter is slightly less.

Figure 6

Result comparison for the proposed algorithm versus the moving average, ${\sigma }^{\mathrm{2}}=\mathrm{4}$ .

Figure 7

Values of γ for the different points of the field, ${\sigma }^{\mathrm{2}}=\mathrm{4}$ .

To investigate more thoroughly this point, we have evaluated how D varies with ${\sigma }^{\mathrm{2}}$ . For a better understanding, the metrics $D_{\mathrm{1}}$ and $D_I$ are also analyzed. We compare the proposed approach against the moving average, considering as reference the values of the metrics calculated on the raw measurements (X). In Figure 8, we report the resulting curves (from (a) to (c), D, $D_I$ , and $D_{\mathrm{1}}$ ): as clearly visible, the proposed approach outperforms the competitor in the whole span of ${\sigma }^{\mathrm{2}}$ . The two algorithms are comparable only in terms of $D_{\mathrm{1}}$ at very high noise level when the measurements are unreliable; hence, no great improvement is possible in terms of denoising compared to a low-pass filtering. However, it is also evident that the Empirical Bayes approach always exhibits the best accuracy in terms of $D_I$ . In other words, while the overall performance obviously degrade as ${\sigma }^{\mathrm{2}}$ increases, the moving average algorithm shows a roughly constant $D_I$ , meaning that even for very accurate measurements (low ${\sigma }^{\mathrm{2}}$ ) such an algorithm is not able to retain the low-probability values. Therefore, it comes as little surprise that when measurements are accurate enough the moving average approach is worse (in terms of D) than simply taking the raw measurements. The proposed approach, conversely, yields best performance by automatically adapting its parameter γ to the different contexts.

Figure 8

Comparison between the raw measurements, proposed Empirical Bayes algorithm, and the competitor algorithm (moving average): from (a) to (c), the metrics D, $D_I$ , and $D_{\mathrm{1}}$ for varying ${\sigma }^{\mathrm{2}}$ .

Finally, the case of multiple “hot spots,” that is, multiple sources of far-from-average spikes is analyzed in Figure 9. In addition to the already observed properties of the proposed approach, which remain valid, here an additional feature can be observed in terms of resilience to masquerading effects. Indeed, when multiple sources of far-from-average events are present, conventional filtering techniques like the moving average do not have enough resolution power even for moderate noise level. This is due to the averaging of points in spatial proximity, so that sources not sufficiently separated become indistinguishable and collapse onto a same blurred blob (as in the top right corner of the field in Figure 9). Conversely, the proposed approach is able to detect all sources and also to maintain the proportionality in their intensity while ensuring denoising.

Figure 9

Result comparison for the proposed algorithm versus the moving average in case of multiple far-from-average events, ${\sigma }^{\mathrm{2}}=\mathrm{1}$ .

5. Conclusions

The problem of improving the accuracy of the measurements collected by a sensor network has been addressed. Aiming at simplicity and cost-effectiveness, which are of utmost importance in application contexts, where sensor networks are often deployed, a low-complexity automatic approach has been proposed, which does not require manual setting of parameters nor recalibration. By following an adaptive (Empirical Bayes) rationale, the algorithm is able to improve the estimation accuracy by leveraging cooperation between nodes. Remarkably, it can provide a “denoised” version of the monitored field without losing accuracy in detecting less probable values. Using a novel performance metric, the capability to both reduce the measurement error and retain highly informative characteristics has been quantified, revealing that the proposed approach can outperform conventional low-pass filtering.