A Volunteer Computing Project for Solving Geoacoustic Inversion Problems

1 Introduction

Desktop Grid computing [11] uses resources of geographically distributed personal computers, connected over general-purpose network. Huge increase in personal computers’ performance has made Desktop Grid computing a powerful tool for solving numerous scientific problems. Volunteer computing [16] is a special kind of desktop grid computing. Its defining characteristic is that it uses computational resources of private persons. Usually, one volunteer computing project is designed to solve one or several closely related hard problems. Volunteer computing is very cheap — to maintain a project one needs only a dedicated server (connected to Internet) and several client applications. Volunteer computing suits well to solve hard problems, which can be decomposed into a set of independent subproblems. In last two decades, a lot of important results in various areas (astronomy, medicine, mathematics, etc.) were obtained in volunteer computing projects.

There exists a number of classes of computationally intense geoacoustic inversion problems, which can be solved using embarrassing parallelism. This is the main reason why we launched the volunteer computing project Acoustics@home, aimed at solving such problems. The project is based on the popular platform BOINC (Berkley Open Infrastructure for Network Compuitng [3]). The notion of geoacoustic inversion refers to a variety of techniques in underwater acoustics that can be used for the reconstruction of water column and bottom parameters from acoustic data [20]. While previously the data for the geoacoustic inversion was mostly obtained using expensive receiver arrays, recently it was shown that a single-hydrophone recording of a broadband pulse signal can be also successfully used for estimating the acoustical parameters of sea bottom [8,9,10]. Due to the lack of the spatial diversity in single-hydrophone measurements, the inversion procedure in this case relies on the so-called waveguide dispersion, i.e., the dependence of arrival times on frequency and mode number (see details, e.g., in [9, 23]). Some interesting recent results [25] also indicate that single-hydrophone dispersion-based inversion method outlined in [9] can be used for the estimation of a sound-speed profile in a shallow-water waveguide. The implementation of this method in practice can be thought of as a solution of a minimization problem in a (very large) discrete search space [26], and every evaluation of the cost function requires numerous solutions of an acoustic spectral problem [23, 26]. Thus, the whole computational burden can be easily divided into a large number of relatively simple independent tasks, offering many opportunities for the application of parallel computing.

The study has the following structure. In Section 2 we outline the dispersion-based geoacoustic inversion scheme following [8, 9, 23, 26]. The synthetic acoustic data used in our work is discussed in Section 3. In Section 4 we describe the project Acoustics@home and show the results of computational experiment, performed within its framework. In Section 5 it is discussed in which way we plan to improve Acoustics@home. In the rest of the paper we discuss related work and draw conclusions.

2 Dispersion-based geoacoustic inversion

The dispersion-based technique for the geoacoustic waveguide parameters estimation [9] introduced by J. Bonnel and his colleagues several years ago is currently a very dynamic field of research in underwater acoustics [8, 10, 23, 25]. In this technique, the dispersion data is used as an input of the inversion algorithm. By the dispersion data we understand the set of functions τ_m, m = 1, 2, … M, where τ_m(f) denotes the arrival time of m-th modal component [19, 23] of a pulse acoustical signal as a function of frequency f (M is the total number of modes that can be filtered from the data). The curves t = τ_m(f) in the two-dimensional time-frequency space are called dispersion curves [20].

The experimental dispersion curves can be obtained from a pulse signal recorded by a single receiver (hydrophone) [8, 10]. Usually, various approaches to the mode filtering are used to extract the dispersion curves from a time series obtained in an experiment. The rise of the dispersion-based inversion schemes popularity can be associated with the so-called warping transform [5, 8] offering a very elegant way to separate the modal components of a pulse signal. At the same time, for any fixed model of a waveguide dispersion curves can be computed theoretically by solving acoustic spectral problem [19, 26]. The mismatch between the experimental and theoretical arrival times indicates to what extent the given model of the waveguide is consistent with the observation results. Normally the waveguide model includes certain number of unknown parameters, and the geoacoustic inversion serves the purpose of estimating their values. The set of parameters corresponding to the minimal mismatch determines the theoretical model that is the most adequate to the experimental data. Hereafter this set is referred to as media parameters estimate based on the given experimental data. Note that probably the most important advantage of the dispersion-based inversion schemes is their ability to provide some information on the medium using single-hydrophone measurements. By contrast, more conventional geoacoustic inversion methods [20] usually require the deployment of a vertical or horizontal array of receivers (which makes the experiment much more expensive and complicated). The lack of spatial diversity of the measurements in this case is compensated for by the frequency diversity.

