Search

We have developed a Liouville partial-differential-equation (PDE)-based method for computing complex-valued eikonals in real phase space in the multivalued sense in attenuating media with frequency-independent qualify factors, where the new method computes the real and imaginary parts of the complex-valued eikonal in two steps by solving Liouville equations in real phase space. Because the earth is composed of attenuating materials, seismic waves usually attenuate so that seismic data processing calls for properly treating the resulting energy losses and phase distortions of wave propagation. In the regime of high-frequency asymptotics, the complex-valued eikonal is one essential ingredient for describing wave propagation in attenuating media because this unique quantity summarizes two wave properties into one function: Its real part describes the wave kinematics and its imaginary part captures the effects of phase dispersion and amplitude attenuation. Because some popular ordinary-differential-equation (ODE)-based ray-tracing methods for computing complex-valued eikonals in real space distribute the eikonal function irregularly in real space, we are motivated to develop PDE-based Eulerian methods for computing such complex-valued eikonals in real space on regular meshes. Therefore, we solved novel paraxial Liouville PDEs in real phase space so that we can compute the real and imaginary parts of the complex-valued eikonal in the multivalued sense on regular meshes. We call the resulting method the Liouville PDE method for complex-valued multivalued eikonals in attenuating media; moreover, this new method provides a unified framework for Eulerianizing several popular approximate real-space ray-tracing methods for complex-valued eikonals, such as viscoacoustic ray tracing, real viscoelastic ray tracing, and real elastic ray tracing. In addition, we also provide Liouville PDE formulations for computing multivalued ray amplitudes in a weakly viscoacoustic medium. Numerical examples, including a synthetic gas-cloud model, demonstrate that our methods yield highly accurate complex-valued eikonals in the multivalued sense.

Nondestructive testing (NDT) methods are widely used in the rail industry to detect and characterize rolling contact fatigue (RCF) defects in railroads, which is very important for railway inspection and maintenance to prevent catastrophic accidents. Existing NDT methods, e.g., ultrasonic testing (UT), magnetic flux leakage (MFL), and eddy current testing (ECT) have been successfully applied in the rail industry, and the state-of-the-art UT method reported recently achieved a high speed of 40 mi/h with a probability of detection (POD) over 80% under 30% false alarms. However, NDT methods still suffer from a bottleneck in that a higher inspection speed causes lower detection sensitivity due to their physical limits, such as negative velocity effect and long sensing time. One of the leading challenges to the rail NDT community is to develop a high-speed high-sensitivity (HSHS) capability that can provide an improved POD of rail defects in high-speed inspection scenarios over 60 mi/h. This article proposes a horizontal U-shaped magnets-based motion-induced eddy current array (MIECA) method to detect rail surface defects with the HSHS capability. The MIECA method deploys a three-axis magnetic sensor array along the rail transverse direction at the middle of the magnets to measure the MIECA signals, which utilizes the wake effect of the diffused motion-induced eddy current (MIEC) caused by the relative high-speed motion between the magnets and the rail track. Finite-element method (FEM) simulations with a wide speed range from 0 to 62.5 mi/h are carried out to investigate the relationships between the MIEC generation, diffusion and magnitude, and the three-axis MIECA signals. The simulation results show that the higher the speed, the greater the magnitude of diffused MIEC, and the greater the peak-to-peak values of three-axis MIECA signals caused by rail surface defects, which shows a great promise for detecting rail surface defects in high-speed scenarios and is superior inspection relative to existing NDT methods in terms of inspection speed, detection sensitivity, and defect characterization capability.

First-arrival traveltime tomography is an essential method for obtaining near-surface velocity models. The adjoint-state first-arrival traveltime tomography is appealing due to its straightforward implementation, low computational cost, and low memory consumption. Because solving the point-source isotropic eikonal equation by either ray tracers or eikonal solvers intrinsically corresponds to emanating discrete rays from the source point, the resulting traveltime gradient is singular at the source point, and we denote such a singular pattern the imprint of ray-illumination. Because the adjoint-state equation propagates traveltime residuals back to the source point according to the negative traveltime gradient, the resulting adjoint state will inherit such an imprint of ray-illumination, leading to singular gradient-descent directions when updating the velocity model in the adjoint-state traveltime tomography. To mitigate this imprint, we solve the adjoint-state equation twice but with different boundary conditions: one being taken to be regular data residuals and the other taken to be ones uniformly, so that we are able to use the latter adjoint state to normalize the regular adjoint state and we further use the normalized quantity to serve as the gradient direction to update the velocity model; we call this process ray-illumination compensation. To overcome the issue of limited aperture, we have developed a spatially varying regularization method to stabilize the new gradient direction. A synthetic example demonstrates that our method is able to mitigate the imprint of ray-illumination, remove the footprint effect near source points, and provide uniform velocity updates along raypaths. A complex example extracted from the Marmousi2 model and a migration example illustrate that the new method accurately recovers the velocity model and that an offset-dependent inversion strategy can further improve the quality of recovered velocity models.

