Towards automatically building starting models for full-waveform inversion using global optimization methods: A PSO approach via DEAP + Devito

TL;DR

This paper presents an automated approach to building initial models for full-waveform inversion by integrating PSO with DEAP and Devito, leveraging parallel computing to reduce computational costs.

Contribution

It introduces a novel combination of PSO, DEAP, and Devito frameworks for efficient macro-model estimation in seismic inversion.

Findings

Efficiently estimates subsurface velocity models using global optimization.

Reduces computational effort through parallelization with Dask.

Automates the initial model building process for FWI.

Abstract

In this work, we illustrate an example of estimating the macro-model of velocities in the subsurface through the use of global optimization methods (GOMs). The optimization problem is solved using DEAP (Distributed Evolutionary Algorithms in Python) and Devito, python frameworks for evolutionary and automated finite difference computations, respectively. We implement a Particle swarm optimization (PSO) with an "elitism strategy" on top of DEAP, leveraging its transparent, simple and coherent environment for implementing of evolutionary algorithms (EAs). The high computational effort, due to the huge number of cost function evaluations (each one demanding a forward modeling step) required by PSO, is alleviated through the use of Devito as well as through parallelization with Dask. The combined use of these frameworks yields not only an efficient way of providing acoustic macro models of the P-wave velocity field (Vp), but also significantly reduces the amount of human effort in fulfilling this task.

Equations

$${\text{v}}_i^{k+1}=\chi \left[{\text{v}}_i^k+c_1{\text{u}}_1^k\otimes ({\text{p}}_i^k-{\text{x}}_i^k)+c_2{\text{u}}_2^k\otimes ({\text{g}}^k-{\text{x}}_i^k)\right]$$

$${\text{x}}_i^{k+1}={\text{x}}_i^k+{\text{v}}_i^{k+1},$$

$$\chi =\frac{2}{\mid 2-\phi -\sqrt{{\phi }^2-4\phi }\mid },\text{where}\phi =c_1+c_2\ge 4.$$

Full text

1 Introduction

Global optimization methods (GOMs) have presented themselves as an alternative to estimate a starting model for Full Waveform Inversion (FWI) - even when using real data (Galuzzi et al.,, 2018). GOMs are an interesting choice since a proper parameterization technique coupled with sufficient computing power allow for a reduction of the time and effort required to build an initial estimate of the velocity model. Simulated annealing (SA) and genetic algorithms (GA) have been used for a long time to solve geophysical inverse problems, they (or a variant) have already been used to estimate a starting model for FWI (Sajeva et al.,, 2016; Datta and Sen,, 2016). In recent years, however, particle swarm optimization (PSO) algorithms have rapidly become an attractive alternative for solving geophysical inverse problems and seem to enjoy an ever increasing popularity. In line with this trend, we propose to use an elitism-based PSO to produce a background model that is used as an input for FWI. The main impediment to PSO’s usage as other GOMs is the high computational cost due to the huge number of cost function evaluations involved. However, this restriction is limited to cases where the forward modeling is computationally expensive. Thanks to advances in high-performance computing, computationally faster forward modeling is now a reality in many important cases in geophysics, including seismic methods. Particularly noteworthy in the seismic modeling context is the emergence of the Devito domain-specific language (DSL) for automated stencil computations (Luporini et al.,, 2018; Louboutin et al.,, 2018). Devito provides a Python-based syntax to easily express Finite differences (FD) approximations of partial differential equations (PDEs) such as the acoustic wave equation. Although the use of GOMs to provide an initial estimate for the velocity model comes at a computational cost, we believe this cost is compensated by the significantly reduced human effort. Combining different packages in the python ecosystem allows one to use tested, performance-optimized code instead of reinventing the wheel. To this end, we use DEAP (Fortin et al.,, 2012) which is a framework written in Python that simplifies the execution of many optimization ideas with parallelization features. For the parallelization of fitness evaluations over multiple nodes, we used Dask. DEAP enabled us to integrate Dask and Devito easily in our application and to incorporate an “elitism strategy”, which was used successfully by another GOM in solving the 2D seismic optimization problem (Sajeva et al.,, 2016). By integrating DEAP, Dask, and Devito, we can nearly automate the determining of a reliable starting model for FWI, since only a few parameters have to be set before start the application. We demonstrate the effectiveness of our approach using synthetic acoustic data of the Marmousi model.