The geoacoustic inversion problem can be transformed into a problem of minimization of a certain mismatch function [26]. The simplest natural choice for such function is the standard mean square fitness function measuring the average squared discrepancy between the theoretical arrivals τmth(f, A) computed for the parameters vector A and experimental arrivals τmexp(f):

The minimization is performed over a certain domain of parameters values which is typically a cuboid in a N_p⁻ dimensional Euclidean space, where N_p is the number of the waveguide parameters being inverted. The cuboid boundaries are determined from certain physical considerations.

3 Waveguide model and the synthetic data

In the present study, the geoacoustic inversion problem is considered in the case of a homogeneous two-dimensional waveguide Ω = {(x, z)|0 ≤ z ≤ H}, where z denotes depth, and x is horizontal coordinate. The waveguide consists of the water column 0 ≤ z ≤ h and a single bottom layer h ≤ z ≤ H. The sound speed profile in the water is continuous function c = c(z) (the dependence of sound speed on z).

For simplicity the sound speed profile is assumed to be independent of x. Although in real ocean this condition is never fulfilled exactly, the variation of sound speed in range x is usually much slower than in depth z. It is often reasonable to assume that the sound speed profile is the same for all x in a certain area of a shallow sea containing the source and the receiver. In our simplified model the bottom is represented by a fluid halfspace with the sound speed c_b and density ρ_b. Although this may also look somewhat non-realistic, the situation where the bottom depth h varies by few meters over a track of several kilometres is quite common in underwater acoustics. In such cases it is reasonable to assume the bottom depth to be constant. A more complicated inversion algorithm where the variations of the latter are taken into account is described in [23].

Our goal in this study is to evaluate the possibility of the sound speed profile c(z) reconstruction from the dispersion data. The results from [25] show that certain approximation of the function c(z) can be obtained by the dispersion-based inversion (although in [25] it is not really perfect even for synthetic data). In the present work we show that at least in principle it can be inverted for very accurately (when inverting for the sound speed profile alone and considering all other parameters known).

For our numerical experiments we prepared a synthetic data by solving the direct propagation problem as described in [24]. A chirp signal with instantaneous frequency varying from f₀ = 20 Hz to f₁ = 300 Hz was emitted by the source at z_s = 10 m. The sound speed profile in the water used in simulations is shown in Figure 1, while the bottom parameters were chosen as follows: c_b = 1700 m/s, ρ_b = 1.7 g/cm³. The time series of the signal was computed at R = 7 km from the source.

Figure 1

Estimated sound speed profiles c(z) for N_c = 1, 2, …, 5 and the true one (black line on all panels). Top panel: N_c = 1 (green line) and N_c = 2 (red line). Middle panel: N_c = 3 (green line) and N_c = 4 (red line). Bottom panel: N_c = 5 (green line) and N_c = 6 (red line).

We applied the warping-based method of the extraction of the dispersion curves from the simulated data. As a result, the dispersion curves τ_m(f) of the five modes were obtained (i.e., in our case M = 5).

4 Computational experiments in Acoustics@home

In order to solve computationally hard geoacoustic inversion problems, we launched the volunteer computing project Acoustics@home [2] on 28 March 2017. It is based on the BOINC platform [3]. In Acoustics@home all daemons (that operate on the server) and client application (it operates on volunteers’ PCs) are based on the CAMBALA MPI-program [26]. Work generator daemon decomposes an original problem into independent subproblems by varying several parameters of a search space. The rest of the parameters are varied by the client application. For each set of the parameters A the value of the mismatch function (1) is calculated by the client application. Program sources of the project are avaliable online [1].

