Improved the Accuracy of Seaﬂoor Topography from Altimetry-Derived Gravity by the Topography Constraint Factor Weight Optimization Method

: Gravity geologic method is one of the important to derive seaﬂoor topography by using altimetry-gravity, and its committed step is gridding of regional gravity anomaly. Hence, we proposed a topography constraint factor weight optimization (TCFWO) method based on ordinary kriging method. This method fully considers the inﬂuence of topography factors on the construction of regional gravity grid besides horizontal distance. The results of regional gravity anomaly models constructed in the Markus-Wake seamount area show that the TCFWO method is better than ordinary kriging method. Then, the above two regional gravity models were applied to invert the seaﬂoor topography. The accuracy of derived topographic models was evaluated by using the shipborne depth data and existing seaﬂoor topography models, including ETOPO1 and V19.1 model. The experimental results show that the accuracy of ST_TCFWO (seaﬂoor topography model inverted by TCFWO method) is better than ST_KR (seaﬂoor topography model inverted by kriging method) and ETOPO1 model. Compared with the ST_KR, the accuracy of the ST_TCFWO has improved about 26%. In addition, the accuracy of seaﬂoor topography is affected by the variation of depth, the distribution of control points and the type of terrain. In different depth layers, the ST_TCFWO has better advantages than ST_KR. In the sparse shipborne measurements area, the accuracy of ST_TCFWO is better than that of V19.1, ETOPO1 and ST_KR. Moreover, compared to other models, ST_TCFWO performs better in ﬂat submarine plain or rugged seamount area.


Introduction
The oceans account for about 71% of the world's total area and are rich in natural resources. Comprehensive seafloor topographic information, on one hand, can provide support for proper development of resources. On the other hand, the topographic information combined with gravity anomaly information ensures navigation safety [1]. Generally, shipborne sounding technology is one of the methods to obtain seafloor topography, which has high accuracy but is time-consuming for large areas. Fortunately, the satellite altimetry can quickly derive global marine gravity information [2], and the topography can be constructed based on the correlation between seafloor topography and gravity information [3]. The method of inversion of seafloor topography by satellite altimetry makes up for the deficiency of traditional sounding technology. Therefore, the exploration of using gravity information to invert seafloor topography has never stopped. With the development of space technology, there are more and more altimetry satellites [4][5][6][7], such as, Geo-sat, ERS-1/2, Jason-1/2/3, Envisat, Cryosat-1/2, HY-2, etc. The accuracy of marine gravity constructed by satellite altimetry has reached 2-3 mGal in most areas [8][9][10][11][12]. Consequently, it promotes the research of seafloor topography inversion [13][14][15][16][17]. Many methods have been proposed for inversion of seabed topography using gravity information, for example, linear response function technique [18], gravity geologic method (GGM) [19], least square configuration method [20] and linear regression analysis [21]. One of them, the GGM was originally proposed to predict the depth of bedrock by Ibrahim and Hinze [22], and then applied to the inversion of seafloor topography. According to the control points of shipborne sounding measurements, the free air gravity anomaly is divided into regional and residual gravity, and then the terrain is inverted based on the correlation between residual gravity anomaly and seafloor topography. The method is carried out in the space domain, which is simple and has high inversion accuracy. Therefore, it has been applied and improved in many references. In 2011, Kim, J.W. et al. [23] proposed the tuning density difference GGM which solved the problem of determining the density contrast. In 2017, Xiang et al. [24] proposed an adaptive triangulated finite element approximation method for non-uniform control points which solved the long wavelength accuracy difference is too large because of the uneven distribution of depth control points. In 2018, Kim and Yun [25] verified the effectiveness of the GGM in the shallow water area. In 2020, Xing [26] used a three-dimensional rectangular model to replace horizontal thin plate model to construct short wavelength gravity anomaly, and integrated prior terrain information by regularization method to improve the accuracy of seafloor topography inversion.
According to GGM principle, grid processing of discrete regional gravity anomaly data is crucial for seafloor topography inversion. However, there are few studies on how to improve the grid accuracy of regional gravity field. Moreover, there are many grid methods [27][28][29][30], so it is necessary to select them. Among them, kriging algorithm has advantages in space prediction and uncertainty analysis. It is widely used in geophysics, geology, aerospace, meteorology, image processing and other fields [31][32][33][34][35]. Hence, we try to improve this method to optimize the regional gravity anomaly model.
Due to the limitation of shipborne sounding technology, the distribution of depth control points is relatively sparse especially in the Southern Ocean [14], which has a great impact on the grid construction of regional gravity anomaly model. Compared with other interpolation methods, the kriging algorithm is based on spatial autocorrelation analysis and the spatial variation analysis of data to obtain the unbiased optimal estimate [36], it has advantages for sparse data interpolation. Moreover, the GGM assumes that the residual gravity anomaly has an approximate linear correlation with the seafloor topography, the non-linear influence is ignored [37]. This leads to the existence of some topographic information in the regional gravity anomalies. Hence, we proposed the topography constraint factor weight optimization (TCFWO) method by introducing a topographic factor to optimize the weight of kriging algorithm for regional gravity data interpolation. In the TCFWO method, the horizontal and vertical variogram models of regional gravity data were constructed, and the anisotropic effect of the horizontal variogram was considered. Subsequently, we used the TCFWO method to carry out relevant experiments on Marcus-Wake seamount area to verify the effectiveness and practicability of the method.

