Underwater Source Localization via Spectral Element Acoustic Field Estimation

Underwater source localization (USL) from a passive array of acoustic sensors is a challenging problem, especially in complex environments characterized by multipath and reverberation effects, irregular seabed geometry, and low signal-to-noise ratio (SNR). This article proposes a recursive Bayesian approach that propagates a spectral-element approximation of the wave equation to model the discretized space–time dynamics of the acoustic field conditioned on the position of the source, and sequentially estimates the field and the position of the radiating source from the acoustic measurements. We pursue a multiple-model approach where each model assumes either the source absence or its presence within a specific spectral element. To handle the high dimension of the large-scale field estimation problem and reduce the computational complexity, the multiple-model filter is implemented using ensemble Kalman filters (EnKFs). Finally, the effectiveness of the proposed multiple-model spectral-element ensemble Kalman filter is demonstrated through simulation experiments in underwater acoustic environments with regular and irregular seabed geometry and via comparison with the standard matched-field processing (MFP) method.


Underwater Source Localization via Spectral
Element Acoustic Field Estimation Graziano A. Manduzio , Nicola Forti , Roberto Sabatini , Giorgio Battistelli , and Luigi Chisci , Senior Member, IEEE Abstract-Underwater source localization (USL) from a passive array of acoustic sensors is a challenging problem, especially in complex environments characterized by multipath and reverberation effects, irregular seabed geometry, and low signal-to-noise ratio (SNR).This article proposes a recursive Bayesian approach that propagates a spectral-element approximation of the wave equation to model the discretized space-time dynamics of the acoustic field conditioned on the position of the source, and sequentially estimates the field and the position of the radiating source from the acoustic measurements.We pursue a multiple-model approach where each model assumes either the source absence or its presence within a specific spectral element.To handle the high dimension of the large-scale field estimation problem and reduce the computational complexity, the multiple-model filter is implemented using ensemble Kalman filters (EnKFs).Finally, the effectiveness of the proposed multiple-model spectral-element ensemble Kalman filter is demonstrated through simulation experiments in underwater acoustic environments with regular and irregular seabed geometry and via comparison with the standard matched-field processing (MFP) method.

NOMENCLATURE R +
Set of nonnegative real numbers.N + Set of nonnegative integers.

G(•; m, P)
Gaussian probability density with mean m and covariance P. a ∼ G(•; m, P) a is Gaussian with mean m and covariance P.
[a; b] a b .

