Staggered-grid Finite-difference Method Research Articles

A time-domain SGFD-FK hybrid method for 2D teleseismic elastic wave modeling and inversion

AbstractFull waveform inversion (FWI) has proved to be a reliable tool for high-resolution imaging of lithospheric structures at various depths down to the upper mantle. However, when the size of the model is large, the computational burden is significative and applications are restricted to low frequencies. To tackle this issue, we developed a new 2D time-domain hybrid method to simulate high-frequency teleseismic body waves propagating through a local heterogeneous elastic Earth model: The frequency-wavenumber (FK) integration method is coupled with the staggered grid finite difference method (SGFD). The FK method is used to compute the wavefield due to obliquely incident plane P and SV waves in a 1D multilayered half space that excites a local heterogeneous region. Inside this region, the velocity-stress staggered grid FD method (SGFD) is used to accurately deal with wave propagation in heterogeneous media. Spurious waves that might be generated at the boundaries of the local region are avoided using convolutional perfectly matched layers (CPML). This new hybrid method inherits the low-memory requirements of the FK method and the accuracy, efficiency and easy implementation of the SGFD. The new hybrid method is benchmarked against the analytical FK method for some canonical models and shows good agreement with analytical solutions. Subsequently, our modeling tool is incorporated into a full waveform inversion algorithm adapted for teleseismic configurations to invert the incident P wave and its coda. The inversion is carried out using a gradient approach that is efficiently implemented via the adjoint-state method. The results suggest that our hybrid method and FWI algorithm represent a valuable tool for 2D forward and inverse regional applications using teleseismic data sets.

Fractional-order velocity-dilatation-rotation viscoelastic wave equation and numerical solution based on a constant-Q model

Viscoelastic wave equations based on the constant- Q (CQ) model can accurately describe the amplitude dissipation and phase distortion of waves in anelastic media. However, only three velocity or displacement components can be obtained directly by solving such equations. Starting from the time-domain second-order displacement viscoelastic wave equation, we derived the decoupled P- and S-wave displacement vector viscoelastic wave equation by using the polarization difference of P- and S-wave propagation in isotropic media. The equation can be transformed into the velocity-dilatation-rotation viscoelastic wave equation containing the first-order temporal derivative and fractional Laplacian operators, which can be solved directly by using the staggered-grid finite-difference and pseudospectral methods. We use the low-rank decomposition method to approximate the derived mixed space-wavenumber domain fractional Laplacian operators for modeling wave propagation in heterogeneous attenuating media. We also demonstrated the precision of our equation by comparing the numerical solutions with the analytical solutions. Furthermore, compared with the conventional velocity-stress viscoelastic wave equation, experimental results demonstrate that our equation can separate the pure P and S waves from the mixed wavefield during wavefield continuation. In addition, it can be separated into an equation containing predominantly an amplitude attenuation or a phase distortion term.

A linear, second-order accurate, positivity-preserving and unconditionally energy stable scheme for the Navier–Stokes–Poisson–Nernst–Planck system

Electrokinetic phenomena involving the transport of ions in fluids generally appears in biological system, hydrodynamics and geodynamics. The underlining mathematical models involve the nonlinear coupling between the fluid motion and ion transportation. Performing efficient, accurate and structure preserving computation for these electrohydrodynamic models is challenging due to their multiphysics and complex nature. The Navier–Stokes–Poisson–Nernst–Planck (NSPNP) system is a widely used electrokinetic model in the existing literature. In this paper, a linear, second-order accurate, positivity-preserving and unconditionally energy stable scheme is developed for the NSPNP system. First, a square root transformation is employed to ensure the positivity of ion concentrations. Next, a specific ordinary differential equation is introduced to deal with the nonlinear coupling terms satisfying the “zero-energy-contribution” feature, which is combined with the scalar auxiliary variable method for nonlinear potentials to achieve a decoupled energy stable scheme. Therefore, only several linearly independent equations need to be solved at each time step. The staggered grid finite difference method is used for the spatial discretizations, and the unconditional energy stability of the fully discrete scheme is established. Finally, the unique solvability of the scheme is proved. Numerical experiments are carried out to demonstrate the accuracy and performance of the proposed scheme, and the effect of electric Coulomb force on the flow and motion of ions are also illustrated.

High-Accuracy Simulation of Rayleigh Waves Using Fractional Viscoelastic Wave Equation

