Xiaohong Meng

3D seismic reverse time migration on GPGPU

Reverse time migration (RTM) is a powerful seismic imaging method for the interpretation of steep-dips and subsalt regions; however, implementation of the RTM method is computationally expensive. In this paper, we present a fast and computationally inexpensive implementation of RTM using a NVIDIA general purpose graphic processing unit (GPGPU) powered with Compute Unified Device Architecture (CUDA). To accomplish this, we introduced a random velocity boundary in the source propagation kernel. By creating a random velocity layer at the left, right, and bottom boundaries, the wave fields that encounter the boundary regions are pseudo-randomized. Reflections off the random layers have minimal coherent correlation in the reverse direction. This process eliminates the need to write the wave fields to a disk, which is important when using a GPU because of the limited bandwidth of the PCI-E that is connected to the CPU and GPU. There are four GPU kernels in the code: shot, receiver, modeling, and imaging. The shot and receiver insertion kernels are simple and are computed using a GPU because the wave fields reside in GPUs memory. The modeling kernel is computed using Micikeviciuss tiling method, which uses shared memory to improve bandwidth usage in 2D and 3D finite difference problems. In the imaging kernel, we also use this tiling method. A Tesla C2050 GPU with 4GB memory and 480 stream processing units was used to test the code. The shot and receiver modeling kernel occupancy achieved 85%, and the imaging kernel occupancy was 100%. This means that the code achieved a good level of optimization. A salt model test verified the correct and effective implementation of the GPU RTM code.

Preferential filtering for gravity anomaly separation

We present the preferential filtering method for gravity anomaly separation based on Green equivalent-layer concept and Wiener filter. Compared to the conventional upward continuation and the preferential continuation, the preferential filtering method has the advantage of no requirement of continuation height. The method was tested both on the synthetic gravity data of a model of multiple rectangular prisms and on the real gravity data from a magnetite area in Jilin Province, China. The results show that the preferential filtering method produced better separation of gravity anomaly than both the conventional low-pass filtering and the upward continuation.

GICUDA: A parallel program for 3D correlation imaging of large scale gravity and gravity gradiometry data on graphics processing units with CUDA

The 3D correlation imaging for gravity and gravity gradiometry data provides a rapid approach to the equivalent estimation of objective bodies with different density contrasts in the subsurface. The subsurface is divided into a 3D regular grid, and then a cross correlation between the observed data and the theoretical gravity anomaly due to a point mass source is calculated at each grid node. The resultant correlation coefficients are adopted to describe the equivalent mass distribution in a quantitate probability sense. However, when the size of the survey data is large, it is still computationally expensive. With the advent of the CUDA, GPUs lead to a new path for parallel computing, which have been widely applied in seismic processing, astronomy, molecular dynamics simulation, fluid mechanics and some other fields. We transfer the main time-consuming program of 3D correlation imaging into GPU device, where the program can be executed in a parallel way. The synthetic and real tests have been performed to validate the correctness of our code on NVIDIA GTX 550. The precision evaluation and performance speedup comparison of the CPU and GPU implementations are illustrated with different sizes of gravity data. When the size of grid nodes and observed data sets is 1024x1024x1 and 1024x1024, the speed up can reach to 81.5 for gravity data and 90.7 for gravity vertical gradient data respectively, thus providing the basis for the rapid interpretation of gravity and gravity gradiometry data.

3D correlation imaging of the vertical gradient of gravity data

We present a new 3D correlation imaging approach for vertical gradient of gravity data for deriving a 3D equivalent mass distribution in the subsurface. In this approach, we divide the subsurface space into a 3D regular grid, and then at each grid node calculate a cross correlation between the vertical gradient of the observed gravity data and the theoretical gravity vertical gradient due to a point mass source. The resultant correlation coefficients are used to describe the equivalent mass distribution in a probability sense. We simulate a geological syncline model intruded by a dike and later broken by two vertical faults. The vertical gradient of gravity anomaly of the model is calculated and used to test the approach. The results demonstrate that the equivalent mass distribution derived by the approach reflects the basic geological structures of the model. We also test the approach on the transformed vertical gradient of real Bouguer gravity data from a geothermal survey area in Northern China. The thermal reservoirs are located in the lower portion of the sedimentary basin. From the resultant equivalent mass distribution, we produce the depth distribution of the bottom interface of the basin and predict possible hidden faults present in the basin.

A correlation-based approach for determining the threshold value of singular value decomposition filtering for potential field data denoising

We present a correlation coefficient analysis (CCA) method for obtaining threshold when using singular value decomposition (SVD) filtering method to reduce noise in potential field data. Before computation of correlation coefficients, SVD is performed on the gridded potential field data with the purpose of obtaining singular values of the data. A sliding window is utilized to truncate the acquired singular values, which allows us to obtain different singular value sequences. The lower limit of the sliding window is generally set to zero and the upper limit of the sliding window is the threshold. Then, we calculate and plot the correlation coefficients associated with the initial sequence and the newly obtained sequences, choosing the inflection point of the plotted correlation coefficients as the threshold. The CCA method offers a quantitative way to determine a threshold, which can be easily implemented by a computer program. We illustrate the method using synthetic datasets and field data from a metallic deposit area in the middle-lower reaches of the Yangtze River in China. The results show that the proposed method is effective and is able to provide an optimal threshold.

