An edge-assisted smooth method for potential field data

Edge detection is one of the most commonly used methods for the interpretation of potential field data, because it can highlight the horizontal inhomogeneous of underground geological bodies (faults, tectonic boundaries, etc.). A variety of edge detection methods have been reported in the literature, most of which are based on the combined transformation results of horizontal and vertical derivatives of the observations. Consequently, these edge detection methods are sensitive to noise. Therefore, noise reduction is desirable ahead of applying edge detection methods. However, the application of conventional filters smears discontinuities in the data to a certain extent, which would inevitably induce unfavourable influence on subsequent edge detection. To solve this problem, a novel edge-preserving smooth method for potential field data is proposed, which is based on the concept of guided filter developed for image processing. The new method substitutes each data point by a combination of a series of coefficients of linear functions. It was tested on synthetic model and real data, and the results showed that it can effectively smooth potential field data while preserving major structural and stratigraphic discontinuities. The obtained data from the new filter contain more obvious features of existing faults, which brings advantageous to further geological interpretations.


Introduction
Gravity and magnetic methods usually play an important role in geological study and mineral resources exploration (Lee et al. 1990;Zhang et al. 2007;Wang et al. 2018). To study underground geological bodies, it is useful to analyse horizontal inhomogeneous of them. In general, the inhomogeneity often occurs along the boundaries of geological bodies, and usually agree well with the gradient belts in gravity and magnetic data. Among various processing methods, the edge detection method can be used to extract these important characteristics. Therefore, edge detection is usually considered as one of the most important processes for the interpretation of potential field data (Fedi & Florio 2001;Li et al. 2013;Guo et al. 2015;Wang et al. 2015;Yan et al. 2019).
Many mathematical algorithms for edge detection have been proposed in the literature, such as the total horizontal derivative (THD) (Cordell & Grauch 1983), analytic signal (AS) (Nabighian 1984), tilt angle (Miller & Singh 1994) and theta map (Wijns et al. 2005). These four methods were the most frequently used edge detection methods, which were based on a combination of different vertical and horizontal derivatives of potential field data. In recent years, various improved methods have also been proposed, for example Wang et al. (2009) proposed a new edge detection method based on the normalised vertical derivative of the THD. Ma & Li (2012) introduced the ratio of the horizontal derivative to the maxima of nearby values to recognise edges. Cooper (2014) suggested a modified AS amplitude method. Wang et al. (2015) proposed the principal component analysis to the curvature gravity gradient tensor to outline edges. Wang et al. (2017) improved Euler inversion and phase congruency to determine lineaments. Yuan et al. (2016) first used the horizontal directional theta to define a new edge detector. Wu et al. (2017) utilised the THD and the modulus of full tensor gravity gradient to recognise edges. Yan et al. (2018) applied the multi-scale edge detection method to marine magnetic survey data to identify gold ore-controlling faults. However, these methods are often very sensitive to noise. If the original potential field data are noise contaminated, the obtained edge detection results will be distorted. In practical applications, noise reduction is usually executed before applying other gradient-based methods, such as Wiener filtering (Pawlowski & Hansen 1990), median filtering (Tukey 1974), statistical filtering (Schau 1980), Laplacian filtering (Kim et al. 2011) and preferential filtering (Guo et al. 2013). These methods work well for noise attenuation, but also remove some useful information, especially abrupt boundary signatures.
This issue also arises in the interpretation of other geophysical data, such as seismic data. Different scholars have proposed a number of algorithms to solve this problem. For example, Luo et al. (2002) applied the Kuwahara multi-window filter as a pre-processing step of seismic data (Kuwahara et al. 1976) to form an edge-preserving smoothing algorithm. Fehmers & Hocker (2003) proposed structure-orientated filtering to remove high frequency noise. The core of this edge-preserving method is a smoothing operation parallel to seismic reflections. Nasher (2006) introduced an edge-preserving noise reduction technique based on which a 3D subsurface is divided into a number of blocks, and values of blocks neighbouring the given one are used to reduce noise. Lu & Lu (2009) proposed an edge-preserving polynomial-fitting method for cases where the traces of seismic events are nonlinear.
However, a few studies in the literature have examined the application of edge-preserving smooth operators to potential field data. In this paper, we are the first to apply the guided filter, which was proposed by He et al. (2013) to process digital images, for edge-preserving smoothing of potential field data. The detailed steps and applications of this method are introduced and tested on designed models. The results indicate that the guided filter could effectively smooth gravity and magnetic data while preserving major structural and stratigraphic discontinuities. Therefore, the subsequent performance based on gradient transformation is significantly improved. Finally, we apply the proposed method to real magnetic data from southwestern Fujian in China.

