Proceedings of the Institute of Acoustics

 

Inverting the sound speed profile from multi-beam echo sounder data and historical measurements – a simulation study

 

Y. Gourret, University of South-Eastern Norway, Horten, Norway
T. Brander, University of South-Eastern Norway, Horten, Norway
K. T. Hjelmervik, University of South-Eastern Norway, Horten, Norway

 

 

1 INTRODUCTION

 

The ocean’s lack of transparency offers plenty of opportunities for illegal trade enterprises¹, terrorism, and covert military operations. Unmanned surface vehicles and autonomous underwater vehicles are possible disruptive technologies that might counter these threats². The lack of human decision makers puts strong requirements on the autonomy of the vessels. Also, the inherent compactness of these vessels imposes strong restrictions on the amount and size of installed hardware. 

 

For manned sonar operations, sound speed probes are used to determine the local sound speed profile, which is then used to estimate the sonar performance. For unmanned vessels, the extra payload for sound speed probes should be avoided. Ocean models provide an alternative source of information on the underwater environment. However, their predictions are still unable to provide the level of detail and accuracy required for small-scale sonar operations in littoral environments. Climatological data³ and other statistically derived oceanographic data^4,5  may also provide information on the present sound speed, but fail to account for local spatio-temporal effects. 

 

Through-the-sensor techniques is an alternative means of assessing the environment. Sound speed inversion from multi-beam echo sounder (MBES) data is addressed in several studies where the main goal is to improve topographic measurements^6,7,8,9. Common for all these studies is that both the sound speed and the topography are assumed unknown in the inversion, and that they focus on the accuracy of the topographic measurements, rather than the quality of the sound speed profile (SSP). Hjelmervik¹⁰  proposed inverting the sound speed from single-beam echo sounder (SBES) data and direct measurements of the sound speed at limited depths in an area with known topography. A major limitation in this work is that an SBES yields information on the average SSP only, not its structure. Keyzer et al.⁹  inverts the SSP from MBES data in a similar manner using an empirical orthogonal function (EOF) representation for the SSPs. The advantage of using MBES data is that beams directed off the vertical will be refracted and thus contain information on the sound speed gradients throughout the depths. However, the cost function used by Keyzer et al. does not include any constraints that ensure that the resulting profile is physical. 

 

We invert the SSP from MBES data in an environment with known topography. The SSPs are represented using EOFs, such that the proposed SSPs can be modelled by a limited number of coefficients¹⁰. An acoustic raytracer models the travel time from the transmitter to the bottom and back to the receiver for proposed SSPs. The cost function compares the modelled and measured travel times for all beam angles, but also includes a priori knowledge on the SSPs through Tikhonov-type regularization¹¹. We assess the validity and performance of our method through simulations, e.g. to explore its sensitivity to errors in the measured travel times. Our work thereby expands on former work^10,9, as we include both MBES data to capture the vertical structure in the sound speed profile as well as constraining the EOF space explored by the method to physically feasible regions through the use of regularization. We also show that the problem is non-unique¹¹ without regularization.

 

2 THEORY

 

Recovering the sound speed is a travel time tomography inverse problem¹², since the travel time depends on the SSP¹³. Each measurement yields a set of travel times corresponding to different beams processed from the received data on the MBES. We assume range-independent sound speed. 

 

As a simple example, consider a piecewise constant SSP, where the ith vertical layer has a constant sound speed of c_i  and a depth of ∆Z_i . Consider a receiving beam directed with an initial angle θ₁ to the vertical at the receiver. Assume a flat bottom with a seafloor depth z_b. The arrival time of a signal arriving in this beam following a reflection off the sea floor, is then given by

 

 

where σ_i is the length of the trajectory through the ith layer. The latter expression follows from the trajectory following a straight path due to the constant sound speed in each layer. By taking into account Snell-Descartes law we then have

 

 

The expression does not depend on the order of the vertical layers. Thus two layers may be switched without impacting the travel times. Hence, the mapping from SSP to travel time is not injective, since there exist several different SSPs that give the same travel times. This ambiguity causes the non-uniqueness of the inverse problem. 

 

Since we have plenty of physical information about the SSP, we can use a priori information to deal with the lack of uniqueness. Not all vertical sound speed structures are physical even if they provide solutions to the inverse problem, and even if they are physically possible, they will have different associated probabilities of occurring. 

 

To structure the historical sound speed measurements, the SSPs are represented by a weighted sum of EOFs following the steps of Preisendorfer¹⁴. The EOFs are determined from a set of historical sound speed profiles from an area of limited geographical extent and a limited part of the year. 

 

