Distributed Least-Squares Estimation of a Remote Chemical Source via Convex Combination in Wireless Sensor Networks

This paper investigates the problem of locating a continuous chemical source using the concentration measurements provided by a wireless sensor network (WSN). Such a problem exists in various applications: eliminating explosives or drugs, detecting the leakage of noxious chemicals, etc. The limited power and bandwidth of WSNs have motivated collaborative in-network processing which is the focus of this paper. We propose a novel distributed least-squares estimation (DLSE) method to solve the chemical source localization (CSL) problem using a WSN. The DLSE method is realized by iteratively conducting convex combination of the locally estimated chemical source locations in a distributed manner. Performance assessments of our method are conducted using both simulations and real experiments. In the experiments, we propose a fitting method to identify both the release rate and the eddy diffusivity. The results show that the proposed DLSE method can overcome the negative interference of local minima and saddle points of the objective function, which would hinder the convergence of local search methods, especially in the case of locating a remote chemical source.

summand, we derive a set of probable source locations by substituting the noisy measurement to the concentration distribution model. An estimate of the source location is locally determined from each of the location sets by empirically maximizing the probability of being the real source location based on the prior information about the source location. Then, the problem of CSL using a WSN can be solved by comprehensively incorporating the information of these local estimates into a global estimate. We consider a convex combination of these local estimates as the global estimate of the source location. The convex combination of these local estimates, which are distributed on individual sensor nodes, is realized in a distributed manner using the distributed average consensus algorithm proposed in [14]. Finally, the above-mentioned process repeats iteratively until the termination condition is satisfied and the global estimate is considered as the prior information about the source location for determining the local estimates at the next iteration. The proposed method was assessed in both simulations and real experiments, including the case of locating a remote chemical source.
The rest of this paper is organized as follows: the chemical concentration distribution model, the sensor-source distance function, and the measurement model are described in Section 2. The DLSE method for CSL using WSN is proposed in Section 3. The results of numerical simulations and real experiments are presented in Sections 4 and 5, respectively. Some concluding remarks are given in Section 6.

Formulation and Preliminaries
In this section, the advection-diffusion chemical concentration distribution model is introduced, and then the distance between the unknown source and the sensor node is represented as a function of their relative angle and the theoretical concentration at the sensor node location. Finally, the concentration measurement model of the sensor nodes in the WSN is presented.

Chemical Concentration Distribution Model
The advection-diffusion model proposed in [15] is adopted in this paper. In the theory of eddy diffusion in the atmosphere, the diffusion rate at a certain point under each pressure level is parameterized by the eddy diffusivity, K cm 2 /s. The concentration of the chemical substance, c g/mL, can reach a steady state in a homogeneous wind field in which the eddy diffusivity is the same at every point [15]. As shown in Figure 1, the continuous chemical source located at point O is releasing at a rate, q g/s. The released chemical substance is advected by a homogeneous wind [16], the speed of which is v cm/s along the positive x-axis. The sensor, which is represented by point P at (x,y,z), lies in the plume caused by the chemical source. The steady concentration at the point (x,y,z) satisfies: The concentration distribution is symmetry about the x-axis due to the homogenous flow, and then Equation (1) becomes: 2 where β is the angle between OP and OX [15]. Coordinate systems used in this paper. (a) 3-dimensional Cartesian coordinate system; the distance between O and P is d cm, the angle between OP and OX is β.
(b) 2-dimensional Cartesian coordinate system; the case that P is located on the XOY plane, and α is the angle from the positive x-axis to PO. (c) the polar coordinate system in which P is taken as the origin, so that the position of O can be determined based on the position of P.
Considering that the concentrations at the source location and the point extremely far away from the source (i.e., d   ) approximately equal  and zero, respectively. Then, the solution of Equation (2) is: In the case studied here, the chemical source and the sensor nodes are situated on an impenetrable floor. To analyze the concentration distribution on the floor, we first suppose that the chemical source is placed at (0,0,h). Since no chemical matters can diffuse through the impenetrable floor, the concentration at P can be calculated using the method of mirror images [17] as follows: . Then, substituting h = 0 and z = 0 into Equation (4) for the concentration in our case, the concentration at P is given by: More generally, when the chemical source is located at another point (x 0 ,y 0 ) in the 2-dimensional coordinate system, the associated distance will be 22 00 To emphasize the importance of the node locations in determining the source location, the location of the chemical source is transferred into the polar coordinate system in which the sensor location P is taken as the origin, as shown in Figure 1c. Then, the coordinate (x 0 ,y 0 ) is expressed as:   is the angle from the positive x-axis to the broken line pointing from the measuring location to the source. The relationship between α and β is shown in Figure 1b. Obviously, =+    and Equation (5) becomes: The advection-diffusion model introduced above was also referred to as Robert's model in [18]. Some of other studies based on this model can be found in [4,17,19,20].

