Wave Sensing, Undersampling

Wavefield retrieval, manifold learning and data reduction

Wavefield retrieval from surface array data

There have been extensive efforts in developing data regularization strategies and interpolation techniques accommodating the limitations of present day acquisition. Acquisition designs away from regular sampling and exploiting frames with respect to which data compress (motivated by compressed sensing) have been considered; however, the underlying theory is subtle and incomplete in the case of solutions to the system of equations describing elastic waves. The multi-scale physics of waves provides us with opportunities to discover new directions in the fields of probing, retrieval and compression.

  • Sparsity, compression, geometry. The wavepackets, like curvelets, for given scale and orientation, are positioned via translation on a lattice. Naturally the points in the lattice have no knowledge of the waves in the data. We developed a nonlinear technique where the points become a function of the data to be decomposed, thus departing from the use of frames, and obtain optimal compression. This technique makes use of block Hankel operators and so-called AAK theory [75,76].Recently, we introduced a method for frequency extrapolation assuming that the data (in time) can be sparsely represented by a sum of Lorentzian functions. This condition restricts, for example, the ‘density’ of reflections in a given time window. The current implementation of this technique makes use of knowledge of the band-pass filter, which requirement needs to be further relaxed. The extrapolation provides insight in the ultimate information contained in finite-frequency data.

    Via compression, we can also recover geometrical information (slopes, curvatures, etc.) connecting points in phase space corresponding with the significant wavepacket coefficients using graph cuts. These appeared to allow conflicting slopes in the data and proved to be a powerful tool in interpolation.

  • Gradients, nonlinear diffusion. Through multi-component data we have implicit information about the spatial derivatives (and, hence, directionality) of the wavefield. We can use these derivatives to form structure tensors (which we can extend to allow for conflicting slopes). We constructed a nonlinear anisotropic diffusion equation using these tensors; the associated process provides effectively a new interpolation concept.

The Dirichlet-to-Neumann map represents the data in the analysis of the inverse boundary value problem for seismic waves. This map models observations when acting on a boundary source, or equivalently simultaneous sources, which can be synthesized from actual sources in multiple ways. We have studied the probing of this Dirichlet-to-Neumann map in the frequency domain, carried out an analysis and obtained conditions and a procedure to enable accurate low-rank approximations. This property lends itself to time-harmonic wavefield recovery with matrix completion techniques.

Classification of seismic phases and manifold learning

The analysis of seismic waves remains a tedious task that seismologists often perform by inspection of the available seismograms. The complexity and noisy nature of the signal combined with a sparse and irregular sampling make this analysis difficult and imprecise. Yet, a detailed interpretation of the geometric contents of these datasets can provide valuable prior information for the solution of corresponding inverse problems.

This state of affair is indicative of a lack of effective and robust algorithms for the computational parsing and interpretation of seismograms. Indeed, the limited frequency content, strong nonlinearity, temporally scattered nature of these signals prove challenging to standard signal processing techniques.

In this context, recent advances in signal processing research show remarkable promise to enable the automatic characterization of the subtle and irregular patterns that are present in seismic waves. The scattering transform, in particular, through its invariance to translation and stability under small temporal distortions, constitutes a compelling new direction for the processing of seismic waves. In a pilot study, we have shown that the scattering transform combined with dimensionality reduction techniques can elucidate the fundamental organization of seismic signals. Our approach leverages ideas from signal processing and machine learning to create an efficient computational pipeline through which individual events and their constitutive wave patterns can be identified and connected across signals. We expect a wide range applications, such as the classification of tremors and identification of repeat earthquakes.

Uncertainty quantification

We developed data analysis tools for uncertainty quantification of discretized linearized inverse problems, based on l2– andl1-regularization [88]. These can be used to reveal the presence of systematic errors or to validate the assumed statistical model. Our methods include bounds on the performance of randomized estimators of large matrices, bounds for the bias, and resampling methods for model validation.

Wave Packets: Nonlinear Approximation, Compression, Sparse Sampling

Multi-scale physics of seismic waves

In building a detailed understanding of the propagation and scattering of waves, that is, their interaction with complex media and boundaries, I have followed a hierarchical description of medium variations or coefficients, starting with a smooth component, a C1,1 component, general interfaces with regularity constraints, and a class of random fluctuations.

  • Multi-scale ray-wave duality. The description of waves up to C1,1 coefficients is naturally obtained introducing the notion of wavepackets. Wavepackets were introduced in harmonic analysis and later identified as curvelets. A wavepacket can be viewed as a localized ‘fat’ plane wave, or a ‘quantum’. Wavepackets follow a scale description (dyadic parabolic decomposition of phase space) and can be used to form a transform pair (generalizing in some sense the notion of beamforming). We obtained a construction of the solution to the wave equation in two steps: in step (i) we smooth the medium on the scales corresponding with the scales of the wavepackets and obtain an approximate construction using geometry (rays on each scale); in step (ii) we incorporate the scattering via a Volterra integral equation. We obtain a concentration of wave packets result; the degree of concentration reveals the smoothness of the medium [65]. The construction also answers precisely questions about how asymptotic representations are related to full wave solutions.

    We generalized the construction and description from scalar waves to elastic waves, and addressed the question of the minimum regularity in Lame´ parameters for which polarized waves can be distinguished, which appears to be C1,1[86].

  • Multi-scale wave-equation transmission tomography. The above construction can be utilized to analyze the physics behind cross-correlation-based wave-equation transmission tomography, that is, the ‘finite-frequency’ analogue of the geodesic ray transform, by taking Fréchet-type derivatives of our solution construction [78], and resolves the seeming insensitivity of the transform along the unperturbed rays via a multi-scale analysis. Indeed, in the case of smooth wavespeeds, and the limit of infinite bandwdith, the geodesic ray transform is recovered [43].
  • Scattering series, interfaces. Reducing the regularity of the medium further, leads us to the introduction of (general) interfaces (conormal singularities). We succeeded to track the order of scattering through the regularity of the wave constituents in the case of C1 (ε > 0) coefficients [108]. In this case, the reflected wave is more regular than the incident wave in a Sobolev sense. Via this analysis, we also obtained the first meaningful introduction of a scattering series for waves (beyond the case of planarly layered media).
  • Random fluctuations, moments, cross correlations. Via random fluctuations we reduce the regularity of the medium further. Assuming a wavepacket decomposition of the source, we studied waves in random media essentially via a description in terms of moments. Through the parabolic scaling underlying the construction of wavepackets, one obtains Itô-Schrödinger diffusion processes which reproduce the moments, in the limit of fine scales. We incorporated deterministic interfaces as well, moderately tilted relative to the orientation of the wavepacket source, and obtained theoretically an extended notion of enhanced backscattering [101]. We found that the backscattering cone does not depend on the tilt.

    We utilized the results pertaining to the Itô-Schrödinger diffusion processes and partly coherent waves to study the retrieval of Green’s functions from ‘field-field’ cross correlations in the strongly cluttered regime. We obtain filters via Wigner transforms which satisfy a set of transport equations, and determine their asymptotic behaviors which contain information about the statistics of the random fluctuations. We found that the retrieval is partial only, contains these (blurring) filters, and depends on the directionality of the source [69,102].

    • We studied the ‘field-field’ correlations from data acquired on the Tibetan plateau when focussing on surface waves. We assessed the retrieval from coda and from ambient noise by comparing the results with direct surface wave observations from earthquakes [54] (including directionality [73]). We carried the comparison over to the estimation of phase velocities and enabled a fusion of the results while enlarging the bandwidth. Here, we employed a semi-classical analysis point of view inserting a small parameter in the variations in medium properties in the surface normal coordinate.
    • In highly discontinuous and random planarly layered media, we developed an understanding of tunneling, hidden equivalent medium averaging, and microlocal attenuation which cannot originate from standard visco-elastic models [1,2,5,8].

