# Research and Construction of a Global Hexagonal Marine Gravity Gradient Reference Map for Navigation

Authors: Huang Yan; Li Shanshan; Zhou Chuyuan; Li Xinxing; Lv Minghao; Fan Diao; Tan Xuli; Wang Aoming
Journal: Geofluids
 (2023)

---

## 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.


---

## Body

## 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 gradientTM and the residual disturbing gravity gradient δT. Finally, according to the characteristics of the hexagonal TM and δT 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 Figure1.Figure 1
Positive icosahedral positioning diagram.In Figure1, 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 Figure1 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.Figure 2
Hexagonal grid generation process.Table 1
Parameter statistics of the icosahedral Snyder equal area aperture 4 hexagon (ISEA4H) grid (with a spherical radius ofR=6378.136km).
LevelLNumber of gridsGrid area (km2)Grid spacing (km)0124.260×1077.530×1031421.217×1073.682×10321623.156×1061.968×10336427.963×1051.002×10342,5621.995×1055.035×102510,2424.991×1042.520×102640,9621.248×1041.260×1027163,8423.120×1036.303×108655,3627.800×1023.152×1092,621,4421.950×1021.576×101010,485,7624.875×107.8791141,943,04212.1613.76712167,772,1623.0401.88413671,088,6420.7600.942142,684,354,5620.1900.471n10×4n+2322n+1L24R⋅arcsin12n32π⋅L2RNote:L represents the hexagonal grid level, and R represents the average radius of the Earth.


