Sequential Geoacoustic Inversion Using an Improved Kalman Particle Filter

Geoacoustic inversion is an efficient method to study the physical properties and structure of ocean bottom while sequential geoacoustic inversion is a challenging task due to the complexity and non-linearity of the underwater environment. In this paper, an ensemble Kalman Particle filter is described to address the sequential geoacoustic inversion problem of range-dependent environment in shallow water. This filter combines the advantages of Particle filter and ensemble Kalman filter so its ability of tracking dynamical geoacoustic parameters is improved. The proposed filtering method is demonstrated with simulated data in a changing oceanic environment and outperforms Particle filter and ensemble Kalman filter. This method is also tested in sea-trial data collected from a shallow-water experiment in the East China Sea. The inverted sound speed in sediment is consistent with in situ measurement and the error between transmission loss predicted by inverted parameters, and the experimental transmission loss is small.


Introduction
The characteristics of seafloor sediments such as sound speed have significant effects on the propagation of underwater acoustic signals. Geoacoustic inversion is an effective method to investigate the properties of sediments without sampling from seafloor directly. A variety of different inversion techniques have been applied in geoacoustic models, including matched field processing [1][2][3], bottom reflection analysis [4], and model dispersion inversion [5,6]. In matched field processing, array-processing algorithms that exploit the full-field structure in the ocean waveguide match the measured field at the array with replicas of the expected field for the location of source and acoustic parameters from propagation models. If the location of the source and the received array are known, searching over a parameter space for unknown geoacoustic parameters can be considered to be an inversion problem, referred to as matched field inversion (MFI) [7]. In general, the complexity of inversion problems increases with the number of unknown parameters and the search for the global minimum becomes very time-consuming. The genetic algorithm (GA) [8] and adaptive simulated annealing algorithm (ASSA) [9] are widely used global optimization tools. The inversion also needs to consider the uncertainty of obtained parameters and relies on Bayesian inference to obtain the multidimensional posterior probability using the Markov Chain Monte Carlo technique [10,11].
Changes in ocean and seafloor parameters will affect the acoustic field and make it "evolve" in space or time. The estimation of geoacoustic parameters of such a varying field can be reformulated as a sequential filter-tracking problem. Kalman filter (KF) framework was first used to process underwater acoustic target tracking and localization and has also been widely used for frequency tracking, estimation of environmental parameters and arrival-time tracking. A review about the application of sequential filters in ocean acoustics is given in Ref. [12]. KF uses two equations to describe a state-space model that varies with the time series, taking into account the relationship between the measured and state parameters as the estimated state parameters evolve with the sequence. The classic Kalman filter is based on the linear noise assumption of a Gaussian distribution, but the interaction between geoacoustic parameters and acoustic field has a high level of non-linearity and non-Gaussianity. The use of classic Kalman filtering method for geoacoustic tracking therefore produces large errors.
As a result of complexity and non-linearity, the extended Kalman filter (EKF) [13,14] and unscented Kalman filter (UKF) [15] are applied in sequential tracking and environment parameter estimation. Meanwhile, the ensemble Kalman filter (EnKF) [16][17][18][19] is a Kalman filter with Monte Carlo (MC) sampling and assumes a large number of state variables and all the probability density functions (PDFs) are Gaussian. Compared to classic KF method, EnKF uses ensemble sampling covariance matrix to replace the covariance matrix when computing Kalman gain, so it would perform better for a high-dimensional parameters problem. These improvements mean that EnKF can be used for non-linear and high-dimensional systems.
Another type of sequential Bayesian filter is Particle filter (PF) [20][21][22][23]. PF estimates system parameters by setting a series of particle points with weights and then updating them according to weights of particles. It does not need to make assumptions about the form of posterior probability, so it can be applied to non-linear and non-Gaussian systems. In addition, it can represent the entire posterior PDF of state variables, which means that PF can better reflect the evolution in non-linear systems. However, PF has a problem with particle degeneracy [12,24]. The weights of most particles will be negligible as iteration proceeds. One method used to address this problem is resampling. During resampling, particles with high weights can be copied many times and include particles with low weights.
Many efforts have been made to combine different types of filters. In Refs. [25,26], EnKF was suggested to produce the proposal distribution of the PF to improve its performance. This combined method is referred to as the ensemble Kalman Particle filter (EnKPF). Because the proposal distribution generated by EnKF contains information on the latest observations, the new particles are closer to observations. This method incorporates the latest observations into weight calculation and reduces the risk of particle degradation. EnKPF has been demonstrated in data assimilation [27] and sound speed profiles (SSP) tracking [28].
Taking into consideration the non-linear and non-Gaussian properties of shallow-water geoacoustic inversion problem, we propose EnKPF method to track a moving source and estimate spatial variation of geoacoustic parameters. The vertical line array (VLA) was fixed in one position with hydrophones placed from sea surface to seafloor and the source transmitted signals continuously at different locations, so the acoustic fields measured from VLA were observed data of inversion. The tracking was initiated with the result of a full environment GA inversion implemented at the first time index. To compare the tracking capability, we used synthetic data to compute mean errors of EnKF, PF and EnKPF in the same dynamic environmental model. Satisfactory results were obtained when sequential inversion tracking was implemented using EnKPF combined with the experimental data from a shallow-water acoustic experiment. The results of sequential inversion were verified by in situ measurement and experimental transmission loss (TL) of signal.
In this paper, the EnKPF is described in the context of the application of sequential geoacoustic inversion and shows superiority in the shallow-water environment model. Section 2 outlines the theory of Bayesian sequential inversion, geoacoustic environment model and gives a brief derivation of EnKPF. The performance of proposed filtering algorithm is compared with EnKF and PF, and implemented with simulation data in Section 3. A sensitivity analysis is also implemented. Section 4 presents sea-trial data in a shallow-water experiment. Section 5 presents some discussions, and a conclusion is presented in Section 6.

