Image analysis and resolution for detection-based synthetic-aperture passive source localization

This paper addresses the problem of determining the locations of stationary sources, from which are radiating unknown electromagnetic or acoustic waves, using data measured by passive receivers moving along known flight paths. Examples of sources of interest are radars and wireless communication systems ² . The ambient medium is assumed to be homogeneous with speed of wave propagation c. Both sources and receivers are assumed to be isotropic.

This paper does not consider imaging in clutter [1], nor does it address imaging through random media [2]. Nor does it consider 'passive radar' [3–6] or 'passive coherent localization' [7] or 'hitchhiker' imaging (e.g. [8, 9]), which are systems that form images of targets using signals from sources of opportunity that reflect from the targets to the sensors. This latter problem ³ might seem similar to the source localization problem, but in fact the two problems differ significantly in their statistical formulation [4, 10].

The passive source localization problem considered here can be broken down into two components: first study (1) the statistical problem of detecting the presence of a common source signal appearing in measurements on separate receivers, and then, given the presence of a common signal, (2) determine the location of the source.

The statistical component (1) has been addressed, for example in [4, 10–14]. Some of this work includes the observation that the detection statistic involves quantities computed from a matrix of receiver cross-ambiguity functions. Two types of signal models for Component (1) are relevant for different types of sources. When the sources are random, the incoherent field approximation [9, 15] can be used, or the covariance could be modeled as in [16]. In this paper, however, we develop the theory relevant for sources transmitting deterministic (but unknown) coded signals such as chirps. Such signals are typically narrowband.

The localization component (2) has been studied from a number of points of view, including triangulation and 'two-step' methods such as [17–21], in which pairwise time-difference-of-arrivals (TDOAs) and frequency-difference-of-arrivals (FDOAs) are found first, and then the source location is obtained by solving the system of nonlinear TDOA and FDOA equations. Approaches to Component (2) alone have also been proposed that involve collecting data over a synthetic aperture [22–24]; these papers did not address Component (1). In addition, [22, 23] both made extra assumptions (orthogonality or separability) on the transmitted waveforms in order to obtain a linear imaging method.

Both components (1) and (2) are addressed simultaneously in the 'direct position determination' approach as in [25–27], in which all the data is used in a maximum-likelihood (ML) estimate for the position. The work [26], in particular, sets up the problem for a synthetic-aperture array of receivers, and observes that the ML estimation involves a matrix of cross-ambiguity functions. These cross-ambiguity functions are quadratic in the data, and consequently standard linear approaches for obtaining image resolution cannot be applied. In particular, [26] does not analyze resolution or other properties of the resulting image.

Both components (1) and (2) have also been studied together in [6, 28], which develop detectors for source localization and also form images by plotting a detection statistic, typically as a function of parameters such as delay and Doppler shift. The idea here [6, 10, 13, 14, 25–28] is that when the hypothesized parameters are correct and consequently the shifted signals are properly aligned, the detection statistic will be larger than in the misaligned case, and the resulting image will show a peak at the correct location. These detection statistics are typically nonlinear functions of the data; consequently a theoretical analysis of the resulting image is difficult.

The problem of passive source localization also occurs in radio astronomy, where it is usually addressed with large sensor arrays. The paper [24] proposes a synthetic-aperture approach using only temporal information; it is thus closely related to a multi-sensor, single-source version of [23],

The only papers to estimate image resolution are [23, 24]; but these papers did not address the issue of noise. The approach of [23] required various extra assumptions to convert the quadratic imaging method to a linear method, to which standard resolution analysis applies. In the present paper, on the other hand, we develop a method for obtaining resolution for the nonlinear synthetic-aperture approach.

The new contributions of the present paper are the following.

• (i)  
  We develop a synthetic-aperture passive source localization (SAPSL) algorithm that, unlike the previous synthetic-aperture approaches [22, 23], is based on a solid statistical foundation. This work answers the questions implicit in [22, 23] regarding whether it is better to start with a cross-correlation or a cross-ambiguity function, and how measurements with widely different noise levels should be handled.
• (ii)  
  We extend the prior statistical approaches [6, 10, 13, 14, 25–28] to a synthetic-aperture formulation. Our approach, which is based on the generalized likelihood ratio (GLR) approach of [10, 13, 14], is different from [6, 25–28] in which the unknown attenuation factors ('gains') are estimated first, before the unknown waveform is estimated. In this paper these gains are assumed to be known; in other words we separate the imaging problem from the 'autofocus' problem. Thus we study the imaging problem in which the sensor positions are known and the system is well-calibrated and synchronized. We treat the autofocus problem in a separate paper.
• (iii)  
  This paper is the first to show how the full SAPSL point-spread function, which is a sum of shifted ambiguity functions, can be expressed as an integral over wavenumber space or ' k -space'. This is nontrivial because the imaging algorithm depends quadratically on the data. This analysis, moreover, can be applied also to other problems (e.g. [29]) in which the image is expressed as a sum of ambiguity functions.
• (iv)  
  This paper is the first to provide a comprehensive resolution analysis for the full nonlinear SAPSL method. Here resolution is defined as the width of the point-spread function. This resolution analysis is an extension of the analysis in [23], which studied a linearized version of the method for the case of broadband waveforms.We provide some examples of the resolution calculation for a simple representative geometry.
• (v)  
  This paper provides an optimal imaging filter. This was done in [23] only under the linearizing approximations; it is not done in any of the other previous work on passive source localization.
• (vi)  
  This paper addresses the issue of redundant information (overlapping regions in k -space) in the data. This is also new.
• (vii)  
  This paper also analyzes the case of multiple sources transmitting different waveforms. This is new and replaces the ad hoc assumptions used in [22, 23].

The paper is organized as follows.

• Section 2 sets out the mathematical model and assumptions.
• Section 3 derives the image formation method, following the development of the detector developed in [13, 14].
• Section 4 analyzes the image formed by the method in section 3; this section contains the new results. Each of the four subsections derives a new result.The first subsection is concerned with resolution: we show how to express the point-spread function in terms of the k -space data collection manifold, and from the data collection manifold, we derive resolution. We provide explicit resolution formulas for an example of a simple representative geometry; the resolution formulas of [23, 24] are special cases.The second subsection shows how to determine an optimal imaging filter; the third subsection discusses the case in which the data collection involves redundant information (overlapping regions in Fourier space); the fourth subsection analyzes the case of multiple sources.
• Section 5 summarizes the conclusions and directions for further work.
• The paper concludes with three appendices containing details of certain calculations and a detailed discussion of figure 4.

We consider sources located at unknown positions $\boldsymbol{x}_j , \, j = 1, 2, \ldots$, each transmitting an (unknown) waveform of the form $\psi_{\boldsymbol{x}_j}(t) = \hbox{Re} \left[ e^{i \omega_j t} a_{\boldsymbol{x}_j} (t) \right]$. Here ω_j denotes the angular carrier frequency for the source j and $a_{\boldsymbol{x}_j}$ is the complex baseband signal transmitted from location x _j . Consequently the overall time-varying source distribution is $\rho(t, \boldsymbol{x}) = \sum_j \psi_{\boldsymbol{x}} (t) \delta(\boldsymbol{x} - \boldsymbol{x}_j)$.

We assume that all the transmitted signals are narrowband, which means that they are each of the form $\hbox{Re} \left[ e^{i \omega_0 t} a (t) \right]$, where the bandwidth of the complex baseband signal a(t) is small relative to the carrier frequency ω₀. Typical examples are the signals transmitted by weather radars and airport surveillance radars. These systems typically transmit at S-band ($\omega_0 \sim$ 2 - 3 GHz), with pulse duration on the order of 1–100 microseconds, with bandwidth of a few MHz, and pulse repetition frequency of about 1 kHz. These signals, as well as most wireless communications signals, are all narrowband.

The transmitted waves are assumed to propagate through a homogeneous medium with speed c. Here we model propagation of electromagnetic waves with the scalar wave equation

$$\begin{equation} \left( \nabla^2 - \frac{\partial^2 }{c^2 \partial t^2} \right) u\left(t, \boldsymbol{x}\right) = -\rho\left(t, \boldsymbol{x}\right) . \end{equation} \tag{ 1 }$$

The waves are received by L sensors moving along known trajectories $\boldsymbol{\gamma}_\ell, \, \ell = 1, 2, \ldots, L$, so that the signal at sensor $\ell$ is then $d_\ell^{RF}(t) = u(t, \boldsymbol{\gamma}_\ell(t))$, which is given by

$$\begin{align} d_\ell^{RF}\left(t \right) & = \int \frac{\delta\left(t - t^{^{\prime}} - |\boldsymbol{x}^{^{\prime}} - \boldsymbol{\gamma}_\ell\left(t \right)| /c\right) }{4 \pi |\boldsymbol{x}^{^{\prime}} - \boldsymbol{\gamma}_\ell\left(t \right)| } \rho\left(t^{^{\prime}}, \boldsymbol{x}^{^{\prime}}\right) dt^{^{\prime}} d\boldsymbol{x}^{^{\prime}} + \hbox{noise}\nonumber\\ & = \sum_j \frac{\psi_{\boldsymbol{x}_j} \left( t - \tau_\ell\left(\boldsymbol{x}_j, t \right) \right) }{4 \pi R_{\ell} \left(\boldsymbol{x}_j, t \right) } + \hbox{noise} = \sum_j G_{\ell}\left(\boldsymbol{x}_j , t \right) \psi_{\boldsymbol{x}_j}\left(t - \tau_{\ell}\left(\boldsymbol{x}_j, t \right)\right)+ \hbox{noise}, \nonumber \end{align} \tag{ 2 }$$

where $R_{\ell}(\boldsymbol{x}, t) = | \boldsymbol{x} - \boldsymbol{\gamma}_\ell(t )|$, $\tau_{\ell}(\boldsymbol{x}, t) = R_{\ell} (\boldsymbol{x}, t) / c$, and, for the isotropic sources considered in this paper, $G_{\ell}(\boldsymbol{x}, t) = (4\pi R_{\ell}(\boldsymbol{x}, t))^{-1}$. We assume that all sensor clocks are synchronized.

The goal is to determine the locations x _j from the data $\{d_\ell (t), \, \ell = 1, 2, \ldots, L \}$. We assume we know the receiver trajectories $\boldsymbol{\gamma}_\ell$ and associated velocities $\dot{\boldsymbol{\gamma}}_\ell$, but the waveforms a_j are unknown, as are the ranges $R_{\ell} (\boldsymbol{x}_j, t)$ and time delays $\tau_{\ell}(\boldsymbol{x}_j, t)$. A model similar to (2) applies to the case of moving sources, but typically then the source velocities are also unknown.

We note that this is a nonstandard type of inverse problem, in that the forward map (2) from $\rho(t, \boldsymbol{x})$ to $\{d_\ell^{RF}(t) \}$ requires knowledge not only of the source positions, but also of the waveforms transmitted from those positions. The linear inverse problem of determining $\rho(t, \boldsymbol{x})$ from $\{d_\ell^{RF} (t) \}$, however, would be severely underdetermined.

We begin by assuming that some method (e.g. [30–32]) has been used to estimate one of the carrier frequencies, say ω₀. If the estimate is slightly wrong, we merely redefine $a_{\boldsymbol{x}_j}$ to include the extra factor of $e^{i (\omega_j - \omega_0)t}$.

We then demodulate [33–35] the signal (2) by the carrier frequency ω₀, obtaining a complex signal model in which all sources transmitting at frequency ω₀ are converted to baseband:

$$\begin{align} d_{\ell}^{BB}\left(t \right) & = \sum_{\left\{j: \omega_j = \omega_0\right\} } G_\ell \left(\boldsymbol{x}_j, t \right) e^{- i \omega_0 \tau_\ell \left(\boldsymbol{x}_j, t \right) } a_{\boldsymbol{x}_j} \left(t - \tau_\ell \left(\boldsymbol{x}_j, t \right)\right) + \hbox{noise}. \end{align} \tag{ 3 }$$

The noise term in (3) is different from that in (2); in particular, not only have we removed the carrier frequency, but we have also lumped the signals from sources transmitting at frequencies other than ω₀ into 'noise'. This is not necessarily the best possible way to treat all possible interfering signals, but this assumption makes the analysis tractable. Exploration of other approaches is left for future work.

Depending on the speed with which the receiving antennas move, the sensor geometry typically changes on the order of seconds, which is a factor of 10⁶ to 10¹⁰ greater than the period of the oscillatory carrier wave, a factor of 10⁴ greater than a typical radar pulse duration, and a factor of 10³ greater than typical pulse repetition intervals. In order to isolate individual pulses or groups of pulses and 'freeze' the geometry, we extract a snippet of the data around the 'slow time' $t = s_m$, say by multiplying by a smooth window (see figure 1). For example, for a radar with a pulse repetition frequency of 1 kHz, a suitable time window might have duration 10⁻³ seconds, so that each window would capture the response from a single pulse ⁴ . A similar windowing operation can be applied to non-pulsed signals, such as communications signals. The windows, taken together, can form a partition of unity, so that no data is lost.

