An Improved Particle Filtering Technique for Source Localization and Sound Speed Field Inversion in Shallow Water

Both source localization and environmental inversions are practical problems for long-standing applications in underwater acoustics. This paper presents an approach of the moving source localization and sound speed field (SSF) inversion in shallow water. The approach is formulated in a state-space model with a state equation for both the source parameters (e.g., source depth, range, and speed) and SSF parameters (first three empirical orthogonal function coefficients, EOFs) and a measurement equation that incorporates underwater acoustic information via a vertical line array (VLA). As a sequential processing algorithm that operates on nonlinear systems with non-Gaussian probability densities, an improved sequential importance resampling type particle filtering (SIR PF) is proposed to counter degeneracy. The improved PF performs tracking of source and SSF parameters simultaneously, and evaluates their uncertainties in the form of time-evolving posterior probability densities (PPDs). The performance of improved PF is illustrated with well-tracked simulations of real-time source localization and time-varying SSF inversion. Moreover, the influence of different particle numbers on PF tracking accuracy and computational cost is also demonstrated. Simulation results show that the high-particle-number PF has an outperform performance. For a given hardware system, the reasonable compromise between accuracy and computational cost is a matter of tradeoff.


I. INTRODUCTION
In underwater acoustic applications, the passive source localization/tracking and ocean parameters inversion are always hot issues of serious concerned [1]- [5]. Matched-field processing (MFP) is a popular approach used for source localization and tracking [6], [7], which has been applied with excellent results both to synthetic and real data. This approach requires the best match between sound fields measured by receive hydrophones and replica fields computed using a numerical propagation model. However, the modeled field never exactly matches the measured data due to ocean medium fluctuation. The mismatch of environment affects the propagation model that results in biased localization errors.
The associate editor coordinating the review of this manuscript and approving it for publication was Jiajia Jiang .
The sound speed profile (SSP), one of the important environmental parameters has a significant influence on determining acoustic waveguide propagation. Although the SSP can be obtained by combining model outputs or in situ measurement, valid approaches have been developed to invert the SSP in the water volume [8], [9], eliminating the longterm observation and large-scale sampling. To reduce the degrees of freedom, the SSP is often represented in terms of EOFs [10]. In a previous paper, we verified the inversion of SSP has a good agreement with the measured data [11]. The discrete adjacent SSPs in each time interval are combined to form the SSF in the whole period. To invert the time-varying SSF success, EOFs are included in the set of parameters to be estimated.
For a moving source, the path between source and receivers is evolving with time. Especially in shallow water areas, VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ physical variability and dynamic processes cause the environmental parameters (including the SSP) in the propagation change in time and space [12]. These characteristics of stepwise variability can be considered as a tracking problem, which is suitable for the sequential Bayesian algorithm to deal with [13]. Therefore, an algorithm that can simultaneously track the source and SSF inversion is of interest. Sequential Bayesian filters such as extended Kalman [14], ensemble Kalman [11], [15], unscented Kalman [16], and PF [17]- [19] have been used successfully for ocean acoustic applications in previous studies. An excellent overview has described that the sequential Bayesian filtering provide a suitable framework for estimating and updating the parameters of a system as data become available, examples include target localization, spatial arrival time tracking, and geoacoustic inversion [20]. These applications are addressed in a statespace framework. The probability density functions (PDFs) of parameters estimated in a dynamic system are usually non-Gaussian, and the measurement equation is highly non-linear. Although the PF is versatile and robust in handling nonlinear and non-Gaussian estimation problems [21], [22], a wellknown problem of degeneracy results in a poor estimation performance. The particle degeneracy problem in classical PF seriously limits its development. An improved algorithm to counter degeneracy is always a research focus. In addition, recent developments based on the real-time monitoring and control framework, oriented prognosis and diagnosis applications are reported [23], [24].
Guiding by the above-mentioned, our interest in this paper is proposing an improved PF tracking scheme to achieve real-time source localization and SSF inversion simultaneously. PF as a Bayesian method estimates the optimum source and environmental parameters based on the PPDs [25]. Here, two kinds of PFs are used to track the parameters of both moving source (depth, range, and speed) and SSF (three EOFs) including their underlying uncertainties in the form of time-evolving PPDs. Then the tracking results are compared between classical and improved SIR PF. It is expected that the improved SIR PF overcomes the degeneracy phenomenon. With the ability enables us to track the real-time localization of the moving source in range and depth, along with the time-varying SSF inversion. As a tracking algorithm, the computational cost is important too. Because the running time depends on the hardware of the computer, the influence of different particle numbers on PF tracking accuracy and computational cost should also be demonstrated.
The remainder of this paper is organized as follows: the theory of the acoustic propagation model, SSF parameterization, and the state-space model are given in Section II, together with the algorithm of the improved SIR PF. Section III illustrates the approach of source localization and SSF inversion. The computational cost and accuracy compromise for different particle numbers of improved PFs is discussed in Section IV. Finally, Section V concludes this work.

