Three-dimensional positioning and magnetic moment direction evaluation method based on magnetic anomaly vector

Abstract

The invention belongs to the field of magnetic exploration, and relates to a three-dimensional positioning and magnetic moment direction evaluation method based on a magnetic anomaly vector, which comprises the following steps of: obtaining uniformly distributed two-dimensional magnetic vector data; carrying out convolution on the two-dimensional magnetic vector data by using a magnetic vector two-dimensional orthogonal basis function which has a preset height and is subjected to window truncation to obtain all orthogonal basis coefficients of a magnetic vector two-dimensional orthogonal basis; calculating the sum of all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector and the sum of squares of all the orthogonal basis coefficients to obtain energy data of a pseudo l1 norm and a pseudo l2 norm of each magnetic abnormal component; performing energy data fusion to obtain inversion energy and total energy; performing two-dimensional constant false alarm rate detection on the total energy to obtain the number of targets and the horizontal position of each target; performing parameter inversion of a vertical position and a magnetic moment direction; the vertical positions and the magnetic moment directions of all the targets are evaluated, and the situation that omission of the targets with the small energy peak value occurs due to the large difference between the energy peaks is avoided.

Description

Technical Field

[0001] The present invention belongs to the field of magnetic exploration, and specifically relates to a three-dimensional positioning and magnetic moment direction evaluation method based on magnetic anomaly vectors. Background Art

[0002] Magnetic anomaly refers to the change in the spatial distribution of the geomagnetic field caused by ferromagnetic targets. Magnetic anomaly detection (MAD) has become an attractive technology for passive detection and positioning of ferromagnetic sources in pipelines, mines and unexploded ordnance due to its concealment, immunity to environmental conditions and non-contact measurement. Although scalar sensors have low noise levels, they are not only expensive, but also affected by the inclination and declination of the geomagnetic field when used for magnetic anomaly positioning and cannot fully utilize the magnetic anomaly information. Vector sensors such as fluxgates can not only fully utilize the magnetic anomaly information but are also inexpensive. Magnetic anomaly positioning and magnetic moment assessment provide valuable technical support for scientific research, resource development and safety monitoring in various fields. They are widely used and have important application value. Accurate magnetic target positioning and magnetic moment direction estimation are crucial for the successful removal of unexploded ordnance and precise navigation of vehicles.

[0003] At present, the location of magnetic anomalies and the evaluation of magnetic moment direction are mainly carried out around dual sensors or even sensor arrays. The increase in the number of sensors not only increases the cost of the magnetometer exponentially, but also increases the difficulty of operation and noise interference. The few existing target position and magnetic moment parameter evaluation methods based on single magnetic sensors are mainly through solving linear equations and nonlinear equations. For linear equation methods, the linearization of nonlinear equations will introduce method errors. For nonlinear equation methods, these methods currently mainly transform them into optimization problems and use intelligent algorithms to search parameters. However, such methods are highly sensitive to the number of parameters and noise, especially in multi-target situations, and also require prior information such as the number of targets. Summary of the invention

[0004] The present invention provides a three-dimensional positioning and magnetic moment direction evaluation method based on magnetic anomaly vector, which solves the problems of high sensitivity to the number of parameters and noise and dependence on the number of targets in a multi-target situation.

[0005] The present invention is achieved in this way.

[0006] A three-dimensional positioning and magnetic moment direction evaluation method based on magnetic anomaly vector, the method comprising:

[0007] The one-dimensional survey line magnetic anomaly data obtained by the three-component magnetometer is preprocessed, interpolated and gridded to obtain uniformly distributed two-dimensional magnetic vector data;

[0008] The two-dimensional magnetic vector data is convolved with a two-dimensional orthogonal basis function of a magnetic vector with a preset height and truncated by a window to obtain all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector;

[0009] The sum of all orthogonal basis coefficients and the sum of squares of all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector are calculated to obtain energy data of pseudo-l1 norm and pseudo-l2 norm of each magnetic anomaly component;

[0010] The energy data of the pseudo-l1 norm and pseudo-l2 norm of each magnetic anomaly component are fused to obtain the inversion energy and the total energy;

[0011] Perform two-dimensional constant false alarm rate detection on the total energy to obtain the number of targets and the horizontal position of each target;

[0012] According to the horizontal position of the target, the inversion model composed of the inversion energy is combined to perform the parameter inversion of the vertical position and magnetic moment direction;

[0013] After evaluating the vertical position and magnetic moment direction of all targets, determine whether the evaluation of all targets is completed. If not, continue to repeat the parameter inversion of the vertical position and magnetic moment direction until the parameter inversion of all targets is completed.

