Wideband passive source localization

This paper develops a mathematical method for determining the locations of multiple transmitters from passive measurements of the signals at two or more receivers. The method applies to the case of emitters transmitting either wideband or narrowband signals.


Introduction
This paper develops a mathematical method for determining the locations of multiple transmitters from passive measurements of the signals at two or more receivers. A typical geometry is shown in figure 1. In this problem, the emitted signals are unknown, and, in particular, the time at which they originate is unknown, so ordinary echolocation and triangulation methods cannot be used.
We consider the case of emitters that are transmitting acoustic or electromagnetic energy for the purposes of sensing or communication. The receivers of such emissions must be able to identify emissions from a particular source. For example, different broadcast radio stations use different frequency bands. Similarly, a cell-phone receiver must be able to identify the signals due to a particular call; in the US, this is done by means of the code division multiple access (CDMA) protocol, in which different calls use the same frequency band but are assigned to orthogonal codes.
On the other hand, there are some emitters that interfere with others, either accidentally, as in the case of emissions from an electric motor, or intentionally, as in the case of jamming. In this paper we do not consider the case of multiple interfering emitters.
The problem of determining locations of emitters has been previously studied; see e.g. [1,3,6,7,9,12,16,18] and references in these works. Most of this work assumes the presence of a single point-like transmitter. The problem for multiple transmitters has been addressed in [10,11] which approach the problem from an imaging point of view. A detection-theoretic approach for multiple transmitters has been recently developed by the authors of [13]. The problem of passive source localization is closely connected to the problem of passive imaging, which has been studied in e.g. [5,17,20].
Efforts to locate unknown sources require decisions regarding two issues: (i) What exactly is being imaged? Much of the previous work assumes that only a single point-like source is present. Other previous work, such as [10,11], assumes a statistical distribution of sources that is incoherent, i.e. delta-correlated. In contrast, we seek to form an image of a non-statistical distribution of sources. (ii) How should we handle the fact that the emitted waveform is unknown? In the case of a single source, the unknown waveform can be eliminated. Much of the work dealing with multiple emitters simply ignores the waveform and instead focuses on the phase. In contrast, we exploit the (approximate) waveform orthogonality that the spectrum users themselves require.
One of the challenges in passive sensing is that the natural data consists of receiver signal cross-correlations. These cross-correlations are quadratic quantities, which means that the usual linear methods do not apply.
In this work, we assume that the transmitters are stationary and that at least one of the receivers is moving. The fact that one receiver is moving enables us to compute time-frequency transforms of the received signals. Specifically, we use the data from two receivers to compute the receiver cross-ambiguity function, which we first relate to the emitter auto-ambiguity function and then backproject to form a spatial image of the emitter locations. This technique combines the known efficacy of time-difference-of-arrival (TDOA) and frequency-differenceof-arrival (FDOA) methods in passive detection, with the ability of synthetic-aperture-like backprojection imaging [4] to coherently focus a weak signal. Signals are received on two receivers, labelled γ 1 (moving) and γ 2 (stationary), and the goal is to find the locations of the sources, here shown in the lower left corner.
We use wideband ambiguity functions [14], to allow for long integration times and emitter waveforms with high frequency diversity. Thus, our method applies to the case of emitters transmitting either wideband or narrowband signals. Relevant work on the wideband ambiguity function includes [14,19,20].
The method developed in this paper does not require the introduction of an emitter density function, and does not rely on the theory of Fourier integral operators. Rather, the method exploits a relation similar to that of [2] to derive a wideband version of a Moyal-like identity [19].

Mathematical model for the data
We consider the case in which the receivers, which are assumed to be pointlike, are moving along paths γ m (t), m = 1, 2, . . .. For example, for a stationary receiver located at z, we have simply γ m (t) = z.
We denote the waveform emitted from a transmitter at location y by p y (t). For example, if the scene consists of a single cell phone at location y 0 emitting the CDMA signal p y 0 (t), then we could take p y (t) = δ(y − y 0 )p y 0 (t). If there is no source at a location y, then p y (t) = 0. We assume that the transmitters are isotropic, so that the signal from each transmitter is received at each receiver. We also assume that the emitted waveforms p y (t) are smooth in t, and we assume that the emissions from different sources (i.e. from different locations) can be disentangled in a sense that will be made more precise below.