Seismic waves in earth media usually undergo attenuation, causing energy losses and phase distortions. In the regime of high-frequency asymptotics, a complex-valued eikonal is an essential ingredient for describing wave propagation in attenuating media, where the real and imaginary parts of the eikonal function capture dispersion effects and amplitude attenuation of seismic waves, respectively. Conventionally, such a complex-valued eikonal is mainly computed either by tracing rays exactly in complex space or by tracing rays approximately in real space so that the resulting eikonal is distributed irregularly in real space. However, seismic data processing methods, such as prestack depth migration and tomography, usually require uniformly distributed complex-valued eikonals. Therefore, we have developed a unified framework to Eulerianize several popular approximate real-space ray-tracing methods for complex-valued eikonals so that the real and imaginary parts of the eikonal function satisfy the classic real-space eikonal equation and a novel real-space advection equation, respectively, and we dub the resulting method the Eulerian partial-differential-equation method. We further develop highly efficient high-order methods to solve these two equations by using the factorization idea and the Lax-Friedrichs weighted essentially nonoscillatory schemes. Numerical examples demonstrate that our method yields highly accurate complex-valued eikonals, analogous to those from ray-tracing methods. Our methods can be useful for migration and tomography in attenuating media.

We propose a data and knowledge driven approach for SPECT by combining a classical iterative algorithm of SPECT with a convolutional neural network. The classical iterative algorithm, such as ART and ML-EM, is employed to provide the model knowledge of SPECT. A modified U-net is then connected to exploit further features of reconstructed images and data sinograms of SPECT. We provide mathematical formulations for the architecture of the proposed networks. The networks are trained by supervised learning using the technique of mini-batch optimization. We apply the trained networks to the problems of simulated lung perfusion imaging and simulated myocardial perfusion imaging, and numerical results demonstrate their effectiveness of reconstructing source images from noisy data measurements.

Starting from Hadamard's method, we develop Babich's ansatz for the frequency-domain point-source elastic wave equations in an inhomogeneous medium in the high-frequency regime. First, we develop a...

We propose an efficient PDE based approach for solving the first-arrival traveltime tomography problem in irregular domains. We consider the mismatch functional based on the $L_1$ or $L_2$ fidelity...

We propose a new fast multiscale Gaussian beam method to solve the three-dimensional elastic wave equation in a bounded domain in the high-frequency regime. We develop a novel multiscale transform ...

We propose a Newton-type Gauss–Seidel Lax–Friedrichs sweeping method to solve the generalized eikonal equation arising from wave propagation in a moving fluid. The Lax–Friedrichs numerical Hamiltonian is used in discretization of the generalized eikonal equation. Different from traditional Lax–Friedrichs sweeping algorithms, we design a novel approach with a line-wise sweeping strategy. In the local solver, the values of traveltime on an entire line are updated simultaneously by Newton’s method. The global solution is then obtained by Gauss–Seidel iterations with line-wise sweepings. We first develop the Newton-based first-order scheme, and on top of that we further develop high-order schemes by applying weighted essentially non-oscillatory (WENO) approximations to derivatives. Extensive 2-D and 3-D numerical examples demonstrate the efficiency and accuracy of the new algorithm. The combination of Newton’s method and Gauss–Seidel iterations improves upon the convergence speed of the original Lax–Friedrichs sweeping algorithm. In addition, the Newton-type sweeping method manipulates data in a vectorized manner so that it can be efficiently implemented in modern programming languages that feature array programming, and the resulting advantages are extremely significant for large-scale 3-D computations.

The need to improve the depth resolution of the magnetic susceptibility model recovered from surface magnetic data is a well-known challenge, and it becomes increasingly important as exploration moves to regions under cover and at great depths. Incorporating borehole magnetic data can be an effective means to achieve increased model resolution at depth. The recently developed level-set method for magnetic inversion provides a novel means for constructing the geometric shape of causative bodies and opens a new avenue for the joint inversion of surface and borehole magnetic data for the purpose of achieving improved depth resolution. We have developed an extension of the algorithm to the joint inversion and find that the level-set algorithm can resolve the configuration and spatial separation of complex magnetic sources using the information in the magnetic data from sparse boreholes. We further examine the use of borehole intersection information in estimating the crucially important susceptibility values within the context of level-set inversion and find that the susceptibility value can also be used as a probing parameter to assess the uncertainty in the spatial extent of the causative bodies. We determine that the modified level-set inversion leads to an effective means to image and delineate magnetic causative bodies with complex structure by combining the information from surface magnetic data, borehole magnetic data, and sparse drillhole intersection data.

We develop efficient algorithms for computing the multidimensional fractional operator $( -\Delta_{x})^{\frac{\alpha}{2}}$ in the form of hypersingular integral in the entire space, where the opera...

