Geolocation of Multiple Noncooperative Emitters Using Received Signal Strength: Sparsity, Resolution, and Detectability

TL;DR

This paper explores the limits of locating multiple RF emitters using only RSS data, focusing on resolution, detectability, and the effectiveness of sparse recovery methods in uncharted environments.

Contribution

It introduces a sparse approximation framework for RF emitter localization using RSS data without prior environment mapping, highlighting resolution and detection limits.

Findings

Resolution depends on sensor geometry and noise levels.

Sparse recovery algorithms can effectively identify multiple emitters.

Detectability threshold varies with emitter power and environmental factors.

Abstract

In this paper we investigate the problem of locating multiple non-cooperative radio frequency (RF) emitters using only received signal strength (RSS) data. We assume that the number of emitters is unknown and that individual emitters cannot be distinguished in the RSS data. Moreover, we assume that the environment in which the data has been collected has not been mapped or "fingerprinted" by the prior collection of RSS data. Our primary interest is the limiting resolution that can be obtained by this type of data, and the lowest power emitters that can be detected, as a function of noise level, sensor geometry, and other variables. We formulate the recovery problem as one of sparse approximation or compressed sensing, and investigate an appropriate recovery algorithm for this setting, and use it to illustrate our conclusions. We also include a reconstruction based on sampled data we collected, to illustrate the reasonableness of our parameter choices and conclusions.

Equations

$$d_j=\sum _{i=1}^N{P}_{ij}=\sum _{i=1}^N{p}_i{\left(\frac{r_0}{r_{ij}}\right)}^n$$

$${\mathbf{d}}_0=\mathbf{\Phi }{\mathbf{p}}_0.$$

$$d_j=p_i{\left(\frac{r_0}{r_{ij}}\right)}^n{e}^{\eta R_i}$$

$$d_j=\sum _{i=1}^N{p}_i{\left(\frac{r_0}{r_{ij}}\right)}^n{e}^{\eta R_{ij}}$$

$$\mu (\mathbf{x},\mathbf{y})=\frac{|\mathbf{x}\cdot \mathbf{y}|}{\Vert \mathbf{x}{\Vert }_2\Vert \mathbf{y}{\Vert }_2}.$$

$$\mu (\mathbf{\Phi })=\underset{i\ne j}{\max }\mu ({\Phi }_i,{\Phi }_j).$$

$$\Vert \mathbf{d}-{\mathbf{d}}_k{\Vert }_2\le C\sqrt{({\mu }_0^2+{\sigma }_0^2)}\Vert \mathbf{d}{\Vert }_2$$

$${\mathbf{b}}_1=\frac{{\mathbf{d}}_1}{\Vert {\mathbf{d}}_1{\Vert }_2}\text{ and }{\mathbf{b}}_2=\frac{{\mathbf{d}}_2}{\Vert {\mathbf{d}}_2{\Vert }_2}$$

$$\left[\begin{array}[]{cc}\mathbf{b}_{1}&\mathbf{b}_{2}\end{array}\right]\left[\begin{array}[]{c}p_{1}\\ p_{2}\end{array}\right]=\mathbf{d}$$

$$\mu ({\mathbf{b}}_1,\mathbf{d})>\mu ({\mathbf{b}}_2,\mathbf{d})$$

$$\mathbf{c}\cdot \mathbf{d}>0$$

$$\mathbf{c}={\mathbf{b}}_1-{\mathbf{b}}_2.$$

$$Q=\sum _{j=1}^M{w}_j{e}^{\eta R_j}$$

$$Q^{+}=\sum _{w_j\ge 0}{w}_j{e}^{\eta R_j}\text{ and }{Q}^{-}=\sum _{w_j<0}(-w_j)e^{\eta R_j}.$$

$$H(x)={\int }_{-\infty }^{\infty }F(x+y)g(y)dy.$$

$$H(x)$$

$$\times \left[\frac{1}{y{\sigma }_{-}\sqrt{2\pi }}{e}^{-\frac{(\ln (y)-{\mu }_{-}{)}^2}{2{\sigma }_{-}^2}}\right]dy.$$

$$\Vert \tilde{\mathbf{d}}-\mathbf{d}\Vert >ϵ.$$

$$\Vert \tilde{\mathbf{d}}-\mathbf{d}\Vert$$

$$P>\frac{2ϵ}{\Vert {\mathbf{\Phi }}_i\Vert }.$$

$$\Vert \mathbf{d}-{\mathbf{d}}_0{\Vert }_2^2$$

$${\mu }_0=e^{{\eta }^2{\sigma }_{dB}^2/2}-1,{\sigma }_0^2=e^{{\eta }^2{\sigma }_{dB}^2}(e^{{\eta }^2{\sigma }_{dB}^2}-1).$$

$$E(\Vert \mathbf{d}-{\mathbf{d}}_0{\Vert }_2^2)$$

$$E\left[{\left(\sum _{j=1}^N\frac{p_j}{r_{i,j}^n}{X}_{i,j}\right)}^2\right]$$

$$=\sum _{j,k=1}^N\frac{p_j{p}_k}{r_{i,j}^n{r}_{i,k}^n}E(X_{i,j}{X}_{i,k})$$

$$=\sum _{j=1}^N\frac{p_j^2}{r_{i,j}^{2n}}E(X_{i,j}^2)$$

$$+\sum _{j,k=1,j\ne k}^N\frac{p_j{p}_k}{r_{i,j}^n{r}_{i,k}^n}E(X_{i,j}{X}_{i,k})$$

$$E(X_{i,j}{X}_{i,k})$$

$$E(X_{i,j}^2)$$

$$E(\Vert \mathbf{d}-{\mathbf{d}}_0{\Vert }_2^2)$$

$$=d_0^{2n}\sum _{i=1}^M(({\mu }_0^2+{\sigma }_0^2)\sum _{j=1}^N\frac{p_j^2}{r_{i,j}^{2n}}$$

$$+{\mu }_0^2\sum _{j,k=1,j\ne k}^N\frac{p_j{p}_k}{r_{i,j}^n{r}_{i,k}^n})$$

$$=d_0^{2n}\sum _{i=1}^M\left({\mu }_0^2\sum _{j,k=1}^N\frac{p_j{p}_k}{r_{i,j}^n{r}_{i,k}^n}+{\sigma }_0^2\sum _{j=1}^N\frac{p_j^2}{r_{i,j}^{2n}}\right)$$

