Automatically Differentiable Higher-Order Parabolic Equation for Real-Time Underwater Sound Speed Profile Sensing

<jats:p>This paper is dedicated to the acoustic inversion of the vertical sound speed profiles (SSPs) in the underwater marine environment. The method of automatic differentiation is applied for the first time in this context. Representing the finite-difference Padé approximation of the propagation operator as a computational graph allows for the analytical computation of the gradient with respect to the SSP directly within the numerical scheme. The availability of the gradient, along with the high computational efficiency of the numerical method used, enables rapid inversion of the SSP based on acoustic measurements from a hydrophone array. It is demonstrated that local optimization methods can be effectively used for real-time sound speed inversion. Comparative analysis with existing methods shows the significant superiority of the proposed method in terms of computation speed.</jats:p>

section_content: Spatial variations in sound speed and the refraction they cause have a significant impact on the propagation of acoustic waves in the underwater environment [1, 2]. Sound is the primary mean of transmitting information in the underwater marine environment, so the sound speed profile (SSP) affects underwater communication [3], navigation [4], monitoring [5], and control of underwater objects. Depending on the specific SSP, the signal can either propagate for hundreds or thousands of kilometers or quickly attenuate near the source [6]. 

There are several ways to estimate the SSP [7] ; these can be divided into direct measurement methods and inversion techniques. The most direct method is to measure the travel time of a high-frequency acoustic wave between a source and a receiver placed at a known distance from each other within a specialized sensor. Another common direct method involves measuring temperature, salinity, and depth-the three main parameters that influence the sound speed in water. Although direct measurement methods are the most accurate, the requirement for the device to be located at the measurement point makes real-time monitoring of the SSP quite challenging. 

The method of ocean acoustic tomography involves determining the SSP by measuring the characteristics of the acoustic wave propagation. Sound from a known source is received by an array of hydrophones that measure the signal's amplitude and phase. By analyzing the deviations that the signal undergoes on its path from the source to the receivers, it is theoretically possible to assess the factors that caused these deviations. 

From a mathematical perspective, the complexity of the tomography problem lies in its ill-posedness according to Hadamard [8, 9]. Typically, the problem has multiple solutions, and additional constraints must be sought and applied to reduce the problem to a single solution. Finding the inverse nonlinear operator is one of the most challenging tasks in modern functional analysis and the theory of inverse problems. 

The most well-studied inversion method is the matched field processing (MFP) technique [10]. MFP involves using a forward model that associates a given SSP and other input data with the acoustic pressure at specific points in space. Based on such a model, it is straightforward to formulate the problem of finding an SSP for which the measured and model-predicted values match. It is important to note that the relationship between SSP and the spatial distribution of acoustic pressure is highly nonlinear. Minimizing this nonlinear functional is a challenging task from both theoretical and computational perspectives. Stochastic optimization methods [11], such as genetic algorithms, simulated annealing, or particle swarm optimization, are commonly applied. To regularize the problem and narrow the search space, the empirical orthogonal function (EOF) method is used, which analyzes historical data on sound speed distribution in a specific region to find an orthogonal basis [12, 13]. However, solving this problem using stochastic optimization methods can take hours. 

MFP is widely used not only for determining SSPs but also for geoacoustic inversion [14, 15], such as estimating seabed parameters. Additionally, MFP is employed to locate the source of acoustic signals [1, 16]. The general principle remains the same across applications: the goal is to adjust the model parameters so that the predicted and measured values align. 

The application of the compressive sensing method [17] [18] [19] reduces the problem to an underdetermined system of linear algebraic equations. This is achieved by utilizing a basis in the space of SSP functions and linearizing the nonlinear operator using the first term of the Taylor series. Given the nonlinear nature of the problem, this approach is only applicable for small variations in the sound speed. 

Inversion methods based on machine learning and deep neural networks are gaining popularity [20] [21] [22] [23] [24]. Like most machine learning models, these inversion methods are difficult to interpret and analyze, and it is difficult to determine the limits of their applicability. Additionally, it is not always feasible to collect a large enough dataset for training. Despite the challenges, the results presented in these works are quite promising and, in many cases, surpass traditional methods. 