Zoom In Zoom Out Reset image size

In the window corresponding to slow time s_m , we make the Taylor expansion

$$\begin{equation} R_{\ell}\left(\boldsymbol{x}_j, t\right) = R_{\ell,m} \left(\boldsymbol{x}_j\right)+ \dot{R}_{\ell,m} \left(\boldsymbol{x}_j\right)\left(t-s_m\right) + \ldots, \end{equation} \tag{ 4 }$$

where $R_{\ell,m} (\boldsymbol{x}_j) = R_{\ell}(\boldsymbol{x}_j, s_m)$, etc. This expansion induces the corresponding expansion in τ:

$$\begin{equation} \tau_{\ell}\left(\boldsymbol{x}_j, t\right) = \tau_{\ell,m} \left(\boldsymbol{x}_j\right)+ \dot{\tau}_{\ell,m} \left(\boldsymbol{x}_j\right)\left(t-s_m\right) + \cdots. \end{equation} \tag{ 5 }$$

We note that $\dot{\tau}$ is of order $v/c$, where v is a sensor velocity, and is consequently small; it can be neglected in the argument of the slowly varying function a_j but must be retained in the exponent where it is multiplied by a large frequency. We introduce the 'slow time' reference points s_m , and then make the 'fast time' change of variables $t \mapsto \tilde{t} = t - s_m$ in (3):

$$\begin{equation} d_{\ell,m}^{BB} \left(\tilde{t}\right) = \sum_j \underbrace{G_{\ell,m}\left(\boldsymbol{x}_j\right) e^{-i \omega_0 \tau_{\ell,m} \left(\boldsymbol{x}_j\right)} }_{\displaystyle g_{\ell,m} \left(\boldsymbol{x}_j\right)} e^{-i \omega_0 \dot{\tau}_{\ell, m} \left(\boldsymbol{x}_j\right) \tilde{t} } a_{\boldsymbol{x}_j} \left(s_m + \tilde{t} - \tau_{\ell, m} \left(\boldsymbol{x}_j\right)\right) + \hbox{noise} , \end{equation} \tag{ 6 }$$

where we have used the fact that a_j is slowly varying to neglect the term $\dot{\tau}_{\ell,m} (\boldsymbol{x}_j)\tilde{t}$ in the argument. Here $g_{\ell,m} (\boldsymbol{x}_j) \doteq G_{\ell,m}(\boldsymbol{x}_j) e^{-i \omega_0 \tau_{\ell, m}(\boldsymbol{x}_j) }$, which for an isotropic source is

$$\begin{equation} g_{\ell,m} \left(\boldsymbol{x}_j\right) \doteq G_{\ell,m}\left(\boldsymbol{x}_j\right) e^{-i \omega_0 \tau_{\ell, m}\left(\boldsymbol{x}_j\right) } = \frac{e^{-i k_0 R_{\ell}\left(\boldsymbol{x}_j, s_m\right) } } {4\pi R_{\ell}\left( \boldsymbol{x}_j, s_m \right)}, \end{equation} \tag{ 7 }$$

where $k_0 = \omega_0 /c$. We note that $g_{\ell, m}(\boldsymbol{x}_j)$ is nonzero and is independent of the fast time t.

In (6), we write $a_{\boldsymbol{x}_j,m}(t) = a_{\boldsymbol{x}_j} (s_m + t)$, write $\nu_{\ell, m} (\boldsymbol{x}_j) = \omega_0 \dot{\tau}_{\ell, m} (\boldsymbol{x}_j)$ for the angular Doppler-shift frequency, and drop the tildes, to obtain the baseband data model

$$\begin{equation} d_{\ell,m}^{BB} \left(t \right) = \sum_j g_{\ell,m} \left(\boldsymbol{x}_j\right) e^{-i \nu_{\ell, m} \left(\boldsymbol{x}_j\right) t } a_{\boldsymbol{x}_j,m} \left( t - \tau_{\ell, m} \left(\boldsymbol{x}_j\right)\right) + \hbox{noise}. \end{equation} \tag{ 8 }$$

The sum over j is only over those sources that are transmitting at carrier frequency ω₀. For sources such as radars, which typically transmit repeated pulses, and windows chosen to correspond to those pulses, the mth baseband waveform $a_{\boldsymbol{x}, m}$ corresponds simply to the mth pulse transmitted from location x .

Our approach is similar to that of [26], except that instead of the maximum likelihood, we use the GLR detection statistic derived in [10, 14]. An image is formed by the following steps:

• (i)  
  For each hypothesized source location y , we time-shift and Doppler shift the data measured at each receiver to remove the propagation effects (time delay and Doppler shift) assuming the source is at y . This process can be thought of as aligning the data.We note that the complex attenuation/gain $g_{\ell, m}$, which has the assumed form (7), is also known since y is assumed known.
• (ii)  
  Once the data is aligned, we apply a statistical test [equations (23), (24), and (28)], the derivation of which follows [10, 14], to determine the degree to which the aligned signals are the same. The statistic we use indicates the likelihood that a common signal is present at the various receivers. We plot the likelihood ratio as the image value at pixel y .
• (iii)  
  We apply the same process to each hypothesized source location (image pixel) y in turn.

We describe these steps in more detail below.

3.1. Alignment of the data for a hypothesized location

For a hypothesized source position y , the time delay and Doppler shifts appearing in (8) would be

$$\begin{align} \tau_{\ell, m} \left(\boldsymbol{y}\right) & = |\boldsymbol{y} - \boldsymbol{\gamma}_\ell\left(s_m\right) |/c \nonumber\\ \nu_{\ell, m} \left(\boldsymbol{y}\right) & = \omega_0 \dot{\tau}_{\ell, m} \left(\boldsymbol{y}\right)/c = - \omega_0 \widehat{\boldsymbol{R}}_{\ell, m} \left(\boldsymbol{y}\right) \cdot \dot{\boldsymbol{\gamma}}_\ell\left(s_m\right)/c . \end{align} \tag{ 9 }$$

Consequently, to test $\{ d_{\ell, m}^{BB} \}$ for the presence of a common signal emanating from location y , we apply the appropriate time advance and Doppler (un-)shift to each $d_{\ell, m}^{BB}$:

$$\begin{align} \eta_{\ell, m} \left(\boldsymbol{y}, t \right) \doteq e^{i \nu_{\ell, m} \left(\boldsymbol{y}\right)\left( t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) }d_{\ell, m}^{BB} \left(t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) . \end{align} \tag{ 10 }$$

In terms of the data model (8), η can be written

$$\begin{align} & \eta_{\ell, m} \left(\boldsymbol{y}, t \right)\nonumber\\ & \quad = e^{i \nu_{\ell, m} \left(\boldsymbol{y}\right)\left( t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) } d_{\ell, m}^{BB} \left(t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) \nonumber\\ & \quad = e^{i \nu_{\ell, m} \left(\boldsymbol{y}\right) \left( t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) } \sum_j g_{\ell,m} \left(\boldsymbol{x}_j\right) e^{-i \nu_{\ell, m} \left(\boldsymbol{x}_j\right) \left(t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) } \ a_{\boldsymbol{x}_j,m} \left( \left(t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) - \tau_{\ell, m} \left(\boldsymbol{x}_j\right) \right) + \hbox{noise} \nonumber\\ & \quad = \sum_j g_{\ell,m} \left(\boldsymbol{x}_j\right) e^{i \Delta \nu_{\ell, m}\left(\boldsymbol{y}, \boldsymbol{x}_j\right) \left( t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) } \ a_{\boldsymbol{x}_j,m} \left( t + \Delta \tau_{\ell, m}\left(\boldsymbol{y}, \boldsymbol{x}_j\right) \right) + \hbox{noise} \end{align} \tag{ 11 }$$

where $\Delta \tau_{\ell, m}(\boldsymbol{y}, \boldsymbol{x}_j)$ and $\Delta \nu_{\ell, m}(\boldsymbol{y}, \boldsymbol{x}_j)$ are the differences

$$\begin{align} \Delta \tau_{\ell, m}\left(\boldsymbol{y}, \boldsymbol{x}_j\right) & = \tau_{\ell, m} \left(\boldsymbol{y}\right)- \tau_{\ell, m} \left(\boldsymbol{x}_j\right) \qquad \hbox{delay difference} \nonumber\\ \Delta \nu_{\ell, m}\left(\boldsymbol{y}, \boldsymbol{x}_j\right) & = \nu_{\ell, m} \left(\boldsymbol{y}\right) - \nu_{\ell, m} \left(\boldsymbol{x}_j\right) \qquad \hbox{frequency difference} . \end{align} \tag{ 12 }$$

Figure 2 shows data, as measured at the sensors after it has been demodulated, time-advanced, and Doppler-shifted, due to a repeating sequence of linear chirps transmitted from the origin. The transmitted chirps have bandwidth of 1.6 kHz, center frequency of 3 GHz, duration 1 second, and a pulse repetition interval of 2 s. The processing windows are 1 second long. These parameters are chosen for visibility in the plots and not for realism. The geometry is shown in figure 3: sensor 1 starts at (−10 km,−10 km, 5 km) and moves at 360 m s⁻¹ in the direction (1,0,0); sensor 2 starts at (0, −10 km, 5 km) and moves at the same speed in the same direction. The top row of figure 2 shows graphs of initial segments of the signals received at each of 6 slow times. The straight lines at slow times 2, 4, and 6 correspond to slow times at which no signal is received. The bottom row shows spectrograms, i.e. windowed Fourier transforms, using sliding windows much shorter than the 1-second processing window. The spectrograms are plotted on a dB scale. The dark vertical stripes in the spectrograms correspond to windows in which no signal is received; these correspond to the horizontal lines at slow times 2, 4, and 6 in the top row.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If there is indeed a source at location y , then the frequency difference and delay difference are zero, and η becomes

$$\begin{align} \eta_{\ell, m} \left(\boldsymbol{y}, t \right) & = g_{\ell, m} \left(\boldsymbol{y}\right)a_{\boldsymbol{y}, m} \left( t \right) \nonumber\\ & \quad + \sum_{\left\{j: \boldsymbol{y} \neq \boldsymbol{x}_j \right\} } g_{\ell,m} \left(\boldsymbol{x}_j\right) \exp \left( i \left(\nu_{\ell, m} \left(\boldsymbol{y}\right) - \nu_{\ell, m} \left(\boldsymbol{x}_j\right)\right) \left( t + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) \right) \nonumber\\ & \quad \times \ a_{\boldsymbol{x}_j ,m}\left( t + \tau_{\ell, m} \left(\boldsymbol{y}\right)- \tau_{\ell, m} \left(\boldsymbol{x}_j\right) \right) + \hbox {noise} , \end{align} \tag{ 13 }$$

where the subscript y corresponds to the appropriate index j for which $\boldsymbol{y} = \boldsymbol{x}_j$. The sum over j for which $\boldsymbol{y} \neq \boldsymbol{x}_j$ is lumped in with the rest of the noise. Thus we assume that when a source is present at y ,

$$\begin{equation} \eta_{\ell, m} \left(\boldsymbol{y}, t \right) = g_{\ell, m} \left(\boldsymbol{y}\right)a_{\boldsymbol{y}, m} \left( t \right) + \hbox {noise} . \end{equation} \tag{ 14 }$$

Here we are using a strategy similar to that used in active synthetic-aperture radar, where one develops an imaging algorithm based on a model for scattering from a single point-like target, even though scenes of interest typically consist of a large number of scatterers [35]. In section 4.4, we address in detail the multiple-source case.

If the flight trajectories are known precisely, then $g_{\ell, m} (\boldsymbol{y})$ of (14) is known from (7) with our choice of y and ω₀. The left side of (14) is obtained from measurements, but $a_{\boldsymbol{y}}$ on the right side is unknown. Thus the problem of determining whether a source is present at location y can be formulated as the problem of determining whether the signals $\eta_{\ell, m}(\boldsymbol{y}, t )$, obtained from measurements at different sensors $\ell$, are proportional to a common (aligned) signal $a_{\boldsymbol{y}, m} (t )$, plus noise.

3.2. Estimation of the waveform

Under the assumption that the hypothesized source is at location y , we first determine an estimate of the transmitted waveform(s) $a_{\boldsymbol{y}, m}$. To do this we make use of the data on all sensors for a particular slow time m, and follow the procedure below. For each m, we assemble these data into the L-dimensional vector-valued functions $\boldsymbol{\eta}_m(t)$, whose $\ell$th element is $\eta_{\ell, m}(t)$, and similarly for g _m and d _m . We use the superscript $^\dagger$ for adjoint. Here we temporarily drop the y subscripts and arguments.

