Mean Sea Surface Model over the Sea of Japan Determined from Multi-Satellite Altimeter Data and Tide Gauge Records

: Mean sea surface (MSS) is an important datum for the study of sea-level changes and charting data, and its accuracy in coastal waters has always been the focus of marine geophysics and oceanography. A new MSS model with a grid of 1 (cid:48) × 1 (cid:48) over the Sea of Japan and its adjacent ocean (named SJAO2020) (25 ◦ N~50 ◦ N, 125 ◦ E~150 ◦ E) was established. It ingested 12 di ﬀ erent satellites altimeter data (including TOPEX / Poseidon, Jason-1 / 2 / 3, ERS-1 / 2, Envisat, GFO, HaiYang-2A, SRL / Altika, Sentinel-3A, Cryosat-2) and 24 tide gauge stations’ records and joint GNSS data. The latter were used to correct the sea surface height within 10 km from the coastline by using the Gaussian inverse distance weighting method in SJAO2020. The di ﬀ erences among SJAO2020, CLS15, and DTU18, as well as the di ﬀ erences between them and the altimeter data of HY-2A, Jason-3, and Sentinel-3A were introduced. By comparing with tide gauge records, satellite altimeter data, and other models (DTU18, DTU15, CLS15, CLS11 and WHU13), it was demonstrated that SJAO2020 produces the smallest errors, and its coastal accuracy is relatively reliable. the along-track SLA in di ﬀ erent wavelength ranges were obtained by using the fast Fourier transform (FFT). The Sentinel-3B data used in this study was completely independent of all models, whether time or satellite ground tracks. the differences among the MSS models (SJAO2020, DTU18, DTU15, CLS15, CLS11, and WHU13), this study obtained the power spectral density of the along-track sea level anomalies (SLA) between Sentinel-3B and these models. The SLA change of SJAO2020 is 11.76% and 59.24% lower than that of DTU18 and DTU15 in the wavelength range of 35~300 km, respectively. Through the STD analysis of the along-track SLA between different satellite data and each model in the wavelength range of 35 to 300 km, the SJAO2020 model has the smallest error and improved several millimeters. This shows that the method proposed in this study is effective for improving the accuracy of the SSH of the MSS model in the offshore region.


Introduction
Mean sea surface (MSS) refers to the mean dynamic sea surface height (SSH) relative to the reference ellipsoid at a certain period and includes the following two pieces of information: the mean dynamic topography (MDT) and the geoid [1]. MSS, as one of the key parameters of geodesy and oceanography, is widely used in ocean gravity calculations [2][3][4], as well as water depth detection, the determination of geoid fluctuations, and the analysis of crustal deformation [5] which is a key issue in environmental science and earth science today.
Multiple global or regional MSS models have been established from multi-satellite altimeter data. In the early 1990s, Marsh et al. [8] established an MSS model named MSS-9012 with a grid of 1/8 • × 1/8 • between 62 • S and~62 • N from Geos-3 and Seasat data. With the decryption of Geosat data and the launch of other altimeter missions, the temporal and spatial resolution of altimeter data have also been significantly improved, and therefore integrating multi-satellite altimeter data to establish an MSS model has entered a period of rapid development. A series of global high-precision MSS models have been established such as the MSS model CLS01 (2 × 2 ) [9], CLS10 (1 × 1 ) [10], CLS11

Satellite Altimeter Data
The satellite altimeter data used in this study are the Leve2+(L2P) SSH products published by Archiving Validation and Interpretation of Satellite Oceanographic Data (AVISO+). L2P products are generated by the 1 Hz mono mission along-track altimeter data processing segment for T/P, Jason-1/2/3, ERS-1/2, Envisat, GFO, Cryosat-2, SRL, HY-2A, and Sentinel-3A missions [28]. These missions' data have been preprocessed (including quality control and editing of data to select valid ocean data) and corrected for various errors (including instrument errors, environmental perturbations, the ocean sea state bias, the tide effect and atmospheric pressure) [28]. The reference ellipsoid used for ERS-1/2, Envisat, GFO, Cryosat-2, SRL, HY-2A, and Sentinel-3A along-track L2P product is the first-order definition of the non-spherical shape of Earth with (same as for TOPEX/Poseidon, Jason-1/2/3 series) an equatorial radius of 6378.1363 km and a flattening coefficient of 1/298.257 [28]. The effects of ocean tide for the L2P products are corrected by the ocean tide model of FES2014 [24,28]. The corrections for the L2P products are detailed in the Along-track L2P Product Handbook [28].

