Inversion and Validation of Improved Marine Gravity Field Recovery in South China Sea by Incorporating HY-2A Altimeter Waveform Data

Abstract

HaiYang-2A (HY-2A, where ‘Haiyang’ means ‘Ocean’ in Chinese) has provided reliable sea surface height observations for gravity with uniform ocean data coverage on a global scale for more than 8 years, particularly with denser across track sampling during the geodetic mission since March 2016. This paper aims at modeling and evaluating the regional marine gravity field at 1′×1′ resolution by incorporating HY-2A altimeter waveform data from 7 complete 168-day cycles in the geodetic mission phase. Initial evaluation indicates that, firstly, the measurements in the geodetic mission stay at a consistent accuracy level with observations at the start-of-life stage according to statistics of discrepancies at crossover points cycle by cycle. Secondly, range precision improvement can be achieved using a two-pass weighted least-squares retracker. Thirdly, a downsampling procedure combined with a low-pass filter is designed for HY-2A 20 Hz data to obtain 5 Hz measurements with enhanced precision. We calculate the 1′×1′ marine gravity field model over the South China Sea area by using the EGM2008 model as a reference field with the remove/restore method. The verifications with published models and shipborne gravimetric data show that HY-2A GM data is capable of improving marine gravity field modeling. Results show slightly higher accuracy than other models with similar input datasets but not including HY-2A. The accuracy is also compared with the latest DTU17 and SIO V27.1 model.

1. Introduction

Accurate knowledge of marine gravity anomalies is essential for a range of scientific and engineering applications, such as earth gravity modeling, oceanic resources exploration, submarine topography inversion, and yielding insight into the structures of the lithosphere [1,2,3,4]. Over the past few decades, satellite altimetry has evolved as an effective tool for recovering regional or global marine gravity fields. The gravitational parameters (e.g., geoid height, vertical deflection, gravity anomaly, gravity gradient,) can be retrieved or deduced from altimeter measurements of sea surface heights or slopes [5,6,7,8,9,10]. The improvement of the marine gravity field model mainly relies on the improvement to range precision and spatial resolution of accumulated multi-satellite altimeter measurements, especially those from geodetic missions (GM) [6,11].

Early studies show that retracking techniques significantly improves the accuracy of the recovered marine gravity field, even in marginal seas and sea-ice covered regions [12,13,14,15,16]. In addition, marine gravity field modeling has gone through huge improvement since the completion of the Jason-1 geodetic mission and accumulation of the CryoSat-2 data with a repeat period around 369-day [11,17]. Recently, geodetic missions (e.g., HaiYang-2A (HY-2A), SARAL/AltiKa Drifting Phase and Jason-2 GM) are expected to improve gravity recovery.

HY-2A is China’s first satellite for exploring the oceanic dynamic environment. HY-2A has been running operationally in orbit for almost 8 years, providing a long series of measurements from exact repetitive mission (~4.5 years) and subsequent geodetic mission (>3 years). HY-2A supplement other altimeter missions with new ground tracks due to its orbit inclination of 99.34° and different revisiting period, and potentially give detective signals over somewhere unreached before. HY-2A started its geodetic mission in March 2016, and the revisiting period extended to 168 days from the original 14 days during exact repeat mission (ERM).

Zhang et al. (2018) has verified that HY-2A provides reliable measurements and shown that the two-pass waveform retracker enhances the range precision by a factor of 1.6 [18]. A previous study, conducted by Zhu et al. (2019) indicated that HY-2A can provide similar marine gravity recovery to other altimeters on the basis of 1 Hz measurements [19]. Here we investigate if the 20 Hz waveform data records can be used to extract valid measurements with a relatively higher rate (e.g., 5 Hz) by adopting the two-pass weighted least-square retracker, such that a refined marine gravity can be derived.

The South China Sea (SCS) (105°~125°E, 4°~26°N) has a relatively complex seafloor topography, and this paper aims to assess the regional marine gravity modeling by incorporating HY-2A GM altimeter waveform data records in to multi-altimeter data sets. The method of constructing a marine gravity model from vertical deflection and calculating vertical deflection from sea surface slopes (SSS) and sea surface height (SSH) measurements will be described in Section 2. Section 3 and Section 4 will investigate the performance of HY2 observations during the GM, as well as the strategies for two-pass waveform retracking, low-pass filtering, eliminating outliers, and marine gravity modeling. Section 5 will evaluate the HY-2A derived model using marine data, including published models and shipborne measurements.