I. INTRODUCTION
U NDERWATER source localization (USL) using a passive array of acoustic sensors aims at detecting a radiating source and estimating its position in space by analyzing the acoustic field measured by an array of hydrophones.Due to the special characteristics and complexity of the underwater environment, USL using a passive array of acoustic sensors is a challenging task, attracting great interest within both the control and signal processing communities, especially in low signal-to-noise ratio (SNR) scenarios where the source can be difficult to detect [1].When the source is in the far-field of the sensor array, i.e., far enough for the plane-wave approximation to be valid, conventional passive sonar systems based on array signal processing can only estimate the direction of arrival (DOA) of an acoustic signal impinging on the array via plane-wave beamforming [1], [2].DOA is determined using the time differences in the sound arrival among the spatially separated hydrophones of the array.This is the main difference from active sonar systems where the absolute time of signal transmission from a projector is known, also allowing estimation of the range from the sonar system to the source of interest.In addition to conventional beamforming, other DOA estimation methods have been developed and applied to USL [3], including subspace-based approaches such as multiple signal classification (MUSIC), estimation of signal parameters via rotational invariance technique (ESPRIT), and their modifications [4].
In low SNR environments, the major practical challenge is to extract the DOA information from acoustic signals that are distorted by different types of noise (e.g., ambient, measurement, and thermal noise).In fact, classical methods based on DOA estimation can only exploit thresholded measurements.Then, for low SNRs, the received signals are more likely to be undetected using conventional DOA methods since the information contained in the measurements may be irreversibly discarded after the thresholding process [5], [6], [7].Moreover, in realistic underwater environments, DOA estimation methods can be highly affected by multipath and reverberation effects, inhomogeneous domains, and irregular seabed geometries, resulting in inaccurate localization.In such scenarios, acoustic propagation and its effects can be exploited to enable passive localization.To this end, the major developments in underwater passive localization since the early 1980s have focused on the inclusion of acoustic propagation modeling into signal processing algorithms [8].When the acoustic propagation model is accounted for in the passive sonar array signal processing, we refer to the procedure as matched-field processing (MFP) [9], [10], [11], [12].MFP is a well-established generalization of plane-wave beamforming for USL where the steering or replica vector is derived from a solution of the wave equation in the frequency domain (i.e., the Helmholtz equation) for a point source, especially common for shallow water applications where multipath propagation can yield useful information to infer the source range and depth.MFP is based on matching (correlating) the acoustic pressure field measured at the hydrophone array with modeled replica fields computed for a given acoustic waveguide environment via a numerical propagation model over a set of test source positions, typically distributed over a uniform search grid covering the given range-depth domain.The maximum in the cross correlation or ambiguity surface gives an estimate of the position of the underwater source, usually taken to be the grid point associated with the highest match.Source localization using common MFP is effective when the source is already detected and stationary, its frequency is known, and there is full knowledge of the propagation model.However, it can be highly degraded or precluded in the case of inaccurate correspondence (mismatch) between the propagation model and the real oceanic waveguide, i.e., for noisy or uncertain source and environmental parameters [13], [14], [15].To improve USL performance, MFP approaches combining information from multiple arrays have recently been studied for shallow-water scenarios [16], [17].It is demonstrated that spatially coherent processing of multiple arrays can yield significant improvement in localization performance over incoherent processing.However, it can also be more susceptible to model mismatch than incoherent processing.
The main challenge of USL is to devise estimators which are less sensitive to noise level and environmental mismatch.Conventional MFP methods rely on numerically efficient propagation models such as normal modes that are not capable of treating inhomogeneous environments with complex geometry characterizing scattering and reverberation problems.More general numerical methods based on direct discretization of the governing full wave equation include the finite difference and finite element method (FEM).In particular, FEM is extremely versatile in terms of geometries and environment properties that can be treated, and boundary conditions and point sources can be directly incorporated in its formulation.However, with increasing frequencies, the FEM mesh requires a growing number of finite elements, thus resulting in increased computational time.The spectral element method (SEM) [18], [19], [20] is a particular kind of FEM that uses high-degree piecewise polynomials as basis functions within quadrangular (in 2-D) or hexahedral (in 3-D) elements.The SEM provides high accuracy with fewer degrees of freedom in comparison to the standard FEM.Moreover, it uses nonuniformly distributed nodes, the so-called Gauss-Lobatto-Legendre (GLL) points, and, when used in conjunction with the GLL quadrature rule, results in a diagonal mass matrix, allowing for reduced complexity, computational time, and memory.
This work exploits recent advancements in numerical simulation of partial differential equation (PDE) systems [18], [19], [20], [21], [22] to approximate the wave equation into a finite-dimensional space-time model of the underwater acoustic field.In addition, this article builds on large-scale field estimation of discretized PDE systems [23], [24] and previous work on source identifiability and estimation in such systems [25], [26].Further related work focused on USL for shallow-water environments and high-frequency signals using a multiray propagation model [27], decentralized detection in underwater sensor networks [28], decentralized USL via generalized likelihood ratio test [29], a self-supervised learning architecture that exploits joint time-frequency processing for USL [30], and acoustic source localization and tracking using a cluster of mobile agents [31].
This article presents a novel multiple-model spectralelement ensemble Kalman filter (MM-SE-EnKF) algorithm for underwater source detection and localization.The proposed method incorporates an SEM-based model of underwater acoustic propagation in the time domain to sequentially estimate the acoustic field directly from acoustic pressure measurements.This allows USL in complex environments characterized by inhomogeneous properties and/or irregular seabed geometries.The ensemble Kalman filter (EnKF) [32] is adopted for computationally efficient acoustic field estimation.Source localization is performed using a multiple-model approach that runs in parallel a bank of field estimators, each conditioned to the source being placed in a given element of the discretized computational domain, plus a null hypothesis accounting for the possible absence of the source.This article improves previous work on USL using FEM-based acoustic field estimation [33] by developing a more computationally efficient SEM-based formulation of the acoustic propagation model, a faster algorithm for multiple-model estimation of the element containing the source position, and by including results in terms of USL in the case of both regular and irregular seabed geometry.In addition, the former scenario includes a comparison of results against standard MFP.Summing up, the main contributions are as follows: 1) we provide a Bayesian formulation for the problem of joint pressure field estimation and source detection/localization directly exploiting the acoustic measurements and show how to directly integrate the SEM within the Bayesian formulation; 2) we develop a novel multiple-model filtering algorithm based on a bank of SE EnKFs for joint pressure field estimation and source detection/localization capable of dealing with arbitrary geometries and inhomogeneous properties; and 3) to improve the computational efficiency, we develop a fast algorithm in which, at each time, only the SE EnKF corresponding to the most likely mode is activated.
The rest of this article is organized as follows.Section II describes the formulation of the considered USL problem.Section III presents an SEM-based approximation and time integration of the generalized wave equation modeling underwater acoustic propagation.In Section IV, a novel MM-SE-EnKF algorithm is proposed to address USL.Section V shows the simulation results, while Section VI concludes this article.

A. Notation
Vectors and matrices are denoted by lower case and, respectively, upper case letters in bold, while scalars are denoted by normal lower case letters.For convenience of the reader, notation and symbols used throughout the article are listed in Nomenclature.

