A high-precision and fast algorithm for three-dimensional magnetic anomalies on undulating terrain

Abstract

This study addresses the challenge of efficiently and accurately computing three-dimensional undulating surface magnetic field data, which has become increasingly difficult with the development of airborne magnetic survey technologies. To solve this problem, we propose a new computational method that is both efficient and highly precise. This method utilizes a planar three-dimensional fast-forward algorithm to calculate multiple planar magnetic fields at different heights, followed by local continuation combined with trilinear interpolation to achieve an accurate conversion from planar magnetic fields to undulating surface magnetic fields. This enables the completion of three-dimensional high-precision fast forward calculations that can then be applied to the inversion algorithm to calculate the distribution of the three-dimensional magnetization rate. Multiple model examples were tested to verify the accuracy and efficiency of the proposed method.

1 Introduction

The magnetic method is a well-developed and widely applied geophysical exploration technique, with extensive use in mineral exploration, engineering geological surveying, geological mapping, and deep structural research (Prutkin and Saleh, 2009; Eppelbaum, 2011; Li and Oldenburg, 2010; David et al., 2019; Kirkby and Duan, 2019). Quantitative interpretation of magnetic anomalies faces challenges in the forward modeling and inversion of three-dimensional magnetic susceptibility models (Bhattacharyya, 1964; Bhattacharyya, 1980; Li and Oldenburg, 1996; Pilkington, 1997; Pilkington, 2008; Li and Sun, 2016; Liu et al., 2018). However, owing to the significant computational resource requirements, the efficiency of three-dimensional inversion is limited, which hampers its widespread application (Portniaguine and Zhdanov, 2000; Davis and Li, 2013; Rezaie et al., 2017).

Recently, researchers have proposed enhanced methodologies and targeted computational strategies to address this issue. When the observation point is situated on a regular grid in the plane, significant reductions in the computational and storage requirements can be achieved through equivalent calculations based on the symmetry of the geometric framework (Shengxian et al., 2018). Building on this foundation, further optimization algorithms have been developed by leveraging the fast Fourier Transform (FFT) to execute convolution operations in the spatial domain via multiplication operations in the frequency domain. These algorithms optimize both spatial computation accuracy and frequency-domain computation speed (Chen and Liu, 2019; Jing et al., 2019; Hogue et al., 2019; Yuan et al., 2022), enabling efficient and accurate three-dimensional forward and inversion computation.

When the observation point is situated on an undulating observation surface, the symmetry between the observation point and field source is disrupted, rendering direct utilization of the aforementioned algorithm unfeasible. To improve computational efficiency, Zheng et al. (2019) enhanced Parker’s formula in the wavenumber domain, achieving rapid and high-precision gravity forward modeling on undulating interfaces, as well as fast and accurate gravity forward simulation of density interface models with topography. In the spatial wavenumber domain, a numerical simulation algorithm for gravity anomalies based on an arbitrary Fourier transform using quadratic interpolation functions was developed; CPU-GPU acceleration was employed to further amplify the advantages of three-dimensional gravity anomaly forward modeling in this domain (Dai et al., 2024). In terms of spatial domain calculations, multiple planar fields encompassing undulating surfaces were computed using efficient algorithms for forward modeling. The forward field values at any given observation points are obtained by using local interpolation techniques (Li et al., 2018; Qiang et al., 2019). Although Li et al. (2018) utilized a cubic spline interpolation method capable of achieving a precise conversion from a planar magnetic field to an undulating surface magnetic field, its non-linear nature poses challenges when directly applied to three-dimensional linear inversion.

This study proposes a high-precision local linear interpolation method to address the aforementioned issues. This method initially employs physically meaningful local extrapolation for interpolation, subsequently enhancing the accuracy of interpolation through trilinear interpolation correction. Integrating this composite linear interpolation approach with a rapid algorithm for the magnetic field plane facilitates the efficient forward and inverse computations of large-scale aeromagnetic data.

2 Methods

2.1 Fast forward modeling algorithm for planar magnetic field

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

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).where vector represents the forward modeling anomaly of the model with a dimension of , is the magnetization intensity according to the forward formula, also with a dimension of , represents the forward operator that connects observation data and physical parameters , known as the forward coefficient matrix, which has a dimension of , represents the fast Fourier transform in the positive direction while represents the fast Fourier transform in the inverse direction, ()denotes the prism forward kernel matrix at depth k, and () 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 elements in matrix represent the density matrix elements of the kth layer in the model while the remaining elements are set to zero. The first elements in matrix correspond to the magnetic field response values at surface observation points for prism with index number in the density matrix of the kth layer (with default starting index as 1). At this time, prism has a density set to the unit magnetization intensity. The remaining elements in matrix are stored symmetrically with respect to index number .

2.2 Fast forward algorithm for magnetic anomaly of undulating surfaces

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

