Compressive Sound Speed Profile Inversion Using Beamforming Results

Abstract

Sound speed profile (SSP) significantly affects acoustic propagation in the ocean. In this work, the SSP is inverted using compressive sensing (CS) combined with beamforming to indicate the direction of arrivals (DOAs). The travel times and the positions of the arrivals can be approximately linearized using their Taylor expansion with the shape function coefficients that parameterize the SSP. The linear relation between the travel times/positions and the shape function coefficients enables CS to reconstruct the SSP. The conventional objective function in CS is modified to simultaneously exploit the information from the travel times and positions. The SSP is estimated using CS with beamforming of ray arrivals in the SWellEx-96 experimental environment, and the performance is evaluated using the correlation coefficient and mean squared error (MSE) between the true and recovered SSPs, respectively. Five hundred synthetic SSPs were generated by randomly choosing the SSP dictionary components, and more than 80 percent of all the cases have correlation coefficients over 0.7 and MSE along depth is insignificant except near the sea surface, which shows the validity of the proposed method.

1. Introduction

It is crucial to estimate the sound speed profile (SSP) of a water column to reliably locate a source or to conduct geoacoustic inversion [1,2,3,4,5]. Many studies have used acoustic measurement data to infer the SSP (i.e., acoustic tomography) [6,7,8,9,10]. Originally, ocean acoustic tomography is conducted using travel times of arrivals measured at a single hydrophone [6]. Its application is limited by ocean noises that mask the arrivals. To mitigate the problem, multiple measurement at sensors in array is used with matched field inversion (MFI) [8,9]. In the SSP inversion process, replicas (synthetic acoustic data at the sensors of array) are generated according to the SSPs using underwater acoustic propagation model [11], and an objective function is defined to evaluate the similarity between the measured data and the replicas. The replicas corresponding to the maximum (or minimum) value from the objective function should then resemble the SSP from the ocean. In general, a genetic algorithm or simulated annealing is applied to solve the optimization problem [2,3]. It estimates the SSP more robustly at the cost of computational burden resulting from iterative use of propagation model. However, its performance relies on the propagation model and it cannot be applied to environments, where the propagation model is restricted owing to complexities of the ocean including geoacoustic parameters (sound speed, density, and attenuation) of sediment. This work suggests a novel inversion scheme for SSP, which employs compressive sensing (CS) with beamforming results (i.e., direction of arrivals (DOAs)). It inverts the SSP with enhanced signal-to-noise ratio (SNR) due to the beamforming and without intensive use of propagation model as well as geoacoustic parameters.

CS is a method that can solve an underdetermined linear system ($\mathbf{y}=\mathbf{A}\mathbf{x}$) by using a sparsity promoting objective function [12,13]. $\mathbf{y}$ is a measurement vector, and $\mathbf{x}$ is an unknown vector that is to be recovered. $\mathbf{A}$ is a sensing matrix that shows the linear relation between $\mathbf{x}$ and $\mathbf{y}$, and these are defined differently according to the given problem. In underwater acoustics, CS has been applied to obtain a high-resolution DOA with a limited number of sensors in an array [14,15,16,17,18], and the applicability of CS has been expanded to source localization [19,20] and ocean parameter inversion [10,21,22].

Recently, CS has been applied to obtain the SSP represented as a linear combination of shape functions in the SSP dictionary, with a sparsity-promoting objective function used to exploit the sparsity of the solution after constructing a linear relation between the (given) measurement and (unknown) parameter [10]. The measurement and the parameter to be recovered in the compressive SSP estimation are the acoustic fields at the sensors of a vertical line array (VLA) and the shape function coefficients, respectively. The acoustic fields are approximately linearized in terms of the shape function coefficients parameterizing the SSP. The SSP is reconstructed using the compressive SSP estimation that replicates the acoustic pressures at the sensors. During the process of replicating the pressures, information of geoacoustic parameters is needed.

In this work, beamforming results from the array are used with CS to infer the SSP. Practically, the DOAs are estimated with various beamforming schemes [11,14,15,16,17,18]. Here, the DOAs are calculated with eigenrays from a ray model. For convenience, a numerical experiment showing the DOAs matching the arrival angles of the eigenrays is illustrated. While the SSP is represented as a linear combination of the shape functions (e.g., empirical orthogonal functions (EOFs)) in the SSP dictionary [1,2,4,8,10,23], the travel time and position are linearized at a specific DOA using their Taylor expansion with the shape function coefficients, respectively. Thus, the proposed method eliminates the need for geoacoustic parameter information. The scheme on the inversion of SSP with the given beamforming results is delineated as follows. Section 2 introduces the detailed process for compressive SSP inversion using the beamforming results. A conventional objective function in CS is modified to simultaneously exploit the information from the position and travel time. Numerical experiments are conducted in the SWellEx-96 environment to examine the performance of the compressive SSP inversion. Finally, Section 4 summarizes this work.

2. Compressive SSP Inversion