We assume that the noise on each sensor $\ell$ is bandlimited complex Gaussian white noise with mean zero and known variance $\sigma^2_\ell$. With the L-dimensional vector notation η and g , and with the notation $\mathcal R = \mathrm{diag} \{\sigma_1^2, \sigma_2^2, \ldots, \sigma_L^2\}$ for the covariance, we write the Gaussian likelihood for η as proportional to

$$\begin{equation} \prod_m \exp \left( - \int \left( \boldsymbol{\eta}_m\left(t\right) - \boldsymbol{g}_m a_m\left(t\right)\right)^\dagger {\mathcal R}^{-1} \left(\boldsymbol{\eta}_m\left(t\right) - \boldsymbol{g}_m a_m\left(t\right) \right) dt\right) . \end{equation} \tag{ 15 }$$

Here for each m, the t integration extends over the support of the mth time window. The integrals in (15) can be obtained by a large-bandwidth limiting process as described in section 6.4 of [36]. More details are given in appendix A. In the following, we temporarily drop the argument t.

To maximize (15), we minimize, for each m, the weighted mean-square error

$$\begin{equation} E\left(a_m\right) = \int \left( \boldsymbol{\eta}_m - \boldsymbol{g} _m a_m\right)^\dagger {\mathcal R}^{-1} \left(\boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right) dt \end{equation} \tag{ 16 }$$

with respect to a_m , which we do by a variational argument. Thus for all perturbations $\tilde{a}_m$, we have

$$\begin{equation} 0 = \left. \frac{\partial}{\partial \epsilon} \right|_{\epsilon = 0} E\left(a_m + \epsilon \tilde{a}_m\right) = -2 \hbox{Re } \int \boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \left(\boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right) \, \tilde{a}_m \, dt \end{equation} \tag{ 17 }$$

which must be valid for all complex $\tilde{a}_m$. Thus we obtain the pseudoinverse expression

$$\begin{equation} a_{\boldsymbol{y}, m}\left(t\right) = \left(\boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{g}_m\right)^{-1} \left( \boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m \right) = \frac{4\pi} {\sum_\ell \frac{1}{R_{\ell}^2 \left(\boldsymbol{y}, s_m\right) \sigma_{\ell}^2} } \left( \sum_\ell \frac{ e^{i k_0 R_\ell\left(\boldsymbol{y}, s_m\right) } \eta_{\ell, m}\left(t\right) }{ R_\ell\left(\boldsymbol{y}, s_m\right) \ \sigma_\ell^2 } \right) . \end{equation} \tag{ 18 }$$

This is the ML estimate for the waveform transmitted from hypothesized location y during the mth slow-time interval.

3.3. Derivation of the GLR detection statistic

This derivation follows [10, 14], using slightly different notation. For the hypothesized point y , we consider two hypotheses, namely H₀, the hypothesis that the data η_m consists only of noise, and H₁, the hypothesis that the data consists of noise plus a signal $\boldsymbol{g}_m a_m$ common to the sensors. The likelihoods under H₀ and H₁ are proportional to

$$\begin{align} H_0: & \qquad \exp \left( - \sum_m\int \boldsymbol{\eta}_m ^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m dt \right) \nonumber\\ H_1:& \qquad \exp \left( - \sum_m \int \left( \boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right)^\dagger {\mathcal R}^{-1} \left(\boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right) dt \right) . \end{align} \tag{ 19 }$$

The log of the likelihood ratio is $\mathcal L (\boldsymbol{y}) = \sum_m \mathcal L_m (\boldsymbol{y})$, where

$$\begin{equation} \mathcal L_m\left(\boldsymbol{y}\right) = \int \left( \boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right)^\dagger {\mathcal R}^{-1} \left(\boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right) dt - \int \boldsymbol{\eta}_m ^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m dt . \end{equation} \tag{ 20 }$$

The next step is to insert (18) and rearrange; the details are given in appendix B. The result is

$$\begin{align} \mathcal L_m\left(\boldsymbol{y}\right) = \frac{ \int \boldsymbol{\beta}_m^\dagger \boldsymbol{M}_m\left(t\right) \boldsymbol{\beta}_m dt } {| \boldsymbol{\beta}_m|^2} , \end{align} \tag{ 21 }$$

where $\boldsymbol{\beta}_m = {\mathcal R}^{-1/2}\boldsymbol{g}_m = (g_{1,m}/\sigma_1, \ldots g_{L,m}/\sigma_L)$, and $| \boldsymbol{\beta}_m|^2 = \boldsymbol{\beta}_m^\dagger \boldsymbol{\beta}_m$, and where $\boldsymbol{M}_m(t)$ is the L × L self-adjoint matrix

$$\begin{align} \boldsymbol{M} _m\left(t\right) & = {\mathcal R}^{-1/2} \boldsymbol{\eta}_m \left(t\right) \boldsymbol{\eta}_m^\dagger \left(t\right) {\mathcal R}^{-1/2}. \end{align} \tag{ 22 }$$

Equation (21) can be written as the sum

$$\begin{equation} \mathcal L_m \left(\boldsymbol{y}\right) = \sum_{\ell, \ell^{^{\prime}}} \mathcal L_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) , \end{equation} \tag{ 23 }$$

where

$$\begin{align} \mathcal L_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) &\doteq \frac{\overline{\beta }_{\ell,m} \left(\boldsymbol{y}\right) \beta_{\ell^{^{\prime}},m} \left(\boldsymbol{y}\right) }{|\boldsymbol{\beta}_m\left(\boldsymbol{y}\right)|^2} \frac{\langle \eta_{\ell, m} , \eta_{\ell^{^{\prime}},m} \rangle } {\sigma_\ell \sigma_{\ell^{^{\prime}}}} = \frac{\overline{g_{\ell, m}\left(\boldsymbol{y}\right)} g_{\ell^{^{\prime}},m} \left(\boldsymbol{y}\right) } { \sigma_\ell^2 \sigma_{\ell^{^{\prime}}}^2 \sum_{\ell^{^{\prime}\!\!^{\prime}}} |g_{\ell^{^{\prime}\!\!^{\prime}}, m} / \sigma_{\ell^{^{\prime}\!\!^{\prime}}}|^2 } \langle \eta_{\ell, m} , \eta_{\ell^{^{\prime}},m} \rangle , \end{align} \tag{ 24 }$$

where $\langle \,f, g \rangle = \int f(t) \overline{g(t)} dt$ denotes the usual L² inner product and the bar denotes complex conjugate. We note that (24) contains factors that automatically weight contributions of the data from sensors $\ell$ and $\ell^{^{\prime}}$ according to the (known) signal-to-noise ratios at those sensors. We note also that $\mathcal L_{\ell, \ell^{^{\prime}}, m} = \overline{\mathcal L_{\ell^{^{\prime}}, \ell, m}}$.

In appendix C we show that the inner product of (24) can be written

$$\begin{equation} \langle \eta_{\ell, m} \left(\boldsymbol{y}, \cdot\right) , \eta_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \cdot\right) \rangle = \displaystyle \mathcal A_{\ell, \ell^{^{\prime}}} \left( T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) , \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) \right) \ e^{-i \nu_{\ell, m} T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right)} \end{equation} \tag{ 25 }$$

where $\mathcal A$ is the receiver cross-ambiguity function

$$\begin{equation} \mathcal A_{\ell, \ell^{^{\prime}}}\left(\tau, \nu\right) = \int d^{BB}_{\ell, m} \left(t + \tau\right)\ \overline{d^{BB}_{\ell^{^{\prime}}, m}\left(t\right) } e^{i \nu t} \ dt \end{equation} \tag{ 26 }$$

and where T and $\mathcal V$ are the time difference of arrival (TDOA) and frequency difference of arrival (FDOA), respectively:

$$\begin{align} T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) & = \tau_{\ell,m} \left(\boldsymbol{y}\right) - \tau_{\ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) \qquad \ \hbox{(TDOA)} \nonumber\\ \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) & = \nu_{\ell , m} \left(\boldsymbol{y}\right) - \nu_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right) \qquad \hbox{(FDOA)} . \end{align} \tag{ 27 }$$

In (26), the bar denotes complex conjugate.

Equation (25) shows that, as we would anticipate from e.g. [22, 26, 37, 38], the image (24) is constructed from receiver cross-ambiguity functions. Consequently, to compute the contribution to the image at each slow time, the receiver cross-ambiguity function can be computed once, and then evaluated at the computed values of TDOA and FDOA for a grid of various hypothesized points y in the image.

The image $\mathcal L(\boldsymbol{y})$ is formed by summing over all slow times m:

$$\begin{equation} \mathcal L\left(\boldsymbol{y}\right) = \sum_m \mathcal L_m\left(\boldsymbol{y}\right) = \sum_{\ell, \ell^{^{\prime}}, m} \mathcal L_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) . \end{equation} \tag{ 28 }$$

We note that this image formation process depends quadratically on the data through (26).

The likelihoods in (28) are derived under the assumption that the signals are correctly aligned, or in other words, that the source is indeed located at the hypothesized location and thus $\eta_{\ell,m}$ has the form (14). When this assumption is satisfied, the sum (28) is non-negative. In the following section, we investigate the case when the signals are misaligned, or in other words when the hypothesis y for the source position is wrong.

For sources that repeatedly transmit the same waveform, as many radar systems do, there is a possibility that a signal that is time-shifted and Doppler-shifted to correspond to a signal from a distant location y might allow the mth pulse at the receiver to align with the $(m-1)$st pulse from the transmitter. This could result in image artifacts that are related to the range ambiguities that can occur in an active radar system and which can be predicted from the pulse repetition frequency of the radar. In active radar, such ambiguities can be handled by design of the antenna beam pattern so that the ambiguities are not illuminated, or by design of the sequence of transmitted waveforms so that targets at ambiguous ranges have negligible reflected energy. Here we assume that the image region is adequately restricted ⁵ so these ambiguities do not occur.

To determine the fidelity of the image, into (24) we insert the signal model (11). The strategy [35] is to analyze the image by examining its Fourier transform and discovering which Fourier components are determined by the measured data. To do this, we can consider the contributions from each sensor pair, and then add those contributions together coherently.

For simplicity we consider the case of only a single source at location x . We note that because of the quadratic nature of the image formation process, the image of multiple sources is not simply a sum of single-source images. Instead, there are cross terms, which we study in section 4.4. A special case of a similar problem was studied in [23].

The contribution to the image at y from a single sensor pair $\ell, \ell^{^{\prime}}$ and a single slow-time index m is $\mathcal L_{\ell, \ell^{^{\prime}}, m}$, which, from (24) and (11), is

$$\begin{align} \mathcal L_{\ell,\ell^{^{\prime}}, m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \frac{\overline{\beta_{\ell, m}\left(\boldsymbol{y}\right)} \beta_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) } {|\boldsymbol{\beta} \left(\boldsymbol{y}\right)|^2 \sigma_{\ell} \sigma_{\ell^{^{\prime}}} } \int g_{\ell ,m} \left(\boldsymbol{x}\right) e^{i \Delta \nu_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}\right) \left( t + \tau_{\ell , m} \left(\boldsymbol{y}\right)\right) } \ a_{\boldsymbol{x},m} \left( t + \Delta\tau_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}\right) \right) \nonumber\\ & \quad \times \overline{g_{\ell^{^{\prime}} ,m} \left(\boldsymbol{x}\right) e^{i \Delta \nu_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right) \left( t + \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right)\right) } \ a_{\boldsymbol{x}, m} \left( t + \Delta\tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right) \right) } dt, \end{align} \tag{ 29 }$$

where the delay and Doppler differences at sensor $\ell$ from points x and y are

$$\begin{align} \Delta\nu_{\ell, m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \nu_{\ell , m} \left(\boldsymbol{y}\right) - \nu_{\ell , m} \left(\boldsymbol{x}\right) \nonumber\\ \Delta \tau_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}\right) & = \tau_{\ell , m} \left(\boldsymbol{y}\right)- \tau_{\ell , m} \left(\boldsymbol{x}\right) . \end{align} \tag{ 30 }$$

In (29), we make the change of variables $t \mapsto t^{^{\prime}} = t + \Delta\tau_{\ell^{\prime},m} (\boldsymbol{y}, \boldsymbol{x})$. This results in

