Introduction

In gravimetric interpretation, determining the parameters of disturbing bodies is usually difficult due to the ambiguity of solutions. This ambiguity is due to Gauss' and Stokes' laws, trying to find the parameters of three-dimensional bodies based on the surface or profile distribution. The parameters that you want to achieve are of course the shape of the body, its density (density contrast) and the depth of the position. Due to the fast reduction of the gravitational influence of the body with depth (inversely proportional to the square of the distance), the body shape is usually assumed in the form of simple solid figures, such as a sphere, vertical and horizontal cylinder, or a prism, and depth is looked for. It should also be noted that the linear dependence of the gravitational effect and the density (density contrast) of the body does not affect the shape of the measurement curve, but only its amplitude.

In the case of microgravity studies, the problem is even more complicated, because there are two factors that make interpretation difficult. Firstly, due to the shallow position of the bodies, the approximation of the disturbing bodies with simple solid figures is in many cases insufficient, and secondly, a much higher error to useful signal ratio is observed for microgravity studies. In this publication, the authors assumed that the simple solid figures approximation is sufficient and focused on the impact of the error that was simulated by a random error.

There are many depth estimation techniques in potential methods, and thus also in the gravity method, performed using different methods (Hinze et al. 2012; Cooper 2011). Generally, they can be divided into two types—graphical techniques and semi-automated methods. Of course, there are other depth estimation solutions in potential methods, but they are of less application.

The graphical techniques are rather historical. They are based on the relationships between the distance of the characteristic elements of the anomaly and the depth of the top or centre of the solid figure. The first one is the half-width method (Nettleton 1940). It is based on the relationship between the characteristic points of the theoretical anomaly distribution—horizontal distance between maximum amplitude and half of its value. The second one is the straight-slope method (Ram Babu et al. 1987). The vertical gravity derivative is often used instead of the gravity anomaly because this method is based on the distance at which the gradient of the anomaly or its derivative is constant. The Smith rules (Smith and Bott 1958; Smith 1959) are one of the graphical methods independent of source geometry. It is based on horizontal gravity derivatives.

The most modern depth estimation methods are the semi-automated methods. They operate on the basis of the FFT—fast Fourier transform, especially to calculate vertical and horizontal derivative of gravity anomaly. All results are an approximation of the depth and depend on the window length, quality of anomaly and geological conditions. The window length is the length of the operator window used to calculate depths and can be movable or fixed. If the window is fixed, the best length of this should be close to half the width of the gravity anomaly. The correctness of depth estimation, of course, depending on the quality of gravity anomaly, it means it depends on the measurement error. The most often used semi-automated methods are Euler deconvolution, Werner deconvolution and statistical spectral techniques. Euler deconvolution (Thompson 1982) and now Extended Euler deconvolution (Mushayandebvu et al. 2001) are based on the structural index (SI) for a couple of simple solid figures but can be used for almost all gravity sources. The solution is obtained by the least-squares inversion of the Euler’s homogeneity equation for assumed window. Werner deconvolution determines the depth of thin dike or sheet, on the understanding that the source of anomaly can be described as a complex of those dikes or sheets (Werner 1953). The solution is obtained by standard matrix inversion for six coefficients (Dimri 1992; Ku and Sharp 1983). This deconvolution is also used for 2D contact depth estimation. The statistical spectral technique consists of spectral analysis of the anomalies; it means the wavelengths of the anomaly analysis. This method can be applied both on profiles and on maps. The average depth is the result of this technique (Spector and Grant 1970; Bhattacharyya and Chan 1977).

In the article, the Werner and Extended Euler deconvolution methods were selected for the analysis of depth estimation in microgravity studies. All depth calculations were made in the Geosoft Oasis montaj application, and their further processing was based on our own programs.

Optimal step in gradient calculations