$$=d_0^{2n}\sum _{i=1}^M\left({\mu }_0^2{\left(\sum _{j=1}^N\frac{p_j}{r_{i,j}^n}\right)}^2+{\sigma }_0^2\sum _{j=1}^N\frac{p_j^2}{r_{i,j}^{2n}}\right)$$

$$={\mu }_0^2\Vert {\mathbf{d}}_0{\Vert }_2^2+{\sigma }_0^2{d}_0^{2n}\sum _{i=1}^M\sum _{j=1}^N\frac{p_j^2}{r_{i,j}^{2n}}$$

Full text

Key Words: source localization, compressed sensing, detection algorithms, signal mapping, sensor networks.

I Introduction

Locating radio frequency (RF) sources from remotely collected RF data is an essential task in many settings, and is commonly referred to as RF localization, or geolocation. Applications are numerous, for example, the localization of subscribers in cell phone or other wireless networks (indoor or outdoor, see [1, 2, 3, 4, 5]). Localizing transmitters in a cognitive radio network ([6, 7, 8]) allows for the more efficient allocation of network resources, for example, frequency bands. Autonomous vehicles may rely on RF localization to augment navigation [6]. In military applications we may be tasked with geolocating RF transmitters that are non-cooperative or evasive [9, 10, 6, 11, 12]. See [1] for a number of other applications.

A variety of techniques for localizing RF emitters from remote data have been developed. Some techniques use range information deduced from the signal time-of-arrival (TOA), time-difference-of-arrival (TDOA), or the received signal strength (RSS), perhaps collected from multiple sensors at spatially diverse locations. Others, such as angle-of-arrival (AOA), rely on directional information collected from sensors. We may or may not have information about the nature of the RF signals, e.g., emitted power or correlation of measured data from distinct sensors. The accuracy of the resulting position estimates depends on uncertainties in the channel models, sensor placement, and precision of the data collected.

The accuracy of RSS as a method for geolocation is known to suffer from multipath and shadowing effects, but RSS-based localization methods have the advantage that sensor design can be low-complexity; complicated timing, synchronization, or other sophisticated hardware is not needed. Thus, the sensors can be relatively low-cost and low-power. The availability of such sensors is particularly important when many sensors are required, or the sensors are required to be battery-powered (e.g., remote or mobile sensors). Since Received Signal Strength Indicator (RSSI) values are available directly from systems implementing standard communication protocols (see for example, [13]), many WLAN applications do not need any additional hardware to implement RSS-based localization algorithms.

In this work we will focus on the problem of geolocating a “small” but unknown number of non-cooperative RF emitters using RSS measurements from multiple sensors dispersed geographically, or alternatively, from a single sensor on a moving platform, or some combination thereof. In particular, we are interested in methods for estimating the best possible resolution one can obtain from RSS-based location estimates.

Localization of RF sources from RSS data has been considered before ([4, 14, 15, 2, 16]) in a variety of scenarios. Some ([17, 18, 2, 19]) have taken the rough approach we use—a compressed sensing view that exploits spatial sparsity in assuming a small number of emitters are present. But many focus on situations in which the emitters are cooperative [14], or only one emitter is present, or emitters can be distinguished in some manner ([15, 20, 18]) in the data. Since RSS-based localization relies on a propagation model relating signal strength and distance to an emitter(s), RSS-based methods suffer if the signal strength model is inaccurate. Hence some prior work ([2, 17, 21, 14]) assumes that the environment has been “fingerprinted,” that is, sensors have been placed in known locations (“anchors”, [22]), and then empirical measurements taken to map the RF environment. This improves the channel model and accuracy of emitter location estimates. Some methods focus on prediction of lower bounds for the variance of location estimates from RSS data [23, 24, 25].

We consider the problem of geolocating multiple non-cooperative RF emitters in an outdoor environment with a known (at least approximately) channel model. We assume that RSS data is collected by multiple RF sensors (possibly mobile, e.g., mounted on UAVs) whose location(s) are known. We specifically focus, for illustrative purposes, on the case in which multiple elevated sensors are used in an unobstructed open-air scenario, with stationary emitters on the ground, though our analysis is not tied to this arrangement. In particular, we assume that the emitters transmit at a common (known) frequency, such that:

• •

The number of emitters is not known, but is “small,” in a sense to be quantified later.

• •

Emitter signals cannot be distinguished by any characteristic in the time or frequency domain. Thus the RSS data collected by any sensor is the “aggregate” power summed over all emitters.

• •

The RF sensors are “limited” in number and have isotropic sensitivity, so no directional information is available.

From such data we seek to recover the number of emitters, the location of each, and possibly the power at which each emitter transmits.

The unique contributions outlined in this paper are to:

• •

Determine the limiting resolution (ability to distinguish two close emitters) from this type of data, as a function of the data noise/uncertainty level, sensor placement, channel attenuation model, and other relevant physical parameters.

• •

Determine the limiting power threshold for an emitter’s “detectability” (the lowest power emitter than can be detected) as a function of the above-mentioned quantities.

• •

Demonstrate that an appropriate algorithm that makes use of the above assumptions can in fact determine the number of emitters present and their locations.

• •

Use this algorithm to illustrate our conclusions on resolution and emitter detectability.

In Sections II and III below we formulate the problem of locating RF emitters from RSS data as one of finding a sparse solution to an underdetermined linear system of equations, and include an appropriate noise model. We then examine the notion of coherence, which plays a central role in our analysis of resolution and emitter detectability, and then briefly examine an appropriate algorithm for solving the resulting system. In Section IV we analyze the resolution that can be obtained with this type of data, and the limits on emitter detectability. The algorithm developed is used to illustrate our conclusions with computational examples. Finally, in Section V we detail data we collected to validate our channel model parameters, and the geolocation of an emitter from measured RSS data.

II Problem Formulation

II-A System Model

Consider the problem of identifying an unknown number of non-cooperative emitters located on the ground using RSS data from airborne sensors. Specifically, let $\Omega\subset\mathbb{R}^{2}$ and assume the emitters lie at three-dimensional coordinates $(x,y,z)$ with $(x,y)\in\Omega$ and $z=0$.