The propagation of Rayleigh waves is usually accompanied by dispersion, which becomes more complex with inherent attenuation. The accurate simulation of Rayleigh waves in attenuation media is crucial for understanding wave mechanisms, layer thickness identification, and parameter inversion. Although the vacuum formalism or stress image method (SIM) combined with the generalized standard linear solid (GSLS) is widely used to implement the numerical simulation of Rayleigh waves in attenuation media, this type of method still has its limitations. First, the GSLS model cannot split the velocity dispersion and amplitude attenuation term, thus limiting its application in the Q-compensated reverse time migration/full waveform inversion. In addition, GSLS-model-based wave equation is usually numerically solved using staggered-grid finite-difference (SGFD) method, which may result in the numerical dispersion due to the harsh stability condition and poses complexity and computational burden. To overcome these issues, we propose a high-accuracy Rayleigh-waves simulation scheme that involves the integration of the fractional viscoelastic wave equation and vacuum formalism. The proposed scheme not only decouples the amplitude attenuation and velocity dispersion but also significantly suppresses the numerical dispersion of Rayleigh waves under the same grid sizes. We first use a homogeneous elastic model to demonstrate the accuracy in comparison with the analytical solutions, and the correctness for a viscoelastic half-space model is verified by comparing the phase velocities with the dispersive images generated by the phase shift transformation. We then simulate several two-dimensional synthetic models to analyze the effectiveness and applicability of the proposed method. The results show that the proposed method uses twice as many spatial step sizes and takes 0.6 times that of the GSLS method (solved by the SGFD method) when achieved at 95% accuracy.

Fast 3D Magnetotelluric Forward Modeling Based on a Regularization Approach in Spherical Coordinate

The earth’s curvature shouldn’t be ignored for the large-scale magnetotelluric (MT) modeling. A forward code of spherical coordinate MT method based on staggered grid finite difference method is developed. In addition, the forward problem is usually solved numerically, mostly based on iterative solution. Krylov subspace method is widely used to solve electromagnetic forward modeling problems in frequency domain. However, as the period increases, the convergence of these methods becomes significantly slower (especially for long-period MT data). We added a scaled grad-div operator to the double curl equation of the electric field to perform the nondivergence condition (CCGD). The results show that CCGD method is more efficient than other methods, which is an efficient method for our code.

3D Forward and Inversion Research on Frequency Domain Electromagnetic Tipper Data

The tipper is one of the main parameters of the electromagnetic exploration data in the frequency domain and has an important role and significant advantage in the exploration under complex geological conditions. In this paper, the 3D forward modeling is carried out based on the staggered grid finite difference method, and the analytical solution of the layered medium is used to verify its correctness; the L-BFGS algorithm is used to perform the 3D inversion of the tipper data. Through the three-dimensional complex model, the forward response characteristics and spatial distribution of the tipper are analyzed. The inversion results show that the inversion using the inversion data has less influence on static displacement and better inversion results.

A simple and high-efficiency viscoacoustic reverse time migration calculated by finite difference

The viscoacoustic wave equation with decoupled amplitude attenuation and phase dispersion terms is convenient for implementing reverse time migration (RTM) in attenuating media. However, the traditional viscoacoustic wave equations are expressed by the fractional Laplacian, which requires computationally expensive fast Fourier transforms to perform numerical solutions. This is not conducive for commercial applications with seismic data in terabyte, particularly in three dimensions. Based on the memory-variable-represented viscoacoustic wave equation, the new dispersion- and dissipation-dominated wavefield extrapolation operators that can be solved using the economical finite-difference method (FDM) are developed. In numerical simulation, the staggered-grid FDM is used to solve the temporal and spatial derivatives. The computational efficiency analysis indicates that, compared to the traditional viscoacoustic wave equation, the new viscoacoustic wave equation has an efficiency improvement of approximately more than two times. The accuracy of simulations of amplitude attenuation and phase dispersion is evaluated in viscoacoustic media with varying Q values. The proposed viscoacoustic wave equation is also extended to Q-compensated RTM. Two synthetic tests and one field data test indicate that the Q-compensated RTM can produce high-resolution images by correcting phase dispersion and compensating amplitude loss.

Higher-order thin layer method (HTLM) based wavefield modeling approach

