1

Matched-Field Processing Performance under the Stochastic and Deterministic Signal Models Yann Le Gall, Student Member, IEEE, Francois-Xavier Socheleau, Member, IEEE and Julien Bonnel, Member, IEEE Abstract Matched-field processing (MFP) is commonly used in underwater acoustics to estimate source position and/or oceanic environmental parameters. Performance prediction of the multisnapshot and multifrequency MFP problem is of critical importance. To this end, two signal models are usually considered: the stochastic model which assumes that the source signal is a stochastic process, and the deterministic model which assumes that the source signal is a deterministic quantity. The Ziv-Zakai bound (ZZB) and the method of interval errors (MIE), that both rely on the computation of a so-called pairwise error probability, proved to be useful tools for MFP performance prediction. However, only the stochastic model has been considered so far. This paper provides a method that allows to compute the pairwise error probability, hence to use the ZZB and MIE, under both the stochastic and deterministic signal models. The proposed approach, based on recent results on quadratic forms in Gaussian variables, unifies the two models under the same formalism. The results are illustrated through the computation of the ZZB and MIE performance analysis. The Bayesian and the hybrid Cram`er-Rao bounds are also given for comparison. Index Terms Underwater acoustics, matched-field processing, maximum likelihood, performance analysis, ZivZakai bound, method of interval errors, Cram`er-Rao bound

Copyright (c) 2014 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to [email protected]. Y. Le Gall and J. Bonnel are with Lab-STICC, ENSTA Bretagne, UMR CNRS 6285, 2 rue Franc¸ois Verny, 29806 Brest cedex 9, France (e-mail: yann.le gall, [email protected]). F.-X. Socheleau was with ENSTA Bretagne, he is now with Institut Mines-Telecom/Telecom Bretagne, UMR CNRS 6285 Lab-STICC, Technople Brest-Iroise, 29238 Brest, France (e-mail: [email protected]).

September 29, 2014

DRAFT

2