### 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.Table 2
Parameter statistics of the H3 grid (with a spherical radius ofR=6378.136km).
LevelLNumber of gridArea of grid (km2)Average side length (km)01224,250,546.84770001,107.7125910001842607,220.9782429418.67600550025,88286,745.8540347158.244655800341,16212,392.264862159.8108579404288,1221,770.323551722.60637940052,016,842252.90336458.544408276614,117,88236.12905213.229482772798,825,1625.16129321.2206297598691,776,1220.73732760.46135468494,842,432,8420.10533250.1743756681033,897,029,8820.01504750.06590780711237,279,209,1620.00214960.024910561121,660,954,464,1220.00030710.0094155261311,626,681,248,8420.00004390.0035598931481,386,768,741,8820.00000630.00134857515569,707,381,193,1620.00000090.000509713Note: the information in this table is taken from the Internet (https://h3geo.org/).


### 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 Figure3(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 levelsFigure 3
Schematic diagrams of four- and seven-aperture division.
(a)
Four apertures(b)
Seven aperturesBased 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 of10°×10° (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.Table 3
Error settings of the inertial components.
Inertial componentsError termError parameterGyroConstant drift0.001∘/hAngle random walk0.001∘/hAccelerometerZero bias50μgSpeed random walk5μg/HzWe 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 Figure4, 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.Figure 4
Schematic diagram of the grid model error.In Figure4, 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 Table4.Table 4
Correspondence between the ISEA4H hexagonal grids and the quadrilateral grids.
GridResolution (level)Distance between adjacent grid midpoint (km)Average area of grid (km2)Number of grids in the test areaHexagon11th level3.76712.16192,76712th level1.8843.040371,05913th level0.9420.7601,484,237Quadrilateral2′3.71113.77190,0001′1.8553.441360,00030″0.9270.8591,440,000Table4 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
Grid model error statistics, where Std indicates the standard deviation.
TrackGridResolutionTrack grid numberMaximum (km)Minimum (km)Average (km)Std (km)1Hexagon11th level3052.13750.00460.99550.5725Quadrilateral2′3612.53760.01311.24420.7167Hexagon12th level6061.09440.00150.49640.2845Quadrilateral1′7211.26680.00300.62050.3584Hexagon13th level12020.55300.00100.24560.1420Quadrilateral30′′14350.63140.00120.30860.17872Hexagon11th level3142.15320.01741.05550.6103Quadrilateral2′3482.53760.00521.22740.7117Hexagon12th level6251.10010.00150.52330.3042Quadrilateral1′6951.26680.00360.61230.3573Hexagon13th level12550.55300.00110.26380.1524Quadrilateral30′′13900.63140.00210.30680.1778Table5 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Δx2+Δy2), which can reflect the effect of the grid system applied to the matching navigation, as shown in Figure 5.Figure 5
Schematic diagram of the grid matching positioning error.In Figure5, 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 Table6.Table 6
Statistics of the grid matching positioning error.
TrackMatching pointsGridResolutionMaximum (km)Minimum (km)Average (km)Std (km)1519Hexagon11th level4.55870.29732.01171.2539Quadrilateral2′5.45600.25692.37411.5116Hexagon12th level3.98880.10451.56010.8121Quadrilateral1′4.71490.12941.88551.0197Hexagon13th level2.44340.03450.91360.4129Quadrilateral30′′3.61270.03111.13820.52972505Hexagon11th level4.74540.31281.91201.3632Quadrilateral2′5.61100.28122.25751.6521Hexagon12th level3.81010.14631.55810.9199Quadrilateral1′4.35250.21301.90241.1242Hexagon13th level2.50120.10380.90040.4008Quadrilateral30′′3.56490.12491.13320.5041Table6 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.


## 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 Figure1.Figure 1
Positive icosahedral positioning diagram.In Figure1, 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 Figure1 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.Figure 2
Hexagonal grid generation process.Table 1
Parameter statistics of the icosahedral Snyder equal area aperture 4 hexagon (ISEA4H) grid (with a spherical radius ofR=6378.136km).
LevelLNumber of gridsGrid area (km2)Grid spacing (km)0124.260×1077.530×1031421.217×1073.682×10321623.156×1061.968×10336427.963×1051.002×10342,5621.995×1055.035×102510,2424.991×1042.520×102640,9621.248×1041.260×1027163,8423.120×1036.303×108655,3627.800×1023.152×1092,621,4421.950×1021.576×101010,485,7624.875×107.8791141,943,04212.1613.76712167,772,1623.0401.88413671,088,6420.7600.942142,684,354,5620.1900.471n10×4n+2322n+1L24R⋅arcsin12n32π⋅L2RNote:L represents the hexagonal grid level, and R represents the average radius of the Earth.


## 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.Table 2
Parameter statistics of the H3 grid (with a spherical radius ofR=6378.136km).
LevelLNumber of gridArea of grid (km2)Average side length (km)01224,250,546.84770001,107.7125910001842607,220.9782429418.67600550025,88286,745.8540347158.244655800341,16212,392.264862159.8108579404288,1221,770.323551722.60637940052,016,842252.90336458.544408276614,117,88236.12905213.229482772798,825,1625.16129321.2206297598691,776,1220.73732760.46135468494,842,432,8420.10533250.1743756681033,897,029,8820.01504750.06590780711237,279,209,1620.00214960.024910561121,660,954,464,1220.00030710.0094155261311,626,681,248,8420.00004390.0035598931481,386,768,741,8820.00000630.00134857515569,707,381,193,1620.00000090.000509713Note: the information in this table is taken from the Internet (https://h3geo.org/).


## 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 Figure3(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 levelsFigure 3
Schematic diagrams of four- and seven-aperture division.
(a)
Four apertures(b)
Seven aperturesBased 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 of10°×10° (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.Table 3
Error settings of the inertial components.
Inertial componentsError termError parameterGyroConstant drift0.001∘/hAngle random walk0.001∘/hAccelerometerZero bias50μgSpeed random walk5μg/HzWe 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 Figure4, 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.Figure 4
Schematic diagram of the grid model error.In Figure4, 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 Table4.Table 4
Correspondence between the ISEA4H hexagonal grids and the quadrilateral grids.
GridResolution (level)Distance between adjacent grid midpoint (km)Average area of grid (km2)Number of grids in the test areaHexagon11th level3.76712.16192,76712th level1.8843.040371,05913th level0.9420.7601,484,237Quadrilateral2′3.71113.77190,0001′1.8553.441360,00030″0.9270.8591,440,000Table4 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
Grid model error statistics, where Std indicates the standard deviation.
TrackGridResolutionTrack grid numberMaximum (km)Minimum (km)Average (km)Std (km)1Hexagon11th level3052.13750.00460.99550.5725Quadrilateral2′3612.53760.01311.24420.7167Hexagon12th level6061.09440.00150.49640.2845Quadrilateral1′7211.26680.00300.62050.3584Hexagon13th level12020.55300.00100.24560.1420Quadrilateral30′′14350.63140.00120.30860.17872Hexagon11th level3142.15320.01741.05550.6103Quadrilateral2′3482.53760.00521.22740.7117Hexagon12th level6251.10010.00150.52330.3042Quadrilateral1′6951.26680.00360.61230.3573Hexagon13th level12550.55300.00110.26380.1524Quadrilateral30′′13900.63140.00210.30680.1778Table5 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Δx2+Δy2), which can reflect the effect of the grid system applied to the matching navigation, as shown in Figure 5.Figure 5
Schematic diagram of the grid matching positioning error.In Figure5, 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 Table6.Table 6
Statistics of the grid matching positioning error.
TrackMatching pointsGridResolutionMaximum (km)Minimum (km)Average (km)Std (km)1519Hexagon11th level4.55870.29732.01171.2539Quadrilateral2′5.45600.25692.37411.5116Hexagon12th level3.98880.10451.56010.8121Quadrilateral1′4.71490.12941.88551.0197Hexagon13th level2.44340.03450.91360.4129Quadrilateral30′′3.61270.03111.13820.52972505Hexagon11th level4.74540.31281.91201.3632Quadrilateral2′5.61100.28122.25751.6521Hexagon12th level3.81010.14631.55810.9199Quadrilateral1′4.35250.21301.90241.1242Hexagon13th level2.50120.10380.90040.4008Quadrilateral30′′3.56490.12491.13320.5041Table6 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.


## 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 Figure4, 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.Figure 4
Schematic diagram of the grid model error.In Figure4, 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 Table4.Table 4
Correspondence between the ISEA4H hexagonal grids and the quadrilateral grids.
GridResolution (level)Distance between adjacent grid midpoint (km)Average area of grid (km2)Number of grids in the test areaHexagon11th level3.76712.16192,76712th level1.8843.040371,05913th level0.9420.7601,484,237Quadrilateral2′3.71113.77190,0001′1.8553.441360,00030″0.9270.8591,440,000Table4 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
Grid model error statistics, where Std indicates the standard deviation.
TrackGridResolutionTrack grid numberMaximum (km)Minimum (km)Average (km)Std (km)1Hexagon11th level3052.13750.00460.99550.5725Quadrilateral2′3612.53760.01311.24420.7167Hexagon12th level6061.09440.00150.49640.2845Quadrilateral1′7211.26680.00300.62050.3584Hexagon13th level12020.55300.00100.24560.1420Quadrilateral30′′14350.63140.00120.30860.17872Hexagon11th level3142.15320.01741.05550.6103Quadrilateral2′3482.53760.00521.22740.7117Hexagon12th level6251.10010.00150.52330.3042Quadrilateral1′6951.26680.00360.61230.3573Hexagon13th level12550.55300.00110.26380.1524Quadrilateral30′′13900.63140.00210.30680.1778Table5 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Δx2+Δy2), which can reflect the effect of the grid system applied to the matching navigation, as shown in Figure 5.Figure 5
Schematic diagram of the grid matching positioning error.In Figure5, 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 Table6.Table 6
Statistics of the grid matching positioning error.
TrackMatching pointsGridResolutionMaximum (km)Minimum (km)Average (km)Std (km)1519Hexagon11th level4.55870.29732.01171.2539Quadrilateral2′5.45600.25692.37411.5116Hexagon12th level3.98880.10451.56010.8121Quadrilateral1′4.71490.12941.88551.0197Hexagon13th level2.44340.03450.91360.4129Quadrilateral30′′3.61270.03111.13820.52972505Hexagon11th level4.74540.31281.91201.3632Quadrilateral2′5.61100.28122.25751.6521Hexagon12th level3.81010.14631.55810.9199Quadrilateral1′4.35250.21301.90241.1242Hexagon13th level2.50120.10380.90040.4008Quadrilateral30′′3.56490.12491.13320.5041Table6 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.Figure 6
Reference map’s construction flowchart.


### 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]
(1)Tαβ=TαβM+δTαβ,α,β=x,y,z.Here,Tαβ represents the disturbing gravity gradient, TαβM represents the model disturbing gravity gradient calculated using the gravitational field model of the Earth, and δTαβ 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 Txx 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]
(2)TxxMR,φ,λ=fMR3∑n=2N∑m=0nC¯∗nmcosmλ+S¯nmsinmλPnmxxsinφ,where TxxM is the horizontal east-west component of disturbed gravity gradient, R is the average radius of the Earth, φ,λ represents the geocentric latitude and longitude of the ocean calculation point, f is the geocentric gravitational constant, M is the mass of the earth, n and m represent the truncation degrees and orders, C¯n,m∗ and S¯n,m are the coefficients of the gravitational field model of the Earth, and P¯nmsinφ 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 Pnmxxsinφ needs to be de singular. The expression after desingularity is as follows [11].
(3)Pnmxxsinφ=12nn−1n+1n+22P¯n,2sinφ−12n+1n+2P¯n,0sinφm=0,14n−2n−1n+2n+3P¯n,3sinφ−143n+1n+2P¯n,1sinφm=1,14n−3n−2n+3n+4P¯n,4sinφ−12n2+3n−2P¯n,2sinφ+142nn−1n+1n+2P¯n,0sinφm=2,14n−m−1n−mn+m+1n+m+2P¯n,m+2sinφ+14n+m−1n+mn−m+1n−m+2P¯n,m−2sinφ−12n2+3n−m2+2P¯n,msinφm>2.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δΔg of the surrounding points. δΔg=Δg−ΔgM, where Δg represents the actual gravity anomaly and ΔgM 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]
(4)δTxx=R4π∬σ1r2cosφ2∂2Sr,ψ∂ψ2∂ψ∂λ2+∂Sr,ψ∂ψ∂2ψ∂λ2−tanφr2∂Sr,ψ∂ψ∂ψ∂φ+1r∂Sr,ψ∂rδΔgdσ,where
(5)∂Sr,ψ∂r=−t2R1−t2D3+4D+1−6D−tcosψ13+6ln1−tcosψ+D2,∂Sr,ψ∂ψ=−t2sinψ2D3+6D−8−31−tcosψ−DDsinψ−3ln1−tcosψ+D2,∂2Sr,ψ∂r2=t2R21−tcosψ31−t2D5−4D3−1+t2D3−10D−18D+2−3tcosψ15+6ln1−tcosψ+D2,∂2Sr,ψ∂ψ2=−t2cosψ2D3+6D−8−31−tcosψ−DDsin2ψ−3ln1−tcosψ+D2+t3sin2ψ6D5+6D3+3D−1D2sin2ψ+3D+1D1−tcosψ+D−31−tcosψ−DDsin2ψ2cosψtsin2ψ+1D2,∂2Sr,ψ∂r∂ψ=t3sinψR31−t2D5+4D3+6D−13−6ln1−tcosψ+D2+6cosψD+1D1−tcosψ+Dt=Rr,D=1−2tcosψ+t21/2ψ=cos−1sinφsinφ′+cosφcosφ′cosλ′−λ.Using equations (4) and (5), the gravity gradient component of the residual disturbing δTxx 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 isx,y,0 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 s0=Rψ0, where ψ0 is the angular distance between the center of the central region and the edge. After gridding the central area, the radius s0 can be approximated to Δx×Δy.If the rectangular coordinatex,y 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 s,α, where s=x′−x2+y′−y2, x and y are the coordinates of the calculation point, and x′ and y′ are the coordinates of the flow point. Then, it is easy to obtain x=scosα, y=ssinα, and dx′dy′=sdsdα. Therefore, the residual disturbing gravity gradient in the central area can be obtained as follows:
(6)δTxx=12π∫02π∫0s0δΔgs23cos2α−1dsdα.The Taylor series expansionδΔg at the calculation point in equation (6) can then be obtained:
(7)δΔg=δΔgP+xg1+yg2+12x2g11+2xyg12+y2g22+⋯,where g1=∂Δg/∂x, g2=∂Δg/∂y, g11=∂2Δg/∂x2, g2=∂2Δg/∂y2, and g12=∂2Δg/∂x∂y.Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form:
(8)δΔg=δΔgP+sg1cosα+g2sinα+s22g11cos2α+2g12sinαcosα+g22sin2α+⋯.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 g1, g2, g11, g12, and g22.Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form.
(9)g1n,m=δΔgn+1,m−δΔgn−1,m2Δx,g2n,m=δΔgn,m+1−δΔgn,m2Δy,g11n,m=g1n+1,m−g1n−1,m2Δx,g12n,m=g1n,m+1−g1n,m−12Δy,g22n,m=g2n,m+1−g2n,m−12Δy.Here,n represents the number of grid rows, n=1,⋯,N, and m represents the number of grid columns, m=1,⋯,M.Substituting equations (8) and (9) into equation (6) gives
(10)δTxx=5g11−g2216s0.The gravity gradient component of the residual disturbingδTxx 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 Section3.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 functionP¯nmcosθ, it is often necessary to create a two-dimensional array PnM to store P¯nmcosθ. However, the Legendre function P¯nmcosθ 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.Figure 7
Schematic diagram of the array dimensionality reduction.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 Figure8. The computing units in each group are serially computed and share memory. Taking a range of 10°×10° 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 1+100×100/2 to store the Legendre functions. Such an array requires 39.5 KB of memory. If it is directly parallelized, at least 371,059×39.5KB 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 I is 100, only 3.86 MB of memory is required, which is achievable for most computers.Figure 8
Schematic diagram of the grouping calculation.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 Figure9. There are 15 (3×5) thread blocks in a single thread grid, each of which contains 9 (3×3) 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 x and y indexes of a thread in the process grid can be calculated. The calculation formulas are shown in
(11)Idx=blockIdx.x∗blockDim.x+threadIdx.x,Idy=blockIdx.y∗blockDim.y+threadidx.y.Figure 9
CUDA thread index diagram.When calculating various gravity gradient maps, the number of grid points is large. Taking a range of10°×10° 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.Figure 10
Schematic diagram of the one-dimensional coding of a hexagonal grid.Using the parallel computing scheme designed in Sections4.1 and 4.2, the rapid construction of a hexagonal marine gravity gradient reference map can be completed using CPU+GPU hybrid programming.