2 Theory

Elitist-mutated PSO: PSO is a population-based stochastic algorithm and is a member of the broad category of swarm intelligence techniques based on the metaphor of social interaction. The PSO algorithm is initialized with a population of random candidate solutions, conceptualized as particles. Each particle $\textbf{x}_{i}=x_{i1},x_{i2},\allowbreak\ldots,x_{iD}$ is assigned a randomized velocity $\textbf{v}_{i}=v_{i1},v_{i2},\ldots,v_{iD}$ and is iteratively moved through the $D$-dimensional problem space. It is attracted towards the location of the best fitness achieved so far by the particle itself $\textbf{p}_{i}=p_{i1},p_{i2},\ldots,p_{iD}$ and by the location of the best fitness achieved so far across the whole population $\textbf{g}_{i}=g_{1},g_{2},\ldots,g_{D}$ ($gbest$-global version of the algorithm). At iteration $k$, the basic PSO algorithm (Clerc,, 1999) can be described in vector notation as follows:

[TABLE]

[TABLE]

In Eq. 1, $\chi$, $c_{1}$, and $c_{2}$ are the control parameters called the constriction factor, cognitive parameter, and social parameter, respectively. The former is a function of $c_{1}$ and $c_{2}$ as reflected in Eq. 3.

[TABLE]

On the other hand, vectors $\textbf{u}_{1}$ and $\textbf{u}_{2}$ are $D$-dimensional vectors of uniformly distributed and independent random numbers in the $[0,1]$ range used to maintain the population diversity ($\otimes$ denotes element-by-element vector multiplication). We used an improved PSO, named EMPSO (Nagesh Kumar and Janga Reddy,, 2007), which introduces an elitist-mutation strategy into the PSO to improve its performance. Pseudocode of the EMPSO algorithm is presented in Fig. 1. In EMPSO, the elitist-mutation step is computed as follows: first, all particles are sorted in ascending order based on their fitness function and the index numbers for the respective particles are obtained; second, the elitist mutation (EM) is performed on $NM$ worst particles and the respective particle position vectors are replaced with the new mutated position vectors, whereas the velocity vectors of these particles are unchanged.

3 Algorithm overview

