Combining checkpointing and data compression for large scale seismic inversion

Seismic inversion and imaging are adjoint-based optimization problems that processes up to terabytes of data, regularly exceeding the memory capacity of available computers. Data compression is an effective strategy to reduce this memory requirement by a certain factor, particularly if some loss in accuracy is acceptable. A popular alternative is checkpointing, where data is stored at selected points in time, and values at other times are recomputed as needed from the last stored state. This allows arbitrarily large adjoint computations with limited memory, at the cost of additional recomputations. In this paper we combine compression and checkpointing for the first time to compute a realistic seismic inversion. The combination of checkpointing and compression allows larger adjoint computations compared to using only compression, and reduces the recomputation overhead significantly compared to using only checkpointing

Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion

We introduce a technique to compute exact anelastic sensitivity kernels in the time domain using parsimonious disk storage. The method is based on a reordering of the time loop of time-domain forward/adjoint wave propagation solvers combined with the use of a memory buffer. It avoids instabilities that occur when time-reversing dissipative wave propagation simulations. The total number of required time steps is unchanged compared to usual acoustic or elastic approaches. The cost is reduced by a factor of 4/3 compared to the case in which anelasticity is partially accounted for by accommodating the effects of physical dispersion. We validate our technique by performing a test in which we compare the $K_\alpha$ sensitivity kernel to the exact kernel obtained by saving the entire forward calculation. This benchmark confirms that our approach is also exact. We illustrate the importance of including full attenuation in the calculation of sensitivity kernels by showing significant differences with physical-dispersion-only kernels

Imagerie de milieux salifères aux échelles crustales et expérimentales par méthodes de migration sismique et méthode de l'adjoint : applications marines

