Multivariate time series (MTS) clustering is an important technique for discovering co-evolving patterns and interpreting group characteristics in many areas including economics, bioinformatics, data science, etc. Although time series clustering has been widely studied in the past decades, no enough attention has been paid to capture time-varying correlation patterns in MTS. In this article, we propose a novel clustering approach for MTS data based on time-varying features. We introduce a time-varying Gaussian Markov Random Fields (T-GMRF) model to describe the correlation structure between MTS variables, and formulate the time-varying feature extraction problem as a convex optimization problem, which can be solved by a T-GMRF learning algorithm based on random block coordinate descent. We further apply a principal component analysis (PCA) based method on GMRF sequences to obtain low-dimensional feature vectors, and adopt a multi-density based clustering approach to form the cluster assignments. We conduct extensive experiments to compare the proposed T-GMRF method with 11 clustering algorithms based on 33 open MTS datasets, which show that T-GMRF significantly outperforms the state-of-the-arts with performance improvement up to 16%-64.5% on a variety of clustering performance metrics. The source codes of T-GMRF are publicly available at GitHub.