Although a communication network may employ a number of different frequencies, it is often the case that those of interest for geolocation are relatively few. Several methods exist to detect and classify signals by their frequency content (see, for example, [26, 27]). We will not focus on this aspect of the problem, but rather assume that a frequency (or narrow range of frequencies) of interest has been identified, and that the emitters of interest are operating at these frequencies.

When the number of emitters is sufficiently small, localizing them is a problem well-suited to formulation in the context of compressed sensing, that of finding a sparse solution to a linear system of equations, where “sparse” means that most components in the relevant solution vector are zero (or close to zero.) Specifically, let $S=\cup_{i=1}^{N}\mathbf{r}_{i}$, where $\mathbf{r}_{i}=(x_{i},y_{i})$, be a subset of $N$ distinct points in $\Omega$; these will be the potential locations of any emitters. These points should be chosen to provide a reasonable sampling of the potential locations of any emitters. For example, if $\Omega$ is a rectangle it may be convenient to define $S$ as the nodes on a finely-spaced rectangular grid. It is not essential that emitters be located precisely at any of the $\mathbf{r}_{i}$.

Suppose there are $M$ sensors that measure the RSS at known $(x,y)$ positions $\mathbf{s}_{j}=(a_{j},b_{j})$, each at a fixed altitude $h$ above the $xy$-plane (we take each at the same altitude only for simplicity; the sensors need not be at a single altitude nor directly above $\Omega$.) The $\mathbf{s}_{j}$ may represent distinct sensors, or a single sensor taking data at different points along a path, or some combination thereof. We assume that the sensors’ antennae are isotropic, though more complicated antenna patterns are easily accommodated in the analysis. The distance $r_{ij}$ from the $j$th sensor to the $i$th point in $S$ is $r_{ij}=\sqrt{\|\mathbf{s}_{j}-\mathbf{r}_{i}\|^{2}_{2}+h^{2}}$ where $\|\cdot\|_{2}$ denotes the usual Euclidean norm in the plane.

One common model for the power $P_{ij}$ received at sensor $j$ from an emitter at position $\mathbf{r}_{i}$ is that $P_{ij}=p_{i}(r_{0}/r_{ij})^{n}$ where $p_{i}\geq 0$ is a reference power measured at distance $r_{0}$ from the emitter $i$ and $n$ is the pathloss exponent that governs the attenuation of the signal power as a function of distance; see [10, 28]. In the ideal case the RSS at sensor $j$ from all emitters is then modeled as

[TABLE]

This assumes receiver antennas are equally sensitive, isotropic, and that the emitters are isotropic and incoherent. If no emitter is present at position $\mathbf{r}_{i}$ then $p_{i}=0$, so if few emitters are present we expect $p_{i}>0$ for only a few indices $i$.

We amalgamate the data $d_{j}$ into a column vector $\mathbf{d}_{0}\in\mathbb{R}^{M}$ and express the ideal RSS data (1) in matrix form,

[TABLE]

Here $\mathbf{\Phi}$ is the measurement matrix, an $M\times N$ matrix with known entry $(r_{0}/r_{ij})^{n}$ in row $i$, column $j$. The vector $\mathbf{p}_{0}\in\mathbb{R}^{N}$ has $i$th entry $p_{i}$, the reference power of the emitter at $\mathbf{r}_{i}$, and is sparse if few emitters are present. Note that the entries of $\mathbf{\Phi}$ are known. The $j$th row of $\mathbf{\Phi}$ embodies the data from the sensor at position $\mathbf{s}_{j}$, and the $i$th column corresponds to a potential emitter location $\mathbf{r}_{i}$. We assume that we can measure the quantity $\mathbf{d}_{0}$, the power received by each sensor. The problem of interest is to recover an estimate of $\mathbf{p}_{0}$ from $\mathbf{d}_{0}$ and $\mathbf{\Phi}$. Of course $\mathbf{d}_{0}$ will be corrupted by noise or other error.

II-B Measurement Noise Model

Departure of measured RSS data from the ideal model above is consider at length in, for example, [36]. We assume that data has been suitably processed to eliminate the effects of so-called “fast-fading” and that the error that remains conforms to the standard log-normal noise model. Specifically, if an emitter with reference power $p_{i}$ is present at location $\mathbf{r}_{i}$, the contribution to the data $d_{j}$ collected at the $j$th sensor from this emitter is of the form

[TABLE]

where $\eta=\ln(10)/10$ and $R_{i}$ is a normal random variable with mean [math] and standard deviation $\sigma_{dB}$. Note that $e^{\eta R_{i}}=10^{R_{i}/10}$. Here $\sigma_{dB}$ is the noise level in dB. Values for $\sigma_{dB}$ vary widely depending on the setting, but the application of interest here (outdoors, a relatively open and obstruction-free area) values from $2$ to $5$ dB are common; see [14] or our data in Section V-A.

For multiple emitters we take

[TABLE]

with the additional assumption that the $R_{ij}$ are independent. The model in (3) is valid when the sensors are sufficiently well-separated.

II-C Underdetermined Systems, Coherence, and Sparse Solutions

Let $\mathbf{d}\in\mathbb{R}^{M}$ denote the noisy data vector with components given by (3). Under the assumption that the number of sensors is much smaller than the number of potential emitter locations ($M<<N$), the system $\mathbf{\Phi}\mathbf{p}=\mathbf{d}$ to be solved for $\mathbf{p}$ (an estimate of $\mathbf{p}_{0}$) is underdetermined, and so almost certainly possesses infinitely many solutions. However, as noted we will make the reasonable assumption that there are few emitters, so that the solution vector $\mathbf{p}_{0}$ is sparse. More specifically, a vector $\mathbf{p}$ is said to be $k$-sparse if $\mathbf{p}$ has at most $k$ nonzero components. Under the assumption that $\mathbf{p}_{0}$ is $k$-sparse for sufficiently small $k$, it is highly likely that a physically relevant solution can be found, although the existence of a unique sparse solution and the ease with which it can be found depend on the measurement matrix $\mathbf{\Phi}$.

One property that $\mathbf{\Phi}$ can possess that leads to favorable recovery results is that of low “mutual coherence.” First, the coherence of vectors $\mathbf{x},\mathbf{y}\in\mathbb{R}^{N}$ is the quantity

