Interior-point methods for PDE-constrained optimization

In applied sciences PDEs model an extensive variety of phenomena. Typically the final goal of simulations is a system which is optimal in a certain sense. For instance optimal control problems identify a control to steer a system towards a desired state. Inverse problems seek PDE parameters which are most consistent with measurements. In these optimization problems PDEs appear as equality constraints. PDE-constrained optimization problems are large-scale and often nonconvex. Their numerical solution leads to large ill-conditioned linear systems. In many practical problems inequality constraints implement technical limitations or prior knowledge. In this thesis interior-point (IP) methods are considered to solve nonconvex large-scale PDE-constrained optimization problems with inequality constraints. To cope with enormous fill-in of direct linear solvers, inexact search directions are allowed in an inexact interior-point (IIP) method. This thesis builds upon the IIP method proposed in [Curtis, Schenk, Wächter, SIAM Journal on Scientific Computing, 2010]. SMART tests cope with the lack of inertia information to control Hessian modification and also specify termination tests for the iterative linear solver. The original IIP method needs to solve two sparse large-scale linear systems in each optimization step. This is improved to only a single linear system solution in most optimization steps. Within this improved IIP framework, two iterative linear solvers are evaluated: A general purpose algebraic multilevel incomplete L D L^T preconditioned SQMR method is applied to PDE-constrained optimization problems for optimal server room cooling in three space dimensions and to compute an ambient temperature for optimal cooling. The results show robustness and efficiency of the IIP method when compared with the exact IP method. These advantages are even more evident for a reduced-space preconditioned (RSP) GMRES solver which takes advantage of the linear system's structure. This RSP-IIP method is studied on the basis of distributed and boundary control problems originating from superconductivity and from two-dimensional and three-dimensional parameter estimation problems in groundwater modeling. The numerical results exhibit the improved efficiency especially for multiple PDE constraints. An inverse medium problem for the Helmholtz equation with pointwise box constraints is solved by IP methods. The ill-posedness of the problem is explored numerically and different regularization strategies are compared. The impact of box constraints and the importance of Hessian modification on the optimization algorithm is demonstrated. A real world seismic imaging problem is solved successfully by the RSP-IIP method

High-resolution imaging beneath the Santorini volcano

Volcanoes are surface expressions of much deeper magmatic systems, inaccessible to direct observation. Constraining the geometry and physical properties of these systems, in particular detecting high melt fraction (magma) reservoirs, is key to managing a volcanic hazard and understanding fundamental processes that lead to the formation of continents. Unfortunately, unambiguous evidence of magma reservoirs has not yet been provided due to the limited resolving power of the geophysical methods used so far. Here, a high-resolution imaging technique called full-waveform inversion was applied to study the magmatic system beneath the Santorini volcanic field, one of the most volcanically and seismically active regions of Europe. Quality-controlled inversion of 3d wide-angle, multi-azimuth ocean-bottom seismic data revealed a previously undetected high melt fraction reservoir 3 km beneath the Kolumbo volcano, a centre of microseismic and hydrothermal activity of the field. To enable the above method to handle land data, two major algorithmic improvements were added to the high-performance inversion code. First, to simulate instrument response of land seismometers, a pressure-velocity conversion has been implemented in a way that ensures reciprocity of the discretised 2nd-order acoustic wave equation. Second, the immersed-boundary method, originally developed for computational fluid dynamics, was implemented to simulate the wave-scattering off the irregular topography of the Santorini caldera. These advancements can be readily used to provide a higher-resolution image of the melt reservoir beneath the Santorini caldera already detected by means of travel-time tomography.Open Acces

Full-Waveform Inversion To 3D Seismic Land Data