I. I NTRODUCTION Matched-field processing (MFP) has received significant attention in underwater acoustics for source localization and ocean environmental parameters estimation. MFP compares the data measured on an array of spatially distributed sensors with replicas of the acoustic field derived from the wave equation in order to infer the desired parameter values [1], [2]. A cost function termed processor is used to determine the values that produce the best match among the scanned parameters. Several MFP processors have been developed on an ad hoc basis or in the framework of maximum likelihood (ML) to address various assumptions on the source signal [3]–[5]. The source signal is either modeled as a random Gaussian process in the so-called stochastic signal model or as a deterministic quantity in the so-called deterministic signal model. Under the stochastic model, the source spectral density is usually assumed known [5] while under the deterministic model, the source signal is most of the time considered as partially or totally unknown [3], [4]. This paper aims at providing a performance analysis tool for MFP under both the stochastic and deterministic signal models. Prediction of MFP performance is a challenging problem. In previous works, MFP performance analysis has been carried out through limiting bounds or method of interval errors (MIE) analysis. Limiting bounds provide absolute bounds on the performance of any estimator for the observation model considered. In underwater acoustics, the Cram`er-Rao bound (CRB) has been used to provide necessary conditions on the signal-to-noise ratio (SNR) and/or the number of snapshots so that the estimation error stays small [6], [7]. It has also been used to predict MFP performance [8]–[10]. However, it is well known that the CRB is tight only at high SNR or for a high number of snapshots, namely in the asymptotic region (see Figure 3), when mostly local errors are made. Below a certain level of SNR or a certain sample size called threshold (see Figure 3), the CRB fails to capture the outliers encountered by the estimator. Fortunately, there exist more relevant bounds than the CRB that account for these global errors. Some of them have been applied to the MFP problem. The Barankin bound [11] was derived by Tabrikian and Krolik under both the stochastic and deterministic signal models in the special case where only one frequency is used [12]. The Weiss-Weinstein bound (WWB) [13] and the Ziv-Zakai bound (ZZB) [14], [15] were derived by Xu et al. under the stochastic signal model [16], [17]. These two later bounds are known to be among the tightest existing bounds and have shown to give tight and reliable error prediction for the MFP problem. Other kinds of bounds based on information theory have also been proposed recently by Meng and Buck [18]. The MIE is an alternative to these lower bounds. The idea behind MIE is to separate the errors into local errors and global errors in order to obtain an approximation of the mean square error (MSE) for a

September 29, 2014

DRAFT

3

given estimator. Thus, MIE does not give a lower bound for any estimator but an approximation of the MSE that captures local errors as well as global errors for a particular estimator. This technique has also proven to be very useful for numerous problems [19], [20] and was recently applied successfully to the MFP maximum likelihood (ML) estimator for the stochastic signal model [21]. For ML MIE analysis and for the evaluation of the ZZB, the main difficulty lies in the computation of the ML pairwise (or two points) error probability. The pairwise error probability is the probability Pe (θ 1 |θ 0 ) of deciding in favor of the parameter θ 1 given the true parameter θ 0 , in the binary hypothesis

test {H1 : θ 1 , H0 : θ 0 } with equally likely hypotheses. In MIE analysis, this probability is used to approximate the probability of outliers associated with each sidelobe peak of the ambiguity function. In the ZZB, the ML pairwise error probability is used to lower bound a suboptimal decision rule, which ultimately provides the bound on the MSE. This paper provides a method that allows to compute the exact MFP ML pairwise error probabilities under the stochastic and the deterministic signal models. The originality of this work primarily lies in the consideration of the deterministic signal model since the stochastic model has already been settled in a different approach [5], [16], [21]. The deterministic model on the other hand has, to our knowledge, never been considered in the full multisnaphot and multifrequency MFP problem, although a solution exists for the simpler monofrequency beamforming problem [19]. Furthermore, the proposed approach unifies the two signal models under the same formalism. In this paper, the results will be illustrated through the computation of the ZZB and MIE performance analysis. The Bayesian and the hybrid CRB are also given for comparison. The paper is organized as follows. The MFP problem, the data models and their corresponding ML estimators are presented in Section II. The MIE and the ZZB are introduced in Section III and IV respectively. The pairwise error probabilities for the considered data models are derived in Section V, and finally Section VI illustrates the theoretical results. Notation: Throughout this paper, lowercase boldface letters denote vectors, e.g., x, and uppercase boldface letters denote matrices, e.g., A. The superscripts

T

and

H

mean transposition and Hermitian

transposition, respectively. Re[·] denotes the real part and Im[·] the imaginary part. We let diag(x) designate a diagonal square matrix whose main diagonal contains the elements of vector x. The N × N identity matrix is denoted by IN and the N × M zero matrix is denoted by 0N,M . The operator |·| means determinant for matrices and modulus for scalars. The operator k·k designates the standard Euclidean norm. The distribution of a jointly proper Gaussian random vector with mean m and covariance matrix R is denoted CN (m, R). Finally, E[·] denotes expectation, Pr (·) denotes probability and δ is the Dirac

delta function. September 29, 2014

DRAFT

4

II. M ATCHED -F IELD

PROBLEM FORMULATION AND DATA MODEL

Consider an array of N sensors (generally a vertical line array) in the oceanic waveguide on which impinge the waveforms from an acoustic source farther in the waveguide. The goal of MFP is to estimate an unknown parameter set θ from the array measurements. This parameter set may be the source position (range, depth and sometimes bearing) for source localization, and/or ocean environmental parameters (bathymetry, sound speed, density...) for tomography and geoacoustic inversion. A typical MFP configuration is presented in Figure 1. This example will be used later to illustrate the theoretical results. For each snapshot l and each frequency fm , the complex array output is modeled by the N × 1 vector yl (fm ) = sl (fm ) · g(fm , θ) + wl (fm ), l = 1, ..., L, m = 1, ..., M

(1)

where, •

L is the number of snapshots l available for each of the M frequencies fm .

•

g(fm , θ) is a complex N × 1 vector representing the transfer function of the medium at frequency fm for the propagation from the source to each of the N receivers. It is usually termed Green’s

function. •

sl (fm ) is a complex scalar representing the source amplitude and phase.

•

wl (fm ) is a complex N × 1 vector representing a circularly-symmetric, zero mean Gaussian noise

that is independent of the source signal. It is assumed to be frequentially and temporally white with 2 (f )C (f ), i.e., positive definite spatial covariance matrix Σw (fm ) = σw m w m

  2 E wl1 (fm1 )wlH2 (fm2 ) = σw (fm1 )Cw (fm1 )δl1 ,l2 δfm1 ,fm2 ,   E wl1 (fm1 )wlT2 (fm2 ) = 0, ∀ l1 , l2 , fm1 , fm2 . (2)

In this paper we assume known covariance matrices. Since the matrices Cw (fm ) are positive definite, their square root decomposition are invertible and can be obtained via the Cholesky decomposition Cw (fm ) = Lw (fm )LH w (fm ). Therefore, the following whitening transform ˜ l (fm ) = L−1 ˜ (fm , θ) = L−1 y w (fm )yl (fm ), g w (fm )g(fm , θ), l = 1, ..., L, m = 1, ..., M

(3)

can be applied, which allows to reduce to the white noise case. In the remainder of this paper we will 2 (f )I to simplify subsequent derivations work in whitened coordinates and assume that Σw (fm ) = σw m N

and analyses. The set of observations is collected in the following vector T T y = [y1T (f1 ) . . . yL (f1 ) . . . y1T (fM ) . . . yL (fM )]T .

September 29, 2014

(4)

DRAFT

5

Two different models for the source signal sl (fm ) are usually considered, namely, the stochastic and the deterministic signal models which are sometimes also referred to as the unconditional and the conditional models: •

(A1) Stochastic signal model: In the stochastic model, the source signal is assumed to be a stationary Gaussian random process with known power spectral density σs2 (f ). Specifically, it is assumed to be temporally and frequentially white, zero mean and circularly-symmetric, i.e.,     T 2 E sl1 (fm1 )sH l2 (fm2 ) = σs (fm1 )δl1 ,l2 δfm1 ,fm2 , E sl1 (fm1 )sl2 (fm2 ) = 0, ∀ l1 , l2 , fm1 , fm2 . (5)

This model is justified by the large number of natural processes that can be modeled as stationary Gaussian random processes. The wideband noise generated by ships and submarines typically fall into this category [22]–[24]. Furthermore, randomness and variability of the ocean medium can be so important that the received signal is fully randomized and follows a Gaussian distribution even when the initial source signal is deterministic. This phenomenon is known as full saturation [25]. •

(A2) Deterministic signal model: In the deterministic signal model, the quantities sl (fm ), l = 1...L, m = 1...M are assumed to be deterministic. Depending upon the context, they are either

considered as unknown or partially known. This model is justified because there exist many source signals that exhibit a deterministic behavior or at least that are not Gaussian. This signal model is representative for instance of tonal components that may arise from submarines, ships [22], [26], [27] or deployed tomographic sources [28] and it may also be suitable for marine mammals vocalizations [29]. These two models are conventionally considered in MFP performance analysis [3]–[8], [10], [12], [16], [17], [21], but it should be noted that the deterministic model is very popular when it comes to ML MFP processor design [3], [4] and parameter estimation in the Bayesian framework [30]–[32]. A. Stochastic Maximum Likelihood estimate Consider the stochastic signal model (A1) defined previously. Given the source/environmental parameter set θ , the conditional probability density function (pdf) of the observation y is p(y|θ) = QM

m=1

1

QL

L M Y Y

l=1 |πΣy (fm , θ)| m=1 l=1

exp(−ylH (fm )Σy −1 (fm , θ)yl (fm )),

(6)

where Σy (fm , θ) is the covariance matrix of yl (fm ) given by 2 Σy (fm , θ) = σs2 (fm )g(fm , θ)gH (fm , θ) + σw (fm )IN .

September 29, 2014

(7)

DRAFT

6

The stochastic Maximum Likelihood (SML) estimate of the desired parameter set θ is obtained by 2 (f ) are assumed known, the SML estimate is maximizing (6) with respect to θ . When σs2 (fm ) and σw m

[5], [21] θˆSM L = arg max θ

"

M X

C1 (fm , θ) +

M X L X

#

C2 (fm , θ) · |¯ gH (fm , θ)yl (fm )|2 ,

m=1 l=1

m=1

where g ¯(fm , θ) is a normalized version of g(fm , θ) defined as g ¯(fm , θ) =

g(fm ,θ) kg(fm ,θ)k ,

and

2 C1 (fm , θ) = −L ln(σs2 (fm )kg(fm , θ)k2 + σw (fm )),

C2 (fm , θ) =

(8)

σs2 (fm )kg(fm , θ)k2 . 2 (f )(σ 2 (f )kg(f , θ)k2 + σ 2 (f )) σw m m s m w m

(9) (10)

Assuming a weak dependence of the norm of the Green’s function, this expression is often approximated by the more simple expression [5], [16], [17] θˆSM L = arg max θ

M X L X

|¯ gH (fm , θ)yl (fm )|2 .

(11)

m=1 l=1

Nevertheless, only the exact expression given by equation (8) will be considered in this paper.

B. Deterministic Maximum Likelihood estimate Consider the deterministic signal model (A2). Given the source/environmental parameter set θ , the conditional probability density function (pdf) of the observation y is   M Y L Y 1 kyl (fm ) − sl (fm )g(fm , θ)k2 p(y|θ) = QM QL exp − . 2 (f ) 2 N σw m m=1 l=1 |πσw (fm )| m=1 l=1

(12)

The deterministic Maximum Likelihood (DML) estimate of the desired parameter set θ is obtained by maximizing (12) with respect to θ , considering the source signal sl (fm ) as a nuisance parameter that is replaced by its ML estimate according to the profile likelihood technique [33]. Various amount of information on the sequence sl (fm ) can be considered, leading to different DML processors. 1) Incoherent DML: In the case where no information is available on the sequence sl (fm ), when 2 (f ) are assumed known, the DML estimate of the desired parameter set θ is [4] σw m

θˆDM L = arg max θ

M X L X |¯ gH (fm , θ)yl (fm )|2

m=1 l=1

2 (f ) σw m

.

(13)

This MFP processor is referred to as incoherent because the multiple snapshots and frequencies are summed incoherently, i.e., with the absolute value that cancels phase information inside the double sum.

September 29, 2014

DRAFT

7

2) Coherent DML: In the case where the relative amplitudes and phases between snapshots and frequencies are available (but the absolute values are unknown), the signal model reduces to sl (fm ) = s and the DML estimate of the desired parameter set θ is [4] 2 M X L X H (f , θ)y (f ) g 1 m m l θˆDM L = arg max PM kg(f ,θ)k2 . 2 m σw (fm ) θ L 2 m=1

σw (fm )

(14)

m=1 l=1

This MFP processor is referred to as coherent because the multiple snapshots and frequencies are summed coherently, i.e., with the absolute value outside the double sum, which allows to take into account phase information. It is mostly suitable for applications with controlled sources like active tomography. These two ML processors are end-members in terms of source information. Other DML processors can be derived for intermediate state of information [4] but only the above DML estimators are in the scope of this paper. It should be noted however that the performance of the other processors are likely to fall between that of the incoherent and coherent processors. III. M ETHOD

OF INTERVAL ERRORS

The MIE provides an approximation of the mean square error (MSE) for a given estimator (see the tutorial treatment by Van Trees and Bell in [34]). The idea behind the MIE is to separate the MSE into a sum of two terms: one for the local errors and the other for the global errors. In MFP or other non linear estimation technique like beamforming, the noise free ML criterion usually exhibits a characteristic mainlobe/sidelobes behavior. The error is then dominated by two kinds of errors: the local errors that concentrate around the mainlobe peak and outliers that concentrate around the sidelobe peaks. Consider the true parameter θ 0 and a discrete set of No parameter points {θ 1 , θ 2 . . . θ No } sampled at the sidelobe maxima of the noise free ML criterion, the conditional MSE of the ML estimator can be approximated as [17], [19]–[21]: i h ˆ − θ 0 )(θ ˆ − θ 0 )T Ey|θ0 (θ ≈

1−

No X

n=1

!

Pe (θ n |θ 0 )

× ΣCRB (θ 0 ) +

No X

Pe (θ n |θ 0 ) × (θ n − θ 0 )(θ n − θ 0 )T , (15)

n=1

where ΣCRB (θ 0 ) is the Cram`er-Rao bound at θ 0 , and Pe (θ n |θ 0 ) is the ML estimator pairwise error probability, i.e. the probability of deciding in favor of the parameter θ n given the true parameter θ 0 in the binary hypothesis test {θ n , θ 0 }. The pairwise error probability Pe (θ n |θ 0 ) is used as an approximation P o of the probability that the estimate falls on the sidelobe n and (1 − N n=1 Pe (θ n |θ 0 )) is used as an approximation of the probability that the estimate falls on the mainlobe (i.e the probability of local September 29, 2014

DRAFT

8

errors). The CRB serves as an MSE approximation of the ML estimator local errors. The CRB for the considered data models are presented in Appendix A. In a Bayesian framework, the global MSE is found by evaluating the expected value of the conditional MSE with respect to the a priori parameters probability density function pθ (θ). This method has shown to give rather accurate MSE prediction in the asymptotic and threshold regions [17], [19]–[21]. IV. Z IV-Z AKAI B OUND The Ziv-Zakai bound (ZZB) originally derived in [14] is a Bayesian bound that provides a bound on the MSE averaged over the a priori parameters probability density function. It is obtained by lower-bounding the probability of outage, i.e. the probability that the estimation error fall above a given threshold. Define
ˆ − θ . In the scalar parameter case the probability of outage Pr |ǫ| ≥ h the estimation error as ǫ = θ(y) 2 is lower bounded by [15]    Z ∞ h Pr |ǫ| ≥ 2 × min(pθ (θ), pθ (θ + h)) · Pmin (θ, θ + h)dθ , ≥V 2 −∞