Deriving the Sensor-Source Distance Function
Based on the concentration distribution model in Equation (7), given the concentration c, flow speed v and the parameters K and q, the distance between the chemical source and the sensor node can be expressed as a function of α by reverse derivation.

  
to Equation (7), we get Then, moving the denominator in the right-hand-side of Equation (8) Since the principle branch of the Lambert W function is single-valued and monotonically increasing, according to Equation (11), d is inversely proportional to c when other variables are kept constant.

Measurement Model
Suppose there is a WSN that consists of n chemical sensor nodes randomly deployed at known global coordinates ( , ), [ where i c represents the theoretical concentration at the location of i x , i e stands for the measurement error which is mainly comprised of the errors caused by the uncontrollable drift of the sensors and the sensors' response to foreign substances due to the poor selectivity of the commonly used chemical sensors [2]. Typically, it can be considered that the error terms for all the sensor nodes at the same sampling time satisfies the normal distribution with a positive mean value [11]. Due to the relatively small variance of the normal distribution of , [ ii ze [6,7]. For a better localization performance, a measurement threshold is set to eliminate the measurements with comparatively low SNRs. The nodes with measurements that are smaller than the threshold th z will be scheduled to hibernate, and the remaining l (l < n) nodes will participate in the localization.

Distributed Least-Squares Estimation Method for CSL Using WSN
In the LSE-based CSL using a WSN, the objective function to be minimized is the sum of squares of the l participating scattered sensor nodes' measurement errors [22], which is as follows: are known measurements and unknown theoretical concentration (since 0 x is unknown during the localization), respectively. By varying the location of 0 x in its domain, i.e., 2 0  x , the distribution of this objective function can be obtained. The logarithm of this distribution is illustrated in Figure 2. There is a global minimum located at the actual source location, while there are some local minima near the global minimum, as well as some saddle points around the local maxima. As shown in Figure 2, there is a global minimum located at the actual source location. It is straightforward that the objective function in Equation (13) would achieve its global minimum, when the l summands, i.e.,  After the estimates of these local estimations are achieved, they are convexly combined to form the global estimate by using an average consensus algorithm. The overall scheme of the proposed method is given in Figure 3, in which the sub-processes are detailed in the following subsections.

Local Estimation of the Source Location
In Section 2, the location of the chemical source has been transformed into the polar coordinate system, in which the location of sensor node is taken as the origin. According to Equation (11), when   the coordinates of the source location can be further expressed as functions of α as follows: For the case that   , the source location is 00 . However, the calculation of Equation (14) involves determining the theoretical concentration c, which is impossible to be obtained due to the unknown source location.
. In the context of probabilistic inference, the points in can be considered as the probable source locations inferred based on , Compared with the domain of 0 x , i.e.,  After i z is substituted for c, Equation (14) becomes an equation set with a single unknown α whose domain is (0,2 ]  . Therefore, as shown in Figure 4, the problem of estimating 0 x can be transformed to another problem of estimating the angle of the vector pointing from i x to 0 x , i.e., 0 i  . Although it is difficult to assign a rigorous theoretical probability of being the actual source location to each point in , if a predictive location of the source exists, it can be empirically considered that the point with the same polar angle as the predictive location has the largest probability of being the actual source location among all the points in , [1, ] i S i l  . Inspired by the process of iterative estimation, the is considered as the prior information (predictive location of the source) for the local estimations at the k-th iteration of localization. Let us define ˆ() i k  as the angle between the vector 0 ( 1) i k  xx and the positive direction of the polar axis. Then, a local estimate of the source location can be considered as: zk. This process of local estimation can be considered as a quasi-maximum a posterior estimation because the maximum posterior probability is determined empirically, rather than theoretically, based on the location of 0 ( 1) k  x . Note that, at the beginning of the localization there is no prior information about the source location. To start the localization, an initial point must be set as the prior information for the first iteration of localization. Since the initial point has little influence on the localization performance (see Section 4.3), it can be randomly initialized.

Combining the Local Estimates to Form the Global Estimate
Section 3.1 presents a framework for the local estimation of the chemical source location based on the prior information about the predictive location of the source. However, it does not detail how this prior information is achieved and how it is accessible to the participating scattered sensor nodes. Moreover, the individual local estimates can only minimize an individual summand of the objective function in Equation (13) with the maximum probability. Thus, to comprehensively consider all the summands of the objective function in Equation (13), these local estimates are combined to form a global estimate, which serves as the prior information about the source location at the next iteration of localization.

Convex Combination of the Local Estimates
A convex combination of the local estimates at the k-th iteration of localization, i.e., 0 () i k x , is considered as the k-th global estimate 0 () k x : (16) where () is calculated using Equation (15).
To set reasonable weights for the local estimates in Equation (16), let us consider the positional relationship between 0 x and the points in , [1, ] i S i l  from which the local estimates are selected.
When the different contaminated measurements , [1, ] i z i l  are substituted for c in Equation (14), the size of the associated ovals formed by the points in different refined sample spaces (i.e., , are different. The oval generated by substituting a larger i z for c in Equation (14) is smaller than that generated using a smaller i z . This is mainly because d will decrease if c increases and other parameters are kept constant according to Equation (11). Moreover, as shown in Figure 5, the points on a smaller oval are generally closer to 0 x than those on a bigger oval. Therefore, it is reasonable to construct a direct proportion between () i wk and () i zk to make the sequence of 0 () k x approach 0 x . The weights () i wk are tentatively set as the proportion of () i zk in the sum of all the l measurements, as in Equation (17), and achieve good performance in the results presented in Sections 4 and 5:

Distributed Average Consensus Algorithm
It is readily seen that there are two individual accumulating operators in Equation (16) and Equation (17). Moreover, the accumulating results, i.e., 0 () k x and the denominator in Equation (17), should be accessible to all the participating scattered sensor nodes. The centralized processing method in which a fusion center gathers the summands from and then transmits the sum back to all the sensor nodes is not power-efficient for WSN, and thus is avoided in this paper. Since the sum can be considered as the product of the amount and the average of the summands, the process of accumulation can be transformed to an averaging process and an additional multiplication between the amount and the associated average as follows: . Then, the distributed average consensus algorithm in [14] can be used to calculate the averaging term in Equation (18), which is involved in the convex combination of local estimates twice, i.e., in calculating Equations (16) and (17). The distributed average consensus algorithm means that only local communications within the neighboring set of each participating scattered sensor node are needed during the averaging procedure [23]. The steps of adopting this algorithm to our specific problem are elaborated as follows: Step 1: Generate the corresponding graph and the Laplacian matrix. Suppose each participating scattered sensor node can only communicate with the sensor nodes within its neighboring range. A graph describing the corresponding network architecture is generated. In the generated graph, there is an edge between two nodes which can communicate with each other. Then, the Laplacian matrix (denoted by L) of the graph is calculated as L = AA T , where A is the incidence matrix of the graph. Since A is a matrix with l rows and g columns, where g is the number of edges in the graph, L is a l-dimensional square matrix. Figure 6. The associated graph of the mesh architecture in Figure 1d. The sensor nodes are denoted as circles, in which the numbers are the indexes of the associated nodes. Communication links within the neighboring ranges are represented by solid lines, and the indexes of these edges are given aside them.
For example, as shown in Figure 6, there are l = 18 vertexes and g = 32 edges in the graph corresponding to the mesh architecture in Figure 1d, so the corresponding Laplacian matrix is a 18-dimensional square matrix.
Step 2: Calculating the modulus matrix. A modulus matrix, which is denoted by M, is calculated as follows: where ( ), [1, ] i L i l   is the i-th largest eigenvalue of L, and I is the identity matrix. Then, the nonzero elements in the i-th row of M, which is denoted by i M , are transmitted to i x . These nonzero elements corresponds to i x and its neighboring nodes. For example, in Figure 5, there are three nonzero elements in 1 M , i.e., 11 12 , MM and 15 M , which correspond to 12 , xx and 5 x , respectively. Since only a few nodes are connected with each participating scattered sensor node, the matrix M is sparse and the number of the elements to be transmitted is small.
Step 3: Performing the iteration of averaging. Since the proposed localization algorithm and the average consensus algorithm used here are both iterative algorithms, it is important to distinguish them from each other. The iteration in the former is called the iteration of localization (counted by k), while the iteration in the latter is called the iteration of averaging (counted by t). In the k-th iteration of localization, two whole processes of iterative averaging (i.e., from the initial state to the convergent state) should be conducted to calculate iteration of averaging. This step actually realizes a matrix multiplication as follows: where () ft is a column vector that consists of ( ), [ iterations of consensus, has been proved in [14].
If the variation of the global local estimation, i.e., performing the iterations of consensus. Thus, some information exchanges between the sink and the participating scattered sensor nodes are engaged in the first two steps. However, the first two steps are conducted only once at the first iteration of localization, i.e., k = 1. Therefore, the communication burdens in these two steps are acceptable. Basically, there should be some routine communications between the sink and the scattered nodes to guarantee the conventional operations of the WSN.

Simulation Results
In this section, the localization performance of the proposed DLSE method is assessed through simulations. First, the basic simulation setup, which is used to generate Figures 2 and 5 and served as the prototype of the following simulations, is described. Then, the process of distributed averaging, which is the premise of distributed implementation of DLSE, is demonstrated. Afterwards, the localization performance of DLSE is compared with the trust-region-reflective algorithm which is recommended for solving centralized nonlinear LSE problem. Finally, the ability of locating a remote chemical source using DLSE is assessed.

Basic Simulation Setup
The simulated WSN was composed of 300 nodes that were randomly distributed over an 80 m × 80 m square. The lower left corner and bottom margin of this square were taken as the origin and abscissa axis of the global Cartesian coordinate system, respectively. In Figures 2 and 5, the chemical source was located at (300,4000) cm, and the chemical substances released from the source was diffused by the flow with a near-surface velocity of v = 70 cm/s. The eddy diffusivity, K, was set as 10 4 cm 2 /s. The unit of the release rate in Robert's original paper [15] is g/s, and thus the associated unit of concentration is g/cm 3 . At 20 °C and standard atmospheric pressure, the mass concentration in g/cm 3 can be transferred to the volume concentration in the most commonly used unit, ppm, by multiplying an exponential 9 (

) 10
M  , where M is the molecular weight of the released substance. The release rate q in our simulations was set as 114 g/s, which can cause the concentration range up to dozens of ppm when M = 46 (i.e., the molecular weight of ethanol) and the sensor-source distance is tens of meters. Referring to [7], the mean and standard deviation of the measurement errors were 10 −5 kg/m 3 and 8 × 10 −6 kg/m 3 . Correspondingly, the concentration threshold was set as 10 −5 kg/m 3 to eliminate the measurements with extremely low SNRs. By setting a presumed chemical source at the center of each of the mesh grids on the node-deployment square, which is denoted as 0  x , a group of , [1, ] i c i l   can be calculated by substituting 0  x for 0 x in Equations (6)  to Equation (13). The logarithmic transform of these values are illustrated in Figure 2. According to Figure 2, a global minimum appears near the actual source location, while there are also several local minima and saddle points at other locations in the square. To generate Figure 5, the measurements of 30 participating scattered sensor nodes, i.e., , [1,30] i zi  were substituted for c in Equation (14).

Demonstrating the Process of Distributed Averaging
As mentioned in Section 3, the convex combination of local estimates is calculated by using a distributed average consensus algorithm, which serves as the prerequisite of the distributed implementation of our method. Thus, before assessing the localization performance of our method, it is necessary to demonstrate the process of distributed averaging.
At the k-th iteration of localization, the coordinates of the local estimates ˆ( ), [ , at different iterations of localization were illustrated in Figure 7. As shown in Figure 7, the distributed averaging of ˆ( ), [1, ] i k i l  x was performed eight times, which were associated with eight iterations of localization. At each time of the averaging, the abscissas maintained by different participating scattered sensor nodes converged to the abscissa of the global estimate, i.e., 0 ( ), [1,8] kk  x , in about ten iterations of consensus. Moreover, the abscissas of these global estimates converged near the abscissa of the actual chemical source which was located at (300,4000) cm.

The Influence of the Initial Point
The trust-region-reflective (TRR) algorithm [24], which is recommended in Matlab, was used to solve the standard centralized nonlinear LSE of the source location and as the benchmark to assess the performance of our method. The TRR algorithm, which is based on the interior-reflective Newton method, is a subspace trust-region method. Trust region approaches approximate the function to be minimized with a simpler function, which can reasonably reflect the behavior of the original function in a neighborhood around the current point, i.e., the trust region. Then, the trust-region sub-problem, which minimizes the simpler substituted function over the trust region, is solved at each iteration of trust region approaches. The trust-region methods are local search methods since the solution can only move in the trust regions. Thus, TRR is a typical local search algorithm, although the structure of the nonlinear LSE problem is exploited in TRR to enhance efficiency.
The influence of initial point on the localization performance of TRR and DLSE were assessed in two different cases, and the results are illustrated in Figure 8. Note that in Figure 8, the initial points are directly connected to the associated final estimates for TRR, while the intermittent global estimates in a single localization were connected to each other for DLSE. As shown in Figure 8, TRR succeeded several times when the initial point was near the source location, otherwise, it would be stagnated at some of the local optima or failed in starting the search. Contrarily, DLSE succeeded from all of the initial points in the two cases, and achieved relatively small localization errors compared to the large size of the node deployment area. Therefore, the locations of initial points hardly influence the localization performance of DLSE. In other words, DLSE can overcome the problem of poor convergence caused by the multiple local minima and saddle points of the objective function in Equation (13), which would obstruct the convergence of TRR and other local search algorithms.

The Performance of Locating a Remote Chemical Source
As shown in Figure 8d, DLSE succeeded in locating a remote chemical source which was located far away from the node deployment area. However, the results presented in Figure 8d was obtained based on limited times of simulation of locating the same source. To get a more convictive assessment of the performance of locating a remote chemical source, the actual chemical source was located at different locations far away from the node deployment area. For each location of the chemical source, 100 times of simulations were conducted, with the measurements collected by the nodes which were randomly re-deployed at the beginning of each simulation. TRR was not considered in these cases due to its poor convergence even when the actual chemical source was located in the node deployment area as shown in Figure 8a. Since the initial point can hardly influence the performance of DLSE, it was set as the centroid of the node deployment area, i.e., (4000,4000) cm, in all of these simulations. The statistical results of these simulations were represented with error bars in Figure 9.
As shown in Figure 9, all of the chemical sources were successfully located with relatively small localization errors compared with the large sensor-source distances. The localization errors became larger with an increase in the distance between the source and the node deployment area. Because the theoretical concentration is inversely proportional to the sensor-source distance when other variables are kept constant according to Equation (11), the SNR of the measurements decreased, which would influence the localization performance, as the distance between the source and the node deployment area increased. However, the performance of any localization method based on concentration measurements might be influenced by the SNR of the measurements. Figure 9. The statistical performance of locating a remote chemical source using DLSE. All the ordinates of these chemical sources were set as 4000 cm, thus, only the abscissas of the sources were displayed to denote the associated sources. The center and half-length of the error bar are the mean and standard deviation of the localization errors, respectively.

Real Experiment Results
In addition to the simulations, our methods were assessed in realistic experiments using a WSN which consists of twenty five real sensor nodes.

The Sensor Nodes and Realistic Environment
The sensor node was constructed based on the C51RF-CC2431 module (Wireless Dragon, Co. Ltd., Chengdu, China). Each of the nodes was equipped with a MiCS-5521 (SGX sensor Technology, Co. Ltd. Neuchatel, Switzerland) gas sensor to measure the chemical concentration. The nodes can communicate with the sink node via the wireless ZigBee protocol. The sensor nodes are able to reliably transmit the real-time sensing voltages of MiCS-5521 to the sink node. Due to the limitation of our currently available nodes, the wind speed was measured by two additional WindSonic anemometers (Gill, Co. Ltd., Hampshire, England). The general picture of the experimental installations is shown in Figure 10. A reduced scale wind tunnel, which is a cuboid with 300 × 400 cm 2 area by 90 cm height, was built up to create an approximately homogeneous wind field. The inlet and outlet of the wind tunnel are two parallel sides of the cuboid. An array of electrical fans was integrated in a flat rack, which was mounted in the wind inlet. Clean air was blown into the tunnel through the wind inlet. More details about the wind tunnel can be found in [25]. A pump and a flow controller were used to blow saturated ethanol vapor from the bubbler into the wind tunnel through a long PVC tube. The density of the saturated ethanol vapor in the bubbler was maintained by a thermostatic bath. The saturated vapor pressure and density can be calculated using Antoine's equation and the Clausius-Clapeyron's equation, respectively. Then, the product of the volume flow rate and the density of the saturated ethanol vapor can be considered as the theoretical release rate of the chemical source.

Sensor Calibration
Before the sensor readings can be used to represent the concentration, the sensors should be calibrated to find an accurate mapping from the sensor readings to the concentration. Due to the controlled homogenous wind field, the concentration distribution and the sensor readings are approximately steady after a period of time [4]. Since this specific work condition approaches the closed sampling space [26], during the calibration process the sensor nodes were enclosed inside a closed glass box whose cubage is 108 L, as shown in Figure 11. The box was placed in the wind tunnel to make the sensors work under the same humidity and temperature like those in the process of localization. Then, multiple calibration points can be created by dropping in certain volumes of absolute ethyl ethanol and blowing them to accelerate their evaporation with convective wind. The relationship between the sensitivity S, i.e., the ratio of the sensing resistance to the baseline resistance, and the concentration C s (in ppm) is well fitted to the transform relationship log log s S a C b    , in the range [10,1000] s C  ppm [27]. Thus, the unknown parameters, i.e., a and b, can be identified by fitting the steady sensing voltages and the associated actual concentrations to the calibration curve. Note that we refer to the final steady concentration, which was calculated by substituting the approximately steady voltage to the identified transform relationship, as the measurement in the rest of this section.

Identification of the Release Rate and the Eddy Diffusivity
According to Equation (8), when   , which means the sensor node lies on the positive abscissa axis of the Cartesian coordinate system in Figure 1b, the concentration of the node is: where qK   can be considered as a constant parameter during the process of localization.
However, the premise is that the value of q and K are identified before the localization. Suppose a standard chemical source, which was releasing at a controlled release rate 0 q = 0.01 g/s, was placed in the wind tunnel and taken as the origin of a Cartesian coordinate system like the one in Figure 1b. All the sensor nodes were placed along the positive abscissa axis. The abscissas of these nodes were set on the premise that the reciprocals of the distances from these nodes to the standard source should range from 0.003 to 0.008 in increments of 0.0005. The relationship of the concentration measurements of these nodes and the reciprocals of the associated distances was fitted into a straight line, which was shown in Figure 12a. Note that in real cases, the source to be located and the standard source may release simultaneously, so the concentration measurements used to identify K should be the increments of the measurements caused by the standard source. With the identified value of K, we can further identify q by measuring the concentration along the centerline of the plume. The first step is to determine the unknown centerline of the plume caused by the source to be located. Because the distribution of the concentration measurements collected along the line perpendicular to the wind direction satisfies the normal law of errors [15], the unknown centerline can be determined as the line passing the maximum concentration point along the downwind direction. To verify this normal distribution, the twenty five sensor nodes were placed along a line vertical to the plume centerline. The measurements of these sensor nodes were plotted in Figure 12b. Since the measurements on the centerline of the plume caused by the source to be located satisfy ,

Experiment Results
After the above-mentioned preliminaries, the performance of DLSE and TRR were assessed with real measurements. Note that, to avoid the complicated programming with the Zigbee protocol, the concentration measurements were transmitted to the workstation and the distributed implementation of DLSE was realized off-line. However, the online distributed implementation of DLSE is applicable after a further study on hardware design and embedded programming, which was not focused in this paper. Four scenarios were created by changing the flow rate of the bubbler from 400 sccm to 1000 sccm in increments of 200 sccm, while keeping the temperature of the thermostatic bath, the wind speed at 60 °C , 70 cm/s, respectively. The concentration threshold is set as 10 ppm to eliminate the measurements with low SNRs. Based on the measurements collected in each of the groups, one hundred of tests with different randomly selected initial points were conducted using TRR and DLSE. The boxplot of the localization errors in these groups were shown in Figure 13. As shown in Figure 13a, although the minimum localization errors obtained using TRR were comparatively small, their overall distribution was highly fragmented and not symmetrical with respect to the median. Thus, the good performance of TRR is not guaranteed if the initial point is not properly chosen. In Figure 13b, DLSE succeeded in all of the tests with considerably high accuracy. The localization errors obtained based on the same group of measurements were generally symmetrically distributed with respect to the associated medians. In addition, higher flow rate of the bubbler, which means larger release rate of the chemical source and more participating scattered sensor nodes, results in smaller mean of localization errors.

Conclusions
In this paper, a distributed chemical source localization method is proposed for a WSN. Aiming to minimize the summands in the objective function of LSE-based CSL using WSNs, i.e., the sum of squared measurement errors, individual local estimations of the source location are conducted on the participating scattered sensor nodes. This scheme avoids gathering the raw measurements and enables the realization of a distributed estimation method. The local estimates are convexly combined into the global estimate, which is used as the prior information for the local estimation at the next iteration of localization. The global estimate includes the information of all the local estimates in different degrees and enables collaboration among all the participating scattered sensor nodes. Our method shows a global convergence property with fast convergence rate, even when the actual chemical source is far away from the deployment area of the WSN. Thus, the initial point of our method can be preset as a random point before the process of measurement.
Our method iteratively incorporates the global information about the source location in a distributed manner so as to achieve a global convergence property. When the distribution models of radioactive substance and acoustic energy are utilized, our method can also be adapted to locate the radioactive and acoustic sources. Future works may cover deriving a theoretical probability of being the real chemical source for the local estimates. When such a probability is used as the weight of the associated local estimate, our method could evolve into a distributed particle filter.