II. PROBLEM FORMULATION
Consider the 2-D infinite oceanic waveguide depicted in Fig. 1, bounded from above by a flat free surface U and from below by a possibly irregular seabed D .Define a Cartesian coordinate system Oξ ζ , with origin O on U and vertical axis ζ oriented toward the seafloor, and denote by p p p ∈ the position vector and by t ∈ R + the time variable.Let ζ = b(ξ ) be the function describing the variable bathymetry, and let α(p p p), c(p p p), and ρ(p p p) be the assumed-known spacedependent damping coefficient, ambient speed of sound, and water density, respectively.Suppose that a perturbation of pressure x(p p p, t; p p p s (t)) is generated in the water column by a point source f located at p p p s (t).The proposed algorithm aims at detecting the presence of the sound emitter and estimating its position p p p s (t) and the induced acoustic field x(p p p, t; p p p s (t)), given discrete time pressure signals y i (t), recorded by N h sensors (hydrophones) at known locations p p p h i , i = 1, . . ., N h .To this end, the propagation of the acoustic perturbation x(p p p, t; p p p s (t)) is assumed governed by the generalized wave equation with initial conditions and boundary constraints The forcing term is null in the absence of source and modeled as f (p p p, t; p p p s (t)) = s(t) δ(p p p − p p p s (t)) otherwise, where s(t) is the source temporal waveform and δ denotes the Dirac delta.In this work, the function s(t) is considered known, while the source position p p p s (t) and the acoustic perturbation x(p p p, t; p p p s (t)) are to be estimated.More specifically, the temporal envelope s(t) is a sinusoidal waveform with known frequency f s This assumption is reasonable for narrowband sources typically encountered in underwater applications.
The above-stated dynamic estimation problem is clearly infinite-dimensional.To make it numerically tractable, a finitedimensional approximation of the solution x(p p p, t; p p p s (t)), based on the spectral-element method, is introduced in Section III.

III. SPECTRAL-ELEMENT DISCRETIZATION AND TIME INTEGRATION A. Computational Domain and Integral Formulation
Problem (1) is solved numerically through the SEM proposed in [18], [19], [20], and [21].The infinite domain is first truncated along the horizontal axis, as illustrated in Fig. 1.The resulting computational domain d is delimited by the piecewise closed line U ∪ R ∪ D ∪ L and spans the physical region of interest Let n n n denote the outward pointing unit normal on the boundary of d .To ensure that the acoustic waves leave the domain without significant spurious reflections, the following radiation condition [1,Ch. 7] is applied to the left and right boundaries, As highlighted by Jensen et al. [1,Ch. 7], this constraint is only exact for a plane wave impinging normally onto a plane boundary.Nonetheless, it significantly reduces the amplitude of the numerical reflections even for moderate angles of incidence.
As a standard finite-element method, SEM is based on an integral formulation of problem (1).Multiplying (1a) by a generic space-dependent test function ψ(p p p), integrating over d , and enforcing the boundary constraints (1c) and the radiation condition (4) yield

B. Computational Mesh
The computational frame d (see Fig. 2) is partitioned into a set of N e nonoverlapping quadrilateral elements e , e = 1, . . ., N e , such that These elements are constructed with curvilinear boundaries to adapt to the irregular bottom.To this end, the following change in coordinates is first introduced: In (7), the parameter D represents a characteristic reference depth.As a result, the complex waveguide The Each element e in the physical plane Oξ ζ is thus mapped to the reference square Equations ( 7) and ( 9) define a local vector-valued invertible function T e : → e that maps the pair (u, v) ∈ to a point p p p ∈ e , i.e., p p p = T e (u, v).
For the sake of clarity, it is worth pointing out that the change in variables here used to transform the irregular physical region into a rectangular domain allows to avoid using mesh generation software and dramatically simplifies the structure of the numerical solver.Nevertheless, such a change in variables is not strictly necessary, and unstructured meshes could be used to handle irregular seabeds.

C. Polynomial Approximation in Spectral Elements
In a traditional finite-element method, low-degree polynomials are used as basis functions for computing the unknown field x(p p p, t).In the SEM, the pressure x(p p p, t) is approximated by a higher degree Lagrangian interpolant.To this purpose, a set of g nodes, called GLL points, is first defined along each direction u, v, within the reference element (see Fig. 2).On a direction w = u, v, the GLL nodes w i , i = 1, . . ., g, are defined as the roots of the equation where P ′ g represents the first derivative of the Legendre polynomial of degree g − 1.Typically, g is between 4 and 10.The GLL points are computed numerically and, as shown in Fig. 2, they form a nonuniform grid of g × g nodes in that is mapped, via T e , to a nonuniform grid e in e e = T e u i u , v i Within a generic element e , the restriction x e (p p p, t) of the pressure field x(p p p, t) to e is then approximated as where x e i u i v (t) is the time-dependent value of x e (p p p, t) at the GLL point p p p = T e (u i u , v i v ) and L i is the gth degree Lagrange interpolating polynomial associated with node i, i.e., By definition, the function L i (w) satisfies the property Boundary nodes are shared by neighbor elements, and the union of all the local grids contains grid points p p p j , j = 1, . . ., N n , which span the whole computational domain d .The triple of indices (e, i u , i v ), with i u , i v = 1, . . ., g, identifies a unique node j of the global grid .As a result, the pressure field x(p p p, t) can be expressed as a linear combination of N n spatially varying basis functions j (p p p), j = 1, . . ., N n , each of which is associated with a grid point j of .Accordingly, ( 13) is rewritten as where x j (t) is the time-dependent value of the pressure field x at node j, and (p p p), x x x(t) ∈ R N n are the column vectors The basis function j (p p p) does not vanish only at points p p p belonging to elements that contain the node j (more than one in the case of boundary nodes).More specifically, they are defined as

