GRV_D_inv: A graphical user interface for 3D forward and inverse modeling of gravity data

Геофизический журнал No 1, Т. 43, 2021 18

cludes the FFT-based methods based on rearrangement of forward algorithm of [Parker, 1972], e.g. [Oldenburg, 1974;Granser, 1986Granser, , 1987Reamer, Ferguson, 1989;Guspi, 1993;Pilkington, 2006;Zhang et al., 2015;Wu, Lin, 2017]. The methods utilize the relationship between the Fourier transforms of the gravity data and the sum of the Fourier transforms of the powers of the depth to the interface. The main advantage of these methods is the rapid calculation of gravity anomalies that leads them to be computer time-efficient practical applications. However, on the other hand, their constraints in practice are that they require a predetermined average depth of the interface and a high-cut filter to ensure convergence of the inversion procedure [Oldenburg, 1974;Pham, Do, 2017;Pham et al., 2018Pham et al., -2020. Based on Parker-Oldenburg method, a number of computer programs are available, for instance, [Nagendra et al., 1996] developed computer programs in FORTRAN language for the interpretation of 2D gravity data, and [Shin et al., 2006] for the analysis of 3D gravity data. [Gomez-Ortiz, Agarwal, 2005] presented a MATLAB code for inverting the gravity anomaly over a 3D density interface by Parker-Oldenburg's algorithm. The code has been used later by many other researchers for interpreting real data from different parts of the world [Ouyed et al., 2010;Kaya, 2010;Hsieh et al., 2010;Van der Meijde et al., 2013;Tugume et al., 2013;Oruç, 2014;Grigoriadis et al., 2015;Gao et al., 2016;Oruç, Sönmez, 2017;Altinoğlu et al., 2018;Nguiya et al., 2019 and among others]. However, a recent paper published by [Gao, Sun, 2019] showed that [Gomez-Ortiz, Agarwal, 2005] misused Parker's formula, thereby leading to incorrect forward and inversion results.
This paper aimes to present an alternative software as a time-efficient gravity inversion tool to determine the 3D depth structure of a density interface related to the observed gravity anomalies. The algorithm linked to the code GRV_D_inv.m is built in Matlab including a graphical interface, which makes it easy for the user to configure the required presets prior an inversion or forward procedure simply without any need of coding experience. The software has been tested for its practical application and accuracy on both synthetic and actual gravity data.
Theory. Describing the model space by a three-dimensional Cartesian coordinate system where its pair-wise perpendicular axes are x, y and z pointing towards east, north and positive downwards, respectively, the 3D expression of the gravity anomaly in the Fourier domain due to a uniform layer of material is given by [Gao, Sun, 2019]: ( 1) where F-1[ ] symbolizes the inverse Fourier transform. Hereby, Eq.
(2) can be used through an iterative inversion procedure to estimate the depth to the undulating density interface. The procedure begins by a preset model of the depth h, for instance zero for all. With this initial starter model, taking the inverse Fourier transform of the first term of Eq. (2) leads the first approximation of the topography interface h. A new set of depth estimate is then calculated by using the h values from the previous step again in equation Eq. (2). Updating for a new depth estimate h continues until a desired fit between two successive depth estimates is achieved whereas the root mean square error (RMS) can be used as the measure of this goodness ( ) where M and N are the numbers of grid size in north and east, respectively, t stands for the iteration step.
To ensure the convergence of the inverse series in Eq. (2), a low-pass filter B(k) is applied during the calculation. The filter is defined in the following equation [Oldenburg, 1974] ( ) ( ) where SH and WH define the roll-off frequen-cies of the filter design. The filter passes all the frequencies lower than WH, attenuates between WH and SH, and cuts off all the frequencies above the SH frequency. An appropriate choice of the roll-off frequencies can be obtained by spectrum analysis of the observed gravity anomaly. Using high-cut filtering results in loss of information's at high frequencies induced by near-surface structures, hence the gravity anomalies due to the inverted depths do not necessarily fit the highfrequency content of the actual anomaly.
Overview of the GRV_D_inv GUI. The GRV_D_inv code is supported by a functional user interface including graphical control functions that allow the user to customize the inversion/forward modeling configurations and setting of the output export formats. When the first time the code is run, it opens the inversion panel as the main graphical interface covering the entire screen. Fig. 1 illustrates the configuration of this panel structure. It includes graphical control items for loading of gridded gravity data, editable cells for settings of the field and filtering parameters, the iteration stop criterion, and a confirmation button that initiates the inversion procedure. The remaining part of the panel contains the input map and its spectrum display area and the display area of the obtained outputs. The interactive [LOAD DATA] menu item at the top left of the panel is used to import the gridded gravity data set which supports the grid formats (*.grd) of the Surfer program (Golden Software). The mesh grid intervals to the east and north are required to be equal and any blank in input is not allowed. After completion of the data loading, the label of the data-loading item is replaced with the file name of the input, and a contour map and a spectrum graph corresponding to the input data are displayed. Next, the set of the density contrast, the mean depth of the density interface, the roll-off frequency parameters of the filtering, and the stopping criterion of the inversion are required. Distances of the map units are expressed in km, the gravity in mGal and the unit of density contrast in g/cm 3 . The set of the SH and WH roll-off frequency values can be done either by entering them manually into the corresponded edit boxes or interactively by mouse controls on the spectrum graph. The iterative procedure stops when the goodness of the fit between two successive depth estimates is below the predefined RMS threshold or when it completes the preset limit number of iterations.
After validating the entries, the process starts with building the initial depth estimates from Eq. (2) and uses them again in Eq.
(2) for their improvement. This operation continues until the matching criterion between two successive depth estimates is achieved. Finally, the code records the gravity response and the inverted basement depths obtained at the stopping iteration step, the residual between the actual and modelled gravity anomalies and the values of RMS after each iteration. Following, the code enables us to display or to export any of the output either in numeric or in image format by using the menu items at the output panel (Fig. 2). The results are exported with a user-defined file left side of the screen. Synthetic data application. The applicability and effectiveness of the above-described program is demonstrated by the inversions of noise-free and noise corrupted gravity data sets synthetically produced from a simulated 3D density interface. Fig. 3 shows the 3D and plan views of the topography of the interface with an average depth of 20 km. The data set is calculated through Eq. 1 on a 256×256 mesh grid of 1 km intervals with using a density contrast of 0.4 g/cm 3 . In the first case, the program is applied to reconstruct the density interface geometry from the noise-free anomaly. Fig. 4, a displays the noise-free gravity anomaly map due to the model in Fig. 3. For the inversion of this data, the high-cut filter frequency parameters were selected as SH=0.035 km −1 and WH=0.025 km −1 . The iterative process performed twenty-fouriterations to fall below the preset allowable RMS of 10 -4 km between two successive interface approximations. The RMS after the initial stage was 0.0132 km and reduced to 9.8098·10 −5 km at the stopping step (Fig. 4, f). Fig. 4, b shows the geometry of the estimated density interface and Fig. 4, d the residual between these and the actual density interface. Here, the absolute depth differences between the actual model and the re-constructed one are mainly less than 0.1 km; hence, the estimated relief well represents the true model. The RMS between them is 0.0291 km. Fig. 4, c is the calculated gravity anomaly due to the inverted density interface shown in Fig. 4, b. By comparing of Figs. 4, a and 4, c, one can observe that the calculated gravity anomalies are also very close in shape to the original gravity anomalies. The residuals between them (Fig. 4, e) range from −0.2171 to 0.1965 mGal and the RMS is 0.0745 mGal.
The sensitivity of the proposed code to noise-corrupted anomaly observations is analyzed in the second example. In this case, the theoretical observations in Fig. 4, a were contaminated with random noise of amplitude 10 % of the original data. Fig. 5, a shows the noise-corrupted gravity data. For inverting this data, the roll-off frequencies of the high-cut filter have been set as SH=0.035 km −1 and WH=0.025 km −1 . The inversion process iterated until the RMS between two successive interface approximations drops below the threshold set 10 −4 km value. In this case, the RMS error decreased from 0.0143 km at the first iteration to 9.8209·10 −5 km at the end of twenty-fifth iteration (Fig. 5, f). The inverted depths and the modeled gravity anomalies after the termination step are given in Figs. 5, b and 5, c, respectively. Although the input data is noisy, it can be seen that the inverted depths (Fig. 5, b) still matches closely with those of the actual model given in Fig. 3 where the residuals between them are in general less than 0.1 km (Fig. 5, d). Consequently, the gravity response of the estimated depth construct (Fig. 5, b) also correlates well with the original model gravity anomaly. The difference between the original model and the calculated gravity anomaly is shown in Fig. 5, e where the RMS error between them is only 0.1230 mGal.
Field Example. The functional applicability of the proposed code is also demonstrated by interpreting a real gravity anomaly from  Inversion results for noise-corrupted data example, (a) theoretical gravity anomalies including 10 % noise (b) the estimated depths from the code application to the noise-corrupted data, (c) calculated gravity anomalies from the inverted interface, (d) the residuals between the actual and computed depths, (e) the residual between the actual and computed gravity anomalies, (f) plot of RMS errors obtained during the on-going iterative process.
Brittany, France, previously analyzed as well by other researchers for estimating the Moho depths of the region, i.e. [Lefort, Agarwal, 2000;Gomez-Ortiz, Agarwal, 2005]. The lo-cation of the study area is shown in Fig. 6. Fig. 7, a shows the residual gravity anomaly ascribed to the Moho interface in Brittany, which obtained from [Gomez-Ortiz, Agarwal, 2005] program on a 51×51 mesh point with square grid intervals of 4 km. The four input parameters are used for the inversion of this data by the present code. They are the density contrast between the lower crust and the upper mantle as 0.4 g/cm 3 and a mean depth of the density interface as 30 km [Gomez-Ortiz, Agarwal, 2005], the roll-off frequency parameters of the high-cut filter chosen as SH=0.012 km −1 and WH=0.01 km −1 [Lefort, Agarwal, 2000]. Here, the inversion scheme performed seventy-one iterations to fall below an allowable RMS error of 2·10 −4 km between two successive depth approximations (Fig. 8). The estimated depths of the Moho interface and the gravity field generated from these depths are given in Figs. 7, b and 7, c, respectively. As can be seen by comparing Fig. 7, a and 7, c, the calculated anomalies fit very closely with the observations. The differences between them are generally less than 1mGal (Fig. 7, d), and the RMS error is only 0.24 mGal. Further, the closeness of the fit between the calculated and observed anomalies is obtained better in this study when compared to those obtained by [Gomez-Ortiz, Agarwal, 2005] with a maximum difference of 1.5 mGal. The depth estimates from the present study ranges from 27.77 km to 32.38 km (Fig. 7, b). For comparison, Fig. 7, e shows the geometry of the density interface estimated by [Gomez-Ortiz, Agarwal, 2005] and Fig. 7, f the interpreted model by [Lefort, Agarwal, 2000] using an inversion method described by [Tsuboi, 1983]. Although the depths obtained from this study are quite similar to those obtained by [Gomez-Ortiz, Agarwal, 2005], a few deviations can be observed, for example, the deeper parts of the interface observed in the southeast of the area show a NW-SE directional trend, which is also figured in the structure of [Lefort, Agarwal, 2000]. By and large, the result obtained by the application of the code to the gravity data of Brittany is comparable with those reported by [Lefort, Agarwal, 2000].
Conclusions. A GUI-based Matlab program GRV_D_inv is presented to estimate the relief of a density interface from gravity anomalies by an iterative procedure. The algorithm linked to the code operates iteratively in the frequency-domain. Hence, the computational speed is fast and therefore it can be a computer time-efficient practical application in processing large data sets. The inversion of a 256×256 gridded data completed within 24 iterations in a PC equipped with a 2.4 GHz central processing unit (CPU) took only 3.9 s, which can be considered as a short duration for this kind of inverse problem. As a further advantage, the code is supported with a practical user interface that includes interactive control functions that makes it easy for the user to configure the settings prior the inversion as well as to manage the extraction of the outputs with optional preference without requiring any coding expertise.
The present code has been tested for its practical application and accuracy on both noise-free and noise corrupted synthetic gravity data produced from a 3D density interface model, as well as on a real gravity data from Brittany, France, where the results were compared with those of previous studies. The depth results inverted from the noisefree data well matched the real structure. Al- Fig. 6. Location of the study area. Fig. 7. The real field example, (a) the residual gravity anomaly from Brittany, France, (b) the estimated Moho relief from practice of the proposed code, (c) the gravity anomaly calculated from the estimated interface, (d) the difference between the actual and the modeled gravity anomalies, (e) the Moho relief reported by [Gomez-Ortiz, Agarwal, 2005], (f) the Moho relief reported by [Lefort, Agarwal, 2000]. though there are some minor deviations between the actual and estimated depths from the noise-corrupted data sample, the results are still in good agreement even the given gravity anomalies are noisy. The obtained depth structure of the Moho in Brittany from the present study is in close agreement with the previously published structures by other studies. As a result, the proposed code has proved effective in interpreting gravity anomalies to determine the depth structure of a density interface and can therefore be considered as a practical tool for modeling purposes in geophysical applications.
Code availability. The GRV_D_invis a free program and it can be downloaded together with synthetic data files and a user manual from the web site [https://github.com/ PhamLT/GRV_D_inv].