Localization of acoustic sources using a decentralized particle filter

This paper addresses the decentralized localization of an acoustic source in a (wireless) sensor network based on the underlying partial differential equation (PDE). The PDE is transformed into a distributed state-space model and augmented by a source model. Inferring the source state amounts to a non-linear non-Gaussian Bayesian estimation problem for whose solution we implement a decentralized particle filter (PF) operating within and across clusters of sensor nodes. The aggregation of the local posterior distributions from all clusters is achieved via an enhanced version of the maximum consensus algorithm. Numerical simulations illustrate the performance of our scheme.


Background and state of the art
In this paper, we use a physics-based model and a Bayesian approach to develop a decentralized particle filter (PF) for acoustic source localization in a sensor network (SN).In a decentralized PF, the processing is done locally at the sensors without using a fusion center.Thereby, the estimated position is known at every sensor in consequence of this decentralized process.
The problem formulation in this paper is motivated by indoor localization of an acoustic source.A hallway is modeled including basic boundary conditions for windows (membranes) and walls.
The source localization problem has been studied, e.g., in [1][2][3], [ [4], p. 4089 ff], [ [5], p. 746 ff] and [6], all of which use a sequential Bayesian estimator [7] to infer the source position states from observations using multiple sensors.These papers build on a state-space transition equation describing the global source state trajectory over time and the measurement equation between these states and the measurements.The underlying model of the physical process is modeled in the measurement equation.A decentralized approach aims at identifying global source states that are common to all decentralized units.Each decentralized unit typically consists of a sensor and a Bayesian estimator associated with the sensor's neighborhood.
A different approach consists of incorporating the partial differential equation (PDE) describing the dynamics of the physical process.In source tracking applications, this implies that the field itself becomes part of the state, which thus is distributed over all space.For instance, the acoustic wave field is described by a hyperbolic PDE for pressure and hence the state vector comprises the spatio-temporal pressure field.This approach is used in (ocean) acoustic models [8,9] and geophysical models [ [4], p. 4089 ff], [10][11][12][13].For localization, the model is augmented with a source model providing a relation between global source states, e.g., position, and distributed field states, i.e., pressure.
Our approach belongs to the realm of the second approach.The novel aspects include the formulation of a source model suitable for distributed processing, the design of a distributed particle for the estimation of the posterior distribution of field and source states, and the development of a modified version of the maximum consensus (MC) algorithm [14] for the maximum a-posteriori (MAP) estimation of the source location.For several loosely connected agents, a consensus algorithm makes the agents to converge to a group decision based on local information.

Contributions and outline
We consider inference problems in space r and time t which are modeled via partial differential equations (PDE) of the form [15,16] L{p(r, t)} = s(r, t), (1a) where L denotes the PDE operator; B the boundary/ initial conditions; p(r, t) the quantity of interest, and s(r, t) the source term.If the PDE parameters, the source term, and the boundary/initial conditions are known, determining p(r, t) is the forward problem.In contrast, inverse problems amount to estimating PDE parameters or states like source locations from measurements of p (r, t).
Discretization of (1) and extending to a stochastic process leads to the state transition equation where k denotes discrete time, x k is the state vector (incorporating samples of p(r, t)), u k is the input vector (here corresponding to the sources), w k is a random noise vector (process noise), and g k is the state transition mapping.The state transition equation is complemented by the measurement equation Here, y k is the observation, v k denotes the measurement noise, and the mapping h k characterizes the measurement.Taken together, (2) and (3) constitute the state-space model, see Section II and [17].
In the Gaussian case, Bayesian estimation based on the state-space model ( 2), (3) leads to various kinds of Kalman filters [1,7,18].Here, the Bayesian estimator builds on the particle filter (PF) framework due to (i) various possible geometries and (ii) the non-linearity of the state-space model.After discussing a centralized PF for source localization and tracking in Section III, we develop a decentralized implementation of the PF by splitting the nodes of the sensor networks (SN) into clusters.The clustered SN architecture entails a corresponding decomposition of the state-space model, and the decentralized PF performs intra-cluster computation and inter-cluster communication on the decomposed state-space model (see Section IV).
The decentralized PF yields local posterior distributions within each cluster.Localization of the acoustic sources amounts to finding the maxima of the global posterior distribution.To this end, we propose a modified maximum consensus algorithm in Section V.After a summary in Section VI, in Section VII, we describe extensive numerical simulations that illustrate the properties and performance of our source localization method.

