Ground-penetrating radar (GPR) is a popular geophysical tool for mapping the underground. High-resolution 3D GPR data carry a large amount of information and can greatly help to interpret complex subsurface geometries. However, such data require a dense collection along closely spaced parallel survey lines, which is time consuming and costly. In many cases, for the sake of efficiency, a choice is made during 3D acquisitions to use a larger spacing between the profile lines, resulting in a dense measurement spacing along the lines but a much coarser one in the across-line direction. Simple interpolation methods are then commonly used to increase the sampling before interpretation, which can work well when the subsurface structures are already well sampled in the across-line direction but can distort such structures when this is not the case. In this work, we address the latter problem using a novel multiple-point geostatistical (MPS) simulation methodology. For a considered 3D GPR dataset with reduced sampling in the across-line direction, we attempt to reconstruct a more densely spaced, high-resolution dataset using a series of 2D conditional stochastic simulations in both the along-line and across-line directions. For these simulations, the existing profile data serve as training images from which complex spatial patterns are quantified and reproduced. To reduce discontinuities in the generated 3D spatial structures caused by independent 2D simulations, the target profile being simulated is chosen randomly, and simulations in the along-line and across-line directions are performed alternately. We show the successful application of our approach to 100 MHz synthetic and 200 MHz field GPR data under multiple decimation scenarios where survey lines are regularly deleted from a dense 3D reference dataset, and the corresponding reconstructions are compared with the original data.