Improving the Reconstruction Accuracy of Marine Gravity Anomaly Encrypted Reference Map Using the New Mean Sea Surface 3-D Correction Method

The reconstruction methods of marine gravity anomaly encrypted reference map are studied in this study, which can improve the accuracy of underwater gravity matching navigation. Firstly, the mean sea surface as the third dimension component was introduced into the semivariance function of the traditional Kriging two-dimensional (2-D) interpolation method for the first time. After that, a new mean sea surface three-dimensional (3-D) correction method was proposed in order to improve the spatial resolution of the marine gravity anomaly reference map and minimize the loss of accuracy. Secondly, encrypting the spatial resolution of the marine gravity anomaly reference map from <inline-formula> <tex-math notation="LaTeX">$4^\prime \times 4^\prime $ </tex-math></inline-formula>, <inline-formula> <tex-math notation="LaTeX">$3^\prime \times 3^\prime $ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$2^\prime \times 2^\prime $ </tex-math></inline-formula> to <inline-formula> <tex-math notation="LaTeX">$1^\prime \times 1^\prime $ </tex-math></inline-formula> as an example, respectively, under the same conditions, the reconstruction accuracy of the new proposed method improved 62.25%, 140.92% and 104.93%, respectively. Thus, compared with the traditional Kriging 2-D interpolation method, we verified that the new proposed method had better reconstruction accuracy. Thirdly, based on the proposed new method, the spatial resolution of the marine gravity anomaly reference map was firstly encrypted from <inline-formula> <tex-math notation="LaTeX">$1^\prime \times 1^\prime $ </tex-math></inline-formula> to <inline-formula> <tex-math notation="LaTeX">$0.25^\prime \times 0.25^\prime $ </tex-math></inline-formula>. After that, transforming the spatial resolution of <inline-formula> <tex-math notation="LaTeX">$0.25^\prime \times 0.25^\prime $ </tex-math></inline-formula> into <inline-formula> <tex-math notation="LaTeX">$1^\prime \times 1^\prime $ </tex-math></inline-formula>, compared with the original model data of <inline-formula> <tex-math notation="LaTeX">$1^\prime \times 1^\prime $ </tex-math></inline-formula>, its the final RMSE is <inline-formula> <tex-math notation="LaTeX">$4.143\times 10^{-1}$ </tex-math></inline-formula> mGal, which further verified the validity and reliability of the new mean sea surface 3-D correction method in application.


