The global navigation satellite system (GNSS) position time series provides essential data for geodynamic and geophysical studies. Interpolation of the GNSS position time series is necessary because missing data will produce inaccurate conclusions made from the studies. The spatio-temporal correlations between GNSS reference stations cannot be considered when using traditional interpolation methods. This paper examines the use of machine learning models to reflect the spatio-temporal correlation among GNSS reference stations. To form the machine learning problem, the time series to be interpolated are treated as output values, and the time series from the remaining GNSS reference stations are used as input data. Specifically, three machine learning algorithms (i.e., the gradient boosting decision tree (GBDT), eXtreme gradient boosting (XGBoost), and random forest (RF)) are utilized to perform interpolation with the time series data from five GNSS reference stations in North China. The results of the interpolation of discrete points indicate that the three machine learning models achieve similar interpolation precision in the Up component, which is 45% better than the traditional cubic spline interpolation precision. The results of the interpolation of continuous missing data indicate that seasonal oscillations caused by thermal expansion effects in summer significantly affect the interpolation precision. Meanwhile, we improved the interpolation precision of the three models by adding data from five stations which have high correlation with the initial five GNSS reference stations. The interpolated time series for the North, East, and Up (NEU) are examined by principal component analysis (PCA), and the results show that the GBDT and RF models perform interpolation better than the XGBoost model.