3D correlation imaging of magnetic total field anomaly and its vertical gradient

We present a new 3D correlation imaging approach for magnetic total field anomaly and its vertical gradient for deriving a 3D equivalent magnetic dipole distribution in the subsurface. In this approach, we divide the subsurface space into a 3D regular grid, and then at each grid node calculate cross correlation between the observed magnetic total field anomaly (or its vertical gradient) and the theoretical magnetic total field anomaly (or its vertical gradient) due to a magnetic dipole. The resultant correlation coefficients are used to describe the equivalent magnetic dipoles distribution in a probabilistic sense. The approach was tested both on the synthetic magnetic data of a model of multiple rectangular prisms and on the real aeromagnetic data from an iron-ore deposit area in the middle–lower reaches of the Yangtze River, China. The results show that the equivalent magnetic dipole distribution derived by the approach basically reflects the subsurface magnetic sources and also illustrate that the approach for the vertical gradient produces a higher resolution of the equivalent magnetic source distribution than that for magnetic total field anomaly alone.

Global correlation imaging of magnetic total field gradients

Firstly we introduce the correlation imaging approach for the x-, y- and z-gradients of a magnetic total field anomaly for deriving the distribution of equivalent magnetic sources of the subsurface. In this approach, the subsurface space is divided into a regular grid, and then a correlation coefficient function is computed at each grid node, based on the cross-correlation between the x-gradient (or y-gradient or z-gradient) of the observed magnetic total field anomaly and the x-gradient (or y-gradient or z-gradient) of the theoretical magnetic total field anomaly due to a magnetic dipole. The resultant correlation coefficient is used to describe the probability of a magnetic dipole occurring at the node. We then define a global correlation coefficient function for comprehensively delineating the probability of an occurrence of a magnetic dipole, which takes, at each node, the maximum positive value of the corresponding correlation coefficient function of the three gradients. We finally test the approach both on synthetic data and real data from a metallic deposit area in the middle-lower reaches of the Yangtze River, China.

Numerical simulation of the electrical properties of shale gas reservoir rock based on digital core

In this paper we study the electrical properties of shale gas reservoir rock by applying the finite element method to digital cores which are built based on an advanced Markov Chain Monte Carlo method and a combination workflow. Study shows that the shale gas reservoir rock has strong anisotropic electrical conductivity because the conductivity is significantly different in both horizontal and vertical directions. The Archie formula is not suitable for application in shale reservoirs. The formation resistivity decreases in two cases; namely (a) with the increase of clay mineral content and the cation exchange capacity of clay, and (b) with the increase of pyrite content. The formation resistivity is not sensitive to the solid organic matter but to the clay and gas in the pores.

Three-dimensional density interface inversion of gravity anomalies in the spectral domain

Based on the Fourier transform, the Parker–Oldenburg algorithm in the frequency domain was extended for the three-dimensional case where the density changes with depth. From this, a gravity interface inversion formula was derived in which the assumed density can be varied laterally and vertically. Iterative convergence is assured by fixing a particular depth as the datum plane below the surface to reduce the interface fluctuation. The results of an example set of synthetic gravity data indicate that the proposed method gives high precision and rapid convergence, with high practical value for the inversion of density interfaces. This method was also used to determine the Moho depth beneath northern China. The results were confirmed by seismic sounding data. Differences between seismic sounding data and inverted depth were insignificant and were in the range of −0.92–1.67 km.

A computationally efficient scheme for the inversion of large scale potential field data

Three dimensional (3D) inversion of potential field data from large scale surveys attempts to recover density or magnetic susceptibility distribution in the subspace for geological interpretation. It is computationally challenging and is not feasible on desktop computers. We propose an integrated scheme to address this problem. We adopt adaptive sampling to compress the dataset, and the cross curve of the data compression ratio and correlation coefficient between the initial and sampled data is used to choose the damping factor for adaptive sampling. Then, the conventional inversion algorithm in model space is transformed to data space, using the identity relationship between different matrices, which greatly reduces the memory requirement. Finally, parallel computation is employed to accelerate calculation of the kernel function. We use the conjugate gradient method to minimize the objective function and a damping factor is introduced to stabilize the iterative process. A wide variety of constraint options are also considered, such as depth weighing, sparseness, and boundary limits. We design a synthetic magnetic model with three prismatic susceptibility causative bodies to demonstrate the effectiveness of the proposed scheme. Tests on synthetic data show that the proposed scheme provides significant reduction in memory and time consumption, and the inversion result is reliable. These advantages hold true for practical field magnetic data from the Hawsons mining area in Australia, verifying the effectiveness of the proposed scheme. An integrated scheme for the inversion of large scale potential field data is proposed.The scheme greatly reduces memory and time consumption compared with a conventional inversion algorithm.Parallelization computation is incorporated to accelerate the inversion process.Aero magnetic data from the Hawsons prospect in Australia is inverted to recover the subsurface magnetic susceptibility.