I. INTRODUCTION
Underwater gravity matching navigation generally uses the gravity feature information such as gravity anomaly or gravity gradient to correct the accumulated errors of inertial navigation system over time, and then achieve the purpose of improving the positioning accuracy of underwater vehicles [1], [2]. At present, the main research direction of improving the accuracy of gravity matching navigation is to solve the key technical problems, such as high-precision gravity measurement systems [3], [4], gravity matching positioning methods [5]- [9], high-precision and high-resolution The associate editor coordinating the review of this manuscript and approving it for publication was Jenny Mahoney. marine global gravity anomaly reference maps [10], [11]. Among them, the reference map possesses unique advantages include stability, anti-interference, passivity and irregular distribution when the gravity information of underwater vehicle is obtained by gravimeter [12]. However, the improvement of accuracy and spatial resolution of the marine gravity anomaly reference map is an inverse relationship (that is, with the increase of the spatial resolution, the accuracy will be reduced). Therefore, it is necessary to find a balance between spatial resolution and accuracy for the current demand of high-spatial resolution and high-precision the reference map, so that it can meet the requirements of high-precision underwater gravity matching navigation. Nevertheless, how to propose an reconstruction method to improve spatial resolution and minimize the accuracy loss is still is an important research direction in the field of underwater navigation.
As a carrier of storing regular grid gravity information, spatial resolution and precision information are included in the marine gravity anomaly reference map, which are the key factors to ensure the location accuracy of underwater gravity matching navigation. Therefore, the research on the reconstruction methods of marine gravity anomaly reference map has attracted a lot of scholars to study. In 1998, Joachim and Bian firstly proposed a fast Fourier transform method to approximate the local marine gravity field [13]. In 2000, Guo et al. suggested a non-uniform B-spline least square approximation method and verified that the method had good accuracy and stability [14]. In 2001, Yang et al. proposed a comprehensive approximation method combining function model and statistical model, and came to a conclusion that the proposed method had better approximation effect than a single method [15]. In 2007, Cheng and Wu proved that Kriging interpolation method can better realized the encrypted of marine gravity anomaly reference map [16], [17]. In 2009, the traditional Kohonen artificial neural network method was improved by sarzeaud et al., and they verified that the improved artificial neural network method and the traditional Kriging two-dimensional (2-D) interpolation method were quite accurate, but the calculation time was shorter [18]. In 2010, combined with the continuous characteristics of the geophysical field, the Koons surface was introduced in the reconstruction of the marine gravity anomaly reference map by Li et al. [19]. In 2011, Zhao et al. improved a method based on the traditional BP (back propagation) neural network, which had higher accuracy and less time compared with the traditional BP method, but the influence of reference point distribution and irregular shape on the reconstruction accuracy was not easy to be determined [20]. In 2012, Tong et al. further developed the 2-D Gaussian spline function method and proved this method also met the matching requirements compared with B-spline function approximation method [21]. In 2015, the compressed sensing theory was taken into account in the reconstruction of airborne gravity profile by Yang et al. [22]. In 2016, Ren et al. conducted experiments of reconstruction of marine gravity anomaly reference map with the current commonly interpolation methods such as inverse distance weighting method, polynomial fitting method, trigonometric network method, radial basis function method, improved Shepard method and Kriging interpolation method, and verified that Kriging interpolation method had smaller error mean square deviation after comparison and analysis [23]. In 2017, A Rossi and Perracchione proposed a radial basis function method for the first time and verified that the method had good accuracy and applicability in encrypting marine gravity anomaly reference map based on simulated experiment [24]. In 2018, Xu et al. proposed a gravity data reconstruction method based on Nonequispaced Fourier transformin framework of compressed sensing theory and verified that the proposed method had better approximation effect compared with the traditional Kriging 2-D interpolation method and the minimum curvature radius method even with few sample data [25]. At present, in other fields, some 3D interpolation techniques are also used to improve some interpolation effects [26]- [29]. Our research team of navigation and detection based on the information of aerospace-aeronauticsmarine integration (http://www.qxslab.cn/ndia/index.html) of Qian Xuesen Laboratory of Space Technology had carried out preliminary exploration research on the theory, method, and key technology of constructing the high-precision and high-resolution marine gravity anomaly reference map for improving underwater gravity matching navigation accuracy based on the GNSS-R altimetry constellation principle [1], [2], [8], [9], [30].
Different from previous studies, based on the traditional Kriging interpolation method, we introduced the third dimension variable (Mean sea surface) to construct a new method. According to the principle of statistics, the new proposed method can improve the prediction accuracy of unknown points by increasing the relevant information to optimize the prediction weight value. Under the same conditions, compared with the traditional Kriging 2-D interpolation method, the proposed new method improves the spatial resolution of the marine gravity anomaly reference map and the accuracy loss is greatly reduced simultaneously.

II. PRINCIPLE OF THE NEW MEAN SEA SURFACE 3-D CORRECTION METHOD (NMSST-DCM)
The traditional Kriging 2-D interpolation method was an optimal interpolation method named after D.G.kriging in France [31], [32]. Because the relative position and relationship between the known point and the point to be estimated was comprehensively considered, the method had unbiased and optimal properties compared with several commonly used interpolation methods [33]. In geological statistics, the traditional Kriging 2-D interpolation method was the most important and basic estimation method, and the spatially distributed data of the method was usually stable when performing encrypted experiments [34]. Additionally, based on the traditional Kriging 2-D interpolation method, the global and local characteristics of the fluctuation of the marine gravity field could be more comprehensively reflected, and this method had become a powerful interpolation technology [35], [36]. We compared the interpolation effect of several commonly used interpolation algorithms. Table 1 shows the statistical results of different interpolation algorithms.
The statistical results in Table 1 show that the Kriging has better interpolation accuracy than other commonly used interpolation algorithms. However, the traditional Kriging 2-D interpolation method still has a great loss of accuracy in the reconstruction of marine gravity anomaly encrypted reference map. Therefore, based on the traditional Kriging 2-D interpolation method, the third component mean sea surface (MSS) is introduced in the semivariance function for the first time, and the new mean sea surface 3-D correction method (NMSST-DCM) is proposed for minimize the accuracy loss. The proposed new method optimizes the weights of samping points by increasing the new variable information, and then uses the new weights to predict the function values of the point to be estimated in space. However, in the reconstruction process of the marine gravity anomaly reference map, the MSS plays a key role in the NMSST-DCM. Generally speaking, there is a certain relationship between geoid and mean sea surface. Moreover, the geoid, as an information quantity to described the characteristics of the earth and the gravitational equipotential plane, was a continuous smooth closed surface which describes the basic shape, structure and a physical gravity field element of the earth [37], [38]. And then, according to the theory of the earth's gravity field, the marine geoid in the sea was a basic concept that describes the shape of the earth's part (especially the sea part), which reflected the inhomogeneity of the earth's internal material density and the abnormal characteristics of the marine gravity field [39], [40]. As the improvement of spatial positioning technology, the research and application fields of altimetry technology have been further expanded and deepened. In general, the marine geoid was mainly obtained by satellite altimetry technology, and this technology provides the marine geoid with many advantages such as uniform elevation reference, high-precision and high-spatial resolution for the sea area [37]. The indirect application of this technology was to carry out the marine gravity acquisition, which mainly consisted of two parts: the acquisition of the marine geoid and the inversion of marine gravity field based on relevant algorithms [43]- [45].
Although in theory, the MSS was not the marine geoid, but their influence mainly came from the of the dynamic sea surface topography [41]. Therefore, we can know that there is also an inseparable relationship between MSS and gravity anomaly in numerical value. Figure 1 illustrates the relationship between marine geoid and MSS. In Figure 1, the Instantaneous Sea Level on the surface of the sea is represented by the blue wavy line, the black dotted line in the middle refers to the MSS, and the wave line below MSS is the N (marine geoid). Moreover, some scholars had produced qualitive result that the sum of part of Sea Surface Topography (SST) and the N was the approximation of MSS, and the difference between the two quantities was mainly the influence of part of Sea Surface Topography [46]. Therefor, we can obtain the mean sea surface P (H ) of a given (measured MSS) point P by constructing a linear combination of the measurements as: where H is the MSS, ζ is the SST, and ζ is part of the SST. However, the SST is usually classified as Sea Surface Height residual according to some factors, such as the complexity of the SST, the large regional error of the global SST model, the small impact of the part SST value ζ ' and the immature local applicability model [46]. Therefore, under the condition of neglecting the dynamic sea surface topography, the value of H variable can be converted into the value of N variable. Finally, the linear relationship between the N and the MSS can be determined as follows: where the H is set as the third dimension variable except longitude and latitude of the proposed new method. According to the principle of statistics, this study introduced the H which was determined by Equation (1) and Equation (2) to increase the information. Thus, the weights of samping points are optimized by increasing the known information, and the gravity anomaly values of unknown points can be better predicted. And then, the size of the rectangular local calculation area randomly selected is defined as a A×B grid, we assume that the corresponding marine gravity field information (including n points in total) is obtained such as longitude, latitude, mean sea surface and gravity anomaly of local sea area in space. In geostatistics, the relationship of spatially distributed data is usually set to stationarity. Thus we assume that the correlation function at any two points only depends on the distance between them and the mathematical expectation is constant and has nothing to do with its location. In this case, we can easily obtained the following system of equation: where the X , Y and H are the longitude, latitude and MSS after standardization, respectively, and G(X , Y , H ) is the marine gravity anomaly at actual location of corresponding point. And then, the h X , h Y and H correspond to the coordinate changes between adjacent grid points on each coordinate axis, respectively, and E(G(X, Y, H )) is the mathematical expectation of gravity anomaly G(X , Y , H ), and K refers to arbitrary constant. If the condition of isotropy is satisfied, the variance relationship can be obtained as: where Var is the sign of variance, D is the Euclidean distance between any two points, respectively, and γ (D) is the semivariance function describing the spatial correlation of random functions. If the condition of isotropy is not satisfied, the semivariance function should be in general form . What's more, as a basic tool of the traditional Kriging 2-D interpolation method, the semivariance function have some advantages such as reflecting the spatial structural change, controling random change and indicating error information of variables, which are the important factors to accurately obtain the prediction point values.
In addition, in order to keep the dimension of the coordinate axis vector of the 3-D coordinate system consistent, this study firstly standardizes the coordinate information of longitude, latitude and mean sea surface of all points. In this case, by making different dimensions consistent, the calculation can be facilitated, and the error is reduced. To standardize the sample data more efficiently, this study adopts the Z-score standardization method which is commonly used and convenient at present. Thus, it is easy to obtain standardized equation system as follows: where x, y and H correspond to the longitude, latitude and mean sea surface coordinate information of the grid points (non-standardized), respectively, and µ and σ are the mean value and the standard deviation of the corresponding sample points (non-standardized), respectively.
In addition, according to the spatial correlation theory of the first law of geography, the semivariance function is determined by the principle that the spatial similarity leads to the attribute similarity. As a result, the spatial similarity can be expressed by the Euclidean distance between two points in space [47], [48]. However, for the traditional Kriging 2-D interpolation method, the Euclidean distance of any two points is obtained as follows: where X and Y are the standardized longitude and latitude value respectively, and d is the Euclidean distance between any two points in space. For the traditional Kriging 2-D interpolation method, the MSS is not introduced in calculation of Euclidean distance. In this case, for the traditional Kriging 2-D interpolation method, the constraint relationship variable d is formed by latitude and longitude, which affects the weights accuracy of samping points. According to the principle of statistics, the more accurate weights will not be obtained. Therefore, due to the accuracy of unknown points mainly depends on the weights, it is the key to improve the accuracy of points to be estimated by optimizing the weights of samping points. In order to optimize the weights of the samping points, the H is introduced as the third variable factor based on the longitude and latitude coordinates. Therefore, the 2-D coordinate system of the traditional Kriging 2-D interpolation method has transform into the 3-D coordinate system [49] of the proposed new method. And then the Euclidean distance equation of arbitrary point P under the new coordinate system can be obtained: where H is the standardized mean sea surface and D is the Euclidean distance of arbitrary two points in 3-D space. Figure 2 illustrates the diagram of transformation from 2-D to 3-D coordinate system. Figure (a) is a 2-D plane interpolation, while Figure (b) is a 3-D coordinate system interpolation. As shown in Figure 2, the change at gravity anomaly of arbitrary point P in the 2-D is from at the prediction point VOLUME 8, 2020 in the 2-D coordinate system to the more accurate gravity anomaly value at the prediction point in the 3-D coordinate system. Due to the errors in the measurement of sampling points, we sure that the gravity anomaly P in the 2-D coordinate system hardly coincides with the value P in the 3-D coordinate system. Moreover, according to the relationship between the prediction point and the expectation of the difference between the known points, the following relationship can be obtained: where G(X P , Y P , H P ) is the gravity anomaly value at the location of (X P ,Y P ,H P ) the carrier, the λ i (i = 1, 2, 3, . . . , n) is the weights of samping points. Moreover, according to Equation (8), we usually do not need to calculate the value of K. And then, the weights are a set of optimal coefficients to meet the minimum difference between the estimated and the actual values of the marine gravity anomaly in space, which play a key role in the accuracy of the prediction of unknown points.
Since the NMSST-DCM has the same characteristics of unbiased and optimal as the traditional Kriging 2-D interpolation method, the weights of samping points can be obtained from the satisfied relationship as follows [50]: Additionally, according to the condition that the sum of Equation (9) is 1, the error variance equation of predicted values can be further calculated as follows: where γ ((X i , Y i , H i ), (X j , Y j , H j )) denotes the value of the semivariance function between the point (X i , Y i , H i ) and the point (X j , Y j , H j ). In order to get the best interpolation result, it is necessary to ensure the minimum error variance under unbiased conditions. Hence, it is also necessary to construct the objective function as follows: where µ P is the Lagrange multiplier at the location of P point. For the calculation of objective function, the partial derivative of the Equation (11) is calculated, and then the partial derivative is set to zero, so the following equation system is introduced: where γ ((X n , Y n , H n ), (X P , Y P , H P )) is the value of the semivariance function between the samping points and corresponding point to be estimated. Therefor, the weights λ i and Lagrange multiplier µ P can be determined according to the matrix Equation (13). If all the variables satisfy the above stationary hypothesis, according to a set of gravity anomaly values given H 2 ),... G(X n , Y n , H n ), we can get the prediction equation as follows: where G(X P , Y P , H P ) is the gravity anomaly value at location of P point in space, G(X i , Y i , H i ) is the marine gravity anomaly value at location i in space. According to Equation (14), the values at the unknown grid points in the prediction area can be determined.

