Distributed Estimation of a Parametric Field: Algorithms and Performance Analysis

# Distributed Estimation of a Parametric Field: Algorithms and Performance Analysis

Salvatore Talarico,  Natalia A. Schmid,  Marwan Alkhweldi, and Matthew C. Valenti,  Copyright (c) 2013 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to pubs-permissions@ieee.org. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Ruixin Niu. This work was supported by the Office of Naval Research under Award No. N000140911189. Portions of this paper were accepted for publication at the IEEE Military Communications Conference (MILCOM),Orlando, FL, Oct. 2012. The authors are with the Department of Computer Science and Electrical Engineering, West Virginia University, Morgantown, WV 26506 USA (e-mail: salvatore.talarico81@gmail.com; Natalia.Schmid@mail.wvu.edu; malkhwel@mix.wvu.edu. Matthew.Valenti@mail.wvu.edu).Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.Digital Object Identifier 10.1109/TSP.2013.2288684
###### Abstract

This paper presents a distributed estimator for a deterministic parametric physical field sensed by a homogeneous sensor network and develops a new transformed expression for the Cramer-Rao lower bound (CRLB) on the variance of distributed estimates. The proposed transformation reduces a multidimensional integral representation of the CRLB to an expression involving an infinite sum. Stochastic models used in this paper assume additive noise in both the observation and transmission channels. Two cases of data transmission are considered. The first case assumes a linear analog modulation of raw observations prior to their transmission to a fusion center. In the second case, each sensor quantizes its observation to levels, and the quantized data are communicated to a fusion center. In both cases, parallel additive white Gaussian channels are assumed. The paper develops an iterative expectation-maximization (EM) algorithm to estimate unknown parameters of a parametric field, and its linearized version is adopted for numerical analysis. The performance of the developed numerical solution is compared to the performance of a simple iterative approach based on Newton’s approximation. While the developed solution has a higher complexity than Newton’s solution, it is more robust with respect to the choice of initial parameters and has a better estimation accuracy. Numerical examples are provided for the case of a field modeled as a Gaussian bell, and illustrate the advantages of using the transformed expression for the CRLB. However, the distributed estimator and the derived CRLB are general and can be applied to any parametric field. The dependence of the mean-square error (MSE) on the number of quantization levels, the number of sensors in the network and the SNR of the observation and transmission channels are analyzed. The variance of the estimates is compared to the derived CRLB.

{keywords}

Wireless sensor network, EM algorithm, maximum-likelihood estimation, Cramer-Rao lower bound, distributed parameter estimation.

## I Introduction

\PARstart

Many military and civilian applications use distributed sensor networks as a platform to perform environmental monitoring, surveillance, detection, tracking, object classification and counting. Typically, such applications impose constraints on power, bandwidth, latency, and/or complexity. Numerous research findings have been reported on these topics over the past two decades (for example, see the reviews in [1, 2]). Some publications are concerned with a specific application [3, 4, 5, 6, 7, 8, 9]; other papers have formulated and solved problems in the area of distributed estimation, detection and tracking [10, 11, 12, 13, 14, 15, 16]. In particular, an important research thrust in the field of distributed estimation is the design of practical distributed algorithms, choosing either to optimize a distributed sensor network with respect to energy consumption during transmission [17, 18, 19] or to impose bandwidth constraints and thus focus on designing an optimal quantization strategy for the distributed network [20, 21, 22, 23]. Alternatively, both constraints could be considered [24, 25]. For instance, Wu et al. [18] studied the problem of minimizing the estimation error, mean square error (MSE), under the constraint of limited power; based on the work in [18], Li and Al-Regib [19] developed an upper bound on the lifetime of a wireless sensor network and then proposed a methodology to increase it; Xian and Luo [20] provided a distributed estimation scheme that deals with low bandwidth channels efficiently; and Cui et al. [25], imposing constraints on both power and bandwidth, proposed a transmission scheme that results in considerable power savings.

Among research groups working on the problem of distributed estimation, there are a few dealing with distributed estimation of a field (a multidimensional function, in general) [26, 27, 28, 29, 30]. In many real-world applications, distributed estimation of a multidimensional function may provide additional information that aids in making a high-fidelity decision or in solving another inference problem. Motivated by this, we contribute to this topic by formulating and solving the problem of a parametric field estimation from sparse noisy sensor measurements using an iterative solution. We also focus on the development of theoretical limits for distributed estimation of the parameters associated with a parametric field. These limits are used to compare the performance of the estimates obtained to the best performance that can be achieved in theory.

