The Global Navigation Satellite System (GNSS) can provide the daily position time series for the geodesy and geophysical studies. However, due to various unpredictable factors, such as receiver failure or bad observation conditions, missing data inevitably exist in GNSS position time series. Most traditional time series analysis methods require the time series should be completed. Therefore, filling the missing data is a valuable step before analyzing the GNSS time series. In this study, a new method named Iteration Empirical Mode Decomposition (Iteration EMD) is proposed to fill the missing data in GNSS position time series. The simulation experiments are performed by randomly removing different missing percentages of the synthetic time series, with the added different types noise. The results show that Iteration EMD approach performs well regardless of high or low missing percentage. When the missing percentage increases from 5 % to 30 % with a step of 5 %, all the Root Mean Square Errors (RMSE) and Mean Absolute Errors (MAE) of Iteration EMD are smaller than Interpolation EMD. The relative improvements at different percentages of Iteration EMD relative to Interpolation EMD are significant, especially for the high missing percentage. The real GNSS position time series of eight stations were selected to further evaluate the performance of Iteration EMD with an average missing percentage 8.15 %. Principal Component Analysis (PCA) was performed on the filled time series, which is used to assess the interpolation performance of Iteration EMD and Interpolation EMD. The results show that Iteration EMD can preserve variance 75.9 % with the first three Principal Components (PC), more than 66.5% of interpolation EMD. Therefore, we can conclude that Iteration EMD is an efficient interpolation method for GNSS position time series, which can make full use of available information in existing time series to fill the missing data.