We use the vector x = [x₁, x₂, ..., x_K]^T, where x_k  are the EOF coefficients, to write the SSP as

 

 

The average SSP, ĉ, may be estimated from historical SSPs, which are first interpolated to a common depth grid with K depth steps. The matrix U  then becomes a K×K  matrix containing the EOFs estimated from the historical SSPs. 

 

The coefficient vector x_n of the nth historical sound speed profile c_n may be found by

 

 

Given a large amount of historical measurements and statistical independence between different coefficients, the probability density functions of each coefficient, f_xk (x), may be derived. Assuming Gaussian behaviour, the distribution for the EOF coefficients x then becomes

 

 

where σ_k is the standard deviation of the kth EOF coefficient and may be estimated from the historical dataset.

 

The number of EOFs corresponds to the number of depth steps K. Most of the sound speed variance is captured by the first few EOFs. The lower order of EOFs represent the bulk of the information¹⁴. We therefore choose N_EOF  as a cutoff for the number of EOFs used and thereby reduce the matrix U to a size of K × N_EOF  and the vector x will only have N_EOF  elements.

 

Assumptions

 

We consider an environment with a flat bottom of known depth. The position and altitude of the MBES is also known. The MBES records N_beam infinitely narrow beams, and the beam directions are perfectly known. The error on travel time is independent and Gaussian distributed with standard deviation σ_τ . This error accounts for the measurement error at the receiver due to the imperfect detection of a received signal embedded in noise, but also for imperfect knowledge of the seafloor depth.

 

In this simulation study we limit the error source to an error in the received travel time. In reality, the results will also be influenced by a plethora of other sources of error for different parameters, such as angular errors and uncertainty in the forward model, but this is out of scope for this study.

 

2.1 Cost Function

 

Let t(x) be a function of the EOF coefficients x = [x₁, x₂, ..., x_K]^T. This corresponds to the modelled travel times which equation (1) is a simplified model for. Here we employ an acoustic raytracer to estimate the travel time for increased accuracy. 

 

The proposed cost function minimizes the difference between measured (T_measure) and the modelled (t (x)) travel times. The cost function also takes into account the historical information by including a regularization term containing the log-likelihood function for the EOF coefficients x. We obtain

 

 

To solve the minimization problem, we use the Gauss-Newton method. The problem is non-linear, but may empirically be solved for small variations. The mapping x 7→ t(x) is continuous for a flat bottom and the angles and SSPs where the rays are not refracted back to the surface, but the differentiability assumed by Gauss-Newton is unknown. The validity of the proposed method is limited to flat or slowly varying seafloors. For complex seafloors small variations in the ray trajectory may result in large variations in travel time. 

 

The value α > 0 is a Tikhonov-type regularization parameter¹¹, that ensures a balance between historical sound speed measurements and the travel time measurements. A large α gives more weight to the historical information, while a low value gives more weight to the acoustic measurements. Generally, selection of the regularization parameter is a difficult task¹¹. Here we use a simple neural network (NN) to determine a good choice for α. 

 

To work, a NN should be trained with many examples. The training dataset used to generate the EOFs may be too small for this training. However, by randomizing the coefficients x_i according to the Gaussian distribution, we can recreate synthetic SSPs cf. eq. (3). This allows the generation of arbitrary amounts of synthetic examples to train our NN. 

 

Figure 1 gives an overview of how the neural network works. We invert the SSP for N_inv different values of α (on a logarithm scale). For each inversion, we get travel time errors and EOF coefficients x, that are normalized respectively by the formula (a) and (b) shown in Fig. 1. We define A_i as the normalized information coming from the inversion, cf. formula (c) in Fig. 1. The neural network is a multilayer perceptron trained to output the SSP error for α_i  with list A_[i−k,i+k]  as input, with k as the size of the operative window. This operation will be repeated for every value of α. The NN requires a fixed number of inputs, so at boundaries we duplicate A_i , if i /∈ [1, N_inv] then A_i = A₁ or A_Ninv . Even if not indicated on the figure, the number of beams is also an input of the NN. 

 

The final α used for the inversion is the α_i  that gives the lowest response from the NN.

 

 

Figure 1: An overview of the neural network. The blue border represents the operative window.

 

3 DATASET

 

The MBES dataset is generated through ray-based simulations. The simulations employ a large set of sound speed profiles measured in the Norwegian trench on the West coast of Norway from 1990 to 2010. The sound speed measurements were taken from a validated dataset¹⁵.

 

