# Distributed Source Localization in Wireless Underground Sensor Networks

## Abstract

Node localization plays an important role in many practical applications of wireless underground sensor networks (WUSNs), such as finding the locations of earthquake epicenters, underground explosions, and microseismic events in mines. It is more difficult to obtain the time-difference-of-arrival (TDOA) measurements in WUSNs than in terrestrial wireless sensor networks because of the unfavorable channel characteristics in the underground environment. The robust Chinese remainder theorem (RCRT) has been shown to be an effective tool for solving the phase ambiguity problem and frequency estimation problem in wireless sensor networks. In this paper, the RCRT is used to robustly estimate TDOA or range difference in WUSNs and therefore improves the ranging accuracy in such networks. After obtaining the range difference, distributed source localization algorithms based on a diffusion strategy are proposed to decrease the communication cost while satisfying the localization accuracy requirement. Simulation results confirm the validity and efficiency of the proposed methods.

Wireless underground sensor networks, source localization, Chinese remainder theorem, time-difference-of-arrival (TDOA).

## I Introduction

Wireless underground sensor networks (WUSNs) are an important extension of terrestrial wireless sensor networks, as they can be used to estimate the location of earthquake epicenters, underground explosions, microseismic activities in mines, etc. Normally, the sensor nodes in WUSNs are buried underground and they exchange information wirelessly via the dispersive underground channel. Some experimental results with WUSNs are reported in [5].

In terrestrial wireless sensor networks, time difference of arrival (TDOA) measurements are widely used for node localization [1]. Because of the physical characteristics of the dispersive underground channel and the heterogeneous network architecture of WUSNs, source localization for WUSNs based on TDOA is more challenging [2]. In a dispersive medium, we cannot directly obtain the range differences between the source and sensors from TDOA measurements since the propagation velocity is a function of frequency; different frequency components will have different propagation delays [3].

Determining the location of sensor nodes is important in many practical applications of wireless underground sensor networks. The objective of a positioning system is to determine accurate node locations with low complexity and communication cost. Localization algorithms for traditional terrestrial wireless sensor networks can be classified into two types: range-based methods and range-free methods. Range-based methods usually have higher location accuracy than range-free ones while demanding additional hardware cost [4].

In [6], a distributed TDOA estimation method that relies only on radio transceivers without other auxiliary measurement equipment was presented. Ultra-wideband (UWB) signaling can be used to accurately achieve time of arrival (TOA) or TDOA measurements, which has the advantages of low-cost and penetrating ability, but also has the weakness of short-range. TDOA based algorithms provide high localization accuracy, and represent a practical method for estimating range differences and source positions in WUSNs. However, this method also faces many challenges. In particular, limited range and directionality constraints decrease the accuracy of range difference estimation. We notice that the TDOA can usually be obtained from the measurement of a signal’s phase which is susceptible to phase ambiguity problems. The Chinese remainder theorem (CRT) offers a closed-form analytical algorithm to calculate a dividend from several of its corresponding divisors and remainders, and can be applied to solve the ambiguity problem here. In our ranging application, the remainders in the CRT are the measured “remainder” wavelengths, the divisors are the measuring wavelengths, and the dividend is the range difference to be estimated. However, directly using the CRT is not feasible due to its over-sensitivity to noise, i.e., a small error in a remainder can lead to very large error in the estimated dividend. To avoid this weakness of the CRT, we propose to use a robust Chinese remainder theorem (RCRT) algorithm to estimate the range difference, in which the dividend can still be reconstructed with only a small error if the errors on remainders are bounded within a certain level [7]. As a result, the range differences or TDOAs can be robustly estimated from noisy measurements in WUSNs by using the RCRT.

After obtaining the range differences using the RCRT, we can estimate the source position based on statistical signal processing methods. For traditional terrestrial wireless sensor networks, TDOA based localization algorithms are normally implemented in a centralized way. In the centralized solution, all nodes relay their TDOA measurements to a fusion center, which uses a conventional localization algorithm to obtain the source position, and then sends the global estimate back to every node. This strategy requires a large amount of energy for communications [9] and has a potential failure point (the central node). Distributed strategies are an attractive alternative, since they are in general more robust, require less communications, and allow for parallel processing. To address this limitation of centralized processing, we propose distributed source localization algorithms using a diffusion strategy in this paper. Diffusion algorithms were proposed in [9, 10, 11] and are applicable to distributed implementations since nodes communicate in an isotropic manner with their one-hop neighbor nodes, and no restrictive topology constraints are imposed. Thus, the algorithms are easier to implement and are also more robust to node and link failures. This approach allows nodes to obtain better estimates than they would without cooperation.

