Multiscale Physics of Waves

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

• Multi-scale ray-wave duality. The description of waves up to C^1,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. with increasing smoothness of the medium. We obtain a concentration of wave packets result; the degree of concentration reveals the smoothness of the medium ¹. 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 Lamé parameters for which polarized waves can be distinguished, which appears to be C^1,1 .

• 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 X-ray transform, by taking Fréchet-type derivatives of our solution construction, and resolves the seeming insensitivity of the transform along the unperturbed rays via a multi-scale analysis. In the case of smooth wavespeeds, and the limit of infinite bandwdith, the geodesic X-ray transform is recovered .

• Scattering series, interfaces. Reducing the regularity of the medium further, leads us to the introduction of (general) interfaces. We succeeded to track the order of scattering through the regularity of the wave constituents in the case of C^1,ε (ε > 0) coefficients. 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. 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 .

  • 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  (including directionality ). 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 cannnot originate from standard visco-elastic models 
• 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. Under certain conditions, we can incorporate and analyze the role of ‘evanescent’ wave constituents . 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 . These consistent paraxial equations enable the development of fast algorithms, 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, and the point of view of classification, characterization and construction of waveforms via Picard-Fuchs equations . 

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 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.

Algorithms

• We designed and implemented a fast discrete almost symmetric wavepacket transform pair in three dimensions . 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. 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.) 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 

• 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 .

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 avail- able), 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 asso- ciated with inverse scattering is still incomplete with so-called I^plclasses appearing in the analysis. To ensure that false discontinuities cannot be generated, we invoked the Bolker condition , which may be violated in particular if the acquisition looses a dimenion. 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 . 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-called curvilinear double-square-root (DSR) condition . 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

• Annihilators, wave-equation reflection tomography. We introduced annihilators (fundamentally extending the notion of differential semblance to allow for the formation of caustics in the unknown background) of the data. Annihilators characterize the range of the single scattering operator 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. We note that a Sobolev gradient appears naturally in the minimization of the annihilators, which is critical in the case of non-smooth (unknown) reflectors (in the contrast). We 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

• Generalized Radon transform (GRT). We have constructed a GRT for inverse scattering incorporating the asymp- totic inverse (parametrix) of the normal operator, in anisotropic media, assuming that the equations for elastic waves are of principal type and carried out an associated resolution analysis. 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. Moreover, upon adapting the order of the Fourier integral operator representing the GRT, we obtained stability and a method for statistical inference of singularities. We developed an approach to precise partial reconstruction from partial data with the GRT and illumination correction using our frame of wavepackets.
  • We applied the GRT to image the D″ 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. We have begun the integration with SKKS coda to enlarge the region of illumination on the one hand and enable the further estimation of anisotropy on the other hand.
  • We also applied the GRT to image the hot spot underneath Hawaii using underside reflections, that is, SS precursors. 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. (We are currently using other phases to extend our study of temperature anomalies associated with a deep plume.)  

• Reverse-time-migration (RTM)-based inverse scattering. We constructed a new inverse scattering transform based on RTM , also in the anisotropic case, 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. In this algorithm it is straightforward to extract information about the local incidence angle and thus generate angle-dependent images. As an alternative, for teleseismic applications, we can invoke our paraxial equations coupled to a decomposition into wave packets, which can be viewed as a replacement of beam migration.
  • We can remove the knowledge of the source in the case of converted waves. In 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 ensuring statistical stability.  

• 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) equa- tion. 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. We also constructed explicitly the corresponding evolution equation (and thus an extended exploding re- flector model) which provided a new insight in the geometry of extended imaging in the presence of caustics via the introduction of extended isochron rays.

(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, 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.