On planetary scale, that is, a ball, we studied a random radially layered medium including random angular fluctuations on distinct scales as well as deterministic discontinuities. We obtained a stochastic analysis in terms of modes and partly coherent waves excited by an interior source. Through cross correlations of the trace (field at the surface) we obtained a (virtual) data set which can be utilized for the imaging of the deterministic discontinuities.

  • Directional wavefield decomposition and propagation, sufficiently smooth coefficients. In the case of sufficient medium regularity one can interchange time locally with a (curvilinear) spatial coordinate (‘pseudodepth’) to describe the propagation of waves in terms of an evolution or one-way wave equation. In certain non-trivial wavespeed profiles this can be accomplised exactly, using spectral theory (modes) and special functions [23]. Under certain conditions, we can incorporate and analyze the role of ‘evanescent’ wave constituents [24]. In general, this technique is desgned for relatively high frequency wave propagation. In appropriate coordinates, we derived a novel paraxial equation and obtained estimates for its solution in terms of an energy norm converging to the full wave solution as the frequency becomes large. This is accomplished, again, via parabolic scaling. The ‘up-down’ coupling can be incorporated through a generalized Bremmer series [11]. These consistent paraxial equations enable the development of fast algorithms [15], and are of current interest in, for example, long-range propagation of and imaging with teleseismic body waves.
  • Anisotropy. We studied medium anisotropy from the point of view of polarized waves and a corresponding local slowness surface characterization (for example, by introducing anellipticity) and expansions [20], and the point of view of classification, characterization and construction of waveforms via Picard-Fuchs equations [7].

The behaviors and phenomena described above provide fundamental insight in developing multi-scale methods to estimate local continuity, anisotropy, regularity etc. of medium properties and substructures from broadband seismic observations, and can lead to a deeper understanding of the coupling of geological and geodynamical processes, (complex) phase boundaries etc. and the interaction (and, hence, probing) with seismic waves. In this context, we obtained some basic results in the case of viscous flow in the mantle wedge induced by a subducting lithospheric plate; we approximated the subduction of the Philippine Sea plate beneath Eurasia, explored the effects of different rheological and thermal constraints, and modelled the evolution of the finite strain ellipses [61].

Algorithms

Following the multi-scale physics of waves a family of new computational algorithms has emerged:

  • We designed and implemented a fast discrete almost symmetric wavepacket transform pair in three dimensions [80]. This transform pair can also be utilized for wave-based signal processing. Waves (as well as coefficient models containing interfaces) can be compressed with this transform.
  • Following the dyadic parabolic decomposition of phase space underlying the construction of wave packets, we developed a fast algorithm (with error estimates) for a class of so-called Fourier integral operators comprising propagators, asymptotic solution operators of the wave equations and (inverse) scattering transforms [92]. This algorithm is based on accurate low-rank separated representations, prolate spheroidal wave functions, and ray-geometric computations. (From a physics of waves point of view, in the case of one-way wave propagators, these representations are related to generalizations of phase screens [34,35].) In the application to the half-wave equation, we obtained a theoretical, multi-scale description of the Gouy phase which affects the extraction of travel time in long-range propagation. We designed a special procedure to incorporate the (numerical) generation of general caustics[104].
  • We estimated the regularity of the kernel of the Volterra operator in the above mentioned construction of weak, full-wave solutions in media of limited smoothness. This regularity determined the optimal choice of quadrature to discretize the integral. Moreover we estimated the approximation of solutions by finite sets or clusters of wave packets to complete the discretization of the Volterra equation and obtain a finite matrix representation [95].

We constructed frames of wave packets produced by parabolic dilation, rotation and translation of (a finite sum of) Gaussians and obtained asymptotics on the analogue of Daubechies’ frame criterion [110]. We also obtained a frame of wave-atom-like packets generated from pure Gaussians; these can be identified as multi-scale Gaussian beams. We developed an associated parametrix construction and algorithm.

Extended Isochron Rays

Microlocal analysis of seismic body waves and linearized inverse problems – Imaging

A separated scale representation of variations in medium properties leads to the introduction of a background model and a contrast or reflectivity. This in turn leads to the introduction and formulation of a separated inverse problem for estimating these, and is tuned to particular (relative high) frequency windows in the data. This formulation encompasses inverse scattering and ‘wave-equation’ tomography and is focussed on the singularities in the data. In this framework, the phases or constituents may be multiply scattered but can interact with the unknown contrast only once.

One key complication in developing proper imaging techniques for transmission and reflection tomography and inverse scattering are the generation of caustics. This is typically associated with regions of low wavespeeds, for example, hot spots and mantle upwellings. A related complication is associated with anisotropy. To avoid misinterpretation of images, it is critical to understand uniqueness of reconstruction (even though full data coverage will not be available), the possible generation of (singular) artifacts (beyond standard resolution analysis), and partial reconstruction (including illumination). One may think of these techniques, which have grown out of traditional tomography and (Kirchhoff) migration, as providing a ‘leading-order’ identification of variations in material properties and structures in Earth’s interior.

Inverse scattering

Even though maximal (5-dimensional) surface data seem to define an overdetermined inverse problem, this is, in fact, not the case in the presence of caustics. In reality, this inverse problem is further complicated by the availability of data (for example, linear arrays) often constituting a lower-dimensional set.

  • Characterization of the normal operator, extensions. The general characterization of the normal operator associated with inverse scattering is still incomplete with so-called Ip,l classes appearing in the analysis. To ensure that false discontinuities cannot be generated, we invoked the Bolker condition [29], which may be violated in particular if the acquisition loses a dimension. Possible mitigation of artifacts caused by such a violation is a current subject of research. We pioneered the introduction of extensions leading to (microlocally) invertible single scattering operators. The propagation of singularities by these extensions are described by canonical transformations with associated evolution equations. With these extensions the imaging of reflectivity in inaccurate background models can be undone; beyond that, we obtained an evolution equation (and an associated Hamiltonian) for extended images under changing background velocities [72]. The extensions play a key role in obtaining estimates for the background model. We constructed explicit extensions via local directional decomposition subject to the introduction of a so-calledcurvilinear double-square-root (DSR) condition [46,74]. Our standard extension involves interior curvilinear offset; however, as an alternative we constructed a (wave-equation) angle transform, which can be related (microlocally) to reflection coefficients [41].
  • Annihilators, wave-equation reflection tomography. We introduced annihilators (extending the notion of differential semblance to allow for the formation of caustics in the unknown background) of the data [55]. Annihilators characterize the range of the single scattering operator [52] and provide a procedure for conditional source-receiver continuation. The range can be used as a criterion for background model recovery, which is under certain conditions unique. We designed an iterative method for this recovery and analyzed the associated sensitivity kernels [57].

    We adapted and applied our methodology using the generalized Radon transform jointly for PP and PS reflections to constrain TI parameters in a sedimentary environment in the Norwegian sector of the North Sea [48].

  • Generalized Radon transform (GRT). We have constructed a GRT for inverse scattering incorporating the asymptotic inverse (parametrix) of the normal operator, in anisotropic media, assuming that the equations for elastic waves are of principal type [22] and carried out an associated resolution analysis [17]. Indeed, for the purpose of imaging Earth’s deep interior, the GRT can be applied to global network data, in principle; one does not need uniform spatial sampling as we demonstrated by invoking quasi-Monte Carlo techniques [10]. Moreover, upon adapting the order of the Fourier integral operator representing the GRT, we obtained stability and a method for statistical inference of singularities [62]. We developed an approach to precise partial reconstruction from partial data with the GRT and illumination correction using a frame of wavepackets [70].
    • We applied the GRT to image the Dll layer beneath Central America using ScS (precursors), and provided a new view of the post-perovskite phase transition near the CMB. With thermodynamic properties of phase transitions in mantle silicates, we interpreted the images and estimated in situ temperatures. Exploiting the presence of multiple phase transitions, we estimated the core heat flux in the coldest region (a site of deep subduction) and in a neighborhood away from it [59].
    • We also applied the GRT to image the hot spot underneath Hawaii using underside reflections, that is, SS precursors [90]. The detection of a thermal plume has been seismically inconclusive. We explained the depths of the discontinuities in the transition zone below the Central Pacific with olivine and garnet transitions in a pyrolitic mantle. The presence of a thermal anomaly west of Hawaii suggest that hot material accumulates near the base of the transition zone before being entrained in flow toward Hawaii.
    • The modern view of Earths lowermost mantle considers a Dll region of enhanced (seismologically inferred) heterogeneity bounded by the core mantle boundary (CMB) and an interface some 150-300 km above it, with the latter attributed to the post-perovskite phase transition. We carried out a deep-Earth exploration with ScS and SKKS wavefields, which probe this region from above and below, respectively, which suggested that this view be modified. Inverse scattering revealed that, beneath Central America and East Asia, two areas known for subduction of oceanic plates deep into Earths mantle, seismic reflectors exist above the conventional Dllregion [107]. The occurrence of multiple interfaces is inconsistent with expectations from a thermal response of a single isochemical post-perovskite transition but can be explained with post-perovskite transitions in differentiated slab materials. Candidate compositions for a seismically detectable post-perovskite transition at pressures less than the CMB include mid-oceanic ridge basalt (MORB) and harzburgite components of differentiated oceanic lithosphere transported to the lowermost mantle by subduction.
  • Reverse-time-continuation (RTM)-based inverse scattering. We constructed a new inverse scattering transform based on RTM [93], also in the anisotropic case [chapter in book 5], which automatically removes the so-called low-frequency artifacts. This transform applies to data generated by a finite set of (point) sources (events). The formulation in curvilinear coordinates is straightforward and an elastic wave-equation driven imaging procedure is obtained.  However, we also designed a fast, approximate algorithm with associated estimates using the dyadic parabolic decomposition of phase space and bypassing time stepping [preprint 1]. In this algorithm it is straightforward to extract information about the local incidence angle and thus generate angle-dependent images.
    • We can remove the knowledge of the source in the case of converted waves. In [99] we obtained an inverse scattering procedure which is effectively bilinear in the data and replaces traditional receiver functions, with the resolving power of RTM-based inverse scattering.
    • We extended our inverse scattering procedure to data from (known) deterministic to (known or unknown, ambient) noise sources [98] ensuring statistical stability. Exploiting mode conversions, we obtained a related result using the coda (only) of an earthquake.

