This study aims to establish an integrated theoretical framework of a time-varying power spectral density (PSD) estimation for non-stationary earthquake ground motion excitations. First, time-frequency analysis methods, including the short-time Fourier transform, complex Morlet wavelet, and generalized harmonic wavelet transforms, are adopted to estimate the time-varying PSD in the presence of multiple samples. Priestley's estimation method is also introduced when the number of samples is limited. Second, by regarding the uniformly modulated and general modulated non-stationary Kanai-Tajimi models as the target seismic spectrum, a comparative study and an evaluation of the accuracy and the convergence of different time-varying PSD estimation approaches are performed. The guidance suggestions of the estimation method and the principle of parameter selection are also proposed herein to facilitate engineering application. Finally, a group of samples from three directions, including EW, NS, and UD of the 45th earthquake recorded in Strong Motion Array in Taiwan, Phase I or SMART-I, are used to reveal the typical transient features of the time-varying PSD and the coherence function among different directions. The results show that the complex Morlet wavelet and the generalized harmonic wavelet transforms exhibit a better performance compared to the short-time Fourier transform. The Priestley method also has advantages to one available sample. The earthquake ground motions of SMART-I exhibit dual non-stationary properties of intensity and frequency, and the three directions show a weak coherence. The research conclusions could also extend the application of spectral estimation theory in engineering and provide references for a further study on the multi-dimension and multi-point non-stationary seismic random response of large-span structures.