Jin-Hai Zhang

Optimized explicit finite-difference schemes for spatial derivatives using maximum norm

We propose an optimized scheme using the maximum norm and the simulated annealing.The maximum norm offers the largest set of possible solutions for solvers to search.The explicit finite-difference operator is greatly improved by our optimized scheme.We use a tight error limitation to make accuracy improvement to be high and solid.Our optimized scheme allows greater saving of computational cost and memory demand. Conventional explicit finite-difference methods have difficulties in handling high-frequency components due to strong numerical dispersions. One can reduce the numerical dispersions by optimizing the constant coefficients of the finite-difference operator. Different from traditional optimized schemes that use the 2-norm and the least squares, we propose to construct the objective functions using the maximum norm and solve the objective functions using the simulated annealing algorithm. Both theoretical analyses and numerical experiments show that our optimized scheme is superior to traditional optimized schemes with regard to the following three aspects. First, it provides us with much more flexibility when designing the objective functions; thus we can use various possible forms and contents to make the objective functions more reasonable. Second, it allows for tighter error limitation, which is shown to be necessary to avoid rapid error accumulations for simulations on large-scale models with long travel times. Finally, it is powerful to obtain the optimized coefficients that are much closer to the theoretical limits, which means greater savings in computational efforts and memory demand.

Accelerating 3D Fourier migration with graphics processing units

Computational cost is a major factor that inhibits the practicalapplicationof3Ddepthmigration.Wehavedevelopeda fastparallelschemetospeedup3Dwave-equationdepthmigration on a parallel computing device, i.e., on graphics processing units GPUs. The third-order optimized generalized-screen propagator is used to take advantage of the builtinsoftwareimplementationofthefastFouriertransform.The propagator is coded as a sequence of kernels that can be called from the computer host for each frequency component. Moving the wavefield extrapolation for each depth leveltotheGPUsallowshandlingalarge3Dvelocitymodel,but this scheme can be speeded up to a limited degree over the CPU implementation because of the low-bandwidth data transfer between host and device. We have created further speedupinthisextrapolationschemebyminimizingthelowbandwidth data transfer, which is done by storing the 3D velocity model and imaged data in the device memory, and reducing half the memory demand by compressing the 3D velocity model and imaged data using integer arrays instead of float arrays. By incorporating a 2D tapered function, timeshiftpropagator,andscalingoftheinverseFouriertransform intoacompactkernel,thecomputationtimeisreducedgreatly. Three-dimensional impulse responses and synthetic data exampleshavedemonstratedthattheGPU-basedFouriermigration typically is 25 to 40 times faster than the CPU-based implementation. It enables us to image complex media using 3D depth migration with little concern for computational cost.Themacrovelocitymodelcanbebuiltinamuchshorter turnaroundtime.

Volcanic history of the Imbrium basin: A close-up view from the lunar rover Yutu.

Significance After the Apollo and Luna missions, which were flown about 40 years ago, the Moon was explored only from orbit. In addition, no samples were returned from the young and high-FeO and TiO2 mare basalt in the northern Imbrium basin. Such samples are important to understand the formation and evolution of the Procellarum KREEP [potassium (K), rare earth elements (REE), and phosphorus (P)] terrain, a key terrain highly enriched in radioactive nuclides. The Chang’e-3 mission carried out the first in situ analyses of chemical and mineral compositions of the lunar soil and ground-based measurements of the lunar regolith and the underlying basalt units at this specific site. The lunar regolith layer recorded the surface processes of the Moon, whereas the basalt units recorded the volcanic eruption history. We report the surface exploration by the lunar rover Yutu that landed on the young lava flow in the northeastern part of the Mare Imbrium, which is the largest basin on the nearside of the Moon and is filled with several basalt units estimated to date from 3.5 to 2.0 Ga. The onboard lunar penetrating radar conducted a 114-m-long profile, which measured a thickness of ∼5 m of the lunar regolith layer and detected three underlying basalt units at depths of 195, 215, and 345 m. The radar measurements suggest underestimation of the global lunar regolith thickness by other methods and reveal a vast volume of the last volcano eruption. The in situ spectral reflectance and elemental analysis of the lunar soil at the landing site suggest that the young basalt could be derived from an ilmenite-rich mantle reservoir and then assimilated by 10–20% of the last residual melt of the lunar magma ocean.

Optimized Chebyshev Fourier migration: A wide-angle dual-domain method for media with strong velocity contrasts