These procedures, in combination with advanced wavefield retrieval techniques, provide a new platform honoring wave dynamics, for example, for high-resolution imaging and improved detection and characterization of interfaces in the upper mantle, including the Moho, the lithosphere-asthenosphere boundary and the 660km discontinuity, as well as low velocity features underneath North America using USArray (including temporary and ‘flexible’ array) data.

  • Downward-continuation-based inverse scattering. The downward continuation concept is based on directional wavefield decomposition and its implementation has its origin in the so-called double-square-root (DSR) equation. We introduced curvilinear coordinates and a pseudodepth via a Riemannian metric. We constructed the normal operators for the above mentioned extended imaging procedures in this downward continuation approach [55]. We also constructed explicitly the corresponding evolution equation (and thus an extended exploding reflector model) which provided a new insight in the geometry of extended imaging in the presence of caustics via the introduction of extended isochron rays [82].

 (Wave-equation) tomography

  • Condition for uniqueness, multiple scattering. In the case of transmission tomography it is well known that the presence of caustics obstructs the unique recovery of velocities. With multiply broken rays and using the scattering relation as the data – with a finite set of favorably located (unknown) reflectors – the linearized inverse problem and (partial) reconstruction allows for more general background models in which certain caustics may form.
  • Anisotropy, shear-wave splitting. We extended our analysis of wave-equation transmisison tomography to shear-wave splitting intensity tomography. We designed, again, an iterative method in a heterogeneous environment [64], and carried out computational experiments for SKS and SKKS phases to demonstrate the possibility of constraining the deformation and flow beneath the Ryukyu arc, Japan. Flow modelling is an essential part of the regularization in estimating local anisotropy.

Downward Continuation Based Extended Imaging

Microlocal analysis of seismic body waves and linearized inverse problems – Imaging

A separated scale representation of variations in medium properties leads to the introduction of a background model and a contrast or reflectivity. This in turn leads to the introduction and formulation of a separated inverse problem for estimating these, and is tuned to particular (relative high) frequency windows in the data. This formulation encompasses inverse scattering and ‘wave-equation’ tomography and is focussed on the singularities in the data. In this framework, the phases or constituents may be multiply scattered but can interact with the unknown contrast only once.

One key complication in developing proper imaging techniques for transmission and reflection tomography and inverse scattering are the generation of caustics. This is typically associated with regions of low wavespeeds, for example, hot spots and mantle upwellings. A related complication is associated with anisotropy. To avoid misinterpretation of images, it is critical to understand uniqueness of reconstruction (even though full data coverage will not be available), the possible generation of (singular) artifacts (beyond standard resolution analysis), and partial reconstruction (including illumination). One may think of these techniques, which have grown out of traditional tomography and (Kirchhoff) migration, as providing a ‘leading-order’ identification of variations in material properties and structures in Earth’s interior.

Inverse scattering

Even though maximal (5-dimensional) surface data seem to define an overdetermined inverse problem, this is, in fact, not the case in the presence of caustics. In reality, this inverse problem is further complicated by the availability of data (for example, linear arrays) often constituting a lower-dimensional set.

  • Characterization of the normal operator, extensions. The general characterization of the normal operator associated with inverse scattering is still incomplete with so-called Ip,l classes appearing in the analysis. To ensure that false discontinuities cannot be generated, we invoked the Bolker condition [29], which may be violated in particular if the acquisition loses a dimension. Possible mitigation of artifacts caused by such a violation is a current subject of research. We pioneered the introduction of extensions leading to (microlocally) invertible single scattering operators. The propagation of singularities by these extensions are described by canonical transformations with associated evolution equations. With these extensions the imaging of reflectivity in inaccurate background models can be undone; beyond that, we obtained an evolution equation (and an associated Hamiltonian) for extended images under changing background velocities [72]. The extensions play a key role in obtaining estimates for the background model. We constructed explicit extensions via local directional decomposition subject to the introduction of a so-calledcurvilinear double-square-root (DSR) condition [46,74]. Our standard extension involves interior curvilinear offset; however, as an alternative we constructed a (wave-equation) angle transform, which can be related (microlocally) to reflection coefficients [41].
  • Annihilators, wave-equation reflection tomography. We introduced annihilators (extending the notion of differential semblance to allow for the formation of caustics in the unknown background) of the data [55]. Annihilators characterize the range of the single scattering operator [52] and provide a procedure for conditional source-receiver continuation. The range can be used as a criterion for background model recovery, which is under certain conditions unique. We designed an iterative method for this recovery and analyzed the associated sensitivity kernels [57].

    We adapted and applied our methodology using the generalized Radon transform jointly for PP and PS reflections to constrain TI parameters in a sedimentary environment in the Norwegian sector of the North Sea [48].

  • Generalized Radon transform (GRT). We have constructed a GRT for inverse scattering incorporating the asymptotic inverse (parametrix) of the normal operator, in anisotropic media, assuming that the equations for elastic waves are of principal type [22] and carried out an associated resolution analysis [17]. Indeed, for the purpose of imaging Earth’s deep interior, the GRT can be applied to global network data, in principle; one does not need uniform spatial sampling as we demonstrated by invoking quasi-Monte Carlo techniques [10]. Moreover, upon adapting the order of the Fourier integral operator representing the GRT, we obtained stability and a method for statistical inference of singularities [62]. We developed an approach to precise partial reconstruction from partial data with the GRT and illumination correction using a frame of wavepackets [70].
    • We applied the GRT to image the Dll layer beneath Central America using ScS (precursors), and provided a new view of the post-perovskite phase transition near the CMB. With thermodynamic properties of phase transitions in mantle silicates, we interpreted the images and estimated in situ temperatures. Exploiting the presence of multiple phase transitions, we estimated the core heat flux in the coldest region (a site of deep subduction) and in a neighborhood away from it [59].
    • We also applied the GRT to image the hot spot underneath Hawaii using underside reflections, that is, SS precursors [90]. The detection of a thermal plume has been seismically inconclusive. We explained the depths of the discontinuities in the transition zone below the Central Pacific with olivine and garnet transitions in a pyrolitic mantle. The presence of a thermal anomaly west of Hawaii suggest that hot material accumulates near the base of the transition zone before being entrained in flow toward Hawaii.
    • The modern view of Earths lowermost mantle considers a Dll region of enhanced (seismologically inferred) heterogeneity bounded by the core mantle boundary (CMB) and an interface some 150-300 km above it, with the latter attributed to the post-perovskite phase transition. We carried out a deep-Earth exploration with ScS and SKKS wavefields, which probe this region from above and below, respectively, which suggested that this view be modified. Inverse scattering revealed that, beneath Central America and East Asia, two areas known for subduction of oceanic plates deep into Earths mantle, seismic reflectors exist above the conventional Dllregion [107]. The occurrence of multiple interfaces is inconsistent with expectations from a thermal response of a single isochemical post-perovskite transition but can be explained with post-perovskite transitions in differentiated slab materials. Candidate compositions for a seismically detectable post-perovskite transition at pressures less than the CMB include mid-oceanic ridge basalt (MORB) and harzburgite components of differentiated oceanic lithosphere transported to the lowermost mantle by subduction.
  • Reverse-time-continuation (RTM)-based inverse scattering. We constructed a new inverse scattering transform based on RTM [93], also in the anisotropic case [chapter in book 5], which automatically removes the so-called low-frequency artifacts. This transform applies to data generated by a finite set of (point) sources (events). The formulation in curvilinear coordinates is straightforward and an elastic wave-equation driven imaging procedure is obtained.  However, we also designed a fast, approximate algorithm with associated estimates using the dyadic parabolic decomposition of phase space and bypassing time stepping [preprint 1]. In this algorithm it is straightforward to extract information about the local incidence angle and thus generate angle-dependent images.
    • We can remove the knowledge of the source in the case of converted waves. In [99] we obtained an inverse scattering procedure which is effectively bilinear in the data and replaces traditional receiver functions, with the resolving power of RTM-based inverse scattering.
    • We extended our inverse scattering procedure to data from (known) deterministic to (known or unknown, ambient) noise sources [98] ensuring statistical stability. Exploiting mode conversions, we obtained a related result using the coda (only) of an earthquake.