In this paper, the problem of distributed estimation of a physical field from sensory data collected by a homogeneous sensor network is stated as a maximum-likelihood estimation problem. The physical field is a deterministic function and has a known spatial distribution parameterized by a set of unknown parameters, such as the location of an object generating the field or the strength and spread of the field in the region occupied by the sensors. It is assumed that white Gaussian noise is added at the sensors and over the transmission channel from the sensors to a fusion center, and that signal processing is performed both locally at the sensor nodes as well as globally at the fusion center.

Two cases of transmission channels are considered. The first case, which we refer to as a amplify-and-forward, assumes that noisy sensor measurements are transmitted as analog signals to a fusion center. The second case, which we refer to as a quantize-and-forward, assumes that each sensor locally quantizes its observation to levels, and the quantized data are digitally modulated and communicated to a fusion center. For both cases, orthogonal additive white Gaussian noise channels are considered and the effect of interference is neglected, since we assume an interference avoidance protocol is used.

An iterative algorithm to estimate unknown parameters is formulated, and a simplified numerical solution involving additional approximations is developed for both cases of the transmission channel. Furthermore, a new transformed expression for the CRLB on the variance of distributed estimates of field parameters is derived. The applied transformation reduces a multidimensional integral to a single infinite series. Although it is an infinite series, only the first few terms need to be evaluated in order to obtain a reasonable accuracy. The developed CRLB is applied to evaluate the performance of the estimates of the field parameters for both types of transmission channels. The developed expectation-maximization (EM) solution and the bound are general. However, this work uses a Gaussian bell function as the parametric field to illustrate the EM approach and the derived bound. Our numerical evaluation shows that the variance of the EM estimates approaches the CRLB as the density of the network increases, provided initial values of the estimates are selected sufficiently close to the true parameters. The accuracy, performance, and the rate of convergence of the developed EM algorithm are compared to those of the Newton-Raphson (NR) algorithm. Numerical results highlight the advantages and disadvantages of the developed EM algorithm.

The rest of the paper is organized as follows. Sec. II describes the network. Sec. III describes the proposed EM solution. Sec. IV develops the CRLB on the variance of parameter estimates. Sec. V provides numerical performance evaluation of the estimator and compares the variance of the estimated parameters to the CRLB. Finally, Sec. VI presents a summary of the results.

## Ii Network Model

Consider a distributed network of homogeneous sensors monitoring the environment for the presence of a substance or an object. Assume that each substance or object is characterized by one or more location parameters and by a spatially distributed physical field generated by them. Examples of physical fields include (1) a magnetic field generated by a dipole and modeled as a function of the inverse cube of the distance to the dipole [31, 32, 33], (2) a radioactive field modeled as a stationary spatially distributed Poisson field with a two-dimensional intensity function decaying according to the inverse-square law [34] or (3) a cloud of pollution or chemical fumes [35] that, if stationary, has a spatial intensity that can often be modeled as a Gaussian bell.

A block diagram of a distributed sensor network used for the estimation of parameters of a physical field is shown in Fig. 1. The network is composed of sensors randomly distributed over a finite area . The network is calibrated, and the relative locations of the sensors are known, which are common assumptions [13, 22]. Sensors act independently of one another and take noisy measurements of a physical field defined over the area A sample of at a location is denoted as . The parametric field is characterized by unknown parameters We use to emphasize this dependence and use for brevity. The sensor noise at different locations, denoted by , is modeled as independent Gaussian random variables with zero mean and known variance . Let , be the noisy samples of the field at the location of distributed sensors. Then is modeled as These observations are transmitted over noisy parallel channels to a fusion center (FC), which estimates the vector parameter We assume that the communication channels do not interfere. The method used to send these observations is chosen accordingly to the application constraints and the channel characteristics, and it is denoted in Fig. 1 by .

If the sensor observations are transmitted as analog samples over the channels without any prior processing (aside from linearly modulating with a suitable carrier frequency and an amplification to assure the transmit power is at a desired level), we refer to it as a amplify-and-forward. Denote by the noisy observations received by the FC and obtained by coherently demodulating the signal sent by the sensors, as shown in Fig. 1. In this case, the observations are given as , , where the noise in the channel is white Gaussian with variance , and and are independent 111The channel gain is normalized to unity and the effect of fading is absorbed into the variance of the transmission channel noise..

