A spectral-spatial domain iterative method for constructing spatial geomagnetic maps

Abstract

The present invention relates to a spectral domain-spatial domain iterative spatial geomagnetic map construction method, by constructing an upward extended geomagnetic map; calculating high-order derivative items based on robust upward extension according to the upward extended geomagnetic map; performing downward iterative extension processing in the spatial domain according to the calculation results of the high-order derivative items; and obtaining a vertical high-precision geomagnetic map by iteratively correcting and processing the downward iterative extension processing results in the spatial domain. The present invention uses the vertical spatial distribution of geomagnetism at each point on the geomagnetic reference plane as a function model, and the vertical geomagnetic space domain high-order derivatives as the parameters to be solved for fitting, converting the downward extended geomagnetic construction into the fitting and extension based on the geomagnetic reference map and the upward extended geomagnetic map, and performing high-order derivative truncation through iterative processing. Therefore, this method can fundamentally avoid the problem of solving high-dimensional ill-conditioned matrices in the downward extension process, thereby obtaining a robust downward extension effect.

Description

Technical Field

[0001] The present invention belongs to the technical field of space geomagnetic maps, and in particular to a spectrum domain-space domain iterative space geomagnetic map construction method. Background Art

[0002] The theory of geomagnetic field continuation mainly involves two parts: upward continuation and downward continuation. The upward continuation is a boundary value problem of the Laplace equation, which can obtain a stable solution; the downward continuation is the inverse problem of the upward continuation and is an ill-posed problem; generally starting from the integral equation of the upward continuation, the original downward continuation operator is derived in the frequency domain and solved by the approximate approximation method. According to the requirements of the technical specifications for aeromagnetic survey, the general aeromagnetic survey operation height is generally 600m-800m, but the aeromagnetic exploration operation height is 100m-200m. Therefore, it is urgent to extend the aerial survey benchmark geomagnetic map downward to the aeromagnetic exploration height.

[0003] According to classical field theory, in a passive space, the upward extension of the field position can be expressed as a boundary value problem of the Laplace equation, as shown in equation (1), and its schematic diagram is shown in Figure 1 shown.

[0004]

[0005] Among them, Γ S , Γ Up , Γ Down Respectively represent the actual observation surface, the upward extension surface and the downward extension surface, Γ ∞ represents the infinite boundary, represents the distant area extending outward along the observation plane, Γ S , Γ far , Γ ∞ Together they form a closed surface, the enclosed area is Ω, and the plane to be extended is within Ω. Since the potential field decays with distance, the infinite boundary Γ ∞ The magnetic value on the surface is generally 0, which is called the natural boundary condition; the prospect area Γ far The contribution to the main area extending upward and downward is very small, and generally only affects the boundary area, which can be regarded as an infinite boundary.

[0006] From formula (1), we can see that once the magnetic field value at the boundary is determined, the magnetic field value at any spatial position in the region Ω can be uniquely determined. In practice, the magnetic field at infinity has little effect on the upward extension area, so the boundary conditions can be simplified. In the processing process, the boundary condition that the boundary at infinity is zero can be ignored, and only the observation surface Γ is retained. S As the only boundary condition. Therefore, the problem of extending the geomagnetic field upward is actually to use the geomagnetic field of the observation plane as the boundary condition, solve the partial differential equation (2) with the help of Green's second formula, calculate the geomagnetic field value in the potential field Ω, and thus obtain the integral formula for upward extension.

[0007]

[0008] Where h is the extension height, and h>0, and u(x,y,-h) is the plane to be extended Γ Up Magnetic field data, u(x,y,0) is the observation plane Γ S From the above formula, it can be seen that the integral equation of the geomagnetic field can be extended upward to obtain a stable solution.

