A three-dimensional fast forward and inversion method and system for magnetic anomaly on undulating surface

Abstract

The present invention discloses a method and system for three-dimensional rapid forward and inversion of magnetic anomalies on an undulating surface, which belongs to the field of geological exploration technology. The method obtains magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and obtains magnetic anomaly values ​​of measuring points on the undulating surface through a local extension method interpolation algorithm based on the magnetic field values ​​of spatial observation grid points on the multiple horizontal observation surfaces. The magnetic anomaly values ​​of measuring points on the undulating surface are reversely weighted to determine the grid point matrix of the spatial observation grid points. Based on the conventional rapid inversion algorithm formula, the conjugate gradient algorithm is used to perform inversion calculation on the magnetic anomaly field to complete the rapid forward and inversion of the magnetic anomaly on the undulating surface. The upward extension algorithm is used to realize local interpolation, which not only has practical physical significance, but also can keep the interpolation results stable and reliable, effectively guaranteeing the overall calculation accuracy, and by making full use of the advantages of the abnormal fast algorithm, the calculation efficiency of three-dimensional inversion of magnetic anomalies on the undulating surface is greatly improved.

Description

Technical Field

[0001] The invention belongs to the technical field of geological exploration, and in particular relates to a three-dimensional rapid forward and inversion method and system for magnetic anomaly on an undulating surface. Background Art

[0002] Potential field (gravity field, geomagnetic field) is a function of spatial position, and its amplitude changes with altitude. Potential field data obtained by ground measurement, ship measurement, aerial measurement and other means are usually located on undulating surfaces. Changes in altitude will complicate the changes in potential field data and increase the difficulty of subsequent data interpretation. In actual aeromagnetic measurements, the observation surface is often an undulating surface, and the number of observation points is large. The traditional three-dimensional inversion algorithm for magnetic anomalies is limited by low computational efficiency and is difficult to carry out inversion calculations of large amounts of data. The magnetic fast forward and inversion algorithm is based on the measurement points being located on a regular plane grid. Therefore, it is also difficult to directly use the fast algorithm to carry out three-dimensional inversion of the magnetic field on the undulating surface. Summary of the invention

[0003] In view of this, the present invention provides a three-dimensional rapid forward and inversion method and system for undulating surface magnetic anomalies, which adopts an upward continuation algorithm. It not only has practical physical significance, but also the interpolation results are stable and reliable, effectively ensuring the overall calculation accuracy of the algorithm to solve the above-mentioned technical problems, and is specifically implemented by the following technical scheme.

[0004] In a first aspect, the present invention provides a three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly, comprising the following steps:

[0005] Acquire magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determine magnetic field values ​​of multiple spatial observation grid points on the horizontal observation surfaces according to the magnetic anomaly fields, wherein the spatial observation grid points are observation points on multiple horizontal observation surfaces, the observation points include horizontal surface measurement points and undulating surface measurement points, the vertical positions of the spatial observation grid points are located on the multiple horizontal observation surfaces, and the horizontal positions of the spatial observation grid points coincide with the horizontal projections of the centers of the underground subdivided grids;

[0006] Performing forward calculation on the observation point according to the magnetic field value and a local extension interpolation algorithm to obtain a magnetic anomaly value of the measurement point on the undulating surface;

[0007] Inversely weighting the magnetic anomaly values ​​of the undulating surface measurement points to determine a grid point matrix of the spatial observation grid points;

[0008] The magnetic anomaly field is inverted and calculated based on the grid point matrix and the conjugate gradient algorithm to complete the rapid forward and inversion of the magnetic anomaly of the undulating surface.

[0009] As a preferred embodiment of the above technical solution, forward calculation is performed on the observation point according to the magnetic field value and the local continuation interpolation algorithm to obtain the magnetic anomaly value of the measurement point on the undulating surface, including:

[0010] The three-dimensional magnetic field fast forward algorithm is used to calculate the magnetic field of the horizontal plane measuring point to obtain the magnetic field forward value of the horizontal observation plane;

[0011] The magnetic field of the measuring point on the undulating surface is calculated based on the magnetic field forward modeling value and the local continuation algorithm to obtain the magnetic anomaly value of the measuring point on the undulating surface.

[0012] As a preferred embodiment of the above technical solution, a three-dimensional magnetic field fast forward algorithm is used to calculate the magnetic field of the horizontal plane measuring point to obtain the magnetic field forward value of the horizontal observation plane, including:

[0013] The underground inversion space is gridded, and a plurality of horizontal magnetic field observation planes are set to obtain a horizontal observation plane, wherein the plurality of horizontal magnetic field observation planes are arranged at equal intervals according to a longitudinal length dz of the underground space grid;