[TABLE]

The Cauchy-Schwarz inequality shows that $0\leq\mu(\mathbf{x},\mathbf{y})\leq 1$, with $\mu(\mathbf{x},\mathbf{y})=0$ when $\mathbf{x}$ and $\mathbf{y}$ are orthogonal and $\mu(\mathbf{x},\mathbf{y})=1$ when one vector is a scalar multiple of another. The mutual coherence of an $M\times N$ matrix $\mathbf{\Phi}$ with columns $\mathbf{\Phi}_{i}$ is the quantity

[TABLE]

Again, $0\leq\mu(\mathbf{\Phi})\leq 1$. If $\mu(\mathbf{\Phi})=1$ then two or more distinct columns of $\mathbf{\Phi}$ are scalar multiples of each other, while $\mu(\mathbf{\Phi})=0$ means $\mathbf{\Phi}$ is an orthogonal matrix, which is impossible in the present situation since $M<N$.

Low coherence matrices are desirable when seeking sparse solutions to a linear system $\mathbf{\Phi}\mathbf{p}=\mathbf{d}$. It can be shown that if $\mu(\mathbf{\Phi})<1/(2k-1)$ then any $k$-sparse solution $\mathbf{p}$ is unique and many compressed sensing algorithms will converge to this solution (see Section 5.1 of [29]). Low mutual coherence also leads to more favorable bounds on the error in the presence of noisy data ([30]).

Unfortunately, for the localization problem described above, low mutual coherence will not hold for any realistic sensor configuration. First, our measurement matrix has entirely positive entries, so no cancelation occurs in the dot product of columns of $\mathbf{\Phi}$; as a result, the pairwise coherence for any two columns is likely to be larger than for a matrix with mixed sign entries. Also, if potential emitter locations $(x_{i},y_{i})$ and $(x_{j},y_{j})$ are closely spaced, then the $i$th and $j$th columns $\mathbf{\Phi}_{i}$ and $\mathbf{\Phi}_{j}$ of the measurement matrix will be nearly identical, and so have high pairwise coherence. Thus, if we work on a fine grid (to obtain higher source resolution) we confront measurement matrices with high mutual coherence. This presents a challenge for the finding the correct sparse solution.

III Algorithm for Sparse Solutions

In this section we briefly detail an algorithm appropriate for finding sparse solutions to the problem at hand. Our goal is not so much to focus on this specific algorithm, but to use it to gain insight into the ill-posedness of this inverse problem, and provide examples that illustrate the analysis for resolution and clearance.

III-A BLOOMP

Finding the sparsest solution to a linear system of equations is, in general, computationally intractable, even if a sparse solution is known to exist [31]. However, a number of efficient computational approaches have been devised that, under the right conditions, find such a sparse solution with high probability. In this section we justify use the algorithm “Band-excluded Locally Optimized Orthogonal Matching Pursuit” (BLOOMP, see [32]) for the present problem, and include an illustrative computational example.

Briefly, the BLOOMP algorithm is a modification of Orthogonal Matching Pursuit (OMP). OMP is a “greedy” algorithm that iteratively builds up a sparse solution to $\mathbf{\Phi}\mathbf{p}=\mathbf{d}$ one nonzero component at a time. Let $\mathbf{p}^{0}=\mathbf{0}$ denote our initial guess at a solution, $\mathbf{p}^{k}$ the $k$th iterate (at most $k$-sparse) in OMP, and $S^{k}=\{i:\mathbf{p}^{k}_{i}\neq 0\}$; $S^{k}$ is called the support set of $\mathbf{p}^{k}$. The set $S^{k}$ indexes those columns of $\mathbf{\Phi}$ that are being used to synthesize the data $\mathbf{d}$. OMP constructs $\mathbf{p}^{k+1}$ by augmenting the support $S^{k}$ with a new index $i_{k}$ chosen so that the residual $\|\mathbf{\Phi}\mathbf{p}^{k+1}-\mathbf{d}\|_{2}$ is minimized. This continues until a maximum sparsity bound or a termination criterion is met. One common stopping criterion takes the form $\|\mathbf{\Phi}\mathbf{p}^{k}-\mathbf{d}\|_{2}\leq C\epsilon$ where $\epsilon$ is comparable to the expected noise level in the data as measured in the Euclidean norm and $C\approx 1$; see [29]. We say more on this in our specific application below in Section III-B and Appendix VIII.

A drawback of OMP is that once an index has been added to the support set $S^{k}$, it is never removed at a later iteration, so sub-optimal early choices cannot be undone. Many modifications to OMP have been proposed to overcome this problem. We have adopted one such modification, BLOOMP [32], because it is particularly suited to “high-coherence” measurement matrices. Like OMP, BLOOMP builds a sparse solution by adding one index at each iteration to the potential support set. In our application this means adding one estimated emitter at each iteration. However, in the BLOOMP algorithm the column in $\mathbf{\Phi}$ corresponding to the emitter added at a given iteration cannot have high coherence with any column of $\mathbf{\Phi}$ corresponding to previously added emitters. Physically, the next estimated emitter location cannot be too close to those already determined to be present—this is the “band exclusion” modification of OMP. Moreover, at each iteration the emitters currently estimated to be present are subject to local adjustments in location and power to better fit the data; this is the “local” optimization portion of the algorithm. The authors in [32] show that in situations such as these—high coherence matrices, but in which the correct solution index support corresponds to columns with lower pairwise coherence, such as well-separated emitters on a finely-spaced grid—the BLOOMP modifications increase the probability of recovering the correct solution support indices, or in our case, the correct emitter number and location(s). We also add a constraint to the algorithm to require that at each iteration the emitter power estimates must remain nonnegative.

III-B Recovery Example

To illustrate, let $\Omega$ be the $50\times 50$ meter region $\{(x,y);0\leq x,y\leq 50\}$ and consider a $50\times 50$ rectangular grid for potential emitter locations, of the form $(x_{i},y_{j})$ where $x_{i}=(i-0.5),y_{j}=(j-0.5)$ for $1\leq i,j\leq 50$, so here $N=50^{2}=2500$. In many settings it is the case that randomness in the construction of the measurement matrix is an asset in using sparsity or compressed sensing recovery algorithms [33]. We thus consider $M=30$ RSS data points collected from an emitter on a “meandering” path above $\Omega$, at altitude $h=10$ meters. The sensor locations are displayed as crosses in Fig. 1. We use pathloss exponent $n=3.5$ in equation (1) (assumed known for now) and noise level $\sigma_{dB}=3$ dB in equation (3). See Section V for data that supports these parameter choices, and for a recovery from measured data.