Due to constraints imposed by practical technology, each sensor may be required to quantize its measurements prior to transmitting them to the FC. Assume a deterministic quantizer with quantization levels at each sensor location. Let be the reproduction points of a quantizer and be the boundaries of the quantization regions. Then the output of the -th quantizer is a random variable taking value with probability

 pk,j(θ) = ∫τj+1τj1√2πσ2kexp(−(t−Gk)22σ2k)dt, (1)

which depends on the unknown parameters of the physical field .

The quantized data are modulated using a linear digital modulation scheme, such as on-off keying (OOK), an M-ary quadrature amplitude modulation (QAM) or a pulse amplitude modulation (PAM), and then transmitted to the FC over noisy parallel channels. We refer to this channel as a quantize-and-forward. To further clarify the modulation and demodulation steps, denote by a bit representation of the quantized observation of the -th sensor in vector form. Then by applying a linear modulation scheme (OOK or QAM, for example) and then coherently demodulating it, the vector of observed data at the FC will be in the form where is a white Gaussian noise vector independent of

## Iii Distributed parameter estimator

Given the noisy measurements and the relative location of the sensors in the network, the task of the FC is to estimate the vector parameter In this work, the maximum-likelihood (ML) estimation approach [36, 37] is used to solve the problem of distributed parameter estimation for both channels. Denote by the log-likelihood function of the vector of the random measurements parameterized by vector , where is the vector of measurements To be more specific, is defined as with being the parameterized probability density function (pdf) of the vector The ML solution is the vector that maximizes the log-likelihood function :

 ^θ = argmaxθ∈Θl(Z:θ), (2)

where is a set of admissible vector parameters.
The necessary conditions to find the maximizer are given by:

 (3)

where denotes the gradient with respect to the vector and the maximizer is the interior point of

### Iii-a Amplify-and-forward channel

For the amplify-and-forward channel, the signals received at the FC are independent but not identically distributed and the log-likelihood function is given as:

 l(Z:θ) = K∑k=1log(fZk(Zk)) (4) = −12K∑k=1[Zk−G(xk,yk:θ)]2(σ2k+η2k)+C1,

where is not a function of .

The condition (3) applied to (4) generates a set of nonlinear equations in The solution to (2) is found numerically by means of Newton’s method [38]. This solution is applied to generate the results in Sec. V-A.

### Iii-B Quantize-and-forward channel

For the quantize-and-forward channel, the joint likelihood function of the independent quantized noisy measurements can be written as

 l(Z:θ) = (5)

where is defined in (1), is a binary representation of the integer and, thus, of , characterized by bits; is not a function of . The ML solution is the solution that maximizes the expression (5), but since applying (3) to (5) results in a set of nonlinear equations, an iterative solution to the problem can be developed: (1) a set of EM iterations [39] are formulated and then (2) a Newton’s linearization is used to solve EM equations for the unknown parameters.

#### Iii-B1 Expectation Maximization Solution

The random variables , are selected as complete data. The complete data log-likelihood, denoted by , is given by

 lcd(R,N)=−12K∑k=1(Rk−Gk)2σ2k+C3, (6)

where is not a function of The measurements form incomplete data and the mapping from complete data space to incomplete data space is given by

Denote by an estimate of the vector obtained at the -th iteration. To update the estimates, the expectation and maximization steps are alternated. During the expectation step, the conditional expectation of the complete data log-likelihood is evaluated as follows:

 T(m+1) = E[−12K∑k=1(Rk−Gk)2σ2k∣∣ ∣∣Z,^θ(m)], (7)

where the expectation is with respect to complete data, given the incomplete data (measurements) and the estimates of the parameters at the -th iteration. During the maximization step, (7) is maximized:

 dT(m+1)dθt=E[K∑k=1(Rk−Gk)σ2k(dGkdθt)∣∣ ∣∣Z,^θ(m)]∣∣ ∣∣^θ(m+1)=0, t=1,…,L. (8)

To find the expectation, it can be noted that the conditional probability density function of given is Gaussian with mean and a diagonal covariance matrix with identical diagonal entries and the pdf of is Gaussian with mean and variance It can be also noted that at the -th iteration the conditional pdf of given implicitly involves the estimates of the parameters obtained at the -th iteration.

Denote by the estimate of the field at the location with the vector of parameters replaced by their estimates Then the following lemma can be formulated.

###### Lemma 3. 1

The expression for the iterative evaluation of the unknown parameters can be written as:

 K∑k=1dG(m+1)kdθtA(G(m)k)−K∑k=1G(m+1)kdG(m+1)kdθtB(G(m)k)=0,t=1,…,L, (9)