II System model
In this section, we develop a state-space model from the PDE of the spatio-temporal acoustic field using the finite difference method (FDM) [15,19] to obtain a discretization in space and time.

A Forward model-spatio-temporal field
In the following, we consider an acoustic problem characterized by the hyperbolic PDE (scalar wave equation) [16,19,20]: Here, p(r, t) denotes pressure, ∂ t is the partial derivative with respect to time, ∇ 2 the Laplace operator, c the sound speed, s(r, t) is the source, and Ω ⊂ ℝ 2 is the 2-D region of interest.Hereafter, let the initial conditions be From three basic boundary conditions we use (4d) and (4f) for modeling a hallway.∂Ω 1 is the transparent part of the boundary of Ω (with normal vector n) modeling an infinite domain for the behind uncovered area.The boundary ∂Ω 2 (disjoint from ∂Ω 1 models windows, whereas ∂Ω 3 (disjoint from ∂Ω 1 and ∂Ω 2 ) models walls.The choice of these boundary conditions indeed affects the resulting state-space model but does not change the general formulation of the decentralized approach.

B Finite difference method
To obtain a space-time-discrete model, the differential operators are approximated by finite differences, see Figure 1.We assume a rectangular region in two dimensions (i.e., r = (x, y)) and use a spatial sampling set given by the finite square lattice L = {(i r , j r ) : i = 1, . . ., I, j = 1, . . ., J}, where Δ r is the spatial sampling interval.For simplicity, we assume identical sampling intervals in both coordinates, but using different sampling intervals for each coordinate is straightforward (Different sampling intervals influence the accuracy of the field approximation only but not the principal features of the decentralized estimator).For simplicity, we assume that there are R sensors whose locations form a subset R of the lattice L.
Similarly, for the second-order temporal derivative, we have Here, k is the discrete time index, and Δ t is the temporal sampling period.It is upper bounded by Δ r /c to ensure numerical stability.The right choice of Δ t is beyond the scope of our paper, so that we refer our reader to [16].

C Forward model
We introduce the auxiliary function q(x, y, t) = ∂ t p(x, y, t) and define the pressure vector p k = vec{P k } with [P k ] ij = p(iΔ r , jΔ r , kΔ t ).The source vector s k and the pressure derivative vector q k are defined similarly.Applying the FDM to (4) then leads to the following linear system of equations: The diagonal matrix F 11 results from the boundary condition (4d).Its diagonal elements are depends on the boundary condition (4f).Similarly, the sparse matrix F 12 stems from (4a) and is given by else.

D Source model
We assume that there are S sources whose positions form a subset S of the discretization lattice L, i.e., s[i, j, k] = s l=1 s 0 [k − k l ]δ(i − i l , j − j l ), where s 0 [k] is a known waveform, but the positions (i l , j l ) and activation times k l are unknown.These unknowns are captured via the integer variables n[i, j, k] that describe, for a lattice point (i, j), the time between the source occurrence and the current time instant k, i.e., for the lth source there is n[i l , j l , k] = max{k -k l , 0}.If there is no source at position (i, j), then n[i, j, k] = 0.
Clearly, the source life span satisfies the state transition equation where S k = {(i l , j l )|k ≥ k l } is the set of sources active at time k.Arranging the variables n[i, j, k] into a vector n k similarly to p k , q k , and s k , we obtain where the elements of δ S k are zero or one depending on whether a source is active at the corresponding position and at time instant k, i.e., Note that the state vector n k has at most S non-zero elements.Using the convention s 0 [0] = 0, the source vector s k in ( 5) is rewritten as thereby linking the state Equation 6 and the forward model ( 5).

E Noise model
So far, no process noise has been considered and specified.Since the source function depends on time and space, these are the only quantities that suffer from noise and are modeled in the following: The temporal noise models the perturbation of a source's life span by an additional term in (6), while this is not possible for the spatial perturbation.This is due to the fact that the position of sources is coded into the sub-vector n k by placing its elements.From a practical perspective, this is done by a time-dependent matrix D k which displaces the elements of a vector to other positions (jitter) according to the mapping between grid and sub-vector n k .Equation 6becomes Here, n' is a random integer perturbation, ⊙ is the Hadamard (element-wise) product, and the lth column of the displacement matrix D k is given by e l+d(l) , with the canonical column unit vector and a random integer jitter d(l) whose probability mass is concentrated about zero.
Because of linearity, ( 9) is rewritten as

F Augmented state-space model
We next combine the state-space model ( 5) with ( 8) and (10) to obtain an augmented state-space model for the extended state vector This gives the state transition equation and Note that non-linearity is inherent in (11).
To complete the state-space model, the measurement equation is introduced.Since the actual observations are given by noisy samples of the pressure field at the sensor positions (i l , j l ) ∈ R, the measurement equation is where v k denotes measurement noise and with e l denoting the lth unit vector.

III Bayesian estimation
Our aim is to perform sequential Bayesian estimation of the state vector n k that characterizes the source positions and activation times.n k is one of the state vectors x k in (11).The data y k is specified in (14).A PF approach [7], i.e., a Monte Carlo approach based on importance sampling, is pursued.This approach exploits that our state-space model ( 11)-( 14) is a hidden Markov model (cf. Figure 2), where (11) implies a state transition distribution f(x k |x k-1 ) and ( 14) leads to a measurement distribution (likelihood function) f(y k |x k ), which both are assumed known in the following.

A Particle filter
To perform Bayesian estimation (e.g., MAP or MMSE) of (part of) the state vector x k given the past observations Using the Bayesian theorem and the fact that y k+1 and y 1:k are statistically independent (due to the Markov chain assumption) given x k+1 , we have which is known as the update step.While the measurement PDF f(y k+1 |x k+1 ) in ( 15) is known, f (x k+1 |y 1:k ) needs to be computed via the so-called prediction step, Here, the transition PDF f (x k+1 |x k ) is known and f (x k | y 1:k ) has been computed in the previous time step.
Since the integral in ( 16) typically is infeasible, it is usually approximated using a Monte Carlo technique known as importance sampling.The approximate sequential computation of the posterior distribution f (x k | y 1:k ) based on importance sampling using the transition PDF f (x k | x k-1 ) as importance (or, proposal) distribution q(x k ) leads to the particle filter.Here, the desired PDFs are approximated in terms of particles, i.e., samples x [l] k and associated weights ω where L is the number of particles.The new samples for the subsequent time instant are generated using the proposal distribution where for the generation of each new particle x . Sampling from q(x k+1 ) can be achieved by generating a noise realization w k and invoking the state transition Equation 11, i.e., k can be computed from the particle x [l] k according to (13).The dependency of the matrices on k issues from spatial noise.
The unnormalized weight for each new particle is where f v (v k ) is the distribution of the measurement noise and we used the measurement Equation 14.For i.

i.d. Gaussian measurement noise with variance σ
Once all unnormalized weights have been obtained, the actual weights are computed via the normalization ω . Particle filters suffer from a general problem termed sample degeneracy, i.e., after sometime only few particles have non-negligible weights.This problem is circumvented using resampling [21].With sampling importance resampling (SIR), new samples are drawn from the distribution and all weights are identical, i.e., ω To obtain initial particles x [l] 0 , samples of the state vector are needed.S random realizations of source positions and activation times are generated according to the prior distributions.Then, we apply the noise-free version of the state-space model (11) where n [l] 0 and u [l] are determined by the realizations of the source parameters (cf.(13) and Section II-D).The random variable k start denotes the time duration between source occurrence and activation of the estimator.

B Source localization
Using (17), the posterior PDF of n k (i.e., the last IJ elements of x k ) is approximated as (Note that n k contains all information about position and activation time of the sources.)The probability P{S k |y 1:k } for sources to be active at the coordinate set S k at time k is obtained via marginali- zation: Here, the function Q : ℝ IJ {0, 1} sets all entries of k to 1 which are unequal to 0. In the case of one source and a SIR PF with w [l] k = 1/L, the probability for a source at position (i, j) at time k is approximately obtained as where L i,j,k is the number of particles for which k ] i+(j−1)I > 0.

