Joint Inversion for Sound Speed Field and Moving Source Localization in Shallow Water

This paper develops a joint approach for time-evolving sound speed field (SSF) inversion and moving source localization in shallow water environment. The SSF is parameterized in terms of the first three empirical orthogonal function (EOF) coefficients. The approach treats both first three EOF coefficients and source parameters (e.g., source depth, range and speed) as state vectors of evolving with time, and a measurement vector that incorporates acoustic information via a vertical line array (VLA), and then the inversion problem is formulated in a state-space model. The processors of the extended Kalman filter (EKF) and ensemble Kalman filter (EnKF) are used to estimate the evolution of those six parameters. Simulation results verify the proposed approach, which enable it to invert the SSF and locate the moving source simultaneously. The root-mean-square-error (RMSE) is employed to evaluate the effectiveness of this proposed approach. The interfile comparison shows that the EnKF outperform the EKF. For the EnKF, the robustness of the approach under the sparse vertical array configuration is verified. Moreover, the impact of the source-VLA deployment on the estimation is also concerned.


Introduction
Ocean parameter inversion and source localization are practical problems with serious concerned and long standing applications in underwater acoustics [1,2]. A common feature of these problems is that underlying environment parameters are estimated from measured acoustic data. The sound speed profile (SSP), one of the important environment parameters extracts from ocean properties, has a significant influence on determining acoustic waveguide propagation [3]. Oceanographic variability and dynamic processes cause the SSP evolve in time and space, especially in shallow water areas [4,5]. The SSP can be obtained by combining model outputs or in situ measurement, however, the temporal and spatial sampling from a large area within a short time and long-term observation is often unrealistic. Numerous approaches have been developed to invert the sound speed in the water volume fast and accurate. To reduce the degrees of freedom, ocean parameters include the SSP are often represented in terms of EOFs [6,7]. The previous studies indicate that the inversion of SSP has a good agreement with the measured data, in both shallow [8] and deep ocean [9]. Besides, the compressive estimation of SSP using EOFs is demonstrated in reference [10]. These characteristics of the spatiotemporal variability can be transformed to an estimating problem.
Since the introduction of the Kalman filter (KF) [11], much more researches on KF family and their extensions have been implemented for practical issues. These filters are suitable approaches to estimate the nonlinear and non-Gaussian problems, and have been used previously in underwater acoustic application [12][13][14]. An excellent overview described that sequential 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 [15]. Yardim et al.

A State-Space Model for Joint Inversion
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. In this section, theory of SSF parameterization is summarized briefly, and then the mathematical model of the moving source tracking is given.