This work adopts CS to recover the fine-scale structure of the SSP since the CS demonstrates a good performance for the SSP inversion with a limited number of observations [10]. However, instead of using acoustic fields that require geoacoustic parameters, the beamforming results from the VLA are used. The DOAs can be obtained using frequency-domain beamforming with the VLA, and the travel times (or relative delay times) of the arrivals can be acquired according to the DOAs using time-domain beamforming. The information for the arrivals regarding the travel times and positions is highly dependent on the SSP, and it is used in the current compressive SSP estimation.

2.1. Formulation of the Compressive SSP Inversion

SSP is represented using a linear combination of the shape functions in the SSP dictionary as below [10,23]:

$$\mathbf{c}\left(\mathbf{x}\right)={\mathbf{c}}_{\mathbf{0}}+\mathbf{Q}\mathbf{x},$$

(1)

where $\mathbf{c}\left(\mathbf{x}\right)\in {\mathbb{R}}^K$ and ${\mathbf{c}}_{\mathbf{0}}\in {\mathbb{R}}^K$ are ocean and reference SSPs, respectively, that are evenly captured at K points along the depth. In this work, the source is located near the array (in the numerical experiment of Section 3, and they are apart by a horizontal range of 2 km), and variations of SSPs along the range are neglected. $\mathbf{Q}\in {\mathbb{R}}^{K\times N}$ is an SSP dictionary composed of N unique shape functions ($\mathbf{q}$${}_{\mathbf{1},\dots ,\mathbf{N}}$$\in {\mathbb{R}}^K$) and $\mathbf{x}\in {\mathbb{R}}^N$ is a set of coefficients that correspond to the N shape functions. The shape functions are defined to express SSP with their linear combination, which enable SSP to be parameterized with the coefficients corresponding to the shape functions. The SSP dictionary is a matrix holding the shape functions as its columns. Recently, a dictionary learning in the form of unsupervised machine learning is applied to obtain the shape functions with SSP data, which are beneficial for CS by representing the SSP sparsely [23]. However, in this work, synthetic EOFs as Gaussian functions and conventional EOFs from principle component analysis are separately used as the N shape functions as those in Bianco and Gerstoft [10].

The travel time and arrival position of a ray are highly dependent on $\mathbf{x}$ (or SSP), and rays with the same arrival angles show distinct trajectories for different $\mathbf{x}$. Figure 1 shows ray trajectories for the ocean SSP (dotted line) and reference SSP (solid line). Consider a ray that emanates from the source to the receiver for the ocean SSP. However, the backpropagated ray from the receiver using the reference SSP deviates from the source position in both the vertical and horizontal directions (Figure 1b,c). When overall SSP differences are small, the ocean SSP can be approximated by the reference SSP using Taylor expansion based on shape function coefficients. The depth and travel time for the ocean SSP at the source range can be denoted as follows:

$$D(\mathbf{x};{\theta }_j)=D(\mathbf{0};{\theta }_j)+\sum _{i=1}^N\frac{\partial D(\mathbf{x};{\theta }_j)}{\partial x_i}{|}_{\mathbf{x}=\mathbf{0}}{x}_i+e_D(\mathbf{x};{\theta }_j),$$

(2)

$$T_D(\mathbf{x};{\theta }_j)=T_D(\mathbf{0};{\theta }_j)+\sum _{i=1}^N\frac{\partial T_D(\mathbf{x};{\theta }_j)}{\partial x_i}{|}_{\mathbf{x}=\mathbf{0}}{x}_i+e_{T_D}(\mathbf{x};{\theta }_j),$$

(3)

where $D(\mathbf{x};{\theta }_j)$ and $T_D(\mathbf{x};{\theta }_j)$ are the depth and travel time corresponding to the arrival for which the DOA is ${\theta }_j$. $x_i$ is an element of $\mathbf{x}$. $e_D$ and $e_{T_D}$ are higher order terms in the Taylor expansion for the depth and travel time, respectively, and are deviations of the linearization from the exact expressions.

The depths and travel times for all arrivals after backpropagation from the center of the VLA (representative position of the array) can be expressed using Equations (2) and (3) with vector and matrix notation as follows:

$$\Delta \mathbf{D}={\mathbf{A}}_{\mathbf{D}}\mathbf{x}+{\mathbf{e}}_{\mathbf{D}},$$

(4)

$$\Delta {\mathbf{T}}_{\mathbf{D}}={\mathbf{A}}_{{\mathbf{T}}_{\mathbf{D}}}\mathbf{x}+{\mathbf{e}}_{{\mathbf{T}}_{\mathbf{D}}}.$$

(5)