$$\begin{align} \mathcal L_{\ell,\ell^{^{\prime}},m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \frac{\overline{\beta_{\ell, m}\left(\boldsymbol{y}\right)} \beta_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \beta_{\ell, m}\left(\boldsymbol{x}\right) \overline{\beta_{\ell^{^{\prime}}, m}\left(\boldsymbol{x}\right)} }{|\boldsymbol{\beta} \left(\boldsymbol{y}\right)|^2 } \nonumber\\ & \quad \times \int e^{i \Delta \nu_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}\right) \left( t^{^{\prime}} + \tau_{\ell , m} \left(\boldsymbol{y}\right) - \Delta\tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right)\right) } \ a_{\boldsymbol{x},m} \left( t^{^{\prime}} + \Delta T_{\ell,\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right) \right) \nonumber\\ & \quad \times \overline{ e^{i \Delta \nu_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right) \left(t^{^{\prime}} + \tau_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) - \Delta\tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right)\right) } \ a_{\boldsymbol{x}, m} \left(t^{^{\prime}}\right) } dt^{^{\prime}} \nonumber\\ & = \frac{\overline{\beta_{\ell, m}\left(\boldsymbol{y}\right)} \beta_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \beta_{\ell, m}\left(\boldsymbol{x}\right) \overline{\beta_{\ell^{^{\prime}}, m}\left(\boldsymbol{x}\right) }}{|\boldsymbol{\beta} \left(\boldsymbol{y}\right)|^2 } e^{i \Phi_{\ell,\ell^{^{\prime}}, m} \left(\boldsymbol{y}, \boldsymbol{x}\right) } \!\mathcal A \left(\! \Delta T_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}, \!\!\boldsymbol{x}) , \Delta \mathcal V_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}, \!\boldsymbol{x}) \!\right) , \end{align} \tag{ 31 }$$

where

$$\begin{equation} \mathcal A \left(\tau, \nu\right) = \int a\left(t+\tau\right) \overline{a\left(t\right)} e^{i \nu t} dt \end{equation} \tag{ 32 }$$

is the classical (narrowband) radar ambiguity function, where

$$\begin{align} \Delta T_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{y}, \boldsymbol{x}\right) & = \Delta\tau_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}\right) - \Delta\tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right) = \tau_{\ell , m} \left(\boldsymbol{y}\right)- \tau_{\ell , m} \left(\boldsymbol{x}\right) - \left[ \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right)- \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{x}\right) \right] \nonumber\\ & = T_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) - T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{x}\right) \qquad \hbox{difference of TDOAs} \nonumber\\ \Delta \mathcal V_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{y}, \boldsymbol{x}\right) & = \Delta \nu_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}\right) - \Delta \nu_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right) = \nu_{\ell , m} \left(\boldsymbol{y}\right) - \nu_{\ell , m} \left(\boldsymbol{x}\right) - \left[ \nu_{\ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) - \nu_{\ell^{^{\prime}} , m} \left(\boldsymbol{x}\right) \right] \nonumber\\ & = \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) - \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{x}\right) \qquad \hbox{difference of FDOAs} \end{align} \tag{ 33 }$$

with

$$\begin{align} \Phi_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \Delta \nu_{\ell,m} \tau_{\ell , m} \left(\boldsymbol{y}\right) - \Delta \nu_{\ell^{^{\prime}},m} \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right) + \Delta\tau_{\ell^{^{\prime}},m} \left(\Delta \nu_{\ell^{^{\prime}},m} - \Delta \nu_{\ell,m} \right)\nonumber \\ & = \left[ \nu_{\ell , m} \left(\boldsymbol{y}\right) - \nu_{\ell , m} \left(\boldsymbol{x}\right) \right] T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) + \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{x}\right) \underbrace{\left[ \mathcal V_{\ell,\ell^{^{\prime}},m} \left(\boldsymbol{y}\right) - \mathcal V_{\ell,\ell^{^{\prime}},m} \left(\boldsymbol{x}\right) \right] }_{\displaystyle \Delta \mathcal V_{\ell,\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}\right) } , \end{align} \tag{ 34 }$$

where more details are given in appendix D.

We note that when $\ell = \ell^{^{\prime}}$, the quantities $\Delta T_{\ell,\ell^{^{\prime}},m}, \Delta \mathcal V_{\ell, \ell^{^{\prime}}, m}$, and $\Phi_{\ell, \ell^{^{\prime}}, m}$ are all zero, and consequently the diagonal elements $\mathcal L_{\ell, \ell}$ depend only on the signal-to-noise ratio. This means that although the detection statistic is affected by the diagonal terms, the image resolution is determined only by the off-diagonal terms of $\mathcal L$.

4.1. Resolution

To determine resolution, we examine the image (28) for y near the true source location x and determine the width of the peak in the image. To do this, we examine each contributing term (31). The strategy, as is common with many other inverse problems, is to express the image in terms of its spatial Fourier components, i.e. its wavenumber space or ' k -space'.

We write $\mathcal A$ in terms of its Fourier transform W, which is the Rihaczek distribution:

$$\begin{align} \mathcal A\left(\tau, \nu\right) = \int e^{i \left(\tau \omega^{^{\prime}} + \nu t^{^{\prime}} \right)} W\left(\omega^{^{\prime}}, t^{^{\prime}}\right) d\omega^{^{\prime}} dt^{^{\prime}} . \end{align} \tag{ 35 }$$

The Rihaczek distribution W is the instantaneous correlation between $a(t^{^{\prime}})$ and its one-term Fourier expansion $A(\omega^{^{\prime}}) e^{i \omega^{^{\prime}} t^{^{\prime}}}$, which we see below in (36).

4.1.1. Support of W.

For a narrowband waveform a, the support of W is determined by its duration and frequency band. We see this from the following calculation:

$$\begin{align} W\left(\omega^{^{\prime}}, t^{^{\prime}}\right) &\propto \int e^{- i \left(\tau \omega^{^{\prime}} + \nu t^{^{\prime}}\right)} \mathcal A\left(\tau, \nu\right) d\tau d\nu \nonumber\\ & = \int e^{- i \left(\tau \omega^{^{\prime}} + \nu t^{^{\prime}}\right)} \int a\left(t + \tau\right) \overline{a\left(t\right)} e^{i \nu t} dt d\tau d\nu \nonumber\\ &\propto \int e^{-i \tau \omega^{^{\prime}}} a\left(t^{^{\prime}} + \tau\right) d\tau \ \overline{a\left(t^{^{\prime}}\right)} = e^{i t^{^{\prime}} \omega^{^{\prime}}} A\left(\omega^{^{\prime}}\right) \overline{a\left(t^{^{\prime}}\right)} , \end{align} \tag{ 36 }$$

where A denotes the Fourier transform of a and where in the last line we have made the substitution $\tau \mapsto \tau^{^{\prime}} = t^{^{\prime}} + \tau$ and carried out the $\tau^{^{\prime}}$ integration. The expression (36) shows that the temporal support of a is the same as the t' support of W, and the frequency band of a (i.e. the support of A) is the same as the $\omega^{^{\prime}}$ support of W.

4.1.2. Fourier components from one sensor pair.

Using (35) in (31) results in

$$\begin{align} \mathcal L_{\ell,\ell^{^{\prime}},m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \frac{\overline{\beta_{\ell, m}\left(\boldsymbol{y}\right)} \beta_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \beta_{\ell, m}\left(\boldsymbol{x}\right) \overline{\beta_{\ell^{^{\prime}}, m}\left(\boldsymbol{x}\right) }} {|\boldsymbol{\beta} \left(\boldsymbol{y}\right)|^2 } e^{i \Phi_{\ell,\ell^{^{\prime}}, m} \left(\boldsymbol{y}, \boldsymbol{x}\right) } \nonumber\\ & \qquad \times \int e^{i \left[ \Delta T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}, \boldsymbol{x}\right) \omega^{^{\prime}} + \Delta \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}, \boldsymbol{x}\right) t^{^{\prime}} \right] } W\left(\omega^{^{\prime}}, t^{^{\prime}}\right) d\omega^{^{\prime}} dt^{^{\prime}}. \end{align} \tag{ 37 }$$

We next expand $\Delta T_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}, \boldsymbol{x})$, $\Delta \mathcal V_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}, \boldsymbol{x})$, and $\Phi_{\ell,\ell^{^{\prime}}, m}$ of (37) in their Taylor series about $\boldsymbol{y} = \boldsymbol{x}$, obtaining

$$\begin{align} \mathcal L_{\ell,\ell^{^{\prime}},m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \frac{\overline{\beta_{\ell, m}\left(\boldsymbol{y}\right)} \beta_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \beta_{\ell, m}\left(\boldsymbol{x}\right) \overline{\beta_{\ell^{^{\prime}}, m}\left(\boldsymbol{x}\right) }} {|\boldsymbol{\beta} \left(\boldsymbol{y}\right)|^2 } \int e^{i \left(\boldsymbol{y} - \boldsymbol{x}\right) \cdot \boldsymbol{\Psi}_{\ell, \ell^{^{\prime}}, m}\left(\omega^{^{\prime}}, t^{^{\prime}}\right)} W\left(\omega^{^{\prime}}, t^{^{\prime}}\right) d\omega^{^{\prime}} dt^{^{\prime}}, \end{align} \tag{ 38 }$$

where

$$\begin{equation} \boxed{ \boldsymbol{\Psi}_{\ell, \ell^{^{\prime}}, m}\left(\omega^{^{\prime}}, t^{^{\prime}}\right) = \boldsymbol{k}^0_{\ell, \ell^{^{\prime}}, m} + \omega^{^{\prime}} \nabla T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) + t^{^{\prime}} \nabla \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right). } \end{equation} \tag{ 39 }$$

Here the spatial gradients of the TDOA and FDOA are given by

$$\begin{align} \nabla\left( \hbox{TDOA}\right) = \nabla T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) & = \frac{1}{c} \left( \widehat{\boldsymbol{R}}_{\ell,m}\left(\boldsymbol{y}\right) - \widehat{\boldsymbol{R}}_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \right) \nonumber\\ \nabla\left( \hbox{FDOA}\right) = \nabla \mathcal V_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) & = - \frac{\omega_0 }{c} \left( \mathcal P_{\ell, m , \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} - \mathcal P_{\ell^{^{\prime}}, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell^{^{\prime}},m} \right) \end{align} \tag{ 40 }$$

where $\boldsymbol{R}_{\ell, m}(\boldsymbol{y}) = \boldsymbol{y} - \boldsymbol{\gamma}_\ell(s_m)$, with the hat denoting unit vector ($\widehat{\boldsymbol{R}} = \boldsymbol{R} / R$), where $\dot{\boldsymbol{\gamma}}_{\ell,m} = \dot{\boldsymbol{\gamma}}_{\ell}(s_m)$ is the sensor velocity vector at slow time m, and where the scaled projection operator $\mathcal P_{\ell, m} = \nabla \widehat{\boldsymbol{R}}_{\ell,m}$ is

$$\begin{equation} \mathcal P_{\ell, m, \boldsymbol{y}} = \frac{\boldsymbol{I} - \widehat{\boldsymbol{R}}_{\ell,m} \left(\boldsymbol{y}\right) \left( \widehat{\boldsymbol{R}}_{\ell,m} \left(\boldsymbol{y}\right) \right)^\dagger}{R_{\ell,m} }. \end{equation} \tag{ 41 }$$

The matrix $\mathcal P_{\ell, m, \boldsymbol{y}}$ projects a spatial vector onto the plane perpendicular to $\widehat{\boldsymbol{R}}_{\ell,m} (\boldsymbol{y})$. The term $\boldsymbol{k}^0_{\ell, \ell^{^{\prime}}, m}$ of (39) results from the expansion of $\Phi_{\ell, \ell^{^{\prime}}, m}$ and is given by

$$\begin{align} \boldsymbol{k}^0_{\ell, \ell^{^{\prime}}, m} & = - \frac{\omega_0}{c}\left[ \tau_{\ell, m}\left(\boldsymbol{y}\right) \mathcal P_{\ell, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} - \tau_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \mathcal P_{\ell^{^{\prime}}, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell^{^{\prime}},m} \right]. \end{align} \tag{ 42 }$$

More details of the calculations are given in appendix D.

We see that the phase of (38) is of the form $(\boldsymbol{y} - \boldsymbol{x}) \cdot \boldsymbol{\Psi}_{\ell, \ell^{^{\prime}}, m}$. Consequently the main contribution to the integral of (38) comes from the point $\boldsymbol{y} = \boldsymbol{x}$.

In (39), (40), and (42), we note that the vector $\mathcal P_{\ell, m} \dot{\boldsymbol{\gamma}}_{\ell,m}$ is orthogonal to $\widehat{\boldsymbol{R}}_{\ell,m}(\boldsymbol{y})$ for every $\ell$ and m. The terms $\boldsymbol{k}^0_{\ell, \ell^{^{\prime}}, m}$ and $t^{^{\prime}} \nabla \mathcal V_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y})$ of (39) can be combined as