The rest of the paper is organized as follows. Section II establishes the mathematical model of the problem, and derives the proposed ranging method based on the RCRT. Section III gives the distributed source localization algorithm based on the diffusion strategy. Simulation results are given in Section IV, and conclusions are drawn in Section V.

## Ii System model and TDOA estimation via the robust CRT

Due to the large attenuations in underground environments, experimental results show that underground to underground (UG2UG) communication is not feasible at the 2.4GHz frequency band [12]. Underground communication becomes practical only at lower frequencies. As reported in [13] and [14], a WUSN system operating at 433MHz with a maximum transmit power of +10dBm usually achieves a communication range of around one meter for UG2UG communication and more than 30m for underground to aboveground (UG2AG) communication. These communication ranges have already exceeded the wavelength of the transmitted signal, which results in phase ambiguity when the distance difference is directly calculated from the phase. In this section, we first propose a method based on the robust CRT to resolve this phase ambiguity when computing the range distance.

Consider a WUSN with sensor nodes at known positions , , receiving a signal from a source at an unknown position through a dispersive medium, as shown in Fig. 1. The distance between the source and the -th sensor is

(1) |

The range difference (RD) between the -th and -th receivers, denoted by , is .

If the medium is non-dispersive, the propagation delay for any frequency is constant. However, in a dispersive medium, the signal propagation velocity is a function of the signal’s frequency, denoted by for frequency ; that is, the propagation delay is frequency dependent. We denote the delay of propagation at frequency for sensor as

(2) |

We assume that the source transmits a sinusoid signal at frequency . The -th sensor receives the signal as

(3) |

where is the noise at sensor , and are independent white Gaussian noise processes. Similarly, the received signal at sensor is

(4) |

By taking the cross-correlation of and , we get

(5) |

It is easy to show that is an asymptotically unbiased estimator of [17], for which it follows that

(6) |

where the range difference is contained in the phase of . We denote the phase of by , i.e.,

(7) |

One can determine from . However, there are two issues of concern with this approach: 1) the value of is folded by and 2) the measurements are noisy, i.e.,

(8) |

where is the quotient (folding integer) and is the noise at frequency .

For issue 1, since , no matter how large the actual is, the RD converted directly from is always within one wavelength of the signal, which is . Therefore, cannot be determined uniquely from a single . We add more constraints to confine the solution space by measuring the phase at different frequencies, . The Chinese remainder theorem (CRT) provides a solution to evaluate a dividend from its remainders. We can regard as the dividend, as the divisor, and as the corresponding remainder, i.e.,

(9) |

CRT tells that a positive integer can be uniquely reconstructed from its remainders modulo positive integers , if , where lcm denotes least common multiple. It seems that the CRT gives a perfect solution to this problem. However, when we consider measurement noise, the traditional CRT [16] is not suitable because it is sensitive to noise. A small error in the remainder can result in a large error in the estimated dividend. Therefore, we adopt a robust CRT method [7]. Theorem 1 in [7] proves that the robust CRT can tolerate an error in that is bounded by , where is the greatest common divisor (GCD) among the divisors. To be more specific, if all the remainder errors are not greater than the error bound , the estimation error of the unknown dividend is upper bounded by , i.e.

(10) |

Based on this robust CRT, we provide the following solution to solve
the problem:

Step 1. Estimate by calculating the phase
of . Denote the estimated phase by , and the corresponding distance converted
from the phase by

(11) |

Note that the CRT is commonly expressed over the ring of integers while the distance is a real number. Therefore, we extend the algorithm from the integers to the reals by introducing a real common factor among the divisors as in [18]. We choose those communication frequencies so that the , have a real common factor of which satisfies

(12) |

where are co-prime integers, i.e., , for , and denotes GCD. According to the prior discussion, we have the following equation:

(13) |

where denotes the measurement errors of
with , and is the error bound.

Step 2. For notational convenience, define the following
symbols:

(14) |

(15) |

Calculate and find the sets,

(16) |

where

(17) |

is the solution space of the quotients.

The set can be regarded as the optimal combination of the
quotients and with which the difference of the estimated
dividends achieves its minimum.

Step 3. Let denote the first element
of the 2-tuples,

(18) |

Calculate the intersection set of :

(19) |

In [7], it was proved that if the error bound is less than , the set contains only the true value of , i.e., . In addition, can also be determined from , that is, if , then for . Therefore, with all the quotients being determined correctly, the error of the estimated is therefore bounded by (10), where is obtained by averaging the estimates corresponding to different wavelengths, i.e.,

(20) |

Remark: One may notice that to find the set in (16) we need to search among the possible values of , which is a 2-D search problem with the order of and requires a high computational complexity. However, the complexity can be decreased by using a fast algorithm proposed in [7] and [8] to a 1-D search problem with the order of only . One can easily apply the fast algorithm to the application in this paper. We do not describe this fast algorithm herein since this is not the focus of this paper.

## Iii Diffusion Algorithms for Source Localization

After obtaining the range differences using the RCRT, a localization algorithm can be used to estimate the source position. Herein, we propose an algorithm for source localization that employs diffusion strategies proposed in [9, 10, 11]. First, the nodes in a cluster send their measurements to the cluster head. Then, the cluster head determines a local estimate using these measurements locally. After obtaining the local estimates, cluster heads exchange their local estimates to achieve diffusion. Before describing the distributed diffusion process, we first discuss the global WLS estimate when all measurements are sent to a fusion center for estimating the source location.

### Iii-a Global Weighted Least Squares Problem

Consider a set of cluster heads and sensor nodes (each cluster head is associated with sensor nodes) spatially distributed over some region with known locations ’s (we consider the 2-D Cartesian coordinate system). The objective of the network is to collectively estimate an unknown deterministic column vector - the location of a source. All the nodes (cluster heads and sensor nodes) in the network can measure the signal transmitted from the source. The sensor nodes transmit their measurements to the corresponding cluster head and each cluster head forms a set of TDOA measurements with itself being the reference. Then, cluster heads can send their TDOA measurements to the fusion center. The TDOA measurements are formed one at a time by comparing the signal from the cluster head and the signal from a sensor node, thus leading to uncorrelated estimates if the estimation period is longer than the typical coherence time of the mobile radio channel. Finally, the fusion center obtains altogether TDOA based range difference measurements ’s in the network.

The scalar model for TDOA based range difference is given by

(21) |

stack all the ’s into a vector and write

where , and . Each can be alternatively denoted as to reflect its location in vector . The corresponding covariance matrix for is denoted as and .

The global weighted least squares estimator for estimating given is thus

(22) |

Assuming the ’s are Gaussian, then the covariance of attains the corresponding CRLB, which is given by [15]

(23) |

where

(24) |

and denotes the true position of the source.

For the TDOA measurements, is a matrix and the elements of are given by

where we assume involves nodes (a sensor node and a cluster head) and .

### Iii-B Local Weighted Least Squares Problem

For each cluster head , it has access to limited data from its neighbors. It can then solve the WLS problem locally as

(25) |

where ’s are the associated weights for node and if is not accessible by node . Let denote the matrix with elements . We require , where denotes the column vector with unit entries.

A local estimate can also be written as

(26) |

where ( is an vector with a unity entry in position and zeros elsewhere) and . can be estimated at .

The covariance matrix associated with is given by . Here contains only local information due to the selection ability of .

The estimation fusion method can then be used to fuse these local estimates in a distributed manner by utilizing the local covariance matrices as proper weights and the estimation fusion method can be shown to achieve the performance of the global estimation. However, it involves covariance estimation and matrix inversion.

### Iii-C Diffusion Algorithm

Besides estimator fusion, we can also use a diffusion algorithm to perform distributed estimation. For each cluster head , at the th time epoch, it exchanges its local estimate with its neighboring cluster heads and updates its local estimate using a diffusion algorithm:

(27) |

where the ’s are the diffusion coefficients. Eq. (27) can be considered as a weighted average of the local estimates in the neighborhood of node . Assume all the local estimates are unbiased. Then in order for to be unbiased, we require where . The diffusion process is repeated until all the local estimates have converged, i.e., , where is a (small positive) design parameter.

