Lightweight high-precision unmanned aerial vehicle aeromagnetic positioning method used in GNSS denial environment

Abstract

The invention provides a lightweight high-precision unmanned aerial vehicle aeromagnetic positioning method used in a GNSS denial environment, belongs to the technical field of high-precision positioning, and aims to solve the problems that the matching precision and the calculation complexity of an existing aeromagnetic matching positioning method are difficult to balance, and the aeromagnetic matching positioning is easy to cause matching failure due to similar results in a wide area. Comprising the steps of 1, synchronously collecting aeromagnetic observation values Mi to construct an aeromagnetic anomaly graph when the unmanned aerial vehicle flies in a specific area; step 2, carrying out gridding and down-sampling on the aeromagnetic anomaly graph and the flight path, and fitting a nearest flight path point of each grid node to a corresponding node to obtain a flight path slice; and step 3, carrying out iterative correlation analysis on the aeromagnetic anomaly graph and the flight path slices until the positioning resolution of an analysis result meets a requirement or the number of positioning iterations exceeds a preset threshold, and obtaining an aeromagnetic positioning result of the unmanned aerial vehicle.

Description

Technical Field

[0001] The present invention relates to a lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environments, belonging to the technical field of high-precision positioning. Background Art

[0002] With the wide application of the Global Navigation Satellite System (GNSS), the positioning and navigation technology of unmanned aerial vehicles (UAVs) has been greatly improved. However, GNSS signal denial poses challenges to navigation tasks in special scenarios such as UAVs. To address this issue, aeromagnetic positioning, as an emerging technology, has gradually become a research hotspot due to its advantages such as strong concealment, high anti-interference ability, and applicability to various weather conditions. Aeromagnetic positioning can achieve relatively fast and high-precision positioning without satellite signal support by measuring and matching the geomagnetic field, demonstrating great application potential in future wars.

[0003] Although aeromagnetic positioning has many advantages, its technology is still in the development stage, especially the wide-area aeromagnetic matching positioning method applied to UAVs. The existing aeromagnetic matching positioning methods still have the following defects: (1) It is difficult to balance the matching accuracy and computational complexity; (2) Aeromagnetic matching positioning is prone to matching failures due to similar results in a wide area. Summary of the Invention

[0004] The present invention aims to solve the problems that it is difficult to balance the matching accuracy and computational complexity of the existing aeromagnetic matching positioning method and that aeromagnetic matching positioning is prone to matching failures due to similar results in a wide area, and further proposes a lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environments.

[0005] The technical solution adopted by the present invention to solve the above problems is as follows: The present invention includes the following steps:

[0006] Step 1: Synchronously collect aeromagnetic observation values M during the flight of the UAV in a specific area i Construct an aeromagnetic anomaly map;

[0007] Step 2: Grid and downsample the aeromagnetic anomaly map and the flight trajectory, and fit the nearest flight trajectory point of each grid node to the corresponding node to obtain a flight trajectory slice;

[0008] Step 3: Perform iterative correlation analysis on the aeromagnetic anomaly map and the flight trajectory slice until the positioning resolution of the analysis result meets the requirements or the positioning iteration times exceed the preset threshold to obtain the UAV aeromagnetic positioning result.

[0009] Further, Step 1 specifically includes:

[0010] Obtain the aeromagnetic observation values M synchronously collected during the flight of the UAV in a specific area iand magnetic measurement related data, based on the collected aeromagnetic observation value M i and magnetic measurement related data, calculate the aeromagnetic anomaly value at the observation position, and construct an aeromagnetic anomaly map. Among them, the magnetic measurement related data includes the Earth's background magnetic field M 0 , magnetic diurnal variation correction M 1 and magnetometer calibration M 2 ;

[0011] The calculation formula for the aeromagnetic anomaly value at the observation position is:

[0012] δM i = M i - M 0 - M 1 - M 2 (1);

[0013] In formula (1), δM i is the aeromagnetic anomaly value at the observation position.

[0014] Furthermore, the gridification of the aeromagnetic anomaly map and the flight trajectory in step 2 specifically includes:

[0015] Step 2.1: Mark the positions without the starting and ending points of the aircraft flight on the aeromagnetic anomaly map, and construct a quadrilateral grid matrix on the aeromagnetic anomaly map based on the grid distance n. Among them, during the construction process, the starting and ending points of the flight form the diagonal points of the quadrilateral matrix;

[0016] Step 2.2: Take the aeromagnetic observation value M i as the aeromagnetic anomaly field intensity of the corresponding point, and the weight coefficients of all points. Calculate the magnetic anomaly field intensity based on the aeromagnetic anomaly field intensity and weight coefficients of the corresponding point

[0017] Step 2.3: Take the corresponding point as a sample point, and use the variogram model to obtain the estimated value of the semi-variance ε i0 between the sample point and the point to be interpolated. On the basis of step 2.2, and under the condition, apply the Lagrange multiplier method to solve the objective function, obtain the Kriging interpolation coefficient, and visualize the interpolation result to complete the gridification of the aeromagnetic anomaly map;

[0018] Step 2.4: Set the flight trajectory as P = {(x i , y i )|i = 1, 2,..., n}, where x i and y i represent the longitude and latitude coordinates of the i-th point, and the resolution of the flight trajectory is the same as that of the aeromagnetic anomaly map;

