Bayesian source localization with uncertain Green’s function Yann Le Gall∗ , Stan E. Dosso†, Franc¸ois-Xavier Socheleau‡ and Julien Bonnel∗ ∗ ENSTA

Bretagne, UMR CNRS 6285 Lab-STICC 2 rue Franc¸ois Verny, 29806 Brest Cedex 9, France Email: yann.le [email protected] † School of Earth and Ocean Sciences, University of Victoria Victoria, British Columbia V8W 3P6, Canada ‡ Institut Mines-Telecom/Telecom Bretagne, UMR CNRS 6285 Lab-STICC Technopˆole Brest-Iroise, 29238 Brest, France Abstract—The localization of an acoustic source in the oceanic waveguide is a difficult task because the oceanic environment is often poorly known. Uncertainty in the environment results in uncertainty in the source position and poor localization results. Hence, localization methods dealing with environmental uncertainty are required. In this paper, a Bayesian approach to source localization is introduced in order to improve robustness and obtain quantitative measures of localization uncertainty. The Green’s function of the waveguide is considered as an uncertain random variable whose probability density accounts for environmental uncertainty. The uncertain distribution over range and depth is then obtained through the integration of the posterior probability density (PPD) over the Green’s function probability density. An efficient integration technique makes the whole localization process computationally efficient. Some results are presented for a simple uncertain Green’s function model to show the ability of the proposed method to give reliable PPDs.

I.

I NTRODUCTION

Matched-Field processing (MFP) has received significant attention in underwater acoustics over the past decades. It allows estimation of source position by comparing data measured on an array of spatially distributed sensors with simulated replicas of the acoustic field derived from the wave equation and a model of the oceanic waveguide [1]. However, poor localization performance may be obtained when there is some mismatch between the true environment and the assumed oceanic waveguide (e.g. errors in sound speed profile, bathymetry or bottom properties) [2]–[5]. To obtain a quantitative measure of localization uncertainty and to improve robustness toward mismatch, environmental uncertainty can be included into the source localization MFP algorithm within a Bayesian framework. Bayes’ rule relates the a priori probability density and the likelihood to the posterior probability density (PPD) which depicts the probability of the possible source locations given the measured data. Environmental uncertainty is generally addressed by incorporating some environmental parameters as additional unknowns in the localization problem [6]–[9]. The analysis is then carried out from an extended PPD involving the localization parameters as well as the unknown environmental parameters. The integration of the PPD over the environmental unknown parameters, known as the marginalization technique, gives the joint marginal probability distribution over source range and 978-1-4799-8736-8/15/$31.00 ©2015 IEEE

depth. However, this approach relies on a parametrization of the waveguide and can be computationally intensive when there are many unknown parameters. In this paper we propose an original Bayesian approach to address the mismatch problem. No environmental parameters are added to the localization problem. The idea is to relax the assumption that the Green’s function of the supposed oceanic waveguide is a deterministic quantity. Environmental uncertainty is handled by considering its rough effect on the Green’s function. The Green’s function is considered as a random vector that accounts for environmental uncertainty. Integrating the PPD over the uncertain Green’s function gives the PPD over source range and depth. This integral is reduced to a one dimensional integral which greatly speeds up the integration and makes the whole process computationally efficient. Another strength of this method is that uncertainty is not limited to specific kinds of parameter mismatch. In the next section the Bayesian MFP theory is reviewed. Then our Bayesian MFP approach is introduced and finally some results are presented. II.

BAYESIAN

MATCHED - FIELD PROCESSING

A. Bayesian inference In a Bayesian context, assumptions about an unknown parameter set θ are made through the study of the conditional probability density function p(θ|y) of θ given the measured data y. This quantity is called the posterior probability density (PPD). It can be obtained using Bayes’ rule which may be written as p(θ|y) ∝ p(y|θ)p(θ), (1) where p(y|θ) is the probability of observing data y given the parameter set θ and p(θ) is the prior information about the unknown parameter set θ. The probability p(y|θ) defines the likelihood function when interpreted as a function of θ. The PPD describes the state of knowledge on θ. Using both data information and prior information it applies all available information to make inferences on θ. The PPD can be analyzed to obtain estimates of the parameters using several Bayesian estimation techniques such as the maximum a posteriori (MAP) estimator or the posterior mean estimator. The MAP is the minimum probability of error estimator when the parameters

