In marine seismic acquisition, surface-related multiples constitute a significant portion of the acquired data. Typically, multiples are removed during early-stage data processing because they can lead to phantom reflectors during migration that might result in erroneous geologic interpretations. However, if properly dealt with, multiples can provide valuable extra information and can complement primaries in illuminating the subsurface. Reverse time migration has limitations in imaging multiples. A computationally efficient inversion procedure is proposed which jointly maps primaries and multiples to the true reflectors and estimates the source function on the fly. As a result, high-quality, mostly artifact-free broad-band images are obtained in which the imprint of the source function is partly removed at a computationally affordable expense compared with the combined costs of wave-equation-based surface-related multiple elimination and reverse time migration. All this is achieved by including the total downgoing wavefields as areal sources in least-squares migration in combination with curvelet-domain sparsity promotion. The proposed method is applied to a shallow-water marine data set from the North Sea which contains abundant short-period surface-related multiples, and the efficacy of the method is shown in eliminating coherent imaging artifacts associated with multiples. The benefits of joint imaging of primaries and multiples are demonstrated compared with imaging these signal components separately.