Traditional free vibration-based forward modeling approaches provide only the theoretical modal solution that lacks information regarding modal amplitudes. Moreover, they ignore the influence of source-receiver layout on the wavefield. Numerical wavefield modeling approaches can generate a dispersion spectrum mimicking the actual field test; however, numerical methods are not efficient in terms of computational time and produce numerical artifacts. Furthermore, they can be contaminated with leaky waves for soil strata with a high Poisson's ratio. This study presents a new wavefield modeling approach for developing full dispersion spectra of stratified media incorporating source-receiver configuration. The proposed method calculates the theoretical modal amplitude using the higher-order thin layer method and determines the vertical displacement responses at specified virtual receiver positions in the frequency-space domain. Afterward, it searches for spectral values at each frequency over a broad range of trial phase velocities to create a full dispersion spectrum. The accuracy and computational efficacy of the present method are demonstrated on different soil profiles adopted from literature, including normally dispersive media, irregularly dispersive media with high-velocity contrast, and a crustal level model. Results indicate the proposed approach can clearly distinguish the relative spectral amplitude of each mode and dominant dispersion trend. Calculated seismograms and dispersion images are compared with the existing 2D staggered grid finite difference method. Although both methods yield similar results, the current method is at least two orders (x102) faster than the staggered grid finite difference method. Unlike numerical approaches, the proposed method does not require additional computational resources to model materials with a high Poisson's ratio. Furthermore, it minimizes the contamination of leaky waves in the dispersion spectrum when modeling a high Poisson's ratio. The proposed method does not produce any numerical artifacts and can be used to generate accurate synthetic seismograms as well as dispersion spectra. In addition, the proposed wavefield modeling approach can be effectively used to develop an accurate and faster inversion scheme by including the entire dispersion spectrum.

Three-Dimensional Ultrasonic Reverse-Time Migration Imaging of Submarine Pipeline Nondestructive Testing in Cylindrical Coordinates

Submarine pipelines are a safe and energy-efficient mode of gas transport. However, due to the complex manufacturing process and harsh operating environment, submarine pipelines are subject to fatigue cracks under long-term cyclic loading. A comprehensive and high-precision characterization strategy for submarine pipelines can effectively prevent potential safety hazards and have significant economic and social repercussions. As a matter of fact, pipeline defects cannot be reliably detected with current traditional 2D methods. On the other hand, in ultrasonic testing, cylindrical geometry increases the complexity of the 3D wave field in the submarine pipeline space and significantly influences the accuracy of the detection results. In this paper, we put forward a novel method for 3D ultrasonic image testing that is suitable for cylindrical coordinates. In order to accurately simulate the ultrasonic signal received from pipelines, we generalize the 3D staggered-grid finite-difference method from Cartesian coordinates to cylindrical ones and simulate the full wave field in the 3D pipeline space. Then, signal processing is performed on the ultrasound simulation records, and 3D reverse-time migration imaging of submarine pipeline defects can be effectively achieved using the reverse-time migration method and cross-correlation imaging conditions. The results obtained from simulations and real field data show that the proposed method provides high-quality 3D imaging of defects in pipelines, taking into account multiple scattering and mode conversion information at the bottom of the defects.

A fast and accurate elastic reverse-time migration method based on decoupled elastic wave equations

The extrapolation and decomposition of different mode waves are very important to efficiently obtain multi-wave depth images for elastic reverse-time migration (ERTM). To achieve this goal successfully, we firstly develop an efficient and accurate staggered-grid finite-difference method (SFDM) based on decoupled elastic wave equations. This developed SFDM uses two sets of difference operator lengths for decoupled elastic wave equations to respectively perform the extrapolation of decoupled P- and S-wavefields and employs the least squares (LS) to minimize the objective functions respectively constructed from the time-space dispersion equations of P- and S-waves to generate their corresponding difference coefficients. This developed method can theoretically be regarded as the optimal SFDM based on least squares (LSSFDM), and it can provide the high-precision separated wavefields of different mode waves simultaneously and reduce the computing time effectively for elastic wave modeling under the given accuracy requirements. The modeling results of different models illustrate that our developed LSSFDM based on decoupled elastic wave equations can generate the decomposed P- and S-wavefields and suppress their numerical dispersions effectively. Secondly, we present a fast and accurate ERTM method based on this developed LSSFDM to successfully achieve the multi-wave imaging results, including PP, PS, SP and SS images. The ERTM examples of different models demonstrate that our ERTM based on this developed LSSFDM can accurately obtain the high-precision depth images of multimode waves and effectively reduce the computing time of wavefield extrapolation.

