A Fast Forward Modeling Method for Variable Dip Angle 3D Magnetic Fields

Abstract

The present invention discloses a fast forward modeling method for a variable dip angle three-dimensional magnetic field, specifically relating to the technical field of geophysical exploration. According to the positions of the observation points in the observation plane of the work area, the underground source space is divided into multiple cuboid models corresponding to the positions of the observation points, the forward calculation formula for the total magnetic field intensity vector of the cuboid model is determined, and based on the position and date of the work area, the magnetic parameter matrix at the center points of each cuboid model and the positions of each observation point is obtained from the International Geomagnetic Reference Field model. Combining with the forward calculation formula for the total magnetic field intensity vector of the cuboid model in the underground source space, a fast algorithm is used to calculate each magnetic field component of each cuboid model for determining the total magnetic field intensity vector, and the total magnetization intensity T of the underground source space in the work area is obtained. The method of the present invention improves the efficiency and accuracy of the fast forward modeling calculation of the variable dip angle three-dimensional magnetic field by integrating the magnetic parameters of each cuboid model and each observation point into the fast algorithm for gravity and magnetic fields, and is applicable to the fast forward modeling of large-area magnetic fields and strong remanent magnetic fields.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical exploration, and particularly relates to a fast forward modeling method for a variable dip angle three-dimensional magnetic field. Background Art

[0002] Magnetic exploration, as a commonly used geophysical exploration method, is mainly used to search for and explore minerals (such as iron ore, lead-zinc ore, copper-tin ore, etc.), conduct geological mapping, study geological structures related to oil and gas, study the earth's tectonics, and conduct military reconnaissance (such as unexploded ordnance, submarines, and mines, etc.). Inversion is an important process for interpreting magnetic targets, mainly to evaluate the contour and shape of unknown magnetic bodies underground or underwater by using the magnetic measurement data of the obtained magnetic targets.

[0003] As the basis of inversion, the efficiency and accuracy of forward modeling calculations directly affect the efficiency and effect of inversion calculations. The gravity and magnetic forward modeling algorithm based on the BTTB (Block-Toeplitz Toeplitz-Block) matrix has the advantages of high precision and high efficiency. Many scholars have carried out relevant research on it. Among them, the fast and high-precision three-dimensional magnetic field forward modeling algorithm based on the BTTB matrix proposed by Yuan Yang et al. can only calculate the case where the geomagnetic field direction and the magnetization intensity direction are fixed, and the forward modeling calculation is greatly affected by the magnetic field direction. For the case of a small area and the magnetic body only contains induced magnetization while ignoring the influence of remanent magnetization, when performing magnetic forward modeling, the magnetic parameters at the center of the work area are generally used to represent the magnetic parameters of the entire work area. Since the magnetic parameters in the work area change little, the forward modeling calculation can achieve a relatively high calculation accuracy. However, when the working area is large or there are strongly remanent magnetic bodies, the magnetic parameters in the work area vary greatly. If fixed magnetic parameters are forcedly used, it will lead to a large error in the forward modeling calculation results.

[0004] Therefore, there is an urgent need to propose a fast forward modeling method for a variable dip angle three-dimensional magnetic field to achieve high-precision and fast forward modeling of a large-area magnetic field or a strongly remanent magnetic field. Summary of the Invention

[0005] In order to solve the above technical problems, the present invention proposes a fast forward modeling method for a variable dip angle three-dimensional magnetic field, which realizes high-precision and fast forward modeling of a large-area magnetic field and a strongly remanent magnetic field, and is beneficial to accurately obtaining the true situation of the three-dimensional magnetic field in the underground space.

[0006] To achieve the above object, the present invention adopts the following technical scheme:

[0007] A fast forward modeling method for a variable dip angle three-dimensional magnetic field specifically includes the following steps:

[0008] Step 1: Obtain the working condition information of the work area, determine the positions of all observation points in the observation plane of the work area, divide the underground field source space of the work area into multiple cuboid models corresponding to the positions of each observation point, and determine the forward calculation formula for the total magnetic field intensity vector of the cuboid models in the underground field source space.

