NHF as an Edge Detector of Potential Field Data and Its Application in the Yili Basin

: Edge enhancement is a frequently used transformation of potential ﬁeld data. Its goal is to sharpen the position of the subsurface structures. Here we propose a new method to enhance the edges of the sources causing the potential anomalies called normalized Harris ﬁlter (NHF), which is based on the Harris ﬁlter and amplitude balance. Three synthetic data sets are used to evaluate the performance of the proposed approach. The presented approach provides a better estimation of the sources’ edges when compared to the other methods. The proposed method is robust to noisy data and can avoid the generation of artiﬁcial edges, thereby reducing the ambiguity of interpretation. The testing on real data set from the Yili basin in Northwestern China demonstrates that the new approach highlights several anomalies not shown in the geological map or other methods. The proposed approach also shows the advantages of gradually enhancing the edges of the deep-seated structure. The results demonstrate that the proposed approach may be a better detector in qualitatively determining the edges of sources causing potential ﬁeld data.


Introduction
Potential field data contains all the signals from geological bodies with different locations and geometries and has a volume effect. Many data processing techniques have been developed to facilitate the interpretation. In recent years, edge enhancement has received increasing attention, benefiting from its simplicity and effectiveness. The main purpose of edge enhancement is to delineate the edges and positions of the subsurface structures. The most popular edge detectors are total horizontal derivative [1] (THD, Equation (1)), tilt angle [2] (TA, Equation (2)), THD of tilt angle [3] (THDT, Equation (3)), theta map [4] (theta, Equation (4)), the tilt angle of the THD [5] (TAHG, Equation (5)), the tilt angle of the analytic signal amplitude [6] (TASA, Equation (6)). The approaches mentioned above have been evaluated in many works [7][8][9][10][11]. In addition to THD, other approaches involve calculating the vertical derivative or second-order derivative of the potential field data, which causes the amplified noise in the data. Therefore, it is necessary to perform denoising before or after the derivative process. with: where f denotes the potential field data. Image processing techniques have been widely used in geophysics. For example, Ref. [12] employed sun shading to extract the features from potential field data. Ref. [13] used the Sobel filter to detect the fault and salt boundary. Ref. [14] used the 2D Gabor filter on the magnetic data to interpret the subsurface lineaments. Ref. [15] used the Canny filter [16] to perform the edge-guided full-waveform inversion. Ref. [17] employed phase congruency to detect lineaments embedded in magnetic data. The Harris filter [18,19] provides strong invariance to rotation, scaling, and noise. It has been widely used in many scientific areas based on image processing techniques, such as camera calibration [20], image stitching [21], and unstructured finite-element meshing [22].
In this work, we employ the Harris filter to detect the edges of the subsurface structures from potential field data. Amplitude balance is applied to equalize the response from the sources at different depths. Depth scales can range from kilometers to meters, depending on the resolution of the data. Source edges are characterized by maxima of the normalized Harris filter (NHF) response. The NHF will be evaluated using three synthetic data sets. The capabilities of NHF to handle noisy data and complex sources are described. Finally, we apply the NHF to real gravity data from the Yili basin.

The Harris Filter
Following [18], the change E in the potential field data f (with a size of m × n) produced by a small shift (u, v) is given by (details are provided in Appendix A) where the local structure matrix M is defined by: The horizontal derivatives A and B, and diagonal derivative C are defined as: and: (12) where n corresponds to the number of neighbors in the four directions (along the row, column, and diagonals) of f within the small shift. The response function R of the Harris filter is defined as: where Det represents the determinant of M and Det(M) = AB − C 2 ; Tr represents the trace of M and Tr(M) = A + B; µ is a user-specified parameter, which defines the sharpness or smoothness of the response function R and µ ∈ [0, 1]. The Harris filter response of the f can be stored in a grid R in the same size as f : A moving window of size 3 × 3 and Equation (13) are employed to compute the elements of the R. The edge/corner areas are determined by the amplitude of the elements in R. If R ij is a local maximum within its neighbors, R ij will be considered as a candidate of edge/corner.

