## 2 Methods
### Forward and Inversion Theory
#### 2.2.2 Inversion Theory  
The linear inversion equation for the known resistivity is  
\[\mathbf{A}\Delta\mathbf{m}=\Delta\mathbf{d}, \tag{9}\]  
where \(\Delta\mathbf{d}\) is the data residual vector, which is the difference between the measured apparent resistivity \(\rho_{a}\) and the simulated apparent resistivity \(\rho_{s}\) in logarithm (\(\Delta d_{i}=ln\rho_{ai}-ln\rho_{si}\), \(i=1,2,\cdots,N\)); \(N\) is the number of measured apparent resistivity data points; \(\Delta\mathbf{m}\) is the correction vector of the model parameters (\(\Delta m_{j}=\Delta ln\rho_{j},j=1,2,\cdots M\)); \(\rho_{j}\) is the resistivity of j-th grid node; \(M\) is the number of the model parameters; and \(\mathbf{A}\) is the partial derivative matrix (\(a_{ij}=\partial ln\rho_{si}/\partial ln\rho_{j}\)). In 3D resistivity inversion, Equation (9) is typically underdetermined and ill-posed. To ensure the stability of the inversion process, prior constraints [44] (including smoothness and background constraints) need to be applied to the geoelectric model. By integrating Equation (9) with these prior constraints, the objective function for resistivity inversion can be constructed in the least-squares sense [45]:  
\[\Phi=\|\Delta\mathbf{d}-\mathbf{A}\Delta\mathbf{m}\|_{2}^{2}+\lambda_{c}\| \mathbf{C}(\mathbf{m}+\Delta\mathbf{m})\|_{2}^{2}+\lambda_{b}\|\mathbf{I}( \mathbf{m}+\Delta\mathbf{m}-\mathbf{m}_{b})\|_{2}^{2}, \tag{10}\]  
where \(\mathbf{m}\) is the model parameter vector; \(\mathbf{m}_{b}\) is the background model representing prior information about subsurface properties, typically derived from geological, geophysical, or historical survey data; \(\mathbf{C}\) is the smoothness matrix [46]; \(\mathbf{I}\) is the identity matrix; and \(\lambda_{c}\) and \(\lambda_{b}\) are regularization factors. In this equation, the first term on the right represents the conventional least-squares data misfit term, quantifying the discrepancy between the observed data and the forward simulated data. The second term imposes smoothness on the inverted resistivity model to reduce unrealistic spatial variations. The third term incorporates prior information through the background model constraint, guiding the solution toward geological plausibility. By differentiating both sides of Equation (10) with respect to \(\Delta\mathbf{m}^{\mathrm{T}}\) and setting it equal to zero, the linear system of equations for resistivity inversion can be obtained after rearranging and transformation:  
\[\left(\mathbf{A}^{\mathrm{T}}\mathbf{A}+\lambda_{c}\mathbf{C}^{\mathrm{T}} \mathbf{C}+\lambda_{b}\mathbf{I}\right)\Delta\mathbf{m}=\mathbf{A}^{\mathrm{ T}}\Delta\mathbf{d}-\lambda_{c}\mathbf{C}^{\mathrm{T}}\mathbf{C}\mathbf{m}+ \lambda_{b}\mathbf{I}(\mathbf{m}_{b}-\mathbf{m}), \tag{11}\]\[\text{If }\mathbf{G}= \begin{vmatrix}\mathbf{A}\\ \sqrt{\lambda_{c}}\mathbf{C}\\ \sqrt{\lambda_{b}}\mathbf{I}\end{vmatrix} ,\mathbf{b}= \begin{vmatrix}\Delta \mathbf{d}\\ -\sqrt{\lambda_{c}}\mathbf{C}\mathbf{m}\\ \sqrt{\lambda_{b}}(\mathbf{m}_{b}-\mathbf{m})\end{vmatrix} ,\text{ then Equation (11) can be simplified to:}\\ \\ \mathbf{G}^{\text{T}}\mathbf{G}\Delta\mathbf{m}=\mathbf{G}^{\text{T}}\mathbf{b }\text{ or }\mathbf{G}\Delta\mathbf{m}=\mathbf{b},\end{vmatrix} \tag{12}\]  
The linear inversion Equation (12) is characterized by being full-rank, ill-posed, and large in scale. The conjugate gradient method has significant advantages in solving such equations and is widely applied in scientific research and engineering computations [47; 48; 49]. In this study, the least-squares conjugate gradient method is used to solve Equation (12) [50]. First, let the initial solution of Equation (12) be \(\Delta\mathbf{m}_{0}=0\); then, the initial residual is, with the initial values of the gradient vector \(\mathbf{g}\) and the conjugate direction vector \(\mathbf{p}\), set to \(\mathbf{g}_{0}=\mathbf{p}_{0}=\mathbf{G}^{\text{T}}\mathbf{r}_{0}=\mathbf{G}^ {\text{T}}\mathbf{b}\).