This paper presents a spectral representation method for simulation of non-stationary ground motion processes on the basis of Priestleys evolutionary spectral theory. Following this method, sample processes can be generated using a cosine series formula. It is shown that, these sample processes accurately reflect the prescribed characteristics of the evolutionary power spectral density function when the number of the terms in the cosine series is large enough; and the ensemble expected value and the ensemble autocorrelation function approach the corresponding target functions, respectively, as the sample size increases; and these sample processes are asymptotically normal as the number of the terms in the series tends to infinity. Finally, a few special cases of the formula are discussed, one of which is non-stationary white noise process, and other one is reduced to the formula for simulation of stationary stochastic processes.