[0009] Step 2: Based on the position and date of the work area, calculate the magnetic parameter matrices at the center points of each cuboid model and the positions of each observation point based on the International Geomagnetic Reference Field model.

[0010] Step 3: Based on the forward calculation formula for the total magnetic field intensity vector of the cuboid models in the underground field source space of the work area, combined with the magnetic parameter matrices at the center points of each cuboid model and the positions of each observation point calculated by the International Geomagnetic Reference Field model, use a fast algorithm to calculate the total magnetic field intensity vector of each cuboid model for determination, and obtain the total magnetization intensity T of the underground field source space of the work area.

[0011] Preferably, in Step 1, obtain the working condition information of the work area, determine the total number P of observation points in the work area and the positions of each observation point, construct a three-dimensional space coordinate system with the ground surface as the observation plane, divide the underground field source space of the work area into p layers along the z direction, divide the underground field source space into m equal intervals along the x direction and n equal intervals along the y direction within each layer. After division, a total of N cuboid models are set in the underground field source space, and the total number N of cuboid models is m×n×p. The positions of each cuboid model in each layer in the horizontal direction correspond one-to-one with the positions of the ground surface observation points, and the total number P of observation points in the work area is m×n.

[0012] In the three-dimensional space coordinate system, obtain the coordinates of each observation point by dividing the grid on the observation plane, and determine that the forward calculation formula for the total magnetic field intensity vector of the cuboid model is a volume integral formula. The forward calculation formula for the total magnetic field intensity vector of the cuboid model is shown in Formula (1):

[0013]

[0014] Where,

[0015]

[0016]

[0017]

[0018]

[0019] r = [(ξ - x) + (η - y) + (ζ - z)] 1 / 2 (1 - 5)

[0020] Where, ΔT is the total magnetic field, (x, y, z) are the coordinates of the grid points on the observation plane; (ξ, η, ζ) are the coordinates of the field source points within the cuboid model; (x0, y0, z0) are the coordinates of the center point of the cuboid model, a is the extension length of the cuboid model in the x direction, b is the extension length of the cuboid model in the y direction, c is the extension length of the cuboid model in the z direction; M is the total magnetization of the cuboid model; μ0 is the magnetic permeability in vacuum; I0 is the dip angle of the total magnetization; D0 is the declination of the total magnetization; I1 is the dip angle of the geomagnetic field direction; D1 is the declination of the geomagnetic field direction; is the unit projection of the magnetization direction of the observation point in the x direction, is the unit projection of the magnetization direction of the observation point in the y direction, is the unit projection of the magnetization direction of the observation point in the z direction; V xx is the second-order derivative of the magnetic potential of the cuboid model in the x direction, V yy is the second-order derivative of the magnetic potential of the cuboid model in the y direction, V zz is the second-order derivative of the magnetic potential of the cuboid model in the z direction, V xy is the second-order derivative of the magnetic potential of the cuboid model in the xy direction, V yx is the second-order derivative of the magnetic potential of the cuboid model in the yx direction, V xz is the second-order derivative of the magnetic potential of the cuboid model in the xz direction, V zx is the second-order derivative of the magnetic potential of the cuboid model in the zx direction, V yz is the second-order derivative of the magnetic potential of the cuboid model in the yz direction, V zy is the second-order derivative of the magnetic potential of the cuboid model in the zy direction; T x is the component of the magnetic field in the x direction, T y is the component of the magnetic field in the y direction, T z is the component of the magnetic field in the z direction; π is the pi; M x is the component of the magnetic susceptibility of the cuboid model in the x direction, M y is the component of the magnetic susceptibility of the cuboid model in the y direction, M z is the component of the magnetic susceptibility of the cuboid model in the z direction; is the unit projection of the magnetization direction of the cuboid model in the x direction, is the unit projection of the magnetization direction of the cuboid model in the y direction, is the unit projection of the magnetization direction of the cuboid model in the z direction; r is the distance between the field source point within the cuboid and the observation point.