Construction of TCFWO Method
The free air gravity anomaly, ∆g f (x, y), obtained by satellite altimetry, which is assumed to divided into two parts regional and residual anomaly according to the shipborne depth measurements in GGM [38]: where x and y represent the latitude and longitude, respectively. ∆g reg (x, y) and ∆g res (x, y) denote as the regional and residual gravity anomalies. In the area of relatively flat terrain change, the relationship between topography and the residual gravity can be assumed by a simple Bouguer plate formula [37]: where E(x, y) is the depth value at point (x, y). G represents the universal gravitational constant (=6.672 × 10 −11 N·m 2 /kg 2 ). D is the deepest depth of the control points within the study area. Both the value of E(x, y) and D are negative from the sea level down. ∆ρ is the difference between bedrock density and seawater density. However, in GGM, ∆ρ needs to be optimized as a parameter to adjust the relationship between residual gravity and topography, its analytical significance exceeds the actual physical meaning [19]. The common methods to determine the density contrast include downward continuation method [38] and iterative method [23], in this paper, the latter method was used. By rearranging the Equation (2), the residual gravity anomalies at the control points, (x i , y i ), (i = 1, 2, . . . , n), can be calculated: where ∆g CP res (x i , y i ) are residual gravity anomalies. E CP (x i , y i ) are the depth of control points. Therefore, the corresponding regional gravity anomalies ∆g CP reg (x i , y i ) can be calculated.
where ∆g CP f (x i , y i ) are the free air gravity anomalies at control points, which can be obtained by bilinear interpolation from the satellite altimetry derived gravity model. Once the residual gravity anomaly model is acquired, the seafloor topography model can be estimated according to Equation (2). However, the regional gravity anomaly model should be constructed first.
Evidentially, the regional gravity anomalies are discrete distribution, hence, it is necessary to grid it. According to the TCFWO method, the construction process of regional gravity field is as follows.
The ∆g CP reg (x i , y i ) at the control points can be calculated according to Equation (4). Then, the unknown value of regional gravity anomaly, ∆g PP reg (x 0 , y 0 ), at the predicted position can be estimated [30] where n is the number of control points. (x 0 , y 0 ) is the predicted position coordinate of unknown attribute value. λ i (i = 1, 2, . . . , n) are weight values assigned to the known points, which determine the results of gridding. Hence, the step to determine weight values is crucial. The estimation of TCFWO method is unbiased and optimal, just like kriging algorithm. Therefore, the estimated value ∆g PP reg (x 0 , y 0 ) satisfies the following conditions.
where ∆g PP * reg (x 0 , y 0 ) is the real value at (x 0 , y 0 ) point. E[·] is the mathematical expectation symbol. Var[·] represents the variance solution symbol. In addition, based on the secondorder stationary hypothesis, the variable data satisfies the following two conditions [39]: Cov ∆g CP reg (x i , y i ), ∆g CP reg x j , y j = E ∆g CP reg (x i , y i )·∆g CP reg x j , y j − K 2 = f (h) (8) where K is any constant. Cov[·] is the covariance calculation symbol of the variable. ∆g CP reg (x i , y i ) and ∆g CP reg x j , y j are the attribute value of two arbitrary points. f (h) is a function only related to the distance h between two variables.
According to the unbiased and error optimal characteristics of TCFWO method for estimation value, a new objective function F(λ, φ) is constructed by combining Equations (5)-(8) and using Lagrange multiplier method. Hence, the problem of determining weight values, λ i , is transformed into the calculation of equation optimization with constraints of where φ is the Lagrange coefficient. γ i0 (d) is the variogram value between the control point (x i , y i ) and the predicted point (x 0 , y 0 ). γ ij (d) is the variogram value between (x i , y i ) and x j , y j . The value of variogram is a statistic describing the spatial correlation of data, which is defined as the variance of the difference between two points. Hence, the γ i0 (d) and γ ij (d) can be calculated by the experimental variogram Equation (10) [30]: where ∆g CP reg (m) and ∆g CP reg (m + d) are the regional gravity anomaly at position m and m + d. d is the Euclidean distance between two points namely the lag distance.
Then, the partial derivative of F(λ, φ) with respect to λ i and φ are calculated, and the results after derivation are zero. Finally, the equations for determining the weight value are as follows [31]: From the Equation (11), the values of variogram are important for weight determination. However, the ordinary kriging algorithm only considers the influence of twodimensional factors on variogram values. In order to calculate and optimize the weight values λ i , this paper introduced the topographic factors, z i (i = 1, 2 . . . , n), as the third variable factor based on longitude and latitude coordinates. Thus, the variogram models of gravity anomaly in horizontal and vertical directions could be constructed, respectively, and added by equivalence weight: where γ L ij and γ Z ij are the experimental variograms in horizontal and vertical directions, respectively. γ L ij , Z ij is the optimized experimental variogram. The L ij and Z ij are the lag distance in horizontal direction and vertical direction, respectively [40]: where (x i , y i , z i ) and x j , y j , z j denote the coordinates of two shipborne measurements control points with known attributes, respectively. The values calculated by Equation (12) are experimental variogram, which are discrete points. Generally, it is necessary to use the least square algorithm to fit the variogram models based on the theoretical variogram function [41] and determine the better models according to the fitting results. Then, using the determined theoretical function model, the variogram values between the corresponding points are calculated. Finally, the formula of TCFWO method for gridding regional gravity anomaly is expressed as follows: whereλ i are the optimized weight values. γ ij L ij , Z ij are the best fitting variograms of γ L ij , Z ij . Other symbols have the same meaning as above. Therefore, based on the above analysis, the process of gridding regional gravity anomaly model and inversion of seafloor topography using TCFWO method is as follows.
Firstly, satellite altimetry derived gravity anomaly at the control points, ∆g CP f (x, y), are divided into two parts residual information, ∆g CP res (x, y), and regional information ∆g CP reg (x, y). There the ∆g CP res (x, y) is determined based on the Equation (3). Secondly, the experimental variogram models of regional gravity anomaly data, ∆g CP reg (x, y), in horizontal and vertical directions are calculated, respectively, according to Equation (12). Based on the theoretical variogram model, the least square algorithm is used to fit the experiment variogram model.
Then, the values of best fitting theoretical variogram are substituted into the TCFWO Equation (14), to realize the grid processing of regional gravity anomaly and construct a model in the region.
Finally, the gridding of regional gravity anomaly model is subtracted from the satellite altimetry derived gravity anomaly model to obtain the residual gravity anomaly model. According to the Equation (2), the topographic model is calculated.

