Multi-Scale Edge Detection Method for Potential Field Data Based on Two-Dimensional Variation Mode Decomposition and Mathematical Morphology

Using observational data to determine the edges of the sources is an important task in the interpretation of potential field data. Extracting the edges of deep and shallow bodies effectively is the key to correctly understanding the underground structure. Based on the good multi-scale decomposition ability of two-dimensional variational mode decomposition (2D-VMD) and the outstanding shape analysis capability of mathematical morphology (MM), a new multi-scale edge detection method for potential field data is proposed. We propose using the variance of this morphological filter as a basis for selecting the optimal structural element (SE) scale. By establishing theoretical models and comparing the results of our method with those of traditional edge detection methods, the proposed method is shown to be effective at detecting edges within potential field data. Taking the Hanmiao area of Chifeng city, Inner Mongolia, China, as an example, 1:50000 aeromagnetic data are processed and analysed by this method. The physical properties of the rocks in the study area are also discussed. The results of theoretical calculations and real data processing show that this method can accurately extract the edges of the sources at different scales. And the real data processing results show that this method is suitable for the identification of structural faults.


I. INTRODUCTION
The geological process is long and complex. Therefore, potential field data are composed of superimposed sources having different depths, shapes, sizes, densities, and magnetism. Shallow-source anomalies with small distribution ranges are called residual anomalies. However, deep-source anomalies with large distribution ranges are called regional anomalies. Separating the potential field from these anomalies and detecting the edges of different sources are important tasks in geophysics.
For many years, scholars have developed a variety of separation methods and edge detection methods for potential field data. Separation methods include matched filtering [1], The associate editor coordinating the review of this manuscript and approving it for publication was Geng-Ming Jiang . the minimum curvature technique [2], [3], analytical continuation [4], and wavelet multi-scale decomposition [5]- [7]. Edge detection methods mainly fall into two categories: mathematical statistics and derivative analysis. The mathematical statistics approaches include Euler deconvolution [8], [9], small subdomain filtering [10], [11], and the normalized standard deviation method (NSTD) [12]. Derivative analysis techniques include the vertical derivative (VDR) [13], tilt angle (Tilt) [14], total horizontal derivative (THDR) [15], and Theta map [16], among others. However, these methods do not completely solve the abovementioned problems, so many scholars have continued to search for new separation and edge detection methods.
N.E. Huang et al. proposed empirical mode decomposition (EMD) in 1998 [17], and Nunes et al. developed bidimensional empirical mode decomposition (BEMD) on the basis of EMD [18]. These methods have been widely used in nonlinear system analysis, earthquake engineering, non-destructive detection, meteorology, biomechanics and many other fields [19], [20]. However, EMD lacks a strict mathematical basis and is prone to problems such as mode confusion, end effects and a low algorithmic efficiency. Dragomiretskiy et al. proposed a new mode decomposition method, variational mode decomposition (VMD), in 2014 and developed 2D-VMD on this basis [21]- [23]. As a non-recursive and self-adaptive variational method, this method is completely different from traditional EMD. VMD, which is based on the Wiener filter and variation method, has a strong mathematical basis and can quickly and efficiently separate several intrinsic mode functions (IMFs) with different frequencies at the same time.
MM is a mathematical tool for analysing image structures and morphologies based on random set theory and integral geometry theory. In 1964, the French scholars Matheron and Serra developed MM. The publication of Image Analysis and Mathematical Morphology in 1982 laid a mature theoretical foundation for this method [24]. MM boasts a fast calculation speed and plays crucial roles in pattern recognition, image filtering, edge detection and signal processing. The application of MM to the geosciences has also developed rapidly. At present, MM is employed mainly in the following regards: seismic signal de-noising [25], seismic crack detection [26], and gravity and magnetic data processing [27].
The common feature of the edge extracting algorithms is the highlighting of source edges by mathematical means. With edge enhancement, source edges are easier to identify than in the original potential field. If the edge can be located automatically on the basis of the above, then the accuracy of interpretation can be improved. For THDR, Theta map and NSTD, the local maximum corresponds roughly to the edge of the geological body, such that the boundary analysis technology proposed by Blakely and Simpson [28] can be used to locate the edge. For cases in which the source edge is not located at the abnormal maximum value, but at the abnormal zero value point or amplitude inflection point (which is also the location of the maximum gradient), such as VDR and Tilt, the Canny operator [29] can be used for edge location. For example, Tilt is enhanced by the vertical first derivative. Although its edge position roughly corresponds to the zero value, nonetheless its gradient is the local maximum, which satisfies the requirement of the Canny operator for image edge localization.
In this paper, 2D-VMD and MM are combined to solve the problem that traditional methods are not effective at extracting the edges of the sources from superimposed field data. The multi-scale decomposition of potential field data is carried out by 2D-VMD, and edge detection operators are constructed by MM to obtain the edges of the sources at different scales. Next, we use the boundary analysis technique proposed by Blakely and Simpson to extract the edge position of each scale. Several models verify the accuracy and stability of this method in multi-scale edge detection tasks. The proposed method is also applied to process real data from the Hanmiao area of Chifeng city, Inner Mongolia. The results show that this method can effectively detect the edges of faults, intrusive rocks and volcanic basins and can provide rich information for the division of geological units and metallogenic prognosis.