These procedures, in combination with advanced wavefield retrieval techniques, provide a new platform honoring wave dynamics, for example, for high-resolution imaging and improved detection and characterization of interfaces in the upper mantle, including the Moho, the lithosphere-asthenosphere boundary and the 660km discontinuity, as well as low velocity features underneath North America using USArray (including temporary and ‘flexible’ array) data.

  • Downward-continuation-based inverse scattering. The downward continuation concept is based on directional wavefield decomposition and its implementation has its origin in the so-called double-square-root (DSR) equation. We introduced curvilinear coordinates and a pseudodepth via a Riemannian metric. We constructed the normal operators for the above mentioned extended imaging procedures in this downward continuation approach [55]. We also constructed explicitly the corresponding evolution equation (and thus an extended exploding reflector model) which provided a new insight in the geometry of extended imaging in the presence of caustics via the introduction of extended isochron rays [82].

 (Wave-equation) tomography

  • Condition for uniqueness, multiple scattering. In the case of transmission tomography it is well known that the presence of caustics obstructs the unique recovery of velocities. With multiply broken rays and using the scattering relation as the data – with a finite set of favorably located (unknown) reflectors – the linearized inverse problem and (partial) reconstruction allows for more general background models in which certain caustics may form.
  • Anisotropy, shear-wave splitting. We extended our analysis of wave-equation transmisison tomography to shear-wave splitting intensity tomography. We designed, again, an iterative method in a heterogeneous environment [64], and carried out computational experiments for SKS and SKKS phases to demonstrate the possibility of constraining the deformation and flow beneath the Ryukyu arc, Japan. Flow modelling is an essential part of the regularization in estimating local anisotropy.

Thermal Anomaly, Transition Zone, Hawaii

Microlocal analysis of seismic body waves and linearized inverse problems – Imaging

A separated scale representation of variations in medium properties leads to the introduction of a background model and a contrast or reflectivity. This in turn leads to the introduction and formulation of a separated inverse problem for estimating these, and is tuned to particular (relative high) frequency windows in the data. This formulation encompasses inverse scattering and ‘wave-equation’ tomography and is focussed on the singularities in the data. In this framework, the phases or constituents may be multiply scattered but can interact with the unknown contrast only once.

One key complication in developing proper imaging techniques for transmission and reflection tomography and inverse scattering are the generation of caustics. This is typically associated with regions of low wavespeeds, for example, hot spots and mantle upwellings. A related complication is associated with anisotropy. To avoid misinterpretation of images, it is critical to understand uniqueness of reconstruction (even though full data coverage will not be available), the possible generation of (singular) artifacts (beyond standard resolution analysis), and partial reconstruction (including illumination). One may think of these techniques, which have grown out of traditional tomography and (Kirchhoff) migration, as providing a ‘leading-order’ identification of variations in material properties and structures in Earth’s interior.

Inverse scattering

Even though maximal (5-dimensional) surface data seem to define an overdetermined inverse problem, this is, in fact, not the case in the presence of caustics. In reality, this inverse problem is further complicated by the availability of data (for example, linear arrays) often constituting a lower-dimensional set.

  • Characterization of the normal operator, extensions. The general characterization of the normal operator associated with inverse scattering is still incomplete with so-called Ip,l classes appearing in the analysis. To ensure that false discontinuities cannot be generated, we invoked the Bolker condition [29], which may be violated in particular if the acquisition loses a dimension. Possible mitigation of artifacts caused by such a violation is a current subject of research. We pioneered the introduction of extensions leading to (microlocally) invertible single scattering operators. The propagation of singularities by these extensions are described by canonical transformations with associated evolution equations. With these extensions the imaging of reflectivity in inaccurate background models can be undone; beyond that, we obtained an evolution equation (and an associated Hamiltonian) for extended images under changing background velocities [72]. The extensions play a key role in obtaining estimates for the background model. We constructed explicit extensions via local directional decomposition subject to the introduction of a so-calledcurvilinear double-square-root (DSR) condition [46,74]. Our standard extension involves interior curvilinear offset; however, as an alternative we constructed a (wave-equation) angle transform, which can be related (microlocally) to reflection coefficients [41].
  • Annihilators, wave-equation reflection tomography. We introduced annihilators (extending the notion of differential semblance to allow for the formation of caustics in the unknown background) of the data [55]. Annihilators characterize the range of the single scattering operator [52] and provide a procedure for conditional source-receiver continuation. The range can be used as a criterion for background model recovery, which is under certain conditions unique. We designed an iterative method for this recovery and analyzed the associated sensitivity kernels [57].

    We adapted and applied our methodology using the generalized Radon transform jointly for PP and PS reflections to constrain TI parameters in a sedimentary environment in the Norwegian sector of the North Sea [48].

  • Generalized Radon transform (GRT). We have constructed a GRT for inverse scattering incorporating the asymptotic inverse (parametrix) of the normal operator, in anisotropic media, assuming that the equations for elastic waves are of principal type [22] and carried out an associated resolution analysis [17]. Indeed, for the purpose of imaging Earth’s deep interior, the GRT can be applied to global network data, in principle; one does not need uniform spatial sampling as we demonstrated by invoking quasi-Monte Carlo techniques [10]. Moreover, upon adapting the order of the Fourier integral operator representing the GRT, we obtained stability and a method for statistical inference of singularities [62]. We developed an approach to precise partial reconstruction from partial data with the GRT and illumination correction using a frame of wavepackets [70].
    • We applied the GRT to image the Dll layer beneath Central America using ScS (precursors), and provided a new view of the post-perovskite phase transition near the CMB. With thermodynamic properties of phase transitions in mantle silicates, we interpreted the images and estimated in situ temperatures. Exploiting the presence of multiple phase transitions, we estimated the core heat flux in the coldest region (a site of deep subduction) and in a neighborhood away from it [59].
    • We also applied the GRT to image the hot spot underneath Hawaii using underside reflections, that is, SS precursors [90]. The detection of a thermal plume has been seismically inconclusive. We explained the depths of the discontinuities in the transition zone below the Central Pacific with olivine and garnet transitions in a pyrolitic mantle. The presence of a thermal anomaly west of Hawaii suggest that hot material accumulates near the base of the transition zone before being entrained in flow toward Hawaii.
    • The modern view of Earths lowermost mantle considers a Dll region of enhanced (seismologically inferred) heterogeneity bounded by the core mantle boundary (CMB) and an interface some 150-300 km above it, with the latter attributed to the post-perovskite phase transition. We carried out a deep-Earth exploration with ScS and SKKS wavefields, which probe this region from above and below, respectively, which suggested that this view be modified. Inverse scattering revealed that, beneath Central America and East Asia, two areas known for subduction of oceanic plates deep into Earths mantle, seismic reflectors exist above the conventional Dllregion [107]. The occurrence of multiple interfaces is inconsistent with expectations from a thermal response of a single isochemical post-perovskite transition but can be explained with post-perovskite transitions in differentiated slab materials. Candidate compositions for a seismically detectable post-perovskite transition at pressures less than the CMB include mid-oceanic ridge basalt (MORB) and harzburgite components of differentiated oceanic lithosphere transported to the lowermost mantle by subduction.
  • Reverse-time-continuation (RTM)-based inverse scattering. We constructed a new inverse scattering transform based on RTM [93], also in the anisotropic case [chapter in book 5], which automatically removes the so-called low-frequency artifacts. This transform applies to data generated by a finite set of (point) sources (events). The formulation in curvilinear coordinates is straightforward and an elastic wave-equation driven imaging procedure is obtained.  However, we also designed a fast, approximate algorithm with associated estimates using the dyadic parabolic decomposition of phase space and bypassing time stepping [preprint 1]. In this algorithm it is straightforward to extract information about the local incidence angle and thus generate angle-dependent images.
    • We can remove the knowledge of the source in the case of converted waves. In [99] we obtained an inverse scattering procedure which is effectively bilinear in the data and replaces traditional receiver functions, with the resolving power of RTM-based inverse scattering.
    • We extended our inverse scattering procedure to data from (known) deterministic to (known or unknown, ambient) noise sources [98] ensuring statistical stability. Exploiting mode conversions, we obtained a related result using the coda (only) of an earthquake.