[0014] The magnetic field value of the horizontal observation surface is calculated by using a three-dimensional magnetic field fast algorithm;

[0015] Calculate an interpolation weight matrix of a local extension of the observation point based on the undulating surface measurement point and the horizontal observation surface;

[0016] The magnetic field value of the measuring point on the undulating surface is obtained according to the magnetic field value of the horizontal observation surface and the interpolation weight matrix.

[0017] As a preferred embodiment of the above technical solution, multiple horizontal magnetic field observation surfaces are arranged at equal intervals according to the longitudinal length of the underground space, including:

[0018] The maximum height of the horizontal observation plane is higher than the vertex of the undulating surface corresponding to the undulating surface measuring point minus one times the longitudinal length dz, and the minimum height of the horizontal observation plane is lower than the lowest point of the undulating surface minus one times the longitudinal length dz.

[0019] As a preferred embodiment of the above technical solution, obtaining the magnetic field value of the measuring point on the undulating surface according to the magnetic field value of the horizontal observation surface and the interpolation weight matrix includes:

[0020] The magnetic field value of the measuring point on the undulating surface is forward calculated using the volume integral formula, and the corresponding expression is:

[0021] (1)

[0022] (2)

[0023] (3)

[0024] Among them, when the magnetic body is a rectangular parallelepiped, is the total magnetic field; are the coordinates of the grid points on the observation plane; is the coordinate of the center point of the cuboid; a, b, c are the extension lengths of the cuboid along the x, y, and z directions; M is the total magnetic susceptibility of the cuboid; is the magnetic permeability in vacuum; is the total magnetic susceptibility inclination; is the total magnetic susceptibility declination; is the inclination angle of the Earth's magnetic field; is the declination angle of the geomagnetic field; They represent the unit projections of the magnetization direction of the observation point in the x, y, and z directions respectively; Represents the second-order derivative of the magnetic potential of the cuboid in the x, y, z directions and their combined directions, ; are the components of the magnetic field in the x, y, and z directions respectively; represents pi; They represent the components of the magnetic susceptibility M of the cuboid in the x, y, and z directions respectively; They represent the unit projections of the magnetization direction of the cuboid in the x, y, and z directions respectively; are the coordinates of the field source point inside the cuboid;

[0025] (4)

[0026] As a preferred embodiment of the above technical solution, when the magnetic inclination and the magnetic declination are constants, the magnetic field forward modeling formula can be simplified as follows:

[0027] (5)

[0028] Among them, V is the forward coefficient matrix for calculating the magnetic field values ​​of multiple horizontal observation surfaces; M is the model magnetic susceptibility matrix; Q represents the interpolation weight matrix.

[0029] As a preferred embodiment of the above technical solution, the magnetic anomaly field is inverted and calculated based on the grid point matrix and the conjugate gradient algorithm to complete the rapid forward and inversion of the undulating surface magnetic anomaly, including:

[0030] The solution to the inverse problem is obtained by minimizing the Tikhonov regularized objective function The corresponding expression is:

[0031] (6)

[0032] in, is the desired susceptibility model, is the regularization factor; is the fitting error objective function of the magnetic field; is the Tikhonov regularization objective function; D is the diagonal data covariance matrix, whose diagonal elements are the estimated data noise variance; is the reference model; W is the model weighting matrix; represents the measured magnetic field vector; T represents the transpose operation;

[0033] According to the inversion calculation, the minimum value of the formula (6) corresponding to the objective function is obtained. According to the gradient of the minimum value of the objective function being 0, the fixed point iteration equation of the model M is obtained in the model space and data space. The corresponding calculation expression is as follows:

[0034] (7)

[0035] Among them, when formula (7) is solved using the conjugate gradient algorithm, its main calculation amount is the matrix The product with vector P, P represents the intermediate vector;

[0036] The calculation form of formula (7) includes and , column vector u and column vector v are intermediate variables and do not have fixed physical meanings.

[0037] As a preferred embodiment of the above technical solution, The vector v of the term is regarded as the model vector M, and its calculation process is regarded as a forward calculation, which is achieved by reverse weighting Calculation of items, including:

[0038] According to the interpolation weight matrix Q of the local extension of the measuring point, the value of the vector u is reversely weighted to multiple horizontal plane measuring points to obtain the pseudo magnetic field vector ;

[0039] The pseudo magnetic field vector The pseudo model vector is calculated by a fast algorithm , thus completing Calculation of items.