II. PRINCIPLE OF THE METHOD A. 2D-VMD ALGORITHM
The 2D-VMD algorithm assumes that each IMF has a different central frequency and limited bandwidth. By updating the bandwidth of each mode, the sum of the estimated bandwidths of all IMFs is minimized, and the central frequency of each mode is finally obtained. The 2D-VMD function is defined as a two-dimensional analytical signal. The decomposed form of the two-dimensional signal f (x) is [21]: To find the optimal solution of the above problem, the constrained variational problem is transformed into an unconstrained variational problem [30]. The penalty parameter α k and the Lagrangian multiplier λ(x) are introduced, and the augmented Lagrangian function is as follows: The alternate direction method of multipliers (ADMM) is used to find the saddle point of the above equation, and the optimal solution of the objective function is obtained by updating u n+1 k , ω n+1 k , and λ n+1 k in an alternating fashion [31]. The expression of u n+1 k is: By using the Parseval/Plancherel Fourier equidistant transformation, equation (3) is transformed to frequency domain, and the solution of the quadratic optimization problem is VOLUME 8, 2020 obtained [32]: where ω ∈ k , k = {ω |ω · ω k ≥ 0 }, and ω n+1 k are the centres of gravity of the power spectra of the different modes.
Step 4: Given the discriminant accuracy ε >0, stop iterat- For the 2D-VMD algorithm, the number of IMFs needs to be set in advance. For different numbers of IMFs, the decomposition results will vary. For a potential field signal, the slope of the radially averaged log power spectrum is inversely proportional to the average buried depth of the source's centre of gravity. Therefore, the larger the slope of the fitting line, the deeper the geological body is buried; the smaller the slope of the fitting line, the shallower the geological body is buried. According to the number of fitting lines with a discernible slope, it is possible to determine how many different equivalent source layers are present and then to obtain the number of IMFs.

B. EDGE ENHANCEMENT ALGORITHM 1) BASIC OPERATION OF MM
The VDR approach has been shown to be an effective edge recognition tool. However, the VDR increases not only the abnormal amplitude but also the interference of noise, which makes the output result unstable. In this paper, the VDR of each IMF is calculated, and then the VDR is processed by MM. The purpose of this is to remove noise and recognize the edge of the field source. MM uses a SE with a certain shape and size to detect the target signal. By continuously translating the SE and matching and modifying the target signal, the desired results, that is, de-noising, image enhancement and feature extraction, can be achieved. The basic operations of MM include dilation, erosion, opening and closing [24].
Assuming that f (x, y) is the VDR of an IMF, g(i, j) is the SE, x and y are the positions of the calculation points, ⊕ denotes dilation, and denotes erosion, then the dilation and erosion operations are defined as: Opening and closing are combined operations comprising both dilation and erosion, which can smooth the image contours. Among them, the operation process of opening is erosion first and then dilation to remove small amounts of image noise, whereas the operation process of closing is dilation first and then erosion to fill miniscule image holes. The opening and closing operations are defined as follows:

2) MULTI-SCALE SE
The scale and shape of the SE have important influences on the calculation results, similar to the filtering window in digital signal processing. The scale of the SE is closely related to the edge width. The de-noising ability of a small-scale SE is weak, but the edge details can be detected. In contrast, largescale SEs have strong de-noising abilities, but the detected edges are relatively rough. Commonly used SEs are disks, squares, diamonds, lines, etc. The specific shape of the SE is usually determined according to the geometric properties of the target geological feature. To ensure a good filtering effect, the SE should be simplified as much as possible. Therefore, after conducting repeated experiments, squares are selected in this study.
Because geological processes are complex, potential field images contain characteristics on different scales. If the image is processed by an SE with a single scale, the geological details will be lost. Therefore, it is necessary to use SEs of varying scale to process the potential field data and obtain detailed information about geological features at different scales. The multi-scale SE is defined as follows [33]: (11) where g denotes the SE, n is the scale parameter, and g dilates itself n − 1 times to obtain ng. In essence, the multi-scale mathematical morphology technique uses SEs to transform an image into a series of filtered images. The shape of the SE does not change, but the scale increases.

3) PROPOSED MORPHOLOGICAL FILTER
Upon improving the morphological operator, the opening and closing operations are combined to achieve de-noising, while dilation and erosion are used for edge detection. A combined filter is constructed by combining the opening and closing operations to eliminate the noise in potential field data that does not match the scale, which constitutes the basis of the edge detection operation [34]: To avoid statistical bias in the use of the opening-closing or closing-opening morphological filters, the two are combined and averaged [34]: On this basis, the gradient operators E 1 (x, y) and E 2 (x, y) are constructed as follows: The traditional MM edge detection operator can eliminate confounding factors from potential field data. However, the recognized edge of the source may deviate from the true position. In this paper, we propose an edge detection method based on MM by improving E 1 (x, y) and E 2 (x, y). We suggest using a filter based on morphological edge enhancement (MEE) to increase the accuracy of the edge detection results.
where E max (x, y) is the maximum value of E 1 (x, y) and E 2 (x, y), E min (x, y) is the minimum value of E 1 (x, y) and E 2 (x, y). E max (x, y) and E min (x, y) can prevent a zero value from being in the denominator. The edges of sources in potential field data can be identified by the maximum value produced by MEE.

