Current non-stationary load models based on the evolutionary power spectral density (EPSD) may lead to overestimation and ambiguity of structural responses. The quasi-stationary harmonizable process with its Wigner-Ville spectrum (WVS) and Loève spectrum, which do not suffer from the deficiencies of EPSD, is suitable for modeling non-stationary loads and analyzing their induced structural responses. In this study, the multi-taper S-transform (MTST) method for estimating WVS and Loève spectrum of multi-variate quasi-stationary harmonizable processes is presented. The analytical biases and variances of the WVS, Loève spectrum, and time-invariant and time-varying coherence estimators from the MTST method are provided under the assumption that the target multi-variate harmonizable process is Gaussian. Using a numerical case of a bivariate harmonizable wind speed process, the superiority and reliability of the MTST method are demonstrated through comparisons with several existing methods for the WVS and Loève spectrum estimations. Finally, the MTST method is applied to two pieces of ground motion acceleration records measured during the Turkey earthquake in 2023.