A wide-angle propagator is essential when imaging complex media with strong lateral velocity contrasts in one-way waveequation migration. We have developed a dual-domain one-way propagator using truncated Chebyshev polynomials and a globally optimized scheme. Our method increases the accuracy of the expanded square-root operator by adding two high-order terms to the traditional split-step Fourier propagator. First, we approximate the square-root operator using Taylor expansion around the reference background velocity. Then, we apply the first-kindChebyshevpolynomialstoeconomizetheresultsofthe Taylor expansion. Finally, we optimize the constant coefficients using the globally optimized scheme, which are fixed and feasible for arbitrary velocity models. Theoretical analysis and numericalexperimentshavedemonstratedthatthemethodhasvery highaccuracyandexceedstheunoptimizedFourierfinite-difference propagator for the entire range of practical velocity contrasts. The accurate propagation angle of the method is always about 60° under the relative error of 1% for complex media with weak, moderate, and even strong lateral velocity contrasts. The method allows us to handle wide-angle propagations and strong lateral velocity contrast simultaneously by using Fourier transform alone. Only four 2D Fourier transforms are required for each step of 3D wavefield extrapolation, and the computing cost is similar to that of the Fourier finite-difference method. Compared with the finite-difference method, our method has no twowaysplittingerrori.e.,azimuthal-anisotropyerrorfor3Dcases and almost no numerical dispersion for coarse grids. In addition, it has strong potential to be accelerated when an enhanced fast Fouriertransformalgorithmemerges.

Comparison of artificial absorbing boundaries for acoustic wave equation modelling

Absorbing boundary conditions are necessary in numerical simulation for reducing the artificial reflections from model boundaries. In this paper, we overview the most important and typical absorbing boundary conditions developed throughout history. We first derive the wave equations of similar methods in unified forms; then, we compare their absorbing performance via theoretical analyses and numerical experiments. The Higdon boundary condition is shown to be the best one among the three main absorbing boundary conditions that are based on a one-way wave equation. The Clayton and Engquist boundary is a special case of the Higdon boundary but has difficulty in dealing with the corner points in implementaion. The Reynolds boundary does not have this problem but its absorbing performance is the poorest among these three methods. The sponge boundary has difficulties in determining the optimal parameters in advance and too many layers are required to achieve a good enough absorbing performance. The hybrid absorbing boundary condition (hybrid ABC) has a better absorbing performance than the Higdon boundary does; however, it is still less efficient for absorbing nearly grazing waves since it is based on the one-way wave equation. In contrast, the perfectly matched layer (PML) can perform much better using a few layers. For example, the 10-layer PML would perform well for absorbing most reflected waves except the nearly grazing incident waves. The 20-layer PML is suggested for most practical applications. For nearly grazing incident waves, convolutional PML shows superiority over the PML when the source is close to the boundary for large-scale models. The Higdon boundary and hybrid ABC are preferred when the computational cost is high and high-level absorbing performance is not required, such as migration and migration velocity analyses, since they are not as sensitive to the amplitude errors as the full waveform inversion. We provide a thorough review of all typical absorbing boundary conditions and derive their equations in a uniform mathematical form. We examine their performance via numerical experiments and qualitatively show their advantages and disadvantages. Finally, we provide some suggestions on choosing different boundary conditions for practical applications.

Third-order symplectic integration method with inverse time dispersion transform for long-term simulation

The symplectic integration method is popular in high-accuracy numerical simulations when discretizing temporal derivatives; however, it still suffers from time-dispersion error when the temporal interval is coarse, especially for long-term simulations and large-scale models. We employ the inverse time dispersion transform (ITDT) to the third-order symplectic integration method to reduce the time-dispersion error. First, we adopt the pseudospectral algorithm for the spatial discretization and the third-order symplectic integration method for the temporal discretization. Then, we apply the ITDT to eliminate time-dispersion error from the synthetic data. As a post-processing method, the ITDT can be easily cascaded in traditional numerical simulations. We implement the ITDT in one typical exiting third-order symplectic scheme and compare its performances with the performances of the conventional second-order scheme and the rapid expansion method. Theoretical analyses and numerical experiments show that the ITDT can significantly reduce the time-dispersion error, especially for long travel times. The implementation of the ITDT requires some additional computations on correcting the time-dispersion error, but it allows us to use the maximum temporal interval under stability conditions; thus, its final computational efficiency would be higher than that of the traditional symplectic integration method for long-term simulations. With the aid of the ITDT, we can obtain much more accurate simulation results but with a lower computational cost. We use the inverse time dispersion transform to the 3rd-order symplectic integration.It can greatly reduce the time-dispersion error, especially for long travel times.It is superior to the conventional 2nd- and high-order finite-difference methods.It allows the maximum temporal interval under stability conditions thus is efficient.We can obtain more accurate simulation results but with a lower computational cost.