4) OPTIMAL SE SCALE
We know that the edge detection results of the MEE filter with different SE scales are inconsistent. Thus, the method for determining the most suitable SE scale is a key problem that we address in this paper. A large number of experiments show that each scale has an optimal SE scale. Variance is an effective tool for measuring the discreteness of data. High variance means that the data are relatively discrete; low variance means that the data are relatively centralized. There is an obvious maximum along the relation curve between the SE scale and the MEE variance (which will be explained later). Here, the scale corresponding to the maximum variance, denoted max (s 2 ), is defined as the optimal scale, and the MEE value corresponding to this scale is called the optimal MEE. The VDR of an IMF is defined as V (x, y), and the SE ng of scale n is selected to calculate the MEE (n) (x, y) of V (x, y). The variance s 2 n of MEE (n) (x, y) is: where MEE (n) (x, y) is the MEE value of scale n, MEE (n) (x, y) is the sample mean of MEE (n) (x, y), and M and N are the number of lines and the number of points, respectively, on each line.

C. PROPOSED MULTI-SCALE EDGE DETECTION PROCESS
The 2D-VMD algorithm has good anti-noise performance, but its ability to filter a pulse signal is not as good as MM. According to IMFs with different scales, using the filtering characteristics of the SE, an edge detection filter based on 2D-VMD and MM is established. The specific implementation process is as follows: Step 1: The radially averaged log power spectrum of potential field data f (x, y) is obtained, and the number k of IMFs is determined according to the number of fitting lines with different slopes in the log power spectrum.
Step 2: To reduce the end effect, a cosine transform is used to extend the boundary of the potential field, and then 2D-VMD is applied to decompose f (x, y) into k IMFs.
Step 3: The VDR of each IMF is calculated, and the SE with the appropriate shape is selected. Then, the relation curve between the SE scale ng and the MEE variance s 2 n is drawn. The scale corresponding to the maximum point on this curve is used as the optimal scale of the edge detection filter for each IMF.
Step 4: MEE is employed for each IMF using the optimal SE scale to obtain source edge information at different scales.
Step 5: Select the Blakely rule to locate the edge automatically. The basic principle of this method is to compare the points around each point of interest. If the point of interest is larger than the surrounding points, then the point will be retained. Each point will be detected in turn through the sliding window and will be retained, or not, through the test standard. If necessary, these points should be connected into a line by using certain rules (the minimum distance between adjacent points and the minimum number of points required). After which point, the multi-scale edge detection of gravity and magnetic data has been completed.

III. MODEL TESTS A. THEORETICAL MODEL
In order to verify the effectiveness of the proposed edge detection method, a synthetic model of seven vertical cubes with different depths, sizes and residual densities is established. The starting point coordinate is 0.0 km, and the ending point coordinate in both directions is 20.0 km. Fig.1 shows the projected distributions of the cubes on a horizontal plan view and a vertical profile. The geometric parameters, densities and magnetism of the cubes are shown in Table 1. Cubes A1 generate gravity anomaly G RegA , cubes B1 and B2 generate gravity anomaly G ResB , cubes C1 and C2 generate gravity anomaly G ResC , and cubes D1 and D2 generate gravity anomaly G ResD.
The model tests are carried out with MATLAB R2016b. The synthetic model is designed with three layers, there is little difference in the buried depth between model C and model D. Therefore, models C and D are regarded as a

B. RESULTS OF 2D-VMD
From the spectrum analysis, the frequency components of the regional field and residual field are different. As a result, the radially averaged log power spectrum can be used to estimate the depth of the source [35]. When the potential field data are superposed by the sources of different frequencies, a number of line segments with different slopes can be fitted from the power spectrum curve. At present, there are many kinds of line fitting methods, such as graphical method, the least square method and so on. The graphic method is relatively large in terms of human factors, and the positions and number of line segments divided by different interpreters will vary. Because the linear least square method is more accurate and objective, we use it in this work. The fitting process entails transforming the nonlinear relationship into a linear relationship. If the line fitting is not accurate, then the relative error will increase significantly. The standard for optimal fitting can be determined as the minimum sum of squared errors [36]. In practical application, where either method is an estimation of the number and depth of the sources, but the spectrum analysis still has a good reference value.
The first synthetic model consists of model B, model C and model D. Through the radially averaged log power spectrum of the synthetic gravity anomaly shown in Fig. 3(a), two frequency bands can be clearly seen, and thus, the number of IMFs is determined to be two. The slopes of the two fitting lines are k1= −5.97 and k2= −0.81, respectively, with a difference of 7.4×. This result shows that there are obvious frequency differences. Then, 2D-VMD is used to decompose the potential field data into two layers, as shown in Fig. 4.
The second synthetic model consists of model B, model C and model D with Gaussian noise and can be divided into three equivalent layers, as shown in Fig. 3(b). The slopes of the fitting lines are k1= −4.78, k2= −1.16 and k3= −0.12.
With the method proposed in this paper, IMF 1 , IMF 2 and IMF 3 are obtained by 2D-VMD through spectral analysis, as shown in Fig. 5. IMF 1 is used to approximate the regional anomaly, IMF 2 is used to approximate the residual anomaly, and IMF 3 is used for noise removal.
The third synthetic model consists of model A, model B, model C and model D and is established to test the multi-scale edge detection capability of the proposed method. As shown in Fig. 3(c), the radially averaged log power spectrum can be divided into three equivalent layers. The slopes of the fitting lines are k1= −12.73, k2= −2.19 and k3= −0.2, respectively. 2D-VMD is used to decompose the potential field data into three layers, as shown in Fig. 6. From these figures, it can be seen that as the number of decomposition layers increases, the frequency aliasing between IMFs will increase. Therefore, we must use MEE or other edge enhancement methods to suppress the noise and enhance the edges of IMFs. This approach will allow us to reduce the frequency aliasing by as much as possible and highlight the edges of sources at different scales.