La découverte de structures géologiques salines est une raison économique importante pour l'exploration dans le monde car elles constituent un piège naturel pour diverses ressources. Cependant, l'imagerie de ces structures est un grand défi. En raison des propriétés du sel, dont les vitesses de propagation sont beaucoup plus élevées que celles des strates adjacentes, les ondes sismiques sont piégées dans ces structures, produisant un grand nombre d'artefacts numériques parasites, tels que des multiples. Cela interfère avec le signal sismique primaire, ce qui empêche de voir clairement ce qui se trouve sous les structures salifères. Parmi toutes les méthodes d'exploration géophysique, la méthode de migration par renversement temporel (RTM), qui fait partie des méthodes qui utilisent la résolution de la forme d'onde sismique complète, est un outil d'imagerie très puissant, même dans les régions à géologie complexe. Dans ce travail, nous utilisons la méthode RTM basée sur l'adjoint, qui consiste essentiellement en trois étapes : la solution de l'équation des ondes, la solution de l'équation des ondes adjointe et la condition d'imagerie, qui consiste en la corrélation des champs d'ondes directs et adjoints. Ce travail peut être divisé en deux cas d'étude : le premier cas consiste en un modèle synthétique bidimensionnel d'un dôme de sel, issu de la migration finale d'une étude réelle dans le Golfe du Mexique. Le second cas consiste en un modèle tridimensionnel expérimental (WAVES), élaboré par le laboratoire LMA de Marseille, qui simule une structure saline, structures sédimentaires environnantes, et un socle. Le modèle a été immergé dans l'eau pour recréer un sondage marin réaliste. Deux types de données différents ont été obtenus dans cette expérience : des données à décalage nul et des données à décalage multiple. Pour résoudre l'équation des ondes impliquée dans la méthode RTM basée sur l'adjoint, nous utilisons des différences finies d'ordre 4 dans les deux cas. De plus, dans le second cas, nous avons utilisé le code UniSolver, qui résout la méthode RTM basée sur l'adjoint en utilisant des différences finies d'ordre 4 et un parallélisme basé sur MPI. Nous avons mis en œuvre les équations viscoélastiques pour simuler l'effet de l'atténuation. Pour cette raison, le schéma "Checkpointing" est introduit pour calculer la condition d'imagerie et assurer la stabilité physique et numérique. Dans le premier cas d'étude, nous analysons la reconstruction de l'image du dôme de sel que produisent différents noyaux de sensibilité. Nous calculons ces noyaux en utilisant différentes paramétrisations (densité - vitesse P), ou (densité - constantes de Lamé) pour une rhéologie acoustique. Nous étudions également comment l'utilisation de différents modèles a priori affecte l'image finale en fonction du type de noyau calculé. En utilisant les résultats obtenus en 2D, nous calculons des noyaux synthétiques tridimensionnels en utilisant une rhéologie élastique. Dans le second cas, nous effectuons d'abord une calibration des propriétés du modèle pour des données à décalage nul, et une fois que les données synthétiques et réelles s'ajustent bien, nous calculons les noyaux tridimensionnels. Nous résolvons le problème direct pour le cas à décalage multiple avec et sans effets d'atténuation. Nous comparons les données synthétiques calculées avec des rhéologies élastiques ou viscoélastiques avec les données réelles. Cela permet ainsi de voir l'impact de l'atténuation dans les signaux. Cela ouvrira la voie à de la RTM et des simulations de la forme d'onde complète viscoélastique dans des contextes tectoniques salifères. Enfin, nous avons implémenté l'algorithme LSRTM pour les données acoustiques synthétiques bidimensionnelles et pour les données viscoélastiques réelles tridimensionnelles, qui est un processus d'inversion itératif, nous avons suivi l'approche du gradient conjugué, en vérifiant que les conditions de Wolfe sont satisfaites.Finding salt geological structures is an important economic reason for exploration in the world because they constitute a natural trap for various resources such as oil, natural gas, water, and also the salt itself can be exploitable. However, the imaging of these structures is a great challenge. Due to the properties of salt, with propagation velocities much higher than the adjacent strata, seismic waves are trapped within these structures, producing a large number of spurious numerical artifacts, such as multiples. This interferes with the primary seismic signal, making it impossible to see clearly what is underneath the salt structures (salt domes for instance). Among all the geophysical exploration methods, the Reverse Time Migration method (RTM), which is part of the methods that solve the complete seismic waveform, is a very powerful imaging tool, even in regions of complex geology. In this work we use the adjoint-based RTM method, which basically consists of three stages: the solution of the wave equation (forward problem), the solution of the adjoint wave equation (adjoint problem), and the imaging condition, which consists in the correlation of the forward and adjoint wavefields. This work can be divided in two cases of study: the first case consists in a two-dimensional synthetic model of a salt dome, taken from the final migration of a real survey in the Gulf of Mexico. The second case consists in an experimental three-dimensional model (WAVES), elaborated by the LMA laboratory in Marseille (France), which simulates a salt structure (with surrounding sedimentary structures), and a basement. The model was immersed in water to recreate a reallistic marine survey. Two different data types were obtained in this experiment: zero-offset and multi-offset data. To compute the adjoint-based RTM method we use fourth-order finite differences in both cases. Furthermore, in the second case we used the UniSolver code, which solves the adjoint-based RTM method using fourth-order finite differences and MPI-based parallelism. It was also necessary to implement the viscoelastic equations to simulate the effect of attenuation. Because of this, the Checkpointing scheme is introduced to calculate the imaging condition and ensures physical and numerical stability in the migration procedure. In the first case study we analyze the recovery of the salt dome image that different sensitivity kernels produce. We calculate these kernels using different parametrizations (density - P velocity), (density - Lamé constants), or (density - P impedance) for an acoustic rheology. We also study how the use of different a priori models affects the final image depending on the kind of kernel computed. Using the results obtained previously in 2D, we calculate synthetic three-dimensional kernels using an elastic rheology. In the second case (the realistic/experimental case), we perform a calibration of the model properties for zero-offset data, and once the synthetic and real data fit well, we calculate the three-dimensional kernels. We compute the forward problem of the multi-offset data with and without attenuation effects. We compare the synthetic data computed with elastic or viscoelastic rheologies with the real data. This will allow to see the impact of attenuation in the signals. This will pave the way to viscoelastic full waveform in salt tectonic context. Finally, we implemented the Least Squares Reverse Time Migration (LSRTM) algorithm for the two-dimensional acoustic synthetic data and for the 3-dimensional viscoelastic real data, which is an iterative inversion process, we have followed the Conjugated Gradient (CG) approach, checking that Wolfe conditions (convergence and curvature of the misfit functions are satisfied

Algorithmic analysis torwards time-domain extended source waveform inversion

Full waveform inversion (FWI) updates the subsurface model from an initial model by comparing observed and synthetic seismograms. Due to high nonlinearity, FWI is easy to be trapped into local minima. Extended domain FWI, including wavefield reconstruction inversion (WRI) and extended source waveform inversion (ESI) are attractive options to mitigate this issue. This paper makes an in-depth analysis for FWI in the extended domain, identifying key challenges and searching for potential remedies torwards practical applications. WRI and ESI are formulated within the same mathematical framework using Lagrangian-based adjoint-state method with a special focus on time-domain formulation using extended sources, while putting connections between classical FWI, WRI and ESI: both WRI and ESI can be viewed as weighted versions of classic FWI. Due to symmetric positive definite Hessian, the conjugate gradient is explored to efficiently solve the normal equation in a matrix free manner, while both time and frequency domain wave equation solvers are feasible. This study finds that the most significant challenge comes from the huge storage demand to store time-domain wavefields through iterations. To resolve this challenge, two possible workaround strategies can be considered, i.e., by extracting sparse frequencial wavefields or by considering time-domain data instead of wavefields for reducing such challenge. We suggest that these options should be explored more intensively for tractable workflows

Modeling for inversion in exploration geophysics

Seismic inversion, and more generally geophysical exploration, aims at better understanding the earth's subsurface, which is one of today's most important challenges. Firstly, it contains natural resources that are critical to our technologies such as water, minerals and oil and gas. Secondly, monitoring the subsurface in the context of CO2 sequestration, earthquake detection and global seismology are of major interests with regard to safety and the environment hazards. However, the technologies to monitor the subsurface or find resources are scientifically extremely challenging. Seismic inversion can be formulated as a mathematical optimization problem that minimizes the difference between field recorded data and numerically modeled synthetic data. The process of solving this optimization problem then requires to numerically model, thousands of times, wave-propagation in large three-dimensional representations of part of the earth subsurface. The mathematical and computational complexity of this problem, therefore, calls for software design that abstracts these requirements and facilitates algorithm and software development. My thesis addresses some of the challenges that arise from these problems; mainly the computational cost and access to the right software for research and development. In the first part, I will discuss a performance metric that improves the current runtime-only benchmarks in exploration geophysics. This metric, the roofline model, first provides insight at the hardware level of the performance of a given implementation relative to the maximum achievable performance. Second, this study demonstrates that the choice of numerical discretization has a major impact on the achievable performance depending on the hardware at hand and shows that a flexible framework with respect to the discretization parameters is necessary. In the second part, I will introduce and describe Devito, a symbolic finite-difference DSL that provides a high-level interface to the definition of partial differential equations (PDE) such as the wave equation. Devito, from the symbolic definition of PDEs, then generates and compiles highly optimized C code on-the-fly to compute the solution of the PDE. The combination of the high-level abstractions and the just-in-time compiler enable research for geophysical exploration and PDE-constrainted optimization based on the paradigm of separation of concerns. This allows researchers to concentrate on their respective field of study while having access to computationally performant solvers with a flexible and easy to use interface to successfully implement complex representations of the physics. The second part of my thesis will be split into two sub-parts; first describing the symbolic application programming interface (API), before describing and benchmarking the just-in-time compiler. I will end my thesis with concluding remarks, the latest developments and a brief description of projects that were enabled by Devito.Ph.D