Xiaofan Li

Scalar Seismic-Wave Equation Modeling by a Multisymplectic Discrete Singular Convolution Differentiator Method

High-precision modeling of seismic-wave propagation in heterogeneous media is very important to seismological investigation. However, such modeling is one of the difficult problems in the seismological research fields. For developing methods of seismic inversion and high-resolution seismic-wave imaging, the modeling problem must be solved as perfectly as possible. Moreover, for long-term computations of seismic waves (e.g., Earth’s free-oscillations modeling and seismic noise-propagation modeling), the capability of seismic modeling methods for long-time simulations is in great demand. In this paper, an alternative method for accurately and efficiently modeling seismic wave fields is presented; it is based on a multisymplectic discrete singular convolution differentiator scheme (MDSCD). This approach uses optimization and truncation to form a localized operator, which preserves the fine structure of the wave field in complex media and avoids noncausal interaction when parameter discontinuities are present in the medium. The approach presented has a structure-preserving property, which is suitable for treating questions of high-precision or long-time numerical simulations. Our numerical results indicate that this method can suppress numerical dispersion and allow for research into long-time numerical simulations of wave fields. These numerical results also show that the MDSCD method can effectively capture the inner interface without any special treatment at the discontinuity.

A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling

We have developed a mixed-grid finite element method (MGFEM) to simulate seismic wave propagation in 2D structurally complex media. This method divides the physical domain into two subdomains. One subdomain covering the major part of the physical domain is divided by regular quadrilateral elements, while the other subdomain uses triangular elements to correctly fit a rugged free surface topography. The local stiffness matrix of any quadrilateral element is identical and matrix-vector production is calculated using an element-by-element technique, which avoids assembling a huge global stiffness matrix. As only a few triangular elements exist in the subdomain containing the rugged free surface topography, the memory requirements for storing the assembled subdomain global stiffness matrix are significantly reduced. To eliminate artificial boundary reflections, the MGFEM is also implemented to solve the system equations of PML absorbing boundary conditions (PML ABC). The accuracy and efficiency of the MGFEM is tested in numerical experiments by comparing it with conventional methods, and numerical comparisons also indicate its tremendous ability to describe rugged surfaces.

A new kind of optimal second-order symplectic scheme for seismic wave simulations

Here we introduce generalized momentum and coordinate to transform seismic wave displacement equations into Hamiltonian system. We define the Lie operators associated with kinetic and potential energy, and construct a new kind of second order symplectic scheme, which is extremely suitable for high efficient and long-term seismic wave simulations. Three sets of optimal coefficients are obtained based on the principle of minimum truncation error. We investigate the stability conditions for elastic wave simulation in homogeneous media. These newly developed symplectic schemes are compared with common symplectic schemes to verify the high precision and efficiency in theory and numerical experiments. One of the schemes presented here is compared with the classical Newmark algorithm and third order symplectic scheme to test the long-term computational ability. The scheme gets the same synthetic surface seismic records and single channel record as third order symplectic scheme in the seismic modeling in the heterogeneous model.

Low-Frequency 3D Wave Propagation Modeling of the 12 May 2008 Mw 7.9 Wenchuan Earthquake

The seismic potential of southern China is associated with the collision between the Indian and the Eurasian plates. This is manifested in the western Sichuan Plateau by several seismically active systems of faults, such as the Longmen Shan. The seismicity observed on the Longmen Shan fault includes recent events with magnitudes of up to 6.5, and the one of 12 May 2008 Mw 7.9 Wenchuan earth- quake. Herewith, as part of an ongoing research program, a recently optimized three- dimensional (3D) seismic wave propagation parallel finite-difference code was used to obtain low-frequency (≤ 0:3 Hz) 3D synthetic seismograms for the Wenchuan earthquake. The code was run on KanBalam (Universidad Nacional Autonoma de Mexico, Mexico) and HECToR (UK National Supercomputing Service) supercompu- ters. The modeling included the U.S. Geological Survey 40 × 315 km 2 kinematic description of the earthquakes rupture, embedded in a 2400 × 1600 × 300 km 3 phys- ical domain, spatially discretized at 1 km in the three directions and a temporal dis- cretization of 0.03 s. The compression and shear wave velocities and densities of the geologic structure used were obtained from recently published geophysical studies performed in the Sichuan region. The synthetic seismograms favorably compare with the observed ones for several station sites of the Seismological and Accelerographic Networks of China, such as MZQ, GYA, and TIY, located at about 90, 500, and 1200 km, respectively, from the epicenter of the Wenchuan event. Moreover, the com- parisons of synthetic displacements with differential radar interferometry (DinSAR) ground deformation imagery, as well as of maximum velocity synthetic patterns with Mercalli modified intensity isoseist of the 2008 Wenchuan earthquake, are acceptable. 3D visualizations of the propagation of the event were also obtained; they show the source rupture directivity effects of the Mw 7.9 Wenchuan event. Our results partially explain the extensive damage observed on the infrastructure and towns located in the neighborhood of the Wenchuan earthquake rupture zone.

Seismic scalar wave equation with variable coefficients modeling by a new convolutional differentiator

