Abstract
A high-precision marine gravity gradient reference map is key to enabling underwater gravity gradient matching navigation. At present, the construction of the reference maps is based on quadrilateral geographic grids. However, quadrilateral grids lead to detriangulation at high latitudes, which limits the global applicability of such maps to underwater gravity navigation. To circumvent the limitations of quadrilateral grids, a hexagonal grid is introduced for constructing the reference map. This paper analyzes the characteristics of the icosahedral Snyder equal area aperture 4 hexagon (ISEA4H) and H3 grid systems and selects an appropriate grid system. In addition, we calculate and analyze the grid model errors and matching positioning errors of hexagonal and quadrilateral grids at the same resolution. The experimental results show that the grid model and matching positioning errors of a hexagonal grid system are more than 14% and 15% less than those of a quadrilateral grid system, respectively, indicating the feasibility and effectiveness of applying hexagonal grids to gravity gradient matching navigation. Given the low construction efficiency of a marine hexagonal grid gravity gradient reference map, we propose an efficient CPU+GPU hybrid parallel scheme. A global total tensor marine hexagonal grid gravity gradient reference map model is then constructed.
1. Introduction
Given China’s strategy to become a maritime power and the proposed 14th five-year plan for future ocean and deep-sea tasks, underwater vehicles have become a development focus of the national marine strategy, which plays an indispensable role in both military and civil fields in China. How to realize the long-term, autonomous, and high-precision positioning and navigation of submarines has become the research focus of multiple scholars. Inertial navigation systems (INSs) are often used in the passive positioning and navigation of underwater vehicles; however, the navigation error of INS systems cannot be independently eliminated and gradually accumulates with time. Therefore, the navigation accuracy of passive positioning and navigation using only an INS system is greatly reduced over long periods of time. This has spurred the development of underwater gravity matching navigation technology. The goal of such technology is to use the gravity sensor carried by underwater vehicles to perform real-time gravity measurements and then to match these measurements with a high-precision gravity map using a relevant algorithm, so as to correct the positioning information and component errors of the INS system. As the second-order differential of the gravity potential, measurements of the gravity gradient are not affected by the acceleration of the underwater vehicle carrier. The gravity gradient has five independent observational components. Therefore, compared with gravity anomaly matching navigation, gravity gradient matching navigation can obtain higher precision positioning and navigation information [1, 2].
The construction of a marine gravity gradient reference map is a prerequisite for the realization of gravity gradient-assisted INS navigation; such a map would be stored in a geographic grid composed of longitude and latitude in the form of discrete points on the basis of a geographic coordinate system [3, 4]. This type of geographic grid has a simple distribution and can easily be used to determine the grid size and index location. A geographic grid system can effectively map and facilitate data management for all types of data structures. However, geographical grids have many limitations because they are closely related to longitude and latitude. For example, the area and shape of geographical grids deform at higher latitudes and gradually degenerate close to the poles. Moreover, the geographical grid structure does not have adjacent consistency and the distance between any quadrilateral grid element with common edges or angles is not equal, which increases the difficulty of spatial analyses. To effectively solve these problems, this paper introduces a hexagonal grid based on polyhedron projection to construct a marine gravity gradient reference map. Compared with geographical grids, hexagonal grids have the following advantages. (1) Hexagonal grids have a good characteristic of equal area. The grid area of a hexagonal grid does not decrease with increasing latitude, and the grid will not degrade in the polar region. (2) Hexagonal grids have good adjacent consistency, and the distance between adjacent grids is relatively fixed. This effectively suppresses errors caused by the inconsistent grid distances of different matching points in the process of gravity matching navigation. (3) Compared with quadrilateral grids, hexagonal grids have greater angular resolution and the particle sense of the discretized boundary is weak. When using a reference map with the same or similar resolution for gravity matching navigation, hexagonal grids provide a higher matching accuracy for tracks. (4) Hexagonal grids are closer to circles. When calculating the disturbing gradient, the influence of the surrounding flow grid points on the center point is the same. Therefore, compared with quadrilateral grids, hexagonal grids can weaken errors caused by asymmetry.
Based on the above advantages, this paper first explores and analyzes the impact of quadrilateral and hexagonal grid geometry errors on matched navigation and positioning. Then, hexagons are selected as the basic grid construction units of the marine gravity gradient reference map. Based on Stokes’ theory, the calculation of the actual disturbing gravity gradient is divided into two parts: the model disturbing gravity gradient and the residual disturbing gravity gradient . Finally, according to the characteristics of the hexagonal and in their calculation processes and using a central processing unit (CPU) parallel method based on OpenMP and a graphics processing unit (GPU) parallel method based on CUDA, a hybrid CPU+GPU parallel programming scheme is designed to enable the efficient construction of a hexagonal marine gravity gradient reference map.
2. Materials and Methods
At present, the most commonly used hexagonal grid systems include the icosahedral Snyder equal area aperture 4 hexagon (ISEA4H) and H3 geospatial indexing systems.
2.1. The ISEA4H System
To understand the ISEA4H system, we start with an icosahedron. An icosahedron has the largest number of faces of the five types of regular polyhedrons (the others being the regular tetrahedron, regular hexahedron, regular octahedron, and regular dodecahedron). Therefore, when it is projected onto a sphere, the deformation of each face element is the smallest. When using an icosahedron to project the Earth, it is also necessary to consider the relative positional relationship between the icosahedron and the Earth. A common positioning method is to align two vertices of the icosahedron with the North and South Poles of the Earth, such that the edge passing through the North Pole is coplanar with the central meridian (0° longitude), as shown in Figure 1.
In Figure 1, point P coincides with the North Pole of the Earth, point Q coincides with the South Pole, and the side PC is in the same plane as the central meridian. According to the properties of an icosahedron, points A, B, C, D, and E are on the same latitude circle, with the latitude being approximately 26° 34 and the longitude interval of each point being 72°.
Expanding the regular icosahedron in Figure 1 into a plane results in the plane area shown in Figure 2(b), which can be used to establish the corresponding positional relationship between the surface of the Earth and the map plane via a projection using the equal area polyhedron orthographic projection method proposed by the American map projection expert Schneider. First, based on the plane triangles in Figure 2(b), the triangular areas are segmented by a hexagonal grid, as shown in Figure 2(c). Then, the Schneider inverse projection transformation is used and the regular icosahedron segmented by hexagons is projected onto the sphere. Finally, the ISEA4H grid is generated. The number of grid cells at each level of the grid system and the grid information are given in Table 1.
2.2. The H3 System
The H3 geospatial index system is a discrete global grid system that is composed of hexagonal tiles with hierarchical indexes [5, 6]. Its positioning method differs from that of ISEA4H; however, it uses a map positioning method similar to that of Dymaxion [7–9] proposed by Fuller. The H3 system gives priority to ensuring that planes expanded in land regions are broken as little as possible. Therefore, a vertex is not located at the pole but at (5.2454°E, 2.3009°N) and the azimuth of the adjacent vertices is fixed at 7.46658°. The H3 system is composed of 16 levels of hexagonal grids. The number of grid cells in each level of the grid system and the grid information are given in Table 2.
2.3. Hexagonal Grid System Selection
Even though the ISEA4H and H3 grid systems both use hexagonal grids as their units, there are some differences between them. (1)The ISEA4H grid system adopts the four-aperture division method; that is, the number of grids between adjacent levels increases by four times (as shown in Figure 3(a)). The H3 system adopts the seven-aperture division method; that is, the number of grids between adjacent levels increases by seven times (as shown in Figure 3(b))(2)The ISEA4H grid system adopts the Schneider equal area projection, and the area of the hexagonal elements in the grid system is basically the same after splitting. To ensure rapid indexing between levels, the H3 system sacrifices the equal area, and its hexagonal grid unit area is not consistent(3)The ISEA4H grid system adopts a regular icosahedral projection. The two vertices of the icosahedron are located at the North and South Poles, and the edge passing through the North Pole is coplanar with the central meridian. Therefore, the ISEA4H grid system is relatively regular and has a stronger correlation with the South and North Poles of the Earth and longitude and latitude. The H3 system adopts a positioning similar to a map. To minimize the splitting of land regions when expanded into a plane, the positioning point deviates from the pole(4)In terms of the coding index and management at all levels, the H3 system is more flexible than the ISEA4H system and can quickly index grids at different levels
(a) Four apertures
(b) Seven apertures
Based on the above comparison of the advantages and disadvantages of the two grid systems and the actual requirements of marine gravity gradient matching navigation, we selected the ISEA4H grid system as the basis to construct the marine gravity gradient reference map. The main reasons are as follows. (1) The ISEA4H grid system has a good characteristic of equal area, and the grid size and shape are relatively consistent; this is conducive to the solution process of the gravity matching navigation algorithm. Conversely, the H3 system grid has a poor characteristic of equal area; therefore, the error caused by the inconsistent grid size needs to be compensated for in the solution process of the gravity matching algorithm. (2) Even though the H3 system has less segmentation in the land regions, this does not have a substantive effect on the construction of a marine gravity gradient reference map and projection positioning of the H3 grid has little correlation with the geographical poles and the central meridian, which is not conducive to solving the gravity matching navigation algorithm. (3) The ISEA4H grid system is a four-aperture grid system; that is, the number of grids at each level increases by four times and the average grid distance is reduced to 1/2. Meanwhile, the H3 system adopts a seven-aperture grid system; that is, the number of grids at each level is increased by seven times and the average grid distance is reduced by 0.378 times. The resolutions of the grids at all levels of the H3 system therefore vary greatly. For example, on the sixth level of the H3 system, the average side length of the grid is 3.23 km and the grid area is 36.13 km2. The number of grids per unit area is equivalent to that of a 3 20 geographic grid, but the number of grids per unit area for the next level, that is, the seventh level, is equivalent to that of a 75 geographic grid. Because of the large variations between the levels, the corresponding resolution does not actually match the quadrilateral grid resolution used in navigation. The 12th level of the ISEA4H grid system has a grid spacing of 1.884 km and a grid area of 3.040 km2. The number of grids per unit area is therefore equivalent to the 1-resolution geographic grid currently used in gravity matching navigation, making the ISEA4H grid system a more suitable grid division method for the construction of a gravity gradient reference map.
2.4. Error Analysis of the Grid Geometry for Matching Navigation and Positioning
Here, we analyze the effectiveness of a reference map based on the hexagonal grid unit for gravity matching navigation. Then, without considering the measurement error of the gravity/gravity gradiometer, the measurement error of the inertial components, or errors in the gravity/gravity gradient reference map, we examine the error caused by the geometric structure of the reference map grid in the process of discretizing the track. Using the inertial navigation solution program and the gravity gradient matching navigation program, we calculate and analyze the hexagonal grid model error and the grid matching positioning error, respectively. The experimental area is a marine area with range of (from 19°N to 29°N and from 41°W to 51°W). The inertial navigation track simulation program is used to simulate the track and measurement data of an underwater vehicle. The inertial navigation sampling rate is 100 Hz. The simulation data include the output of the inertial components (the angular and velocity increments), the gravity gradient observation data, and the real track of the underwater vehicle. The error settings of the inertial components are listed in Table 3.
We used the simulation program to generate two simulated tracks as experimental objects; the navigation speed was set to 10 nautical miles per hour, the navigation time of track 1 was 93,453 s (approximately 25.96 h), and the navigation time of track 2 was 90,984 s (approximately 25.27 h).
2.4.1. Grid Model Error Analysis
The grid model error is the distance from the grid midpoint to the inertial navigation track when the discrete grid is used to approach the inertial navigation track. As shown in Figure 4, when the track is determined, this error is related solely to the grid model. It is the inherent error that occurs when using a grid to generate the track and cannot be compensated for by the algorithm.
In Figure 4, the dark blue line represents the inertial navigation track and the light blue grid represents the grid through which the track passes. These grids can be used as a discrete representation of the track. The distance from the midpoint of the grid to the track is the grid model error. Obviously, a smaller grid model error corresponds to the simulation of the discrete grid sequence being closer to the track and the track restoration degree being higher.
For the two simulated tracks, the 11th, 12th, and 13th levels of the ISEA4H hexagonal grid system were used for the calculations. According to the number of hexagonal grids per unit area and the midpoint distance of the adjacent grids at the different levels, the equivalent quadrilateral grid was selected for comparison. The error of the hexagonal grid model was then evaluated via an analysis of the calculation results. The correspondence between the grids is shown in Table 4.
Table 4 indicates that the 11th level of the hexagonal grid roughly corresponds to a geographical quadrilateral grid with 2 resolution; the 12th level corresponds to a 1-resolution geographic quadrilateral grid; and the 13th level corresponds to a 30-resolution geographic quadrilateral grid. The inertial navigation solutions were found for the two tracks, and on this basis, the grid model errors of the hexagonal and quadrilateral grids at different levels and resolutions were calculated, compared, and analyzed. The calculation results are shown in Table 5.
Table 5 shows that, compared with those of the corresponding resolution quadrilateral grid, the grid model errors of the hexagonal grid system at different levels are reduced to a certain extent.
For track 1, compared with that of the 2-resolution quadrilateral grid, the average grid model error of the 11th level of the hexagonal grid is reduced by 24.87% and the standard deviation is reduced by 20.12%. Compared with that of the 1-resolution quadrilateral grid, the average grid model error of the 12th level of the hexagonal grid is reduced by 20.00% and the standard deviation is reduced by 20.62%. Compared with that of the 30-resolution quadrilateral grid, the average grid model error for the 13th level of the hexagonal grid is reduced by 20.41% and the standard deviation is reduced by 20.54%.
For track 2, compared with that of the 2-resolution quadrilateral grid, the average grid model error of the 11th level of the hexagonal is reduced by 14.01% and the standard deviation is reduced by 14.25%. Compared with the 1-resolution quadrilateral grid, the average grid model error of the 23th level of the hexagonal grid is reduced by 14.54% and the standard deviation is reduced by 14.86%. Compared with the 30-resolution quadrilateral grid, the average grid model error of the 13th level of the hexagonal grid is reduced by 14.02% and the standard deviation is reduced by 14.29%.
For tracks 1 and 2, the adoption of a hexagonal grid reduces the average value and standard deviation of the grid model error by more than 14% and the grid model error accuracy of track 1 is improved by more than 20%. It can be seen therefore that a hexagonal grid has more advantages for track discretization than a quadrilateral geographic grid. In addition, for the same track, the accuracy improvement between each level of the hexagonal grid and its corresponding resolution quadrilateral gird is small, indicating that the grid model error for the same track is related to only the grid selection not the level or resolution.
2.4.2. Grid Matching Positioning Error Analysis
The grid midpoint of the matching point is usually selected as the matching track position when using the gravity gradient matching algorithm for matching navigation. The grid matching positioning error is the distance between the real track position point and the grid midpoint of the grid where the matching track located at the corresponding time (i.e., the distance between the two points ), which can reflect the effect of the grid system applied to the matching navigation, as shown in Figure 5.
In Figure 5, the dark blue line represents the matching track, the pink line represents the real track, the black point represents the grid midpoint of the matching track, the red points represent the real track points, and the distance from the real track points to the grid midpoints of the matching track at the corresponding time is the grid matching positioning error. A smaller grid matching positioning error corresponds to a higher gravity gradient matching navigation accuracy.
The two simulated tracks were matched via gravity gradient filtering, and the navigation speed was set to 10 nautical miles per hour. The measurements were made every 3 min. Taking the matched tracks as the experimental objects, the matching positioning errors of the two tracks using the different grid systems were calculated. The results are shown in Table 6.
Table 6 shows that, compared with those of corresponding resolution quadrilateral grids, the grid matching positioning errors of the hexagonal grid system at different levels and are reduced to a certain extent.
For track 1, compared with that of the 2-resolution quadrilateral grid, the average grid matching positioning error of the 11th level of the hexagonal grid is reduced by 15.26% and the standard deviation is reduced by 17.05%. Compared with that of the 1-resolution quadrilateral grid, the average grid matching positioning error of the 12th level of the hexagonal grid is reduced by 17.25% and the standard deviation is reduced by 20.36%. Compared with that of the 30-resolution quadrilateral grid, the average grid matching positioning error of the 13th level of the hexagonal grid is reduced by 19.73% and the standard deviation is reduced by 22.06%.
For track 2, compared with that of the 2-resolution quadrilateral grid, the average grid matching positioning error of the 11th level of the hexagonal grid is reduced by 15.30% and the standard deviation is reduced by 17.48%. Compared with the 1-resolution quadrilateral grid, the average grid matching positioning error of the 12th level of the hexagonal grid is reduced by 18.09% and the standard deviation is reduced by 18.17%. Compared with the 30-resolution quadrilateral grid, the average grid matching positioning error of the 13th level of the hexagonal grid is reduced by 20.54% and the standard deviation is reduced by 20.50%.
For track 1, different level hexagonal grids are used for the simulation matching. The average grid matching positioning error obtained by matching with the corresponding resolution quadrilateral grid is reduced by more than 15.26% and the standard deviation is reduced by more than 17.05%. For track 2, the average grid matching positioning error is reduced by more than 15.30% and the standard deviation is reduced by more than 17.48%. It can be seen that the hexagonal grid system has more advantages than the quadrilateral geographic grid. With increasing grid resolution, the advantages of the hexagonal grid for matching and positioning gradually increase. Compared with that of quadrilateral geographic grids at the same resolution, the standard deviation of the grid matching positioning error of the hexagonal grids at all levels is reduced by more than 15%, which can effectively suppress the matching positioning error caused by the discretization of the grid model.
3. Fast Construction Method for a Hexagonal Marine Gravity Gradient Reference Map
Because satellite gravity gradient measurements only have access to low-frequency gradient signals, the current global marine gravity gradient grid numerical model used to assist inertial navigation is generally based on satellite altimetry gravity anomalies and calculated using remove–restore technology according to Stokes’ theory [4, 10]. The gravity gradient data obtained via this method can be regarded as new observational information to some extent. However, the calculation process is complex and time-consuming and it is difficult to efficiently complete the construction requirements for a large-scale, or even global, hexagonal marine gravity gradient grid model. Therefore, it is also necessary to design an effective scheme to parallelize the construction process and improve its calculation efficiency. Figure 6 shows the calculation process of reference map’s construction.
3.1. Calculation Process for a Disturbing Gravity Gradient
According to the Stokes remove–restore method, the disturbing gravity gradient can be expressed as [26]
Here, represents the disturbing gravity gradient, represents the model disturbing gravity gradient calculated using the gravitational field model of the Earth, and represents the residual disturbing gravity gradient. Because of the symmetry of the gravity gradient, the sum of the main diagonal elements is 0 and six components (five independent components and one as the check condition) are usually selected for calculation when constructing the marine gravity gradient reference map. The calculation process and the method used for each component are similar. As an example, we introduce the calculation process for the disturbing gravity gradient.
3.1.1. Calculation of the Model Disturbing Gravity Gradient
The spherical harmonic series calculation of the ocean model disturbing gravity gradient can be expressed as [11] where is the horizontal east-west component of disturbed gravity gradient, is the average radius of the Earth, represents the geocentric latitude and longitude of the ocean calculation point, is the geocentric gravitational constant, is the mass of the earth, and represent the truncation degrees and orders, and are the coefficients of the gravitational field model of the Earth, and represents the associated Legendre function. Because the latitude of the hexagonal grid in a single latitude circle is not completely consistent, the exchange integration order method cannot be used to improve the calculation efficiency for the model disturbing gradient calculation based on the hexagonal grid.
Because the ocean gravity gradient needs to be constructed in any region of the world, we need to consider that the recurrence of the quadratic differential of the Legendre function in equation (2) will be singular near the pole. Therefore, the recurrence formula of needs to be de singular. The expression after desingularity is as follows [11].
Using equations (2) and (3), the model disturbing gravity gradient component of an arbitrary point on the globe can be calculated.
3.1.2. Calculation of the Residual Disturbing Gravity Gradient
Because the gravitational field model of the Earth is a spherical harmonic series truncated to a finite order, it can only represent the gravitational field information of the Earth to a finite order, which will inevitably produce a truncation error. The purpose of the remove–restore method is to remove the long-wave component in the disturbing gravity gradient and calculate the remaining part using the residual gravity anomaly of the surrounding points. , where represents the actual gravity anomaly and represents the model gravity anomaly calculated using the gravitational field model of the Earth. In this paper, the DTU18 model data retrieved from satellite altimetry data were selected as the actual gravity anomaly values. Because of the singularity in the calculation of the central area, the calculation process needs to be divided into a central area and an inner area.
(1) Calculation of the Residual Disturbing Gravity Gradient in the Inner Area. The influence of the residual gravity anomaly on the disturbing gravity gradient is [11] where
Using equations (4) and (5), the gravity gradient component of the residual disturbing at any point in the world can be calculated.
(2) Singular Integral Processing in the Central Region. It can be seen from equations (2) and (4) that the Stokes kernel function and its derivative are used in the calculation process; the Stokes kernel function is equal to the square of the distance between the flow point and the calculation point divided by the radial diameter. Therefore, when the flow point approaches the calculation point, its derivative will abnormally increase, resulting in a singularity, especially when calculating the radial disturbing gravity gradient [4, 12]. Therefore, the original formula cannot be used to calculate the disturbing gravity gradient in the central area. We use the Taylor series expansion method to solve the singularity problem in the central region.
First, the regional geoid can be approximately regarded as a plane. Then, the coordinate of the calculation point on the geoid is and the central area is regarded as a circular area with the calculation point at the center and a distance of 1–2 as the radius. Then, the radius is , where is the angular distance between the center of the central region and the edge. After gridding the central area, the radius can be approximated to .
If the rectangular coordinate is transformed into the polar coordinate system with the calculation point as the origin, the coordinate of the flow point in the polar coordinate system is , where , and are the coordinates of the calculation point, and and are the coordinates of the flow point. Then, it is easy to obtain , , and . Therefore, the residual disturbing gravity gradient in the central area can be obtained as follows:
The Taylor series expansion at the calculation point in equation (6) can then be obtained: where , , , , and .
Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form:
In equation (8), only the linear term (assuming that the higher-order terms are small) is considered and the residual gravity anomaly grid data can be used to obtain the coefficient parameters , , , , and .
Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form.
Here, represents the number of grid rows, , and represents the number of grid columns, .
Substituting equations (8) and (9) into equation (6) gives
The gravity gradient component of the residual disturbing in the central region at any point in the world can then be calculated using equation (10).
3.2. The CPU+GPU Hybrid Parallel Scheme
It can be seen from Section 3.1 that, because the latitude of the hexagonal grid in a single latitude circle is not completely consistent, when using the hexagon as the grid unit to construct the marine gravity gradient reference map, it is impossible to use methods such as exchanging the integral order to reduce the time-consuming calculation process of the model disturbing gradient. Moreover, the calculation of the residual disturbing gradient increases rapidly with increasing resolution, calculation range, and integration radius. Therefore, it is necessary to accelerate the construction of the hexagonal marine gravity gradient reference map using the CPU+GPU hybrid parallel algorithm to parallelize the process. According to the calculation process described in Section 3 and the calculation characteristics of the model disturbing and residual disturbing gradients, the calculation of the model disturbing gradient is allocated to the CPU and that of the residual disturbing gradient is assigned to the GPU.
3.2.1. Model Perturbation Gradient CPU Parallel Algorithm
CPUs have a strong overall regulation ability, and the memory in the CPU is easier to expand to create a large memory space for calculations. Currently, the number of CPU threads on a personal computer (PC) can reach 16 or more and the threads of a multi-CPU server can exceed 100. Even though the number of threads on a CPU is often less than that on a GPU, the single thread computing power of a CPU is much stronger than that of a GPU. The single-thread processing frequency of a CPU can generally reach more than 2 GHz. Therefore, CPUs are suitable for overall parallel tasks with large memory consumption and relatively concentrated computing. OpenMP is a parallel algorithm designed and developed for multicore/shared memory processors [13, 14]. Using OpenMP does not necessitate large modifications to the original serial program and only requires making simple improvements to avoid memory read/write conflicts during multithreaded computing. Therefore, it is suitable for parallelizing the model disturbing gradient calculation, which has a large amount of recursive computations. After determining the parallel computational mode, it is necessary to create parallel units. The parallel computing unit is the smallest unit used to allocate tasks to each thread during parallel computing, and the memory space allocated to each computing unit should be independent of the other units. Because of the differing latitudes of the hexagonal grids in the same latitude circle, the method of exchanging the integral order cannot be used to reduce the calculation amount; accordingly, the calculations can only be done one by one using a separate grid as the calculation unit. Therefore, it was decided to take a grid point as the parallel computing unit.
The division of the parallel computing units in the above manner can only be implemented when calculating the small-scale and low-order model disturbing gradient. However, when calculating large-scale and high-order models, such a division will lead to excessive memory consumption, such that most PCs will have difficulty completing the calculation. To reduce the dependence of the parallel computation on the memory, this paper proposes the “array dimensionality reduction + thread grouping” method, which can reduce the memory occupied by the parallel computation and weaken the dependence of the parallel computation on memory.
Array dimensionality reduction: when calculating the Legendre function , it is often necessary to create a two-dimensional array to store . However, the Legendre function is not a complete two-dimensional array but rather an array in the form of a triangular matrix. Directly using a two-dimensional array to store the function will waste half of the memory; therefore, a dimensionality reduction is performed, as shown in Figure 7. We can save half the memory by removing all the zeros in the two-dimensional array and then arranging the Legendre function stored in the lower triangle into a one-dimensional array from left to right for each row.
Thread grouping: in general, high-performance computers have more CPU threads and more memory, while ordinary PCs have fewer threads and less memory. Generally speaking, the number of CPU threads on a computer is often positively correlated with the size of its memory. Therefore, to make the proposed method applicable to as many types of computers as possible, this paper proposes grouping the parallel computing units according to the number of threads used in the parallel computation, as shown in Figure 8. The computing units in each group are serially computed and share memory. Taking a range of and the 12th level of the hexagonal grid in Table 3 as an example, the number of parallel computing units is equal to the number of grids: 371,059. If the computing order is 100, each parallel unit needs to create at least a double array with a size of to store the Legendre functions. Such an array requires 39.5 KB of memory. If it is directly parallelized, at least of memory is required, that is, 13.98 GB. With increasing order, the memory consumption will increase exponentially, which is not practical for ordinary PCs. According to the above “thread grouping” scheme, when the number of available threads is 100, only 3.86 MB of memory is required, which is achievable for most computers.
Using the above “array dimensionality reduction + thread grouping” method, the parallel calculation of the model disturbing gradient can be completed on a multicore processor computer.
3.2.2. Residual Disturbing Gradient GPU Parallel Algorithm
CUDA is a heterogeneous parallel platform developed by NVIDIA for its GPU products [15]. With the further development of hardware technology, the number of GPU cores is increasing [24, 25]. Taking rtx3080 as an example, the number of CUDA cores has reached 5,120, with a price of less than 10,000 yuan, which greatly reduces the threshold for parallel computing. Ordinary computers equipped with this type of graphics card can use the CUDA platform to complete parallel computing tasks. Even though GPUs have a large number of threads, their disadvantages are more prominent than those of CPUs. First, the single-thread GPU computing ability is weak and GPUs have high requirements for the amount of computation of a parallel computing unit. If the amount of computation of the parallel computing unit is too large, the GPU will crash. Second, the GPU video memory is smaller than the computer memory and cannot be expanded; therefore, it is suitable for operations with less initial data but more computing units. The calculation process of the residual disturbing gradient conforms to the above characteristics; therefore, a GPU was selected to accelerate its calculation process in parallel.
The CUDA programming language organizes the GPU threads in the form of a thread grid, as shown in Figure 9. There are 15 () thread blocks in a single thread grid, each of which contains 9 () threads. To accurately index these threads, four built-in variables are introduced: threadIdx, representing the index in the thread block; blockIdx, representing the index in the thread grid of the thread block; and blockDim and gridDim, representing the dimensions (i.e., the length and width) of the thread grid and the thread block, respectively [16–18]. Using these variables, the and indexes of a thread in the process grid can be calculated. The calculation formulas are shown in
When calculating various gravity gradient maps, the number of grid points is large. Taking a range of and the 12th level of the hexagonal grid as an example, there are 371,059 grid points and the edges of the hexagonal grid are irregular. Directly using the two-dimensional grid for indexing not only increases the time consumption of the thread scheduling process [16, 19] but also is not conducive to thread indexing at the edge of the grid. Therefore, it is necessary to reduce the dimension of the thread grid; however, after dimension reduction, the original adjacent grids will be separated, resulting in a reduction of the overall spatial clustering. Therefore, spatial clustering of the one-dimensional grid needs to be ensured as much as possible during the dimension reduction process. Accordingly, this paper uses the one-dimensional thread indexing method to divide the tasks of the GPU threads, taking the calculation of the residual disturbing gravity gradient of any grid point as a calculation unit and sorting all of the calculation units. To ensure the spatial clustering of the grid data in the one-dimensional indexing process, the regional hexagonal grid is sorted accordingly, as shown in Figure 10, and the grid in the left panel is sorted into rows (or columns). Each row (or column) is connected end to end into the form shown in the right panel, the grid of which is recoded. Taking the red hexagonal grid in the 3rd row of the 15th column in the figure as an example, its grid code after dimension reduction sorting is 150.
Using the parallel computing scheme designed in Sections 4.1 and 4.2, the rapid construction of a hexagonal marine gravity gradient reference map can be completed using CPU+GPU hybrid programming.
4. Experimental Verification and Analysis
To verify the effectiveness of the proposed method, the algorithm was experimentally tested on the following hardware and software experimental platforms. The hardware environment was a high-performance desktop computer with an Intel®Xeon(R) Gold 6130 @ 2.10 GHz processor having 22 cores and 44 threads. The GPU was a Tesla V100 16G with 64.00 GB of memory. The software environment was a 64-bit Windows 10 operating system, with a VS2013+CUDA9 2 platform.
The selection range of the test area (the range of the marine grid) was from 19°N to 29°N and from 41°W to 51°W. The deepest sea depth in this area is 6,872 m, the shallowest sea depth is 0 m, and the regional average elevation is −3,657.22 m. The maximum value of the gravity anomaly in this area is 89.8 mGal, the minimum value is −83.3 mGal, and the average gravity anomaly is 4.4 mGal. The specific relationship between the regional gravity anomaly and the submarine topography is shown in Figure 11.
(a) DTU18 gravity anomaly
(b) Submarine topography
The GPU is sensitive to the selection of the data types in the calculation, and the calculation efficiency of different data types varies greatly [16, 18]. The efficiency of using single-precision floating-point data for calculations is four times that of using double-precision data. Therefore, when using GPUs for parallel calculations, single-precision data should be used as much as possible to ensure the accuracy and improve the calculation efficiency. However, in the process of calculating the residual disturbing gravity gradient tensor, it is necessary to use sine and cosine functions and the square operation. The influence of the truncation error caused by using only single-precision data on the calculation results cannot be ignored. Therefore, this section uses both single-precision floating-point data and double-precision type data and uses the GPU to calculate the total tensor of the residual disturbing gravity gradient in the experimental area. We take the serial calculation results of the CPU using the double-precision data as the “true value” and perform a comparative analysis with the GPU calculation results, so as to select the data type of the parallel calculation and verify the reliability of the proposed parallel calculation method. By comparing the GPU calculation results with the CPU calculation results, the accuracy of the parallel calculation results of the different data types on the GPU can be obtained, as shown in Table 7.
Table 7 shows that the significant-digit limitation causes the truncation error generated by the calculation results using single-precision data to have a significant impact on the results. The maximum error is , the minimum error is , the average error is , and the root mean square (RMS) error is . The absolute value of the error calculated with double-precision data is less than or equal to , and the average difference and RMS value are less than . The experimental results show that, even though the single-precision data can effectively improve the calculation efficiency, this precision has difficulty meeting the calculation requirements of the disturbing gravity gradient tensor. Therefore, we chose to use double-precision data to construct the gravity gradient reference map.
To further verify the reliability of the mathematical model of the parallel scheme and the correctness of the programming software, according to the constraint condition , the above parallel calculation scheme was used to calculate the disturbing gravity gradient tensor and evaluate the internal consistency of the calculation results. The experimental results are shown in Table 8.
It can be seen from Table 8 that the maximum value of is , the minimum value is , the standard deviation is , and the RMS error is . The calculation results fluctuate randomly around 0 and meet the internal consistency conditions of the disturbing gravity gradient tensor. This, to a certain extent, indicates the reliability of the mathematical model of the proposed parallel scheme and the correctness of programming software.
4.1. Efficiency Analysis of the Parallel Algorithm
In the previous section, the data type selection, calculation reliability, and calculation accuracy of the parallel computing were analyzed. This section analyzes the efficiency of the proposed parallel computing scheme on the basis of the results found in the previous section.
4.1.1. Efficiency Analysis of Different Parallel Schemes
To fully verify the effectiveness of the proposed parallel scheme, four calculation schemes are examined in this section: (1) a serial computing scheme, (2) a pure CPU parallel computing scheme, (3) a pure GPU parallel computing scheme, and (4) the proposed CPU+GPU hybrid parallel scheme. These four calculation schemes were used to calculate the disturbing gravity gradient tensor in the test area shown in Table 9. The 12th level of the ISEA4H hexagonal grid was used for the calculation, and the EIGEN-6C4 gravity field model was selected to calculate the disturbing gravity gradient tensor of the model, with the truncation order calculated to 2,160 [20–23]. The actual gravity anomaly of the DTU18 data retrieved via satellite altimetry data was selected, and the actual gravity anomaly of the hexagonal grid was interpolated. The central area incorporated the grid area 1 from the midpoint of the central grid, and the integral radius was in the range of 1°. The calculation results are shown in Figure 12, and the calculation time statistics are shown in Table 9.
(a)
(b)
(c)
(d)
(e)
Table 9 indicates that the total tensor of the disturbing gravity gradient in the experimental area is calculated by the serial calculation in 120.35 h or more than 5 days. The pure CPU parallel computing scheme takes 3.18 h, for an acceleration ratio of 40.33 times. The pure GPU parallel computing scheme takes 6.74 h and has an acceleration ratio of 19.03 times. The CPU+GPU hybrid parallel computing scheme takes 2.99 h and has an acceleration ratio of 42.90 times. The acceleration effect of the hybrid parallel computing scheme is better than that of the CPU or GPU single system parallel computing schemes. When only the CPU is used for the parallel calculation, the acceleration ratio of calculating the disturbing gravity gradient of the model is 40.38 times and the acceleration ratio of calculating the residual disturbing gravity gradient is 39.24 times. Using only the GPU for the parallel calculation results in the acceleration ratio of calculating the disturbing gravity gradient of the model being 17.86 times and the acceleration ratio of calculating the residual disturbing gravity gradient being 1126.38 times. It can be seen that the proposed CPU+GPU hybrid parallel computing scheme can fully exploit the available computing power and effectively improve the computing efficiency according to the CPU and GPU characteristics.
4.1.2. Calculation Efficiency Analysis for Different Resolutions
To further verify the universal effectiveness of the proposed CPU+GPU parallel algorithm, here we use ISEA4H hexagonal grid data at different levels to conduct numerical experiments, calculate the marine area within the test area, and calculate the model disturbing gravity gradient to order 2,160, with an integral radius of 1°. The hexagonal grid levels used are the 11th, 12th, and 13th levels. The calculation time and acceleration are shown in Table 10.
It can be seen from Table 10 that, with increasing hexagonal grid level, the calculation time increases; therefore, a large-scale high-resolution hexagonal gravity gradient reference map is increasingly difficult to obtain in a short time. The proposed CPU+GPU parallel algorithm was used to calculate the marine gravity gradient reference map within the test area. The acceleration ratio of the hexagonal gravity gradient reference map at the different levels can reach more than 40 times and up to 50.92 times. It takes 0.70 h to construct the hexagonal marine gravity gradient reference map at the 11th level with a range of , while it takes 28.45 h for a traditional serial calculation. In parallel, it takes 2.98 h to construct of the 12th-level hexagonal marine gravity gradient reference map, while the serial method takes 120.35 h or more than 5 days. In parallel, it takes 11.84 h to construct the 13th-level hexagonal marine gravity gradient reference map, while the serial method takes 602.95 h or more than 25 days. This demonstrates that the proposed parallel scheme is generally effective for the construction of hexagonal disturbing gravity gradient reference maps at different levels and can effectively improve the construction speed of gravity gradient reference maps at various resolutions.
Table 10 shows that, for increasing hexagonal grid level (resolution), the calculation ratio of the proposed hybrid parallel scheme increases. The main reason for this is that the calculation time of the disturbing gravity gradient part of the different level models is different from that of the residual disturbing gravity gradient part. Taking level 12 as an example, even though the calculation ratio using the GPU to calculate the residual disturbing gravity gradient can reach more than 1,000 times, the overall acceleration ratio is less than 43 times. The main reason is that, when using a hexagonal grid to calculate the disturbing gravity gradient, because the latitude of the hexagonal grid in the same latitude circle is not completely consistent, it is impossible to use methods such as exchanging the integration order to accelerate the calculation process of the model disturbing gradient. In the process of calculating the disturbing gravity gradient using the 12th-level hexagonal grid, the calculation time of the gravity gradient part of the model accounts for 93.8% of the total calculation time, while that of the residual gravity gradient part accounts for 6.2% of the total. Therefore, the acceleration effect of the calculation of the model part is dominant. With increasing level (resolution), the calculation amount of the disturbing gradient part of the model increases linearly, while that of the residual disturbing gradient part increases exponentially; accordingly, the proportion of the calculation acceleration effect of the residual part gradually increases. Therefore, the acceleration effect of the CPU+GPU hybrid parallel scheme gradually increases with increasing level (resolution), making it more suitable than other schemes to building a high-resolution marine gravity gradient reference map.
5. Construction of the Global Marine Hexagonal Gravity Gradient Reference Map
Based on the proposed CPU+GPU hybrid parallel scheme, we used the ISEA4H hexagonal grid system to rapidly construct a global, full-tensor, 12th-level hexagonal marine gravity gradient reference map. The EIGEN-6C4 gravity field model was used for the calculation, and the truncation was calculated to order 2160. Using the satellite altimetry DTU18 model gravity anomaly, the hexagonal grid data were formed first and then the residual disturbing gravity gradient was calculated. An integral radius of 1° was selected for the calculation in the inner area, and the central area selected the central grid point and the six surrounding grid points, totaling seven grid areas. The calculation results are shown in Table 11 and Figure 13, and the internal consistency of the gravity gradient tensor is shown in Table 12.
(a)
(b)
(c)
(d)
(e)
The results in Tables 11 and 12 show that the maximum absolute value of the global full-tensor hexagonal marine gravity gradient reference map model is less than and that the standard deviation is less than , which meets the internal conditions of the disturbing gravity gradient tensor. The global marine gravity gradient value is small. The absolute value of the gravity gradient of more than 80% of the grid points is less than , the average value is less than , the median is near , and the gravity gradient is related to the seabed topography. The variation in the gravity gradient value is obvious in areas with severe topographic fluctuations.
6. Discussion and Conclusions
This paper first introduced the hexagonal grid system as a grid unit for constructing a marine gravity gradient reference map, compared the advantages and disadvantages of the ISEA4H and H3 grid systems, and selected the ISEA4H hexagonal grid system according to the characteristics of gravity gradient matching navigation. Then, via calculations, statistics, and analyses of the grid model error and grid matching positioning error, the feasibility and effectiveness of the hexagonal grid for gravity gradient matching navigation were verified. Consequently, based on Stokes’ remove–restore theory and according to the characteristics of the hexagonal grid disturbing gravity gradient calculation, we proposed an efficient CPU+GPU hybrid parallel scheme, which can effectively improve the construction efficiency of a hexagonal marine gravity gradient reference map. The experimental calculations and analyses show the following. (1)Compared with the H3 hexagonal grid system, the ISEA4H hexagonal grid system has a better characteristic of equal area; that is, the number of grids between adjacent levels changes less. In addition, the grid levels of the ISEA4H system better correlate with quadrilateral geographic grids, which is conducive to the construction of reference maps based on existing data. Therefore, the ISEA4H grid system was selected to construct the hexagonal marine gravity gradient reference map(2)On the premise of only considering errors caused by the geometric structure of the reference grid in the process of discretizing a track, the grid model error and matching positioning error of a hexagonal grid and quadrilateral geographic grid with the same resolution were analyzed for the process of matching navigation. The experimental results show that the standard deviation of the grid model error of a hexagonal grid is more than 14% less than that of a quadrilateral grid with the same resolution. The standard deviation of the grid matching positioning error is reduced by more than 15% for a hexagonal grid, which effectively proves the feasibility and effectiveness of hexagonal grids in gravity gradient matching navigation(3)Based on the hexagonal grid, we constructed a reference map of the marine gravity gradient using the Stokes remove–restore technique. To address the low efficiency of calculating the disturbing gravity gradient on a hexagonal grid, a CPU+GPU hybrid parallel computing scheme was proposed. This scheme can make full use of the computing power of the available equipment and effectively improve the computing efficiency. Compared with a serial algorithm, the calculation acceleration ratio of the parallel computing scheme can reach more than 40 times, which is beneficial for the construction of high-resolution gravity gradient reference maps. A global full-tensor hexagonal marine gravity gradient reference map was constructed under the following accuracy conditions: an absolute error of the parallel computing of less than or equal to , a maximum absolute value of less than , and a standard deviation of less than
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
Huang Yan and Li Shanshan were responsible for conceptualization; Li Xinxing and Huang Yan were responsible for methodology and software; Fan Diao and Tan Xuli were responsible for validation; Wang Aoming was responsible for data curation; Huang Yan was responsible for writing—original draft preparation; Lv Minghao was responsible for writing—review and editing; Li Shanshan was responsible for funding acquisition. All authors have read and agreed to the published version of the manuscript.
Acknowledgments
This research was funded by the National Natural Science Foundation of China, grant numbers 42174007 and 42174013.