Sequential Geoacoustic Tacking Modeling
The sequential tracking problem consists of two dynamic equations: (1) a state equation that describes the evolution of geoacoustic parameters and the movement of source; (2) a measurement equation that relates geoacoustic parameters and source location at step k to acoustic field. These are given as where f (·) is a known state function of state vector x k−1 , and h(·) is the function relating the environmental and source parameters x k to the acoustic field y k . v k and w k are state noise and measurement noise, respectively. Because geoacoustic tracking includes two parts (geoacoustic sequential inversion and moving source localization), the state equation is separated into two blocks. The first block deals with the source tracking. Considering a constant-speed movement [22], source location is determined by three parameters: source depth z, range r and corresponding speed v. They are grouped as s k = [z k , r k , v k ] T . Then the state equation for source tracking is presented by and the matrix form is where F s is the matrix of state function in the case of source tracking, B s is scaling matrix and ∆t is the time interval between measurements. The assumption of constant speed is difficult to satisfy in the actual environment. Hence the random variables (v z and v a s ) representing the variation in source depth and acceleration are introduced with the assumption of a Gaussian distribution. The second block of state equation is used for geoacoustic inversion. The model assumes that the rate of change of parameters is slow for each step. If the rate of change in geoacoustic parameters is rapid (i.e., a sudden change in water depth), then the corresponding state noise matrix v m k can be set with a related high covariance Q k to continue tracking successfully.
The full state vector, including both source tracking and geoacoustic inversion, is x T k = [s k , m k ] T and the full state equation can be given as [22] s m The measurement process is implemented by measurement equations with the forward model, in most commonly normal mode method or parabolic equation method in a shallow-water environment. In this paper, the measurement vector y k is acoustic pressure field data received in a VLA in which hydrophones are spaced over the whole water depth. The measurement error w k contains the error of model and receiver signal-to-noise ratio (SNR) in the array. It is easy to understand that the measurement error will increase in an environment with a low SNR.
In this paper, geoacoustic inversion likelihood formulation with an unknown source signal is derived from the classical geoacoustic inversion likelihood with the Bartlett power objective function [2]. The noise for frequency f j at step k is usually assumed as additive complex Gaussian noise, which is written as w k ( f j ) = y k ( f j ) − a k d(x k , f j ) [22]. The corresponding multifrequency likelihood function is [29] L = where φ j (x k ) is the Bartlett objective function and tr is trace operation.
is the cross spectral density. n h and n f are number of hydrophones and frequency points used and ν j is the noise variance in frequency f j , y k ( f j ) is the observed acoustic field at frequency f j , a k is source amplitude and d(x k , f j ) is the computed acoustic field from forward model. We can use the maximum likelihood estimator to obtain unknown noise variance ν j as [29]