It is important to note that acoustic measurements are not the only way to invert the SSP. For instance, in the works [25] [26] [27], inversion models based on satellite observation data have been developed. The accumulated satellite data on the ocean surface and underwater environment [28] enables the construction of complex regression models. 

Given the advancement of machine learning tools, such as automatic differentiation [29] and gradient-based optimization, there arises the idea of applying these tools not to poorly interpretable deep neural networks, but to the existing physical models of wave propagation and numerical methods. In this work, a fully differentiable higher-order numerical scheme based on the parabolic equation (PE) method [30] has been developed. The ability to perform automatic differentiation with respect to SSP parameters enables the use of gradient-based methods for rapid optimization problem solving. 

A somewhat related approach is proposed in [31], where a convolutional generative network is trained on historical data and used to subsequently predict the SSP in real time. 

The paper is organized as follows. Section 2 presents the formal problem statement. Section 3 is dedicated to solving the corresponding forward problem using the PE method. In Section 4, the original problem is reduced to a functional minimization problem, and various solution approaches are discussed. The analysis of numerical modeling results is provided in Section 5. The conclusion outlines several questions for the further research.
section_name: Introduction

section_content: There is an acoustic wave source with known parameters, coordinates, and depth. At a known distance from the source, an array of hydrophones is located, measuring the acoustic signals from the source at different time points. The seabed topography and parameters are also considered known. The amplitude and phase of the signals at the hydrophones vary over time, mainly due to the SSP variations. We assume that the SSP changes smoothly over time. The complex acoustic pressure v i j is measured at points x j, z j, j ∈ [1, N]. Measurements are taken at known discrete time intervals t i. 

The task is to determine the current SSP at any given moment in time t i based on the hydrophone data v i j. This means constructing the time-varying SSP ci (z) as a function of both depth and time. The method should account for smooth temporal variations in the SSP and provide an updated estimate of the sound speed at each time step based on the observed changes in the acoustic pressure v i j. A schematic description of the problem setup is shown in Figure 1.
section_name: Problem Statement

section_content: At the core of almost any method for solving an inverse problem lies the repeated solution of the corresponding forward problem and the search for an optimal solution. In our case, the forward problem involves computing the spatial distribution of acoustic pressure in the underwater environment with known environmental parameters, including the SSP and seabed characteristics. Additionally, the parameters of the acoustic wave source are known. The problem, as formulated, has a unique solution and can be efficiently solved using various numerical methods. The most effective methods include the PE method [30] and the normal modes method [1]. In this work, we use the one-way Helmholtz equation method [32, 33], which is a generalization of the PE.
section_name: Forward Problem and Its Solution

section_content: The complex-valued acoustic field ψ(x, z) at frequency f satisfies the two-dimensional Helmholtz equation

is the wave number dependent on depth z. • c(z) is the SSP, which is the primary variable of interest in this research. While the sound speed may also vary in the horizontal direction x, these variations are generally weak and do not have as significant impact as depth-dependent variations. In the forward problem, c(z) is assumed to be known. • ρ(x, z) is density, a parameter of the seabed that may vary with both depth and range. 

In water, it is typically assumed to be 1.
section_name: Mathematical Formulation of the Forward Problem

section_content: Q is a complex constant responsible for the amplitude and phase of the emitting source. As it will be clear later, solution ψ is linear with respect to coefficient Q, so for simplicity, we will seek the solution of Equation (1) in the case where Q = 1. 

• ψ 0 (z) is a known function that defines the directivity and depth of the emitting source. 

In this work we use a Gauss source, defined as follows:

where z 0 is the antenna depth, θ bw is the beam width, θ e is the elevation angle, and

The problem of determining the acoustic field ψ(x, z) at any point in space is wellposed, meaning that it has a unique solution given the specified environmental parameters.
section_name: •

section_content: Let us rewrite the Helmholtz Equation ( 1 ) in the form of the sequential application of two operators, responsible for wave propagation in the left and right directions, respectively:

By discarding the second term, and thus, ignoring backscattering, we obtain a one-way approximation of the Helmholtz equation:

The method used in this work is a step-by-step evolutionary solution of the one-way Helmholtz Equation (3) along the positive direction of the x-axis. The step-by-step solution is computed simultaneously in both the water column and the seabed. To maintain a finite computational domain, the transparent boundary condition method is applied [34] [35] [36] [37]. 

