The underground space was assumed to be divided into N M N M = m × n × p prism grid cells, where m and n represent the number of grid divisions in the north and east, respectively, and p represents the number of grid divisions in the vertical direction. The observed magnetic field on the surface consisted of N Δ T N Δ T = m × n observation points, each corresponding to a horizontal position on the prism cell (Figure 1).

According to the theory of magnetic field forward modeling, the forward process can be represented by a filtering formula (Equation 1) that balances the spatial domain accuracy and frequency domain speed in the forward calculation (Jing et al., 2019). Δ T = VM = F − 1 ∑ k = 1 p F V ∼ k · F M ∼ k , 1 ≤ k ≤ p , (1) where vector Δ T represents the forward modeling anomaly of the model with a dimension of N Δ T , M is the magnetization intensity according to the forward formula, also with a dimension of N M , V represents the forward operator that connects observation data Δ T and physical parameters M , known as the forward coefficient matrix, which has a dimension of N Δ T × N M , F represents the fast Fourier transform in the positive direction while F − 1 represents the fast Fourier transform in the inverse direction, V ∼ k ( 2 m − 1 × 2 n − 1 )denotes the prism forward kernel matrix at depth k, and M ∼ k ( 2 m − 1 × 2 n − 1 ) refers to an extended matrix obtained by adding zeros to the magnetization intensity matrix at depth k. The symbol “∙” denotes the Hadamard product, representing element-wise multiplication between matrices of the same order. The derivation of Equation 1 is detailed in Supplementary Appendix A.

The first m × n elements in matrix M ∼ k represent the density matrix elements of the kth layer in the model while the remaining elements are set to zero. The first m × n elements in matrix V ∼ k correspond to the magnetic field response values at surface observation points for prism M ∼ m , n , k with index number m , n in the density matrix of the kth layer (with default starting index as 1). At this time, prism M ∼ m , n , k has a density set to the unit magnetization intensity. The remaining elements in matrix V ∼ k are stored symmetrically with respect to index number M ∼ m , n , k .

To calculate the magnetic field at any point on the undulating surface, we adopted a systematic approach calculates the magnetic field distribution on multiple horizontal planes by using a fast algorithm and deriving observation point through a local interpolation technique. There are many local interpolation methods, such as local upward continuation, inverse distance-weighted interpolation, trilinear interpolation, Kriging interpolation, and cubic spline. If only the accuracy of forward modeling is considered, the nonlinear interpolation method may have a higher calculation accuracy. However, when forward calculation is used for linear inversion calculations, the interpolation method should be linear interpolation. In this study, a composite interpolation method that combines local continuation and trilinear interpolation was adopted. The specific steps were as follows. First, multiple horizontal planes (Figure 2) were selected as the calculation basis according to the terrain characteristics of the undulating observation surface. Multiple horizontal observation surfaces were arranged at equal intervals according to the longitudinal length of the underground space dissection grid. Such a setting enables the reuse of a large number of forward coefficient matrices, thereby saving considerable calculation time and storage space (Li et al., 2018).

An efficient three-dimensional magnetic field fast forward algorithm (Equation 1) was then used to perform the magnetic field forward calculation to obtain accurate magnetic field distribution data on these planes. The magnetic field calculation formula (Equation 2) for multiple horizontal planes is as follows: Δ T s = V s M (2) where vector Δ T s represents the magnetic fields on the multiple horizontal observation planes and its dimension is N Δ T s × 1 , i.e., the total number of observation points on multiple horizontal observation planes. Furthermore, V s represents the forward coefficient matrix for calculating the magnetic field values of multiple horizontal observation planes, whose dimension is N Δ T s × N M .