[0014] Furthermore, all coefficients of the two-dimensional orthogonal basis of magnetic vectors are obtained, including:

[0015] The two-dimensional orthogonal basis function of magnetic vector is expressed as follows:

[0016]

[0017] where f Bkl , k = x, y, z, l = 1, 2, 3 represents the orthogonal basis function of the preset height; x, y are the horizontal coordinates; the prior information z is the preset height, and the two-dimensional orthogonal basis function of each magnetic vector meets the following orthogonality and standardization requirements:

[0018]

[0019] A window of preset size is used to intercept all two-dimensional orthogonal basis functions of magnetic vectors, so that each magnetic anomaly component is convolved with each two-dimensional orthogonal basis function of each magnetic vector to obtain the corresponding orthogonal basis coefficient a. Bkl (m,n);

[0020] The convolution operation is:

[0021]

[0022] Where m0 is the modulus of the magnetic moment vector; μ0 is the magnetic permeability of vacuum; B kis the two-dimensional magnetic anomaly vector data; I and J represent the sliding window size, which is used to control the boundary of discrete numerical integration; Δx and Δy are the discrete integration intervals along the x-axis and y-axis respectively; m, n represent the horizontal and vertical coordinates of the orthogonal basis coefficients; x m:i ,y n:j Represents the horizontal and vertical coordinates of the magnetic anomaly component; x i ,y j Represents the horizontal and vertical coordinates of the two-dimensional orthogonal basis function of the magnetic vector.

[0023] Furthermore, the inversion energy is: E Inv =∑ k<x,y,z E k1 , E Inv Represents the inversion energy, the total energy is: E tot =∑ k<x,y,z E k2 , E tot represents the total energy;

[0024] E k1 Indicates magnetic anomaly B k The pseudo l1 norm energy of the component, E k2 Indicates magnetic anomaly B k Pseudo-l2-norm energy of the component.

[0025] Further, k is the category of magnetic anomaly component, a Bkl with a Bkl (m,n) have the same meaning.

[0026] Furthermore, a two-dimensional constant false alarm rate detection is performed on the total energy to obtain the number of targets and the horizontal position of each target. The calculation process includes:

[0027] N R =N TR +N GR

[0028] N C =N TC +N GC

[0029] N=(2N R +1)(2N C +1)-(2N GR +1)(2N GC +1)

[0030]

[0031]

[0032] if E tot >T(x m ,yn )→X s ,Y s = argmax E tot

[0033] Among them, N TR and N TC are the number of rows and columns of training units respectively; N GR and N GC are the number of rows and columns of protection units respectively; N is the number of training units; P FA is the probability of false alarm; E(x m:i ,y n:j ) indicates that the coordinate is x m:i ,y n:j energy, T is the dynamic threshold; θ is the threshold factor; N R is the sum of the number of rows of training units and protection units; N C is the sum of the number of columns of training units and protection units; x m ,y n The horizontal and vertical coordinates represent the dynamic threshold; X s ,Y s is the horizontal position of the target obtained, and the subscript s is the serial number of the evaluation target.

[0034] Furthermore, according to the horizontal position of the target, the inversion model composed of the inversion energy is combined to perform parameter inversion of the vertical position and magnetic moment direction, specifically including:

[0035] The inversion model of inversion energy composition is:

[0036]

[0037] Z s ,I s ,D s =argmin LS s

[0038] Among them, LS s is the loss function value, Z s is the vertical position of the inversion, I s is the inverted magnetic moment inclination, D s is the inverted magnetic moment deflection, and the subscript s is the serial number of the evaluation target; is the estimated inversion energy; P and Q are used to control the size of the matching window, p and q are index subscripts, x s:p ,y s:q are the coordinates of the inverted energy and the estimated inverted energy;

[0039] The horizontal position of the target is substituted into the inversion model composed of inversion energy to perform parameter inversion of the vertical position and magnetic moment direction.

[0040] Compared with the prior art, the present invention has the following beneficial effects:

[0041] Compared with the fixed threshold, the method of the present invention uses two-dimensional constant false alarm rate detection without setting any threshold, which largely avoids the situation that the target with small energy peak value is missed due to the large difference between energy peaks. The pseudo-l2 norm of the orthogonal basis coefficient is used to construct the total energy, strengthen the noise suppression ability, and improve the horizontal positioning accuracy; the pseudo-l1 norm of the orthogonal basis coefficient is used to construct the inversion energy, which avoids the spatial symmetry of the total energy, eliminates the multi-solution of the optimization inversion, and improves the inversion accuracy. Compared with the existing method of directly using magnetic anomalies to construct the optimization function and perform parameter optimization inversion, the method of the present invention not only does not require prior information on the number of magnetic targets, does not need to iterate to confirm the number of magnetic targets, has a small number of optimization parameters, reduces the sensitivity to noise, and improves the accuracy and efficiency of the inversion. Compared with the existing two-dimensional orthogonal basis method of the total magnetic field, the method of the present invention is insensitive to the direction of the magnetic moment, has rich inversion results, high inversion accuracy and low false alarm rate. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 A flowchart of a method provided by an embodiment of the present invention;

[0043] Figure 2 The two-dimensional orthogonal basis surface of magnetic vectors and its contour map provided by the embodiment of the present invention, (a) f1 = f Bx1 , (b) f2=f By1 , (c)f3=f Bx2 =f By2 ,(d)f4=f Bz2 =f By3 , (e)f5=f Bx3 =f Bz3 ,(f)f6=f Bz1 ;

[0044] Figure 3 Two-dimensional magnetic anomaly (a)B generated by simulation according to an embodiment of the present invention x Quantity, (b) B y Components and (c) B z Quantity;

[0045] Figure 4 (a) Energy B of the three components of magnetic anomaly in the simulation experiment provided by the embodiment of the present invention x2 (b) Energy B y2 and (c) energy B z2 ;

[0046] Figure 5 The magnetic anomaly total energy surface and its contour lines of the simulation experiment provided by the embodiment of the present invention;

[0047] Figure 6 A two-dimensional constant false alarm rate detection result diagram of the total energy of magnetic anomaly in a simulation experiment provided by an embodiment of the present invention;

[0048] Figure 7 Two-dimensional magnetic anomaly of field data provided by an embodiment of the present invention (a) B x Quantity, (b) B y Components and (c) B z Quantity;

[0049] Figure 8 (a) Energy B of the three components of magnetic anomaly in field data provided by the embodiment of the present invention x2 (b) Energy B y2 ) and (c) energy B z2 ;

[0050] Fig. 9 The magnetic anomaly total energy surface and contour lines of field data provided by the embodiment of the present invention;

[0051] Fig.10 An embodiment of the present invention provides a two-dimensional constant false alarm rate detection result diagram of the total energy of magnetic anomaly for field data. DETAILED DESCRIPTION

[0052] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0053] See also Figure 1 As shown, a three-dimensional positioning and magnetic moment direction evaluation method based on magnetic anomaly vector, the method comprises:

[0054] The one-dimensional survey line magnetic anomaly data obtained by the three-component magnetometer is preprocessed, interpolated and gridded to obtain uniformly distributed two-dimensional magnetic vector data;

[0055] The two-dimensional magnetic vector data is convolved with a two-dimensional orthogonal basis function of a magnetic vector with a preset height and truncated by a window to obtain all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector;

[0056] The sum of all orthogonal basis coefficients and the sum of squares of all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector are calculated to obtain energy data of pseudo-l1 norm and pseudo-l2 norm of each magnetic anomaly component;

[0057] Fusing the energy of each magnetic anomaly component to obtain inversion energy and total energy;

[0058] Perform two-dimensional constant false alarm rate detection on the total energy to obtain the number of targets and the horizontal position of each target;

[0059] According to the horizontal position of the target, the inversion model composed of the inversion energy is combined to perform the parameter inversion of the vertical position and magnetic moment direction;

[0060] After evaluating the vertical position and magnetic moment direction of all targets, determine whether the evaluation of all targets is completed. If not, continue to repeat the parameter inversion of the vertical position and magnetic moment direction until the parameter inversion of all targets is completed.

[0061] Among them, in the embodiment of the present invention, a strapdown MAD device composed of a fluxgate sensor and an inertial navigation system is used to obtain one-dimensional survey line magnetic anomaly data, attitude data and coordinate data. Specifically, after determining the survey area, the device is moved along several parallel survey lines for continuous measurement, and the strapdown MAD device collects magnetic anomaly vectors, attitude and coordinate data in a continuous sampling manner (the sampling frequency is set to 10Hz). Then, the obtained one-dimensional survey line magnetic anomaly data is subjected to sensor error correction, misalignment error correction, attitude restoration and detrending operations, and the interpolation means selects a linear interpolation method; after interpolation and gridding, the intervals of the data in both directions along the survey line and the vertical survey line are 0.1m. Obtain two-dimensional magnetic vector data that does not contain the geomagnetic field in the geographic coordinate system.