[0021] Preferably, in the step 2, by inputting the position at the center point of the cuboid model and the positions of the observation points into the International Geomagnetic Reference Field Model, the magnetic parameter matrix at the center point of each cuboid model in the underground field source space is calculated using the International Geomagnetic Reference Field Model and the magnetic parameter matrix at the positions of the observation points

[0022] Preferably, in the step 3, when the cuboid model and the observation points in the underground field source space have different inclination angles and declination angles, the total magnetic field intensity vector of the cuboid model is expressed as:

[0023]

[0024] In the formula, is a diagonal matrix with the matrix size of P×P, including the diagonal matrix diagonal matrix and the diagonal matrix is a two-dimensional matrix with the matrix size of N×N, including the diagonal matrix diagonal matrix and the diagonal matrix m is the magnetic susceptibility model matrix with the matrix size of N×1; the magnetic susceptibility forward modeling coefficient matrix V rs , with the matrix size of P×N, including the kernel matrix V xx 、kernel matrix V xy 、kernel matrix V xz 、kernel matrix V yx 、kernel matrix V yy 、kernel matrix V yz 、kernel matrix V zx 、kernel matrix V zy and the kernel matrix V zz ;

[0025] When the cuboid model and the observation points in the underground field source space have the same inclination angle and declination angle, combining the identity matrix I d and I m , the total magnetic field intensity vector of the cuboid model is expressed as:

[0026] ΔT = Vm (3)

[0027] Among them,

[0028]

[0029] In the formula, q x 、q y 、q z are all scalars related to the inclination angle.

[0030] Preferably, in the step 3, the magnetic susceptibility forward modeling coefficient matrix Vrs Each core matrix in is a BTTB matrix with a large number of repeated matrix elements. Based on a fast algorithm, the total magnetic field intensity vector of the cuboid model is calculated. Using the three-dimensional matrix rs where r, s ∈ {x, y, z}, the size of the three-dimensional matrix is (2m - 1) × (2n - 1) × p. Then, using the three-dimensional matrix with size m × n × p to replace the vector m, the three-dimensional matrix with size m × n × p to replace the two-dimensional matrix the three-dimensional matrix with size m × n × p to replace the two-dimensional matrix After that, the total magnetic field intensity vector of each cuboid model in the underground field source space is calculated, which specifically includes the following steps:

[0031] Step 3.1, when the cuboid model in the underground field source space has different dip angles and declination angles from the observation point, is equivalent to By multiplying the corresponding elements in the same-order matrices in the total magnetic field intensity vector calculation formula of the cuboid model, we get:

[0032]

[0033] In the formula, m t is the calculated magnetic susceptibility component, t ∈ {x, y, z}, and calculating the magnetic susceptibility component m t includes the matrix m x , the matrix m y and the matrix m z ;

[0034] Step 3.2, substitute T rs = V rs m t in formula (5) with where r, s, t ∈ {x, y, z}, k is the layer number, * is the convolution calculation operator. Based on the fast algorithm, convolution calculation is performed through the forward and inverse fast Fourier transforms to quickly determine and we get:

[0035]

[0036] In the formula, F is the forward fast Fourier transform, and F -1 is the inverse fast Fourier transform;

[0037] Step 3.3, substitute in formula (6) with End the rapid forward modeling of the variable dip angle three-dimensional magnetic field in the underground source space of the work area, and determine the total magnetic field intensity T of the underground source space in the work area.

[0038] The beneficial technical effects brought by the present invention:

[0039] The present invention proposes a rapid forward modeling method for variable dip angle three-dimensional magnetic fields, which solves the problem that the existing forward modeling method for variable dip angle three-dimensional magnetic fields has low calculation efficiency, resulting in difficulty in carrying out rapid calculations of large-scale data. The method of the present invention divides the underground source space into multiple cuboid models, separately sets magnetic parameters for each cuboid model and observation points, and improves the calculation efficiency and calculation accuracy of the rapid forward modeling of variable dip angle three-dimensional magnetic fields by quickly and accurately obtaining the total magnetic field intensity of each cuboid model, laying a foundation for the detection and evaluation of unknown magnetic bodies underground or underwater. Description of the drawings

