Optimal Space Interpolation Method for Continuous Marine Vertical Datum Based on WGS-84 Ellipsoid

Global Navigation Satellite System (GNSS)-based survey technology improves the accuracy and usability of location information by enabling the rapid updating of different information and the efficient use of integration through an ellipsoid. The global marine sector also uses the GNSS to produce rapid and accurate hydrographic surveys and to expand investment and introduce integrated use, such as the S-10X electronic chart. In this study, the conversion of a vertical-based information provision system through a tidal benchmark that is currently implemented to an area-based information provision system was examined. We analyzed a range that allows area modeling using a tidal benchmark by calculating the height from an ellipsoid to the mean sea level (MSL) using the tidal benchmark. An experiment to determine the optimal spatial interpolation was conducted by comparing and analyzing the external verification and performance results. The results of this study are expected to be actively utilized in the manufacture of a continuous marine vertical datum system based on WGS-84.


Introduction
One of the most important current issues today in hydrographic surveys and information production is the use of the ellipsoidal height as the vertical standard in surveys. Recently, the development of Global Navigation Satellite System (GNSS)-based satellite survey technology has enabled the determination of the ellipsoidal height through improved geodetic and survey techniques in various fields, improving the accuracy and processing time of data. Currently, the vertical datum system of each country is divided into land and sea, which makes it difficult to link, integrate, manage, and utilize national spatial information. (1) To address this problem, countries such as the U.S., Britain, and Australia (2)(3)(4)(5) are making continuous efforts to integrate and utilize a dual vertical datum system between land and sea using the ellipsoidal height and geoid. In 2014, the International Federation of Surveyors (FIG) used two types of data from land and hydrographic surveys, and the International Hydrographic Organization (IHO) announced in 2016 that they would use the ellipsoidal height as an official guideline. (6,7) As such, the global trend is to establish a vertical datum system using the ellipsoidal height. (8) For this purpose, in this study, we defined a continuous marine vertical datum as shown in Fig. 1 and conducted an investigation on switching from the vertical datum information system with points through the current tidal benchmark to an area-type information system. Using the tidal benchmark, the height from the ellipsoid to the mean sea level (MSL) was calculated and the optimal spatial interpolation method for reflecting the characteristics of the terrain was determined. (9,10) A spatial data distribution having continuous features is expressed using area-based information by applying a range of spatial interpolation methods based on observed 3D point data. Curtarelli and co-workers performed the mapping of the water depth of the Amazon hydraulic reservoir by spatial interpolation methods, and through cross-validation, they found the conventional kriging interpolation method to be the most suitable. (11,12) Georgas et al. found spline interpolation with barriers to be the most suitable method by comparing actually observed values with a vertical datum after applying interpolation for a vertical datum of the Hudson river and constructing its tide level model. (13) Kim et al. constructed a tide level model showing the relationship between the approximate lowest low water (ALLW) and the regional MSL by spatial interpolation after extracting the semi-range sum of four largeness tide values (Z 0 ) in the tidal benchmark performance. They deduced, through cross-validation, that spline interpolation would be the most suitable method. (14) Jeong et al. constructed a submarine topography model through spatial interpolation using water depth data for the Jeju area. They deduced that kriging interpolation would be the most suitable method and considered that the point density features of the collected data affected its accuracy. (15)

