## 2 Inversion Methodology  
Generally, the model domain is discretized into many cells, whose physical properties are the model parameters. According to the linear relationship between model parameters and measured data, the inverse problem of potential field data has the matrix form as  
\[\mathbf{d_{obs}}=\mathbf{G}\mathbf{m}, \tag{1}\]  
where \(\mathbf{d_{obs}}\in R^{n}\) represents the measured data vector, \(\mathbf{m}\in R^{m}\) is the unknown parameters, and \(\mathbf{G}\in R^{n\times m}\) represents the kernel matrix. \(n\) represents the number of measured data, and \(m\) represents the number of the unknown parameters. In general, \(m\) is much larger than \(n\). Therefore, the inversion of potential field data belongs to an under-determined problem. Through the singular value decomposition (SVD) of the matrix \(\mathbf{G}\) or the least squares solution, Equation (1) can be used to calculate the inversion result directly. However, in this way, the inversion result may not conform to reality. These ill-posed problems can be solved by introducing a regularization term [4]. Then, the objective function of this inversion has the form  
\[\phi=\|\mathbf{W_{d}}(\mathbf{G}\mathbf{m}-\mathbf{d_{obs}})\|_{2}^{2}+\mu^{2} \|\mathbf{W_{m}}\mathbf{Z}\mathbf{m}\|_{2}^{2}, \tag{2}\]  
In Equation (2), \(\|\mathbf{W_{d}}(\mathbf{G}\mathbf{m}-\mathbf{d_{obs}})\|_{2}^{2}\) is the data misfit, \(\|\mathbf{W_{m}}\mathbf{Z}\mathbf{m}\|_{2}^{2}\) is the regularization term, and \(\mu\) is the regularization parameter, which is used to balance these two terms. In the data misfit, \(\mathbf{W_{d}}\) represents a data weighting matrix. In the regularization term, \(\mathbf{W_{m}}\in R^{4m\times m}\) represents a smooth constraint matrix and its matrix form is \(\mathbf{W_{m}}=[\mathbf{W_{s}};\,\mathbf{W_{s}};\,\mathbf{W_{y}};\,\mathbf{W_{ z}}]\), where \(\mathbf{W_{s}}\in R^{m\times m}\), \(\mathbf{W_{s}}\in R^{m\times m}\), \(\mathbf{W_{y}}\in R^{m\times m}\), and \(\mathbf{W_{z}}\in R^{m\times m}\) are different component matrices, respectively [11; 12]; and \(\mathbf{Z}\) represents a depth-weighting matrix [11; 12]. Generally, the physical bound constraint is incorporated into this inversion for obtaining a geologically plausible inversion result. Then, the inverse problem of potential field data becomes the following constrained minimization problem  
\[ \begin{array}{l}\text{minimize}:\phi=\|\mathbf{W_{d}}(\mathbf{G}\mathbf{m} -\mathbf{d_{obs}})\|_{2}^{2}+\mu^{2}\|\mathbf{W_{m}}\mathbf{Z}\mathbf{m}\|_{2 }^{2}\\ \text{subjectto}:\mathbf{m_{min}}\leq\mathbf{m}\leq\mathbf{m_{max}}\end{array} , \tag{3}\]  
where \(\mathbf{m_{min}}\) and \(\mathbf{m_{max}}\) are vectors consisting of the lower and upper physical bounds on the unknown model values. Because the matrix \(\mathbf{Z}\) is a diagonal matrix, the following transformations can be performed easily:  
\[\mathbf{h}=\mathbf{W_{d}}\mathbf{G}\mathbf{Z}^{-1},\mathbf{r}=\mathbf{W_{d}} \mathbf{d_{obs}},\text{ and }\mathbf{y}=\mathbf{Z}\mathbf{m}. \tag{4}\]  
Then, Equation (2) is rewritten as  
\[\phi=\|\mathbf{h}\mathbf{y}-\mathbf{r}\|_{2}^{2}+\mu^{2}\|\mathbf{W_{m}} \mathbf{y}\|_{2}^{2}. \tag{5}\]  
The constrained minimization problem (3) becomes the following new problem  
\[ \begin{array}{l}\text{minimize}:\ \phi=\|\mathbf{h}\mathbf{y}-\mathbf{r}\|_{2}^{ 2}+\mu^{2}\|\mathbf{W_{m}}\mathbf{y}\|_{2}^{2}\\ \text{subjectto}:\ \mathbf{y_{min}}\leq\mathbf{y}\leq\mathbf{y_{max}}\end{array} , \tag{6}\]  
where \(\mathbf{y_{min}}=\mathbf{Z}\mathbf{m_{min}}\) and \(\mathbf{y_{max}}=\mathbf{Z}\mathbf{m_{max}}\) are vectors of the lower and upper bounds for \(\mathbf{y}\). The logarithmic barrier method is adopted to incorporate the physical bound constraint into the objective function [13; 14; 15].