where

 A(G(m)k) = M∑j=1exp(−(zk−bj)T(zk−bj)2η2k)f(m)Zk(zk)(2πη2k)α/2⎛⎜ ⎜⎝√σ2k2πe−(τj−G(m)k)22σ2k (10) −√σ2k2πe−(τj+1−G(m)k)22σ2k+G(m)kΔT(m)(j,k)⎞⎟ ⎟⎠,
 B(G(m)k)=M∑j=1exp(−(zk−bj)T(zk−bj)2η2k)f(m)Zk(zk)(2πη2k)α/2ΔT(m)(j,k), (11)

with

 ΔT(m)(j,k)=Q⎛⎝τj−G(m)kσk⎞⎠−Q⎛⎝τj+1−G(m)kσk⎞⎠, (12)
 f(m)Zk(zk)=∫f(m)Zk|R(zk|r)f(m)Rk(r)dr. (13)

The expression is used to denote the Q-function, that is given by

 Q(x)=1√2π∫∞xexp(−t22)dt. (14)

Proof: The details of the derivation are moved to the Appendix A.

#### Iii-B2 Linearization

The equations (9) are nonlinear in and have to be solved numerically for each iteration. To simplify the solution, the expression in (9) is linearized by means of Newton’s method. Denote by the vector form of the first derivative of the left side in (9) over and let be the Jacobian of the left side in (9) . The index indicates the iteration of the Newton’s solution. Then solves the following linearized equation:

 J(θ(m+1)n)[θ(m+1)n+1−θ(m+1)n]=−F(θ(m+1)n). (15)

#### Iii-B3 Alternative solution

For a parametric field estimation with the network model described in Sec. II the Newton-Raphson (NR) algorithm [40] can be used as an alternative to the EM algorithm. The EM and NR algorithms have a similar computational complexity that becomes imperceptible as the network becomes sparse and/or as the number of quantization levels decreases. The NR algorithm displays a faster convergence compared to the EM algorithm. However, the proposed EM solution converges more consistently to the local maximum of the likelihood function [41] and guarantees a better accuracy. Sec. V-C provides an illustration, where the EM and NR algorithms are compared both in terms of rate and consistency of convergence as well as in terms of precision of the final results.

## Iv Cramer-Rao Lower Bound

In Sec. III, parameters of the physical field were estimated for the amplify-and-forward and the quantize-and-forward channel by means of a Newton’s method and a linearized version of EM algorithm, respectively. In order to evaluate the efficiency of the estimator for both channels the CRLB on the variances of unknown parameters is evaluated.

In this section, we assume that the parameter estimates are unbiased. Denoted by the covariance matrix of the estimated parameters the CRLB on [36, 37] is given as

 Σθ≥I−1(θ), (16)

where denotes the Fisher information matrix with the entry at the location given by

 Is,t = −E[∂2l(Z:θ)∂θs∂θt]. (17)

Below we detail the derivation of the CRLB for the amplify-and-forward and the quantize-and-forward channel.

### Iv-a Amplify-and-forward channel

In the case of amplify-and-forward channel, substituting (4) into (17), and simplifying all the terms that do not depend on , yields

 Is,t=12K∑k=11(σ2k+η2k)E[∂2∂θs∂θt(Zk−G(xk,yk:θ))2]. (18)

After taking the partial derivatives and noting that , (18) becomes

 Is,t=−K∑k=11(σ2k+η2k)∂G(xk,yk:θ)∂θs∂G(xk,yk:θ)∂θt. (19)

### Iv-B Quantize-and-forward channel

In the case of quantize-and-forward channel, substituting (5) into (17) yields

 Is,t=−E[∂2∂θs∂θtK∑k=1logxk(θ)], (20)

where is given by the following expression

 (21)

After taking the partial derivatives in (20) we have

 Is,t=−K∑k=1E[1xk(θ)∂2xk(θ)∂θs∂θt−1x2k(θ)∂xk(θ)∂θt∂xk(θ)∂θs]. (22)

Substituting (21) into (22) yields

 (23)

where

 (24) Φk,j,i=E⎡⎢ ⎢ ⎢ ⎢ ⎢⎣exp(−(Zk−bj)T(Zk−bj)+(Zk−bi)T(Zk−bi)2η2k)[∑Mv=1pk,v(θ)exp(−(Zk−bv)T(Zk−bv)2η2k)]2⎤⎥ ⎥ ⎥ ⎥ ⎥⎦. (25)