The stepwise solution of the one-way Helmholtz Equation (3) in operator form is written as follows:

where

The finite-difference numerical solution will be sought on a uniform grid in both the x and z dimensions, with steps ∆x and ∆z, respectively. We use the notation u n j (z) = u(n∆x, j∆z), where u n j (z) represents the solution at the discrete points of the grid. For the numerical solution, we use a fourth-order finite-difference approximation of the transversal operator L [38],

where the modified second difference operator is defined as follows:

The next step in constructing the numerical scheme is to approximate propagation operator (5). For this, we treat operator P as a function of operator L and expand it in the vicinity of L ≈ 0 using a rational Padé approximation of order [p/q] [33,39]:

Then, the action of the propagation operator (5) at each step is written as the following system of equations, which is solved sequentially from top to bottom:

Each row of system (7) represents a tridiagonal system of linear algebraic equations, which can be solved in linear time using the Thomas algorithm (tridiagonal matrix algorithm). It is important to note that the entire numerical scheme consists solely of solving tridiagonal systems of equations, which are efficiently solved by sequential application of elementary operations. 

The optimal choice of the computational grid parameters ∆x, ∆z, β, and approximation order [p/q] is discussed in Refs. [40] [41] [42]. The optimal values of the computational grid parameters are highly dependent on the input data of the problem; specifically, the propagation environment and the characteristics of the acoustic source. A typical choice for the Padé approximation order is [7/8]. 

It is worth noting that the frequently used Claerbout's wide-angle (WA) PE, solved using the Crank-Nicolson method [38], is a special case of the above-described scheme when using a [1/1] Padé approximation. Thus, the scheme considered here provides a significantly higher order of accuracy.
section_name: One-Way Helmholtz Equation and Its Numerical Solution

section_content: Let the acoustic pressure values measured at the hydrophones be denoted as the vector

Let us define the operator G : L 2 [0, H] → C N, which maps SSP function c(z) to the vector of acoustic pressure values at the measurement points

The inverse problem consists of finding the inverse operator G -1, which, based on the measurement data, would return the SSP c(z). Operator G is nonlinear, which significantly complicates the task of finding its inverse [8]. We define the inverse operator G -1 through the minimization of the residual (discrepancy) between the measured and predicted acoustic pressure values:

In practice, the amplitude and phase of the source Q are often unknown. Therefore, the minimization problem (8) should be reformulated to account for this uncertainty:

which is equivalent [43] to

Rewriting this as a minimization problem gives

where γ is a regularization parameter [44]. In all subsequent experiments, we set γ = 1. It is important to note that since the primary influence on propagation characteristics comes from the spatial variations in the sound speed rather than its absolute value, regularization is applied with respect to the norm of the transversal derivative of the SSP. 

In this work, we assume that both the amplitude and phase of the signal at the receivers are known. The modifications to the loss function for the case where only the amplitude of the received signal is known are provided in [45]. 

We will seek function c(z) in the form of a piecewise linear function, defined on a predefined uniform grid 0

Note that the discretization grid for the SSP does not necessarily have to coincide with the computational grid used for solving the forward problem. In this case, the minimization problem (8) simplifies somewhat, as we now seek the minimum in an M-dimensional space:

In addition, physical constraints should be imposed on the sound speed values, such as

Nevertheless, the minimization problem remains extremely challenging due to the high dimensionality and the complex landscape of the objective function [46]. Additionally, the value of objective function J(c, v) for each c depends on the numerical solution of the Helmholtz Equation (1), making the evaluation of the objective function a computationally expensive operation.
section_name: Methods for Solving the Inverse Problem

section_content: Thus, to solve the given problem, we need to minimize the objective functional J(c, v). Optimization methods (for finding minimum of an objective function) are generally divided into two categories: local and global [11, 47]. Global methods include genetic algorithms [48], the simulated annealing method [49], and particle swarm optimization method. These methods typically treat the loss function as a black box, requiring only the ability to compute its value at any point within the domain. For example, SAGA 5. 4 software package [50] implements environmental parameter inversion, including SSP inversion, based on genetic algorithms. 

However, global optimization methods are not feasible for real-time applications due to their slow convergence. Global methods are typically used for acoustic inversion of seabed parameters [14], as these parameters do not change over time and the inversion speed is not as critical in such cases. 