EnKPF as a Sequential Tracking Processor
As discussed before, EnKF and PF have a wide range of applications. EnKF can be considered to be a Kalman variant for a system in which the state is composed of a large number of variables. Its basic idea is to use a series of ensembles to approximate PDF with the assumption that all ensembles and PDFs are Gaussian. EnKF improves the ability to handle high-dimensional problems, but working with a large covariance matrix in a high-dimensional system is computationally inefficient for KF framework. PF is an implementation of Bayesian filtering by sequential MC sampling to estimate posterior PDF p(x k |y 1:k ) in non-linear problems. The principle of PF is to represent the posterior PDF of a state vector by a cloud of random particles with associated weights. So, PF performs well in non-linear systems. A key issue of PF design is choosing a good proposal density. The optimum choice of proposal density is the posterior PDF itself, which minimizes the variance of importance weights, but it is hard to sample particles from this distribution. A common choice is to sample from the transit density of state vector [30], but it makes sample impoverishment [12]. Because transit density does not use the newly available observations, leading to a "blind" prediction. The generated particles may be too far away from observations to properly represent the posterior PDF of state variable and may be more severe in the case of small observation error. Then it becomes vulnerable to outliers, especially when the likelihood L is peaked relative to prior.
An improved method is EnKPF, which uses a PF procedure to update particles and then uses EnKF analysis PDF as the proposal density. Because the proposal density includes information of latest observations, proposed particles are much closer to observations than those from original PF. Therefore, more particles can obtain non-negligible weights, reducing the chance of particle degeneracy and impoverishment. The insertion of EnKF into PF gives a better approximation of the proposal density and increases performance of filter.
In this study, we generate an ensemble {x m i,k , m = 1, 2, ..., M} in every particle {x i,k , i = 1, 2, ..., N p } and implement an EnKF process to create a Gaussian distribution as where i and m are indexes of particle and ensemble, respectively. The proposed particlex i,k is sampled from a Gaussian distribution defined by meanx i,k and covariance matrixP i,k . After this step, the new particle is closer to likelihood because a new observation is used in the update step of EnKF. The density of transition from x i,k−1 to x i,k is assumed to be a Gaussian distribution, which is given as Therefore, the full formulation of particle weight can be rewritten as [28] After resampling, new particles are drawn from the approximation to density p(x k |Y k ) [12], so their weights are equal to 1/N p , then state estimation can be written aŝ