Since the expectation in (24) is with respect to the pdf of which is given by

 fZk(zk)=M∑v=1pk,v(θ)(2πη2k)α/2exp(−(zk−bv)T(zk−bv)2η2k), (26)

the expression (24) integrates to 1, that is,

###### Lemma 4. 1

The entry of the Fisher information matrix (16) can be transformed in (27) (displayed at the top of this page) where is given by (40).

Proof: The details of the proof are presented in the Appendix B.

Thus the entry of the Fisher information matrix can be implemented using two approaches: (i) by evaluating numerically (20) and (ii) by truncating the infinite sum in (27) to terms such that a compromise between computational complexity and accuracy is achieved. We involve Simpson’s method [42] as an alternative approach to evaluate It provides a good tradeoff between accuracy and speed. The results from the two different implementations coincide even when very few terms are used in the series representation (27). This is demonstrated at the end of Sec. V-B that presents a comparison of the two implementations.

## V Numerical Analysis

In this section, the performance of the distributed ML estimator is evaluated for both types of channels. Numerical results are presented for the case of the field modeled as a Gaussian bell. However, the distributed ML estimator in Sec.III-A and Sec.III-B and the bounds in (19) and (27) are general and can be applied to estimates of any parametric field

Our numerical analysis assumes that a distributed network of sensors is formed by deploying them uniformly at random over a finite area of size , where the location of each sensor is noted. The Gaussian field used in our simulations is

 G(xk,yk)=hexp[−(xk−xc)22ρ2x−(yk−yc)22ρ2y], (28)

where is the “strength” of the field, and determine the “spread” of the bell in the x and y direction, respectively, and is the position of the object generating the field. In the numerical examples, we set , The location parameters of the field are set to and . For numerical illustration we assume that , , , and are all unknown parameters that have to be estimated; i.e., in our experiments the unknown vector-parameter is . The size of the network is varied from to . A Gaussian field is sampled at the location of the -th sensor, , and a sample of randomly generated Gaussian noise with mean zero and variance is added to each field measurement. Note that our experiments assume i.i.d. noise samples, and, for simplicity, and for all sensors. The noise variance is selected such that the total signal-to-noise ratio (SNR) of the local observations defined as

 SNRO=∫∫AG2(x,y:θ)dxdyAσ2 (29)

takes a predetermined value. For the amplify-and-forward channel, the variance of the noise in the transmission channels is selected such that the total SNR in the channels defined as

 SNRC=∫∫AG2(x,y:θ)dxdyAη2+σ2η2 (30)

is also set to a predefined value. For the quantize-and-forward channel, each sensor observation is quantized to one of levels using a uniform deterministic quantizer. parallel white Gaussian noise channels add samples of noise with variance selected to set the total SNR, that is defined as

 SNRC=∫∫AE[q2(R(x,y:θ))]dxdyAη2, (31)

to a specific value. Note that (31) converges asymptotically to (30) when the number of quantization levels tends to infinity.

It is assumed that the FC observes the noisy quantized samples of the field. The function in (31) is a quantized version of Note that due to the symmetry of the experimental set up and due to the statistical averaging, the results for and and furthermore the results for and are very similar. Therefore, to preserve the space, convergence of solutions of iterative algorithms is demonstrated for the case of , and only.

### V-a Amplify-and-forward channel

First, the convergence of the ML estimator numerically evaluated by means of Newton’s method is illustrated. The estimated values of the x-location, the standard deviation and the strength of the field are plotted in Fig. 2(a) as a function of the number of iterations. The functions in Fig. 2(a) are parameterized by different values of the SNR in the transmission and observation channels. This illustration is based on a single realization of the distributed network with , when the initial values are picked to be for the x-location, for the standard deviation and for the strength of the field. Note that for the case of dB, the estimated values converge to the real values after 14 iterations. For the other two cases ( dB and dB) the increasing discrepancy between the estimated and the true values is due to a lower sensor density in the network and also due to increasing variances of observation and transmission noise. Note how large is the deviation of the asymptotic value of , , and from the real values when both and are set to dB.

To further analyze the estimation performance, the squared error (SE) between the estimated and true location parameters is evaluated. The SE is defined as

 SE=L∑i=1(^θi−θi)2

and the mean square error (MSE) is evaluated numerically by means of Monte Carlo simulations.

In this and the following subsections, we involve a box plot to illustrate the dependencies of the MSE on several parameters. The central mark in each box is the median. The edges of a box present the th and th percentiles. The dashed vertical lines mark the data that extend beyond the two percentiles, but not considered as outliers. The outliers are then plotted individually and marked with a “+” sign.