The components of $\Delta \mathbf{D}\in {\mathbb{R}}^M$ and $\Delta {\mathbf{T}}_{\mathbf{D}}\in {\mathbb{R}}^M$ are the differences between the depths ($D(\mathbf{x};{\theta }_j)-D(\mathbf{0};{\theta }_j)$) and travel times ($T_D(\mathbf{x};{\theta }_j)-T_D(\mathbf{0};{\theta }_j)$) for the ocean and reference SSPs at a specific DOA, respectively, where M is the number of arrivals. Note that $D(\mathbf{x};{\theta }_j)$ is the same as the source depth regardless of ${\theta }_j$ since the backpropagated arrivals for the ocean SSP converge at the source position. $\Delta \mathbf{D}$ is referred to as vertical position difference and is related to ray information on position. The elements of ${\mathbf{A}}_{\mathbf{D},{\mathbf{T}}_{\mathbf{D}}}\in {\mathbb{R}}^{\mathbf{M}\times \mathbf{N}}$ are first partial derivative terms with respect to $x_i$ in Equations (2) and (3), respectively. ${\mathbf{e}}_{\mathbf{D},{\mathbf{T}}_{\mathbf{D}}}$ consist of $e_{D,T_D}$.

Similarly, the range ($R(\mathbf{x};{\theta }_j)$) and travel time ($T_R(\mathbf{x};{\theta }_j)$) for a specific DOA (${\theta }_j$) at the source depth can be linearized using the Taylor expansion and can be denoted in vector and matrix form as follows:

$$\Delta \mathbf{R}={\mathbf{A}}_{\mathbf{R}}\mathbf{x}+{\mathbf{e}}_{\mathbf{R}},$$

(6)

$$\Delta {\mathbf{T}}_{\mathbf{R}}={\mathbf{A}}_{{\mathbf{T}}_{\mathbf{R}}}\mathbf{x}+{\mathbf{e}}_{{\mathbf{T}}_{\mathbf{R}}}.$$

(7)

The range difference ($R(\mathbf{x};{\theta }_j)-R(\mathbf{0};{\theta }_j)$) and horizontal travel time difference ($T_R(\mathbf{x};{\theta }_j)-T_R(\mathbf{0};{\theta }_j)$) are components of $\Delta \mathbf{R}\in {\mathbb{R}}^M$ and $\Delta {\mathbf{T}}_{\mathbf{R}}\in {\mathbb{R}}^M$, respectively. $R(\mathbf{x};{\theta }_j)$ is equivalent to the source range. $T_R(\mathbf{x};{\theta }_j)$ is the same as $T_D(\mathbf{x};{\theta }_j)$ due to their identical travel paths for the ocean SSP. $\Delta \mathbf{R}$ is referred to as horizontal position difference and is related to ray information on position. Elements of ${\mathbf{A}}_{\mathbf{R}}\in {\mathbb{R}}^{\mathbf{M}\times \mathbf{N}}$ and ${\mathbf{A}}_{{\mathbf{T}}_{\mathbf{R}}}\in {\mathbb{R}}^{\mathbf{M}\times \mathbf{N}}$ are $\frac{\partial R(\mathbf{x};{\theta }_j)}{\partial x_i}{|}_{\mathbf{x}=\mathbf{0}}$ and $\frac{\partial T_R(\mathbf{x};{\theta }_j)}{\partial x_i}{|}_{\mathbf{x}=\mathbf{0}}$, respectively. Higher order terms $e_{R,T_R}(\mathbf{x};{\theta }_j)$ in the Taylor expansion constitute ${\mathbf{e}}_{\mathbf{R},{\mathbf{T}}_{\mathbf{R}}}$.

The linear relationship between the ray information and shape function coefficients are constructed as Equations (4)–(7), which enable the CS to be applied to estimate the ocean SSP. However, the linear relationship has the inherent error owing to higher order terms in the Taylor series. The validity of the approximated expression is examined numerically according to the DOAs before reconstruction of the SSP using CS.

The exact difference in the vertical travel time is calculated using a conventional ray model [11] according to the coefficients of the first EOF in SWellEx-96; the SSP is represented by $\mathbf{c}={\mathbf{c}}_{\mathbf{0}}+\mathbf{q}{}_1{x}_1$, where $\mathbf{q}{}_1$ and $x_1$ are respectively the first EOF among 26 EOFs formed by CTD (Conductivity, Temperature, Depth) casts from SWellEx-96 and the corresponding coefficient. It is compared to the approximated value from the linear expression, as in Equation (3). The numerical environment is the same as that in the paper [10] except for the number and position of the receiver, as shown in Figure 2. The center of the VLA is used as a representative position of the receiver. Figure 3 shows comparisons of results for arrivals with low (DOA of $-7.5^{\circ }$, Figure 3a), moderate (DOA of $13.7^{\circ }$, Figure 3c), and high arrival angles (DOA of ${56.3}^{\circ }$, Figure 3e) and their corresponding ray trajectories for the reference SSP (i.e., $x_1=0$). Here, instead of using beamforming, DOAs are calculated with eigenrays from the ray model, which connect the source to the receiver. The rays are backpropagated from the receiver with their arrival angles. The first order derivatives in the linear approximation are calculated numerically according to the DOAs, which correspond to gradients at 0 in Figure 3a,c,e, respectively. A ray trajectory for a low grazing angle tends to have a turning point in the water column as shown in Figure 3b, which is sensitive to the variation in the SSP. Since the ray information depends nonlinearly on the shape function coefficient, the linear approximation is inappropriate at a low grazing angle (Figure 3a). The ray for the moderate angle spans over full water column, and the turning is a result of the reflection by the boundaries of the waveguide (Figure 3d). Its trajectory is less sensitive to the SSP, and its information can be linearized for a small fluctuation in the SSP as seen in Figure 3c. Nonlinearity is shown with the increase in absolute value of the shape function coefficient (i.e., increment of $\left|x_1\right|$), as expected. The ray trajectory at a high grazing angle is insensitive to the variation in the SSP, and no notable difference is observed in the results of the comparison (Figure 3e).