One possible choice for the weights is to consider the degree of connectivity, which is

(28) |

where denotes the cardinality of cluster head ’s neighbors (also cluster heads) and denotes the set of neighboring cluster heads of head . Such choice has been observed to yield good results for the diffusion algorithm in general [9].

However, this method does not consider the reliability of different local estimates. The reliability of a local estimate is reflected in its associated covariance matrix. But estimating the covariance matrix is not an easy task for a nonlinear weighted least-squares estimator (WLSE). Here we propose a method for setting appropriate ’s which reflects the reliability of different local estimates to a certain extent without requiring use of the covariance matrices.

Since each local estimate contains certain errors, when sorting them in respective dimensions (each dimension of is treated separately), the middle ones are more reliable. Therefore, for cluster head , at the th time epoch, it first finds the median along respective dimensions among its received local estimates ’s. Then for each obtained local estimate, head calculates if . Otherwise, = 0. is a parameter which controls how rapidly the weight decays as grows. The diffusion coefficient is then set to

(29) |

Here we explicitly indicate the time epoch of and in the superscript. It can be seen that the larger the deviation between a local estimate and , the smaller the weight assigned to this local estimate and vice versa. Obviously, .

Another method of interest is to use an optimization technique. That is, to set ’s such that the trace of the covariance matrix of is minimized. We start by examining the first time epoch . Then, we have

(30) |

where we have used (26).

The covariance matrix of is thus given by

(31) |

and its trace is

(32) |

where and . To be unbiased, we also require .

Therefore, to find the optimal in the sense of minimizing the diffusion covariance matrix, we need to solve an optimization problem:

(33) |

This problem is convex if and thus can be solved fast and efficiently.

After determining at time epoch , is updated as and will also be updated accordingly. Then the above optimization process can be performed iteratively until estimates converge.

This optimization method can enhance the distributed fusion performance at the cost of slightly higher computational complexity compared with simply setting according to (28).

To decrease the computational cost, the optimization process can be implemented only once for the first diffusion at each node, and the weight remains the same for the latter diffusions. Since after one diffusion, the updated local estimates becomes much less dissimilar and thus the weights will have much less influence on the followed diffusions.

## Iv Simulation results

In this section, simulation results are given. We first study the ranging performance using the RCRT in Subsection A. Then, the localization accuracy is introduced in Subsection B.

### Iv-a Ranging Performance

Assume that the signals are transmitted at 3 frequencies, i.e., =3, =80, =, and the corresponding three dividers are which represent the wavelengths . According to the RCRT, the maximum estimation distance is =32640mm. Each trial of simulation generates random integers , which are uniformly distributed over . And there are 1000 trials for each SNR. The result is shown in Fig. 2.

Then we consider the effect of on the performance. We set , respectively and fix =. The result is shown in Fig. 2 with 10000 trials for each SNR. According to the result, changing does not have significant influence on the relative error of distance estimation.

Next, we compare the performance under different values of with constant and . In the simulation, we fix =50, =3, and let be , , , respectively. Fig.3 demonstrates that the estimation error increases with . The results from Figs. 2 and 3 can be explained by equations (11) and (12): the phase measurement error is amplified by and , i.e., the error of is , where denotes the phase measurement error of . The larger is, the larger error results. However, does not affect the performance because the robustness of the algorithm also linearly increases with which cancels the performance deterioration of the increased . (The error tolerance of the algorithm is )

In addition, we consider the scenario in which different sets of are compared under the constraint of constant maximum range . In Fig. 4, We choose =, , , respectively. The simulation results suggest that the performance is better if the differences between the are smaller.

Fig. 5 demonstrates that using more wavelengths results in better ranging performance. In this simulation, we fix the maximum estimation distance , and vary . We consider the cases when =2,3,4, respectively, with =50 and =3465mm. are set to , , , respectively. Simulation results demonstrate that our RCRT based ranging scheme can estimate the range differences with high accuracy.

### Iv-B Comparison Results for Distributed Source Localization

We now show simulation results for source localization. The algorithms compared are: 1) Global: the global weighted least square estimator (III-A), 2) Diff (con): the diffusion algorithm with coefficients set by considering connectivity (28), 3) Diff (wei): the diffusion algorithm with coefficients set by weighting (29), 4) Diff (opt): the diffusion algorithm with coefficients set by optimization (33) and 5) Local: simple average of all the local estimates, i.e., .