IV Decentralized scheme
The particle filter developed in the previous section is centralized in nature since it requires all pressure measurements and the observation modalities described by the globally assembled likelihood function and operates on the full state vector x k in a fusion center.Additionally, the computed estimates are inherently unknown on the individual sensor nodes.In a SN context, such constraints are undesirable since they imply a large communication overhead to collect the measured data, a high computational effort due to the high-dimensional state vector, a feedback to the sensor nodes to spread the estimates, and a central knowledge of measurement noise.Therefore, a decentralized scheme that distributes the data collection and computational costs among several clusters of sensor nodes is developed.This is achieved by splitting the state-space model ( 11), ( 14) into lower-dimensional sub-models (each corresponding to a cluster), cf. with [22,23].Due to the sparsity of the state-space matrices F and Γ, these sub-models are only loosely coupled, thus a decentralized PF that requires little communication between the clusters can be developed.

A SN clusters and partitioned state-space model
We start with partitioning the region of interest Ω into M disjoint subregions Ω (m) .The sampling lattice corresponding to each subregion is given by with its boundary nodes ∂L (m) , see Figure 3.The sensors within each subregion form clusters, denoted by To each subregion, we associ- ate a subset of elements of the state vector x k given by where and the superscsript (m) refers to region m.Except for F 12 , all of the blocks in the state-space matrices F k and Γ k are diagonal or zero (cf.( 12)).Thus, there is no coupling between the sub-vectors p in the adjacent subregions (in fact, this coupling is limited to samples at the boundaries of the subregions).The same applies for the sub-vectors n (m) k due to the spatial noise.This gives L (1)   L (2)   ∂L (1)   N (1) ⊂ ∂L (2)    This coupling Equation 25 is only possible for the time-independent part of these matrices.However, for uncorrelated noise between clusters, the time-dependent part, i.e., D k , is calculated separately according to Section II-E on every cluster at each time step, see below.
The coupling terms between neighboring subregions are given by with and, analogously, with Here, N (m) is the set of subregions adjacent to Ω (m) , and (m,m ) 12 is obtained from F 12 by extracting the rows and columns corresponding to L (m) and L (m ) .The off- diagonals of F 12 are extremely sparsely populated; in fact, (26)  means that pressure samples at subregion boundaries are exchanged between neighboring clusters in order to compute the finite differences.
Boundary conditions do not play a role in the decomposition step as long as (i) they do not depend on adjacent neighbors and (ii) their numerical solution fits into (5).In the first situation, an additional term (m,m )