[0019] Step 2.5: Construct a grid space consistent with the aeromagnetic anomaly using the same grid distance n as the aeromagnetic anomaly map Figure 1 and the aeromagnetic anomaly

[0020] The calculation formula for the aeromagnetic anomaly field intensity is as follows:

[0021]

[0022] In formula (2), λ i is the weight coefficient of the corresponding point;

[0023] The calculation formula for the weight coefficient is:

[0024]

[0025] In formula (3), J is the optimal coefficient set that minimizes the difference between the estimated aeromagnetic anomaly field intensity and the true value M 0 , and ω i is the Kriging interpolation coefficient;

[0026] The calculation formula for the Kriging interpolation coefficient is:

[0027]

[0028] In formula (4), ε ij (i,j = 1,2,...,n) is the semivariance between sample points, and [ω 1 ω 2 … ω n - φ] is the Kriging interpolation coefficient, and ε i0 (i = 1,2,...,n) is the semivariance between the sample point and the point to be interpolated.

[0029] Furthermore, the steps of downsampling the aeromagnetic anomaly map and the flight trajectory in step 2 include:

[0030] For any flight trajectory point p i =(x i ,y i ), find the nearest grid node g j =(x j ,y j ), and fit the nearest flight trajectory point of each grid node to the corresponding node to obtain a flight trajectory slice;

[0031] The expression for the nearest grid node of the flight trajectory point is:

[0032]

[0033] The expression for fitting the most recent flight trajectory points to the corresponding nodes is:

[0034]

[0035] Furthermore, step 3 specifically includes:

[0036] Step 3.1: Perform a convolution operation on the aeromagnetic anomaly map using a CNN network. Set the k-th convolutional kernel in the CNN network as W k , and obtain the index set p of the non-zero positions in the convolutional kernel k ;

[0037] Step 3.2: Calculate the mean absolute difference between the points on the aeromagnetic anomaly map and the corresponding flight trajectory slices at the non-zero positions of the convolutional kernel

[0038] Step 3.3: Take any area of the aeromagnetic anomaly map as a slice of the aeromagnetic anomaly field intensity, and calculate the gradient intensity of the corresponding points in the slice of the aeromagnetic anomaly field intensity. Among them, the value of the point in the slice of the aeromagnetic anomaly field intensity is M n , M n represents the aeromagnetic anomaly field intensity at each position along the flight trajectory, and the value of the point in the slice of the aeromagnetic anomaly gradient intensity is δM n , δM n represents the aeromagnetic anomaly gradient intensity at each position along the flight trajectory;

[0039] Step 3.4: Regard the slice of the aeromagnetic anomaly field intensity and the slice of the aeromagnetic anomaly gradient intensity as convolutional kernels, and perform pyramid iterative matching analysis on the slice of the aeromagnetic anomaly field intensity and the slice of the aeromagnetic anomaly gradient intensity with the corresponding feature maps until the positioning resolution of the analysis result meets the requirements or the positioning iteration times exceed the preset threshold to obtain the UAV aeromagnetic positioning result;

[0040] The index set p of the non-zero positions in the convolutional kernel k is calculated by the formula:

[0041] p k = {(m,n)|W k ≠ 0} (7);

[0042] In formula (7), |p k | is the number of elements in the set, that is, the number of non-zero elements, and (m,n) is the point on the aeromagnetic anomaly map;

[0043] The mean absolute difference is calculated by the formula:

[0044]

[0045] In formula (8), X is the input feature map, and (i,j) is the point on the flight trajectory slice;

[0046] Midpoint δM of the aeromagnetic anomaly gradient intensity slice n The calculation formula is as follows:

[0047] dif(·) = δM n = M n+1 - M n (9).

[0048] Furthermore, step 3.4 specifically includes:

[0049] Step 3.4.1: Take the aeromagnetic anomaly map M1 as the input, conduct a correlation analysis of the aeromagnetic anomaly field intensity with the flight trajectory slice T1, calculate the MAD value of each point according to formula (8), and select 5 - 10 points with the highest correlation as the candidate result P1;

[0050] Step 3.4.2: Obtain the corresponding candidate region R1 according to the candidate result P1, conduct a correlation analysis of the aeromagnetic anomaly gradient intensity and the flight trajectory slice T1 for the candidate region R1, calculate the MAD value of each point according to formula (8), and select the point with the highest correlation and its corresponding map as the input M2 for the next iteration;

[0051] Step 3.4.3: Repeat steps 3.4.1 - 3.4.2, take the point with the highest correlation and its corresponding map output in each iteration as the input for the next iteration until the positioning resolution of the output result of this iteration meets the requirements or the positioning iteration times exceed the preset threshold to obtain the UAV aeromagnetic positioning result;

[0052] The calculation formula for the candidate result P1 is:

[0053]

[0054] The calculation formula for the candidate region R1 is:

[0055] R1(u, v) = M1(x j - 5δx + u, y j - 5δy + v)(11);

[0056] In formula (11), u ∈ [0, 10δx], v ∈ [0, 10δy], (x j , y j ) is the point in P1, δx is the length of the flight trajectory T1 in the east - west direction, and δy is the length of the flight trajectory T1 in the north - south direction;

[0057] The calculation formula for the input M2 of the next iteration is:

[0058]

[0059] The beneficial effects of the present invention are as follows:

[0060] (1) The multi-scale aeromagnetic matching positioning method based on pyramid iteration proposed by the present invention performs matching positioning on maps with gradually increasing resolution. Compared with directly performing matching on high-resolution maps, while ensuring the positioning accuracy, it reduces the computational complexity and overcomes the problem that it is difficult to balance the matching accuracy and computational complexity in existing aeromagnetic matching positioning methods.

[0061] (2) The multi-factor correlation analysis matching method proposed by the present invention adds aeromagnetic anomaly gradient matching on the basis of conventional aeromagnetic anomaly matching. By screening similar results of aeromagnetic anomalies through the added factors, it improves the success rate of matching positioning and overcomes the problem that existing aeromagnetic matching positioning is prone to matching failure due to similar results within a wide area. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a schematic flowchart of the lightweight high-precision unmanned aerial vehicle aeromagnetic positioning method provided by the present invention for GNSS-denied environments;

[0063] Figure 2 It is a schematic diagram of the nearest neighbor interpolation fitting flight trajectory provided by the present invention;

[0064] Figure 3 It is a schematic diagram of the correlation analysis principle of the aeromagnetic anomaly map and the flight trajectory slice provided by the present invention;

[0065] Figure 4 It is a schematic diagram of the matching error of the deviation pair magnetic anomaly matching positioning method provided by the present invention;

[0066] Figure 5 It is a schematic diagram of the aeromagnetic anomaly field strength slice and the aeromagnetic anomaly gradient strength slice provided by the present invention;

[0067] Figure 6 It is a flow block diagram of the pyramid iteration matching method of correlation analysis provided by the present invention;

[0068] Figure 7 It is a schematic diagram of the aeromagnetic anomaly map and the flight trajectory used in the simulation experiment provided by the present invention, Figure 7 where (a) is the aeromagnetic anomaly map and (b) is the flight trajectory map;

[0069] Figure 8 It is a schematic diagram of the aeromagnetic anomaly field strength and the aeromagnetic anomaly field gradient of the flight trajectory with a resolution of 10 meters obtained from the simulation experiment provided by the present invention, Figure 8 where (a) is the aeromagnetic anomaly field strength of slice 1, (b) is the aeromagnetic anomaly gradient of slice 1, (c) is the aeromagnetic anomaly field strength of slice 2, and (d) is the aeromagnetic anomaly gradient of slice 2;

[0070] Figure 9 Schematic diagram of the aeromagnetic anomaly positioning result under the 10-meter resolution map obtained from the simulation experiment provided by the present invention;

[0071] Figure 10 Schematic diagram of the aeromagnetic anomaly positioning results under the 5-meter and 3-meter resolutions obtained from the simulation experiment provided by the present invention, Figure 10 in which (a) is the schematic diagram of the aeromagnetic anomaly positioning result under the 5-meter resolution, and (b) is the schematic diagram of the aeromagnetic anomaly positioning result under the 3-meter resolution;

[0072] Figure 11 Schematic diagram of the longitude and latitude positioning results and true values of the aeromagnetic anomalies of the slices provided by the present invention, Figure 11 in which (a) is the map of the latitude positioning result of the aeromagnetic anomaly of the slice, (b) is the enlarged view of a partial area of the map of the latitude positioning result of the aeromagnetic anomaly of the slice, (c) is the map of the longitude positioning result of the aeromagnetic anomaly of the slice, and (d) is the enlarged view of a partial area of the map of the longitude positioning result of the aeromagnetic anomaly of the slice. Detailed implementation manners

[0073] Detailed implementation manner 1: With reference to Figure 1 explain this implementation manner. As Figure 1 shown, the steps of the lightweight high-precision UAV aeromagnetic positioning method for GNSS-denied environments described in this implementation manner include:

[0074] S1: Construct an aeromagnetic anomaly map;

[0075] The aeromagnetic anomaly map is obtained by synchronously collecting aeromagnetic observation values M i when the UAV flies in a specific area. These observation values are combined with calculations including the Earth's background magnetic field M 0 , magnetic diurnal variation correction M 1 and magnetometer calibration M 2 . The calculation method is as follows:

[0076] δM i = M i - M 0 - M 1 - M 2 (1);

[0077] In formula (1), δM i is the aeromagnetic anomaly value at the observation position. The Earth's background magnetic field can be calculated using the WMM model or the IGRF model. The purpose of this implementation manner is to introduce the matching method, so the process of obtaining aeromagnetic anomaly data is omitted. By using the aeromagnetic anomaly field for matching, the work process is simplified. In fact, the aeromagnetic anomaly map can be superimposed on the geomagnetic map, and then the same aeromagnetic matching method can be applied for positioning.

[0078] S2: Grid and downsample the aeromagnetic anomaly map and the UAV flight trajectory;

[0079] The flight of a UAV usually covers a vast area, and for large-scale map positioning in aeromagnetic anomaly map matching without an initial position, a large amount of calculation is required. In most cases, the onboard platform cannot handle such a huge computational load, and it is necessary to reduce the resolution of the magnetic anomaly map to lighten the burden. Therefore, the aeromagnetic anomaly map needs to be gridded at an appropriate resolution.