Elastic full-wave field simulation in 3D tunnel space with the variable staggered-grid finite-difference method in cylindrical coordinates

Tunnel advance detection technology is an important method for determining the structure of complex geological bodies in front of the tunnel face. Among the tunnel advance detection technologies, the seismic method is one of the most accurate methods with long detection distances. In seismic tunnel advance detection, the cylindrical configuration aggravates the complexity of the wave field in the 3D tunnel space and significantly influences the accuracy of the detection results. Thus, it is crucial to simulate an accurate seismic full-wave field of the tunnel space and to understand the propagation and wave-field characteristics of individual seismic waves for seismic tunnel advance detection. Usually, in 3D Cartesian coordinates, the tunnel wall is approximated with a staircase boundary, but it is not sufficiently accurate in shape and generates numerical dispersion in the simulation, especially in the presence of surface waves. Therefore, we developed a variable staggered-grid finite-difference method in cylindrical coordinates to simulate the elastic full-wave field in a 3D tunnel space. Properly setting free-surface boundary conditions solves the propagation of surface waves on the tunnel wall and face. The interference of the instability and discontinuity of the pole axis in the seismic wave field simulation was eliminated using our method. According to practical engineering situations, models containing three types of geological bodies in front of the tunnel face were designed. Their seismic wave propagation and wave field characteristics were analyzed. Compared with the simulation results in Cartesian coordinates, the results in cylindrical coordinates show that numerical dispersion is negligible and conclude that a higher signal-to-noise ratio and more accurate seismic wave field can be simulated with cylindrical coordinates in the 3D tunnel space. Our simulation method provides theoretical and practical guidance for analyzing and interpreting seismic wave fields in tunnel advance detection.

Efficient temporal high-order staggered-grid scheme with a dispersion-relation-preserving method for the scalar wave modeling

Staggered-grid finite-difference (FD) method is widely used to solve the wave equation for the numerical seismic modeling, and it is a key step of the advanced seismic imaging and inversion problem. However, the conventional FD method is prone to instability and dispersion error due to the insufficient approximation accuracy. In this work, we propose an efficient temporal high-order finite-difference (FD) scheme with the cross-rhombus stencil. Compared with the standard cross-rhombus method, the new method has less computational cost due to we simplify the FD scheme. Moreover, the dispersion relation of the new method is easy to obtain the dispersion-relation-preserving (DRP) FD coefficients, which can significantly alleviate the spatial and temporal dispersion errors. Dispersion and stability analyses indicate that the new scheme has better performance in seismic modeling than the conventional method, and numerical experiments also indicate that the new scheme can effectively mitigate dispersion error and improve the numerical accuracy.

Modeling seismic wave propagation in the Loess Plateau using a viscoacoustic wave equation with explicitly expressed quality factor

The thick Quaternary loess on the Loess Plateau of China produces strong seismic attenuation, resulting in weak reflections from subsurface exploration targets. Accurately simulating seismic wavefield in the Loess Plateau is important for guiding subsequent data processing and interpretation. We present a 2D/3D wavefield simulation method for the Loess Plateau using a viscoacoustic wave equation with explicitly expressed quality factor. To take into account the effect of irregular surface, we utilize a vertically deformed grid to represent the topography, and solve the viscoacoustic wave equation in a regular computational domain that conforms to topographic surface. Grid deformation introduces the partial derivatives such as ∂vx/∂z and ∂vy/∂z in the wave equation, which is difficult to be accurately computed using traditional staggered-grid finite-difference method. To mitigate this issue, a finite-difference scheme based on a fully staggered-grid is adopted to solve the viscoacoustic wave equation. Numerical experiments for a simple layer model and 2D/3D realistic Loess Plateau models demonstrate the feasibility and adaptability of the proposed method. The 3D modeling results show comparable amplitude and waveform characteristics to the field data acquired from the Chinese Loess Plateau, suggesting a good performance of the proposed modeling method.

Three-dimensional elastic reverse-time migration using a high-order temporal and spatial staggered-grid finite-difference scheme

