Exploration Geophysics - Volume 49, Issue 5, 2018

Measurements of natural time-fluctuating magnetic and electric fields at the surface of Earth produce magnetotelluric data, most of which can be expressed as quantities invariant to the directions of field measurement. Such invariants give information on geological structure.

A decomposition of the magnetotelluric tensor is described in terms of quantities which are invariant to the rotation of observing axes, and which also are distinct measures of the 1D, 2D or 3D characteristics of the tensor and so may be useful in dimensionality analysis. When the in-phase and quadrature parts of the tensor are analysed separately there are two invariants which gauge 1D structure, two invariants which gauge 2D structure, and three invariants which gauge 3D structure. A matrix method similar to singular value decomposition is used to determine many of the invariants, and their display is then possible on Mohr diagrams. A particular set of invariants proposed some seventeen years ago is revised to yield an improved set. Several possibilities for the seventh invariant are canvassed, and illustrated by examples from field data. Low values of Δβ, the invariant now preferred for ‘the 7th’, may indicate a particular simplification of otherwise complicated three-dimensional structure.

This paper shows that a method used to identify a peak of horizontal-to-vertical spectral ratios (HVSRs) of microtremors, which was proposed in 2001, is unreliable, and clarifies the reason for the unreliability. Because of the resemblance between the distribution of normalised HVSRs and the F distribution, we frequently fail to detect the deviation between them.

The Albarello test is a statistical test proposed in 2001 with the aim of quickly measuring microtremors (e.g. 5 min) to estimate horizontal-to-vertical spectral ratios (HVSRs). This test is used to identify the deviation between the empirical distribution of normalised (amplitude-corrected) HVSRs and the F distribution so as to identify a physically meaningful peak from numerous peaks, including local peaks inevitably generated by chance when a short-time record is analysed. A HVSR peak is regarded to be meaningful when a significant deviation is detected between the two distributions. However, a paper published in 2005 empirically checked the reliability of this test itself and concluded that the test was unreliable. In this study, we investigate the cause of the unreliability and evaluate the applicable condition. First, we show theoretically that the empirical distribution of normalised HVSRs is generally expected to resemble the F distribution. Then, we show that these two distributions are very similar to each other on the basis of both simulated and observed microtremors. This means that we tend to fail to identify the deviation (commit a type II error) when the sample size is not adequate. To be specific, a sample size on the order of 10²–10⁴ is required for a Kolmogorov–Smirnov test to avoid a type II error. Such a sample size corresponds to records on the order of 10⁰–10² h in length for a representative parameter set for the spectral analyses. The original purpose of this test cannot be achieved because such long records are required to identify the deviation. It may be useful from a practical perspective to consider that normalised HVSRs can approximate the F distribution when developing methods to evaluate the statistical property of HVSRs.

Microseismic monitoring is used to optimise shale gas production or enhanced geothermal stimulation. The technical tools for microseismic monitoring, which is a passive seismic method, are similar to those used in earthquake detection, but differ in that the target area is much smaller than areas affected by earthquakes. Therefore, it is important to use an accurate velocity model. However, such models require conducting an additional survey, which can be both expensive and time-consuming. Many microseismic monitoring studies have used an approximated velocity model constructed from well logging data to reduce these additional costs. In this study, we used a simple approximated model in which velocity increases linearly with depth and creates an accurate velocity model, eliminating the need for an additional survey. We analytically derived formulas for seismic ray traveltime and inverted the velocity gradient using the Gauss–Newton method. Using a numerical example, we verified that the proposed algorithm accurately describes the long-wavelength trend of the true velocity model in a negligibly short time. We performed a Monte Carlo simulation to evaluate the effects of traveltime picking errors. The simulation results indicated that the proposed algorithm provides a reasonable solution under the probable uncertainty of traveltime picking. Finally, we verified that our algorithm was not sensitive to the initial velocity gradient through inversion tests using various initial values. Thus, the numerical example and analysis confirm that the proposed algorithm is efficient and robust.

In this study, we suggest a simple algorithm to estimate 1D velocity gradient for microseismic monitoring. The proposed algorithm is based on the analytically-derived ray formulas for a linearly increasing velocity model and the Gauss–Newton method. Numerical examples show that the proposed algorithm is robust to picking errors and initial guess.