II. THEORY
In this section, the theory of the acoustic propagation model as well as the SSF parameterization is summarized, and then the state-space model for tracking is given, along with the algorithm of improved SIR PF.

A. ACOUSTIC PROPAGATION MODEL
The proposed approach in this work is modeled on a rangedependent shallow water environment. For ranges greater than several depths in shallow water, the normal mode theory is suitable for reconstructing the complicated acoustic field [15], [26]. The acoustic pressure field propagation equation from a point source of frequency ω in cylindrical coordinates (r, z) can be expressed as where p (r, z; ω) stands for the acoustic pressure, c (r, z) is the sound speed and s represents the source. In the range-independent condition c (r, z) = c (z), we can get the normal mode equation where ψ (z) and k represent the vertical modes and the horizontal wave number, respectively. Depending on the boundary conditions, the pressure field is proportional to the sum of the normal modes, multiplied by the Hankel functions (first kind) of the range where H (1) 0 is Hankel function of the first kind. Considering the range-dependent medium as a set of L range-independent vertical segments, the field of each segment is solved by standard normal mode equation. The local modes are coupled through the boundary conditions at the interfaces. Assuming that there are M normal modes, the pressure field in the l th segment can be expressed as where the coefficients a l m and b l m are determined by source and boundary conditions.H 1 andH 2 are the following ratios of Hankel functions of the first and second kind respectivelỹ The continuity condition of pressure at the lth interface is written as Following the derivation by Evans [27], the couple mode formulation requires the solution of a whole family of normal mode problems, one for each segment, followed by the solution of a large banded, block linear system. To obtain numerical convergence in cases of continuously varying properties, the range segments are frequently required to be less than a wavelength, which leads to extremely long computation times. It noted that the numerical calculation of the coupled normal mode is very complicated and time-consuming. To simplify the calculation, adiabatic normal mode theory is used to neglect the energy coupling between modes, i.e. the modes are adiabatic from the current distance (vertical segment) to the next adjacent distance. In the range-dependent environment, the approximate solution of (1) can be derived as during the simulation in this work, using the asymptotic approximation to the Hankel function, the acoustic pressure is calculated in the form of where ρ is the density, z is the source depth, r is the source range, ψ m (·) is the eigenfunction for m th mode and k rm is the eigenvalue.

B. SSF PARAMETERIZATION
In this work, the inversion of SSF is carried out in the unit of SSP. An independent SSP can be measured in each time interval. Combining these discrete adjacent intervals, the approximate distributions of sound speed variation in continuous time are obtained. EOF is a reasonable way to fully parameterize the SSP [28]. It can be obtained from a direct measurement of the database and are efficient in reducing the degrees of freedom of the estimated parameter vector.
Assuming that there are N SSP measurements, the mean SSP is calculated as where performing an eigendecomposition of a covariance matrix R provides the required EOFs (eigenvectors) f n with corresponding eigenvalues defined by λ n . Therefore, any SSP in the ocean can be approximated with the first k eigenvectors as where α k and f k are called the kth EOF coefficient and EOF, respectively. In this way, any sound-speed estimates should be constrained to lie in the space spanned by the EOFs. To handle the range-dependent environment, that is, dividing the range between the source and the VLA into L segments using K EOFs, each segment has only one SSP, which is expressed as where g l (r) is defined as a gate function For the l th segment, the sound speed is calculated in the form of where α k,l is the k th EOF coefficient for the profile in the l th segment. In this way, the SSF is recovered by a total of l SSPs, and α k,l is the corresponding SSF parameter.