Full-waveform inversion (FWI) is a technique that seeks to find a high-resolution high-fidelity model of the Earth's subsurface that is capable of matching individual seismic waveforms, within an original raw field dataset, trace by trace. The method begins from a best-guess starting model, which is then iteratively improved using a sequence of linearized local inversions to solve a fully non-linear problem. In principle, FWI can be used to recover any physical property that has an influence upon the seismic wavefield, but in practice the technique has been used predominantly to recover P-wave velocity, and this is the route that is followed here. Full-waveform tomographic techniques seek to determine a highly resolved quantitative model of the sub-surface that will ultimately be able to explain the entire seismic wavefield including those phases that conventional processing and migration seek to remove such as refracted arrivals. Although the underlying theory of FWI is well established, its practical application to 3D land data, and especially to seismic data that have been acquired using vibrators, in a form that is effective and robust, is still a subject of intense research. In this study, 2D and 3D FWI techniques have been applied to a vibrator dataset from onshore Oman. Both the raw dataset and the subsurface model cause difficulties for FWI. In particular, the data are noisy, have weak early arrivals, are strongly elastic, and especially are lacking in low-frequency content. The Earth model appears to contain shallow low-velocity layers, and these compromise the use of first-arrival travel-time tomography for the generation of a starting velocity model. The 2D results show good recovery of the shallow part of the velocity models. The results show a low-velocity layer that extends across the velocity model, but lacking in a high-resolution image due to the absence of the third dimension. The seismograms of the final inversion models give a good comparison with the field data and produce a reasonably high correlation coefficient compared to the starting model. An inversion scheme has been developed in this study in which only data from the shorter offsets are initially inverted since these represent the subset of the data that is not cycle skipped. The offset range is then gradually extended as the model improves. The final 3D model contains a strongly developed low-velocity layer in the shallow section. The results from this inversion appear to match p-wave logs from a shallow drill hole, better flatten the gathers, and better stack and migrate the reflection data. The inversion scheme is generic, and should have applications to other similar difficult datasets

Seismic full-waveform inversion of the crust-mantle structure beneath China and adjacent regions

