Modelling and analysis of radio tomography

Radio tomographic imaging (RTI) has recently been proposed for tracking object location via radio waves without requiring the objects to transmit or receive radio signals. The position is extracted by inferring which voxels are obstructing a subset of radio links in a dense wireless sensor network. This paper proposes a refined model for signal attenuation in RTI based on measured data, which is used to provide analytic support for previous qualitative observations. We also provide an analytic method for choosing the weighting of the regularization term, and investigate methods for dealing with negative observations caused by noise.


I. INTRODUCTION
Ubiquitous position awareness is a recurrent theme in many disparate areas of signal processing. One may wish to know one's own position (i.e. navigation), the position of resources, or the position of other entities. In imagebased positioning, the object must not be obstructed and illumination conditions must be good. Radio Frequency (RF) based techniques are a popular alternative. If the object to be tracked has a radio transmitter, then the problem is called source localization [1], [2], [3]. However, in many applications, the object in question may not wish to carry a transmitter (e.g. patients in a nursing home) or may not use a known radio device (e.g. in detection of intruders in law enforcement and military applications). To this end, Wilson and Patwari recently proposed Radio Tomographic Imaging (RTI) as a means of locating an object via radio waves [4], [5].
In RTI, each sensor in a Wireless Sensor Network (WSN) repeatedly sends packets to all other sensors. If an object is physically obstructing a given link, then the link's Received Signal Strength (RSS) will drop relative to a pre-calibrated value. By observing link attenuations, it is possible to determine which voxels are occupied.
This paper provides measured data to refine the RTI model, which is used to provide analytic support for previous qualitative observations. We also provide an The authors are funded in part by the Office of Naval Research. The views expressed in this paper are those of the authors, and do not reflect the official policy or position of the United States Air Force, Department of Defense, or the U.S. Government. This document has been approved for public release; distribution unlimited. analytic method for choosing the weighting of the regularization term, and develop a method for dealing with negative observations caused by noise.

II. WILSON AND PATWARI'S MODEL
This section reviews the system model of [4]. Consider a WSN consisting of K sensors. The convex hull of the WSN is divided into N voxels, which for simplicity are cubes of size δ × δ × δ m 3 . The WSN is fully connected, hence there are M = K(K − 1)/2 unique links. Fig. 1 shows an example, and videos are available at [6].
A calibration step is performed to determine the baseline RSS. Later, the differences in RSS y ∈ R M×1 due to differences in voxel occupancies x ∈ R N ×1 are with weighting matrix W ∈ R M×N and measurement noise n ∈ R M×1 . The RSS in dBm on link m is denoted y m , the occupancy of voxel n is x n , and a link's attenuation is a sum over voxel occupancy times the link attenuation (weight) the voxel would cause if the voxel were obstructed. We now review models for W, n, x (stochastic and deterministic, in turn), and x from [4].
Wilson and Patwari [4] define the weighting W by w m,n = 1/ d(m), d 1 (m, n)+d 2 (m, n) < d(m)+λ 0, else (2) where d(m) is the length of link m, d 1 (m, n) and d 2 (m, n) are the distances from the center of voxel n to the two endpoints of link m, and the tunable parameter λ defines an ellipse with the endpoints of link m as the foci. Increasing λ will increase the minor axis of the ellipse. Typical values of λ are 0.1 to 0.01 feet.
In [4] measured data was used to fit a Gaussian mixture model to the noise, with two zero-mean Gaussians with weights of 0.548 and 0.452 and corresponding standard deviations of 0.971 dBm and 3.003 dBm.
The model in [4] for x is presented in two different ways in different sections of the paper. For the error bound (Cramer-Rao Lower Bound (CRLB)) derivation, x is modeled as zero-mean Gaussian random field with covariance matrix C x , with where d (n 1 , n 2 ) is the distance between voxels n 1 and n 2 , σ 2 x = 0.1 dB 2 , and δ c = 1.3m is the correlation "space constant." In the experimental results section of [4], a cylindrical human model is used, with where χ(n) is the position vector of voxel n (with occupancy x n ) and c o and r o are the center point and radius of the obstruction. Neither (3) or (4) was used to help estimate x.
To estimate x, [4] used a least-squares solution with an additive Tikhonov regularization term, where D d computes the derivative in dimension d and α is a user-determined parameter indicating the relative emphasis on the regularization term. Aside from the regularization term, this is also the Maximum Likelihood (ML) estimate if the noise is Gaussian. Without the regularization term, the matrix W T W is not full rank.
In [4], α was experimentally determined; in Section IV, we will derive the theoretical value that should be used.