C. A STATE SPACE-MODEL FOR TRACKING
Tracking of the moving source and SSF inversion requires two dynamic equations: A state equation that models the movement of the source and the evolution of SSF, and a measurement equation that relates the environment and source localization at step k to the acoustic field. Two equations define a state-space model are given as where f (·) is a state transition function of the state vector x k−1 , and h (·) is the nonlinear function that relates the x k to y k . State and measurement noise covariances Q k and R k are The state vector in (16) includes the parameters of source and SSF is given by merging these two blocks as x T k = s T k α T k .s and α denote the parameters of source and SSF respectively. Using a constant velocity (CV) track model for the moving source, the depth Z , horizontal range R, and radial speed v are formed as the first block where t is the update time increment (i.e., the time interval between two measurements), v d and v a are the random variables representing the variation in source depth and acceleration, respectively. As a multi-dimensional tracking problem, the computational complexity and estimation accuracy often depend on the number of estimated variables. In practice, the number of EOFs representation can be reduced. It is known that the first three EOFs are capable of representing SSP [11], the second block is given by In (17), the measurement equation is used to characterize the signal of acoustic pressure, which is received by a VLA. According to the normal mode theory in Section II.A, the synthetic acoustic pressure can be expressed as a general nonlinear function of the source depth Z , range R, sound speed c and sea bottom boundary condition ( BCs) where y k is the acoustic pressure measured at k th time frame, h is the normal mode propagation model. Considering the SSF is parameterized in terms of the first three EOFs, the resulting measurement equation becomes where α = [α 1 α 2 α 3 ] T . As mentioned above, the state equation describes the source localization and evolution of the SSF, and the measurement equation relates the state vectors and ocean environment to the acoustic field. In this way, the state-space model is used as a sequential method enable to track the evolving state.

D. PARTICLE FILTERING ALGORITHM
The state-space model allows straightforward implementation of sequential filtering algorithm. However, the posterior PDFs of source parameters and EOFs are usually non-Gaussian, and the high nonlinearity of the measurement function h (·) requires nonlinear filtering methods. As discussed in Section I, the PF is appropriate. Normally, degeneracy can be a problem for the classical PF algorithm that results in poor performance, especially for the low process noise system. The degeneracy will cause adverse effects: the posterior probability is only represented by a few particles with larger weights, and the contribution of most particles to PPD is close to 0. In addition, it will cause a waste of computational resources. A lot of computational resources are wasted on the particles that make little contribution to state estimation.
To prevent the degeneracy phenomenon, an improved-SIR PF is proposed, which uses the likelihood function of optimization for guiding the importance sampling particles to the high likelihood regions to ensure the resampling process.
Following (17), the additive complex Gaussian noise for step k, the likelihood function is of the form where n h is the number of hydrophone of VLA, and w j is the noise variance. A single iteration at step k of the SIR PF is summarized as follows.
Predict: Using a set of particles from previous state based on a given state model described in (16).
Update: Update the predictions in the previous step using the acoustic pressure data y k at the current step k, evaluate the normalized weight w i k of each particle x i k|k−1 from its likelihood function where p (x k |x k−1 , y k ) is the likelihood function calculated by (23), we have p y k |x i k|k−1 = L x i k|k−1 . The posterior PDF p (x k |x k−1 , y k ) is expressed by using the weight w i k as Resample: New set of particles x j k|k , is drawn from the discrete approximation density p (x k |x k−1 , y k ) obtained at the update stage. All particle weights are now equal to 1/N p . The stage of resampling integrates both the predictions from previous values and the information from the new measurement.