These procedures, in combination with advanced wavefield retrieval techniques, provide a new platform honoring wave dynamics, for example, for high-resolution imaging and improved detection and characterization of interfaces in the upper mantle, including the Moho, the lithosphere-asthenosphere boundary and the 660km discontinuity, as well as low velocity features underneath North America using USArray (including temporary and ‘flexible’ array) data.

  • Downward-continuation-based inverse scattering. The downward continuation concept is based on directional wavefield decomposition and its implementation has its origin in the so-called double-square-root (DSR) equation. We introduced curvilinear coordinates and a pseudodepth via a Riemannian metric. We constructed the normal operators for the above mentioned extended imaging procedures in this downward continuation approach [55]. We also constructed explicitly the corresponding evolution equation (and thus an extended exploding reflector model) which provided a new insight in the geometry of extended imaging in the presence of caustics via the introduction of extended isochron rays [82].

 (Wave-equation) tomography

  • Condition for uniqueness, multiple scattering. In the case of transmission tomography it is well known that the presence of caustics obstructs the unique recovery of velocities. With multiply broken rays and using the scattering relation as the data – with a finite set of favorably located (unknown) reflectors – the linearized inverse problem and (partial) reconstruction allows for more general background models in which certain caustics may form.
  • Anisotropy, shear-wave splitting. We extended our analysis of wave-equation transmisison tomography to shear-wave splitting intensity tomography. We designed, again, an iterative method in a heterogeneous environment [64], and carried out computational experiments for SKS and SKKS phases to demonstrate the possibility of constraining the deformation and flow beneath the Ryukyu arc, Japan. Flow modelling is an essential part of the regularization in estimating local anisotropy.

Core-Mantle Boundary, Post-Perovskite Phase Transition and D”

Microlocal analysis of seismic body waves and linearized inverse problems – Imaging

A separated scale representation of variations in medium properties leads to the introduction of a background model and a contrast or reflectivity. This in turn leads to the introduction and formulation of a separated inverse problem for estimating these, and is tuned to particular (relative high) frequency windows in the data. This formulation encompasses inverse scattering and ‘wave-equation’ tomography and is focussed on the singularities in the data. In this framework, the phases or constituents may be multiply scattered but can interact with the unknown contrast only once.

One key complication in developing proper imaging techniques for transmission and reflection tomography and inverse scattering are the generation of caustics. This is typically associated with regions of low wavespeeds, for example, hot spots and mantle upwellings. A related complication is associated with anisotropy. To avoid misinterpretation of images, it is critical to understand uniqueness of reconstruction (even though full data coverage will not be available), the possible generation of (singular) artifacts (beyond standard resolution analysis), and partial reconstruction (including illumination). One may think of these techniques, which have grown out of traditional tomography and (Kirchhoff) migration, as providing a ‘leading-order’ identification of variations in material properties and structures in Earth’s interior.

Inverse scattering

Even though maximal (5-dimensional) surface data seem to define an overdetermined inverse problem, this is, in fact, not the case in the presence of caustics. In reality, this inverse problem is further complicated by the availability of data (for example, linear arrays) often constituting a lower-dimensional set.

  • Characterization of the normal operator, extensions. The general characterization of the normal operator associated with inverse scattering is still incomplete with so-called Ip,l classes appearing in the analysis. To ensure that false discontinuities cannot be generated, we invoked the Bolker condition [29], which may be violated in particular if the acquisition loses a dimension. Possible mitigation of artifacts caused by such a violation is a current subject of research. We pioneered the introduction of extensions leading to (microlocally) invertible single scattering operators. The propagation of singularities by these extensions are described by canonical transformations with associated evolution equations. With these extensions the imaging of reflectivity in inaccurate background models can be undone; beyond that, we obtained an evolution equation (and an associated Hamiltonian) for extended images under changing background velocities [72]. The extensions play a key role in obtaining estimates for the background model. We constructed explicit extensions via local directional decomposition subject to the introduction of a so-calledcurvilinear double-square-root (DSR) condition [46,74]. Our standard extension involves interior curvilinear offset; however, as an alternative we constructed a (wave-equation) angle transform, which can be related (microlocally) to reflection coefficients [41].
  • Annihilators, wave-equation reflection tomography. We introduced annihilators (extending the notion of differential semblance to allow for the formation of caustics in the unknown background) of the data [55]. Annihilators characterize the range of the single scattering operator [52] and provide a procedure for conditional source-receiver continuation. The range can be used as a criterion for background model recovery, which is under certain conditions unique. We designed an iterative method for this recovery and analyzed the associated sensitivity kernels [57].

    We adapted and applied our methodology using the generalized Radon transform jointly for PP and PS reflections to constrain TI parameters in a sedimentary environment in the Norwegian sector of the North Sea [48].

  • Generalized Radon transform (GRT). We have constructed a GRT for inverse scattering incorporating the asymptotic inverse (parametrix) of the normal operator, in anisotropic media, assuming that the equations for elastic waves are of principal type [22] and carried out an associated resolution analysis [17]. Indeed, for the purpose of imaging Earth’s deep interior, the GRT can be applied to global network data, in principle; one does not need uniform spatial sampling as we demonstrated by invoking quasi-Monte Carlo techniques [10]. Moreover, upon adapting the order of the Fourier integral operator representing the GRT, we obtained stability and a method for statistical inference of singularities [62]. We developed an approach to precise partial reconstruction from partial data with the GRT and illumination correction using a frame of wavepackets [70].
    • We applied the GRT to image the Dll layer beneath Central America using ScS (precursors), and provided a new view of the post-perovskite phase transition near the CMB. With thermodynamic properties of phase transitions in mantle silicates, we interpreted the images and estimated in situ temperatures. Exploiting the presence of multiple phase transitions, we estimated the core heat flux in the coldest region (a site of deep subduction) and in a neighborhood away from it [59].
    • We also applied the GRT to image the hot spot underneath Hawaii using underside reflections, that is, SS precursors [90]. The detection of a thermal plume has been seismically inconclusive. We explained the depths of the discontinuities in the transition zone below the Central Pacific with olivine and garnet transitions in a pyrolitic mantle. The presence of a thermal anomaly west of Hawaii suggest that hot material accumulates near the base of the transition zone before being entrained in flow toward Hawaii.
    • The modern view of Earths lowermost mantle considers a Dll region of enhanced (seismologically inferred) heterogeneity bounded by the core mantle boundary (CMB) and an interface some 150-300 km above it, with the latter attributed to the post-perovskite phase transition. We carried out a deep-Earth exploration with ScS and SKKS wavefields, which probe this region from above and below, respectively, which suggested that this view be modified. Inverse scattering revealed that, beneath Central America and East Asia, two areas known for subduction of oceanic plates deep into Earths mantle, seismic reflectors exist above the conventional Dllregion [107]. The occurrence of multiple interfaces is inconsistent with expectations from a thermal response of a single isochemical post-perovskite transition but can be explained with post-perovskite transitions in differentiated slab materials. Candidate compositions for a seismically detectable post-perovskite transition at pressures less than the CMB include mid-oceanic ridge basalt (MORB) and harzburgite components of differentiated oceanic lithosphere transported to the lowermost mantle by subduction.
  • Reverse-time-continuation (RTM)-based inverse scattering. We constructed a new inverse scattering transform based on RTM [93], also in the anisotropic case [chapter in book 5], which automatically removes the so-called low-frequency artifacts. This transform applies to data generated by a finite set of (point) sources (events). The formulation in curvilinear coordinates is straightforward and an elastic wave-equation driven imaging procedure is obtained.  However, we also designed a fast, approximate algorithm with associated estimates using the dyadic parabolic decomposition of phase space and bypassing time stepping [preprint 1]. In this algorithm it is straightforward to extract information about the local incidence angle and thus generate angle-dependent images.
    • We can remove the knowledge of the source in the case of converted waves. In [99] we obtained an inverse scattering procedure which is effectively bilinear in the data and replaces traditional receiver functions, with the resolving power of RTM-based inverse scattering.
    • We extended our inverse scattering procedure to data from (known) deterministic to (known or unknown, ambient) noise sources [98] ensuring statistical stability. Exploiting mode conversions, we obtained a related result using the coda (only) of an earthquake.