D. Time Integration
To compute N n unknowns x j , j = 1, . . ., N n , expansion ( 18) is first introduced in the weak form (5).Then, by replacing the generic test function ψ with the basis functions, the following system of N n second-order-in-time ordinary differential equations is obtained: In (21), M, D, K ∈ R N n ×N n are, respectively, the mass, damping, stiffness matrices, and Integrals in (22) are evaluated through the GLL quadrature method [34].In conjunction with the use of Lagrangian interpolants based on GLL nodes, this choice leads by construction to a diagonal mass matrix, which results into a drastic reduction of both the complexity of the numerical method and the computational time.No costly matrix inversion algorithm is indeed needed to compute the solution in time.This is one of the main differences between the SEM used in this work and more classical FEMs.Time discretization of ( 21) is achieved via an explicit Newmark scheme [18].Let t be the time integration step.Then, the approximations x x x k+1 and ẋ x x k+1 at instant t k+1 = (k +1) t, k ∈ N + , of x x x and its time derivative ẋ x x are computed as with the initial conditions Using (21) for ẍ x x k and ẍ x x k+1 finally yields the following discrete-time linear system: where . Equation ( 25) are solved sequentially: the discretized acoustic pressure x x x k+1 is first computed via (25a) and is then inserted in (25b) to calculate its derivative ẋ x x k+1 .To conclude, at each time instant t k = k t, k ∈ N + , the pressure signals x x x h k ≜ col{x(p p p h i , t k )} N h i=1 measured by the hydrophones at locations p p p h i , i = 1, . . ., N h , are given by where the elements of the N h × N n matrix C C C turn out to be

IV. SOURCE DETECTION AND LOCALIZATION A. MM Approach
Relying on the SEM of the previous section, the key idea for source detection and localization pursued in this work is to consider multiple hypotheses, one corresponding to the absence of the source and the others corresponding to a source located in a generic element of the SE mesh.In this way, recalling that the source presence/location only affects the forcing term f f f k = f f f (t k ) in (25), it is possible to associate an appropriate model to each hypothesis and adopt a multiple-model filtering approach [35,Ch. 11] to ascertain which model/hypothesis is more likely with the available sensor measurements and accordingly decide whether the source is present (detection) and, in such a case, estimate its position (localization).
To be more precise, let N e ⊂ {1, . . ., N n } denote the subset of nodes of element e ∈ {1, . . ., N e } and let e = 0 refer to the null model/hyphotesis accounting for the absence of the source in the monitored area.Then, for model/hypothesis e ∈ {1, . . ., N e }, the forcing term f f f k = f f f e k turns out to be where is the restriction of (p p p) to element e for e = 1, . . ., N e (source located in e ), while for e = 0 (nonexisting source) f f f e k is null.Accordingly, the model matched to hypothesis e ∈ {0, 1, . . ., N e } consists of the following state equations for acoustic pressure x x x k and its derivative ẋ x x k : Furthermore, whenever e ̸ = 0 (existing source located in e ), the dynamical model for the source position p p p s k is p p p s k+1 = p p p s k + w w w 3,k .
In (30), w w w 1,k , w w w 2,k , w w w 3,k are process disturbances that account for uncertainty on pressure, pressure derivative, and source position, respectively.Such process disturbances can be used to account for different sources of uncertainty including parameter uncertainties, time and space discretization errors, etc.Note that in (30c), the position of the source is assumed to vary slowly with time and, hence, to follow a discrete-time random walk.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Then, each of the above models involves the augmented state In view of ( 26), the measurement equation turns out to be independent of e equal to where y k y k y k is the vector of all the measurements collected from the sensor array and v v v k accounts for measurement noise.
The multiple-model filter for source detection and localization runs a bank of N e + 1 spectral-element nonlinear filters associated with the above discussed hypotheses (e = 0 for source absence and e = 1, . . ., N e for its presence in spectral element e).To this end, each filter e must propagate in timepredicted (filtered) estimates ẑ ẑ ẑe k|k−1 (ẑ ẑ ẑe k|k ) of the augmented state z z z k defined in (31) according to hypothesis e, along with the probability that the current mode e k is equal to hypothesis e. Transitions among hypotheses are modeled, as usual, by means of a homogeneous Markov chain with constant transition probabilities Remark 1: The choice of the above transition probabilities could be related to the ratio r = T /τ min of the filter sampling interval T to the minimum sojourn time τ min = ( L/v max ) of the source within an SE, v max denoting the maximum source speed.In fact such ratio, satisfying r < 1 for typical values of ( L , v max , T ) can be regarded as a probability that the source leaves an element e ′ to enter a neighboring element e; hence, r should be equally divided among all neighbors e of e ′ to define π ee ′ and the residual probability be assigned to permanence in e ′ , i.e. π e ′ e ′ = 1 − r .In particular, if the source is known to be motionless (infinite sojourn time), the suggested choice is π ee ′ = 1 for e = e ′ or π ee ′ = 0 otherwise.
Remark 2: Tracking of a quickly moving source could be accomplished in the multiple model (MM) framework by: 1) adopting suitable higher order kinematic models in place of the random-walk model (30c); 2) suitably specifying the transition probabilities in (34); and 3) resorting to an interacting MM (IMM) approach [35].It is worth to point out, however, that the standard IMM filter involves mixing of the MM states and needs therefore to be tailored to the source tracking case where the source states of the various models, constrained to belong to different spectral elements, cannot clearly be mixed.The extension of this approach to the case of a moving source will be the objective of future work.