Finally, a high-precision local interpolation algorithm was adopted. The idea of the local interpolation algorithm (Equation 3) is as follows. First, the limited planar measuring points below the undulating surface measuring points were utilized to calculate the magnetic field values of the undulating surface measuring points and the eight corner points around the measuring points through upward continuation. Subsequently, based on the difference between the corner values calculated by continuation and those calculated by forward modeling, trilinear interpolation was used to calculate the magnetic field correction value at the undulating surface measuring point. Finally, the continuation and correction values were summed to obtain the magnetic field value of the undulating surface. The local interpolation formula can be expressed as follows: Δ T c = Q Δ T s = Q 1 + P J − Q 2 Δ T s (3) where vector Δ T c represents the magnetic field on the undulating surface, whose dimension is N Δ T c × 1 , where N Δ T c is the number of observation points on the undulating surface; Q is the interpolation weighting matrix for calculating Δ T c based on the forward field Δ T s of multiple horizontal layer magnetic fields, whose dimension is N Δ T c × N Δ T s , where N Δ T s is the number of observation points of forward field Δ T s of multiple horizontal layer magnetic fields. Both Q 1 and Q 2 are local upward continuation matrices, among which matrix Q 1 calculates the magnetic field of the undulating surface by continuation from the magnetic field on the horizontal plane, whose dimension is N Δ T c × N Δ T s . Matrix Q 2 is used to calculate the magnetic field values of the eight observation points on the horizontal planes around the observation point on the undulating surface by continuation from the magnetic field on the horizontal plane, whose dimension is 8 N Δ T c × N Δ T s . Matrix J is the observation point extraction matrix, which is used to extract the magnetic field values of the eight observation points around the observation point on the undulating surface from matrix Δ T s , whose dimension is 8 N Δ T c × N Δ T s . Matrix P is the trilinear interpolation weighting matrix, which is used to perform weighted summation of the magnetic field values of the eight observation points around the observation point on the undulating surface, whose dimension is N Δ T c × 8 N Δ T c .

The basic formula (Equation 4) for extending from the original spatial point ξ , η , ζ to the calculation point x , y , z is as follows (Zeng, 2005): U x , y , z = ζ − z 2 × ∑ i = 1 m ∑ j = 1 n U ξ i , η j , ζ i j ξ i − x 2 + η j − y 2 + ζ i j i − z 2 3 2 Δ ξ Δ η (4) where Δ ξ and Δ η are the point distances in the x and y directions, respectively.

To ensure a good balance between the accuracy and efficiency of the continuation, the continuation height was generally limited to between one and two times the horizontal grid spacing. The number of grid points used for extension for x and y directions is the same and it is denoted by ns ranging from 32 × 32 to 128 × 128 points. The maximum height of the plane should be higher than the vertex of the undulating surface, and the minimum height of the horizontal plane should be lower than the lowest point of the undulating surface minus one time (Figure 2). This design ensures that the interpolation calculation can fully cover the complex terrain of the undulating surface while ensuring the accuracy and reliability of the calculation results.

Trilinear interpolation was used for error correction. Trilinear interpolation is a linear interpolation method of the tensor-product grid of three-dimensional discrete sampling data (Figure 3). This method linearly approximates the value of the point (x, y, z) on the local rectangular prism through the data points on the grid, and the interpolation formula (Equation 5) as follows: Δ T = Δ T 000 1 − x d 1 − y d 1 − z d + Δ T 010 x d 1 − y d 1 − z d + Δ T 001 1 − x d 1 − y d z d + Δ T 011 x d 1 − y d z d + Δ T 100 1 − x d y d 1 − z d + Δ T 110 x d y d 1 − z d + Δ T 101 1 − x d y d z d + Δ T 111 x d y d z d (5) where Δ T i j k i , j , k ∈ 0 , 1 represents the values of the eight corner points used for linear interpolation; x d , y d ,and z d are intermediate variables, specifically x d = x − x 0 / x 1 − x 0 , y d = y − y 0 / y 1 − y 0 , and z d = z − z 0 / z 1 − z 0 ; x , y , z represents the coordinates of the interpolation point; and x 0 , y 0 , z 0 , x 1 , y 1 , z 1 represents the coordinate range of the eight corner points.

Δ T i j k i , j , k ∈ 0 , 1 represents the values of the eight corner points used for linear interpolation. The red dot represents the interpolation point, and the plane it appears on indicates an interpolation sequence. The interpolation sequence does not affect the interpolation results.

In the calculation, the interpolation weighting matrices, i.e., Q 1 , Q 2 , J , and P , are all sparse matrices. After removing the zero-valued elements in the above matrices, the storage space can be significantly reduced. Among them, the dimension of matrix Q 1 can change from N Δ T c × N Δ T s to N Δ T c × N n s × n s . n s × n s represents the number of observation points used for local continuation. During calculation, the value of p ranging from 32 to 128 can meet the calculation accuracy requirements; the dimension of the matrix Q 2 can change from N Δ T c × N Δ T s to 8 × N Δ T s ; the dimension of the matrix J can change from 8 N Δ T c × N Δ T s to N Δ T c × N n s × n s ; the dimension of the matrix P can change from N Δ T c × 8 N Δ T c to N Δ T c × 8 . By reducing the dimension of the above matrices, the fast calculation can be completed with only a small amount of computation cost.