III. RESULTS AND DISCUSSION
In this section, simulation studies are carried out to demonstrate the performance of the proposed approach, consisting of the following contents: tracking source and SSF parameters, evaluating how the PF tracks their uncertainties as PPDs, and the results of real-time source localization and time-varying SSF inversion.

A. TRACKING PARAMETERS OF THE MOVING SOURCE AND SSF
To synthesize the acoustic pressure data used in the measurement equation, the environment model used is shown in Figure 1, the simulation setup involving a fixed VLA and moving source in a range-dependent shallow water environment. The water depth is 106 m, and the modeled bottom is homogeneous, with sound speed 1610 m/s, density 1.7 g/cm 3 , and attenuation rate 0.1dB per wavelength. These values correspond to average values deduced from the core measurements of AXIAEX, ECS, 2001 Experiment [29]. The VLA is configured with 16-element equispaced 4 m that spanned the 15-75 m water depths. The initial position of the source is at 40 m in depth and 1 km in range. A frequency of 400 Hz is selected. The source has 180 steps (k = 1, · · · , 180) movement at a speed of 2.5 m/s, resulting in a total track length of 1 hour during which the source moves from 1 to 10 km in the horizontal range.
All the simulation parameters, environmental constants, six parameters as state variables and their start values x 0 , prior standard deviations P 1/2 0 , state noise Q k were selected as given in Table 1. Selection of suitable noise covariances for the PF is essential to tune the filter. For the convenience of calculation, we assumed that both additive noise terms Q k and R k were represented by the Gaussian PDFs. The adiabatic normal mode propagation model was used to compute the acoustic field by inputting the environmental and source parameters into KRAKEN [30].
According to the state-space model formulated in Section II.C, the source depth, range, speed and three EOFs are taken as the state variables. Two types of PFs are used to track these parameters simultaneously. Higher approximation accuracy can be obtained for the posterior distribution with a larger number of particles. To ensure the tracking accuracy, the number of particles is set to N p = 500. All the particles are propagated through the state equation, updated via the measurement equation and resampled, and then to predict the evolution of state variables.
Here the same setting, both classical SIR PF and improved SIR PF are used to track the source and SSF parameters. The tracking results of two types PF are given in Figure 2    The evolutions of the 1-D marginal posterior densities for improved SIR PF are given in Figure 3. Compared with the classical SIR PF, it overcomes the particle degeneracy. The plots in the left column belong to the source parameters. The depth and range of source are tracked accurately with small deviations from the true values. Variation in the radial source velocity is also well tracked with the constant velocity evolution model. The three SSF parameter tracks are in the right column. The improved SIR PF also does a good job of tracking the EOFs. Even though SSPs in the water column are expected to be evolving slowly most of the time, there are rapid perturbations of sound speed in some regions of the true environment. The perturbations result in larger variations of the EOFs, which exceed the statistical values assumed for the random walk, and resulting in a model mismatch.
To continue tracking the EOFs through the perturbation, the PF will have to incorporate a state noise term v k with a high covariance Q k . Therefore, the evolution of EOFs in some periods is tracked with deviations from the true trajectories because of the noise v k . Moreover, it is seen that the tracking of EOF3 is not as good as the other two. The possible reason is that EOF3 plays a minor position compared with the first two on dominating the SSP characterization, which leads to a minor variation of the sound speed (in (13)) and little pressure variation information in synthetic acoustic pressure data, so the PF cannot fully capture the variation of EOF3.