[0040] As a preferred embodiment of the above technical solution, the fast forward algorithm is accelerated by using the FFT algorithm, and the observation point in the FFT algorithm coincides with the center point of the underground subdivision grid in the horizontal projection.

[0041] In a second aspect, the present invention further provides a three-dimensional rapid forward and inversion system for undulating surface magnetic anomaly, which is applied to the above-mentioned three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly, comprising:

[0042] A measuring point determination module, used to obtain magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determine the magnetic field values ​​of the spatial observation grid points on the multiple horizontal observation surfaces according to the magnetic anomaly fields, wherein the spatial observation grid points are observation points on the multiple horizontal observation surfaces, the observation points include horizontal surface measurement points and undulating surface measurement points, the vertical positions of the spatial observation grid points are located on the multiple horizontal observation surfaces, and the horizontal positions of the spatial observation grid points coincide with the horizontal projections of the centers of the underground subdivision grids;

[0043] A forward calculation module, used for performing forward calculation on the observation point according to the magnetic field value and a local extension interpolation algorithm to obtain a magnetic anomaly value of the measurement point on the undulating surface;

[0044] A reverse weighting module, used for reverse weighting the magnetic anomaly values ​​of the undulating surface measurement points to determine the grid point matrix of the spatial observation grid points;

[0045] The inversion calculation module is used to perform inversion calculation on the magnetic anomaly field based on the grid point matrix and the conjugate gradient algorithm to complete the rapid forward and inversion of the undulating surface magnetic anomaly.

[0046] The present invention provides a three-dimensional fast forward and inversion method and system for magnetic anomaly of an undulating surface. The method obtains magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determines magnetic field values ​​of multiple spatial observation grid points on the horizontal observation surfaces according to the magnetic anomaly fields. The observation points are forward calculated according to the magnetic field values ​​and a local continuation interpolation algorithm to obtain magnetic anomaly values ​​of the undulating surface measurement points. The magnetic anomaly values ​​of the undulating surface measurement points are reversely weighted to determine a grid point matrix of the spatial observation grid points. The magnetic anomaly fields are inverted based on the grid point matrix and a conjugate gradient algorithm to complete the fast forward and inversion of the magnetic anomaly of the undulating surface. The upward continuation algorithm is used to realize local interpolation, which not only has practical physical significance, but also can keep the interpolation results stable and reliable, effectively guarantees the overall calculation accuracy, and greatly improves the calculation efficiency of the three-dimensional inversion of the undulating surface magnetic anomaly by making full use of the advantages of the abnormal fast algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

[0048] Figure 1 A flow chart of the three-dimensional rapid forward and inversion method for magnetic anomaly on an undulating surface provided by the present invention;

[0049] Figure 2A coordinate schematic diagram of a magnetic cuboid provided by the present invention;

[0050] Figure 3 A schematic diagram of the corresponding relationship between the model unit grid and the measuring point grid provided by the present invention;

[0051] Figure 4 A schematic diagram of the relationship between the undulating observation surface and the horizontal observation surface provided by the present invention;

[0052] Figure 5 A schematic diagram of the distribution of arbitrary measuring points on the undulating surface provided by the present invention;

[0053] Figure 6 A comparison chart of the calculation time of the traditional algorithm and the fast algorithm provided by the present invention;

[0054] Figure 7 A schematic diagram of a model provided by the present invention;

[0055] Figure 8 A schematic diagram for comparing the calculation accuracy provided by the present invention;

[0056] Fig. 9 A schematic diagram of a magnetic body model provided by the present invention;

[0057] Fig.10 A three-dimensional schematic diagram of the undulating observation surface and measuring points provided by the present invention;

[0058] Fig.11 The forward field diagram of the model magnetic field theory provided by the present invention;

[0059] Fig.12 A schematic diagram of the three-dimensional inversion result of magnetic susceptibility provided by the present invention;

[0060] Fig.13 A slice diagram of the three-dimensional inversion result of magnetic susceptibility provided by the present invention;

[0061] Fig.14 This is a structural block diagram of the three-dimensional rapid forward and inversion system for undulating surface magnetic anomalies provided by the present invention. DETAILED DESCRIPTION

[0062] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and cannot be understood as limiting the present invention.

[0063] See also Figure 1 In a first aspect, the present invention provides a method for rapid three-dimensional forward and inversion of undulating surface magnetic anomalies, comprising the following steps:

[0064] S1: Acquire magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determine the magnetic field values ​​of multiple spatial observation grid points on the horizontal observation surfaces according to the magnetic anomaly fields, wherein the spatial observation grid points are observation points on multiple horizontal observation surfaces, the observation points include horizontal surface measurement points and undulating surface measurement points, the vertical positions of the spatial observation grid points are located on the multiple horizontal observation surfaces, and the horizontal positions of the spatial observation grid points coincide with the horizontal projections of the centers of the underground subdivision grids;

[0065] S2: forward calculating the observation point according to the magnetic field value and the local extension interpolation algorithm to obtain the magnetic anomaly value of the measurement point on the undulating surface;

[0066] S3: performing inverse weighting on the magnetic anomaly values ​​of the undulating surface measurement points to determine a grid point matrix of the spatial observation grid points;

[0067] S4: Based on the grid point matrix and the conjugate gradient algorithm, the magnetic anomaly field is inverted and calculated to complete the rapid forward and inversion of the undulating surface magnetic anomaly.

[0068] See also Figure 2 , Figure 3 and Figure 4 In this embodiment, forward calculation is performed on the observation point according to the magnetic field value and the local extension interpolation algorithm to obtain the magnetic anomaly value of the undulating surface measurement point, including: using a three-dimensional magnetic field fast forward algorithm to perform magnetic field calculation on the horizontal surface measurement point to obtain the magnetic field forward value of the horizontal observation surface; based on the magnetic field forward value and the local extension algorithm, the magnetic field is calculated on the undulating surface measurement point to obtain the magnetic anomaly value of the undulating surface measurement point. Using a three-dimensional magnetic field fast forward algorithm to perform magnetic field calculation on the horizontal surface measurement point to obtain the magnetic field forward value of the horizontal observation surface includes: gridding the underground inversion space, and setting multiple horizontal magnetic field observation surfaces to obtain a horizontal observation surface, wherein the multiple horizontal magnetic field observation surfaces are arranged at equal intervals according to the longitudinal length dz of the underground space division; using a three-dimensional magnetic field fast algorithm to calculate the magnetic field value of the horizontal observation surface; based on the undulating surface measurement point and the local extension interpolation weight matrix of the observation point calculation observation point on the horizontal observation surface; according to the magnetic field value of the horizontal observation surface and the interpolation weight matrix, the magnetic field value of the undulating surface measurement point is obtained. Multiple horizontal magnetic field observation surfaces are arranged at equal intervals according to the longitudinal length of the underground space division, including: the maximum height of the horizontal surface is higher than the vertex of the undulating surface corresponding to the undulating surface measuring point minus one times the longitudinal length dz, and the lowest height of the horizontal observation surface is lower than the lowest point of the undulating surface minus one times the longitudinal length dz.

[0069] See also Fig.10It should be noted that the magnetic body can be a cuboid or other three-dimensional body. For the magnetic field calculation of any observation point on the undulating surface, the three-dimensional magnetic field fast forward algorithm is used to calculate the forward values ​​of multiple horizontal layer magnetic fields, and then the magnetic field value of the observation point is obtained by the local extension algorithm. Since the three-dimensional magnetic field fast forward algorithm uses FFT for acceleration, it has high computational efficiency, but the algorithm requires that the position of the observation point is relatively strict and needs to coincide with the center point of the underground subdivision grid in the horizontal projection. For the setting of multiple horizontal observation surfaces, in order to facilitate calculation, they are generally arranged at equal intervals according to the longitudinal length dz of the underground space subdivision grid, where the maximum height of the horizontal observation surface is higher than the vertex of the undulating surface minus one dz; the lowest height of the horizontal observation surface is lower than the lowest point of the undulating surface minus one dz. The data interpolation in the forward algorithm adopts the algorithm of local upward extension. In order to take into account the extension accuracy and extension efficiency, the extension height is generally limited to between 1 times dz and 2 times dz. The grid range used for extension is 64×64 measurement points to 128×128 measurement points.

[0070] It should be understood that by obtaining the magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determining the magnetic field values ​​of multiple spatial observation grid points on the horizontal observation surfaces based on the magnetic anomaly fields, forward modeling is performed on the observation points based on the magnetic field values ​​and the local continuation interpolation algorithm to obtain the magnetic anomaly values ​​of the undulating surface measurement points, the magnetic anomaly values ​​of the undulating surface measurement points are reversely weighted to determine the grid point matrix of the spatial observation grid points, and the magnetic anomaly field is inverted based on the grid point matrix and the conjugate gradient algorithm to complete the rapid forward and inversion of the magnetic anomaly of the undulating surface. The upward continuation algorithm is used to achieve local interpolation, which not only has practical physical significance, but also can keep the interpolation results stable and reliable, effectively ensuring the overall calculation accuracy, and by making full use of the advantages of the abnormally fast algorithm, the calculation efficiency of the three-dimensional inversion of the undulating surface magnetic anomaly is greatly improved.