have discrete probability distributions whereas the posterior mean is the minimum mean-squared error estimator (MMSE). The MAP estimate is θˆ = arg max p(θ|y), θ

and the posterior mean estimate is  θ¯ = θp(θ|y)dθ.

(2)

(3)

Deeper analysis of the PPD can provide information about the reliability of the estimate and other likely solutions to the inverse problem. Thus, in a Bayesian framework the inverse problem is not limited to a point estimate with no measure of confidence. B. Likelihood functions The statistical distribution of the data given the unknown parameters determines the likelihood function. In MFP, the data are the output of an N -element array measuring the waveforms from an acoustic source at a distance in the waveguide. 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

the Cholesky decomposition Cw (fm ) = Lw (fm )LH w (fm ). Therefore, the whitening transform ˜ l (fm ) = L−1 y w (fm )yl (fm ), ˜ (fm , θ) = L−1 g w (fm )g(fm , θ), (6) can be applied if the matrices Cw (fm ) are known, which allows reduction to the white noise case. To simplify subsequent derivations and analyses, in the remainder of this paper we will work in whitened space and consider that 2 Σw (fm ) = σw (fm )IN . The source signal was considered as an unknown quantity s(fm ) = sl (fm ) with a Rayleigh distributed amplitude and a uniform phase in the first development of Bayesian MFP by Richardson and Nolte [6]. However, a deterministic model where the source term sl (fm ) is regarded as a deterministic quantity appears to be more flexible and prevails in more recent MFP applications [7]–[12]. This deterministic source signal model will be used in what follows. Thus, the likelihood may be expressed as p(y|θ) = M

m=1

L M   m=1 l=1

(4)

where •

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

•

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

•

g(fm , θ) is a complex N × 1 vector corresponding to the transfer function or Green’s function of the waveguide at frequency fm for the propagation between the source and each of the N receivers of the array.

•

wl (fm ) is a complex N × 1 vector representing circularly-symmetric, zero mean Gaussian noise independent of the source signal. It is assumed to be white in time and in frequency with positive definite spatial 2 covariance matrix Σw (fm ) = σw (fm )Cw (fm ). This noise model is suited to the ambient noise which is always present in realistic applications. Its interpretation is often broadened to modeling errors like environmental mismatch [10], [11]. However, since modeling errors do not generally behave as a temporally independent additive Gaussian noise, this interpretation may be challenged, in particular when these errors are large. In this paper, an uncertain Green’s function is considered to address the mismatch problem.

The whole dataset y is the concatenation of all snapshots and frequencies T T y = [y1T (f1 ) . . . yL (f1 ) . . . y1T (fM ) . . . yL (fM )]T . (5)

The matrices Cw (fm ) are positive definite, their square root decomposition are invertible and can be obtained through

L

1

2 N l=1 |πσw (fm )|

×

  yl (fm ) − sl (fm )g(fm , θ)2 . exp − 2 (f ) σw m

(7)

In realistic applications of MFP the source terms sl (fm ) 2 and the noise variances σw (fm ) are generally not known and must be treated as nuisance parameter in the Bayesian inversion process. A few techniques have been investigated to handle these nuisance parameters [10], [11]. A convenient and reliable method involves replacing the source terms and the noise variances by their maximum likelihood estimate [7]– [12]. Depending upon the information available on the noise variance, two likelihood functions may be used. When the noise variances are known and the source terms are unknown, the likelihood is expressed as [13, Eq. (3.3)]   M  φm (θ) 2 −LN p(y|θ) ∝ , (8) (πσw (fm )) exp −L 2 σw (fm ) m=1 and when both the noise variances the source terms are unknown, the likelihood is expressed as [13, Eq. (3.11)] p(y|θ) ∝

M 

exp (−LN loge φm (θ)) ,

(9)

m=1

where φm (θ) = tr(Rm ) −

gH (fm , θ)Rm g(fm , θ) , gH (fm , θ)g(fm , θ) L 1 yl (fm )ylH (fm ). Rm = L

(10)

l=1

III.

BAYESIAN

SOURCE LOCALIZATION WITH UNCERTAIN G REEN ’ S FUNCTION