B. THE REAL-TIME SOURCE LOCALIZATION AND TIME-VERYING SSF INVERSION
The improved SIR PF is used throughout the following work. To simplify the expression, all the PFs mentioned below refer to improved SIR PF.
In the above contents, we discussed the PFs tracking scheme, and gave the time-evolving marginal posterior PDFs of the six parameters. It is also of interest to observe the temporally track at each step k. A vertical slice from each of the evolving PPDs in Figure 3 is given in Figure 4 for t = 19. This snapshot shows source depth, range, EOF1 and EOF2 are well-tracked at this time with narrow and sharp densities, whereas the probability density of EOF3 is dispersed and flat in the search interval, which indicates the increase of tracking uncertainty. This result is consistent with that in Figure 3. In addition, the probability density of source speed is concentrated, but the peak deviates from the true value, resulting in unreliable tracking result at this time point.
As a sequential filtering algorithm, the PF enables temporal tracking of source and SSF parameters and their underlying probability densities, which allows real-time locating of the moving source and updating of the SSF.
For the moving source localization, we investigate the PF performance by computing and plotting the posterior source depth-range probabilities. The evolution of the 2-D PPD of the source at three time steps is shown in Figure 5. The results at 7, 29 and 54 min point of the track are given in Table 2. The PF makes only 0.1-0.2 m, 1 m and about 0.01 m/s error at the three mid-tracks for the source depth, range and speed, respectively. Even after the full track length of 60 min, the error terms stay at 0.04 m, 1.74 m and 0.003 m/s. Therefore, this real-time localization of PF enables us to track the source accurately in range and depth with small errors.
As the tracking time goes on, the PF utilizes both current data and previous parameter estimated result to update the parameter at the next step. This enables us to obtain  histograms of particles that represent the marginal posterior PDFs of the parameters estimated in any period time. Selecting a 12-minute tracking duration, time-evolving marginal posterior PDFs for three EOF from 7 to 19 min are given in Figure 6. Vertical slices of the PDFs are given at every 2 min intervals. The first two EOFs are well estimated through the track and the posterior PDFs correspond to narrow Gaussian PDFs, and most of the peaks overlap the true values at the corresponding time interval. Even though the relatively large uncertainty for EOF3 is observed, it plays a minor position in dominating the SSF characterization. Therefore, the tracking results of EOFs are reliable. The inversion of SSF is performed next.
As described in Section II.B, the time-varying SSF is recovered by a total of EOFs estimated in the whole tracking period. Figure 7 shows pseudo-color maps of the SSF inversion result and relative error between inverted SSF and the true one. Moreover, the errors near the thermocline are large compared with other regions because sound speed changes fast in this region. This fast-changing feature results in larger variations of the EOFs which resulting in a model mismatch at that time.

IV. COMPUTATIONAL COST AND ACCURACY COMPROMISE
Desired accuracy in tracking results is one of the major factors in selecting the number of particles. As a tracking algorithm, the computational cost is important too. In this section, the influence of different particle numbers on PF tracking accuracy and computational cost is demonstrated.
Each of those six parameters is tracked by PFs using 50, 100, 200 and 500 particles designated by PF-50, PF-100, PF-200 and PF-500, respectively. The tracking results of PFs are given in Figure 8 along with the true trajectories. Only PF-50 has a large tracking error when the true value changing abruptly. Too few particles resulting in sample impoverishment, thereby losing diversity and leading the PF to divergence. Noting that increasing the particle number improves the PF tracking performance obviously.
For the convenience of calculation and comparison, the inverted SSF is taken as a gathered parameter. The performance of four PFs for source parameters and inverted SSF are evaluated by the root-mean-square-errors (RMSEs) between the true values and tracking results, are shown in Figure 9.  It is confirmed that both PF-500 and PF-200 outperform the PF-50 and PF-100 over the whole tracking period. To quantify the average tracking performance for each filter, TARMSE is defined as the time-averaged root mean square error calculated for the interval (k = 180).
Increasing the number of particles initially provides large performance improvement in a PF, see Table 3. Even though PF-500 outperforms the PFs in terms of TARMSEs for all the parameters, it is also necessary to compare the computational cost of each filter. Because the computational cost depends on the hardware of the computer, and the number of operations a filter has to perform at each step. In the following, we employ the computational operation O (·) to evaluate the computational cost.  According to Section II.C, we assume the number of state vector is S, and the size of measurement vector is M . The computational operation of sound propagation model is O (h). After ignoring the lower order term, a total computational operation of the PF is evaluated as O (N P (h + Sµ + M + S)), where N P is the particle number and µ represents the random number complexity. From the equation, the computational cost of PF is critically affected by N P .
For a normal desktop computer (AMD Ryzen 9 3900x 12-Core processor with 32-GB memory) running the programs, the computational time is shown in Table 4. The number of operations required for propagation model calculations is significant relative to the filtering algorithm. It takes a lot of time to calculate the sound propagation model in the measurement equation. In order to save the computational cost, it is necessary to reduce the N P , that leads to the compromising accuracy. To assess the degree of compromise, equation (26) is used to calculate the accuracy compromise of a filter with respect to the PF-500.  The compromise over PF-500 percentages and computational times of the PFs are given in Table 4. We can see that the computational time doubles as the N P increases. PF-50 and PF-100 have a high degree of accuracy compromise over PF-500. However, PF-200 achieves 2.7%, 29.9%, 10.1% and 10.6% accuracy compromise respective for source depth, range, speed and SSF, but saves nearly 2/3 of computational time.
An upper limit for particle number is determined by the maximum number that can be processed with limited computational resources, which is important especially for real-time filters. However, after an accuracy-dependent particle number is reached, the performance stays relatively flat. Increasing the computational cost provides only marginal benefits [20]. In practice, based on the hardware employed to run the algorithms, a reasonable compromise between accuracy and computational cost is always the matter of tradeoff.

