2.2 Inverse Problem
Source inversion is an ill-posed, non-unique problem that can be solved by incorporating model regularization into the inverse problem. According to the SP inversion method proposed by Ahmed and Jardani (2013), the spatial distribution of potential anomalies caused by the underground current density vector is obtained by M:
\(\psi\left(N\right)=\int_{\Omega}{G\left(N,M\right)}\mathbf{j}_{s}\left(M\right)\text{dV}\text{\ \ \ }\)(2.8)
Where \(\mathbf{j}_{s}(M)\) is the current density source at a set of points M, \(\psi(N)\) is the potential at a set of \(N\) electrodes, and\(G\left(N,M\right)\) represents the measured natural potential data at point \(N\) and is currently at point \(M\) Kernel matrix. The calculation of the \(G\left(N,M\right)\) kernel function depends on the selection of two parameters: 1) the conductivity value of the medium; 2) the number of discrete elements \(M\). The kernel matrix can be solved by using the finite element method (Trujillo-Barreto et al., 2004).
\(\phi_{d}=\left\|W_{d}(G\omega-\psi^{\text{obs}})\right\|_{2}^{2}\)(2.9)
\(\psi^{\text{obs}}=[\psi_{1}^{\text{obs}},\psi_{2}^{\text{obs}},\psi_{3}^{\text{obs}},\ldots{,\psi}_{N}^{\text{obs}}]^{T}\)(2.10)
\(W_{d}=diag[\frac{1}{\delta_{1}},\frac{1}{\delta_{2}},\ldots,\frac{1}{\delta_{N}}]\)(2.11)
Where \(\left\|\right\|_{2}\) is the L2 norm, \(G\) represents the kernel matrix of N × 2M (corresponding to the self-potential measured at each site), and \(\omega\) is the current density \(j_{s}\) related to the potential and the kernel matrix \(G\).\(\psi^{\text{obs}}\) denotes a vector of N elements corresponding to a set of self-potential data, \(W_{d}\) is an N × N diagonal weighted square matrix, and \(\delta_{1}\) denotes the square of the deviation (Linde et al., 2007). Because the inversion leads to ill-posedness of the solution, a regularization term, and a depth weighting function \(\phi_{m}\) are added to the analysis process to reduce the problem of overfitting. Establish the global objective function term \(\phi_{T}\) (Menke, 1984):
\(\phi_{T}=\phi_{d}+\phi_{m}=\left\|W_{d}(G\omega-\psi^{\text{obs}})\right\|_{2}^{2}+\lambda^{2}\left\|W_{m}(\omega-\omega_{0})\right\|_{2}^{2}\)(2.12)
\(W_{m}=LW_{z}={(z_{0}+z)}^{-\frac{\beta}{2}}\ \)denotes a depth-weighted matrix of N × 2M (Li & Oldenburg, 1998),\(z_{0}\ \)denotes the observation height of the model unit, and \(L\)is the smoothness based on the first derivative of \(\omega\), \(\beta\)is a constant term between 1 and 4. λ denotes the regularized constraint term 0 <λ <∞, and the objective function of the above formula is in the standard form:
\({\phi_{T}=\left\|\overset{\overline{}}{G}\overset{\overline{}}{\omega}-{\overset{\overline{}}{\psi}}^{\text{obs}}\right\|}_{2}^{2}+\lambda\left\|\overset{\overline{}}{\omega}-{\overset{\overline{}}{\omega}}_{0}\right\|_{2}^{2}\)(2.13)
Where\(\overset{\overline{}}{G}\)\(\overset{\overline{}}{\omega}\)\({\overset{\overline{}}{\psi}}^{\text{obs}}\)\(\text{and\ }L_{P}\)are matrices deduced by the Elédn’s algorithm based on QR decomposition (Elédn, 1977). The solution to finding the minimum value of the objective function \(\phi_{T}\) (Hansen, 1998):
\({min(\left\|\overset{\overline{}}{G}\overset{\overline{}}{\omega}-{\overset{\overline{}}{\psi}}^{\text{obs}}\right\|}_{2}^{2}+{\lambda\left\|\overset{\overline{}}{\omega}-{\overset{\overline{}}{\omega}}_{0}\right\|}_{2}^{2})\)(2.14)
\(\hat{\omega}\left(\lambda\right)=[{\overset{\overline{}}{G}}^{T}\overset{\overline{}}{G}+\lambda E]^{-1}({\overset{\overline{}}{G}}^{T}{\overset{\overline{}}{\psi}}^{\text{obs}}+\lambda{\overset{\overline{}}{\omega}}_{0})\)(2.15)
When the permeability is known, the current density distribution can be calculated for the prior model of the objective function, and in the case where there is no prior information, we assume the model\({\overset{\overline{}}{\omega}}_{0}=0\):
\(\hat{\omega}\left(\lambda\right)=[{\overset{\overline{}}{G}}^{T}\overset{\overline{}}{G}+\lambda E]^{-1}{\overset{\overline{}}{G}}^{T}{\overset{\overline{}}{\psi}}^{\text{obs}}\)(2.16)
The singular value decomposition (SVD) method is used\(\overset{\overline{}}{G}=\sum_{i=1}{u\Lambda v^{T},=min(N,M)}\)and redefine \(\hat{\omega}\left(\lambda\right)\) such that:
\(\hat{\omega}\left(\lambda\right)=\frac{\Lambda_{i}}{\Lambda_{i}^{2}+\lambda}u_{i}^{T}v_{i}{\overset{\overline{}}{\psi}}^{\text{obs}}\)(2.17)
Where \(\Lambda\) is a singular diagonal matrix, \(\Lambda_{i}\) denotes a singular value component on the diagonal of the matrix and \(uv^{T}\)is a singular vector. In terms of calculation efficiency, the equation after SVD can better reflect the information of the main components, which effectively reduces the amount of calculation. Besides, the selection of the regularization parameter λ is also crucial. This paper uses the L-curve method (logarithmic-logarithmic intersection graph of\(\phi_{d}\) and \(\phi_{m}\)) to define the optimal value of the regularization parameter (Hansen, 1998).