Estimating the Extreme Values in Gridding Data using Gauss Elimination Method: A New Approach in Potential Field Processing

The position of a maximum point of a function depends on its coefficients and order. The maximum horizontal gradient method is a popular method that greatly contributes to the detection of maximum points and approximation of geological structures edges. By adopting a mathematical logic, Blakely and Simpson established a quadratic function based on the characteristic of three points of a straight line in the fundamental directions. However, for potential field data like gravity and magnetic data, the coefficients of a quadratic function in each direction are not only dependent on the values of three points on a straight line, but also, they depend on the values of the surrounding points. This article proposes an algorithm which can detect maximum points more effectively in order to delineate geological structures boundaries from potential field data. The proposed algorithm uses a 3×3 neighborhood data grid to establish a two-variables function and to determine its coefficients by applying the Gaussian elimination method. After the two-variables function has been established, the algorithm estimates any extreme points and their positions from a set of four single-variable functions which correspond to the horizontal, vertical and the two diagonal directions by the cases x = 0, y = 0, y = -x and y = x of the main function. Finally, the conditions to detect the maximum point from the extreme points are defined. The validity of the algorithm was demonstrated on synthetic datasets generated by two different model structures. A real data application of the method has also been realized by estimating the geological boundaries by gravity data in the Vietnam’s continental shelf. The results obtained from the synthetic applications of the algorithm proved that it can determine more maximum points as compared to Blakely and Simpson method, and as a result, in all the test cases, it has drawn the real boundaries of the model structures more accurately. The application results of the method on real data revealed new boundary delineations in the research area, interpreted to be faults or fractures which lies between deep trench in the East Vietnam Sea.


Introduction
Information of the potential field source boundaries plays a significant role in mineral resource exploration.Thus, determining the boundaries of geological structures by evaluating potential field data is a frequently requested goal in shallow and deep sub-surface investigations (Ahmed et al., 2017).At present, there are many different methods which applied to estimate the geological boundaries from analyzing potential field data.These methods have provided detailed studies ranging from the theoretical basis to the numerical model, as well as its application to real data.It is worthy to mention first the method of the normalized full gradient of gravity anomaly (NFG) of Berezkin (1967).The original method was later applied, improved and developed by different authors (Aghajani et al., 2009;Oruc and Keskinsezer, 2008;Oruc, 2012;Karsli and Bayrak, 2010;Ekinci et al., 2013;Ekinci andYiğitbaş, 2012, 2015;Sheng and Xiaohong, 2015).
A second type of methods is the horizontal gradient method (Cordell, 1979), the total horizontal gradient (Cordell and Grauch, 1985) and the analytical signal amplitude (Nabighian, 1972;Roest et al., 1992).However, the performance of these filters to clearly detect the edges of deep sources is known to be limited (Eldosouky et al., 2020;Pham et al., 2020aPham et al., , 2020b)).In order to enhance the detection of the boundaries of deeper structures, many other methods have been proposed, based on the horizontal and vertical derivatives of the potential field.For instance, the enhanced analytical signal method of Hsu et al. (1996), the magnitude of the analytic signal of horizontal gradients of Bournas and Bake (2001), the enhanced horizontal gradient method introduced by Fedi and Florio (2001), the multiscale edge detection method presented by Cella et al. (2009) and the method using the derivatives of the directional analytic signals of Beiki (2010).Besides, many other studies also mentioned that the amplitude-based filters in general are limited in simultaneously determining source edges of structures being at different depths due to the effect of varying amplitudes of the potential field (Eldosouky et al., 2020;Pham et al., 2020aPham et al., , 2020b)).To overcome this problem, a number of phase-based filters have been developed over the years that their basis in determining the source edges can provide the advantage of simultaneously balance of structures amplitude at different depths (Miller and Singh, 1994;Wijns et al., 2005;Cooper and Cowan, 2008;Verduzco et al., 2004;Ferreira et al., 2013;Chen et al., 2017;Nasuti and Nasuti, 2018;Pham et al., 2019).However, although these methods are very effective in this aspect, they may bring misleading edge information which may go considerable beyond the true edges (Zhou et al., 2013;Ma et al., 2014;Eldosouky et al., 2020).On the other hand, using the derivatives of the field of low amplitude anomalies which produced by deep structures can also be applied to delineate structure edges, especially when the maxima of the transformed field are quantitatively processed to locate the peak values of the algorithms.The maximum horizontal gradient method of Blakely and Simpson (1986) based on interpolating a parabola through three successive data points and the curvature method of Phillips et al. (2007) based on fitting a restricted quadratic function to the nine points of a moving window are the most common in this regard.
In this work, only the method of Blakely and Simpson (1986) and the one of Phillips et al. (2007) are discussed in detail.Regarding the method of Blakely and Simpson, the main limitation is that it can only detect the maxima points that lies on four directions along the horizontal, vertical, and two diagonals and is the maximum point of a quadratic function that is established by three points on a straight line.Therefore, the accuracy of the approximated buried geological structure boundaries may not be efficient enough.In the method of Phillips et al. (2007), two-variables function was established to detect the maximum point.Therefore, as an advantage, the maximum point does not necessarily lie along the four principal directions, but also it can lye anywhere within the frame of the nine data points.Nevertheless, the determination of the coefficients of the quadratic surface of this method has unstable accuracy because the number of equations is nine where as the number of variables (the number of coefficients) are six, thus, it is not a unique root.
In the current research, we proposed a new algorithm which originated from a two-variables function, then employing the Gauss elimination method in order to improve the detection of the maxima points.The efficiency of the proposed algorithm was demonstrated on synthetic data samples where the results have been compared with results which obtained by other common methods.The introduced algorithm has also been applied on a real Bouguer anomaly gravity data of Vietnam continental shelf as a practical application.