[0040] Figure 1 It is a schematic diagram of the cuboid model in this embodiment.

[0041] Figure 2 It is a schematic diagram of the corresponding relationship between the cuboid model and the observation points in this embodiment.

[0042] Figure 3 It is a schematic diagram of the magnetic body model in this embodiment.

[0043] Figure 4 It is the magnetic parameter distribution map of the underground space in this embodiment; Figure 4 In it, (a) is the magnetic dip angle distribution map, and (b) is the magnetic declination distribution map.

[0044] Figure 5 It is the calculation result of the total magnetic field intensity vector of the cuboid model in this embodiment; Figure 5 In it, (a) is the forward modeling result of ΔT of the cuboid model, and (b) is the theoretical polarization result of the cuboid model. Specific implementation manners

[0045] The present invention will be further described in detail below in conjunction with the drawings and specific implementation manners:

[0046] Embodiment 1

[0047] The present invention proposes a rapid forward modeling method for variable dip angle three-dimensional magnetic fields, which specifically includes the following steps:

[0048] Step 1, obtain the working condition information of the work area, determine the positions of all observation points in the observation plane of the work area, divide the underground source space of the work area into multiple cuboid models corresponding to the positions of each observation point, and determine the forward calculation formula of the total magnetic field intensity vector of the cuboid models in the underground source space.

[0049] In this embodiment, the working condition information of the work area is obtained, the total number P of observation points in the work area and the positions of each observation point are determined, a three-dimensional space coordinate system is constructed with the ground surface as the observation plane, the underground field source space in the work area is divided into p layers along the z direction, and the underground field source space in each layer is equally divided into m parts along the x direction and n parts along the y direction at equal intervals, so as to evenly divide the underground field source space into a plurality of cuboid models, as Figure 1 shown. After division, a total of N cuboid models are set in the underground field source space, and the total number N of cuboid models is m×n×p.

[0050] Since both the underground field source space and the observation points are in discrete form during actual calculation, in order to ensure the resolution of the forward calculation using the cuboid model, the positions of each cuboid model in each layer in the horizontal direction correspond one by one to the positions of the surface observation points, as Figure 2 shown, that is, the total number P of observation points in the work area is m×n.

[0051] In the three-dimensional space coordinate system, the coordinates of each observation point are obtained by dividing the grid on the observation plane, and the forward calculation formula for the total magnetic field intensity vector of the cuboid model is the volume integral formula. The forward calculation formula for the total magnetic field intensity vector of the cuboid model is as shown in formula (1):

[0052]

[0053] Among them,

[0054]

[0055]

[0056]

[0057]

[0058] r = [(ξ - x)+(η - y)+(ζ - z)] 1 / 2 (1-5)

[0059] In the formula, ΔT is the total magnetic field, (x, y, z) are the coordinates of the grid points on the observation plane; (ξ, η, ζ) are the coordinates of the field source points in the cuboid model; (x0, y0, z0) are the coordinates of the center point of the cuboid model, a is the extension length of the cuboid model along the x direction, b is the extension length of the cuboid model along the y direction, c is the extension length of the cuboid model along the z direction; M is the total magnetization intensity of the cuboid model; μ0 is the magnetic permeability in vacuum; I0 is the dip angle of the total magnetization intensity; D0 is the declination of the total magnetization intensity; I1 is the dip angle of the geomagnetic field direction; D1 is the declination of the geomagnetic field direction; is the unit projection of the magnetization direction of the observation point in the x direction, is the unit projection of the magnetization direction of the observation point in the y direction, is the unit projection of the magnetization direction of the observation point in the z direction; V xx is the second-order derivative of the magnetic potential of the cuboid model in the x direction, V yy is the second-order derivative of the magnetic potential of the cuboid model in the y direction, V zz is the second-order derivative of the magnetic potential of the cuboid model in the z direction, V xy is the second-order derivative of the magnetic potential of the cuboid model in the xy direction, V yx is the second-order derivative of the magnetic potential of the cuboid model in the yx direction, V xz is the second-order derivative of the magnetic potential of the cuboid model in the xz direction, V zx is the second-order derivative of the magnetic potential of the cuboid model in the zx direction, V yz is the second-order derivative of the magnetic potential of the cuboid model in the yz direction, V zy is the second-order derivative of the magnetic potential of the cuboid model in the zy direction; T x is the component of the magnetic field in the x direction, T y is the component of the magnetic field in the y direction, T z is the component of the magnetic field in the z direction; π is the circumference ratio; M x is the component of the magnetic susceptibility of the cuboid model in the x direction, M y is the component of the magnetic susceptibility of the cuboid model in the y direction, M z is the component of the magnetic susceptibility of the cuboid model in the z direction; is the unit projection of the magnetization direction of the cuboid model in the x direction, is the unit projection of the magnetization direction of the cuboid model in the y direction, is the unit projection of the magnetization direction of the cuboid model in the z direction; r is the distance between the field source point and the observation point inside the cuboid.