## 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]
(1)Tαβ=TαβM+δTαβ,α,β=x,y,z.Here,Tαβ represents the disturbing gravity gradient, TαβM represents the model disturbing gravity gradient calculated using the gravitational field model of the Earth, and δTαβ 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 Txx 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]
(2)TxxMR,φ,λ=fMR3∑n=2N∑m=0nC¯∗nmcosmλ+S¯nmsinmλPnmxxsinφ,where TxxM is the horizontal east-west component of disturbed gravity gradient, R is the average radius of the Earth, φ,λ represents the geocentric latitude and longitude of the ocean calculation point, f is the geocentric gravitational constant, M is the mass of the earth, n and m represent the truncation degrees and orders, C¯n,m∗ and S¯n,m are the coefficients of the gravitational field model of the Earth, and P¯nmsinφ 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 Pnmxxsinφ needs to be de singular. The expression after desingularity is as follows [11].
(3)Pnmxxsinφ=12nn−1n+1n+22P¯n,2sinφ−12n+1n+2P¯n,0sinφm=0,14n−2n−1n+2n+3P¯n,3sinφ−143n+1n+2P¯n,1sinφm=1,14n−3n−2n+3n+4P¯n,4sinφ−12n2+3n−2P¯n,2sinφ+142nn−1n+1n+2P¯n,0sinφm=2,14n−m−1n−mn+m+1n+m+2P¯n,m+2sinφ+14n+m−1n+mn−m+1n−m+2P¯n,m−2sinφ−12n2+3n−m2+2P¯n,msinφm>2.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δΔg of the surrounding points. δΔg=Δg−ΔgM, where Δg represents the actual gravity anomaly and ΔgM 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]
(4)δTxx=R4π∬σ1r2cosφ2∂2Sr,ψ∂ψ2∂ψ∂λ2+∂Sr,ψ∂ψ∂2ψ∂λ2−tanφr2∂Sr,ψ∂ψ∂ψ∂φ+1r∂Sr,ψ∂rδΔgdσ,where
(5)∂Sr,ψ∂r=−t2R1−t2D3+4D+1−6D−tcosψ13+6ln1−tcosψ+D2,∂Sr,ψ∂ψ=−t2sinψ2D3+6D−8−31−tcosψ−DDsinψ−3ln1−tcosψ+D2,∂2Sr,ψ∂r2=t2R21−tcosψ31−t2D5−4D3−1+t2D3−10D−18D+2−3tcosψ15+6ln1−tcosψ+D2,∂2Sr,ψ∂ψ2=−t2cosψ2D3+6D−8−31−tcosψ−DDsin2ψ−3ln1−tcosψ+D2+t3sin2ψ6D5+6D3+3D−1D2sin2ψ+3D+1D1−tcosψ+D−31−tcosψ−DDsin2ψ2cosψtsin2ψ+1D2,∂2Sr,ψ∂r∂ψ=t3sinψR31−t2D5+4D3+6D−13−6ln1−tcosψ+D2+6cosψD+1D1−tcosψ+Dt=Rr,D=1−2tcosψ+t21/2ψ=cos−1sinφsinφ′+cosφcosφ′cosλ′−λ.Using equations (4) and (5), the gravity gradient component of the residual disturbing δTxx 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 isx,y,0 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 s0=Rψ0, where ψ0 is the angular distance between the center of the central region and the edge. After gridding the central area, the radius s0 can be approximated to Δx×Δy.If the rectangular coordinatex,y 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 s,α, where s=x′−x2+y′−y2, x and y are the coordinates of the calculation point, and x′ and y′ are the coordinates of the flow point. Then, it is easy to obtain x=scosα, y=ssinα, and dx′dy′=sdsdα. Therefore, the residual disturbing gravity gradient in the central area can be obtained as follows:
(6)δTxx=12π∫02π∫0s0δΔgs23cos2α−1dsdα.The Taylor series expansionδΔg at the calculation point in equation (6) can then be obtained:
(7)δΔg=δΔgP+xg1+yg2+12x2g11+2xyg12+y2g22+⋯,where g1=∂Δg/∂x, g2=∂Δg/∂y, g11=∂2Δg/∂x2, g2=∂2Δg/∂y2, and g12=∂2Δg/∂x∂y.Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form:
(8)δΔg=δΔgP+sg1cosα+g2sinα+s22g11cos2α+2g12sinαcosα+g22sin2α+⋯.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 g1, g2, g11, g12, and g22.Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form.
(9)g1n,m=δΔgn+1,m−δΔgn−1,m2Δx,g2n,m=δΔgn,m+1−δΔgn,m2Δy,g11n,m=g1n+1,m−g1n−1,m2Δx,g12n,m=g1n,m+1−g1n,m−12Δy,g22n,m=g2n,m+1−g2n,m−12Δy.Here,n represents the number of grid rows, n=1,⋯,N, and m represents the number of grid columns, m=1,⋯,M.Substituting equations (8) and (9) into equation (6) gives
(10)δTxx=5g11−g2216s0.The gravity gradient component of the residual disturbingδTxx in the central region at any point in the world can then be calculated using equation (10).