$$\begin{align} - \frac{\omega_0 }{c} \left[ \frac{\Pi_{\ell, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} \left(\tau_{\ell, m}\left(\boldsymbol{y}\right) + t^{^{\prime}} \right) } {R_{\ell, m}\left(\boldsymbol{y}\right)} - \frac{\Pi_{\ell^{^{\prime}}, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell^{^{\prime}},m} \left(\tau_{\ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) + t^{^{\prime}} \right) } {R_{\ell^{^{\prime}}, m}\left(\boldsymbol{y}\right)} \right] \end{align} \tag{ 43 }$$

where $\Pi_{\ell, m, \boldsymbol{y}} = \boldsymbol{I} - \widehat{\boldsymbol{R}}_{\ell, m}(\boldsymbol{y}) \widehat{\boldsymbol{R}}_{\ell, m}(\boldsymbol{y}) ^T$ is the operator that projects a vector onto the plane perpendicular to $\widehat{\boldsymbol{R}}_{\ell, m}(\boldsymbol{y})$, so that $\mathcal P = \Pi / R$. Thus each term in (43) is the ratio of the range distance R to the distance traveled perpendicular to R ; this is the ratio of two sides of a right triangle (see figure 3) and is therefore approximately the angular aperture (in radians) covered by the sensor over the mth data collection interval. Equation (43) is the difference of the corresponding angles for the sensors $\ell$ and $\ell^{^{\prime}}$.

In (38), we make the change of variables

$$\begin{align} \left(\omega^{^{\prime}}, t^{^{\prime}}\right) \mapsto \boldsymbol{k} & = \boldsymbol{\Psi}_{\ell, \ell^{^{\prime}}, m}\left(\omega^{^{\prime}}, t^{^{\prime}}\right). \end{align} \tag{ 44 }$$

The contribution to the image from sensor pair $(\ell, \ell^{^{\prime}})$ and slow time m can then be written

$$\begin{align} \mathcal L_{\ell, \ell^{^{\prime}},m} \left(\boldsymbol{y}\right) & = \frac{\overline{\beta_{\ell, m}\left(\boldsymbol{y}\right)} \beta_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \beta_{\ell, m}\left(\boldsymbol{x}\right) \overline{\beta_{\ell^{^{\prime}}, m}\left(\boldsymbol{x}\right) }} {|\boldsymbol{\beta} \left(\boldsymbol{y}\right)|^2 } \nonumber\\ & \qquad \times \int \chi_{\Omega_{\ell, \ell^{^{\prime}}, m}}\left(\boldsymbol{k}\right) e^{i \left(\boldsymbol{y} - \boldsymbol{x}\right) \cdot \boldsymbol{k}} \underbrace{W\left(\omega^{^{\prime}}\left(\boldsymbol{k}\right), t^{^{\prime}}\left(\boldsymbol{k}\right)\right) }_{\displaystyle W\left(\boldsymbol{k}\right) } J_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{k}\right) d\boldsymbol{k}, \end{align} \tag{ 45 }$$

where the Jacobian $J_{\ell, \ell^{^{\prime}}, m} = | \partial (\omega^{^{\prime}}, t^{^{\prime}})/ (\partial \boldsymbol{k})|$ is often called the 'Beylkin determinant' [39] and where $\chi_{{\Omega}_{\ell, \ell^{^{\prime}}, m}}(\boldsymbol{k})$ is the characteristic function of the set

$$\begin{equation} \Omega_{\ell, \ell^{^{\prime}}, m} = \left\{ \boldsymbol{k} = \boldsymbol{\Psi}_{\ell, \ell^{^{\prime}}, m}\left(\omega^{^{\prime}}, t^{^{\prime}} \right) : \left(\omega^{^{\prime}}, t^{^{\prime}}\right) \in \hbox{support of } W \right\} . \end{equation} \tag{ 46 }$$

The set $\Omega_{\ell, \ell^{^{\prime}}, m}$ is the set of image Fourier components, of the image near $\boldsymbol{y} = \boldsymbol{x}$, determined by measurements on sensor pair $\ell, \ell^{^{\prime}}$ at slow time index m. An example of a set $\Omega_{\ell, \ell^{^{\prime}}, m} (\boldsymbol{k})$ is shown in figure 4.

Zoom In Zoom Out Reset image size

The union of these sets,

$$\begin{equation} \Omega_{\boldsymbol{y}} = \cup_{\ell, \ell^{^{\prime}}, m} \Omega_{\ell, \ell^{^{\prime}}, m}, \end{equation} \tag{ 47 }$$

we call the data collection manifold (DCM) [35]. As we see below, resolution at y is mainly determined by the size of this DCM, which is determined by the set of values of Ψ of (39).

4.1.3. Combining data from all sensors.

The full likelihood image is the coherent sum of the images from each sensor pair and every slow time. In the neighborhood of the image point $\boldsymbol{y} = \boldsymbol{x}$, the image can be written:

$$\begin{align} \mathcal L\left(\boldsymbol{y}\right) = \int_{\Omega_{\boldsymbol{y}}} \Xi\left(\boldsymbol{y}, \boldsymbol{x}, \boldsymbol{k}\right) e^{i \left(\boldsymbol{y} - \boldsymbol{x}\right) \cdot \boldsymbol{k}} d\boldsymbol{k}, \end{align} \tag{ 48 }$$

where the amplitude Ξ is

$$\begin{equation} \Xi\left(\boldsymbol{y}, \boldsymbol{x}, \boldsymbol{k}\right) = \sum_{\ell, \ell^{^{\prime}}, m} \frac{\overline{\beta_{\ell, m}\left(\boldsymbol{y}\right)} \beta_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \beta_{\ell, m}\left(\boldsymbol{x}\right) \overline{\beta_{\ell^{^{\prime}}, m}\left(\boldsymbol{x}\right) }} {|\boldsymbol{\beta} \left(\boldsymbol{y}\right)|^2 } \chi_{\Omega_{\ell, \ell^{^{\prime}}, m}}\left(\boldsymbol{k}\right) W\left(\boldsymbol{k}\right) J_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{k}\right). \end{equation} \tag{ 49 }$$

4.1.4. Resolution formulas.

The resolution at y is determined by the data collection manifold $\Omega_{\boldsymbol{y}}$ because the image (48) at y is the Fourier transform (of Ξ) over the region $\Omega_{\boldsymbol{y}}$ in Fourier ( k -)space. For a source distribution consisting of a single isotropic point source at y , we determine the resolution of its image (the point-spread function) by noting that

$$\begin{equation} \int_a^b e^{i k x} dk \propto \hbox{sinc} \left( \frac{b-a}{2} \ x \right) \end{equation} \tag{ 50 }$$

has a main-lobe peak-to-null width of $x = 2\pi/(b-a)$. Consequently, if $\Omega_{\boldsymbol{y}}$ extends, in direction $\hat{\boldsymbol{e}}$, from $\boldsymbol{k} \cdot \hat{\boldsymbol{e}} = a$ to $\boldsymbol{k} \cdot \hat{\boldsymbol{e}} = b$, then the point-spread function will have peak-to-null resolution in direction $\hat{\boldsymbol{e}}$ of approximately $2\pi/(b-a)$.

For the geometry shown in figure 3, we refer to the direction from the sensors to the top of the page as the 'range' direction, and the direction of travel as the 'cross-range' direction. We see in figure 4 that the direction $\widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1$ is roughly in the cross-range direction.

Figure 4 is a sketch of the corresponding values of Ψ, as defined in (39). A detailed discussion of figure 4 appears in appendix E. The main point is that for a single view (i.e. a single slow time s_m ) in this geometry, bandwidth of the transmitted waveform roughly corresponds to k -space extent in the horizontal (cross-range) direction, and waveform duration roughly corresponds to k -space extent in the vertical (range) direction. With multiple views, the k -space region rotates counter-clockwise about the origin by the total viewing aperture angle α, which we define as the angular aperture swept out by the midpoint between the sensors.

The resolution for the geometry of figure 4 is obtained as follows.

• Resolution from the TDOA term of (39) (relevant for short waveforms with significant bandwidth). This analysis is the same as that in [23, 24].
  • 1  
    Cross-range resolution. From the first (TDOA) term of (39), we see that the extent in direction $\widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1$ is determined by $\omega^{^{\prime}} | \widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1| / c$, as $\omega^{^{\prime}}$ ranges over the angular frequency band of the source. Consequently the single-look peak-to-null resolution in the cross-range direction is $c/( B | \widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1| )$, where B is the source bandwidth in Hertz.The cross-range resolution is not appreciably affected by a small synthetic aperture (see figure 6).
  • 2  
    Range resolution. As the sensors move along the flight path, the vector $\widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1$ rotates. With the small-angle approximation (for the aperture α), the vertical (range) extent is approximately $2 \pi f_0 \alpha/c$, which gives rise to resolution of $c/(f_0 \alpha)$.
• Resolution from the FDOA term of (39) (relevant for waveforms of significant duration).
  • Range resolution. The FDOA term of (39) is a vector that is roughly perpendicular to $\widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1$. Its extent, as discussed below (43), is roughly $\omega_0 (\theta_2 - \theta_1)/c$, where $\theta_2 - \theta_1$ is the difference of angular apertures covered by the sensors. Consequently the single-look peak-to-null resolution in the range direction is approximately $c/[f_0(\theta_2 - \theta_1)]$, where f₀ is the center frequency in Hertz.Again, as the sensors move along the flight path to angle α, the single-look contributions to the DCM rotate to the angle α. With the small-angle approximation for α, the total vertical (range) extent of the DCM is roughly $\omega_0[ \alpha + (\theta_2 - \theta_1)]/c$. Consequently the range resolution is $c/[f_0 (\alpha + \theta_2 - \theta_1)]$.
  • Cross-range resolution. The cross-range resolution deriving from the FDOA term arises mainly from the dwell and the synthetic aperture. The horizontal (cross-range) extent is approximately $\alpha \omega_0 (\theta_2 - \theta_1)/c$, which is the single-look vertical extent multiplied by $\sin \alpha \approx \alpha$. Consequently the cross-range resolution is approximately $c /[f_0 \alpha (\theta_2 - \theta_1)]$.

These formulas are listed in table 1. Data collection geometry significantly different from that of figure 3 may result in different formulas. For example, one stationary sensor and one moving sensor might provide desirable resolution properties.

Table 1. Approximate resolution for geometry of figure 3. B is the bandwidth in Hertz and α is the total viewing aperture.

	Range	Cross-range
Single look	$\displaystyle \frac{ c }{f_0 (\theta_2 - \theta_1) }$	$\displaystyle \frac{c}{B | \widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1|}$
Aperture α	min$\left\{\displaystyle \frac{ c }{f_0 (\alpha + \theta_2 - \theta_1) }, \frac{c}{\alpha f_0} \right\}$	min $\left\{ \displaystyle \frac{c}{B | \widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1|} , \displaystyle \frac{ c }{f_0 \alpha (\theta_2 - \theta_1) } \right\}$

4.1.5. Resolution examples

4.1.5.1. Example: an airport surveillance radar or weather radar.

We consider the resolution for a source consisting of an airport surveillance radar operating at the carrier frequency 2.8 GHz, and transmitting pulses with bandwidth $B =$ 10 MHz.

The cross-range (TDOA-derived) resolution is roughly

$$\begin{equation} \frac{c}{B} = \frac{3 \times 10^8 \,\hbox{m}\,\hbox{s}^{-1}}{10^7 /\hbox{s} } = 30 \,\hbox{m}. \end{equation} \tag{ 51 }$$

where we have assumed that the sensor geometry is such that $| \widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1| \approx 1$. Broader-band sources provide correspondingly better resolution; sensors too close together give worse resolution.

The range (FDOA-derived) resolution, on the other hand, for data collection geometry that results in a one-degree aperture $\theta_2 - \theta_1$, is approximately

$$\begin{equation} \frac{ c }{f_0 \left(\theta_2 - \theta_1\right) } = \frac{3 \times 10^8 \,\hbox{m}\,\hbox{s}^{-1} }{2.8 \times 10^9 /\hbox{s} \times 1^\circ} \times \frac{180^\circ}{\pi} \approx 7 \,\hbox{m}. \end{equation} \tag{ 52 }$$

An airport surveillance radar may have a narrow beam and may be rotating, so a one-degree aperture could be difficult to obtain.

4.1.5.2. Example: an X-band radar.

For an X-band radar with a 100 MHz bandwidth and center frequency of $f_0 = 10^{10}$, the cross-range resolution is approximately

$$\begin{equation} \frac{c}{B} = \frac{3 \times 10^8 \hbox{m s}^{-1} }{10^8 /\hbox{sec} } = 3 \,\hbox{m} \end{equation} \tag{ 53 }$$

