Image source detection for geoacoustic inversion by the Teager-Kaiser energy operator^a) Open Access

The image source method (ISM) is a recently developed technique for geoacoustic inversion¹ which provides the sound speed profile locally with good resolution and with low computational cost compared to classical approaches (see, e.g., Refs. 2 or 3). This local inversion strategy requires only one measurement of a signal emitted by a broadband source and recorded by an array of hydrophones. The effectiveness of the ISM has been shown on synthetic and real data.^1,4 However ISM is conditioned by some parameter settings. For further improvement of the computational time and the reduction of false positive sources, a new time arrival estimation strategy is introduced. This strategy is based on the Teager-Kaiser energy operator (TKEO)⁵ which is a differential operator well-dedicated for the detection of signal peaks such as those reflected by the seafloor. Using the arrival times, estimated by the TKEO, arrival angles are calculated by a triangulation method. Once the times and the angles of arrival are estimated, the velocity profile is calculated using the same inversion strategy as described in Ref. 1.

ISM inversion is based on the analysis of an acoustic wave emitted by a point source and reflected by a layered seafloor¹ where only one array recording of the seabed reflected path for one shot of the source is used [Fig. 1(a)]. Under the Born approximation, the reflected signal can be modeled as a sum of contributions coming from image sources which are symmetric to the real source with respect to the seabed layers [Fig. 1(b)]. The location of these image sources is related to the sound-speed profile; it is therefore possible to invert this structure by detecting the image sources.

The first step of ISM is the detection and location of the image sources. In the original work, two maps are reconstructed from the recorded signals: The migration map and semblance map.¹ The combination of these maps provides a location of the image sources which is accurate enough for inversion. This method is later called ISM-MS. The maps are processed in a homogeneous medium, water, so the image sources are not aligned on the source's nadir because the refraction is not taken into account [Fig. 1(b)]. Nevertheless, times (⁠$ti, j$⁠) and angles (⁠$θi, j$⁠) of arrival from each image source $i$ to each hydrophone $j$ of the array are correct.

At the resolution of the signal which is roughly the wavelength, the number of layers is directly given by the number of detected image sources. The first image source, corresponding to travel only in water, can be used to obtain the system's geometry.⁴ Then, the arrival time and angle of the second image source give the sound speed and thickness of the first layer using Snell-Descartes laws.⁴ In a recursive way, these quantities are then obtained for all the layers.

In this section, we show that the $ti, j$ and $θi, j$ parameters can be determined in a simple way with a reduced computational time and with good precision using a strategy based on the TKEO.⁵ This method is later called ISM-TK.

TKEO is a nonlinear tracking energy operator and its output to a given signal $x(t)$ is the actual physical energy required to generate $x(t)$⁠.⁵ It has mainly been used for its demodulation properties of signals⁶ but rarely for its self-properties.⁷ In this work the TKEO output is exploited for detection purposes. In continuous time, it is defined as

For the discrete time signal $x(n)$⁠, the TKEO can be approximated as follows:^6,7

Equation (2) shows that the TKEO computes a running estimate of the signal energy at each instant that takes into account the signal's strength at its immediate neighbors. An important aspect of the TKEO is that it amplifies discontinuities and sudden changes in amplitude while the smooth transitions between the samples are reduced. This property is mainly attributed to the use of second order derivative and is exploited here for the estimation of arrival time by identification of peaks of TKEO output.

The TKEO-based detection is illustrated on simulated and real data (Sec. IV). Synthetic signals are computed by a numerical evaluation of the Sommerfeld integral with the Born approximation in a stratified fluid. The signals are corrupted with white Gaussian noise of SNR $≃$ 3 dB and deconvoluted by a cross-correlation with the source signal. The peaks of the TKEO output (Fig. 2) correspond to arrival times but can also be attributed to noise. To overcome this problem, a threshold, T, is used to isolate the most prominent peaks that are considered as the reflexion of the waves on the surface. More precisely, this threshold is introduced to minimize the missing of true peaks, while reducing the number of false detected of peaks introduced by noise, within a reasonable limit. For this, the T value is taken as a scaled version of the average of the means of the TKEO output over all the hydrophones,

where

where $xk$ is the signal received by the kth hydrophone, $K$ is the number of hydrophones, $N$ is the number of samples, and $C$ is a scaling factor. This scaling value depends essentially on the input signal-to-noise ratio (SNR). For the present synthetic signal the true peaks are successfully detected with $C$ set to 0.9 (Fig. 2). Simulation results also show that the TKEO-based approach is also effective for synthetic signals with linear chirp and impulsion sources, horizontal and vertical array, and at low SNR. For a more robust estimation of arrival time, the identified peaks of hydrophones are aligned to those of the first hydrophone. This step ensures that each peak corresponds to the arrival time of the associated image source. A good resolution of arrival time values is required for a better estimation of the parameters. This method is quite simple, but still sensitive to errors of arrival time estimation. Consequently, for better resolution of peak position an oversampling of the deconvoluted signal is necessary. In practice oversampling of noisy signals can induce fluctuations on peaks which in turn complicate their detection. To remedy this problem, a smoothing strategy based on the Savitzky-Golay filter is used. This filter performs a local polynomial regression to determine the smoothed value for each data point.⁸ This smoothing is superior to the moving average because it preserves the characteristics of the data such as height and width of peaks, which are usually lost by moving average.