## 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]
(2)TxxMR,φ,λ=fMR3∑n=2N∑m=0nC¯∗nmcosmλ+S¯nmsinmλPnmxxsinφ,where TxxM is the horizontal east-west component of disturbed gravity gradient, R is the average radius of the Earth, φ,λ represents the geocentric latitude and longitude of the ocean calculation point, f is the geocentric gravitational constant, M is the mass of the earth, n and m represent the truncation degrees and orders, C¯n,m∗ and S¯n,m are the coefficients of the gravitational field model of the Earth, and P¯nmsinφ 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 Pnmxxsinφ needs to be de singular. The expression after desingularity is as follows [11].
(3)Pnmxxsinφ=12nn−1n+1n+22P¯n,2sinφ−12n+1n+2P¯n,0sinφm=0,14n−2n−1n+2n+3P¯n,3sinφ−143n+1n+2P¯n,1sinφm=1,14n−3n−2n+3n+4P¯n,4sinφ−12n2+3n−2P¯n,2sinφ+142nn−1n+1n+2P¯n,0sinφm=2,14n−m−1n−mn+m+1n+m+2P¯n,m+2sinφ+14n+m−1n+mn−m+1n−m+2P¯n,m−2sinφ−12n2+3n−m2+2P¯n,msinφm>2.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δΔg of the surrounding points. δΔg=Δg−ΔgM, where Δg represents the actual gravity anomaly and ΔgM 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]
(4)δTxx=R4π∬σ1r2cosφ2∂2Sr,ψ∂ψ2∂ψ∂λ2+∂Sr,ψ∂ψ∂2ψ∂λ2−tanφr2∂Sr,ψ∂ψ∂ψ∂φ+1r∂Sr,ψ∂rδΔgdσ,where
(5)∂Sr,ψ∂r=−t2R1−t2D3+4D+1−6D−tcosψ13+6ln1−tcosψ+D2,∂Sr,ψ∂ψ=−t2sinψ2D3+6D−8−31−tcosψ−DDsinψ−3ln1−tcosψ+D2,∂2Sr,ψ∂r2=t2R21−tcosψ31−t2D5−4D3−1+t2D3−10D−18D+2−3tcosψ15+6ln1−tcosψ+D2,∂2Sr,ψ∂ψ2=−t2cosψ2D3+6D−8−31−tcosψ−DDsin2ψ−3ln1−tcosψ+D2+t3sin2ψ6D5+6D3+3D−1D2sin2ψ+3D+1D1−tcosψ+D−31−tcosψ−DDsin2ψ2cosψtsin2ψ+1D2,∂2Sr,ψ∂r∂ψ=t3sinψR31−t2D5+4D3+6D−13−6ln1−tcosψ+D2+6cosψD+1D1−tcosψ+Dt=Rr,D=1−2tcosψ+t21/2ψ=cos−1sinφsinφ′+cosφcosφ′cosλ′−λ.Using equations (4) and (5), the gravity gradient component of the residual disturbing δTxx 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 isx,y,0 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 s0=Rψ0, where ψ0 is the angular distance between the center of the central region and the edge. After gridding the central area, the radius s0 can be approximated to Δx×Δy.If the rectangular coordinatex,y 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 s,α, where s=x′−x2+y′−y2, x and y are the coordinates of the calculation point, and x′ and y′ are the coordinates of the flow point. Then, it is easy to obtain x=scosα, y=ssinα, and dx′dy′=sdsdα. Therefore, the residual disturbing gravity gradient in the central area can be obtained as follows:
(6)δTxx=12π∫02π∫0s0δΔgs23cos2α−1dsdα.The Taylor series expansionδΔg at the calculation point in equation (6) can then be obtained:
(7)δΔg=δΔgP+xg1+yg2+12x2g11+2xyg12+y2g22+⋯,where g1=∂Δg/∂x, g2=∂Δg/∂y, g11=∂2Δg/∂x2, g2=∂2Δg/∂y2, and g12=∂2Δg/∂x∂y.Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form:
(8)δΔg=δΔgP+sg1cosα+g2sinα+s22g11cos2α+2g12sinαcosα+g22sin2α+⋯.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 g1, g2, g11, g12, and g22.Using the above conversion relationship between rectangular and polar coordinates, equation (7) can be converted into polar coordinates in the following form.
(9)g1n,m=δΔgn+1,m−δΔgn−1,m2Δx,g2n,m=δΔgn,m+1−δΔgn,m2Δy,g11n,m=g1n+1,m−g1n−1,m2Δx,g12n,m=g1n,m+1−g1n,m−12Δy,g22n,m=g2n,m+1−g2n,m−12Δy.Here,n represents the number of grid rows, n=1,⋯,N, and m represents the number of grid columns, m=1,⋯,M.Substituting equations (8) and (9) into equation (6) gives
(10)δTxx=5g11−g2216s0.The gravity gradient component of the residual disturbingδTxx 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 Section3.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 functionP¯nmcosθ, it is often necessary to create a two-dimensional array PnM to store P¯nmcosθ. However, the Legendre function P¯nmcosθ 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.Figure 7
Schematic diagram of the array dimensionality reduction.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 Figure8. The computing units in each group are serially computed and share memory. Taking a range of 10°×10° 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 1+100×100/2 to store the Legendre functions. Such an array requires 39.5 KB of memory. If it is directly parallelized, at least 371,059×39.5KB 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 I is 100, only 3.86 MB of memory is required, which is achievable for most computers.Figure 8
Schematic diagram of the grouping calculation.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 Figure9. There are 15 (3×5) thread blocks in a single thread grid, each of which contains 9 (3×3) 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 x and y indexes of a thread in the process grid can be calculated. The calculation formulas are shown in
(11)Idx=blockIdx.x∗blockDim.x+threadIdx.x,Idy=blockIdx.y∗blockDim.y+threadidx.y.Figure 9
CUDA thread index diagram.When calculating various gravity gradient maps, the number of grid points is large. Taking a range of10°×10° 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.Figure 10
Schematic diagram of the one-dimensional coding of a hexagonal grid.Using the parallel computing scheme designed in Sections4.1 and 4.2, the rapid construction of a hexagonal marine gravity gradient reference map can be completed using CPU+GPU hybrid programming.