As it was mentioned above, the goal of our computational experiments is the estimation of the sound-speed profile (SSP) c = c(z) in the water column (see Section 3). Hence, the set of the inversion parameters A consists of the sound speed values c₁, c₂, …, c_Nc on a uniform grid of the depth values z₀, z₁, z₂, …, z_Nc. The grid is introduced in the water column in such a way that z₀ = 0, z_Nc = h (z_i+1 − z_i = Δz = const). We also assume that c_Nc = c_Nc−1, and therefore N_c parameters should be inverted. For each z_j corresponding c_j in our experiment can take the values in a discrete set of even integers from 1450 m/s to 1510 m/s (31 possible values in total). Although it means that the estimation error can be ~ 2 m/s just due to the search space discretization, such accuracy is sufficient for many applications in underwater acoustics.

Clearly, the resolution of the SSP c(z) depends on the number of nodes N_c. The more nodes we can afford, the better is the possible accuracy of the SSP estimation. We launched our dispersion-based geoacoustic inversion algorithm for N_c = 1, 2, …, 6. The inversion procedure for N_c = 1, 2, …, 4 was successfully completed via CAMBALA on the computational cluster “Academician V.M. Matrosov” of Irkutsk Supercomputer Center of SB RAS [21]. In our experiments we used 1 node of the cluster, which is equipped with two 18-core CPUs Intel Xeon E5-2695 (36 CPU cores in total) and 128 gigabytes of RAM.

The search spaces for N_c = 5 and N_c = 6 turned out to be too large for a cluster. That is why we launched them in Acoustics@home on 12 September 2017. For N_c = 5 the search space was divided into 29 791 workunits, which were made by varying the values of c₁, c₂, c₃. The client application had to vary the values of c₄, c₅ for each workunit, so 961 points had to be processed within it. For N_c = 6 the search space was divided into 923 521 workunits (by varying the values of c₁, c₂, c₃, c₄), each of them consisted of 961 points (the values of c₅, c₆ were varied). On average it took about 1 hour to process one workunit on 1 CPU core. Both experiments were successfully completed, it took about 2 months to perform the hardest of them.

Recently Acoustics@home has reached the average performance of 7 teraflops, 2544 PCs of 1547 volunteers have already participated in the project. Let us remind, that during our experiments for N_c = 1, 2, …, 4 we used 36 cluster’s CPU cores. In the considered experiments for N_c = 5 and N_c = 6, the project’s performance was comparable to that of a computational cluster equipped with 500 CPU cores of the mentioned kind. Details of the conducted experiments are presented in Table 1. The estimated sound speed profiles are presented in Figure 1.

Let us briefly discuss the inversion results from the acoustical point of view. It should have been expected that in the case N_c = 1 we obtained some value of the sound speed that is close to its average over the whole water column. In the case N_c = 2 the approximation is still very inaccurate, yet it already catches the main features of the true profile. Indeed, to certain extent it reproduces the thermocline near the surface and the homogeneous layer near the bottom.

From the middle and the bottom panels of Figure 1 it is clear that starting from N_c = 3 we already obtain a relatively accurate estimate of the sound speed profile. On average the approximations corresponding to N_c = 4 and N_c = 5 are almost identically accurate. Thus, the grid of a few points is in fact sufficient to estimate the hydrological properties of the shallow-water waveguide via geoacoustic inversion. The results for N_c = 6 also indicate that the increase of the number of nodes does not necessarily improves the inversion accuracy. This can be considered an additional justification for using trans-D inversion technique [25], where the preference in given to the sound speed profile approximations with fewer nodes.

One of the most efficient ways to boost the performance of a BOINC-based project is to organize a competition on the BOINCstats site [7]. In the course of such competitions teams try to earn as many project’s credits as possible. In September 2017 the “Ocean of fantasy” competition was launched to help Acoustics@home perform the aforementioned experiments faster. As a result, 26% of workunits for N_c = 6 were processed within a week (let us remind, that the whole experiment took 58 days). At the end the competition, the project’s has reached the record performance of 18 teraflops. It should be noted, that competitions of the mentioned type can lead to high load of a project’s server. In order to minimize the server load, an improved task scheduling can be used (e.g., [22]). However, in our case the default task scheduling coped with the competition.

5 Preliminary experiments on GPU