B Decentralized particle filter
For the decentralized PF, we need to distribute the sampling (particle generation) step and the weight computation step.Based on the local particles weights, each cluster can then compute posterior source probabilities in a similar manner as in Section III-B.
1) Particle Generation: Sub-particles x [l,m] k within cluster R (m) are generated according to (25), cf. also (18), is a randomly chosen previous particle and n' , respectively.In order to compute the latter, only elements of x [l,m ] k that correspond to pressure samples from the boundaries of adjacent subregions are exchanged, and in the event of source hopping from one to another cluster, a message is sent.
2) Weights: Assuming independent measurement noise in the individual subregions, i.e., k ) , the weight update ( 19) is computed in each cluster as where the partial weights are computed within each cluster and then are shared among all clusters to obtain the final unnormalized weight [24] and [25] are treating the issue of computation of the global factorizable likelihood by means of distributed protocols.If these take longer than the time span between two estimator iterations, the particle filter converts to a particle predictor.
3) (Re)sampling: A remaining problem with the decentralized PF is that the sampling (particle generation) step (30) requires that the clusters pick local particles .The same problem occurs for the resampling procedure.Since a central random number generator whose output is distributed to each cluster incurs a large communication overhead, we propose to use identical pseudo-number generators in all clusters and initialize those with the same seed, thereby ensuring that all clusters perform the same (re)sampling (cf.with [24] and [26]).