Methodology
According to the mathematical theory of local linear models, a specific point of a function is linearly related to the points around it (He et al. 2013). Therefore, a complex function can usually be represented by a linear combination of numerous local functions. When the value of a specific point is required, we can calculate the values of all linear functions containing the point, and their average value is the final result.
In general, the expression for removing noise is as follows: where Res(x, y) and Ori(x, y) are the expected result and the original data, respectively, N(x, y) is noise, and x and y are the coordinates of the data points. The guided filter method regards the result as a linear model of the guided dataset I(x, y). The expression for it is as follows: where w x 0 ,y 0 is the calculated window, the coordinates of the central are (x 0 , y 0 ). a x 0 ,y 0 and b x 0 ,y 0 are the coefficients of each window. If we take the gradient of both sides of equation (2), the following can be obtained: where a x 0 ,y 0 represents the degree of edge preservation in the result. It is well known that the gradient of the anomaly plays an important role in edge detection. Equation (3) means that the gradient of the result is similar to that of the guided data. That is to say, the guided filter has an edge-preserving function. The guided data of the initial iteration can be selected as a constant dataset, such as a constant 1, and each subsequent iteration can use the result of the previous iteration as the guided filter. To minimise the difference between the result and the original anomaly while filtering out noise, the least squares algorithm is used to fit equation (1): where is the regularisation parameter that prevents the obtained a x 0 ,y 0 from being too large. According to the principle of the least squares algorithm, equation (4) has an optimal solution as follows: Here, all values are related to window w x 0 ,y 0 , |w| is the number of points in each window, x 0 ,y 0 and x 0 ,y 0 2 are, respectively, the mean value and the variance of I(x, y), ∀x, y ∈ w x 0 ,y 0 and Ori x 0 ,y 0 is the mean value of Ori(x, y), ∀x, y ∈ w x 0 ,y 0 . When calculating the coefficients of every window, each point is used simultaneously by multiple windows. That is to say, each point is described by multiple linear model functions. When we want to obtain the result for each point, we need to only calculate the average coefficients of all linear functions containing this point. The expression is as follows: We summarise the flow chart of the guided filter method as Figure 1. We now directly analyse the edge-preserving function of the guided filter through a simple case in which the guided anomaly was equal to the original anomaly; that is to say, I ≡ Ori(x, y). Equations (5) can then be rewritten as follows: The regularisation parameter is generally a positive number. According to the change in the gradient of the guided anomaly, the two following cases should be considered: (i) the guided anomaly changes significantly in window w x 0 ,y 0 . In this case, x 0 ,y 0 2 >> and so a x 0 ,y 0 ≈ 1, b x 0 ,y 0 ≈ 0. (ii) The guided image changes slightly in window w x 0 ,y 0 . In this case, x 0 ,y 0 2 << and so a x 0 ,y 0 ≈ 0, b x 0 ,y 0 ≈ x 0 ,y 0 . Once the parameters of each window have been obtained, equation (6) is used to obtain the final result. When the value of the central point in the window changes significantly compared to the value around it, the result remains unchanged, which means that a(x, y) ≈ 1, b(x, y) ≈ 0, and Res(x, y) ≈ Ori(x, y). When the value of the central of the window changes slightly compared with the values of points around it, the result is the average of the values of these points, which means that a(x, y) ≈ 0, b(x, y) ≈ x,y , Res(x, y) ≈ x,y . In other words, the result can retain the corresponding gradient information when it changes significantly, and can also smooth the anomaly when the gradient changes slightly.