Satellite Altimeter Data
The satellite altimeter data used in this study are the Leve2+(L2P) SSH products published by Archiving Validation and Interpretation of Satellite Oceanographic Data (AVISO+). L2P products are generated by the 1 Hz mono mission along-track altimeter data processing segment for T/P, Jason-1/2/3, ERS-1/2, Envisat, GFO, Cryosat-2, SRL, HY-2A, and Sentinel-3A missions [28]. These missions' data have been preprocessed (including quality control and editing of data to select valid ocean data) and corrected for various errors (including instrument errors, environmental perturbations, the ocean sea state bias, the tide effect and atmospheric pressure) [28]. The reference ellipsoid used for ERS-1/2, Envisat, GFO, Cryosat-2, SRL, HY-2A, and Sentinel-3A along-track L2P product is the first-order definition of the non-spherical shape of Earth with (same as for TOPEX/Poseidon, Jason-1/2/3 series) an equatorial radius of 6378.1363 km and a flattening coefficient of 1/298.257 [28]. The effects of ocean tide for the L2P products are corrected by the ocean tide model of FES2014 [24,28]. The corrections for the L2P products are detailed in the Along-track L2P Product Handbook [28].
To establish a high-precision MSS model, the ERM data of T/P, Jason-1, Jason-2, Jason-3, ERS-1, ERS-2, GFO, Envisat, HY-2A, SRL, and Sentinel-3A missions are used in this study. To improve the spatial resolution of the MSS model, the GM data of ERS-1/GM, Jason-1/GM, Cryosat-2/LRM (low-resolution mode), SRL/DP (drifting phase), and HY-2A/GM missions are also used. The data used are shown in Table 1. In Table 1, TOPEX/A, Jason-1/A, Jason-2/A, Jason-3/A, and Envisat/A are the ERM data before the orbital transfer of each satellite, while TOPEX/B, Jason-1/B, and Envisat/B are the ERM data after the orbital transfer of each satellite. All ERM data in Table 1 are selected from full-year observations to minimize the temporal oceanic variability in the MSS model.

Tide Gauge Records
The accuracy of the MSS model established solely from satellite altimeter data is seriously affected due to the poor quality of the coastal altimeter data. Tide gauge records have long-term characteristics, relative stability, and high coastal accuracy, and therefore they can be used to improve the coastal accuracy of the MSS model established solely from satellite altimeter data [22]. Thirty-five tide gauges along the coast of Japan on the Permanent Service for Mean Sea Level (PSMSL) have continuous annual records from 1993 to 2018. Among them, only 28 tide gauges (shown in Figure 1) have continuous GNSS joint measurements. Considering that the tide gauge records and GNSS data are consistent in time and space, 24 tide gauges, the black solid circle in Figure 1, are selected to improve the coastal accuracy of the MSS model established only from satellite altimeter data. The remaining four tide gauges, the red solid circle in Figure 1, are used for the validation of coastal accuracy. The annual tide gauge records are derived from the revised local reference (RLR) sea level data released by the PSMSL [29]. The detailed information on 28 tide gauge stations is shown in Appendix A.
The missing values of the annual tide gauge records downloaded from the PSMSL are filled with extreme values of −99999. The missing rate of the 28 tide gauges records is 1.236%. The method of singular spectrum analysis (SSA) iterative interpolation is used to complete the data [30,31].

GNSS Data
The sea level measured by satellite altimeter is relative to the reference ellipsoid. However, the sea level observed from the tide gauges is relative to a certain benchmark. Therefore, there are differences between these two surfaces. Fortunately, the ellipsoidal height of the benchmark can be obtained by GNSS (equipped on the tide gauges) observation, which can be used to unify the sea level obtained by the tide gauges to the reference ellipsoid. Figure 2 shows the relationship between the sea level measured from the satellite altimeter relative to the reference ellipsoid, the sea level measured from the tide gauges relative to a certain benchmark, and the height of the benchmark measured from GNSS relative to the reference ellipsoid.
The GNSS data downloaded from Systèmed' Observation duNiveau des Eaux Littorales (SONEL) is the last solution, named ULR6, completed by the University of La Rochelle (ULR) with GAMIT/GLOBK software. Its reference ellipsoid and frame are GRS80 and ITRF08, respectively, and the baseline processing strategy has been detailed in Santamaria-Gomez et al. [32,33]. For information about GNSS stations, refer to Appendix B. GNSS time series data have made corresponding data corrections [32,33] for earthquakes and other emergencies, and the data over three-times larger than the standard deviation have been eliminated by an iterative outlier detection. Then, the eliminated data are filled by the method of SSA iterative interpolation [30,31]. Since the tide gauge records are the annual records, the GNSS daily data for each year are added and averaged as the GNSS annual data for that year for maintaining the consistency of the time scale between the tide gauge records and GNSS measurements.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 21 the annual records, the GNSS daily data for each year are added and averaged as the GNSS annual data for that year for maintaining the consistency of the time scale between the tide gauge records and GNSS measurements. Figure 2. The relationship between the sea level measured from the satellite altimeter relative to the reference ellipsoid, the sea level measured from the tide gauges relative to a certain benchmark (revised local reference, RLR), and the height of the benchmark measured from GNSS relative to the reference ellipsoid.

Methodology
In this study, the MSS model solely established from multi-satellite altimeter data was established based on the following steps: data selection and preprocessing, collinear adjustment of ERM data, removal of the temporal oceanic variability of GM data, crossover adjustment, and gridding. The 19-year moving average method was used for establishing the MSS model from multisatellite altimeter data, and then the tide gauge records and joint GNSS data were used to correct the SSH of the MSS model in the offshore area (e.g., 10 km from the coastline).