This paper describes the development of a deterministic prestack inversion for estimating elastic and petrophysical parameters. The Gassmann equation is used to construct the relationship between the seismic data and petrophysical parameters. Seismic facies constraints were introduced to improve the accuracy. The very fast-simulated annealing method is used to quickly find the optimal solutions.

Comprehensive utilisation of elastic and petrophysical parameters to predict favourable reservoirs can help reduce the possibility of misidentification of resources. Existing simultaneous inversion methods for estimating the elastic and petrophysical parameters are typically based on either the Gassmann equation, with which these parameters are inverted from prestack seismic data through stochastic optimisation methods, or Wyllie’s modified equation, with which these parameters are inverted from poststack seismic data using deterministic optimisation methods. The purpose of this work is to develop a strategy for estimating the elastic and petrophysical parameters based on the Gassmann equation using deterministic prestack inversion. We employ the Gassmann equation to construct the relationship between the prestack seismic data and petrophysical parameters. We treat the joint posterior probability of elastic and petrophysical parameters as the objective function under a Bayesian framework. Given the macroscopic geological background and the poor-quality prestack seismic data, seismic facies regularisation constraints were introduced to improve the robustness and accuracy of the inversion. The very fast-simulated annealing method is used to quickly find the optimal solutions for the elastic and petrophysical parameters. Based on a model test and the application of real data demonstrates that the proposed inversion method has high accuracy and strong reliability.

This paper examines the capability of a multi-static ground penetrating radar (GPR) system ‘Yakumo’ to estimate a vertical profile of dielectric constants of stratified geological layers. Yakumo is equipped with eight transmitting and eight receiving antennas. By using Yakumo, we propose to obtain an eight-channel common mid-point (CMP) gather from the 64-channel multi-static radar dataset acquired at each point on the survey line. From each CMP gather, a vertical profile of interval velocities can be estimated by velocity spectrum analysis. By a laboratory sandbox experiment, we validate that the thickness of an 83 cm sand layer is estimated with an error of less than 4%. A field survey carried out on tsunami deposits demonstrates that our GPR method can provide the spatial distribution of layer velocities of the sedimentary structure and that the results are comparable with the geological logs at three points and a conventional common-offset GPR profile.

Velocity analysis is applied to explore the capability of a multi-static ground penetrating radar (GPR) system to estimate a vertical profile of dielectric constants of stratified tsunami layers. Both laboratory and field experiments demonstrate that multi-static GPR can accurately delineate the layered geological structures.

The Marine Vibrator Joint Industry Project (MVJIP), sponsored by Shell, Total and ExxonMobil, has been ongoing for four years. During this time significant progress has been made towards developing fully commercial marine vibrators. We begin by describing the history of, and need for, marine vibrator development, including the successes and shortcomings of the technology development before the MVJIP. Following this, we discuss the motivation for initiating the MVJIP, including the specifications required and initial phases of the project. Next, we move on to the vendors that are developing the different options. For each of the three vendors, we shall discuss their background with the technology, the type of technology being developed and their results to date. After discussing the vendors, we shall move on to the impact the MVJIP has had on the commercial marine vibrator industry. Finally, we will discuss our experiences to date and conclude with further work and the way forward for the project.

The Marine Vibrator Joint Industry Project (MVJIP), sponsored by Shell, Total and ExxonMobil, has been ongoing for four years. This paper discusses the history, motivation, technology and latest developments from the project.

We formulate the amplitude versus angle (AVA) inversion in terms of a Markov Chain Monte Carlo (MCMC) algorithm and apply it for reservoir characterisation and litho-fluid facies prediction in offshore Nile Delta. A linear empirical rock physics model is used to link the petrophysical properties (porosity, water saturation and shaliness) to the elastic attributes (P-wave velocity, S-wave velocity and density), whereas the exact Zoeppritz equations are used to convert the elastic properties into AVA responses. The exact Zoeppritz equations allow us to take advantage of the long offset seismic acquisition and thus to consider a wide range of incidence angles in the inversion. The proposed algorithm reliably estimates the non-uniqueness of the solution that is the uncertainties affecting the estimated subsurface characteristics (both in terms of litho-fluid facies and petrophysical properties), taking into consideration the uncertainties in the prior information, the uncertainties in the estimated rock-physics model and the errors affecting the observed AVA responses. A blind test, based on available well log information, demonstrates the applicability of the proposed method and the reliability of the results. In addition, comparisons between the results provided by the implemented MCMC algorithm with those yielded by a linear AVA inversion and an analytical approach to facies prediction, show the benefits introduced by wide-angle reflections in better constraining the inverted parameters and in attenuating the noise in the predicted subsurface models.

