## 2 Theory and Methodology  
### Apparent Resistivity Calculated by the Magnetic Field Gradient  
The schematic diagram of GAFEM detection system is shown in Figure 1; the electric source is located in the origin along the \(x\)-axis, and the \(z\)-axis is positive downward, with the origin at the center of the electric source. The received magnetic sensor is in the air, lifted by an unmanned aerial vehicle, \((x_{r},y_{r},z_{r})\) is the coordinate of the receiver point. The FEM response of the magnetic field response \(H_{z}\) of the receiver point in a uniform half-space can be expressed as [2]  
\[ \begin{array}{l}H_{z}=\frac{IL}{2\pi t}\cdot\frac{v_{r}}{R}\int_{0}^{\infty} \frac{\lambda^{2}}{\lambda+\mu_{1}}e^{\lambda z_{r}}J_{1}(\lambda R)d\lambda\\ \phantom{H_{z}=}=\frac{3i\lambda\mu_{1}\mu_{2}}{2\pi\omega rR^{2}}[1-e^{-ik_{1 }R}(1+ik_{1}R-\frac{1}{3}k_{1}^{2}R^{2})],\end{array}  \tag{1}\]  
where \(R=(x_{r}{}^{2}+y_{r}{}^{2}+z_{r}{}^{2})^{1/2}\), \(J_{1}\) is the first-order Bessel function, \(u_{1}=(\lambda^{2}-k_{1}{}^{2})^{1/2}\), \(k_{1}{}^{2}=\omega^{2}\mu\varepsilon-i\omega\mu\sigma,f\) is the transmitting frequency, the angular frequency \(\omega=2\pi f\) is sufficiently low that displacement currents are negligible, i.e., \(k_{1}{}^{2}=-i\omega\mu\sigma\), \(\lambda\) is the [PERSON] integral variable, \(\rho\) is the resistivity of the earth, \(\sigma=1/\rho\), \(\mu\), and \(\varepsilon\) are the conductivity, magnetic permeability, and electric permittivity of the earth, respectively. \(IL\) is the transmitting electric moment, \(i\) is the imaginary unit. Figure 1: Schematic diagram of ground-airborne electromagnetic detection. According to [PERSON] et al. [9], the space magnetic gradient anomaly can recognize the underground electrical anomalies effectively, the spatial gradient of \(H_{z}\) at the receiver point can be calculated by:  
\[H_{z}^{\Delta x}(x_{i},y_{j},z_{k},f) =\frac{\big{|}H_{z}(x_{i},y_{j},z_{k},f)\big{|}-\big{|}H_{z}(x_{i-1 },y_{j},z_{k},f)\big{|}}{\Delta L_{x}}, \tag{2}\] \[H_{z}^{\Delta y}(x_{i},y_{j},z_{k},f) =\frac{\big{|}H_{z}(x_{i},y_{j},z_{k},f)\big{|}-\big{|}H_{z}(x_{i}, y_{j-1},z_{k},f)\big{|}}{\Delta L_{y}}. \tag{3}\]  
where \(\Delta L_{x}=\left\lvert x_{i}-x_{i-1}\right\rvert\) and \(\Delta L_{y}=\left\lvert y_{j}-y_{j-1}\right\rvert\) represents the distance between the two adjacent points in the \(x\) direction and \(y\) direction, respectively. Subscripts \(i,j,k\) represent the number of the receiver point coordinates in \(x,y,z\) direction, respectively. Equations (1)-(3) reveal that the primary influencing factors of the magnetic field gradient include the coordinates of the receiver point, transmitting frequency, and subsurface electric parameter \(\rho\). To investigate the relationship between the magnetic field gradient and the resistivity, we calculate the magnetic field gradient for a changing resistivity of the homogeneous earth as an example, the receiver point and transmitting frequency are set as (3000, 5000, \(-\)30) m, \(f=64\) Hz, as shown in Figure 2.  
As shown in Figure 2, as the resistivity increases, both the magnetic gradients in the \(x\) and \(y\) directions exhibit a consistent upward trend, i.e., at a certain observation point, there exists a one-to-one mapping relationship between the magnetic field gradient and resistivity. Due to the implicit relationship between the magnetic field gradient and resistivity, Equations (2) and (3) are equations with unknown resistivity, not a solution to the equation. Therefore, the optimal value of apparent resistivity can be obtained by solving Equations (2) and (3) using the iteration-matching algorithm. For a known value \(H_{z}^{\Delta x}\) or \(H_{z}^{\Delta y}\) at a given survey point and a given frequency, an initial range of resistivity value is provided, and the resistivity \(\rho\) is iteratively calculated using the bisection method until the error between the calculated and known magnetic field gradients is less than a predefined threshold (e.g., 1%).