# A Stable and High-Precision Downward Continuation Method of Magnetic Data
## 2 Methodology of the Proposed Method
 \tag{5}\]  
where \(\zeta(r)\) is the depth weighted operator, \(k(r)\) is the unsolved model parameter; \(k_{0}\) is the reference model; and \(\alpha_{s}\), \(\alpha_{x}\), \(\alpha_{y}\), and \(\alpha_{z}\) are weighted coefficients, respectively. The \(\alpha_{s}\ll 1\), \(\alpha_{z}=\alpha_{x}=\alpha_{y}=1\) empirically. The matrix form for model objective equation is  
\[ \begin{split}\Phi_{k}=\mathbf{k}\mathbf{W}_{k}^{T}\mathbf{W}_{k} \mathbf{k}\\ \mathbf{W}_{k}^{T}\mathbf{W}_{k}=\alpha_{s}(\mathbf{w}_{s}^{T} \mathbf{w}_{s})+\alpha_{z}(\mathbf{w}_{z}^{T}\mathbf{w}_{z})+\alpha_{x}( \mathbf{w}_{x}^{T}\mathbf{w}_{x})+\alpha_{y}(\mathbf{w}_{y}^{T}\mathbf{w}_{y}) \end{split}  \tag{6}\]  
Naturally, the matrix form for the data misfit equation between model forward data and observed data is  
\[\Phi_{d}=(\mathbf{G}\mathbf{k}-\mathbf{d}_{obs})^{T}(\mathbf{G}\mathbf{k}- \mathbf{d}_{obs}) \tag{7}\]  
The model susceptibility parameters could be calculated by Equation (2); however, the results obtained by directly calculating this equation may not conform to reality. In case of that, we introduce the Tikhonov regularization to solve this problem [21]. The objective function to be solved is formed as Equation (6):  
\[ \begin{split}&\text{minimize}:\Phi=\|\mathbf{G}\mathbf{k}-\mathbf{d}_{ obs}\|_{2}^{2}+\lambda^{2}\|\mathbf{W}_{k}\mathbf{Z}\mathbf{k}\|_{2}^{2}\\ &\text{subject to}:\mathbf{k}_{min}\leq\mathbf{k}\leq\mathbf{k}_{ max}\end{split}  \tag{8}\]  
where \(\|\mathbf{G}\mathbf{k}-\mathbf{d}_{obs}\|_{2}^{2}\) is the data misfit, \(\lambda^{2}\|\mathbf{W}_{k}\mathbf{Z}\mathbf{k}\|_{2}^{2}\) is the regularization term, and \(\lambda\) is the regularization parameter, which is applied to prevent overfitting. \(\mathbf{W}_{k}\in\mathbf{R}^{4m\times m}\) and \(\mathbf{W}_{k}=[\mathbf{W}_{s};\mathbf{W}_{x};\mathbf{W}_{y};\mathbf{W}_{z}]\) represent the smooth constraint matrix, \(\mathbf{Z}\) represents a depth-weighting matrix, and \(\mathbf{k}_{min}\) and \(\mathbf{k}_{max}\) are vectors of the lower and upper physical bounds, respectively, of the unknown model values [21]. Finding an optimal \(\lambda\) is particularly important; in this paper, the unbiased predictive risk estimator method (UPRE) is used to calculate \(\lambda\)[22]:  
\[ \begin{split} UPRE(\lambda)=&\|\mathbf{d}_{obs}- \mathbf{G}(\mathbf{G}^{T}\mathbf{G}+\lambda^{2}\mathbf{W}_{k}^{T}\mathbf{W}_{ k})^{-1}\mathbf{G}^{T}\mathbf{d}_{obs}\|_{2}^{2}\\ +& 2 trace\big{[}\mathbf{G}(\mathbf{G}^{T}\mathbf{G}+ \lambda^{2}\mathbf{W}_{k}^{T}\mathbf{W}_{k})^{-1}\mathbf{G}^{T}\big{]}- \mathbf{N}\end{split}  \tag{10}\]  
The conjugate gradient method (CG) is adopted to solve \(\mathbf{G}(\mathbf{G}^{T}\mathbf{G}+\lambda^{2}\mathbf{W}_{k}^{T}\mathbf{W}_{ k})^{-1}\mathbf{G}^{T}\mathbf{d}_{obs}\). In order to find the trace of the matrix in the denominator, the random estimation method is applied as Equation (11)[23]:  
\[trace(\mathbf{A})=\mathbf{u}^{T}\mathbf{A}\mathbf{u} \tag{11}\]  
where \(\mathbf{u}\) is the random vector with 50% probabilities of \(-1\) and \(1\), respectively. The Equation (10) transform to Equation (12):  
\[ \begin{split} UPRE(\lambda)&\approx\|\mathbf{d}_{ obs}-\mathbf{G}(\mathbf{G}^{T}\mathbf{G}+\lambda^{2}\mathbf{W}_{k}^{T}\mathbf{W}_{ k})^{-1}\mathbf{G}^{T}\mathbf{d}_{obs}\|_{2}^{2}\\ +& 2\mathbf{u}^{T}\mathbf{G}(\mathbf{G}^{T}\mathbf{G} \mathbf{+}\lambda^{2}\mathbf{W}_{k}^{T}\mathbf{W}_{k})^{-1}\mathbf{G}^{T} \mathbf{u}-\mathbf{N}\end{split}  \tag{12}\]  
The value range of \(\lambda\) in this paper is \(\ln^{1}\sim\ln^{12}\), which are put into Equation (12) to find \(\lambda\), which minimizes the function result, which is the optimal regularization parameter. After solving the regulation problem, we focus on how to make the iterate process fast and to obtain a reasonable result directly.