Reverse time migration of controlled-order multiples (RTM-CM) can significantly remove strong crosstalk artifacts in multiples imaging. However, the computational cost is demanding if all pairs of consecutive-order multiples are considered, and an RTM-CM image retains residual crosstalk. To address these issues, we have developed a phase-encoding-based migration approach. This approach first phase-encodes different-order multiples by randomly modifying time shifts and polarity reversals, and then blends the encoded data into supergathers. As a result of using supergathers, migrations of different-order multiples are simultaneously accomplished, which considerably improves computational efficiency. As the simultaneous migration inherits the advantages of RTM-CM and provides good gradients for least-squares migration (LSM), we introduced an iterative LSM to gradually attenuate residual crosstalks in our approach. Numerical experiments on synthetic and field data sets revealed that our approach can significantly enhance imaging quality by suppressing crosstalks, increasing spatial resolutions, and reducing computational costs.