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).