Amplitude Balance
Although the local maxima of R can delineate edges of the source, their amplitudes vary greatly. We employ a multi-step normalization technique to balance the amplitude of the local maxima of R. Thus; all edges are featured by a fixed value of 1.
(i) Pick out the local maxima of R We adopt the method of [23] to pick out the local maxima of R. The method compares each R i,j (except the marginal rows and columns, i ∈ [1, m] and j ∈ [1, n]) and its eight nearest neighbors in four directions (along rows, columns, and two diagonals). If R i,j is greater than its neighbors in at least two directions, it is regarded as a candidate of local maximum [24].
(ii) Threshold of the local maxima Some local maxima can be named as invalid, while they produce artifacts (false edges) in the solutions. Therefore, threshold processing is necessary to suppress the local maxima related to noise. If a local maximum is less than the threshold threshold = λ·max(R) (15) it will be disqualified from the candidate of the local maximum. In theory, λ is limited to between 0 and 1. For different data sets, the optimal value of λ may be different.
(iii) Calculate the upper envelope surface of the R We employ the local maxima above the threshold to interpolate the upper envelope surface, and the natural neighborhood method [25] is used. For the upper envelope surface boundary, we increase their values by 0.1 times the maximum of R. This will most likely guarantee that the upper envelope surface is always above R.
(iv) Normalize R using upper envelope surface The normalized Harris filter (NHF) response is the ratio of R to the upper envelope surface, and its maxima (with a fixed value of 1) represent the source edges.

Synthetic Data
This section applies the proposed approach to three synthetic data sets to evaluate its performance. All derivatives are calculated in the frequency domain. Three scenarios are considered as follows:

The Single Parallelepiped
We have generated a single parallelepiped located at a depth of 20 m to assess the feasibility of the proposed approach. Figure 1a shows the noise-free gravity anomaly (with a grid interval of 10 m) produced by the parallelepiped with a density of 1.0 g/cm 3 , where the black lines indicate the actual edges.

Three Prisms
To further examine the proposed approach's capability, we employ a model consisting of three prisms located at different depths [5].    Figure 1a with different empirical parameters µ. All corners and edges are characterized by local maxima. Fur-thermore, a small µ value allows the Harris filter to detect only the sharpest corners of the source (Figure 1b-d), and as the µ value increases, edges of the source are detected. When µ equals 0.1 and 1.0, the edges of the source can be clearly observed (Figure 1f,g). Figure 1h shows the NHF response. It demonstrates that the NHF response is amplitude balanced, and maxima characterize the edges. In general, when the application of interest is edge detection, high µ is preferred. Based on the above results and the knowledge from [18], we set the value of µ in the following examples to 1, unless otherwise noted.