The numerical experiment demonstrates that the valid region for the approximation is dependent on the DOA, whose accuracy increases with an increase in the DOA. While the ray information at a high grazing angle can be linearized as in Equations (4)–(7), it is unavailable due to the large reflection loss by multiple reflections (Figure 3f). Thus, in this work, DOAs of less than ${30}^{\circ }$ are used for the CS SSP inversion. Rays having turning points within the water column are excluded since the approximated expressions are improper for these rays. The results of the inversion with moderate DOAs will be shown in Section 3.

2.2. Objective Function in the Compressive SSP Inversion

In a standard approach, the CS is applied individually to each ray information to recover the ocean SSP as follows:

$$\widehat{\mathbf{x}}\left(\mu \right)=\underset{\mathbf{x}\in {\mathbb{R}}^N}{argmin}{\parallel \mathbf{y}-\mathbf{A}\mathbf{x}\parallel }_2+\mu {\parallel \mathbf{x}\parallel }_1.$$

(8)

Here, a measurement vector $\mathbf{y}$ is the difference in Equations (4)–(6), or (7) corresponding to $\Delta \mathbf{D}$, $\Delta {\mathbf{T}}_{\mathbf{D}}$, $\Delta \mathbf{R}$, or $\Delta {\mathbf{T}}_{\mathbf{R}}$, respectively. The corresponding sensing matrix $\mathbf{A}$ is ${\mathbf{A}}_{\mathbf{D}}$, ${\mathbf{A}}_{{\mathbf{T}}_{\mathbf{D}}}$, ${\mathbf{A}}_{\mathbf{R}}$, or ${\mathbf{A}}_{{\mathbf{T}}_{\mathbf{R}}}$, respectively. $l_0$-norm of $\mathbf{x}$ which counts the number of non-zero elements in $\mathbf{x}$ is originally used in CS to impose the sparsity on the recovery solution $\widehat{\mathbf{x}}$. $l_0$-norm of $\mathbf{x}$ is relaxed to $l_1$-norm of $\mathbf{x}$ with the assumption of a sufficiently low coherence of the sensing matrix for the spare solution, which enables convex optimization. Regularization parameter $\mu$ controls the sparsity of the solution. When $\mu$ is small, the sparsity is less important, and this should result in overfitting of $\widehat{\mathbf{x}}$ to the measurement vector. The opposite occurs as $\mu$ increases.

The SSP can be obtained from the mean of the solution for the ray information after a separate application of CS. Note that an unknown parameter $\mathbf{x}$ is common for all linear transformation. To exploit the property of the commonality, the conventional optimization in CS as Equation (8) is modified as follows:

$$\widehat{\mathbf{x}}\left(\mu \right)=\underset{\mathbf{x}\in {\mathbb{R}}^N}{argmin}\parallel \Delta \mathbf{D}-{\mathbf{A}}_{\mathbf{D}}{\mathbf{x}\parallel }_2+{\parallel \Delta {\mathbf{T}}_{\mathbf{D}}-{\mathbf{A}}_{{\mathbf{T}}_{\mathbf{D}}}\mathbf{x}\parallel }_2+{\parallel \Delta \mathbf{R}-{\mathbf{A}}_{\mathbf{R}}\mathbf{x}\parallel }_2+{\parallel \Delta {\mathbf{T}}_{\mathbf{R}}-{\mathbf{A}}_{{\mathbf{T}}_{\mathbf{R}}}\mathbf{x}\parallel }_2+\mu {\parallel \mathbf{x}\parallel }_1.$$

(9)

Equation (9) is a convex problem, which can be solved with the CVX tool [24]. An optimal value from Equation (9) is supposed to be simultaneously fitted to the measurement vectors. However, since the ray travel distance is roughly proportional to the corresponding travel time multiplied by the speed of sound (3 orders of magnitude), a scale of inherent error of the higher order term in the approximation is different according to the ray information, and the optimization of Equation (9) is inclined to reduce the errors in the depth and range more than those in the vertical and horizontal travel times. In order for each ray information to equally influence the optimization process, each $l_2$-norm in Equation (9) is normalized by the $l_2$-norm of the corresponding measurement vector as follows:

$$\widehat{\mathbf{x}}\left(\mu \right)=\underset{\mathbf{x}\in {\mathbb{R}}^N}{argmin}\frac{\parallel \Delta \mathbf{D}-{\mathbf{A}}_{\mathbf{D}}{\mathbf{x}\parallel }_2}{{\parallel \Delta \mathbf{D}\parallel }_2}+\frac{\parallel \Delta {\mathbf{T}}_{\mathbf{D}}-{\mathbf{A}}_{{\mathbf{T}}_{\mathbf{D}}}{\mathbf{x}\parallel }_2}{\parallel \Delta {\mathbf{T}}_{\mathbf{D}}{\parallel }_2}+\frac{\parallel \Delta \mathbf{R}-{\mathbf{A}}_{\mathbf{R}}{\mathbf{x}\parallel }_2}{{\parallel \Delta \mathbf{R}\parallel }_2}+\frac{\parallel \Delta {\mathbf{T}}_{\mathbf{R}}-{\mathbf{A}}_{{\mathbf{T}}_{\mathbf{R}}}{\mathbf{x}\parallel }_2}{\parallel \Delta {\mathbf{T}}_{\mathbf{R}}{\parallel }_2}+\mu {\parallel \mathbf{x}\parallel }_1.$$

(10)

Here, $\mathbf{y}/{\parallel \mathbf{y}\parallel }_2={(\mathbf{A}/\parallel \mathbf{y}\parallel }_2)\mathbf{x}$ is solved instead of the original problem $\mathbf{y}=\mathbf{A}\mathbf{x}$, which does not affect the solution from the original problem. Thus, Equation (10) is used for reconstructing the SSP to exploit the ray information simultaneously.

3. Numerical Experiments: SWellEx-96

As shown in Figure 2, the numerical environment is identical to that of Bianco and Gerstoft [10], and similar numerical experiments are conducted to examine the performance of the compressive SSP estimation. The sediment geoacoustic parameters including sound speed, density, and attenuation denoted in Figure 2 are unnecessary using present compressive SSP inversion. The average of 26 CTD casts from SWellEx-96 is used as the reference SSP, and the ocean SSP to be recovered by the CS is synthetically generated using a shape function, such as Gaussian functions or EOFs using the 26 CTD casts. Geometry for experimental site is known a priori for the inversion, which includes source and array positions and bottom bathymetry.

The SSP is inverted using the formulation in Section 2, which requires beamforming results that indicate the incident angles and travel times of the arrivals. While the travel times can be replaced by the relative delay times from a specific arrival (e.g., the earliest arrival) with slight modification of Equations (5) and (7), the travel times are used in numerical experiments. Figure 4 shows vertical arrival structures at an array for the reference SSP and their corresponding DOAs. The synthetic acoustic signals are generated with the conventional ray model [11] for the array composed of evenly spaced 10 sensors from 10 m to 150 m. Matched filter is applied to the acoustic signals with the transmitted signal of 0.5–2 kHz LFM waveform to show the arrival structures clearly. After applying a simple delay-and-sum beamforming [11], DOAs are revealed in a time domain, which have much higher intensities compared to regions around them. It is helpful to refine the DOAs to exploit frequency-domain beamforming based on advanced techniques [14,15,16,17,18] and high-resolution channel impulse response [21] for the center of the array simultaneously with the delay-and-sum beamforming result.

Practically, the DOAs should be used with the compressive SSP inversion. However, in this work, the inversion scheme after acquiring the DOAs is focused and examined. Eigenrays connecting the source to the receiver are shown in Figure 5a. Their arrival angles are overlaid on the DOAs from the beamforming according to their travel times, which are marked with the circles in Figure 5b. In particular, later arrival angles match well with the DOAs, and thus, instead of the DOA estimation via beamforming, the arrival angles and travel times of eigenrays are used for a convenient numerical experiment to examine the compressive SSP inversion while they are not perfectly identical to those from the beamforming; actually, DOAs and corresponding travel times from Figure 4 have inevitable deviations from the actual values owing to plane-wave approximation in the beamforming and uncertainties in experimental environments and the proposed scheme will be applied to more realistic situations after examining its performance with exact arrival angles and travel times from the ray model. On the other hand, indistinctive DOAs observed around 1.35 s correspond to refracted rays (dashed lines in Figure 5a) and are discarded during the inversion.

The steps for the inversion is summarized as follows: first, the beamforming is conducted for detection of DOAs (in this work, DOAs are replaced with the arrival angles and travel times of eigenrays from the ray model). Then, the detected arrivals are respectively backpropagated with their arrival angles from the center of the array for the reference SSP. The arrivals inappropriate for the approximation (e.g., the refracted rays in Figure 5a) explained in Section 2 are removed and the first derivatives in the sensing matrices are calculated for the remaining arrivals. Finally, the ocean SSP is inverted (or the deviation of SSP from the reference SSP is estimated) with information of the remaining arrivals using CS.

Note that information of arrivals is used for the inversion. When the ocean noises are dominant over acoustic signals from the source, they mask the arrivals after beamforming and the proposed scheme cannot be applied. Thus, in the following simulations, SNR sufficient to identify the arrivals is assumed.

3.1. SSP Dictionary: Gaussian Shape Functions