C. RATIONALITY ANALYSIS OF THE STRUCTURAL ELEMENT 1) DETERMINATION THE MEE VARIANCE
The purpose of MEE is to detect the edges of the sources. The MEE value fluctuates greatly near such edges. This method considers the SE scale. If the scale is different, the MEE value will be different, and thus, the MEE variance will be significantly different. The more appropriate the scale is, the more obvious the edge will be, which directly reflects an increase in variance. Therefore, variance can effectively distinguish the MEE calculation effect.
The proposed method is used to analyse IMF 1 and IMF 2 in Fig. 4. First, the MEE values of IMF 1 and IMF 2 corresponding to different SE scales are obtained, and the variances of the different MEE values are calculated on this basis. These parameters can be analysed according to the relation curve between the SE scale and MEE variance. As seen from Figs. 7(a) and 7(b), the variance increases at   first and then decreases with an increase in the SE scale; the scales corresponding to the maximum values are 18 and 3. Therefore, the optimal sizes of the SEs for IMF 1 and IMF 2 are 18 × 18 and 3 × 3, respectively. This is also equivalent to their window sizes. The same method is used to analyse IMF 1 and IMF 2 in Fig. 5. According to the relation curves in Figs. 7(c) and 7(d), the optimal sizes of the SEs for IMF 1 and IMF 2 are determined to be 20 × 20 and 3 × 3, respectively. Similarly, as shown in Figs. 7(e) to 7(g), the optimal size of the SE for IMFs in Fig. 6 are 46 × 46, 14 × 14 and 3 × 3.

2) SE SELECTION
To verify the influence of the SE scale on the MEE value, the IMFs are analysed by using different SE scales.     Figure 4; (c) and (d) IMF 1 and IMF 2 of the gravity anomaly with Gaussian noise in Figure 5; (e), (f) and (g) IMF 2 , IMF 3 and IMF 1 of the gravity anomaly in Figure 6.
resolving power of small-scale SEs for such sources is stronger than that of large-scale SEs. As the SE scale increases, the edge recognition ability of MEE decreases gradually. Therefore, small-scale SEs are more suitable for this kind of shallow-buried models.
To further prove the influence of different SE scales on the edge extraction results, the IMF 1 in Fig. 6(a) is taken as an example, and the SEs of different sizes are used for edge enhancement. The comparison is shown in Figs. 8(q) -8(t). It can be seen that for a large-scale IMF, as the SE scale is increased, the filtering effect of the SE scale on small-scale geological bodies is continuously enhanced. When the SE scale reaches the optimal SE scale in Fig. 7(e), the edge increment effect is optimal.
On the basis of the above analysis, the profiles in Figs Fig. 2(a), and the results are shown in Fig. 9. As seen from these figures, when the SE scale is small, the width of the MEE high-value area is narrow. There are also large amplitudes at the edge of the shallow-buried model. As the size of the SE increases, the MEE amplitudes at the edges of model B1 and model B2 with deeper burial depths are obviously large, as shown in Fig. 9(d). In the same case, as shown in Fig.9(p), for model A in Fig. 6, when the size of the SE increases to 46 × 46, the IMF can achieve the best edge enhancement results. This shows that using optimal SE scales can enable the edges of the sources with different depths and sizes  Figure 4; (e)-(h) IMF 2 of the gravity anomaly in Figure 4; (i)-(l) IMF 1 of the gravity anomaly with Gaussian noise in Figure 5; (m)-(p) IMF 2 of the gravity anomaly with Gaussian noise in Figure 5; (q)-(t) IMF 1 of the gravity anomaly in Figure 6 (the black boxes are the locations of the models).
to be better recognized. Moreover, the detection of the residual field is best suited to small-scale SEs, while the detection of the regional field is best suited to large-scale SEs. The above comparative analysis thoroughly confirms that the edge detection method proposed in this paper not only has high accuracy but also has a stable calculation process.

D. EDGE DETECTION EFFECT OF DIFFERENT METHODS
In this section, to verify the effectiveness of the proposed method, we compare the effect of different potential field separation methods and different edge detection methods on edge detection accuracy. Therefore, 2D-VMD, BEMD and wavelet multi-scale decomposition are selected to separate the gravity anomaly data of the synthetic models. The edge detection effect of the separation field is compared by using MEE, Tilt and Theta map. For MEE and Theta map, we use the Blakely rule to locate the edges, and then, we overlay the results of significance level N=2 on the anomaly map with blue lines. We also use the Canny operator to locate the edges of Tilt. With this method, the edges are also overlaid on the anomaly map with blue lines. Then, by comparing the distance between each detected edge position and its corresponding true edge position, the edge detection accuracy of the different methods is verified.