[0062] The two-dimensional magnetic vector data is convolved using the two-dimensional orthogonal basis function of the magnetic vector with a preset height and truncated by a window to obtain all the orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector. Including:

[0063] The two-dimensional orthogonal basis function of magnetic vector is expressed as follows:

[0064]

[0065] where f Bkl , k = x, y, z, l = 1, 2, 3 represents the orthogonal basis function of the preset height; x, y are the horizontal coordinates; the prior information z is the preset height, and the two-dimensional orthogonal basis function of each magnetic vector meets the following orthogonality and standardization requirements:

[0066]

[0067] A window of preset size is used to intercept all two-dimensional orthogonal basis functions of magnetic vectors, so that each magnetic anomaly component is convolved with each two-dimensional orthogonal basis function of each magnetic vector to obtain the corresponding orthogonal basis coefficient a. Bkl (m,n);

[0068] In this embodiment, the truncation window is set to a square with a side length of 5m, that is, the side length of the sliding window is set to 5m (2I=2J=5m), and the window is used to truncate all two-dimensional orthogonal basis functions of magnetic vectors.

[0069] The convolution operation is:

[0070]

[0071] Where m0 is the modulus of the magnetic moment vector; μ0 is the magnetic permeability of vacuum; B k is the two-dimensional magnetic anomaly vector data; I and J represent the sliding window size, which is used to control the boundary of discrete numerical integration; Δx and Δy are the discrete integration intervals along the x-axis and y-axis respectively; m, n represent the horizontal and vertical coordinates of the orthogonal basis coefficients; x m:i ,y n:j Represents the horizontal and vertical coordinates of the magnetic anomaly component; x i ,y j Represents the horizontal and vertical coordinates of the two-dimensional orthogonal basis function of the magnetic vector.

[0072] In this embodiment, the preset height z=1m, and the surface and contour map of the two-dimensional orthogonal basis function (OBF) of the magnetic vector with the range of x (west-east) and y (north-south) both being 10m are calculated, as shown in FIG. Figure 2 shown. Figure 2 The OBF f1 in (a) is the basis function f Bx1 The surface diagram of Figure 2 (b) OBF f2 is the basis function f By1 The surface diagram of Figure 2 The OBF f3 in (c) is the two-dimensional orthogonal basis function f of the magnetic vector Bx2 and f By2 The surface diagram of Figure 2 The OBF f4 in (d) is the two-dimensional orthogonal basis function f of the magnetic vector Bz2 and f By3 The surface diagram of Figure 2 (e)OBF f5 in is the two-dimensional orthogonal basis function f of magnetic vector Bx3 and f Bz3 The surface diagram of Figure 2 The OBF f6 of (f) is the two-dimensional orthogonal basis function f of the magnetic vector Bz1 Surface map of . Figure 2 (a) and Figure 2 (b) shows that the rotation of x and y conforms to the law of two-dimensional orthogonal basis function of magnetic vector; Figure 2 (d) and Figure 2 (e) in the figure also shows the rotation of x and y, which conforms to the law of two-dimensional orthogonal basis functions of magnetic vectors.

[0073] The sum of all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector and the sum of the squares of all orthogonal basis coefficients are calculated to obtain the energy data of the pseudo-l1 norm and pseudo-l2 norm of each magnetic anomaly component. The energy data includes: k is the category of magnetic anomaly component, Ek1 Indicates magnetic anomaly B k The pseudo l1 norm energy of the component, E k2 Indicates magnetic anomaly B k Pseudo-l2-norm energy of the component.

[0074] The energy of each magnetic anomaly component is integrated to obtain the inversion energy and the total energy. In essence, the sum of all orthogonal basis coefficients (pseudo-l1 norm) and the sum of the squares of all orthogonal basis coefficients (pseudo-l2 norm) are calculated to obtain the inversion energy E. Inv and the total energy E tot , and finally perform energy normalization. The normalization of total energy is conducive to comparative analysis during multi-target inversion; the normalization of inversion energy can eliminate a magnetic moment modulus parameter, thus significantly increasing the inversion accuracy.

[0075] E Inv =∑ k<x,y,z E k1

[0076]

[0077] Perform two-dimensional constant false alarm rate detection on the total energy to obtain the number of targets and the horizontal position of each target. The calculation process includes:

[0078] N R =N TR +N GR

[0079] N C =N TC +N GC

[0080] N=(2N R +1)(2N C +1)-(2N GR +1)(2N GC +1)

[0081]

[0082] if E tot >T(x m ,y n )→X s ,Y s = argmax E tot