The dependence of the MSE on the number of sensors, in the distributed network for the case of dB is shown as a set of box plots in Fig. 3(a). The dependence of the MSE on the SNR of the observation and transmission channels, when the number of sensors is fixed to , is displayed as a set of box plots in Fig. 4(a). Each box in Fig. 3(a) and Fig. 4(a) is generated using Monte Carlo realizations of the network and ML runs. To take a closer look at the distribution of outliers, we define the probability of outliers as a probability that SE exceeds a positive valued threshold Denote by the probability of outliers at threshold Then mathematically is defined as

 PO(τ)=P[SE>τ]. (32)

We vary the value of the threshold and display the percentage of outliers as a function of in Fig. 5(a). Note the large percentage of outliers for small values of This corresponds to the case when one or more of the five parameters (x-location, y-location, the standard deviation and and the strength of the field ) did not converge to its true value.

The convergence of the Newton’s iterations to the true values of the location parameter is analyzed with respect to a choice of initial guess for all five unknown parameters (). In particular the initial values for the iterative solution due to Newton’s method are selected randomly for each new Monte Carlo realization and for each unknown parameter within eight regions, that we indicate with We have enumerated these regions starting from the one closest to the true values. In particular the possible initial values for the unknown parameter, that can be chosen inside the -th region, are in the interval , where is the true value for the parameter. Fig. 6(a) shows a box plot of the MSE between the estimated and true parameters of the field displayed as a function of the regions within which the initial values for the Newton’s method are selected at random. Each box plot is obtained by using Monte Carlo realizations of the network and ML runs. The number of sensors is fixed and equal to The variance and are chosen such that the SNR for the observation and transmission channels is each equal to 15 dB. Fig. 6(a) can be interpreted as a sensitivity analysis of the implemented ML solution due to Newton’s method. As Fig. 6(a) demonstrates, it is still possible to estimate the vector of unknown parameters with a relatively low value of MSE, even when the initial values are selected far apart from the true values, but after a certain region the Newton’s method is not able to converge anymore and the MSE starts to increase rapidly. This effect is due to the presence of multiple local maxima in the log-likelihood function given by (4).

Finally, Fig. 7(a) shows the variance of the estimated parameters , , and obtained numerically by means of the Newton’s method. The results are averaged over Monte Carlo realizations. The plots are compared to the CRLB displayed as a function of the number of sensors when dB. The results in Fig. 7(a) demonstrate that for the amplify-and-forward channel, the empirical variance of the estimated parameters converges to the values obtained by means of the CRLB, underlining that with only a few distributed measurements () the Newton’s iterative solution is efficient.

### V-B Quantize-and-forward channel

Fig. 2(b) shows the estimated values obtained by the ML estimator in (9) for the -location, the standard deviation and the strength of the field as a function of the number of EM iterations. Different numbers of quantization levels have been considered. For this plot, a single realization of a sparse distributed network composed of is considered, the initial values are picked to be for the -location, for the standard deviation and for the strength of the field. The SNR for the observation and transmission channel is fixed to dB. Fig. 2(b) illustrates the convergence of the EM algorithm, pointing out the role of the quantization levels: when the EM algorithm converges to the true values after only few iterations. For the other two cases, the EM algorithm does not converge as fast, and there is a larger discrepancy between the estimated and the true values of the parameters , and . This discrepancy grows as the number of quantization levels decreases, as it is expected. In addition to rough quantization, it is affected by the low sensor density in the network and by the noise in observation and transmission channels.

In order to analyze the performance of the estimator, in the same way it has been done for the amplify-and-forward channel, the MSE between the estimated and true parameters is used as a performance metric. Fig. 3(b) shows the dependence of the MSE on the number of sensors, in the distributed network for the case of dB. Fig. 4(b) demonstrates the dependence of MSE on the varying values of SNR in the observation and transmission channels. The results are shown for the case of Fig. 4(b) is generated considering a sparse network composed of sensors. All three plots are generated using Monte Carlo realizations of the network and EM runs. Fig. 4(b) indicates that the noise in observation channel (expressed as ) prevails over the noise in transmission channel (expressed as ) in terms of its effect on the estimation error (the average SE and its variance).

