Source localization and tracking in a dispersive medium using wireless sensor network

In this paper, we address the issue of collaborative information processing for diffusive source localization and tracking using wireless sensor networks capable of sensing in dispersive medium/environment. We first determine the space-time concentration distribution of the dispersion from physical modeling and mathematical formulations of an underwater oil spill scenario, considering the effect of laminar water velocity as an external force. For static diffusive source localization, we propose two parametric estimation techniques based on maximum-likelihood (ML) and best linear unbiased estimator for the special case of our physical dispersion model. We prove the consistency and asymptotic normality of the obtained ML solution when the number of sensor nodes and samples approach infinity, and derive the Cramér-Rao lower bound on its performance. We also propose a particle filter-based target tracking scheme for moving diffusive source and derive the posterior Cramér-Rao lower bound for the moving source state estimates as a theoretical performance bound. The performance of the proposed schemes are shown through numerical simulations and compared with the derived theoretical bounds.


Introduction
The release of liquid petroleum hydrocarbon into the ocean or coastal water due to human activity has attracted tremendous attention because of its environmental, biological, and economical impact. Recent BP oil disaster in the Gulf of Mexico is a perfect example of how spill stemmed from a sea-floor oil gusher can severely damage the marine and wildlife habitats as well as the Gulf 's fishing and tourism industries. Research in modeling and predicting the extent of such oil spill can assist in planning and emergency decision-making, thereby reducing the threats and hazardous effects on the environment as well as the economic cost. Considering the fact that this is a diffusive source estimation and tracking problem, such research can in general be applicable in many other similar contexts such as homeland security, environmental and industrial monitoring, pollution control, servers, and data center temperature monitoring as well [1][2][3][4][5][6][7][8]. For example, http://asp.eurasipjournals.com/content/2013/1/147 make them suitable candidates for the set of applications involving monitoring of diffusion phenomena that we are interested in.
Source or target localization using distributed sensor arrays is an area of active research interest for a considerable period of time [19,20]. In the past, detection and localization problems of diffusive sources in WSN have been a topic of interest, specially in the case of chemical/biological threat detection. Interesting research in this context can be found in [3,4,9,10], where biochemical concentration distribution in space and time for different types of diffusive sources, diffusion models, and/or sensor networks is estimated. For instance, remotely localizing a gas or odor source using mobile robot was proposed in [3] by fitting the gas distribution model to the gas sensor response at the sensor locations. However, the mobile sensor dynamics model therein was obtained empirically, which does not allow for dynamic environment and moving diffusive source. In [4], a maximum-likelihood (ML) estimator was developed for localizing vapor-emitting sources, and its asymptotic normality of the obtained ML estimator was proved when the signal-to-noise-ratio (SNR) approaches infinity. Many other estimation techniques have also been used in diffusive source parameters estimation literature [9,10,[21][22][23]. In particular, Bayesian estimation has been applied in [9,21] in a sequential manner, which is not suitable in many practical scenarios where faster estimation and immediate actions based on the estimation are top priorities. A real-time maximumlikelihood estimation method was proposed in [23] for estimating diffusive source parameters, where consistency and asymptotic efficiency of the obtained estimator were proved when the density of sensors becomes infinite. In [24], the problem of impulsive diffusive source localization was solved assuming the spatial sensor measurements at any sensor location as a scaled and shifted version of a common prototype function, leading to solving a set of linear equations. However, the physical diffusion models used in [23,24] are oversimplified with the diffusive sources assumed to be impulsive or instantaneous in nature.
Although research has been done in tracking and/or estimating time-varying parameter estimation in general [25][26][27][28], to the best of our knowledge, very few attempts have been made in time-varying diffusive source parameter estimation. Some of these methods cannot be applied directly into our time-varying parameter estimation model since, e.g., for a moving source, the concentration at the current time is affected by all past values of source position. Therefore, time-cumulation effects on the concentrations (i.e., observations) must be taken into account to estimate time-varying parameters. Among previous works, a parametric moving path model for a diffusive moving source was discussed in [10], where the moving source path was approximated using finite number basis functions. Tracking performance in this case depends on the smoothness of the source trajectory, prior information about the moving source trajectory, and choosing a suitable finite set of basis functions. In [29], a novel recursive algorithm was proposed to track the intensity of a diffusive point source, but the source location was considered as an unknown static value.
The aforementioned limitations may be overcome by developing or exploiting state-of-the-art Bayesian-based location tracking methods suitable for handling highly nonlinear diffusion processes. In the Bayesian approach, the key is to construct the posterior probability density function (PDF) of the underlying state vector based on all available information. For linear and Gaussian state dynamics and observation models, the optimal minimum mean squared error (MMSE) solution is tractable and is given by the well-known Kalman filter [30]. However, for most of the real-world scenarios, dynamic state estimation problems are nonlinear and non-Gaussian, and obtaining optimal closed-form solution is not tractable under the Bayesian approach. In these cases, suboptimal approaches such as extended Kalman filter and Gaussian-sum filter [31] are used with certain approximations. These suboptimal algorithms become inefficient for highly nonlinear and non-Gaussian systems. In these cases, numerical techniques based on sequential Monte Carlo methods are used to achieve better performance for highly nonlinear systems. To that end, the idea of particle filtering was introduced in [32] as an effective method of representing PDF in terms of a set of random sampling.
In this paper, our main objectives are to efficiently estimate and track diffusive source location using a wireless network of chemical sensor capable of sensing in diffusive environment. To cater to the objectives, we formulate and derive a physical model for the space-time substance dispersion mechanism of an underwater diffusive source. The modeling and the proposed solution methods can also be extended to other important diffusion phenomena involving biochemical contaminant materials as well. We propose and implement ML and best linear unbiased estimator (BLUE)-based parameter estimation techniques for a static diffusive source which is continuously emitting substance [33]. In the previous literature, such as in [4], the asymptotic normality of the obtained ML estimator was proved when the SNR approaches infinity. We prove both the consistency and asymptotical normality of our obtained ML-based solution when the number of sensor nodes and time samples go to infinity, thus allowing for the option of tweaking these two parameters. We derive the Cramér-Rao lower bound (CRLB) as a theoretical performance bound for a special case of our obtained physical dispersion model. We also propose a particle filter (PF)based target tracking method for moving diffusive source. http://asp.eurasipjournals.com/content/2013/1/147 To the best of our knowledge, moving diffusive source tracking using particle filtering approach has not been attempted before. The posterior Cramér-Rao lower bound (PCRLB) for the moving source state estimates is also derived as a theoretical performance bound [34].
The remainder of this paper is organized as follows: Sections 2 and 3 discuss, respectively, modeling of an underwater oil spill scenario and measurement model for static diffusive source localization using sensor network. The proposed statistical methods for static diffusive source localization and corresponding theoretical performance bound are discussed in Section 4. Section 5 presents the proposed particle filter-based method for moving diffusive source tracking with theoretical performance bound analysis in detail. Section 6 shows the validity and effectiveness of our proposed methods for diffusive source localization and tracking through numerical simulations. Finally, Section 7 concludes the paper by summarizing our results.

Physical model for dispersion
We first derive the physical models for the space-time substance dispersion mechanisms from a diffusive source and then transform the obtained dispersion model to a statistical measurement model. The transport model of a substance from a diffusive source can be obtained by solving the corresponding diffusion equation. Diffusion equation describes the dispersion of particles from a region of high concentration to regions of lower concentration due to random molecular motion. Let us denote the concentration of the diffused substance at a position r = [x, y, z] T and at time t as c(r, t). Ignoring the effects of external forces for a source-free volume and for space-invariant diffusivity constant κ, the concentration of a dispersed substance follows the following diffusion equation [35]: To solve the above differential equation, appropriate boundary and initial conditions are required. We first compute the concentration for a stationary impulse point source of unit mass to obtain Green's function. The obtained result is then extended for a continuous source by integrating the source-release rate with the Green's function. Denoting the Green's function of the impulse source as c G (r, t), the concentration of a continuous point source with mass release rate μ(t) and initial release time t I can then be given by the following integral: For parametric estimation case, it is to be noted that from the concentration measurements taken by the sensors, we can first estimate the source parameters of interest and then predict its cloud evolution in space and time by inserting the estimated parameters into (1).
Although the main focus of this paper are diffusive source localization and tracking, we introduce a special diffusion phenomenon, i.e., an underwater oil spill, to demonstrate how to model and solve for a practical diffusion phenomenon and also to motivate the practical importance of the problem we are discussing. As shown in Figure 1, we model an underwater oil spill as a diffusion occurring in a two-layer semi-infinite medium (i.e., water and air). We assume that the oil spilling source is located at the bottom (i.e., river/sea bed) at a location r 0 = [x 0 , y 0 , z 0 ] T . The depth of water level is 0 ≤ z ≤ L with diffusivity κ w and concentration c w . The same quantities for air (z > L) are denoted as κ a and c a respectively. Along the z-axis, we need to solve the following differential equations: Considering only point impulse source located at z = z 0 , where 0 ≤ z 0 ≤ L and impermeable boundary at z = 0, we have the following initial condition: Boundary condition (2) implies the continuity of concentration at the interface z = L. Condition (3) represents the fact that there is no accumulation of diffusing substance at z = L. Finally, the third boundary condition in (4) reflects the assumption that the medium at z = 0 is impermeable. By applying the concept of Laplace transform on the above system of partial differential equations, we can obtain the solution to the spatio-http://asp.eurasipjournals.com/content/2013/1/147 temporal concentration distribution (omitting the details in [10,35,36]): As can be seen, the concentration curve can be considered to be the superimposed curve resulting from each successive reflection (from the surface layer) being superimposed on the original curve. In practice, if κ w κ a , then ρ → 1. Therefore we have, Considering the laminar water velocity working along the X-Y -plane as an external force, we have v = [v x , v y , 0] T . The diffusion equations along the x and y axes will include additional advection term [35]: For X-Y -plane, there is no boundary condition and the initial condition is given as Using the concept of Fourier transform for solving partial differential equations, we can solve for the following concentration distribution along x and y axes [37]: Based on our assumptions on initial and boundary conditions and for rectangular parallelepiped space, the Green's function solution for 3-spatial-variable case is the product of the solutions of the three single spatial-variable cases with stationary impulse point source [10,35]. Therefore, the Green's function c G (r, t) for the space-time http://asp.eurasipjournals.com/content/2013/1/147 concentration distribution can be obtained as the product of the solutions in (5), (6), and (7): Considering the source mass release rate to be constant μ(t) = μ, the final solution for concentration of oil diffusion in water for stationary continuous source with mass rate of μ(t) can be obtained from (1): where Derivation to (9) is given in Appendix 1. For the sake of simplicity from here on, we denote the diffusivity constant κ w = κ.

Moving diffusive source
For a moving diffusive source-emitting substance continuously in a semi-infinite medium similar to our case, space-time concentration distribution can be obtained using the concept of convolution integral from the Green's function solution corresponding to stationary impulsive source. In this case, substance concentration at any time instant is affected by all the past values of source position and release rate. Therefore, time-cumulation effect on the concentrations has to be considered to obtain complete physical model. For a moving diffusive source continuously releasing substance at a mass rate μ(t), the space-time concentration distribution in a semi-infinite medium can be obtained for a given Green's function c G (r, t) using the following integral: where T represents the source moving path. The advantage of solving the physical diffusion model corresponding to a moving diffusive source using (10) is that the initial, boundary, and other necessary conditions can be taken into account to solve for the stationary case in the first step before extending it to the moving source case.

Measurement and system models for static diffusive source localization
We consider a WSN consisting of a fusion center (FC) and N spatially distributed biochemical static sensor nodes capable of sensing in dispersive environment. For practical consideration, we assume that the N distributed sensors are located in a rectangular volume in space such that It is also assumed that the source-to-sensor distances are much higher than the source and sensor dimensions. Each sensor node takes measurements at times t k ; ∀k ∈ {1, 2, . . . , T}, where T is the total number of time samples. Assuming that the physical model discussed before is the underlying dispersion mechanism, we may obtain a measurement model for a sensor at a position r j and at time t k as y(r j , is the concentration of interest, b is a bias term, and e(r j , t k ) ∼ N (0, σ 2 ) is the sensor noise assumed to be independent in both time and space. For the sake of brevity, it can be rewritten in the simplified form as where y j,k = y(r j , t k ), e j,k = e(r j , t k ), c j,k (θ) = c(r j , t k ), θ ∈ R n×1 is the unknown source and medium parameter vector that we are interested to estimate, and b is the bias or clutter term representing the sensor's response to foreign substances that may be present in a diffusive field of interest. The bias term is assumed to be space and timeinvariant such that the foreign substances interfering with the actual measurements are in steady state. If we want to localize a static diffusive source, then only [x 0 , y 0 , z 0 ] are the parameters of interest. It is to be noted that some of the parameters, such as the diffusivity constant κ, bias term b, and noise variance σ 2 , can be measured at the calibration stage, thereby reducing the cost of computation during the detection/estimation phase. We assume that the sensor nodes are in sleep mode until they are activated by some central control (i.e., FC) due to a possible release of a substance from a diffusive source. The activated sensor nodes take measurements of the substance's concentration at time instants t k s and then return to sleep mode. For N number of nodes in a WSN and with each node taking T number of time samples of the substance concentrations at their respective locations, let y ∈ R NT×1 be the measurement vector received at the FC.

Static diffusive source localization
In this section, we use the maximum-likelihood estimator (MLE) and the BLUE to estimate the location of an underwater diffusive source diffusing oil into water. For simplicity of exposition, we consider a special case of our obtained physical model when an oil spill occurs in an infinite (L → ∞) underwater medium. In this case, the concentration at any position r j at time t k is reduced to the following expression [4,35]: where erfc(.) is the complementary error function.

Maximum-likelihood-based source localization
From the measurement model discussed in Section 3, the conditional PDF of the measurements taken by the jth node at time t k is p(y j,k |θ ) ∼ N (c j,k (θ )+b, σ 2 ). Hence, the log-likelihood function formed at the FC can be written as The log-likelihood equations are obtained by for u = 1, 2, 3, where θ u is the uth element of θ , and Since the system of equations in (14) is nonlinear, there is no closed-form solution to it. We can obtain an ML estimation of the source location using any suitable nonlinear optimization technique. In this case, (14) is solved using simplex search algorithm [38].
The CRLB provides a lower limit on the mean squared estimation error of an unbiased estimator of nonrandom parameter [30]. CRLB in this case can be obtained as CRLB ≥ I −1 θ , where I θ ∈ R 3×3 is the Fisher information matrix (FIM) formed at the FC. The u-vth element of the FIM can be found as where (16) was obtained using the independence assumption of observations in space and in time.
A sequence of estimatorsθ n to an unknown parameter vector θ is said to be consistent if the sequence converges in probability to θ, i.e., lim n→∞θ n = θ , where n is the sample size [30]. It is desirable to have a consistent MLE as consistency ensures that for large data sets, the MLE will converge to the true parameter. The obtained MLE to our source localization problem is consistent when the Once consistency for the obtained MLE is established, the next important thing is to check the asymptotic normality. An asymptotically normal estimator is a consistent estimator whose distribution around the true parameter θ approaches a normal distribution with standard deviation shrinking in proportion to 1/ √ n as the sample size n grows, i.e., √ nI θ θ n − θ −→ N 0, I −1 , where I θ and I http://asp.eurasipjournals.com/content/2013/1/147 are the Fisher information and identity matrices, respectively [30]. It ensures that the estimator not only converges to the unknown parameter, but it converges fast enough at a rate 1/ √ n. We address this issue with the following theorem on asymptotic normality. is an open subset of , the following is true: in distribution where the (u, v)th element of the matrixĪ θ is given by Proof. See Appendix 3.

Best linear unbiased estimator-based source localization
The advantages of using the BLUE for static diffusive source localization are that there are no constraints on the PDF and also knowing only the mean and covariance of the measurements is enough. However, observations have to be linear for the BLUE algorithm. In this subsection, we assume that the distributed sensing nodes are capable of estimating their respective distances from the source using BLUE.
Since the complementary error function can be approximated as erfc(z) ≈ 1 − 2 √ π z, our observation model for jth node at the kth time instant can be linearized in terms of the inverse of the source-to-node distances from (11) and (12): . Since all the parameters are known except for the diffusive source location, we can writeỹ j,k = y j,k − a k = hd inv j +e j,k . Therefore, the observation vector formed at the jth node can be written as where h is a column vector of all hs and e j = [e i1 , e i2 , . . . , e iT ] T . Since e j,k ∼ N 0, σ 2 for ∀j, k and measurement noise is assumed to be independent and identically distributed across space and time, the covariance matrix ofỹ j is˜ j = diag σ 2 , σ 2 , . . . , σ 2 ∈ R T×T . The optimal BLUE estimator formed at jth node is given bŷ After the distributed nodes estimate their respective distancesd j = |r j − r 0 | from the source using BLUE, all nodes sendd j s to the FC for further processing. It is to be noted that the source-to-node distance estimation can also be performed at the FC. The signal received at the FC from the jth node can be expressed as f j =d j + w j , where w j is normally distributed with mean 0 and variance σ 2 m . For N number of nodes, the data vector available at the FC can be written as (20) To solve for the source location from (20), simplex search algorithm [38] has been used.

State dynamics model
For the simplicity of exposition and computation, we consider the problem of tracking a diffusive source moving in a 2D X-Y -plane. The assumption can be easily extended to the 3D case without any loss of generality. Let us denote by s k = x s,k y s,kẋs,kẏs,k T , the state vector associated with the moving source at time t k , where the first two elements represent the source position in 2D and the next two elements represent the speed of the moving source, respectively. We assume linear dynamic model for the source state vector: for k = 1, 2, . . ., with the initial known distribution p(s 0 ) for s k , where F is a 4 × 4 matrix that models the state kinematics [39]: where T s is the time difference between two consecutive measurements. The noise vector u k is assumed to be 0 mean Gaussian with covariance matrix Q [39]: which models the acceleration terms in the spatial directions, and σ 2 u is the variance of the process noise.

Observation model
In case of a moving diffusive source continuously emitting diffusing substance in 2D, we may obtain a measurement model for a sensor at a position r j,k and at time t k as where z j,k is the jth node's observation at time t k ; c j,k c(r j,k , t k ) = t k t I μ(τ )c G r j,k − r s (τ ), t − τ is the substance concentration at jth node location at time t k ; moving diffusive source location at time t k is r s,k =s k = x s,k , y s,k T ; location of jth node at time t k is r j,k = x j,k , y j,k T ; and ν j,k ∼ N 0, σ 2 ν is the sensor measurement noise assumed to be independent in both time and space. Note that for static sensor node locations, we use r j,k = r j = x j , y j T by dropping the time index, since node locations do not change over time. By assuming the additive white Gaussian noise channel for the sake of simplicity, the received signal at the FC from the jth node at time t k can be written as where j,k is the received noise which is assumed to be Gaussian with mean 0, variance σ 2 and e j,k = j,k + ν j,k and σ 2 = σ 2 ν + σ 2 . We denote y j,1:k as the measurement vector from jth node up to the time t k , and y c,1:k {y 1,1:k , y 2,1:k , . . . , y N,1:k } T as the collection of all measurements at the FC from N-distributed sensor nodes.
In a realistic moving source scenario, the instantaneous velocity is restricted by some practical upper limit. Hence, for lower sampling time T s , we can assume that the moving diffusive source moves in a linear fashion between two observations with an average velocity determined by the source locations r s,k and r s,k+1 . For 2D moving diffusive source tracking with no external force in action, the Green's function can be obtained from (6) and (7) as Therefore, for a continuous moving diffusive source with constant mass rate μ(t) = μ, observations taken by the jth node at kth time instant can be written as where

Target tracking using particle filters
In Bayesian belief update, to estimate state vector s k at time instant k, we need to construct posterior distribution p s k |y c,1:k with initial PDF p(s 0 ). The Bayesian belief update is done in two stages: prediction and update. Prediction. Considering that p s k−1 |y c,1:k−1 is available at time k, the PDF p s k |y c,1:k−1 can be obtained as [40] p s k |y c,1:k−1 = p s k |s k−1 p s k−1 |y c,1:k−1 ds k−1 .
Update. If observations y c,1:k are available at time instant k, the posterior distribution to estimate the state vector s k is given by [40] p s k |y c,1:k = p y c,k |s k p s k |y c,1:k−1 p y c,k |y c,1:k−1 .
Since the observation model is highly nonlinear, the analytical solution for the optimal estimator is not tractable in our case. Hence, we use sequential Monte Carlo method to approximate the posterior PDF (27) with particle filters [32]. Let us denote X k = s i k , w i k P i=1 to be the random measure that characterizes the posterior PDF p s k |y c,1:k , where P is the number of particles. Then p s k |y c,1: is the Dirac delta function. The state vector estimate at time t k can be obtained as http://asp.eurasipjournals.com/content/2013/1/147 predicted stateŝ k+1|k and the corresponding covariance matrix U k+1|k can be obtained from the state dynamics in (21) asŝ k+1|k = Fŝ k|k and U k+1|k = FU k|k F T + Q.

Posterior Cramer-Rao lower bound analysis
Analogous to the CRLB, the PCRLB provides a lower bound for the mean squared error of random parameter estimation [34]. Let us define the joint probability distribution of S k and y c,1:k for an arbitrary k is p S k , y c,1:k = p k , where y c,1:k is the observation vector formed at the FC at kth time instant and S k = (s 0 , s 1 , . . . , s k ). Following (26), the concentration at any time k + 1 for any node j can be written as c(r j , t k+1 ) c j,k+1 = ζ j,0:1 +ζ j,1:2 +. . .+ζ j,k−1:k +ζ j,k:k+1 .
Based on the assumed observation model in (25), the log-likelihood function L k+1 = log p y c,k+1 |s k+1 , S k at (k + 1)th time instant formed at the FC is given by Let I(S k ) ∈ R 4k×4k be the information matrix derived from the joint distribution p k . We wish to solve for the information submatrix for estimating s k , denoted by I k .
The following theorem gives a two-step recipe for computing I k .
Theorem 3. The sequence {I k+1 } of the posterior information submatrices for estimating state vectors s k+1 can be computed as follows: and = ∇ ∇ T with ∇ being the Laplacian operator.
Proof. See Appendix 4. Note that the information submatrix computation in (28) involves computation of the inverse of a matrix of size 4k × 4k. This is because of the output y j,k+1 at the jth node at (k+1)th time instant being a function of all the previous states S k+1 .

Simulation results
In the following, we show the performances of our proposed models and schemes through numerical simulations.

Simulations for the physical model in Section 2
We show the space-time concentration distribution of a static continuous point source (oil spill source) located at the bottom of a sea at r 0 based on the physical diffusion model formulated in Section 2. The parameters used for this simulation are oil release rate μ = 10 3 kg/s, diffusivity constant of oil in saline water κ = 25 m 2 /s, initial release time t I = 0 s and laminar water velocity v = [50, 50, 0] m/s. The oil spill source is assumed to be located at r 0 = [0, 0, 0] T and the depth of water is taken to be L = 100 m from the sea bed. Figure 2 shows the spatial concentration distribution for two different time instants t = 1 and t = 100 s. It can be seen from Figure 2 that as the oil source is located at the origin, the concentration is high near the origin at t = 1 s. By the time it is 100 s, oil has diffused over larger distance from the source. It is interesting to see that since laminar water flow is assumed to be only active in the positive x and y directions, concentration increases more along the positive X-Y -plane with the increase in time.

Static diffusive source localization
Here, we show the simulation results in estimating the location of a static diffusive source using the proposed MLE and BLUE-based methods from the concentration observations taken by the sensing nodes. For the sake of simplicity, we consider a 2D diffusive field volume of = [−50, 50] ×[−50, 50] m 2 . We assume that the sensors are placed in a uniform 2D grid such that the distance between adjacent sensors along the same ordinate is approximately 14.3 m. Parameters used for simulations are number of nodes N = 64, r 0 = [0, 0] T , μ = 1, 000 kg/s, b = 10 −4 kg/m 2 , t I = 0 s, and  Figure 3 shows the normalized mean squared error (MSE) and CRLB (in dB) with the increase in the number of nodes and samples. The normalized MSE and CRLB are obtained by dividing each with the diffusive field volume. As one would expect, the estimation error decreases as more distributed nodes and samples are considered for estimation purpose. Since it is a 2D location estimation problem, we need at least three nodes to determine the source location correctly. It is interesting to note that the estimation performance is slightly better than the CRLB in some cases. This is due to the fact that the ML estimator in this case is biased (suggested from simulation), and thus can outperform the CRLB by trading variance for bias. In this particular case, the continuous diffusive source can be localized with a resolution of less than 12 cm.
The estimated source location using the BLUE estimator is shown in Figure 4 as a function of the number of nodes and time samples for different values of σ 2 m s. As one would expect, the overall performance obtained from the BLUE estimator is not as good as that from the MLE due to the linear approximation applied on the observation http://asp.eurasipjournals.com/content/2013/1/147 model in (17). However, the performance of the BLUE estimator-based localization improves as the number of nodes and/or time samples increases. This is because for N, T → ∞, the complementary error function in (12) tends to be equal to 1, causing the linearization to have almost no effect on the approximation.

Moving diffusive source tracking
In this subsection, we analyze the performance of our proposed moving diffusive source tracking scheme. We use the same sensor network setup as described in Section 6.2. The initial source state vector is assumed to be Gaussian with mean μ = [0, 0, 0, 0] T and covariance matrix 0 = diag [0.01, 0.01, 0.01, 0.01] T . The intensity of the state process noise is σ 2 u = 0.1. The sampling time is assumed to be T s = 0.5 s, and the total number of random realizations used for simulations is 50. The tracking is performed for 30 s and the number of particles in the PF is N p = 1, 000. The rest of the parameters is same as in Section 6.2. The performance measure is taken as the root mean squared error (RMSE) of the moving source position esti- The RMSE is compared with the square root of the PCRLB components of the position error, PCRLB k ≈ . Figures 5 and 6 show the tracking performances of the proposed tracking scheme using particle filter for gridbased and random node deployment strategies, respectively. It can be seen that the target trajectory can be tracked with better accuracy in Figure 6 compared to that in Figure 5. Figures 5B and 6B show the RMSEs on the tracking performances for the aforementioned two node deployment strategies respectively. The obtained RMSE with the random node deployment case is better and closer to the derived PCRLB than those of the grid-based node deployment case. This is because for a fixed node density, the expected nearest neighbor node distance (from the source) in case of random node deployment is less than the inter-node spacing in grid-based node deployment, which in our case is 14.3 m. The random node deployment is specially suitable when there is no pre-designed infrastructure for sensor network and also when the diffusive field is hazardous for human deployment.
It is of interest also to investigate the performance of the proposed target tracking method when the sampling time T s is varying. Figure 7 shows the effect of sampling time T s on the tracking performances of the proposed moving diffusive source tracking scheme using grid-based node deployment strategy, keeping all the other parameters same as mentioned before. As one would expect,   the tracking performance decreases with the increase of sampling time T s . This is because for higher values of T s , the process noise will increase according to (23). Since we are also assuming that the movement of the diffusive source is almost linear between two successive time instants, the lower T s will result in better accuracy of the proposed tracking scheme.

Conclusion
In this paper, we obtained spatio-temporal distribution of the substance concentration by solving physical diffusion model for an underwater oil spill scenario which considers laminar water velocity as an external force. The obtained mathematical model was found to be capable of modeling satisfactorily the underlying physical diffusion phenomenon. We proposed two parametric estimation methods based on MLE and BLUE for determining static diffusive source location using wireless sensor network. We also obtained the CRLB as theoretical performance bound for source localization. It was observed that though the MLE performs better than the BLUE-based diffusive source localization method, the latter shows satisfactory performance trend for large number of sensing nodes and time samples. We also proposed a particle filter-based target tracking method for moving diffusive source-emitting substance continuously into the dispersive medium. PCRLB corresponding to the moving diffusive source tracking was obtained as a theoretical performance measure and was compared with the simulation results. The effect of sampling time on the moving source tracking was also investigated. The performance of the proposed estimation and tracking methods are shown to be excellent using numerical simulations. In future research, we plan to combine our obtained analytical results with non-model-based numerical techniques to make them applicable for more realistic and complex scenarios. ∂x 0 are continuous functions of x 0 , by using Cauchy-Schwarz inequality, we can obtain ≤ TS with S being some positive real value. For practical consideration, assuming 0 ≤ |x j −x 0 | |r j −r 0 | 2 ≤ P, 0 < 1 |r j −r 0 | ≤ M 1 , and 0 < 1 |r j −r 0 | ≤ M 2 , K N,T x 0 ; x 0 can be written as If d N,T = N 3 T 3 > 0 for N ≥ 1, T ≥ 1, then we can claim that lim N,T→∞ For the proof of (33), assuming 0 ≤ |x j −x 0 | 2 |r j −r 0 | 2 ≤ Q 1 and T k=1 1 t k −t I < TQ 2 with Q 1 and Q 2 being some positive real numbers, we obtain the following from (15): ilarly, for y 0 and z 0 , we can also claim that the MLE to the diffusive source localization problem is consistent when the number of sensor nodes and time samples go to infinity.

Appendix 3 Proof of Theorem 2
To prove the asymptotic normality of the MLE, we define j,k y j,k |θ = log p y j,k |θ ,˙ j,k,u y j,k |θ = ∂ ∂θ u j,k y j,k |θ , and¨ j,k,u,v y j,k |θ = ∂ 2 ∂θ u ∂θ v j,k y j,k |θ . Below, we verify the necessary conditions mentioned in [42] for our obtained MLE to be asymptotically normal.
From practical point of view, there is no loss in generality in assuming that is an open subset of . Also because the obtained MLE to source localization is consistent, it is also consistent even when θ 0 ∈ • ⊂ . Thus conditions N1 and N2 are satisfied.
From the notations defined above, since ∂c j,k (θ) ∂θ u and ∂ 2 c j,k (θ) ∂θ u ∂θ v exist for u, v = 1, 2, 3, it can be easily verified that both˙ j,k,u y j,k |θ and j,k,u,v y j,k |θ exist almost surely. Therefore N3 is satisfied.
Since θ ∈ and¨ j,k y j,k |θ is a continuous mapping of θ, we can claim that¨ j,k y j,k |θ is indeed uniformly continuous on θ in j and k [41]. Also, because¨ j,k y j,k |θ : y j,k → R is a continuous function of y j,k with y j,k being Lebesgue measurable function,¨ j,k y j,k |θ is also a measurable function of y j,k and condition N4 is satisfied. To satisfy N5, it is easy to verify that E ˙ j,k,u y j,k |θ = 0 for all j, k and u. Since p(y j,k |θ ) ∼ N c j,k (θ) + b, σ 2 and p(y j,k |θ ) is continuous and Lebesgue measurable in y j,k , ∂ 2 ∂θ u ∂θ v p(y j,k |θ )dy j,k = ∂ 2 ∂θ u ∂θ v p(y j,k |θ )dy j,k is valid for all j, k, u, and v, and thus N6 is satisfied. To prove condition N9 since¨ j,k,u,v y j,k |θ is a uniformly continuous function of θ (shown in condition N4), for any > 0, there exists one δ > 0 such that ¨ j,k,u,v y j,k |θ −¨ j,k,u,v y j,k |θ 0 < δ, ∀ ||θ − θ 0 || < . is a finite real number. Therefore, the obtained MLE of the diffusive source location is asymptotically normal when the number of sensor nodes and time samples go to infinity.

Appendix 4 Proof of Theorem 3
For p (s 0 ) ∼ N (μ 0 , 0 ), the initial condition for the FIM is I(S 0 ) = E − s 0 s 0 log p (s 0 ) = −1 0 . Decomposing S 1 as S 1 = [s T 0 , s T 1 ] T , I(S 1 ) can be written as Since error is independent across space and time, using the concept from block matrix inversion, the information submatrix that provides the mean square error estimate of s 1 is given by where D 1 = E − s 1 s 1 L 1 +Q −1 , R 1 = E − s 0 s 0 L 1 +F T Q −1 F, and B 1 = E − s 0 s 1 L 1 −F T Q −1 . Similarly, decomposing S 2 as S 2 = [s T 0 , s T 1 , s T 2 ] T , the FIM I(S 2 ) can be written as follows: The information submatrix I 2 can be found as an inverse of the right-lower 4 × 4 submatrix of [I(S 2 )] −1 : By extending the above procedure and decomposing S k+1 = s T 0 , s T 1 , . . .
The information submatrix I k+1 can be generalized as an inverse of the right-lower 4 × 4 submatrix of I(S k+1 ) s k+1 L k+1 , and R k+1 is defined in (29