Background on the wave equation
We assume that the field u emanating from the source satisfies the scalar wave equation The free-space Green's function is which satisfies

The model for received data
The solution of (1) is obtained by convolving the right side of (1) with the Green's function (2): The data d m (t) measured by receiver m is u(t, γ m (t)): A receiver antenna beam pattern could easily be included in (5); this would limit the region of integration in y at each t. In (5), by assuming that time is the same in all reference frames, we are neglecting any effects of special relativity.

Variable counts and time scales
We note that in the case in which the sources lie on a known surface, the unknown source p y (t) is a function of three variables, namely y 1 , y 2 , and t, whereas the signal d m (t) is a function of only one variable. However, this problem involves multiple time scales. First, the speed of wave propagation is 3 · 10 9 m s −1 , which is much faster than the speed at which the receivers move, which is typically subsonic, less than about 300 m s −1 . Consequently we introduce 'fast' and 'slow' time scales with the windowing technique below.
Second, the emitted waveforms are likely to be highly oscillatory: cell phones, for example, operate at roughly 1 GHz =10 9 Hz. We exploit the oscillatory nature of the signals by performing a type of time-frequency analysis in the fast time scale. In this manner we convert the original one-dimensional signal into a function of three variables, namely slow time (s), fast time (t or τ) and fast frequency or scale factor σ.

Ambiguity functions
For a single delta-like source at y 0 , it can be shown 3 that over short time intervals, the signal where a is the Doppler scale factor 4 (also known as the Doppler stretch factor), and b is related to the time delay required for the signal to propagate from y 0 to the receiver. The time delay and Doppler scale factor cannot be determined from the data (5), because the transmitted signal is unknown. We can, however, measure the signal at two different receivers, and compare the measurements.
In fact, the usual procedure for detecting the presence of a known expected signal in a received signal contaminated by noise is to use a matched filter [15], which cross-correlates the received signal with the expected signal. In our case, we do not know the transmitted signal, but we do have two received signals, one measured at each receiver. Although both may be contaminated by noise, if both receivers record a signal due to the same transmission, we expect that the noise-free part of the signal recorded at one receiver should be a time-delayed, Doppler-scaled version of the noise-free part of the signal recorded at the other receiver. Consequently we first use the signal at one receiver as a matched filter for the signal at the other receiver. In other words, we compute the (windowed) wideband cross-ambiguity function of the received signals: Here σ is the Doppler scale factor due to the fact that the receivers are in relative motion, and τ is the time delay due to the fact that the receivers are at different distances from the emitter. In (6), the time integration is limited automatically by the support of the window function h.
The time window h is included in (6) in order to separate the different time scales. We refer to the time t as the 'fast' time, corresponding to the time scale on which the waveforms and wave propagation takes place, and s as the 'slow' time. The time delay and Doppler scale factor change with the geometry, which changes only over the slow time scale. The length of the time window should be chosen so that the geometry, and therefore the time delay and Doppler scale, do not change signficantly over the duration of the time window.
Using (5) in (6), we obtain where R m,y (t) = |γ m (t) − y| is the range of the emitter from the receiver platform m. We note that the data (7) depends quadratically on the source terms p z (t) and p y (t) about which we wish to obtain information.
In order to identify the Doppler effect, we perform a Taylor expansion of R about the window center s: where the dot denotes differentiation with respect to the argument of R. We write Using (8) in (7) results in 5 where in the last line we have introduced notation for the Doppler scale factor as a m,y = (1 −Ṙ m,y (s)/c) (and temporarily suppressed the s dependence).
In (9) we make the change of variables t → t ≡ ta l,z − r l,z (s)/c, which has inverse to obtain where we have neglected the Jacobian because it is of the form 1 + O(Ṙ/c). We assume in this work that if multiple emitters are present, each emitter must be able to distinguish its own signal from that of other emitters. Consequently we assume that the emitter waveforms are approximately orthogonal in the sense that where h denotes a windowing function and χ y,s denotes the windowed (wideband) emitter auto-ambiguity function for the waveform p y . We see that (11) is also of the form of a windowed cross-ambiguity function; we can put it in the form of the left side of (12) by noting that the argument of p * y can be written With assumption (12) about orthogonality of the transmitted waveforms, we write (11) as where we have neglected the order Ṙ /c terms in the argument of the smoothly varying window h (see appendix B). In (14), we have used the notation The significance of (14) is that the receiver cross-ambiguity function χ m,l,s , which we measure, is related to the emitter auto-ambiguity function χ y,∆ s,l , about which we wish to obtain information, evaluated at the arguments involving the expressions (15). The arguments of (14) can be interpreted as follows. We note that in the case of a single emitter, the integration over y in (14) disappears, and the transmitter waveform auto-ambiguity function on the right side has a maximum when σ and τ are related to y by The first condition of (16) is known as the frequency difference of arrival (FDOA) condition, and the second is known as the time difference of arrival (TDOA) condition. The set of points y determined by the TDOA condition τ = constant and FDOA condition σ = constant are each algebraic curves (hyperbolas in the case when one receiver is stationary). When the geometry is favorable, these curves intersect in a unique point y, which is the desired source location. Such a favorable situation is indicated in figure 2. The TDOA and FDOA curves depend on the receiver positions and velocities; in unfavorable geometry, for example, the TDOA and FDOA curves could be nearly tangent. The case when both receivers are moving results in FDOA curves that can be much more complicated; see figure 3. This process can be viewed in the following alternative way. In the single-emitter case, the values of σ and τ at the cross-ambiguity maximum give the TDOA and FDOA of the emitter. In other words, the ambiguity function locates the emitter in the TDOA-FDOA coordinate system, and under favorable conditions it is possible to transform the TDOA-FDOA coordinates into spatial coordinates ( y).
Difficulties arise when there are multiple emitters present or when the cross-ambiguity function is ridge-like, so that the precise location of the maximum may be difficult to determine. For these circumstances, we perform an additional matched filtering step: we use the receiver auto-ambiguity function as a matched filter ('second filter') for the cross-ambiguity function.