We divided the data into two sub-datasets; a training and a testing dataset. Both are geographically cropped according to a bounding box defined by [59.849^◦, 62.092^◦] in latitude and [2.924^◦, 4.990^◦] in longitude (South-West coast of Norway). The profiles are also cropped to a depth of 300 m. The training data (151 SSPs) were measured in April in the years from 1990 and 2000, while the testing data (111 SSPs) were measured in April in the years between 2001 and 2010. The objective is to invert the test SSPs based on historical information provided by the training data.

 

 

Figure 2: Left: Dataset area in orange, training data in blue, testing data in red. Right: Training dataset information, average SSP and 5 first EOFs.

 

4 RESULTS AND DISCUSSION

 

The proposed method was applied on a simulated raytrace based on the test dataset. Figure 3 presents the statistics for the errors of the inverted SSPs along with a single representative example of an inverted SSP. The main bulk of errors lies around 0.7 – 0.9 m/s on average over the entire depth. By comparison, the mean SSP error between the test dataset and the average SSP ĉ (see eq. 3) is 3.22 m/s.

 

Tab. 1 shows the results of a sensitivity analysis for different parameters. The central column (green) lists the parameters used to produce the results shown in Fig. 3. The table presents the average RMS error of the inverted SSPs of the testing dataset for different parameters. This permits us to have an overview of the impact of every parameter. For comparison the average error obtained from simply comparing the average SSP in the training data set to all the test SSPs is 3.22 m/s, and 2.55 m/s when using the average SSP from the test data set. Lacking the historical measurements, a natural choice of SSP to use during a sonar operation could be a climatological estimate from for instance the World Ocean Atlas¹⁶. This would result in an error of 2.62 m/s, which is still twice as high as the maximum error observed in the parameter variations shown in Tab. 1.

 

 

Table 1: RMS error for different variations of parameters used in the overall method. The central column, in green, corresponds to the results shown in fig. 3. The spatial error σ_x is here given by σ_x = c _t / 2 .

 

The first parameter considered is the amount of beams times the amount of pings. This is a combination of the amount of beams along the swath and the amount of pings accumulated before inverting the sound speed. This also includes extra measurements accumulated by exploiting that the bottom might be detected at different time samples within the same beam. However, it represents the number of independent measurements. An increase here results in a decreasing SSP RMS error. This is expected since an increased amount of travel time measurements should result in an improved estimate of the SSP, in much the same way as repeated, independent measurements of a random value will increase the accuracy of an estimate for that value.

 

The second parameter is the swath angle width of the MBES. For narrow swath angle, the sonar measurements are all directly below the ship. This only provides information on the average SSP, as a wave travelling directly downwards in a range-independent environment remains unrefracted and has a travel time dependent on the average sound speed and seafloor depth only. By increasing the opening angle, we increase the incidence angle of the sound wave to the water layers at the side of the MBES. This increases the refraction effects providing us more information for the inversion of the SSP. As we observe on Tab. 1, the accuracy of the inversion increases with the opening angle. The achieved gain from widening the opening angle is limited to the spectrum of angles that result in bottom interaction. Depending on the sound speed profile, angles close to the horizontal can be reflected by water layers, and thus there might be no bottom interaction.

 

The third parameter is the number of EOFs used for the inversion. As we can see, the SSP error decreases until six EOFs are used, but for seven EOFs, the SSP error increases. Thus, increasing the number of EOFs further does not increase the quality of the inversion. The number of EOFs should be chosen carefully. The quality of the EOFs depends on the number of SSPs used to obtain them. Higher- order EOFs require more SSPs to be determined accurately, and they are more likely to change over the years than lower-order EOFs.

 

The last parameter is the range accuracy of the MBES. Smaller spatial error leads to greater accuracy of sound speed recovery, as expected.

 

 

Figure 3: Left: Histogram of the RMS error for the inverted SSPs for the test dataset. The average SSP RMS error of 0.83 m/s is shown as a vertical dashed line. Right: One example of a true and inverted SSP with an error of 0.83 m/s.

 

5 SUMMARY

 

We showed with an elementary example that the inverse problem of recovering an SSP from travel time measurements is generally non-unique, leading to an ill-posed inverse problem. However, by using EOFs from historical data and a Tikhonov-type regularization we showed that we are able to recover SSPs with an accuracy of 0.83 m/s for an MBES system with a maximum error in its spatial estimates of 1 cm. Larger spatial errors translate into increased errors in the SSP up to 1.36 m/s for spatial errors of 10 cm, which still outperforms the climatological alternative from the World Ocean Atlas¹⁶ by a factor close to 2.

 