FIGURE 2

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:where vector represents the magnetic fields on the multiple horizontal observation planes and its dimension is , i.e., the total number of observation points on multiple horizontal observation planes. Furthermore, represents the forward coefficient matrix for calculating the magnetic field values of multiple horizontal observation planes, whose dimension is .

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:where vector represents the magnetic field on the undulating surface, whose dimension is , where is the number of observation points on the undulating surface; is the interpolation weighting matrix for calculating based on the forward field of multiple horizontal layer magnetic fields, whose dimension is , where is the number of observation points of forward field of multiple horizontal layer magnetic fields. Both and are local upward continuation matrices, among which matrix calculates the magnetic field of the undulating surface by continuation from the magnetic field on the horizontal plane, whose dimension is . Matrix 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 . Matrix 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 , whose dimension is . Matrix 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 .

The basic formula (Equation 4) for extending from the original spatial point to the calculation point is as follows (Zeng, 2005):where and are the point distances in the and 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:where represents the values of the eight corner points used for linear interpolation; , ,and are intermediate variables, specifically , , and ; represents the coordinates of the interpolation point; and represents the coordinate range of the eight corner points.

FIGURE 3

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., , , , and , 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 can change from to . 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 can change from to ; the dimension of the matrix can change from to ; the dimension of the matrix can change from to . By reducing the dimension of the above matrices, the fast calculation can be completed with only a small amount of computation cost.

2.3 Fast algorithm for magnetic anomaly of undulating surface

Based on the magnetic forward model in Equations 2, 3, the Tikhonov (Tikhonov and Arsenin, 1977) regularization objective function, , was constructed as follows:where is the expected susceptibility model, is the regularization factor; is the fitting difference objective function of the magnetic field, is the Tikhonov regularization objective function, is the diagonal data covariance matrix, where the elements on its diagonal are the estimated data noise variance, is the reference model, is the model weighting matrix, and 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 in the model space can be obtained as follows:

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

3 Theoretical model experiments

Theoretical model tests were conducted for computation and testing on a computer with a processor speed of 2.20 GHz and a memory of 64 GB. This was the mainstream computer configuration for notebook workstations.

3.1 Efficiency comparison of forward modeling calculation

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.

FIGURE 4

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.

FIGURE 5

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.

3.2 Accuracy analysis of forward modeling computation

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.

FIGURE 6

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:where represents the forward value of the magnetic field theory and 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).

FIGURE 7

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

FIGURE 8

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

3.3 The model study

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.

FIGURE 9

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.

FIGURE 10

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.

FIGURE 11

FIGURE 12

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.

4 Conclusion

The algorithm proposed in this study has high potential for parallelism. The execution efficiency of this algorithm should improve with the introduction of parallel-computing technology. It should be noted that this method may not work well if there is significant remanent magnetization or the case that the magnetization is very high or anisotropy.

References

• 1

  BhattacharyyaB. K. (1964). Magnetic anomalies due to prism‐shaped bodies with arbitrary polarization. Geophysics29 (4), 517–531. 10.1190/1.1439386

  • CrossRef
  • Google Scholar

• 2

  BhattacharyyaB. K. (1980). A generalized multibody model for inversion of magnetic anomalies. Geophysics45, 255–270. 10.1190/1.1441081

  • CrossRef
  • Google Scholar

• 3

  ChenL.LiuL. (2019). Fast and accurate forward modelling of gravity field using prismatic grids. Geophys. J. Int.2019 (2), 1062–1071. 10.1093/gji/ggy480

  • CrossRef
  • Google Scholar

• 4

  DaiS.ZhuD.YingZ.KunL.ChenQ.JiaxuanL.et al (2024). Three dimensional fast forward modeling of gravity anomalies under arbitrary undulating terrain. Chin. J. Geophys. (In Chinese)67 (2), 768–780. 10.6038/Cjg2023q0605

  • CrossRef
  • Google Scholar

• 5

  DavidW. M.MejerH. T.ArneD. (2019). Inference of unexploded ordnance (UXO) by probabilistic inversion of magnetic data. Geophysical Journal International1 (1). 10.1093/gji/ggz421

  • CrossRef
  • Google Scholar

• 6

  DavisK.LiY. (2013). Efficient 3D inversion of magnetic data via octree-mesh discretization, space-filling curves, and wavelets. Geophysics78 (5), J61–J73. 10.1190/GEO2012-0192.1

  • CrossRef
  • Google Scholar

• 7

  EppelbaumL. V. (2011). Study of magnetic anomalies over archaeological targets in urban environments. Physics and Chemistry of the Earth36 (16), 1318–1330. 10.1016/j.pce.2011.02.005

  • CrossRef
  • Google Scholar

• 8

  HogueJ. D.RenautR. A.VatankhahS. (2019). A tutorial and open source software for the efficient evaluation of gravity and magnetic kernels. Computers and Geosciences144, 104575. 10.1016/j.cageo.2020.104575

  • CrossRef
  • Google Scholar