These procedures, in combination with advanced wavefield retrieval techniques, provide a new platform honoring wave dynamics, for example, for high-resolution imaging and improved detection and characterization of interfaces in the upper mantle, including the Moho, the lithosphere-asthenosphere boundary and the 660km discontinuity, as well as low velocity features underneath North America using USArray (including temporary and ‘flexible’ array) data.

  • Downward-continuation-based inverse scattering. The downward continuation concept is based on directional wavefield decomposition and its implementation has its origin in the so-called double-square-root (DSR) equation. We introduced curvilinear coordinates and a pseudodepth via a Riemannian metric. We constructed the normal operators for the above mentioned extended imaging procedures in this downward continuation approach [55]. We also constructed explicitly the corresponding evolution equation (and thus an extended exploding reflector model) which provided a new insight in the geometry of extended imaging in the presence of caustics via the introduction of extended isochron rays [82].

 (Wave-equation) tomography

  • Condition for uniqueness, multiple scattering. In the case of transmission tomography it is well known that the presence of caustics obstructs the unique recovery of velocities. With multiply broken rays and using the scattering relation as the data – with a finite set of favorably located (unknown) reflectors – the linearized inverse problem and (partial) reconstruction allows for more general background models in which certain caustics may form.
  • Anisotropy, shear-wave splitting. We extended our analysis of wave-equation transmisison tomography to shear-wave splitting intensity tomography. We designed, again, an iterative method in a heterogeneous environment [64], and carried out computational experiments for SKS and SKKS phases to demonstrate the possibility of constraining the deformation and flow beneath the Ryukyu arc, Japan. Flow modelling is an essential part of the regularization in estimating local anisotropy.

Scattering, Random Media and Interferometry

Multi-scale physics of seismic waves

In building a detailed understanding of the propagation and scattering of waves, that is, their interaction with complex media and boundaries, I have followed a hierarchical description of medium variations or coefficients, starting with a smooth component, a C1,1 component, general interfaces with regularity constraints, and a class of random fluctuations.

  • Multi-scale ray-wave duality. The description of waves up to C1,1 coefficients is naturally obtained introducing the notion of wavepackets. Wavepackets were introduced in harmonic analysis and later identified as curvelets. A wavepacket can be viewed as a localized ‘fat’ plane wave, or a ‘quantum’. Wavepackets follow a scale description (dyadic parabolic decomposition of phase space) and can be used to form a transform pair (generalizing in some sense the notion of beamforming). We obtained a construction of the solution to the wave equation in two steps: in step (i) we smooth the medium on the scales corresponding with the scales of the wavepackets and obtain an approximate construction using geometry (rays on each scale); in step (ii) we incorporate the scattering via a Volterra integral equation. We obtain a concentration of wave packets result; the degree of concentration reveals the smoothness of the medium [65]. The construction also answers precisely questions about how asymptotic representations are related to full wave solutions.

    We generalized the construction and description from scalar waves to elastic waves, and addressed the question of the minimum regularity in Lame´ parameters for which polarized waves can be distinguished, which appears to be C1,1[86].

  • Multi-scale wave-equation transmission tomography. The above construction can be utilized to analyze the physics behind cross-correlation-based wave-equation transmission tomography, that is, the ‘finite-frequency’ analogue of the geodesic ray transform, by taking Fréchet-type derivatives of our solution construction [78], and resolves the seeming insensitivity of the transform along the unperturbed rays via a multi-scale analysis. Indeed, in the case of smooth wavespeeds, and the limit of infinite bandwdith, the geodesic ray transform is recovered [43].
  • Scattering series, interfaces. Reducing the regularity of the medium further, leads us to the introduction of (general) interfaces (conormal singularities). We succeeded to track the order of scattering through the regularity of the wave constituents in the case of C1 (ε > 0) coefficients [108]. In this case, the reflected wave is more regular than the incident wave in a Sobolev sense. Via this analysis, we also obtained the first meaningful introduction of a scattering series for waves (beyond the case of planarly layered media).
  • Random fluctuations, moments, cross correlations. Via random fluctuations we reduce the regularity of the medium further. Assuming a wavepacket decomposition of the source, we studied waves in random media essentially via a description in terms of moments. Through the parabolic scaling underlying the construction of wavepackets, one obtains Itô-Schrödinger diffusion processes which reproduce the moments, in the limit of fine scales. We incorporated deterministic interfaces as well, moderately tilted relative to the orientation of the wavepacket source, and obtained theoretically an extended notion of enhanced backscattering [101]. We found that the backscattering cone does not depend on the tilt.

    We utilized the results pertaining to the Itô-Schrödinger diffusion processes and partly coherent waves to study the retrieval of Green’s functions from ‘field-field’ cross correlations in the strongly cluttered regime. We obtain filters via Wigner transforms which satisfy a set of transport equations, and determine their asymptotic behaviors which contain information about the statistics of the random fluctuations. We found that the retrieval is partial only, contains these (blurring) filters, and depends on the directionality of the source [69,102].

    • We studied the ‘field-field’ correlations from data acquired on the Tibetan plateau when focussing on surface waves. We assessed the retrieval from coda and from ambient noise by comparing the results with direct surface wave observations from earthquakes [54] (including directionality [73]). We carried the comparison over to the estimation of phase velocities and enabled a fusion of the results while enlarging the bandwidth. Here, we employed a semi-classical analysis point of view inserting a small parameter in the variations in medium properties in the surface normal coordinate.
    • In highly discontinuous and random planarly layered media, we developed an understanding of tunneling, hidden equivalent medium averaging, and microlocal attenuation which cannot originate from standard visco-elastic models [1,2,5,8].

