## 2 Method  
### Mathematical Method  
#### 2.1.1 Fractal Surface Variance Statistics Method  
The calculation method of the fractal dimension was first proposed for topographic surfaces [24]. In this study, the method is improved, modified and introduced in the field of gravity. The gravity anomaly (abbreviated to \(\Delta g\)) of the study area is treated as the self-similar topographic surface. This property is inherited from the fractal nature of geology itself. Therefore, the gravity anomaly \(\Delta g(i,j)\) is defined as a function with distance, and its expected variance between the point \((i,j)\) and the other point \((p,q)\) is presented by Equation (1):  
\[E\left\{\left[\Delta g(i,j)-\Delta g(p,q)\right]^{2}\right\}=kr^{2H}, \tag{1}\]  
where \(r\) is the horizontal distance between points \((i,j)\) and \((p,q)\), and \(H\) equals 3-D. \(D\) is the fractal dimension, and \(k\) is the constant related with the square deviation. This is a power law model for the variogram in two dimensions. The natural logarithm is taken on both sides of Equation (1), then  
\[\ln E\Big{\{}[\Delta g(i,j)-\Delta g(p,q)]^{2}\Big{\}}=\ln k+2H\ln r \tag{2}\]  
Equation (2) is further simplified to  
\[\ln E\Big{\{}[\Delta g(i,j)-\Delta g(p,q)]^{2}\Big{\}}=\ln k+(3-D)\ln\Big{\{} \left[\Delta x(i-p)\right]^{2}-\left[\Delta y(j-q)\right]^{2}\Big{\}}, \tag{3}\]  
where \(\Delta x\) and \(\Delta y\) are the grid interval of the \(x\)-axis and \(y\)-axis, respectively. Based on Equation (2) and the log-log coordinate, there is a linear region in the unscaled range when fractal analysis is performed. This property illustrates that the linear region for a gravity anomaly contains self-similarity and develops a 1-, 2- and 3-dimensional predictive model of correlated random data. Therefore, the gradient \(b\) of the linear region could be obtained by the calculation of the fractal dimension, such that we have  
\[D=3-b/2\ \text{and}\ 2\ <\ D<\ 3, \tag{4}\]  
The gradient of the fractal surface variance statistics filter is a vector quantity. 1.2 Conventional Method of the Directional Gradient Filter and the Second Order Derivate Filters of Laplacian 1 and 2  
The gradient of function \(\Delta g(x,y)\) is the vector, and its module is defined as  
\[G[\Delta g(x,y)]=\sqrt{\left(\frac{\partial\Delta g}{\partial y}\right)^{2}+ \left(\frac{\partial\Delta g}{\partial x}\right)^{2}}, \tag{5}\]  
where \(G[\Delta g(x,y)]\) is the module of the gradient, and it is called the gradient in this study. The sliding window of the \(3\times 3\) matrix is characterized with directions, which include east, southeast, south and southwest (Figure 1). They are the first-order derivatives. \(xx\) and \(yy\) second-order derivatives of Laplacian 1 and 2 are introduced to discern the whole boundaries of the gravity model and to make the comparison. The eastwest gradient filter follows Equation (6) and Figure 1b:  
\[ \begin{array}{l}G_{east}[\Delta g(x,y)]\ =|-\Delta g(i-1,j-1)+\Delta g(i-1,j)+ \Delta g(i-1,j+1)\\ -\Delta g(i,j-1)-2\Delta g(i,j)+\Delta g(i,j+1)+\Delta g(i+1,j-1)+\Delta g(i+1, j)-\Delta g(i+1,j+1)|\end{array}  \tag{6}\]  
The southeast gradient filter follows Equation (7) and Figure 1c:  
\[ \begin{array}{l}G_{south}[\Delta g(x,y)]\ =\ |-\Delta g(i-1,j-1)-\Delta g(i-1,j)- \Delta g(i-1,j+1)\\ -\Delta g(i,j-1)-2\Delta g(i,j)+\Delta g(i,j+1)+\Delta g(i+1,j-1)+\Delta g(i+1, j)+\Delta g(i+1,j+1)|\end{array}  \tag{7}\]  
The southnorth gradient filter follows Equation (8) and Figure 1d:  
\[ \begin{array}{l}G_{south}[\Delta g(x,y)]=|-\Delta g(i-1,j-1)-\Delta g(i-1, j)-\Delta g(i-1,j+1)\\ +\Delta g(i,j-1)-2\Delta g(i,j)+\Delta g(i,j+1)+\Delta g(i+1,j-1)+\Delta g(i+1, j)+\Delta g(i+1,j+1)|\end{array}  \tag{8}\]  
The southwest gradient filter follows Equation (9) and Figure 1e:  
\[ \begin{array}{l}G_{southness}[\Delta g(x,y)]\ =\ |\Delta g(i-1,j-1)-\Delta g(i-1, j)-\Delta g(i-1,j+1)\\ \Delta g(i,j-1)-2\Delta g(i,j)-\Delta g(i,j+1)+\Delta g(i+1,j-1)+\Delta g(i+1, j)+\Delta g(i+1,j+1)|\end{array}