Receiver auto-ambiguity function
When m = l in (14), the function χ l,l,s (τ , σ) is the receiver auto-ambiguity function: where Thus the receiver auto-ambiguity function (which we measure) is also a shifted, weighted sum of transmitter auto-ambiguity functions (about which we wish to obtain information). When only a single emitter is present, the right side of (17) attains its maximum at σ = 1, τ = T l,l,y,s (1) = 0.

Image formation
The idea is to use a shifted version of the receiver auto-ambiguity function (17) as a matched filter for the cross-ambiguity function (14). In order to determine the correct shift, we need to find the correct arguments to use in the auto-ambiguity function. In particular, we want to choose τ and σ so that at location y = x = z, the arguments of p * z in Thus we should choose σ so that and choose τ so that Solving (22) for τ , we obtain To form an image, we use the shifted auto-ambiguity function χ * l,l,s (τ (τ , x), σ (σ, x)) as a matched filter for the cross-ambiguity function χ m,l,s (τ , σ) and then integrate over 'slow time' s: where the shifts τ (τ , σ, x) and σ (σ, x) are given by (23) and (21), respectively.

Image analysis
To understand what information is provided by the image (24), we relate (24) to the transmitted waveforms. To do this, we substitute (17) and (11) into (24). After some calculations (see appendix C) , we obtain the approximate result where P y denotes the Fourier transform of p y . The t integral of (25) is the (windowed) power of the emitter waveform p y . The ω integral of (25) has its maximum when the arguments of P y match and when the exponent is zero. These conditions are 1 = a m,x a l,y a m,y a l, which are both satisfied at x = y. The integration over s provides a coherent sum similar to that of synthetic-aperture imaging. Thus we expect the image to show a peak at the correct spatial location, and the intensity of the image is determined by the signal energy.
We can understand (25) by considering each slow time s separately. In particular, we can think of the (windowed) source energy density (the t integral) as the quantity we want to image. The ω integral is then the point-spread function for that s. These individual images are then combined coherently as s ranges over the time interval for which we have measurements.
We see that each slow-time point-spread function depends not only on the receiver trajectories, but also on the waveform transmitted. We leave a systematic exploration of this pointspread function for the future. In the following section we consider a very simple special case.