Since the late 1970s, seismic tomography has been emerging as the pre-eminent tool for imaging the Earth’s interior from the meter to the global scale. Significant recent advances in seismic data acquisition, high-performance computing, and modern numerical methods have drastically progressed tomographic methods. Today it is technically feasible to accurately simulate seismic wave propagation through realistically heterogeneous Earth models across a range of scales. When seismic waves propagate inside the Earth and encounter structural heterogeneities with a certain scale, wave propagation speed changes, reflection, and scattering phenomena occur, and interconversions between compressional and shear waves happen. The combined effect of multiple heterogeneities produces a highly complicated wavefield recorded in the form of three-component (vertical, radial, transverse) seismograms. The ultimate objective of seismic imaging is to utilize the full information from waveforms recorded at seismic stations distributed around the globe in a broad frequency range to characterize detailed tomographic images of Earth’s interior by fitting synthetic seismograms to recorded seismograms. The full-waveform inversion technique based on adjoint and spectral-element methods can be employed to maximumly exploit the information contained in these seismic wavefield complexities to determine the fine-scale structural heterogeneities from which they originated across various orders of magnitude in frequency and wavelength. Over the past two decades, the three-dimensional crust and mantle structure beneath the broad Asian region has attracted much attention in seismic studies due to its complicated and unique geological setting involving active lithospheric deformation, intracontinental rifting, intraplate seismotectonics, volcanism and magmatism, continent-continent collision, oceanic plate deep subduction, and mantle dynamics. To enhance our understanding of the subsurface behavior of cold subducting slabs and hot mantle flows and their dynamic relation to the tectonic evolution of the overriding plates as mentioned above, we construct the first-generation full-waveform tomographic model SinoScope 1.0 for the crust-mantle structure beneath China and adjacent regions. The three-component seismograms from 410 earthquakes recorded at 2,427 seismographic stations are employed in iterative gradient-based inversions for three successively broadened period bands of 70 - 120 s, 50 - 120 s, and 30 - 120 s. Synthetic seismograms were computed using GPU-accelerated spectral-element simulations of seismic wave propagation in 3-D anelastic models, and Fréchet derivatives were calculated based on an adjoint-state method facilitated by a checkpointing algorithm. The inversion involved 352 iterations, which required 18,600 wavefield simulations. The simulations were performed on the Xeon Platinum ‘SuperMUC-NG’ at the Leibniz Supercomputing Center and the Xeon E5 ‘Piz Daint’ at the Swiss National Supercomputing Center with a total cost of ~8 million CPU hours. SinoScope 1.0 is described in terms of isotropic P-wave velocity (Vp), horizontally and vertically polarized S-wave velocities (Vsh and Vsv), and mass density (ρ) variations, which are independently constrained with the same data set coupled with a stochastic L-BFGS quasi-Newton optimization scheme. It systematically reduced differences between observed and synthetic full-length seismograms. We performed a detailed resolution analysis by repairing input random parametric perturbations, indicating that resolution lengths can approach the half-propagated wavelength within the well-covered areas. SinoScope 1.0 exhibits strong lateral heterogeneities in the crust and uppermost mantle, and features correlate well with geological observations, such as sedimentary basins, Holocene volcanoes, Tibetan Plateau, Philippine Sea Plate, orogenic belts, and subduction zones. Estimating lithospheric thickness from seismic velocity reductions at depth reveals significant variations underneath the different tectonic units:~50 km in Northeast and North China, ~70 km in South China, ~90 km in the South China Sea, Philippine Sea Plate, and Caroline Plate, and 150-250 km in the Indian Plate. The thickest (200-250 km) cratonic roots are present beneath the Sichuan, Ordos, and Tarim basins. The large-scale lithospheric deformation is imaged as low-velocity zones from the Himalayan orogen to the Baikal rift system in central Asia. A thin horizontal layer of ~100-200 km depth extent below the lithosphere points toward the existence of the asthenosphere beneath East and Southeast Asia, with heterogeneous anisotropy indicative of channel flows. This provides independent, high-resolution evidence for the low-viscosity asthenosphere that partially decouples plates from mantle flow beneath and allows plate tectonics to work above. Beneath the Indo-Australian Plate, we observe distinct low-velocity anomalies from a depth of ~200 km to the bottom of the mantle transition zone (MTZ), continuously extending northward below western China from the lower MTZ down to the top of the lower mantle. Furthermore, we observe an enhanced image of well-known slabs along strongly curved subduction zones, including Kurile, Japan, Izu-Bonin, Mariana, Ryukyu, Philippines, Indonesia, and Burma. Broad high-velocity bodies persist from the lower MTZ to 1000 km depth or deeper beneath the north of the Indo-Australian Plate. They might be pieces of the ancient Tethyan slab sinking down to the lower mantle before the Indo-Australian Plate and Eurasian Plate collision. The deep geodynamic processes controlling the large-scale tectonic activity of the broad Asian region are very complicated and not yet well understood, which is a source of much debate. In this thesis, the main focus is on deciphering the three-dimensional seismic structure and dynamics of the lithosphere and mantle beneath China and adjacent regions with the help of the high-resolution full-waveform tomographic model. More importantly, in the subsequent works of geodynamic inversion, it provides the improved seismological basis for a sequential reconstruction of the late Mesozoic and Cenozoic plate motion history of the broad Asian region and the present-day mantle heterogeneity state estimates that, in turn, allow one to track material back in time from any given sampling location through retrodicting past mantle states

Seismic Waves

The importance of seismic wave research lies not only in our ability to understand and predict earthquakes and tsunamis, it also reveals information on the Earth's composition and features in much the same way as it led to the discovery of Mohorovicic's discontinuity. As our theoretical understanding of the physics behind seismic waves has grown, physical and numerical modeling have greatly advanced and now augment applied seismology for better prediction and engineering practices. This has led to some novel applications such as using artificially-induced shocks for exploration of the Earth's subsurface and seismic stimulation for increasing the productivity of oil wells. This book demonstrates the latest techniques and advances in seismic wave analysis from theoretical approach, data acquisition and interpretation, to analyses and numerical simulations, as well as research applications. A review process was conducted in cooperation with sincere support by Drs. Hiroshi Takenaka, Yoshio Murai, Jun Matsushima, and Genti Toyokuni