Materials and Methods
The first step in locating the maximum is to solve within each 3×3window of grid data for the coefficients of a two-variable function of a fitting surface passing through the nine data points of the window.The two-variable function that we focus on is expressed as: Here we propose using the Gauss elimination method to determine the coefficients (ai) of the function (1).For each 3×3 window of grid data, we get 9 equations: { a 1 x 1 2 + a 2 y 1 2 + a 3 x 1 2 y 1 2 + a 4 x 1 2 y 1 + a 5 x 1 y 1 2 + a 6 x 1 y 1 + a 7 x 1 + a 8 y 1 + a 9 = g 1 a 1 x 2 2 + a 2 y 2 2 + a 3 x 2 2 y 2 2 + a 4 x 2 2 y 2 + a 5 x 2 y 2 2 + a 6 x 2 y 2 + a 7 x 2 + a 8 y 2 + a 9 = g 2 ⋮ a 1 x 9 2 + a 2 y 9 2 + a 3 x 9 2 y 9 2 + a 4 x 9 2 y 9 + a 5 x 9 y 9 2 + a 6 x 9 y 9 + a 7 x 9 + a 8 y 9 + a 9 = g 9 (2) Where 1 to  9 and  1 to  9 are the coordinates of the nine data.
where gmax is the value of the function g(x,y) at point M. A counter N is increased by one for each satisfied case.The parameter N as the significance level of the maximum.We found that N ≥ 2 produces the most useful maps.Table 1 summarizes a comparison of the differences in the conditions used in this paper with those of Blakely and Simpson (1986).Base on the theoretical basis presented above, a MATLAB code has been built to detect the extreme points and the maximum points of a gridded data.Application samples of the algorithm linked to the code has been given in the following sections.

Parameter Blakely and Simpson (1986) Current study
The number of equations 4 4 The order of polynomial 2 2 (for case: x=0, y=0); and 4 (for case: y=-x, y=x) Determining the coefficients of function Based on 3 points Based on 9 points Detecting the maximum point The value at a central grid location exceeds that of the outer two points (straight line) The value at the extreme point exceeds that of the outer two points-the extreme point (straight line)

Synthetic Data Applications
In this section, we analyze the performance of the proposed algorithm on two synthetic gravity data produced with different model configurations.

Model 1
The first model consists of two prism models in same depths with different density contrasts, and one side of their edges are close to each other.Fig. 2a shows the perspective view of the model structures (labeled M1-A and M1-B), and the numerical parameters of each are listed in Table 2.The synthetic gravity anomalies due to the prisms were calculated using the formula given by Talwani and Ewing (1960) on a regular grid of 101×101 observation points with a spacing of 0.1 km.Fig. 2b and Fig. 3a show the noise-free gravity anomaly due to the models and its total horizontal gradients used for the maxima detections, respectively.The maximum detection results obtained by the method of Blakely and Simpson (1986) and the results obtained with the present algorithm are shown in Fig. 3b comparatively.From the results obtained, it can be seen that the number of maximums obtained by the proposed algorithm (red markers in Fig. 3b) for both the structures models is higher than the number obtained according to the Blakely and Simpson's method (blue markers in Fig. 3b) and gives more detail the actual edges of the models M1-A and M1-B.Besides, the locations of maxima estimated by using the Blakely and Simpson's algorithm for structure M1-A (lower density contrast source) seems limited in delineating the boundaries which are close to the edges of M1-B.

Model 2
The second model is more complex, and composed of five prisms (M2-A, M2-B, M2-C, M2-D and M2-E) with different density contrasts.The three-dimensional view of the model is illustrated in Fig. 4a.The parameters of the sources are listed in Table 3 and Fig. 4b, the gravity anomaly of the model calculated on a 101×101 grid with an interval of 0.1 km in x and y directions.The total horizontal gradient map of the data in Fig. 4b is shown in Fig. 5a.Fig. 5b shows the maximum locations obtained from the proposed algorithm and the Blakely and Simpson method.As can be seen from this figure, both methods are effective in detecting the peaks of the total horizontal gradient, however, the proposed method can detect the peaks over the west edge of the M2-C in more detail.In addition, the proposed method also detects more maximum locations over the edges of other sources than the Blakely and Simpson method.Clearly, the use of the proposed method proved useful to bring the maximum locations dealing with source edges in more detail.

Real Data Application
Once demonstrated the capabilities of the proposed algorithm to delineate the boundaries of different model bodies from the analysis of their synthetic gravity anomalies, an application example using a real gravity data is presented here.The study area comprises the continental shelf of Vietnam and adjacent areas, extending between latitudes 8°N and 18°N, and longitudes 108°E and 116°E (Fig. 6).The aim of the study here is to detect the maxima points of the total horizontal gradient amplitude function of the second order vertical derivatives of the gravity field, and approximate the geological boundaries in the study area.The Bouguer anomaly data of the region have been obtained from the available free-air gravity anomaly (Sandwell and Smith, 2009;Sandwell et al., 2013Sandwell et al., , 2014) ) and the seabed topography (Smith and Sandwell, 1997) by using the algorithm of Parker (1972)   It is a well-known fact that the gravity field includes the effect of both deep and local shallow sources.Thus, the data has been processed with an upward continuation operator for different altitudes (10, 15, 20, 25, 30 km) ( Reynolds, 1997) in order to remove the residual anomalies corresponding to both the seabed topography variations and the local geological structures.Next, the first vertical derivatives of each the upward continued gravity field have been calculated.The resulting grids have been used to calculate both the horizontal gradient and the horizontal gradient amplitude.Finally, from these horizontal gradient amplitude grids, the maxima points have been detected by the methods (using n>=2 for both cases).The maxima points detected in the study area at different upward continuation altitudes based on the proposed algorithm are shown in Fig. 8.In order to compare in detail, the results of the geological boundaries obtained from the two different approaches, a small zone of 2º×2º (112º-114º E, 11º-13º N) has been selected (Fig. 9).Yellow points represent the location of the geological boundaries obtained by the present approach whereas red points are the result of the Blakely and Simpson's approach.Cyan arrows represent the horizontal gradient vectors, the length of vectoris multiplied with 2 in Fig. 8 and with 10 in Fig. 9.A comparison of the results at different upward continuation altitudes show that the Blakely and Simpson's approach provides a lower amount of data points defining the geological boundaries.Moreover, their horizontal gradient vectors have small amplitude and exhibit a more confusing direction.The geological boundaries can be approximated by linking the maxima points together into lines.This linking needs to be combined with the direction of their horizontal gradient vectors.If these directions are towards a center, it may be interpreted as an uplift; if it shows directions from a center towards outside, it may be interpreted as a depression.All the results achieved are an interpretation, and must be confirmed with joint analysis of other geological information and geophysical data.
Fig. 9f displays an overlapped map of all the boundaries detected by both the methods for the whole range of different upward continuation altitudes.As it can be seen, many boundaries have been detected by the two different approaches.However, one boundary (represented by the green line) has only been detected by this new proposed approach.This boundary is introduced as a spreading axis or spreading ridges of the South China Sea (or the East Vietnam Sea), but isn't described detail (Zhao et al., 2020;Wang et al., 2020;Zhang et al., 2018;Ding et al., 2018;Pubellier et al., 2017;Claude et al., 2016).In some paper, this boundary isn't described (Dung et al., 2013;Hiep N., 2005) and, consequently, it reveals a new result for the study area.In order to verify the new boundary and obtain information about its nature, more geophysical and geological works should be carried out.