1) MODEL ANALYSIS
We take the first synthetic model, which consists of model B, model C and model D, as shown in Table 1, as an example to compare the edge detection effect of the different methods. First, 2D-VMD is used to separate the gravity anomaly data of this synthetic model. As shown in Fig. 4, each IMF is more or less mixed with information from other scale sources. The regional field anomaly caused by model B is significantly larger than the residual field anomaly caused by model C and model D, so there is still a small amount of regional field anomaly residue in small-scale IMF 2 .
Next, we compare the edge enhancement effects of MEE and Tilt on the same IMFs. In Fig. 10(a), the MEE of IMF 2 mainly enhances the edges of model C and model D, whereas in Fig. 10  located by using the Canny operator is more accurate than that determined by the zero value. However, Tilt can amplify frequency aliasing while highlighting weak anomalies. Fig. 10(d) shows that while the Tilt method highlights the edges of deep sources, it also identifies the edges of shallow sources.
To compare the influence of different potential field separation methods on the edge detection results, BEMD, which also incorporates multi-scale decomposition, is introduced. A total of two components, IMF 1 and RES, are obtained by BEMD. IMF 1 is used to approximate the residual anomaly, and RES is used to approximate the regional anomaly. The edges of IMF 1 and RES are detected by the Theta map method. As shown in Figs. 10(e) and 10(f), the results display obvious frequency aliasing, and there are false edges.
Furthermore, the anti-noise effects of the different methods are compared. Gaussian noise is added to the gravity anomaly of the synthetic model, which consists of model B, model C and model D, and the edge detection is carried out by each of the above methods. The decomposition results of 2D-VMD for the gravity anomaly of this synthetic model are shown in Fig. 5. Figs. 11(a) and 11(b) show that the edges of IMF 2 and IMF 1 are detected by MEE, and the detection effect of the proposed method is good. Figs. 11(c) and 11(d) show the edge detection result obtained by the Tilt method. To eliminate the influence of noise, the data should be filtered before performing edge detection. In this paper, upward continuation is used for filtering, and the upward continuation height is determined to be 0.5 km through experiments. It can be seen from these two figures that the Tilt method are still greatly affected by noise after filtering. We use the Canny operator's non-maximum suppression ability and set the threshold to 0.6 in MATLAB to reduce the edge error. After setting the threshold, the detected edge is more accurate. However, false edges are still identified, which will significantly interfere with the recognition of true edges.
Five IMFs and one residue (RES) are obtained by BEMD. The highest-frequency signal, IMF 1 , is treated as noise, and IMF 2 is used to approximate residual anomaly G Res , while the sum of IMF 3 , IMF 4 , IMF 5 and RES is used to approximate regional anomaly G Reg . The Theta map method is also used to recognize the edges of G Res and G Reg . As shown in Figs. 11(e) and 11(f), BEMD also has a frequency aliasing problem, and the Theta map edge detection algorithm is sensitive to noise.
Through the analysis of the above two synthetic models, we can see that the edge detection effect of MEE is better than that of Tilt for the same IMFs. In addition, the edge detection effect based on 2D-VMD and MEE is better than that based with BEMD and Theta map. In the above model tests, there are only two layers. Therefore, we need to build a three-layer model for testing to illustrate that the proposed method is a multi-scale edge detection method. The third synthetic model consists of model A, model B, model C and model D. According to the power spectrum, we divide the gravity anomaly of this model into three layers, and then decompose the anomaly into three IMFs by VOLUME 8, 2020  2D-VMD. This process has been described in detail earlier, and the results are shown in Fig. 6. Through MEE variance, the optimal sizes of the SEs for IMF1, IMF2 and IMF3 are 46 × 46, 14 × 14 and 3 × 3, respectively. Using MEE to detect the edge of IMFs has a good effect, as shown in Figs. 12(a) -12(c). For comparison, we use wavelet multi-scale decomposition to decompose the gravity anomaly of this model. After repeated experiments, we take n = 4 to decompose the gravity anomaly. The sum of the first-order and secondorder wavelet details is called the shallow anomaly. The sum of the third-order and fourth-order wavelet details is called the deep anomaly. The fourth-order wavelet approximation is taken as the regional field anomaly. We still use the same SE sizes as above, in addition to using MEE for edge detection with the results of the wavelet multi-scale decomposition. The results of the edge detection are shown in Figs. 12(d) -12(f). It can be seen that the methods for potential field separation are different; that is, even if the same edge detection method is used, the results are different. Compared with wavelet multi-scale decomposition, 2D-VMD is more suitable for multi-scale edge detection with MEE.
To further compare the edge detection results, the main profiles are extracted. The profile data are derived from Fig.10 and Fig. 11, and the results are shown in Fig. 13 and Fig. 14. Comparing the results shown in these figures, the method proposed in this paper is less affected by noise and can recognize the main edges of sources at multiple scales. Moreover, the proposed method can detect the edges of a model with different depths. When the potential field data contain noise, there is no need for filtering before using the approach presented herein. By contrast, it can be seen that this method can more accurately recognize the true edges of the source without increasing the errors at the edges.

2) COMPARATIVE OF EDGE DETECTION ACCURACY
In the above experiments, the edge detection results of the different methods are qualitatively compared. Next, we perform a quantitative accuracy comparison of the edges that have been detected by the different methods, as shown in Figs. 10 -12. We use the true position of each model in the profiles as a reference. The plan positions of these profiles are shown in Fig. 2. The location of data in Table 2 is shown in Figs. 13 and 14. First, we calculate the distance between the extracted edge position and the true edge position. The edge detection accuracy of the different methods is compared based on the distance between these values. We take the edge detection results of model B2 in the first synthetic model in Table 2 as an example. In the profile x=7.7 km, the first true edge position of model B2 is at 14 km. As shown in Fig. 10(b), the edge position detected by 2D-VMD and MEE is at 13.9 km. The distance between these values is 0.1 km. As shown in Fig. 10(d), the edge position detected by 2D-VMD and Tilt is at 13.7 km. Therefore, the distance of this value from the true edge is 0.3 km. Moreover, Tilt will produce false edges. It can also be seen from a comparison of the other models in Table 2 that based on the separation field of 2D-VMD, the edge detection accuracy of MEE is better than that of Tilt.
As also shown in Table 2, the edge detection effect of the proposed method is significantly better than that of BEMD and Theta map in Fig. 10(e). For data with Gaussian noise, MEE is used for the edge detection of model B1. The distance differences between the two detected edge positions and the true edge positions are both 0.3 km, whereas those of model B2 are 0.2 km and 0.3 km. This result shows that the MEE method can still accurately estimate the edge position of abnormal bodies after the addition of noise. Comparing the edge detection results before and after the addition of Gaussian noise, the difference between them is small. This shows that the method of MEE based on SEs scale optimization to improve the edge accuracy has a strong anti-noise ability.    From Table 3, it can also be seen that when the source is three layers, for model A1, in the profile x=7.7 km, the distance differences between the edge positions detected by the proposed method and the true edge position are 0.3 km and 0.3 km. The accuracy of edge detection is better than that of wavelet multiscale decomposition and MEE.
Through the above comparison, we know that the edge detection accuracy of the proposed method is not the highest in all positions. But on the whole, compared with other methods, the proposed method has a higher edge detection accuracy. Therefore, the proposed method can accurately extract the edges of the sources at different scales.