2. Marine Gravity Modeling from Altimeter Measurements

The process of calculating vertical deflection from SSH can effectively restrain the effect of radial orbit error and other long-wavelength corrections. The use of vertical deflection makes the crossover adjustment procedure unnecessary and gets comparable results [20]. The relationship between vertical deflections and gravity anomalies can be deduced by the Laplace equation without any complex kernel function involved. The method is briefly summarized in Section 2.

Marine Gravity Anomaly by Vertical Deflection

Reference to a previous marine gravity recovery research [6], we also assume a flat earth approximation. The disturbing potential $T$ satisfies the Laplace equation in the given local planar coordinate system $(x,y,z)$. The gravity disturbances $\delta g$ can be expressed as the radial gradient of $T$ under such planar approximation. The two directional components of vertical deflections ($\xi$ and $\eta$) can be described as the horizontal gradients of $T$ according to Bruns formula. The relationship between gravity and vertical deflection [6] can be established as,

$$\frac{\partial \delta g}{\partial z}=-{\gamma }_0\left(\frac{\partial \xi }{\partial x}+\frac{\partial \eta }{\partial y}\right)$$

(1)

where ${\gamma }_0$ is normal gravity. Noting $\delta g(K,z)$, $\xi (K)$, and $\eta (K)$ represents the 2-dimensional Fourier Transform of $\delta g(x,y,z)$, $\xi (x,y)$, and $\eta (x,y)$, where the wavenumber $K$ satisfies the following conditions: $K=(k_x,k_y), k_x=1/{\lambda }_x, k_y=1/{\lambda }_y, \left|K\right|=\sqrt{k_x^2+k_y^2}$. Here, $k_x$ and $k_y$ are wavenumbers in the $(x,y)$ directions, while ${\lambda }_x$ and ${\lambda }_y$ represent corresponding wavelengths. According to the principles of Fourier transformation for the derivatives, the relationship between gravity and vertical deflection (Equation (1)) can be converted to the following equation.

$$\frac{\partial \delta g(K,z)}{\partial z}=-i2\mathsf{\pi }{\gamma }_0\left[k_x\xi (K)+k_y\eta (K)\right],i=\sqrt{-1}$$

(2)

The relationship between gravity disturbance on the ground and at elevation $z$ can be established based on the solution of Laplace equation in the wavenumber domain [6]. When the $z$ value is set to zero, the derivative of gravity along the radial direction can be written as Equation (3).

$$\frac{\partial \delta g(K,z)}{\partial z}|{}_{z=0}=-2\mathsf{\pi }\left|K\right|\cdot \delta g(K,0)e^{-2\mathsf{\pi }\left|K\right|z}|{}_{z=0}=-2\mathsf{\pi }\left|K\right|\delta g(K,0)$$

(3)

The right parts of Equations (2) and (3) are theoretically equivalent. Therefore,

$$\delta g(K,0)=\frac{i}{\left|K\right|}{\gamma }_0\left[k_x\xi (K)+k_y\eta (K)\right]$$

(4)

By introducing the operator $F^{-1}$ as the inverse Fourier transform,

$$\delta g(x,y,0)=F^{-1}\left[\delta g(K,0)\right]=\frac{i}{\left|K\right|}{\gamma }_0\left[F^{-1}\left[k_x\xi (K)\right]+F^{-1}\left[k_y\eta (K)\right]\right]$$

(5)

That is an expression of Fourier transform in order to calculate gravity from two components of vertical deflection [6]. Taking the difference between gravity disturbance and gravity anomaly into account, gravity anomaly is calculated,

$$\Delta g(x,y)=\delta g(x,y)-\frac{2}{R}T(x,y)=\delta g(x,y)-2\frac{{\gamma }_0}{R}N(x,y)$$

(6)

where, $T$ is the disturbing potential, $R$ is the average radius of Earth, and $N$ is the geoid height, which can be provided by geo-potential models. The last term can be further simplified as $-0.3076N$, where the unit of $N$ is meter. Marine gravity anomalies are practically calculated by combining Equations (5) and (6) once the vertical deflections are computed from satellite altimetry measurements.

To calculate vertical deflection, we use the derived along-track gradient information between consecutive observations, divided by the on-orbit velocity for the satellite [6]. These gradients can be calculated by differentiating each profile with respect to time using a first difference formula [7]. The latitudinal and longitudinal velocity of a satellite at nadir points could be calculated by using Equations (7) and (8) under the assumption of a circular orbit with respect to a flat earth locally [6].