[0071] Optionally, obtaining the magnetic field value of the measuring point on the undulating surface according to the magnetic field value of the horizontal plane and the interpolation weight matrix includes:

[0072] The magnetic field value of the measuring point on the undulating surface is forward calculated using the volume integral formula, and the corresponding expression is:

[0073] (1)

[0074] (2)

[0075] (3)

[0076] Among them, when the magnetic body is a rectangular parallelepiped, is the total magnetic field; are the coordinates of the grid points on the observation plane; is the coordinate of the center point of the cuboid; a, b, c are the extension lengths of the cuboid along the x, y, and z directions; M is the total magnetic susceptibility of the cuboid; is the magnetic permeability in vacuum; is the total magnetic susceptibility inclination; is the total magnetic susceptibility declination; is the inclination angle of the Earth's magnetic field; is the declination angle of the geomagnetic field; They represent the unit projections of the magnetization direction of the observation point in the x, y, and z directions respectively; Represents the second-order derivative of the magnetic potential of the cuboid in the x, y, z directions and their combined directions, ; are the components of the magnetic field in the x, y, and z directions respectively;

[0077] (4)

[0078] in, represents pi; They represent the components of the magnetic susceptibility M of the cuboid in the x, y, and z directions respectively; They represent the unit projections of the magnetization direction of the cuboid in the x, y, and z directions respectively; are the coordinates of the field source point inside the cuboid.

[0079] In this embodiment, when the magnetic inclination and the magnetic declination are constants, the magnetic field forward modeling formula can be simplified as follows:

[0080] (5)

[0081] Among them, V is the forward coefficient matrix for calculating the magnetic field values ​​of multiple horizontal observation surfaces; M is the model magnetic susceptibility matrix; Q represents the interpolation weight matrix. In other words, the specific calculation process is: first, the underground inversion space is reasonably gridded, and then multiple horizontal magnetic field observation surfaces are set; based on the three-dimensional magnetic field fast algorithm, the magnetic field values ​​of multiple horizontal surfaces are quickly calculated. ; Based on the spatial position of the measuring point on the undulating surface and the horizontal observation surface, the interpolation weight matrix Q for the local extension of the measuring point is calculated and stored; according to the magnetic field values ​​of multiple horizontal planes And the interpolation weight matrix Q, calculate the magnetic field value of the measuring point on the undulating surface .

[0082] Optionally, performing inversion calculation on the magnetic anomaly field based on the grid point matrix and the conjugate gradient algorithm to complete rapid forward and inversion of the undulating surface magnetic anomaly includes:

[0083] The solution to the inverse problem is obtained by minimizing the Tikhonov regularized objective function The corresponding expression is:

[0084] (6)

[0085] in, is the desired susceptibility model, is the regularization factor; is the fitting error objective function of the magnetic field; is the Tikhonov regularization objective function; D is the diagonal data covariance matrix, whose diagonal elements are the estimated data noise variance; is the reference model; W is the model weighting matrix; represents the measured magnetic field vector; T represents the transpose operation;

[0086] According to the inversion calculation, the minimum value of the formula (6) corresponding to the objective function is obtained. According to the gradient of the minimum value of the objective function being 0, the fixed point iteration equation of the model M is obtained in the model space and data space. The corresponding calculation expression is as follows:

[0087] (7)

[0088] Among them, when formula (7) is solved using the conjugate gradient algorithm, its main calculation amount is the matrix The product with the vector P, P represents the intermediate vector;

[0089] The calculation form of formula (7) includes and , column vector u and column vector v are intermediate variables and do not have fixed physical meanings.

[0090] In this embodiment, P represents a temporary vector. It has different meanings in different calculation steps, but the vector size is the same. It is mainly used to express the vector in the multiplication of the matrix and the vector, and has the same function as the vectors u and v. The vector v of the term is regarded as the model vector M, and its calculation process is regarded as a forward calculation, which is achieved by reverse weighting Calculation of items, including:

[0091] According to the interpolation weight matrix Q of the local extension of the measuring point, the value of the vector u is reversely weighted to multiple horizontal plane measuring points to obtain the pseudo magnetic field vector ;

[0092] The pseudo magnetic field vector The pseudo model vector is calculated by a fast algorithm , thus completing Calculation of items.