Inner product preconditioned trust-region methods for frequency-domain full waveform inversion

Full waveform inversion is a seismic imaging method which requires to solve a large-scale minimization problem, typically through local optimization techniques. Most local optimization methods can basically be built up from two choices: the update direction and the strategy to control its length. In the context of full waveform inversion, this strategy is very often a line search. We here propose to use instead a trust-region method, in combination with non-standard inner products which act as preconditioners. More specifically, a line search and several trust-region variants of the steepest descent, the limited memory BFGS algorithm and the inexact Newton method are presented and compared. A strong emphasis is given to the inner product choice. For example, its link with preconditioning the update direction and its implication in the trust-region constraint are highlighted. A first numerical test is performed on a 2D synthetic model then a second configuration, containing two close reflectors, is studied. The latter configuration is known to be challenging because of multiple reflections. Based on these two case studies, the importance of an appropriate inner product choice is highlighted and the best trust-region method is selected and compared to the line search method. In particular we were able to demonstrate that using an appropriate inner product greatly improves the convergence of all the presented methods and that inexact Newton methods should be combined with trust-region methods to increase their convergence speed

Construção de um modelo inicial em profundidade usando métodos robustos de análise de velocidade por migração em tempo

Orientadores: Joerg Dietrich Wilhelm Schleicher, Maria Amélia Novais SchleicherTese (doutorado) - Universidade Estadual de Campinas, Faculdade de Engenharia Mecânica e Instituto de GeociênciasResumo: A necessidade de investigar regiões formadas por estruturas geológicas complexas tem motivado o desenvolvimento de métodos de imageamento que atuem no domínio da profundidade. Exemplos notáveis são as técnicas de migração pré-empilhamento em profundidade (PSDM, do inglês "prestack depth migration") e a tomografia de onda completa (FWT, do inglês "full-waveform tomography"). No entanto, a aplicação desses métodos enfrenta ao menos dois desafios: eles requerem (1) um modelo de velocidade (inicial) preciso, e (2) elevado poder computacional. Por outro lado, a migração em tempo provou ser um processo robusto e muito rápido, tornando-se rotineiramente empregado para o imageamento sísmico. Além disso, a construção de modelos de velocidade em tempo é um processo bem compreendido. Portanto, é altamente desejável usar as técnicas de conversão tempo-profundidade para construir, a partir desses modelos de velocidade no domi?nio do tempo, modelos de velocidade iniciais para te?cnicas que operam em profundidade. Neste trabalho, investigamos a aplicabilidade de um fluxo de trabalho formado por alguns recém-desenvolvidos métodos (semi-) automáticos de análise de velocidade de migração em tempo (MVA, do inglês "migration velocity analysis"), capazes de gerar modelos de velocidade e imagens migradas no tempo sem precisar de informações a priori, seguido por uma técnica robusta de conversão tempo-profundidade. Discutimos as vantagens e limitações desse fluxo de trabalho e suas perspectivas para se tornar uma ferramenta plenamente automática, capaz de gerar modelos de velocidade sísmica para o uso subsequente em métodos de FWT. Nos nossos testes em diferentes versões dos dados Marmousi, o procedimento proposto produziu modelos de velocidade iniciais suficientemente precisos para uma FWT sob condições quase ideais. Começando no modelo de velocidade do domínio do tempo convertido para profundidade, a FWT convergiu para um modelo final com qualidade comparável a quando feito a partir de uma versão suavizada do modelo de velocidade verdadeiro. Isso indica que a correta informação sobre a velocidade de fundo pode ser extraída com sucesso pela MVA automática no domínio do tempo mesmo em meios onde a migração em tempo não pode fornecer imagens sísmicas satisfatórias. Como resultado, esta tese não só contribui para o desenvolvimento de um fluxo de trabalho para a construção de modelos de velocidade iniciais para a FWT, mas também apresenta várias aplicações inovadorasAbstract: The need to investigate regions with complex geology has encouraged the development of imaging methods that act in the depth domain. Notable examples are prestack depth migration (PSDM) and full-waveform tomography (FWT). However, the application of these techniques faces at least two challenges: they require (1) an accurate (initial) velocity model and (2) massive computation power. In contrast, time migration has proven to be a fast and robust process, making it routinely used for seismic imaging. Moreover, time-domain velocity-model building is a well-understood process. Therefore, it is highly desirable to use time-to-depth conversion to construct starting models for depth-imaging techniques from these time-domain velocity models. In this work, we investigate the applicability of a workflow consisting of some recent (semi-) automatic time migration-velocity-analysis (MVA) methods, which can generate a velocity model and a time-migrated image without a priori information, followed by a robust time-to-depth conversion technique. We discuss advantages and limitations of this workflow and its perspectives to become a fully automatic tool, capable of generating initial seismic depth velocity models for subsequent FWT methods. In our tests on different versions of the Marmousi data, the proposed procedure produced sufficiently accurate initial models for an FWT under nearly ideal conditions. Starting at the depth-converted time-domain model, FWT converged to a final model of comparable quality as when starting at a smoothed version of the true velocity model. This indicates that correct background velocity information can be successfully extracted from automatic time-domain MVA even in media where time-migration cannot provide satisfactory seismic images. In effect, this thesis not only contributes to the development of a workflow for the construction of initial velocity-models for FWT but also presents several innovative applications thereofDoutoradoReservatórios e GestãoDoutor em Ciências e Engenharia de Petróle