Simulation
This section presents the simulation of geoacoustic parameters and acoustic source tracking in a changing environment. The environment model is shown in Figure 1 and SSP in water is assumed to be known and range-independent. In the environment model, signals are received by a VLA with 20 hydrophones varying in depth from 5 m to 100 m and source is moving away from the VLA at a constant speed. The environment parameters, including water depth (d water ), sediment thickness (h sed ), sound speed (c sed ), density (ρ sed ), attenuation (α sed ), and geoacoustic parameters in the basement are inverted sequentially. The range (r s ), depth (z s ), and radial speed (v s ) of source are also treated as unknown parameters. Almost all the unknown parameters vary slowly with time. Corresponding to Equation (1), the state vector x k consists of environmental and source parameters, y k is the received complex acoustic pressure field data at each hydrophones and h(·) is forward model of acoustic propagation relating the x k to y k . A sensitivity analysis of inversion parameters is carried out for different frequencies at k = 0 before the implementation of simulation, shown in Figure 2. The method used for sensitivity analysis is to compute the objective function while varying one parameter and fixing other parameters at their true values. The normalized objective function is given as [21] It is important to select an appropriate frequency for inversion, because the sensitivity curves of different frequencies are diverse. Different colored curves in Figure 2 indicate objective function values for different frequencies at 150, 250 and 300 Hz. For certain parameters, different frequencies would show significantly different sensitivities and when deviating from true value, the object function for more sensitive frequency would fall faster. Lower frequency signal penetrates the seafloor more easily and can be used to estimate the parameters of bottom layer, while the higher frequency signal is more sensitive to sediment layer and location of the sound source. However, the higher frequency signal is more difficult to detect at a longer range as a result of greater attenuation. The sensitivity plots show significant differences at different frequencies of the same parameter. For example, parameters of the sediment layer are more sensitive than parameters in the bottom layer, because sensitivity curves show that the density and attenuation of the basement have very flat curves for all frequencies considered. For this reason, in the bottom layer, only the sound speed (c base ) is set to be inverted and other parameters are considered to be constant in simulation. So the inverted parameter space consists of 9 parameters: the water depth (d water ), sediment thickness (h sed ), sound speed (c sed ), density (ρ sed ), attenuation (α sed ), sound speed (c base ) in basement layer, the source-receiver range (r s ), depth (z s ), and radial speed (v s ).  The measurement error w k in the simulation is computed by the array SNR formula to synthesize observed acoustic data [2], In this paper, the array SNR is 30 dB. A frequency of 250 Hz was selected, and 100 measurements were made at 6 s intervals (∆T = 6 s) to give a total tracking time of 10 min. The source traveled from 2 to 5 km in that time. The sediment thickness remained constant in the first 6 min and rapidly increased to 25 m after 8 min, then stayed at this depth until the end of simulation. Parameters of sediment layer such as sound speed, density and attenuation were also range-dependent and changed with the movement of source. Although parameters of the bottom layer were found to be insensitive in the inversion, it is set to a time-varying parameter to test the performance of filters under this condition.
The prior result in k = 0 was assumed to be known and each parameter was consistent with a Gaussian distribution whose covariance matrix is P 0 . The state noise Q k affected the scope of evolution with time. The greater the state noise and prior deviation, the more likely the filter was to track rapid changes in the environment and the greater risk of divergence, and vice versa. Both the length of ensemble and particle number were set to 200 to compare the tracking performance among PF, EnKF and EnKPF techniques. Because EnKPF uses EnKF as proposal distribution of each particle in the PF, its computational cost would increase obviously. These parameters and their prior deviation P 1/2 0 , state noise deviation Q 1/2 k and other variables are given in Table 1. The evolution of parameters is shown in Figure 3a with true trajectories given by solid lines. The tracking results of EnKF, PF, and EnKPF are shown as different colors, respectively. Table 2 gives the time-averaged root mean square error (RMSE) of three filters. Overall, EnKPF performed best, followed by PF, and EnKF showed the poorest performance. Water depth is a sensitive parameter, so all three filters tracked rapid changes well. EnKPF and PF closely match source depth and range of true values in the tracking. In the constant-speed simulation environment, tracking error from one of the source range or velocity would cause the mismatch of the other. The source velocity is not a geometric parameter that directly affects acoustic field like the source depth, but acts on distance together with the sampling interval. Considering the evolution trajectory of source speed, PF is more oscillating than EnKPF, which causes the error of distance to be larger. Though EnKF performs better in source speed, mismatches in other parameters (sediment thickness, sediment sound speed) accumulate errors and cause the divergence in tracking of source range. Reference [18] indicated that if PDFs become non-Gaussian, the performance of Kalman filter will be significantly reduced, which means that when SNR decreases or range increases, EnKF may lose track of location of the source. The slight oscillations occur in the tracking of the sediment layer (sound speed and thickness) which become stronger varying with time. EnKPF and PF control the oscillation and rapidly track changing trajectory well, while RMSE of EnKPF is less than the latter. In addition, the lower sensitivity of density and attenuation in sediment layer (ρ sed , α sed ) leads to divergence, especially for EnKF. For the considered frequency, sound speed in the bottom layer shows insensitivity in the Figure 2, so these filters cannot fully track it. This phenomenon is similar to previous research [21]. To overcome this problem, the number of particles and ensemble need be larger, while increasing the computational costs. Furthermore, to test the performance of EnKPF initialized in a random model, we make sediment and basement parameters deviate from their true values in Figure 3a. The tracking trajectories shows that only relative sensitive parameters (h sed and c sed ) can return their true values in time.
To compare the performance of PF and EnKPF, the same environment was tracked using both PF and EnKPF with different number of particles, shown in Figure 3b. Their time-averaged RMSE are also listed in Table 2. With increase of the number of particles, the time-averaged RMSE of most parameters are decreased. It is worth noting that when number of particles is small, EnKFP is superior to PF. On the other hand, though EnKPF shows the superiority in geoacoustics, this does not mean that it applies equally to other issues. In some systems, such as geophysical data assimilation, the dimension of state vector would be larger than geoacoustic inversion problem. For this situation, the selection of particle number should not simply be increased but carefully considered.    Figure 3.