Collinear Adjustment of Exact Repeat Mission (ERM) Data
After precise data editing and preprocessing, the SSHs obtained by satellite altimeter had higher accuracy. However, the instantaneous SSHs still had large fluctuations and contained temporal oceanic variability signals. A collinear adjustment of the ERM data effectively eliminated the temporal oceanic variability of the SSHs with a period shorter than the period of the collinear tracks used and weaken the sea level anomalies (SLA) caused by large-scale ocean anomalies (such as El Niño and La Niña) in a specific period.
One of the collinear tracks with maximum observations was selected as the reference track. After the reference tracks were determined, the SSH of each point of the collinear tracks corresponding to the point of the reference track was computed by collinear adjustment. The method of collinear adjustment used in our study was the same as that described in Jiang et al. [16] and Jin et al. [17] and it was also used to calculate the mean along-track SSH.
In the process of calculating the mean sea surface height (MSSH), the following steps were implemented [25]: (ⅰ) If the difference between SSH and MSSH was larger than 1.0 m, the data were deleted, and the new MSSH was recomputed. (ⅱ) Additionally, it should be ensured that the adjusted data can at least eliminate the time-varying effects of the year, therefore, when the observations participating in the collinear adjustment were less than one year, the point was eliminated. . The relationship between the sea level measured from the satellite altimeter relative to the reference ellipsoid, the sea level measured from the tide gauges relative to a certain benchmark (revised local reference, RLR), and the height of the benchmark measured from GNSS relative to the reference ellipsoid.

Methodology
In this study, the MSS model solely established from multi-satellite altimeter data was established based on the following steps: data selection and preprocessing, collinear adjustment of ERM data, removal of the temporal oceanic variability of GM data, crossover adjustment, and gridding. The 19-year moving average method was used for establishing the MSS model from multi-satellite altimeter data, and then the tide gauge records and joint GNSS data were used to correct the SSH of the MSS model in the offshore area (e.g., 10 km from the coastline).

Collinear Adjustment of Exact Repeat Mission (ERM) Data
After precise data editing and preprocessing, the SSHs obtained by satellite altimeter had higher accuracy. However, the instantaneous SSHs still had large fluctuations and contained temporal oceanic variability signals. A collinear adjustment of the ERM data effectively eliminated the temporal oceanic variability of the SSHs with a period shorter than the period of the collinear tracks used and weaken the sea level anomalies (SLA) caused by large-scale ocean anomalies (such as El Niño and La Niña) in a specific period.
One of the collinear tracks with maximum observations was selected as the reference track. After the reference tracks were determined, the SSH of each point of the collinear tracks corresponding to the point of the reference track was computed by collinear adjustment. The method of collinear adjustment used in our study was the same as that described in Jiang et al. [16] and Jin et al. [17] and it was also used to calculate the mean along-track SSH.
In the process of calculating the mean sea surface height (MSSH), the following steps were implemented [25]: (i) If the difference between SSH and MSSH was larger than 1.0 m, the data were deleted, and the new MSSH was recomputed. (ii) Additionally, it should be ensured that the adjusted data can at least eliminate the time-varying effects of the year, therefore, when the observations participating in the collinear adjustment were less than one year, the point was eliminated.

Temporal Oceanic Variability Corrections of GM Data
Since the GM data did not have the characteristics of repeated cycles, it was impossible to remove the ocean variability by collinear adjustment method. Fortunately, for the period of GM data (e.g., ERS-1/GM span 1994-1995), the ocean variability has been simultaneously observed by ERM data (e.g., T/P span 1994-1995). Therefore, the ocean variability of GM data could be corrected by the ocean variability of the ERM data which was considered as a reference at the spatial and temporal positions of GM data. Altimeter data used for temporal ocean variability corrections of GM data are listed in Table 2. The method of optimal analysis (OA) [10,11] was used in our study. OA was used to interpolate the SLA of one or more missions considered to be a reference at the spatial and temporal position of the satellite that would be corrected for ocean variability. Table 2. Corresponding data used for sea level variation corrections of geodetic mission (GM) data. The altimeter data of T/P, Jason-1, Jason-2, and Jason-3 are acknowledged as having the highest orbit and measuring accuracy. Therefore, the mean along-track SSH of uninterrupted joint T/P + Jason-1 + Jason-2 + Jason-3 (hereafter T/P series) of each group data (shown in Appendix B) was the fundament to calculate the SLA of T/P, Jason-2 and Jason-3. The former was used to correct the ocean variability of a given SSH of GM data.