The percentage of outliers due to divergence of the EM algorithm is illustrated in Fig. 5(b). A sparse network composed of sensors is considered for and three different combinations of and are analyzed. The three combinations are (1) dB dB, (2) dB, and (3) dB dB. Fig. 5(b) highlights that the observation channel not only affects more the performance in terms of SE than the transmission channel, but it also affects the convergence of the EM algorithm in a more tangible way, increasing the probability of outliers. This is a consequence of quantizing very noisy measurements.

In order to analyze the robustness of the EM algorithm to the initial values of estimates, eight regions , inside which the initial values are chosen randomly, are considered. In particular the possible initial values for the unknown parameter, that can be chosen inside the -th region, are in the interval , where is the true value for the parameter. Fig. 6(b) shows a box plot of the MSE between the estimated and true parameters of the field displayed as a function of the regions within which initial values of parameters are selected randomly. Each box is generated by using Monte Carlo realizations of the network and EM runs. Here the network is composed of sensors, the number of quantization levels is set to , and the variances and are selected such that the SNR for the observation and transmission channels is equal to 15 dB. Fig. 6(b) demonstrates that the EM algorithm is not very sensitive to the choice of the initial value of estimates. Note that up to the -th region the median value of MSE is reasonably low. The abrupt increment of the MSE is caused by the presence of multiple local maxima in the log-likelihood function given by (5).

Fig. 9 shows the effect of the number of quantization levels on the MSE, when the variance and are chosen such that the SNR for both observation and transmission channels is set to dB. As expected the MSE decreases when the number of quantization levels is increased. As the number of quantization levels at local sensors increases, the percentage of this performance improvement decreases and tends to converge to the case when raw observations are transmitted. Nevertheless the distributed parameter estimation system achieves an acceptable performance in terms of MSE even when and for a sparse network composed of sensors. This emphasizes on the energy efficiency of the proposed parameter estimation framework that could lead to a higher lifetime of distributed sensors in the network. In other words, local sensors do not need to waste a lot of energy to send high-resolution quantized observations to achieve an acceptable estimation performance in terms of the MSE.

Fig. 7(b) shows the variance of the estimated parameter , and , obtained by means of the EM algorithm for the case of The results are averaged over 1000 Monte Carlo realizations. The plots are compared to the CRLB displayed as a function of the number of sensors. Fig. 7(b) demonstrates that for the quantize-and-forward channel, the EM solution is efficient. The convergence rate of variance of the estimated parameters to the value provided by the CRLB is slightly slower compared to the similar plots in the case of the amplify-and-forward channel.

Lastly, Fig. 9 shows the dependence of the truncated CRLB on the number of terms retained in the infinite series (27). The results of truncation are compared to those obtained by means of numerical integration using the Simpson’s method. The plots are obtained for the case of , and dB. Fig. 9 shows that by truncating the infinite series in (27), it is possible to achieve the same values of the CRLB as those achieved by the Simpson’s method.

### V-C Comparison Between EM and NR

Fig. 11 shows the estimated values for the -location, the standard deviation and the strength of the field as a function of the number of iterations for both the EM and the NR algorithms. For this plot, a single realization of a sparse sensor network composed of is considered, the initial values are picked to be for the -location, for the standard deviation and for the strength of the field for both algorithms. The SNR in the observation and transmission channels is fixed to dB, and the number of quantization levels is fixed to . Fig. 11 shows that when both algorithms converge, the NR algorithm, as stated in Sec. III-B3, requires fewer iterations to reach the final values compared to the EM algorithm, underlining the advantage of the NR algorithm in terms of convergence rate.

In order to compare the two algorithms in terms of their accuracy, the MSE between the estimated and the true parameters is evaluated. Fig. 11 shows the MSE for both EM and NR algorithms. For this example, the network is composed of sensors, the SNRs are set to dB and . Fig. 11 shows that the EM algorithm has lower MSE compared to the NR. To further analyze the performance of the algorithms, Fig. 12 displays the percentage of outliers as a function of the threshold (see equation (32)) for both algorithms and for three different combinations of and : (1) dB, (2) dB, and (3) dB.

## Vi Summary