Conclusions
An algorithm that can detect the maxima points and approximate the geological structure boundaries by analyzing potential field data has been proposed.The algorithm uses a 3×3 neighborhood data grid to establish a two-variables function, whose coefficients are established by the Gauss elimination method.The function of two variables is established for nine data points (a 3×3 data grid), that are examined to obtain the extreme points from 4 functions of one variable that corresponds with 4 specific cases of the function (cases: x=0, y=0, y=-x and y=x).Based on the theoretical basis described in this work, a MATLAB code was built to perform the algorithm proposed.By this, different synthetic models have been tested to demonstrate the validity and robustness of the algorithm as well the linked code.From the results obtained by applying both synthetic models and real data; it has been shown that the maxima points defining the location of the geological boundaries detected by the new approach are more abundant than Blakely's approach.Moreover, using the maxima points detected by the proposed approach have delineated more accurately the edges of the subsurface structures.By using this algorithm with real data, a new geological boundary has been defined on the research area.This feature has been interpreted as a large fracture or fault lines which lies between a deep trench of the East Vietnam Sea.

Fig. 1 .
Fig. 1.Location of grid intersections to detect a maximum point near gi,j

Fig. 2 .
Fig. 2. Model-1 (a) the three-dimensional view of the prisms construct; (b) the gravity anomaly due to the models (the white lines are the borders of the actual prism models in plain view)

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3. Model-1 results; (a) the total horizontal gradient of data in Fig. 2b; (b) maxima locations obtained by using algorithm of Blakely and Simpson (1986) (blue marker) and by the present algorithm (red marker) as described in Pedraza et al. (2014).Both datasets have a resolution of 1 minute in east and north directions.The resulting Bouguer anomaly map is shown in Fig. 7 comprising gravity values ranging from -35 m Gal to 325 m Gal.

Fig. 6 .
Fig. 6.Location map of the study area for the real case application of the proposed method

Fig. 7 .
Fig. 7. Bouguer anomaly map of the study area

Fig. 9 .
Fig. 9. Edge detection results in the sub-region shown in Fig.8 obtained from varying upward continuations of the field data; (a) 10 km; (b) 15 km; (c) 20 km; (d) 25 km and (e) 30 km.Here the results obtained from Blakely and Simpson method and the proposed algorithm are indicated with red and yellow markers, respectively; (f) shows the combined results from varying continuation fields by the present algorithm Let point M is the maxima point then M need satisfied condition:

Table 1 .
Blakely and Simpson (1986)ns used in this study with those byBlakely and Simpson (1986)

Table 2 .
Geometric and density contrast parameters used in the first model case shown in Fig.2