[0060] Step 2: According to the location and date of the work area, calculate the magnetic parameter matrix at the center point of each cuboid model and at the position of each observation point based on the International Geomagnetic Reference Field model.

[0061] In this embodiment, the International Geomagnetic Reference Field model is the prior art for those skilled in the art. By inputting the position of the center point of each cuboid model and the positions of each observation point into the International Geomagnetic Reference Field model, the magnetic parameter matrix at the center point of each cuboid model in the underground field source space and the magnetic parameter matrix at the position of each observation point

[0062] Step 3: Based on the forward calculation formula of the total magnetic field intensity vector of the cuboid model in the subsurface source space of the work area, combined with the magnetic parameter matrix at the center points of each cuboid model and the positions of each observation point calculated by the International Geomagnetic Reference Field model, use the fast algorithm to calculate the total magnetic field intensity vector of each cuboid model for determining the total magnetic field intensity vector, and obtain the total magnetization intensity T of the subsurface source space of the work area.

[0063] When the cuboid model in the subsurface source space has different dip angles and declination angles from the observation point, the total magnetic field intensity vector of the cuboid model is expressed as:

[0064]

[0065] In the formula, is a diagonal matrix with matrix size P×P, including the diagonal matrix diagonal matrix and diagonal matrix is a two-dimensional matrix with matrix size N×N, including the diagonal matrix diagonal matrix and diagonal matrix m is the magnetic susceptibility model matrix with matrix size N×1; the magnetic susceptibility forward coefficient matrix V rs , with matrix size P×N, includes the kernel matrix V xx 、kernel matrix V xy 、kernel matrix V xz 、kernel matrix V yx 、kernel matrix V yy 、kernel matrix V yz 、kernel matrix V zx 、kernel matrix V zy and kernel matrix V zz .

[0066] When the cuboid model in the subsurface source space has the same dip angle and declination angle as the observation point, combined with the identity matrix I d and I m , the total magnetic field intensity vector of the cuboid model is expressed as:

[0067] ΔT = Vm (3)

[0068] where,

[0069]

[0070] In the formula, q x 、q y 、q z are all scalars related to the dip angle.

[0071] In the formula (2), the magnetic susceptibility forward coefficient matrix V rsEach core matrix in it is a BTTB matrix, also known as a block Toeplitz matrix. There are a large number of repeated elements in such matrices. Therefore, only some elements in the matrix need to be calculated to obtain all the elements contained in the matrix, thus significantly shortening the time for calculating the core matrix in the fast forward modeling of three-dimensional magnetic fields with variable dip angles.

[0072] Calculate the total magnetic field intensity vector of the cuboid model based on the fast algorithm, and use the three-dimensional matrix to replace the susceptibility forward coefficient matrix V rs , where r, s ∈ {x, y, z}, and the three-dimensional matrix has a size of (2m - 1) × (2n - 1) × p. Then, use the three-dimensional matrix with a size of m × n × p to replace the vector m, and use the three-dimensional matrix with a size of m × n × p to replace the two-dimensional matrix The three-dimensional matrix with a size of m × n × p is used to replace the two-dimensional matrix Using small three-dimensional or two-dimensional matrices to represent large two-dimensional matrices or vectors significantly reduces the memory requirements during the forward modeling calculation. Calculate the total magnetic field intensity vector of each cuboid model in the underground field source space, which specifically includes the following steps:

[0073] Step 3.1, when the cuboid model in the underground field source space has different dip angles and declination angles from the observation point, make equivalent to By multiplying the corresponding elements in the matrices of the same order in the total magnetic field intensity vector calculation formula of the cuboid model, we get:

[0074]

[0075] In the formula, m t is the calculated susceptibility component, t ∈ {x, y, z}, and calculating the susceptibility component m t includes the matrix m x , the matrix m y and the matrix m z .

[0076] Step 3.2, make T rs = V rs m t equivalent to where r, s, t ∈ {x, y, z}, k is the layer number, and * is the convolution calculation operator. Based on the fast algorithm, perform the convolution calculation through the forward and inverse fast Fourier transforms to quickly determine and obtain:

[0077]

[0078] Where F is the fast Fourier transform, F -1 is the inverse fast Fourier transform;

[0079] Step 3.3: Substitute Equivalent to The variable-angle three-dimensional magnetic field rapid forward modeling of the underground field source space in the work area is completed to determine the total magnetic field intensity T of the underground field source space in the work area.

[0080] In the process of three-dimensional forward modeling of variable-angle magnetic field, the method of the present invention does not need to explicitly generate a kernel matrix, which greatly saves memory space. At the same time, due to the use of fast Fourier transform, the computational efficiency of the three-dimensional forward modeling of variable-angle magnetic field is greatly improved, making the computational efficiency comparable to that of frequency domain calculation.

[0081] Example 2

[0082] In order to verify the effectiveness of the fast forward modeling method of the variable inclination three-dimensional magnetic field proposed by the present invention, three rectangular magnetic bodies were constructed in this embodiment, such as Figure 3 As shown, the fast forward modeling method of the three-dimensional magnetic field with variable inclination angle described in Example 1 is used to perform forward calculations.

[0083] In this embodiment, three rectangular-shaped magnetic bodies are set in the underground space. The spatial position and magnetic susceptibility of each magnetic body are shown in Table 1. A three-dimensional space coordinate system is established in the underground space and the underground space is divided. The underground space is divided into 150 layers along the z direction. In each layer, the underground field source space is divided into 300m parts along the x direction and 600 parts along the y direction. After the division, a total of 27 million rectangular models are set in the underground field source space. In this embodiment, the magnetic field inclination angle of the underground space varies from -20° to 39.9°, and the magnetic declination angle varies from 0° to 11.98°. Figure 4 shown.

[0084] Table 1 Spatial position and magnetic susceptibility of various magnetic bodies in underground space

[0085]

[0086] In this embodiment, a computer with a main frequency of 2.20 GHz and a memory of 64 GB is used for forward simulation calculation to test the feasibility of the method of the present invention. It takes about 297 seconds to calculate the kernel matrix, and about 52 seconds to calculate the forward calculation based on the kernel matrix. The ΔT forward modeling results and theoretical extreme results of the rectangular parallelepiped model are shown in Figure 2. Figure 5 As shown ( Figure 5 Shown are the calculation results of the total magnetic field intensity vector of the rectangular model; in the figure, (a) is the ΔT forward result of the rectangular model, and (b) is the theoretical polar result of the rectangular model).

[0087] It can be seen from Figure 5 that the characteristics of the theoretical pole field have a good correspondence with the magnetic body. Since the magnetic dip angle of the theoretical pole field is 90°, and the total magnetic field changes with the magnetic dip angle and magnetic declination, the magnetic field characteristics show obvious differences. The dip directions of magnetic body 1 and magnetic body 2 are symmetrically opposite with respect to the equator (magnetic dip angle is zero), so the forward fields of the two are basically symmetric, and magnetic body 3 shows the magnetic anomaly characteristics of a lower latitude.