Sound Speed Field Parameterization
EOF is a reasonable way to fully parameterize the SSP [22], which can be obtained from a direct measurement of database and are efficient in reducing the degrees of freedom of estimated parameter vector.
Assuming that there are N SSP measurements c 1 (z i ), c 2 (z i ), · · · c N (z i ) i = 1, 2, · · · , M, where M is the number of discrete points in depth. The mean sound speed is given by c(z) = 1 N N i=1 c i (z), and the sound speed perturbation is expressed as: where c i z j denotes the sound speed value at depth j of the i-th SSP. The sample covariance matrix R is computed to characterize the uncertainty of the sound speed: where T is the transpose. Using the orthogonal decomposition, we have R = M n=1 λ n f n f T n , where λ 1 > λ 2 > · · · > λ M are the respective eigenvalues of eigenvectors f n . Due to the eigenvalues diminishing rapidly, only a few EOFs are able to dominate the SSP characterization. Therefore, each SSP in the ocean can be approximated with the first K eigenvectors as: where α k and f k (z) are called the kth EOF coefficient and EOF, respectively. As stated above, inversion for SSF is a multi-dimensional optimization problem. The computational complexity and inversion accuracy usually depend on the number of inverted variables. Therefore, a limited number of EOFs are used to invert the SSF. Assuming the number of EOFs to be retained K is based on the amount of total energy contained in the associated eigenvalues { λ i } . Defining the ratio of partial-total energy as a percentage: where N is the total number of eigenvalues, which enable us to select K. In order to demonstrate the effectiveness of inversion method, data measured in the 2001 ASIAEX, ECS experiment were used [23]. The sound speeds were obtained from the onsite conductivity-temperature-depth (CTD) measurements. They were sampled at 106 discrete depth points from the surface to the bottom, equally spaced by 1 m. Figure 1a shows 54 samples of sound speed perturbations, which were calculated by the SSP data measurements, the fluctuation amplitude was more than 11 m/s between 20 m and 60 m in depth. The variability of SSP was modeled using EOF by calculating from Equations (2) and (4), and the first three EOFs are plotted in Figure 1b. In this case, Figure 1c shows that choosing K = 3 could represent approximately 90% of the energy disturbance. Then, the first three EOFs were chosen to invert the pth SSP in 54 measured samples, the EOF coefficient was calculated as , it is seen that the inverted SSP was in good agreement with the measured SSP in Figure 1d. The depth-integrated RMSE was employed to evaluate the discrepancy between the measured SSP and inverted SSP. After calculation, the RMSE of inverted SSP by EOF was 0.37 m/s. As the result, the first three EOFs were capable of represent SSP. As introduced above, the inversion problem of SSP corresponds to the projection of first three EOF coefficients. The proposed scheme developed in the paper was to reconstruct the time-evolving SSF. The SSF evolution of vertical slice was then discretized into a finite number of adjacent rectangular intervals, each of which describes its characteristics with an independent SSP. Dividing the evolution into L intervals and using the first three EOFs. For the th l interval, the sound speed was calculated in the form of: where , k l α is the th k EOF coefficient of the th l profile. In this way, any sound speed estimates should be constrained to lie in the evolution of the EOFs. That is, the SSF is recovered by a total of l SSPs, and the coefficients , k l α constitute the SSF parameters.

The Mathematical Model of the Joint Inversion
The joint inversion for SSF and source localization requires two dynamic equations: A state equation that models the evolution of SSF and the movement of source, and a measurement equation that relates the environment and source location at step k to the acoustic field. As introduced above, the inversion problem of SSP corresponds to the projection of first three EOF coefficients. The proposed scheme developed in the paper was to reconstruct the time-evolving SSF. The SSF evolution of vertical slice was then discretized into a finite number of adjacent rectangular intervals, each of which describes its characteristics with an independent SSP. Dividing the evolution into L intervals and using the first three EOFs. For the lth interval, the sound speed was calculated in the form of: where α k,l is the kth EOF coefficient of the lth profile. In this way, any sound speed estimates should be constrained to lie in the evolution of the EOFs. That is, the SSF is recovered by a total of l SSPs, and the coefficients α k,l constitute the SSF parameters.

The Mathematical Model of the Joint Inversion
The joint inversion for SSF and source localization requires two dynamic equations: A state equation that models the evolution of SSF and the movement of source, and a measurement equation that relates the environment and source location at step k to the acoustic field.
It is known that ocean sound speed can be parameterized in terms of the first three EOF coefficients in Section 2.1. In this work, we have: w k is the state noise vector. For a moving source, the depth Z, horizontal range R, and radial speed v are to be estimated: where s denotes the source parameters, ∆t is the update time increment (i.e., the time interval between two measurements), and v d and v a are the random variables representing the variation in source depth and acceleration, respectively. Then the state vector that includes the parameters of source and SSF is given by merging these two blocks. Defining the state vector as x T k = s T k α T k . The state and measurement equations are considered as: x where f(·) is a state transition function of the state vector x k−1 , and h(·) is the normal mode propagation model, and v k and w k are the process/state and measurement noise vectors, respectively. Here, y k is the acoustic pressure measured at kth time frame, the measurement equation was used to characterize the signal of acoustic pressure, which was received by a VLA. According to the normal mode theory, the synthetic acoustic pressure can be expressed as a general nonlinear function of the horizontal range R, depth Z and sound speed c.
As mentioned above, the state equation describes the evolution of the SSF and source location, 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 estimate the evolving state.

Nonlinear Kalman Filter as an Inversion Processor
In dynamical evolution of time-evolving SSF and moving source localization, the processor of the sequential method is used to predict and estimate the joint inversion results. For the nonlinearity of measurement function h(·), two kinds of nonlinear filtering algorithms are given as follows.