At first, the ocean SSP is generated by adding two Gaussian functions both with the same standard deviation as shown in Figure 6a. The peak positions of the Gaussian functions are 63.3 m and 126.7 m, respectively, and the peak values are +2 m/s and as −2 m/s, respectively, as in the paper [10].

The SSP dictionary is composed of Gaussian functions with the same peak/width and different peak positions. The peak is 1, the standard deviation related to the width is 5, and the peak position is from 10 m to 150 m every 1.5 m.

The sound speed fluctuation from the reference SSP is estimated using Equation (10) with the SSP dictionary composed of 94 shape functions ($N=94$). The inversion uses 12 rays ($M=12$). The rays with DOAs less than ${30}^{\circ }$ are used, and the rays for which the turning points are within the water column are excluded. Thus, 12 is the number of eigenrays from the ray model satisfying the aforementioned conditions. Practically, the number of rays used for the inversion should be determined with the DOA estimation as Figure 4b. $\mu$ is chosen as $0.01$. Note that the two Gaussian functions used to generate the ocean SSP do not match those in the SSP dictionary, and perfect recovery of the ocean SSP is impossible. The activated components of $\mathbf{x}$ to reconstruct the ocean SSP and the sound fluctuation by the activated components are shown in Figure 7. The solid and dotted lines are respectively the true and estimated values. There is a notable difference because of the basis mismatch and the inherent error terms in the approximation. The correlation coefficient between the solid and dotted lines in Figure 7b is $0.43$, which is calculated to evaluate the similarity between the true and estimated SSPs. However, components near the true values (Figure 7a) are selected for the estimation, which results in a modest correspondence between the true and estimated sound speed fluctuations in terms of their structures. The two peak positions in the estimated SSP are shifted upward from those in the true SSP (Figure 7b). Figure 8a,b show backpropagated ray trajectories from the receiver for the reference SSP and the estimated SSP, respectively. Backpropagated rays for the reference SSP do not converge exactly at the source position and one ray highlighted by the circles misses the source position greatly (Figure 8a). While the estimated SSP does not perfectly match the true ocean SSP, the backpropagated rays for the estimated SSP converge to the source position, including the ray highlighted by the circles, as shown in Figure 8b (see the inset parts of Figure 8).

3.2. SSP Dictionary: EOFs

A more realistic ocean SSP is generated by using EOFs based on the 26 CTD casts from SWellEx-96. The first, third, and seventh components out of 26 EOFs composing the SSP dictionary ($N=26$) are used to produce the ocean SSP (Figure 6b), and their coefficients are +6, −6, and +6, respectively, as in the paper [10]. The SSP fluctuates over the whole depths, and some arrivals deviate far from the source vertically at the source range, unlike in the previous case when the ocean SSP is replaced by the reference SSP. The arrivals for which the DOAs are less than ${30}^{\circ }$ are used to recover the SSP. An arrival with a turning point within the water column or arrival for which the vertical difference from the source at the source range is more than 10 m is excluded for the inversion. The number of arrivals satisfying the conditions and $\mu$ used for the inversion are 11 ($M=11$) and $0.01$, respectively. Figure 9 shows the inversion results and true values. Six components are used to reconstruct the SSP, and the largest three of these correspond to the actual components (Figure 9a). The inherent error terms in the linear approximation are attributed to an imperfect reconstruction of the SSP. However, the recovered sound speed captures the major features of the true SSP, as shown in Figure 9b, and the correlation coefficient between the sound speed fluctuations is $0.81$.

Figure 10 shows the ray trajectories corresponding to the reference and estimated SSPs. Differences of ray trajectories between the reference and estimated SSP are more significant than the previous case. Rays are not focused on the source position, as shown in Figure 10a. However, when the SSP is rectified by compressive SSP estimation (Figure 10b), 14 arrivals out of 15 arrivals, of which the DOAs are less than ${30}^{\circ }$, come closer to the source. Most of the backpropagated rays for the reconstructed SSP converge at the source position, and the rays highlighted with circles are notably enhanced by the recovered SSP.

By comparing Figure 7 and Figure 9 with Figure 2 and Figure 4 in Bianco and Gerstoft [10], the proposed scheme shows almost equivalent performance with their method. Note that the proposed scheme does not require geoacoustic parameters of sediment and replicas at VLA for SSP inversion. Thus, when the geoacoustic parameters are unknown or the proper propagation model is unavailable for an experimental environment, the proposed scheme can be applied with identified arrivals from beamforming to estimate an SSP.

Now, consider a situation where the DOAs are estimated with errors. To investigate the effects of errors on the performance of compressive SSP inversion, Equation (3) is slightly modified with an additional Taylor expansion for DOA as follows:

$$\begin{array}{cc}\hfill T_D(\mathbf{x},{\theta }_{tr})=& T_D(\mathbf{0},{\theta }_{est})+\sum _{i=1}^N\frac{\partial T_D(\mathbf{0},{\theta }_{est})}{\partial x_i}{x}_i+{ϵ}_{T_D}(\mathbf{x},{\theta }_{est})\hfill \\ & +\left[\frac{T_D(\mathbf{0},{\theta }_{est})}{\partial \theta }+\sum _{i=1}^N\frac{{\partial }^2{T}_D(\mathbf{0},{\theta }_{est})}{\partial x_i\partial \theta }{x}_i\right]({\theta }_{tr}-{\theta }_{est})+\dots ,\hfill \end{array}$$

