Abstract: The frequencies estimation of two-dimensional multiple harmonics in additive noise is considered in this paper, by a robust iterative algorithm on the basis of a two-stage estimation. It can be proved that the algorithm needs only three iterative steps to converge and the final estimator attains the same convergence rate as Least Squares Estimator (LSE)，while the variance of the estimator attains the Cramer-Rao low bound. Since the statistics based estimation makes full use of the inner character of the harmonic model and the noise distribution, thus only three iterations are needed for the algorithm to work, while the robustness and efficiency in computation is also retained. Moreover, the estimator after three iterations is proved to be unbiased and consistent. We prove that the algorithm in 9 for the mono harmonic model can be generalized to the condition of multiple harmonics. Moreover, the parallel strategy for the frequency pair estimation can avoid the influence of the former frequency pair estimated on the later frequency pair to be estimated. Finally, the unbiasedness and consistency are verified via some simulation results, as well as for the high-accuracy under the condition of middle sample size.