[0080] S201: The method proposed in the present invention must ensure that the resolutions of the map and the trajectory are the same. During the gridding process of the aeromagnetic anomaly map, first, the positions of the starting point and the ending point need to be determined, and then a grid matrix is constructed based on the grid distance, where the ending point should be the diagonal point of the quadrilateral formed with the starting point. Next, calculate the aeromagnetic anomaly field intensity at each point in the matrix. In the interpolation of the geophysical field, the Kriging interpolation method is usually used, and the estimation method of the magnetic anomaly field intensity at a specific point is as follows:

[0081]

[0082] In formula (2), λ i is the weight coefficient corresponding to the point, and the weight coefficient is the optimal coefficient set that minimizes the difference between the estimated aeromagnetic anomaly field intensity and the true value M 0 and can be expressed as:

[0083]

[0084] In formula (3), J is the optimal coefficient set that minimizes the difference between the estimated aeromagnetic anomaly field intensity and the true value M 0 ;

[0085] S202: Take the corresponding point as a sample point, and use the variogram model to obtain the estimated value of the semi-variance ε i0 between the sample point and the point to be interpolated. On the basis of S201, and under , apply the Lagrange multiplier method to solve the objective function, obtain the Kriging interpolation coefficient, and visualize the interpolation result to complete the gridding of the aeromagnetic anomaly map:

[0086]

[0087] In formula (4), ε ij (i, j = 1, 2,..., n) is the semi-variance between sample points, and the semi-variance ε ij (i, j = 1, 2,..., n) between sample points is known, [ω 1 ω2 …ω n -φ] is the Kriging interpolation coefficient, ε i0 (i = 1, 2, ..., n) is the semi-variance between the sample point and the point to be interpolated, ε i0 It is usually estimated using a variogram model (such as a linear function model, exponential model, spherical model, Gaussian model, or cubic spline model).

[0088] S203: After the aeromagnetic anomaly map is gridded, the trajectory needs to be refitted using the nearest neighbor interpolation method at the same resolution. The main process is as follows Figure 2 shown. Let the flight trajectory be P = {(x i , y i ) | i = 1, 2, ..., n}, where x i and y i represent the latitude and longitude coordinates of the i-th point. A grid space similar to the map is constructed using the same grid distance n During the trajectory fitting process, for any flight trajectory point p i = (x i , y i ), the process of finding the nearest grid node g j = (x j , y j ) can be expressed as:

[0089]

[0090] Similarly, for all grid nodes g j that can be fitted, the nearest flight trajectory point of each grid node is fitted to that node. This process can be expressed as:

[0091]

[0092] S3: Perform iterative correlation analysis on the aeromagnetic anomaly map and the flight trajectory slices to obtain the UAV aeromagnetic positioning result;

[0093] After constructing the aeromagnetic anomaly map and the flight trajectory slices, the approximate position of the flight trajectory can be obtained by analyzing the correlation between the aeromagnetic anomaly map and the trajectory slices, thereby determining the position at each moment. The flow chart of this method is as Figure 3 shown. This method is inspired by the convolution operation between the convolution kernel and the feature map in the convolutional neural network. By comparing the similarity between the convolution kernel and the feature map, when applied to positioning, the position of the convolution kernel on the map can be determined. However, when implementing this method in aeromagnetic positioning, the following process needs to be optimized:

[0094] S301: Locate the non - zero positions in the convolution kernel: Use the CNN network to perform convolution operations on the aeromagnetic anomaly map. In traditional convolution operations, all positions of the convolution kernel are involved in the calculation. However, in this embodiment, only the non - zero positions in the convolution kernel are selected for operation. Set the k - th convolution kernel in the CNN network as W k , and obtain the index set p of the non - zero positions in the convolution kernel k :

[0095] p k ={(m,n)|W k ≠0}(7);

[0096] In formula (7), |p k | is the number of elements in the set, that is, the number of non - zero elements, and (m,n) are the points on the aeromagnetic anomaly map;

[0097] S302: Obtain the correlation: In aeromagnetic positioning, the correlation measurement problem is a one - dimensional data matching process. Whether it is direct data association or feature association, it is essentially a concept of distance measurement and belongs to the category of data association. Commonly used data association criteria include Product Correlation (PROD), Normalized Product Correlation (NPROD), Mean Absolute Difference (MAD), Mean Square Deviation (MSD), and Hausdorff distance measurement. Among them, PROD pays more attention to the matching of data change characteristics, but has lower accuracy in numerical matching (or when there is no scale factor error). The Hausdorff distance measurement is greatly affected by noise, and a single outlier may cause the entire slice matching to fail. Therefore, MAD and MSD are more suitable for the aeromagnetic matching process. The MSD and MAD matching methods can obtain good matching results by minimizing the residual offset in the case of no noise and limited data length. This embodiment selects the MAD matching method, as follows:

[0098]

[0099] In formula (8), X is the input feature map, and (i,j) are the points on the flight trajectory slice;

[0100] S303: Feature selection: In S302, similarity measurement is performed. However, the feature selection for correlation measurement will affect the accuracy and effectiveness of positioning. The commonly used matching feature is the aeromagnetic anomaly field intensity. However, the matching based on the aeromagnetic anomaly field strength may be affected by regions with similar aeromagnetic anomaly values. Aeromagnetic anomalies are affected by the distribution of the earth's mineral resources, terrain shape, and other ground interference factors. Therefore, their distribution is uneven, but changes in an irregular and gradual manner. For example Figure 4As shown, curve A represents the flight path, and curve B represents the matching result based on deviation data. If there are areas with similar aeromagnetic anomaly values, the matching process may show a smaller MAD value at curve B, which means that the position of curve B has a higher correlation, resulting in its alignment to that position. This may cause significant errors in the matching of this part. If there are irregular or inconsistent matching errors in multiple consecutive parts, the matching method may fail, posing a serious risk to the drone. Using the gradient as a matching factor can effectively filter out smaller types of deviation interference.

