Remanent magnetization complicates the inversion of magnetic data due to altering the strength and direction of the magnetization vector. To deal with the problem, we have developed a 2D sequential inversion method of successively recovering the distributions of magnetization intensity and susceptibility and estimating the magnetization direction. The magnetization intensity distribution is first recovered from the total magnitude anomaly that is frequently transformed from the observed total-field anomaly and is invariant with the magnetization direction. With the recovered magnetization intensity distribution, we forward model the total-field anomalies caused by sources with different magnetization directions and calculate the correlations between the observed and predicted data. The orientation when the correlations attain a peak of maximum is defined as the optimal magnetization direction. Finally, the estimated magnetization direction helps to recover the susceptibility distribution by further inverting for total-field data. This method was tested by use of synthetic data and field data of two iron-ore deposits involving significant remanence, and all tests returned favorable results. The method obtaining the magnetization intensity and susceptibility distributions and an averaged magnetization direction made full use of the amplitude and phase information of magnetic anomalies, and it was more applicable for scenarios with a homogeneous magnetization direction.