B. Ensemble Kalman Filter
The critical point of this approach is that the state z z z k in (31) tends to be of very high dimension since it is in the order of 2N n .In fact, from (17) it turns out that where N e = N ξ N ζ is the number of elements of the mesh and g is the number of nodes along each side of the SE, and hence is typically very large especially when the size of the computational domain d increases and/or the mesh resolution L gets finer and/or g increases.In this respect, the cubic computational complexity and quadratic memory complexity (with respect to the state dimension n) of any Kalman-like filter propagating in time the state covariance matrix would become infeasible.To this end, the so-called ensemble Kalman filter (EnKF) has been devised [32], [36] to efficiently cope with problems (e.g., in meteorology, geoscience, oceanography) that involve a large number of state variables typically arising from spatial discretizations of PDEs.The idea of EnKF is to replace the n × n state covariance matrix (where n is the dimension of the state vector) with an ensemble consisting of randomly sampled q ≪ n state vectors and to compute the required statistics for estimation (mean, cross-covariance between state and measurement, measurement covariance) as sample averages over the ensemble.This allows to drop the computational complexity from O(n 3 ) to O(n 2 q) and the memory complexity from O(n 2 ) to O(nq) [38,Sec. 4.2], with significant savings provided that the ensemble size q is much smaller than the state dimension n.
The pseudocode of a multiple-model spectral-element EnKF (MM-SE-EnKF) for acoustic source detection and localization is provided in Algorithm 1.Such algorithm consists of the following steps.1) Prediction: First, for each element hypothesis e and member of the ensemble i, the acoustic pressure field is predicted by the space-time discretized acoustic wave propagation model (30a) and (30b), and, if e ̸ = 0, the source position is also predicted by the random-walk model (30c); then, mode probabilities are predicted by the Markov chain transition model (34).2) Correction: The augmented states, for any e and i, and mode probabilities, for any e, are corrected on the basis of the available measurements from hydrophones.

3) Source Detection and Localization:
The hypothesis e with highest probability is selected and, if e ̸ = 0 (detected source), the source is localized according to the corrected source position relative to such hypothesis.
MM-SE-EnKF runs a bank of N e + 1 SE-EnKFs each associated with a possible hypothesis (mode) e ∈ {0, 1, . . ., N e } and propagating in time predicted and filtered augmented state estimates ẑ z z e k|k−1 and, respectively, ẑ z z e k|k and the mode probability µ e k .Each filter can run independently of the others performing, at each time k, the prediction step (lines 4-25 of the pseudocode) followed by the correction step (lines 26-47) and by the source detection and localization step (lines 48-58).Note that each mode-matched SE-EnKF uses an ensemble of size q to compute the required joint statistics for the augmented state z z z and measurement y y y (lines 22, 33, and 34).
Note also that in each EnKF, the model uncertainties are incorporated through the use of ensemble simulations to generate multiple realizations of the state vector in the prediction step to approximate the underlying probability distribution.Such model uncertainties are quantified in terms of the variances σ 2 w 1 , σ 2 w 2 , σ 2 w 3 of the process disturbances Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.w w w 1,k , w w w 2,k , w w w 3,k which determine the spread of the ensemble simulations.

C. Constraints in Prediction and Correction Steps
Recalling that filter e, for e > 0, is associated with the hypothesis that the source is located in e , then the source position samples p p p s,i,e k|k−1 and p p p s,i,e k|k must be enforced to be in e .
This can be accomplished in several different ways.For the sake of computational simplicity, in this article we adopt the following procedure.
1) In the prediction step, p p p s,i,e k|k−1 is resampled until it falls inside e (see Fig. 3).2) In the correction step, any sample p p p s,i,e k|k falling outside e is projected to the closest point on the boundary of e (see Fig. 4).The reason for this different handling of constraint violation in prediction and correction steps can be explained as follows.For source location prediction p p p s,i,e k|k−1 = p p p s,i,e k−1|k−1 + w w w i,e 3,k−1 starting from p p p s,i,e k−1|k−1 ∈ e and provided that σ w 3 is sufficiently small compared with the element size L, it is quite unlikely that p p p s,i,e k|k−1 ̸ ∈ e for several consecutive draws of w w w i,e 3,k−1 ∼ G(•; 0 0 0, σ 2 w 3 I I I ); hence, in this case, resampling represents a simple and effective solution.On the other hand, correction (line 37) tends to push p p p s,i,e k|k outside e for elements e that are far away from the true source position; in such a case, resampling would be highly inefficient while projection onto the element boundary is a simple and efficient solution.