and the range resolution, again for a data collection resulting in a one-degree aperture, is roughly

$$\begin{equation} \frac{ c }{f_0 \left(\theta_2 - \theta_1\right) } = \frac{3 \times 10^8 \hbox{m s}^{-1} }{10^{10} /\hbox{s} \times 1^\circ} \times \frac{180^\circ}{\pi} \approx 1.8 \,\hbox{m}. \end{equation} \tag{ 54 }$$

4.1.5.3. Example: 16 kHz bandwidth, long-duration waveform.

Communications waveforms are often long-duration and have a narrow bandwidth. In figure 5, we show the DCM and PSF for a source located at the origin, transmitting pulses of duration .08 second with (angular) center frequency $\omega_0 = 10^{7}$ radians s⁻¹ (≈ 1.6 MHz) and bandwidth of 10⁵ radians s⁻¹ (≈ 16 kHz), with sensors 50 m to the south [initially located at (−20 m, −50 m) and (20 m, −50 m)] traveling at speed 10 m s⁻¹ along the x axis, with the geometry as shown in figures 3 and 4. (These parameters are chosen for ease of computation and for visibility of the relevant features.)

Zoom In Zoom Out Reset image size

The top left plot in figure 5 is the DCM for a single look from the initial position; immediately below is the characteristic function of the DCM, which is the integrand of (57), and below that is the corresponding normalized PSF plotted with the vertical axis on a dB scale. The top right plot is the DCM for 10 looks, at slow times 0, 0.1111, 0.2222, 0.3333, 0.4444, 0.5556, 0.6667, 0.7778, 0.8889, and 1.0, as the sensors move along their trajectories. We see that the successive single-look contributions to the DCM rotate and shift counter-clockwise with successive looks. Immediately below the 10-look DCM is the characteristic function of that DCM, and below that is the corresponding normalized PSF plotted with the vertical axis on a dB scale. The units for the k_x and k_y axes are 1/m.

4.1.5.4. Example: 250 MHz bandwidth, short-duration waveform.

Figure 6 shows a DCM for a short-pulse waveform with bandwidth of 250 MHz around a center frequency of 2 GHz, again for the same geometry as shown in figure 3. (Again the parameters are chosen for visibility.) The different colors show the contributions for each of the 40 different looks (i.e. slow time samples). Here the contributions from different looks have minimal overlap.

Zoom In Zoom Out Reset image size

The resolution analysis here is similar to that of [23]. It is also similar to the analysis of standard synthetic-aperture radar (SAR) [35] except that the $2\widehat{\boldsymbol{R}}$ in SAR is replaced here by $\widehat{\boldsymbol{R}}_2 - \widehat{\boldsymbol{R}}_1$, which is roughly perpendicular to the range direction $\widehat{\boldsymbol{R}}$.

4.2. Optimal filtering of the image

In this subsection, we follow the general strategy of [35, 40].

To the detection-based image (48), we apply a filter with kernel Q:

$$\begin{align} \mathcal I\left(\boldsymbol{p}\right) & = \int \!\!\! \int e^{i \left(\boldsymbol{p} - \boldsymbol{y}\right) \cdot \boldsymbol{\xi}} Q\left(\boldsymbol{p}, \boldsymbol{\xi}\right) d\boldsymbol{\xi} \ \mathcal L \left(\boldsymbol{y}\right) d\boldsymbol{y} \nonumber\\ & = \int \!\!\! \int e^{i \left(\boldsymbol{p} - \boldsymbol{y}\right) \cdot \boldsymbol{\xi}} Q\left(\boldsymbol{p}, \boldsymbol{\xi}\right) d\boldsymbol{\xi} \int_\Omega \Xi\left(\boldsymbol{y}, \boldsymbol{x}, \boldsymbol{k}\right) e^{i \left(\boldsymbol{y} - \boldsymbol{x}\right) \cdot \boldsymbol{k}} d\boldsymbol{k} \ d\boldsymbol{y} \nonumber\\ &\propto \int_\Omega e^{i \left(\boldsymbol{p} - \boldsymbol{x}\right) \cdot \boldsymbol{k}} Q\left(\boldsymbol{p}, \boldsymbol{k}\right) \Xi\left(\boldsymbol{x}, \boldsymbol{x}, \boldsymbol{k}\right) d\boldsymbol{k}, \end{align} \tag{ 55 }$$

where we have used $\int \exp(-i \boldsymbol{y} \cdot (\boldsymbol{\xi} - \boldsymbol{k})) d\boldsymbol{y} \propto \delta(\boldsymbol{\xi} - \boldsymbol{k})$. We see that to make the image as close as possible to the Dirac delta function $\delta(\boldsymbol{p}\! - \! \boldsymbol{x}) \propto \int \exp(i (\boldsymbol{p} - \boldsymbol{x}) \cdot \boldsymbol{k}) d\boldsymbol{k}$, we should choose Q to be the filter

$$\begin{equation} Q\left(\boldsymbol{p}, \boldsymbol{k}\right) = \frac{\chi_\Omega\left(\boldsymbol{k}\right) }{\Xi \left(\boldsymbol{p}, \boldsymbol{p}, \boldsymbol{k}\right) } = \frac{\chi_\Omega\left(\boldsymbol{k}\right) }{\sum_{\ell, \ell^{^{\prime}}, m} \frac{g_{\ell ,m} \left(\boldsymbol{p}\right) \overline{g_{\ell^{^{\prime}} ,m} \left(\boldsymbol{p}\right) } } {\sigma_\ell^2 \sigma_\ell^{^{\prime}} 2} \chi_{\Omega_{\ell, \ell^{^{\prime}}, m}}\left(\boldsymbol{k}\right) W\left(\boldsymbol{k}\right) J_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{k}\right) }. \end{equation} \tag{ 56 }$$

Here the $\chi_\Omega$ in the numerator is 1 in Ω and 0 otherwise. We note that this filter depends, in a slowly varying manner, on the image point p through the functions $g_{\ell, m} ,\ \chi_{\Omega_{\ell,\ell^{^{\prime}},m}}$ and the Jacobian $J_{\ell, \ell^{^{\prime}},m}$. Typically the Jacobian involves a factor of $1/\omega$ that produces image sharpening.

With filter (56), the image (55) becomes, to leading order,

$$\begin{equation} \mathcal I^Q\left(\boldsymbol{p}\right) = \int_{\Omega_{\boldsymbol{y}}} e^{i \left(\boldsymbol{p} - \boldsymbol{x}\right) \cdot \boldsymbol{k}} d\boldsymbol{k} \end{equation} \tag{ 57 }$$

which is an approximation to a delta function located at $\boldsymbol{p} = \boldsymbol{x}$.

Thus we see that source localization is reduced to geometry and Fourier analysis. The flight trajectories, source location, and transmitted waveforms determine the set Ω, and the Fourier transform of the characteristic function of the set Ω determines the image of a single source. The connection between sets in Fourier space and associated images has been studied in [41, 42].

For k -space regions of overlap for multiple sensors and/or multiple looks m, (56) provides automatic downweighting. This is studied further in the following subsection.

4.3. Redundant information

We see in figure 7 that even a single pair of sensors can provide information that is redundant, in the sense that the k -space regions from different looks can overlap. The analysis above suggests that these overlapping regions should be downweighted by means of the filtering operation (55); but how important is this really? In figure 7 we compare the weighted PSF with the unweighted PSF for the 10-look, 1.6 kHz-bandwidth case. The weighted PSF (bottom right ) is nearly indistinguishable from the unweighted PSF (middle row, right). A close examination shows that the weighted PSF is slightly better than the unweighted one.

Zoom In Zoom Out Reset image size

The conclusion is that, at least in this case, it is the size of the k -space region that is most important for resolution, not the weighting. (And this is consistent with other numerical experiments, not shown.) In other words, data that overlaps in k -space can often be simply added together coherently, with little danger that the overlapping regions will result in appreciably degraded point-spread functions. The more overlap, of course, the more important weighting could become.

4.4. Image analysis for multiple sources

Although the image formation process was derived under the assumption that only a single source is present, in fact there may be multiple sources. This section investigates the behavior of the imaging algorithm when multiple sources are present. To do this, we insert the signal model (11) into (24):

$$\begin{align} \mathcal L_{\ell,\ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) & = \frac{\overline{g_{\ell, m}\left(\boldsymbol{y}\right)} g_{\ell^{^{\prime}},m} \left(\boldsymbol{y}\right) } { \sum_{\ell^{^{\prime}\!\!^{\prime}}} |g_{\ell^{^{\prime}\!\!^{\prime}}, m} / \sigma_{\ell^{^{\prime}\!\!^{\prime}}}|^2 \sigma_\ell^2 \sigma_{\ell^{^{\prime}}}^2 } \nonumber\\ & \quad \times \int \sum_j g_{\ell ,m} \left(\boldsymbol{x}_j\right) e^{i \Delta \nu_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}_j\right) \left( t + \tau_{\ell , m} \left(\boldsymbol{y}\right)\right) } \ a_{\boldsymbol{x}_j,m} \left( t + \Delta\tau_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}_j\right) \right)\nonumber\\ & \quad \times \sum_{j^{^{\prime}}} \overline{g_{\ell^{^{\prime}} ,m} \left(\boldsymbol{x}_{j^{^{\prime}}} \right) e^{i \Delta \nu_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}\right) \left( t + \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right)\right) } \ a_{\boldsymbol{x}_{j^{^{\prime}}}, m} \left( t + \Delta\tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}\right) \right) } dt\nonumber \end{align} \tag{ 58 }$$

where as before

$$\begin{align} \Delta\nu_{\ell, m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \nu_{\ell , m} \left(\boldsymbol{y}\right) - \nu_{\ell , m} \left(\boldsymbol{x}\right) \nonumber\\ \Delta \tau_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}\right) & = \tau_{\ell , m} \left(\boldsymbol{y}\right)- \tau_{\ell , m} \left(\boldsymbol{x}\right). \end{align} \tag{ 59 }$$

In (58), we make the change of variables $t \mapsto t^{^{\prime}} = t + \Delta\tau_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}})$. This results in

$$\begin{align} \Lambda_{\ell,\ell^{^{\prime}}, m}(\boldsymbol{y}) & = \sum_{j,j^{^{\prime}}} \frac{g_{\ell ,m} (\boldsymbol{x}_j) \overline{g_{\ell^{^{\prime}} ,m} (\boldsymbol{x}_{j^{^{\prime}}} ) } \overline{g_{\ell, m}(\boldsymbol{y})} g_{\ell^{^{\prime}},m} (\boldsymbol{y}) } { \sum_{\ell^{^{\prime}\!\!^{\prime}}} |g_{\ell^{^{\prime}\!\!^{\prime}}, m} / \sigma_{\ell^{^{\prime}\!\!^{\prime}}}|^2 \sigma_\ell^2 \sigma_{\ell^{^{\prime}}}^2 } \int e^{i \Delta \nu_{\ell,m} (\boldsymbol{y}, \boldsymbol{x}_j) ( t ^{^{\prime}} - \Delta\tau_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) + \tau_{\ell , m} (\boldsymbol{y})) } \nonumber\\ &\quad \times \ a_{\boldsymbol{x}_j,m} \big( t^{^{\prime}} - \Delta\tau_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) + \Delta\tau_{\ell,m} (\boldsymbol{y}, \boldsymbol{x}_j) \big) \nonumber\\ &\quad \times e^{-i \Delta \nu_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) ( t^{^{\prime}} - \Delta\tau_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) + \tau_{\ell^{^{\prime}} , m} (\boldsymbol{y})) } \ \overline{a_{\boldsymbol{x}_{j^{^{\prime}}}, m} \big( t^{^{\prime}} ) \big) } dt^{^{\prime}} \nonumber\\ & = \sum_{j,j^{^{\prime}}} \frac{g_{\ell ,m} (\boldsymbol{x}_j) \overline{g_{\ell^{^{\prime}} ,m} (\boldsymbol{x}_{j^{^{\prime}}} ) } \overline{g_{\ell, m}(\boldsymbol{y})} g_{\ell^{^{\prime}},m} (\boldsymbol{y}) } { \sum_{\ell^{^{\prime}\!\!^{\prime}}} |g_{\ell^{^{\prime}\!\!^{\prime}}, m} / \sigma_{\ell^{^{\prime}\!\!^{\prime}}}|^2 \sigma_\ell^2 \sigma_{\ell^{^{\prime}}}^2 } e^{i \Phi_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}, \boldsymbol{x}_j, \boldsymbol{x}_{j^{^{\prime}}}) } \int e^{i \left[ \Delta \nu_{\ell,m} (\boldsymbol{y}, \boldsymbol{x}_j) - \Delta \nu_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) \right] t^{^{\prime}} } \nonumber\\ &\quad \times \ a_{\boldsymbol{x}_j,m} \big( t^{^{\prime}} - \Delta\tau_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) + \Delta\tau_{\ell,m} (\boldsymbol{y}, \boldsymbol{x}_j) \big) \ \overline{a_{\boldsymbol{x}_{j^{^{\prime}}}, m} \big( t^{^{\prime}} ) \big) } dt^{^{\prime}} \nonumber\\ & = \sum_{j,j^{^{\prime}}} \frac{g_{\ell ,m} (\boldsymbol{x}_j) \overline{g_{\ell^{^{\prime}} ,m} (\boldsymbol{x}_{j^{^{\prime}}} ) } \overline{g_{\ell, m}(\boldsymbol{y})} g_{\ell^{^{\prime}},m} (\boldsymbol{y}) } { \sum_{\ell^{^{\prime}\!\!^{\prime}}} |g_{\ell^{^{\prime}\!\!^{\prime}}, m} / \sigma_{\ell^{^{\prime}\!\!^{\prime}}}|^2 \sigma_\ell^2 \sigma_{\ell^{^{\prime}}}^2 } e^{i \Phi_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}, \boldsymbol{x}_j, \boldsymbol{x}_{j^{^{\prime}}}) } \nonumber\\ &\quad \times \mathcal A_{\boldsymbol{x}_j, \boldsymbol{x}_{j^{^{\prime}}} } \left(\Delta\tau_{\ell,m} (\boldsymbol{y}, \boldsymbol{x}_j) - \Delta\tau_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) \ , \ \Delta \nu_{\ell,m} (\boldsymbol{y}, \boldsymbol{x}_j) - \Delta \nu_{\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}) \right) \end{align} \tag{ 60 }$$