(16)

where pθ (θ) is the a priori parameters probability density function, V is the valley filling function defined by V{f (h)} = max f (h + ξ) and Pmin (θ, θ + h) is any tractable lower bound for the probability of error ξ≥0

of the suboptimal decision scheme Decide H0 : θ Decide H1 : θ + h

h ˆ if θ(y) ≤θ+ 2 h ˆ if θ(y) >θ+ 2

(17)

associated to the binary hypothesis test {H0 : θ, H1 : θ + h} with equally likely hypotheses. The MSE bound is then obtained through the following equality [15]   Z  2 1 ∞ h Pr |ǫ| ≥ hdh, E ǫ = 2 0 2 so that

  Ey,θ ǫ2 ≥

Z

∞

V 0

Z

∞

(18)



min(pθ (θ), pθ (θ + h)) · Pmin (θ, θ + h)dθ hdh. −∞

(19)

The right hand side of expression (19) is the Ziv-Zakai bound. It has been extended by Bell et al. to handle the N -dimensional vector parameter case [15]  Z Z ∞  T min(pθ (θ), pθ (θ + ∆)) · Pmin (θ, θ + ∆)dθ h dh, V max a Σa ≥ T 0

September 29, 2014

∆:a ∆=h

(20)

RN

DRAFT

9

  where Σ = Ey,θ ǫǫT is the error covariance matrix and a is an arbitrary vector. When a is the unit

vector with a one in the nth position, expression (20) yields a bound on the MSE of the nth component of θ . As can be seen in expressions (19) and (20), the ZZB involves two integrations. Most of the time, and it is the case here, these integrals must be evaluated numerically. As for the probability Pmin (θ, θ + ∆), it can be expressed as Pmin (θ, θ + ∆) =

1 [Pe (θ|θ + ∆) + Pe (θ + ∆|θ)] , 2

(21)

where Pe (θ +∆|θ) and Pe (θ|θ +∆) are pairwise error probabilities. The computation of these probabilities represents the central part of the ZZB derivation. In the absence of nuisance parameters, these pairwise error probabilities are given by the error probabilities obtained from the Likelihood ratio test (LRT), which is optimal for testing {H1 : θ 1 , H0 : θ 0 }. When dealing with nuisance parameters, these pairwise error probabilities can be given by the error probabilities obtained from the Generalized Likelihood ratio test (GLRT) [35]. While the GLRT optimality is not guaranteed, it provides good results in evaluating the ZZB for performance analysis of the ML processors. V. M AXIMUM L IKELIHOOD

PAIRWISE ERROR PROBABILITY

The pairwise error probability Pe (θ 1 |θ 0 ) associated with the ML estimator is the kernel of the ZZB and MIE analysis, and is therefore a key of MFP performance analysis. In this section we will derive simple expressions for the ML pairwise error probabilities of the stochastic and the deterministic signal models presented herein. These expressions are derived in the framework of recent results on the cumulative distribution function (CDF) of quadratic forms in Gaussian random variables [36]. In fact, it can be seen from the expressions (8), (13) and (14) that the ML estimators presented here are all quadratic criteria in the observation vector y which is assumed to follow a Gaussian distribution. Therefore, the computation of the pairwise error probabilities Pe (θ 1 |θ 0 ) involves the comparison of quadratic forms in Gaussian variables. In the following subsection, we briefly review the theoretical results on quadratic forms in Gaussian variables that will be used to derive the pairwise error probability expressions. A. Cumulative distribution function of quadratic forms in Gaussian random variables Let A be a deterministic Hermitian symmetric matrix, and denote by A = VΛVH the eigendecomposition of A, Λ = diag(λ1 , λ2 , . . . , λP ), where λp are the eigenvalues of A, and V = [v1 , v2 , . . . , vP ],

September 29, 2014

DRAFT

10

where vp are the eigenvectors of A. For a Gaussian random vector x ∼ CN (µ, IP ), the CDF of the quadratic form xH Ax can be written in terms of the eigenvalues and eigenvectors as [36] Z ∞ q(jω+β) 1 e−c e H Pr (x Ax ≤ q) = dω, QP 2π −∞ jω + β p=1 (1 + λp (jω + β)) for some β > 0 such that 1 + βλp > 0 for all p, and with c =

λp |vpH µ|2 (jω+β) p=1 (1+λp (jω+β)) .

PP

(22)

Note that the CDF

is then characterized by a simple one-dimensional integral which can be evaluated numerically. For our MFP problem, as will be seen in the next subsections, the following quadratic form Q is of particular interest: Q=

L M X X

λ1m |z1m,l |2 + λ2m |z2m,l |2 ,

(23)

m=1 l=1

where z1m,l ∼ CN (µ1l,m , 1), z2m,l ∼ CN (µ2l,m , 1), l = 1...L, m = 1...M are independent random variables. Based on the above result, the CDF of Q can be expressed as Z ∞ q(jω+β) 1 e−c e Pr (Q ≤ q) = dω, QM 2π −∞ jω + β m=1 (1 + λ1m (jω + β))L (1 + λ2m (jω + β))L

with

c=

M X L X λ1m |µ1l,m |2 (jω + β)

m=1 l=1

+

(1 + λ1m (jω + β))

(24)

λ2m |µ2l,m |2 (jω + β) , (1 + λ2m (jω + β))

(25)

for some β > 0 such that 1 + βλ1m > 0, 1 + βλ2m > 0. These results are now used to derive the pairwise error probability expressions. B. Stochastic signal model The derivation and the computation of the pairwise error probability Pe (θ 1 |θ 0 ) for the MFP SML estimator of equation (8) has already been addressed by Xu et al. [21]. Nevertheless, a different approach that unifies the stochastic and deterministic models under the same formalism is presented here in the light of the above results on the CDF of quadratic forms in Gaussian variables. The pairwise error probability Pe (θ 1 |θ 0 ) of the SML estimator is given by Pe (θ 1 |θ 0 ) = Pr

M X

C1 (fm , θ 1 ) +

L M X X

C2 (fm , θ 1 ) · |¯ gH (fm , θ 1 )yl (fm )|2

m=1 l=1

m=1

−

M X

C1 (fm , θ 0 ) −

M X L X

C2 (fm , θ 0 ) · |¯ g (fm , θ 0 )yl (fm )| ≥ 0 , (26)

M X L X

Qm,l ,

H

2

m=1 l=1

m=1

!

where yl (fm ) is defined as in expression (1) for θ = θ 0 . Define the following quadratic form Q=

(27)

m=1 l=1 September 29, 2014

DRAFT

11

with h Qm,l = ylH (fm ) C2 (fm , θ 1 )¯ g(fm , θ 1 )¯ gH (fm , θ 1 ) Define also the scalar q=

M X

i − C2 (fm , θ 0 )¯ g(fm , θ 0 )¯ gH (fm , θ 0 ) yl (fm ). (28)

C1 (fm , θ 0 ) −

M X

C1 (fm , θ 1 ).

(29)

m=1

m=1

The pairwise error probability Pe (θ 1 |θ 0 ) can then be written as Pe (θ 1 |θ 0 ) = Pr (Q ≥ q) .

(30)

It is shown in Appendix C-A that the quadratic forms Qm,l can be expressed as Qm,l = λ1m |z1m,l |2 + λ2m |z2m,l |2 ,

(31)

where z1m,l ∼ CN (0, 1), z2m,l ∼ CN (0, 1) are independent random variables and λ1m , λ2m are the two non-zero eigenvalues of h g(fm , θ 1 )¯ gH (fm , θ 1 ) Σy 1/2 (fm , θ 0 ) C2 (fm , θ 1 )¯

i − C2 (fm , θ 0 )¯ g(fm , θ 0 )¯ gH (fm , θ 0 ) Σy 1/2 (fm , θ 0 ). (32)

Thus, using equations (23), (24), (25) and the above results, it can be shown that the pairwise error probability has the following expression Z ∞ q(jω+β) 1 1 e dω, Pe (θ 1 |θ 0 ) = 1 − Q L L 2π −∞ jω + β M m=1 (1 + λ1m (jω + β)) (1 + λ2m (jω + β))

(33)

for some β > 0 such that 1 + βλ1m > 0, 1 + βλ2m > 0.

The derivation of the eigenvalues λ1m , λ2m has been achieved by Xu et al. [21]. For the sake of completeness, it is recalled here in Appendix C-A. Ultimately, the eigenvalues can be expressed in terms of simple quantities. Denote ρ(fm ) the correlation of the Green’s functions g(fm , θ 1 ) and g(fm , θ 0 ) ρ(fm ) =

|gH (fm , θ 1 )g(fm , θ 0 )| . kg(fm , θ 1 )kkg(fm , θ 0 )k

(34)

Denote also γi (fm ) the SNR at the receiver for the frequency fm when the true parameter is θ i γi (fm ) =

September 29, 2014

σs2 (fm ) kg(fm , θ i )k2 . 2 (f ) σw m

(35)