Three-dimensional (3D) elastic reverse-time migration (ERTM) can image the subsurface 3D seismic structures, and it is an important tool for the Earth’s interior imaging. A common simulation kernel used in 3D ERTM is the current staggered-grid finite-difference (SGFD) method of the first-order elastic wave equation. However, the mere second-order accuracy in time of the current SGFD method can bring non-negligible time dispersion, which reduces the simulation accuracy and further leads to the distortion of the imaging results. This paper proposes a vector-based 3D ERTM using the high-order accuracy SGFD method in time to obtain high-accuracy images. This approach is a new high-resolution ERTM workflow that improves the imaging accuracy of conventional ERTM from numerical simulation. The proposed ERTM workflow is established on a quasi-stress–velocity wave equation and its vector wavefield decomposition form. Advanced SGFD schemes and their corresponding coefficients with fourth-order temporal accuracy solve the quasi-linear wave equation system. The normalized dot product imaging condition produces high-quality images using high-accuracy vector wavefields solved using the SGFD method. Through the numerical examples, we test the simulation efficiency and analyze how temporal accuracy in numerical simulations affects migration imaging quality. We include that the proposed method obtains highly accurate images.

An optimized three-dimensional time-space domain staggered-grid finite-difference method

Numerical simulation of three-dimensional (3D) seismic wavefields forms the basis of the research on the migration methods of 3D seismic data based on wave equations. Because the simulation precision of wavefield extrapolation determines the imaging accuracy to a certain extent, it is very important to study how to enhance the forward modeling precision of 3D seismic wavefields. Thus, we build on an optimized 3D staggered-grid finite-difference (SFD) method with high simulation precision based on two-dimensional (2D) seismic modeling. Since it generates the corresponding difference coefficients by utilizing the least square (LS) method to minimize the objective function constructed by the time-space domain dispersion relation of the 3D acoustic wave equation, our optimized time-space domain LS-based 3D SFD method can effectively enhance the modeling precision of the 3D seismic wavefields in theory compared with the 3D SFD methods based on the Taylor-series expansion (TE), especially for the large wavenumber range. Examining the numerical dispersion, algorithm stability and computational cost, we compare our optimized time-space domain LS-based 3D SFD method with three conventional TE-based and LS-based 3D SFD methods to illustrate and demonstrate its effectiveness and feasibility. The numerical examples from different 3D models suggest that our optimized time-space domain LS-based 3D SFD method can generate less numerical dispersion and higher simulation accuracy for 3D seismic wavefields than three other conventional 3D SFD methods, but its stability condition is stricter and its computational cost is slightly higher.

Decoupled acoustoelasticity equations and stable FD modeling for wave propagation in prestressed media

Seismic exploration of deep oil/gas reservoirs involves wave propagation in high-pressure media. Acoustoelasticity with the third-order elastic constants can account for nonlinear strain responses due to stresses of finite magnitude. Since current seismic exploration depends primarily on P waves, further development in prospecting technology to handle deep reflection data will benefit from the ability to decouple the elastic waves from high-pressure deep formations. We decouple the first-order velocity-stress acoustoelasticity equations for wave propagation under isotropic (confining) and anisotropic (uniaxial and pure shear) prestress conditions. A rotated staggered-grid finite-difference (RSG-FD) method is used to solve the decoupled acoustoelasticity equations with an unsplit convolutional perfectly matched layer (CPML) absorbing boundary. Numerical examples demonstrate the significant impact of prestressed conditions on seismic responses in both phase and amplitude. The stress-induced seismic anisotropy with orthotropic characteristics is strongly related to the orientation of prestresses. We compare the coupled and the decoupled wavefields to validate the proposed decoupling method. Investigation of the influence of prestresses on the acoustoelastic decoupling illustrates that the isotropic prestress does not affect the decoupling, but the anisotropic prestresses significantly affect the decoupling accuracy.

Propagation Characteristics of Seismic Waves in Shallow Sea Sedimentary Layers

Zhao, S.; Gu, D.; Guo, C.; Wang, X., and Zhang, J., 2023. Propagation characteristics of seismic waves in shallow sea sedimentary layers. Journal of Coastal Research, 39(1), 32–44. Charlotte (North Carolina), ISSN 0749-0208. To explore new techniques for detecting and identifying targets in water, the propagation characteristics of seismic waves in shallow sea were deeply studied. Wave theory and staggered grid finite difference method were employed to analyze the propagation characteristics of seismic waves in a layer of underwater medium with a sedimentary layer. The results show that a Scholte wave is prone to disperse in an underwater medium with a sedimentary layer. The sediment thickness and acoustic parameters have great influence on the propagation characteristics of seismic waves. Moreover, the seismic wave signals were mainly concentrated below 20 Hz by collecting the vibration velocity of underwater interface. Furthermore, Scholte wave signals in underwater experiments were determined by utilizing four filtering methods and Scholte polarization characteristics, and the Butterworth filtering method can better preserve the integrity of vibration velocity signal. The theoretical, numerical, and experimental polarization characteristic results are basically consistent, which indicate that Scholte wave can be well extracted by the method of polarization characteristics.