where

$$\begin{equation} \mathcal A_{\boldsymbol{x} , \boldsymbol{x}^{^{\prime}}} \left(\tau, \nu\right) = \int a_{\boldsymbol{x}} \left(t+\tau\right) \overline{a_{\boldsymbol{x}^{^{\prime}}} \left(t\right)} e^{i \nu t} dt \end{equation} \tag{ 61 }$$

is the narrowband cross-ambiguity function, and where

$$\begin{align} \Phi_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{y}, \boldsymbol{x}_j, \boldsymbol{x}_{j^{^{\prime}}}\right) & = \Delta \nu_{\ell,m} \left(\boldsymbol{y}, \boldsymbol{x}_j\right) \left( - \Delta\tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}\right) + \tau_{\ell , m} \left(\boldsymbol{y}\right)\right) \nonumber\\ & \qquad - \Delta \nu_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}\right) \left( - \Delta\tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}, \boldsymbol{x}_{j^{^{\prime}}}\right) + \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right)\right). \end{align} \tag{ 62 }$$

Thus if the waveforms emitted at x _j and $\boldsymbol{x}_{j^{^{\prime}}}$ have a small cross-ambiguity function, the contribution to the image of the associated cross terms will be small.

This paper has developed a statistics-based synthetic-aperture approach for forming images that show the locations of stationary sources radiating (electromagnetic or acoustic) waves, using measurements from multiple moving passive receivers. This approach is based on the detection theory of [10, 13, 14]. The resulting imaging method is nonlinear; consequently the standard linear theory cannot be used to analyze the imaging performance.

This paper, however, has also developed a method to determine the image resolution obtained from this nonlinear approach. The resolution is obtained from the data collection manifold in k -space, which gives the set of scene Fourier coefficients that can be obtained from the data. This resolution depends not only on the geometry, but also on the radiated waveforms, which are initially unknown. The resolution from multiple views is obtained by taking the union of the associated regions in k -space.

This paper also derived the optimal filtering for the image. This optimal filtering automatically downweights regions of k -space for which redundant measurements are obtained.

Finally, this paper also addressed the case of multiple sources; the image artifacts, which are due to undesirable cross terms in the nonlinear image formation process, are determined by the cross-ambiguity function of the associated waveforms. Thus, for example, a synthetic-aperture image formed of sources transmitting orthogonal waveforms will show only negligible artifacts.

The analysis in this paper may enable source-localization systems to find sources that radiate a variety of signals, including those that might radiate long-duration signals in an extremely narrow frequency band. In addition, the resolution analysis may be useful for the design and trajectory planning of systems to achieve a desired location accuracy.

We leave many related problems for the future, including addressing sources that do not radiate isotropically, source waveforms that do not satisfy the narrowband assumption, artifacts due to multiple sources radiating correlated or identical waveforms, different choices of slow times and window lengths, interference from multiple radiating sources, different noise models, use of polarimetric data, and determination of various system parameters that may be poorly known, such as the sensor positions, unknown clock offsets, and noise variances. Certain of these issues may be especially important in different specific applications. We also leave for the future the development of alternative approaches that might exploit sparsity of the sources.

We are grateful to Drs James A Given, L Todd McWhorter, and Ivars Kirsteins for many valuable discussions. This work was supported in part by the Office of Naval Research under Contract N00014-21-1-2145 and by the Air Force Office of Scientific Research ⁶ under the Contract FA9550-21-1-0169.

No new data were created or analysed in this study.

A.1. Approach of [36]

A.2. Estimation of the waveform a_m from discrete samples of $\eta_{\ell,m}$

For each slow time m and each hypothesized location y , we consider the time-sampled version of (14), which we write as

$$\begin{equation} z_{\ell,j} = g_\ell a_j + n_{\ell,j} , \qquad j = 1, 2, \ldots , J = TB, \end{equation} \tag{ A.1 }$$

where $z_{\ell,j} = \eta_{\ell, m}(\boldsymbol{y}, j/B)$, $g_\ell = g_{\ell, m}(\boldsymbol{y})$, $a_j = \ a_{\boldsymbol{y},m}(j/B)$, and $n_{\ell,j} = n_{\ell}(j/B) \sim \mathbb{C}\mathcal N( 0, \sigma_\ell^2)$ with $n_{\ell}(t)$ the noise on the $\ell$th receiver. We assemble all the elements $z_{\ell,j}$ first into a column vector $\boldsymbol{z}_j = (z_{1,j}, \ldots , z_{L,j})^T$ and then into a matrix $\boldsymbol{Z} = (\boldsymbol{z}_1, \ldots, \boldsymbol{z}_J)$ and assemble the elements $g_\ell$ into a column vector g . Then the probability density function for Z is

$$\begin{equation} p\left(\boldsymbol{Z}\right) = \frac{1}{\left(\pi \sigma^2\right)^J } e^{- \sum_j \left(\boldsymbol{z}_j - \boldsymbol{g} a_j \right)^\dagger \mathcal R^{-1} \left(\boldsymbol{z}_j - \boldsymbol{g} a_j \right)} \end{equation} \tag{ A.2 }$$

which is a more conventional discretized version of the mth factor of (19). To find a , we maximize (A.2) with respect to the variables $\boldsymbol{a} = (a_1, \ldots, a_J)$, obtaining the pseudo-inverse $\boldsymbol{a} = (\boldsymbol{g}^\dagger \mathcal R^{-1} \boldsymbol{g})^{-1} \boldsymbol{g}^\dagger \mathcal R^{-1} \boldsymbol{z}$, which is a discrete version of (18). Substituting this result into (A.2) gives

$$\begin{equation} p\left(\boldsymbol{Z}\right) = \frac{1}{\left(\pi \sigma^2\right)^J } e^{- \sum_j \left[\left( I - P_{\boldsymbol{g}} \right)\boldsymbol{z}_j \right]^\dagger \mathcal R^{-1} \left( I - P_{\boldsymbol{g}} \right)\boldsymbol{z}_j } \end{equation} \tag{ A.3 }$$

where $P_{\boldsymbol{g}} \boldsymbol{z} = \boldsymbol{g}(\boldsymbol{g}^\dagger \mathcal R^{-1} \boldsymbol{g})^{-1} \boldsymbol{g}^\dagger \mathcal R^{-1} \boldsymbol{z}$ and I denotes the L × L identity matrix. A discrete version of the log-likelihood ratio (20) is then

$$\begin{align} & - \mathrm{tr} \sum_j \left[\left( I - P_{\boldsymbol{g}} \right)\boldsymbol{z}_j \right]^\dagger \mathcal R^{-1} \left( I - P_{\boldsymbol{g}} \right)\boldsymbol{z}_j + \mathrm{tr} \sum_j \boldsymbol{z}_j^\dagger \mathcal R^{-1} \boldsymbol{z}_j \nonumber\\ & \quad = - \sum_j \mathrm{tr} \left[ \left( I - P_{\boldsymbol{g}} \right) \mathcal R^{-1/2} \boldsymbol{z}_j \boldsymbol{z}_j^\dagger \mathcal R^{-1/2} \right] + \sum_j \mathrm{tr} \left[ \mathcal R^{-1/2} \boldsymbol{z}_j \boldsymbol{z}_j^\dagger \mathcal R^{-1/2} \right] \nonumber\\ & \quad = \sum_j \mathrm{tr} \left[ P_{\boldsymbol{g}} \mathcal R^{-1/2} \boldsymbol{z}_j \boldsymbol{z}_j^\dagger \mathcal R^{-1/2} \right]. \end{align} \tag{ A.4 }$$

At this stage we can take the limit $B \rightarrow \infty$ to obtain (21) or (B.4).

With the estimate (18) for the waveform, the first term of (20) is

$$\begin{align} &\int \left( \boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right)^\dagger {\mathcal R}^{-1} \left(\boldsymbol{\eta}_m - \boldsymbol{g}_m a_m\right) dt\nonumber\\ &\quad = \int \left[ \boldsymbol{\eta}_m - \boldsymbol{g}_m \left(\boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{g}_m \right)^{-1} \left( \boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m \right) \right]^\dagger {\mathcal R}^{-1} \left[\boldsymbol{\eta}_m - \boldsymbol{g}_m \left(\boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{g}_m\right)^{-1} \left( \boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m \right) \right] dt \nonumber\\ &\quad = \int \left[ {\mathcal R}^{-1/2} \boldsymbol{\eta}_m - {\mathcal R}^{-1/2}\boldsymbol{g}_m \left(\boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{g}_m \right)^{-1} \left( \boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m \right) \right]^\dagger \nonumber\\ & \qquad \times \left[{\mathcal R}^{-1/2}\boldsymbol{\eta}_m - {\mathcal R}^{-1/2}\boldsymbol{g} \left(\boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{g}_m\right)^{-1} \left( \boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m \right) \right] dt\nonumber\\ & \quad = \int | \left( {\boldsymbol{I}} - \boldsymbol{P}_m \right) {\mathcal R}^{-1/2} \boldsymbol{\eta}_m |^2 dt, \end{align} \tag{ B.1 }$$

where $| \cdot |$ denotes the usual L-dimensional Euclidean norm, where

$$\begin{equation} \boldsymbol{P}_m = {\mathcal R}^{-1/2}\boldsymbol{g}_m \left(\boldsymbol{g}_m^\dagger {\mathcal R}^{-1} \boldsymbol{g}_m \right)^{-1} \left( \boldsymbol{g}_m^\dagger {\mathcal R}^{-1/2} \right) = \boldsymbol{\beta}_m \left(\boldsymbol{\beta}_m^\dagger \boldsymbol{\beta}_m\right)^{-1} \boldsymbol{\beta}_m^\dagger = \frac{\boldsymbol{\beta}_m \boldsymbol{\beta}_m^\dagger} {| \boldsymbol{\beta}_m|^2} , \end{equation} \tag{ B.2 }$$

where $\boldsymbol{\beta}_m = {\mathcal R}^{-1/2}\boldsymbol{g}_m = (g_{1,m}/\sigma_1, \ldots g_{L,m}/\sigma_L)$, and $| \boldsymbol{\beta}_m|^2 = \boldsymbol{\beta}_m^\dagger \boldsymbol{\beta}_m$. We note that $\boldsymbol{\beta}_m$, P , and $\mathcal R$ are independent of t. Equation (B.1) can also be found in [10, 14].

Since $\boldsymbol{I} - \boldsymbol{P}_m$ is an orthogonal projection and is thus self-adjoint and idempotent, the integrand of (B.1) can be written [10, 14] as