The root mean square error (RMSE) is used as a performance metric, which is defined as . The Cramer-Rao Lower Bound (CRLB) on the RMSE based on the entire data (equals to for Gaussian measurement noise) is also presented as a benchmark. Each simulated point is averaged over 200 runs.

The simulated network consists of cluster heads that are regularly deployed at grid points. The distance between neighboring cluster heads is set to 50 (the units are meters, here and below). Each cluster head has associated sensor nodes which are distributed uniformly around the corresponding cluster head. The source location is fixed at in all the simulations. The TDOA measurements are generated according to (21) with ’s being Gaussian noises. is set to with . is set to 1. The initial point for the WLSE is always set as the center of the deployment area.

Here we consider the scenario in which each cluster head exchanges its local measurements with its neighboring cluster heads to perform a local estimate, and then exchanges its local estimate with its neighboring cluster heads to perform diffusion until convergence. is set to

(34) |

where

(35) |

Here, denotes the weight assigned with respect to cluster heads and , and represents the weight assigned to the th measurement used by the cluster head .

First, we fix and vary the value of , which changes from 5 to 50 with a step size of 5. The RMSEs of respective algorithms are shown in Fig. 7. It can be observed that Global is the best which can attain the CRLB, local is the worst and Diff (con) is always better than local. In general, Diff (opt) is better than Diff (wei), and Diff (wei) is better than Diff (con). The performance improvement of Diff (opt) and Diff (wei) compared with Diff (con) comes from the consideration of the reliability of the local estimates. Though the diffusion algorithms are always worse than Global, the performance differences are not significant. As grows, all the algorithms perform better.

Fig. 8 shows the corresponding average CPU times of respective algorithms except Diff (opt) whose CPU time is typically 10 times that of Global due to its numerical optimization nature (the same below). It can be seen that Diff (con) is very time efficient with a CPU time almost the same as local. Diff (wei) consumes a little more CPU time than Diff (con), while the CPU time of Global is much larger. This demonstrates the advantage of the diffusion algorithm in terms of CPU time and computational complexity. Furthermore, we can say the diffusion algorithm has lower communication cost and computation complexity than the centralized solution, while the localization accuracy is close to the centralized method.

Fig. 9 shows the average number of iterations before convergence of the three diffusion algorithms. We can observe that Diff (opt) requires many fewer iterations compared with Diff (con) and Diff (wei). Diff (wei) needs slightly more iterations than Diff (con), which may explain our observation of a slightly longer CPU time consumed by Diff (wei).

We now fix and examine the effects of the number of cluster heads on the performance of the algorithms. We vary from 4 to 36. The RMSEs of the algorithms are shown in Fig. 10. It is clear that adding more clusters (and thus sensor nodes) will not necessarily improve the estimation performance. This is known as the geometric effect of the localization problem. Since the source location is outside the convex hull formed by respective added clusters, the corresponding local estimates are not good enough. Diff (con) and Diff (wei) thus also show performance degradation. However, Diff (wei) is much better than Diff (con) when more clusters are added while the sensor nodes associated with each cluster head is fixed. The performance of Global is almost unchanged as becomes large. Due to the consideration of the estimation covariance matrices, Diff (opt) shows very good performance. However, it will consume much more CPU time. The corresponding average CPU times and average number of iterations are shown in Figs. 11 and 12 respectively. Similarly, Diff (con) is the most time efficient and Diff (opt) requires the smallest number of iterations.

Then we examine the effect of on the performance of the algorithms. We set and . The RMSEs of respective algorithms are shown in Fig. 13. The relative performance of respective algorithms is the same as before. As enlarges, all of them show performance degradation. The average CPU times and average number of iterations among different algorithms have the same relationship as before and thus the figures are not shown here.

Finally, we examine the choice of on the performance of Diff (wei). We set , and . We generate 200 realizations of the overall TDOA measurement vector, then store and use them for all the corresponding simulations. The corresponding RMSEs are shown in Fig. 14. It can be seen that an optimal exists which results in the minimum RMSE for Diff (wei). However, in a large range of the choice of , Diff (wei) can generate better results than Diff (con).

