An Improved Coastal Marine Gravity Field Based on the Mean Sea Surface Height Constraint Factor Method

Abstract

Construction of a high spatial resolution and high precision marine gravity field in coastal areas is constrained by the low quality and sparse coverage of altimetry data, except for limited shipborne and airborne gravity surveys. To address this problem, a mean sea surface height constraint factor (MSSHCF) method based on the ordinary kriging method and the remove-restore technique is proposed from the perspective of interpolation. In this method, the data is standardized during the interpolation process to reduce the error and mean sea surface as variables related to the marine gravity field are added to the semi-variance function in ordinary kriging to obtain a marine gravity field with a spatial resolution of 1′ × 1′. Validation experiments show that the MSSHCF method more closely agrees with the referenced SS V28, DTU17 global marine gravity models than the ordinary kriging method. Our results were further validated against shipborne data; the accuracy of the MSSHCF method is 0.13 and 0.33 mGal higher than that of the ordinary kriging method in two experimental areas. The effects of ocean depth and offshore distance on the results were also assessed. These results show that the proposed method is more accurate than the ordinary kriging method, when the distance and depth varied. Therefore, our study demonstrates that the MSSHCF method is an innovative and feasible tool for extracting gravity fields along coastal, beach, and island areas.

1. Introduction

Accurate knowledge of the coastal gravity field is important when studying coastal ecosystem processes, as well as facilitating other off shore activities [1]. Gravity in coastal areas can be observed close to the coastline by shipborne gravity measurement, but this is costly and time-consuming [2]. Airborne methods can provide uniform measurements, but this approach might not be applicable due to the high cost [3,4,5,6]. Currently, there are two ways to improve the marine gravity field, including increasing the number of measurements and improving their quality [7]. A detailed marine gravity map of the Southern Ocean and Antarctic margin was obtained by adding high quality Geosat/ERM satellite data [8]. Meanwhile, the overall accuracy of the regional marine gravity field was improved to 5.5 mGal by using high density Geosat/GM altimeter satellite data [9]. The fusion of different satellite altimeters, especially GM altimeter satellite data types can further improve the accuracy of the marine gravity field to some extent [10,11,12,13,14,15]. Versions of the global marine gravity field calculated by fusing Cryosat-2, Envisat, and Jason-1 data achieved accuracy greater than 2 mGal in local areas [16]. In addition to the data fusion strategy, waveform retracking of satellite altimetry data, especially in the immediate vicinity of the coasts, can also increase the number and quality of measurements. At present, satellite altimeter technology leads to sparse data and low quality in coastal areas [17]. However, the use of an improved threshold re-tracker for Geosat/GM data improved the accuracy of the marine gravity field by approximately 11% within 10 km of Taiwan island [18]. Furthermore, a “Double Retracking” technique was used to process ERS-1 geodetic mission waveform data and improved the accuracy of the marine gravity field by 20% overall and by 40–50% in coastal areas [19]. In another study, Cryosat-2, Envisat, and Jason-1 satellite altimeter waveform maps were retracked for a significant improvement in the accuracy of the marine gravity field [13]. Other work has also shown that the use of the proposed fitting (also called retracking) strategy (ALES+) can improve the determination of mean sea surface at high latitudes [20], and the combination of different waveform retracking methods can also improve the accuracy of the marine gravity field to some extent, especially in coastal areas [21,22,23,24].

The interpolation method is needed along the coastal area because of a lack of real measurements. For example, a marine gravity field that meets the requirements can be generated using the Shepard interpolation method [25]. Moreover, the kriging interpolation technique is widely applied to produce geophysical maps, including gravity, a magnetic field, and sea depth topographic maps [26,27,28,29]. However, even for the same spatial resolution, varying interpolation methods make a significant difference in the accuracy of gridded marine gravity fields from the same data source [30]. Furthermore, when commonly used interpolation methods were compared, including inverse distance weighted (IDW), radial basis function (RBF), Shepard, and ordinary kriging, kriging has the highest accuracy but with an unacceptable accuracy loss in areas with sparse data density [31,32]. From the perspective of statistical analysis, adding data correlated with the predicted variables can improve the weight of the kriging predicted points [33]. Therefore, in order to improve accuracy of the gridded results, some researchers began to add variables to optimize the weights for the prediction points when executing kriging interpolation [34,35]. Using the Earth’s ellipsoid as a reference, the mean sea surface includes the geoid and the dynamic seafloor topography, commonly approximated from the marine geoid, a factor related to earth gravity estimated by geodesists and geophysicists [36]. Likewise, the mean height of the sea surface can play a role in reconstructing the marine gravity field [33]. In view of the sparse data density along coastal areas, this paper will present a way to improve the derived marine gravity field along the coast, using a mean sea surface height constraint factor (MSSHCF).

Unlike previous research, this study investigates the impact of interpolation methods on construction of high spatial resolution marine gravity field in coastal regions. The goal of our study is to develop an improved interpolation method that enhances the quality of the coastal marine gravity field. To address the quality of coastal gravity, we propose a mean sea surface height constraint factor (MSSHCF) method based on the ordinary kriging interpolation and a remove-restore technique. The proposed method introduces mean sea surface height factor as the distance in the vertical direction to improve the accuracy of the results at unknown points by optimizing the weight information of the control points and then improves the accuracy of the coastal marine gravity field. In this study, we take two local coastal areas as an example in the South China Sea (SCS) to validate the approach.

2. Materials and Methods