(11)

where ${\theta }_{tr}$ and ${\theta }_{est}$ are respectively true arrival angle and estimated arrival angle including an error by array tilt, etc. Note that extra error terms emerge, which are proportional to the error in DOA. The same error analysis can be conducted on the other ray information used for the inversion.

Figure 11 shows comparison results between true and estimated values in terms of sound speed difference along depth from the reference SSP when exact DOAs are contaminated with uniformly distributed error angles between $-{0.05}^{\circ }$ and $+{0.05}^{\circ }$. The same SSP in Figure 9 is used. First, the ocean SSP is estimated with compressive SSP inversion neglecting the errors in DOAs as before (Figure 11a). As expected, the CS performance deteriorates owing to the extra error terms. However, the estimated SSP still captures major sound speed variation near sea surface and corresponding ray trajectories tend to merge at the source position as in Figure 10b.

A simple method to mitigate the error effects on the inversion performance is to estimate the errors in DOAs as well as SSP (or $\mathbf{x}$) during the inversion process; ${\theta }_{tr}-{\theta }_{est}$ is treated as an unknown parameter. The terms related to the errors in DOAs are added to Equations (4)–(7), respectively. Here, the first derivative with respect to $\theta$ is calculated numerically as the first derivative with respect to $x_i$ whereas the second derivative terms are ignored. The measurement vectors ($\Delta \mathbf{D}$, $\Delta {\mathbf{T}}_{\mathbf{D}}$, $\Delta \mathbf{R}$, and $\Delta {\mathbf{T}}_{\mathbf{R}}$) are still linearly approximated with unknown parameters related to the SSP and errors in DOAs. The sparsity condition is imposed on $\mathbf{x}$ as before. On the other hand, $l_2$-norm of the errors in DOAs multiplied by regularization parameter is added to Equation (10) and it prevents the errors from growing too large during the inversion. The CVX tool [24] is used for finding the unknown parameters. By the simple consideration of the errors in DOAs, the CS inversion performance is improved as shown in Figure 11b; the correlation coefficient between true and estimated sound speed increases from $0.4$ to $0.64$.

The errors in DOAs are increased to be between $-{0.5}^{\circ }$ and $+{0.5}^{\circ }$. The modified inversion scheme is applied to estimate the DOA errors as well as the SSP. Backprogated rays deviate more from the source position as the increment of the DOA errors and the number of available rays for the inversion is reduced to 5 ($M=5$) in the simulation for Figure 12 showing estimation results. While the errors are included in DOAs and the number of eigenrays used for the inversion is decreased, the estimated SSP still shows a good agreement with the true SSP (Figure 12a). The error effects on the SSP inversion is diminished with reliable error estimation for DOAs (Figure 12b). It is beneficial for reconstructing the SSP using beamforming results.

To enhance the CS inversion performance, the errors in DOAs should be suppressed via refining beamforming results through advanced schemes [14,15,16,17,18,21]. Furthermore, the sparsity in the solution can be exploited with probability distribution of the errors as well as its probability distribution (i.e., sparse Bayesian learning) [25]. However, here, the focus is on the inversion process after accurate DOA detection and the advanced schemes on alleviating the error effects are beyond the scope of this work.

3.3. Examination of the Compressive SSP Estimation

To examine the performance of the compressive SSP estimation, 500 synthetic SSPs are generated using the EOFs. Three randomly-selected components out of the 26 EOFs are activated for each SSP, and a half-normal distribution is used for the random selection to consider relative importance of the EOFs. The amplitude of the activated component is +6 or −6.

Three experiments are conducted with different objective functions for the compressive SSP inversion, with 500 SSPs used for each experiment. In the first experiment, part of the information related to the positions of the rays (vertical and horizontal position differences) is used for the objective function. Then, all of the ray information is exploited for the objective functions of the second and the third experiments without and with normalization, respectively.

Similarity between the true and estimated SSPs is evaluated with their correlation coefficient. Figure 13 shows a portion of the cases corresponding to the correlation coefficients on the x-axis. No errors are in arrival angles and travel times. As in the previous cases, the rays with DOAs less than ${30}^{\circ }$ are used for the inversion, and the ray with the turning point in the water column or vertically far from the source at the source range is excluded. The number of available eigenrays for the inversion is different according to the randomly generated SSP. The average number of rays used for the inversion is $10.5$. A proper $\mu$ has to be chosen for each SSP to provide better performance for the compressive SSP estimation. However, for convenience, $0.005$, $0.05$, and $0.01$ are used for all cases in the first, second, and third experiments, respectively. The most frequent cases in the first experiment have correlation coefficients between $0.8$ and $0.9$ (Figure 13a). The cases with correlation coefficients of more than $0.7$ occur 302 times out of the 500 cases during the first experiment. The compressive inversion using limited information regarding the position shows a fairly good performance, and this performance can be improved by exploiting all ray information, as shown in Figure 13b, where normalization is not applied (Equation (9)). The correlation coefficients of most frequent cases in the second experiment are between $0.9$ and 1, and the number of cases for which the correlation coefficients are over $0.7$ increases slightly (316 times). In the second experiment, the information of the position is more important than that of the travel times in the compressive SSP estimation. To use all information equally, normalization is conducted on the objective function, as in Equation (10). Figure 13c shows the result from the third experiment, which demonstrates much better performance than the other two experiments. The cases with correlation coefficients over $0.9$ appear more than twice compared to those in the second experiment. More than 80 percent of the cases have a correlation coefficient over $0.7$.