Local optimization methods typically require knowledge of the gradient of the loss function

and an initial guess c0 (z). The gradient helps them find the minimum more quickly and accurately, but only if the minimum is located near the initial guess. 

There are also hybrid algorithms that combine the strengths of both global and local optimization methods. In [51], a combination of a genetic algorithm and the Gauss-Newton method was used for acoustic inversion of layered medium parameters. The synthesis of the simulated annealing and simplex methods for seabed parameter inversion was investigated in [52]. 

In this work, it is assumed that changes in the SSP occur smoothly over time, meaning that each subsequent change is small compared to the current state. As a result, it is relevant to use local optimization methods based on gradient descent. Since the profile evolves gradually, these methods can effectively track the small variations and provide accurate updates with each time step.
section_name: Local and Global Optimization

section_content: Local optimization methods utilize not only the values of the objective function but also its gradient. This imposes certain constraints on the objective function-it must be differentiable, and it can no longer be treated as a "black box". Although the gradient can be computed using finite-difference methods, this approach is typically inefficient in practice. The advantage of using the gradient can be fully exploited if it is computed analytically. However, this is often a non-trivial task. For example, in our case, the optimization takes place over an infinite-dimensional function space, and the loss function depends on the solution of the Helmholtz equation. 

Nevertheless, several successful attempts have been made to compute the gradient of the objective function using the adjoint equation method, followed by its application to invert the SSP and seabed parameters within the PE method [44, [53] [54] [55]. One drawback of these methods is that the analytical expressions for the gradient are often complex and computationally intensive. As a result, in all these works, a simple finite-difference approximation of the PE was used, based on the Crank-Nicolson method and corresponding to a Padé approximation of order [1/1]. 

As we will see in further computational experiments, using a higher-order approximation can significantly increase the computational efficiency. Moreover, as expected, the convergence of this local method is highly dependent on the initial guess. Specifically, it converges to the correct solution only when the initial approximation is sufficiently close to the true profile. However, practical applications of such local optimizers have not been widely explored until now.
section_name: Gradient Using the Adjoint Method

section_content: In this work, an alternative to manually computing gradient expressions of the PE is proposed. The numerical scheme is represented as a computation graph, over which automatic differentiation is performed [29]. Indeed, the considered numerical scheme (7) consists of a large, but fixed and predefined, number of elementary operations. The gradient does not need to be derived or programmed manually; it is computed analytically using the computation graph of the numerical scheme. It should be noted that the construction of the computation graph requires a predetermined finite number of iterations, which means that various adaptive algorithms cannot be differentiated. 

The development of deep learning methods has introduced several highly efficient and accessible software frameworks for automatic differentiation, such as TensorFlow 2. 18, PyTorch 2. 5, and JAX 0. 4. While these frameworks are primarily designed for models in the form of deep neural networks, they can also be used for other finite-dimensional models. From a computational perspective, determining the coefficients of a neural network is not significantly different from determining the coefficients c 0,. . ., c M of the SSP. 

Programming models in frameworks like TensorFlow and PyTorch requires a very specific representation as a computation graph, typically designed for deep neural networks. An alternative approach is offered by the JAX framework [56]. JAX allows users to program models in a more familiar format, using an interface that closely mirrors standard numerical modeling frameworks like NumPy 2. 1 and SciPy 1. 14. The model is automatically translated into a computation graph suitable for automatic differentiation. 

The JAX framework is already being actively used to build differentiable numerical schemes and solve similar problems across various fields. The JAX-FEM library [57] is an implementation of the finite element method (FEM) in pure Python 3. By using JAX, the library enables automatic differentiation of FEM-based numerical schemes with respect to any parameter, simplifying the solution of inverse problems, particularly in computational mechanics. Similarly, the XLB library [58] provides differentiable algorithms for solving fluid dynamics problems. The Diffrax library [59] offers a JAX-based implementation of well-known methods for solving first-order differential equations, such as the Runge-Kutta method. Numerical solutions are differentiable with respect to any predefined parameters, making them suitable for use in gradient-based optimization methods. This library also supports stochastic differential equations. 

It is also important to mention the pioneering work in this field [60], which demonstrated the connection between first-order differential equations and ResNet architectures. In this work, the authors were the first to propose an algorithm for training a differential equation by differentiating its solution with respect to certain parameters.
section_name: Automatic Differentiation of the Numerical Scheme

section_content: Although having access to the gradient of the numerical scheme simplifies the optimization problem, finding a solution within an acceptable time frame is still not feasible without prior information. In this work, we assume that the SSP at time t 0 is known. 

We introduce an additional dependence of the inverse operator on the SSP from the previous step ci-1 (z) (denoted as G -1 ci-1 v), which helps the algorithm to converge faster. Given that the problem is solved in real time and has continuous access to the acoustic pressure measurements, this requirement can be easily satisfied. The only challenge is the "cold start" problem-how to initialize c0 (z) on the very first iteration. In this work, we propose initializing c0 (z) using data from direct measurements. If direct measurements are not available, the problem can be approached with global optimization methods. The initial solution using global optimization methods can take a long time to find (maybe several hours), but then, by simultaneously saving measurements over time, it is possible to quickly catch up to real time using the proposed method, since it searches for the optimum much faster than the SSP actually changes. 

By using the mesh optimizer proposed in [40], the grid shapes can be precisely computed in advance, eliminating the need for adaptive algorithms, which would be impossible to differentiate. 

In this work, the operator G -1 ci-1 is implemented using the L-BFGS-B algorithm [61], with the initial guess set to ci-1. The implementation from the SciPy library [62] is used. The gradients and function values computed using JAX can be easily utilized in external libraries. Specifically, the gradient of the numerical scheme in our implementation is seamlessly passed to the SciPy optimizer. The pseudocode for the proposed method is shown in Listing 1. 

Listing 1. Proposed real-time method for the SSP inversion based on the acoustic pressure measurements. c0 (z) ←ssp_direct_measure(t 0 ) // direct measurement of the SSP at time t 0 for i = 1. . . // measurement of acoustic pressure at time

, near ci-1 // return the current SSP and continue the inversion process yield ci (z)
section_name: The Proposed Algorithm

section_content: The proposed method has been implemented as part of the open-source library Py-WaveProp [63], developed by the author of the present research. The examples presented in this section were implemented using the mentioned library. 

To demonstrate the functionality of the proposed algorithm, we replace the actual pressure measurements on the hydrophones with data generated from forward modeling. Using the pre-known dynamics of the sound speed changes c i (z), i = 1,. . ., we compute the acoustic pressure values at the hydrophones at some time points t i. Gaussian noise is added to the computed values at the hydrophones, and these noisy values are used to solve the inverse problem. We then compare the true SSP with the inverted one. The pseudocode of the simulation process is shown in Listing 2. 

Listing 2. Pseudocode of the simulation process. 

// simulate data on the hydrophones using forward modeling

This simulation approach allows us to validate the performance of the inversion algorithm by generating synthetic data that mimics real-world conditions, including noise, and observing how well the method can recover the true SSP. 

As an example, we model the sensing of a bilinear SSP in shallow water. The depth is constant and equal to 200 m. The acoustic wave source is located at a depth of 50 m and emits a monochromatic signal. Frequencies of 200 Hz and 500 Hz are considered. A vertical array of eight uniformly distributed hydrophones (5-75 m) measures the signal from the source. The array is placed at distances of 5 and 10 km from the source. The signal-to-noise (SNR) ratio on the receivers is simulated equal to 30 dB. 

The initial profile is a constant SSP of 1500 m/s. First, we increase the positive sound speed gradient, which encourages multiple reflections of waves from the sea surface. Then, we gradually add a negative gradient near the surface, which causes the acoustic waveguide to sink to a certain depth. A total of 50 sequential measurements and SSP inversions are performed. The dynamics of the SSP changes are shown in Figure 2. This figure also shows the inversion results for various frequencies and hydrophone array locations. For the readers convenience, enlarged images of SSPs for several values of t are also shown in Figure 3. The pointwise difference between the original and inverted SSPs c(z)c(z) is shown in Figure 4. 

Given that propagation of acoustic waves is more influenced by the spatial variation of sound speed rather than its absolute value, it makes sense to evaluate the error between the true and inverted profiles in terms of the norm of the difference of the gradients (spatial changes) rather than the absolute sound speed values:

The error dynamics over time are shown in Figure 5. It can be observed from the inversion results that at a frequency of 200 Hz, the method performed excellently at distances of 5 and 10 km from the source. Increasing the frequency to 500 Hz, which decreases the wavelength and increases the distance from the source to the receiver in terms of wavelengths, already introduces difficulties for the inversion. Specifically, at a distance of 5 km, the algorithm still managed to perform the inversion, but at 10 km, it diverged. It should also be noted that increasing the number of measurements did not help to achieve convergence. Evidently, in this case, the data from hydrophones at a single distance are no longer sufficient to uniquely determine the SSP. In situations where the algorithm converges, the error remains bounded over time, indicating the stability of the proposed algorithm. 

Let us now turn to the question of the convergence speed of the proposed algorithm. In Figure 6, the number of function evaluations and gradient computations at each time step is shown. It can be seen that this number generally does not exceed 500. However, in the most challenging case, the number of calls increased sharply, yet this still did not help the algorithm to converge. As previously noted, earlier gradient-based inversion methods were based on the Crank-Nicolson method for the WA PE models. Since the Crank-Nicolson method is a special case of the rational Padé approximation used in this work, we can easily compare the methods even without having their specific software implementations. In Table 1, the grid cell sizes used and the average inversion time are shown. The optimal grid parameters were selected using the algorithm developed in [40]. It is evident that the proposed method is approximately 20 times faster than the previously used methods. Using a higher-order scheme makes it possible to perform computations much faster on a more sparse grid. Moreover, it is clear that the inversion time is directly dependent on both the frequency and the distance from the source to the receivers. The greater the distance from the source to the receivers in wavelengths, the longer it takes for the inversion algorithm to converge. 

The effect of divergence as the receivers move away from the source has been noted for local methods before [44, 64]. The complexity of the problem at large distances lies in the complication of the objective function in the vicinity of the initial approximation. This is confirmed by the increase in convergence time with increasing distance. Thus, the local method hits a local optimum. Ref. [44] notes that horizontally arranged hydrophone arrays are less susceptible to this effect, but they are significantly more difficult to deploy and operate. Two-dimensional distributions of the acoustic pressure at different time points are shown in Figure 7. These distributions were computed using forward modeling for both the true and inverted SSPs. It can be observed that, at the same time points, the distributions are nearly indistinguishable, which attests to the accuracy of the obtained solution. Furthermore, the waveguide effects, which emerge at various moments, are clearly visible. These waveguide effects are what the hydrophones detect, and they form the basis for the inversion process. Finally, we compare the accuracy of the proposed algorithm depending on the SNR value at the receivers. Figure 8 depicts absolute error of the inversion for four values of SNR. Obviously, decreasing the SNR leads to increasing the inversion error. However, even for an SNR value of 10 dB, the error remains bounded over time, indicating that the proposed algorithm is robust to receiver noise.
section_name: Numerical Results and Discussion

section_content: The proposed method significantly outperforms, on average, the existing approaches in terms of inversion speed, making it suitable for real-time applications on low-powered computational equipment. 

The method leverages modern machine learning tools while remaining physically and computationally well-founded and interpretable. Unlike many machine learning approaches, the model does not require the prior collection of large amounts of statistical data about the environment. It only requires an initial estimate of the SSP. 

A key area for future research is addressing the "cold start" problem. This will involve developing modifications that incorporate elements of global optimization to determine the initial profile without the need for direct measurements. Combining global and local optimization could also improve the method's convergence, particularly when the source and receivers are located far apart. 

Gradient-based local search methods have already been tested on real data and showed good results there [53]. The method proposed in the present research performs the same process, but much faster. Nevertheless, a comprehensive analysis of the method on experimental data is planned for future work. Another important direction is the uncertainty quantification of the solution obtained using the proposed method.
section_name: Conclusions

section_content: Funding: This study was supported by the Russian Science Foundation grant No. 23-71-01069. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable.

section_name: Conflicts of Interest:

section_content: Data Availability Statement: Data is contained within the article.

section_content: The author declares no conflicts of interest.
section_name: Conflicts of Interest:

section_name: Abbreviations

section_content: The following abbreviations are used in this manuscript:
section_name: Abbreviations