DRAFT

12

The eigenvalues λ1m , λ2m , and the scalar q are expressed as follows " 1 λ1m ,2m = − × γ0 (fm ) − γ1 (fm ) + γ0 (fm )γ1 (fm )(1 − ρ2 (fm )) 2(1 + γ1 (fm )) ±

h

2 γ0 (fm ) + γ1 (fm ) + γ0 (fm )γ1 (fm )(1 − ρ2 (fm ))

(36)

# i1/2 − 4 × γ0 (fm )γ1 (fm )ρ2 (fm ) , q=L

M X

m=1

ln

1 + γ1 (fm ) . 1 + γ0 (fm )

(37)

Therefore, the computation of the exact SML pairwise error probability only requires the calculation of simple quantities and the numerical evaluation of a one-dimensional integral.

C. Deterministic signal model Just like the SML, the DML processors of equations (13) and (14) can be seen as quadratic forms in Gaussian variables. In this case however, the Gaussian variables are noncentral and to the best of our knowledge no closed-form expression exists for the pairwise error probability. This makes the computation of the exact pairwise error probability tedious. We provide a solution using the results presented in Section V-A. 1) Incoherent DML: The pairwise error probability Pe (θ 1 |θ 0 ) of the DML estimator of equation (13) is given by Pe (θ 1 |θ 0 ) = Pr

M X L X |¯ gH (fm , θ 1 )yl (fm )|2 2 (f ) σw m

m=1 l=1

−

Define the following quadratic from

Q=

L M X X

M X L X |¯ gH (fm , θ 0 )yl (fm )|2

m=1 l=1

2 (f ) σw m

!

≥0 .

Qm,l ,

(38)

(39)

m=1 l=1

where

hg ¯(fm , θ 1 )¯ gH (fm , θ 1 ) g ¯(fm , θ 0 )¯ gH (fm , θ 0 ) i yl (fm ). Qm,l = ylH (fm ) − 2 (f ) 2 (f ) σw σw m m

(40)

Pe (θ 1 |θ 0 ) = Pr (Q ≥ 0) .

(41)

The pairwise error probability Pe (θ 1 |θ 0 ) can then be written

As shown in Appendix C-B, the quadratic forms Qm,l can be expressed as Qm,l = λ1m |z1m,l |2 + λ2m |z2m,l |2 , September 29, 2014

(42) DRAFT

13

where z1m,l ∼ CN (µ1m,l , 1), z2m,l ∼ CN (µ2m,l , 1) are independent random variables. The mean values are µ1m,l =

sl (fm ) H σw (fm ) u1m gn (fm , θ 0 )

and µ2m,l =

sl (fm ) H σw (fm ) u2m gn (fm , θ 0 ),

where u1m , u2m are the

eigenvectors associated with the two non-zero eigenvalues λ1m , λ2m of g ¯(fm , θ 1 )¯ gH (fm , θ 1 ) − g ¯(fm , θ 0 )¯ gH (fm , θ 0 ).

(43)

Thus, using equations (23), (24), (25) and the above results, it can be shown that the pairwise error probability has the following expression Z ∞ 1 e−c 1 Pe (θ 1 |θ 0 ) = 1 − dω, Q L L 2π −∞ jω + β M m=1 (1 + λ1m (jω + β)) (1 + λ2m (jω + β)) c=

L M X X λ1m |µ1l,m |2 (jω + β)

m=1 l=1

(1 + λ1m (jω + β))

+

λ2m |µ2l,m |2 (jω + β) , (1 + λ2m (jω + β))

(44)

(45)

for some β > 0 such that 1 + βλ1m > 0, 1 + βλ2m > 0. The derivation of the eigenvectors u1m , u2m and the eigenvalues λ1m , λ2m is given in Appendix C-B. Ultimately, the eigenvalues and the mean values can be expressed in terms of the correlation function ρ(fm ) defined in (34) and the SNR at the receiver for the snapshot l and the frequency fm defined as γ(l, fm ) =

|sl (fm )|2 kg(fm , θ 0 )k2 . 2 (f ) σw m

The eigenvalues λ1m , λ2m , and the mean values µ1m,l , µ2m,l are expressed as follows p λ1m ,2m = ± 1 − ρ2 (fm ) r γ(l, fm )(1 − λ1m ,2m ) µ1m,l ,2m,l = ∓ . 2

(46)

(47) (48)

2) Coherent DML: The pairwise error probability Pe (θ 1 |θ 0 ) of the DML estimator of equation (14) is given by 

Pe (θ 1 |θ 0 ) = Pr  PM L m=1

1 kg(fm ,θ1 )k2 2 (f ) σw m

2 M X L X gH (fm , θ 1 )yl (fm ) 2 (f ) σw m m=1 l=1

− PM L m=1

1

kg(fm ,θ0 )k2 2 (f ) σw m

 M L X X gH (f , θ )y (f ) 2 m 0 l m ≥ 0 . (49) 2 (f ) σw m m=1 l=1

Denote x the modified set of observations  T T (f ) T (f ) T yL yL y1T (fM ) y1 (f1 ) 1 M ... ... ... , x= σw (f1 ) σw (f1 ) σw (fM ) σw (fM )

and g(θ) the corresponding modified set of Green’s functions  T T 1 gT (f1 , θ) gT (fM , θ) gT (fM , θ) g (f1 , θ) g(θ) = q P ... ... ... . kg(fm ,θ)k2 σw (f1 ) σw (f1 ) σw (fM ) σw (fM ) L M m=1 σ2 (fm )

(50)

(51)

w

September 29, 2014

DRAFT

14

We define the following quadratic from h i Q = xH g(θ 1 )gH (θ 1 ) − g(θ 0 )gH (θ 0 ) x.

(52)

The pairwise error probability Pe (θ 1 |θ 0 ) can then be written

Pe (θ 1 |θ 0 ) = Pr (Q ≥ 0) .

(53)

As shown in Appendix C-B, the quadratic form Q can be expressed as Q = λ1 |z1 |2 + λ2 |z2 |2 ,

(54)

where z1 ∼ CN (µ1 , 1), z2 ∼ CN (µ2 , 1) are independent random variables whose mean values µ1 and µ2 are functions of the eigenvectors u1 , u2 associated with the two non-zero eigenvalues λ1 , λ2 of g(θ 1 )gH (θ 1 ) − g(θ 0 )gH (θ 0 ).

(55)

Using equations (23), (24), (25) and the above results, it can be shown that the pairwise error probability has the following expression Pe (θ 1 |θ 0 ) = 1 −

1 2π

c=

Z

∞

−∞

e−c 1 dω, jω + β (1 + λ1 (jω + β))(1 + λ2 (jω + β))

λ2 |µ2 |2 (jω + β) λ1 |µ1 |2 (jω + β) + , (1 + λ1 (jω + β)) (1 + λ2 (jω + β))

(56) (57)

for some β > 0 such that 1 + βλ1 > 0, 1 + βλ2 > 0. The derivation of the eigenvectors u1 , u2 and the eigenvalues λ1 , λ2 is given in Appendix C-B. Ultimately, the eigenvalues and the mean values can be expressed as follows q λ1,2 = ± 1 − |gH (θ 1 )g(θ 0 )|2 , s P (1 − λ1,2 )L M m=1 γ(fm ) , µ1,2 = ∓ 2

(58) (59)

where γ(fm ) is the SNR at the receiver for the frequency fm γ(fm ) =

|s|2 kg(fm , θ 0 )k2 . 2 (f ) σw m

(60)

Similarly to the stochastic signal model, the computation of the exact DML pairwise error probabilities in the incoherent and coherent cases only requires the calculation of simple quantities and the numerical evaluation of a one-dimensional integral.

September 29, 2014

DRAFT

15

VI. A PPLICATION A. Computational issues 1) Evaluation of the pairwise error probability: The computation of the pairwise error probability Pe (θ 1 |θ 0 ) derived in this paper requires the numerical evaluation of an infinite one-dimensional integral

(see equations (33), (44) and (56)). In general, the integrand is highly concentrated around its maximum at ω = 0 and oscillates slightly while it decreases to zero. A typical example is represented in Figure 2. The numerical integration of this kind of integrands can be done with standard numerical integration algorithms. In our simulations we used the quadgk function of Matlab which showed good speed as well as good accuracy. Furthermore, the evaluation of the integral requires the assignment of a value to the almost free parameter β . This parameter must be chosen so that 0 < β <

1 . max(−λ1m ,2m )

The

value of β has an influence on the maximum value of the integrand at ω = 0 and on the decay rate of its oscillations. In practice, values of β too close to the upper limit lead to instabilities in the integral evaluation owing to a slow decay of the oscillation, whereas a choice of β included in the interval h i 1 1 , 1000×max(−λ1 ,2 ) 5×max(−λ1 ,2 ) enabled integration without difficulties. Also, the value of β that m

m

m

m

minimize the integrand at ω = 0 appears to be a very good choice to obtain a fast decay with a moderate maximum value at ω = 0. 2) MIE analysis: A critical step in MIE performance analysis is the selection of the sampling points θn of equation (15). This process can be done manually but this can be quite tedious especially if multiple