V Decentralized source localization
The PF yields the posterior PDF of the sources' position and life span.To obtain the current MAP position estimates the maximum and the maximizing state of the posterior PDF Ps(i, j, k) in ( 23) must be found.In the decentralized scheme, each cluster disposes only of the local posterior PDF for the state sub-vector x (m) k .To find the global maximizing state, each cluster determines the local maximizing state and afterward the clusters use a distributed consensus protocol to determine the global maximum.For simplicity, this procedure is here developed for one source.
For the centralized PF, the posterior probability for a source to be active at time k at position (i, j) is given by (23).In the decentralized case, each cluster determines a similar probability according to ] i+(j−1)I > 0. Since the probabilities P (m) s (i, j, k) have disjoint support, the maximization underlying the MAP estimates (32) is While the local maxima with regard to L (m) can be determined within each cluster, the global maximization with regard to m requires communication between the clusters.Since sharing the local maxima among all clusters via broadcast transmissions requires a large coordinated transmission, we compute the global maximum via the maximum consensus (MC) algorithm [14].For the MC algorithm, we assume that only neighboring clusters communicate with each other.Thus, each cluster sends to the adjacent clusters a message which contains the local maximum and the position for which the local maximum is achieved.In the subsequent steps, each cluster compares the incoming "maximum" messages with their current estimate of the global position and retain the most likely and its associated position.In the next iteration, this message will be sent to the neighboring clusters.
Denote the current estimate of the maximum P k,max for cluster m by P(m)  k,max and let ( î(m) k , ˆj(m) k ) be the associated position estimate (initially, P(m)  k,max = P (m) k,max ).In our MC algorithm, termed argumentum-maximi consensus (AMC), at time instant k, each cluster performs the following steps: 1) Send a message containing the estimates P(m) k,max and ( î(m) k , ˆj(m) k ) to the neighbor clusters N (m) .2) Receive corresponding messages from the neighbor cluster, if a neighbor m ∈ N (m) remains silent, then P(m ) k,max = P(m ) k−1,max .3) Update the maximum probability and position as k+1,max = P(m) k,max to go 1), otherwise go to 2).When the maximum is fixed, all clusters converge to the true maximum after some iterations (depending on the diameter of the cluster communication graph).Here, the position of the maximum moves as the distributed PF evolves and the AMC will then allow the clusters to jointly track the maximum.

A Dimensions and trade-offs
Since we are estimating the 2-D position and activation time for each of the S sources, the number of unknowns equals 3S.This is relevant for the choice of the number of particles, cf.[4].For the calculation of the forward model (state transition), however, the dimension of the state vector x k is relevant which equals 3IJ.In the decentralized case, the computational complexity of the forward model is distributed across all clusters.
We now face the behavior of a high number of clusters.Generally, the volume of a polytope (cluster) L (m)  with edge lengths i .Generally, the dimension per cluster of the equation system to be calculated is 3|L (m) | which, in comparison, equals in the centralized case 3|L|.
In our 2-D problem, let the lattice L be partitioned into M = M i M j clusters of same size, M i clusters in idirection and M j clusters in j-direction.Then, e 1 = I/M i and e 2 = J/M j .Furthermore, the volume ∞, then the dimension of the equation system, which specifies the amount of computation, becomes in O(1/M) [27].Thus the compu- tational effort per cluster decreases when the number of clusters increases.On the other hand, an increasing number of clusters leads to a larger number of boundaries and hence to a larger communication overhead (i.e., message exchange between adjacent clusters).
Algorithm 1: Global initialization generate priors X 0 ; / / E q u a t i o n (20) decompose X 0 to {X (m) 0 } ; // Equation (24) choose seed s 0 (Section IV-B3); for m = 1 to M parallel do DD-SIR-PF( X (m) 0 , s 0 ) of cluster m; Algorithm 2: DD-SIR-PF(): Decentralized distributed SIR particle filter of cluster m input : X (m) 0 , s 0 k 1; wait while no signal sensed and no wake-up call; send wake-up call to other clusters; while estimating do observe:

B Communication between clusters
The variables that are broadcast by cluster m are summarized by the set The first subset , signifies a message about sources which migrate across boundaries from one cluster to another.Every message includes the new location and the current time duration since the occurrence of the sources.The last two terms stem from the AMC algorithm where Ŝ(m . Note that the cardinality of (34) which is a measure of the amount of transmission per cluster is given by the sum Here, the μ (i,m) k messages are disregarded.The amount of transmission in the decentralized case to adjacent neighbors for Note that there is no approximation compared to the centralized method and thus neither source coding nor approximations reducing the weight communication have been considered.For the communication of the weights, either the graph needs to be fully connected or the clusters need to act as relay.A summary is drawn in Table 1.

C Algorithm
The algorithm of the decentralized and distributed SIR PF together with the AMC is drawn in Algorithms 1-4.Compare it with that one in [28] and note that the forloop can be parallelized.
The joint setup of the computational nodes is shown in Algorithm 1 which consists of the calculation of the priors and the synchronization of the pseudo-random generator.Subsequently, each individual PF is launched (Algorithm 2).Two important sub-routines are plotted in their own tableaus: • Algorithm 3 calculates particles and sends messages when a source jumps over to another cluster.

All
*Source migration denotes the information that a source changes from one cluster to another.
• Algorithm 4 adds states from the neighbor clusters according to (25) and calculates the overall weight (31).
Algorithm 3: SI(): sample importance part Input: Algorithm 4: modify(): contribution of the neighbors.T (m) is a mapping from neighbors' pressure substates to the own sub-states with T (m) / / E u ation (31) normalize Ŵk ;

VII Simulations
In this section, we present simulations illustrating the performance of the proposed Algorithms 1-4.The configuration used in the simulations is shown in Figure 4 with parameters in Table 2 (N {μ, σ 2 } denotes the Gaussian distribution with mean μ and variance s 2 ).In particular, we used M = 5 subregions Ω (m) corresponding to 5 clusters each with 2 sensors.We considered a single source located in Ω (3) at the lattice point (i 0 , j 0 ) = (25,25); it is modeled by choosing the source function as s 0 [n] = s 0 (nΔ t ) where s 0 (t) is a time-shifted Ricker wavelet.A Ricker wavelet [29] is defined by the negative second derivative of a Gaussian function such that Here, ν is approximately the peak frequency.A Ricker wavelet shifted by 16.7 ms with ν = 60 Hz is used, i.e. s 0 (t) = ricker(t -16.7 ms), see Figure 5.The acoustic pressure field is simulated using the FDM introduced in Section II.A snapshot of the field at time k = 160 is shown in Figure 6.
The parameters used in the decentralized PF are summarized in Table 3 (U{a, b} represents a discrete uniform PDF with support [a, b]).For the fixed source position, we used a discrete uniform distribution on the 50 × 50 lattice.The spatio-temporal noise and the observation noise are drawn from a Gaussian distribution.The PF is initialized at time k = 0, and the source is assumed to become active at time instant k < 0. The maximum value of the random variable k start is a prior and is proportional to the maximal possible time duration between source arise and first detection (cf.(20)).Larger values of k start necessitate a larger number of particles to cover the time interval [-k start , 0] and thus to achieve the same approximation accuracy.

A Estimation of posterior PDF
For the centralized PF, Figure 7a shows an example of the posterior PDF P s (i, j, k) for the source position obtained with the centralized particle filter at time instant k = 160 (cf.( 23)).For comparison, Figure 7b shows the result obtained with the decentralized PF, i.e., the composition 5 m=1 P (m) s (i, j, k) of the local posterior PDF obtained by each cluster.It is seen that the centralized and the decentralized PF obtain similar results, and both yield a posterior PDF which is well concentrated about the true position (i 0 , j 0 ) = (25, 25) of the source.
Figure 8a, b shows the MAP and MMSE of the source's i coordinate and j coordinate, respectively.The   (25,25).Hence, in this specific case, the MMSE estimates outperform the MAP estimates for small k.After a certain number of PF iterations (around k > 6), however, the MAP estimates match the true source position better than the MMSE estimates.The variance of P s (i, j, k) for any given k (which can be interpreted as MMSE) is shown in Figure 9 and corroborates that for small-to-medium k, the i coordinate estimate is more reliable; this can be attributed to the specific sensor arrangements which favors better i-resolution (cf. Figure 4).

B Decentralized MAP source localization
This subsection illustrates the decentralized source localization using the AMC algorithm proposed in Section V (simulation setup unchanged).Recall that with AMC, each cluster has estimates P(m) k,max of the MAP probability and ( î(m) k , ˆj(m) k ) of the associated position.Figure 10a shows the local MAP probabilities P (m) k,max (cf.(33)) for all five clusters; clearly, only the third cluster builds up a distinguished maximum over time, which indicates that the source is located within Ω (3) .
All clusters track the global MAP probability, Figures 10b and 11, and eventually agree on the source position provided by cluster 3 whose behavior over time resembles the global estimates using the centralized PF (cf. Figure 8a, b).
After about 6 iterations, the PF achieves a localization accuracy on the order of the lattice spacing Δ r .These estimates could be further improved (with higher computational complexity) by refining the discretization lattice and increasing the number of particles.

VIII Conclusions
We proposed a scheme for the localization of multiple acoustic sources in a sensor network (SN).The method uses an augmented non-linear non-Gaussian state-space      obtained by the individual clusters using the AMC algorithm.In the latter case, the maximum time difference between two lines equals two.

Figure 1
Figure1The FDM model showing the discretization lattice, boundaries, sources, and sensors.L is the discretization lattice, while Ω denotes the area.

Figure 2
Figure 2 Hidden Markov model representation of the state-space model.
similarly for the sub-vector q (m) k .Coupling between state vectors from different regions, induced by the non-diagonal structure of F 12 , is between the sub-vectors q (m) k in one subregion and the sub-vectors p (m) k
contains only few non-zero terms corresponding to adjacent pressure samples and the change of sources from one to another cluster.D (m,m ) k is generated from every cluster m' such that the composition of all submatrices D (m) k and D (m,m ) k equals D k .From a practical perspective, elements of D (m) kare calculated separately on every cluster by means of spatial noise with additional triggering of a message to neighbor clusters whenever a source hop (migration) from one cluster to another is detected (this takes over the purpose of D

x
[l,m] k m = 1, ..., M, that correspond to the same global particle x [l] k .This choice is made at random according to the weights ω [l] k m) k collects all pressure sub-state particles on the boundary.The third, μ (i,m) k

Figure 4
Figure 4 Simulation setup comprising sensors, a single source, and SN cluster structure.

Figure 5
Figure 5 Ricker wavelet shifted by 16.7 ms with ν = 60 Hz (a) in the time domain and (b) its Fourier transform.

Figure 6
Figure 6 Pressure field from finite difference modeling after 41,750 time steps (corresponding to estimation time k = 160).

Figure 7 2 jFigure 9
Figure 7 Posterior source position PDF P s (i, j, k) at time k = 160 obtained with (a) centralized and (b) decentralized PF.

Figure 8
Figure8MAP and MMSE estimate of the i and j coordinate of the source (note that the lines of the centralized and decentralized MMSE estimations are close together).

Figure 10
Figure 10 Local (a) MAP probabilities P (m) k,max in comparison with (b) MAP probability estimates P(m) k,max

Figure 11
Figure 11 Source coordinate estimates of the individual clusters.

Table 1
Necessary message exchange

Table 2
Parameters for simulated hallway

Table 3
PF parameters