Virtual Green's functions obtained by seismic interferometry (SI) can provide valuable reflectivity data that can complement tomographic inversion schemes. However, virtual reflections are affected by illumination irregularities that are typical of earthquake-induced wavefields recorded by the receiver array. As a consequence, irregular source distributions, scattering (in case of suboptimal illumination), and complex source mechanisms can significantly disturb the retrieval of Green's function approximations by conventional SI methods. We introduce SI by full-field multidimensional deconvolution (MDD) for elastodynamic wavefields as an alternative method to obtain body wave Green's functions under those typical circumstances. The advantage of this method compared to other MDD methods is that the kernel of its governing equation is exact. This alternative formulation of the kernel pertains to several advantages: the solution is less sensitive to artifacts and utilizes the free-surface multiples in the data to estimate primary reflections. Moreover, the point spread function of the full-field MDD method corrects more affectively for irregular illumination because it also addresses irregularities caused by scattering inside the medium. In order to compare full-field MDD to existing SI methods, we model synthetic earthquake recordings in a subduction zone setting using an elastodynamic finite-difference scheme with double couples of different orientations and peak frequencies. Our results show that SI by cross correlation suffers most under these circumstances. Higher-quality reflections are obtained by the MDD methods, of which full-field MDD involves the most stable inversion, and its results are least contaminated by artifacts.@en