[0009] However, when solving the downward continuation problem, the downward continuation surface is below the known measurement plane (such as Figure 1 (shown) is an extension in the direction of the field source, which cannot meet the condition that the downward extension space is passive. The prior art points out that the method of describing the downward extension problem in mathematical language is to regard it as the inverse problem of the upward extension problem, that is, the upward extension is to calculate the regional potential field value through the boundary potential field observation value, while the downward extension is to calculate the boundary potential field value through the regional potential field observation value, so it is an ill-posed problem.

[0010] like Figure 1 As shown, by establishing a Cartesian coordinate system, the Z axis is downward as the positive direction. Assuming that the observation plane is known to be Z = 0, the observation plane u(x, y, 0) is extended downward to a height of h, and the potential field value u(x, y, h) of the downward extended surface to be calculated can be obtained as:

[0011]

[0012] By performing spatial Fourier transform on the above equation, it can be converted to the spectral domain for solution, that is,

[0013]

[0014] It can be seen from the above formula that the potential field value U of the downward extension surface h (k x ,k y ) can be solved by discretization into a linear system of equations:

[0015] Y=Aθ(5)

[0016] Where Y = U0(k x ,k y ), θ=U h (k x ,k y ). x ,k y is the wave number at the discretized grid points.

[0017] The downward continuation problem is transformed into the problem of solving θ through equation (5). In theory, θ can be solved by least squares estimation, but if A is a high-dimensional ill-conditioned matrix, the solution of θ is not feasible in practice. For the downward continuation problem, since the discretized grid is high-dimensional, the matrix A is also high-dimensional. At the same time, the ill-conditioning of A can be determined by the condition number. That is,

[0018]

[0019] Among them, ·2 represents the 2-norm of the matrix, σ max is the maximum eigenvalue of the matrix, σ min is the smallest eigenvalue of the matrix. In addition, the downward extension operator It is equivalent to a high-pass filter, which causes the high-frequency signal components such as target information and noise in the geomagnetic field on the observation surface to diverge, resulting in instability in the downward extension process.

[0020] From the above analysis, we can see that the downward continuation problem belongs to the problem of solving high-dimensional ill-conditioned matrix equations. In order to solve high-dimensional ill-conditioned linear equations, the traditional downward continuation generally adopts the Landweber iteration method, that is, gradually approximating the optimal solution of the objective function through iteration, reducing the solution error caused by the inversion of the ill-conditioned matrix. First, the solution problem of (5) is defined as an optimization problem:

[0021] θ=argminJ(θ),

[0022]

[0023] Then the optimal solution is found through iteration. The iteration format of Landweber iteration method is:

[0024] θ n+1 =θ n +sA T (Y-Aθ n )(8)

[0025] Here, s is the step size.

[0026] The Landweber iteration method can effectively improve the instability problem of solving high-dimensional ill-conditioned matrix equations. It can be understood as "connecting in series" a low-pass filter that changes with the number of iterations on the basis of the original downward extension. However, as the number of iterations increases, the low-pass frequency band shifts to the right and becomes an all-pass filter, thereby losing the low-pass effect and causing divergence in the later stages of iteration. Therefore, the convergence of the iterative method is greatly affected by the number of iterations. On the one hand, it may require an infinite number of iterations to converge; on the other hand, there may be a semi-convergence phenomenon, that is, the stable solution can be approached at the beginning of the iteration, but as the number of iterations increases, the iterative solution tends to diverge. This limitation is mainly because the iterative method only uses "successive approximation" to solve, avoiding direct matrix inversion. This method reduces the influence of matrix ill-conditioning, but does not fundamentally overcome the ill-conditioning problem of the matrix equation to be solved (Equation (5)). Summary of the invention

[0027] The purpose of the present invention is to overcome the shortcomings of the prior art and propose a spectral domain-spatial domain iterative spatial geomagnetic map construction method. Based on the spectral domain results of robust upward extension, the first law of geography, that is, the smoothness of the spatial geomagnetic field distribution, is utilized to fit the geomagnetic spatial distribution through high-order derivative terms in the spatial domain to achieve downward extension of the geomagnetic map.

[0028] The present invention solves the technical problem by adopting the following technical solutions:

[0029] A spectral domain-spatial domain iterative spatial geomagnetic map construction method comprises the following steps:

[0030] Step 1, constructing an upward extended geomagnetic map;

[0031] Step 2, calculating high-order derivatives based on robust upward continuation according to the upward extended geomagnetic map;

[0032] Step 3: Perform downward iterative extension processing in the spatial domain according to the calculation results of the high-order derivative terms;

[0033] Step 4: Iteratively extend the processing results downward in the spatial domain through iterative correction processing to obtain a vertical high-precision geomagnetic map.

[0034] Moreover, the specific implementation method of step 1 is: the integral formula for upward extension is:

[0035]

[0036] Where h is the extension height, and h>0, and u(x,y,-h) is the plane to be extended Γ Up Magnetic field data, u(x,y,0) is the observation plane Γ S Magnetic field data, where ξ is the dimension of the data observed along the x direction, and η is the dimension of the data observed along the y direction;

[0037] Performing Fourier transform on the upward extension equation, we get its spectral domain expression as follows:

[0038]

[0039] because

[0040]

[0041] in k x , k y Respectively represent the spectral domain variables corresponding to the spatial domain variables x and y:

[0042]

[0043] Among them, L x , L y Represent the fundamental wavelength in the x and y directions respectively:

[0044] L x =MΔx,L y =NΔy

[0045] Wherein, Δx and Δy represent the point distance and line distance of the observation data respectively, the number of sampling points along the x direction is M, and the number of sampling points along the y direction is N;

[0046] When M and N are odd numbers,

[0047]

[0048] When M and N are even numbers,

[0049]

[0050] The spectral domain calculation formula for the extension of the geomagnetic map plane u0 upward to height h is:

[0051] U -h (u,v)=F[u0(x,y)]e -2πfh

[0052] Perform inverse Fourier transform on it and calculate the geomagnetic map of the upward extension plane:

[0053] u -h (x,y)=F -1 {F[u0(x,y)]e -2πfh}.

[0054] Moreover, the specific implementation method of step 2 is: according to the upward continuation equation, the following is obtained:

[0055]

[0056] Among them, the p plane is the measurement plane, the o plane is the flight plane of the magnetic detection aircraft, the s plane is the sea level, and q1~q n are n planes extended upward based on the p plane. The coordinate system is located on the p plane, downward is positive, and the height from the p plane to the o plane is h, and the height from the p plane to q1~q n The heights of the surfaces are -△z, -2△z, -3△z, ..., -n△z, and the height from the p surface to the s surface is l;

[0057] Combined with the geomagnetic map constructed by the known plane p, we can get different levels q1~q n Magnetic map of:

[0058]

[0059] Different levels q1~q n The geomagnetic map of the known reference layer can be extended upward to a geomagnetic map with a height of nΔz. n Based on the plane geomagnetic map, we further analyze the different levels q1~q at the plane single coordinate point. n The geomagnetic map is expanded by Taylor series in the space domain:

[0060]

[0061] in, is the plane geomagnetic map obtained by upward extension, n represents the number of upward extension layers, u p is the plane geomagnetic map of the observation plane, Δz is the height difference of the upward extension plane; Taylor series expansion describes the vertical space geomagnetic field distribution of a single coordinate point on the plane at different heights. For a single coordinate point on the observation plane p and the upward extension plane q1...q n For a single coordinate point above, a magnetic function with respect to height z is formed in the vertical direction.

[0062] Get the different heights of the first point q1...q n The equations for geomagnetic data are:

[0063]

[0064] Writing the system of equations in matrix form:

[0065] Y=AX

[0066] mmm

[0067] in,

[0068] Calculated using the least squares method:

[0069] X m =(Am T A m ) -1 A m T Y m

[0070] Moreover, the specific implementation method of step 3 is:

[0071] At the downward extension plane o and the upward extension plane q, respectively:

[0072]

[0073]

[0074] Where Δz po =z p -z o represents the height difference between plane p and plane o; Δz pq =z p -z q represents the height difference of plane p relative to plane q; Represents the high-order derivative in the spatial domain; take Δz qp =Δz po , add and sort to get:

[0075]

[0076]

[0077] Among them, R 2n (z po ) is the remainder of the downward extension correction number.

[0078] Moreover, the specific implementation method of step 4 is: iteratively process the height difference Δz in the spatial domain, perform upward extension by setting a small distance, and perform downward extension processing, and then use the plane as a known plane to continue downward extension processing, and iterate until the required downward height is reached;

[0079] The downward extended Δz height plane is corrected by extending it upward to a known plane, and the correction difference is judged by the threshold. If the difference is higher than the threshold, it is considered that the order of the high-order derivative term is too low, so an order high-order derivative term is added and the above calculation is repeated until the difference is less than the threshold, and the calculation accuracy is considered to meet the requirements.

[0080] The advantages and positive effects of the present invention are:

[0081] The present invention constructs an upward extended geomagnetic map; calculates high-order derivatives based on robust upward extension according to the upward extended geomagnetic map; performs downward iterative extension processing in the space domain according to the calculation results of the high-order derivatives; and obtains a vertical high-precision geomagnetic map by iteratively correcting the results of the downward iterative extension processing in the space domain. The present invention uses the vertical spatial distribution of geomagnetism at each point on the geomagnetic reference plane as a function model, and the vertical high-order derivatives of the geomagnetic space domain as the parameters to be solved for fitting, transforms the downward extended geomagnetic construction into the fitting and extension based on the geomagnetic reference map and the upward extended geomagnetic map, and performs high-order derivative truncation through iterative processing. Therefore, this method can fundamentally avoid the problem of solving high-dimensional ill-conditioned matrices in the downward extension process, thereby obtaining a robust downward extension effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1 This is a three-dimensional schematic diagram of the extended passive area and boundary of the Earth's magnetic field.

[0083] Figure 2 It is a flow chart of the iterative downward extension method combining the spatial domain and the spectral domain of the present invention;

[0084] Figure 3 It is a schematic diagram of constructing a geomagnetic map based on the upward extension of the spectrum domain according to the present invention;

[0085] Figure 4 It is a schematic diagram of the extension process of the present invention;

[0086] Figure 5 It is a schematic diagram of an iterative downward extension algorithm combining the spatial domain and the spectral domain of the present invention;

[0087] Figure 6 The conventional geomagnetic background images at different altitudes of the embodiment of the present invention;

[0088] Figure 7 This is a downward extension result and error distribution diagram based on different iterative methods of conventional geomagnetic maps in an embodiment of the present invention. DETAILED DESCRIPTION

[0089] The present invention is further described in detail below with reference to the accompanying drawings.

[0090] A spectral domain-spatial domain iterative spatial geomagnetic map construction method, such as Figure 2 As shown in the figure, based on the spectrum domain results of robust upward extension, the first law of geography, that is, the smoothness of the spatial geomagnetic field distribution, is used to fit the geomagnetic spatial distribution through the high-order derivative terms in the spatial domain to achieve the downward extension of the geomagnetic map. The following steps are included:

[0091] Step 1: Construct an upward extended geomagnetic map.

[0092] Performing Fourier transform on the upward extension equation (2), we get its spectral domain expression:

[0093]

[0094] because

[0095]

[0096] in k x , k y They represent the spectral domain variables corresponding to the spatial domain variables x and y, ξ is the dimension of the observed data along the x direction, and η is the dimension of the observed data along the y direction;

[0097]

[0098] Among them, L x , L y Respectively represent the fundamental wavelengths in the x and y directions, and they are:

[0099] L x =MΔx,L y =NΔy(12)

[0100] Wherein, Δx and Δy represent the point distance and line distance of the observation data respectively, the number of sampling points along the x direction is M, and the number of sampling points along the y direction is N.

[0101] When M and N are odd numbers,

[0102]

[0103] When M and N are even numbers,

[0104]

[0105] From this, the spectral domain calculation formula for the upward extension height h of the geomagnetic map plane u0 is obtained as follows:

[0106] U -h (u,v)=F[u0(x,y)]e -2πfh (15)

[0107] By performing an inverse Fourier transform, the geomagnetic map of the upward extension plane can be calculated:

[0108] u -h (x,y)=F -1 {F[u0(x,y)]e -2πfh} (16)

[0109] Step 2: Calculate high-order derivatives based on robust upward extension according to the upward extended geomagnetic map.

[0110] By performing inverse Fourier transform on equation (15), the geomagnetic map of the upward extension plane can be calculated:

[0111] u -h (x,y)=F -1 {F[u0(x,y)]e -2πfh} (16)

[0112] From the above process, it can be seen that the calculation process of the upward extension of the geomagnetic field does not require inverse calculation. The stable upward extension geomagnetic map can be obtained through fast Fourier transform and inverse transform. In this way, the geomagnetic information at multiple heights can be obtained. The calculation method is as follows: Figure 3 As shown. It is known that the p plane is the measurement plane, the o plane is the flight plane of the magnetic detection aircraft, the s plane is the sea level, q1~q n are n planes extended upward based on the p plane. The coordinate system is located on the p plane, downward is positive, and the height from the p plane to the o plane is h, and the height from the p plane to q1~q n The heights of the surfaces are -△z, -2△z, -3△z, ..., -n△z. The height from the p surface to the s surface is l. According to the upward extension equation, we can get:

[0113]

[0114] Using formula (17) and combining the geomagnetic map constructed with the known plane p, we can get the different levels q1~q n The geomagnetic map of is given by:

[0115]

[0116] This formula can be used to extend the geomagnetic map of the known reference layer upward to a geomagnetic map with a height of nΔz. n Based on the plane geomagnetic map, we further expand the Taylor series of equation (18) in the space domain at a single coordinate point on the plane to obtain:

[0117]

[0118] in, is the plane geomagnetic map obtained by upward extension, n represents the number of upward extension layers, u p is the plane geomagnetic map of the observation plane, and Δz is the height difference of the upward extension plane. Equation (19) describes the vertical spatial geomagnetic field distribution of a single coordinate point on the plane at different heights, that is, for a single coordinate point on the observation plane p and the upward extension plane q1...q n For a single coordinate point above, a magnetic function with respect to height z is formed in the vertical direction.

[0119] In formula (19), u pand Δz are both known quantities. It is the high-order gradient term of the vertical geomagnetic field in the spatial domain and is an unknown quantity. n These n upward extension planes can construct n linear equations to solve the high-order gradient terms of the vertical geomagnetic field in the spatial domain. Take point as an example (such as Figure 3 As shown), according to formula (19), the different heights of the first point q1...q n The geomagnetic data equations are as follows:

[0120]

[0121] Write (20) in the form of a matrix equation:

[0122] Y m =A m X m (twenty one)

[0123] in,

[0124] By using the least squares method to solve equation (21), we can obtain:

[0125] X m =(A m T A m ) -1 A m T Y m (twenty two)

[0126] Step 3: According to the calculation results of the high-order derivative terms, perform downward iterative extension processing in the spatial domain.

[0127] Equations (20) and (21) describe the process of solving the spatial domain high-order derivatives of a single plane coordinate point. By traversing all coordinate points on the observation plane, the spatial domain high-order derivatives of all points can be obtained. On the basis of obtaining the spatial domain high-order derivatives of all points, according to the Morit potential field analytical extension theory, the downward extension plane o and the upward extension plane q are respectively:

[0128]

[0129] Where Δz po =z p -z o represents the height difference between plane p and plane o; Δz pq =z p -z q represents the height difference of plane p relative to plane q; is the high-order derivative in the spatial domain obtained by (22).

[0130] Take Δz qp =Δz po , add and rearrange equation (23) to obtain:

[0131]

[0132] Among them, R 2n (z po ) is the remainder of the downward extension correction number, which can be calculated by formula (22). Figure 4 As shown, plane o is assumed to be the downward extension plane. The definitions of plane p and plane q are consistent with the upward extension part.

[0133] At the downward extension plane o and the upward extension plane q, respectively:

[0134]

[0135] Where Δz po =z p -z o represents the height difference between plane p and plane o; Δz pq =z p -z q represents the height difference of plane p relative to plane q; is the high-order derivative in the spatial domain obtained by (22).

[0136] Take Δz qp =Δz po , add and rearrange equation (23) to obtain:

[0137]

[0138] Among them, R 2n (z po ) is the remainder of the downward extension correction number and can be calculated by formula (22).

[0139] Step 4: Iteratively extend the processing results downward in the spatial domain through iterative correction processing to obtain a vertical high-precision geomagnetic map.

[0140] Equations (20) to (24) describe the principle of the downward extension method. It should be pointed out that, in theory, the high-order derivatives of the vertical geomagnetic field need to be taken to positive infinity, that is, t→+∞. At the same time, the more planes are extended upward, the more conducive it is to ensure the well-posedness of equation (22). However, in fact, the matrix dimension can be reduced by truncating the high-order derivatives, and the matrix dimension can be designed by the spectral domain-spatial domain iteration method, that is, the values ​​of n and t are determined. The schematic diagram of the specific iterative algorithm is shown in the figure. Figure 5 shown.

[0141] First, the iterative processing is performed on the height difference Δz in the spatial domain. Since the iterative downward extension algorithm combining the spatial domain and the spectral domain is based on the first law of geography (i.e., the smoothness of the spatial distribution of the geomagnetic field), and this law is more suitable when the height difference Δz is small. By setting a small distance for upward extension (such as Figure 5 (a) shown in the figure), use equation (24) to perform downward extension, and then use the plane as a known plane to continue downward extension (as shown in the figure). Figure 5 As shown, (b) its rectangular area moves to the left by Δz) compared to (a) the rectangular area, and the above (a) and (b) processes are iterated until the required downward height is reached.