Crossover Adjustment
Long-wavelength sea level changes of satellite observations can be weakened by collinear adjustment and temporal oceanic variability correction, such as the radial orbit error and the temporal variability of SSH. However, the residual radial orbit error, the short-wavelength signal of temporal oceanic variability and geophysical correction residuals, are still the main influences on the determination of MSS. The crossover adjustment is based on the difference between two observations at the same point to integrate different satellite altimeter data (including ERM data and GM data) or to determine corrections to measurements [25,34,35].
The classical crossover adjustment regards the radial orbit error as one of the dominant sources of errors affecting altimeter data and that error can be sufficiently modelled by either a time-or a distance-dependent polynomial [36][37][38]. With the improvement of the precision orbit determination technology, the radial orbit error of the new generation of satellite altimeter data has been effectively controlled. Therefore, the radial orbit error is no longer the main factor that affects the accuracy of the altimeter data but is a comprehensive effect of the same magnitude as other errors, such as short-wave ocean time-varying signals and geophysical correction residuals.
In this study, the posterior compensation method was used [34,35]. First, the conditional adjustment was used to adjust the crossover observation equation, and then a new error model was used to filter and predict the SSH along the track.
When the condition adjustment is carried out, the SSH at any point along the track can be expressed as: where h is the SSH observation, h 0 is the true value of h, and ∆ is the observation error (systematic error and random error). According to Equation (1), at the crossover point of the track i and j, the conditional equation can be established as: where d ij is the SSH difference at the crossover point. For a track network with multiple crossover points, the matrix form of the crossover point condition equation can be written as: where A is the coefficient matrix, which consists of 1 and −1; V is the correction vector of the observation error, and W is the difference vector of crossover points. The least squares solution of Equation (3) is: and the corresponding cofactor matrix is: where P is the weight matrix of SSH observations. If each observation point on the track is an independent observation, we can derive the following: where p a ij and p d ij represent the weight factors of observations at the crossover point of the track i and j. In the crossover adjustment of single satellite, p a ij = p d ij = 1/2. In the crossover adjustment of multi-satellite, matrix P is constructed by 1/ √ 2 times STD of the single-satellite crossover differences. The comprehensive effect of residuals changes is very complicated, including parts that change linearly and parts that change periodically, and more are parts that have more complicated changes. According to Wagner [36] and Rummel [37], the traditional error model was extended to the following mixed polynomial model with respect to the observed time as the independent variable [35]. It can be expressed as where t is the observed time of SSH; a 0 , a 1 , b 1 , and c i (i = 1, 2, · · · , m) are the undetermined coefficients to be evaluated, and m is a positive integer which is determined with the length of the track. Here, m is proposed to be 1~2 for a short track, 3~5 for a middle-long track, and 6~8 for a long track by experience [35]; w represents the angular frequency corresponding to the periodic change of the error, which can generally be expressed as: where T 0 and T 1 represent the corresponding observation times at the beginning and end of the track, respectively. After the conditional adjustment, V can be regarded as a kind of virtual observation vector. Equation (7) is used as the error model, and then an error equation can be established at the crossover point. This error equation is the following: where v is the virtual observations and δ is the correction of the virtual observations. The matrix form of Equation (9) is as follows: where l is the correction vector of the virtual observations; B is a matrix of known coefficients; X is the vector of the undetermined coefficient; and V is the virtual observation vector. The least squares solution of Equation (10) is as follows:X where P V is the weight matrix of virtual observations. The estimated parameter vectorX is put into Equation (7). According to the observation time of the along-track SSH of the track, the residuals of SSH systematic error can be calculated and corrected by using Equation (7).

Least-Squares Collocation Technique for Gridding
Most of the orbital errors and residual temporal oceanic variability errors are weakened by the crossover adjustment. However, the differences at the crossover point after adjustment indicate that there is still residual error, especially the GM observations. The least-squares collocation (LSC) technique for gridding can effectively use the prior information of the observations to solve the optimal estimate of the interpolated value [39,40], that is, different weights are given according to the residual error after adjustment to reduce the impact on the accuracy of the SSH.
When using the LSC to grid the SSH, the signal must have a zero-mean characteristic [17,25]. Therefore, the reference MSS model must first be removed from the satellite altimeter data after crossover adjustment. In this study, the CLS15 MSS model in the L2P data product was selected as the reference MSS model. Then, the residual SSH was gridded. Before gridding, the average of the residual SSH should be subtracted to satisfy the zero-mean characteristic. Finally, the MSS model was obtained by adding the grid value to the average and restoring the removed reference MSS model. Suppose a certain observation vector y contains a zero mean signal t and a zero-mean noise vector v and is expressed as: where the self-covariance matrices of t and v are C tt and C vv , respectively; there is no correlation between t and v, that is, C tv = 0. Using the LSC technique, for any zero-mean signal s in the data distribution, the fitted value [17] is:ŝ where C st is the cross-covariance between signal s and signal t, and again, there is no correlation between s and v. If the self-covariance C ss of the signal s is known, the estimation error of s can be expressed as: When the fitting point and the observation point coincide, the observation point is considered to have no error, that is, C st = C tt = C ss , C vv = 0. From Equations (13) and (14), the estimated valueŝ = t and the error Eŝŝ = 0 satisfy the general interpolation method. Therefore, when the cross-covariance between the a priori information and the signal as well as the error is accurately known, LSC can effectively use the a priori information of the observations to obtain accurate interpolation values.
Because the satellite altimeter observations are very large and the data around each network point are densely distributed, they are insensitive to the accuracy of the covariance function when determining the MSSH of each network point [41]. Therefore, a second-order Markov process was used to describe the one-dimensional covariance function [42]. The process can be expressed as [43]: where d is the two-dimensional distance between the observation point and the grid point; C 0 is the local variance parameter, which can be expressed by the variance of all observations in the local range involved in the grid. α = 0.595ξ, α is the parameter to be estimated and ξ is the correlation length (where a 50% correlation is obtained) [44], and the value is 70 km [45,46]. The single-satellite crossover differences accuracy of 1/ √ 2 times after the crossover adjustment is introduced into the LSC as the noise of the corresponding satellite data.

