Virtual-source surface wave responses can be retrieved using the crosscorrelation (CC) of wavefields observed at two receivers. Higher mode surface waves cannot be properly retrieved when there is a lack of subsurface sources that excite these wavefields, as is often the case. In this paper, we present a multidimensional-deconvolution (MDD) scheme that is based on an approximate convolution theorem. The scheme introduces an additional processing step in which the CC result is deconvolved by a so-called point-spread tensor. The involved point-spread functions capture the imprint of the lack of subsurface sources and possible anelastic effects, and quantify the associated spatial and temporal smearing of the virtualsource components that leads to the poor surfacewave retrieval. The functions can be calculated from the same wavefields as used in the CC method. For a 2-D example that is representative of the envisaged applications, we show that the deconvolution partially corrects for the smearing. The retrieved virtual-source response only has some amplitude error in the ideal situation of having the depth of the required vertical array equal to the depth penetration of the surface waves. The error is due to ignored cross-mode terms in the approximate convolution theorem. Shorter arrays are also possible. In the limit case of only a single surface receiver, the retrieved virtual-source response is still more accurate than the CC result. The MDD scheme is valid for horizontally layered media that are laterally invariant, and includes exclusively multicomponent point-force responses (rather than their spatial derivatives) and multicomponent observations. The improved retrieval of multimode surface waves can facilitate dispersion analyses in shallow-subsurface inversion problems and monitoring, and surface wave removal algorithms.