Three emitters with unit power at reference distance $r_{0}=1$ meter are placed at $(24.0,41.0)$, $(19.3,20.1)$, and $(36.4,12.8)$ (hereafter referred to as emitters $1,2,$ and $3$). Note that these are not themselves grid points; nonetheless, one would hope to recover emitter estimates that correspond to nearby grid points. We then simulate noisy data $\mathbf{d}$ using equation (3) and perform a reconstruction from $\mathbf{d}$ using the BLOOMP algorithm, to recover an estimate of the emitter number, location(s), and power(s). This process of generating noise and reconstructing is repeated $500$ times, each with a different noise realization. The number of emitters is not assumed a priori.

One can show (see Appendix VIII) that for a modest noise level $\sigma_{dB}\leq 5\textrm{ dB}$ the expected value of $\|\mathbf{d}-\mathbf{d}_{0}\|^{2}_{2}$ is bounded by and comparable to the quantity $\epsilon=(\mu_{0}^{2}+\sigma_{0}^{2})\|\mathbf{d}_{0}\|_{2}^{2}$ where $\mu_{0}=e^{\eta^{2}\sigma_{dB}^{2}/2}-1$ and $\sigma_{0}^{2}=e^{\eta^{2}\sigma_{dB}^{2}}(e^{\eta^{2}\sigma_{dB}^{2}}-1)$ (recall $\eta=\ln(10)/10$). Of course we expect the noiseless data $\mathbf{d}_{0}$ is unknown, but the noisy data $\mathbf{d}$ provides a reasonable estimate. We thus terminate the iteration when the fit to the data is comparable to (or a bit smaller than) this noise level, specifically, when

[TABLE]

where $\mathbf{d}_{k}=\mathbf{\Phi}\mathbf{p}^{k}$ denotes the estimated data at the $k$th iteration of BLOOMP and $C$ is a constant less than $1$ (we use $C=1/4$). For high noise levels the random variable $\|\mathbf{d}-\mathbf{d}_{k}\|_{2}$ is more highly skewed to the right, and so $E(\|\mathbf{d}-\mathbf{d}_{k}\|_{2})$ may be somewhat smaller than $\sqrt{\|\mathbf{d}-\mathbf{d}_{k}\|_{2}^{2}}$, hence a value of $C$ somewhat less than $1$ can be helpful to prevent the iterative algorithm from terminating too early.

The results of these 500 simulated cases are shown in Fig. 1. The image is an average of the recovered power at each grid location, coded so [math] recovered power is white, $1$ or higher is black. The sensor locations are illustrated as crosses and the true position of each emitter is represented by a star. The average estimated power for each of the three emitters is $1.7424,0.8830,$ and $1.1712$ for emitters $1,2,$ and $3$.

The grey areas indicating positive power recovery clustered around the three emitters, which are reasonably well resolved. The spread of each cluster gives an indication of the resolution one can achieve with this sensor configuration and noise level. An analysis of this resolution is the focus of the next section.

The pathloss exponent $n$ would appear to be a rather critical value in estimating the number and position of the emitters, but we find that this is not the case. Specifically, an incorrect pathloss exponent has little effect on the recovery of the emitter count and locations, but does significantly affect the estimated power of each emitter. As an illustration, in Fig. 2 is shown a recovery with exactly the same parameters as Fig. 1, but with the (erroneous) assumption of a pathloss exponent of $2.5$ (whereas $n=3.5$ was used to generate the data). The average estimated power is $0.1368,0.0595,$ and $0.0612$ for emitters $1,2,$ and $3$, respectively, considerably off from the correct values of $1$ for each. Nonetheless, the number and location are quite accurate.

IV Analysis of Resolution and Detection Limits

The goal in this section is to develop a method for quantifying the local resolution one can obtain at any fixed potential emitter location from RSS data for a given noise level and sensor configuration, and to provide a bound on the weakest emitters that can be reliably detected.

IV-A Resolution Analysis

Suppose an emitter lies at one of two potential locations, say $\mathbf{q}_{1}=(x_{1},y_{1})$ or $\mathbf{q}_{2}=(x_{2},y_{2})$. We collect noisy RSS data from $M$ sensors. The goal is to determine at which location the emitter actually lies, with sufficiently high probability (to be specified). If this can be done we will say the two potential locations are “resolvable.”

Let $\mathbf{d}_{k}\in\mathbb{R}^{M}$ denote the noiseless RSS data we would collect from an emitter at location $\mathbf{q}_{k}$, where $k=1$ or $k=2$. This data vector is assumed to obey the model (1), with a single nonzero power location. For convenience we define normalized data vectors

[TABLE]

so $\|\mathbf{b}_{k}\|_{2}=1$ for $k=1,2$. Note that the reference power $p_{k}$ will not matter in either case.

Suppose we collect noisy data $\mathbf{d}\in\mathbb{R}^{M}$ from the sensors, stemming from an emitter at location $\mathbf{q}_{1}$; the components of $\mathbf{d}$ are given by (3) (with only a single nonzero summand). The goal is to use $\mathbf{d}$ to correctly assign the emitter to location $\mathbf{q}_{1}$, as opposed to $\mathbf{q}_{2}$. If we formulate this as a compressed sensing problem as above in equation (2) then we obtain linear system

[TABLE]

(equivalently, $p_{1}\mathbf{b}_{1}+p_{2}\mathbf{b}_{2}=\mathbf{d}$) in which the sensing matrix $\mathbf{\Phi}$ is $M\times 2$ with unit norm columns $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$. We seek a $1$-sparse solution to this system. In this very simple case OMP or any standard sparse solver (e.g. basis pursuit) will provide a $1$-sparse solution consisting of a multiple of that column of $\mathbf{\Phi}$ which has the highest coherence with the data $\mathbf{d}$, with power estimate $\tilde{p}_{k}=\mathbf{b}_{k}\cdot\mathbf{d}$ for either $k=1$ or $k=2$. That is, the emitter is correctly assigned to location $\mathbf{b}_{1}$ if