$${\dot{\phi }}_c({\phi }_c)={\omega }_s{(1-\frac{{\cos }^{2}I}{{\cos }^2{\phi }_c})}^{1/2}$$

(7)

$$\dot{\lambda }={\omega }_s\frac{\cos I}{{\cos }^2{\phi }_c}-{{\omega }^{\prime }}_e{{\omega }^{\prime }}_e={\omega }_e-{\omega }_n$$

(8)

where, ${\omega }_s$ is orbit frequency, $I$ is orbit inclination, ${\phi }_c$ is geocentric latitude, ${\omega }_e$ is earth rotation rate, and ${\omega }_n$ is orbit precession frequency.

3. Data Analysis and Waveform Retracking

3.1. HY-2A Data Analysis

Through China’s National Satellite Ocean Application Center (NSOAC, http://www.nsoas.gov.cn/) data distribution, we obtained 7 complete cycles of the HY-2A Sensor Data Records (SDR) in standard NETCDF format. The general information of the obtained HY-2A GM altimeter waveform data is listed in

Table 1. The general information of acquired HaiYang-2A (HY-2A) GM altimeter waveform data sets.

Cycle	Time	Number of Available Passes (*)	Available Passes over SCS
GM C0001	20160324-20160908	4332 (298)	271
GM C0002	20160908-20170223	4373 (257)	309
GM C0003	20170223-20170810	4360 (270)	301
GM C0004	20170810-20180125	4399 (231)	303
GM C0005	20180125-20180712	4386 (244)	309
GM C0006	20180712-20181227	4333 (297)	313
GM C0007	20181227-20190613	4278 (352)	294

* In parentheses, the number of unprocessed passes is given. A small proportion of data cannot be normally processed to generate standard data products.

.

The geographic distribution of the 7 cycles of HY-2A ground tracks is shown in Figure 1. The across track distance between adjacent passes drops sharply from ~180 km in the ERM phase to ~16 km during the GM phase. The ground tracks are repeated within 3.0 km range. The overall amount of unprocessed passes correspond to 6% which is slightly higher than for the exact repetitive mission. The resulted number of valid passes in the SCS area for each cycle is listed in the last column of Table 1.

The SSH is computed by correcting for path delay and geophysical environment effects [21]. Most corrections are supplied in the level 2 product. However, several corrections have to be treated individually. No available corrections for the dry tropospheric path delay, inverted barometer effect, and pole tide are provided in the HY-2A GM SDR products. The effect of geocentric pole tide height can be ignored for gravity field retrieval due to its small effect on SSS [22] The dry tropospheric correction and inverted barometer correction are recalculated using a 6-hour interval real-time pressure model distributed by the National Centers for Environmental Prediction (NCEP) [23,24]. We adopt the Goddard Ocean Tide (GOT) V4.7 model for computing the geocentric ocean tidal effect. Detailed information of applied corrections used to estimate SSH is listed in

Table 2. Corrections used to estimate the sea surface height (SSH) for HY-2A geodetic mission (* means default item in the Sensor Data Records (SDR)).

Corrections	Used Models or Methods
Dry troposphere correction	NCEP
Wet troposphere correction	ECMWF (*)
Ionospheric correction	GIM (*)
Sea state bias	NSOAS empirical solution (*)
Ocean tide	GOT V4.7
Solid earth tide	Cartwright and Tayler tables (*)
Pole tide	Not applied
Inverted barometer correction	NCEP

.

HY-2A has run for more than 8 years and its orbit has gone through several maneuvers. Initially, we investigate the accuracy of HY-2A GM measurements using crossover differences. Smaller crossover differences between SSHs measured from ascending and descending passes indicate better stability and accuracy, which may decay over time at the later stage of satellite operations. Hence, we compute the crossover differences of HY-2A data for each GM cycle. The statistic values are listed in

Table 3. Statistics of the SSH discrepancy between ascending passes and descending passes at crossover points over SCS for HY-2A GM data. NC_A and NC_B denote the number of crossover points after and before deleting outliers, respectively. Unit: m.

Mission	Cycle	No Time Limit				Measurements within 1 Day
		NC_A/NC_B	Min	Max	STD	NC_A/NC_B	Min	Max	STD
HY-2A GM	1	2156/2198	−0.884	0.971	0.232	65/67	−0.332	0.482	0.156
	2	4242/4300	−0.976	0.951	0.227	81/82	−0.976	0.267	0.185
	3	2820/2986	−0.992	0.994	0.311	65/72	−0.865	0.794	0.218
	4	3364/3396	−0.999	0.958	0.253	67/67	−0.228	0.258	0.104
	5	2948/3012	−0.973	0.975	0.270	69/70	−0.643	0.287	0.151
	6	4173/4229	−0.984	0.967	0.244	77/78	−0.436	0.867	0.187
	7	2834/2869	−0.998	0.968	0.248	68/68	−0.756	0.968	0.204
	1-7	155121/158253	−0.999	0.999	0.278	498/510	−0.976	0.968	0.175

. A small proportion (~1%) of anomaly points with large discrepancies was eliminated before the statistical evaluations. Although we considered the effects by ocean tide and dynamic atmospheric correction, the natural variability of ocean surface and errors in data would inevitably affect the crossover discrepancies. As a result, both crossover differences with and without time limits were considered during the statistical analysis.

Compared to crossover statistic values in earlier HY-2A studies [18,19] for the ERM, the latest GM measurements stay at a consistent accuracy level with respect to ERM observations. A noticeable decrease in standard deviation (STD) values is observed for measurements within 1 day compared to those with no time limit. This indicates that the time-varying effect of the ocean surface is still an important signal in observations.

3.2. HY-2A Two-Pass Waveform Retracking

Waveform retracking is an effective tool to enhance range precision and improving marine gravity recovery. Over the ocean surface, the typical shapes of HY-2A raw waveforms are well-modelled using the Brown type model with five parameters: arrival time, rise time, amplitude, square of the antenna mis-pointing angle, and thermal noise [25].

Considering the inherent relationship between the errors in retracked estimates of arrival time and rise time, the two-pass waveform retracker estimates arrival time with higher precision by fixing the rise time parameter in a two-step model-fitting procedure [12,26,27]. In the first step, the waveforms are fitted by the three-parameter Brown model (arrival time, rise time, and amplitude). The thermal noise and square of the antenna mis-pointing angle are treated as constant parameters but need suitable pre-determination. In the second step, the rise time parameter is smoothed along-track and being fixed, before retracking the waveforms again using a two-parameter Brown model (arrival time and amplitude) [18]. Either the former three parameters or the latter two parameters will be solved during a weighted least squaring procedure.

$$M(t)=\frac{A}{2}\left[1+erf(\frac{t-t_0}{\sqrt{2\sigma }})\right]\exp \left[-\alpha (t-t_0)\right]$$

(9)

$${\chi }^2={\displaystyle {\sum }_{i=1}^{Num}\left[\frac{P_i-M(t_i-t_0,{\sigma }_h,A)}{W_i}\right]}$$

(10)

$$W_i=\frac{P_i+P_0}{\sqrt{K}}$$

(11)

The power $M(t)$ in HY-2A waveforms can be modeled by Equation (9), where $A$ is the amplitude, $t_0$ is the arrival time, $t$ is time for each sample gate, $\sigma$ is the effective pulse length, and $\alpha$ is the trailing edge decay factor relate to the antenna mis-pointing angle [18]. The effective pulse length $\sigma$ can be calculated from ${\sigma }^2={\sigma }_h^2+{\sigma }_p^2$, where ${\sigma }_h$ is the broadening of the pulse due to the interaction of ocean surfaces with certain significant wave heights (SWH) ${\sigma }_h=SWH/2c$. Here $c$ is the speed of the radar signals. The ${\sigma }_p=0.513{\tau }_p$ is the STD of the length of the outgoing pulse where ${\tau }_p=1/B$ and $B$ represents the chirp bandwidth. The chi-square measure of the misfit between the returned power samples $P_i$ and the model estimate $M(t_i)$ is described as Equation (10), where $Num$ is the number of gates (time samples), and $W_i$ is the uncertainty assigned to each sample and given by Equation (11), where $K=96$ is the number of statistically independent return echoes that were averaged to produce waveform records, and $P_0$ is the power offset value reflecting the thermal noise level of the instrument. The parameters to be estimated are respectively $[A,t_0,{\sigma }_h]$ and $[A,t_0]$ during the two steps.

Our previous study [18] resulted in a set of optimal constant parameters $[\alpha ,P_0]$ for HY-2A ERM waveform records. The numerical test for HY-2A GM passes indicates that the pre-determined $\alpha$ value of 0.0105 remains optimal by considering the minimum misfit STD values between normalized waveform and modeled waveform. We also inherit the power offset value $P_0$ (5500) from the ERM and a threshold value (0.015) with respect to the cumulative sum of waveform power [18].

We indicate successfully fitted waveform records using a maximum of 10 iterations. Successful fitting is defined for HY-2A waveforms if the estimated amplitudes are within the range of 70–130% of maximum echo power, while the misfit STD should not exceed 0.8 for a normalized waveform. We also remove waveforms having SWH outside of the range 0–10 m to exclude observations over extreme sea state conditions.

To evaluate the improvement in HY-2A GM data accuracy by the two-pass waveform retracker, we performed a statistical analysis on the retracked SSH measurements. The noise was estimated as the standard deviation of the SSHs with respect to the mean difference from the EGM2008 geoid plus dynamic ocean topography model computed over a 1-s interval. The dynamic ocean topography model was estimated together with EGM2008 and complete to degree and order (d/o) 180 (denote as DOT2008A_n180). This involved the retracked 20 estimates corresponding to the HY-2A sampling rate. Figure 2 (top panel) shows HY-2A GM noise levels for a sample pass. The noise level increases with increasing SWH. The solid smoothed curves (Figure 2, bottom panel) are medians of these estimates in 0.5-m SWH bins. Similar to the HY-2A ERM data [18], we find that the noise level for two-step retracked data (blue) is 1.6 times lower than that for one-step retracked data (green) and retracked results from SDR (red).

3.3. HY-2A SSS Calculation

After the two-pass waveform retracking procedure, we obtained 20 Hz retracked observations. Early studies showed that the high-sampling rate measurements retain more short-wavelength gravity information than 1 Hz observations [28], and the shortest wavelength that can be resolved in conventional radar altimeter data is about 10 km [29]. The spacing of the altimeter at 20 Hz sampling rate is ~350 m and the 20 Hz measurements are noisy. We resample the data to ~5 Hz or about 1.4 km using an along-track low-pass filter before computing the along-track SSS. The weight function during this process is calculated by a low-pass finite impulse response (FIR) filter. We apply the Parks-McClellan optimal equi-ripple FIR filter designed with a 0.5 gain corresponding to the wavelength of 6.7 km [30]. The filter weights and its impulse response in the frequency domain are designed according to the time-space intervals of HY-2A measurements (Figure 3). The filter rolls off nearly linearly over a band from 10 km to 5 km (Figure 3, lower). The corresponding moving window involves a sequence of 99 HY-2A 20 Hz measurements, which corresponds to a spatial range of about 32 km.

Based on the 20 Hz retracked observations, along-track 5 Hz SSHs are firstly computed by a downsampling procedure with weight coefficients of low-pass filter involved. After that, the along-track 5 Hz SSS were computed, along with the along-track gradient values of EGM2008 for editing. Editing was performed using a threshold of the slope differences (~30 microradian) which was set at three times the STD between observations and model values to eliminate outliers. The histogram for a sample pass (C0001P4632) is shown in Figure 4. Approximately 0.88% of data was removed. However, the along-track difference procedure in computing SSS amplified high-frequency noise signal again. Hence, we used another Parks-McClellan low-pass filter similarly designed with a sequence of 30 along-track 5 Hz SSS dataset involved to obtain filtered SSS [30]. In case of missing data, we adopted linear interpolation to enforce equal spacing and an edge window of ~2 s to taper the ends of the profiles.

4. Regional Marine Gravity for the SCS

The 7 cycles of HY-2A GM is augmented with a multi-satellite altimeter dataset to derive the 1 min marine gravity. This multi-satellite altimeter dataset, derived from eight radar altimeter missions (Geosat 1985–1989, ERS-1 1991–1998, T/P 1992–2006, ERS-2 1995–2011, Jason-1 2001–2013, Envisat 2002–2012, CryoSat-2 2010-2015, SARAL/AltiKa ERM 2013-2015), is preprocessed to along-track 5 Hz SSS by two-pass waveform retracker and similar down-sampling [30]. Except for HY-2A GM, CryoSat-2, and SARAL/AltiKa, the preprocessed datasets for other missions are generated by Scripps Institution of Oceanography (SIO) [31]. The dataset information used in our study is listed in

Table 4. General information of multi-satellite altimeter dataset for constructing regional marine gravity field model over SCS.

Mission	No. of 5 Hz SSS	No. of Outliers	Eliminating Rate (%)	Stacked Track Used
HY-2A GM	5680168	76799	1.35	No
Geosat GM	3107788	49314	1.59	No
Geosat ERM	117686	559	0.47	Yes
ERS GM	3195774	53377	1.67	No
ERS ERM	238514	1204	0.50	Yes
Envisat	11396374	101255	0.89	No
SARAL ERM	2172606	23636	1.09	No
Jason-1/GM	2612280	31209	1.19	No
Jason-1/ERM	23025417	263610	1.14	No
T/P	63509	147	0.23	Yes
T/P Tandem	647969	13225	2.04	No
CryoSat-2/LRM	6616436	65418	0.99	No
CryoSat-2/SAR	1104892	27060	2.45	No

. In order to evaluate the impact of HY-2A, we compare with the result by Zhang [30] who used the identical multi-satellite altimeter dataset over the SCS area without the HY-2A’s GM.

The DOT2008A_n180 and EGM2008 geoid (d/o 2190) were added and generated a SSS grid. Then the grid was interpolated and subtracted from the filtered along-track SSS dataset. The along-track residual slopes were divided by corresponding velocity at the observation position to calculate residual vertical deflections (Equations (7) and (8)). All along-track residual vertical deflections from 7 GM cycles were treated separately and we did not try to compute an average track during the process. This is because the ground track distance varies up to 3 km. This is also suitable for other altimeter data marked as “No” for stacked track used in the last column of Table 4. Based on the along-track residual vertical deflection, we firstly calculate them into gridding points and also calculates STDs at given gridding points. The STD value is used to assign a different weight for each data observation. Then we combined gridding residual DOV from different data observation to calculate north and east components at given gridding points. This gridding procedure is executed by biharmonic spline approach with a tension parameter of 0.25 described in [28,32]. Then the residual gravity anomaly grid was inverted based on Equations (1) and (2). Finally, gravity anomaly values of the EGM2008 model were restored to derive the 1′×1′ marine gravity field model over the SCS area (Figure 5).

5. Evaluation and Discussion

We performed comparisons and statistics of misfit between the new result (marked as HY-2A/GM derived result) and published altimetric marine gravity models and marine gravity measurements.

5.1. Verification with Marine Gravity Models

A series of global marine gravity models has been released by SIO and from the Technical University of Denmark (DTU) with the latest versions being V27.1 and DTU17, respectively. The V23.1 and DTU13 were introduced as they adopted a similar multi-satellite altimeter dataset to Zhang [30]. The DTU10 was included to represent marine gravity fields derived solely from the first generation GM data Geosat and ERS-1. Finally, Zhu’s result [19] was also included, as they adopted HY-2A 1 Hz data using standard retracking. All the models provide values at 1′×1′ gridding points and the difference was found using the Generic Mapping Tools (GMT) [33].

The STDs of the difference between the result with and without HY-2A/GM are shown in

Table 5. Statistics of difference between gravity with and without HY-2A/GM and other altimetric marine gravity models, and earth gravity models as well. Unit: mGal.

Model	Zhang’s Result (No HY-2A/GM)				HY-2A/GM Derived Result
	Min	Max	Mean	STD	Min	Max	Mean	STD
DTU10	−75.95	76.00	0.02	5.03	−79.15	74.22	0.02	4.95
DTU13	−72.65	70.35	0.02	4.83	−74.80	67.95	0.01	4.75
DTU17	−74.05	63.58	0.01	4.63	−77.25	60.77	0.01	4.55
V23.1	−81.70	198.00	−0.04	1.57	−72.60	198.00	−0.05	1.43
V27.1	−81.00	96.40	−0.07	2.47	−79.40	96.40	−0.08	2.44
Zhu	−85.361	101.860	0.23	5.03	−76.26	101.86	0.23	4.93
EGM2008	−85.29	87.37	0.03	4.95	−76.19	84.37	0.02	4.86
XGM2019	−81.99	93.14	0.29	6.46	−72.88	93.14	0.28	6.41
GOCO06S	−184.81	229.18	0.07	25.04	−183.07	229.18	0.06	25.00

. Discrepancies with altimetric marine gravity models are in the range of 1.4~4.9 mGal. Clearly, the HY-2A/GM derived result shows best consistency with SIO models due to similar usage of along-track vertical deflections for computing gravity. The smallest STD (~1.4 mGal) is found for the V23.1 model. This is attributed to largely similar dataset and method. The HY-2A/GM derived result has a smaller STD of discrepancy with the most recent DTU17 model than its two predecessors, DTU10 and DTU13. HY-2A/GM derived result fits all other models superior to the results by Zhang not using the HY-2A/GM data. This validation implies that HY-2A GM data is capable of improving marine gravity. The comparisons with different earth gravity models also indicate a slight enhancement in consistency with HY-2A/GM measurements involved. In addition, the much larger STD values show that the satellite-based model (GOCO06S, d/o 300) has large omission errors and is not suitable as reference model in our remove-restore procedure. Meanwhile, we wanted to compare the XGM2019 (d/o 2190) and EGM2008 (d/o 2190) at gridding points to show which model is more consistent with altimetric derived models. Compared to EGM2008, the XGM2019 has larger STDs with our derived results.

The differences between HY-2A/GM derived result and V23.1, V27.1 at gridding points over the SCS area are shown in Figure 6. The HY-2A/GM derived result is in better agreement with V23.1, while the differences at coastal regions, are significantly larger due to complex environment and contaminated echoes. The V27.1 model adopts more observations than V23.1, including 30 more months of CryoSat-2, 21 months of SARAL/AltiKa and 12 months of Jason-2 GM data. These additional data resolves more short wavelength gravity features and leads to slightly larger differences when comparing with our HY-2A/GM derived result.

5.2. Validation with Shipboard Gravity Measurements

We compare with marine gravity anomaly measurements along a number of surveys. This comparison provides an independent assessment of model accuracy but is limited to areas where shipboard data are available. The National Geophysical Data Center (NGDC) has 14 survey lines in the SCS area. Their geographic distributions are shown in Figure 7 (left). These marine data are edited to eliminate systematic bias and outliers through comparison with EGM2008. This leads to removal of 3.17% of the data [30]. Figure 7 right panel shows a zoom-in on one of the 14 survey lines in this investigation.

We identified and computed the cross-over statistics for all the location where the survey lines intersect with themselves to get the internal accuracy of the marine gravity. The cross-over differences and locations were acquired using a sliding-window searching method and discrepancies were calculated by linear interpolation from two adjacent observations [34]. The right panel of Figure 7 displays the cross-over locations with red circles for survey line 08020071. The crossover statistics and the number of crossovers are listed in the last column of

Table 6. Validation of altimetric gravity model with NGDC marine gravity along 14 survey lines in the SCS area with cross-over information. Number of measurements (NB_M) and number of crossover points (NB_C) are given in parentheses. Unit: mGal.

Cruise (NB_M)	HY-2A/GM Derived Result	Zhang [30]	Zhu [19]	V23.1	V27.1	DTU10	DTU17	EGM2008 (Reference)	Crossover Discrepancy (NB_C)
01010258 (33861)	8.46	8.47	8.90	8.51	8.40	8.97	8.72	8.94	3.97 (248)
01010259 (2505)	6.99	7.00	7.19	7.21	6.87	7.09	6.98	7.08	9.59 (46)
01010260 (28028)	8.00	8.03	8.11	8.02	7.97	8.12	7.78	8.10	4.22 (131)
08020049 (18068)	10.91	11.05	12.82	11.66	10.98	13.14	11.62	13.49	8.97 (109)
08020071 (45474)	6.91	6.96	6.48	6.87	6.79	6.47	5.98	6.53	6.72 (263)
08020072 (13068)	10.56	10.61	13.43	10.57	10.80	13.67	11.33	13.90	4.50 (68)
67010140 (9532)	5.38	5.39	5.33	5.45	5.27	5.23	5.12	5.29	9.04 (315)
67010141 (23843)	6.19	6.24	6.50	6.24	6.12	6.40	6.16	6.49	4.61 (75)
67010142 (15985)	9.28	9.35	9.47	9.69	8.81	9.55	8.63	9.72	7.60 (171)
67010146 (21896)	6.73	6.73	6.47	6.70	6.54	6.49	6.22	6.51	7.38 (139)
j1020003 (1860)	7.65	7.69	6.78	7.74	7.49	6.84	6.57	6.76	5.85 (86)
j1020026 (1132)	7.04	7.10	6.96	6.99	6.76	7.31	6.78	7.30	4.99 (178)
j1030028 (4206)	7.42	8.00	6.82	6.94	6.22	7.14	5.94	7.22	5.42 (248)
j2010003 (1151)	8.26	8.33	7.62	8.06	7.92	7.57	7.27	7.62	1.07 (92)
Average	7.86	7.92	8.18	7.95	7.75	8.23	7.63	8.31	6.17
Accuracy	4.88	4.96	5.38	5.02	4.69	5.45	4.49	5.57	--

. They range between 1.07 mGal (j2010003) and 9.59 mGal (01010259).

Comparisons between each model and marine gravity measurements are listed in Table 6. The averaged values are computed by weighting the STD for each cruise with the number of measurements involved in the verification (marked as ‘Average’ in Table 6). A smaller STD value generally represents better consistency between marine gravity result and shipborne measurements. More than 92% of cruises (13 out of 14) indicate that the HY-2A/GM derived result (column 2) has smaller STD values than Zhang’s result (column 3). This highlights the importance of adopting HY-2A/GM dataset since it is the only difference between the two results. The comparison between V23.1 and V27.1, as well as DTU10 and DTU17, prove a noticeable enhancement by using more geodetic altimeter data as done by SIO and DTU, separately. The averaged STD values show that the HY-2A/GM derived result has a slightly worse performance than DTU17 and V27.1, and a better performance than others. Among these statistics results from 14 selected cruises, 8 and 4 cruises have best agreement with DTU17 and V27.1, respectively. The smallest STD for each cruise is highlighted by bold font in Table 6. By contrast, the HY-2A/GM derived result has 2 smallest STD, indicating better performance than Zhu’s result (0), V23.1 (0), DTU10 (0), and Zhang’s result (0).

If we consider the crossover verification results as the synthesis of uncertainty of shipborne data, we computed a general accuracy estimate ${\Delta }_{result}$ for altimetric gravity according to the following equation.

$${\Delta }_{result}=\sqrt{{\sigma }_{diff}^2-{\sigma }_{ship}^2}$$

(12)

where ${\sigma }_{diff}$ is average STD with altimetric gravity models and ${\sigma }_{ship}$ is ship marine crossover discrepancy. Consequently, we deduced an accuracy of around 4.5~5.5 mGal for these altimetric gravity models. The specific values (marked as ‘Accuracy’) are listed in the last row in Table 6. In many cases, these are better than marine gravity.

6. Conclusions

For the SCS, satellite altimetry provides the most comprehensive images of the marine gravity field with accuracies better than available shipboard gravity data. Two-pass waveform retracker, which was known for optimal gravity recovery, was adopted for HY-2A. This paper attempted to evaluate the improved marine gravity field recovery over the SCS area (105°~125°E, 4°~26°N) by incorporating HY-2A GM retracked waveform data.

Little or no degradation between ERM and GM was found. Pre-determined parameters from the ERM mission remain optimal for the GM mission while the editing criteria need a slight adjustment. The retracked 20 Hz measurements were downsampled to 5 Hz by using a low-pass FIR filter to reduce the noise level. Residual vertical deflections were obtained through adding range and geophysical corrections, differencing, eliminating outliers, low-pass filtering, and removing sea surface topography and geoid height. Then we incorporated multi-satellite altimeter dataset of filtered SSS from the Geosat, ERS-1, Envisat, Jason-1, CryoSat-2, SARAL/AltiKa and HY-2A. Finally, regional marine gravity anomalies at 1′×1′ resolution were estimated based on the remove/restore procedure.

Two approaches were used to evaluate the accuracy of the new gravity model. We compared with global gravity models and subsequently with 14 NGDC marine gravity cruises. Our findings imply that the HY-2A/GM derived result is consistent with other published models. When HY-2A/GM measurements were incorporated, better consistency with shipborne gravity data was achieved for 13 out of 14 survey lines. The HY-2A/GM derived result is still slightly less accurate than the most recent global marine gravity fields (V27.1 and DTU17), demonstrating the importance of incorporating the latest CryoSat-2, SARAL/AltiKa, and Jason-2 GM data into the solution. Taking into account the uncertainty of marine gravity data from the averaged STD values by calculating squared difference, we got an accuracy of around 4.5~5.5 mGal for these altimetric gravity models.

In the near future, more GM observations from satellites with diverse track orientations and dense spatial distributions play complementary roles in calculating vertical deflections and provide a great opportunity to explore the finer structure of marine gravity fields. HY-2B, a follow-on mission to HY-2A, was successfully launched on 24 October 2018 and providing the continuation global and regional SSHs estimates. The third flight unit of the HY-2 program (HY-2C) is planned to launch in the near future. The potential altimeter satellite network can provide more measurements for marine gravity recovery. The next big step in gravity field improvement may come with the higher range precision and finer resolution provided by the Surface Water and Ocean Topography (SWOT) mission, which is scheduled for launch in 2021 [35].