Selected methods of calculating the depth to the disturbing body (Werner and Extended Euler) in their calculations contain horizontal and vertical gradients of gravity. Undoubtedly, the calculation of both gradients is largely influenced by the error in determining the value of gravity at a measurement station which was presented for the horizontal gradient by Porzucek (2010). In this paper, the authors decided that the correctness of the calculation of the depth to the body generating the anomaly would be analysed for theoretical data disturbed by random errors with a maximum error equal to 5% and 10% of the anomaly amplitude. The analysis of many different bodies for both methods showed that the horizontal cylinder would be the best body to analyse.

The horizontal cylinder selected for testing had parameters that could reflect a horizontal gallery or a karst void hollowed by a flowing underground river. The parameters of the test cylinder were: radius 1.7 m, depth 4 m and density contrast   − 2.55 g/cm3. For the parameters adopted in this way, the amplitude of the calculated anomaly was 773 m·s−2, i.e. within the limits of microgravimetric surveys.

Porzucek (2010) noted that the calculation of the horizontal gradient is obviously influenced by the measurement error, but also by the distance between the computational points. In the first stage of the analysis, an attempt was made to determine the optimal calculation step based on the calculation of the horizontal gradient. In order to thoroughly analyse the solutions, an anomaly was generated by taking a calculation step of 0.5 m, and then a random error in the range of ± 5% and ± 10% of the anomaly amplitude was added to the anomaly value. The calculated values were used to create data sets for analysis with a calculation step of 0.5 m, 1 m, 1.5 m, 2 m, 2.5 m, 3 m, 4 m, 5 m, 6 m and 7.5 m. Figures 1 and 2 show the results of the horizontal gradient calculation for the selected steps for 5% and 10% error, respectively.

Fig. 1
figure 1

Horizontal gradient theoretical (black) and calculated (grey) for horizontal cylinder anomaly with random error 5%

Full size image
Fig. 2
figure 2

Horizontal gradient theoretical (black) and calculated (grey) for horizontal cylinder anomaly with random error 10%

Full size image

The analysis of the results showed that the calculated distributions are similar for both errors. For small calculation steps, the mapping of the theoretical horizontal gradient distribution is strongly distorted, while for too large steps, the calculated distribution poorly reflects the horizontal gradient distribution. Obviously, the computational step depends on the horizontal extent of the anomaly, and therefore, the step is related to the half-width anomaly (x1/2), which is the distance between the maximum of the anomaly and the place where the anomaly reaches half of its maximum value. From the obtained distributions, it can be assumed that the measuring step should be between 0.5x1/2 and x1/2. Additional tests on other bodies (sphere, vertical cylinder) showed that a step in this interval can be considered correct. In further calculations used in the analysis of the correctness of the determination of the depth with the Werner and Extended Euler methods, a step equal to 0.75x1/2 was adopted.

The methodology of depth calculation

During the test of depth calculations using the Werner method, many solutions were obtained for theoretical data disturbed by the error (Fig. 3). This was to some extent the result of using multiple sizes of calculation windows. Therefore, in the next stage of the analysis, an attempt was made to limit the number of solutions. After many attempts, an algorithm was developed that allowed to limit the number of solutions to one.

Fig. 3
figure 3

Sample of solutions of Werner deconvolution for a horizontal cylinder with its axis at depth 4 m (random error 5%)

Full size image

In the first stage, the algorithm limited solutions due to their location on the profile. Due to the introduced error, the range of solutions is large and they cover the entire length of the profile. In actual measurements, it is usually not possible to provide a specific location of the characteristic point of the disturbing body (e.g. the centre of a sphere); therefore, it was decided to adopt a relatively wide range of possible body positions for profiles. Namely, it was assumed that the solution should be located no further than x1/2 from the anomaly centre; the remaining solutions were removed from the solution set.

In the next stage, the mean value of the depth from the remaining solutions was calculated. Then, the difference between the calculated values and the mean value was calculated. The solution error was determined as a percentage of the calculated mean depth. If there was no value greater than the assumed error in the set of differences, the mean value was considered a solution. Otherwise, a solution with the maximum deviation from the mean was sought and removed from the set of solutions. The process was then repeated until a clear solution was obtained.