$$\begin{align} | \left( {\boldsymbol{I}} - \boldsymbol{P}_m \right) {\mathcal R}^{-1/2} \boldsymbol{\eta}_m |^2 & = \left[ \left({\boldsymbol{I}} - \boldsymbol{P}_m \right) {\mathcal R}^{-1/2} \boldsymbol{\eta}_m \right]^\dagger \left({\boldsymbol{I}} - \boldsymbol{P}_m \right) {\mathcal R}^{-1/2} \boldsymbol{\eta}_m \nonumber\\ & = \hbox{tr } \left[ \left( {\mathcal R}^{-1/2} \boldsymbol{\eta}_m \right)^\dagger \left({\boldsymbol{I}} - \boldsymbol{P}_m \right) {\mathcal R}^{-1/2} \boldsymbol{\eta}_m \right] \nonumber\\ & = \hbox{tr } \left[ \left({\boldsymbol{I}} - \boldsymbol{P}_m \right) \underbrace{{\mathcal R}^{-1/2} \boldsymbol{\eta}_m \boldsymbol{\eta}_m^\dagger {\mathcal R}^{-1/2} }_{\displaystyle \boldsymbol{M}_m } \right], \end{align} \tag{ B.3 }$$

where in the third line we have used the cyclic property of trace.

Using (B.3) in (20), and writing $\boldsymbol{\eta}_m^\dagger {\mathcal R}^{-1} \boldsymbol{\eta}_m = \hbox{tr } \left[ {\mathcal R}^{-1/2} \boldsymbol{\eta}_m \boldsymbol{\eta}_m^\dagger {\mathcal R}^{-1/2} \right]$, we obtain

$$\begin{align} \mathcal L_m\left(\boldsymbol{y}\right) & = \int \hbox{tr } \left[ \boldsymbol{P}_m \boldsymbol{M}_m\left(t\right) \right] dt = \int \frac{\hbox{tr } \left[ \boldsymbol{\beta}_m \boldsymbol{\beta}_m^\dagger \boldsymbol{M}_m\left(t\right) \right] } {| \boldsymbol{\beta}_m|^2} dt = \frac{\hbox{tr } \left[ \int \boldsymbol{\beta}_m^\dagger \boldsymbol{M}_m\left(t\right) \boldsymbol{\beta}_m dt \right] } {| \boldsymbol{\beta}_m|^2} , \end{align} \tag{ B.4 }$$

where $\boldsymbol{M}_m(t)$ is the L × L self-adjoint matrix

$$\begin{align} \boldsymbol{M} _m\left(t\right) & = {\mathcal R}^{-1/2} \boldsymbol{\eta}_m \left(t\right) \boldsymbol{\eta}_m^\dagger \left(t\right) {\mathcal R}^{-1/2} . \end{align} \tag{ B.5 }$$

We note that in the final expression of (B.4), we can drop the trace because the quantity inside the square brackets is a scalar.

In this section we write the inner product of (24) directly in terms of receiver data.

$$\begin{align} \langle \eta_{\ell, m} , \eta_{\ell^{^{\prime}}, m} \rangle & = \int e^{i \nu_{\ell, m} (\boldsymbol{y})( t + \tau_{\ell, m} (\boldsymbol{y})) } d_{\ell, m}^{BB} (t + \tau_{\ell, m} (\boldsymbol{y})) \nonumber\\ & \qquad \times e^{-i \nu_{\ell^{^{\prime}}, m} (\boldsymbol{y}) ( t + \tau_{\ell^{^{\prime}}, m} (\boldsymbol{y})) } \overline{d_{\ell^{^{\prime}}, m}^{BB} (t + \tau_{\ell^{^{\prime}}, m} (\boldsymbol{y})) }\ dt. \end{align} \tag{ C.1 }$$

We make the change of variables $t \mapsto t^{^{\prime}} = t + \tau_{\ell^{^{\prime}}, m} (\boldsymbol{y})$, which is inverted as $t = t^{^{\prime}} - \tau_{\ell^{^{\prime}}, m} (\boldsymbol{y})$:

$$\begin{align} \langle \eta_{\ell, m} , \eta_{\ell^{^{\prime}}, m} \rangle & = \int d_{\ell, m}^{BB} \left(t^{^{\prime}} - \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) \overline{d_{\ell^{^{\prime}}, m}^{BB} \left(t^{^{\prime}}\right) } \nonumber\\ &\qquad \times e^{i \left[\nu_{\ell, m} \left(\boldsymbol{y}\right)\left( t^{^{\prime}} - \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) + \tau_{\ell, m} \left(\boldsymbol{y}\right)\right) - \nu_{\ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) t^{^{\prime}} \right] } dt^{^{\prime}} \nonumber\\ & = \underbrace{\int d_{\ell, m}^{BB} \left(t^{^{\prime}} + T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) \right) \overline{d_{\ell^{^{\prime}}, m}^{BB} \left(t^{^{\prime}}\right) } e^{i t^{^{\prime}} \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) } dt^{^{\prime}} }_{\displaystyle \mathcal A_{\ell, \ell^{^{\prime}}} \left( T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) , \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) \right) } \ e^{-i \nu_{\ell, m} T_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right)} \end{align} \tag{ C.2 }$$

where $\mathcal A$ is the cross-ambiguity function of (26) and where T and $\mathcal V$ are the TDOA and FDOA as defined in (27).

The expression for Φ is obtained as

$$\begin{align} \Phi_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}, \boldsymbol{x}\right) & = \Delta \nu_{\ell,m} \tau_{\ell , m} \left(\boldsymbol{y}\right) - \Delta \nu_{\ell^{^{\prime}},m} \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right) + \Delta\tau_{\ell^{^{\prime}},m} \left(\Delta \nu_{\ell^{^{\prime}},m} - \Delta \nu_{\ell,m} \right) \nonumber\\ & = \Delta \nu_{\ell,m}\tau_{\ell , m} \left(\boldsymbol{y}\right) - \Delta \nu_{\ell^{^{\prime}},m} \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right) + \left[ \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right)- \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{x}\right) \right] \left(\Delta \nu_{\ell^{^{\prime}},m} - \Delta \nu_{\ell,m} \right) \nonumber\\ & = \Delta \nu_{\ell,m} \left[ \tau_{\ell , m} \left(\boldsymbol{y}\right)- \left[ \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right)- \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{x}\right) \right] \right] \nonumber\\ & \qquad + \Delta \nu_{\ell^{^{\prime}},m} \left[ - \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right) + \left[ \tau_{\ell^{^{\prime}} , m} \left(\boldsymbol{y}\right)- \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{x}\right) \right] \right] \nonumber\\ & = \left[ \nu_{\ell , m} \left(\boldsymbol{y}\right)-\nu_{\ell, m} \left(\boldsymbol{x}\right)\right]\left[ T_{\ell,\ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) + \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{x}\right) \right] - \left[ \nu_{\ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) - \nu_{\ell^{^{\prime}} , m} \left(\boldsymbol{x}\right) \right] \tau_{\ell^{^{\prime}}, m} \left(\boldsymbol{x}\right) \nonumber\\ & = \left[ \nu_{\ell , m} \left(\boldsymbol{y}\right) - \nu_{\ell , m} \left(\boldsymbol{x}\right) \right] T_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}) \nonumber\\ &\qquad + \tau_{\ell^{^{\prime}}, m} (\boldsymbol{x}) \left[ \nu_{\ell , m} (\boldsymbol{y}) - \nu_{\ell , m} (\boldsymbol{x}) - \left[ \nu_{\ell^{^{\prime}}, m} (\boldsymbol{y}) - \nu_{\ell^{^{\prime}} , m} (\boldsymbol{x}) \right] \right] \nonumber\\ & = \left[\nu_{\ell, m}(\boldsymbol{y})-\nu_{\ell, m} (\boldsymbol{x})\right]T_{\ell, \ell^{^{\prime}},m}(\boldsymbol{y})+ \tau_{\ell^{^{\prime}}, m} (\boldsymbol{x}) \underbrace{\left[ \mathcal V_{\ell,\ell^{^{\prime}},m} (\boldsymbol{y}) - \mathcal V_{\ell,\ell^{^{\prime}},m} (\boldsymbol{x}) \right] }_{\displaystyle \Delta \mathcal V_{\ell,\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x}) } . \end{align} \tag{ D.1 }$$

Using Taylor expansions for $\nu_{\ell , m} (\boldsymbol{y}) - \nu_{\ell , m} (\boldsymbol{x})$ and $\Delta \mathcal V_{\ell,\ell^{^{\prime}},m} (\boldsymbol{y}, \boldsymbol{x})$, we obtain

$$\begin{align} \Phi_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}, \boldsymbol{x}) & = T_{\ell, \ell^{^{\prime}}, m} (\boldsymbol{y}) (\boldsymbol{y} - \boldsymbol{x}) \cdot \nabla \nu_{\ell, m} \nonumber\\ & \qquad + \left[(\tau_{\ell^{^{\prime}},m}(\boldsymbol{y}) + (\boldsymbol{y} - \boldsymbol{x}) \cdot \nabla \tau_{\ell^{^{\prime}},m}(\boldsymbol{y}) + \cdots\right] \left[ (\boldsymbol{y} - \boldsymbol{x}) \cdot \nabla \mathcal V_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}) + \cdots \right] \nonumber\\ & = (\boldsymbol{y} - \boldsymbol{x}) \cdot \boldsymbol{k}^0_{\ell, \ell^{^{\prime}}, m}(\boldsymbol{y}) + \cdots , \end{align} \tag{ D.2 }$$

where, having dropped terms that are quadratic in $|\boldsymbol{y} - \boldsymbol{x}|$, we have

$$\begin{align} \boldsymbol{k}^0_{\ell, \ell^{^{\prime}}, m}&\left(\boldsymbol{y}\right) = T_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) \nabla \nu_{\ell, m} + \tau_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \nabla \mathcal V_{\ell, \ell^{^{\prime}}, m}\left(\boldsymbol{y}\right) \nonumber\\ & = \frac{\omega_0}{c}\left[ - T_{\ell, \ell^{^{\prime}}, m} \left(\boldsymbol{y}\right) \mathcal P_{\ell, m,\boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} +\tau_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \left( -\mathcal P_{\ell, m , \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} + \mathcal P_{\ell^{^{\prime}}, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell^{^{\prime}},m} \right) \right] \nonumber\\ & = \frac{\omega_0}{c} \left[ \left( \tau_{\ell^{^{\prime}},m} \left(\boldsymbol{y}\right) - \tau_{\ell, m}\left(\boldsymbol{y}\right) \right) \mathcal P_{\ell, m,\boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} +\tau_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \left( -\mathcal P_{\ell, m , \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} + \mathcal P_{\ell^{^{\prime}}, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell^{^{\prime}},m} \right) \right] \nonumber\\ & = \frac{\omega_0}{c}\left[ -\tau_{\ell, m}\left(\boldsymbol{y}\right) \mathcal P_{\ell, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell,m} + \tau_{\ell^{^{\prime}},m}\left(\boldsymbol{y}\right) \mathcal P_{\ell^{^{\prime}}, m, \boldsymbol{y}} \dot{\boldsymbol{\gamma}}_{\ell^{^{\prime}},m} \right] , \end{align} \tag{ D.3 }$$

Figure 4 is a sketch of the values of Ψ, as seen from the second equality of (39). First, the unit vectors $\widehat{\boldsymbol{R}}$ from the top diagram are copied to the bottom left diagram, and and their difference drawn. That gives the direction vector for the TDOA term of (39). Then that vector is multiplied by $\omega^{^{\prime}}/c$, where $\omega^{^{\prime}}$ varies over the bandwidth of the waveform. The roughly horizontal arrows at the bottom of the bottom right plot of figure 4 shows the difference of the unit vectors $\widehat{\boldsymbol{R}}$ multiplied by different lengths to represent the interval of $\omega^{^{\prime}}/c$ factors.

The FDOA term of (39) is more complicated. First, in (43), the scaled projections $\mathcal P$ are written in terms of an actual orthogonal projection Π divided by R. When Π operates on a vector, it projects that vector onto the plane perpendicular to $\widehat{\boldsymbol{R}}$. The vector that each Π is operating on is the corresponding sensor velocity vector $\dot{\boldsymbol{\gamma}}$. Those projections are shown in the top diagram of figure 4, and are copied into the bottom left diagram. Then in (43), those vectors are multiplied by a scalar $\omega_0 (\tau + t^{^{\prime}})/(cR)$ where t' ranges over the duration of the waveform. The bottom left plot of figure 4 shows those projection vectors multiplied by a scalar factor. Then (43) is the difference of those vectors for the same value of t'. Those difference vectors are represented as the downward pointing arrows in the bottom left diagram of figure 4. Those arrows all have the same direction but different lengths, ranging from a certain minimum to a certain maximum. Then, in (39), the negatives of those vectors are added to the first term. This is represented by the upward pointing arrows in the bottom right plot of figure 4: to each of the multiples of the difference of the unit vectors from the first term, we add an upward-pointing vector with a certain interval of lengths. That is represented by the blue shaded region, which is the DCM for this single look