Methodology
In this study, to construct a new ellipsoid-based marine vertical datum, the Gyeong-gi Bay area on the west coast of Korea ( Fig. 2) was selected as the target region in our experiment because of its large tidal difference, the complicated coastal line owing to the presence of bays and islands of various sizes, and its diverse tidal and oceanic currents. If a satisfactory spatial interpolation can be achieved for the west coast, the method used should also be satisfactory for the less complicated south and east coasts.
The height and location criteria of the continuous marine vertical datum of the WGS84 ellipsoid were used and the grid size was set to 5′′ (approximately 160 × 120 m 2 ). The reason for this is that if the grid size is set to 10′′, the resolution will be reduced to a value only suitable for large-scale areas.
The height from the ellipsoid to the MSL was calculated using the information provided by the tidal benchmark performance table for each region obtained by subtracting the MSL elevation from the ellipsoidal height. On this basis, information on tidal benchmark properties, such as the height of the MSL from the ellipsoid, was entered using the ArcGIS and Suffer tools.  Then, basic data, such as coastlines, were extracted from an electronic navigational chart, a polygon file was created, and the range that can be modeled was analyzed through a correlation analysis by the least squares collocation (LSC) method. Each spatial interpolation method was then used to select its parameters.
An experimental verification that comprised an external verification and an evaluation of the performance for verification purposes was performed, and the optimal spatial interpolation method for establishing an ellipsoidal-based vertical datum was determined after the accuracy verification. A correlation distance analysis was performed using 133 performances of the tidal benchmark in 2015 for the west coast where the ellipsoidal height performance had been obtained by the Korea Hydrographic and Oceanographic Agency (KHOA). The height (reference MSL) from the ellipsoid to the MSL was calculated from the tidal benchmark to obtain the ellipsoidal height, and GPS observations for more than 4 h and GPS analysis for Bernese GNSS software were carried out.
It is assumed that when determining the correlation distance using the LSC technique, the differential distribution above the potential depends on the maximum dimension used for the spherical harmonization of a particular potential model. However, the development of specific potential models does not include data on the target area, and if actual observations and potential models used in the LSC method do not fit well, they will change very irregularly. This is because the first order of the magnitude distribution above the potential is very closely related to the error of the potential model coefficients. Therefore, it can be assumed that the differential distribution has a proportional relationship with the error, degree, and variance of the potential model as follows.
where α is a scale factor that must be determined by gravitational field measurement. The differential distribution above the potential can also be determined through a proportional relationship, but it is necessary to determine the covariance function for the potential difference. To address these problems, a differential distribution model, which consists of a functional relationship between the number of dimensions and the number of variances is used. Although various differential distribution models exist, in this study, the following secondary Markov covariance distribution model was used as an analytical covariance function between the reference MSLs for LSC interpolation.
( ) C εε is the correlation distance, C 0 is the variance at 0, and α is the scale factor defining the relationship between the error displacement of the geometric geoid and the error distribution of the gravitational geoid. (16) These are determined automatically and the correlation is determined empirically by the user. Typically, a covariance model for potential anomalies is used as a common Tscherning/Rapp model, which is known to provide the best results for determining the differential distribution of potential anomalies worldwide. However, for the covariance function analyzed in this study, the secondary Markov model was found to be more suitable than the Tscherning/Rapp model owing to the ideal air distribution model between two critical anomalies. Reference 17 provides more details of experimental procedures. Figure 4 shows the variations in the MSL and EGM2008 geoid height over the Korean sea area; they are highly similar regardless of the distance used. Figure 5 shows empirical and analytical covariances of the MSL on the west coast determined through the secondary Markov distribution model. The correlation analysis by the LSC method requires the determination of the relationship between a particular plane and the reference MSL.
Our analysis shows a slightly larger error in the initial variance, which means that the performance of the tidal benchmark in coastal regions does not match that of the geoid height. It was found that the actual covariance and the covariance obtained using the analytical model were relatively consistent and within the range of approximately 0.10 to 0.16° (approximately 16 km). Thus, the scope of surface modeling using the datum level point performance of the west coast region was determined to be up to 16 km.
As spatial interpolation generates discontinuous point data in the form of continuous area data, the spatial interpolation parameter shall be established differently on the basis of the distribution and conditions of the point data to obtain more accurate area data. Therefore, before starting this experiment, an experiment to select the optimal parameter by a different spatial interpolation method [inverse distance weighted (IDW), spline, or kriging interpolation, or spline interpolation with barriers] was performed. For the parameter obtained by each spatial interpolation method, RMSE was deduced and the optimal parameter was selected by comparing the differences between forecast and observed values for every point through crossvalidation. Cross-validation is a method of verifying accuracy through the difference from actually observed values after obtaining forecast values for all the target points while excluding experimental target points present in a certain area one by one.

IDW interpolation
The parameter of IDW interpolation is the power index (distance index), and a parameter selection test was performed for eight power indices (0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.0). The test results are shown in Table 1, and 3.0 was selected as the final parameter. The smallest RMSE means that the difference between the model and actual values is small.