Genetic full waveform inversion to characterise fractures

Active seismic methodologies provide a non-invasive tool to remotely characterise the physical properties of fractures at a wide range of scales, and have a positive impact in helping to solve rock engineering problems in a variety of geo-industrial applications. With current advances in seismic processing tools, such as full waveform inversion (FWI), and accurate models of seismic wave interaction with fractures, seismic characterisation of fractures can be tackled by utilising the entire seismic wavefield recorded at the receiver locations. A two-step strategy, using the genetic algorithm (GA) for global optimisation and the Neighbourhood Algorithm (NA) for evaluating uncertainties, was developed to simultaneously estimate the fracture properties (both fracture specific stiffness and equivalent fracture stiffness) and the background material properties directly from seismic waveforms. The optimisation involves minimising the difference between the observed (measured) and forward-modelled full waveforms through the finite difference code WAVE3D. The development, named Genetic Algorithm Full-Waveform Fracture Inversion (GAFWFI), looks beyond conventional seismic methods which focus on characterising fracture-induced anisotropy, by reducing the need to manually condition the data (e.g. manual picking of seismic phases), and by providing a robust means to explore multiple solutions. The development also allows the gap between different representations of fracturing to be bridged within a comprehensive method which can employ both discrete fracture and effective fracture models. GA-FWFI is tested initially on synthetic ultrasonic experiments with parallel fractures. Results confirm that the method can effectively invert for physical properties such as fracture stiffness, location, background material properties, while the posterior probability density (PPD) show that inversions are very well constrained. GA-FWFI is then applied to waveforms from a laboratory experiment investigating fracture slip and again results show high degree of accuracy. GA-FWFI is then utilised to unveil the coupling between discrete fracture networks (DFNs) and their equivalent fracture zone properties. The results reveal that the transition from a medium with open cracks to one with welded interfaces leads to the equivalent media having the equivalent medium stiffness non-linearly related to the crack specific stiffness. An attribute χ is proposed which helps guide the interpretation of a cracked medium by giving a range of likely values for crack size and crack stiffness. This work paves the way for novel strategies to seismically characterise fractures

Marchenko-Lippmann-Schwinger inversion