We present a butterfly-compressed representation of the Hadamard-Babich (HB) ansatz for the Green's function of the high-frequency Helmholtz equation in smooth inhomogeneous media. For a computational domain discretized with $N_v$ discretization cells, the proposed algorithm first solves and tabulates the phase and HB coefficients via eikonal and transport equations with observation points and point sources located at the Chebyshev nodes using a set of much coarser computation grids, and then butterfly compresses the resulting HB interactions from all $N_v$ cell centers to each other. The overall CPU time and memory requirement scale as $O(N_v\log^2N_v)$ for any bounded 2D domains with arbitrary excitation sources. A direct extension of this scheme to bounded 3D domains yields an $O(N_v^{4/3})$ CPU complexity, which can be further reduced to quasi-linear complexities with proposed remedies. The scheme can also efficiently handle scattering problems involving inclusions in inhomogeneous media. Although the current construction of our HB integrator does not accommodate caustics, the resulting HB integrator itself can be applied to certain sources, such as concave-shaped sources, to produce caustic effects. Compared to finite-difference frequency-domain (FDFD) methods, the proposed HB integrator is free of numerical dispersion and requires fewer discretization points per wavelength. As a result, it can solve wave-propagation problems well beyond the capability of existing solvers. Remarkably, the proposed scheme can accurately model wave propagation in 2D domains with 640 wavelengths per direction and in 3D domains with 54 wavelengths per direction on a state-the-art supercomputer at Lawrence Berkeley National Laboratory.

We present efficient numerical methods for solving a class of nonlinear Schrodinger equations involving a nonlocal potential. Such a nonlocal potential is governed by Gaussian convolution of the intensity modeling nonlocal mutual interactions among particles. The method extends the Fast Huygens Sweeping Method (FHSM) that we developed in Leung et al. (Methods Appl Anal 21(1):31–66, 2014) for the linear Schrodinger equation in the semi-classical regime to the nonlinear case with nonlocal potentials. To apply the methodology of FHSM effectively, we propose two schemes by using the Lie’s and the Strang’s operator splitting, respectively, so that one can handle the nonlinear nonlocal interaction term using the fast Fourier transform. The resulting algorithm can then enjoy the same computational complexity as in the linear case. Numerical examples demonstrate that the two operator splitting schemes achieve the expected first-order and second-order accuracy, respectively. We will also give one-, two- and three-dimensional examples to demonstrate the efficiency of the proposed algorithm.

We develop a deep learning approach to extract ray directions at discrete locations by analyzing highly oscillatory wave fields. A deep neural network is trained on a set of local plane-wave fields to predict ray directions at discrete locations. The resulting deep neural network is then applied to a reduced-frequency Helmholtz solution to extract the directions, which are further incorporated into a ray-based interior-penalty discontinuous Galerkin (IPDG) method to solve the Helmholtz equations at higher frequencies. In this way, we observe no apparent pollution effects in the resulting Helmholtz solutions in inhomogeneous media. Our 2D and 3D numerical results show that the proposed scheme is very efficient and yields highly accurate solutions., 30 pages

In the present article we extend to the three-dimensional elliptic Monge–Ampere equation the method discussed in Glowinski et al. (J Sci Comput 79:1–47, 2019) for the numerical solution of its two-dimensional variant. As in Glowinski et al. (2019) we take advantage of an equivalent divergence formulation of the Monge–Ampere equation, involving the cofactor matrix of the Hessian of the solution. We associate with the above divergence formulation an initial value problem, well suited to time discretization by operator splitting and space approximation by low order mixed finite element methods. An important ingredient of our methodology is forcing the positive semi-definiteness of the approximate Hessian by a hard thresholding eigenvalue projection. The resulting method is robust and easy to implement. It can handle problems with smooth and non-smooth solutions on domains with curved boundary. Using piecewise affine approximations for the solution and its six second-order derivatives, one can achieve second-order convergence rates for problems with smooth solutions.

We have developed a novel Kantorovich-Rubinstein (KR) norm-based misfit function to measure the mismatch between gravity-gradient data for the inverse gradiometry problem. Under the assumption that an anomalous mass body has an unknown compact support with a prescribed constant value of density contrast, we implicitly parameterize the unknown mass body by a level-set function. Because the geometry of an underlying anomalous mass body may experience various changes during inversion in terms of level-set evolution, the classic least-squares ([Formula: see text]-norm-based) and the [Formula: see text]-norm-based misfit functions for governing the level-set evolution may potentially induce local minima if an initial guess of the level-set function is far from that of the target model. The KR norm from the optimal transport theory computes the data misfit by comparing the modeled data and the measured data in a global manner, leading to better resolution of the differences between the inverted model and the target model. Combining the KR norm with the level-set method yields a new effective methodology that is not only able to mitigate local minima but is also robust against random noise for the inverse gradiometry problem. Numerical experiments further demonstrate that the new KR norm-based misfit function is able to recover deep dipping flanks of SEG/EAGE salt models even at extremely low signal-to-noise ratios. The new methodology can be readily applied to gravity and magnetic data as well.

We develop a novel numerical microlocal analysis (NMLA) method using fast Gaussian wave packet transforms. Our new NMLA method extracts ray directions at discrete locations by analyzing highly osci...

We propose an $L^2$ regularized Kantorovich--Rubinstein (KR) metric as a novel misfit function for the level-set based inverse gravity problems without assuming mass conservation. We study the well...

