Inversion of Gravity Anomalies Using Primal-Dual Interior Point Methods

Structural inversion of gravity datasets based on the use of density anomalies to derive robust images of the subsurface (delineating lithology and their boundaries) constitutes a fundamental non-invasive tool for geological exploration. The use of gravity data to estimate and interpret the substructure based on its density properties have proven efficient; however, the inherent non-uniqueness associated with most non-seismic geophysical datasets make this the ideal scenario for the use of recently developed robust constrained optimization techniques. We present a constrained optimization approach for a least squares inversion problem to characterize 2-Dimensional Earth density structure models based on gravity anomalies. The formulation inverts Bouguer gravity anomalies using polygons along with Primal-Dual Interior-Point methods for the optimization of the results, which include equality and inequality physical and structural constraints. We validate our results using synthetic density crustal structure models with varying complexity and illustrate the behavior of the algorithm using different initial density structure models and increasing noise levels in the observations. Moreover, we apply the approach to the Southern Rio Grande Rift (SRGR) region using previously obtained receiver function results as constraints for the inverted density profiles. We produce constrained crustal models that characterize the SRGR showing a shallower Moho (30 km) in the region, which is thicker than previously suggested. Based on its validation and implementation, we conclude that the algorithm using Primal-Dual Interior-Point methods is robust and always honors the geophysical constraints. Advantages of using this approach for structural inversion of gravity data are the incorporation of a priori information related to the model parameters (coming from actual physical properties of the subsurface) and the reduction of the solution space contingent on these boundary conditions.