Nineteen-Year Moving Average Method
Almost all altimeter data from 1993 to 2018 with a period of 26 years were used in the establishment of the MSS model. To further eliminate the influence of the residual errors of tide model error on the altimeter data (in Table 1), these data were put into 8 groups over 19-year-long moving windows, starting in 1993 and shifted by one year. The 8 sets of data were independently established an MSS model, respectively, and eight MSS models were obtained. Then, the final MSS model was determined by calculating the average at each grid point of these eight MSS models.
Also, the SSHs (span 1993-2018) obtained from the corrected tide gauge records and the joint GNSS data were grouped with the same method used in altimeter data. In this way, eight groups SSHs of each tide gauge station (a total of 28 tide gauges) at different periods were obtained. Then, the average value of the 8 SSHs of the tide gauges was calculated for the SSH of that tide gauge, and used to correct the SSHs of the MSS model in the offshore region.
The calculating process is described as follows [25]: where SSH i,s is the SSH at the grid point i in the SJAO2020A (see Section 3.1.3) model or at the i-th of the 28 tide gauges, and ssh i,j is the SSH at the grid point i in each of the eight models or the SSH of the eight groups SSH of the i-th tide gauge. The mean along-track SSH of uninterrupted joint T/P series over the period of each group was used as fundament, for example, the mean along-track SSH of uninterrupted joint T/P series during 1992-2011 was the fundament for the first MSS model. The SJAO2020A model was derived from the average of eight MSS models. Therefore, the fundament for SJAO2020A model was the average of the fundaments for the eight MSS models, i.e., the mean along-track SSH of uninterrupted joint T/P series during 1993-2018 could be considered to be the fundament for the SJAO2020A model.

The Method for Improving the Coastal Accuracy of the Model
First, inverse barometer correction was performed on the tide gauge records, after preprocessing described earlier. Then, the SSHs (span 1993-2018) obtained from the corrected tide gauge records and the joint GNSS data were adjusted to have the same reference ellipsoid and frame as T/P. Finally, the SSH of each tide gauge station was obtained with the 19-year moving average method. These SSHs were used to correct the SSHs of the grid points of the SJAO2020A model within 10 km from each tide gauge station. Regarding why the data of 10 km away from the coastline was selected to be corrected, please see Appendix C for specific reasons.
The method of Gaussian inverse distance weighting was used in the process of the corrections. The correction value ∆h of each grid point is expressed as: where H tide is the tide gauge SSH, H ssh is the SSH of the grid point of the SJAO2020A model within 10 km from each tide gauge station, and p is a weighting factor weighted by the Gaussian inverse distance, which can be expressed as: where α is the Gaussian distance smoothing factor, which is 10 km in this study and d is the spherical distance from the grid point to the tide gauge station. The calculation equation is: where a is the semimajor axis of the T/P reference ellipsoid; y t and x t are the latitude and longitude of the tide gauge station; y s and x s are the latitude and longitude of the grid point; and ∆y = (y t − y s )/2, and ∆x = (x t − x s )/2. Finally, in the study area, the correction of the coastal (10 km from the coastline) grid points of the MSS model was realized by the spline interpolation method [47] of adjustable tensor continuous curvature in the Generic Mapping Tools (GMT) software [48].

Correction of Temporal Oceanic Variability
To validate the effect of correction of temporal oceanic variability, the crossover differences of SSH before and after temporal oceanic variability correction are separately counted. The statistical results are shown in Table 3.  Table 3 shows that the accuracy of the ERM data is greatly improved after collinear adjustment. This shows that the temporal oceanic variability signal has a significant influence on the MSS. Collinear processing can weaken the effect of temporal oceanic variability on ERM data and improve (better than 10 cm) the calculation accuracy of the SSH. The T/P series, as the spatial-temporal fundamental, has the highest accuracy of crossover differences. According to Table 2, temporal oceanic variability correction is performed on the GM data. The effect of temporal oceanic variability correction is shown in Table 3, in which the crossover difference is improved by approximately 7 cm.

The Results of the Crossover Adjustment
In the crossover adjustment of multi-satellite, the mean along-track SSH of uninterrupted joint T/P series is used as fundament. Satellite observations with high orbit accuracy improve that with low orbit accuracy by combining all satellite data with the method of crossover adjustment. The combined effects of the above radial orbit errors and other errors are further eliminated, realizing the unity and coordination of the combined processing of multiple satellite altimeter data.
It can be seen from Tables 3 and 4 that after the self-crossover adjustment between the satellites, the STD of the crossover difference of SSH in ERM data is at the level of 1~4 cm and the GM data is approximately 9 cm, which effectively reduces the residual radial orbit errors of each satellite. The STD of the crossover difference after the joint crossover adjustment of each satellite is 0.0811 m, realizing the improvement of satellite observations with high accuracy to the satellite observations with low accuracy. The T/P series data are used as fundament, and the STD of the self-crossover difference is 0.0068 m, which is smallest as compared with other satellites. Therefore, it is still in a dominant position in the joint crossover adjustment, and the coordination between multiple satellites is achieved, while the benchmark is unchanged.