MIE analyses are performed over the a priori parameters probability density function pθ (θ) to obtain the global MSE. In our simulation we select the sampling points θn automatically by searching for the No = 12 highest local maxima on the noise free ML criterion (except for the mainlobe). The choice of No = 12 is based on the observed sidelobe structure of the noise free ML criterion. It should be noted

that the results were found to be little sensitive to the choice of No in the asymptotic and threshold regions (No = 8 or No = 20 lead to similar results) which may be explained by the results in [37]. Furthermore, the Cram`er-Rao bounds are computed using the results presented in Appendix A. 3) Evaluation of the ZZB: The two integrals of the ZZB are evaluated numerically on a predefined grid. The integration steps must be chosen with care. The first integration over the parameter θ in equations (19) or (20) does not need a coarse grid since the probability Pmin (θ, θ + h) is mildly dependent of the parameter θ , a large integration step can greatly improve the calculation time without damaging the accuracy. The second integration over h however may require a small integration step to fully capture errors in the asymptotic region. Typical integration steps are given in the following subsection.

September 29, 2014

DRAFT

16

B. Illustrations and discussions The theoretical results are illustrated for MFP source localization in a shallow water environment. The MFP configuration is presented in Figure 1. The waveguide is a Pekeris waveguide [38] which consists in an isospeed water layer overlaying a semi-infinite fluid basement. This waveguide is classically used to model shallow water environments. The Pekeris waveguide is chosen with the following parameters: depth of the water column D = 100 m, sound speed in the water cw = 1500 m.s−1, density in the water ρw = 1000 kg.m−3 , sound speed in the bottom cb = 1800 m.s−1 and density in the bottom ρw = 2000 kg.m−3 . The receiving array is a vertical array of N = 12 elements linearly spaced between z1 = 5 m and zN = 95 m. The source depth is assumed known and is fixed to zs = 30 m, whereas the

range is unknown with a uniform prior distribution on the interval r = [4000, 6000] m. We assume that we have L = 5 snapshots and M = 8 frequencies logarithmically spaced from 50 − 500 Hz. The Green’s functions are computed using normal mode theory [39]. The ZZB and the MIE are implemented with the results developed in Section V. The two numerical integrations of the ZZB in equation (19) are done with integration steps of δθ = 80 m and δh = 0.25 m for the stochastic and deterministic incoherent model, and with δh = 0.125 m for the deterministic coherent model. The MIE analyses are averaged over Nl = 25 locations on the uniform prior in order to approximate the global MSE. The number of location Nl was chosen so that higher values of Nl do not change significantly the result. The Bayesian CRB (BCRB) and the hybrid CRB (HCRB) presented in Appendix B are also evaluated. Monte-Carlo simulations averaged over Nl = 25 locations on the uniform prior are carried with Nc = 5000 iterations for each location and each SNR. For the stochastic model the SNR is defined by equation (35) while for the deterministic models it is defined by equations (46) and (60). Note that we take the same SNR for all frequencies and snapshots. The results are presented in Figure 3 for the deterministic incoherent model, in Figure 4 for the deterministic coherent model and in Figure 5 for the stochastic model. They are consistent with the general expectations from the ZZB and MIE analysis. The threshold phenomenon is well captured as well as the asymptotic Cram`er-Rao region. Further observations can be made from these simulations. First, one can see that the coherent model greatly departs from the other models. As might be expected, the performance is much better in the coherent case which gives an overview of the benefit that can be earned by incorporating more information on the source signal. The deterministic incoherent model and the stochastic model on the other hand provide close results. The threshold in the deterministic case is only 1.5 dB lower than in the stochastic case and the MSE are of the same order. This might be explained

September 29, 2014

DRAFT

17

by the fact that the underlying assumptions behind these two models are not so different, in both cases the signal is poorly known. VII. C ONCLUSIONS In this paper, we used recent results on the CDF of quadratic form in Gaussian variables to derive simple expressions for the pairwise error probabilities of the MFP ML estimators under both the stochastic and deterministic signal models. These expressions can be used to analyze MFP performance through the computation of the ZZB or through MIE analysis. The consistency of the results is assessed here with the computation of the ZZB and MIE performance analysis under the two signal models. It offers good prospects for many underwater acoustics applications. For instance, the results presented here can be easily adapted to handle a hybrid stochastic/deterministic signal, i.e., with both stochastic and deterministic frequency components. They can also be transposed to the Matched-Mode processing problem [40], [41]. Furthermore, our results could possibly be used to incorporate background mismatch in the analysis [42], [43] and analyze the impact of mismatched Green’s function on MFP. VIII. ACKNOWLEDGMENT This work was funded by the French Government Defense procurement agency (Direction G´en´erale de l’Armement). The authors are very grateful to the editor and the reviewers for their insightful comments that contributed to the improvement of the paper. A PPENDIX A C RAM E` R -R AO

BOUNDS

The CRB provides a lower bound on the variance of any unbiased estimator. The CRB states that i h ˆ − θ 0 )(θ ˆ − θ 0 )T ≥ ΣCRB (θ 0 ) = JD −1 (θ 0 ), (61) Ey|θ0 (θ

where ΣCRB (θ 0 ) is the CRB at θ 0 and JD (θ) is the local Fisher information matrix defined as: "   # ∂ ln p(y|θ) ∂ ln p(y|θ) H JD (θ) = −Ey|θ · . ∂θ ∂θ

(62)

The computation of the CRB requires the computation of the Fisher information matrix. In the following, expressions are given for the Fisher information matrix under the data models considered in this paper.

September 29, 2014

DRAFT

18

A. Stochastic model The expression of the local Fisher information matrix for the stochastic data model has been derived in [8]: [JD ]i,j (θ) = L

M X

m=1

"

 σs4 ǫ(fm , θ) · Re[d2 (fm , θ)li,j (fm , θ) − li (fm , θ)ljH (fm , θ)]

#  + ǫ(fm , θ)Re[li (fm , θ)]Re[lj (fm , θ)] , (63)

where d2 (fm , θ) = ǫ(fm , θ) =

kg(fm , θ)k2 , 2 (f ) σw m 2

1+

1

σs2 d2 (fm , θ)

(64) ,

∂g(fm , θ) , ∂θi   1 ∂g(fm , θ) H ∂g(fm , θ) li,j (fm , θ) = 2 . σw (fm ) ∂θi ∂θj li (fm , θ) =

2 (f ) σw m

gH (fm , θ)

(65) (66) (67)

The Green’s function does not have closed-form expressions. Hence, the above quantities and the Fisher information matrix inverse have to be evaluated numerically. B. Deterministic model In the deterministic model, the parameters to be estimated can be decomposed into the set of parameters p = [θ, α] where θ represents the Nθ parameters of interest and α represents the Nα nuisance parameters.

The nuisance parameters are complex valued and can be decomposed into α = [αre , αim ] where αre stands for the real part of the parameters and αim stands for the imaginary part. The local Fisher information matrix can then be written as the following block matrix:   Jθ,αim Jθ Jθ,αre     H  JD (p) =  Jθ,αre Jαre Jαre ,αim  .   H H Jαim Jθ,αim Jαre ,αim

September 29, 2014

(68)

DRAFT

19

1) Incoherent model: In the deterministic incoherent model, there are 2 × L × M nuisance parameters: h i αre = Re [s1 (f1 ) . . . s1 (fM ) . . . sL (f1 ) . . . sL (fM )] , (69)

and

im

α

h i = Im [s1 (f1 ) . . . s1 (fM ) . . . sL (f1 ) . . . sL (fM )] .

It can be shown from equations (12) and (62) that "   #  L X M X ∂sl (fm )g(fm , θ) H ∂sl (fm )g(fm , θ) 1 . Re 2 [JD ]i,j (p) = 2 σw (fm ) ∂pi ∂pj

(70)

(71)

l=1 m=1

Thus

[Jθ ]i,j = 2

L X M X

l=1 m=1

[Jθ,αre ]i,m+M (l−1)

"

|sl (fm )|2 Re 2 (f ) σw m

"



∂g(fm , θ) ∂θi

sl (fm )H = 2 × Re 2 (f ) σw m



H 

∂g(fm , θ) ∂θi

∂g(fm , θ) ∂θj

H

#

, ∀(i, j) ∈ {1, 2 . . . Nθ }2 ,

(72)

#

g(fm , θ) ,

∀(i, m, l) ∈ {1, 2 . . . Nθ } × {1, 2 . . . M } × {1, 2 . . . L}, (73)

[Jθ,αim ]i,m+M (l−1)

"

sl (fm )H = −2 × Im 2 (f ) σw m



∂g(fm , θ) ∂θi

H

#

g(fm , θ) ,

∀(i, m, l) ∈ {1, 2 . . . Nθ } × {1, 2 . . . M } × {1, 2 . . . L}, (74) [Jαre ,αim ] = 0(L×M )×(L×M ) ,

(75)

Jαre = Jαim = 2 × diag (v) ,

(76)

and

where v is the following vector concatenated L times:   kg(f1 , θ)k2 kg(fM , θ)k2 ... . 2 (f ) 2 (f ) σw σw 1 M

(77)

2) Coherent model: In the deterministic coherent model, there are only two nuisance parameters: αre = Re[s] and αim = Im[s]. The expression (71) is still true but with sl (fm ) = s. Thus, "   #  M X ∂g(fm , θ) H ∂g(fm , θ) |s|2 , ∀(i, j) ∈ {1, 2 . . . Nθ }2 , Re 2 [Jθ ]i,j = 2 × L σw (fm ) ∂θi ∂θj m=1 " #   M X sH ∂g(fm , θ) H Re 2 g(fm , θ) , ∀i ∈ {1, 2 . . . Nθ }, [Jθ,αre ]i = 2 × L σw (fm ) ∂θi

(78)

(79)

m=1

September 29, 2014

DRAFT

20

[Jθ,αim ]i = −2 × L

M X

m=1

"

sH Im 2 σw (fm )



∂g(fm , θ) ∂θi

H

#

g(fm , θ) , ∀i ∈ {1, 2 . . . Nθ },

(80)

[Jαre ,αim ] = 0,

Jαre = Jαim

(81)

M X kg(fm , θ)k2 =2×L . 2 (f ) σw m

(82)

m=1

Again, the above quantities and the Fisher information matrix inverse have to be evaluated numerically to obtain the CRB. The block structure of the Fisher information matrix can be taken into account in order to obtain directly the CRB on the parameters of interest. Using a blockwise inversion formula and the previous results, we get:   h i [JD −1 (p)]θ = Jθ − Jθ,αre Jθ,αim · 

J−1 αre

0(L×M )×(L×M )

0(L×M )×(L×M )

J−1 αre

  ·

Jθ,αre Jθ,αim

−1 

.

(83)

A PPENDIX B BAYESIAN C RAM E` R -R AO