We discuss in this article a novel method for the numerical solution of the two-dimensional elliptic Monge–Ampere equation. Our methodology relies on the combination of a time-discretization by operator-splitting with a mixed finite element based space approximation where one employs the same finite-dimensional spaces to approximate the unknown function and its three second order derivatives. A key ingredient of our approach is the reformulation of the Monge–Ampere equation as a nonlinear elliptic equation in divergence form, involving the cofactor matrix of the Hessian of the unknown function. With the above elliptic equation we associate an initial value problem that we discretize by operator-splitting. To enforce the pointwise positivity of the approximate Hessian we employ a hard thresholding based projection method. As shown by our numerical experiments, the resulting methodology is robust and can handle a large variety of triangulations ranging from uniform on rectangles to unstructured on domains with curved boundaries. For those cases where the solution is smooth and isotropic enough, we suggest also a two-stage method to improve the computational efficiency, the second stage being reminiscent of a Newton-like method. The methodology discussed in this article is able to handle domains with curved boundaries and unstructured meshes, using piecewise affine continuous approximations, while preserving optimal, or nearly optimal, convergence orders for the approximation error.

For scattering problems of time-harmonic waves, the boundary integral equation (BIE) methods are highly competitive, since they are formulated on lower-dimension boundaries or interfaces, and can automatically satisfy outgoing radiation conditions. For scattering problems in a layered medium, standard BIE methods based on the Green's function of the background medium must evaluate the expensive Sommefeld integrals. Alternative BIE methods based on the free-space Green's function give rise to integral equations on unbounded interfaces which are not easy to truncate, since the wave fields on these interfaces decay very slowly. We develop a BIE method based on the perfectly matched layer (PML) technique. The PMLs are widely used to suppress outgoing waves in numerical methods that directly discretize the physical space. Our PML-based BIE method uses the Green's function of the PML-transformed free space to define the boundary integral operators. The method is efficient, since the Green's function of the PML-transformed free space is easy to evaluate and the PMLs are very effective in truncating the unbounded interfaces. Numerical examples are presented to validate our method and demonstrate its accuracy., Comment: 37 pages, 14 figures, 1 table

Understanding effective Hamiltonians quantitatively is essential for the homogenization of Hamilton--Jacobi equations. We propose in this article a simple efficient operator-splitting method for computing effective Hamiltonians when the Hamiltonian is either convex or nonconvex in the gradient variable. To speed up our Lie scheme--based operator-splitting method, we further propose a cascadic initialization strategy so that the steady state of the underlying time-dependent Hamilton--Jacobi equation can be reached more rapidly. Extensive numerical examples demonstrate the efficiency and accuracy of the new algorithm.

Starting from Hadamard's method, we extend Babich's ansatz to the frequency-domain point-source (FDPS) Maxwell's equations in an inhomogeneous medium in the high-frequency regime. First, we develop...

Microwave-induced thermoacoustic imaging (TAI) is a hybrid imaging technique that combines electromagnetic radiation and ultrasonic waves to achieve high imaging contrast and submillimeter spatial resolution. These characteristics make TAI a good candidate to detect material anomalies that change the material electric properties without a noticeable variation in material density. Conventional pulsed TAI systems work by sending a single short pulse to the imaged target and then detecting the generated pressure signal; therefore, a very high peak power microwave pulse or data averaging is needed to produce images with a high signal-to-noise ratio (SNR). In this paper, we propose to enhance the SNR of pulsed TAI systems by using non-coherent pulse compression. In this approach, a predefined pulse coded signal is used to illuminate the imaged sample and the received pressure signal is cross correlated with a template that is related to the power profile of the excitation signal. The proposed approach can be easily deployed to pulsed TAI systems without the need for major system modifications to the RF source because it only requires a timing circuit to control the triggering time of the RF pulses. In this paper, we demonstrate experimentally that the proposed approach highly improves the SNR of TAI signals and images and can be used to reduce the acquisition time by lowering the number of data averaging or reduce the required peak power from RF sources.

We propose a numerically efficient algorithm for simulating the multi-color optical self-focusing phenomena in nematic liquid crystals. The propagation of the nematicon is modeled by a parabolic wave equation coupled with a nonlinear elliptic partial differential equation governing the angle between the crystal and the direction of propagation. Numerically, the paraxial parabolic wave equation is solved by a fast Huygens sweeping method, while the nonlinear elliptic PDE is handled by the alternating direction explicit (ADE) method. The overall algorithm is shown to be numerically efficient for computing high frequency beam propagations.

This paper is the first attempt to systematically study properties of the effective Hamiltonian $$\overline{H}$$ arising in the periodic homogenization of some coercive but nonconvex Hamilton–Jacobi equations. Firstly, we introduce a new and robust decomposition method to obtain min–max formulas for a class of nonconvex $$\overline{H}$$ . Secondly, we analytically and numerically investigate other related interesting phenomena, such as “quasi-convexification” and breakdown of symmetry, of $$\overline{H}$$ from other typical nonconvex Hamiltonians. Finally, in the appendix, we show that our new method and those a priori formulas from the periodic setting can be used to obtain stochastic homogenization for the same class of nonconvex Hamilton–Jacobi equations. Some conjectures and problems are also proposed.

We have developed a multiple level-set method for inverting magnetic data produced by weak induced magnetization only. The method is designed to deal with a specific class of 3D magnetic inverse problems in which the magnetic susceptibility is known and the objective of the inversion is to find the boundary or geometric shape of the causative bodies. We adopt the conceptual representation of the subsurface geologic structure by a set of magnetic bodies, each having a uniform magnetic susceptibility embedded in a nonmagnetic background. This representation enables us to reformulate the magnetic inverse problem into a domain inverse problem for those unknown domains defining the supports of the magnetic causative bodies. Because each body may take on a variety of shapes, and we may not know the number of bodies a priori either, we use multiple level-set functions to parameterize these domains so that the domain inverse problem can be further reduced to an optimization problem of multiple level-set functions. To efficiently compute gradients of the nonlinear functional arising from the multiple level-set formulation, we take advantage of the rapid decay of the magnetic kernels with distance to significantly speed up the matrix-vector multiplications in the minimization process. We apply the new method to the synthetic and field data sets and determine its effectiveness.