[0093] It should be noted that the fast forward algorithm uses the FFT algorithm for acceleration, and the observation point in the FFT algorithm coincides with the center point of the underground subdivision grid in the horizontal projection. The solution to the inverse problem is generally achieved by minimizing the Tikhonov regularization objective function. Regularization means that in linear algebra theory, ill-posed problems are usually defined by a set of linear algebraic equations, and this set of equations usually comes from ill-posed inverse problems with large condition numbers. Large condition numbers mean that rounding errors or other errors will seriously affect the results of the problem. FFT (Fast Fourier Transformation) is a fast Fourier transform, which is a fast algorithm for discrete Fourier transform. It is obtained by improving the discrete Fourier transform algorithm based on the odd, even, imaginary, and real characteristics of discrete Fourier transform, thereby effectively improving the algorithm's computational efficiency.

[0094] Specifically, in order to effectively explain the calculation process of the above-mentioned rapid forward and inversion of the undulating surface magnetic anomaly and propose a theoretical model test, the theoretical model test was tested and calculated on a computer with a main frequency of 2.20GHz and a memory of 64GB. The computer configuration is a mainstream notebook workstation configuration. The theoretical model test includes forward calculation efficiency comparison, forward calculation accuracy analysis and inversion calculation effect analysis. The specific process is as follows:

[0095] Comparison of forward modeling efficiency: In order to compare the computational efficiency of the algorithms, models of different sizes were set up for comparative testing. The size of the prism used to divide the underground space is 2m×2m×2m, and the number of model grid divisions ranges from 9216 to 45278208, which is about 4913 times the difference between the two. The lowest point of the undulating surface where the measuring point is located is 3m, the highest point is set to 11m, the maximum height difference is 8m, and the spacing between multiple horizontal observation surfaces is 2m. The traditional algorithm and the fast algorithm are used to perform forward modeling on this series of models. The range of data points used for local extension in the fast algorithm is 64×64 measuring points. The number of observation points for the forward calculation is consistent with the number of single-layer grid divisions, and is randomly distributed with the horizontal center point of the grid as the center and a radius of one times the grid spacing (such as Figure 5 ). When using traditional algorithms for forward modeling, the required calculation time exceeds the practical scope and is difficult to obtain through complete actual calculation. In order to effectively estimate the calculation time, the average time required for single-point forward modeling is obtained by calculating the limited observation points, and then multiplied by all observation points to obtain the estimated time for the model to complete a forward modeling. The calculation results show (such as Figure 6 ) The fast algorithm can effectively improve the computational efficiency, and as the number of mesh divisions increases, the acceleration effect of the fast algorithm becomes more obvious (as shown in Table 1), that is, the larger the speedup ratio, the maximum speedup ratio in the model example reaches 29167 times.

[0096] Table 1 Computational efficiency analysis of different model segmentation numbers

[0097]

[0098] Forward calculation accuracy analysis: The model of this experiment consists of a rectangular magnetic body (such as Figure 7 ), the spatial position and magnetic susceptibility of the magnetic body are shown in Table 2. The underground space is divided into 64×64×18, with a total of 73728 cubic grids. The side length of the cube is 2m, the magnetic field inclination is 60º, the magnetic declination is -9º, and the fluctuation range of the magnetic field observation surface is 3.0m~11.0m. Figure 8 The calculation results show that (as shown in (b) Figure 8 ), the fast algorithm has a high calculation accuracy, and the maximum absolute error is only about 1.0nT. According to the relative error calculation formula (8), the data relative error of the model is only 0.28%. The expression of formula (8) is:

[0099] (8)

[0100] in, represents the theoretical forward value of magnetic field; Indicates the forward modeling value of the magnetic field based on the fast algorithm; Figure 8 (a) and (b) show the calculation accuracy of the traditional algorithm for magnetic bodies. Figure 8 (c) and (d) show the fast algorithm calculation accuracy results of magnetic bodies.

[0101] Table 2 Parameters of magnetic body model

[0102]

[0103] Inversion calculation effect analysis: This test model consists of five rectangular magnetic bodies (such as Fig. 9 ), the spatial position and magnetic susceptibility of the magnetic body are shown in Table 3. The underground space is divided into 200×200×36, with a total of 1.44 million cubic grids. The side length of the cube is 2m, the magnetic field inclination is 60.0º, and the magnetic declination is -9.0º. The forward field of this model combination is shown in Fig.11 shown.

[0104] Table 3 Parameters of magnetic body model

[0105]