## 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 functionP¯nmcosθ, it is often necessary to create a two-dimensional array PnM to store P¯nmcosθ. However, the Legendre function P¯nmcosθ 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.Figure 7
Schematic diagram of the array dimensionality reduction.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 Figure8. The computing units in each group are serially computed and share memory. Taking a range of 10°×10° 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 1+100×100/2 to store the Legendre functions. Such an array requires 39.5 KB of memory. If it is directly parallelized, at least 371,059×39.5KB 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 I is 100, only 3.86 MB of memory is required, which is achievable for most computers.Figure 8
Schematic diagram of the grouping calculation.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 Figure9. There are 15 (3×5) thread blocks in a single thread grid, each of which contains 9 (3×3) 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 x and y indexes of a thread in the process grid can be calculated. The calculation formulas are shown in
(11)Idx=blockIdx.x∗blockDim.x+threadIdx.x,Idy=blockIdx.y∗blockDim.y+threadidx.y.Figure 9
CUDA thread index diagram.When calculating various gravity gradient maps, the number of grid points is large. Taking a range of10°×10° 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.Figure 10
Schematic diagram of the one-dimensional coding of a hexagonal grid.Using the parallel computing scheme designed in Sections4.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.Figure 11
Physical characteristic map of the calculation area.
(a)
DTU18 gravity anomaly(b)
Submarine topographyThe 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
Accuracy analysis of performing parallel computing using different data types, where RMS indicates the root mean square.
Data typeMaximum differenceEMinimum differenceEAverage differenceERMSEDouble precision data10-6-10-6<10-12<10-12Single precision data2.72-4.660.180.37Table7 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 2.72E, the minimum error is −4.66E, the average error is 0.18E, and the root mean square (RMS) error is 0.37E. The absolute value of the error calculated with double-precision data is less than or equal to 10−6E, and the average difference and RMS value are less than 10−12E. 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 conditionTxx+Tyy+Tzz=0, 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.Table 8
Numerical characteristics ofTxx+Tyy+Tzz, where T is the disturbing gravity gradient tensor.
MaximumEMinimumEAverageEStdERMSETxx+Tyy+Tzz0.105-0.127<10-40.00750.0091It can be seen from Table8 that the maximum value of Txx+Tyy+Tzz is 0.105E, the minimum value is −0.127E, the standard deviation is 0.0075E, and the RMS error is 0.0091E. 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 Table9. 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.Table 9
Time consumption of the different parallel schemes.
Parallel schemeModelTMcalculation timeModelTM calculation acceleration ratioResidualδT calculation timeResidualδT calculation acceleration ratioTotal timeTotal acceleration ratioSerial computing120.35 h17.92 h1128.27 h1Pure CPU parallelism2.98 h40.3812.11 min39.243.18 h40.33Pure GPU parallelism6.73 h17.8625.32 s1126.386.74 h19.03GPU+CPU hybrid parallel2.98 h40.3825.32 s1126.382.99 h42.90Figure 12
Gravity gradient reference map of the calculation area.
(a)
Tyy(b)
Tzz(c)
Txy(d)
Tzx(e)
TzyTable9 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 the10°×10° 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.Table 10
Calculation time for different resolutions.
LevelSerial computing time (h)Parallel computing time consuming (h)Acceleration ratio11th level28.450.7040.6412th level120.352.9842.9013th level602.9511.8450.92It can be seen from Table10 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 10°×10°, 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.Table10 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.