For source localization, the unknown parameter set θ represents the range and depth of the source. Uncertainty in the environment is generally handled by adding environmental parameters to θ at the cost of greatly increased computational

time [6]–[9]. In this paper no environmental parameters are added to θ. Environmental uncertainty is handled by considering its rough effect on the Green’s function as detailed in this section. For simplicity only monofrequency MFP is considered in what follows, although the results can be easily extended to the multifrequency problem. Also the dependence of the Green’s function g on θ is omitted for clarity. Suppose that the Green’s function of the assumed waveguide is a random vector with a prior density function p(g) that accounts for environmental uncertainty, the PPD of the localization parameters can be obtained through marginalization over p(g)  p(θ|y) ∝ p(y|g)p(g)dg p(θ). (11) A. Uncertain Green’s function model The prior p(g) aims at characterizing the a priori knowledge on the possible mismatch by its effect on the Green’s function of the waveguide. A simple and flexible way to do so is to consider that the Green’s function is a zero-mean Gaussian random vector characterized by its spatial covariance matrix Σg . The covariance matrix Σg determines how the true Green’s function of the waveguide departs from the assumed Green’s function when the true oceanic waveguide departs from the assumed oceanic waveguide. Various amount of information may be considered. The least informative model considers Σg = diag(σg2 , 1, · · · , 1) in a basis where σg2 is a variance in the direction of the assumed Green’s function. A high value of σg2 means that the true Green’s function of the waveguide does not statistically depart much from the assumed the Green’s function while a low value means that it statistically departs a lot. This simple model aims solely at characterizing the strength of the mismatch on the Green’s function. More informative models may consider expected directions and correlations leading to more complicated covariance matrices Σg . Furthermore, the Green’s function is assumed to be constant in time considering that the environment should not vary much for the L snapshots considered. B. Green’s function marginalization The Green’s function g is a N × 1 vector, hence marginalization of the Green’s function over the PPD in Eq. (11) involves an N -dimensional integral. Monte-Carlo integration of such a high dimensional integral is prohibitively costly in general. In this paper, we address this problem by reducing this high dimensional integral to a one-dimensional integral. Denote u the quantity u=

1 gH Rg , tr(R) gH g

(12)

where R is defined in Eq. (10). Since φ in Eq. (10) can be written as φ = tr(R)(1 − u), the likelihoods in Eqs. (8) and (9) can be expressed as functions of u instead of g. The multidimensional integral can therefore be replaced by a simple onedimensional integral to obtain the PPD over the localization parameters:  1 p(y|u)p(u)du p(θ), (13) p(θ|y) ∝ 0

as the random variable u satisfies 0 ≤ u ≤ 1 1 . The probability density p(u) is not known analytically but can be computed using results on quadratic form in Gaussian variables as the quantity u is a ratio of quadratic form in Gaussian variables. One may express the cumulative distribution function (CDF) FU (u) of u as:   1 gH Rg FU (u) = Pr < u , (14) tr(R) gH g    R (15) − uIN g < 0 . = Pr gH tr(R) 1/2

Define h such that g = Σg h, where Σg is the covariance matrix of g introduced in the previous subsection, this operation is a whitening transform and h ∼ CN (0, IN ). The CDF FU (u) can be written as 

(16) FU (u) = Pr hH Au h < 0 , where Au =

Σ1/2 g



R − uIN Σ1/2 g . tr(R)

(17)

Thus, FU (u) can also be seen as the CDF at zero of the quadratic form hH Au h. Let {λui , i = 1, .., P } be the positive eigenvalues of the matrix Au and assume that they are all different, the CDF can be computed as follows [14] FU (u) = 1 −

P  i=1

N 

−1 λN ui

.

(18)

(λui − λul )

l=i

of the positive eigenvalues of the matrix Obviously, none R are repeated, however we don’t know if this − uI N tr(R) is enough to ensure that all the positive eigenvalues of Au are different. On the other hand, this seems quite reasonable and has always been the case in our simulations. In any case, if some positive eigenvalues are not distinct a more general expression may be used [14]. The probability density p(u) is computed by the numerical differentiation of the CDF FU (u), and the integrand in Eq. (13) is evaluated at a finite set of points to perform the numerical integration. Integration on a regular grid with about 100 sampling points has been shown to give accurate results. IV.

R ESULTS