We develop new efficient algorithms for a class of inverse problems of gravimetry to recover an anomalous volume mass distribution (measure) in the sense that we design fast local level-set methods to simultaneously reconstruct both unknown domain and varying density of the anomalous measure from modulus of gravity force rather than from gravity force itself. The equivalent-source principle of gravitational potential forces us to consider only measures of the form \begin{document}$ \mu = f\,\chi_{D} $\end{document} , where \begin{document}$ f $\end{document} is a density function and \begin{document}$ D $\end{document} is a domain inside a closed set in \begin{document}$ \bf{R}^n $\end{document} . Accordingly, various constraints are imposed upon both the density function and the domain so that well-posedness theories can be developed for the corresponding inverse problems, such as the domain inverse problem, the density inverse problem, and the domain-density inverse problem. Starting from uniqueness theorems for the domain-density inverse problem, we derive a new gradient from the misfit functional to enforce the directional-independence constraint of the density function and we further introduce a new labeling function into the level-set method to enforce the geometrical constraint of the corresponding domain; consequently, we are able to recover simultaneously both unknown domain and varying density from given modulus of gravity force. Our fast level-set method is built upon localizing the level-set evolution around a narrow band near the zero level-set and upon accelerating numerical modeling by novel low-rank matrix multiplication. Numerical results demonstrate that uniqueness theorems are crucial for solving the inverse problem of gravimetry and will be impactful on gravity prospecting. To the best of our knowledge, our inversion algorithm is the first of such for the domain-density inverse problem since it is based upon the conditional well-posedness theory of the inverse problem.

We have developed a new level-set-based structural parameterization for joint inversion of gravity and traveltime data, so that density contrast and seismic slowness are simultaneously recovered in the inverse problem. Because density contrast and slowness are different model parameters of the same survey domain, we assume that they are similar in structure in terms of how each property changes and where the interface is located, so that we are able to use a level-set function to parameterize the common interface shared by these two model parameters. The level-set parameterization makes it easy to maintain the structural similarity between the two geophysical properties. The inversion of gravity and traveltime data is carried out by minimizing a joint data-fitting function with respect to density contrast, slowness, as well as the level-set function. An adjoint state method is used to compute the traveltime gradient efficiently. We have tested our algorithm on various synthetic examples, including a 2D ovoid model and the 2D SEG/EAGE salt model. The results show that the joint-inversion algorithm effectively improves recovery of subsurface features. To the best of our knowledge, this is the first time that the level-set method was used to structurally link density contrast and slowness distribution systematically.

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that the Helmholtz equation can be viewed as an evolution equation in one of the spatial directions. With such applications in mind, starting from Babich's expansion, we develop a new high-order asymptotic method, which we dub the fast Huygens sweeping method, for solving point-source Helmholtz equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of this method is that we develop a new Eulerian approach to compute the asymptotics, i.e. the traveltime function and amplitude coefficients that arise in Babich's expansion, yielding a locally valid solution, which is accurate close enough to the source. The second novelty is that we utilize the Huygens-Kirchhoff integral to integrate many locally valid wavefields to construct globally valid wavefields. This automatically treats caustics and yields uniformly accurate solutions both near the source and remote from it. The third novelty is that the butterfly algorithm is adapted to accelerate the Huygens-Kirchhoff summation, achieving nearly optimal complexity O ( N log ź N ) , where N is the number of mesh points; the complexity prefactor depends on the desired accuracy and is independent of the frequency. To reduce the storage of the resulting tables of asymptotics in Babich's expansion, we use the multivariable Chebyshev series expansion to compress each table by encoding the information into a small number of coefficients.The new method enjoys the following desired features. First, it precomputes the asymptotics in Babich's expansion, such as traveltime and amplitudes. Second, it takes care of caustics automatically. Third, it can compute the point-source Helmholtz solution for many different sources at many frequencies simultaneously. Fourth, for a specified number of points per wavelength, it can construct the wavefield in nearly optimal complexity in terms of the total number of mesh points, where the prefactor of the complexity only depends on the specified accuracy and is independent of frequency. Both two-dimensional and three-dimensional numerical experiments have been carried out to illustrate the performance, efficiency, and accuracy of the method.

