Free-surface-related multiples can provide extra illumination of the subsurface and thus can be usefully included in migration procedures. However, most multiple migration approaches require separation of primaries and free-surface-related multiples or at least prediction of multiples in advance, which is time consuming and prone to errors. The data-to-data migration (DDM) method migrates free-surface-related multiples by forward and backward propagating the recorded full data (containing primaries and free-surface-related multiples). For DDM, there is no need to predict or separate multiples, but the migration results suffer from the crosstalk generated by crosscorrelations of undesired seismic events, e.g., primaries and second-order free-surface-related multiples. We have developed least-squares DDM (LSDDM) for marine data to eliminate the crosstalk generated by DDM. In each iteration, the forward-propagated primaries and free-surface-related multiples are crosscorrelated with the backward-propagated primary and free-surface-related multiple residuals to form the reflectivity gradient. We use a three-layer model and the Marmousi model for numerical tests. The results validate that LSDDM can provide a migrated image with higher signal-to-noise ratio and more balanced amplitudes than DDM. The LSDDM approach might be valuable for general subsurface imaging for marine seismic data when the migration velocity is accurate, and the acquired data have sufficient recording time.