Spline interpolation
The parameter of spline interpolation is mainly determined to be of the regularized or tension type and a selection test for parameters of 0.1, 0.2, 0.3, and 0.4 was performed. The test results are shown in Table 2, and as the final parameter, 0.4 was selected from the regularized type. The smallest RMSE means that the difference between the model and actual values is small.

Kriging interpolation
In the case of kriging interpolation, the parameter is mainly determined to be one of the conventional and two universal types, which are spherical, circular, exponential, Gaussian, and linear, and linear with linear drift and linear with quadratic drift, respectively. The results are shown in Table 3, and the universal-type linear with linear drift parameter was selected. The smallest RMSE means that the difference between the model and actual values is small.

Spline interpolation with barriers (minimum curvature)
In spline interpolation with barriers, the parameter is the smoothing factor, which ranges from 0 to 1. This factor determines how smoothly spatial interpolation is performed. We performed a parameter selection test for parameter values of 0.6, 0.7, 0.8, 0.9, and 1.0. The test results are shown in Table 4, and 1.0 was selected as the smoothing factor because the smallest RMSE means that the difference between the model value and the actual value is small.

Spatial Interpolation Experiment
In this study, we attempted to express the spatial distribution of the regional MSL through surface modeling based on spatial interpolation using the height from the ellipsoid of the tidal benchmark to the MSL.
As the test procedure, the height value (ellipsoid height -MSL) was first estimated from the ellipsoid height and the MSL, which are information provided by each tidal benchmark performance table for the test target region. On the basis of this result, a parameter selection test was performed for each spatial interpolation method (IDW, spline, or kriging, or spline interpolation with barriers) after generating a polygon file by extracting basic data such as the coastline from an electronic navigational chart subsequent to entering attribute data of the tidal benchmark, including the height value of the MSL from the altitude/longitude coordinate, and the ellipsoid by using the ArcGIS tool and Suffer tool. Reference 18 provides more details of experimental procedures. Subsequently, a spatial interpolation test was carried out with two kinds of experimental verification: external verification and comparison of observations. Then the optimal spatial interpolation method was obtained after verifying the accuracy, and an ellipsoid-based marine vertical datum (MSL) was constructed.

Spatial interpolation experiment and comparative validation
The test carried out to determine the optimum method of spatial interpolation for constructing the ellipsoid height-based marine vertical datum was performed for each spatial interpolation method (IDW, spline, or kriging interpolation, or spline interpolation with barriers) with the ArcGIS tool by estimating the height from the ellipsoid to the MSL on the basis of information provided by 67 tidal benchmarks for the target region in the test.
Using 60 of the 67 tidal benchmarks, where those of Gungpyeong port, Deokjeokdo bukri, Eoeundol port, Incheon port, Jumun port, Palmido, and Pungdo port were excluded, spatial interpolation was performed and external validation was carried out. The external validation was performed to verify the accuracy of the model using points not utilized at the time of modeling.
For validation, Gungpyeong port, Deokjeokdo bukri, Eoeundol port, Incheon port, Jumun port, Palmido, and Pungdo port were selected because these locations are strongly affected by oceanic and tidal currents If there is little difference between the forecast and observed values, the reliability of the spatial interpolation can be ensured. There are a few tidal benchmarks in their surroundings considered to be sufficiently accurate for deducing forecast values by minimizing the effect of the surroundings. Tables 5 and 6 show the difference between the observed and forecast values for each location obtained by each spatial interpolation. In the case of IDW interpolation, the difference was between 0.86 and 28.11 and the RMSE was 14.90 cm. In the case of spline interpolation, the difference was between 0.29 and 19.02 cm and the RMSE was 9.57 cm. In the case of kriging interpolation, the difference was between 0.89 and 17.60 and the RMSE was 8.22 cm. In the case of spline interpolation with barriers, the difference was between 0.83 and 19.63 cm and the RMSE was 8.21 cm.