BOUNDS

When the parameter set θ to be estimated is assumed to be random with a probability density function pθ (θ), the Bayesian CRB (BCRB) provides a bound for any estimator on the MSE averaged over the

prior pθ (θ) [44]. When the parameter set p = [θ, α] to be estimated has both a random part θ and a deterministic part α the bound can be provided by the hybrid CRB (HCRB) [45], [46]. The BCRB applies to the stochastic data model while the HCRB applies to the deterministic data model.

A. Stochastic model The BCRB states that [44] i h ˆ − θ)(θ ˆ − θ)T ≥ ΣBCRB = JB −1 , Ey,θ (θ

(84)

where ΣBCRB is the BCRB and JB is the Bayesian information matrix. The Bayesian information matrix is related to the Fisher information matrix JD (θ) through the relation:

where JA

JB = Eθ [JD (θ)] + JA , (85) h 2 i ln pθ (θ) is the a priori information matrix, [JA ]i,j = Eθ ∂ ∂θ . The BCRB for the stochastic i ∂θj

model can therefore be computed using the results in Appendix A-A.

September 29, 2014

DRAFT

21

B. Deterministic model Denote, as in Appendix A-B, p = [θ, α] the parameter set to be estimated. The parameters of interest θ are random while the nuisance parameters α are deterministic. The HCRB has to be used instead of

the BCRB. The HCRB states that [45], [46]   Ey,θ|α0 (ˆ p − p)(ˆ p − p)T ≥ ΣHCRB (α0 ) = JH −1 (α0 ),

(86)

where ΣHCRB (α0 ) is the HCRB for the set of nuisance parameters α0 and JH (α) is the Hybrid information matrix. When θ and α are independent, the Hybrid information matrix is related to the Fisher information matrix JD (p) through the relation:  JH (α) = Eθ|α [JD (p)] + 

JA

0Nθ ×Nα

0Nα ×Nθ

0Nα ×Nα



.

(87)

The HCRB for the deterministic model can therefore be computed using the results in Appendix A-B. A PPENDIX C E IGENVECTORS

AND EIGENVALUES

A. Stochastic model Consider the quadratic form Qm,l defined in Section V-B. Let xm,l be defined such that yl (fm ) = Σy 1/2 xm,l . This operation is actually a whitening transform and xm,l ∼ CN (0, IN ). The quadratic form Qm,l can then be written

h 1/2 C2 (fm , θ 1 )¯ g(fm , θ 1 )¯ gH (fm , θ 1 ) Σ Qm,l = xH m,l y

i − C2 (fm , θ 0 )¯ g(fm , θ 0 )¯ gH (fm , θ 0 ) Σy 1/2 xm,l . (88)

Denote UΛUH the eigendecomposition of h i g(fm , θ 1 )¯ gH (fm , θ 1 ) − C2 (fm , θ 0 )¯ g(fm , θ 0 )¯ gH (fm , θ 0 ) Σy 1/2 . Σy 1/2 C2 (fm , θ 1 )¯

(89)

It follows that

Qm,l = zH m,l Λzm,l ,

(90)

where zm,l = UH xm,l . Since UH U = IN , we have zm,l ∼ CN (0, IN ). Furthermore, the matrix defined in (89) is rank two, hence Λ = diag(λ1m , λ2m , 0...0) and the quadratic form Qm,l reduce to Qm,l = λ1m |z1m,l |2 + λ2m |z2m,l |2 ,

(91)

where z1m,l ∼ CN (0, 1), z2m,l ∼ CN (0, 1) are independent random variables. September 29, 2014

DRAFT

22

The matrix defined in expression (89) shares the same eigenvalues with h i C2 (fm , θ 1 )¯ g(fm , θ 1 )¯ gH (fm , θ 1 ) − C2 (fm , θ 0 )¯ g(fm , θ 0 )¯ gH (fm , θ 0 ) Σy .

(92)

Since it is a rank two matrix, the eigenvectors corresponding to the two non-zero eigenvalues are linear combinations of g ¯(fm , θ 1 ) and g ¯(fm , θ 0 ). Therefore h i C2 (fm , θ 1 )¯ g(fm , θ 1 )¯ gH (fm , θ 1 ) − C2 (fm , θ 0 )¯ g(fm , θ 0 )¯ gH (fm , θ 0 ) Σy h i h i · g ¯(fm , θ 0 ) + β¯ g(fm , θ 1 ) = λ g ¯(fm , θ 0 ) + β¯ g(fm , θ 1 ) . (93)

Replacing Σy by its expression in (7) gives an expression of the form a¯ g(fm , θ 0 ) + b¯ g(fm , θ 1 ) = 0. Since g ¯(fm , θ 0 ) and g ¯(fm , θ 1 ) are not collinear, both a and b must be zero. Therefore, we get two equations with two unknowns β, λ, and finally obtain the following equation for λ (1 + γ1 (fm ))λ2 + (γ0 (fm ) − γ1 (fm ) + γ0 (fm )γ1 (fm )(1 − ρ2 (fm )))λ

− γ0 (fm )γ1 (fm )(1 − ρ2 (fm )) = 0, (94)

where γi (fm ) and ρ(fm ) are defined in Section V-B. Solving for this equation gives the two eigenvalues in (36).

B. Deterministic model 1) Incoherent model: Consider the quadratic form Qm,l defined in Section V-C. Let xm,l be defined such that xm,l =

yl (fm ) σw (fm ) ,

we have xm,l ∼ CN ( σswl (f(fmm)) g(fm , θ 0 ), IN ). The quadratic form Qm,l can then

be written h i H H g ¯ (f , θ )¯ g (f , θ ) − g ¯ (f , θ )¯ g (f , θ ) Qm,l = xH m 1 m 1 m 0 m 0 xm,l . m,l

(95)

g ¯(fm , θ 1 )¯ gH (fm , θ 1 ) − g ¯(fm , θ 0 )¯ gH (fm , θ 0 ).

(96)

Qm,l = zH m,l Λzm,l ,

(97)

Let UΛUH denote the eigendecomposition of

It follows that

where zm,l = UH xm,l . Since UH U = IN , we have zm,l ∼ CN (UH σswl (f(fmm)) g ¯(fm , θ 0 ), IN ). Furthermore, the matrix defined in expression (96) is rank two, hence Λ = diag(λ1m , λ2m , 0...0) and the quadratic form Qm,l reduce to Qm,l = λ1m |z1m,l |2 + λ2m |z2m,l |2 , September 29, 2014

(98) DRAFT

23

where z1m,l ∼ CN (µ1m,l , 1), z2m,l ∼ CN (µ2m,l , 1) are independent random variables with mean values µ1m,l =

sl (fm ) H ¯(fm , θ 0 ) σw (fm ) u1m g

and µ2m,l =

sl (fm ) H ¯(fm , θ 0 ). σw (fm ) u2m g