We have developed a level-set method for the inverse gravimetry problem of imaging salt structures with density contrast reversal. Under such a circumstance, a part of the salt structure contributes two completely opposite anomalies that counteract with each other, making it unobservable to the gravity data. As a consequence, this amplifies the inherent nonuniqueness of the inverse gravimetry problem so that it is much more challenging to recover the whole salt structure from the gravity data. To alleviate the severe nonuniqueness, it is reasonable to assume that the density contrast between the salt structure and the surrounding sedimentary host depends upon the depth only and is known a priori. Consequently, the original inverse gravity problem reduces to a domain inverse problem, where the supporting domain of the salt body becomes the only unknown. We have used a level-set function to parametrize the boundary of the salt body so that we reformulated the domain inverse problem into a nonlinear optimization problem for the level-set function, which was further solved for by a gradient descent method. Both 2D and 3D experiments on the SEG/EAGE salt model were carried out to demonstrate the effectiveness and efficiency of the new method. The algorithm was able to recover dipping flanks of the salt model, and it only took 40 min in a 2.5 GHz CPU to invert for a 3D model of 97,000 unknowns.

Because angle-domain common-image gathers (ADCIGs) from reverse time migration (RTM) are capable of obtaining the correct illumination of a subsurface geologic structure, they provide more reliable information for velocity model building, amplitude-variation versus angle analysis, and attribute interpretation. The approaches for generating ADCIGs mainly consist of two types: (1) indirect approaches that convert extended image gathers into ADCIGs and (2) direct approaches that first obtain propagating angles of wavefronts and then map the imaging result to the angle domain. In practice, however, generation of ADCIGs usually incurs high computational cost, poor resolution, and other drawbacks. To generate efficient ADCIGs using RTM methods, we have introduced a novel approach to obtain polarization vectors — directions of particle motion — from the Cauchy wavefield (CWF) and an efficient localized plane-wave decomposition algorithm to implement the angle-domain imaging condition. The CWF is a wavefield constructed from the Cauchy condition of the wave equation at any given time, and it only contains negative frequencies of the original wavefield so that the polarization vector is obtained from the local CWF in the wavenumber domain. With polarization vectors at our disposal, we have further developed an efficient localized plane-wave decomposition algorithm to implement the angle-domain imaging condition. Numerical examples have indicated that the new approach is able to handle complex wave phenomenon and has advantages in illuminating subsurface structure.

Wave propagation in an isotropic acoustic medium occupied by a moving fluid is governed by an anisotropic eikonal equation. Since this anisotropic eikonal equation is associated with an inhomogeneous Hamiltonian, most of existing anisotropic eikonal solvers are either inapplicable or of unpredictable behavior in convergence. Realizing that this anisotropic eikonal equation is defined by a sum of two well-understood first-order differential operators, we propose novel operator-splitting based fast sweeping methods to solve this generalized eikonal equation. We develop various operator-splitting methods relying on the Peaceman--Rachford scheme, the Douglas--Rachford scheme, the $\theta$-scheme, and the regularized $\theta$-scheme. After applying the operator-splitting strategy, each splitting step corresponds to a much simpler Hamilton--Jacobi equation so that we can apply the Lax--Friedrichs sweeping method to solve these splitted equations efficiently and easily. Two- and three-dimensional examples demons...

We propose a novel Babich-like ansatz consisting of an infinite series of dyadic coefficients (three-by-three matrices) and spherical Hankel functions for solving point-source Maxwell's equations in an inhomogeneous medium so as to produce the so-called dyadic Green's function. Using properties of spherical Hankel functions, we derive governing equations for the unknown asymptotics of the ansatz including the traveltime function and dyadic coefficients. By proposing matching conditions at the point source, we rigorously derive asymptotic behaviors of these geometrical-optics ingredients near the source so that their initial data at the source point are well-defined. To verify the feasibility of the proposed ansatz, we truncate the ansatz to keep only the first two terms, and we further develop partial differential equation--based Eulerian approaches to compute the resulting asymptotic solutions. Since the system of governing equations for each dyadic coefficient is strongly coupled, we introduce auxiliary variables to transform these strongly coupled systems into decoupled scalar equations. Furthermore, we develop high-order Lax--Friedrichs weighted essentially nonoscillatory schemes for computing these auxiliary variables so that the Green's function can be constructed. Numerical examples demonstrate that our new ansatz yields a uniform asymptotic solution in the region of space containing a point source but no other caustics.

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that Maxwell's equations may be viewed as an evolution equation in one of the spatial directions. With such applications in mind, we propose a new Eulerian geometrical-optics method, dubbed the fast Huygens sweeping method, for computing Green's functions of Maxwell's equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of the fast Huygens sweeping method is that a new dyadic-tensor-type geometrical-optics ansatz is proposed for Green's functions which is able to utilize some unique features of Maxwell's equations. The second novelty is that the Huygens--Kirchhoff secondary source principle is used to integrate many locally valid asymptotic solutions to yield a globally valid asymptotic solution so that caustics associated with the usual geometrical-optics ansatz can be treated automatically. The third novelty is that a butterfly algorithm is adapted to carry out the matrix-vector products induced by the Huygens--Kirchhoff integration in $O(N\log N)$ operations, where $N$ is the total number of mesh points, and the proportionality constant depends on the desired accuracy and is independent of the frequency parameter. To reduce the storage of the resulting traveltime and amplitude tables, we compress each table into a linear combination of tensor-product based multivariate Chebyshev polynomials so that the information of each table is encoded into a small number of Chebyshev coefficients. The new method enjoys the following desired features: (1) it precomputes a set of local traveltime and amplitude tables; (2) it automatically takes care of caustics; (3) it constructs Green's functions of Maxwell's equations for arbitrary frequencies and for many point sources; (4) for a specified number of points per wavelength it constructs each Green's function in nearly optimal complexity $O(N\log N)$ in terms of the total number of mesh points $N$, where the prefactor of the complexity depends only on the specified accuracy and is independent of the frequency parameter. Three-dimensional numerical experiments are presented to demonstrate the performance and accuracy of the new method.