E. INFLUENCE OF THE MAGNETIZATION DIRECTION ON MEE
To study the influence of magnetic inclination on the MEE algorithm, model B in Table 1 is taken as an example to apply MEE with an SE size of 19 × 19, and the profiles with x=7.7 km are extracted, as shown in Fig. 15. Fig. 15 Fig. 15(f).
As seen from these figures, MEE can effectively enhance the edges of the source. However, in the case of magnetic inclination, there is a deviation between the recognized edge and the true edge. When the inclination is 90 • , the recognized edge corresponds to the true edge. Therefore, the proposed method can directly recognize the edges of gravity anomalies. However, in the interpretation of the magnetic anomalies, due to the influence of the magnetization direction, the magnetic anomalies need to be reduced to the pole (RTP) first.

IV. AEROMAGNETIC EDGE DETECTION IN THE HANMIAO AREA, INNER MONGOLIA, CHINA A. GEOLOGICAL BACKGROUND
The Hanmiao area is located in the eastern segment of the Central Asian Orogenic Belt, which is sandwiched between the North China plate and Siberian plate [39], as shown in Fig. 16. The west side of the area is a folded fault belt on the eastern slope of the Da Hinggan Mountains, and the east side of the area is the Songliao basin. The study area has a two-layer structure composed of late Paleozoic basement and Yanshanian caprock with developed Cenozoic, Mesozoic and Paleozoic strata.
The intrusive rocks in the Hanmiao area are mainly the products of Early Cretaceous magmatism. The diverse lithological composition is related to the origin of metal The tectonic traces in the study area originated in the middle Permian. Since the Mesozoic, the regional tectonic framework has undergone substantial changes due to intense tectonic movements and extensive magmatism, and a series of transtensional faults and graben-style basins have formed.The overall structures in the area are distributed NE-SW, and the Mesozoic volcanic rocks that outcrop in large areas are controlled by NE-directed tectonism. The E-W faults are basement faults. Due to several subsequent tectonic events, surface tectonic traces are not evident.

B. AEROMAGNETIC DATA
Aeromagnetic data in the study area with a resolution of 1:50000 are selected as an example to verify the effectiveness of the proposed method. The data area is bounded by 119 • 45'00"E-120 • 00'00"E, 44 • 30'00"N-44 • 40' 00"N; the magnetic inclination is 62.35 • , and the magnetic declination is −8.75 • . The aeromagnetic data are RTP first and then interpolated to a 100 m×100 m grid, as shown in Fig.17(a).
The NE-SW anomaly in the central part of the Hanmiao area reflects the distribution of Mesozoic volcanic basins. This anomaly is narrow in the NE and wide in the SW with a length of 10 km. The main anomaly area is near the Wote and Bayanchagan volcanic basins.
In the northwestern extent of the study area near Manhatu, there is an aeromagnetic anomaly belt distributed NNE-SSW with an intensity of 250 nT. In the southern part of the area near Baiyinhuo, the magnetic anomaly is between 100 and 300 nT with a maximum of 400 nT. The positive magnetic anomaly in this area is large, and the amplitude of the anomaly is obviously higher than that of the background field. Through a field verification, it can be seen that the lithology is mainly acidic volcanic rocks from the Jurassic Manketou'ebo Formation.

C. 2D-VMD RESULT AND DEPTH ESTIMATION
First, the RTP aeromagnetic anomaly data f (x,y) are used to calculate the radially averaged log power spectrum, as shown in Fig. 18(a). According to the characteristics of this curve, f (x,y) can be divided into three components. Therefore, we use 2D-VMD to decompose f (x,y) into IMF 1 , IMF 2 and IMF 3 , as shown in Figs. 17(b) -17(d). Then, the IMFs are used to estimate the depths of the equivalent source layers, as shown in Figs. 18(b) -18(d).
From Fig. 17(b), we can see that IMF 1 is a large-scale component, reflecting the characteristics of the magnetic anomaly; the depth of the equivalent source layer is 2920 m. The anomaly area is large, and the variation in the anomaly is relatively gentle, reflecting the characteristics of the regional magnetic anomaly in the study area. We speculate that the sources distributed NW-SE in the study area are the main cause of the magnetic anomaly.
As shown in Figs. 17(c) and 17(d), IMF 2 and IMF 3 are mainly sub-high-frequency and high-frequency residual magnetic anomalies, and the scale of IMF 2 is smaller than that of IMF 1 . The anomaly characteristics of IMF 2 are as follows: the gradient belts with NE and NW trends surround many circular magnetic anomalies of different sizes. The positive anomalies are distributed mainly in the NW and SE of the study area, which is consistent with the large areas of igneous rocks on the surface. The depth of the equivalent source layer of IMF 2 is approximately 830 m, and the local fractures basically extend to this point. IMF 3 is a small-scale component, reflecting the anomaly characteristics of the magnetic anomalies with an equivalent source layer depth of 710 m; this anomaly is caused mainly by shallow geological bodies and local fractures.