This paper has presented a distributed ML estimation procedure based on an iterative solution for estimating a parametric physical field. Furthermore, it has also detailed the derivation of a transformed expression for the CRLB on the variance of distributed estimates of field parameters by a homogeneous sensor network. The model of the network assumed independent sensor measurements and transmission over a noisy environment. Two channel models were considered: (1) the measurements from the sensors were sent directly to the FC by means of a linear analog modulation; (2) each sensor quantized its measurement to levels and the quantized and encoded data were communicated to the FC over parallel additive white Gaussian channels. The stability of the distributed parameter estimator has been analyzed for both models and also its robustness to the initial values of estimates has been considered. The results have shown that for the quantize-and-forward channel the SNR of the observation channel dominates the SNR of the transmission channel in terms of the values of MSE and it also affects the convergence of the EM algorithm in a more tangible way, increasing the probability of outliers. The developed CRLBs were further compared to the numerical values of the variance of estimates. The results showed that the estimates are nearly efficient for small and become efficient for both channel models when the density of the sensor network increases. We have demonstrated numerically that the derived EM algorithm generates more accurate estimates compared to the NR method. Finally, we have shown that truncating the infinite sum in the developed CRLB to only few first terms produces highly accurate results.

As a future work, the parametric field will be replaced by a mixture model or nonparametric field (for generality), mimicking the case of unknown number of multiple objects generating a cumulative physical field sensed by a distributed sensor network within an area The models for transmission channels will involve fading and shadowing effects. Efforts to eliminate the FC and make the sensor network decentralized will be made.

## Appendix A

This section provides details leading to the equation (9). Consider the -th term under the sum in (8):

 E[(Rk−Gk)∂Gk∂θt∣∣∣zk,^θ(m)]
 =∫+∞−∞(rk−Gk)∂Gk∂θtexp(−(rk−G(m)k)22σ2k)f(m)Zk(zk)√2πσ2k
 =M∑j=1∫τj+1τj(rk−Gk)∂Gk∂θtexp(−(rk−G(m)k)22σ2k)f(m)Zk(zk)√2πσ2k
 ×exp(−(zk−bj)T(zk−bj)2η2k)(2πη2k)α/2drk
 =M∑j=1exp(−(zk−bj)T(zk−bj)2η2k)f(m)Zk(zk)(2πη2k)α/2∂Gk∂θt
 ×∫τj+1τj(rk−Gk)exp(−(rk−G(m)k)22σ2k)√2πσ2kdrk.

Note that the difference in the last integral can be rewritten as Then

 ×⎧⎪ ⎪⎨⎪ ⎪⎩1√2πσ2k∫τj+1τjexp⎛⎝−(rk−G(m)k)22σ2k⎞⎠d(rk−G(m)k)22
 +(G(m)k−Gk)1√2πσ2k∫τj+1τjexp⎛⎝−(rk−G(m)k)22σ2k⎞⎠drk⎫⎪ ⎪⎬⎪ ⎪⎭.

Replacing the last integral with a difference of two Q-functions we obtain:

 K∑k=1E[(Rk−Gk)∂Gk∂θt∣∣∣zk,^θ(m)]=K∑k=1M∑j=1exp(−(zk−bj)T(zk−bj)2η2k)f(m)Zk(zk)(2πη2k)α/2
 ×∂Gk∂θt⎧⎪ ⎪⎨⎪ ⎪⎩σ2k√2πσ2k⎧⎨⎩exp⎛⎝−(τj−G(m)k)22σ2k⎞⎠
 −exp⎛⎝−(τj+1−G(m)k)22σ2k⎞⎠⎫⎬⎭+(G(m)k−Gk)
 ⎧⎨⎩Q⎛⎝τj−G(m)kσk⎞⎠−Q⎛⎝τj+1−G(m)kσk⎞⎠⎫⎬⎭⎫⎬⎭∣∣ ∣∣Gk=G(m+1)k=0.

## Appendix B

This section provides details leading to a closed form expression for the expectation in (25).

After some basic manipulations the right side of (25) can be represented in the following form:

 Φk,j,i = ∫...∫Rα1(2πη2k)α/2 (33)

The denominator in (33) can be further expressed as a series (see page 15 in [43]):

 ((x(θ)−1)+1)−1 = ∞∑n=0(−1)n(x(θ)−1)n. (34)

After replacing the denominator of the integrand in (33) with the expression in the right hand side of (34), we have

 Φk,j,i = ∞∑n=0(−1)n∫...∫Rα1(2πη2k)α/2 (35) ×exp(−(zk−bj)T(zk−bj)+(zk−bi)T(zk−bi)2η2k) ×[M∑v=1pk,v(θ)exp(−(zk−bv)T(zk−bv)2η2k)−1]ndzk.

Applying the binomial theorem [43]

 (a+b)n = n∑m=0(nm)an−mbm, (36)

to the power of term in (35), yields:

 Φk,j,i = ∞∑n=0(−1)nn∑m=0{(−1)m(nm) (37) ×∫...∫