Surface loading effects related to atmospheric, hydrological, non-tidal ocean, are one of the principal sources of the seasonal oscillations in GNSS time series, and it should be taken into account for improving GNSS accuracy. In this study, the daily vertical time series of 9 GNSS stations at Hong Kong was used to investigate the surface loading (sum of atmospheric loading, hydrological loading, non-tidal ocean loading (AHNL)) contributors of seasonal oscillations in GNSS observations. This paper reveals a correlation between the AHNL deformation and the GNSS vertical time series, with an average correlation coefficient of 0.5. The GNSS vertical time series and the corresponding AHNL deformation at all stations exhibit identical amplitudes and phases. The average root mean square (RMS) reduction is 15% at all stations after removing the AHNL deformation from the GNSS vertical time series, implying that AHNL may contribute to non-linear fluctuations in GNSS observations at Hong Kong. Furthermore, the independent component analysis (ICA) method was performed to extract periodic signals from the GNSS time series. ICA method can effectively separate the seasonal signals related to AHNL, the seasonal signals show a strong correlation with AHNL deformations, with Lin correlation coefficients above 0.6. Finally, we carried out cross wavelet transform (XWT) method to quantitatively express the annual phase relationship between GNSS vertical time series and AHNL deformation. The XWT result shows AHNL mainly contribute to the annual oscillation in GNSS observations.