On planetary scale, that is, a ball, we studied a random radially layered medium including random angular fluctuations on distinct scales as well as deterministic discontinuities. We obtained a stochastic analysis in terms of modes and partly coherent waves excited by an interior source. Through cross correlations of the trace (field at the surface) we obtained a (virtual) data set which can be utilized for the imaging of the deterministic discontinuities.

  • Directional wavefield decomposition and propagation, sufficiently smooth coefficients. In the case of sufficient medium regularity one can interchange time locally with a (curvilinear) spatial coordinate (‘pseudodepth’) to describe the propagation of waves in terms of an evolution or one-way wave equation. In certain non-trivial wavespeed profiles this can be accomplised exactly, using spectral theory (modes) and special functions [23]. Under certain conditions, we can incorporate and analyze the role of ‘evanescent’ wave constituents [24]. In general, this technique is desgned for relatively high frequency wave propagation. In appropriate coordinates, we derived a novel paraxial equation and obtained estimates for its solution in terms of an energy norm converging to the full wave solution as the frequency becomes large. This is accomplished, again, via parabolic scaling. The ‘up-down’ coupling can be incorporated through a generalized Bremmer series [11]. These consistent paraxial equations enable the development of fast algorithms [15], and are of current interest in, for example, long-range propagation of and imaging with teleseismic body waves.
  • Anisotropy. We studied medium anisotropy from the point of view of polarized waves and a corresponding local slowness surface characterization (for example, by introducing anellipticity) and expansions [20], and the point of view of classification, characterization and construction of waveforms via Picard-Fuchs equations [7].

The behaviors and phenomena described above provide fundamental insight in developing multi-scale methods to estimate local continuity, anisotropy, regularity etc. of medium properties and substructures from broadband seismic observations, and can lead to a deeper understanding of the coupling of geological and geodynamical processes, (complex) phase boundaries etc. and the interaction (and, hence, probing) with seismic waves. In this context, we obtained some basic results in the case of viscous flow in the mantle wedge induced by a subducting lithospheric plate; we approximated the subduction of the Philippine Sea plate beneath Eurasia, explored the effects of different rheological and thermal constraints, and modelled the evolution of the finite strain ellipses [61].

Algorithms

Following the multi-scale physics of waves a family of new computational algorithms has emerged:

  • We designed and implemented a fast discrete almost symmetric wavepacket transform pair in three dimensions [80]. This transform pair can also be utilized for wave-based signal processing. Waves (as well as coefficient models containing interfaces) can be compressed with this transform.
  • Following the dyadic parabolic decomposition of phase space underlying the construction of wave packets, we developed a fast algorithm (with error estimates) for a class of so-called Fourier integral operators comprising propagators, asymptotic solution operators of the wave equations and (inverse) scattering transforms [92]. This algorithm is based on accurate low-rank separated representations, prolate spheroidal wave functions, and ray-geometric computations. (From a physics of waves point of view, in the case of one-way wave propagators, these representations are related to generalizations of phase screens [34,35].) In the application to the half-wave equation, we obtained a theoretical, multi-scale description of the Gouy phase which affects the extraction of travel time in long-range propagation. We designed a special procedure to incorporate the (numerical) generation of general caustics[104].
  • We estimated the regularity of the kernel of the Volterra operator in the above mentioned construction of weak, full-wave solutions in media of limited smoothness. This regularity determined the optimal choice of quadrature to discretize the integral. Moreover we estimated the approximation of solutions by finite sets or clusters of wave packets to complete the discretization of the Volterra equation and obtain a finite matrix representation [95].

We constructed frames of wave packets produced by parabolic dilation, rotation and translation of (a finite sum of) Gaussians and obtained asymptotics on the analogue of Daubechies’ frame criterion [110]. We also obtained a frame of wave-atom-like packets generated from pure Gaussians; these can be identified as multi-scale Gaussian beams. We developed an associated parametrix construction and algorithm.

Anisotropic surface-wave tomography – Tibetan plateau

Free oscillations, surface waves and spectral theory

I revisited the construction and analysis of the relevant system of elasto-gravitational equations describing the free oscillations of the earth. The emphasis here is to reconcile (non-smooth, typically Lebesgue measurable) variations in material properties and the geodynamical processes in Earth’s interior with this construction. We adapted the weak formulation, in particular, with a view to solid-fluid boundaries, introduced appropriate function spaces and proved well-posedness using the theory of Hille-Yoshida and applying the Lumer-Phillip theorem. We also obtained energy estimates.

Semi-classical analysis of surface waves. We carried out an analysis of surface waves in a medium which is stratified near its boundary at some scale comparable to the wave length. Such a medium can also be thought of as a surficial layer (which can be thick) overlying a half space. We described how the (dispersive) propagation of such waves is governed by effective Hamiltonians on the boundary, and showed that the system is displayed by a space-adiabatic behavior. We obtained pseudodifferential surface wave equations and, via parametrix constructions, expansions of surface-wave amplitudes. These play a role, for example, in surface-wave tomography.

(Our analysis applies to the study of surface waves in Earth’s ‘near’ surface in the scaling regime mentioned above. The existence of such waves, that is, propagating wave solutions which decay exponentially away from the boundary of a homogeneous (elastic) half-space was first noted by Rayleigh. Rayleigh and (‘transverse’) Love waves can be identified with Earth’s free oscillation triples nSl and nTl with n « l/4. Love was the first to argue that surface-wave dispersion is responsible for the oscillatory character of the main shock of an earthquake tremor, following the ‘primary’ and ‘secondary’ arrivals.)

The Hamiltonians are pseudodifferential operators with non-trivial dispersion relations, that is, principal symbols which are non-homogeneous. We derived this with techniques from semi-classical analysis. Each Hamiltonian is identified with an eigenvalue of a locally one-dimensional Schrödinger problem on the one hand, and generates a flow identified with surface-wave bicharacteristics on the other hand. We showed that the eigenvalues exist under certain assumptions for the local profiles of coefficients in the boundary normal direction. The surface waves correspond with the discrete spectrum of the Schro¨dinger operator. The semi-classical asymptotics of the discrete spectrum provides a way to study the uniqueness and stability of the corresponding locally one-dimensional spectral inverse problems.

We used our results from semi-classical analysis to further develop a method, related to so-called eikonal tomography, for the recovery of anisotropic elastic parameters [preprint 9]. We first obtained a reconstruction of surface-wave bicharacteristics in the presence of azimuthal anisotropy enabled by appropriate geometric regularization from ambient noise generated dense array data. We then addressed the one-dimensional spectral inverse problems in the piecewise constant coefficients case.  We applied our method to Rayleigh wave signals extracted from cross-correlations of ambient noise recorded at a seismograph network to constrain heterogeneity and anisotropy in the crust of SE Tibet.

Carbonite Rock and Fracture

Fast algorithms, massively parallel computing

Direct structured solvers and time-harmonic waves

In the invere boundary value problem, to probe the Dirichlet-to-Neumann map, one needs to solve the system of equations describing time-harmonic waves for a large number of right-hand sides. We use (compact) FD (with mass lumping) and CG-FEM methods for the discretization. A typical region comprises a three-dimensional structured or unstructured tetrahedral mesh consisting of about 109 elements.

