## 2 Mixed L\({}_{\rm p}\) Norm Regularization for 3-D Magnetic Inversion  
### Methodology  
The mixed L\({}_{\rm p}\) norm regularization algorithm employed in this paper was proposed by [PERSON] and [PERSON] (2019) [30]. It replaces the regularization term in the traditional Tikhonov regularization (L\({}_{2}\) norm) objective function with the norms \(p_{x,y,z,s}\in\{0,\,1,\,2\}\). These norms control the regularization term and play a crucial role in determining the final model's features. Typically, the objective function \(\phi(m)\) of the inverse problem can be expressed by the following formulation:  
\[\underset{m}{min}\phi(m)=\phi_{d}+\beta\phi_{m} \tag{1}\]  
where \(\phi_{d}\) denotes the data fitting term, \(\phi_{m}\) represents the regularization term, and \(\beta\) is the regularization parameter that balances the two. Assuming that the noise contained in the observed data follows an independent Gaussian distribution, the data fitting term \(\phi_{d}\) is typically expressed in the following mathematical form:  
\[\phi_{d}=\sum_{i=1}^{N}\left(\frac{d_{i}^{pred}-d_{i}^{obs}}{\sigma_{i}}\right) ^{2} \tag{2}\]  
where \(d_{i}^{obs}\) represents observed data, \(d_{i}^{pred}\) represents calculated magnetic anomalies caused by the inversion model, and \(\sigma_{i}\) represents their uncertainties. Among these, \(\phi_{m}\) is the regularization term of the objective function, which plays a critical role in determining the characteristics of the final model. In the mixed L\({}_{\rm p}\) norm inversion, \(p_{x,y,z,s}\) can take the independent values of 0, 1, or 2. The regularization term is as follows:  
\[\phi_{m}=\sum_{i=s,x,y,z}\alpha_{i}\int_{V}\left|f_{j}(m)\right|^{p_{i}}dV \tag{3}\]  
where:  
\[f_{s}=m,f_{x}=\frac{\partial m}{\partial x},f_{y}=\frac{\partial m}{\partial y },f_{z}=\frac{\partial m}{\partial z} \tag{4}\]  
The discrete form of the formula is expressed as follows:  
\[\phi_{m}=\alpha_{s}\left\|V_{s}m\right\|_{p_{s}}^{p_{s}}+\sum_{j=x,y,z}\alpha _{j}\left\|V_{j}G_{j}m\right\|_{p_{j}}^{p_{j}} \tag{5}\]  
where \(V_{s}\), \(V_{x}\), \(V_{y}\) and \(V_{z}\) are diagonal matrices corresponding to the \"volumes\" of the cells used in discretizing the problem, and \(G_{x}\), \(G_{y}\), and \(G_{z}\) are the discrete gradient operators. Byemploying various mixed norms for model regularization, the mixed L\({}_{\text{p}}\) norm inversion method yields a range of solutions. This approach highlights the inherent multiplicity of solutions in the inversion process and mitigates the risk of over-interpretation that may arise from considering only a single solution. The data experiments with the method mentioned above indicate that when sparse norms (L\({}_{0}\), L\({}_{1}\)) are applied to the model vector \(\mathbf{m}\), an issue of dependence on the upper limit of physical properties arises. Conversely, when sparse norms are applied to the model gradient \(\
abla\mathbf{m}\), a \"staircase\" effect is observed. The former leads to the physical-property values of the anomaly body in the inversion approaching the upper limit, while the latter causes the anomaly body to exhibit horizontal or vertical boundaries. To address these issues, [PERSON] and [PERSON] (2021) [28] further refined this method based on the work of [PERSON] and [PERSON] (2019) [30]. The regularization term in the objective function of the mixed L\({}_{\text{p}}\) norm inversion consists of four components, as follows:  
\[\phi_{m}=\alpha_{s}\int|f_{s}(m)|^{p}\,dv+\sum_{j=x,y,z}\alpha_{j}\int|f_{j}(m )|^{q}\,dv \tag{6}\]  
The discrete form of the formula is as follows:  
\[\phi_{m}=\alpha_{s}\|W_{s}m\|_{p}^{p}+\sum_{j=x,y,z}\alpha_{j}\Big{\|}W_{j} \
abla_{j}m_{j}\Big{\|}_{q}^{q} \tag{7}\]  
where \(f\) is defined as in Equation (4). The parameters \(\alpha_{s}\) and \(\alpha_{j}\) are weighting coefficients that control sparsity and smoothness, while \(p\) and \(q\) are the norm exponents that regulate the smoothness of the model's small-scale features and gradients, respectively.