A. THE DATA AND MODEl
In order to compare the accuracy between the proposed new method and the traditional Kriging 2-D interpolation method, zone of the south China sea (longitude 114 • E −115 • E, latitude 15 • N −16 • N) was selected for test. Moreover, according to the analysis of the characteristic parameters of gravity field such as roughness, standard deviation of gravity anomaly, standard deviation of slope and entropy of gravity anomaly difference, the adaptability of this area was better. Figure 3 illustrats Satellite image and partial enlarged map of the study area. The DTU13 global marine gravity field model (1 ×1 ) and Mean Sea Surface (MSS) model (1 × 1 ) from Andersen team of Danish University of Science and Technology (https://www.space.DTU.dk/) were used as original calculation data. In addition, this study obtained gravity anomaly sample data (3600 values) whose range were from −4.040 × 10 1 to 1.103 × 10 2 mGal, and then, the MSS sample data (3600 values) were from 1.765 × 10 1 to 2.422 × 10 1 m. Figure 4(a) and (b) depict 2-D and 3-D marine gravity anomaly reference maps with a spatial resolution of 1 ×1 , respectively, which reflect the change trend of marine gravity fielid in the selected region, such as the more intense regionin the north and the gentle region in the south.
However, because the earth is not a medium density equilibrium sphere and the different types of medium distribution in different regions lead to the difference between topography and gravity anomaly values, the local variation trend in Figure 4 is not completely consistent with the remote sensing image in Figure 3.
In this study, the original data of marine gravity anomaly and MSS with spatial resolution of 1 ×1 in the selected area were sparse processed. Thus, the the global marine gravity anomaly and MSS test data with spatial resolution of 2 ×2 and 4 ×4 were obtained, respectively.

B. THE EVALUATION OF RESULTS
For the experimental verification, the reconstructed marine gravity anomaly encrypted reference maps are compared with DTU13 model in the same area, and then the root mean square error (RMSE) of the difference between the two points is taken as the accuracy evaluation standard, and the relationship is as follows: (15) where G(X i , X i , H i ) and G(X i , X i , H i ) are the predicted values and original model values of gravity anomaly at grid points, respectively.

C. THE SELECTION OF THE BEST APPROXIMATION SEMIVARIANCE FUNCTION MODEL
Since the function relationship between Euclidean distance and semivariance function can not be determined directly, it is necessary to determine change trend according to the spatial VOLUME 8, 2020  distribution of gravity anomaly values at the grid points in the selected region. Thus, in the reconstruction experiments of marine gravity anomaly encrypted reference maps, some appropriate and simplified semivariance models are selected such as linear, quadratic, exponential and logarithmic and other functional relationship. For the reconstruction accuracy of the marine gravity anomaly encrypted reference maps, the appropriate semivariance function model with the best approximation effect is very important for the NMSST-DCM and the traditional Kriging 2-D interpolation method. In order to select the semivariance function model better, this study encrypts the spatial resolution of the marine gravity anomaly reference maps as an example from 2 ×2 to 1 ×1 .
Additionally, making a comparative analysis of the seven semivariance function models such as Exponential model - [53]. The reconstruction accuracy results of the marine gravity anomaly encrypted reference maps are shown in Table 2 and Table 3. Table 2 shows the reconstruction results of the traditional Kriging 2-D interpolation method based on seven semivariance function models. According to the statistical results of Table 1, based on previous traditional method, the seven semivariance function models selected had obvious differences in  the reconstruction accuracy of marine gravity anomaly reference maps. When maximum, minimum and mean errors did not change much, the traditional Kriging 2-D interpolation method based on G-EXPM and SPLM had better fitting effect and the reconstruction accuracy was 2.744 × 10 −1 mGal and 2.771 × 10 −1 mGal, respectively. Table 3 shows the reconstruction statistical results of the NMSST-DCM based on seven semivariance function models. According to Table 3, when conducting the reconstruction experiments, the new proposed method based on the G-EXPM and the CUBM performed better, but the error of maximum, mean and minimum values of the NMSST-DCM based on the G-EXPM were relatively less. Moreover, the reconstruction accuracy of the new proposed method based on the G-EXPM reaches 1.339 × 10 −1 mGal, which was superior to other semivariance function models.
Additionally, according to the statistical results of Table 2 and Table 3, the NMSST-DCM based on seven semivariance function models had obvious changes in the maximum, minimum and mean values of the reconstructed marine gravity anomaly reference maps compared with the traditional Kriging 2-D interpolation method. In Table 2 and  Table 3, the reconstruction accuracy of the proposed new method had been improved to varying degrees compared VOLUME 8, 2020 with the traditional Kriging 2-D interpolation method. The approximation result was the the new proposed method based on the G-EXPM, whose reconstruction accuracy reached better.
Based on seven semivariance function models, Figure 5 shows the reconstruction accuracy comparison histogram of the marine gravity anomaly encrypted reference maps by the NMSST-DCM and the traditional Kriging 2-D interpolation method. The gray rectangle box represented the reconstruction accuracy based on the NMSST-DCM, the black rectangle box illustrated the reconstruction accuracy based on the traditional Kriging 2-D interpolation method, and the arrow showed the percentage of accuracy improvement of the the new proposed method compared with the traditional Kriging 2-D interpolation method. By showing statistical results of Figure 5, when the spatial resolution of the marine gravity anomaly encrypted reference maps were reconstructed from 2 ×2 to 1 ×1 , the NMSST-DCM improved the reconstruction accuracy in different degrees compared with the traditional Kriging 2-D interpolation method. Comparing with the traditional Kriging 2-D interpolation method, the result indicated that the reconstruction accuracy of the NMSST-DCM based on the CUBM and the G-EXPM was improved by 145.93% and 104.93%, respectively. However, for the final approximation effect, the comprehensive performance of the NMSST-DCM based on the G-EXPM was better. Therefore, the the new proposed method based on the G-EXPM was the best for the prediction of marine gravity anomaly values in space compared with other semivariance function models.

D. THE VERIFICATION OF THE NMSST-DCM
In addition, in order to further verify the validity and reliability of the NMSST-DCM based on the G-EXPM, the reconstruction experiments of a variety of marine gravity anomaly reference maps were conducted with different spatial resolutions. The experimental process ranged from 4 ×4 to 2 ×2 and 1 ×1 , 3 ×3 to 1 ×1 and 2 ×2 to 1 ×1 . Figure 6 depicts marine gravity anomaly reference maps of the reconstruction experiments based on the traditional Kriging 2-D interpolation method. At the same time, Figure 7 shows the error approximation diagrams of the 3-D corresponding to Figure 6. Firstly, according to the comparison of (a), (b), (c) and (d) in Figure 6, we proved that there were more regional differences such as northwest region and other individual regions for reconstruction of the different spatial resolutions. Then, compared with the error approximation diagrams (a), (b), (c) and (d) of the 3-D in Figure 7, we verified that the approximation error gradually reduced as the improvement of the spatial resolution of the original model data. For example, the error approximation ranged about from −30 to 15 mGal in Figure 7 (a), and that of Figure 7 (d) was about from −5 to 3 mGal.
At the same time, under the same conditions and sea area, the reconstruction experiments of the marine gravity anomaly reference maps were conduced based on the NMSST-DCM. Figure 8 illustrates marine gravity anomaly encrypted reference maps of the reconstruction experiments based on the proposed new method. Figure 9 depicts the error approximation diagrams of the 3-D corresponding to   Figure 6 and Figure 8 indicated some regional shapes of different sizes and linear trend. Moreover, the boundary errors also had a great influence on the results which were shown in Figure 6 and Figure 8. Certainly, if the boundary errors of marine gravity anomaly reference maps were not taken into account, the accuracy changes would be smaller between the proposed new method and the traditional Kriging 2-D interpolation method. However, in order to fairly compare both methods, we believed that all points had to be considered in this study. In other words, the advantages of NMSST-DCM became more obvious in the boundary area. Morever, comparing the error approximation Figure 7 based on the traditional Kriging 2-D interpolation method, Figure 9 indicated that the error approximation absolute values of the marine gravity anomaly based on the NMSST-DCM were smaller and the overall accuracy was higher. Obviously, the error approximation range of Figure 7 (a) was about from −30 to 15 mGal, while that of the corresponding Figure 9 (a) was about −4 to 12 mGal, and the error approximation ranged about from −5 to 3 mGal in Figure 7 (d), while that of the corresponding Figure 9 (d) was about −1.5 to 1.5 mGal. Table 4 lists the statistical reconstruction results from Figure 6 and Figure 8. According to Table 4, for reconstruction experiments of the same spatial resolution based on the NMSST-DCM and the traditional Kriging 2-D interpolation method, the errors including maximum and minimum errors were gradually reduced, and the higher the spatial resolution of the original marine gravity anomaly reference maps, the higher accuracy of the reconstruction. Moreover, the absolute value of the mean error was also gradually reduced with the increase of spatial resolution, for example, 2 ×2 to 1 ×1 . And then, comparing the NMSST-DCM and the traditional Kriging 2-D interpolation method, we could prove that the proposed new method was better than the traditional Kriging 2-D interpolation method. For example, compared with the traditional Kriging 2-D interpolation method, when the spatial resolution of the marine gravity anomaly reference maps were encrypted from 4 ×4 , 3 ×3 and 2 ×2 to 1 × 1 , respectively, the overall accuracys based on the NMSST-DCM improved 62.25%, 140.92% and 104.93%, respectively. Although the loss of reconstruction accuracy was all decreasing, the range of the reduction was different. As a result, compared with the traditional Kriging 2-D interpolation method, the proposed new method did not show a strict increasing relationship. But, it still had a great improvement in reducing the accuracy loss. Therefore, compared with the traditional Kriging 2-D interpolation method, the proposed new method could reduce the loss of reconstruction accuracy to a large extent. In other words, we further proved the effectiveness and reliability of the NMSST-DCM in reducing the accuracy loss.

IV. THE DISCUSSION AND APPLICATION OF THE NMSST-DCM A. THE ACQUISITION OF THE SIMULATED MSS DATA
To date, mean sea surface of the best spatial resolution (i.e. 1 ×1 ) is published from the authoritative Andersen team of Danish University of Science and Technology, which greatly limits the reconstruction of marine gravity anomaly reference map of higher spatial resolution. Due to the higher spatial resolution demand of underwater gravity matching navigation, this paper will continue to study higher spatial resolution encrypted marine gravity anomaly reference map in depth. In this study, based on traditional Kriging 2-D interpolation method, the MSS data with spatial resolution of 4 ×4 , 3 ×3 and 2 ×2 was interpolated as 1 ×1 simulated MSS data, respectively. Figure 10 shows the MSS reference maps and the corresponding error approximation scatter maps based on the traditional Kriging 2-D interpolation method, respectively. In Figure 10   with a small range of MSS changes. There were some local differences, which were mainly concentrated in the northeast and southwest. In Figure 10 (b), (d) and (f), we could conclude that the interpolation accuracy is 5.070 ×10 −2 m, 3.300 × 10 −2 m and 1.560 × 10 −2 m, respectively. As the improvement of the spatial resolution of the original MSS data, the accuracy loss of the simulated MSS was also decreasing, which further indicated the stability of the simulated MSS used in NMSST-DCM.

B. THE APPLICATION VERIFICATION OF THE NMSST-DCM
In order to verify the effect of reducing the loss of reconstruction accuracy by NMSST-DCM, the simulated MSS data was used for experiments in the selected area. Under the same conditions, conducting the experiments of reconstruction, such as 4 ×4 to 2 ×2 and 1 ×1 , 3 ×3 to 1 ×1 , 2 ×2 to 1 ×1 . Figure 11 depicts the marine gravity anomaly encrypted reference maps based on the NMSST-DCM using the simulated MSS data. In addition, comparing Figure 11 with the same spatial resolution marine gravity anomaly encrypted reference maps Figure 6, we verified regional size differences and linear trend still existed, such as the northeast and northwest regions. Figure 12 illustrates the error approximation diagrams of the 3-D corresponding to Figure 11. By comparing Figure 12 with Figure 7, we verified that the error approximations also showed different ranges of accuracy variation. For example, the error approximations in Figure 7 (a) and Figure 7 (d) were about from −30 to 15 mGal and −5 to 3 mGal, respectively. But the error approximations in Figure 12 (a) and Figure 12 (d) ranged about from −4 to 12 mGal and −2 to 3 mGal, respectively. Therefore, based on the NMSST-DCM using simulated MSS, the advantage of the proposed new method became also clearer than the traditional Kriging 2-D interpolation method in the boundary area when all the points were taken into account. Table 5 provides the statistical experimental results of two methods. According to Table 5, compared with results of the  traditional Kriging 2-D interpolation method, in the same sea area and under the same conditions, we proved that the absolute values of the maximum and minimum errors were reduced to different degrees and mean errors changed little. When encrypting the spatial resolution of the marine gravity anomaly reference maps from 4 ×4 , 3 ×3 and 2 ×2   to 1 ×1 , as an example respectively, Figure 13 provides improvement percentage histogram. According to Figure 13, compared the NMSST-DCM using simulated MSS data with the traditional Kriging 2-D interpolation method, the overall reconstruction accuracy improved 62.42%, 132.90% and 64.11%, respectively. Although the loss of reconstruction accuracy based on the proposed method using simulated MSS data increased compared with the NMSST-DCM using actual MSS data, it still reduced the accuracy loss to a large extent compared with the traditional Kriging 2-D interpolation method. As a result, compared with the traditional Kriging 2-D interpolation method, the proposed new method still reduced the loss of reconstruction accuracy of encrypting marine gravity anomaly reference maps to a large extent. Therefore, based on the NMSST-DCM using simulated MSS data, we verified that the proposed new method had far-reaching research significance for the further reconstruction of higher spatial resolution marine gravity anomaly reference maps.

C. FURTHER APPLICATION OF CONSTRUCTING MARINE GRAVITY ANOMALY REFERENCE MAP
In addition, in order to verify the effectiveness of the proposed new method in reducing the loss of reconstruction accuracy, this study encrypted the marine gravity anomaly reference map with spatial resolution of 1 ×1 to 0.25 ×0.25 as an example. The Figure 14 (a) and (b) depict the 2-D and 3-D marine gravity anomaly reference maps of spatial resolution of 0.25 ×0.25 , respectively. Nevertheless, in order to have a qualitative evaluation of the reconstructed results, the nearest neighbor distance method [54] with faster operation speed and more convenient was selected for restoring test in this paper. Therefore, compared with the original model data of 1 ×1 , when the marine gravity anomaly encrypted reference map with spatial resolution of 0.25 ×0.25 was restored to 1 ×1 based on the nearest distance method, its RMSE was 4.143 × 10 −1 mGal, and its 3σ error was 1.243 × 10 0 mGal (3σ is about 99.74%). Since the new procedure of high spatial resolution to low spatial resolution, the new error from the interpolation process was introduced and the actual RMSE was less than 4.143 × 10 −1 mGal. As a result, the NMSST-DCM provides the feasibility and effectiveness for the future reconstruction of higher spatial resolution marine gravity anomaly reference maps. VOLUME 8, 2020

V. CONCLUSIONS
The global high-precision and high-spatial resolution marine gravity anomaly reference map is one of the key elements to improve the accuracy of the underwater vehicle gravity matching navigation, which determines the navigation security and the precise positioning of the underwater vehicle. However, in the field of marine gravity field, the global marine gravity field model (the highest spatial resolution of 1 ×1 ) published could not meet the high-precision navigation requirements of underwater vehicles [55], [56]. Therefore, a new mean sea surface 3-D correction method (NMSST-DCM) was proposed to improve the accuracy of gravity matching navigation.
(1) In this study, we constructed the new mean sea surface 3-D correction method (NMSST-DCM). Firstly, based on the principle of the traditional Kriging 2-D interpolation method, we verified that the best approximation semivariance function model was Generalized Exponential model from seven semivariance function models. Then, the third variable MSS (i.e. by increasing the variable information) was introduced to construct the 3-D coordinate system, and the NMSST-DCM was proposed.
(2) Under the excellent-suitability area of the South China Sea, the results showed that the NMSST-DCM could improve the spatial resolution of the marine gravity anomaly reference map and reduced the accuracy loss. Among them, compared with the traditional Kriging 2-D interpolation method, when encrypting the spatial resolution of the marine gravity anomaly reference maps from 4 ×4 , 3 ×3 and 2 ×2 to 1 ×1 as an example, respectively, the loss of reconstruction accuracy based on the NMSST-DCM was significantly reduced. The reconstruction accuracy of the NMSST-DCM improved 62.25%, 140.92% and 104.93%,respectively.
(3) Exploration and research: due to the lack of MSS with spatial resolution better than 1 ×1 , the MSS with spatial resolution of 0.25 ×0.25 was obtained by using the traditional Kriging 2-D interpolation method. Moreover, compared with the original model data of 1 ×1 , the NMSST-DCM encrypted the spatial resolution of marine gravity anomaly reference map from 1 ×1 to 0.25 ×0.25 , and then used the nearest neighbor distance method to restore the spatial resolution of 1 ×1 , its final RMSE is 4.143×10 −1 mGal. Thus, the proposed new method not only obtained the higher spatial resolution marine gravity anomaly reference map, but also reduced the accuracy loss greatly, which provided an effective basis for improving the gravity matching navigation accuracy of underwater vehicles. In other words, the study provides a new method support for constructing a higher spatial resolution marine gravity anomaly reference map.