Three Prisms
To further examine the proposed approach's capability, we employ a model consisting of three prisms located at different depths [5].  Figure 2a shows the magnetic anomalies of three prisms. We chose seven frequently used approaches to compare the results. They are THD [1], tilt angle [2], sun shading [12], THDT [3], theta map [4], TAHG [5], and tilt-ASA [6]. To reduce the long-tailed histogram of THD, we plot the square root of THD in this and the following sections. We assess the proposed approach's performance to detect edges from anomalies with different amplitudes and compare the results with other approaches. Figure 2b-h show THD, tilt angle, THDT, theta map, TAHG, tilt-ASA, and the sun shading response of the data in Figure 2a. Figure 3 presents profile data (the location is shown in Figure 2a with a red line) to provide insightful details. THD is presented by maxima at the edges of the sources, and the strength of the THD is gradually attenuated with depth. The tilt angle equalizes the response from the anomalies at different depths, and the zero values of the tilt angle can be used to trace the edges [26]. However, there are some erroneous edges between different prisms. THDT delineates the edges by maxima, but the amplitude of THDT decreases as the source's depth increases. Also, some wrong edges pollute the map. Theta map produces the edges of sources by maxima; however, for deep sources, the maxima are shifted out from the actual edges, making the sources appear bigger than they are. Similar to THDT, erroneous edges appeared in the results of the theta map. The TAHG balances the response of sources from different depths and shows maxima at the edges of We assess the proposed approach's performance to detect edges from anomalies with different amplitudes and compare the results with other approaches. a red line) to provide insightful details. THD is presented by maxima at the edges of the sources, and the strength of the THD is gradually attenuated with depth. The tilt angle equalizes the response from the anomalies at different depths, and the zero values of the tilt angle can be used to trace the edges [26]. However, there are some erroneous edges between different prisms. THDT delineates the edges by maxima, but the amplitude of THDT decreases as the source's depth increases. Also, some wrong edges pollute the map. Theta map produces the edges of sources by maxima; however, for deep sources, the maxima are shifted out from the actual edges, making the sources appear bigger than they are. Similar to THDT, erroneous edges appeared in the results of the theta map. The TAHG balances the response of sources from different depths and shows maxima at the edges of the sources, as expected for TAHG that uses second-order derivatives. Like the tilt-angle, the tilt-ASA displays the edges with zero values [6]. The sun shading presents the minima at the edges of the sources, and its strength gradually decreases with depth. In the implementation of NHF, we set λ to 10 −3 . It is found that the NHF response clearly and focused characterizes edges of sources at different depths. implementation of NHF, we set λ to 10 −3 . It is found that the NHF response clearly and focused characterizes edges of sources at different depths. Noise is an unavoidable factor in interpreting real field data. Figure 4 shows the robustness of different approaches to noisy data. Figure 4a shows the magnetic anomaly in Figure 2a by adding independent Gaussian noise. Noise is added to data with the following equation:  Noise is an unavoidable factor in interpreting real field data. Figure 4 shows the robustness of different approaches to noisy data. Figure 4a shows the magnetic anomaly in Figure 2a by adding independent Gaussian noise. Noise is added to data with the following equation: with: where d n i is the noisy datum; d i is the noise-free datum; normrnd is the Matlab built-in function, which generates a random number from the normal distribution with mean parameter mu and standard deviation parameter std; d min is the minimum value of noisefree data; d max is the maximum value of noise-free data. The percent is set to 2% in the noisy data test.
where n i d is the noisy datum; di is the noise-free datum; normrnd is the Matlab built-in function, which generates a random number from the normal distribution with mean parameter mu and standard deviation parameter std; dmin is the minimum value of noise-free data; dmax is the maximum value of noise-free data. The percent is set to 2% in the noisy data test. The black rectangles outline the true edges. Figure 4 displays the results of the noisy data test. Figure 5 shows a depth section of the synthetic model and responses of different edge detectors. In general, all results are affected by the noise, and as the source depth increases, this effect becomes more serious. Specifically, Figure 4b shows that THD is less sensitive to noise. The tilt angle map equalizes the anomaly amplitudes, and it is less sensitive to noise than THDT, theta map, TAHG, and tilt-ASA. The edges recognized by the THDT are very blurred. The theta map can only roughly identify the edges of the sources. TAHG and tilt-ASA are particularly sensitive to noise and can almost only reflect the edge of a shallow source. Sun shading is less sensitive to the noise, and the recognition effect deteriorates as the source depth increases. Figure 4i shows the NHF response (λ = 1.2 × 10 −2 ), and its maxima outline the sources' edges. It shows that the proposed approach is capable of handling noisy data.   Figure 5 shows a depth section of the synthetic model and responses of different edge detectors. In general, all results are affected by the noise, and as the source depth increases, this effect becomes more serious. Specifically, Figure 4b shows that THD is less sensitive to noise. The tilt angle map equalizes the anomaly amplitudes, and it is less sensitive to noise than THDT, theta map, TAHG, and tilt-ASA. The edges recognized by the THDT are very blurred. The theta map can only roughly identify the edges of the sources. TAHG and tilt-ASA are particularly sensitive to noise and can almost only reflect the edge of a shallow source. Sun shading is less sensitive to the noise, and the recognition effect deteriorates as the source depth increases. Figure 4i shows the NHF response (λ = 1.2 × 10 −2 ), and its maxima outline the sources' edges. It shows that the proposed approach is capable of handling noisy data. To illustrate the choice of parameters, we change λ for noisy data tests. Figure 6 shows the NHF response of the data in Figure 4a with different parameters. A small λ value can cause erroneous edges related to noise, and a large λ value can suppress the noise effect. In this case, an λ of 1.2 × 10 −2 produces satisfactory results (Figure 4i). This is a benefit of the proposed approach because it can be adapted to handle noise data. To illustrate the choice of parameters, we change λ for noisy data tests. Figure 6 shows the NHF response of the data in Figure 4a with different parameters. A small λ value can cause erroneous edges related to noise, and a large λ value can suppress the noise effect. In this case, an λ of 1.2 × 10 −2 produces satisfactory results (Figure 4i). This is a benefit of the proposed approach because it can be adapted to handle noise data.