D. Computational Complexity of MM-SE-EnKF
It is first worth pointing out that the computational burden of MM-SE-EnKF is mostly due to the N e + 1 mode-matched filters, as the remaining part (lines 48-56) involves a negligible amount of computation.Moreover, such filters can clearly run in parallel so that they can be assigned to N p processors roughly reducing the overall computational load of a factor N p .
Hence, hereafter the focus is on the computational complexity of the single SE-EnKF (lines 6-46).Recalling that ≥ N e and that the ensemble size q must be chosen such that N h ≤ q ≪ N n for the sake of computational reduction as well as to guarantee that the N h × N h matrix P P P e y,k computed in line 33 of Algorithm 1 be invertible, the most burdensome task of SE-EnKF is for pressure and pressure derivative prediction which takes in the order of N 2 n q operations.Conversely, the correction step has O(N n q N h ) computational complexity required for the computation of matrix P P P e zy,k (line 34).Summing up, the overall computational complexity of MM-SE-EnKF turns out to be O(N 2  n N e q).Assigning the various SE-EnKFs to O(N e ) parallel processors, it is clearly possible to reduce the processing time for a single recursion of MM-SE-EnKF to O(N 2 n q).

E. Fast MM-SE-EnKF
To further reduce the computational load and thus allow working at faster sampling rates, a computationally cheaper version of MM-SE-EnKF referred to as fast MM-SE-EnKF (FMM-SE-EnKF) has been devised (see the pseudocode of Algorithm 2) to deal with a motionless acoustic source possibly located in the surveillance area d .
FMM-SE-EnKF switches between two operating phases.1) An initial source detection phase that aims to detect the presence of the source and single out the element e where it is located.2) A source localization phase that aims to accurately localize the source within e .The source detection phase must clearly take into account all possible hypotheses e ∈ {0, 1, . . ., N e } but, for the sake of computational efficiency, each mode-matched filter only performs prediction (but not correction with the available sensor measurements) of the pressure field estimate (lines 5 and 6 of Algorithm 2) and exploits the measurements only to update the hypothesis probabilities µ e k (lines 7-9).This allows to drop the computational complexity from O(N 2 n N e q) to O(N 2 n N e ) as no ensemble is required to compute the sample statistics for augmented state correction.In fact, prediction (25a) and (25b) involves multiplication of off-line computed N n × N n matrices by N n × 1 vectors for N e + 1 models.Conversely, the subsequent source localization phase only needs to run a single SE-EnKF for the selected hypothesis êk resulting from the source detection phase; in this phase, the computational complexity is just O(N 2 n q), the complexity of a single EnKF with state dimension 2N n + 2 = O(N n ) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

33:
P P P e y,k = (q − 1) −1 E E E e y,k (E E E e y,k ) T + R R R

34:
P P P e zy,k = (q − 1) −1 E E E e z,k (E E E e y,k ) T

35:
L L L e k = P P P e zy,k (P P P e y,k ) −1

36:
for i = 1, 2, . . ., q do 37: 38: if e ̸ = 0 then only for mode êk , run the SE-EnKF as described 25: in Algorithm 1 26: end if and ensemble size q.The switching the initial detection phase to the localization phase can be performed either after a prescribed setup time or whenever the selected hypothesis êk has remained unchanged for a sufficiently high number k s of consecutive recursions in the detection phase (see Algorithm 2 where the counter k is initialized to 1).Note that Algorithm 2, wherein the selected hypothesis is held fixed after the switching, is appropriate for a static (i.e., motionless and always  emitting) source but should clearly be reconceived for the dynamic source case which, however, is outside the scope of this article and left for possible future work.

V. SIMULATION RESULTS
We considered two scenarios to characterize the capabilities of the FMM-SE-EnKF algorithm.The first configuration, shown in Fig. 5, is used to assess the filter performance.The FMM-SE-EnKF algorithm is additionally compared with the minimum variance, distortionless filter MFP (MVDF-MFP) algorithm proposed by [10].It is worth emphasizing that current implementations of MFP are based on normal modes [1], [37] and are, by construction, applicable only in depth-dependent media with flat horizontal boundaries.On the contrary, the present algorithm can be used in more general space-time-varying underwater environments involving complex boundaries.The second scenario, shown in Fig. 6, is therefore designed to show the FMM-SE-EnKF algorithm's performance in the case of a variable seabed.