Cold and Wet Mantle Transition Zone Beneath the Mediterranean Sea Revealed by the Electrical Image

A three-dimensional electrical conductivity model of the mantle beneath North Africa and Southern Europe is obtained by geomagnetic depth sounding. C-responses are estimated from geomagnetic data observed at observatories in and around the region and converted to the electrical structure of the mantle transition zone. The limited-memory quasi-Newton method is chosen to minimize the nonlinear objective function of inversion, while the forward modeling relies on a staggered-grid finite difference method in the spherical coordinate system. The data misfit term of the inversion function is measured using the L1-norm in order to suppress the response instability caused by the significant noise contained in the observed data. In order to adequately correct the ocean effect in observatories near the coast, a fixed shell comprised of ocean and land is incorporated in inversion. A banded zone with reduced conductivity is present in the three-dimensional model, primarily seen in the lower mantle transition zone and lower mantle beneath the Mediterranean Sea. Combining laboratory-measured conductivity models, we propose that subducted slabs causing reduced temperature and a water reservoir in the mantle transition zone should be responsible for the observed electrical model.

Dispersive staggered grid finite difference modelling of Rayleigh waves

Non-destructive surface wave-based ground exploration techniques such as multichannel analysis of surface wave (MASW) require a rigorous numerical model to predict the soil properties accurately. The staggered grid finite difference (SGFD) method is widely used to simulate elastic wave propagation due to its stability and ease of implementation. However, in the presence of low-velocity layers, the SGFD method requires computationally expensive high spatial grid resolution to accurately model the surface waves. On the other hand, velocity dependent grid approaches become unstable for the highly irregular soil profile. In order to overcome these challenges, a vertically non-uniform dispersive staggered grid scheme is proposed that significantly reduces the computing time without sacrificing accuracy. The proposed non-uniform discretization strategy is based on the Rayleigh wave dispersion relation, where the penetration depth of each frequency is assumed to be roughly one and half of the wavelength. Based on the above assumption, an exponential relationship between grid size and grid number is developed. The proposed method calculates continuously increasing grid length while ensuring grid points at the exact layer interface. Five different soil profiles: a homogenous elastic half-space, an irregularly dispersive crustal level profile, a low velocity layer at a deeper depth model, a real-world model, and a lateral heterogeneous model are used to demonstrate the efficacy of the proposed approach. For all the cases, the results of the proposed dispersive grid match those of the uniform grid with a correlation coefficient value (R) > 0.99. Compared to the uniform and velocity dependent grid approaches, the proposed dispersive grid scheme requires ½ to 1/4th of elements and CPU time. Furthermore, unlike the velocity dependent grid scheme, the proposed method does not become unstable for a highly irregular profile. Results indicate that the proposed method is superior to the traditional uniform and velocity dependent grid approaches for modelling elastic waves in layered media.

Wavelet-Based Three-Dimensional Inversion for Geomagnetic Depth Sounding

The complexity of Earth’s structure poses a challenge to the multiscale detection capability of geophysics. In this paper, we present a new wavelet-based three-dimensional inversion method for geomagnetic depth sounding. This method is based on wavelet functions to transfer model parameters in the space domain into the wavelet domain. The model is represented by wavelet coefficients containing both large- and fine-scale information, enabling wavelet-based inversion to describe multiscale anomalies. L1-norm measurement is applied to measure the model roughness to accomplish the sparsity constraint in the wavelet domain. Meanwhile, a staggered-grid finite difference method in a spherical coordinate system is used to calculate the forward responses, and the limited-memory quasi-Newton method is applied to seek the solution of the inversion objective function. Inversion tests of synthetic data for multiscale models show that wavelet-based inversion is stable and has multiresolution. Although higher-order wavelets can lead to finer results, our tests present that a db6 wavelet is suitable for geomagnetic depth sounding inversion. The db6 inversion results of responses at 129 geomagnetic observatories around the world reveal a higher-resolution image of the mantle.