Studying seismic wavefields in the Earths interior requires an accurate calculation of wave propagation using accurate and efficient numerical techniques. In this paper, we present an alternative method for accurately and efficiently modeling seismic wavefields using a convolutional generalized orthogonal polynomial differentiator. Our approach uses optimization and truncation to form a localized operator. This preserves the fine structure of the wavefield in complex media and avoids non-causal interaction when parameter discontinuities are present in the medium. We demonstrate this approach for scalar wavefield modeling in heterogeneous media and conclude that the method could be readily extended to elastic wavefield calculations. Our numerical results indicate that this method can suppress numerical dispersion and allow for the study of wavefields in heterogeneous structures. The results hold promise not only for future seismic studies, but also for any field that requires high-precision numerical solution of partial differential equation with variable coefficients.

Finite-element implementation of reverse-time migration for anisotropic media

Reverse-time migration for post-stack seismic data in anisotropic media is implemented using the finite-element method. As an accurate digital method, the finite-element method is flexible for dealing with complicated geological structures, inner and man-made boundaries despite its intensive computation. Applying it in reverse-time migration may produce accurate images for anisotropic media. To eliminate man-made boundary reflections, the absorbing boundary condition for anisotropic elastic waves is also studied. An efficient and stable absorbing boundary scheme is presented combining a discrete transparent boundary condition with an attenuation boundary condition.

A symplectic method for structure-preserving modelling of damped acoustic waves

In this paper, a symplectic method for structure-preserving modelling of the damped acoustic wave equation is introduced. The equation is traditionally solved using non-symplectic schemes. However, these schemes corrupt some intrinsic properties of the equation such as the conservation of both precision and the damping property in long-term calculations. In the method presented, an explicit second-order symplectic scheme is used for the time discretization, whereas physical space is discretized by the discrete singular convolution differentiator. The performance of the proposed scheme has been tested and verified using numerical simulations of the attenuating scalar seismic-wave equation. Scalar seismic wave-field modelling experiments on a heterogeneous medium with both damping and high-parameter contrasts demonstrate the superior performance of the approach presented for suppression of numerical dispersion. Long-term computational experiments display the remarkable capability of the approach presented for long-time simulations of damped acoustic wave equations. Promising numerical results suggest that the approach is suitable for high-precision and long-time numerical simulations of wave equations with damping terms, as it has a structure-preserving property for the damping term.

Pre-stack full wavefield inversion for elastic parameters of TI media

Pre-stack full wavefield inversion for the elastic parameters of transversely isotropical media is implemented. The Jacobian matrix is derived directly with the finite element method, just like the full wavefield forward modelling. An absorbing boundary scheme combining Liaos transparent boundary condition with Sarmas attenuation boundary condition is applied to the forward modelling and Jacobian calculation. The input data are the complete ground-recorded wavefields containing full kinematic and dynamic information for the seismic waves. Inversion with such data is desirable as it should improve the accuracy of the estimated parameters and also reduce data pre-processing, such as wavefield identification and separation. A scheme called energy grading inversion is presented to deal with the instability caused by the large energy difference between different arrivals in the input data. With this method, parameters in the shallow areas, which mainly affect wave patterns with strong energy, converge before those of deeper media. Thus, the number of unknowns in each inversion step is reduced, and the stability and reliability of the inversion process is greatly improved. As a result, the scheme is helpful to reduce the non-uniqueness in the inversion. Two synthetic examples show that the inversion system is reliable and accurate even though initial models deviate significantly from the actual models. Also, the system can accurately invert for transversely isotropic model parameters even with the introduction of strong random noise.

2D efficient ray tracing with a modified shortest path method

The computation effort of ray tracing with the shortest path method (SPM) is strongly dependent on the number of the discretized nodes in a model and the number of ray directions emanating from a secondary source node. In the traditional SPM, a secondary source emanates rays to all the surrounding nodes. Obviously, most of them are not minimal traveltime raypaths. As a result, the efficiency of SPM can be greatly improved if some measures are taken to avoid those unnecessary computations. In the current study, we apply the traveltime information of neighbouring nodes and the incident rays to determine the effective target propagation directions of secondary source nodes in 2D case. Generally, the effective propagation directions are narrow bands with few surrounding nodes. Thus, most unnecessary ray directions of secondary source nodes are avoided. 2D model tests show that the computational speed of the improved method is about several to tens of times of that of the traditional SPM with the increase of network nodes and cells.

Ray tracing with the improved shortest path method

The shortest path ray tracing is a kind of efficient and flexible method to calculate global minimum traveltimes and raypaths. However, the derived zigzag raypaths can generate errors in both traveltimes and loci. Increasing the number of nodes and propagation directions from a secondary source can improve the situation. But it also increases much computation. Here, several measures are presented to improve the efficiency and accuracy of the shortest path method (SPM). First, Snell’s law is considered when searching wave propagation directions from a secondary source. Thus, a large number of unnecessary directions are discarded. Second, when tracing reflections, a measure is applied to keep from tracing areas below interfaces. Third, traveltimes and raypath loci are corrected with a nonlinear interpolation scheme. In addition, a special SPM, named interface point method, is presented. It can tracing blocky models with high efficiency and accuracy.