Establishment of the Model
The LSC method is a statistical estimation method based on the observation covariance information that takes full account of the statistical correlation between the data; the smoothing function is better than Shepard, continuous curvature tension spline, and other analytical methods, and therefore it is more suitable to grid satellite altimeter along-track data after crossover adjustment [17]. Eight MSS models are obtained from the eight groups of altimeter data in Table A1. The MSS model over the Sea of Japan and its adjacent ocean was determined from multi-satellite altimeter data (named SJAO2020A) by taking the average at the grid point of these eight MSS modes. To validate the accuracy of the SJAO2020A model, it is compared with DTU18, DTU15, Whu13, CLS15, and CLS11 within different distances from the coastline, respectively. The STD of the differences of the SSH between SJAO2020A and DTU18, DTU15, Whu13, CLS15, and CLS11, within 10 km, 10-20 km, 20-30 km, 30-40 km, and 40-50 km from the coastline, are shown in Figure 3.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 21 SSH between SJAO2020A and DTU18, DTU15, Whu13, CLS15, and CLS11, within 10 km, 10-20 km, 20-30 km, 30-40 km, and 40-50 km from the coastline, are shown in Figure 3. It can be seen from Figure 3 that the STDs of the differences of the SSH between SJAO2020A and CLS15 are always the smallest, and the largest is that of CLS11. The CLS11 model was established using only the satellite data from 1993 to 2009, and the GM data were only used from the ERS-1 satellite from 1994 to 1995. The STD of the difference between SJAO2020A and the other four models decreases with increasing distance from the coastline, and the largest decrease is within the range of 10~20 km. After 40 km, the STD is relatively stable, at approximately 0.025 m, indicating that the accuracy of SJAO2020A is relatively stable and reliable. However, these MSS models are very different in coastal areas, which is caused by the low coastal accuracy of satellite altimeter data. Therefore, improving the coastal accuracy of SSH is extremely important for improving the overall accuracy of the MSS model. In this study, the coastal accuracy of SJAO2020A was improved using 24 tide gauge stations and joint GNSS around the Japanese coastline.

Improvement of Model Coastal Accuracy
The Gaussian inverse distance weighting method was used to correct the SSH of grid points within 10 km from the 24 tide stations. Then, the adjustable tensor continuous curvature spline interpolation method was used to correct all coastal grid points of SJAO2020A. The final MSS model (reference to T/P reference ellipsoid) over the Sea of Japan and its adjacent ocean (named SJAO2020) was established, as shown in Figure 4. Figure 4 shows that the SSH change over the Sea of Japan and its adjacent ocean areas is complicated. The overall SSH of the Philippine plate is higher than that of the Pacific plate. Among them, the SSH of the Izu-Ogasawara Trench, Japan Trench, and Thousand Islands Trench is low. The lowest point is the Thousand Islands Trench, which is −2.826 m. The highest point at SSH is 52.9643 m, near the Ogasawara Islands.
To verify the coastal accuracy of SJAO2020 corrected by the tide gauge station, the SSHs of DTU18, CLS15, SJAO2020, and SJAO2020A at the position of the tide gauge station (the red solid circle in Figure 1) were interpolated and compared with the actual measured SSHs of the tide gauge stations. The statistical results are shown in Table 5. It can be seen from Figure 3 that the STDs of the differences of the SSH between SJAO2020A and CLS15 are always the smallest, and the largest is that of CLS11. The CLS11 model was established using only the satellite data from 1993 to 2009, and the GM data were only used from the ERS-1 satellite from 1994 to 1995. The STD of the difference between SJAO2020A and the other four models decreases with increasing distance from the coastline, and the largest decrease is within the range of 10~20 km. After 40 km, the STD is relatively stable, at approximately 0.025 m, indicating that the accuracy of SJAO2020A is relatively stable and reliable. However, these MSS models are very different in coastal areas, which is caused by the low coastal accuracy of satellite altimeter data. Therefore, improving the coastal accuracy of SSH is extremely important for improving the overall accuracy of the MSS model. In this study, the coastal accuracy of SJAO2020A was improved using 24 tide gauge stations and joint GNSS around the Japanese coastline.

Improvement of Model Coastal Accuracy
The Gaussian inverse distance weighting method was used to correct the SSH of grid points within 10 km from the 24 tide stations. Then, the adjustable tensor continuous curvature spline interpolation method was used to correct all coastal grid points of SJAO2020A. The final MSS model (reference to T/P reference ellipsoid) over the Sea of Japan and its adjacent ocean (named SJAO2020) was established, as shown in Figure 4. Figure 4 shows that the SSH change over the Sea of Japan and its adjacent ocean areas is complicated. The overall SSH of the Philippine plate is higher than that of the Pacific plate. Among them, the SSH of the Izu-Ogasawara Trench, Japan Trench, and Thousand Islands Trench is low. The lowest point is the Thousand Islands Trench, which is −2.826 m. The highest point at SSH is 52.9643 m, near the Ogasawara Islands.
To verify the coastal accuracy of SJAO2020 corrected by the tide gauge station, the SSHs of DTU18, CLS15, SJAO2020, and SJAO2020A at the position of the tide gauge station (the red solid circle in Figure 1) were interpolated and compared with the actual measured SSHs of the tide gauge stations. The statistical results are shown in Table 5. represent the SSH at the tide gauge station calculated by SJAO2020A and SJAO2020, respectively. DTU18 and CLS15 are a few centimeters away from the sea level at tide gauge stations Ogi, Ito II, and Kushimoto, and are about 4 dm away from the sea level at tide gauge station Tajiri, while the sea level difference between SJAO2020a and tide gauge station Tajiri is about 4 cm. The SSH difference between SJAO2020 and the four tide gauge stations is approximately 2 cm, which improves the SSH difference between SJAO2020A and the tide station. This indicates that the Gaussian inverse distance weighting method based on tide gauge stations and the joint GNSS can effectively improve the coastal accuracy of the SJAO2020 MSS model.