We formulate the amplitude versus angle (AVA) inversion in terms of a Markov Chain Monte Carlo algorithm and apply it for reservoir characterisation and litho-fluid facies prediction in offshore Nile Delta. A blind test, based on available well log information, demonstrates the applicability of the proposed method and the reliability of the results.

Channel detection plays a significant role in seismic interpretation. Shearlet transform as a multi-scale and multi-directional transformation is capable of detecting anisotropic singularities. We applied complex-valued shearlet edge measure to synthetic and real seismic time-slices from the South Caspian Sea. The proposed algorithm outperformed both Sobel and Canny edge detectors.

Channels are important sedimentary features in hydrocarbon plays either as targets for drilling or geohazards that should be avoided, depending on burial depth and fluid-fill. Either way, for well design purposes it is important to image channels before drilling. Shearlet transform, as a multi-scale and multi-directional transformation, is capable of detecting anisotropic singularities in two and higher dimensional data. In this study, the complex-valued shearlet-based edge measure was implemented for the aim of channel boundary detection. The method was applied to synthetic seismic time-slices containing channels with different signal-to-noise ratios as well as a real time-slice from the South Caspian Sea. The performance of the shearlet-based algorithm was compared both qualitatively and quantitatively with well known gradient-based edge detectors such as Sobel and Canny, resulting in successfully localising edges and detecting less false positives.

We present a compensation model of aircraft effects for full magnetic gradient tensor data acquired in a towed bird. Training flight tests were designed so that aircraft compensation parameters could be estimated. The feasibility of the compensation method was verified by modelling and by a flight simulation test.

The effect of magnetic interference from the helicopter on full gradient tensor measurements acquired in a towed bird is substantial and must be corrected if airborne data are to be usable. During the actual flight process, the helicopter’s and bird’s attitudes, as well as the relative position between the helicopter and bird, change continuously. Thus the traditional method of compensation for aircraft effects is not suitable for this mode of measurement. For a particularly long towline, magnetic interference from helicopter can be reduced, but the errors introduced by substantial variations in altitude and orientation of the bird may be even greater than the magnetic interference from the helicopter. We have developed a compensation model from the perspective of forward modelling to correct full gradient tensor data measured in a towed bird for the magnetic effects of the helicopter, taking into account variations in attitude of the bird and helicopter. We designed training flight projects that allow compensation parameters to be estimated. The feasibility of the compensation method was verified by modelling and simulated flight tests. Finally, this method was applied to real data and the data quality was improved.

An airborne fluxgate magnetic tensor gradiometer is built on fluxgates to measure directional derivatives of the magnetic field. It has been used to carry out many geophysical exploration programs quickly and efficiently. However, two key issues greatly reduce the data quality of a tensor gradiometer. One is that the fluxgate magnetic tensor gradiometer suffers various errors, such as scale drift, non-orthogonality, misalignment, zero offset, dynamic, and nonlinear errors of individual fluxgates, along with differences between characteristics of fluxgates. The other is that manoeuvring an aircraft flying in the geomagnetic field can generate magnetic interference effects on a tensor gradiometer. Regarding the common airborne fluxgate magnetic tensor gradiometer that has a cross-shaped structure, we have proposed the magnetic interference model of the aircraft and the error model of a single fluxgate. Then we have seamlessly combined these two models into the unified calibration model of a tensor gradiometer by a recursive method. Finally, we have simultaneously determined the correction coefficients and the magnetic properties of the aircraft by a calibration flight at high altitude in an area of low magnetic gradient. We have evaluated the performance of the proposed method through simulation and actual flight results using a microlight aircraft. The root-mean-square noise of each component has reached the level of less than 1.5 nT/m, and the improvement ratios are from 4096 to 17444 in terms of the measured field data of tensor components. The proposed method reduces the reliance of the installation on the aircraft and can easily be applied to other tensor gradiometers, such as airborne superconducting magnetic tensor gradiometers.