EKF Model
A common Kalman filter applied to the nonlinear systems was EKF. Since KF is an exact linear tracking filter, its extended version works by converting the system into a form over which the KF can operate. In an EKF, this is done by locally linearizing the state and measurement equations using the first order Taylor series expansion of the nonlinear transformations (such as the normal mode propagation model h(·)). After linearization, the state and measurement equations given in (17) and (18) to: x where F k and H k are Jacobian matrices of f and h, respectively. Table 1 summarizes the principle of EKF, applied to the proposed state-space model. In order for EKF to perform well, nonlinearities should be small and the underlying densities should also be close to Gaussian.

EnKF Model
The EnKF is a KF that includes Monte Carlo sampling, which enables us to deal with moderate nonlinear and non-Gaussian problems. The basic principle is initializing a set of system state samples as a background ensemble, and then using the latest observation information to update each background ensemble member in current time by Kalman filter to obtain the analysis ensemble. The sampling ensemble is transferred by the system model to get the background ensemble of the next time. Each member is propagated through the state transition and measurement functions and the statistics of the estimates are propagated without the use of any linearization step. An EnKF algorithm is summarized in Table 2.

Principles
Definitions X k = x k,i , i= 1, 2, · · · , N Background ensemble member: Where i is the ensemble member index and k is the time frame x k Background ensemble mean value at k X a k = x k,i a , i= 1, 2, · · · , N Analysis ensemble member x a k Analysis ensemble mean value at k z k,i = y k + ε i Ensemble of observation Predict: Kalman gain:

Estimation of EKF and EnKF for SSF and Moving Source Localization
To synthesize the acoustic pressure data used in the measurement equation, the simulated environment is shown in Figure 2. The water depth was 106 m according to the measurement. The bottom was modeled as a homogeneous fluid half space, which consisted of the silt clay sediment (sound speed of 1610 m/s, density of 1.7 g/cm 3 , and attenuation rate of 0.1dB/λ). These values of the geoacoustic parameters correspond to average values deduced from the core measurements of AXIAEX, ECS, 2001 [23]. The VLA was configured with 16-element equispaced 4 m that spanned the 15-75 m water depths. The source was placed at 40 m in depth and 1 km in range with a frequency of 400 Hz. Moreover, the source had 180 steps (k = 1, · · · , 180) movements away from VLA with the speed of 2.5 m/s. As a result, the continuous SSF was discretized in 180 intervals with the corresponding SSPs.  For the EnKF, the ensemble number is related to the precision and computational complexity. In order to balance the trade-off between the computational complexity and estimation accuracy, the ensemble number was set to = 50 N . We assumed that both additive noise terms k v (in Equation (8)) and k w (in Equation (9)) were represented by the Gaussian probability density functions. The adiabatic normal mode propagation model was used to compute the acoustic field by inputting the environmental and source parameters (except source speed) into KRAKENC [24]. As the state-space model described in Section 2.2, the first three EOFs coefficients and source parameters (depth, range, speed) were taken as the state vectors, and the acoustic pressures received by VLA were served as the measurement vectors. Under the proposed joint approach, the evolutions of EOFs coefficients were estimated to reconstruct the SSF, and the source parameters were tracking over time to give the localization. That is, this approach needs to estimate the evolution of those six parameters simultaneously. The evolution of the six parameters is presented in Figure 3, where the estimation results were compared with the true values. We found that two filters were able to estimate the parameters. For the EKF, the processor shows good tracking performance at the beginning. However, there were divergences between the estimation results and true trajectories over time. It was seen from the curves of the first three EOFs coefficients in Figure 3a and source depth in Figure 3b, the large divergences appeared around the 55-60th min. As it describes, we found that the estimation of source depth would be affected when the EOF coefficients were not tracked well. The one reason for this phenomenon was explained that the acoustic pressure field was more sensitive to the depth of the source in normal mode theory. The other one was related to the basic assumption of the random walk in the state-space model. Noting that there were rapid perturbations of sound speed in some regions of the true environment, the above assumption might be explained that the perturbations result in For the EnKF, the ensemble number is related to the precision and computational complexity. In order to balance the trade-off between the computational complexity and estimation accuracy, the ensemble number was set to N= 50. We assumed that both additive noise terms v k (in Equation (8)) and w k (in Equation (9)) were represented by the Gaussian probability density functions. The adiabatic normal mode propagation model was used to compute the acoustic field by inputting the environmental and source parameters (except source speed) into KRAKENC [24]. As the state-space model described in Section 2.2, the first three EOFs coefficients and source parameters (depth, range, speed) were taken as the state vectors, and the acoustic pressures received by VLA were served as the measurement vectors. Under the proposed joint approach, the evolutions of EOFs coefficients were estimated to reconstruct the SSF, and the source parameters were tracking over time to give the localization. That is, this approach needs to estimate the evolution of those six parameters simultaneously.
The evolution of the six parameters is presented in Figure 3, where the estimation results were compared with the true values. We found that two filters were able to estimate the parameters. For the EKF, the processor shows good tracking performance at the beginning. However, there were divergences between the estimation results and true trajectories over time. It was seen from the curves of the first three EOFs coefficients in Figure 3a and source depth in Figure 3b, the large divergences appeared around the 55-60th min. As it describes, we found that the estimation of source depth would be affected when the EOF coefficients were not tracked well. The one reason for this phenomenon was explained that the acoustic pressure field was more sensitive to the depth of the source in normal mode theory. The other one was related to the basic assumption of the random walk in the state-space model. Noting that there were rapid perturbations of sound speed in some regions of the true environment, the above assumption might be explained that the perturbations result in larger variations of the EOF coefficients, which exceed the statistical values assumed for the random walk, thus resulting in a model mismatch. In addition, Figure 3a shows tracking of the EOF3 was worse than EOF1 and EOF2 coefficients. The possible reason was attributed to the third EOF playing a minor position compared with the first two orders on dominating the SSP perturbation, which leads to a minor variation of the sound speed in Equation (3). In the measurement equation, a little pressure variation information will affect the processor to track. By comparison, the estimation results confirmed that the EnKF was superior to the EKF. The EnKF was capable of estimating all parameters with a stable and narrower fluctuant result, which almost coincided with the true value completely.  According to the above estimation results of EOF coefficients, both two filters enabled us to reconstruct the true SSF, the inversion of SSF is shown in Figure 4. Figure 4b,d show pseudo-color maps of the relative error between estimated SSF and true SSF. Moreover, the errors near the thermocline were larger than the other regions. The inverted SSF of an EKF as Figure 4a,b indicates that the processor was correctly tuned at the beginning of estimation, but the sound speed mismatched because of the disturbance appeared around the 55-60th min, that was corresponding with the evolution of EOFs coefficients in Figure 3a. From Figure 4c,d, the EnKF shows a steady estimated result, especially in the region where intensive perturbation of sound speed existed. Meanwhile, the relative errors stay near-zero over the whole evolution period. According to the above estimation results of EOF coefficients, both two filters enabled us to reconstruct the true SSF, the inversion of SSF is shown in Figure 4. Figure 4b,d show pseudo-color maps of the relative error between estimated SSF and true SSF. Moreover, the errors near the thermocline were larger than the other regions. The inverted SSF of an EKF as Figure 4a,b indicates that the processor was correctly tuned at the beginning of estimation, but the sound speed mismatched because of the disturbance appeared around the 55-60th min, that was corresponding with the evolution of EOFs coefficients in Figure 3a. From Figure 4c,d, the EnKF shows a steady estimated result, especially in the region where intensive perturbation of sound speed existed. Meanwhile, the relative errors stay near-zero over the whole evolution period.
For both two filters, the inversion performance of SSF and source parameters was evaluated by root-mean-square-errors (RMSE) between the true values and estimated results. The RMSE curves of SSF and source parameters were calculated as shown in Figure 5, for comparison. We found that both processors still converged during an hour, but the RMSE values of EKF were much larger than EnKF. The values of SSF were especially enlarged (the maximum RMSE value was close to 7 m/s for the EKF) around the 55-60th min, which could be predicted by the tracking results in Figure 4. After calculation, the time-average RMSE values of SSF, source depth, range and speed for EKF were 0.6675 m/s, 0.0134 m/s, 0.8736 m/s and 0.0013 m/s, respectively. The EnKF, in contrast, maintained a low stable value of RMSE, which time-average values remained near-zero clearly. Therefore, it could be concluded that the EnKF outperformed EKF of this joint approach. reconstruct the true SSF, the inversion of SSF is shown in Figure 4. Figure 4b,d show pseudo-color maps of the relative error between estimated SSF and true SSF. Moreover, the errors near the thermocline were larger than the other regions. The inverted SSF of an EKF as Figure 4a,b indicates that the processor was correctly tuned at the beginning of estimation, but the sound speed mismatched because of the disturbance appeared around the 55-60th min, that was corresponding with the evolution of EOFs coefficients in Figure 3a. From Figure 4c,d, the EnKF shows a steady estimated result, especially in the region where intensive perturbation of sound speed existed. Meanwhile, the relative errors stay near-zero over the whole evolution period.  For both two filters, the inversion performance of SSF and source parameters was evaluated by root-mean-square-errors (RMSE) between the true values and estimated results. The RMSE curves of SSF and source parameters were calculated as shown in Figure 5, for comparison. We found that both processors still converged during an hour, but the RMSE values of EKF were much larger than EnKF. The values of SSF were especially enlarged (the maximum RMSE value was close to 7 m/s for the EKF) around the 55-60th min, which could be predicted by the tracking results in Figure 4. After calculation, the time-average RMSE values of SSF, source depth, range and speed for EKF were 0.6675 m/s, 0.0134 m/s, 0.8736 m/s and 0.0013 m/s, respectively. The EnKF, in contrast, maintained a low stable value of RMSE, which time-average values remained near-zero clearly. Therefore, it could be concluded that the EnKF outperformed EKF of this joint approach.
The main reason for this different estimated capacity of two filters can be summarized as follows: In an EKF, the covariance matrices propagated through the nonlinear model with the Jacobian of the transition and measurement functions. Such a first-order approximation, which ignores the higherorder terms, could result in large estimation errors, which could cause the filter to diverge (the gain tends to zero). In the inversion scheme, high nonlinearity of the measurement function requires highdimensional nonlinear methods. While, the EnKF represents the underlying distributions with a large ensemble of model realizations, and the statistics of estimates are propagated without the use of any linearization step. It is a method for nonlinear transformation of the mean value and covariance to be practical for the realistic high-dimensional system [25].