This algorithm was tested on an anomaly from a horizontal cylinder with parameters identical to the previous chapter. Similarly, the data were also distorted by two random errors with a maximum value of 5% and 10% of the anomaly amplitude. In order to obtain meaningful results for each error, 100 draws were made and the results are presented in Fig. 4.

Fig. 4
figure 4

Solutions of Werner deconvolution for a horizontal cylinder with the axis at depth 4 m for data with random error (100 drawings)

Full size image

Figure 4a shows the distribution of all solutions for all 100 data sets for a 5% error. It can be easily seen that the vast majority of solutions fall within the range ( − x1/2, x1/2.), i.e. for the examined horizontal cylinder ( − 4, 4), which is better seen in Fig. 4b. There are also many solutions that differ significantly from the theoretical depth. The distributions for the 10% error are not included, but they are largely similar to those presented.

For both data sets, the above-described algorithm was used, while the calculations assumed four variants of the depth error (εz), i.e. 5%, 10%, 15% and 20%, and the results are shown in Fig. 4.

Figure 4c–j shows the depths obtained for both data sets after applying the algorithm proposed above for the above four variants. The best and most stable results were obtained for the depth error (εz) of 5%. The greater the error value, the greater the scatter of the results, which is obviously logical—incorrect solutions are less likely to be rejected. It is worth noting that the algorithm is not perfect and there are solutions that are not correct. A typical example of this is the depth value of 5 m found in Fig. 4c. The analysis of this case showed that the Werner method gave solutions in which there were two clusters of solutions with significantly different values and the algorithm chose a bad solution. Nevertheless, the goal was to obtain a single solution, and as you can see, the vast majority of solutions are correct.

Theoretical examples

The developed conditions for selecting the measurement step as well as the methodology for calculating the depth have been tested on several theoretical examples of solid figures that best described the most common anomalous bodies in the microgravity method (Table 1).

Table 1 Parameters of used solid figures
Full size table

For analysis were chosen:

  • Sphere 1 as the approximation of void in karst,

  • Sphere 2 as the approximation of an unconsolidated zone in the rock mass,

  • Vertical cylinder 1 as the approximation of a mine shaft

  • Vertical cylinder 2 as the approximation of an unconsolidated zone in the rock mass,

  • Horizontal cylinder as the approximation of a tunnel, a gallery etc.

  • Prism 1as the approximation of unknown basement;

  • Prism 2 as the approximation of an unconsolidated zone above a collapsed mine drift

  • Contact as the approximation of a fault.

All solid figure parameters were selected so that the amplitude of the anomaly was about 0.08 mGal.

For all described solid figures, distributions of microgravity anomalies were obtained. Two cases of the obtained distributions were considered for each a solid figure: the first one with a random error of 5% of an anomaly amplitude value and the second—10% of an anomaly amplitude value. For the identifiable of the results, in the case of analysis of the distributions for 5% error, the name of the solid figures is marked with the symbol “a”, for 10% error with the symbol “b”.

Each of the anomalies was used to calculate the depth using the Werner and Extended Euler method for two depth errors (εz): 5% and 10% of a depth. In the Werner method, there are two possibilities of choosing the body: a dike and a contact. The tests showed that a contact body gives acceptable results only for faults, and for this reason, a dike was used to calculate the remaining solids. The authors, of course, are aware that the resulting depths will be affected by a mismatched body shape error for some figures, but it will be interesting to see how far the results differ from theoretical. In the Extended Euler method, the structural indices corresponding to given bodies were used (0 for contact, 1 for dike and 2 for centre point). All solution results are presented in Table 2.

Table 2 Results of the depth calculations using Werner and Extended Euler methods
Full size table

For a sphere, a theoretical depth is the depth to the centre of it, for a horizontal cylinder depth to its axis, for a vertical cylinder 1 a depth to the top of it and for the rest of solid figures a depth to an average value of its top and bottom (as similar to 3D bodies). For a contact as a half-width (x1/2) of the anomaly was taken a distance between the middle of the anomaly and a place, where the value of anomaly is a quarter of its amplitude.