Since the matrix defined in expression (96) is a rank two matrix, the eigenvectors corresponding to the two non-zero eigenvalues are linear combinations of g ¯(fm , θ 1 ) and g ¯(fm , θ 0 ). Therefore h i h i g ¯(fm , θ 1 )¯ gH (fm , θ 1 ) − g ¯(fm , θ 0 )¯ gH (fm , θ 0 ) · g ¯(fm , θ 0 ) + β¯ g(fm , θ 1 ) h i =λ g ¯(fm , θ 0 ) + β¯ g(fm , θ 1 ) . (99)

Developing the above expression gives an expression of the form a¯ g(fm , θ 0 ) + b¯ g(fm , θ 1 ) = 0. Since g ¯(fm , θ 0 ) and ¯ g(fm , θ 1 ) are not collinear, both a and b must be zero. Therefore, we get the two following

equations

Solving for λ gives

  −λβ + β + ¯ gH (fm , θ 1 )¯ g(fm , θ 0 ) = 0  1 + β¯ gH (fm , θ 0 )¯ g(fm , θ 1 ) + λ = 0 λ1m ,2m = ±

and

p 1 − ρ2 (fm ),

(100)

g ¯H (fm , θ 1 )¯ g(fm , θ 0 ) . λ−1

(101)

g ¯(fm , θ 0 ) + β1m ,2m g ¯(fm , θ 1 ) , k¯ g(fm , θ 0 ) + β1m ,2m g ¯(fm , θ 1 )k

(102)

β1m ,2m =

Thus, we also have the eigenvectors u1m ,2m =

and after some straightforward manipulation we get the mean values r γ(l, fm )(1 − λ1m ,2m ) µ1m,l ,2m,l = −sign(λ1m ,2m ) . 2

(103)

2) Coherent model: Consider the vector x and the quadratic form Q defined Section V-C. The vector x follows x ∼ CN (mx , IN ) where  T s · gT (f1 , θ 0 ) s · gT (f1 , θ 0 ) s · gT (fM , θ 0 ) s · gT (fM , θ 0 ) mx = ... ... ... . σw (f1 ) σw (f1 ) σw (fM ) σw (fM )

(104)

Denote UΛUH the eigendecomposition of g(θ 1 )gH (θ 1 ) − g(θ 0 )gH (θ 0 ).

(105)

Q = zH Λz,

(106)

It follows that

September 29, 2014

DRAFT

24

where z = UH x. Since UH U = IN , z ∼ CN (UH mx , IN ). Furthermore, the matrix defined in expression (105) is rank two, hence Λ = diag(λ1 , λ2 , 0 . . . 0) and the quadratic form Q reduce to Q = λ1 |z1 |2 + λ2 |z2 |2 ,

(107)

where z1 ∼ CN (µ1 , 1), z2 ∼ CN (µ2 , 1) are independent random variables with mean values µ1 = uH 1 mx and µ2 = uH 2 mx . Since the matrix defined in expression (105) is a rank two matrix, the eigenvectors corresponding to the two non-zero eigenvalues are linear combinations of g(θ 1 ) and g(θ 0 ). Therefore h i h i h i g(θ 1 )gH (θ 1 ) − g(θ 0 )gH (θ 0 ) · g(θ 0 ) + βg(θ 1 ) = λ g(θ 0 ) + βg(θ 1 ) .

(108)

Developing the above expression gives an expression of the form ag(θ 0 ) + bg(θ 1 ) = 0. Since g(θ 0 ) and g(θ 1 ) are not collinear, both a and b must be zero. Therefore, we get the two following equations   −λβ + βkg(θ )k2 + gH (θ )g(θ ) = 0 1 1 0  2 H kg(θ )k + βg (θ )g(θ ) + λ = 0 0

0

1

Since kg(θ)k2 = 1, solving for λ gives

λ1,2

q = ± 1 − |gH (θ 1 )g(θ 0 )|2 ,

(109)

and β1,2 = −

gH (θ 1 )g(θ 0 ) . 1 − λ1,2

(110)

Thus, we also have the eigenvectors u1,2 =

g(θ 0 ) + β1,2 g(θ 1 ) , kg(θ 0 ) + β1,2 g(θ 1 )k

and after some straightforward manipulation we get the mean values s P (1 − λ1,2 )L M m=1 γ(fm ) µ1,2 = −sign(λ1,2 ) . 2

(111)

(112)

R EFERENCES [1] A. Tolstoy, Matched field processing for underwater acoustics, vol. 52, World Scientific Singapore, 1993. [2] A.B. Baggeroer, W.A. Kuperman, and P.N. Mikhalevsky, “An overview of matched field methods in ocean acoustics,” IEEE J. Ocean. Eng., vol. 18, no. 4, pp. 401–424, 1993. [3] C.F. Mecklenbr¨auker and P. Gerstoft, “Objective functions for ocean acoustic inversion derived by likelihood methods,” Journal of Computational Acoustics, vol. 8, no. 02, pp. 259–270, 2000. [4] S.E. Dosso and M.J. Wilmut, “Maximum-likelihood and other processors for incoherent and coherent matched-field localization,” The Journal of the Acoustical Society of America, vol. 132, pp. 2273–2285, 2012.

September 29, 2014

DRAFT

25

[5] W. Xu, Performance bounds on matched-field methods for source localization and estimation of ocean environmental parameters, Ph.D. thesis, MIT, Cambridge, MA, 2001. [6] A. Thode, M. Zanolin, E. Naftali, I. Ingram, P. Ratilal, and N.C. Makris, “Necessary conditions for a maximum likelihood estimate to become asymptotically unbiased and attain the cramer–rao lower bound. ii. range and depth localization of a sound source in an ocean waveguide,” The Journal of the Acoustical Society of America, vol. 112, pp. 1890–1910, 2002. [7] M. Zanolin, I. Ingram, A. Thode, and N.C. Makris, “Asymptotic accuracy of geoacoustic inversions,” The Journal of the Acoustical Society of America, vol. 116, pp. 2031–2042, 2004. [8] A.B. Baggeroer and H. Schmidt, “Cramer-rao bounds for matched field tomography and ocean acoustic tomography,” in Proc. ICASSP. IEEE, 1995, vol. 5, pp. 2763–2766. [9] J. Tabrikian and H. Messer, “Three-dimensional source localization in a waveguide,” IEEE Trans. Signal Process., vol. 44, no. 1, pp. 1–13, 1996. [10] W. Xu and J. Li, “Study of statistical signal models in low-frequency underwater acoustic applications,” in Proc. OCEANS 11. IEEE, 2011, pp. 1–5. [11] E.W. Barankin, “Locally best unbiased estimates,” The Annals of Mathematical Statistics, pp. 477–501, 1949. [12] J. Tabrikian and J.L. Krolik, “Barankin bounds for source localization in an uncertain ocean environment,” IEEE Trans. Signal Process., vol. 47, no. 11, pp. 2917–2927, 1999. [13] E. Weinstein and A.J. Weiss, “A general class of lower bounds in parameter estimation,” IEEE Trans. Inf. Theory, vol. 34, no. 2, pp. 338–342, 1988. [14] J. Ziv and M. Zakai, “Some lower bounds on signal parameter estimation,” IEEE Trans. Inf. Theory, vol. 15, no. 3, pp. 386–391, 1969. [15] K.L. Bell, Y. Steinberg, Y. Ephraim, and H.L. Van Trees, “Extended ziv-zakai lower bound for vector parameter estimation,” IEEE Trans. Inf. Theory, vol. 43, no. 2, pp. 624–637, 1997. [16] W. Xu, A.B. Baggeroer, and C.D. Richmond, “Bayesian bounds for matched-field parameter estimation,” IEEE Trans. Signal Process., vol. 52, no. 12, pp. 3293–3305, 2004. [17] W. Xu, A.B. Baggeroer, and H. Schmidt, “Performance analysis for matched-field source localization: Simulations and experimental results,” IEEE J. Ocean. Eng., vol. 31, no. 2, pp. 325–344, 2006. [18] T. Meng and J.R. Buck, “Rate distortion bounds on passive sonar performance,” IEEE Trans. Signal Process., vol. 58, no. 1, pp. 326–336, 2010. [19] F. Athley, “Threshold region performance of maximum likelihood direction of arrival estimators,” IEEE Trans. Signal Process., vol. 53, no. 4, pp. 1359–1373, 2005. [20] C.D. Richmond, “Mean-squared error and threshold snr prediction of maximum-likelihood signal parameter estimation with estimated colored noise covariances,” IEEE Trans. Inf. Theory, vol. 52, no. 5, pp. 2146–2164, 2006. [21] W. Xu, Z. Xiao, and L. Yu, “Performance analysis of matched-field source localization under spatially correlated noise field,” IEEE J. Ocean. Eng., vol. 36, no. 2, pp. 273–284, 2011. [22] V. Bush, J.B. Conant, and J.T. Tate, Principles and applications of underwater sound, Defense Technical Information Center, 1946. [23] M. Nicholas, J.S. Perkins, G.J. Orris, L.T. Fialkowski, and G.J. Heard, “Environmental inversion and matched-field tracking with a surface ship and an l-shaped receiver array,” The Journal of the Acoustical Society of America, vol. 116, pp. 2891– 2901, 2004. [24] S.A. Stotts, R.A. Koch, S.M. Joshi, V.T. Nguyen, V.W. Ferreri, and D.P. Knobles, “Geoacoustic inversions of horizontal