[0101] Since similar aeromagnetic anomaly values may affect the positioning result, in order to overcome this problem in this embodiment, multi-factor correlation analysis matching is performed by combining the aeromagnetic anomaly field intensity and the aeromagnetic anomaly gradient, as Figure 5 shown. In the aeromagnetic anomaly field intensity slice shown on the left of 5, the value of each point (e.g., from M 1 to M 5 ) represents the aeromagnetic anomaly field intensity at each position along the flight trajectory. In Figure 5 the aeromagnetic anomaly gradient intensity slice shown on the right, the value of each point is defined as:

[0102] dif(·) = δM n = M n+1 - M n (9);

[0103] Similar to a convolutional neural network, these two types of slices can be regarded as convolutional kernels. During the correlation analysis process, each slice is analyzed with its corresponding feature map to obtain the degree of correlation, thereby determining the result of aeromagnetic positioning. In summary, the multi-factor correlation analysis matching method proposed in this embodiment, on the basis of conventional aeromagnetic anomaly matching, adds aeromagnetic anomaly gradient matching, screens similar results of aeromagnetic anomalies through the added factors, improves the success rate of matching positioning, and overcomes the problem that existing aeromagnetic matching positioning is prone to matching failure due to similar results within a wide area.

[0104] S304: Regard the aeromagnetic anomaly field intensity slice and the aeromagnetic anomaly gradient intensity slice as convolutional kernels, and perform pyramid iterative matching analysis on the aeromagnetic anomaly field intensity slice and the aeromagnetic anomaly gradient intensity slice with the corresponding feature maps until the positioning resolution of the analysis result meets the requirements or the positioning iteration times exceed the preset threshold to obtain the drone aeromagnetic positioning result;

[0105] The core of aeromagnetic anomaly matching is to extract a data segment from a larger dataset. However, there is a contradiction between low computational cost and high positioning accuracy. High precision requires a large amount of computation, while reducing computation will lower the precision. Inspired by two-dimensional image matching methods, this paper proposes a pyramid iterative matching method for aeromagnetic positioning. This algorithm balances global search and local positioning by evolving from the traversal process. The process of the pyramid iterative matching method is as Figure 6 shown.

[0106] S30401: Different from the common global scaling and feature extraction methods in two-dimensional image processing, the pyramid iterative matching method proposed in this embodiment processes the aeromagnetic anomaly map at only one scale at a time. After the pre-acquired aeromagnetic anomaly map M1 is input into this method, it will first perform a correlation analysis with the flight trajectory slice (such as slice T1), and this slice is downsampled to 10 meters. The MAD value of each point on the map is used as the basis for correlation evaluation, and the MAD value of each point is calculated according to formula (8). Select the 5 to 10 points with the highest correlation (the smallest value) as the candidate results. The above process can be described as follows:

[0107]

[0108] 30302: The lengths of the flight trajectory T1 in the east-west direction and the north-south direction can be represented as δx and δy respectively. Therefore, for the set P1, the candidate region R1 that can be selected can be represented as:

[0109] R1(u,v) = M1(x j -5δx + u,y j -5δy + v)(11);

[0110] In formula (11), u ∈ [0,10δx], v ∈ [0,10δy], (x j ,y j ) is the point in P1, δx is the length of the flight trajectory T1 in the east-west direction, and δy is the length of the flight trajectory T1 in the north-south direction;

[0111] S30303: Obtain the corresponding candidate region R1 according to the candidate result P1, perform a correlation analysis on the candidate region R1 between the aeromagnetic anomaly gradient intensity and the flight trajectory slice T1, calculate the MAD value of each point according to formula (8), and select the point with the highest correlation and its corresponding map as the input M2 for the next iteration, which can be expressed as:

[0112]

[0113] S30304: Repeat S30301 - S30303. Use the point with the highest correlation in the output of each iteration and its corresponding map as the input for the next iteration until the positioning resolution of the output of this iteration meets the requirements or the number of positioning iterations exceeds the preset threshold, and obtain the UAV aeromagnetic positioning result.

[0114] In summary, the proposed multi - scale aeromagnetic matching positioning method based on pyramid iteration performs matching positioning on maps with gradually increasing resolution. Compared with directly performing matching on high - resolution maps, it ensures positioning accuracy while reducing computational complexity, overcoming the problem that it is difficult to balance matching accuracy and computational complexity in existing aeromagnetic matching positioning methods.

[0115] Specific Embodiment 2: To verify the technical effects of the lightweight and high - precision UAV aeromagnetic positioning method proposed in Specific Embodiment 1 for GNSS - denied environments, the following simulation experiments were conducted in this embodiment:

[0116] 1. Data and experimental settings:

[0117] This specific embodiment is mainly based on the flight simulation software Ardupilot and the aeromagnetic anomaly data and maps collected during the flight of UAVs in the MAGNAV project. The aeromagnetic anomaly map used in the experiment is near Ontario, Canada, covering a latitude range from 44.8°N to 45.9°N and a longitude range from 76.1°W to 77.4°W.