Based on the magnetic forward model in Equations 2, 3, the Tikhonov (Tikhonov and Arsenin, 1977) regularization objective function, ρ φ α , was constructed as follows: ρ φ α = D Δ T c − Q V s M φ 2 2 + α W M φ − M φ 0 2 2 , (6) where M φ is the expected susceptibility model, α is the regularization factor; D Δ T c − Q V s M φ 2 2 is the fitting difference objective function of the magnetic field, W M φ − M φ 0 2 2 is the Tikhonov regularization objective function, D is the diagonal data covariance matrix, where the elements on its diagonal are the estimated data noise variance, M φ 0 is the reference model, W is the model weighting matrix, and Δ T c represents the measured magnetic field value.

The inverse calculation determines the minimum solution value of the objective function (Equation 6). The gradient at the minimum value of the objective function is 0. The expression for the fixed-point iteration equation of model M φ in the model space can be obtained as follows: M φ = M φ 0 + V s T Q T D T DQ V s + α W T W − 1 V s T Q T D T D Δ T c − Q V s M φ 0 . (7)

In this study, the three-dimensional inversion result of the magnetic susceptibility was obtained by applying the conjugate gradient algorithm to solve Equation 7.

To compare the computational efficiency of the algorithms, models of different scales were used for comparative tests. The size of the prisms used to subdivide the underground space was 2 × 2 × 2 m. The number of subdivide cells of the model ranged from 9,216 to 45, 278, 208 (Table 1), and the latter is 4,913 times greater than the former. The lowest point of the undulating surface was 3 m, the highest point was set at 11 m, the maximum height difference was 8 m, and the spacing between the horizontal observation surfaces was 2 m.

For this series of models, forward calculations were performed using both traditional and fast algorithms. The traditional method uses the superposition principle to calculate the forward values of each observation point one by one and does not introduce any acceleration techniques. In the fast algorithm, the range of data points used for local continuation was set as 64 × 64 observation points. The number of computation points for the forward calculation was consistent with that of single-layer observation surface and they were randomly distributed with the horizontal center point of the observation grid as the center and a grid spacing as the radius (Figure 4). As the calculation time required for the forward calculation using the traditional algorithm exceeded the practical range, obtaining solutions through complete actual calculations was difficult. To effectively estimate the calculation time, the average time required for a single-point forward calculation was obtained by calculating a limited number of observation points and then multiplying by all of the observation points to obtain the estimated time for completing one forward calculation of this model.

The calculation results (Figure 5) show that the fast algorithm effectively improved the calculation efficiency. With an increase in the number of mesh divisions (Table 1), the acceleration effect of the fast algorithm became more notable. In other words, the speedup ratio increased and the maximum speedup ratio in the model example reached 29,167.

We note that the acceleration effect shown in the above calculations was not fixed. The speedup ratio was affected by several factors. One of the most important factors was the ratio between the degree of undulation of the terrain and the distance between the planes. If the elevation difference of the curved surface was greater and the plane spacing was smaller, more planes were required to cover the undulating curved surface, resulting in an increase in the number of calculations and a decrease in the speedup ratio.

The purpose of this experiment was to compare the computational accuracy of three different interpolation methods: one is an interpolation method that only adopts local continuation; the second uses trilinear interpolation to correct all data based on the local continuation interpolation calculation; and the third uses trilinear interpolation to correct only the data at the boundary based on the local continuation interpolation calculation. We compared and analyzed the influence of different data ranges used for local continuation on the computational accuracy of these interpolation methods.

The model used in the test experiment consisted of four cuboidal magnetic bodies (Figure 6a). Table 2 lists the spatial positions and magnetic susceptibilities of the magnetic bodies. The side length of each model cube was 100 m and the elevation range of the underground space was −5,600 to 0 m. The magnetic field dip angle was 60° and the magnetic declination was –9°. The elevation range of the magnetic field undulation observation surface was 205.00 to 1,004.97 m (Figure 6b). The observation points were randomly on an undulating surface with average lateral distance of 100 m in each direction, totaling 39,204 observation points.

The calculation accuracy of the forward field values of the model (Wu, 2018) was evaluated using the relative mean square error formula (Equation 8) as follows: E err = ∑ i Δ T i c − Δ T i 2 ∑ i Δ T i 2 (8) where Δ T i represents the forward value of the magnetic field theory and Δ T i c represents the forward value of the magnetic field based on the fast algorithm.