Accuracy Assessment of SJAO2020
Usually, reliability and accuracy are validated by comparing with mean along-track altimetry data and other models [1]. The difference between MSS models depends on the dataset used for calculation and the data processing method [11]. To evaluate the model error of SJAO2020 at different wavelength scales and to better quantify the difference between it and DTU18, DTU15, CLS15, CLS11, and WHU13, based on the mean along-track ERM data of the one-year uninterrupted Sentinel-3B satellite in 2019, the along-track SLAs between the Sentinel-3B and SJAO2020, DTU18, DTU15, CLS15, CLS11, and WHU13 models were obtained. The power spectral density (PSD) of the along-track SLA in different wavelength ranges were obtained by using the fast Fourier transform (FFT). The Sentinel-3B data used in this study was completely independent of all models, whether time or satellite ground tracks.  In Table 5, H DTU15 and H CLS15 represent the SSH at the tide gauge station calculated by DTU18 and CLS15 model, respectively. H SJAO2020A and H SJAO2020 represent the SSH at the tide gauge station calculated by SJAO2020A and SJAO2020, respectively. DTU18 and CLS15 are a few centimeters away from the sea level at tide gauge stations Ogi, Ito II, and Kushimoto, and are about 4 dm away from the sea level at tide gauge station Tajiri, while the sea level difference between SJAO2020a and tide gauge station Tajiri is about 4 cm. The SSH difference between SJAO2020 and the four tide gauge stations is approximately 2 cm, which improves the SSH difference between SJAO2020A and the tide station. This indicates that the Gaussian inverse distance weighting method based on tide gauge stations and the joint GNSS can effectively improve the coastal accuracy of the SJAO2020 MSS model.

Accuracy Assessment of SJAO2020
Usually, reliability and accuracy are validated by comparing with mean along-track altimetry data and other models [1]. The difference between MSS models depends on the dataset used for calculation and the data processing method [11]. To evaluate the model error of SJAO2020 at different wavelength scales and to better quantify the difference between it and DTU18, DTU15, CLS15, CLS11, and WHU13, based on the mean along-track ERM data of the one-year uninterrupted Sentinel-3B satellite in 2019, the along-track SLAs between the Sentinel-3B and SJAO2020, DTU18, DTU15, CLS15, CLS11, and WHU13 models were obtained. The power spectral density (PSD) of the along-track SLA in different wavelength ranges were obtained by using the fast Fourier transform (FFT). The Sentinel-3B data used in this study was completely independent of all models, whether time or satellite ground tracks.
As shown in Figure 5a, the error of each model is obviously different for wavelengths longer than 300 km. This is due to the independence of each model and the systematic error between each model included in the SLA time series. Among them, the systematic error between SJAO2020 and DTU18, as well as CLS15, is relatively small due to the types and quantities of SSH data used by the three are the closest. The error of each model drops rapidly at wavelengths of 250~300 km. At wavelengths of 35~250 km, the SJAO2020 error is significantly improved as compared with other models. In the smaller wavelength range of 12~18 km, the errors of each model have a small increase, of which satellite altimeter noise is dominant [12]. Figure 5b shows the ratio of the along-track PSD of DTU15, DTU18, and CLS15 in Figure 5a to the along-track PSD of SJAO2020. Compared with DTU18 and DTU15, the SLA change of SJAO2020 is reduced by 11.76% and 59.24%, respectively in the wavelength range of 35~300 km, and the improvement is greatest when the wavelength is 66 km. The SLA change of SJAO2020 is reduced by 8.91% as compared with CLS15 in the wavelength range of 35~300 km.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 21 As shown in Figure 5a, the error of each model is obviously different for wavelengths longer than 300 km. This is due to the independence of each model and the systematic error between each model included in the SLA time series. Among them, the systematic error between SJAO2020 and DTU18, as well as CLS15, is relatively small due to the types and quantities of SSH data used by the three are the closest. The error of each model drops rapidly at wavelengths of 250~300 km. At wavelengths of 35~250 km, the SJAO2020 error is significantly improved as compared with other models. In the smaller wavelength range of 12~18 km, the errors of each model have a small increase, of which satellite altimeter noise is dominant [12]. Figure 5b shows the ratio of the along-track PSD of DTU15, DTU18, and CLS15 in Figure 5a to the along-track PSD of SJAO2020. Compared with DTU18 and DTU15, the SLA change of SJAO2020 is reduced by 11.76% and 59.24%, respectively in the wavelength range of 35~300 km, and the improvement is greatest when the wavelength is 66 km. The SLA change of SJAO2020 is reduced by 8.91% as compared with CLS15 in the wavelength range of 35~300 km. As seen from the above, satellite altimeter data are also an effective method for evaluating the MSS model. The SLA of satellite altimeter data and five models were obtained separately. Among them, the satellite altimeter data include collinear data of the T/P series, ERS-1, HY-2A, and Sentinel-3B, as well as GM data of Jason-2 and SRL (due to technical problems, the SRL satellite drifted in repetitive periodic tracks in March 2015; compared with historical repetitive tracks, the maximum drift amount is up to 10 km). The root mean square error of each SLA in the wavelength range of 35 to 300 km was calculated by Chebyshev bandpass filtering. The statistical results are shown in Table  6. As seen from the above, satellite altimeter data are also an effective method for evaluating the MSS model. The SLA of satellite altimeter data and five models were obtained separately. Among them, the satellite altimeter data include collinear data of the T/P series, ERS-1, HY-2A, and Sentinel-3B, as well as GM data of Jason-2 and SRL (due to technical problems, the SRL satellite drifted in repetitive periodic tracks in March 2015; compared with historical repetitive tracks, the maximum drift amount is up to 10 km). The root mean square error of each SLA in the wavelength range of 35 to 300 km was calculated by Chebyshev bandpass filtering. The statistical results are shown in Table 6. As shown in Table 6, the STD of the SLA given by SJAO2020 and its benchmark (T/P series) data in the wavelength range of 35 to 300 km is the smallest as compared with the other five models, indicating that the model is the most stable and the data processing results are more reliable. The STD of the along-track SLA given by between HY-2A and SJAO2020 is significantly smaller than that of the other five models. This is because HY-2A data are only used in SJAO2020, indicating that the accuracy of the MSS model at the ground track of the satellite can be improved by adding more high accuracy satellite observations. Sentinel-3B, SRL, and Jason-2/GM are not used in the establishment of the SJAO2020, CLS15, DTU18, DTU15, WHU13, and CLS11 models, and the STDs of the three-satellite data along-track SLA are sequentially reduced. This shows that the accuracy of the SJAO2020 model is the highest. The accuracy of CLS15 is equivalent to that of DTU18. As compared with DTU15, the accuracy of DTU18 has been improved, and the accuracy of both models is better than WHU13. However, the ERS-1 along-track data are consistent with CLS11, only second to SJAO2020, and better than CLS15, DTU18, DTU15, and WHU13. This may be caused by the establishment of the CLS11 model using only the GM data of one ERS-1 satellite. The model has a strong correlation with the ERS-1 satellite.