In this article, we report the results we obtained when investigating the numerical solution of some nonlinear eigenvalue problems for the Monge-Ampère operator v → det D2v. The methodology we employ relies on the following ingredients: (i) a divergence formulation of the eigenvalue problems under consideration. (ii) The time discretization by operator-splitting of an initial value problem (a kind of gradient flow) associated with each eigenvalue problem. (iii) A finite element approximation relying on spaces of continuous piecewise affine functions. To validate the above methodology, we applied it to the solution of problems with known exact solutions: The results we obtained suggest convergence to the exact solution when the space discretization step h → 0. We considered also test problems with no known exact solutions.

The usual geometrical-optics expansion of the solution for the Helmholtz equation of a point source in an inhomogeneous medium yields two equations: an eikonal equation for the traveltime function, and a transport equation for the amplitude function. However, two difficulties arise immediately: one is how to initialize the amplitude at the point source as the wavefield is singular there; the other is that in even-dimension spaces the usual geometrical-optics expansion does not yield a uniform asymptotic approximation close to the source. Babich (USSR Comput Math Math Phys 5(5):247---251, 1965) developed a Hankel-based asymptotic expansion which can overcome these two difficulties with ease. Starting from Babich's expansion, we develop high-order Eulerian asymptotics for Helmholtz equations in inhomogeneous media. Both the eikonal and transport equations are solved by high-order Lax---Friedrichs weighted non-oscillatory (WENO) schemes. We also prove that fifth-order Lax---Friedrichs WENO schemes for eikonal equations are convergent when the eikonal is smooth. Numerical examples demonstrate that new Eulerian high-order asymptotic methods are uniformly accurate in the neighborhood of the source and away from it.

Accurate simulation of seismic waves is of critical importance in a variety of geophysical applications. Based on recent works on staggered discontinuous Galerkin methods, we have developed a new method for the simulations of seismic waves, which has energy conservation and extremely low grid dispersion, so that it naturally provided accurate numerical simulations of wave propagation useful for geophysical applications and was a generalization of classical staggered-grid finite-difference methods. Moreover, it could handle with ease irregular surface topography and discontinuities in the subsurface models. Our new method discretized the velocity and the stress tensor on this staggered grid, with continuity imposed on different parts of the mesh. The symmetry of the stress tensor was enforced by the Lagrange multiplier technique. The resulting method was an explicit scheme, requiring the solutions of a block diagonal system and a local saddle point system in each time step, and it was, therefore, very efficient. To tailor our scheme to Rayleigh waves, we developed a mortar formulation of our method. Specifically, a fine mesh was used near the free surface and a coarse mesh was used in the rest of the domain. The two meshes were in general not matching, and the continuity of the velocity at the interface was enforced by a Lagrange multiplier. The resulting method was also efficient in time marching. We also developed a stability analysis of the scheme and an explicit bound for the time step size. In addition, we evaluated some numerical results and found that our method was able to preserve the wave energy and accurately computed the Rayleigh waves. Moreover, the mortar formulation gave a significant speed up compared with the use of a uniform fine mesh, and provided an efficient tool for the simulation of Rayleigh waves.

Multiarrival Green’s functions are essential in seismic modeling, migration, and inversion. Huygens-Kirchhoff (HK) integrals provide a bridge to integrate locally valid first-arrival Green’s functions into a globally valid multiarrival Green’s function. We have designed robust and accurate finite-difference methods to compute first-arrival traveltimes and amplitudes, so that first-arrival Green’s functions can be constructed rapidly. We adapted a fast butterfly algorithm to evaluate discretized HK integrals. The resulting fast Huygens’ sweeping method has the following unique features: (1) it precomputes a set of local traveltime and amplitude tables, (2) it automatically takes care of caustics, (3) it constructs Green’s functions of the Helmholtz equation for arbitrary frequencies and for many point sources, and (4) for a fixed number of points per wavelength, it constructs each Green’s function in nearly optimal complexity [Formula: see text] in terms of the total number of mesh points [Formula: see text], where the prefactor of the complexity only depends on the specified accuracy, and is independent of the frequency. The 2D and 3D examples revealed the performance of the method.

We propose a new methodology for carrying out eikonal based traveltime tomography arising from important applications such as seismic imaging and medical imaging. The new method formulates the traveltime tomography problem as a variational problem for a certain cost functional explicitly with respect to both traveltime and sound speed. Furthermore, the cost functional is penalized to enforce the nonlinear equality constraint associated with the underlying eikonal equation, bihar- monically regularized with respect to traveltime, and harmonically regularized with respect to sound speed. To overcome the difficulty associated with the inherent nonlinearity of the eikonal equation, the Euler-Lagrange equation of the penalized-regularized variational problem is reformulated into an equivalent, mixed optimality system. This mixed system is associated with an initial value problem which is solved by an operator-splitting based solution method, and the splitting approach effec- tively reduces the optimality system into three nonlinear subproblems and three linear subproblems. Moreover, the nonlinear subproblems can be solved pointwise, while the linear subproblems can be reduced to linear second-order elliptic problems. Numerical experiments show that the new method can carry out traveltime tomography successfully and recover sound speeds efficiently.