III. REVISED MODEL
The RSS measurements in Fig. 2 were recorded using commercially available Sun SPOT sensor motes, placed at a height of 1 meter. RSS was recorded for open and obstructed links, and the values in Fig. 2  The model presented in Section II has several issues: • The scaling of the attenuation values is inconsistent with experimental data in [4] as well as the data in Fig. 2. Using (3) or (4) with (2) leads to link attenuations of small fractions of a dB, whereas a human will cause an attenuation of 5 to 10 dB [4], and Fig. 2 shows an average attenuation of 5 dB. • The weights in (2) decrease with increasing link length, which is not justified by any physical model. (2) interferes with extracting analytical intuition from the CRLB analysis. In Fig. 2, two models are considered: that of (2), proportional to 1/ √ d; and a constant attenuation, independent of link length, The two models are numerically similar, and the difference between the two is much less than the noise. The 1/ √ d fit looks slightly better, but that is skewed by the fact that at large distances, the RSS was already small, so if it happened to drop a lot, no packets were observed.
Statistically, if multiple models are equivalent fits, it is generally preferable to use the simplest model with fewest parameters [7], [8]. This is sometimes called Occam's razor or the law of parsimony. The residual error from the constant fit is shown in the inset of Fig. 2, and the residual from the 1/ √ d fit (not shown) had an almost identical histogram. Thus, since the constant model is simpler and makes more sense from a physical perspective, we prefer its use. This simpler weighting matrix allows us to extract more analytically intuition. As an example, consider a spherical obstruction where A is attenuation (dB) per voxel of obstruction. A link passing through the center of the obstruction would be attenuated by To start, we can use the CRLB on x from [4], since even though the structure of our W is different, the dependence of y on x is functionally the same. Thus, where f v (v) is the Probability Density Function (PDF) of any sample of the noise. For the noise model of [4], γ = 0.548; for Gaussian noise, γ = σ −2 n . Since x depends on θ, we can extend the CRLB on x to a CRLB on θ [9] by left-and right-multiplying the Fisher Information Matrix (FIM) by the matrix ∂x T /∂θ, Using the revised model of (9), This form holds regardless of the model of x.
Our CRLB is distinct from the one in [4] for two reasons. First, in [4], the CRLB was a pixel-by-pixel bound on estimating the amount of obstruction; whereas our bound is on estimating an object's position, size, and attenuation. Second, our bound is based on the revised data model which makes it possible to extract more analytic intuition from the bound.
Consider the matrix W T W. Its [n 1 , n 2 ] element is the number of links that touch both voxel n 1 and voxel n 2 . As such, it is diagonally dominant, with diagonal elements equal to the numbers of links crossing the corresponding voxels. Row n of matrix P is the sensitivity of voxel n to changes in parameter values. For position and size estimation, this is largest in voxels along the perimeter of the obstruction; and for attenuation estimation, this is largest in the center voxel of the obstruction.
Combining these interpretations of W T W and P leads to rules of thumb for arranging a sensor network to improve estimation performance. Estimation of an object's position and size improves as more links cross the object's perimeter. Moreover, this holds for each dimension -estimation of the position in dimension d gets better as more links cross the extremes of the object in dimension d. Attenuation may be useful to estimate in aid of object classification; in such a case, estimation gets better as more links cross the center of the object.

IV. SOLUTION TECHNIQUES
It is known that regularization is equivalent to assuming a Bayesian prior on x. In this section, we use that fact to derive a theoretical value for α. We also investigate methods for dealing with aphysical situations in which some entries in y are negative, which can provide negative estimates for the obstructionx.
Given a Bayesian prior f (x), then instead of a regularized ML estimate, we can compute a Maximum A Posteriori (MAP) estimate. Specifically, if y|x ∼ N Wx, σ 2 n I M and x ∼ N (0, Σ) as in (3), then = arg min This is equivalent to the regularized ML solution if αQ = σ 2 n Σ −1 . For (3), considering just the x dimension for simplicity, (18) where c = e −δ/δc = 0.9260 for the parameters in [4]. This is very nearly Note that σ 2 x /σ 2 n is a Signal to Noise Ratio (SNR) and δ c /δ is the number of voxels of separation required for 63% decorrelation. Using the parameter values from [4]   in (19) suggests α ≈ 0.3 for weak fading up to α ≈ 30 for strong fading. The experimental optimum in [4] was α ≈ 5, with good values in the range [2,20]. Thus, rather than empirically searching for a good α, it can be chosen based on statistics of the environment. Now consider the problem of negative entries in y. This can occur due to inadequate calibration of the empty network, or due to large amounts of noise. Fig. 3 shows the search space for a toy problem of 2 pixels side by side, with 2 nodes above and 2 nodes below, leading to 4 links total. If a link has a large enough negative value, the contours of the cost surface of (5) are as in Fig. 3. If we simply compute (a), the solution of (6), and replace the negative value ofx 1 by zero, we get a poor solution, (b). However, if we modify the Gaussian prior on x to be a truncated Gaussian allowing only positive x values, that effectively constrains the search space. The minimum in the constrained space is a much better solution, (c). Alternatively, we could force y to be positive, leading to a comparable solution, (d); however, it is more palatable to constrain the search space than to modify the observed data. The true x is (e).
In order to find the constrained solution without performing a global search, note that the constrained solution occurs where the negative values of the erroneous solution are set to zero. Thus, a possible procedure to approximate method (c) is: 1) Computex via (6).
2) Find the negative elements ofx and omit the corresponding columns of W and of each D d . 3) Repeat (6) using the reduced matrices. 4) Iterate through steps 2-3 if necessary. The iteration is needed because with more than 2 voxels, the cost surface can be more complex than depicted in Fig. 3. The first iteration usually removes most of the negativity, but a few elements in the re-solved (6) may become negative in the next iteration.
To test the utility of this approach, data was taken on Hospital Point at the US Naval Academy in an 18-node roadside surveillance network. The nodes used commercially available ZigBee transceivers with omnidirectional antennas, and the obstruction was a large van. The left side of Fig. 4 depicts the results of method (b), just setting negative values ofx to zero; and the right side shows the proposed method (c) as described above.