[0088] Thus, it is verified that the forward modeling of the present invention is based on a high-precision fast algorithm, realizing the fast forward modeling of the variable dip three-dimensional magnetic field, improving the calculation efficiency and calculation accuracy of the variable dip three-dimensional magnetic field forward modeling, being applicable to large-area magnetic fields and strong remanent magnetic fields, and being beneficial to the detection and evaluation of unknown magnetic bodies underground or underwater.

[0089] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the essence of the present invention should also fall within the protection scope of the present invention.

Claims

1. A fast forward modeling method for variable dip angle three-dimensional magnetic fields, characterized in that, The specific steps include: Step 1, obtaining the working condition information of the work area, determining the positions of all observation points in the observation plane of the work area, dividing the underground field source space of the work area into multiple rectangular parallelepiped models corresponding to the positions of the observation points, and determining the forward calculation formula of the total magnetic field intensity vector of the rectangular parallelepiped model in the underground field source space; Step 2: Based on the location and date of the work area, calculate the magnetic parameter matrix at the center points of each cuboid model and at the positions of each observation point based on the International Geomagnetic Reference Field model Step 3, based on the forward calculation formula of the total magnetic field intensity vector of the rectangular model in the underground field source space of the work area, combined with the magnetic parameter matrix at the center point of each rectangular model and each observation point calculated by the International Geomagnetic Reference Field Model, use a fast algorithm to calculate each rectangular model for determining the total magnetic field intensity vector, and obtain the total magnetization intensity T of the underground field source space of the work area; In step 2, by inputting the positions at the center points of the cuboid models and the positions of the observation points into the International Geomagnetic Reference Field (IGRF) model, the magnetic parameter matrices at the center points of each cuboid model in the underground field source space are calculated using the IGRF model and the magnetic parameter matrices at the positions of the observation points In step 3, when the rectangular parallelepiped model and the observation point in the ground field source space have different inclinations and declinations, the total magnetic field intensity vector of the rectangular parallelepiped model is expressed as: In the formula, is a diagonal matrix with the matrix size of P×P, including diagonal matrix diagonal matrix and diagonal matrix is a two-dimensional matrix with the matrix size of N×N, including diagonal matrix diagonal matrix and diagonal matrix m is the magnetic susceptibility model matrix with the matrix size of N×1; the forward magnetic susceptibility coefficient matrix V rs , with the matrix size of P×N, includes kernel matrix V xx , kernel matrix V xy , kernel matrix V xz , kernel matrix V yx , kernel matrix V yy , kernel matrix V yz , kernel matrix V zx , kernel matrix V zy and kernel matrix V zz ; When the cuboid model and the observation point in the local subsurface source space have the same dip angle and declination angle, combined with the identity matrix I d and I m , the total magnetic field intensity vector of the cuboid model is expressed as: ΔT=Vm (3) in, where q x , q y , q z are all scalars related to the dip angle.