Fast image inpainting using exponential- threshold POCS plus conjugate gradient

Abstract Image inpainting can remove unwanted objects and reconstruct the missing or damaged portions of an image. The projection onto convex sets (POCS) is a classical method used in image inpainting. However, the traditional POCS converges slowly due to the linear error threshold. We propose an exponential-threshold scheme, which greatly improves the convergence of the POCS. Although the exponential-threshold POCS can recover the image in about 20 iterations, it cannot reconstruct the image details very well even with hundreds of iterations. Thus, we append the non-local restoration to the exponential-threshold POCS to further refine the image details, and then we solve this objective function using the conjugate gradient. Numerical experiments show that for each iteration, the exponential-threshold POCS and the conjugate gradient have very similar computational efficiencies. For an image with various topologies of the missing areas, our scheme can recover missing pixels simultaneously and obtain a satisfied inpainting result in only 20 iterations of the exponential-threshold POCS and 20 iterations of the conjugate gradient. The proposed method can excellently restore damaged photographs and remove superimposed text. This method has less computational cost than the conjugate gradient and has a higher resolution than the POCS.

Reducing two-way splitting error of FFD method in dual domains

The Fourier finite-difference (FFD) method is very popular in seismic depth migration. But its straightforward 3D extension creates two-way splitting error due to ignoring the cross terms of spatial partial derivatives. Traditional correction schemes, either in the spatial domain by the implicit finite-difference method or in the wavenumber domain by phase compensation, lead to substantially increased computational costs or numerical difficulties for strong velocity contrasts. We propose compensating the two-way splitting error in dual domains, alternately in the spatial and wavenumber domains via Fourier transform. First, we organize the expanded square-root operator in terms of two-way splitting FFD plus the usually ignored cross terms. Second, we select a group of optimized coefficients to maximize the accuracy of propagation in both inline and crossline directions without yet considering the diagonal directions. Finally, we further optimize the constant coefficient of the compensation part to further improve the overall accuracy of the operator. In implementation, the compensation terms are similar to the high-order corrections of the generalized-screen method, but their functions are to compensate the two-way splitting error rather than the expansion error. Numerical experiments show that optimized one-term compensation can achieve nearly perfect circular impulse responses and the propagation angle with less than 1% error for all azimuths is improved up to 60 � from 35 � . Compared with traditional single-domain methods, our scheme can handle lateral velocity variations (even for strong velocity contrasts) much more easily with only one additional Fourier transform based on the two-way splitting FFD method, which helps retain the computational efficiency.

Restoration of clipped seismic waveforms using projection onto convex sets method

The seismic waveforms would be clipped when the amplitude exceeds the upper-limit dynamic range of seismometer. Clipped waveforms are typically assumed not useful and seldom used in waveform-based research. Here, we assume the clipped components of the waveform share the same frequency content with the un-clipped components. We leverage this similarity to convert clipped waveforms to true waveforms by iteratively reconstructing the frequency spectrum using the projection onto convex sets method. Using artificially clipped data we find that statistically the restoration error is ~1% and ~5% when clipped at 70% and 40% peak amplitude, respectively. We verify our method using real data recorded at co-located seismometers that have different gain controls, one set to record large amplitudes on scale and the other set to record low amplitudes on scale. Using our restoration method we recover 87 out of 93 clipped broadband records from the 2013 Mw6.6 Lushan earthquake. Estimating that we recover 20 clipped waveforms for each M5.0+ earthquake, so for the ~1,500 M5.0+ events that occur each year we could restore ~30,000 clipped waveforms each year, which would greatly enhance useable waveform data archives. These restored waveform data would also improve the azimuthal station coverage and spatial footprint.

Unsplit complex frequency shifted perfectly matched layer for second-order wave equation using auxiliary differential equations

The complex frequency shifted perfectly matched layer (CFS-PML) can improve the absorbing performance of PML for nearly grazing incident waves. However, traditional PML and CFS-PML are based on first-order wave equations; thus, they are not suitable for second-order wave equation. In this paper, an implementation of CFS-PML for second-order wave equation is presented using auxiliary differential equations. This method is free of both convolution calculations and third-order temporal derivatives. As an unsplit CFS-PML, it can reduce the nearly grazing incidence. Numerical experiments show that it has better absorption than typical PML implementations based on second-order wave equation.