We have developed a local level-set method for inverting 3D gravity-gradient data. To alleviate the inherent nonuniqueness of the inverse gradiometry problem, we assumed that a homogeneous density contrast distribution with the value of the density contrast specified a priori was supported on an unknown bounded domain [Formula: see text] so that we may convert the original inverse problem into a domain inverse problem. Because the unknown domain [Formula: see text] may take a variety of shapes, we parametrized the domain [Formula: see text] by a level-set function implicitly so that the domain inverse problem was reduced to a nonlinear optimization problem for the level-set function. Because the convergence of the level-set algorithm relied heavily on initializing the level-set function to enclose the gravity center of a source body, we applied a weighted [Formula: see text]-regularization method to locate such a gravity center so that the level-set function can be properly initialized. To rapidly compute the gradient of the nonlinear functional arising in the level-set formulation, we made use of the fact that the Laplacian kernel in the gravity force relation decayed rapidly off the diagonal so that matrix-vector multiplications for evaluating the gradient can be accelerated significantly. We conducted extensive numerical experiments to test the performance and effectiveness of the new method.

We propose a hybrid approach to solve the high-frequency Helmholtz equation with point source terms in smooth heterogeneous media. The method is based on the ray-based finite element method (ray-FEM) [29] , whose original version can not handle the singularity close to point sources accurately. This pitfall is addressed by combining the ray-FEM, which is used to compute the smooth far-field of the solution accurately, with a high-order asymptotic expansion close to the point source, which is used to properly capture the singularity of the solution in the near-field. The method requires a fixed number of grid points per wavelength to accurately represent the wave field with an asymptotic convergence rate of O ( ω − 1 / 2 ) , where ω is the frequency parameter in the Helmholtz equation. In addition, a fast sweeping-type preconditioner is used to solve the resulting linear system. We present numerical examples in 2D to show both accuracy and efficiency of our method as the frequency increases. In particular, we provide numerical evidence of the convergence rate, and we show empirically that the overall complexity is O ( ω 2 ) up to a poly-logarithmic factor.

We propose a multiple level set method for inverting three-dimensional magnetic data induced by magnetization only. To alleviate inherent non-uniqueness of the inverse magnetic problem, we assume that the subsurface geological structure consists of several uniform magnetic mass distributions surrounded by homogeneous non-magnetic background such as soil, where each magnetic mass distribution has a known constant susceptibility and is supported on an unknown sub-domain. This assumption enables us to reformulate the original inverse magnetic problem into a domain inverse problem for those unknown domains defining the supports of those magnetic mass distributions. Since each uniform mass distribution may take a variety of shapes, we use multiple level-set functions to parameterize these domains so that the domain inverse problem can be further reduced to an optimization problem for multiple level-set functions. To compute rapidly gradients of the nonlinear functional arising in the multiple-level-set formulation, we utilize the fact that the kernel function in the field-susceptibility relation decays rapidly off the diagonal so that matrix-vector multiplications for evaluating the gradients can be speeded up significantly. Numerical experiments are carried out to illustrate the effectiveness of the new method.

We present a ray-based finite element method for the high-frequency Helmholtz equation in smooth media, whose basis is learned adaptively from the medium and source. The method requires a fixed number of grid points per wavelength to represent the wave field; moreover, it achieves an asymptotic convergence rate of $$\mathcal {O}(\omega ^{-\frac{1}{2}})$$ O ( ω - 1 2 ) , where $$\omega $$ ω is the frequency parameter in the Helmholtz equation. The local basis is motivated by the geometric optics ansatz and is composed of polynomials modulated by plane waves propagating in a few dominant ray directions. The ray directions are learned by processing a low-frequency wave field that probes the medium with the same source. Once the local ray directions are extracted, they are incorporated into the local basis to solve the high-frequency Helmholtz equation. This process can be continued to further improve the approximations for both local ray directions and high-frequency wave fields iteratively. Finally, a fast solver is developed for solving the resulting linear system with an empirical complexity $$\mathcal {O}(\omega ^d)$$ O ( ω d ) up to a poly-logarithmic factor. Numerical examples in 2D are presented to corroborate the claims.

The numerical approximation of high-frequency wave propagation in inhomogeneous media is a challenging problem. In particular, computing high-frequency solutions by direct simulations requires several points per wavelength for stability and usually requires many points per wavelength for a satisfactory accuracy. In this paper, we propose a new method for the acoustic wave equation in inhomogeneous media in the time domain to achieve superior accuracy and stability without using a large number of unknowns. The method is based on a discontinuous Galerkin discretization together with carefully chosen basis functions. To obtain the basis functions, we use the idea from geometrical optics and construct the basis functions by using the leading order term in the asymptotic expansion. Also, we use a wavefront tracking method and a dimension reduction procedure to obtain dominant rays in each cell. We show numerically that the accuracy of the numerical solutions computed by our method is significantly higher than that computed by the IPDG method using polynomials. Moreover, the relative errors of our method grow only moderately as the frequency increases.