2. The fast forward modeling method for variable dip angle three-dimensional magnetic field according to claim 1, wherein In the step 1, the working condition information of the work area is obtained, the total number P of observation points in the work area and the position of each observation point are determined, a three-dimensional space coordinate system is constructed with the surface as the observation plane, the underground field source space of the work area is divided into p layers along the z direction, and the underground field source space is divided into m parts with equal intervals along the x direction and n parts with equal intervals along the y direction in each layer. After the division, a total of N rectangular parallelepiped models are arranged in the underground field source space, and the total number N of rectangular parallelepiped models is m×n×p. The position of each rectangular parallelepiped model in each layer in the horizontal direction corresponds one-to-one to the position of the surface observation point, and the total number P of observation points in the work area is m×n; In the three-dimensional space coordinate system, the coordinates of each observation point are obtained by dividing the grid on the observation plane, and the forward calculation formula of the total magnetic field intensity vector of the rectangular model is determined to be a volume integral formula. The forward calculation formula of the total magnetic field intensity vector of the rectangular model is shown in formula (1): in, r = [(ξ - x) + (η - y) + (ζ - z)] 1 / 2 (1 - 5) Where, ΔT is the total magnetic field, (x, y, z) are the coordinates of the grid points on the observation plane; (ξ, η, ζ) are the coordinates of the field source points inside the cuboid model; (x0, y0, z0) are the coordinates of the center point of the cuboid model, a is the extension length of the cuboid model in the x direction, b is the extension length of the cuboid model in the y direction, c is the extension length of the cuboid model in the z direction; M is the total magnetization intensity of the cuboid model; μ0 is the magnetic permeability in vacuum; I0 is the dip angle of the total magnetization intensity; D0 is the declination of the total magnetization intensity; I1 is the dip angle of the geomagnetic field direction; D1 is the declination of the geomagnetic field direction; is the unit projection of the magnetization direction of the observation point in the x direction, is the unit projection of the magnetization direction of the observation point in the y direction, is the unit projection of the magnetization direction of the observation point in the z direction; V xx is the second-order derivative of the magnetic force potential of the cuboid model in the x direction, V yy is the second-order derivative of the magnetic force potential of the cuboid model in the y direction, V zz is the second-order derivative of the magnetic force potential of the cuboid model in the z direction, V xy is the second-order derivative of the magnetic force potential of the cuboid model in the xy direction, V yx is the second-order derivative of the magnetic force potential of the cuboid model in the yx direction, V xz is the second-order derivative of the magnetic force potential of the cuboid model in the xz direction, V zx is the second-order derivative of the magnetic force potential of the cuboid model in the zx direction, V yz is the second-order derivative of the magnetic force potential of the cuboid model in the yz direction, V zy is the second-order derivative of the magnetic force potential of the cuboid model in the zy direction; T x is the component of the magnetic field in the x direction, T y is the component of the magnetic field in the y direction, T z is the component of the magnetic field in the z direction; π is the pi; M x is the component of the magnetic susceptibility of the cuboid model in the x direction, M y is the component of the magnetic susceptibility of the cuboid model in the y direction, M z is the component of the magnetic susceptibility of the cuboid model in the z direction; is the unit projection of the magnetization direction of the cuboid model in the x direction, is the unit projection of the magnetization direction of the cuboid model in the y direction, is the unit projection of the magnetization direction of the cuboid model in the z direction; r is the distance between the field source point inside the cuboid and the observation point.

3. The fast forward modeling method for variable dip angle three-dimensional magnetic field according to claim 1, characterized in that In step 3, the susceptibility forward coefficient matrix V rs each kernel matrix in is a BTTB matrix, with a large number of repeated matrix elements. Based on a fast algorithm, the total magnetic field intensity vector of the cuboid model is calculated, and a three-dimensional matrix is used to replace the susceptibility forward coefficient matrix V rs , where r, s ∈ {x, y, z}, and the size of the three-dimensional matrix is (2m - 1) × (2n - 1) × p. Then, a three-dimensional matrix with a size of m × n × p is used to replace the vector m, a three-dimensional matrix with a size of m × n × p is used to replace the two-dimensional matrix After a three-dimensional matrix with a size of m × n × p is used to replace the two-dimensional matrix , the total magnetic field intensity vector of each cuboid model in the underground field source space is calculated, which specifically includes the following steps: Step 3.1, when the cuboid model and the observation point in the local subsurface source space have different dip angles and declination angles, equivalent to by multiplying the corresponding elements in the same-order matrices in the total magnetic field intensity vector calculation formula of the cuboid model, we get: where m t is for calculating the magnetic susceptibility component, t ∈ {x, y, z}, and calculating the magnetic susceptibility component m t includes matrix m x , matrix m y and matrix m z ; Step 3.2, substitute T in formula (5) with rs = V rs m t which is equivalent to where r, s, t ∈ {x, y, z}, k is the layer number, * is the convolution calculation operator, and based on the fast algorithm, convolution calculation is performed through the forward and inverse fast Fourier transforms to quickly determine to obtain: where F is the forward fast Fourier transform, and F -1 is the inverse fast Fourier transform; Step 3.3, substitute the equivalent to End the fast 3D magnetic field forward modeling with variable dip angles for the underground source space in the work area, and determine the total magnetic field intensity T of the underground source space in the work area.

Images