D. MULTI-SCALE EDGE DETECTION RESULTS
In this case, the SE shape is a square. As Fig. 19 shows, through the relation curve between the SE scale and MEE variance, the optimal SE sizes for IMF 1  The results show that the edge of the RTP aeromagnetic anomaly is in good agreement with the boundary of the measured geological body. In detail, this method is able to recognize relatively complete edges, and the edges are relatively convergent. In general, the edges of IMF 1 basically reflect the large-scale, deep-cutting and concealed geological tectonic boundaries throughout the region. The edges of IMF 2 and IMF 3 reflect some shallow tectonic boundaries that are small in scale and do not extend to great depth. These linear traces reflect the trend of tectonic events and their changes at different depths. The blocks delineated by these edges represent subsurface geological bodies with different physical properties.
We use the Tilt method to detect the edges of the same IMFs and then compare the results obtained with Tilt against those of the proposed method. Here, the Canny operator is used for edge location, the results of which are shown in Figs. 20(d) -20(f). Comparing the results of the two methods, it can be seen that Tilt, as a balanced edge enhancement algorithm, can enhance weak anomalies and highlight deep and shallow sources. However, in this case, Tilt is based on the separated field data, which is different from the existing Tilt method, which is based on the original noise data. Theoretical models have verified that each IMF inevitably retains other scale information. Tilt is sensitive to noise and can detect more edges than are actually present.
The edge extracted by MEE is less than that of Tilt, which raises the question as to why. After analysis, we find that in the proposed method, whether the edge can be detected depends on the SE size. If the IMF is a high-frequency component, then a smaller SE scale is needed. The frequency of the large, low and gentle anomaly is obviously different from that of this IMF, so it is considered as confusion. As a result, the small-scale SE usually cannot detect the edge of low and gentle anomaly. Therefore, for a complete comparison, it is also necessary to determine whether the proposed method detect low and gentle anomalies or indistinct weak anomalies. Previously, we have proved through the model that the proposed method uses 2D-VMD to separate the potential field data and highlight weak anomalies of different scales in this process. Therefore, low-frequency IMFs use large-scale SEs to eliminate high-frequency interference and enhance the edges of large-scale sources. In contrast, high-frequency IMFs use small-scale SEs to enhance the edges of small-scale sources.   For example, the proposed method for IMF 1 edge detection shows that the NW trending fault in the central part of the study area (numbered F5 in Fig. 22) is a linear fracture, as shown in Fig. 20(a). However, a series of circles is present in Fig. 20(d), which makes interpretation difficult. This fault shows a low and gentle anomaly in IMF 2 . The SE size used by this IMF is small, which can filter this anomaly. Therefore, the fault cannot be detected in the edge of IMF 2 extracted by MEE, as shown in Fig. 20(b). Correspondingly, in Fig. 20(e), the obvious NW trending fault can be seen after the edge detection of IMF 2 by Tilt. This fault is characterized by high frequency in IMF 3 , so it can be detected by both MEE and Tilt.
Therefore, MEE extracts the main edges of each scale based on 2D-VMD, separates the edges of the different scales, and weakens the frequency aliasing between IMFs by setting different SE scales. Correspondingly, the advantage of the Tilt method is that it can recognize the edge of a weak anomaly better. In addition, MEE has a better edge recognition ability for linear fracture, while Tilt has a stronger ability for the edge recognition of geological bodies. Pilkington and Keating [40] compared several edge detection algorithms and concluded that no single edge detection method can extract all edge features. Therefore, the coordinated use of multiple methods can effectively improve the reliability of edge detection.