This Bayesian source localization method is applied to a mismatched localization problem in a shallow water waveguide. The oceanic environment is illustrated in Figure 1. It has an uncertain water column sound speed profile set by the sound speeds c1 (at surface), c2 (at 10 m depth), c3 (at 50 m depth) and c4 (at seabed, depth D), and the seabed consists of an uncertain bottom halfspace. The uncertain parameters bounds, the assumed waveguide and the true environment for the two test case considered are given in Table I. The range is unknown with a uniform prior distribution on the interval rs = [500, 10000] m and the depth is unknown with 1 The positive semi-definite matrices R and ggH obey the Cauchy-Schwarz inequality with the trace operator as the inner product which leads to a proof that u ≤ 1.

a uniform prior distribution on the interval zs = [2, 95] m. The receiver array is an N = 15 element vertical line array sampling the water column on a regular grid between z = 2 m and z = 86 m. The computation of the PPD is carried out for the likelihood with unknown variance in Eq. (9) on a numerical grid involving depth and range increments of 2 m and 50 m. Only the f = 300 Hz frequency is used to carry out the inversion. The number of snapshots considered is L = 10. The SNR is 3 dB for every snapshots and defined 2 2 m ,θ0 ) as SNR = |sl (fmN)|σg(f . At this SNR, the MFP source 2 w (fm ) localization works well when there is no mismatch however this is not a realistic assumption. Figure 2(a) and 3(a) show the PPD of range and depth when inference is made as in Section II without knowing that there is mismatch, respectively for test case 1 and test case 2. Figure 2(b) and 3(b) show the PPD of range and depth after marginalization over the least informative Green’s function model introduced in Section III-A, respectively for test case 1 and test case 2. The parameter σg2 = 5 was roughly evaluated from a quick sensitivity analysis of the acoustic field on the N = 15 element vertical line array for the environmental parameter bounds given in Table I.

 
  











  




 

   



Fig. 1: Schematic diagram of the ocean environment for the mismatched source localization examples.

Parameters rs (km) zs (m) c1 (m/s) c2 (m/s) c3 (m/s) c4 (m/s) D(m) cb1 (m/s) ρb1 (g/cm3 ) αb1 (dB/λ)

Bounds [0.5, 10] [2, 95] [1515, 1525] [1514, 1522] [1510, 1516] [1508, 1512] [98, 102] [1550, 1750] [1.2, 2.2] [0.05, 0.5]

Assumed waveguide − − 1520 1518 1513 1510 100 1650 1.7 0.275

Test case 1 1 49 1523.5 1521.5 1512.5 1511.5 100.1 1770.7 2.06 0.052

Test case 2 6 30 1521.4 1519.8 1511.9 1511.5 98.6 1630.5 1.84 0.341

TABLE I: Simulation parameters for the mismatched localization examples. The two cases depict typical results. In Figure 2 both PPDs suggests almost certainly (i.e. with a probability close to one) the right source position. In Figure 3, while the PPD in Figure 3(a) suggests almost certainly a wrong position for the source, the PPD in Figure 3(b) takes into account the uncertainty on the Green’s function and suggests other positions including the true one. These two examples shows the ability of the method to give more reliable PPDs and thus to avoid making incorrect assumptions on the source position. When the localization is uncertain because of the mismatch, the PPD is also uncertain and does not suggest a

wrong position with a strong probability. Furthermore, more informative uncertain Green’s function covariance matrices Σg than the simple matrix used here may also be used to further enhance robustness toward mismatch. V.

C ONCLUSION

In this paper we presented an original Bayesian localization method to handle uncertainty in the waveguide environment. No environmental parameters are added to the localization problem, instead the Green’s function of the waveguide is considered as an uncertain Gaussian random variable whose probability density accounts for environmental uncertainty. The uncertain distribution over range and depth is then obtained through the integration of the PPD over the Green’s function probability density. An efficient integration technique makes the whole localization process computationally efficient. A very simple model of uncertain Green’s function where only the strength of the mismatch effect on the Green’s function is known a priori has been considered. This model does not change the location of the maximums of the PPD but avoid suggesting a wrong position with a strong probability which results generally in more reliable PPDs. More informative models of uncertain Green’s function may be considered. For instance, physical insight of the mismatch effect on the Green’s function can be obtained through the modal representation of the wavefield. Using the fact that different groups of modes are affected by different kind of mismatch may allow to construct an uncertain Green’s function model that capture some important features of mismatch. Adding some information may enable to suppress wrong positions from the PPD and further enhance robustness.

(a) no mismatch integration

20

0.8

40

0.6

60

0.4

depth (m)

0.2

80

20

0.8

40

0.6

60

0.4

0

0.2 0.4 0.6 0.8

0.2

80 0.2 0.4 0.6 0.8

depth (m)

(a) no mismatch integration

0.8 0.6 0.4 0.2

0.8 0.6 0.4 0.2 2

4 6 range (km)

8

10

2

(b) mismatch integration

4 6 range (km)

8

10

(b) mismatch integration 0.08 0.6

40

0.4

60

20 depth (m)

depth (m)

20

0.06 40 0.04 60

0.2

0.02

80 0.02 0.06 0.1

0.2 0.4 0.6

80

0.1 0.08 0.06 0.04 0.02

0.6 0.4 0.2 2

4 6 range (km)

8

10

Fig. 2: Test case 1 : joint marginal PPD over source range and depth together with the associated one dimensional marginals, (a) no integration over the uncertain Green’s function model, (b) integration over the uncertain Green’s function model. The circle indicates the MAP estimate, the triangle and the dashed lines indicate the true position.

2

4 6 range (km)

8

10

Fig. 3: Test case 2 : joint marginal PPD over source range and depth together with the associated one dimensional marginals, (a) no integration over the uncertain Green’s function model, (b) integration over the uncertain Green’s function model. The circle indicates the MAP estimate, the triangle and the dashed lines indicate the true position.

ACKNOWLEDGMENT This work was funded by the French Government Defense procurement agency (Direction G´en´erale de l’Armement). The authors would like to thank the GdR ISIS (Groupement de Recherche en Information Signal Image et viSion) for its support. R EFERENCES [1]

[2]

[3]

[4]

[5]

[6]

[7] [8]

[9]

[10]

[11]

[12]

[13]

[14]

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. D.R. Del Balzo, C. Feuillade, and M.M. Rowe, “Effects of water-depth mismatch on matched-field localization in shallow water,” The Journal of the Acoustical Society of America, vol. 83, no. 6, pp. 2180–2185, 1988. C. Feuillade, D.R. Del Balzo, and Mary M. Rowe, “Environmental mismatch in shallow-water matched-field processing: Geoacoustic parameter variability,” The Journal of the Acoustical Society of America, vol. 85, no. 6, pp. 2354–2364, 1989. E.C. Shang and Y.Y. Wang, “Environmental mismatching effects on source localization processing in mode space,” The Journal of the Acoustical Society of America, vol. 89, no. 5, pp. 2285–2290, 1991. 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. A.M. Richardson and L.W. Nolte, “A posteriori probability source localization in an uncertain sound speed, deep ocean environment,” The Journal of the Acoustical Society of America, vol. 89, no. 5, pp. 2280– 2284, 1991. S.E. Dosso, “Environmental uncertainty in ocean acoustic source localization,” Inverse Problems, vol. 19, no. 2, pp. 419–431, 2003. S.E. Dosso and M.J. Wilmut, “Bayesian focalization: Quantifying source localization with environmental uncertainty,” The Journal of the Acoustical Society of America, vol. 121, no. 5, pp. 2567–2574, 2007. S.E. Dosso and M.J. Wilmut, “Bayesian multiple-source localization in an uncertain ocean environment,” The Journal of the Acoustical Society of America, vol. 129, no. 6, pp. 3577–3589, 2011. P. Gerstoft and C.F. Mecklenbrauker, “Ocean acoustic inversion with estimation of a posteriori probability distributions,” The Journal of the Acoustical Society of America, vol. 104, no. 2, pp. 808–819, 1998. C.F. Huang, P. Gerstoft, and W.S. Hodgkiss, “Uncertainty analysis in matched-field geoacoustic inversions,” The Journal of the Acoustical Society of America, vol. 119, no. 1, pp. 197–207, 2006. Y. Le Gall, F-X. Socheleau, and J. Bonnel, “Matched-field processing performance under the stochastic and deterministic signal models,” Signal Processing, IEEE Transactions on, vol. 62, no. 22, pp. 5825– 5838, 2014. 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. 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.