## V Conclusion

We have presented energy efficient localization schemes that can achieve high localization accuracy in wireless underground sensor networks. These distributed localization algorithms require low computational complexity and energy consumption based on a diffusion strategy. An accurate RCRT based ranging scheme using TDOA to determine range differences between sensors and source that does not require time synchronization is also proposed. It has been shown via simulation results that the proposed localization algorithms achieve excellent localization accuracy with lower communication cost. In future work, we plan to implement our localization scheme in a testbed and verify its performance with an actual WUSN.

### References

- J. Luo, H. V. Shukla, and J.-P. Hubaux, “Non-interactive location surveying for sensor neworks with mobility-differentiated ToA,” in Proc. IEEE INFOCOM, Barcelona, Spain, Apr. 2006, pp. 1-12.
- Z. Sun and I. F. Akyildiz, “Connectivity in wireless underground sensor networks,” in Proc. SECON, Boston, Massachusetts, USA, Jun. 2010, pp. 1-9.
- Y.T. Chan, T.K.C. Lo, and H.C. So, “Passive range-difference estimation in a dispersive medium,” IEEE Transactions on Signal Processing, vol. 58, no. 3, pp.1433-1439, March 2010.
- C. Wang, J. Chen, and Y. Sun, “Sensor network localization using kernel spectral regression,” Wireless Communications and Mobile Computing (Wiley), vol. 10, pp. 1045-1054, June 2010.
- A.R. Silva and M.C. Vuran, “Communication with aboveground devices in wireless underground sensor networks: An empirical study,” in Proc. IEEE ICC’10, Cape Town, South Africa, May 2010.
- W. Zhang, Q. Yin, and W. Wang, “Distributed TDoA estimation for wireless sensor networks,” in Proc. ICASSP’10, Dallas, Texas, Mar. 2010, pp. 2862-2865.
- X. Li, H. Liang, and X. Xia, “A robust Chinese remainder theorem with its applications in frequency estimation from undersampled waveforms,” IEEE Transactions on Signal Processing, vol. 57, no. 11, pp. 4314-4322, Nov. 2009.
- G. Li, J. Xu, Y.-N. Peng, and X.-G. Xia, “An efficient implementation of robust phase-unwrapping algorithm,” IEEE Signal Process. Lett., vol. 14, pp. 393-396, Jun. 2007.
- F. S. Cattivelli, C. G. Lopes, and A. H. Sayed, “Diffusion recursive leastsquares for distributed estimation over adaptive networks,” IEEE Trans. on Signal Process., vol. 56, no. 5, pp. 1865-1877, May 2008.
- F. Cattivelli and A. H. Sayed, “Diffusion LMS strategies for distributed estimation,” IEEE Trans. on Signal Process., vol. 58, no. 3, pp. 1035-1048, March 2010.
- C. G. Lopes and A. H. Sayed, “Incremental adaptive strategies over distributed networks,” IEEE Trans. on Signal Process., vol. 55, no. 8, pp. 4064-4077, August 2007.
- E. Stuntebeck , D. Pompili, and T. Melodia, “Underground wireless sensor networks using commodity terrestrial motes,” in Proc. IEEE SECON’06, Reston, USA, Sep. 2006, pp. 112-114.
- A. R. Silva and M. C. Vuran, “Empirical evaluation of wireless underground-to-underground communication in wireless underground sensor networks,” in Proc. DCOSS’09, Marina Del Rey, CA, Jun. 2009, pp. 231-244.
- A. R. Silva and M. C. Vuran, “Development of a testbed for wireless underground sensor networks,” EURASIP Journal on Wireless Communications and Networking, vol. 2010, Article ID 620307, 14 pages, 2010.
- S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Englewood Cliffs, NJ: Prentice-Hall, 1993.
- T. Koshy, Elementary Number Theory with Applications, 2nd ed., Burlington, MA: Elsevier, 2007, pp. 295-302.
- J.G. Proakis and M. Salehi, Digital Communications, 5th ed., New York, NY: McGraw-Hill, 2008.
- W. Wang and X.-G. Xia, “A closed-form robust Chinese remainder theorem and its performance analysis,” IEEE Transactions on Signal Processing, vol.58, no.11, pp.5655-5666, Nov. 2010.