E. GEOLOGICAL INTERPRETATION
According to the analysis of the geophysical data from the Hanmiao area, the magnetic bodies are mainly various kinds of intermediate-felsic volcanic rocks and intrusive rocks.  Shale, sandstone and other surrounding rocks are mostly weakly magnetic or non-magnetic. Therefore, the multi-scale edge detection results can be used to delineate the distributions of volcanic rocks and intrusive rocks. Each anomaly has its own scale and characteristics that can be represented by the multi-scale detection algorithm. The method proposed by Blackly and Simpson is used to detect edge positions. The local maximum points detected by each IMF are connected into lines, and the lines with lengths of 2 km or greater are retained, thereby completing the multi-scale edge detection process for a set of magnetic data. In this study, the multiscale edge detection results are directly added to the geological map in Fig. 21. Lines of different colours are used to represent edges: pink, green and blue represent the edges of IMF 1 , IMF 2 and IMF 3 , respectively. By using the similarity among the edges at different scales, we can obtain the penetration depth and strike of the faults. The pink lines indicate faults with a deep penetration depth, while the green and blue lines indicate faults with relatively shallow penetration depths. The closer together these lines are, the larger the dip angles of the faults. Set of lines can reflect the fault characteristics, such as the strike direction, dip direction and cutting depth. Using the multi-scale edge detection results in combination with remote sensing and field geological surveys, a comprehensive geological interpretation of the study area is shown in Fig. 22. F 5 is called the Oerge-Chaganhada fault, which is revealed very clearly in the IMF 1 and IMF 2 edge detection results. The Oerge-Chaganhada fault is a very important metallogenic structure that is reflected in the edge detection results at different scales and may be a deep fault developed in this region. In terms of physical properties, the Paleozoic belt of uplift in the northwestern part of the study area is close to the volcanic basin. Among the rocks in that area, the average magnetic susceptibility of the intrusive granodiorite is 0.013 SI, while that of the Permian Linxi Formation is 0.00015 SI. The edge between the granodiorite and Linxi Formation is very obvious in the edge recognition results from the aeromagnetic data. The magnetic susceptibilities of the volcanic rocks and intrusive rocks in the Wote-Bayanchagan basin are approximately 0.006∼0.012 SI, while those of the Permian sedimentary strata along the periphery of the basin are approximately 0.00007∼0.00027 SI; these susceptibilities differ on the order of 10 2 . By comparing the geological boundaries of the Wote-Bayanchagan basin with the recognized magnetic edges, it can be seen that the ring-shaped edges are the positions of volcanic rocks, which is consistent with the field geological survey.
For the Abuga rock mass, the detected magnetic edges are also in good agreement with the geological boundaries. The lithology of this rock mass is mainly adamellite, the magnetic susceptibilities of which are approximately 0.00027∼0.00435 SI. The lithologies of the surrounding rocks are mainly metamorphic sandstone and metamorphic siltstone; the magnetic susceptibilities of these rocks are approximately 0.00007∼0.00027 SI. The physical properties of these two masses are very different. The Aorege-Gulbalkin basin is dominated by meso-basaltic volcanic rocks, the magnetic susceptibilities of which range from 0.005 SI to 0.0248 SI with an average of 0.01244 SI. These magnetic properties are caused mainly by andesite.
There are also some circular or elliptical edges along the periphery of the Mesozoic volcanic basin in the Hanmiao area. These edges may correspond to deep sources. For example, the Baiyinhuo area has recognized edges covering a large area, but only a small part of a tuff outcrop is exposed at the surface. This indicates that the actual distribution range of volcanic rocks in the Manketou'ebo Formation may be quite large under the Quaternary sediments. We also employ magnetic data to recognize the edges within the Manghatu area, and the results similarly show that the granodiorite in this area has obviously been limonitized, and there is hidden magnetite at depth.
Interpreting the edges of small-scale IMFs is more difficult than the edges of large-scale IMFs. However, through field verification, the edges of large-scale IMFs are mainly related to regional geological bodies and faults, while the edges of small-scale IMFs are mostly related to small geological bodies and shallow linear fracture. The ring and radial faults near the volcanic apparatus are also related to the edges of the small-scale IMFs. Among them, the edge identified by IMF 2 is mainly a NE trending fault, which is most obvious in the Wote area. Examples include the F7, fault verified in the field, and the F6 fault, verified by remote sensing. In addition, the edges identified by IMF 2 are mostly related to the volcanic apparatus. In Gulbalkin and Baiyinhuo, the ringshaped area is the location of the volcanic apparatus. A series of faults near Baiyinhuo, including the verified F15 fault, are also obvious at this scale.
The edge identified by IMF 3 is mainly a series of NW trending faults with dense distribution and shallow depth. The verified NW trending faults such as F1, F5 and F13 have been well recognized at the edges of this scale. The F9 fault that has been verified in the Abuga area is also well identified in this scale edge. In particular, the Oerge-Chaganhada fault (F5) has a high intensity and small range in IMF 3 . This feature is well reflected after edge enhancement with an SE size of 4 × 4. This IMF also identified ring edges in Gulbalkin, Bayanchagan, Baiyinhuo and other volcanic rock areas. A series of radial fractures are distributed around the ring. Some areas are confirmed to be volcanic apparatus after verification, and other unverified areas with similar characteristics are also presumed to be volcanic apparatus.

V. CONCLUSION
In this paper, 2D-VMD is applied to the multi-scale decomposition of potential field data for the first time. This method can effectively decompose potential field data into k IMFs. On this basis, we propose an edge detection algorithm, MEE, which is based on MM. MEE can effectively enhance the edges of IMFs by setting different SE scales.
Using the MEE variance to quantitatively determine the optimal SE scale is helpful to calculate the MEE values more effectively. Through model tests, compared with other potential field separation methods and edge detection methods, 2D-VMD and MEE have superior applicability. When positive and negative abnormal bodies exist in the model at the same time, the proposed method does not generate false extra edges and has a stronger anti-noise capability. Hence, the presented approach is suitable for both gravity and magnetic data. In this paper, the edge detection accuracy of several different methods is compared. The experimental results show that the proposed method can better solve the edge detection problem, with improved accuracy and anti-noise performance.
The proposed method is applied to the data processing and interpretation of 1:50000 RTP aeromagnetic data from the Hanmiao area of Chifeng city, Inner Mongolia, China. The results show that the edges detected by utilizing aeromagnetic data from the middle of the study area are consistent with the geological boundaries of Mesozoic volcanic rocks. The granodiorite in the northwestern part of the study area is also clearly revealed in the edge detection results, and the latest field geological survey shows concealed magnetite at depth. Both these theoretical model tests and this application example show that the novel multi-scale edge detection method based on 2D-VMD and MEE can effectively extract edges from potential field data. MEE can suppress frequency aliasing by changing the SE scale. This technique is a powerful tool for identifying regional tectonics and delineating concealed geological bodies. The proposed method has a good reference value for enriching multi-scale edge detection methods.