[0142] Secondly, to ensure that the truncation of high-order derivatives is reasonable, it is necessary to correct the downward extended Δz height plane by extending it upward to a known plane, and make a judgment based on the corrected difference and the threshold. If the difference is higher than the threshold, it is considered that the order of the high-order derivative term is too low, so an order high-order derivative term is added and the above calculation is repeated until the difference is less than the threshold, and the calculation accuracy is considered to meet the requirements.

[0143] According to the above-mentioned spectral domain-spatial domain iterative spatial geomagnetic map construction method, the construction and calculation of the spatial geomagnetic map based on simulation data have verified the effect of the present invention.

[0144] In order to verify the feasibility of the algorithm, the present invention adopts a conventional model for verification, wherein the model parameters are shown in Table 1, the burial depth is set to 200m underground, and the flight altitude Z=200m ( Figure 6 (a)) and Z = 500m( Figure 6 (b)) The geomagnetic map of the height plane is taken as the true value of the geomagnetic map. The geomagnetic map of the Z = 500m plane is extended downward by 300m using the iterative downward extension method combining Landweber iteration and spatial domain-spectral domain, and the extended geomagnetic map is compared with the true value geomagnetic map of Z = 200m to obtain the downward extension error. The results are shown in Figure 7 shown.

[0145] Table 41 Parameters of conventional geomagnetic background map simulation model