Mean squared errors (MSE) of estimated SSPs from true values are also calculated according to objective functions and the number of available rays (Figure 14). As expected, inversion using all ray information after normalization outperforms the others. Surprisingly, the reduction in the number of available rays insignificantly affects on the performance of our inversion scheme except peaks near the sea surface when it is estimated with MSE.

In this work, underdetermined linear system is solved with CS. For better CS performance, it is desirable to have more arrivals, making fat sensing matrices thinner. However, to investigate the effect of arrival number on its performance, arrivals with angles less than ${20}^{\circ }$ are used to estimate the ocean SSP while satisfying other necessary conditions. The average number of rays used for the inversion is $5.5$. Inversion performances for the three experiments with different objective functions deteriorate in terms of similarity (Figure 15). Number of cases with correlation coefficients over $0.7$ drops from 302 to 252 cases for the first experiment, and from 316 to 298 cases for the second experiment, and from 411 to 381 cases for the third experiment. On the other hand, the reduction of available ray number has the slight effect on the performance when it is evaluated with MSE (Figure 14). Regardless of the estimation methods on the performance, the inversion exploiting all ray information equally demonstrates better results.

The inversion scheme delineated in this work cannot be applied when DOAs are unavailable. However, the inversion using all ray information shows a good performance in the numerical environment with relatively few arrivals.

3.4. Application: SSP in SWellEx-96

In previous simulations, incident grazing angles and travel times of arrivals from ray model are used and they are ideal for the CS SSP estimation. Practically, to obtain them, a beamforming is applied. DOAs and corresponding travel times from the beamforming are not exactly the same as the actual values owing to plane-wave approximation in the beamforming. Furthermore, the difference is increased by uncertainties in experimental environment including source/array position, array tilt, and bottom depth, etc.

To demonstrate the feasibility of the proposed scheme to real application, peak values of beamforming results as in Figure 4b are used to recover one of SSPs in SWellEx-96 as shown in Figure 16a; the SSP has the largest deviation from the mean SSP. Eleven arrivals less than ${30}^{\circ }$ are detected and two arrivals are filtered out after backpropagation. The estimated DOAs and corresponding travel times of nine arrivals whose backpropagated rays are within 10 m from the source position are used for the inversion. As shown in Figure 16b, the sound speed difference along depth evaluated with the modified scheme shows a good agreement with the true values while it overestimates major SSP variation near the sea surface. To account for the uncertainties in experimental environments, random errors between $-2^{\circ }$ and $+2^{\circ }$ are added to the DOAs, which make the backpropagated rays deviate far from the source at the source range. Here, to acquire the arrivals used for the SSP inversion, the criterion increases from 10 m to 20 m. The number of available arrivals is reduced to six and their backpropagated rays are located farther from the source position. Still, the estimated SSP are in good agreement with the true values by capturing the major variation in SSP near the sea surface and following the overall trend.

4. Conclusions

In this work, the sound speed profile (SSP) is represented by a linear combination of shape functions in an SSP dictionary. The range/depth and horizontal/vertical travel time of arrival are linearly approximated from the reference SSP using Taylor series with the shape function coefficients at the corresponding arrival angle. Compressive sensing (CS) is applied to estimate the ocean SSP with the linear approximations for all direction of arrivals (DOAs). Thus, the inversion scheme proposed in this work can be used when the DOAs are resolved. The objective function is normalized to equally consider all the ray information.

The performance of present compressive SSP estimation is examined for the SWellEx-96 experimental environment. For simplicity, the beamforming results needed for the inversion are calculated using a conventional ray model by treating arrival angles from the ray model as DOAs, which are assumed to be identifiable. Five hundred sound speed profiles are generated synthetically using empirical orthogonal functions. Three numerical experiments are conducted using different objective functions. In the first experiment, the information on the position of the arrival is used. All ray information on the position and travel times is used without normalization and with normalization in the second and third experiments, respectively. The performance of the compressive SSP estimation is enhanced by simultaneous usage of position and travel time of ray. In particular, in the third experiment, over 80 percent of 500 cases show a correlation coefficient exceeding $0.7$. The performance is degraded with decrease in the number of arrivals when it is estimated with the correlation coefficient. However, the inversion scheme using all ray information shows a reliable result with relatively few arrivals.