Thermal expansion of GPS monuments and nearby bedrock could result in vertical changes in the coordinate time series of GNSS reference stations. In this paper, an improved method was developed to compute the magnitude of vertical variations caused by thermal expansion. Firstly, we calculated the effect on GPS monument and bedrock caused by thermal expansion based on land surface temperature data of GNSS reference stations and thermal expansion model. Secondly, we estimated the circular frequencies, amplitudes and phases using the method of least squares fitting instead of the current method which estimated only the amplitudes and phases information. Finally, we studied the periodic characteristics of the vertical variations caused by our modified thermal expansion model. Through analyzing the results of 9 representative IGS stations, we concluded that thermal expansion of GPS monuments and nearby bedrock could result in vertical variations of GNSS stations. The maximum variations could reach up to 0.57 mm and 1.85 mm at these stations respectively. The vertical variation caused by thermal expansion exhibited both annual and semiannual characteristics, which could explain 11.2% and 3.3% of the total annual and semi-annual variations in the up component of the coordinate time series respectively, and the magnitudes became larger with the increasing of their latitudes. Meanwhile, the amplitudes of the annual variations were much larger than that of the semi-annual variations. Meanwhile, some other small period (about 51 days) was also detected at some of these stations. In addition, we chose 107 IGS reference stations and computed the annual amplitudes and phases caused by thermal expansion of all these stations based on the method aforesaid. The results show that the maximum annual amplitude can reach to 3.3 mm, and their magnitudes show positive correlation with their latitudes prominently.