The Bishop Model
We use the Bishop model's magnetic anomaly to illustrate the proposed approach's performance in the case of complex sources. Various authors have used the Bishop model to test approaches of estimating source depths using magnetic data [27][28][29]. Figure 7a displays the magnetic anomalies for a basement susceptibility varying from 1.26 × 10 −2 to 1.0 × 10 −1 SI units and an ambient geomagnetic field with strength 50,000 nT, declination 0 • , and inclination 90 • .

The Bishop Model
We use the Bishop model's magnetic anomaly to illustrate the proposed approach's performance in the case of complex sources. Various authors have used the Bishop model to test approaches of estimating source depths using magnetic data [27][28][29]. Figure 7a displays the magnetic anomalies for a basement susceptibility varying from 1.26 × 10 -2 to 1.0 × 10 -1 SI units and an ambient geomagnetic field with strength 50,000 nT, declination 0°, and inclination 90°.

The Bishop Model
We use the Bishop model's magnetic anomaly to illustrate the proposed approach's performance in the case of complex sources. Various authors have used the Bishop model to test approaches of estimating source depths using magnetic data [27][28][29]. Figure 7a displays the magnetic anomalies for a basement susceptibility varying from 1.26 × 10 -2 to 1.0 × 10 -1 SI units and an ambient geomagnetic field with strength 50,000 nT, declination 0°, and inclination 90°. The THD, tilt angle, THDT, theta map, TAHG, tilt-ASA, and sun shading are shown in Figure 7b-h, respectively. It can be seen that all approaches generate the main features of shallow sources. Among them, TAHG provides improved structural details yet contains very little information for deep structures (the southeast part of the figure). We can also find that tilt-ASA shows some striped artifacts in the east and south parts of the figure. Figure 7i-l show the NHF responses with different choices of λ. As λ decreases, we can enhance our understanding of the edges of subsurface structures. Specifically, the NHF response generated by the large λ corresponds to the edges of the main structures, and the small λ can extract the edges associated with the minor structures. The THD, tilt angle, THDT, theta map, TAHG, tilt-ASA, and sun shading are shown in Figure 7b-h, respectively. It can be seen that all approaches generate the main features of shallow sources. Among them, TAHG provides improved structural details yet contains very little information for deep structures (the southeast part of the figure). We can also find that tilt-ASA shows some striped artifacts in the east and south parts of the figure. Figure 7i-l show the NHF responses with different choices of λ. As λ decreases, we can enhance our understanding of the edges of subsurface structures. Specifically, the NHF response generated by the large λ corresponds to the edges of the main structures, and the small λ can extract the edges associated with the minor structures.

Field Gravity Data from Yili Basin
We will now discuss applying the proposed approach to gravity data collected from Xinjiang, northwestern China (Figure 8a). The study area is located in the middle portion of the Yili Basin, centered on the Baishidun sub-high, covering an area of about 1200 km 2 . Tectonic units in the study area include Awulela high, Baishidun sub-high, Gongliu sub-sag, Southern ramp region, and Wusunshan uplift. The geological structure map of the study area is given in Figure 8b. The red lines indicate the structures inferred from the historical geophysical data, while the blue lines present the structures deduced from newly collected gravity data [30]. There is no clear understanding of the structural edges of the study area.
Sunshading using an elevation angle of 90°. (i) NHF response of the data in (a) (with λ = 1.0 × 10 −3 ). (j) NHF response of the data in (a) (with λ = 1.0 × 10 −6 ). (k) NHF response of the data in (a) (with λ = 1.0 × 10 −9 ). (l) NHF response of the data in (a) (with λ = 0). The THD, tilt angle, THDT, theta map, TAHG, tilt-ASA, and sun shading are shown in Figure 7b-h, respectively. It can be seen that all approaches generate the main features of shallow sources. Among them, TAHG provides improved structural details yet contains very little information for deep structures (the southeast part of the figure). We can also find that tilt-ASA shows some striped artifacts in the east and south parts of the figure. Figure 7i-l show the NHF responses with different choices of λ. As λ decreases, we can enhance our understanding of the edges of subsurface structures. Specifically, the NHF response generated by the large λ corresponds to the edges of the main structures, and the small λ can extract the edges associated with the minor structures.