[0146]

[0147] The calculation results show that:

[0148] (1) The iterative downward extension method combining the spatial domain and the spectral domain can obtain a good extension effect for the main area of ​​the magnetic source. Except for the oscillation of the boundary, the extension result is almost completely consistent with the plane true value.

[0149] (2)(2) The Gibbs effect in the boundary area of ​​the Landweber method is serious, but the high-frequency components of the main body of the magnetic source are smoothed, resulting in oscillation distortion in the area surrounding the magnetic source.

[0150] The extension error distribution can be intuitively seen from the error matrix:

[0151] (1) The Landweber iteration method produces a continuation error greater than 100 nT in the center of the magnetic source and the adjacent area;

[0152] (2)(2) The error of the iterative downward extension method combining the spatial domain and the spectral domain is much smaller in the area with small gradient and gentle change (i.e., medium and low frequency components). At the same time, the calculated MAE and RMSE are lower than those of the Landweber iteration method (as shown in Table 2). The mean absolute error MAE is reduced by 63.1%, and the root mean square error RMSE is reduced by 83.7%, which illustrates the effectiveness of the iterative downward extension method combining the spatial domain and the spectral domain.

[0153] Table 2 Comparison of the results of downward extension of conventional geomagnetic maps by different iterative methods

[0154]