Seismic wave reflections recorded at the Earth’s surface provide a rich source of information about the structure of the subsurface. These reflections occur due to changes in the material properties of the Earth; in the acoustic approximation, these are the density of the Earth and the velocity of seismic waves travelling through it. Therefore, there is a physical relationship between the material properties of the Earth and the reflected seismic waves that we observe at the surface. This relationship is non-linear, due to the highly scattering nature of the Earth, and to our inability to accurately reproduce these scattered waves with the low resolution velocity models that are usually available to us. Typically, we linearize the scattering problem by assuming that the waves are singly-scattered, requiring multiple reflections to be removed from recorded data at great effort and with varying degrees of success. This assumption is called the Born approximation. The equation that describes the relationship between the Earth’s properties and the fully-scattering reflection data is called the Lippmann-Schwinger equation, and this equation is linear if the full scattering wavefield inside the Earth could be known. The development of Marchenko methods makes such wavefields possible to estimate using only the surface reflection data and an estimate of the direct wave from the surface to each point in the Earth. Substituting the results from a Marchenko method into the Lippmann-Schwinger equation results in a linear equation that includes all orders of scattering. The aim of this thesis is to determine whether higher orders of scattering improve the linear inverse problem from data to velocities, by comparing linearized inversion under the Born approximation to the inversion of the linear Lippmann-Schwinger equation. This thesis begins by deriving the linear Lippmann-Schwinger and Born inverse problems, and reviewing the theoretical basis for Marchenko methods. By deriving the derivative of the full scattering Green’s function with respect to the model parameters of the Earth, the gradient direction for a new type of least-squares full waveform inversion called Marchenko-Lippmann-Schwinger full waveform inversion is defined that uses all orders of scattering. By recreating the analytical 1D Born inversion of a boxcar perturbation by Beydoun and Tarantola (1988), it is shown that high frequency-sampling density is required to correctly estimate the amplitude of the velocity perturbation. More importantly, even when the scattered wavefield is defined to be singly-scattering and the velocity model perturbation can be found without matrix inversion, Born inversion cannot reproduce the true velocity structure exactly. When the results of analytical inversion are compared to inversions where the inverse matrices have been explicitly calculated, the analytical inversion is found to be superior. All three matrix inversion methods are found to be extremely ill-posed. With regularisation, it is possible to accurately determine the edges of the perturbation, but not the amplitude. Moving from a boxcar perturbation with a homogeneous starting velocity to a many-layered 1D model and a smooth representation of this model as the starting point, it is found that the inversion solution is highly dependent on the starting model. By optimising an iterative inversion in both the model and data domains, it is found that optimising the velocity model misfit does not guarantee improvement in the resulting data misfit, and vice versa. Comparing unregularised inversion to inversions with Tikhonov damping or smoothing applied to the kernel matrix, it is found that strong Tikhonov damping results in the most accurate velocity models. From the consistent under-performance of Lippmann-Schwinger inversion when using Marchenko-derived Green’s functions compared to inversions carried out with true Green’s functions, it is concluded that the fallibility of Marchenko methods results in inferior inversion results. Born and Lippmann-Schwinger inversion are tested on a 2D syncline model. Due to computational limitations, using all sources and receivers in the inversion required limiting the number of frequencies to 5. Without regularisation, the model update is uninterpretable due to the presence of strong oscillations across the model. With strong Tikhonov damping, the model updates obtained are poorly scaled, have low resolution, and low amplitude oscillatory noise remains. By replacing the inversion of all sources simultaneously with single source inversions, it is possible to reinstate all frequencies within our limited computational resources. These single source model updates can be stacked similarly to migration images to improve the overall model update. As predicted by the 1D analytical inversion, restoring the full frequency bandwidth eliminates the oscillatory noise from the inverse solution. With or without regularisation, Born and Lippmann-Schwinger inversion results are found to be nearly identical. When Marchenko-derived Green’s functions are introduced, the inversion results are worse than either the Born inversion or the Lippmann-Schwinger inversion without Marchenko methods. On this basis, one concludes that the inclusion of higher order scattering does not improve the outcome of solving the linear inverse scattering problem using currently available methods. Nevertheless, some recent developments in the methods used to solve the Marchenko equation hold some promise for improving solutions in future