September 29, 2014

DRAFT

26

and vertical line array acoustic data from a surface ship source of opportunity,” IEEE J. Ocean. Eng., vol. 35, no. 1, pp. 79–102, 2010. [25] A. Fredricks, J.A. Colosi, J.F. Lynch, G. Gawarkiewicz, C.-S. Chiu, and Ph. Abbot, “Analysis of multipath scintillations from long range acoustic transmissions on the new england continental slope and shelf,” The Journal of the Acoustical Society of America, vol. 117, pp. 1038–1057, 2005. [26] D. Tollefsen and S.E. Dosso, “Bayesian geoacoustic inversion of ship noise on a horizontal array,” The Journal of the Acoustical Society of America, vol. 124, pp. 788–795, 2008. [27] A.V. van Leijen, J.P. Hermand, and M. Meyer, “Geoacoustic inversion in the north-eastern caribbean using a hydrographic survey vessel as a sound source of opportunity,” Journal of Marine Systems, vol. 78, pp. S333–S338, 2009. [28] C.-F. Huang and W.S. Hodgkiss, “Matched-field geoacoustic inversion of low-frequency source tow data from the asiaex east china sea experiment,” IEEE J. Ocean. Eng., vol. 29, no. 4, pp. 952–963, 2004. [29] A.M. Thode, G.L. D’Spain, and W.A. Kuperman, “Matched-field processing, geoacoustic inversion, and source signature recovery of blue whale vocalizations,” The Journal of the Acoustical Society of America, vol. 107, pp. 1286–1300, 2000. [30] S.E. Dosso, “Quantifying uncertainty in geoacoustic inversion. i. a fast gibbs sampler approach,” The Journal of the Acoustical Society of America, vol. 111, pp. 129–142, 2002. [31] S.E. Dosso and P.L. Nielsen, “Quantifying uncertainty in geoacoustic inversion. ii. application to broadband, shallow-water data,” The Journal of the Acoustical Society of America, vol. 111, pp. 143–159, 2002. [32] S.E. Dosso and M.J. Wilmut, “Quantifying data information content in geoacoustic inversion,” IEEE J. Ocean. Eng., vol. 27, no. 2, pp. 296–304, 2002. [33] Y. Pawitan, In all likelihood: statistical modelling and inference using likelihood, Oxford University Press, 2001. [34] H.L. Van Trees and K.L. Bell, Bayesian bounds for parameter estimation and nonlinear filtering/tracking, Wiley-IEEE Press, 2007. [35] N. Decarli and D. Dardari, “Ziv-zakai bound for time delay estimation of unknown deterministic signals,” in Proc. ICASSP, 2014. [36] T.Y. Al-Naffouri and B. Hassibi, “On the distribution of indefinite quadratic forms in gaussian random variables,” in Proc. Int. Symp. Inf. Theory (ISIT’09). IEEE, 2009, pp. 1744–1748. [37] D. Meng, W. Xu, and M. Xia, “Understanding the method of interval errors from the information theory perspective,” in Proc. ICASSP, 2009, pp. 2117–2120. [38] C.L. Pekeris, Theory of propagation of explosive sound in shallow water, vol. 27, Geol. Soc. Am. Mem., 1948. [39] F.B. Jensen, W.A. Kuperman, M.B. Porter, and H. Schmidt, Computational ocean accoustics, 2nd ed, New York: American Institute of Physics, 2011. [40] J. Bonnel, G. Le Touz´e, B. Nicolas, and J. Mars, “Physics-based time-frequency representations for underwater acoustics: Power class utilization with waveguide-invariant approximation,” IEEE Signal Process. Mag., vol. 30, no. 6, pp. 120–129, 2013. [41] G. Le Touz´e, J. Torras, B. Nicolas, and J. Mars, “Source localization on a single hydrophone,” in Proc. OCEANS 2008. IEEE, 2008, pp. 1–6. [42] W. Xu, A.B. Baggeroer, and K.L. Bell, “A bound on mean-square estimation error with background parameter mismatch,” IEEE Trans. Inf. Theory, vol. 50, no. 4, pp. 621–632, 2004. [43] C.D. Richmond, “Capon algorithm mean-squared error threshold snr prediction and probability of resolution,” IEEE Trans. Signal Process., vol. 53, no. 8, pp. 2748–2764, 2005.

September 29, 2014

DRAFT

27

[44] H.L. Van Trees, Detection, Estimation, and Modulation Theory Part I, Wiley, 1968. [45] Y. Rockah and P.M. Schultheiss, “Array shape calibration using sources in unknown locations–part i: Far-field sources,” IEEE Trans. Acoust., Speech, Signal Process., vol. 35, no. 3, pp. 286–299, 1987. [46] H. Messer, “The hybrid cramer-rao lower bound-from practice to theory,” in in Proc. of IEEE Workshop on Sensor Array and Multi-channel Processing (SAM). IEEE, 2006, pp. 304–307.

Yann Le Gall (S’12) received the Dipl.Ing. degree in signal processing from Grenoble Institut National Polytechnique (Grenoble INP), Grenoble, France, in 2012. PLACE PHOTO HERE

He is currently working toward the Ph.D. degree at Lab-STICC (UMR 6285), ENSTA Bretagne in Brest, France. His research interests in the field of signal processing and underwater acoustics include performance prediction, source localization and geoacoustic inversion.

Francois-Xavier Socheleau (S’08-M’12) graduated in electrical engineering from ESEO, Angers, France, in 2001 and received the Ph.D. degree from Telecom Bretagne, Brest, France, in 2011. PLACE PHOTO HERE

From 2001 to 2004, he was a Research Engineer at Thales Communications, France, where he worked on electronic warfare systems. From 2005 to 2007, he was employed at Navman Wireless (New Zealand/U.K.) as an R&D Engineer. From 2008 to 2011, he worked for Thales Underwater Systems, France. In November 2011, he joined ENSTA Bretagne as an Assistant Professor. Since September 2014, he has been an

Associate Professor at Telecom Bretagne, France. His research interests lie in the field of signal processing and underwater acoustics.

Julien Bonnel (S’08 - M’11) received the Ph.D. degree in signal processing from Grenoble Institut National Polytechnique (Grenoble INP), Grenoble, France, in 2010. PLACE PHOTO HERE

Since 2010, he has been an Assistant Professor at Lab-STICC (UMR 6285), ENSTA Bretagne in Brest, France. His research in signal processing and underwater acoustics include time-frequency analysis, source detection/localization, geoacoustic inversion, acoustical tomography, passive acoustic monitoring, and bioacoustics. He is a Member of the the Acoustical Society of America.

September 29, 2014

DRAFT

28

Fig. 1.

Matched-field processing configuration.

300

Real part of the integrand

250 200 150 100 50 0 −50 −100 −0.2

Fig. 2.

−0.15

−0.1

−0.05

0

ω

0.05

0.1

0.15

0.2

Example of integrand for the pairwise error probability computation. Only the real part of the integrand is represented

here since the imaginary part does not contribute to the integral (the result is real).

September 29, 2014

DRAFT

29

Deterministic incoherent model

Global Mean−Square Error (dB)

70

Monte Carlo HCRB ZZB MIE

60

50

40

30

20

10

0 No information region −10 −20

−15

−10

Threshold region −5

0

Asymptotic region 5

10

15

SNR(dB) Fig. 3.

Hybrid Cram`er-Rao bound (HCRB), Ziv-Zakai bound (ZZB) and method of interval errors (MIE) for source range

estimation under the deterministic incoherent model. Dash dot line: HCRB; Solid line: ZZB; Dashed line: MIE; *: Monte-Carlo simulation.

Deterministic coherent model 70

Monte Carlo HCRB ZZB MIE

Global Mean−Square Error (dB)

60 50 40 30 20 10 0 −10 −20 −30 −25

−20

−15

−10

−5

0

5

10

SNR(dB) Fig. 4.

Hybrid Cram`er-Rao bound (HCRB), Ziv-Zakai bound (ZZB) and method of interval errors (MIE) for source range

estimation under the deterministic coherent model. Dash dot line: HCRB; Solid line: ZZB; Dashed line: MIE; *: Monte-Carlo simulation.

September 29, 2014

DRAFT

30

Stochastic model

Global Mean−Square Error (dB)

70

Monte Carlo BCRB ZZB MIE

60

50

40

30

20

10

0

−10 −20

−15

−10

−5

0

5

10

15

SNR(dB) Fig. 5.

Bayesian Cram`er-Rao bound (BCRB), Ziv-Zakai bound (ZZB) and method of interval errors (MIE) for source range

estimation under the stochastic model. Dash dot line: BCRB; Solid line: ZZB; Dashed line: MIE; *: Monte-Carlo simulation.

September 29, 2014

DRAFT