Validation by comparison with observed values
Spatial interpolation was performed using 67 tidal benchmarks and verified by comparison with the observed values (Fig. 6). The validation performance was estimated through actual tide observation and GNSS surveying at locations north of Pungdo and south of Incheon Grand Bridge, where no tidal benchmarks are observed. These locations were selected on the basis of the judgment that if the difference between the forecast and observed values is minimal at locations strongly affected by oceanic and tidal currents, the reliability of spatial interpolation can be ensured; these locations were suitable for validation because tidal benchmarks are evenly distributed around them. The tidal values at these two points are annually revised using yearly data for standard ports (Incheon and Yeongheungdo) after observation for 30 days and nights. (19) Tables 7 and 8 show the difference between the observed and forecast values for each location obtained by each spatial interpolation. In the cases of IDW, spline, and kriging interpolations, and spline interpolation with barriers, the RMSE values were 13.01, 9.12, 9.36, and 9.85 cm, respectively.

Analysis of validation results
The spatial interpolation methods were applied to the test target area and two types of comparison and validation were conducted (external validation and comparison with observed values). The previously analyzed results are summarized in Table 9. As a result of analyzing the general results, we found that the RMSE values for IDW, spline, and kriging interpolations and spline interpolation with barriers are 14.503, 9.470, 8.490, and 8.599 cm, respectively.
Kriging interpolation gave a smaller RMSE than spline interpolations with barriers by 0.109 cm and it may be considered to be more suitable. However, the problem with kriging interpolation is that if a barrierlike coastline impeding physical flow is present in the space to be interpolated, interpolation cannot be performed. If the coastline is not considered, such as when interpolation can be interpolated to land and thus affect the forecast values, barriers such as coastlines must be considered and applied. Considering these features, it was concluded that spline interpolation with barriers using the minimum curvature technique (20) was the most suitable spatial interpolation method in view of the fact that spatial interpolation considering a barrierlike coastline can be performed.
In addition, Table 10 shows the minimum standards of channel surveying specified by IHO. (21) The maximum allowable total vertical uncertainty (TVU) at a confidence level of 95% satisfies a minimum of 26 cm for the special grade, which is the highest grade, and spline interpolation with barriers using the minimum curvature technique, in which spatial interpolation considering the barrierlike coastline may be performed, was considered to be the most suitable spatial interpolation method.

Conversion to each marine vertical datum
Through the previous test, the ellipsoid-based MSL height from the ellipsoid of each tidal benchmark to MSL was constructed by spline interpolation with barriers. In a tide model, when adding and subtracting as much as the semi-range of four largeness tide values (Z 0 ), the marine vertical datum of MSL can be converted to the marine vertical datum of approximate highest high water (AHHW) or the marine vertical datum of ALLW. Among the various tide models, the tidebed system of KHOA, which is applicable to all the waters of Korea, was utilized, and by extracting a range identical to that of the test target region, the conversions to AHHW and ALLW were performed. Figure 7 shows Z 0 and the marine vertical datum of AHHW, MSL, and ALLW from the ellipsoid.

Analysis of ellipsoid-based marine vertical datum
After constructing the ellipsoid-based MSL by spline interpolation with barriers, as a result of comparative analysis with actual observed values to determine whether the conversion to AHHW and ALLW was performed well, we obtained the results in Table 11.
For AHHW and ALLW, the RMSE values were 3.235 and 3.529 cm, respectively, both of which satisfy the allowable error specified by IHO. Therefore, it was confirmed that the conversion to AHHW and ALLW was performed well.

Conclusions
The objectives of this study were to improve the precision of geological surveys using GNSS and to create a continuous marine vertical standard for rapid information updates and the efficient use of integrated information through the representation of locations using ellipsoids. As a result, the method of spatial interpolation was employed to define and produce a continuous marine vertical datum, and verification tests were used to determine and justify  the optimal method of spatial interpolation. We consider that the results can be actively utilized to develop a continuous marine vertical datum based on the WGS-84 ellipsoid in the future. However, in the case of an open sea area where tidal benchmark information is insufficient, it is necessary to analyze physical oceanic information obtained through satellite altimeter data or GNSS buoys and to construct an ellipsoid-based continuous marine vertical datum for the open sea.