We introduced massively parallel structured direct solvers for this purpose using Hierarchically Semi-Separable matrices, transforming the way time-harmonic waves are computed approaching linear complexity. We exploit the fact that data are acquired with limited accuracy and use techniques to control finite accuracy computation. Our methods involve multiple layers of hierarchical tree structures and advanced graph and sparse matrix techniques.

  • Hierarchically Semi-Separable (HSS) matrices, structured multifrontal solvers. HSS matrix algorithms have been emerging in the construction of the superfast direct solvers for both dense and sparse linear systems. We developed a set of new parallel algorithms for the key HSS operations that are used for solving such large systems. These include the parallel rank-revealing QR factorization, the HSS constructions with hierarchical compression, the ULV HSS factorization, and the HSS solutions [105]. We coupled these to the multifrontal method, using a matrix re-ordering following nested dissection based domain decomposition, and designed and implemented massively parallel algorithms [89,94,100]. The rank patterns of the off-diagonal blocks are precisely such that we obtain a storage requirement which is nearly O(N) (N signifies the size of the matrix) and a computational complexity of the factorization which falls just below O(N 4/3) in dimension three.
  • The method in [105] follows a geometric approach with fits discretizations based on rectangular grids or structured meshes. We have now completed the development of an algebraic approach, which is amenable to the use of unstructured tetrahedral meshes (arising in the conditional Lipschitz stability estimates for the inverse boundary value problem and multi-level schemes).

    The direct structured solvers also have applications in certain problems of computational geodynamics; the unstructured tetrahedral meshes – with multi-scale refinement and deformation – provide an interface between these and (seismic) inverse problems.

    We developed distributed memory parallel randomized HSS algorithms, which avoid the storage of any dense matrix. The parallelism in the HSS tree structure is exploited by overlapping the computations on each separate subtree using different subgroups of processors. Within each subgroup, a 2D process grid is created for the block-cyclic matrix distribution of each HSS generator. If r is the rank bound of the off-diagonal blocks and n is the matrix size, the parallel run time for randomized HSS algorithms is essentially O(r) with O(rn) processors.

  • Nonlinear eigenvalue problems and normal modes. We anticipate that the structured matrix approach will lead to a superfast algorithm to solve nonlinear eigenvalue problems (three-dimensional anisotropic, rotating, self-gravitating Earth models) in prescribed frequency intervals. This would make it feasible to compute directly normal modes in a three-dimensional heterogeneous earth.

The structured multifrontal factorization techniques apply to numerically solving certain integral equations, such as the Lippman-Schwinger equation and, via a Neumann series, the integral equation driving (time-harmonic) surface-related multiple elimination.

Time-domain solvers

In the time domain, we developed a novel Discontinuous Galerkin (DG) method for acousto-elastic waves allowing general anisotropy. In my program, the use of the DG method is primarily motivated by the introduction of unstructured tetrahedral meshes (and polyhedral interfaces) through the relevant hyperbolic inverse problems subject to conditions of well-posedness. This method yields a block-diagonal mass matrix.

We started from a coupled first-order elastic strain-velocity and acoustic velocity-pressure formulation. We introduced a novel penalty flux and avoided diagonalization. Our penalty flux accommodates fluid-solid boundaries. We obtained a convergence result. In the implementation, we made use of a fast implicit-explicit time integration approach so as to handle a local multi-scale refinement of the mesh, again, arising in iterative schemes for hyperbolic inverse problems. (We note that the structured multifrontal factorization using HSS matrices directly applies to the fast computation of the propagator, that is, a complex exponential map, associated the half scalar-wave equation.)

We have designed a parallel computational platform from the laboratory (microstructures such as fracture systems) and exploration (geological structures, salt tectonics) to global (planetary) scale.

Model structures, flexible updating and inverse problems

It seems to be natural to integrate (the modelling of) structural geological and geodynamical processes with domain partitioning (meshing) and associated sparse multi-scale model representations for which the mentioned Lipschitz stability estimates hold. To constrain the deformation during updating in iterative methods, we exploit connections with deforming meshes (and mesh refinement). Techniques based on topological completion formulas obtained by a homology criterion are applicable.

Discontinuous Galerkin Method: Rough surface of a carbonate rock sample.

Fast algorithms, massively parallel computing

Direct structured solvers and time-harmonic waves

In the invere boundary value problem, to probe the Dirichlet-to-Neumann map, one needs to solve the system of equations describing time-harmonic waves for a large number of right-hand sides. We use (compact) FD (with mass lumping) and CG-FEM methods for the discretization. A typical region comprises a three-dimensional structured or unstructured tetrahedral mesh consisting of about 109 elements.

We introduced massively parallel structured direct solvers for this purpose using Hierarchically Semi-Separable matrices, transforming the way time-harmonic waves are computed approaching linear complexity. We exploit the fact that data are acquired with limited accuracy and use techniques to control finite accuracy computation. Our methods involve multiple layers of hierarchical tree structures and advanced graph and sparse matrix techniques.

  • Hierarchically Semi-Separable (HSS) matrices, structured multifrontal solvers. HSS matrix algorithms have been emerging in the construction of the superfast direct solvers for both dense and sparse linear systems. We developed a set of new parallel algorithms for the key HSS operations that are used for solving such large systems. These include the parallel rank-revealing QR factorization, the HSS constructions with hierarchical compression, the ULV HSS factorization, and the HSS solutions [105]. We coupled these to the multifrontal method, using a matrix re-ordering following nested dissection based domain decomposition, and designed and implemented massively parallel algorithms [89,94,100]. The rank patterns of the off-diagonal blocks are precisely such that we obtain a storage requirement which is nearly O(N) (N signifies the size of the matrix) and a computational complexity of the factorization which falls just below O(N 4/3) in dimension three.
  • The method in [105] follows a geometric approach with fits discretizations based on rectangular grids or structured meshes. We have now completed the development of an algebraic approach, which is amenable to the use of unstructured tetrahedral meshes (arising in the conditional Lipschitz stability estimates for the inverse boundary value problem and multi-level schemes).

    The direct structured solvers also have applications in certain problems of computational geodynamics; the unstructured tetrahedral meshes – with multi-scale refinement and deformation – provide an interface between these and (seismic) inverse problems.

    We developed distributed memory parallel randomized HSS algorithms, which avoid the storage of any dense matrix. The parallelism in the HSS tree structure is exploited by overlapping the computations on each separate subtree using different subgroups of processors. Within each subgroup, a 2D process grid is created for the block-cyclic matrix distribution of each HSS generator. If r is the rank bound of the off-diagonal blocks and n is the matrix size, the parallel run time for randomized HSS algorithms is essentially O(r) with O(rn) processors.

  • Nonlinear eigenvalue problems and normal modes. We anticipate that the structured matrix approach will lead to a superfast algorithm to solve nonlinear eigenvalue problems (three-dimensional anisotropic, rotating, self-gravitating Earth models) in prescribed frequency intervals. This would make it feasible to compute directly normal modes in a three-dimensional heterogeneous earth.

The structured multifrontal factorization techniques apply to numerically solving certain integral equations, such as the Lippman-Schwinger equation and, via a Neumann series, the integral equation driving (time-harmonic) surface-related multiple elimination.

Time-domain solvers

In the time domain, we developed a novel Discontinuous Galerkin (DG) method for acousto-elastic waves allowing general anisotropy. In my program, the use of the DG method is primarily motivated by the introduction of unstructured tetrahedral meshes (and polyhedral interfaces) through the relevant hyperbolic inverse problems subject to conditions of well-posedness. This method yields a block-diagonal mass matrix.

We started from a coupled first-order elastic strain-velocity and acoustic velocity-pressure formulation. We introduced a novel penalty flux and avoided diagonalization. Our penalty flux accommodates fluid-solid boundaries. We obtained a convergence result. In the implementation, we made use of a fast implicit-explicit time integration approach so as to handle a local multi-scale refinement of the mesh, again, arising in iterative schemes for hyperbolic inverse problems. (We note that the structured multifrontal factorization using HSS matrices directly applies to the fast computation of the propagator, that is, a complex exponential map, associated the half scalar-wave equation.)

We have designed a parallel computational platform from the laboratory (microstructures such as fracture systems) and exploration (geological structures, salt tectonics) to global (planetary) scale.

Model structures, flexible updating and inverse problems

It seems to be natural to integrate (the modelling of) structural geological and geodynamical processes with domain partitioning (meshing) and associated sparse multi-scale model representations for which the mentioned Lipschitz stability estimates hold. To constrain the deformation during updating in iterative methods, we exploit connections with deforming meshes (and mesh refinement). Techniques based on topological completion formulas obtained by a homology criterion are applicable.