[TABLE]

and incorrectly to location $\mathbf{p}_{2}$ otherwise. Condition (8) is quite natural—the emitter is assigned to a location according to which vector $\mathbf{b}_{1}$ or $\mathbf{b}_{2}$ best matches the collected data $\mathbf{d}$ after optimal scaling for power. This notion of resolution is not wedded to a compressed sensing approach to the problem, nor any particular algorithm.

Equation (8) is equivalent to

[TABLE]

where

[TABLE]

Equations (9) and (10) can be written equivalently as $Q>0$ where

[TABLE]

with $w_{j}=c_{j}/r_{1j}^{n}$ and where $c_{j}$ denotes the $j$th component of $\mathbf{c}$. For a given sensor configuration the $w_{j}$ are known. We want to compute the probability $P(Q>0)$, so that we correctly assign the emitter to location $\mathbf{q}_{1}$. It should be noted that we will have $0.5\leq P(Q>0)\leq 1$, with $P(Q>0)=1$ as the best case—the emitters are certainly resolvable—and $P(Q>0)=0.5$ as the worst case, in which resolving the emitter locations becomes a “coin toss.”

The random variable $Q$ is a signed linear combination of lognormal random variables (the $w_{i}$ are generally of mixed sign). The next section is devoted to accurately approximating the probability $P(Q>0)$ in an easily computable fashion.

IV-B Approximating a Signed Sum of Lognormal Random Variables

Though $Q$ has coefficients of mixed sign, we first consider the case in which all coefficients are positive. The distribution of such a sum of lognormals is a well-studied problem, though such a sum has no closed-form density function. However, it has long been noted that such a sum is itself approximately lognormal, and so can be characterized as being of the form $e^{N(\mu,\sigma^{2})}$ for suitable $\mu$ and $\sigma$ (Note that $\mu$ stands for the mean of the noise distribution here, not mutual coherence).

In [34] the authors provide a simple and effective method for fitting $\mu$ and $\sigma$ to such a sum. The individual lognormals in the sums they consider are of the form $e^{N(\mu_{i},\sigma_{i}^{2})}$ with varying $\mu_{i}$ and $\sigma_{i}$, and are assumed independent. For a linear combination of the form (11) with weights $w_{i}$ that are positive, the weighted sum in $Q$ is easily adapted to this setting, by absorbing the $w_{i}$ into the $R_{i}$ (we can shift the mean of $R_{i}$ by $\ln(w_{i})$). If we split the sum defining $Q$ into a piece with positive weights and a piece with negative weights, we can write $Q=Q^{+}-Q^{-}$ where

[TABLE]

The method of [34] provides a lognormal random variable approximation for $Q^{+}$ in the form $e^{R}$ where $R=N(\mu_{+},\sigma_{+})$, by determining an appropriate mean and variance $\mu_{+}$ and $\sigma_{+}^{2}$. A similar approximation is made to obtain $\mu_{-}$ and $\sigma_{-}$ for $Q^{-}$.

The probability density function (pdf) and cumulative density function (cdf) for the lognormal random variable are well-known. Moreover, if a random variable $X$ has cdf $F(x)$ and random variable $Y$ has pdf $g(x)$ then the cumulative distribution function $H(x)$ of $X-Y$ is given by

[TABLE]

Then, for example, $P(X-Y>0)$ is given by $1-H(0)$. In the present case the cdf $H$ for $Q=Q^{+}-Q^{-}$ can be expressed as

[TABLE]

The $\max(0,-x)$ lower limit cuts off the integral as soon as the cdf or pdf of either random variable equals zero. The value we are interested in is $P(Q>0)=1-H(0)$, and this can be computed easily from (13).

The overall procedure is as follows: Given potential emitter locations $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$, we compute $\mathbf{c}$ as in (10) and set $w_{j}=c_{j}/r_{1j}^{n}$ with $r_{1j}$ as the distance from location $\mathbf{q}_{1}$ to the $j$th sensor. We then use the procedure in [34]) to estimate $\mu_{+},\sigma_{+},\mu_{-},$ and $\sigma_{-}$ for $Q^{+}$ and $Q^{-}$ and compute $P(Q>0)$ using (13). If $P(Q>0)$ exceeds some threshold probability $p_{min}$ we will say the emitter location $\mathbf{q}_{1}$ is resolvable from location $\mathbf{q}_{2}$.

To illustrate the accuracy of the approximation, Fig. 3 shows the quantity $P(Q>0)$ computed by this procedure versus the simulated probability of correctly resolving the emitters locations for a variety of sensor counts and noise levels. In each base we use $\mathbf{q}_{1}=(24.5,41.5)$ and $\mathbf{q}_{2}=(19.5,20.5)$ with sensors at random $(x,y)$ locations in $0<x,y<50$ at altitude $h=10$. We generate $10^{4}$ realizations of synthetic noisy data $\mathbf{d}$ for a sensor at location $\mathbf{q}_{1}$ and assign it to location $\mathbf{q}_{1}$ if $\mu(\mathbf{d},\mathbf{b}_{1})>\mu(\mathbf{d},\mathbf{b}_{2})$, location $\mathbf{q}_{2}$ otherwise. The pathloss exponent is $3.5$.

As an example of how this can be used to quantify local resolution, consider the three-emitter configuration of Fig. 1, with the same noise level and other parameters. What local resolution might we expect near the emitter at location $(24,41)$? Let $\mathbf{p}_{1}=(24,41)$ and $\mathbf{p}_{2}=(x,y)$ for $0<x,y<50$, so that $P(Q>0)$ as computed above is a function of $(x,y)$. In Fig. 4 we show a contour plot of this function.

The red oval delineates, for comparison, the Cramer-Rao lower bounds on the uncertainty in estimating the location of the emitter (discussed below).

To illustrate the validity of the resolution analysis, in Fig. 5 is shown a situation similar to that of Fig. 1, but in which the emitter at position $(19.3,20.1)$ has been moved to $(19.0,36.0)$, which is only 7 meters away from the emitter at position at $(24,41)$. The newly moved emitter lies outside the Cramer-Rao bounds, on about the $P(Q>0)=0.85$ contour. The emitters are not as reliably resolved.