Sparse Array Configuration
In the above simulation, the initial configuration was a 16-element VLA. However, it is of practical interest to study whether the joint approach was robust when using a few hydrophones or a shorter array. The EnKF of better performance was chosen as a processor. The different configurations tested are given in Table 3. Keeping the other parameters settings consistent with the Figure 2, and reducing the number of VLA to 8-element in Test B (few hydrophones) and C (shorter The main reason for this different estimated capacity of two filters can be summarized as follows: In an EKF, the covariance matrices propagated through the nonlinear model with the Jacobian of the transition and measurement functions. Such a first-order approximation, which ignores the higher-order terms, could result in large estimation errors, which could cause the filter to diverge (the gain tends to zero). In the inversion scheme, high nonlinearity of the measurement function requires high-dimensional nonlinear methods. While, the EnKF represents the underlying distributions with a large ensemble of model realizations, and the statistics of estimates are propagated without the use of any linearization step. It is a method for nonlinear transformation of the mean value and covariance to be practical for the realistic high-dimensional system [25].

Sparse Array Configuration
In the above simulation, the initial configuration was a 16-element VLA. However, it is of practical interest to study whether the joint approach was robust when using a few hydrophones or a shorter array. The EnKF of better performance was chosen as a processor. The different configurations tested are given in Table 3. Keeping the other parameters settings consistent with the Figure 2, and reducing the number of VLA to 8-element in Test B (few hydrophones) and C (shorter array) to compare with Test A. The evolutions of six parameters within an hour were estimated for three different tests, and the inversion for SSF under sparse array configuration is shown in Figure 6. The existence of fast-changing sound speed near the thermocline made a larger disturbance for EOF coefficients, which led to increased error. The inversion performances are summarized by the time-evolving RMSE in Figure 7. Although the estimation of EnKF was always converged, the RMSEs of Test B were larger than the other tests, seen from the curves of both SSF and source parameters. Moreover, we found the values of RMSE were enlarged around 40 and 55th min, that was due to the EOF coefficients (i.e., SSP) that were impacted with sound speed perturbation, especially in Test B and C. The evolutions of six parameters within an hour were estimated for three different tests, and the inversion for SSF under sparse array configuration is shown in Figure 6. The existence of fastchanging sound speed near the thermocline made a larger disturbance for EOF coefficients, which led to increased error. The inversion performances are summarized by the time-evolving RMSE in Figure 7. Although the estimation of EnKF was always converged, the RMSEs of Test B were larger than the other tests, seen from the curves of both SSF and source parameters. Moreover, we found the values of RMSE were enlarged around 40 and 55th min, that was due to the EOF coefficients (i.e., SSP) that were impacted with sound speed perturbation, especially in Test B and C.    . Figure 7. RMSE of the EnKF estimation results obtained with different array configurations (in Table  3 Table 4 gives the time-averaged RMSEs of the three tests. As expected, the joint inversion performance mainly depended on the number of effective measurements. For a same given number  Table 3 Table 4 gives the time-averaged RMSEs of the three tests. As expected, the joint inversion performance mainly depended on the number of effective measurements. For a same given number of measurements, the approach was more suitable for the case involving a smaller spacing in depth (see Test B versus C). Test B gave the least accurate estimation, while Test A, characterized by the maximum number of measurements, outperformed the other configurations. Compared with 16-VLA, the RMSEs of the 8-VLA in the same time dimension increased slightly. Such differences might be explained that both the sparse elements and large spacing reduced the rate of spatial sampling, which led to an increasing error. To sum up, EnKF was a superior processor to invert the SSF and locate the moving source. Simulation results show the approach was robust even under a short array or lack of hydrophone elements.

Impact of the Source-VLA Deployment
With the gradual development of underwater vehicles in practical applications, the position of the source will be more flexible. As it is known from the normal mode theory, the source in different depths excites relative modes to represent the acoustic field, which may have an impact on the performance of the joint approach. The influence of the source-receiver deployment on the estimation is necessary to concern. To evaluate the effective of the joint approach, the source was placed at 10/96 m in depth, and moved away from the VLA in a horizontal direction. The VLA with 16 and 8 elements were deployed spanning upper/middle/bottom part of the water, respectively. After simulation, the time-average RMSEs of EnKF for SSF, source depth, range and speed estimation are listed in Table 5. For both two depths of source, the excellent estimation results were obtained when the VLAs were deployed in the middle depth of the sea, respectively. It was found that the source-VLA was placed in a cross-laying configuration, which had a better effect than the relative horizontal placement. Such differences in filter performance might be explained in terms of the optimal position of the hydrophones with respect to the number of the latter used in the VLA. Furthermore, the cross-laying source-VLA deployment would get a better collection of the acoustic measurements information with a higher resolution and larger sampling measurement.

Conclusions
In this paper, a joint approach of SSF inversion and moving source localization was proposed. The flexibility of the moving source could increase the sampling density of acoustic pressure field, and enrich the information of environmental parameters in the ocean. The parameters of SSF (e.g., first three EOF coefficients) and source (e.g., source depth, range and speed) were inverted simultaneously by assimilating measurements of the acoustic filed on a VLA into the processors of EKF and EnKF. The inversion scheme was verified by the realistic environment simulated from the ASIAEX ECS 2001 experiment measured data. The RMSEs were used to evaluate the approach in this work, and the comparison showed that the EnKF outperformed the EKF. For the EnKF, the robustness of the proposed approach was verified under the sparse array configuration. Moreover, the impact of the source-VLA deployment on the estimation was also concerned.
Source localization and geoacoustic parameters inversion are still hot issues in the research field of underwater acoustics. The method we proposed was useful for studying the inversion of dynamic ocean parameters and even providing a reference on the deployment of underwater observing system. A joint approach allows real-time updating of the SSF and accurate locating of the source, which could improve the accuracy than the traditional ocean acoustic tomography (OAT) method in a large scale, and provide high-resolution measurement of the ocean environment. The theory presented in this paper was based on the motion of source, which had no effect on the surrounding environment assumption. For the next work, radiated noises of the moving source should be taken into account, incorporating the Doppler shift effects to better match the measurement equation, an acoustic model will be derived based on the normal mode propagation and waveguide Doppler.