A. Scenario With Regular Seabed
This first scenario, inspired by a benchmark used in [39], involves a motionless source emitting a signal of frequency f s = 100 Hz for a total simulation time of 10 s and time integration step t = 0.23 ms.The source is located within a shallow water environment with an isospeed water column (c = 1500 m/s, ρ = 1000 Kg/m 3 , α = 4 • 10 −6 s/m 2 ) of 90 m depth and 1500 m length located at ξ = 1212.5m and ζ = 57.5 m.Fig. 5 shows the source position p p p s and the hydrophone positions p p p h i , for i = 1, . . ., N h , and provides a comparison between mesh elements e,l used by the ground-truth simulator with side length l = L g , and by the filter with side length l = L f within the domain d .We set the mesh size ratio (MSR), defined as the ratio MSR ≜ L f / L g between the side length of the mesh element used by the filter and the one used by the groundtruth simulator, to MSR = 0.5.Mesh parameters of the ground-truth simulator and the filter are reported in Table I.Filter parameters have been set as reported in Table II.The filter collects measurements from a vertical array (VA) of N h = 20 hydrophones located 1 Km away from the source in an environment with an SNR of 10 dB.Given the standard deviation of simulated measurement noise σ y , the SNR can be defined as where: is the acoustic pressure on the hth sensor at time k; and N is the total number of time integration steps.In this scenario, the source starts to emit at time t = 0 s while the data assimilation process is activated when the array starts to collect data, i.e., at time t 0 = 0.8 s when the acoustic wave has already reached the VA of hydrophones.Performance of the filter, in terms of detection and localization, is shown in Fig. 7.It can be seen how no detection is unveiled until t * = 1.516 s, when the array measurements begin to correctly match the predicted measurements ŷ y y e k generated from the model associated with e s .The mode in operation before t * is êk = 0, i.e., the nosource mode, while êk = e s , k ≥ t * , i.e., the true mode in operation e s is correctly detected starting from time t * .Once the source mode is estimated, the localization process within the source element is activated.Furthermore, in Fig. 7 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the distance between the estimated position p p p s k|k and the true source position p p p s is plotted versus time, compared with the distance obtained with the MVDF-MFP for the same dataset.It can be seen how the proposed algorithm yields better localization accuracy than the MVDF-MFP over the whole simulation time.In Fig. 8, we see that when the local EnKF is activated and the estimation is performed within the selected element e s , the initial distance corresponding to the center of the element becomes even smaller and the source localization performance further improves.Indeed, the position samples p p p s,i,e k|k (displayed in red) initially distributed all over the source element e s , and their average p p p s k|k (displayed in blue) gets closer to the true source position p p p s (displayed in magenta) with a final offset of about 3 m.Fig. 9 provides several plots comparing the estimated and ground-truth acoustic fields over the domain d at different time instants.Such fields can be compared with the ambiguity surface A(p p p s ) provided by MVDF-MFP in Fig. 10 using a grid of N s = 13 800 possible source position hypotheses.Note that the maximum value of A(•), surrounded by a red box, coincides with the estimated source position provided by MVDF-MFP.
Furthermore, for the same environmental configuration but with source frequency f s = 30 Hz, L g = 15 m, and L f = 30 m, we evaluated the performance of FMM-SE-EnKF in terms of time-averaged root mean square error of the estimated where k * is the discrete time instant corresponding to t * , and RMSE pos,k is the error at time k defined as with p p p s k|k,r being the estimated source position at time k for the r th Monte Carlo run, and M = 50 the total number of Monte Carlo runs.   the source.It can be seen that source localization turns out to be satisfactory at all levels of SNR.Furthermore, Table III also provides position RMSE for a VA with a varying number N h ∈ {10, 20, 30, 40, 50} of hydrophones, SNR = 10 dB, and distance from the source d hs = 1000 m.In this case, it can be seen how performance clearly improves as the number of sensors increases.Note that the reported simulations also show how the space between sensors affects performance.In fact, in the performed Monte Carlo simulations, the N h sensors were regularly spaced to span the entire 90 [m] depth of the surveiled water column so that for each value of N h , a different space among sensors of d * hh = 90/N h [m] is considered.We remark that the choice of varying sensor number N h and sensor spacing d hh while keeping constant the product N h d hh is motivated by the considerations that: decreasing the intersensor distance d hh for fixed N h would reduce the depth span of the vector array with consequent deterioration of localization performance; increasing d hh above d * hh would actually reduce the number of effective sensors Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.As a final remark, we point out that the model used in the filter is different from the one used to generate the data (see the different parameters of the two meshes in Table I), thus showing that the proposed approach is able to provide satisfactory performance even under model mismatch.