We proposed a mean sea surface height constraint factor (MSSHCF) method based on the ordinary kriging interpolation method. Ordinary kriging provides a more comprehensive picture of the global and local characteristics of marine gravity field fluctuations than other interpolation methods. In geo-statistics, the ordinary kriging method is a widely used estimation method, and it is usually stable in terms of spatial distribution data when gridding discrete points [37]. Furthermore, the ordinary kriging interpolation method can predict values for un-measured points using the weight information generated by the semi-variance function from the measured points. The introduction of correlation variables can improve the weighting information and then further improve the accuracy of the points to be measured. Thus, we proposed an MSSHCF method based on the ordinary kriging interpolation method to improve the accuracy of the coastal points to be measured.

2.1. The Mean Sea Surface Height Constraints Factor (MSSHCF) Method

The methods for deriving a marine gravity field using satellite altimeter data can be broadly classified into two types: least squares configurations (LSC) and deflection of vertical (DOV). Unlike the LSC method, the DOV method can suppress long wavelength radial orbital errors below the noise level of the altimeter, using a simple gradient along the orbital geodesic [28]. We derived the marine gravity field by the Inverse Vening–Meinesz (IVM) formula based on north and east components of the deflection of vertical [38]. In order to reduce the introduction of errors while ensuring an adequate amount of gridded marine gravity data, the marine gravity field with a spatial resolution of 2-arc minute was derived for the selected region. In contrast to getting 1- and 4-arc minute marine gravity directly, this approach ensures a balance between the accuracy of the starting gravity field and the amount of computational data required for interpolation.

A high spatial resolution marine gravity field is constructed based on the principles of ordinary kriging interpolation method [29]:

$$G_p(x_0,y_0)={\displaystyle \sum _{i=1}^n{a}_i(x_0,y_0)}\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_i,y_i)+{\mathit{g}}_{\mathit{\text{EGM}} 2008}\mathrm{}(x_0,y_0)$$

(1)

where $\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_i,y_i)$ is the final residual marine gravity anomaly obtained using DOV. The term ${\mathit{g}}_{\mathit{\text{EGM}} 2008}$ denotes the reference gravity anomaly based on the EGM2008 model. The value $a_i(x_0,y_0)$ is the weight factor of known control points in the grid. The predicted values $G_p(x_0,y_0)$ in the ordinary kriging method also satisfy isotropic conditions [27]:

$$\{\begin{array}{l}\mathrm{E}[\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_0,y_0)-\mathsf{\Delta }{\tilde{g}}_{\mathit{\text{res}}}(x_0,y_0)]=0\\ \mathrm{Var}[\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_0,y_0)-\mathsf{\Delta }{\tilde{g}}_{\mathit{\text{res}}}(x_0,y_0)]=\min \end{array}$$

(2)

where Var[.] is the variance sign and Equation (2) to ensure the minimum; E[.] is the sign of the mathematical expectation of the variable. The term $\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_0,y_0)$ refers to the inverse gravity field residuals and $\mathsf{\Delta }{\tilde{g}}_{\mathit{\text{res}}}(x_0,y_0)$ refers to the true point gravity field residuals at the point $(x_0,y_0)$. The semi-covariance function is calculated as follows:

$$\gamma (d)=\frac{1}{2}\mathrm{E}[{(\mathsf{\Delta }{\tilde{g}}_{\mathit{\text{res}}}(r)-\mathsf{\Delta }{\tilde{g}}_{\mathit{\text{res}}}(r+d))}^2]$$

(3)

where d is the Euclidean distance between the point $(x_0,y_0)$ to be predicted and the current point $(x_i,y_i)$. The r is the position of the starting point, and $\gamma$ is the semi-covariance function describing the spatial correlation of the random function. Equation (2) can be obtained based on the expected relationship between the predicted and known points:

$$\begin{array}{l}\mathrm{E}[\mathsf{\Delta }{\tilde{g}}_{\mathit{\text{res}}}(x_0,y_0)-\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_0,y_0)]=\mathrm{E}[{\displaystyle \sum _{i=1}^n{a}_i}\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_i,y_i)-\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_0,y_0)]\\ =\mathrm{M}({\displaystyle \sum _i^n{a}_i}-1)\end{array}$$

(4)

where the M is an arbitrary constant. Based on the relationship between $(x_0,y_0)$ and $(x_i,y_i)$, and according to Lagrange’s theorem, the weight coefficients $a_i$ (i = 1, 2, 3, …, n) can be obtained according to the following equation:

$$\{\begin{array}{l}{\displaystyle \sum _{i=1}^n{a}_i{\gamma }_{\mathit{\text{ij}}}(d)+\mu ={\gamma }_{\mathit{\text{j}}\mathrm{\text{0}}}(d)}\\ {\displaystyle \sum _{i=1}^n{a}_i=1}\end{array}$$

(5)

With Equation (5), the weight coefficients in the ordinary kriging method can be determined. When the mean sea surface is introduced as a constraint, the distance is redefined and calculated as follows:

$$\{\begin{array}{l}{S}_{\mathit{\text{ij}}}=\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}\\ D_{\mathit{\text{ij}}}=|h_i-h_j|,(i,j=0,1,2,3,...,n)\end{array}$$

(6)

where $S_{\mathit{\text{ij}}}$ and $D_{\mathit{\text{ij}}}$ are the distance variables and the variables in the vertical direction of the perpendicular face, respectively. The definition of the distance variable and the semi-variance function in the spatial domain redefined according to Equation (6) can be regained as follows [30]:

$$\{\begin{array}{l}{\gamma }_{\mathit{\text{ij}}}(D_{\mathit{\text{ij}}},S_{\mathit{\text{ij}}})=\frac{1}{2}({\gamma }_{\mathit{\text{ij}}}(D_{\mathit{\text{ij}}})+{\gamma }_{\mathit{\text{ij}}}(S_{\mathit{\text{ij}}}))\\ {\displaystyle \sum _{i=1}^n{\widehat{a}}_i{\gamma }_{\mathit{\text{ij}}}(D_{\mathit{\text{ij}}},S_{\mathit{\text{ij}}})+\mu ={\gamma }_{\mathit{\text{j}}\mathrm{\text{0}}}(D_{\mathit{\text{ij}}},S_{\mathit{\text{ij}}})}\\ {\displaystyle \sum _{i=1}^n{\widehat{\lambda }}_i=1}\end{array}$$

(7)

where ${\widehat{\lambda }}_i$ is the redetermined optimized weight factor. The values ${\gamma }_{\mathit{\text{ij}}}(D_{\mathit{\text{ij}}})$ and ${\gamma }_{\mathit{\text{ij}}}(S_{\mathit{\text{ij}}})$ are the semi-covariance functions obtained for the distance variables in Equation (6), respectively. The optimized weighting coefficients ${\widehat{\lambda }}_i$ can be obtained by Equation (7). The values predicted by the MSSHCF method can be obtained from the following equation:

$$G^{MSSFCW}(x_0,y_0)={\displaystyle \sum _{i=1}^n{\widehat{\lambda }}_i}\mathsf{\Delta }{\mathit{g}}_{\mathit{\text{res}}}(x_i,y_i)+{\mathit{g}}_{\mathit{\text{EGM}} 2008}(x_0,y_0) (i j=1, 2, 3, \dots , n)$$

(8)

Therefore, the inversion method of the marine gravity field is divided into two main steps: the first one is to calculate the low spatial resolution marine gravity field using the deflection of vertical method and the second step is to interpolate low spatial resolution marine gravity field using the ordinary kriging and the proposed method, i.e., to obtain a high spatial resolution marine gravity field. Before the second step is performed, the data involved in the reconstruction is simultaneously normalized to reduce errors.

2.2. Data Preparation and Pre-Processing

Two local areas of the South China Sea (2° × 2°) are arbitrarily selected as the main experimental area A (20°~22°N, 113°E~115°E) and B (21°~23°N, 116°E~118°E). The selected areas are not only adjacent to the coast but also have a large amount of shipborne data. The satellite altimeter data sourced from Jason-1/GM is used to inverse the marine gravity field. Because the data distribution of Jason-1/GM is denser, this data provides more than three times as much geodetic missions altimetric sea surface height observations than previously available [19]. The study areas and satellite distribution trajectories are shown in Figure 1.

For the satellite altimetry data in area A and area B, we processed the data to obtain many reliable data sources before deriving the marine gravity field. The specific correction models and methods are shown in

Table 1. Correction methods and models for obtaining sea surface height from Jason-1/GM.

Correction	Modification Model or Method
Reprocessed method	ALES
Dry troposphere correction	European Center for Medium Range Weather Forecasting
Wet troposphere correction	European Center for Medium Range Weather Forecasting
Ionospheric correction	CNES
Marine state bias	CNES
Derived barometer height correction	European Center for Medium Range Weather Forecasting
High frequency fluctuations of the sea surface topography	LEGOS/CLS/CNES
Marine tide	LEGOS/CNES
Load tide	GOT4.10
Pole tide	Aviso

.

In Table 1, we have chosen different correction models for different errors, including the path and environmental impact errors introduced in the process of sea surface height acquisition. The ALES waveform retracking method was used to process the altimetry data to obtain more reliable data. (https://openadb.dgfi.tum.de/en/data_access/ (accessed on 27 May 2022)).

2.3. Dynamic Sea Surface Topography, Mean Sea Surface and Reference Gravity Field Model

The dynamic sea surface topography we selected was the DTU15 MDT model; it combines the gravity model EIGEN-6C4 with DTU15MSS mean sea surface model from Technical University of Denmark (DTU). Furthermore, the DTU18MSS of 1-arc minute for mean sea surface has been developed from various data sources, which contained 3 years Sentinel-3 A and the improved 7 years Cryosat-2 LRM records with a data-corrected model, which is currently the more commonly used model [39] (i.e., using FES2014 as the ocean tide model) (https://www.space.dtu.dk/ (accessed on 27 May 2022)).

The reference gravity field model plays a key role in the remove-restore technique. In our research, we used the EGM2008 global gravity field model released by the National Geospatial Intelligence Agency of the United States as a reference. The spherical harmonic expansion order of this model reached 2159, and the extended parameters extended to 2190 are also provided, which resolved the Earth’s gravity field [40] (http://icgem.gfz-potsdam.de (accessed on 27 May 2022)).

2.4. The Gravity Data Source for Comparative Verification

The gravity data sources for comparative verification come from two forms: global marine gravity field models and shipborne gravity data. The two marine gravity field models are DTU17 published by the Technical University of Denmark (DTU) and SS V28 published by the Scripps Institution of Oceanography (SIO), respectively. Both accurate gravity field models fused data from multiple satellite altimeters to provide the global marine gravity field values with a spatial resolution of 1′ × 1′. The shipborne gravity data comes from national geophysical data center (NGDC) [22].

3. Results

After obtaining a reliable source of satellite altimetry data, the 2-arc minute residual marine gravity field was obtained using the deflection of vertical method, and the 1-arc minute marine gravity field was obtained using the MSSHCF method, the ordinary kriging method and the EGM2008 reference gravity field model. A flow diagram showing how the 1-arc minute marine gravity reference map was acquired as shown in Figure 2. Figure 2 shows the two operations to obtain the 1-arc minute marine gravity field, including data processing, DOV, and the use of interpolation methods to obtain the final marine gravity field.

As shown, the first step consisted of data processing and the acquisition of an initial 2-arc marine gravity field by the same conditions. The second step ensured that the control conditions remained unchanged and for the normalized data to be interpolated using the proposed method and the ordinary kriging method. The results from geostatistical methods were tested by comparing the calculated values with the overlap points of the model and shipborne values. The root mean square error (RMSE) was used to assess the final results. RMSE is defined as:

$$RMSE=sqrt(\frac{1}{n}{\displaystyle \sum _i^n{(G^{MSSHCF}(X_i,X_i,N_i^{’})-G(X_i,X_i,N_i^{’}))}^2}),i=1,2,3.....n$$

(9)

where $G$ and $G^{MSSHCF}$ are the measured and predicted values, and n is the number of samples.

In the previous section, we detailed the specific process and methods used to obtain the 1′ × 1′ marine gravity field. As discussed, the MSSHCF method introduces the mean sea surface height as the weight optimization factor based on the ordinary kriging interpolation method, and the optimized weight information is used to predict unmeasured values. To compare the performance of the coastal gravity field from the ordinary kriging interpolation method and the MSSHCF method, we assessed the coastal marine gravity field using global marine gravity field models (SS V28, DTU17). The obtained final results for areas A and B and SS V28 and DTU17 marine gravity reference maps in the corresponding areas are shown in Figure 3 and Figure 4. Moreover, the corresponding final statistical results are shown in

Table 2. Statistical results of the comparison between the marine gravity models and the final results for areas A and B (units: mGal).

Region	Method	Model	Max	Min	Mean	Std
A	The ordinary kriging interpolation method VS	SS	28.32	−20.2	−0.83	3.25
		DTU	9.68	−7.15	−0.45	2.01
		DR	43.24	−33.74	−1.27	9.3
		EGM2008	9.25	−5.58	−0.35	1.79
		SS	28.61	−19.36	−0.66	2.78
	The MSSHCF method VS	DTU	8.34	−5.65	0.28	1.48
		DR	36.13	−33.87	−1.09	8.78
		EGM2008	7.88	−6.37	−0.18	1.29
		SS	15.53	−13.48	0.383	3.86
B	The ordinary kriging interpolation method VS	DTU	12.73	−12.43	−0.41	2.44
		DR	61.04	−41.21	0.8	9.9
		EGM2008	10.29	−12.45	0.28	1.9
		SS	14.58	−11.44	0.3	3.5
	The MSSHCF method VS	DTU	9.5	−6.41	0.33	2.04
		DR	58.74	−34.45	0.72	9.44
		EGM2008	7.79	−5.83	0.2	1.44

Note: DR denotes the direct result with a spatial resolution of 1′ × 1′ by IVM.

. In the black and red boxes of the marine gravity field in Figure 3 and Figure 4, we can find inconsistencies in the gravity distribution.

In Figure 3a, the distribution intervals for the DTU17 model are −27.39 to 20.80, and in Figure 3c the distribution intervals for SS V28 model are −46.63 to 21.61 mGal. Figure 3b shows a marine gravity reference map by ordinary kriging interpolation method, and it contains 14,331 points which range from −31.15 to 21.60 mGal. Figure 3d depicts the marine gravity reference map by the MSSHCF method, again containing the same amount of data with values ranging from −28.32 to 21.76 mGal. For region A, the gravity field obtained is closer to DTU17, with a larger difference in the minimum value compared with the SS V28, as can also be seen from the distribution of SS V28 in Figure 3b, but the trend of the gravity field remains unchanged. The distribution intervals from Figure 4a, c is −17.80 to 42.41 and −15.14 to 41.02 mGal, respectively. Figure 4b, d show marine gravity reference maps by the ordinary kriging and MSSHCF method, and it contains 14,224 points, which range from −15.03 to 40.49 and −15.08 to 40.55 mGal, respectively. For region B, in Figure 4, the inversion results do not vary much from the two models SS V28, DTU17, but there are still some different gravity distributions in offshore and complex regions. The final statistical results including DTU17, SS V28, EGM2008, and direct results (DR) by IVM are shown in Table 2.

Table 2 presents a comparison of the results for the tested gravity field models; we can draw the following inferences from this data. The statistical results of areas A and B in Table 2 show that when compared with the gravity field models EGM2008, DTU17, and SS V28, the Standard deviation (Std) of the MSSHCF method is more accurate than the ordinary kriging interpolation method. Furthermore, compared to the DR in region A, the Std for the ordinary Kriging interpolation and the MSSHCF methods were 9.30 and 8.78 mGal, respectively. Similarly, in region B, the Std is 9.90 for the ordinary Kriging interpolation method and 9.44 mGal for the MSSHCF method. Therefore, compared to the ordinary kriging interpolation and the MSSHCF methods, the accuracy of the gravity field obtained by using IVM directly is lower due to the sparse and low-quality data.

Scatter plots of the obtained gravity anomaly versus the global marine gravity models in areas A and B are shown in Figure 5 and Figure 6. We can see that the inversion results of the MSSHCF method and the ordinary kriging method have different coefficients of agreement with the global gravity field models.

As shown in the scatter plots for Figure 5 of area A and Figure 6 of area B, closer agreement is found between the MSSHCF method and the gravity field models, which verifies that the MSSHCF method outperforms the ordinary kriging method. Thus, the final results of MSSHCF method have a high coherence coefficient with the same global gravity field models in two regions. In contrast to the ordinary kriging method, the MSSHCF method settings were: a nested variogram consisting of two variogram function, one for the horizontal and another for the vertical direction. Under the same conditions, the weighting factor of the point to be measured can be optimized by increasing the amount of the relationship related to the oceanic gravity field. by introducing variables in the vertical direction as constraint on the gravity field. Thus, a more realistic map of the gravity distribution is possible to derive in these regions by the MSSHCF method.

4. Discussion

4.1. Verification from Coastal Gravity Field Models

The gravity anomaly residual maps for the corresponding locations are shown in Figure 7 and Figure 8. In the black and red boxes of the marine gravity field in Figure 7 and Figure 8, we can find some inconsistencies in the gravity residual distribution. The influence is also due precisely to the sparse amount of coastal data and the low reliability, compared to gravity field models incorporating multiple sources of data. In addition, regions with a large distribution of gravity anomaly values can be significantly different. Related research showed that there are implications for satellite altimeter inversion of the marine gravity field in some areas with significant variation in the seafloor topography [22]. In addition to the density and reliability of the satellite altimeter data, topographical factors may also lead to inconsistency in the trend of the distribution of values during the interpolation process [32]. In two figures, there are significant differences between the marine gravity fields obtained by the proposed method and the ordinary kriging method. However, excluding these data deficiencies and the topography, we find that the MSSHCF method improves the gravity field along the coast more than the ordinary kriging method.

To fully assess the performance of the two methods in coastal regions A and B, we compared the two methods with the DTU17, SS V28, and EGM2008 gravity field models and DR with varying offshore distances. The statistical results are shown in

Table 3. The statistical results of the two methods at different coastal distances, area A (units: mGal).

Distance/km	Model Comparison Results	Max	Min	Mean	Std
	The ordinary kriging	4.07	−9.02	−1.52	2.59
	The MSSHCF VS DTU17	3.72	−7.51	−1.03	1.92
0~20	The ordinary kriging	3.8	−8.06	−2.03	2.46
	The MSSHCF VS EGM2008	4.55	−7.88	−1.53	1.99
	The ordinary kriging	15.84	−33.89	−13.93	9.97
	The MSSHCF VS DR	13.98	−33.57	−13.44	9.47
	The ordinary kriging	18.9	−26.38	−0.19	6.45
	The MSSHCF VS SS V28	18.22	−26.25	0.3	5.71
	The ordinary kriging	6.3	−5.64	0.9	2.13
	The MSSHCF VS DTU17	4.09	−4.92	0.57	1.58
20~50	The ordinary kriging	4.77	−7.34	0.58	1.83
	The MSSHCF VS EGM2008	3.53	−4.96	0.31	1.19
	The ordinary kriging	23.4	−31.36	0.8	8.9
	The MSSHCF VS DR	21.51	−28.28	0.53	8.25
	The ordinary kriging	11.95	−9.29	1.39	3.96
	The MSSHCF VS SS V28	11.04	−8.57	1.07	3.47
	The ordinary kriging	7.12	−2.3	1.7	1.66
	The MSSHCF VS DTU17	4.47	−2.07	1.07	1.16
50~100	The ordinary kriging	5.57	−1.48	1.69	1.46
	The MSSHCF VS EGM2008	5.23	−0.95	1.04	1.04
	The ordinary kriging	33.72	−10.7	8.38	7.6
	The MSSHCF VS DR	33.87	−9.86	7.72	7.14
	The ordinary kriging	8.32	−4.14	1.79	2.42
	The MSSHCF VS SS V28	6.6	−4.14	1.17	1.93
	The ordinary kriging	5.12	−5.26	−0.08	1.56
	The MSSHCF VS DTU17	3.53	−3.45	−0.07	1.16
100~200	The ordinary kriging	3.0128	−4.1208	−0.1792	1.28
	The MSSHCF VS EGM2008	2.7119	−2.9765	−0.1255	0.9
	The ordinary kriging	29.24	−16.08	−0.31	6.53
	The MSSHCF VS DR	29.07	−14.3	−0.26	6.14
	The ordinary kriging	7.82	−8.29	0.17	2.37
	The MSSHCF VS SS V28	7.38	−6.68	0.18	2.01

and

Table 4. The statistical results of the two methods at different coastal distances, area B (units: mGal).

Distance/km	Model Comparison Results	Max	Min	Mean	Std
	The ordinary kriging	12.53	−5.2	0.39	4.16
	The MSSHCF VS DTU17	7.81	−5.36	0.36	3.62
0~20	The ordinary kriging	12.1	−6.03	0.32	4.57
	The MSSHCF VS EGM2008	5.17	−4.34	−0.07	2.49
	The ordinary kriging	38.81	−29.58	−1.98	17.11
	The MSSHCF VS DR	31.87	−29.31	−2.37	15.15
	The ordinary kriging	14.61	−9.35	3.07	6.02
	The MSSHCF VS SS V28	14.11	−7.4	2.95	5.33
	The ordinary kriging	5.93	−12.31	−1.78	3.9
	The MSSHCF VS DTU17	8.29	−9.29	−1.25	3
20~50	The ordinary kriging	12.4	−10.17	−1.59	3.95
	The MSSHCF VS EGM2008	5.52	−7.79	−1.25	2.5
	The ordinary kriging	41.16	−61.04	−10.16	18.31
	The MSSHCF VS DR	34.45	−58.74	−9.82	16.98
	The ordinary kriging	11.76	−15.39	−2.82	5.01
	The MSSHCF VS SS V28	6.99	−16.43	−3.01	4.67
	The ordinary kriging	5.58	−8.8	−0.29	2.19
	The MSSHCF VS DTU17	4.54	−6.86	−0.3	1.96
50~100	The ordinary kriging	3.6	−6.87	0.1	1.52
	The MSSHCF VS EGM2008	2.54	−4.93	0.08	1.24
	The ordinary kriging	21.34	−32.71	1.81	8.02
	The MSSHCF VS DR	20.83	−30.77	1.8	7.73
	The ordinary kriging	7.92	−10.22	−0.62	3.52
	The MSSHCF VS SS V28	7.06	−9.77	−0.64	3.3
	The ordinary kriging	5.12	−5.83	−0.27	1.82
	The MSSHCF VS DTU17	4.7	−5.63	−0.24	1.65
100~200	The ordinary kriging	2.61	−2.01	−0.26	1.04
	The MSSHCF VS EGM2008	2.53	−2.4	−0.23	0.87
	The ordinary kriging	13.45	−10.24	−0.8	4.99
	The MSSHCF VS DR	13.57	−10.02	−0.78	4.81
	The ordinary kriging	10.38	−7.33	0.69	2.97
	The MSSHCF VS SS V28	10.48	−6.99	0.71	2.82

.

The statistical results from area A presented in Table 3 show that the proposed method improved accuracy more than the ordinary kriging method in 0 to 200 km of distance interval. The statistical results for area B in Table 4, were similar as for area A. The results in Table 3 and Table 4 also show that in the range of distances from 0 to 200 km along the coast, the accuracy of the results when directly using IVM was low when compared to results from the interpolation methods. Among these results, in the range of 0~100 km, directly using IVM had the greatest limitation on accuracy. The absolute Std results of the MSSHCF and the ordinary method after subtracting the marine gravity models are shown in Figure 9.

Figure 9a,b denote regions A and B, respectively. From the variation in Figure 9a,b, we find that as the distance changes, the inversion results are more and more accurate when compared with the gravity field models. Furthermore, when compared with the ordinary kriging method, the proposed method improves accuracy of the results more when the offshore distance is close than when the offshore distance is far. For offshore distances greater than 100 km, although the proposed method has a lower degree of improvement compared with the ordinary kriging method, it is more accurate than the ordinary method for both A and B regions. Thus, it is also demonstrated that the MSSHCF method has greater selectivity than the ordinary interpolation method in coastal areas.

4.2. Effect of Sea Depth on Coastal Gravity Field

Different sea depth in the nearshore areas A and B were considered, and the statistical results of MSSHCF method and ordinary kriging method compared with the DTU17, SS V28, and EGM2008 gravity field models and DR at sea depths are shown in

Table 5. The statistical results of the two methods at different sea depth, area A (units: mGal).

Sea Depth	Model Comparison Results	Max	Min	Mean	Std
	The ordinary kriging VS SS V28	20.31	−28.61	1.04	3.41
	The MSSHCF VS SS V28	19.36	−28.59	0.8	2.95
0~100	The ordinary kriging	5.58	−9.25	0.43	1.9
	The MSSHCF VS EGM2008	6.37	−7.88	0.21	1.38
	The ordinary kriging	33.74	−43.24	1.33	9.92
	The MSSHCF VS DR	33.87	−36.13	1.11	9.37
	The ordinary kriging VS DTU17	7.38	−9.2	0.57	2.11
	The MSSHCF VS DTU17	5.65	−8.34	0.33	1.56
	The ordinary kriging VS SS V28	7.82	−9.26	0.36	2.48
	The MSSHCF VS SS V28	7.77	−7.65	0.31	2.09
100~200	The ordinary kriging	3.01	−4.12	0.1	1.4
	The MSSHCF VS EGM2008	3.2	−3.41	0.08	1
	The ordinary kriging	16.85	−16.12	1.29	7
	The MSSHCF VS DR	16.9	−15.23	1.26	6.56
	The ordinary kriging VS DTU17	4.37	−5.54	0.18	1.6
	The MSSHCF VS DTU17	3.65	−3.7	0.13	1.18
	The ordinary kriging VS SS V28	6.46	−5.37	−0.48	2.18
	The MSSHCF VS SS V28	6.01	−5.21	−0.54	1.97
200~500	The ordinary kriging	2.78	−1.3	0.22	0.67
	The MSSHCF VS EGM2008	1.87	−0.65	0.19	0.47
	The ordinary kriging	16.14	−9.25	1.23	4.1
	The MSSHCF VS DR	14.84	−8.5	1.2	3.88
	The ordinary kriging VS DTU17	3.79	−3.49	0.01	1.2
	The MSSHCF VS DTU17	2.37	−3.06	−0.05	0.98
	The ordinary kriging VS SS V28	−0.49	−2.5	−1.38	0.58
500~1000	The MSSHCF VS SS V28	−0.03	−2.23	−1.17	0.64
	The ordinary kriging	−0.12	−0.69	−0.42	0.18
	The MSSHCF VS EGM2008	−0.07	−0.46	−0.21	0.11
	The ordinary kriging	−1.54	−3.05	−2.3	0.36
	The MSSHCF VS DR	−1.49	−2.57	−2.1	0.32
	The ordinary kriging VS DTU17	−0.63	−2.98	−1.67	0.69
	The MSSHCF VS DTU17	−0.51	−2.75	−1.45	0.71

and

Table 6. The statistical results of the two methods at different sea depth, area B (units: mGal).

Sea Depth	Model Comparison Results	Max	Min	Mean	Std
	The ordinary kriging VS DTU17	12.73	−12.43	−1.03	2.63
	The MSSHCF VS DTU17	8.29	−12.31	−0.43	1.23
0~100 m	The ordinary kriging	12.44	−10.29	−0.77	2.27
	The MSSHCF VS EGM2008	5.83	−7.79	−0.57	1.67
	The ordinary kriging	41.21	−61.04	−3.41	11.92
	The MSSHCF VS DR	34.45	−58.74	−3.21	11.34
	The ordinary kriging VS SS V28	13.48	−15.53	−1.08	4.14
	The MSSHCF VS SS V28	16.43	−10.48	−0.47	3
	The ordinary kriging VS DTU17	5.24	−4.18	−0.69	1.55
	The MSSHCF VS DTU17	3.8	−3.94	−0.23	1.18
100~200 m	The ordinary kriging	3.82	−1.98	−0.75	0.82
	The MSSHCF VS EGM2008	2.53	−2.24	−0.59	0.72
	The ordinary kriging	13.16	−10.67	−2.95	4.19
	The MSSHCF VS DR	13.46	−10.98	−2.8	4.09
	The ordinary kriging VS SS V28	7.13	−10.27	−0.56	2.58
	The MSSHCF VS SS V28	6.99	−9.64	−0.09	2.15
	The ordinary kriging VS DTU17	5.96	−5.83	0.37	2.01
	The MSSHCF VS DTU17	5.29	−4.75	0.17	1.6
200~500 m	The ordinary kriging	2.66	−2.45	0.39	1.07
	The MSSHCF VS EGM2008	2.81	−2.17	0.34	0.93
	The ordinary kriging	14.13	−17.5	2.83	5.4
	The MSSHCF VS DR	14.86	−18.02	2.78	5.25
	The ordinary kriging VS SS V28	9.92	−11.89	0.14	3.43
	The MSSHCF VS SS V28	9.77	−11.45	−0.06	3.03
	The ordinary kriging VS DTU17	5.54	−3.1	1.22	1.76
	The MSSHCF VS DTU17	4.13	−5.12	0.45	1.37
500~1000 m	The ordinary kriging	2.62	−1.97	0.91	0.93
	The MSSHCF VS EGM2008	2.88	−1.96	0.69	0.77
	The ordinary kriging	14.33	−6.61	4.97	4.39
	The MSSHCF VS DR	14.27	−6.38	4.75	4.22
	The ordinary kriging VS SS V28	12.27	−5.32	2.53	2.82
	The MSSHCF VS SS V28	10.14	−5.24	1.76	2.41

.

Table 5 and Table 6 show the comparative results for our research results and different gravity field models with the variations of sea depths. From the statistical results of the two tables, it can be found that with the change of sea depth, the proposed method is more accurate than the ordinary kriging method, and the accuracy of the direct results by using IVM is the least accurate. The absolute Std the results of the MSSHCF and the ordinary method after subtracting the marine gravity models at different sea depth are shown in Figure 10.

Figure 10a,b denote regions A and B, respectively. From the variation in Figure 10a, we find that accuracy of the MSSHCF method is more accurate when compared with the gravity field models at the depth of 0~500 m. Although there is little change of accuracy between 500~1000 m, the overall accuracy is still greatly improved. For Figure 10b, as the sea depth changes, the inversion results based on the MSSHCF method are the most accurate when compared with the gravity field models. Thus, compared with the ordinary method, the MSSHCF method can still optimize the coastal marine gravity field at different sea depths. We infer that this might be because the mean sea surface in the MSSHCF method is approximately equal to the ocean geodetic level, the introduction of the mean sea surface improves the prediction accuracy of the un-measured values.

4.3. Verification against NGDC Shipboard Gravity

So far, our marine gravity field results were compared only with the current gravity field models EGM2008, SS V28, and DTU17, which was only an internal compliance test. Therefore, in external compliance tests our marine gravity field results from NMSSCF and the ordinary kriging were verified against shipborne gravity in areas A and B. The distribution of shipborne gravity in the study area is shown in Figure 11.

As shown in Figure 11, there is a relatively dense shipborne data distribution in both selected areas. There are systematic errors and outliers in the shipborne data due to the influence of sea circulation and various errors. Before using the shipboard gravity data, we processed the shipborne gravity data based on the reference Earth gravity field model EGM2008, which was used as a filter. The specific criteria were as follows: the measurement points with inconsistent values greater than 20 mGal compared with the reference gravity field model EGM2008 were eliminated for shipborne gravity data. After the systematic errors and outliers were removed, area A contained 1621 gravity data points and area B contained 3383 data points. Considering the sparsity of the shipborne data, we sampled the nearest inverse gravity data from the grid where the shipborne data was located for validation. Our marine gravity field results and the different gravity models were compared separately to shipborne gravity; the descriptive statistics are shown in

Table 7. The Shipborne gravity comparison between the inversion results and marine gravity models in the coastal area.

Area	Source of Results		Max	Min	Mean	Std
A	The MSSHCF method	VS shipborne gravity	13.69	−15.94	−4.50	5.39
	The ordinary kriging		14.16	−18.57	−4.35	5.52
	SS V28		15.74	−19.74	−5.19	5.61
	DR		30.74	−33.74	−5.84	9.98
	EGM2008		12.99	−12.85	−4.71	5.47
	DTU17		14.97	−14.65	−4.92	5.36
B	The MSSHCF method	VS shipborne gravity	14.15	−15.90	−1.07	5.99
	The ordinary kriging		17.80	−16.35	−1.05	6.32
	SS V28		17.55	−17.11	−1.20	5.91
	DR		27.66	−17.42	−0.51	6.82
	EGM2008		13.00	−12.99	−1.12	5.84
	DTU 17		15.12	−17.36	−1.32	5.83

.

From the results of the statistical analysis presented in Table 7, we can see that the Std of the results obtained by the MSSHCF method, the ordinary kriging method and directly by using IVM compared with the shipborne gravity are 5.39, 5.52, and 9.98 mGal in area A, respectively. When comparing the proposed MSSHCF method to the shipborne gravity, the Std of SS V28, EGM2008, and DTU17 are 5.61, 5.47, and 5.36 mGal, respectively. In Table 7, comparisons of our inversion results and shipborne gravity data show that the Std for the proposed method is almost equal to that for DTU17 and smaller than the Std for SS V28. Hence, the Std for the proposed method is smaller than that of the ordinary method. In area B, compared with the shipborne data, the Std for the proposed MSSHCF method was 5.99 mGal and for the ordinary kriging method the Std was 6.32 mGal. Compared with the ordinary kriging method, the overall accuracy of the proposed MSSHCF method in the two coastal areas was improved by 0.13 for areas A and by 0.33 mGal for area B.

We compared all gravity results with the shipborne gravity not only to verify the reasonableness of the results obtained by the MSSHCF method but also to show that it has some advantages over the ordinary kriging method. We do not deny the great advantages of global marine gravity model SS V28. However, in coastal areas with sparse data, more accuracy will be lost when the using inversion method to directly obtain the high spatial resolution marine gravity field. Therefore, it may be because the SS V28 model directly uses satellite altimetry data to derive the high spatial resolution marine gravity field in coastal areas.

Scatter plots of the gravity anomalies and shipborne gravity are shown in Figure 12. As shown in Figure 12, The overall agreement R in area A was 83.94% and the ordinary kriging method and 84.55% for the proposed MSSHCF method. In area B, there are 87.06% and 87.38% for two methods. In both areas, the agreement between the gravity field obtained by the proposed MSSHCF method and the shipborne gravity data was greater than the ordinary kriging method. Again, the results from the proposed MSSHCF method were in general agreement with the gravity field models (SS V28 and DTU17).

Using the shipborne gravity for evaluation, although the increased accuracy in the MSSHCF and the ordinary kriging method is comparatively small, it can’t be insignificant. Compared with shipborne gravity, the Std of the gravity field obtained by Zhu after fusing satellite data, such as Jason-1/2/GM, ERS-1/GM, HY-2/GM, and Cryosat-2 was 5.17 mGal [22]. However, the Std of the marine gravity field obtained by Zhu only using Jason-1/GM data is 5.79 mGal, and the overall improved accuracy is only approximately 0.6 mGal [22]. Compared with area A, although the Std of the inversion results in area B is larger than the Std of the gravity field models, this depends on the fact that we use data from only one satellite altimeter in order to validate the proposed method. However, the result also shows that the accuracy of the marine gravity field from the MSSHCF method is more accurate than the ordinary kriging method, which is closer to SS V28.

We combined the results from the proposed method with a comprehensive comparison and analysis of the shipborne coastal gravity field and the results are plausible. We only used Jason-1/GM data and did not consider other types of GM data, such as Geosat/GM, ERS-1/GM, HY-2/GM, and Cryosat-2 satellite altimeter data. However, we only verified that the proposed method is reliable for the inversion of the marine gravity field in the coastal regions, in future work more data will be considered. Our proposed approach may not be sufficiently favorable compared to the optimization of the spatial resolution and accuracy of the gravity field in coastal areas using shipborne or airborne gravity data. Numerical experiments show that the effects introduced by the proposed method appear important for coastal marine gravity optimization. Furthermore, the effects of ocean depth and offshore distance are explored in a localized 2° × 2° coastal area, which may not be rich enough to fully represent the optimization of the proposed method over the ordinary method. Here we accepted this imperfect data processing because our validation results show that the accuracy of the proposed method is more accurate with varying ocean depth and distance offshore. So, when compared with the ordinary kriging method, the proposed method has a certain degree of improvement to nearshore areas.

5. Conclusions

We investigated the effect of the interpolation methods on accuracy when inverting the high spatial resolution marine gravity fields, especially in coastal areas where altimetry data is sparse. We presented an approach based on the ordinary kriging interpolation method, which uses the mean sea surface height as the constraint component in the vertical direction to improve the weight information of the ordinary kriging method and then improved the accuracy of the un-measured values. We used the Jason-1/GM dataset with dense spatial distribution to conduct experiments in the South China Sea region. This data was subjected to error correction and waveform retracking to construct a reliable data set for the experiment [41]. In order to fully verify the validity of the method, we arbitrarily selected two coastal sub-regions with relatively dense ship data distribution; both sea and a small portion of land were included in two areas. We used the DOV method to derive the 2-arc minute residual gravity field of the experimental areas, and applied the combined MSSHCF method and the ordinary kriging method to obtain a 1-arc minute marine gravity field.

The following conclusions can be drawn. Compared to the ordinary kriging method, the experimental results show that the MSSHCF method can better assist in the construction of marine gravity fields at high spatial resolution in coastal regions. In order to fully validate the experimental results, they were evaluated in internal and external tests. The internal test used the gravity field models including the EGM2008, SS V28 and DTU17 and the external test used the shipborne gravity data. These results showed that the accuracy of the coastal marine gravity field produced by the proposed MSSHCF was improved when compared with the ordinary kriging method in two coastal areas. In addition, as the distance of the measurement values from the nearshore and its ocean depth varies, the results from the MSSHCF were more accurate than the ordinary kriging method. When the final results were evaluated with shipborne gravity, accuracy of the MSSHCF method was closer to DTU17 in area A and closer to SS V28 in area B. Thus, this is further evidence that the MSSHCF can deliver accurate coastal marine gravity field estimates in areas with sparse coastal data distribution. The validation experiments showed that the results from the proposed MSSHCF method reduces the loss of accuracy in coastal areas to some extent. This study demonstrates that the MSSHCF method is an innovative and feasible tool for the study of coastal, beach, and island areas where data are sparse and the shipborne and airborne data is lacking, and it demonstrates the validity of the proposed method in the South China Sea region.

Exploration: our results are consistent with some previous studies, where the introduction of variables related to the predicted quantity can optimize the weighting factors of the un-measured points. In order to verify the reliability of the proposed method, the single altimeter satellite data is experimented on in this paper. Thus, the fusion of data from multiple satellites to derive the marine gravity field may also have some influence on the experimental results, which will also be our next work. In addition, in the field of underwater gravity matching navigation, the proposed method can also be well applied to the reconstruction of a high spatial resolution marine gravity field map, and it has good practical significance and popularization value.