Overview of the Experimental Area and Data Preparation
The Marcus-Wake seamount group in the Western Pacific Ocean (156.00 • E-164.47 • E, 17.88 • N-26.26 • N) was selected as the study area. This area has experienced strong volcanic activity, with complex topography and rich natural resources [42]. The high precision seafloor topography is of great significance to the exploitation and utilization of natural resources and marine geological research in the region. The experimental data sources include the shipborne depth measurements data and the satellite altimetry gravity anomaly model, and existing seafloor topographic models.

Shipborne Depth Measurements Data
The distribution depth measurements provided by the National Geophysical Data Center (NGDC) (http://www.ngdc.noaa.gov/mgg, accessed on 2 September 2020) are shown in Figure 1a, there are 53,863 shipborne measurements. One of every five points was selected to form the check data set. There are 10,777 check points, as shown by the blue star dots. The check data set was used to determine the density contrast and evaluate the results of seafloor topography. There are 43,086 remaining points used as control data set for seafloor topography inversion, as shown in red dots. Among them, the maximum depth is −6190.00 m, the minimum depth is −1038.00 m, and the average is −4782.68 m, which suggests the region's terrain fluctuates sharply.

Satellite Altimetry Derived Gravity Anomaly Data
The satellite altimetry derived gravity anomaly model was provided by the Scripps Institution of Oceanography (SIO) (https://topex.ucsd.edu, accessed on 12 July 2020). As shown in Figure 1b, the version is V28.1 with a resolution of 1 × 1 . The maximum gravity value in the area is 254.80 mGal, the minimum is −63.00 mGal, and the average is −4.69 mGal. In addition, it also includes the seafloor topography model ETOPO1 (https://www.ngdc. noaa.gov, accessed on 2 November 2020) and V19.1 (https://topex.ucsd.edu, accessed on 12 July 2020) for accuracy comparison, which were analyzed in Section 4.2 below.

Geostatistical Analysis of Regional Gravity Anomaly Data
Using TCFWO method to grid gravity anomaly model, geostatistical analysis is necessary to master the spatial stationarity and variability of variables comprehensively. Quantile-quantile graphs and variogram models are generally better methods.

Second-Order Stationary Analysis
The TCFWO method is based on kriging algorithm, so the experimental data should satisfy the second-order stationary hypothesis. However, in practice, such data is difficult to exist, and there is no good method to judge the second-order stability. Generally, If the data approximately satisfies the normal distribution, we consider that the data is second-order stationary. Therefore, before gridding the regional gravity anomaly data, it is necessary to test the approximate normality of the data.
The normal distribution test is often analyzed by quantile-quantile diagram, and the regional gravity anomaly data distribution drawn with blue solid dots based on the idea of quantile diagram is shown in Figure 2. According to the figure, the sample data are basically distributed near the red straight line which is normal distribution line. Therefore, we considered that the regional gravity anomaly data approximately satisfies the normal distribution. Indeed, if the data does not normal distribution, it needs to be transformed, such as Box-cox and logarithmic transformation. Quantile figure of regional gravity anomaly data. The red line represents a normally distributed line. The blue solid dots represent the regional gravity anomaly data distribution.

Analysis of Spatial Variability Characteristics
It is crucial to master the variability characteristics of data spatial distribution for grid processing. In geo-statistics, the variogram model is used to analyze it, which mainly involves the calculation and fitting of the variogram model and the anisotropic analysis of data. The anisotropy of data means that its variogram is not identical in different directions, that is, the value of range and sill of variogram model are not identical. Instead, if the variogram is the same in different directions, it satisfies isotropy.

Construction of variogram model in the horizontal direction
According to the distribution of the regional gravity anomaly data in the horizontal direction, the Equation (13) was used to calculate the distance of points in the horizontal direction as the lag distance. In the calculation of variogram, the east-west direction is represented by 0 and 180 degrees. Therefore, in order to comprehensively analyze the variation properties in different directions, the variogram values in the four directions of 0 • , 45 • , 90 • and 135 • were calculated by Equation (12), respectively. The experimental variograms are shown in Figure 3a. From the curve, the range and sill value are different in different directions. Hence, the regional gravity anomaly data is anisotropic in the horizontal direction.
In this paper, the experimental variogram of horizontal direction was established by using the mean values of variogram models in the four directions. The distribution of experimental samples is shown as the black star in Figure 3b. Based on the theoretical variogram models, such as exponential function model, spherical function model and gaussian function model, the least square algorithm was used to fit the experimental variogram curves. The results of the three theoretical variograms fitting are shown in the red triangle line, blue square line and black solid line in Figure 3b, respectively. The statistics of the formula parameters and fitting degree are shown in Table 1.  Note. R 2 namely R-squared represents the degree of fitting, the greater the value, the higher the degree of fitting and the value range is 0~1.
From the Figure 3b and Table 1, the results of the spherical model and the gaussian model are basically the same, and the R 2 is 0.45 and 0.44, respectively. While the result of the exponential model is 0.59, which is better than the other two function models. Therefore, the exponential model was used to calculate the value of the variogram.

Construction of variogram model in vertical direction
Firstly, variogram values and the lag distance in vertical direction were calculated based on Equations (12) and (13), respectively. Then, according to the horizontal variogram fitting method, the variogram fitting curve in the vertical direction is shown in Figure 4 and the fitting results are shown in Table 2.  According to Figure 4 and Table 2, the fitting results of the three theoretical variation function models are similar. By comparison the value of R 2 , the values of spherical model and Gaussian model are 0.94 and 0.93, respectively. However, the value of the exponential model is 0.97, which shows that the fitting accuracy is better. Therefore, the exponential model was used in the construction of the variogram model in vertical direction.

Accuracy Assessment of Regional Gravity Anomaly Model
Combined with Equations (2)-(4) the discrete regional gravity anomalies of shipborne depth control points position in Section 3.1 can be calculated. In order to verify the applicability of TCFWO method, results of this method and kriging algorithm in gridding regional gravity anomalies were compared. Two methods were used to build GR_TCFWO ( Figure 5a) and GR_KR (Figure 5b) models with 1 × 1 resolution, respectively, and the leave-one-out cross-validation method was used to evaluate model's accuracy [43]. The process of cross-validation method is to first remove one point in the regional gravity anomaly data set, and then use the remaining points to predict the gravity value at that location. Then add the excluded point back to the dataset and remove another point. Do this 1000 times for samples in the dataset and compare the differences between the estimated and true values of the selected samples. Taking the estimated values of the model as the abscissa and the true values as the ordinate, the constructed scatter diagram and its regression fitting curve are shown in Figure 6. The statistical results of the difference between the two values are shown in Table 3b.   From the Figure 5, the two models are basically consistent in the whole, but there are obvious differences in some areas with sparse control points. According to the statistical results in Table 3a, the values range of GR_TCFWO model is larger than GR_KR model, but the standard deviation of GR_TCFWO model is slightly lower than GR_KR model. Therefore, the TCFWO method increased the value range of the regional gravity anomaly model, which may have revealed more details, and the stability of the model was slightly improved. From the Figure 6a,b, the results showed that the GR_KR values were typically more scatter compared to GR_TCFWO values, which indicated a lower correlation with sample values. Therefore, TCFWO method could construct a higher accuracy model than kriging method.
From the Table 3b showed that the standard deviation (STD) of cross-validation results of GR_TCFWO model was 3.62 mGal, the R 2 of liner regression is 0.98. By comparison, the STD of GR_KR model is 6.08 mGal, the R 2 is 0.92. The results showed that the accuracy of model constructed by the TCFWO method was about 40% higher than that of the ordinary kriging algorithm. Hence, the TCFWO method was more effective in gridding. Subsequently, in order to further prove the advantages of the TCFWO method, two gravity anomaly models were used to invert seafloor topography.

Inversion and Comparison of Seafloor Topography
The regional gravity anomaly models established by the above two methods provide basic data for seafloor topography inversion. According to the principle of GGM, the quality of the regional gravity anomaly model is important to the accuracy of seafloor topography inversion. Therefore, it is necessary to analyze the inversion results.

Determination of Density Contrast
Density contrast, ∆ρ, is the difference between bedrock density and seawater density, which is important to adjust the relationship between gravity and seafloor topography. Generally, the density of the global ocean bedrock is 2.70 g/cm 3 (1 g/cm 3 = 10 3 kg/m 3 ) which is the mean value between 2.67 g/cm 3 and 2.73 g/cm 3 [19], and the average density of the seawater is 1.03 g/cm 3 , so the density constant is 1.67 g/cm 3 . However, it is not suitable to use the average value for the inversion of local topography, the density constant should be accurately determined. The commonly used methods are the iterative method and the downward continuation method, among which the iterative method is widely used, since it is simple and convenient when there are more control points.
Based on the iterative density constant determination method proposed by Kim, J.W. [23], 43,086 control points and 10,777 check points were selected, and their distribution was as show Figure 1a. Then, the control points were used to estimate the seafloor topography model under different density constants, and bilinear interpolation was used to obtain the predicted topographic values at the location of the check points. By comparing the correlation coefficient and the difference standard deviation between the predicted values and the actual values, we constructed the curve Figure 7. According to Figure 7a,b the correlation coefficient increases and the STD with ∆ρ increases in a certain range. When the density contrast is greater than a certain constant, the correlation coefficient decreases and the STD increases. Through the above-mentioned calculation and fitting of the two curves, the density constant was about 1.11 g/cm 3 , the correlation coefficient curve reached the maximum value and STD was the minimum value. Therefore, in this paper the optimal density contrast 1.11 g/cm 3 was used to estimate seafloor topography.

Accuracy Comparison and Analysis of Seafloor Topography Models
In order to further illuminate the advantages of the TCFWO method compared with kriging method, the following focuses on the analysis of the inverted seafloor topography models. Firstly, the satellite altimetry gravity anomaly model V28.1 was used to subtract GR_TCFWO and GR_KR models constructed in Section 3.3, respectively, and the residual gravity anomaly model was obtained. Then, the corresponding seafloor topography model ST_TCFWO (Figure 8a) and ST_KR (Figure 8b) were obtained by the Equation (2). Finally, the accuracy of inversion models was compared and analyzed by using the shipborne measurements data and existing models such as the ETOPO1 (Figure 8c) and V19.1 (Figure 8d). From the Figure 8, the geomorphic features of the Marcus-Wake seamount group area shown by each model are basically the same in general. However, there are some differences in the local values, which can be seen from the statistical results of the average and STD in Table 4. Obviously, the STD of difference between ST_TCFWO model and ST_KR model was 42.50 m, which was the smallest. In comparison, the STD of difference between ST_TCFWO model and V19.1 model is 213.85 m and ETOPO1 model is 208.83 m, respectively. However, the STD of the difference between the ST_KR model and the general models are 221.75 m and 210.17 m, respectively. Therefore, the ST_TCFWO model is closer to the international general model.
In addition, the two international models V19.1 and ETOPO1 are also different. The mean value of difference between them is −9.4 m, and the STD is 202.32 m. The reasons may be related to their construction methods. ETOPO1 model was mainly based on the seafloor topography of V8.2 inverted by satellite altimetry, and then integrated the shipborne depth data for interpolation and smoothing [44]. However, V19.1 was predicted by satellite altimetry and sparse shipboard bathymetry. Compared with V8.2, V19.1 model has higher accuracy and resolution. Through the above comparison, there are some differences among the four models. Moreover, both the international general seafloor topographic models and the models predicted by gravity in this paper can't completely guaranteed to be correct, especially in the area where there are no shipborne depth control points. Therefore, this paper analyzed the accuracy of the seafloor topographic model by using the external check points, namely four models were interpolated to 10,777 shipborne depth measurements inspection points ( Figure 1) using bilinear interpolation, and the difference between model values and checking values were compared, the statistical results of are shown in Table 5. According to the Table 5, by comparing the STD, the model accuracy in the study area from high to low is V19.1 model, ST_TCFWO model, ST_KR model and ETOPO1 model. The STD values are 117.86 m, 129.14 m, 173.51 m and 191.71m, respectively. The correlation coefficient between the model and the shipborne sounding measurements inspection data, the accuracy can also be reflected. In contrast, the accuracy of V19.1 model is better than ST_TCFWO model. Analyzing the reasons, on one hand, the inversion band of V19.1 model was controlled in the range of 15 km to 160 km, and the shipborne measured data were used in the long wavelength part of topography [14]. On the other hand, it needs to be further explored whether the shipborne check points for model accuracy evaluation were also applied to the V19.1 model construction. However, it is obvious that the accuracy of ST_TCFWO model is better than ETOPO1 and ST_KR model. Compared with the ST_KR model (STD = 173.51 m), the accuracy of the ST_TCFWO model (STD = 129.14 m) was improved about 26%. Therefore, it is reflected that the TCFWO method has obvious advantages compared with the ordinary kriging algorithm in application of inversion seafloor topography. To compare the difference of each model in different depth, the variation of residual values between model values and checking values with depth was drawn, as shown in Figure 9. The results show that with the increase of depth, the residual values of the four models tend to zero value. Comparing the Figure 9a with Figure 9b,c, the results of ST_TCFWO model are more aggregated, which indicates the inversion advantage of TCFWO method. In contrast, the data distribution of V19.1 model is closest to the value of zero, which indicates that the model values are more consistent with the checking values overall. It is also consistent with the conclusion of the above analysis. To further analyze the accuracy of the model at different depths, the depth variation of the check points was taken as the standard. Taking a depth layer every 1000 m, and the above residual values of models in each different depth layers was shown in Table 6. With the increase of depth, the checking accuracy of ST_KR model improved. The verification accuracy of ST_TCFWO model, ETOPO1 model and V19.1 model first decreases and then increases with the increase of depth, and all reached the maximum STD in the depth range of −2000 m to −3000 m. The maximum STD of each model was 235.16 m, 420.47 m and 249.57 m, respectively. Therefore, the stability of ST_TCFWO model is basically consistent with international models. Furthermore, in the range of −1000~−3000 m and more than −6000 m, ST_TCFWO model is better than the other three models. However, in the range of −3000 m to −6000 m (which is also the range where most of the check points gather), the accuracy of models from high to low was V19.1, ST_TCFWO, ST_KR and ETOPO1. This is a very interesting result. In this depth range, the accuracy of V19.1 model is better. The specific reasons need further discussion in future. However, it is obvious that ST_TCFWO model has better accuracy than ST_KR model in all ranges, indicating an advantage of TCFWO method.
Then, as shown in Figure Table 7. We defined relative precision as the ratio of the difference values to the shipborne depth checking values. RESTD in the table is the standard deviation of relative precision. In this paper, the mean value of RESTD of each model was taken as the criterion to judge the large difference. If the relative accuracy is greater than the mean of RESTD value (8.26%), it is considered as a large difference point. The distribution of large difference points was shown in the red dots in Figure 10. We can see that the red dots are basically distributed near the ridge and seamount where the topographic changes dramatically. Among them, ETOPO1 has the largest number of red dots and V19.1 has the smallest number. The distribution of red points in ST_TCFWO model is like that in ST_KR. However, in some regions, the number of dots in ST_TCFWO is significantly less than ST_KR. From the distribution of red dots in the map, the accuracy of the model is affected by the drastic change of terrain. We think that the reasons may be as follows: (1) The accuracy of gravity anomaly obtained by satellite altimetry is reduced due to the sharp relief of topography. (2) Most of the seamount in the study area were formed during the mid-Cretaceous, resulting in a large amount of sediment accumulation. Therefore, the correlation between gravity anomaly and seafloor topography is reduced.  Note. min. = minimum value; max. = maximum value; STD = standard deviation. RESTD= the standard deviation of relative precision. Rate is the standard deviation change of comparison with the overall statistical results in Table 5 which represents the stability of the model.
From the Table 7, the extreme, mean values and STD of ST_TCFWO local model are smaller than that of three other models, with the highest accuracy among the four models. Comparing with the local of ST_KR model (STD = 207.59 m), the accuracy of ST_TCFWO local model was improved by about 26.35%. Moreover, comparing Tables 5 and 7 above, due to the reduction of control points, the accuracy of the four models is significantly reduced and the accuracy change rate is different. The STD of local of ST_TCFWO model is 152.88 m, compared with ST_TCFWO model (STD = 129.14 m), the accuracy change rate is 18.38%. The STD change rate of ST_KR is 19.64%. More obvious is the V19.1 model, with a change rate of 35.47%. Obviously, the ST_TCFWO model is better than ST_KR in both accuracy and stability, and is less affected by the distribution of control points.
Subsequently, the accuracy of four local seafloor topographic models in different areas was compared. We delineate two areas A and B with different geomorphic types for analysis, as shown in Figure 10. The control points in area A and B are rare and distributed discretely. In addition, the terrain of area A is relatively flat, mostly abyssal plain landform. In area B, there are seamounts, ridge and other landforms, and the terrain changes dramatically. The shipborne measurements check points in the two areas are used to evaluate the accuracy of models. The results of difference between the model values and checking values were shown in Table 8. From the Table 8, in the area A, the accuracy of V19.1 model was the best. The second was ST_TCFWO model, its STD was 46.62 m. The accuracy ST_KR model was the lowest, and the STD was 66.36 m. In the area B, the accuracy of ST_TCFWO model was better than the other three models. However, the accuracy of the four models was not ideal. This result was consistent with the distribution of red relative precision points in Figure 10. It is worth affirming that ST_TCFWO model has better accuracy than ST_KR model no matter in flat submarine plain or rugged seamount area, which shows the superiority of TCFWO method.

Conclusions
When GGM is used to invert seafloor topography, due to the distribution of control points, it is necessary to grid the regional gravity anomaly model by using interpolation method, thereby, to improve the accuracy of seafloor topography. A study area of Marcus-Wake seamount group was selected to evaluate the reliability of the TCFWO method. Then, TCFWO and kriging methods were used to grid the regional gravity anomaly, and then the corresponding seafloor topography models were inversed. Finally, the inversion models were compared with shipborne depth measurements data and existing models, including ETOPO1 and V19.1, respectively. The following useful conclusions were drawn.

1.
The cross-validation results show that the STD of the GR_TCFWO regional gravity anomaly model gridded by TCFWO method was 3.63 mGal, which was more accurate than the GR_KR model constructed by the ordinary kriging algorithm nearly 40%. Hence, TCFWO method has optimized the interpolation weight of kriging algorithm by introducing topography constraint factor, and has obtain a better grid model. The GR_TCFWO model is a better data for seafloor topography inversion.

2.
The accuracy of the models was evaluated by using the shipborne sounding data. The STD of ST_TCFWO model was 129.14 m, and its accuracy was slightly lower than that of V19.1 model and better than the ST_KR and ETOPO1 models. Comparing with the ST_KR, and the accuracy of ST_TCFWO model was improved by about 26%. The accuracy of the models was affected by the change of depth. In each depth layer, the accuracy of ST_TCFWO model was better than that of ST_KR model. Moreover, ST_TCFWO model showed better advantages than other models in the depth layer of −1000 m~−3000 m and more than −6000 m.

3.
The accuracy of the model was affected by the distribution of control points, but the influence on ST_TCFWO model was relatively small. In the area with sparse shipborne measured control points, the accuracy of ST_TCFWO model was better than ST_KR, ETOPO1 and V19.1 model. By analyzing the model in different topographic features, the accuracy of the model was greatly affected by the terrain changes and the points with low relative accuracy were mainly distributed in the locations where the terrain changes violently. However, ST_TCFWO model has better accuracy than ST_KR model no matter in flat submarine plain or rugged seamount and ridge area, which shows the superiority of TCFWO method.