This article aims at estimating the multi-variable commensurate fractional order systems based on the subspace identification method. By means of the properties of modulating function and integration by parts, the fractional order derivatives of input and output signals are converted into differentiation of modulating function. A new method via the nuclear norm and the Frobenius norm is proposed to handle the influence of measurement noise, avoiding the use of instrumental variable. At last, a simple simulation example is provided to verify the effectiveness of the proposed method.