Recently, frequency-based least-squares (LS) estimators have found wide application in identifying aircraft flutter parameters. However, the frequency methods are often known to suffer from numerical difficulties when identifying a continuous-time model, espe-cially, of broader frequency or higher order. In this article, a numerically robust LS estimator based on vector orthogonal polynomial is proposed to solve the numerical problem of multivariable systems and applied to the flutter testing. The key idea of this method is to represent the frequency response function (FRF) matrix by a right matrix fraction description (RMFD) model, and expand the numerator and denominator polynomial matrices on a vector onhogonal basis. As a result, a perfect numerical condition (numerical condition equals 1) can be obtained for linear LS estimator. Finally, this method is verified by flutter test of a wing model in a wind tunnel and real flight flutter test of an aircraft. The results are compared to those with notably LMS PolyMAX, which is not troubled by the numerical problem as it is established in z domain (e.g. derived from a discrete-time model). The verification has evidenced that this method, apart from overcoming the numerical problem, yields the results comparable to those acquired with LMS PolyMAX, or even considerably better at some frequency bands.