[0155] It should be emphasized that the embodiments described in the present invention are illustrative rather than restrictive. Therefore, the present invention includes but is not limited to the embodiments described in the specific implementation manner. Any other implementation manners derived by those skilled in the art based on the technical solution of the present invention also fall within the scope of protection of the present invention.

Claims

1. A spectral domain-spatial domain iterative spatial geomagnetic map construction method, characterized by: The following steps are involved: Step 1, constructing an upward extended geomagnetic map; Step 2, calculating high-order derivatives based on robust upward continuation according to the upward extended geomagnetic map; Step 3: Perform downward iterative extension processing in the spatial domain according to the calculation results of the high-order derivative terms; Step 4: Iteratively extend the processing results downward in the spatial domain through iterative correction processing to obtain a vertical high-precision geomagnetic map.

2. The method for constructing a spectral domain-spatial domain iterative spatial geomagnetic map according to claim 1, characterized in that: The specific implementation method of step 1 is: the integral formula for upward extension is: Where h is the extension height, and h>0, and u(x,y,-h) is the plane to be extended Γ Up Magnetic field data, u(x,y,0) is the observation plane Γ S Magnetic field data, where ξ is the dimension of the data observed along the x direction, and η is the dimension of the data observed along the y direction; Performing Fourier transform on the upward extension equation, we get its spectral domain expression as follows: because in k x , k y Respectively represent the spectral domain variables corresponding to the spatial domain variables x and y: Among them, L x , L y Represent the fundamental wavelength in the x and y directions respectively: L x =MΔx,L y =NΔy Wherein, Δx and Δy represent the point distance and line distance of the observation data respectively, the number of sampling points along the x direction is M, and the number of sampling points along the y direction is N; When M and N are odd numbers, When M and N are even numbers, The spectral domain calculation formula for the extension of the geomagnetic map plane u0 upward to height h is: U -h (u,v)=F[u0(x,y)]e -2πfh Perform inverse Fourier transform on it and calculate the geomagnetic map of the upward extension plane: u -h (x,y)=F -1 {F[u0(x,y)]e -2πfh }。 3. The method for constructing a spectral domain-spatial domain iterative spatial geomagnetic map according to claim 2, characterized in that: The specific implementation method of step 2 is: according to the upward continuation equation, we get: Among them, the p plane is the measurement plane, the o plane is the flight plane of the magnetic detection aircraft, the s plane is the sea level, and q1~q n are n planes extended upward based on the p plane. The coordinate system is located on the p plane, downward is positive, and the height from the p plane to the o plane is h, and the height from the p plane to q1~q n The heights of the surfaces are -△z, -2△z, -3△z, ..., -n△z, and the height from the p surface to the s surface is l; Combined with the geomagnetic map constructed by the known plane p, we can get different levels q1~q n Magnetic map of: Different levels q1~q n The geomagnetic map of the known reference layer can be extended upward to a geomagnetic map with a height of nΔz. n Based on the plane geomagnetic map, we further analyze the different levels q1~q at the plane single coordinate point. n The geomagnetic map is expanded by Taylor series in the space domain: Among them, u qn is the plane geomagnetic map obtained by upward extension, n represents the number of upward extension layers, u p is the plane geomagnetic map of the observation plane, Δz is the height difference of the upward extension plane; Taylor series expansion describes the vertical space geomagnetic field distribution of a single coordinate point on the plane at different heights. For a single coordinate point on the observation plane p and the upward extension plane q1...q n For a single coordinate point above, a magnetic function with respect to height z is formed in the vertical direction. Get the different heights of the first point q1...q n The equations for geomagnetic data are: Writing the system of equations in matrix form: AND m =A m X m in, Calculated using the least squares method: X m (A) m T TO m ) -1 TO m T AND m 4. The method for constructing a spectral domain-spatial domain iterative spatial geomagnetic map according to claim 1, characterized in that: The specific implementation method of step 3 is: At the downward extension plane o and the upward extension plane q, respectively: Where Δz po =z p -z o represents the height difference between plane p and plane o; Δz pq =z p -z q represents the height difference of plane p relative to plane q; Represents the high-order derivative in the spatial domain; take Δz qp =Δz po , add and sort to get: Among them, R 2n (z po ) is the remainder of the downward extension correction number.

5. The method for constructing a spectral-spatial domain iterative spatial geomagnetic map according to claim 1, characterized in that: The specific implementation method of step 4 is: iteratively process the height difference Δz in the spatial domain, perform upward extension by setting a small distance, and perform downward extension processing, and then use the plane as a known plane to continue downward extension processing, and iterate until the required downward height is reached; The downward extended Δz height plane is corrected by extending it upward to a known plane, and the correction difference is judged by the threshold. If the difference is higher than the threshold, it is considered that the order of the high-order derivative term is too low, so an order high-order derivative term is added and the above calculation is repeated until the difference is less than the threshold, and the calculation accuracy is considered to meet the requirements.

Images