For better visibility of the results, depths were shown as a bar chart of percent of theoretical depth (Fig. 5).

Fig. 5
figure 5

Results of depth calculations using Werner and Extended Euler methods: a an error in relation to the theoretical depth dt; b an absolute value of an error

Full size image

The analysis of the solutions presented in Fig. 5 shows that for both errors disturbing the anomalies, i.e 5% and 10% of the anomaly amplitude, the results are similar, which confirms the solutions previously obtained for the horizontal cylinder. The values of the depth solutions depend on the estimation method. As we can see in Fig. 5a, solutions of the Werner method give generally lower values of the depth, and the Extended Euler method higher values.

As expected, the best solutions were obtained for the horizontal cylinder and contact. It should be noted that despite the significant error distorting the theoretical data (5% and 10%), the results got with both methods are very good. The calculated depth deviates from the theoretical no more than 10%. The performed calculations show that in the case of contact, increasing the fault throw reduces the correctness of the results.

The solutions for a shallow prism 1 are more divergent than for a deeper prism 2 (Fig. 5b). For the prism 1, the Euler method gives similar depths for both error values (approx. 20% of dt), while the Werner method gives more divergent results (approx. 40% of dt for the prism 1a and 20% for 1b). For a prism 2 simulating an unconsolidated zone over a collapsed mine gallery, both methods give similar values, namely 30% theoretical. This is most likely due to the fact that the shape of the solid figures can be approximated to some extent with a horizontal cylinder (or a dike).

Acceptable results were also obtained for both spheres by the Werner deconvolution method (Fig. 5b). An error of the solutions is generally in the range of 15–40% of dt. For the sphere 2a, the error is much greater, due to the imperfection of the proposed algorithm. The solution analysis showed that two solution clusters were recovered and the algorithm selected the wrong one. Choosing the second cluster would give a solution error of less than 10% of dt. The Extended Euler method gives better solutions for sphere 2, when the disturbing body is deeper, which makes the width of the anomaly larger.

The weakest solution was obtained for the vertical cylinder 1 for the Werner method, for which the calculated depths are even more than three times greater than theoretical. Also, the Extended Euler method gives weaker results, and the calculated error is overstated by 40–50%. Surprisingly, good results were obtained for the vertical cylinder 2 for the Euler method, where the depth calculation error was approx. 10%. The Werner method obtained better results for more distorted theoretical values (random error 10%). The analysis of the solutions for the 5% random error showed that they do not group into any clusters, so the final solution is close to the middle value of all solutions.

Conclusions

Finding the correct depth to the disturbing body in the gravity method is usually not easy, and in the microgravity method even more difficult. The article analyses the possibility of using the Werner and Extended Euler deconvolution methods for microgravity profile data to find bodies that may reflect disturbing bodies found in microgravity surveys. Most of these bodies are not two-dimensional, which of course introduces an additional error in determining their depth.

After testing, it was found that the calculation step should be in the range of 0.5–1 of half-width of an anomaly. Since the Werner deconvolution method gave many, sometimes very divergent, solutions, an algorithm was developed that allowed us to obtain one solution that was assumed to be correct.

The results carried out for various errors clearly indicate that the solutions obtained for the horizontal cylinder and the fault (with a small throw) are very good. Worse results were obtained for the sphere and the prism, but they also seem to be acceptable. Namely, the absolute error in calculating the depth is relatively large, 20–40% percent, but due to the shallow depth of the quantitatively disturbing bodies, these values are not significant. The worst results were obtained for the vertical cylinder 1 reflecting the mine shaft. The results of Werner's deconvolution have a large error, which means that this method should not be used for this type of body. The causes should be looked for in the small horizontal range of the anomaly, and thus large changes in the gradients resulting from the imposed error.

The obtained results make it possible to use Werner and Extended Euler deconvolution in microgravity bodies close to a horizontal cylinder and a fault, using the proposed algorithm. The results for three-dimensional bodies are burdened with a much larger relative error, but quantitatively they are also acceptable.