In this paper, the correction coefficients of a fluxgate magnetic tensor gradiometer and the magnetic properties of an aircraft are determined. A recursive method is used to combine the magnetic interference model and the error model into a unified calibration model. The method is significant as there will be a greater use of airborne magnetic tensor gradiometers using a wide range of aircraft.

Gravity data processing (reduction) generally utilises the best-available formulae. New and improved formulae have been introduced over time and the resulting newly processed gravity will not match the old. Additionally, mistakes made in the gravity reduction process, as well as the incompatibilities between various equations, will inevitably lead to errors in the final product. This can mean that overlapping gravity surveys are often incompatible, leading to incorrect geological interpretations. In this paper I demonstrate the magnitude of change that results when different information is introduced at various stages of the gravity reduction process. I have focussed on differences relating to calibration factors, time zones and time changes, height, geodetic datums, gravity datums and the equations involved therein. The differences range from below the level of detection (0.01 mGal) to over 16.0 mGal.

The results not only highlight the need to be diligent and thorough in processing gravity data, but also how it is necessary to document the steps taken when processing data. Without proper documentation, gravity surveys cannot be reprocessed should an error be identified.

This paper demonstrates the magnitude of change that results when different information (relating to calibration factors, time zones and time changes, height, geodetic datums and gravity datums) is introduced at various stages of the gravity reduction process. The differences range from below 0.01 mGal to over 16.0 mGal.

Laser Doppler interferometers were previously employed to detect ultrasonic waves propagating in different directions and to estimate elastic anisotropy from these measurements. Our numerical simulations and laboratory measurements show that the recorded wavefield contains converted PS-waves, which need to be taken into consideration to obtain robust estimates of anisotropy.

Ultrasonic measurements using laser Doppler interferometry (LDI) have been reported to provide robust estimates of elastic anisotropy of rock samples. In this approach, an ultrasonic wave is emitted by a piezo-electric source and detected by the LDI, which can be configured to measure three components of the particle velocity in a very small area (~1 mm²) of the sample. Repeating these measurements for a dense array of points on the sample’s surface gives a distribution of traveltimes and polarisation fields on the surface. Anisotropy is then obtained by inverting these fields using analytical expressions or numerical algorithms for computing phase and group velocities. The existing implementation of this approach involves the inversion of direct compressional (P) and shear (S) wave arrivals only. A previous study showed that this approach produces stable results if only a small range of source–receiver offsets is included in the inversion. This limitation resulted in a relatively large uncertainty of the result. This uncertainty can be reduced by inverting the entire traveltime field. To this end, we numerically simulate the wavefield in the sample. Analysis of the computed wavefield reveals the presence of P- and S-waves as well as a critically refracted converted PS-wave. Hence, the inversion of the entire traveltime field must include these three waves. We implement this inversion using global minimisation of the traveltime misfit function, coupled with numerical computation of ray velocities. Application of this algorithm to laboratory LDI measurements on a transversely isotropic phenolic sample provides stable anisotropy estimates consistent with previous studies.

The ancient Roman city of Sagalassos (Turkey) is covered by layers of eroded soil that have preserved many secrets waiting to be revealed. A geophysical campaign was planned to highlight the buried structure. Geophysics revealed evidence of a clay quarry, a number of tombs related to the Byzantine period and defensive walls.

The archaeological site of Sagalassos is a very important settlement located in a magnificent mountain landscape, 7 km north from a village named Ağlasun (province of Burdur, south-west Turkey). Since 1990, the University of Leuven (Belgium) has carried out an interdisciplinary archaeological research program that studies >1000 years of uninterrupted human occupation in Sagalassos, concerning all historical aspects of daily life from architecture, to trade and its mechanisms and environmental conditions. The ancient Roman city is covered by layers of eroded soil that has preserved many secrets waiting to be revealed. A geophysical campaign was planned along the south facing terraces of the mountain slopes to highlight the structure of the city that remains covered in soil. Site conditions (high slope, high grass, several obstacles) and the need to investigate to depths greater than 20 m influenced the choice of geophysical methods; we chose to use both passive and active electrical resistivity tomography. Three different areas, labelled Area 1, Area 2 and Area 3, were investigated, with results revealing information about the location, depth, size and extent of buried archaeological features. Of particular interest is the presence of: (i) a deep depression in Area 1, thought to be a clay quarry; (ii) a number of tombs related to the Byzantine period in Area 2; and (iii) defensive walls in Area 3.