Model one
We established a theoretical magnetic susceptibility model that contains two vertical prisms to verify the performance of the proposed method. The inducing geomagnetic field had a strength of 5 × 10 4 nT, with an inclination of 90°and a declination of 0°. Figure 2a shows the 3D view of the model. The range of the observation surface was 10 × 10 km, with a spacing of 0.1 km. The top and bottom of the prism on the left were 0.3 and 0.8 km from the centre at (2.5, 5.25), and the intensity of magnetisation was 1.5 A m −1 . The top and bottom of the prism on the right were 0.6 and 1.5 km from the centre at (6, 5), and the intensity of magnetisation was 2 A m −1 . Figure 2b shows the synthetic magnetic anomalies produced by the model, in which the dotted lines in all panels represent the real locations of the two prisms. Random noise with a zero mean and a standard deviation of 50 nT was added to the synthetic data to form the noisy data in figure 2c, and this was used for the smoothing tests. As there are many edge detection methods, the most commonly used THD method was randomly selected to demonstrate the necessity of smoothing. We first calculated the THD of the noisy data, and the result is displayed in figure 2d. It can be seen that little useful information on the real boundaries could be identified in this case. This test indicated that the conventional derivative-based methods amplified the components of noise. Thus, smoothing usually should be applied before the computation of the derivative-based transforms.
To test the effect of the radius of the window on the guided filter, windows with radii of 1, 2 and 3 were used to filter the noise-contaminated magnetic anomaly with model one, and the results are displayed in figure 3. It can be noted that with the increase in the radius of the window, most effective anomalies were filtered out, but windows with smaller radii achieved a better effect. However, how to choose the suitable appropriate radius is still a problem that deserves further study.
We then applied the proposed smoothing method to the noisy data obtained. To compare the influence of the  filtering results on the subsequent boundary enhancement processing based on derivative transformation, the most common boundary recognition method, THD processing, was used on the result. Figure 4 shows the results of filtering and the corresponding THD results by using the guided filter, Gaussian filter and Wiener filter with different numbers of iterations. It is clear from all the results in figure 4 noise was filtered out to some degree compared with the initial noisy data. After five iterations, the result of the guided filter was still contaminated by unfiltered noise as evidenced by the wavy disturbance in figure 4a1, and the results of the Gaussian filter and the Wiener filter were much smoother. This indicated that the guided filter required more iterations. After 10 iterations, all results had improved significantly compared with the initial noisy data, and the guided filter yielded better results in terms of preserving the amplitudes of smaller anomalies than the Gaussian and the Wiener filters. The THD results in figure 4a2 and d2 show that the maximum value of the THD method could accurately characterise the edges of the prisms even with a large number of iterations. In contrast, the Gaussian filter did well to filter out noise after fewer iterations, but also filtered out information on the horizontal boundary of the prism on a small scale, as shown in figure 4b2 and e2. This happened because conventional filters usually smoothed the sharp gradient in the original data. Similarly, the Wiener filter also filtered noise to a certain extent as shown in figure 4c1 and f1, but the anomaly still contained noise. It is evident that the enhanced linear information obtained using THD was not as smooth as those of the other two methods in figure 4c2 and f2. Overall, the guided filter preserved gradient information much better than the other two filters, and the effect of the reduction was small, especially for cases of a large number of iterations. That is to say, the proposed guided filter had a superior edge-preserving impact, and multiple guided filters could be used to remove noise while preserving information on small structures. However, the amplitudes in the case of all results were smaller than the real case. In other words, amplitude reduction was unavoidable for all filters. This feature was especially serious near areas with large gradients, namely along the real boundaries. The removal of gradient information posed difficulties for further gradient-based transforms.
116 Figure 4. Filtering results of (a1) guided filter, (b1) Gaussian filter, (c1) Wiener filter with five iterations for model one. THD results of (a2) guided filter, (b2) Gaussian filter, (c2) Wiener filter with five iterations for model one. Filtering results of (d1) guided filter, (e1) Gaussian filter, (f1) Wiener filter with 10 iterations for model one. THD results of (d2) guided filter, (e2) Gaussian filter and (f2) Wiener filter with 10 iterations for model one.  Therefore, developing a filter with an edge-preserving effect is important. This is a significant feature of this study.
We extracted the profiles at northing = 5 km from the results shown in figure 4 to compare the results of the three filters at the same level. Figure 5 shows values along this profile in different cases. In all curves, positions with larger slopes represented the boundaries of the geological bodies. It is significant that the slopes along the real boundaries in the results of the guided filter were much larger than those of the Gaussian filter. Although the Wiener filter yielded slopes similar to those of the guided filter, it was obvious that its results were not smooth. This indicated that the guided filter had better edge-preserving effects than the Gaussian and Wiener filters. Moreover, the amplitudes of the curves of the guided filter were much closer to the theoretical curves, indicating that amplitude reduction due to the guided filter was also much less pronounced than in the Gaussian and Wiener filters. This effect was far more evident for a shallow, small prism. Figure 6 shows the root mean square (RMS) errors between the results using different numbers of iterations on the three filters and the theoretical data. It is explicit that the RMS error of the filters increased with the number of iterations, but that of the guided filter was significantly smaller than those of the other two filters. Therefore, a greater number of iterations could be used to filter out noise  Figure 8. Filtering results of (a1) guided filter, (b1) Gaussian filter, (c1) Wiener filter with five iterations for model two. THD results of (a2) guided filter, (b2) Gaussian filter, (c2) Wiener filter with five iterations for model two. Filtering results of (d1) guided filter, (e1) Gaussian filter, (f1) Wiener filter with 10 iterations for model two. THD results of (d2) guided filter, (e2) Gaussian filter and (f2) Wiener filter with 10 iterations for model two. while retaining edge information. That is to say, the guided filter was more efficient at suppressing noise, smoothing anomalies and preserving edges at the same time. According to figure 5, the amplitudes were smaller than the real case and with increasing numbers of iterations the reduction in amplitude became more serious. Therefore, the RMS curves presented a trend of gradual increase. Consequently, how to reasonably select the number of iterations is also a problem worthy of further study.