This analysis makes it clear that, for a given noise level (and other parameters) the resolution obtainable with RSS data is limited, and can be quantified. In particular, in a compressed sensing approach there are little improvements in resolution to be obtained by using too fine of a grid.

IV-C Comparison to Cramer-Rao Bounds

Other authors (e.g., [24, 25, 23]) have examined statistical bounds, for example, Cramer-Rao bounds, on the minimum variance that can obtained by using RSS data to estimate the distance to or position of an emitter. Such a bound provides a natural way to quantify resolution. However, as noted in [24, 25], the Cramer-Rao bounds in this setting cannot be attained by any unbiased estimator, and so are too optimistic. (It should also be noted that our estimates are almost certainly biased.)

To illustrate and compare with the current analysis, we consider a single emitter of unknown power $p_{0}$ at true location $(22,41)$, altitude zero, with the $30$ sensor locations as used in Figs. 1, 4, and 2, pathloss exponent $3.5$, and noise level $3$ dB. Following the computations of Section 3.2 in [25] we establish a Cramer-Rao lower bound on the minimum covariance of any unbiased estimator of the emitter location and power. The red elliptical region in Fig. 4 is a 95 percent confidence region with respect to the spatial variables for an emitter with $p_{0}=1$, though the bounds do not depend on the unknown power $p_{0}$. The lower bound on the variance of any unbiased estimate of $p_{0}$ is $0.16$.

IV-D Detectability and Clearance

In this section we consider the problem of when we can be reasonably certain that we have detected all the emitters above a given power threshold in a region of interest; this could be the entire region $\Omega$ or some subregion thereof.

For a given configuration of $M$ sensors, $N$ potential emitter locations and corresponding $M\times N$ measurement matrix $\mathbf{\Phi}$, suppose that $\mathbf{p}_{0}\in\mathbb{R}^{N}$ embodies the true emitter power vector. The noise-free data $\mathbf{d}_{0}\in\mathbb{R}^{M}$ is given by (2); let $\mathbf{d}\in\mathbb{R}^{M}$ be the collected (noisy) data vector. Suppose that $\mathbf{p}_{r}$ is an estimate of $\mathbf{p}_{0}$ based on the data $\mathbf{d}$, computed using BLOOMP or any other recovery algorithm. We assume, however, that the algorithm produces an estimate $\mathbf{p}_{r}$ for which an error bound of the form $\|\mathbf{\Phi}\mathbf{p}_{r}-\mathbf{d}\|\leq\epsilon$ holds, for some tolerance $\epsilon$, where $\|\cdot\|$ can denote any norm, e.g., the $L^{2}$ or supremum norm. Typically $\epsilon$ is comparable to the expected noise level in the data in the appropriate norm.

Now suppose that a single additional emitter were present at location $\mathbf{r}_{i}$, with power $P$. Let $\tilde{\mathbf{p}}=\mathbf{p}_{r}+P\mathbf{e}_{i}$ denote resulting power vector ($\mathbf{e}_{i}$ is the $i$th standard basis vector). This would yield data $\tilde{\mathbf{d}}=\mathbf{\Phi}\tilde{\mathbf{p}}=\mathbf{\Phi}\mathbf{p}_{r}+P\mathbf{\Phi}_{i}.$ We will consider the additional emitter at $\mathbf{r}_{i}$ to be detectible if

[TABLE]

That is, the presence of this additional emitter would yield reconstructed data $\tilde{\mathbf{d}}$ that is inconsistent with the measured data at the given tolerance level. But we do not require that the reconstructed emitter power configuration $\mathbf{p}_{r}$ be accurate, in that $\|\mathbf{p}_{r}-\mathbf{p}_{0}\|$ need not be small.

The value of $P$ that assures $\|\tilde{\mathbf{d}}-\mathbf{d}\|>\epsilon$ holds can be estimated. We have, using the reverse triangle inequality

[TABLE]

Inequality (14) must hold if $P\|\mathbf{\Phi}_{i}\|-\epsilon>\epsilon$ or

[TABLE]

The threshold on the right in (15) depends on the precision to which we fit the measured data, i.e., the noise level in the data, the norm we use, and on $\mathbf{\Phi}$. By taking the maximum of the right side of (15) over all locations $\mathbf{r}_{i}$ in a given region $\Omega^{\prime}\subseteq\Omega$ we obtain a threshold of the weakest emitters that can be reliably identified in $\Omega^{\prime}$. If a lower threshold is desired, it would be necessary to alter the number and/or placement of sensors. Inequality (15) quantifies what is required. Of course the estimates leading to (15) are likely pessimistic—an emitter may well be detected below this power threshold—but it does provide a rough lower bound for emitter detectability.

To illustrate, again consider the setting of Fig. 1. Let us consider the power threshold for detectability of the emitter at location $(19.3,20.1)$. The closest grid location is $\mathbf{r}_{1020}=(19.5,20.5)$ (that is, index location $i=1020$ in our indexing scheme). We iterate BLOOMP until $\|\tilde{\mathbf{d}}-\mathbf{d}\|_{2}\leq 2\times 10^{-4}$ and compute $\|\mathbf{\Phi}_{1020}\|_{2}\approx 3.56\times 10^{-4}$, leading to a power bound $P\approx 1.1$ for the emitter in this location. As is obvious in Fig. 1, the emitter is clearly detectable at power level $1$. However, under the same conditions but with power level $0.5$ the result is as shown in Fig. 6. At power level $0.25$ the emitter becomes essentially invisible.

V Illustration With Measured Data

In this section we briefly detail an experiment we, with the aid of our students, performed to collect actual RSS data under the conditions that were only simulated above. Our goal here is not to reproduce the resolution or clearance analysis with experimentation, but rather to estimate realistic noise and pathloss parameters under relatively ideal conditions. However, we do perform a reconstruction for a single emitter, and illustrate the effect of using an erroneous pathloss exponent.

V-A Measurement of RSS in Open Air