[0118] In this area, a flight simulation is carried out using the fixed - wing UAV simulation program in Ardupilot. The flight trajectory involves multi - direction flights at an altitude of about 500 meters above the ground, with each direction covering more than 2000 meters to ensure there is sufficient data to support the aeromagnetic matching and positioning method proposed in this paper. After obtaining the flight trajectory, the aeromagnetic anomaly map is interpolated according to the coordinates of each position to obtain the aeromagnetic anomaly field strength at these points. Since the focus of this invention is on the implementation of the aeromagnetic anomaly matching and positioning method, the process of using an airborne magnetometer and IMU for aeromagnetic compensation to obtain the aeromagnetic anomaly field at each moment during the UAV flight is omitted. The aeromagnetic anomaly map and flight trajectory used in the experiment are as shown in Figure 7 (a) and Figure 7 (b). According to the data provided by the map, some areas are obtained by extrapolation. Due to its inaccuracy, it is not recommended to use the data in these areas, and the data in these areas is not selected in this experiment.

[0119] According to the size of the experimental space, in the first round of experiments, the resolution of the aeromagnetic anomaly map and the trajectory was set to 10 meters. To ensure that the matching trajectory is long enough for point-by-point matching on the map, the path length of each slice in this round was set to 50 points, that is, the distance of each trajectory slice is about 200 meters. In addition, there is no overlap between slices.

[0120] According to the method proposed by the present invention, in each subsequent round, the resolution is doubled, that is, 5 meters and 3 meters (rounded up), and a total of three rounds of experiments are carried out to determine the matching accuracy of each round. At these resolutions, the path length of the slice is the same as that at 10-meter resolution, which means that the total length of each path is the same, but the number of matching points increases. By gradually increasing the resolution through multiple iterations, that is, the pyramid-based iterative matching method, the positioning accuracy can be improved while avoiding the high computational load caused by directly using high-resolution maps and flight trajectory matching.

[0121] 2. Data and experimental settings:

[0122] ① Gridding and downsampling of aeromagnetic anomaly maps and trajectories:

[0123] In the first round of matching, the MIT map was gridded with a resolution of the gridding matrix of 10 meters, which was determined according to the size of the map and the expected computational load. Starting from the first point of the flight trajectory, each point was fitted to the nearest 10-meter resolution grid. Then, for each grid point, the nearest point on the flight trajectory was searched, and the aeromagnetic anomaly field strength at that point was used as the aeromagnetic anomaly field strength of the grid point. This fitting process continued until the values of 50 grid points were obtained, and these values were regarded as the first slice segment. This process was repeated, and the path was segmented. To verify the method proposed in this paper, 20 consecutive slices were selected for matching. Based on these slices, the difference between adjacent points was obtained as the aeromagnetic anomaly field gradient. The magnetic anomaly field strength and magnetic anomaly field gradient of the flight trajectory at 10-meter resolution for slice 1 and slice 2 are as Figure 8 shown.

[0124] ② Cross-correlation analysis of slices and aeromagnetic anomaly grid maps:

[0125] The correlation between the aeromagnetic anomaly field strength slices of 20 flight trajectories and the map was analyzed, and the mean absolute deviation (MAD) of each point was calculated. Five points with the smallest MAD were searched for each flight trajectory. At the positions of these points, the correlation was further analyzed using the aeromagnetic anomaly field gradient slices, and the MAD was calculated again to find the point with the smallest MAD. The results on the 10-meter resolution map are as Figure 9 shown.

[0126] At Figure 9In it, the horizontal axis represents the slice number, and the vertical axis represents the MAD. For each slice, the five points with the smallest MAD are selected and marked with squares. Then, the aeromagnetic anomaly field gradient is used for further screening, and the screening results are marked with dots. From the matching results of the 20 slices in the figure, it can be seen that 75% of the slices are successfully matched, and the best matching points of most slices are located at the positions with the highest correlation of the aeromagnetic anomaly field. The remaining points can also be matched through the aeromagnetic anomaly field gradient. However, slices 6, 11, 12, 13, and 15 failed to be successfully matched, mainly because these segments are in the turning stage and the area covered by the trajectory is limited. From the perspective of features, these segments have fewer significant features compared with other segments, resulting in matching failures.

[0127] ③ Pyramid iterative matching:

[0128] Repeat the above experiment. Centered on the position of the successfully matched points, extract a 2000-meter map within a range approximately ten times the track length and grid it at a resolution of 5 meters, and further match it with the track slices. Based on the 5-meter matching results, further grid the map at a resolution of 3 meters within a range approximately five times the track length and extract a 1000-meter map. Finally, high-precision positioning is achieved. These two resolutions are respectively half and one-fourth of the 10-meter resolution, that is, the resolutions are increased by two times and four times respectively. The aeromagnetic anomaly positioning results at 5 meters and 3 meters resolutions are as Figure 10 shown.

[0129] From Figure 10 it can be seen that except for slices 6, 11, 12, 13, and 15 that failed to be matched at 10-meter resolution, the MAD of the matching results of other slices is closer to 0, which means better matching performance. In addition, the positions of the best matching results of the aeromagnetic anomaly field strength and the aeromagnetic anomaly field gradient of almost all slices are almost the same, which also proves the uniqueness of the matching results within a small area.