[0106] The constraints of 3D magnetic inversion are mainly depth weighted constraint, focus weighted constraint and physical property constraint. The physical property constraint sets the maximum magnetic susceptibility value to 0.03SI. Through multiple iterative calculations, the inversion results (such as Fig.12 ), and theoretical models (e.g. Fig. 9 ) comparison shows that the spatial positions and spatial forms of the five magnetic anomalies are well consistent. In addition, the inversion results are further sliced ​​horizontally and vertically (such as Fig.13 ), the slice results show that the magnetic susceptibility distribution of the three-dimensional inversion results is mainly located within the spatial range of the theoretical model (rectangular black frame), indicating that the overall consistency with the theoretical model is high, among which, Fig.13 (a) in the figure represents a horizontal slice of the model. Fig.13 (b) and (c) represent vertical slices of the model.

[0107] See also Fig.14 The present invention also provides a three-dimensional rapid forward and inversion system for undulating surface magnetic anomaly, which is applied to the above-mentioned three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly, comprising:

[0108] A measuring point determination module, used to obtain magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determine the magnetic field values ​​of the spatial observation grid points on the multiple horizontal observation surfaces according to the magnetic anomaly fields, wherein the spatial observation grid points are observation points on the multiple horizontal observation surfaces, the observation points include horizontal surface measurement points and undulating surface measurement points, the vertical positions of the spatial observation grid points are located on the multiple horizontal surfaces, and the horizontal positions of the spatial observation grid points coincide with the horizontal projection of the center of the underground subdivision grid;

[0109] A forward calculation module, used for performing forward calculation on the observation point according to the magnetic field value and a local extension interpolation algorithm to obtain a magnetic anomaly value of the measurement point on the undulating surface;

[0110] A reverse weighting module, used for reverse weighting the magnetic anomaly values ​​of the undulating surface measurement points to determine the grid point matrix of the spatial observation grid points;

[0111] The inversion calculation module is used to perform inversion calculation on the magnetic anomaly field based on the grid point matrix and the conjugate gradient algorithm to complete the rapid forward and inversion of the undulating surface magnetic anomaly.

[0112] In this embodiment, the process of improving the computational efficiency of three-dimensional inversion of undulating surface magnetic anomalies has the following characteristics: (1) Due to the local interpolation of this algorithm, the upward continuation method is adopted, which not only has practical physical significance, but also the interpolation results are stable and reliable, effectively ensuring the overall computational accuracy of this algorithm; (2) The traditional three-dimensional inversion of undulating surface magnetic anomalies requires large amounts of computational resources, while the algorithm proposed in this application greatly improves the computational efficiency of three-dimensional inversion of undulating surface magnetic anomalies by making full use of the advantages of the magnetic anomaly fast algorithm.

[0113] In all examples shown and described herein, any specific values ​​should be interpreted as merely exemplary and not as limiting, and thus other examples of the exemplary embodiments may have different values.

[0114] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, further definition and explanation thereof is not required in subsequent drawings.

[0115] The above-mentioned embodiments only express several implementation methods of the present invention, and the description thereof is relatively specific and detailed, but it cannot be understood as limiting the scope of the present invention. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly, characterized in that: The following steps are involved: Acquire magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determine magnetic field values ​​of multiple spatial observation grid points on the horizontal observation surfaces according to the magnetic anomaly fields, wherein the spatial observation grid points are observation points on multiple horizontal observation surfaces, the observation points include horizontal surface measurement points and undulating surface measurement points, the vertical positions of the spatial observation grid points are located on the multiple horizontal observation surfaces, and the horizontal positions of the spatial observation grid points coincide with the horizontal projections of the centers of the underground subdivided grids; Performing forward calculation on the observation point according to the magnetic field value and a local extension interpolation algorithm to obtain a magnetic anomaly value of the measurement point on the undulating surface; Inversely weighting the magnetic anomaly values ​​of the undulating surface measurement points to determine a grid point matrix of the spatial observation grid points; Based on the grid point matrix and the conjugate gradient algorithm, the magnetic anomaly field is inverted and calculated to complete the rapid forward and inversion of the undulating surface magnetic anomaly; The observation point is forward calculated according to the magnetic field value and the local extension interpolation algorithm to obtain the magnetic anomaly value of the measurement point on the undulating surface, including: The three-dimensional magnetic field fast forward algorithm is used to calculate the magnetic field of the horizontal plane measuring point to obtain the magnetic field forward value of the horizontal observation plane; Based on the magnetic field forward modeling value and the local continuation algorithm, the magnetic field of the measuring point on the undulating surface is calculated to obtain the magnetic anomaly value of the measuring point on the undulating surface; The details are as follows: The underground inversion space is gridded, and a plurality of horizontal magnetic field observation planes are set to obtain a horizontal observation plane, wherein the plurality of horizontal magnetic field observation planes are arranged at equal intervals according to a longitudinal length dz of the underground space grid; The magnetic field value of the horizontal observation surface is calculated by using a three-dimensional magnetic field fast algorithm; Calculate an interpolation weight matrix of a local extension of the observation point based on the undulating surface measurement point and the horizontal observation surface; Obtaining the magnetic field value of the measuring point on the undulating surface according to the magnetic field value of the horizontal observation surface and the interpolation weight matrix of the local extension; Multiple horizontal magnetic field observation surfaces are arranged at equal intervals according to the longitudinal length of the underground space, including: The maximum height of the horizontal observation plane is higher than the vertex of the undulating surface corresponding to the undulating surface measuring point minus one times the longitudinal length dz, and the minimum height of the horizontal observation plane is lower than the lowest point of the undulating surface minus one times the longitudinal length dz.