Experiment Data and Implementation of Filtering
The 1000-point EnKPF was used on acoustic field data collected in a shallow-water acoustic propagation experiment in the East China Sea in winter, 2014. The topography is shown in Figure 4a. Many papers have been published on geoacoustic environment in the East China Sea, and some conclusions were obtained [7,[31][32][33]. A VLA consisting of 32 hydrophones was deployed in water with a depth of 106.8 m. The first hydrophone was deployed at 24.5 m and the last was placed at a depth of 102 m. During the experiment, an air-gun source was towed at about 2.75 m/s at a depth of 10 m and its trajectory was recorded by GPS, as shown by the red line in Figure 4a. In the same figure, the in situ site is marked as a downward triangle. This in situ system can penetrate about 1 m above sea floor then compute the sound speed of surface sediment by acoustic probes. Because the topography is flat, the straight trajectory has almost constant water depth when ship is moving. The average SSP is given in Figure 4b. Because the experiment was conducted in a windless winter, the variation of SSP was small in both in time and space. Most of the energy of the air-gun signal was limited to a relatively low frequency band and frequencies of 41, 54, 83, 117, 128, 149, 181 and 231 Hz were used to calculate the acoustic field in the normal mode [34], because the water depth is an almost constant in experiment, the slightly range-dependent effect could be ignored.
The first signal was recorded by VLA at about 10:40 UTC, then the ship towed an air-gun away from VLA at a roughly constant speed. Figure 5 shows the signal received at VLA from air-gun source; hydrophone 13 and 32 failed during deployment. Under ideal conditions, VLA would receive pulse signals from the air-gun every 120 s and the distance between two signals transmitted by source would be about 330 m. The ship stopped repeatedly during this period to conduct other experiments, so signals transmitted in these gaps were not recorded.
Before implementing the sequential inversion, signal with index k = 0 was used to carry out a matched field inversion to obtain prior information about geoacoustic parameters, and results were used as initial vector in sequential inversion, which are listed in Table 3. Another important aspect is choosing an appropriate model parameterization, which can be solved by Bayesian information criterion (BIC) [11]. BIC balances the contradiction of lower misfit and fewer free parameters where E(m) is misfit value of maximum a posteriori (MAP) model, M is the number of parameters and N is the size of data set. The parameterization with the smallest BIC value is selected as the most appropriate model. The results of model selection is given by 3 figures: Figure 6a shows the misfit values decrease from 1 layer to 3 layer while number of parameters increase in Figure 6b,c shows the 2 layer model has the minimum BIC values.    The measurement error is related to the SNR, although, in reality, the measurement noise is already included in observed acoustic field with no additional consideration. The SNR decreased as distance increased, leading to a lower reliability of the source location and geoacoustic inversion. Another factor influencing tracking results is state noise Q k . The selection of Q k is a trade-off between how much we trust the time evolution model and how noisy an estimate we accept [22]. If the range-dependence is strong or water depth changes fast, Q k is considered to be larger. Otherwise set it to small values. Some research gives references about the selection of P 0 and Q k in different environments [19,[21][22][23]. The P 0 and Q k in our work are given in Table 3. The tracking results, evolution of the marginal posterior probability distribution (PPD) of geoacoustic parameters and source parameters are given in Figure 7. It can be seen that the source parameters are tracked closely. The source depth z s and velocity v s remain roughly constant during tracking. The source range r s is tracked well in the first 50 min, but decreasing SNR accumulates errors in the tracking. Moreover, the discontinuity in r s was due to several stops and drifts by ship. The water depth d water changed as the ship moved, but a rough estimation of the depth can be obtained from Figure 4a, which shows that water depth changed from about 110 m to less than 100 m. In addition, the exact depth was measured by a bathymeter and is drawn as a blue solid line in Figure 7. The evolution of d water with time is consistent with the depth measured onboard. Sediment thickness shows a clear fluctuation between 19 and 10 m in the first 10 min, then drops to about 5 m until the end of the experiment. The sound speed in sediment layer increases from 1670 to 1700 m/s; the nearest in situ site, which is labeled with the downward triangle in Figure 4a, shows that sediment sound speed is 1675.78-1686.08 m/s with an mean of 1680.72 m/s. The density of sediment layer ρ sed increases very slowly with the increase in c sed through sequential filter. Figure 8 gives a snapshot of the evolution of marginal PPDs at t = 76 min in Figure 7, where dashed lines are estimated value for each parameter in that time period. Moreover, the measured source range and water depth are drawn as solid lines in the second and last chart, showing that estimated source range and water depth are close to true values. Comparison of predicted TL computed by inverted parameters with the experimental TL can be used to demonstrate the reliability of inversion results, as shown in Figure 9. The received depths (RD) at the array are 37 and 77 m, the frequencies are 80, 160 and 500 Hz. The predicted TL is computed using the inverted geoacoustic parameters at 65 min. For the experimental TL, each frequency is calculated as the 1/3 octave average, the same as when computing the predicted TL. It is clear that the predicted TL and experimental TL are consistent for most distances. The root mean square error (RMSE) between predicted TL and experimental TL is less than 3 dB.

Discussion
The selection of P 0 and Q k has an important influence on filter convergence. Much research has provided varying combinations of P 0 and Q k in different tracking or filtering questions. In general, the design of these values needs to consider variance from k − 1 to k and the range of value of corresponding parameters. It is hard to choose an appropriate value when a parameter first changes slowly then rapidly. In addition, the range-dependence of real seafloor parameters needs to be considered carefully as well. For our work, though parameters of the bottom change spatially. The slope is still small (water depth is about 107 m at 0 km and 103 m at 30 km). In such a situation, the P 0 and Q k is relatively small to make filter perform stable. If the correlation between spatial/time with parameters increases and the structure of bottom changes significantly, the filter must be modified such as increasing the length of ensemble and particles, or adjusting Q k . It would be worked but increase the computational cost and have greater risk of divergence.
Recently, an adaptive particle filter has been developed to try to solve this problem and was applied to geoacoustic inversion of a horizontal towed array [35]. Compared to typical sequential filters, it uses an adaptive method to update Q k online rather than an assumed diagonal matrix, which ensures tracking performance in the time-varying environment. When dealing with strong range-dependent inversion problems such as continental shelf break [23], the adaptive method may perform better.

Conclusions
In this paper, the sequential inversion of geoacoustic parameters and the tracking of acoustic source locations are achieved using measured acoustic fields in the VLA with ensemble Kalman Particle filter. The simulated environment and East China Sea shallow-water experiment were used to test the efficiency of inversion schemes.
Sequential inversion and source tracking in underwater environments is a complex and non-linear problem that can be handled by sequential Bayesian filters but PF and EnKF have defects such as "blind" prediction and Gaussian assumption for PDFs. EnKPF combines the PF and EnKF to overcome their drawbacks because EnKPF employs the PPD of EnKF as proposal distribution of the PF to improve the performance of the PF. For this reason, the non-linearity between environmental parameters and measured acoustic field is typically encountered. Simulation results showed that the EnKPF outperformed the PF and the EnKF over the whole tracking period. For most of the geoacoustic parameters in sea bottom and sound source parameters, EnKPF showed good tracking performance. In addition, air-gun data from a shallow-ocean experiments was used to invert geoacoustic parameters with source location by using EnKPF. Sequential inversion results showed that source parameters and water depth are consistent with measurement values with small error. Furthermore, sediment sound speed was consistent with the in situ measurement. The RMSE between predicted and experimental TL was less than 3 dB, and the predicted TL was computed by environmental model with parameters inverted by EnKPF.
It should be pointed out that EnKPF needs more computational costs to achieve the best performance because its computational complexity is increased. If the system dimension increases, the filter performance and the running time is always a matter of trade-offs. This method needs further exploration.