[0130] Based on the above experimental results, this embodiment plots the true position of the UAV flight trajectory, the matching results at 10-meter resolution, the matching results at 5-meter resolution, and the matching results at 3-meter resolution in the same figure, as Figure 11 (a)- Figure 11 (d) shown.

[0131] To display the positioning results, this embodiment calculates the MAD of the Euclidean distance between the best matching point and the true position, as shown in Table 1. The calculation formula of MAD is as follows:

[0132]

[0133] In formula (13), n represents the number of points after downsampling of the trajectory segment, and d represents the Euclidean distance between the actual position of each point in the trajectory and the corresponding point in the matching result. The smaller the MAD, the better the matching performance.

[0134] Table 1

[0135]

[0136]

[0137] As can be seen from Table 1, the positioning accuracy is the highest at a resolution of 3 meters, and the positioning accuracy of the vast majority of slices has improved during the matching process from a resolution of 10 meters to 3 meters. This shows that the method proposed in this paper significantly improves the positioning accuracy through multiple iterations, and also proves that there are no interfering regions in the map with similar aeromagnetic anomaly field strengths, such as Figure 9 and Figure 10 the several candidate regions shown.

[0138] Currently, based on other types of sensors in a GNSS-denied environment, at a similar ground height (AGL), their accuracy is usually about 1 meter, and the least accurate is within 10 meters. The method proposed in this paper can achieve a comparable positioning accuracy when the data is accurate enough and has the potential for practical application on drones, especially when GNSS signals are interfered.

[0139] In terms of computational optimization, a statistical analysis was performed on the average computational load of all slices during the three rounds of positioning processes and compared with directly using a 3-meter resolution map for positioning. In this paper, the number of floating-point operations (FLOPs) is used as a metric for computational complexity. The smaller the FLOPs, the lower the complexity, and it is more suitable for low-power platforms with lower performance. The results are shown in Table 2.

[0140] Table 2

[0141]

[0142]

[0143] From the above experimental results, it can be seen that using the method proposed by the present invention, the average computational load of a single slice is 20.08 GFLOPs, while the computational load of directly using a 3-meter resolution map for matching and positioning is 1.28 TFLOPs. Compared with directly using high-resolution slices for large-scale map positioning, the method proposed by the present invention reduces the computational load by more than 98%. It should be noted that the map size and flight trajectory slice length used for positioning by the present invention are the same at 5-meter and 3-meter resolutions. Therefore, the computational load at 3-meter resolution is significantly higher than that at 5-meter resolution. In theory, using a smaller map and flight trajectory slices for 3-meter resolution positioning can further reduce the computational load. However, in order to compare the positioning accuracy, the sizes of these two parameters are not changed in this paper. This result shows that the method proposed by the present invention can effectively reduce the computational load of matching and positioning.

[0144] The above are only the preferred embodiments of the present invention, and do not impose any form of limitation on the present invention. Although the present invention has been disclosed above with the preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to the above-disclosed technical content to form equivalent embodiments with equivalent changes within the scope of the technical solution of the present invention. However, as long as it does not depart from the content of the technical solution of the present invention and is based on the technical essence of the present invention, any simple modification, equivalent replacement, and improvement made to the above embodiments still fall within the protection scope of the technical solution of the present invention.

Claims

1. A lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environments, characterized in that: The steps of the method for lightweight and high-precision aeromagnetic positioning of unmanned aerial vehicles in a GNSS-denied environment include: Step 1: Collect aeromagnetic observations M synchronously while the drone is flying over a specific area i Construct aeromagnetic anomaly maps; Step 2: Grid and downsample the aeromagnetic anomaly map and the UAV flight trajectory, and fit the nearest flight trajectory point of each grid node to the corresponding node to obtain the flight trajectory slice; Step 3: Perform iterative correlation analysis on the aeromagnetic anomaly map and flight trajectory slices until the positioning resolution of the analysis result meets the requirements or the number of positioning iterations exceeds the preset threshold, and obtain the UAV aeromagnetic positioning result.

2. The lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environments according to claim 1, characterized in that: Step 1 specifically includes: Get the aeromagnetic observation value M collected synchronously when the UAV flies in a specific area i As well as magnetic survey related data, based on the collected aeromagnetic observations M i The aeromagnetic anomaly value of the observation position is calculated from the magnetic survey related data, and the aeromagnetic anomaly map is constructed. The magnetic survey related data include the earth's background magnetic field M0, magnetic diurnal variation correction M1 and magnetometer calibration M2; The calculation formula of the aeromagnetic anomaly value at the observation location is: δM i =M i -M0-M1-M2(1); In formula (1), δM i is the aeromagnetic anomaly value of the observation location.