This paper investigates the use of induced polarisation (IP) and resistivity methods to characterise bitumen deposits in western Iran. Specifically we explore the Barreh-palang bitumen deposit in Gilan-gharb, Kermanshah province, Iran. These deposits are mainly found in vertical structures and, in some cases, complicated fault and fracture structures. Additionally, the hydrogeophysical environment adds complexity due to elevated groundwater levels that are of moderate to high salinity. Both resistivity and IP data were recorded using a combined resistivity sounding and profiling (CRSP) array with potential electrode spacing of 5 m in four profiles. Several anomalous responses were detected. For several of these responses, the bitumen-rich zone was difficult to distinguish from other responses without the IP data, implying that resistivity alone is not always sufficient for detection of bitumen, especially in complex geological situations. Interestingly, we also show that the presence of saline water was useful to assist in detecting the target in some situations.

A hydrogeophysically and geologically complex bitumen deposit in western Iran was investigated by means of resistivity and induced polarisation (IP) methods. The results show the high potential of IP data to improve bitumen detection. It was also proven that saline water can sometimes be useful in target identification in mineral explorations.

In laboratories, core cannot always be cut exactly along the axis of symmetry (normal to the bedding plane), which leads to a minor core cutting rotation error. This paper investigates the sensitivity of Thomson’s anisotropic parameters, epsilon (ε), gamma (γ) and delta (δ), to the cutting error.

In laboratories, core cannot always be cut exactly along the axis of symmetry (normal to the bedding plane), which leads to a minor core cutting rotation (CCR) error. The presence of a small CCR error can give rise to an error in computation of anisotropic parameters, which will result in erroneous P-wave and S-wave velocities (VP, VSV and VSH). In this study, we test the sensitivity of Thomson’s anisotropic parameters, epsilon (ε), gamma (γ) and delta (δ), to the CCR error. In the vertical transverse isotropy (VTI) system where no rotation exists, values of ε, γ and δ will not vary with azimuth. Similarly, in VTI the values of VP, VSV and VSH, measured at orientations with respect to the axis of symmetry, will also not vary with azimuth. We modelled and analysed the error generated, in anisotropic parameters and phase velocities, due to rotation error of 5° on VTI ANNIE model data, and the Niobrara Shale core sample measurements. For these two cases, we obtained values of P-wave and S-wave velocities along with anisotropy parameters ε, γ and δ, at 0°, 45° and 90° orientations with respect to the axis of symmetry, before and after a rotation of 5°, for different azimuth directions. This study shows that, when the CCR error exists, ε, γ and δ vary with azimuth. For Niobrara Shale samples, we observed that, among all three Thomsen’s parameters δ is the most sensitive parameter to the CCR error; when we varied azimuth, we observed a sign change of δ from positive to negative. For ε and γ, variation in azimuth lead to only slight changes in values without any change of sign. The CCR error affects VP and Vs measurements the most at 45° orientation, and the least at 90° orientation. The maximum error occurred at 0° azimuth, whereas the minimum error occurs at 90° azimuth. This analysis suggests that, if core cutting rotation exists, phase velocities should be measured at 90° azimuth for accurate results.

Distinguishing high-quality shale effectively is a difficult exploration problem. We propose a constrained spectral inversion method based on compressive sensing to handle this challenge. This method can help us improve the resolution and continuity of the profiles and enhance our ability to discover high-quality shale.

Multi-wave exploration has become one of the main means of unconventional shale gas exploration in the Sichuan Basin, China. How to effectively improve the resolution of the shale gas reservoir layer and distinguish high-quality shale has become one of the difficulties for shale gas exploration. The spectral inversion technology can effectively break the resolution limit of the conventional technology and greatly improve the resolution of the thin shale layer; however, it also needs enough prior information to overcome the problem of multiple solutions in the inversion. The compressed sensing (CS) theory has the ability to reconstruct complete data using incomplete data can use less data information to improve the accuracy of spectral inversion. A novel constrained spectral inversion method based on CS is presented to handle these situations. The CS technique is applied to the objective function of spectral inversion, which can improve the accuracy of the spectral inversion algorithm and create profiles with a higher resolution and greater continuity. Applications through theoretical and real data can illustrate very high performance of the presented algorithm.