DEAP: DEAP (https://github.com/DEAP/deap) is a Python-based evolutionary computation framework. We chose it because it provides many useful features out of the box, and it is current, actively maintained and well documented. Moreover, DEAP is highly versatile, whereby most central members of its class hierarchy, such as individuals and operators, are fully customizable with user defined implementations. We used DEAP’s own PSO and added a function to implement elitism. At first, we create a new population. For each particle in the population, we calculate the fitness. There are three main loops. The outer loop is repeated for every generation until the pre-defined maximum number of generations is reached. At termination, several statistics and the final population is saved to a log file. The second loop iterates over all particles. This loop is distributed over multiple nodes in a cluster using Dask. Each particle’s fitness is calculated using Devito following which the particle is updated as per Eqs. 1 and 2. The innermost loop is executed for a prescribed number of particles ($NM$ in Figure 1) if the EM strategy is enabled.

Devito: Devito (https://github.com/opesci/devito) is a DSL embedded in Python and is specifically designed for finite differences in the context of seismic modeling and inversion. It offers a portable framework for the automated generation of finite-difference code from a symbolic description of PDEs. It allows the description of arbitrary time-dependent PDEs as symbolic Python expressions, from which optimized C code implementing a full time-stepping modeling loop is automatically generated, compiled and executed from the application environment. The Devito compiler introduces multiple performance optimizations when it turns the symbolic PDE representation into stencil code.

Parallelization: DEAP provides an easy way to evaluate individuals in a population on several cores in parallel. The user need merely provide an implementation of a map function. Here we used the distributed map function from Dask. Dask is an open source Python library that provides easy interfaces to scale python code across a large cluster. Although Dask handles arbitrary task graphs, here we only exploit its map functionality. (https://dask.org/).

Configuration values of the algorithm are specified in a JSON-file. The file contains forward modeling parameters such as the model size, vertical and horizontal space sampling or maximum frequency, PSO parameters and variables to be employed in the model parameterization (See section below). It is a dictionary of dictionaries. The values of cgrid dictionary are used to define the number of unknown medium parameters, whereas optimum values of pso dictionary variables will have influence on the optimization process and the output. The number of variables we have to set up is limited to the number of parameters of these two dictionaries. An example of the JSON file is shown in Figure 2. The lambda key in the pso dictionary is used to define the restriction range limits $[-\textbf{v}_{max},\textbf{v}_{max}]$ of the particles velocity according to $\textbf{v}_{max}=(\textbf{x}_{max}-\textbf{x}_{min})\times\lambda$, while the scale key from cgrid dictionary is the number of parameters scaling factor.

showstringspaces = false, keywords = false,true, alsoletter = 0123456789., morestring = [s]””, stringstyle = , MoreSelectCharTable =\lst@DefSaveDef‘:\colon@json\processColon@json, basicstyle = , keywordstyle = **, numbers=left, numberstyle=, stepnumber=1, numbersep=8pt, backgroundcolor=, **

4 Testing approach (EMPSO + local FWI)

In this section, we show the results obtained when the complete workflow (global

• local search methods) was applied for retrieving a cropped Marmousi model ($285\times 369$ samples, with vertical and horizontal space sampling of 8 m and 25 m, respectively). For the forward modeling of the EMPSO algorithm, we use the FD method (optimized FD stencil code generated with Devito), with an accuracy of 2nd order in time and 8th order in space. The acquisition geometry consisted of 50 sources and 369 receivers, one receiver for each sample on the horizontal axis. We generated synthetic data using a Ricker wavelet with a maximum frequency of 17 Hz, and set the sampling and recording times to 0.8 ms and 4 s, respectively. This dataset was then filtered (below 6 Hz), and a Ricker wavelet with a maximum frequency of 6 Hz was used to compute the modeled data. To evaluate the misfit, we use the $l_{2}$ norm. Proceeding as in Sajeva et al., (2017), we use a simple 1D $V_{p}$ model (which together with the water bottom depth constitute the prior information) with velocities linearly increasing with depth from 1500 to 3500 m/s. This model is used to center the EMPSO inversion ranges and also to build the irregular EMPSO grid by following predefined resolution criteria. For full details on such criteria, see Sajeva et al., (2017). The resulting grid (black dots) and the linear 1D model are shown in Figure 3-a. This grid has 176 nodes. These nodes are bilinearly interpolated to the FD grid for the forward-modeling following what has come to be called a “two grid strategy”. The ranges for the $V_{p}$ values during the EMPSO inversion are shown in Figure 3-b. We defined the minimum and maximum limits for the first and last level of depth as a percentage of the velocity value of the grid nodes at these levels. The limits for the intermediate levels of depth are defined by the lines passing through the maximum and minimum points of the shallower and deeper levels.

For the EMPSO we used 360 particles, 100 iterations, $NM$=25%, $P_{em}$=0.3 and the gbest topology (Kennedy,, 1999). In EMPSO, the EM step begins from 10th iteration (10% of the maximum number of iterations) and the coefficients of cognitive ($c_{1}$) and social($c_{2}$) acceleration were set to 1.2 and 2.9, respectively. In the optimization process, if the model variables violate their upper or lower bounds, they are artificially brought back into the search space, i.e. a box constraint. In the EMPSO inversion, we performed 36,000 model evaluations and the final best-fitting model is used as a starting point for a local full-waveform inversion.

While Devito could be used in the full workflow (global + local inversions), our initial work has focused on showing its potential to reduce the time required by the GOMs solving the wave equations a large number of times. This enables the use of GOMs in a reasonable amount of execution time. The implementation of a complete workflow for FWI built with Devito is a larger project that we intend to carry out in the near future. The lessons learnt from this distributed GOM implementation shall be useful in that implementation. For this work, instead of using Devito to implement a local FWI, an in-house FWI code is used to do the job. Our implemented descent-based FWI algorithm uses the steepest-descent method and a multiscale approach (performing thirty iterations for each frequency band with maximum frequencies of 4.6, 11.5, 18.4, 25.3, 32.2 and 39 Hz). The line search along the gradient search directions uses the Barzilai-Borwein (BB) formula for an initial step length (Barzilai and Borwein,, 1988). When required, it applies a backtracking line search method to update the step length. The forward problem in FWI is formulated in the time domain and solved using an FD method having an accuracy of 2nd order in time and 16th order in space (Devito was not used for this), with a time step of 0.8 ms to ensure stability. The recording time and sampling grid (dx and dy) were equal to those used by EMPSO.

The experiments were run on the ÒGÚN Supercomputer at SENAI CIMATEC, which uses an Ethernet interconnection. Each compute node used contains 192 GB of RAM and two sockets, where each socket has an Intel Xeon Gold 6148 CPU at 2.4 GHz.

5 Results

Figure 4 ilustrates EMPSO results for three random trials. Figure 5 shows the comparison of the convergence rate between standard PSO and EMPSO. For both algorithms, the error gradually decreases over time. It should be noted, however, that EMPSO gives better fitness values over different trials than standard PSO. EMPSO does not seem to experience long periods of stagnation as PSO (apparent from the staircase pattern in fitness curves). Figure 6 shows the final models after descent-based FWI using as starting models the velocity estimates retrieved by EMPSO. The correct Marmousi model is shown repeatedly in Figures 4-a and 6-a for ease of comparison.

In the test scenario outlined in the last section, Devito took less than one second to generate a single shot gather with 2 threads. It means that the fitness function evaluation of a particle (50 shots) took approximately one minute. In the best of cases, with an equal number of available parallel workers and tasks (fitness evaluations), one iteration would be completed in this time, which would lead to an extremely fast overall processing time.

In general, the velocity models obtained by FWI using EMPSO outputs as starting points are very similiar to the true velocity model of Figures 4-a and 6-a. Putting it another way, the EMPSO allows the recovery of the low wavenumber components in the background model to avoid the cycle-skipping problem, which leads to good final results at the end of the entire workflow.

6 Conclusions

In order to determine a promising starting model for FWI, we have incorporated a DSL, Devito, designed optimally for solving wave equations into DEAP an Evolutionary Computation Framework, through Dask, a parallelization framework. With the help of DEAP, we have devised an elitis-based PSO that finds solutions to the seismic optimization problem. The fitness evaluation (the most computationally demanding step) is parallelized over a cluster using Dask and uses Devito to generate highly optimized C code for solving wave equations. The Devito, Dask, and DEAP mixture offers an optimised and almost automatic procedure to estimate a starting velocity model for FWI. We test our approach on a 2D acoustic FWI benchmark problem, namely the Marmousi model. The similarity between the true and the final models obtained by gradient-based inversions that started from EMPSO models, makes the proposed approach well suited for the hard task of finding a good first-guess model for FWI.

7 Acknowledgements

Computational resources and services used in this work were provided by the High Performance Computing (HPC) and Research Support Group of SENAI CIMATEC, Salvador, Brazil. The authors also gratefully acknowledge support from Shell Brasil through the PSO-FWI project at SENAI CIMATEC and the strategic importance of the support given by ANP through the R&D levy regulation. This work was supported in part by Intel Parallel Computing Centre at Imperial College London and EPSRC EP/R029423/1.

\append

The source of the bibliography