3. The lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environments according to claim 1, characterized in that: The gridding of the aeromagnetic anomaly map and flight trajectory in step 2 specifically includes: Step 2.1: Mark the locations of the flight start and end points without aircraft in the aeromagnetic anomaly map, and construct a quadrilateral grid matrix in the aeromagnetic anomaly map based on the grid distance n, wherein the flight start and the flight end point in the grid matrix form the diagonal points of the quadrilateral matrix during the construction process; Step 2.2: Convert the aeromagnetic observation value M i As the aeromagnetic anomaly field intensity of the corresponding point and the weight coefficient of all points, the magnetic anomaly field intensity is calculated based on the aeromagnetic anomaly field intensity of the corresponding point and the weight coefficient Step 2.3: Take the corresponding points as sample points and use the variogram model to obtain the semivariance ε between the sample points and the points to be interpolated i0 The estimated value of The Lagrange multiplier method is used to solve the objective function, obtain the Kriging interpolation coefficients, and visualize the interpolation results to complete the gridding of the aeromagnetic anomaly map. Step 2.4: Set the flight trajectory to P = {(x i ,y i )|i=1,2,...,n}, where x i and i Represents the longitude and latitude coordinates of the i-th point, and the flight trajectory and resolution are the same as the resolution of the aeromagnetic anomaly map; Step 2.5: Use the same grid distance n as the aeromagnetic anomaly map to construct a grid space consistent with the aeromagnetic anomaly map Magnetic anomaly field strength The calculation formula is: In formula (2), λ i is the weight coefficient of the corresponding point; The calculation formula of the weight coefficient is: In formula (3), J is the estimated aeromagnetic anomaly field strength The optimal coefficient set that minimizes the difference with the true value M0, ω i is the Kriging interpolation coefficient; The calculation formula of Kriging interpolation coefficient is: In formula (4), ε ij (i,j=1,2,...,n) is the semivariance between sample points, [ω1ω2…ω n -φ] is the Kriging interpolation coefficient, ε i0 (i=1,2,...,n) is the semivariance between the sample point and the point to be interpolated.

4. The lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environment according to claim 1, characterized in that: The steps of downsampling the aeromagnetic anomaly map and flight trajectory in step 2 include: For any flight trajectory point p i =(x i ,y i ), find the nearest grid node g j =(x j ,y j ), fit the nearest flight trajectory point of each grid node to the corresponding node to obtain the flight trajectory slice; The expression of the nearest grid node of a flight trajectory point is: The expression for fitting the nearest flight trajectory point to the corresponding node is:

5. The lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environment according to claim 1, characterized in that: Step 3 specifically includes: Step 3.1: Use the CNN network to perform convolution operations on the aeromagnetic anomaly map, and set the kth convolution kernel in the CNN network to W. k , get the index set p of the non-zero positions in the convolution kernel k ; Step 3.2: At the non-zero position of the convolution kernel, calculate the mean absolute difference between the point of the aeromagnetic anomaly map and the corresponding flight trajectory slice Step 3.3: Take any area of ​​the aeromagnetic anomaly map as an aeromagnetic anomaly field intensity slice, and calculate the gradient intensity of the corresponding point in the aeromagnetic anomaly field intensity slice, where the value of the midpoint of the aeromagnetic anomaly field intensity slice is M n , M n It represents the aeromagnetic anomaly field intensity at each position along the flight track. The value of the midpoint of the aeromagnetic anomaly gradient intensity slice is δM n , δM n It indicates the aeromagnetic anomaly gradient intensity at each location along the flight track; Step 3.4: The aeromagnetic anomaly field intensity slice and the aeromagnetic anomaly gradient intensity slice are regarded as convolution kernels, and the aeromagnetic anomaly field intensity slice and the aeromagnetic anomaly gradient intensity slice are subjected to pyramid iterative matching analysis with the corresponding feature map until the positioning resolution of the analysis result meets the requirements or the number of positioning iterations exceeds the preset threshold, and the UAV aeromagnetic positioning result is obtained; The index set p of the non-zero positions in the convolution kernel k The calculation formula is: p k ={(m,n)|W k ≠0} (7); In formula (7), |p k | is the number of elements in the set, that is, the number of non-zero elements, and (m,n) is the point on the aeromagnetic anomaly map; Mean absolute difference The calculation formula is: In formula (8), X is the input feature map, (i, j) is the point on the flight trajectory slice; Aeromagnetic anomaly gradient intensity slice midpoint δM n The calculation formula is: dif(·)=δM n =M n+1 -M n (9)。 6. The lightweight and high-precision UAV aeromagnetic positioning method for GNSS-denied environments according to claim 5, characterized in that: Step 3.4 specifically includes: Step 3.4.1: Take the aeromagnetic anomaly map M1 as input, and perform correlation analysis on the aeromagnetic anomaly field intensity with the flight track slice T1. Calculate the MAD value of each point according to formula (8), and select the 5 to 10 points with the highest correlation as candidate results P1; Step 3.4.2: Obtain the corresponding candidate region R1 according to the candidate result P1, perform correlation analysis between the aeromagnetic anomaly gradient intensity and the flight track slice T1 on the candidate region R1, calculate the MAD value of each point according to formula (8), and select the point with the highest correlation and its corresponding map as the input M2 for the next iteration; Step 3.4.3: Repeat steps 3.4.1-3.4.2, and use the most correlated point and its corresponding map output from each iteration as the input for the next iteration, until the positioning resolution of the output result of this iteration meets the requirements or the number of positioning iterations exceeds the preset threshold, and the UAV aeromagnetic positioning result is obtained; The calculation formula of candidate result P1 is: The calculation formula of candidate region R1 is: R1(u,v)=M1(x j -5δx+u,y j -5δy+v)(11); In formula (11), u∈[0,10δx], v∈[0,10δy], (x j ,y j ) is a point in P1, δx is the length of the flight trajectory T1 in the east-west direction, and δy is the length of the flight trajectory T1 in the north-south direction; The calculation formula for the input M2 of the next iteration is:

Images