[0083] Among them, N TR and N TC are the number of rows and columns of training units respectively; N GR and N GC are the number of rows and columns of protection units respectively; N is the number of training units; PFA is the probability of false alarm; E(x m:i ,y n:j ) indicates that the coordinate is x m:i ,y n:j energy, T is the dynamic threshold; θ is the threshold factor; N R is the sum of the number of rows of training units and protection units; N C is the sum of the number of columns of training units and protection units; x m ,y n The horizontal and vertical coordinates represent the dynamic threshold; X s ,Y s is the horizontal position of the target obtained, and the subscript s is the evaluation target number.

[0084] The dynamic threshold of the entire measurement area is obtained by performing a two-dimensional constant false alarm rate test on the normalized total energy. Then, the total energy and the dynamic threshold are compared. The areas with total energy higher than the dynamic threshold are regarded as areas containing targets, and the number of these areas is used as the evaluation value of the number of targets. The horizontal coordinate of the energy maximum in each area is used as the evaluation value of the horizontal position of the corresponding target.

[0085] The number of rows in the training unit N TR and the number of columns N TC Set to 40 (4m) points; the number of rows of protection units is N GR and the number of columns N GC Set to 10 (1m) points; N is the number of training units; set the probability of false alarm P FA is 0.01; E(x m:i ,y n:j ) indicates that the coordinate is x m:i ,y n:j energy, T represents the dynamic threshold; X s ,Y s is the target horizontal position obtained, and the subscript s is the serial number of the evaluation target.

[0086] According to the horizontal position of the target, the inversion model composed of the inversion energy is used to invert the parameters of the vertical position and magnetic moment direction, including:

[0087] The inversion model of the inversion energy composition is:

[0088]

[0089] Z s ,I s ,D s =argmin LS s

[0090] Among them, LS s is the loss function value, Z sis the vertical position of the inversion, I s is the inverted magnetic moment inclination, D s is the inverted magnetic moment deflection, and the subscript s is the serial number of the evaluation target; is the estimated inversion energy; P and Q are used to control the size of the matching window, p and q are index subscripts, x s:p ,y s:q are the coordinates of the inverted energy and the estimated inverted energy;

[0091] The horizontal position of the target is substituted into the inversion model composed of inversion energy to perform parameter inversion of the vertical position and magnetic moment direction.

[0092] In order to improve the inversion accuracy and calculation efficiency, it is necessary to select a matching window of appropriate size when performing energy matching. The selection principle aims to maximize the inclusion of energy information of a specific target while minimizing interference from other targets. In this embodiment, a square window with a side length of 3m (2P=2Q=3m) is selected as the matching window, and the horizontal position of the target is used as the center position of the matching window. In a multi-target scenario, the horizontal position of each target is also used as the center position of the matching window. When the estimated inversion energy is closest to the actual inversion energy, that is, when the loss function value is the smallest, the inversion results of the vertical position and magnetic moment direction parameters are obtained.

[0093] In order to verify the effectiveness and practicability of the present invention, multi-objective multi-parameter inversion simulation and field experiments were carried out as demonstrations.

[0094] Simulation experiment:

[0095] The measurement line spacing is set to 0.5m, 30m long, with a total of 41 measurement lines, forming a measurement area of ​​20×30m 2 ; The noise level is Gaussian white noise with a standard deviation of 20nT; the detection plane is 1m above the ground. Taking the southwest corner of the survey area as the coordinate origin, 6 magnetic dipole targets are randomly generated in the survey area. The preset parameters (X, Y, Z, I, D) of targets 1 to 6 are (4m, 4m, 0.3m, -30°, -170°), (10m, 14m, 0.5m, 58°, 86°), (12m, 6m, 0.5m, 70°, 176°), (14m, 5m, 0.6m, 60°, 0°), (20m, 14m, 0.4m, 30°, -3°) and (24m, 8m, 0.2m, 0°, 45°). In order to verify that the method of the present invention is not affected by the magnetic moment modulus value and to ensure that the generated magnetic anomaly data has almost the same signal-to-noise ratio, the magnetic moment modulus value of the calculation end is set to 1A〃m 2 The estimated end moment modulus range is set to [1-100A"m 2 ].

[0096] According to the above conditions, after linear interpolation and gridding, the grid size of the data is 0.1m, and the three components of the calculated magnetic anomaly (Bx, By, Bz) are generated as follows: Figure 3 As shown, x (west-east) and y (north-south), can be obtained from Figure 3 Middle (a), Figure 3 (b) It is found that the magnetic anomaly signal is almost submerged by the noise. Figure 3 In (c), a weak magnetic anomaly signal can be observed in the Bz component, but its shape is also severely distorted by noise.

[0097] The side length of the sliding window is set to 5m, and the corresponding magnetic anomaly component is used to perform a two-dimensional convolution operation with the two-dimensional orthogonal basis function of the magnetic vector of the corresponding component to obtain the corresponding orthogonal basis coefficients. Then, the pseudo-l1 norm and pseudo-l2 norm of the orthogonal basis coefficients of the two-dimensional orthogonal basis of each magnetic vector are calculated to obtain the energy data of the corresponding component. The pseudo-l2 norm energy data of the three components of the magnetic anomaly are shown in Figure 4 As shown, it can be seen that compared with the magnetic anomaly component, the signal-to-noise ratio in the energy domain is significantly improved, and the magnetic target is very easy to detect. Figure 4 (a) and Figure 4 (b) Energy E x2 and energy E y2 The size of the energy E y2 The two targets with smaller energy in the middle target have energy E x2 In contrast, Figure 4 (c) Energy E z2 The energy peak values ​​of each target in the . The peak values ​​of the three components are used to locate the horizontal position, and it is found that the energy E z2 Higher precision positioning can also be achieved. However, the present invention only shows the positioning method with the highest precision, using the fusion energy of the three component energies to determine the horizontal position. The pseudo-l1 norm energy map of the three components of magnetic anomaly is similar to the pseudo-l2 norm energy map and is no longer shown.

[0098] Next, energy fusion is performed to calculate the energy of the three pseudo l1 norm components (E x1 、E y1 and E z1 ) and obtain the inversion energy E Inv , calculate the three pseudo l2 norm component energies (E x2 、E y2 and E z2 ) and obtain the total energy E tot After normalization, the total energy surface and its contour lines are as follows: Figure 5 As shown, it can be seen that energy fusion reduces the difference between the energy maxima of the magnetic dipoles and improves the regularity of their energy shapes, thereby reducing the false alarm rate and improving the accuracy of horizontal positioning. Figure 5Six energy peaks can be easily identified, including two adjacent peaks.

[0099] Then the total energy is tested with a two-dimensional constant false alarm rate. The test results are as follows: Figure 6 As shown, it can be seen that there are 6 energy peaks exceeding the dynamic threshold, so the number of detected targets is 6, and then the horizontal coordinates of the energy peaks of these 6 regions are obtained, that is, the horizontal positioning of the 6 targets is achieved. The horizontal positioning results (X, Y) of these 6 targets are (4.00m, 4.00m), (10.10m, 14.00m), (12.10m, 6.00m), (13.90m, 5.10m), (20.00m, 14.00m), (24.10m, 8.00m).

[0100] Substitute the horizontal position results of the above targets into the inversion model composed of inversion energy in turn, select a 3m square window for the matching window, and then use the intelligent algorithm to perform parameter inversion with the number of targets as the number of times, and use the marine predator intelligent optimization algorithm to perform 6 parameter inversions. The inversion results of the three-dimensional parameters (Z / m, I / °, D / °) of the six targets are (0.32m, -39.71°, -191.88°), (0.32m, 63.76°, 98.21°), (0.65m, 59.22°, 195.26°), (0.68m, 49.10°, -26.04°), (0.43m, 35.40°, -7.59°) and (0.22m, -0.82°, 41.55°). In this 5D parameter inversion simulation, the maximum inversion errors corresponding to these five parameters (X / m, Y / m, Z / m, I / °, D / °) are (0.10m, 0.10m, 0.18m, 10.9°, 26.04°). The estimated errors of the magnetic moment inclination and declination of the fourth target are the largest, which are 10.9° and 26.04° respectively. Considering the nearby targets, such a high noise level and the fact that no filtering operation is included, the inversion accuracy is still satisfactory.

[0101] Field experiment:

[0102] In the field experiment, the survey line was east-west, 20 m long, with a spacing of 0.5 m between survey lines, and a survey area of ​​20 × 20 m. 2 The detection plane is 1.1m above the ground. Taking the ground at the southwest corner of the survey area as the coordinate origin, three iron pipes of the same size and material are randomly placed in the survey area. The outer diameter of the iron pipe is 10cm, the thickness is 3mm, and the length is 30cm. Their three-dimensional coordinates are (7.80m, 9.62m, 0.05m), (13.76m, 12.96m, -0.05m), and (13.95m, 5.95m, 0.05m).

[0103] First, a strapdown MAD device consisting of a fluxgate sensor and an inertial navigation system (INS) is used to obtain one-dimensional line magnetic anomaly data, attitude data, and position data. Then, the Euler angle data of the INS is used to perform sensor error correction, misalignment error correction, attitude restoration, and detrending operations on the magnetic anomaly vector data to obtain one-dimensional line magnetic anomaly data that does not contain the geomagnetic field in the geographic coordinate system. Then, the one-dimensional line magnetic anomaly data is linearly interpolated and gridded to obtain two-dimensional magnetic vector data, that is, two-dimensional magnetic anomaly three-component data.

[0104] After the above processing, the two-dimensional magnetic vector data surface map and contour map with a grid size of 0.1m are obtained as shown in Figure 7 As shown. It can be seen that the noise is large. Figure 7 (a) Figure 7 The three components of the magnetic anomaly in (b) are almost all submerged by the noise, except Figure 7 (c) E z The component can distinguish weak magnetic anomaly signals.

[0105] The side length of the sliding window is set to 5m, so that each magnetic anomaly component and the two-dimensional orthogonal basis function of the magnetic vector in each group are subjected to a two-dimensional convolution operation to obtain the corresponding orthogonal basis coefficients, and then the pseudo-l1 norm and pseudo-l2 norm of the orthogonal basis coefficients of the two-dimensional orthogonal basis of each magnetic vector are calculated to obtain the energy data of the corresponding component. The pseudo-l2 norm energy of the three components of the magnetic anomaly is as follows: Figure 8 As shown, Figure 8 (a) Figure 8 (b) and Figure 8 In (c), we can see that the signal-to-noise ratio of the signal in the energy domain is high, and the signal is easy to detect.

[0106] Then, energy fusion is performed to calculate the three pseudo l1 norm component energies (E x1 、E y1 and E z1 ) and obtain the inversion energy E Inv , calculate the three pseudo l2 norm component energies (E x2 、E y2 and E z2 ) and obtain the total energy E tot After normalization, the total energy surface and its contour map are as follows: Fig. 9 As shown, it can be seen that the signal-to-noise ratio of the total energy is higher, and the magnetic target is more easily detected.

[0107] Then the total energy is subjected to a two-dimensional constant false alarm rate test, and the test results are as follows: Fig.10As shown, we can see that there are 3 energy peaks exceeding the dynamic threshold, so the number of detected targets is 3, and then the peak coordinates of these three peaks are obtained, that is, the horizontal positioning of the three targets is achieved. The horizontal positioning (X, Y) results of the three targets are (7.70m, 9.25m), (13.45m, 13.05m), and (14.00m, 5.85m) respectively.

[0108] Substitute the above horizontal positioning results into the inversion model in turn, select a square window with a side length of 3m as the matching window, and then use the marine predator algorithm to perform three independent parameter optimization inversions. The inversion results of the three-dimensional parameters (Z / m, I / °, D / °) of the three targets obtained are (0.14m, 47.57°, -68.54°), (-0.24m, 43.50°, -147.95°) and (-0.12m, 28.39°, 68.26°). The inversion errors of the three-dimensional positions of the three targets in this field experiment are (0.10m, 0.37m, 0.09m), (0.31m, 0.09m, 0.19m) and (0.05m, 0.10m, 0.17m). Considering such large colored noise conditions, these results are acceptable, which verifies the effectiveness of the method of the present invention and its powerful noise suppression ability.

[0109] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A three-dimensional positioning and magnetic moment direction evaluation method based on magnetic anomaly vector, characterized in that: The method includes: The one-dimensional survey line magnetic anomaly data obtained by the three-component magnetometer is preprocessed, interpolated and gridded to obtain uniformly distributed two-dimensional magnetic vector data; The two-dimensional magnetic vector data is convolved with a two-dimensional orthogonal basis function of a magnetic vector with a preset height and truncated by a window to obtain all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector; The sum of all orthogonal basis coefficients and the sum of squares of all orthogonal basis coefficients of the two-dimensional orthogonal basis of the magnetic vector are calculated to obtain energy data of pseudo-l1 norm and pseudo-l2 norm of each magnetic anomaly component; The energy data of the pseudo-l1 norm and pseudo-l2 norm of each magnetic anomaly component are fused to obtain the inversion energy and the total energy; Perform two-dimensional constant false alarm rate detection on the total energy to obtain the number of targets and the horizontal position of each target; According to the horizontal position of the target, the inversion model composed of the inversion energy is combined to perform the parameter inversion of the vertical position and magnetic moment direction; After evaluating the vertical position and magnetic moment direction of all targets, determine whether the evaluation of all targets is completed. If not, continue to repeat the parameter inversion of the vertical position and magnetic moment direction until the parameter inversion of all targets is completed.

2. The method for three-dimensional positioning and magnetic moment direction evaluation based on magnetic anomaly vector according to claim 1 is characterized in that: Get all the coefficients of the two-dimensional orthogonal basis of magnetic vectors, including: The two-dimensional orthogonal basis function of magnetic vector is expressed as follows: where f Bkl , k = x, y, z, l = 1, 2, 3 represents the orthogonal basis function of the preset height; x, y are the horizontal coordinates; the prior information z is the preset height, and the two-dimensional orthogonal basis function of each magnetic vector meets the following orthogonality and standardization requirements: A window of preset size is used to intercept all two-dimensional orthogonal basis functions of magnetic vectors, so that each magnetic anomaly component is convolved with each two-dimensional orthogonal basis function of each magnetic vector to obtain the corresponding orthogonal basis coefficient a. Bkl (m,n); The convolution operation is: Where m0 is the modulus of the magnetic moment vector; μ0 is the magnetic permeability of vacuum; B k is the two-dimensional magnetic anomaly vector data; I and J represent the sliding window size, which is used to control the boundary of discrete numerical integration; Δx and Δy are the discrete integration intervals along the x-axis and y-axis respectively; m, n represent the horizontal and vertical coordinates of the orthogonal basis coefficients; x m:i ,y n:j Represents the horizontal and vertical coordinates of the magnetic anomaly component; x i ,y j Represents the horizontal and vertical coordinates of the two-dimensional orthogonal basis function of the magnetic vector.

3. The method for three-dimensional positioning and magnetic moment direction evaluation based on magnetic anomaly vector according to claim 2 is characterized in that: The inversion energy is: E Inv =∑ k<x,y,z E k1 , E Inv Represents the inversion energy, the total energy is: E tot =∑ k<x,y,z E k2 , E tot represents the total energy; E k1 Indicates magnetic anomaly B k The pseudo l1 norm energy of the component, E k2 Indicates magnetic anomaly B k Pseudo-l2-norm energy of the component.

4. The method for three-dimensional positioning and magnetic moment direction evaluation based on magnetic anomaly vector according to claim 3 is characterized in that: k is the category of magnetic anomaly component, a Bkl with a Bkl (m,n) have the same meaning.

5. The method for three-dimensional positioning and magnetic moment direction evaluation based on magnetic anomaly vector according to claim 4, characterized in that: Perform two-dimensional constant false alarm rate detection on the total energy to obtain the number of targets and the horizontal position of each target. The calculation process includes: N R =N TR +N GR N C =N TC +N GC <h2 style=";text-align:left;direction:ltr">N=(2N<h2 style=";text-align:left;direction:ltr"> R <h2 style=";text-align:left;direction:ltr"> +1)(2N<h2 style=";text-align:left;direction:ltr"> C <h2 style=";text-align:left;direction:ltr"> +1)-(2N<h2 style=";text-align:left;direction:ltr"> GR <h2 style=";text-align:left;direction:ltr"> +1)(2N<h2 style=";text-align:left;direction:ltr"> GC <h2 style=";text-align:left;direction:ltr"> +1) if E tot >T(x m ,y n )→X s ,Y s =argmax E tot Among them, N TR and N TC are the number of rows and columns of training units respectively; N GR and N GC are the number of rows and columns of protection units respectively; N is the number of training units; P FA is the probability of false alarm; E(x m:i ,y n:j ) indicates that the coordinate is x m:i ,y n:j energy, T is the dynamic threshold; θ is the threshold factor; N R is the sum of the number of rows of training units and protection units; N C is the sum of the number of columns of training units and protection units; x m ,y n The horizontal and vertical coordinates represent the dynamic threshold; X s ,Y s is the horizontal position of the target obtained, and the subscript s is the serial number of the evaluation target.

6. The method for three-dimensional positioning and magnetic moment direction evaluation based on magnetic anomaly vector according to claim 5, characterized in that: According to the horizontal position of the target, the inversion model composed of the inversion energy is used to invert the parameters of the vertical position and magnetic moment direction, including: The inversion model of inversion energy composition is: Z s ,I s ,D s =argmin LS s Among them, LS s is the loss function value, Z s is the vertical position of the inversion, I s is the inverted magnetic moment inclination, D s is the inverted magnetic moment deflection, and the subscript s is the serial number of the evaluation target; is the estimated inversion energy; P and Q are used to control the size of the matching window, p and q are index subscripts, x s:p ,y s:q are the coordinates of the inverted energy and the estimated inverted energy; The horizontal position of the target is substituted into the inversion model composed of inversion energy to perform parameter inversion of the vertical position and magnetic moment direction.

Images