In the future, we are planning to improve Acoustics@home in two respects. Firstly, it is intended to employ advanced black-box optimization algorithms. Secondly, a GPU-based version of the computing application will be implemented. In the following subsection, a new blackbox optimization algorithm is proposed, then we describe its GPU implementation, and present the experimental results obtained by this application on one geoacoustic inversion problem.

5.1 New Black-box Optimization Algorithm

The problem (1) belongs to the class of black-box optimization problems which are very common in practice. In such problems derivatives are either unavailable or very hard to compute. Thus, solution methods should not rely on derivatives. In this paper, in order to minimize the objective function (1) we use uniform mesh to obtain a finite search space. This search space is processed by exhaustive search (see Section 4). In [27] instead of exhaustive search we applied iterative hill climbing algorithm. In the future, we are planning to apply the combination of global and local search techniques. The global search techniques will be used to diversify the search, by starting the local search from several initial points. The local search method will be used to minimize the value of the objective function. We plan to consider several local search techniques: Hooke-Jeeves method [17], pseudo-gradient approach [15] and a variety of coordinate descent techniques.

In black-box optimization problems derivatives are either unavailable or very hard to compute. Thus, solution methods should not rely on derivatives. There are a lot of local search techniques which can be used for this purpose: Hooke-Jeeves method [17], pseudo-gradient approaches [15] and a variety of coordinate descent techniques. We developed a simple modification of the coordinate descent method presented in Algorithm 1

Algorithm 1

ASN Search

The proposed algorithm (ASN search) uses asymmetric neighborhood adaptation. It tries to improve the objective function f(x) by shifting by some value δ_i along each coordinate vector e_i in both directions. In the case of success (failure) the respective step size is increased (decreased) by multiplying it on α > 1(β < 1). When approaching minimum values of δ_i usually decrease. The method terminates when the maximal component of the vector δ becomes less that the predefined ε > 0.

The choice of initial δ and ε significantly affects the performance. For the problem under consideration the best performance were obtained with δ_i = 10⁻¹, i = 1 …, n, ε = 10⁻⁴.

6 Related work

The computing clusters are often used in the practical applications of geoacoustic inversion algorithms (see, e.g., [14] and references therein). While sometimes one can compromise by using heuristic optimization techniques in end-user applications, the problems of development and validation of inversion algorithms anyway set very strong demand for the high-performance computational tools.

In the previous works [26, 27] we solved on a computing cluster the problems similar to that, considered in the present paper. The restrictions on available computational resources forced us to launch Acoustics@home. But we will continue to use the cluster in order to solve simple scenarios, or to test new versions of the computing application.

In [18] general purpose GPUs were used to accelerate underwater acoustic propagation modeling.

7 Conclusions

In the present study we describe the volunteer computing project Acoustics@home that was set up to handle computationally intense geoacoustic inversion problems. With the help of our volunteers we conducted a series of experiments aiming to investigate the capabilities of the dispersion-based inversion method in producing the estimate of the sound speed profile in shallow-water environment. In our view, it seems interesting to find out, how many grid points are required to resolve the function c(z) and how accurate the inversion result can be. In [25] a similar question was studied by using a sophisticated trans-D Bayesian inversion approach. Our results are obtained in a much more simple and straightforward way, and they are also probably easier to understand. The accuracy of the sound speed profile inversion is quite encouraging, and it seems that in principle function c(z) can be estimated from the data almost exactly. At least, it is clear that there is no limitation neither on the side of the warping-based mode filtering algorithm, nor on the side of the global minimization problem (i.e., it features only one global extremum of the fitness function near the true set of parameters values). It is also important to note that increasing the number of nodes in sound speed profile approximation does not necessarily improve the inversion accuracy. At some point adding new dimensions of the search space results in the larger deviation of the inverted sound speed profile from the true one.

The present study highlights great opportunities that volunteer computing can bring into the field of computational underwater acoustics. In the future, we are aiming at developing a powerful computational framework that can be used on demand by any member of ocean acoustics research community who needs to conduct some very intense calculations in order to solve certain direct or inverse problems. To this end, we are planning to expand our code by adding some modules capable of solving propagation problems in inhomogeneous 2D and 3D waveguides and also implement other inversion algorithm (not necessarily relying on single-hydrophone data).