This paper develops an efficient solution for superresolution two-dimensional (2D) harmonic retrieval from multiple measurement vectors (MMV). Given the sample covariance matrix constructed from the MMV, a gridless compressed sensing approach is proposed based on the atomic norm minimization (ANM). In the approach, our key step is to perform a redundancy reduction (RR) transformation that effectively reduces the large problem size at hand, without loss of useful frequency information. For uncorrelated sources, the transformed 2D covariance matrices in the RR domain retain a salient structure, which permits a sparse representation over a matrix-form atom set with decoupled 1D frequency components. Accordingly, the decoupled ANM (D-ANM) framework can be applied for superresolution 2D frequency estimation. Moreover, the resulting RRenabled D-ANM technique, termed RR-D-ANM, further allows an efficient relaxation under certain conditions, which leads to low computational complexity of the same order as the 1D case. Simulation results verify the advantages of our solutions over benchmark methods, in terms of higher computational efficiency and detectability for 2D harmonic retrieval.