Approximate analysis of (25) in a special case
In this section we carry out an approximate analysis of (25) in the case when the emitted waveform is the duration-T continuous-wave (CW) pulse where the rect function is 1 when its argument is less than 1/2 and 0 otherwise. The Fourier transform of p is where sinc x = (sin x)/x. The main lobe of (28) has width 2|ω − ω 0 | = 4π/T , which implies that for large duration T, the CW pulse (27) has a spectrum that is concentrated around ω = ω 0 .
Replacing ω by ω 0 in the denominator of the ω integration of (25) results in the ω integral being proportional to the frequency-domain expression for the (unwindowed) wideband ambiguity function [14]. For a narrowband signal, the wideband ambiguity function reduces to the narrowband ambiguity function (see appendix A): The narrowband ambiguity function for (27) is [8] When T is large, this ambiguity function is a long, thin ridge in the τ direction; in other words, the long-duration CW pulse corresponds to good resolution in the Doppler direction but poor resolution in the delay direction.
If the geometry is such that we can consider the delay (TDOA) resolution and Doppler (FDOA) resolution separately (for example, in the geometry of figure 2), then for the peak-tonull resolution, we have roughly |τ (x, y)| T and a m,x a l,y a m,y a l, These conditions give the 'thickness' of the TDOA and FDOA curves, respectively. From (26), for the geometry shown in figure 2, where one receiver is stationary and the second receiver's trajectory is perpendicular to the direction to the first receiver, the FDOA condition of (32) translates into where λ 0 = 2πc/ω 0 is the wavelength corresponding to angular frequency ω 0 . Since R m,x ≈ |γ m (s) − x|/R, (33) translates into resolution in the direction of the flight velocity vector: For the example considered in appendix B, in which R = 10 4 m, |γ| = 100 m s −1 and λ 0 = .3 m, we find FDOA resolution of ∆R ≈ 5T m. Similarly, the TDOA conditions of (26) and (32) translate into which translates roughly into |x − y| cT . When T is large, this resolution is poor. Consequently we expect the TDOA localization to be accomplished by viewing the emitter from different locations along the synthetic aperture. A flight path encircling the emitter would then result in resolution (34).

Conclusions
We have developed an imaging algorithm for showing the locations of multiple transmitters that could be transmitting broadband, difficult-to-detect waveforms. The algorithm proceeds by first computing the receiver cross-ambiguity function and the receiver auto-ambiguity function. The receiver auto-ambiguity function, appropriately shifted, is then used as a matched filter for the receiver cross-ambiguity function in an image formation algorithm similar to standard synthetic-aperture imaging. We leave for the future the following interesting questions.
• How can the resolution of images be predicted in the case of more general geometry and more general emitted signals? Figure 3 suggests that the general case will involve not only signal processing but also algebraic geometry. • How robust is this method in the presence of additive noise and uncertainties in receiver position and velocity? • Can an approach to emitter localization based on sparsity be developed?
• Is it possible to develop a similar theory that would apply to interfering emitters?

Acknowledgment
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-14-1-0185 6 . We also thank the Office of Naval Research for support of NRL in-house research which contributed to this effort. We would also like to thank the anonymous reviewers for their comments, which helped us to improve the exposition.

Appendix B. Example of window length
For a situation in which the receiver is 10 km from the emitter and flies in a straight line, perpend icular to the direction to the emitter, at a speed of 100 m s −1 , we carry out the following calculation to find an appropriate window length. To make this less than, say, one-tenth of a full cycle, or less than 2π/10 radians, we should have (t − s) 2 6 · 10 −2 or (t − s) .2 s. A window of length 2 · .2 = .4 s will encompass approximately 5 · 10 8 cycles of the carrier wave. A window of this length corresponds to a largest Fourier frequency of 1/.4 = 5/2. In this example, Ṙ /c ≈ 3 · 10 −7 . The term in the window argument neglected in (14) therefore corresponds to the fraction of a full cycle of the highest-frequency component of the window. 6 Consequently the US Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the Air Force Research Laboratory or the US Government.