Artificially determining the model order is necessary for time domain modal identification, which affects the identification efficiency and accuracy. An automatic identification method is proposed to determine the best estimation of the model order. Based on the covariance-driven stochastic subspace identification method (SSI-COV), the singular value decomposition of the Toeplitz matrix is carried out to determine the range of the lower triangular matrix. The similarity coefficient and distance function are introduced to cluster the modes, and the poles of the false modes are removed to obtain the clustering stabilization diagram. The model order is taken as the number of clustering centers of the fuzzy c-means (FCM) algorithm, and all the poles of the stabilization diagram are clustered to get the modal parameters of each order. The effectiveness of the proposed method is verified by adopting a thin-walled cylindrical column, a triangular-prism truss in simulation, and a flexible beam in experiment. The robustness of the method is illustrated by adding white noise to the simulated response data. The advantage of the clustering stabilization diagram is demonstrated by comparing with the traditional stabilization diagram.