ACKNOWLEDGMENTS

 

This research has been supported by the European Commission within the context of the project SMAUG (Smart Maritime and Underwater Guardian), funded under EU Horizon Europe Grant Agreement 101121129.

 

We would like to thank Ailin Brakstad for preparing the sound speed data used in this study.

 

REFERENCES

 

1. B. Ramirez and R. J. Bunker, “Narco-Submarines: Specially fabricated vessels used for drug smuggling purposes,” tech. rep., The Foreign Military Studies Office, Fort Leavenworth, Kansas, USA, Jan. 2015.

2. S. Savitz, I. Blickstein, P. Buryk, R. W. Button, P. DeLuca, J. Dryden, J. Mastbaum, J. Osburg, P. Padilla, A. Potter, C. C. Price, L. Thrall, S. K. Woodward, R. J. Yardley, and J. M. Yurchak, “U.S. Navy employment options for unmanned surface vehicles (USVs),” tech. rep., RAND Corporation, Santa Monica, CA, USA, 2013.

3. M. R. Carnes, “Description and evaluation of GDEM-V 3.0,” Tech. Rep. NRL/MR/7330–09-9165, Naval Research Laboratory, Feb. 2009.

4. K. T. Hjelmervik and K. Hjelmervik, “Estimating temperature and salinity profiles using empirical orthogonal functions and clustering on historical measurements,” Ocean Dynamics, vol. 63, pp. 809–821, July 2013.

5. E. Pauthenet, L. Bachelot, K. Balem, G. Maze, A.-M. Tréguier, F. Roquet, R. Fablet, and P. Tandeo, “Four-dimensional temperature, salinity and mixed-layer depth in the Gulf Stream, reconstructed from remote-sensing and in situ observations with neural networks,” Ocean Science, vol. 18, pp. 1221–1244, Aug. 2022. Publisher: Copernicus GmbH.

6. S. Jin, W. Sun, J. Bao, M. Liu, and Y. Cui, “Sound velocity profile (SVP) inversion through correcting the terrain distortion,” The International Hydrographic Review, 2015.

7. C. Didier, E. Jaouad, G. Gaspard, and L. Michel, “Real-time correction of sound refraction errors in bathymetric measurements using multiswath multibeam echosounder,” in OCEANS 2019 - Marseille, pp. 1–7, June 2019.

8. T. H. Mohammadloo, M. Snellen, W. Renoud, J. Beaudoin, and D. G. Simons, “Correcting multibeam echosounder bathymetric measurements for errors induced by inaccurate water column sound speeds,” IEEE Access, vol. 7, pp. 122052–122068, 2019.

9. L. Keyzer, T. H. Mohammadloo, M. Snellen, J. Pietrzak, C. Katsman, Y. Afrasteh, H. Guarneri, M. Verlaan, R. Klees, and C. Slobbe, “Inversion of sound speed profiles from MBES measurements using Differential Evolution,” Proceedings of Meetings on Acoustics, vol. 44, p. 070035, June 2021. Publisher: Acoustical Society of America.

10. K. T. Hjelmervik, “Inverting the water column sound speed,” in Proceedings of the 10th European Conference on Underwater Acoustics (ECUA 2010), Istanbul, Turkey, July 2010.

11. J. Mueller and S. Siltanen, Linear and nonlinear inverse problems with practical applications, Computational Science and Engineering series, Philadelphia: Society for Industrial and Applied Mathematics, 2012.

12. P. Stefanov, G. Uhlmann, A. Vasy, and H. Zhou, “Travel time tomography,” Acta Mathematica Sinica, English Series, vol. 35, pp. 1085–1114, June 2019.

13. F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt, Computational Ocean Acoustics, Modern Acoustics and Signal Processing, New York, NY: Springer, 2011.

14. R. W. Preisendorfer, Principal component analysis in meteorology and oceanography, no. 17 in Developments in Atmospheric Science, Amsterdam; New York, NY, USA: Elsevier; Elsevier Science Pub. Co., 1988.

15. A. Brakstad, G. Gebbie, K. Våge, E. Jeansson, and S. R. Ólafsdóttir, “Formation and pathways of dense water in the Nordic Seas based on a regional inversion,” Progress in Oceanography, vol. 212, p. 102981, Mar. 2023.

16. M. Locarnini, A. Mishonov, O. Baranova, T. Boyer, M. Zweng, H. Garcia, J. Reagan, D. Seidov, K. Weathers, C. Paver, and I. Smolyar, “World Ocean Atlas 2018, Volume 1: Temperature,” Jan. 2018.