B. Scenario With Irregular Seabed
This second scenario involves a motionless source located at ξ = 552.5 m and ζ = 98.5 m (very close to the seabed) emitting a signal of frequency f s = 100 Hz for a simulation time of 10 s with step integration time t = 0.23 ms, within a rectangular domain d with a depth of 100 m, range of 1500 m, and irregular seabed.The environment scenario, depicted in Fig. 6, presents an isospeed configuration (ρ = 1000 kg/m 3 , c = 1500 m/s, and α = 4 × 10 −6 s/m 2 ).Table I provides the parameters of the mesh for the ground-truth simulator and the filter.The parameters used for the FMM-SE-EnKF algorithm are reported in Table II.Also in this case, as for the previous scenario, we chose the size of the square mesh element of the filter such that MSR = 0.5.In Fig. 11 the performance of the filter, in terms of both source detection and localization, is illustrated.It can be seen how data assimilation starts at t 0 = 0.8 s, and at time t * = 1.119 s, the filter is able to detect the signal and localization is satisfactorily performed as the mode in operation êk matches the source element e s after the detection of the source signal.Once the source element e has been detected, the filter runs a local EnKF for finer localization within e s .Furthermore, Fig. 11 plots, versus time, the distance between true, p p p s , and estimated, p p p s k|k , source locations.It can be observed how the local run of the EnKF enhances localization accuracy starting from an initial position coinciding with the element center.This behavior is verified in Fig. 12, where we can see the position samples p p p s,i,e k|k for the source mode e s moving, within the element ê = e s during the data assimilation process, toward the true source position p p p s .Finally, Fig. 13 provides some snapshots of the estimated field x x x k|k compared with the simulated field x x x k at different time instants.

VI. CONCLUSION
In this work, we proposed a fast multiple model spectral element ensemble Kalman filter (FMM-SE-EnKF), a new MM algorithm specifically devised for acoustic USL in shallow water.Such algorithm combines an SEM-based propagation model and Bayesian field estimation for large-scale systems.Unlike approaches based on MFP, the proposed FMM-SE-EnKF can cope with more general space-timevarying underwater environments.In this work, it has been demonstrated how this algorithm can effectively perform detection and localization in the case of irregular seabed.Future work will concern extensions to: 1) moving sources; 2) space-time-varying speed of sound and/or density; and 3) 3-D environments.

Fig. 1 .
Fig. 1.Sketch of the source localization problem in an underwater acoustic environment.
rectangular domain ¯ d is then divided into N e equal square elements ¯ e with side L, each of which corresponds to a curvilinear element e in the physical plane Oξ ζ .Let N ξ and N ζ be the numbers of elements along the ξ and ζ directions, respectively, so that N e = N ξ • N ζ .Let ( ξe , ζe ) be the coordinates of the eth element's center in the plane O ξ ζ , and let u, v ∈ [−1, 1] be local variables such that

Fig. 3 .
Fig. 3. Resampling for constraint enforcement in the prediction step.

Fig. 4 .
Fig. 4. Projection for constraint enforcement in the correction step.

k 25 :
E E E e y,k = [y y y k − ŷ y y 1,e k , . . ., y y y k − ŷ

Fig. 5 .
Fig. 5. Environment configuration of the waveguide with regular seabed within the domain Ω d .Elements Ω e,l with different side lengths l in meters, used for the ground truth and filter mesh, are compared.

Fig. 6 .
Fig. 6.Environment configuration of the waveguide with irregular seabed within the domain Ω d .

Fig. 7 .
Fig. 7. Mode detection trend and distance between estimated source position p p p s k|k and true source position p p p s , versus distance provided by MVDR-MFP algorithm for the scenario with regular seabed.

Fig. 8 .
Fig. 8. (a)-(c) Position samples ps,i,e k|k (in red), estimated position ps,e k|k (in blue), and true source position p s (in magenta) within the source element domain Ω e s at different time steps for the simulated scenario with regular seabed.

Fig. 9 .
Fig. 9. Comparison of ground-truth simulator acoustic field x x x k (upper plot) and estimated acoustic field x x x k|k (lower plot) by the Fast MM-SE-EnKF algorithm at different times, for the scenario with regular seabed.(a) t = 0.87 s.(b) t = 1.54 s.(c) t = 7.48 s.

Fig. 10 .Fig. 11 .
Fig. 10.Ambiguity surface A(p p p s ) for the scenario with regular seabed provided by the MVDF-MFP algorithm where the maximum value is surrounded by a red box.

Fig. 12 .
Fig. 12. (a)-(c) Position samples ps,i,e k|k (in red), estimated position ps,e k|k (in blue), and true source position p s (in magenta) within the source element domain Ω e s at different time steps for the simulated scenario with irregular seabed.

Fig.
Fig. Comparison of ground-truth simulator acoustic field x x x k (upper plot) and estimated acoustic field x x x k|k (lower plot) by the Fast MM-SE-EnKF algorithm at different times, for the scenario with irregular seabed.leaving some of them out of the water column of interest.Finally, Table III shows comparison of the position RMSE for different distances d hs of the array from the source, with fixed SNR = 10 dB and number of sensors N h = 20, showing how localization performance significantly improves when the sensor array gets closer to the source.As a final remark, we point out that the model used in the filter is different from the one used to generate the data (see the different parameters of the two meshes in TableI), , σ w 2 , σ w 3 , σ v std.dev. of pressure, pressure derivative, source position, measurement noise.

TABLE I SPECTRAL
-ELEMENT MESH PARAMETERS FOR BOTH GROUND-TRUTH SIMULATOR AND FILTER, IN SCENARIOS WITH REGULAR AND IRREGULAR SEABED

TABLE III RMSE
pos VERSUS SNR, N h , d hs Table III reports the position RMSE performance for five different levels of SNR, a VA of N h = 20 hydrophones located at a range distance d hs = 1000 m from