To compare the influence of different data ranges in the local continuation on the accuracy, the number of local continuation points selected in the experiment ranged from 32 × 32 to 158 × 158. The calculation results (Figure 7) show that the forward calculation accuracy of the three interpolation methods improved with an increase in ns. When only the interpolation method of local continuation was adopted, the relative error value first rapidly decreased from 2% to 0.6% (Local continuation interpolation) and converged. When the trilinear interpolation correction was used to correct all of the data, the relative error was basically stable at approximately 0.5% (Local continuation interpolation + Trilinear 1). When trilinear interpolation was used to correct only the boundary data, the relative error rapidly decreased from 1.6% and finally stabilized at approximately 0.2% (Local continuation interpolation + Trilinear 2).

Based on the forward modeling error distribution map (Figure 8), when the value of ns was relatively small, the forward modeling errors of the three interpolation methods were mainly distributed in the regions with higher magnetic field anomaly values (Figures 8a–c). With an increase in ns, the forward modeling error using only the continuation interpolation method was effectively reduced in the middle area of the data (Figures 8a,d,g), but the error at the data edge did not significantly improve, whereas the forward modeling error when using trilinear interpolation to correct all of the data still remained large in the regions with higher magnetic field anomaly values; however, the error at the data edge was relatively small (Figures 8b,e,h). When using trilinear interpolation to correct only the boundary data for forward modeling errors, it combines the advantages of the above two methods, gradually reduces the forward modeling errors in both the central area and edge, and the error distribution is uniform, ultimately achieving the optimal forward modeling accuracy (Figures 8c,f,i).

The experiments show that the local interpolation correction method adopted in this study achieved relatively high calculation accuracy except at the edges.

In this experiment, a forward-field simulation of the measured magnetic field was performed using the model for the accuracy analysis, as described in Section 2.2. The subsurface space was divided into 200 × 200 × 56, totaling 2.24 million prisms. The prisms were set as cubes with side lengths of 100 m. The magnetic field inclination was 60.0° and the magnetic declination was −9.0°. The forward field of this model combination mainly consisted of four local anomalies with distorted isolines (Figure 6c).

The constraint methods for three-dimensional magnetic inversion mainly adopt depth weighting, focusing weighting, and physical property constraints. The physical property constraint sets the maximum magnetic susceptibility to 0.03 SI. The inversion result (Figure 9) was obtained through multiple calculation iterations. Compared with the theoretical model (Figure 6a), the spatial positions and shapes of the four magnetic anomalies had a good degree of coincidence.

The inversion results were further sliced in the horizontal (z = −1,000 m) and vertical directions (Figure 10). The slicing results show that the susceptibility distribution of the 3-D inversion results was mainly located within the spatial range of the theoretical model (blue-dashed box). The overall agreement with the theoretical model was high. This indicates that the forward and inversion algorithms for the undulating terrain proposed in this study have high inversion accuracy.

To simulate realistic data application (Figure 11a), a region with typical undulating characteristics was selected from the actual terrain and appropriately modified to be the undulating observation surface of the magnetic field, with an elevation range of 420.8–1,402.7 m above sea level (Figure 12a). A sloping plate-shaped magnetic body was set as the target body (Table 3). The body sloped northward with a dip angle of 45°. The dip angle of the magnetic field was 60.0° and the declination angle was −9.0°. The underground space was divided into 70 × 120 × 30, totaling 252,000 prisms, which were set as cubes with a side length of 100 m. The lateral average distance between the observation points on the curved surface was approximately 100 m.

Based on the above model, forward magnetic field calculations were separately performed on curved surfaces and on a plane at the average elevation of the curved surface (822.3 m). The results show the forward modeling on the curved surface was significantly affected by the undulating observation surface compared with the smoother variations of the magnetic field on the plane (Figure 12c). Additionally, random noise (±2 nT) was introduced to evaluate the robustness of the method (Figure 11c).

The three-dimensional inversion adopted the same constraint strategy as in the previous model. The inversion result (Figures 11b,c) was obtained through multiple calculation iterations. The inversion results show a good agreement between the spatial geometry of the inverted magnetic anomalies and the theoretical model (Figure 11a), regardless whether the magnetic field contains noise or not. Furthermore, vertical slices of the inversion results were obtained (Figures 12d–g). The slice results show that the susceptibility distribution of the three-dimensional inversion result was mainly located within the spatial range of the theoretical model (blue-dashed box), regardless whether the noise data is contained. The overall coincidence with the theoretical model was high. This further indicates that the forward and inversion algorithms for undulating terrain proposed in this study also had high inversion accuracy in simulating actual situations.