Field Gravity Data from Yili Basin
We will now discuss applying the proposed approach to gravity data collected from Xinjiang, northwestern China (Figure 8a). The study area is located in the middle portion of the Yili Basin, centered on the Baishidun sub-high, covering an area of about 1200 km 2 . Tectonic units in the study area include Awulela high, Baishidun sub-high, Gongliu subsag, Southern ramp region, and Wusunshan uplift. The geological structure map of the study area is given in Figure 8b. The red lines indicate the structures inferred from the historical geophysical data, while the blue lines present the structures deduced from newly collected gravity data [30]. There is no clear understanding of the structural edges of the study area. Simplified geological structure of the study area (after [30]), the red lines and blue lines indicate two possible structural divisions. The gravity data has been gridded in 500 m intervals. Figure 9a shows the Bouguer anomaly of the study area. Similar to the previous sections, we compare the proposed approach with other approaches, and the results are shown in Figure 9b-h. From Figure 9, we can see that these approaches enhance the main structures. However, it shows that the TAHG and tilt-ASA carry severe striped artifacts because they all use high-order derivatives, and the calculation of high-order derivatives is easily affected by the conflicting information contained in the actual data [32]. The proposed approach (with λ = 10 −5 ) produces focused results ( Figure 9i) and clarifies the structures.
anomaly of the study area. Similar to the previous sections, we compare the proposed approach with other approaches, and the results are shown in Figure 9b-h. From Figure 9, we can see that these approaches enhance the main structures. However, it shows that the TAHG and tilt-ASA carry severe striped artifacts because they all use high-order derivatives, and the calculation of high-order derivatives is easily affected by the conflicting information contained in the actual data [32]. The proposed approach (with λ = 10 −5 ) produces focused results (Figure 9i) and clarifies the structures.  When comparing the NHF response (Figure 9i) with known geological structures (Figure 8b), edges of the Awulela high, the Wusunshan uplift, the Gongliu sub-sag, and the Southern ramp region are consistent with the blue lines in Figure 8b. Some minor structures (white solid ellipses in Figure 9i) in the Wusunshan uplift and Baishidun sub-high are compatible with red lines in Figure 8b. These results suggest that edge enhancement can bridge the gap between different prior information. Structures in the Baishidun sub-high (white dashed ellipses in Figure 9i) are identified in the NHF response but not marked on the geological map ( Figure 8b); this means that edge enhancement can help the interpreter enrich the knowledge of subsurface structures.

Conclusions
In this study, the applicability of NHF to potential field data is shown to enhance the edges of the subsurface structures. The edges are identified by fixed value 1 in the NHF response. The NHF involves an adjustable parameter λ. The choice of λ depends on the quality of the potential field data. When data quality is poor, a higher λ is required to obtain an interpretable result (Figures 4 and 6). When the data quality is acceptable, λ can be adjusted to enhance the understanding of edges from different structures (Figure 7). In this case, NHF can be regarded as a gain-controlled edge detector.
The proposed approach has been tested using three synthetic data sets. For the single parallelepiped example, we have demonstrated that setting the empirical parameter µ to 1 can improve the sensitivity of NHF to the source edges. For three prisms example, the effectiveness of the NHF for sources at different depths is illustrated. Furthermore, the results show that a larger λ can aid the NHF against noise. For magnetic data caused by the Bishop model, we present the capability of NHF as a gain control edge detector. In the case of complex sources, NHF can enhance our understanding of the edges of different sources. Finally, the proposed approach was applied to the gravity data collected in the Yili Basin. The results are compared with previous studies, and it can be seen that edge enhancement is a valuable tool for characterizing structural edges.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The data underlying this article will be shared on reasonable request to the corresponding author. The code associated with this paper can be found at https://github. com/GeoGoku/EdgeEnhancement (accessed on 3 June 2021). The users can redistribute it and/or modify it under the terms of the GNU General Public License.