In this work, $θ$ estimation is given by a triangulation method. The calculated arrival times $tij$ (Sec. III B) are converted to distances $dij$ while keeping the water sound-speed $c0$⁠. Since the sound-speed profile is unknown, $dij$ is the distance between hydrophone $j$ and an equivalent image source, located in the same direction as the arrival angle of the real image source [Fig. 1(b)]. But unlike the real image source, the equivalent image source locations depend on the hydrophone's position [Fig. 3(a)]. Instead of seeing a single source, the array perceives an extended source. This could be interpreted as a blurring process due to the sound-speed profile. To obtain an accurate arrival angle for inversion, the intersections located near the nadir of the sources of the pairs of distances $dij$ and $dik$ between image source $i$ and hydrophones $j$ and $k$ are first computed. This produces $M(M−1)/2$ locations with $M$ being the number of hydrophones. The image source's position is taken as the median of all these locations (Fig. 3). The median of the hydrophones' positions is similarly calculated. Finally, for inversion, the distance $di$⁠, the arrival time $ti$⁠, and the arrival angle $θi$ of image source $i$ are computed from these two median locations.

The configuration used is similar to that of the data acquired in the Clutter 09 experiment⁴ where an omnidirectional source and a horizontal array of 32 hydrophones with a spacing of 1.05 m are towed by an Autonomous Underwater Vehicle 12 m above the seafloor. The source transmits a linearly frequency-modulated pulse, in the band 1600 to 3500 Hz, of 1 s duration. The bottom is composed of two layers of homogeneous thickness and of constant sound speed (Table I), over a basement. To accurately estimate the arrival time, the signal emitted by the source must be a zero-phase pulse. Thus, the same preprocessing is applied to find the emitted signal for the zero-phase impulse response. The estimated geoacoustic parameters from the signals by the ISM-TK and ISM-MS methods are reported in Table I. These results show that, overall, the ISM-TK method performs better in terms of estimation of thickness and sound speed parameters compared to the ISM-MS method. Furthermore, the computation time is reduced by a factor of 10.

The ISM-TK method is tested on data acquired near Elba Island in Italy as part of the SCARAB (Scattering And ReverberAtion from the sea Bottom) experiment series.⁹ The data are acquired from site 2 where the water depth is 150 m and from side-scan sonar data. The seabed is flat and featureless; bottom slopes are 0.3 or less. The source is 200 m away from a vertical array of 64 m in length and moored to the seafloor. This array is composed of 15 hydrophones where the lowermost one is around 12 m above the seafloor. Due to ocean currents, the array can be distorted. Thus, array inclination must be taken into account and more specifically the experiment's geometry. This problem is illustrated by Fig. 4 where position error is estimated at 18 to 19 m. The array's correction step is performed with the two arrival times detected with the TKEO on the direct path and the first reflected signal. In this experiment, the signal of hydrophone 5 was not exploitable and thus it was discarded. Estimation results of the ISM-TK, the ISM-MS, and Holland and Osler methods⁹ are compared and are reported in Table II. Since the SNR of the real data is unknown, the $C$ value is set to 1. The ground truth, provided by coring, is only available for the first 15 m of the sediment; comparisons are restricted to this depth. Since the Holland and Osler method compares favorably against the ground truth,⁹ it is used as the reference. Overall, our results confirm the findings of Holland and Osler (Table II). However, only four layers are identified by the ISM-TK method. Also, the first layer has not been identified due to the complexity of the sedimentary stratification of the first meters. However, from the computational point of view, the ISM-TK method is very fast compared to two of the other methods.

In this letter, a new arrival time estimation strategy based on the TKEO is introduced. Arrival times corresponding to each image source are identified as time locations of prominent peaks detected directly on the recorded acoustic signals by the TKEO. The conjunction ISM and TKEO (ISM-TK) is illustrated on both synthetic and real data. Reported results on synthetic data show that ISM-TK performs better than IMS-MS in terms of thickness and sound speed parameter estimations and with significantly reduced computational cost by a factor varying from 10 to 50 depending on the medium complexity. The effectiveness of ISM-TK is also illustrated on real data and the obtained inversion results are, overall, in good agreement with the ground truth. Based on the TKEO, which is an instantaneous operator, the computation burden of the inversion is largely reduced which makes the ISM-TK a good candidate for real time inversion.

The authors would like to thank Professor Holland from Penn State University for making real data available.