To identify multiple degree-of-freedom vibrating structures with local nonlinearities, a two-stage time domain approach based on the subspace method is proposed in this study. The locally nonlinear system to-be-identified is divided into an underlying linear part described by the FRF and a local nonlinear part described by nonlinear coefficients. The identification process of the proposed approach is not the same as the existing single-stage method. It identifies the underlying linear system before the local nonlinearities. The proposed approach reduces the dimensions of matrices which are used to calculate the system's state-space model and also gives a more appropriate estimation of the order of the underlying linear system with the utilization of classical spectrum estimation techniques. Both numerical and experimental examples are given to verify the performance of the method. Results show that the method is more accurate and reliable than the single-stage method, especially in a noisy environment.