2. The three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly according to claim 1 is characterized in that: When the magnetic inclination and declination are constants, the magnetic field forward modeling formula can be simplified as: （5） Among them, V is the forward coefficient matrix for calculating the magnetic field values ​​of multiple horizontal observation surfaces; M is the model magnetic susceptibility matrix; Q represents the interpolation weight matrix of local extension.

3. The three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly according to claim 2 is characterized in that: The magnetic anomaly field is inverted and calculated based on the grid point matrix and the conjugate gradient algorithm to complete the rapid forward and inversion of the undulating surface magnetic anomaly, including: The solution to the inverse problem is obtained by minimizing the Tikhonov regularized objective function The corresponding expression is: （6） in, is the desired susceptibility model, is the regularization factor; is the fitting error objective function of the magnetic field; is the Tikhonov regularization objective function; D is the diagonal data covariance matrix, whose diagonal elements are the estimated data noise variance; is the reference model; W is the model weighting matrix; represents the measured magnetic field vector; T represents the transpose operation; According to the inversion calculation, the minimum value of the formula (6) corresponding to the objective function is obtained. According to the gradient of the minimum value of the objective function being 0, the fixed point iteration equation of the model M is obtained in the model space and data space. The corresponding calculation expression is as follows: （7） Among them, when formula (7) is solved using the conjugate gradient algorithm, its main calculation amount is the matrix The product with the vector P, P represents the intermediate vector; The calculation form of formula (7) includes and , column vector u and column vector v are intermediate variables and do not have fixed physical meanings.

4. The three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly according to claim 3 is characterized in that: Will The vector v of the term is regarded as the model magnetic susceptibility matrix M, and its calculation process is regarded as a forward calculation, which is achieved by reverse weighting Calculation of items, including: According to the interpolation weight matrix Q of the local extension of the measuring point, the value of the vector u is reversely weighted to multiple horizontal plane measuring points to obtain the pseudo magnetic field vector ; The pseudo magnetic field vector The pseudo model vector is calculated by a fast algorithm , thus completing Calculation of items.

5. The three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly according to claim 4 is characterized in that: The fast forward algorithm uses the FFT algorithm for acceleration, and the observation point in the FFT algorithm coincides with the center point of the underground subdivision grid in the horizontal projection.

6. A three-dimensional rapid forward and inversion system for undulating surface magnetic anomaly, applied to the three-dimensional rapid forward and inversion method for undulating surface magnetic anomaly according to any one of claims 1 to 5, characterized in that: include: A measuring point determination module, used to obtain magnetic anomaly fields of multiple horizontal observation surfaces covering the undulating surface, and determine the magnetic field values ​​of the spatial observation grid points on the multiple horizontal observation surfaces according to the magnetic anomaly fields, wherein the spatial observation grid points are observation points on the multiple horizontal observation surfaces, the observation points include horizontal surface measurement points and undulating surface measurement points, the vertical positions of the spatial observation grid points are located on the multiple horizontal observation surfaces, and the horizontal positions of the spatial observation grid points coincide with the horizontal projections of the centers of the underground subdivision grids; A forward calculation module, used for performing forward calculation on the observation point according to the magnetic field value and a local extension interpolation algorithm to obtain a magnetic anomaly value of the measurement point on the undulating surface; A reverse weighting module, used for reverse weighting the magnetic anomaly values ​​of the undulating surface measurement points to determine the grid point matrix of the spatial observation grid points; The inversion calculation module is used to perform inversion calculation on the magnetic anomaly field based on the grid point matrix and the conjugate gradient algorithm to complete the rapid forward and inversion of the undulating surface magnetic anomaly.

Images