Nowadays techniques and technologies used in geophysical exploration focused on finding additional natural resources besides those of fossil origin, has new improvements based on theories applied so far only in researching field, such as spectral, and geostatistical analysis in order to estimate the origin of basal magnetic sources known as the depth of Curie´s isotherm; very important issue on geothermal deposits.
Fourier analysis applied to potential fields allows the depth estimation of interfaces into the Earth's crust with high contrasts of magnetic susceptibility or density. This method, applied to the magnetic field, uses the relationship among the magnetic anomalies' power spectrum, depth and size of such sources. (Spector and Grant, 1970;Bhattacharyya and Leu, 1975;Blakely, 1995 andRuiz andIntrocaso, 2004).
From spectral and geostatistical analysis proposed by Okubo and Bhattacharyya (1985;1970) it is possible to estimate the depth of basal magnetic source associated with of Curie's isotherm depth whose temperature ranges between of 580 and 600 °C. (Blakely, 1988(Blakely, , 1995Tanaka et Al., 1999 andFrost andShive, 1986).
This study is based on the aeromagnetic data acquired by the Mexican Geological Survey and processed to estimate the depths of the magnetic sources in the Tres Vírgenes area, Baja California Sur, Mexico.
The area of interest was divided in seven windows with squares of 64 Kilometers; knowing that the maximum depth that can be estimated depends on the length of the 2π window. Main magnetic field was included in each window.
In order to estimate depths, the Spector and Grant (1970) methodology was applied; such process assumes that a set of vertical prisms are the sources of the anomalies on the aeromagnetic charts. These prisms define a distribution of magnetic sources, computing its radial power spectra, obtaining its square root, and its natural logarithm it is possible to associate its second slope of the resulting function with the magnetic structure´s top depth defined as Zt.
On the other hand, the geostatistical analysis was done following the methodology published by Bhattacharyya, Leu (1975, 1977 and Okubo (1985) where the radial power spectrum rise to the square root and divided by the absolute value of the wave number. After that the natural logarithm was obtained in order to find the centroid depth (Z0) of the magnetic structure. The basal depth of the magnetic structures defined as Zb was calculated by the ratio (2Z0 -Zt ).
The average depth in the area of interest in Tres Virgenes; ranges between 5 to 6 kilometers. It is possible to correlate such results with heat flow values and estimate the current potential of the Tres Virgenes geothermal field and its surroundings.
Curie´s isotherm, Tres Virgenes, radial power spectrum STT55-02 In the fields of earthquake disaster prevention and energy resource exploration, understanding the properties of seismic anisotropy is important for obtaining various information about subsurface, for example, its structure, regional stress field and selective orientation of crack opening direction. In some cases, like near surface and reservoir rocks, formations are estimated to be strongly anisotropic for elastic wave propagation. To extract anisotropic information from seismic exploration data, many researches have conducted elastic anisotropy studies. However, most of these studies are based on the assumption of weakly anisotropic media (Thomsen, 1986). Our previous studies showed that strong anisotropic media in the subsurface significantly influence the seismic waveforms especially on the PS converted waves. In the present study, we apply the reverse time migration (RTM) to the PS converted waves to determine the depth of anisotropic layer. To extract PS converted waves from observed data, we also develop a novel wave separation method. We demonstrate the effectiveness of our method using a numerical experiment.
Our numerical result shows that our method can image layer boundary between isotropic and anisotropic layers which generates strong PS converted waves.

地震波速度異方性、PS変換波
Seismic anisotropy, PS converted waves STT55-03 Finite difference is the most widely used method for seismic wavefield modelling. However, most finite-difference implementations discretise the Earth model over a fixed grid interval. This can lead to irregular model geometries being represented by 'staircase' discretisations, and potentially causes mispositioning of interfaces within the media. This misrepresentation is a major disadvantage to finite difference methods, especially if there exist strong and sharp contrasts in the physical properties along an interface. The discretisation of undulated seabed bathymetry is a common example of such misrepresentation of the physical properties in finite-difference grids, as the seabed is often a particularly sharp interface owing to the rapid and considerable change in material properties between fluid seawater and solid rock. There are two issues typically involved with seabed modelling using finite difference methods: firstly, the travel times of reflections from the sea bottom are inaccurate as a consequence of its spatial mispositioning; secondly, artificial diffractions are generated by the staircase representation of dipping seabed bathymetry. In this paper, we propose a new method that provides a solution to these two issues by positioning sharp interfaces at fractional grid locations. To achieve this, the velocity model is first sampled in a model grid that allows the centre of the seabed to be positioned at grid points, before being interpolated vertically onto a regular modelling grid using the windowed sinc function. This procedure allows for undulated seabed bathymetry to be represented with improved accuracy during modelling. Numerical tests demonstrate that this method generates reflections with accurate travel times and effectively suppresses artificial diffractions.
finite difference, seabed bathymetry, integer grid, fractional grid, staircase, sinc interpolation In the present study, we use a stretched stencil of FD operator not to increase the bandwidth in the impedance matrix. Coefficients of the stencil are determined by a minimization process. We investigate the accuracy of our scheme using dispersion analysis and numerical experiment. They show that the proposed scheme can improve not only accuracy but also efficiency compared to the conventional 9-point scheme proposed by Jo et al. (1996). Seismic inversion is a highly nonlinear and ill-posed problem. In a typical seismic survey, aperture size, source frequency bandwidth, and source and receiver distribution are limited. Therefore, the inversion results usually suffer from strong artifacts and reduced resolution.
In this study, we propose a novel full waveform inversion (FWI) method with a nonlocal similarity prior and adaptive sparsity-promoting regularization in the gradient domain. First, we invoke a self-similarity between patches extracted from a single velocity model. In other words, these patches have a certain level of redundancy. To exploit this property, we adopt a regularization term that finds and stacks together the k most similar patches into a 3D cube for each patch from the velocity model, and then require that the cube be sparse under a 3D orthogonal sparsifying transform. A second regularization term is a generalization of total variation with learning-based dictionaries, which provide a good estimation of the gradient field. Since the gradients emphasize high-frequency components of the velocity model, a better reconstruction of the gradients means a more accurate reconstruction of high spatial frequencies, hence higher resolution.
The proposed optimization problem can be solved by the Alternating Direction Multiplier Method  to monitor the seismic activity and the process of earthquake generation including the stress accumulation.
To elucidate earthquake generation and preparation process, it is necessary to investigate how the stress could be accumulated not only in deeper part but also in the shallow sediments, what the role of interstitial fluid could be in the stress accumulation processes, etc. There are some conventional methods to measure these physical properties, such as borehole strainmeter, borehole breakouts or borehole dynamic tests. However, these methods have some difficulties from the viewpoints of technical and/or cost. Therefore, we need to have some other methods to see the state and time variation of the stress in the subseafloor.
In this study, we applied seismic interferometry technique to ambient noise records observed by DONET seafloor seismometers to obtain time dependent seismic velocity as a proxy of stress state below seafloor.
We calculated zero-offset horizontal and vertical pseudo shot records from every 1 hour ambient noise records in two years continuous data since July 2014, which were observed by three components of DONET seismometers. Obtained pseudo shot records are then stacked every 240 hours (10 days). Clear events are visible in both vertical and horizontal components of the pseudo shot records, which could be reflected P-and S-waves from the bottom of the shallow sediment layer, respectively. Then we applied the stretching interpolation method to the obtained pseudo shot records. Results confirmed that time variation of P-and S-wave velocity in the shallow sediment layers are stable with its velocity changes of less than 1 percent without large seismic events.
On 1 April 2016, a large offshore earthquake, the off Mie earthquake (Mw = 6.0), occurred. DONET KME18 seismometer is located in the vicinity of the source region with horizontal distance of less than 10 km, and clearly observed co-and post-seismic events including aftershocks and tremors. In pseudo shot records obtained from the KME18 seismometer, large S-wave velocity change of 5 percent maximum, was observed as co-and post-seismic velocity change with duration of approximately 20 days, although velocity change of P-wave was less than 1 percent in the same duration. Results confirmed that change of pore fluid migration and crack distributions, which could be caused by strong motion of the large earthquake, are observable as a proxy of stress distribution and changes below each seismometer by using ambient noise records.
We plan to perform further analysis and discussion including S-wave anisotropy analysis and more quantitative discussion of the relationship between seismic velocity, its anisotropy, pore-fluid pressure and stress for this approach to be used as simple stress monitoring tool to understand the mega-thrust earthquake preparation and generation cycle. In reflection seismic survey, relatively low resolution information can be obtained over a wide range, whereas the borehole logging data contains a high precision data in a narrow range. It is important to estimate the physical properties by integrating these two kinds of information. Recently, anisotropic structure in the subsurface has drawn attention because of its impact on the seismic survey and interpretation. The velocity of seismic wave varies due to natural fractures, sediment structures and selective crystal structures of subsurface materials. Understanding anisotropic properties is important to estimate dynamic properties such as deformation and stress of rock inside the earth and to estimate the composition of minerals.

全波形逆解析によるS波異方性の検出
In this research, a numerical experiment is conducted on the sonic logging model which contains the anisotropic layer. A cross dipole, which consists of two orthogonal dipole sources, is used as a transmitter.
For the anisotropic model, a transversely isotropic medium with a horizontal axis of symmetry with five independent elastic elements is set. According to previous studies, it is suggested that the received waveform obtained by the cross dipole measurement changes with the degree of anisotropy of the rock.
Distinctive feature can be observed when shear wave passes through the anisotropic layer and splits into two polarized waves: fast shear and slow shear. We attempt to estimate the elastic parameters using the shear waves. Full waveform inversion (FWI) is applied to estimate the elastic parameters directly. FWI is one of the techniques to estimate the ground physical properties with high accuracy and resolution by updating the model parameters so as to minimize the residual between the observed and the calculated waveforms. As a result, it was suggested that highly accurate estimation of elastic parameters is achieved by applying FWI to the layer where anisotropy exists. Although the resolution depends on the receiver interval in conventional slowness time coherence method, it is suggested that FWI can estimate ground physical properties with high resolution.