An experiment to collect RSS data from a single emitter using $15$ sensors was conducted in the open-air on a flat grass-covered field of $90\times 120$ meters with no overhead obstructions. The transmitter was placed at location $(12.4,17.5)$ meters relative to an origin on a cartesian grid, at a height of $70$ cm. Fifteen different receivers were scattered within a $50\times 50$ meter square area to collect RSS samples at 15 different positions $\mathbf{s}_{j}$, corresponding to 15 different distances $r_{1j}$. The height of the receivers was $50$ cm. The locations of the receivers and emitter are plotted in Fig. 7. This isn’t precisely the configuration we simulated, but the analysis is easily adapted to any emitter/sensor geometry.

The transmitter emitted a continuous-wave, unmodulated signal, centered at 925 MHz (in the ISM band) using a Software Defined Radio (SDR) transceiver (USRP E100, Ettus Research). An omnidirectional vertical dipole antenna was used for the transmitter (VERT900, Ettus Research). The transmitted signal was sampled at each sensor position at a rate of $1.152$ Msamples/s, for a duration of one second, using an SDR radio (receiver only) device with a USB interface (R820T NESDR Mini, Noo Electric). The RTL-SDR has the capability to tune over the range 25 MHz to 1.75 GHz, producing raw, 8-bit IQ data samples, at a programable, baseband sampling rate of up to 2.8MHz [35]. However, the data acquisition sampling rate was set lower to ensure the accuracy of the rate. The gain was set to $32.8$ for each of the receivers, which was tuned so that the receiver closest to the transmitter ($\approx 6.5$ meters away) did not experience saturation. Without automatic gain control, we found the useful dynamic range of the RTL-SDR is around $45$ dB. The receivers used an omnidirectional vertical dipole antenna, approximately $14$ cm in length with an MCX connection.

The raw IQ data were processed using the procedures recommended in [36]. The RSS was calculated by first applying a Chebyshev Type I IIR filter of order $1$ to remove most of the fast-fading variations The RSS values are normalized to the value received at the sensor with the shortest distance to the transmitter (about six meters). A least-squares fit to the log-normal distance trend is used to estimate the path-length exponent, $n\approx 3.45$ for our data. The standard deviation of the log normalized uncertainty term (long-term fading uncertainty) was computed from the variation from the fitted data, $\sigma_{dB}\approx 1.86$ dB. The decimated RSS values and fit are plotted in Figure 8.

V-B A Sample Reconstruction from Data

A reconstruction is performed using the BLOOMP algorithm to reconstruct an estimate of the transmitter power vector, $\mathbf{p}$, from the measured RSS vector, $\mathbf{d}$. The index value of the power vector, $p_{i}$, which are not estimated to be zero will indicate the location of a detected emitter. Although our measured data is known to be received from only one emitter, we allow the algorithm to iterate as many as 12 times, corresponding to a reconstruction with the potential to predict as many as 12 emitter locations. The measurement matrix, $\mathbf{\Phi}$, is calculated by assuming a one-slope propagation model ([36] and our equation (3)), with the reference distance set to 1, the path-loss coefficient set to 3.45, the received power at the reference distance, $K$, is set 1 and all of the RSS values are normalized to this reference power. Band exclusion is applied to the modified BLOOMP algorithm according to [32] with a exclusion parameter set to 0.98. The BLOOMP reconstruction algorithm is terminated after one iteration when the residuals are less than $0.5$ times the expected uncertainty using a log normalized standard deviation of 2dB. In Fig 2, the estimated location of the emitter is plotted as the magenta square. The estimate falls on the closest grid point at (12.50m,15.50m) which is 2.06 meters away from the position of the true emitter’s location (ro = 12.41m, 17.56). It should be noted that the true emitter is not placed on a grid point.

VI Conclusion

We have formulated the problem of geolocating multiple non-cooperative RF emitters in a given region using low-capability sensors as a problem in compressed sensing, and use this formulation to develop methods for examining the limits on resolution and emitter detectability as a function of the data noise level, sensor number and configuration, as well as other relevant variables, for example, the pathloss exponent. We have also implemented an algorithm suitable for actually recovering emitter number and location from simulated data. We also demonstrate the recovery of a single emitter using measured data.

Several natural extensions and refinements of this technique suggest themselves. The model can be easily adapted to directional sensor antennae, and sensors (or emitters) at nonconstant altitude. Also of interest, but more challenging, is the problem of locating anisotropic, intermittent, or moving emitters, and operating in an environment in which sensor positions’ themselves are not known and must be estimated.

VII Acknowledgements

The authors would like to thank the many Rose-Hulman undergraduate and masters students who have worked with us on this project (and continue to do so). Dr. Walter would also like to thank the many colleagues at the Air Force Research Laboratory Sensors Directorate who have supported her.

VIII Appendix A

As noted in Section III-B, we iterate the BLOOMP algorithm until the fit squared residual is comparable to $E(\|\mathbf{d}-\mathbf{d_{0}}\|_{2}^{2})$. The latter quantity can be estimated from $\mathbf{d}$ (the measured data) and the noise level $\sigma_{dB}$.

From the noise-free model (1) and (3) we compute

[TABLE]

where $X_{i,j}=e^{\eta R_{i,j}/10}-1$ with $R_{i,j}$ normal with mean zero, variance $\sigma^{2}_{dB}$. The random variable $X_{i,j}$ is lognormal with mean and variance given by

[TABLE]

Since the expected value is linear,

[TABLE]

A little algebra shows that

[TABLE]

Since the $X_{i,j}$ are independent we have

[TABLE]

Then some mundane algebra shows that

[TABLE]

This provides the basis for the termination criterion (6) (replacing $\mathbf{d}_{0}$ with $\mathbf{d}$.)

Note however that $E(\|\mathbf{d}-\mathbf{d}_{0}\|_{2}^{2})\geq E(\|\mathbf{d}-\mathbf{d}_{0}\|_{2})^{2}$, so the termination criterion (6) may result in under-fitting the data. When $\sigma_{dB}\leq 5$ the quantities $E(\|\mathbf{d}-\mathbf{d}_{0}\|_{2}^{2})$ and $E(\|\mathbf{d}-\mathbf{d}_{0}\|_{2})^{2}$ are comparable in magnitude, but for larger noise levels the random variable $\|\mathbf{d}-\mathbf{d}_{0}\|_{2}$ is skewed heavily higher, to the right. In such a case a smaller value of $C$ in (6) is appropriate.