## 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 Table9. 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.Table 9
Time consumption of the different parallel schemes.
Parallel schemeModelTMcalculation timeModelTM calculation acceleration ratioResidualδT calculation timeResidualδT calculation acceleration ratioTotal timeTotal acceleration ratioSerial computing120.35 h17.92 h1128.27 h1Pure CPU parallelism2.98 h40.3812.11 min39.243.18 h40.33Pure GPU parallelism6.73 h17.8625.32 s1126.386.74 h19.03GPU+CPU hybrid parallel2.98 h40.3825.32 s1126.382.99 h42.90Figure 12
Gravity gradient reference map of the calculation area.
(a)
Tyy(b)
Tzz(c)
Txy(d)
Tzx(e)
TzyTable9 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 the10°×10° 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.Table 10
Calculation time for different resolutions.
LevelSerial computing time (h)Parallel computing time consuming (h)Acceleration ratio11th level28.450.7040.6412th level120.352.9842.9013th level602.9511.8450.92It can be seen from Table10 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 10°×10°, 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.Table10 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.


## 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 Table9. 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.Table 9
Time consumption of the different parallel schemes.
Parallel schemeModelTMcalculation timeModelTM calculation acceleration ratioResidualδT calculation timeResidualδT calculation acceleration ratioTotal timeTotal acceleration ratioSerial computing120.35 h17.92 h1128.27 h1Pure CPU parallelism2.98 h40.3812.11 min39.243.18 h40.33Pure GPU parallelism6.73 h17.8625.32 s1126.386.74 h19.03GPU+CPU hybrid parallel2.98 h40.3825.32 s1126.382.99 h42.90Figure 12
Gravity gradient reference map of the calculation area.
(a)
Tyy(b)
Tzz(c)
Txy(d)
Tzx(e)
TzyTable9 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 the10°×10° 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.Table 10
Calculation time for different resolutions.
LevelSerial computing time (h)Parallel computing time consuming (h)Acceleration ratio11th level28.450.7040.6412th level120.352.9842.9013th level602.9511.8450.92It can be seen from Table10 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 10°×10°, 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.Table10 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 Table11 and Figure 13, and the internal consistency of the gravity gradient tensor is shown in Table 12.Table 11
Numerical characteristics of the global marine gravity gradient.
MaximumEMinimumEAverageEMedianETxx131.43-118.97<10-40.0725Tyy102.45-106.05<10-40.0584Tzz221.83-117.77<10-40.0100Txy78.66-74.72<10-40.0000Tzx112.31-111.50<10-40.0037Tzy103.02-107.41<10-4-0.0064Figure 13
Global hexagonal ocean gravity gradient reference map model.
(a)
Tyy(b)
Tzz(c)
Txy(d)
Txz(e)
TyzTable 12
Numerical characteristics of the globalTxx+Tyy+Tzz.
Maximum (E)Minimum (E)Average (E)Std (E)RMS (E)Txx+Tyy+Tzz0.175-0.168<10-40.000120.00037The results in Tables11 and 12 show that the maximum absolute value of the global full-tensor hexagonal marine gravity gradient reference map model is less than 0.18E and that the standard deviation is less than 0.0002E, 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 10E, the average value is less than 10−6E, the median is near 0E, 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 to10−6E, a maximum absolute value of less than 0.18E, and a standard deviation of less than 0.0002E



---

*Source: 4141572-2023-03-23.xml*