V. CONCLUSION
In this paper, an approach of moving source localization and SSF inversion in shallow water environment has been addressed and illustrated. The parameters of source (e.g., VOLUME 8, 2020 source depth, range and speed) and SSF (e.g., first three EOFs) were tracked simultaneously by assimilating measurements of the acoustic filed on a VLA into classical and improved SIR PFs. To counter degeneracy, an improved SIR PF was demonstrated and used throughout this work. The improved SIR PF enabled providing real-time, continuously updated tracking of the parameters and their underlying uncertainties in the form of a time-evolving PPD. The capabilities of the approach were illustrated on well-tracked simulation for source localization together with time-varying SSF inversion. The improved SIR PF proved to be a promising algorithm in the nonlinear, non-Gaussian underwater acoustic tracking problem. Furthermore, the influence of different particle numbers on PF tracking accuracy and computational cost was also demonstrated. It should be pointed out that although increasing the number of particles provides large performance improvement, its computational cost is greatly increased. Therefore, the choice of accuracy and computational time for a given hardware system is always a matter of tradeoff.
Source localization and environmental parameters inversion are practical problems for long-standing applications in underwater acoustics. The approach we proposed has a useful value for studying the inversion of dynamic ocean parameters and target localization/tracking in more complex environments. In practice, it provides a reference solution to real-time monitoring and control problems in marine science and ocean engineering. He is currently a Professor with the School of Marine Science and Technology, NPU. His current research interests include underwater acoustical signal processing, nonlinear dynamical modeling, and target tracking.
JINYING YE received the B.S. degree in thermal power engineering, the M.S. degree in naval architecture and ocean engineering, and the Ph.D. degree in aerospace propulsion theory and engineering from Northwestern Polytechnical University (NPU), Xi'an, China, in 2010, 2013, and 2018, respectively.
Since 2019, he has been an Assistant Researcher with the School of Astronautics, NPU. His current research interests include rocket-based combinedcycle technology and variable geometry combustor technology.
KUNDE YANG (Member, IEEE) received the B.S., M.S., and Ph.D. degrees in underwater acoustic engineering from Northwestern Polytechnical University (NPU), Xi'an, China, in 1996China, in , 1999China, in , and 2003 He was a Visiting Scholar with the School of Earth and Ocean Sciences, University of Victoria, Victoria, BC, Canada, from 2006 to 2007. Since 2003, he has been with the School of Marine Science and Technology, NPU, where he is currently a Full Professor and the Director of the Department of Acoustic and Information Engineering. His research interests include ocean acoustic modeling, signal processing, and microwave propagation.