Conclusions
A new MSS model named SJAO2020A with a grid of 1 × 1 over the Sea of Japan and its adjacent ocean was established with the 19-year moving average method by combining satellite altimeter data from 1993 to 2018. Different from the latest MSS models CLS15 and DTU18, the measured data of the latest altimetry satellites HY-2A, Jason-3, and Sentinel-3A are ingested in SJAO2020A. To improve the coastal accuracy of SJAO2020A, 24 tide gauges and the joint GNSS stations along the coast of Japan are used to correct the sea surface height within 10 km from the coastline by using the Gaussian inverse distance weighting method. Then, the SJAO2020 model with higher coastal accuracy is obtained. The difference between SJAO2020 and the four tide gauge stations (the red solid circle in Figure 1) is approximately 2 cm, and SJAO2020 coastal accuracy is better than that of CLS15 and DTU18.
To better quantify the differences among the MSS models (SJAO2020, DTU18, DTU15, CLS15, CLS11, and WHU13), this study obtained the power spectral density of the along-track sea level anomalies (SLA) between Sentinel-3B and these models. The SLA change of SJAO2020 is 11.76% and 59.24% lower than that of DTU18 and DTU15 in the wavelength range of 35~300 km, respectively. Through the STD analysis of the along-track SLA between different satellite data and each model in the wavelength range of 35 to 300 km, the SJAO2020 model has the smallest error and improved several millimeters. This shows that the method proposed in this study is effective for improving the accuracy of the SSH of the MSS model in the offshore region.   . 816-517). AVISO+ for providing the along-track L2P products and the waveform data of Jason-1, which can be obtained by downloading freely from AVISO+ official website (ftp://ftp-access.aviso.altimetry.fr). The tide gauge records are available online (https://www.psmsl.org/) and the GPS data are available online (https://www.sonel.org).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix C
Altimeter waveforms are usually contaminated due to land, island, sea reef, sea ice, seabed terrain, etc. If the extracted ranges from these corrupted waveforms are used, sea levels calculated from these ranges are also incorrect [20]. However, for the seas of Japan and its adjacent ocean, how far away from the coastline is the altimeter data that is incorrect and cannot be used? To solve this question, six arcs ( Figure A1) of Jason-1 satellite from sea to land or from land to sea are selected from the study area. The return power of each arc within 30 km from the coastline is shown in Figure A2. Altimeter waveforms are usually contaminated due to land, island, sea reef, sea ice, seabed terrain, etc. If the extracted ranges from these corrupted waveforms are used, sea levels calculated from these ranges are also incorrect [20]. However, for the seas of Japan and its adjacent ocean, how far away from the coastline is the altimeter data that is incorrect and cannot be used? To solve this question, six arcs ( Figure A1) of Jason-1 satellite from sea to land or from land to sea are selected from the study area. The return power of each arc within 30 km from the coastline is shown in Figure A2.  As shown in Figure A2 a, c-e, the waveforms of Jason-1 began to be contaminated approximately 10 km from the coastline. The closer the waveform is to the coastline, the more serious the waveform is contaminated and the greater the accuracy of sea level observations. There is no obvious waveform contamination of arcs b and f in cycle 18, but the two arcs only have waveform data after 9.361 km and 6.596 km from the coastline. Therefore, the tide gauge stations are mainly used to correct the SSH of the MSS model 10 km away from the coastline. As shown in Figure A2a,c-e, the waveforms of Jason-1 began to be contaminated approximately 10 km from the coastline. The closer the waveform is to the coastline, the more serious the waveform is contaminated and the greater the accuracy of sea level observations. There is no obvious waveform contamination of arcs b and f in cycle 18, but the two arcs only have waveform data after 9.361 km and 6.596 km from the coastline. Therefore, the tide gauge stations are mainly used to correct the SSH of the MSS model 10 km away from the coastline.