Model two
Further, we established a more complex theoretical magnetic susceptibility model containing four vertical prisms to validate the proposed method. The inducing geomagnetic field had a strength of 5 × 10 4 nT, with an inclination of 90°and a declination of 0°. Figure 7a shows the 3D view of model two. The range of the observation surface and subdivisions of the grid were identical to those in model one. The top and bottom of the model were 0.5 and 1.5 km, and the intensity of magnetisation was 2 A m −1 . Figure 7b also shows the theoretical anomalies generated by model two, in which the dotted lines also represented the real boundaries of the model. Random noise with a zero mean and a standard deviation of 50 nT was added to form the noisy data in figure 7c. The THD result is displayed in figure 7d, in which it was impossible to detect information on the boundary entirely.
As in case of model one, we applied the three filtering methods to filter noise with different numbers of iterations. Figure 8 shows the results of filtering and the corresponding THD results. As can be seen from all the results in figure 8, the noise was filtered out to a certain degree in comparison with the initial noise-contaminated anomalies. After five iterations, the results of the guided filter and the Gaussian filter still contained some noise, while those of the Wiener filter contained more noise. This was clear from the wavy disturbances in figure 8a1, b1 and c1. Meanwhile, the extreme value in the result shown in figure 8a2 clearly depicts the boundary of the model. As the amplitude of the result in figure 8b2 was small, it was difficult to accurately identify the boundary. Furthermore, the result in figure 8c2 clearly shows noise in the THD result, and the maximum value could not be aligned accurately to depict the boundary of the model. After 10 iterations, although the wave disturbance of the guided filter in figure 8d1 showed that some noise persisted, the amplitude of the results of the Gaussian filter in figure 8e1 was significantly reduced, whereas the contour line of the results of the Wiener filter in figure 8f1 was still disturbed. Similarly, compared with figure 8e2, the maximum value in figure 8d2 still characterised the edges of the model. Moreover, compared with figure 8f2, the enhanced linear information in figure 8d2 was smoother and more accurate, and hence was easier to recognise. Above all, the guided filter preserved gradient information much better, which thus provided a better foundation for the data used in subsequent edge detection based on derivatives. In particular, after more iterations, the guided filter filtered out noise and maintained higher amplitudes of gradient anomalies much better than the other two filters. It is worth mentioning that the guided filter could filter out noise effectively in the case that the actual anomaly was smoothed slightly, but various subsequent edge detection methods were needed to identify the edge information. Therefore, we reached the conclusion that the guided filter has an edge-preserving function.

