Maximum-likelihood spatial spectrum estimation in dynamic environments with a short maneuverable array
Jonathan L. Odoma) and Jeffrey L. Krolikb) Dept. of Electrical and Computer Engineering, Duke Univ., PO Box 90291 Durham NC 27708 Jeffrey S. Rogersc) Acoustics Division, Naval Research Lab, Washington, DC 20375 (Dated: October 31, 2012)
J. Odom et al.: Spatial spectrum estimation with short arrays
1
Abstract This work concerns the development of field directionality mapping algorithms for short acoustic arrays on mobile maneuverable platforms that avoid the left/right ambiguities and endfire resolution degradation common to longer non-maneuverable line arrays. In this paper, it is shown that short maneuverable arrays can achieve a high fraction of usable bearing space for target detection in interference-dominated scenarios, despite their lower array gain against diffuse background noise. Two narrowband techniques are presented which use the ExpectationMaximization maximum likelihood (ML) algorithm under different models of the time-varying field directionality. The first, Derivative Based Maximum Likelihood (DBML), uses a deterministic model while the second, Recursive Bayes Maximum Likelihood (RBML), uses a stochastic model for the time-varying spatial spectrum. In addition, a broadband extension is introduced that incorporates temporal spectral knowledge to suppress ambiguities when the average sensor array spacing is greater than a half-wavelength. Dynamic multi-source simulations demonstrate the ability of a short, maneuvering array to reduce array ambiguities and spatial grating lobes in an interference dominated environment. Monte Carlo evaluation of receiver operating characteristics is used to evaluate the improvement in source detection achieved by the proposed methods versus conventional broadband beamforming.
PACS numbers: 43.60.Jn, 43.30.Wi, 43.60.Fg
2
I. INTRODUCTION
Passive sonar detection in shallow-water littoral areas is typically limited by the presence of surface shipping interference and/or transmission losses due to signal interaction with the bottom. In principle, both of these challenges can be addressed by performing field directionality mapping (or spatial spectrum estimation) using large acoustic arrays, either towed or static. For example, the well-known bearing-time-record (BTR) is essentially a time-varying spatial spectral estimate used for discrimination of sources of interest from interference, as well as signal gain against diffuse noise. However, a difficulty in littoral areas is that the achievable spatial gain is often insufficient to overcome high bottom-induced signal transmission losses. In such cases, detection of weak targets requires moving acoustic sensors physically closer to targets of interest. This has motivated large N system concepts wherein a dense network of omni-directional sensors is used in order to reduce transmission losses to potential targets.1 While, in theory, large N approaches, could provide satisfactory detection against diffuse noise, the concept fails to address the challenge of interference rejection in busy littoral areas. More recently, the use of underwater unmanned vehicles (UUVs) has been proposed.2 While the mobility of these UUV platforms offers the potential to position or distribute them to locations that minimize transmission losses, the need to keep their array apertures short seriously limits their ability to perform field-directionality mapping using conventional beamforming methods. This paper, therefore, presents field directionality mapping algorithms specifically designed for use with short acoustic arrays which provide discrimination of interferers while being on platforms that have the mobility to minimize transmission loss. Besides conventional beamforming approaches, a variety of approaches for field directionality mapping with long towed-arrays have been previously developed.3–8 These methods a)
Electronic address:
[email protected] Electronic address:
[email protected] c) Electronic address:
[email protected] b)
3
often assume stationarity of the acoustic field, which is appropriate for background noise mapping.3 Or alternatively, methods exploit the temporal coherence of the source to achieve a synthetic aperture with improved spatial resolution.9–12 Motion of the sources and platform limit the time window where the acoustic field is stationary.13 This motivates maximum likelihood methods using structured covariance models and/or low-rank approximations that can be estimated with less training data.14–16 More recently, an approach for time-varying spatial spectrum estimation with long maneuverable towed arrays using the EM algorithm was presented by Rogers and Krolik.17 In this paper, the work of Rogers and Krolik is extended for field directionality mapping with short maneuvering arrays. The additional challenge with short arrays is that their inherently poor spatial resolution limits the usable bearing space fraction (UBSF) available for weak target detection, i.e. performance is significantly degraded towards endfire. In this paper, we present methods to overcome this limitation by exploiting the greater mobility and maneuverability offered by the small platforms used to deploy such arrays. Although the use of short arrays limits array gain against diffuse noise, this may be offset by distributing several small maneuverable vehicles over the surveillance area which may offer a favorable cost-benefit trade-off in bottom-loss-limited littoral settings. Two models for the time-varying spatial spectrum are considered in this paper: 1) as unknown deterministic parameters, and 2) as a time-varying stochastic process. The deterministic model is inspired by radar literature on derivative-based updating (DBU) for space-time adaptive processing.18 Under this model, the spatial spectrum and its temporal gradient are treated as deterministic unknowns to be estimated. The stochastic model is a generalization of previous work19 by proposing an explicit time-varying random model of the spectrum. Unlike previous approaches, however, this paper incorporates constraints on the temporal spectrum of the source that are shown to enhance broadband spatial spectral estimation. The remainder of the paper is organized as follows. Section II provides the statistical model for data received by an array. Section III introduces two time-varying narrowband 4
spatial spectrum estimates. In Section IV, the narrowband performance of each estimate is compared using simulated data from an interference-dominated environment with a moving source. Simulations are presented with the assumption that a single snapshot is available for each of a series of array locations consistent with data obtained from a small maneuvering array. Broadband spectrum parameterization is considered in Section V, and Cram`er-Rao bounds are used to compare the traditional and proposed broadband approaches. This motivates a new broadband estimator with simulation results in Section VI. A simulation demonstrating the broadband estimate working in the presence of spatial grating lobes is provided to demonstrate the utility of incorporating temporal spectral knowledge.
II. RECEIVED DATA MODELS This section introduces narrowband and broadband temporal frequency domain models for data received by a maneuvering array. Assume an array of M sensors collects bandlimited discrete-time data for an observation window indexed by n, referred to as a time instance. The observation time multiplied by the received signal bandwidth is assumed to be large. During data collection, the field is assumed static, and the array platform and sources move between time instances. Consider a narrowband received signal model at angular frequency ω with Q sources for the nth time instance of the form, x(n, ω) =
Q X
dq (n, ω)sq (n, ω) + η(n, ω)
q=1
where
the
received
data
of
each
element
forms
the
vector
x(n, ω)
=
[ x1 (n, ω), x2 (n, ω), · · · , xM (n, ω) ]T . The source wavenumber is defined by
cos θ(n) sin φ(n) ω , k(n, ω) = sin θ(n) sin φ(n) c cos φ(n)
where c is the phase speed (m/s) of the medium. Bearing and elevation, θ and φ respectively, are relative to a constant direction (e.g. North) such that array motion does not rotate the 5
coordinate system, and the origin translates with the phase center of the array. The steering vector for the qth source and mth sensor element is given by [dq (n, ω)]m = exp (−jkT q (n, ω)rm (n)), where rm (n) is the position of the sensor relative to the phase center of the array. The received signal vector is expressed x(n, ω) = D(n, ω)s(n, ω) + η(n, ω) where D(n, ω) = [d1 (n, ω), d2 (n, ω), · · · , dQ (n, ω)].
(1)
All random variables are circular
symmetric complex normal distributed but are referred to as complex normal for convenience. The signal vector is complex normal distributed, s(n, ω) ∼ CN (0, Σ(n, ω)) and E[s(n1 , ω)s(n2 , ω)] = Σ(n1 , ω)δ(n1 − n2 ). The signals are independent of the noise, distributed η(n, ω) ∼ CN (0, Q). Note the snapshots are assumed to be uncorrelated across time samples, as is typical in passive sonar applications. The received signal is therefore complex normal x(n, ω) ∼ CN (0, R(n, ω)) with corresponding covariance matrix R(n, ω) = D(n, ω)Σ(n, ω)DH(n, ω) + Q.
(2)
Broadband signals are modeled by Fourier synthesis, using a statistical model first published by Bangs.21 The broadband source signal vector is formed by stacking B narrowband vectors. T ¯ s(n) = sT (n, ω1 ), sT (n, ω2 ), · · · , sT (n, ωB )
Assuming a large time-bandwidth product, the Fourier coefficients of each frequency are un¯ correlated such that s(n, ω) ∼ CN (0, Σ(n)) has block diagonal covariance matrix expressed as
0 Σ(n, ω1 ) . . ¯ .. Σ(n) = 0 Σ(n, ωB )
The frequency domain model is used in order to exploit prior knowledge of source temporal ¯ spectra, which constrains the structure of Σ(n), without full knowledge of the source signals. 6
Similarly, stack the received data vector and steering vectors to form: T x ¯(n) = xT (n, ω1 ), xT (n, ω2 ), · · · , xT (n, ωB )
T ¯ D(n) = DT (n, ω1 ), DT (n, ω2 ), · · · , DT (n, ωB )
which are assumed to be distributed according to a complex normal distribution, x ¯(n) ∼ ¯ ¯ ¯ Σ(n) ¯ D ¯ H (n) + ση2 I or CN (0, R(n)), with block covariance matrix R(n) = D(n) 0 R(n, ω1 ) . . ¯ .. R(n) = 0 R(n, ωB )
(3)
III. NARROWBAND SPATIAL SPECTRUM ESTIMATORS Spatial spectrum estimation is formulated here as a covariance matrix estimation problem, where the source signal covariance is the unknown parameter of interest. In this section, only narrowband data is considered. Assume the angular space around the array platform is a grid of Q points. Discretization provides the bearing angles for a bearing-time-record (BTR) or time-varying field directionality map. Ambient sea noise can be modeled as a sum of uncorrelated signals assuming the field is homogeneous.20 In this case, the source covariance matrix, Σ, is diagonal and represents power at each of the Q grid points. Time-varying field directionality and platform dynamics are physically constrained and captured by either a deterministic or stochastic model of the underlying covariance matrix, as described in Section III.A and III.B respectively. The deterministic method provides more parameters at higher computational cost. Note that unlike methods where tracking is performed post-detection, here the spatial spectrum dynamics are directly modeled in the pre-detection processing.?
A. Derivative Based Maximum Likelihood: DBML The first approach models the time-varying spatial spectrum deterministically through a set of known basis functions. Dependence on frequency is suppressed for notational con7
venience throughout narrowband derivations. The ML estimate of the field directionality at snapshot n is given in terms of the probability density function (pdf), f , of the received data for the previous N snapshots is written as ˆ n = arg max f (x(n), . . . , x(n − N + 1)|Σ(n)) . Σ
(4)
Σ(n)
In the radar literature,18 derivative-based updating (DBU) methods compensate for the time-varying change of returns across the array (and therefore the clutter covariance matrix) using a first-order model of the adaptive weight vector. In this paper, DBU is used to approximate slowly time-varying field directionality at time instance k near n using a linear model given by ˙ Σ(k) = Σo (n) + [k − (n − N + 1)]Σ(n),
(5)
which applies over a sliding window where n − N + 1 ≤ k ≤ n and Σo represents the spatial spectrum at the start of the window. Let k˜ = k − (n − N + 1) be a shifted version of k. By substituting (5) into (2) , the approximation for the received signal covariance matrix is, ˜ ˙ Rx (k) ≈ D(k) Σo (n) + (k − n)Σ(n) DH (k) + Q
(6)
˙ as unknown nonrandom parameters for k = [n, n − 1, . . . , n − N + 1]. Treating Σo and Σ leads to a maximum likelihood estimate given by ˆ˙ = arg max f x | Σ (n), Σ(n) ˆ n, Σ ˙ Σ n n o ˙ Σo (n),Σ(n)
where X(n) = [x(n), x(n − 1), · · · , x(n − N + 1)]. Since brute-force maximization of the likelihood function is computationally infeasible, the iterative Expectation-Maximization (EM) technique is used to approximate the solution.22 The EM algorithm requires splitting data into incomplete, or observable data, and complete, or unobservable data, where there is a many-to-one relation from the complete data to the incomplete data. The first step is to take the expectation of the likelihood containing the complete data given the incomplete data and a previous parameter estimate. The 8
second step is maximizing the result with respect to the unknown parameters. As described in previous array processing work using EM algorithms,14,17,23 define the complete data as the hidden source and noise signals, CN (k) = {s(k), η(k), · · · , s(k − N + 1), η(k − N + 1)} while the incomplete data is X(n). The likelihood function with the complete data is given ˙ = f (CN (n)|Σo , Σ)
n+N Y−1
˙ G(0, Σo (n) + (k˜ − n)Σ(n))G(0, ση2 I)
(7)
˜ k=n
in terms of the shift variable k˜ from (6) and G(µ, Σ) represents the circular symmetric complex normal density function. Assuming each source is uncorrelated, the expectation of (7) with respect to s, η given a previous iteration is h i old ˙ old ˙ ECN (n) ln f (CN (n)|Σo (n), Σ(n)) | Σo , Σ , X(n) ∝ Q n+N X−1 X ˜ k=n
˙ − ln([Σ(n)]qq + (k˜ − n)[Σ(n)] qq )
q=1
− h i H ˜ old ˙ old ˜ ˜ where hq (k) = E s(k)s (k) | Σo , Σ , X(n)
qq
˜ hq (k) (8) ˙ [Σ(n)]qq + (k˜ − n)[Σ(n)] qq
is the element (q, q) of the matrix given by
h i ˙ old, X(n) = Σold + (k˜ − n)Σ ˙ old − E s(k)sH (k) | Σold , Σ o o ˙ old ˜ − n)Σ ˙ old G(k) Σold + (k˜ − n)Σ Σold + ( k o o ˆ and the observed data appears in the estimate through the sample covariance matrix R(k) = X(n)X(n)H in G(k), which is defined in (9). The estimated covariance matrix using the previous estimate is defined as K(k), given in (10). To perform the M-step, the maximization ˙ of (8) with respect to Σ(n) and Σ(n) is preformed numerically using the interior point algorithm.24 The field directionality map, Σ(n), is constrained to be greater than 0 at every location. h i −1 ˆ G(k) = DH (k) K−1 (k) − K−1 (k)R(k)K (k) D(k)
(9)
˙ old H K(k) = D(k)(Σold o + k Σ )D (k) + Q.
(10)
9
B. Recursive Bayes Maximum Likelihood: RBML Using a Bayesian approach and stochastic model for the field directionality, assume the spectrum is a Markov random process. The Q × Q diagonal matrix Σ is assumed to be approximated by Σ(n + 1) = Σ(n) + ∆,
(11)
during the snapshot window from n − N + 1 to n where ∆ represents zero-mean state noise on the power from each angle snapshot window of n − N + 1 to n. A more general model can be used, but this simple model leads to a computationally tractable solution for estimating the time-varying frequency-domain spatial spectrum. Uncertainty in the spatial spectrum dynamics is described by the diagonal state noise ∆. This provides a probabilistic relation between sequential spectra, referred to as the transition pdf f (Σ(k) | Σ(k − 1)) .
(12)
Note that Σ(n) in (11) is only dependent on the previous time-step Σ(n − 1). This allows the joint pdf to be simplified into a product written in terms of (12) given by f (Σ(n), · · · , Σ(n − N + 1)) =
n Y
f (Σ(k)|Σ(k − 1)),
(13)
k=n−N +1
which will be lead to recursion in the estimate. The RBML estimate is defined in terms of the joint distribution of the data and the field directionality according to: ˆ n = arg max f (X(n), Σ(n)). Σ
(14)
Σ(n)
Letting all possible sequences Σ(k) for k = [n, · · · , n − N + 1] be defined as S, the joint distribution in (14) can be written as f (X(n), Σ(n)) =
Z
S
f X(n), Σ(n), · · · , Σ(n − N + 1)
dΣ(n − 1) · · · dΣ(n − N + 1), (15)
but is not computationally feasible to maximize for two reasons: the complex relation between X(n) and Σ(n) as well as the large size of S. 10
To solve (14) and (15), the EM algorithm is again used to simplify the integration over S, as suggested by Rogers and Krolik17 by invoking the forward method from Markov chains25 (pp. 112-113). The work presented in this section considers the general form of the transition pdf in (12) while previous work17 only considered f (Σ(k)|Σ(k − 1)) = δ(Σ(k) − Σ(k − 1)). In order to use the EM algorithm, the joint pdf from (15) is written in terms of the complete data CN and manipulated using Bayes’s theorem, f (a, b) = f (a|b)f (b), to express the joint pdf of the completed data and field directionality: f (CN (n), Σ(n − N + 1), · · · , Σ(n)) = f (CN (n)|Σ(n) · · · Σ(n − N + 1))f (Σ(n) · · · Σ(n − N + 1)). (16) To simplify the joint pdf consider the likelihood of the spatial spectrum given by: f (CN (n)|Σ(n), · · · , Σ(n − N + 1)) = n Y
f (s(k), η(k)|Σ(k)). (17)
k=n−N +1
Using (13) and (17), the joint pdf in (16) is simplified into a recursive relation f (CN (n),Σ(n), · · · , Σ(n − N + 1)) = =
n Y
f (s(k), η(k)|Σ(k)) · f (Σ(k)|Σ(k − 1)) (18)
k=n−N +1
=f (s(n), η(n)|Σ(n))f (CN −1 (n − 1), Σ(n − 1)) f (Σ(n)|Σ(n − 1)) with initial conditions f (C1 (n − N + 1), Σ(n − N + 1)) assuming a uniform distribution over Σ(n − N). The joint pdf from (15) is then written in terms of the complete data and the recursion relation, (18), in the form ( Z
f (CN −1 (k − 1), Σ(k − 1))
f (CN (k), Σ(k)) =
Sk−1
)
f (Σ(k)|Σ(k − 1))dΣ(k − 1) f (s(k), η(k)|Σ(k)) for k = n − N, · · · , n (19) 11
with integration over Sk−1 , all possible values of Σ(k − 1). Note that the form of the integration is similar to an expectation of f (CN (k − 1), Σ(k − 1). A second order approximation is derived in Appendix A. The natural logarithm of the joint pdf in (19) for k = n is now given by, ln f (CN (n), Σ(n)) =
n X
− ln det (πΣ(n))−
k=n−N +1
sH (k)Σ−1 (n)s(k) +
n−1 X
lk , (20)
k=n−N +1
where the second order terms, lk , are defined in (A3).
The expectation and maximization steps of the EM algorithm are now derived by taking the expectation of the joint pdf of the complete data and spectrum given the incomplete data, X(n), and the previous estimate Σold . Note the only terms in (20) dependent on X(n) are the signal terms, s. Let the source signal terms from the expectation step be defined hq (k) , E |sq (k)|2 | Σold , X(n) . Taking the expectation of (20) with respect to the complete data results in
EC ln f (CN (n), Σ(n)|Σold , X(n)) =
Q X q=1
(
− N ln[Σ(n)]qq
) n n−1 X X (N − 1)rq2 rq2 hq (k) hq (k) − , (21) − + 2[Σ(n)]qq [Σ(n)]qq k=n−N +1 2[Σ(n)]2qq k=n−N +1
where rq2 is the variance of [∆]qq in (11). The sum over Q sources is a result of assuming P incoherent sources. Let hq = nk=n−N +1 hq (k) denote the conditional expectation of the qth
source power summed over the observation window of N snapshots. The solution to the maximization step of the EM algorithm is found by setting the derivative of (21) to zero and solving for Σ(n). Consider the qth source and refer to the result from (21) as z where the derivative is expressed 0.5(N − dz −N = + d[Σ(n)]qq [Σ(n)]qq
1)rq2
+
n X
hq (k)
k=n−N +1 [Σ(n)]2qq
rq2 − [Σ(n)]3qq 12
n−1 X
k=n−N +1
hq (k). (22)
Since (22) only depends on q, the maximization of Q parameters in (21) reduces to Q single parameter optimizations. The cubic equation has 2 non-trivial solutions with only one solution assuming the spatial spectrum variance is small, given by new ˆ [Σ(n)] qq =
hq N −1 2 + r 2N 4N q +
1 2N
N −1 2 rq + hq 2
2
− 4Nrq2
n−1 X
k=n−N +1
hq (k)
! 21
(23)
The RBML estimate is given by (23) and at each numerical iteration the solution is set to be positive and greater than zero. This provides a generalization of the time-varying model assumed by Rogers and Krolik.17,19 The result from their previous work can be obtained by setting rq = 0 for q = [1, · · · , Q] in the solution (23). The previous solution simplifies to hq /N, which is the time-average of the expected values of the signal covariance, hq (k). Note the term hq (k) is the qth element of the vector given by E |s(k)|2 |Σold, x = Σold − Σold G(k)Σold
(24)
where G(k) is defined in (9) and estimated received data covariance estimate from the previous EM step is K(k) = D(k)Σold DH (k) + Q.
(25)
Lanterman considers Bayesian estimates using a combination of a likelihood and non-uniform prior (often improper and referred to as regularization)23 , which frequently results in numerical optimizations with high computational cost, as in the DBML method of Section III.A.
IV. NARROWBAND SIMULATION Numerical simulation is used to demonstrate the performance of the two spatial spectrum estimators developed in Section III. The goal is to demonstrate that capturing array dynamics allows a much smaller maneuvering uniform linear array (ULA) to provide field 13
directionality maps comparable to larger arrays by enhancing endfire directions and suppressing backlobes. Assuming white noise Q = ση2 I, the signal-to-noise (SNR) ratio is defined as [Σ]qq /ση2 for the qth source grid point, using received levels before array gain. Non-white noise is considered in broadband case. Consider a scenario with 4 interferers at positions relative to North at -90◦, -14◦, 135◦ each with 10 dB SNR and -135◦ with 3 dB SNR. The target begins at 60◦ and transitions to 15◦ with a maximum instantaneous bearing rate of -0.26◦ per second with SNR of 3 dB. The cartoon bearing-time-record (BTR) illustrates the scenario in Figure 1. BTRs provide a method of viewing the angular map of source directions as it changes over time. 300
3 dB
10 dB
10 dB 3 dB
10 dB Target
250
Interference
Time, sec
200
150
100
50
0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
FIG. 1. Illustrative BTR for all simulations For the narrowband case, each BTR is computed from -180◦ to 180◦ with 4◦ spacing and snapshot length of 1 sec. The estimates are computed every 2 seconds, resulting in an overlap of
N −2 N
for N > 2. The simulated data is generated according to (1) assuming
a narrowband frequency of 750 Hz and speed of sound of 1560 m/s. Conventional beamHˆ ˆ forming is accomplished by computing [Σ(k)] qq = [D(k) R(k)D(k)]qq at each time instance
independently using a rectangular window. The system used to show upper bound or “clairvoyant” performance in this example is a filled 660 element circular array with inter-element spacing of λ/2 and diameter ≈ 110λ. Circular arrays provide uniform angular resolution. This is compared against a short 7 element ULA maneuvering around the same circle at 14
2 m/s for 5 minutes with heading increments of +1◦ per second. The circle is not fully completed in the simulation, and the circular array is matched to the path covered by the short maneuverable array. Unlike expensive, large surveillance arrays, short maneuverable arrays can be dispersed throughout an environment of interest. Assuming a short array can be placed at 1/10 the distance to sources, the received SNR would increase by 20 dB under spherical spreading. An equivalent comparison used in this paper is to reduce the array gain of the 600-element array by 20 dB. This normalizes the array such that the large and small array achieve approximately the same output SNR. The large filled circular array is assumed to have a nearly perfect estimate of the received data covariance matrix for conventional beamforming, assuming 100 snapshots at each time instance. The resulting BTR, shown in Figure 2, displays each source clearly with a noise floor prescribed by the uniform weights. Previous work has compared a large array with several short fixed arrays at scattered locations using the same total number of elements using optimum placement.2
300
0 −2
250
−4 −6
Time, sec
200
−8 −10
150
−12 −14
100
−16 −18
50
−20 −22
0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
FIG. 2. BTR showing power estimates, in dB, from narrowband conventional beamformer output for clairvoyant circular array.
Conventional beamforming is optimal in spatially diffuse noise limited environments with performance determined by the number of sensors and array weighting schemes. When using few sensors with many sources, however, the interference often limits the detection of 15
weak targets. This occurs in the second array configuration, which utilizes 7 sensors in a ULA at λ/2 spacing while maneuvering over the shape of the circular array. The symmetric geometry of a ULA causes ambiguities that appear as backlobes on the BTR. Conventional beamforming at each time instance is unable to distinguish true sources from backlobes, which appear as diagonal stripes or extremely high bearing rate sources with circular platform maneuvers, shown in Figure 3(a). Poor endfire resolution of the ULA is visible as thick lines progressing from -180◦ to 180◦ and also from 0◦ wrapping around to 0◦ over time. These problems can be mitigated with time-varying spatial spectrum estimation methods described in Section III. The DBML and RBML algorithms perform 10 iterations every 2 sec, using only a single snapshot for each location. The algorithms approach convergence in this case, but the size of the sliding snapshot window, N, used by each algorithm is limited by the target motion. It is assumed that a single snapshot is available for each time instance, n so that the number of snapshots used by the estimate is same as the size of the sliding window, N. The DBML estimate, using the previous N = 9 snapshots with 7/9 overlap, displays the interferer and target tracks clearly without backlobes, as seen in Figure 3(b). For the first 50 seconds, the target is transitioning across the backlobe from the source at 135◦ and across the end-fire of the array creating uncertainty in the target power estimate. A similar ambiguity occurs again after the 200 second mark when the target passes through endfire again. The RBML estimate uses N = 1 with no overlap and rq = 0, as was used for a recent online technique originally developed for maneuvering long towed arrays.19 The maneuverability of the short array and DBML/RBML allow significant suppression of the array ambiguities and increased endfire resolution. Comparing Figure 2 with Figures 3(b) and 3(c), the results show that the usable bearing space of short mobile arrays can approximate that of much larger arrays. Differences in array gain can be addressed by increasing the density of short arrays over the surveillance area or by reducing transmission loss through location (e.g. depth), thus exploiting the maneuverability of the small arrays.
16
300
0
300
0
−5
−5
250
250 −10 200
−15 −20
150 −25 100
−30
Time, sec
Time, sec
200
−10 −15 −20 150 −25 100
−30
−35 50
−35 50
−40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
−45
−40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
(a)Conventional BTR
150
−45
(b)DBML BTR
300
0 −5
250 −10
Time, sec
200
−15 −20
150 −25 100
−30 −35
50 −40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
−45
(c)RBML BTR FIG. 3. BTR showing narrowband power estimates, in dB, with a maneuvering short array using (a) conventional beamforming, (b) DBML estimate, and (c) RBML estimate.
V. BROADBAND SPATIAL SPECTRUM ESTIMATE WITH PARAMETER REDUCTION While the previous section only considers narrowband sources, this section introduces a broadband estimate and analyzes the effect of temporal spectrum knowledge. A priori information of the temporal spectrum has the potential to increase detection and localization performance.26 In this paper, temporal spectrum knowledge is incorporated into the model ¯ and is used to reduce the number of estimated of the broadband signal covariance matrix, Σ, 17
parameters. In this section, the Cram´er-Rao bound (CRB) is derived to demonstrate the utility of exploiting a priori knowledge of the broadband source spectra. The bound on the covariance of an unbiased estimator is given by cov(ψ) ≥ J−1 or for the scalar case var(ˆ σ 2 ) ≥ J −1 . Traditionally, the parameter of interest is average (or total) power as a function of direction without a priori knowledge of the temporal spectra. Considering P the average power parameter defined by ψ = σ˙ = B1 B b=1 σ(ωb ) and σ(ωb ) = diag (Σ(ωb )),
the CRB is given by (26) assuming each frequency bin is independent. Note that the Q parameters of interest are a function of Q × B parameters. The bound is derived in Appendix B. J−1 ˙ c (σ)
B 1 X −1 J (σ, ωb ) = 2 B b=1
(26)
Now consider the case where it is assumed a priori that the temporal spectra are flat. The parameter of interest is defined ψ = σ where σ = σ(ωb ) ∀ b = [1, . . . , B]. With this parameterization, the number of estimated parameters is Q. The CRB is given in (27) and derived in Appendix B. J−1 B (σ)
=
" B X b=1
J(σ, ωb )
#−1
(27)
The bound in (26) will be referred to as the full ML approach since it estimates each parameter then averages across frequency. The bound from (27) will be referred to as the reduced ML approach since it performs parameter reduction before estimation. Comparing the full ML approach with the reduced ML case, it can be seen that there will only be a significant difference when J(ωb ) varies between frequencies. A good example of large variations in J(ω) is when spatial aliasing occurs and grating lobes create incorrect power estimates in ambiguous directions. Comparing the forms of (26) and (27), uncertainty in the spectral shape parameters will cause the grating lobe frequency terms in (26) to dominate the full estimate while the reduced estimate is less sensitive to broadband frequency dependent parameter uncertainty. Consider an example with a 10 element ULA with minimum 3λ/2 spacing. A 10 dB SNR source at 10◦ transmits in the 500-750 Hz range. The inter-element spacing, ≈ 3m, results in 2 ambiguous regions over 18
0 Reduced ML ≈ Full ML Std. Deviation/Mean in dB
−5
−10
−15
−20
−25
20
40
60
80 100 120 Bearing of Source, deg
140
160
180
FIG. 4. CRB of power estimate in aliasing region with 10 dB source at 10◦ as function of second source (0 dB) using 10 element ULA array with 3λ/2 spacing 180◦ which can be seen as jumps in the bound shown in Figure 4. The full parameter bound shown in Figure 4 uses an approximation discussed in Appendix B for numerical stability. The reduced parameter bound remains significantly lower than the full method, reducing the impact of spatial aliasing. This motivates incorporating parameter reduction in broadband extensions of the spatial spectrum estimates. Now the RBML estimate is extended for broadband sources using temporal spectral knowledge. In order to focus on the broadband effects, the derivations presented will consider the case where the variance of the spatial spectrum is small, rq = 0 for q = [1, · · · , Q] in (23), without loss of generality since any dynamic model can be incorporated, resulting in different likelihood functions and priors. The derivation for the RBML estimate in Section III.B is extended where the expectation in (21) is replaced with expectation step zB , h i ¯ Σ(n)| ¯ ˆ old , x ¯ n given by EC¯ ln f (C, Σ zB = −N
Q B X X
ln[Σ(n, ωb )]qq
b=1 q=1
−
n X
k=n−N +1
Q B X X hq (k, ωb ) [Σ(n, ωb )]qq q=1 b=1
!
(28)
where hq (k, ωb) , E |sq (k, ωb )|2 | Σold , X(n) . For simplicity, assume the spectrum is flat 19
such that Σ(n) = Σ(n, ωb ) for b = [1, · · · , B]. The estimate is now an optimization directly over Σ(n), which is the reduced broadband parameter to be estimated. The maximization step is obtained by taking the derivative of (28), shown in (29), and setting it equal to 0. n B X X hq (k, ωb ) ∂zB −NB = + . ∂[Σ(n)]qq [Σ(n)]qq k=n−N +1 b=1 ([Σ(n)]qq )2
(29)
This reduced broadband solution using the previous estimate and data across all frequencies has a form given by ˆ new ˆ old − Σ ˆ old 1 Σ =Σ n NB
n X
B X
ˆ old . (G(n, ωb ))Σ
(30)
k=n−N +1 b=1
where the frequency and time dependent terms are given by ˆ −1 )D G(n, ωb ) = DH (K−1 − K−1 RK ˆ old DH + Q. K(n, ωb ) = DΣ The reduced broadband solution can now be contrasted with the full broadband solution that uses no a priori temporal spectrum knowledge. The full broadband estimate averages narrowband estimates and is given by (31). Note that the frequency dependent variables are combined differently in the two estimates. The reduced estimate, (30), averages during an EM iteration while the full estimate, (31), averages after EM iterations converge. Incorporating temporal spectral information causes the reduced estimate to converge to a broadband joint solution. B X ˆn = 1 ˆ n (ωb ) Σ Σ B b=1
(31)
VI. BROADBAND SIMULATION WITH UNDER-SAMPLED ARRAY A simulation is used to demonstrate performance differences between a full ML approach, (31), and the reduced ML approach, (30) with under-sampled arrays subject to spatial aliasing. The simulation from Section IV is extended with true BTR given in Figure 1, now 20
over a 550-750 Hz bandwidth with 2 Hz bins and snapshot length of 1 sec. The array interelement spacing is 3λ/2 for the maximum frequency, ≈ 3m. For a uniform linear array, it is possible to have as many as 5 additional incorrect ambiguities from 2 grating lobes and a total of 3 backlobes. The clairvoyant circular array is reduced to 220 elements. The SNR is normalized for each array such that post-array gain is held constant as before. Conventional beamforming is used on the unrealisable circular array with incoherent frequency averaging resulting in relatively higher sidelobes, seen in Figure 5. The circular array does not suffer, however, from ambiguities due to spatial grating lobes.
300
0 −2
250
−4 −6
Time, sec
200
−8 −10
150
−12 −14
100
−16 −18
50
−20 −22
0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
FIG. 5. BTR showing power estimates, in dB, from broadband conventional beamformer output for under-sampled clairvoyant circular array.
A 7 element ULA array maneuvers over the path of the circular array as before. Note the length of the mobile platform is increased by a factor of 3 due to the larger spacing. Using the ULA, conventional beamforming output contains ambiguities as a result of spatial aliasing and array geometry, which vary as a function of platform orientation. This creates structured patterns appearing as maneuvering sources or swirls in the background of the BTR, shown in Figure 6(a). The full ML technique estimates the covariance matrix at each frequency, which contains spatial grating lobes, across all time and then averages the estimates across frequencies. This reduces the effects of spatial aliasing but results in 21
sporadic sources appearing in ambiguous regions, as shown in Figure 6(b). It is possible to introduce regularization, for example by setting rq in (23) to a non-zero value and forming a broadband solution, but this will also reduce the estimated power of high bearing rate sources. A trade-off exists in this case but is not explored in this paper. The reduced ML method provides significant suppression of spatial grating lobes, shown in Figure 6(c). The performance improvement from full ML to reduced ML is significant, with only a small ambiguity seen near 100◦ . There is also improvement from the narrowband case in Figure 3(c) to the broadband cases, although an increase in resolution is due to the larger length of the array in the broadband simulation. In order to understand the effects of surface wave noise on the field directionality map estimate, consider the case where the background noise covariance, Q, is generated from noise generated by the ocean surface uniformly with downward directivity given by cos2 (β), where β is the angle measured from a downward-facing axis perpendicular to the ocean surface. The resulting spatial correlation of the surface wave noise at the array can be found in closed form in (70) of Ref. 20 where m = 1. For example, consider the case where ˜ such that Q ˜ is the surface generated noise with diagonal equal to the sensor Q = ση2 I + Q noise. The resulting BTR has an increased noise floor due to the increase noise level as well as the correlation, as shown in Figure 7, but the performance does not otherwise change. This, the RBML algorithm is robust to such type of noise. Additionally, consider robustness to temporal spectrum mismatch. Consider the case where each source only transmits power over the middle 20% of the assumed 550-750 Hz band. Mismatch between the assumed source spectra and the received data, as well as the decrease in total signal energy, results in decreased performance. The algorithm is unable reduce spatial grating lobes to the levels achieved without mismatch by comparing the BTRs shown in Figure 6(c) and Figure 8. However, each of the five sources are clearly visible even though 80% of the assumed band is missing. A detection algorithm is used to provide an additional comparison among the different methods in the presence of spatial aliasing. This is also used to demonstrate the potential 22
300
0
300
0
−5
−5
250
250 −10 200
−15 −20
150 −25 100
−30
Time, sec
Time, sec
200
−10 −15 −20 150 −25 100
−30
−35 50
−35 50
−40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
−45
−40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
(a)Conventional BTR
150
−45
(b)Full ML BTR
300
0 −5
250 −10
Time, sec
200
−15 −20
150 −25 100
−30 −35
50 −40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
−45
(c)Reduced ML BTR FIG. 6. BTRs showing broadband power estimates, in dB, with an under-sampled maneuvering short array using (a) conventional beamforming, (b) full ML, and (c) reduced ML.
for low SNR target detection using an array mounted on a small UUV. The null hypothesis (H0 ) is the target-free interference dominated environment, and the alternate hypothesis (H1 ) is a target with given target track in same environment. Since ambiguities result ˜ detector is used that is robust to in false targets that cross the target track, an M-of-N ˜ = 100 observations of 1 second intervals, a detection occurs when spurious peaks.27 Using N M = 80 or more power (BTR) estimates exceed a given threshold along an assume target 23
300
0 −5
250 −10
Time, sec
200
−15 −20
150 −25 100
−30 −35
50 −40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
−45
FIG. 7. Broadband BTR demonstrating the performance of the reduced ML technique in the presence of surface wave noise.
300
0 −5
250 −10
Time, sec
200
−15 −20
150 −25 100
−30 −35
50 −40 0
−150
−100 −50 0 50 100 Bearing Relative to North, degrees
150
−45
FIG. 8. Broadband BTR using the reduced ML method with mismatched source signals. True source bandwidth 20% of assumed band.
track. The true target track is used and the threshold is varied to determine the receiver operating characteristics (ROC) in terms of probability of detection, Pd , and probability of false alarm, Pf . In order to consider low SNR target detection, A target SNR of -25 dB is used, and interference levels are unchanged. Using 100 realizations, the result shown in Figure 9. The conventional delay-and-sum beamforming is unable to distinguish the target 24
from the interferers. The full ML method provides some improvement over conventional, but the reduced ML method provides the best performance.
1
0.8
P
d
0.6
0.4
Reduced ML Full ML Conventional
0.2
0
0.2
0.4
Pf
0.6
0.8
1
˜ detector for under-sampled array FIG. 9. ROC curve with target SNR -25 dB using M-of-N
VII. CONCLUSIONS This paper introduces several narrowband and broadband time-varying spatial spectrum estimates that exploit array maneuvers in order to improve endfire resolution and suppress ambiguities. A traditional broadband extension is compared with a reduced broadband parameter estimation technique that incorporates knowledge of temporal spectral shape. CRB analysis yields a comparison of broadband techniques and the effects of temporal spectrum knowledge. An important result is the ability to mitigate the effects of spatial aliasing by exploiting prior temporal spectral knowledge. Broadband simulations demonstrate the ability of maneuvering short arrays to mitigate the effects of spatial aliasing. This leaves future work to define the mapping limits of under-sampled arrays. Optimal array maneuvers or multiple arrays are not considered here but would provide insight into the ability of distributed short arrays to provide similar mapping performance of much larger arrays. 25
APPENDIX A: SECOND ORDER EXPECTATION APPROXIMATION Consider a function g(x) with parameter x. The first three terms in the Taylors series expansion around µ are g(x) = g(µ) + g ′(µ)(x − µ) + 21 g ′′ (µ)(x − µ)2 . Assuming x has a mean µ and variance σ 2 , then E[g(x)] ≈ g(µ) +
σ2 ′′ g (µ). 2
This is essentially a local approximation
of g(x) assuming the pdf is tight around x near µ. The multivariate equivalent for x with mean µ and covariance R is E[g(x)] ≈ g(µ) + 12 Tr[RHµ g(x)], where Tr is trace and Hx is the Hessian with respect to x evaluated at µ. The approximation applied to the term in brackets from (19) is analogous to assuming the field does not vary quickly between time instances. The joint pdf from (19) is written as
f (CN (n), Σ(n)) = α + f (s(n), η(n)|Σ(n)) n−1 Y
1 G(s(k); 0, Σ(n)) + Tr [RHΣn G(s(k); 0, Σ(n))] , (A1) 2 k=n−N +1 where α is not a function of Σ(n) or X(n). The covariance matrix of ∆ is represented by R. Let the diagonal of the covariance be denoted [R]qq = rq2 for the qth element. Taking the natural log of (A1) results in
ln f (CN (n), Σ(n)) ∝ n X
ln G(s(k); 0, Σ(n)) +
k=n−N +1
n−1 X
lk (A2)
k=n−N +1
where
1 lk = ln 1 + Tr S Σ−1 (n)s(k)sH (k)Σ−1 (n) − Σ−1 (n) 2 Q rq2 1 X rq2 |sq (k)|2 − . ≈ 2 q=1 ([Σ]qq )2 (n) [Σ]qq
(A3) (A4)
When the spectrum changes slowly, the standard deviation rq and the resulting lk are small and allows for the approximation of ln given by (A4). 26
´ APPENDIX B: CRAMER-RAO BOUNDS The Cram´er-Rao bound for source parameters, ψ, using the model from Section V is dependent on the general broadband Fisher information matrix (FIM) given by (B1). [J(ψ]pq
¯ ¯ ∂ R(ω) ∂ R(ω) −1 −1 ¯ ¯ = Tr R (ω) R (ω) ∂ψp ∂ψq
(B1)
1. Bound without a priori knowledge First consider the parameter vector without a priori temporal spectrum knowledge, 2 2 σ ¯ = σ12 (ω1 ), · · · , σQ (ω1 ), · · · , σ12 (ωB ), · · · , σQ (ωB ) . The order of parameters is sources
followed by frequencies. Using this parameter vector, the broadband FIM, (B1), can be written in terms of the narrowband FIM, with form given by
0 J(ω1 ) .. . J= . 0 J(ωB ) The parameter of interest in this case is average power defined by σ˙ =
1 B
PB
b=1
σ(ωb ). Thus
the full ML method estimates each frequency and then averages the estimates. In this case, ¯ σ] ψ = [σ, ˙ so that the FIM contains an additional parameter and is denoted Jc (¯ σ , σ), ˙ given by
J Jσ¯ σ˙ Jc (¯ σ , σ) ˙ = . Jσ˙ σ¯ Jσ˙ The inverse can be written in the form [Jc (¯ σ , σ] ˙ −1 = HH J−1 H where [H]ab = σ ˙
∂σb2 28 . ∂σa2
(p.
230) The effect of H in matrix form is to perform the addition on each estimate. Thus the bound can be simplified into form given by (B2). J−1 ˙ = c (σ)
B 1 X −1 J (σ, ωb ). B 2 b=1
27
(B2)
2. Bound with a priori knowledge Alternatively, the reduced ML method assumes a priori knowledge such that a single 2 parameter describes each source and results in a parameter vector σ = σ12 , · · · , σQ . Since
data from different frequencies is independent, the FIM from (B1) can be written in terms of σ and given by (B3). Note the broadband FIM, JB (σ), can be written in terms of the narrowband FIM, as given by (B4). The broadband bound contains only Q parameters and has final form given by (B5). [JB (σ)]pq =
B X b=1
∂R(ωb ) −1 ∂R(ωb ) Tr R (ωb ) R (ωb ) . ∂σp2 ∂σq2 −1
B X
JB (σ) =
J(σ, ωb)
(B3)
(B4)
b=1
J−1 B (σ) =
" B X
J(σ, ωb )
b=1
#−1
(B5)
3. Remark on numerical computation of CRB ˜ Without loss For numerical stability of (B2), it is useful to consider an arbitrary FIM J. of generality on the total number of parameters, the partitioned case for a single variable is given by
˜= J
J11 JT 21 J21 J22
˜ are non-negative and J22 is positive semidefinite. Note that the where all elements of J bound on the parameter is −1 −1 J11 = J11 − JT 21 J22 J21
−1 J11 ≥ J11
˜−1 ]qq ≥ {[J] ˜ qq }−1 [J
(B6)
where the new bound is looser than the CRB but provides an approximation at least as low as the CRB. 28
ACKNOWLEDGMENTS This work was supported by ONR under grant N000140810947.
REFERENCES 1
T. Wettergren, “Performance of search via track-before-detect for distributed sensor networks”, IEEE Trans. Aerosp. Electron. Syst. 44, 314–325 (2008).
2
I. Bilik and J. Krolik, “Comparative theoretical performance of maneuverable unmanned vehicles versus conventional towed-arrays for passive sonar”, in Conf. Proc. IEEE OCEANS Europe, 1–5 (2007).
3
R. Wagstaff, “Iterative technique for ambient-noise horizontal-directionality estimation from towed line-array data”, J. Acoust. Soc. Am. 63, 863–869 (1978).
4
W. Allensworth, C. Kennedy, B. Newhall, and I. Schurman, “Twinline array development and performance in a shallow-water littoral environment”, Johns Hopkins APL Tech. Dig. 16, 222–232 (1995).
5
J. Sheinvald, M. Wax, and A. Meiss, “Localization of multiple sources with moving arrays”, IEEE Trans. Signal Process. 46, 2736 –2743 (1998).
6
M. V. Greening and J. E. Perkins, “Adaptive beamforming for nonstationary arrays”, J. Acoust. Soc. Am. 112, 2872–2881 (2002).
7
P. Gerstoft, W. Hodgkiss, W. Kuperman, H. Song, M. Siderius, and P. Nielsen, “Adaptive beamforming of a towed array during a turn”, IEEE J. Ocean. Eng. 28, 44–54 (2003).
8
H. Tao, G. Hickman, and J. Krolik, “Field directionality synthesis from multiple array orientations: A least squares approach”, in Proc. IEEE OCEANS, 1–6 (2008).
9
N.-C. Yen and W. Carey, “Application of synthetic-aperture processing to towed-array data”, J. Acoust. Soc. Am. 86, 754–765 (1989).
10
S. Stergiopoulos and H. Urban, “A new passive synthetic aperture technique for towed arrays”, IEEE J. Ocean. Eng. 17, 16–25 (1992).
11
Y. Lee, W. Lee, and A. Lee, “Synthetic stationary array processing for a towed array in 29
maneuver”, J. Acoust. Soc. Am. 120, 3049–3049 (2006). 12
E. J. Sullivan, J. D. Holmes, W. M. Carey, and J. F. Lynch, “Broadband passive synthetic aperture: Experimental results”, J. Acoust. Soc. Am. 120, EL49–EL54 (2006).
13
A. Baggeroer and H. Cox, “Passive sonar limits upon nulling multiple moving ships with large aperture arrays”, in Conf. Signals, Systems, and Computers (ASILOMAR), volume 1, 103–108 (1999).
14
M. Miller and D. Fuhrmann, “Maximum-likelihood narrow-band direction finding and the em algorithm”, IEEE Trans. Acoustics, Speech and Signal Process. 38, 1560–1577 (1990).
15
D. Abraham and N. Owsley, “Beamforming with dominant mode rejection”, in Conf. Proc. IEEE OCEANS, 470–475 (1990).
16
A. Kraay and A. Baggeroer, “A physically constrained maximum-likelihood method for snapshot-deficient adaptive array processing”, IEEE Trans. Signal Process. 55, 4048–4063 (2007).
17
J. S. Rogers and J. L. Krolik, “Time-varying spatial spectrum estimation with a maneuverable towed array”, J. Acoust. Soc. Am. 128, 3543–3553 (2010).
18
S. Hayward, “Adaptive beamforming for rapidly moving arrays”, in Proc. CIE Int. Conf. Radar, 480–483 (1996).
19
J. Rogers and J. Krolik, “An online method for time-varying spatial spectrum estimation using a towed acoustic array”, in Conf. Signals, Systems, and Computers (ASILOMAR), 1837–1841 (2010).
20
H. Cox, “Spatial correlation in arbitrary noise fields with application to ambient sea noise”, J. Acoust. Soc. Am. 54, 1289–1301 (1973).
21
W. Bangs, “Array processing with generalized beam-formers”, Ph.D. thesis, Yale Univ. (1971).
22
A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm”, J. Royal Statistical Society. Series B (Methodological) 39, 1–38 (1977). 30
23
A. Lanterman, “Statistical radar imaging of diffuse and specular targets using an expectation-maximization algorithm”, in Proceedings-SPIE of the International Society for Optical Engineering, 20–31 (2000).
24
R. A. Waltz, J. L. Morales, N. J., and O. D., “An interior algorithm for nonlinear optimization that combines line search and trust region steps”, Mathematical Programming 107, 391–408 (2006).
25
C. W. Therrien, Discrete Random Signals and Statistical Signal Processing (PrenticeHall, Upper Saddle River, NJ), 112–113 (1992).
26
H. Messer, “The potential performance gain in using spectral information in passive detection/localization of wideband sources”, IEEE Trans. Signal Process. 43, 2964–2974 (1995).
27
B. D. Carlson, E. D. Evans, and S. L. Wilson, “Search radar detection and track with the hough transform. iii. detection performance with binary integration”, IEEE Trans. Aerosp. Electron. Syst. 30, 116–125 (1994).
28
L. L. Scharf, Statistical Signal Processing (Addison-Wesley, Reading, MA), 230 (1991).
31
LIST OF FIGURES FIG. 1
Illustrative BTR for all simulations . . . . . . . . . . . . . . . . . . . . . . .
FIG. 2
BTR showing power estimates, in dB, from narrowband conventional beamformer output for clairvoyant circular array. . . . . . . . . . . . . . . . . . .
FIG. 3
14
15
BTR showing narrowband power estimates, in dB, with a maneuvering short array using (a) conventional beamforming, (b) DBML estimate, and (c) RBML estimate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
FIG. 4
CRB of power estimate in aliasing region with 10 dB source at 10◦ as function of second source (0 dB) using 10 element ULA array with 3λ/2 spacing . . .
FIG. 5
19
BTR showing power estimates, in dB, from broadband conventional beamformer output for under-sampled clairvoyant circular array. . . . . . . . . . .
FIG. 6
17
21
BTRs showing broadband power estimates, in dB, with an under-sampled maneuvering short array using (a) conventional beamforming, (b) full ML, and (c) reduced ML. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
FIG. 7
Broadband BTR demonstrating the performance of the reduced ML technique in the presence of surface wave noise. . . . . . . . . . . . . . . . . . .
FIG. 8
24
Broadband BTR using the reduced ML method with mismatched source signals. True source bandwidth 20% of assumed band. . . . . . . . . . . . . . .
FIG. 9
23
24
˜ detector for under-sampled ROC curve with target SNR -25 dB using M-of-N array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
25