• 9

  JingL.YaoC. L.YangY. B.XuM. L.ZhangG. Z.JiR. Y. (2019). Optimization algorithm for rapid 3D gravity inversion. Applied Geophysics16 (4), 507–518. 10.1007/s11770-019-0781-2

  • CrossRef
  • Google Scholar

• 10

  KirkbyA.DuanJ. (2019). Crustal structure of the eastern arunta region, central Australia, from magnetotelluric, seismic, and magnetic data. Journal of Geophysical Research Solid Earth124 (8), 9395–9414. 10.1029/2018JB016223

  • CrossRef
  • Google Scholar

• 11

  LiK.Long-WeiC.ChenQ.DaiS. K.ZhangQ. J.ZhaoD. D.et al (2018). Fast 3D forward modeling of the magnetic field and gradient tensor on an undulated surface. Applied Geophysics015 (003), 500–512. 10.1007/s11770-018-0690-9

  • CrossRef
  • Google Scholar

• 12

  LiY.OldenburgD. W. (2010). Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method. Geophysical Journal of the Royal Astronomical Society152 (2), 251–265. 10.1046/j.1365-246X.2003.01766.x

  • CrossRef
  • Google Scholar

• 13

  LiY.SunJ. (2016). 3D magnetization inversion using fuzzy c-means clustering with application to geology differentiation. Geophysics81 (5), J61–J78. 10.1190/geo2015-0636.1

  • CrossRef
  • Google Scholar

• 14

  LiY. G.OldenburgD. W. (1996). 3-D inversion of magnetic data. Geophysics61 (2), 394–408. 10.1190/1.1443968

  • CrossRef
  • Google Scholar

• 15

  LiuS.MaurizioF.XiangyunH.JamaledinB.BangshunW.ZhangD.et al (2018). Extracting induced and remanent magnetizations from magnetic data modeling. Journal of geophysical research. Solid earth123 (11), 9290–9309. 10.1029/2017jb015364

  • CrossRef
  • Google Scholar

• 16

  PilkingtonM. (1997). 3-D magnetic imaging using conjugate gradients. Seg Technical Program Expanded Abstracts15 (1), 1415–1418. 10.1190/1.1826377

  • CrossRef
  • Google Scholar

• 17

  PilkingtonM. (2008). 3D magnetic data-space inversion with sparseness constraints. Geophysics74 (1), L7–L15. 10.1190/1.3026538

  • CrossRef
  • Google Scholar

• 18

  PortniaguineO.ZhdanovM. S. (2000). 3-D magnetic inversion with data compression and image focusing. Geophysics67 (5), 1532–1541. 10.1190/1.1512749

  • CrossRef
  • Google Scholar

• 19

  PrutkinI.SalehA. (2009). Gravity and magnetic data inversion for 3D topography of the moho discontinuity in the Northern Red Sea area, Egypt. Journal of Geodynamics47 (5), 237–245. 10.1016/j.jog.2008.12.001

  • CrossRef
  • Google Scholar

• 20

  QiangJ.ZhangW.LuK.ChenL.ZhuY.HuS.et al (2019). A fast forward algorithm for three-dimensional magnetic anomaly on undulating terrain. Journal of Applied Geophysics166, 33–41. 10.1016/j.jappgeo.2019.04.009

  • CrossRef
  • Google Scholar

• 21

  RezaieM.MoradzadehA.KalatehA. N. (2017). Fast 3D inversion of gravity data using solution space priorconditioned lanczos bidiagonalization. Journal of Applied Geophysics136, 42–50. 10.1016/j.jappgeo.2016.10.019

  • CrossRef
  • Google Scholar

• 22

  ShengxianL.CenterC.SurveyC. G. (2018). A self-constrained 3D inversion and efficient solution of gravity data based on cross-correlation coefficient. Journal of Jilin University(Earth Science Edition). 10.13278/j.cnki.jjuese.20170167

  • CrossRef
  • Google Scholar

• 23

  TikhonovA. N.ArseninV. Y. (1977). Solutions of Ill-posed problems. Winston, New York.

  • Google Scholar

• 24

  WuL. (2018). Efficient modelling of gravity effects due to topographic masses using the Gauss–FFT method. Geophysical Journal International205 (1), 160–178. 10.1093/gji/ggw010

  • CrossRef
  • Google Scholar

• 25

  YuanY.CuiY. A.ChenB.ZhaoG. D.LiuJ. X.GuoR. W.et al (2022). Fast and high accuracy 3D magnetic anomaly forward modeling based on BTTB matrix. Chin. J. Geophys. (in Chinese)65 (3), 1107–1124. 10.6038/cjg2022P0126

  • CrossRef
  • Google Scholar

• 26

  ZengH. L. (2005). Gravity fields and gravity exploration. Beijing: Geology Press.

  • Google Scholar

• 27

  ZhengX.LiuJ.ChenB.ZhaoG.ChenL.GuoR.et al (2019). A fast and accurate gravity forward method for undulating structure model with surface observation. Computing techniques for geophysical and geochemical exploration.

  • Google Scholar