Application to measured magnetic data
To compare the difference in process between the filters on actual data, we used aeromagnetic data from the Dapai polymetallic deposit in Fujian Province of Southwest China 120 Figure 10. Filtering results of (c1) guided filter, (d1) Gaussian filter, (e1) Wiener filter and the corresponding THD results of (c2) guided filter, (d2) Gaussian filter, (e2) Wiener filter. Different components between the filtering result and the corresponding THD result of (c3) guided filter, (d3) Gaussian filter and (e3) Wiener filter.
as an example to verify the performance of the proposed method. Figure 9 shows the simplified regional geological structure map of the deposit. The study area was located in the Tethys tectonic belt of the Continental margin system around the Pacific Ocean. According to the literature, the mining area mainly contained monzonitic granite porphyry. It provides mineralised hydrothermal fluids for such polymetals as iron, lead and zinc in the deposition medium (Lin 2011;Zhang et al. 2018). The aeromagnetic data were the results of 1:10 000 areal magnetic data measured by using an AS350B3 helicopter equipped with an airborne magnetic measurement system. It was also one of the important scientific achievements of our group. The direction of the survey line was NW and the distance is 100 m between two lines. In this paper, the magnetic data used consisted of 301 rows and 251 columns, with a spacing of 0.02 km. Figure 10 parts a and b, respectively, show the results of the reduction-to-pole magnetic data and the result of its THD processing. The 121 white line represents the geological inference fault. It was clear that owing to the large amount of noise in the data, the result of THD processing did not enhance the known faults. Therefore, it was necessary to filter the anomalies before taking the derivative transformation.
Using the same method of processing as for the previously considered models, we applied the guided filter, Gaussian filter and Wiener filter using five iterations to the anomalies, and the results are shown in figure 10c1, 10d1 and 10e1. We also calculated differences between the original anomaly and the results of the three filters, and these are shown in figure 10c3, 10d3 and 10e3. Different results used the same colour scheme, making it easier to distinguish the differences. It is obvious that these results removed some linear structural information. However, more information has been filtered out in figure 10d3, as indicated from the range of the black ellipses. In the meantime, little information shows in the range of the black ellipses of figure 10e3, but we can guess there was still some noise in the result by analysing the models. This linear structural information plays an important role in fault recognition by using derivative-based transforms.
Following this, THD was applied to the results of figure 10c2, 10d2 and 10e2, in which the black curves show the positions of the faults drawn according to the result of the guided filter and the blue curves show those indicated on the geological structure map. It can be seen that most faults identified by the guided filter were consistent with the geological map. This can provide a basis for further geological interpretation in figure 10c2. However, part of the gradient information was blurred in the result from the Gaussian filter, which reduced the ability to recognise faults from the THD result in figure 10d2. Although the gradient information was retained by the Wiener filter, it did not clearly filter out noise and a serious boundary effect also appeared when the same number of iterations was used for the other filters. In general, all three filtering methods reliably filtered noise out of the anomalies. The difference was that the guided filter retained more gradient information while filtering out noise, thus providing important data for the subsequent processing based on the derivative transformation. This improved the results of the latter, which was beneficial for subsequent geological interpretation.

Conclusion
In this paper, a new method of removing noise while preserving edges in potential field data is proposed. The availability of the proposed filter is verified by synthetic model and real data. All of the tests verify the effectiveness of the guided filter in terms of edge-preserving. Of course, as the guided filter cannot be used to detect edges directly, we propose that it be used to smooth anomalies and simultaneously cooperate with other boundary detection enhancement methods.
It also should be noted that the number of iterations should be selected carefully and this is worth studying further.