Receiver function is important for imaging crustal and upper mantle discontinuities. However, sparsely scattered stations would introduce imaging artifacts or even lead to misinterpretation on complex structures. We regularize 3D teleseismic wavefield to reduce the artifacts using radial basis function interpolation. First, we evaluated the feasibility of the wavefield regularization with several typical models by synthetic data. The results demonstrate the high reliability of our method on recovering local 2D and 3D structures, even when seismic stations are intentionally missed 95% on a uniform fine grid. Then, we applied this method to sparsely deployed stations in Northeast China. The waveforms reconstructed from surrounding stations show good consistency with the observed waveforms; furthermore, the imaging results using the regularized data are highly comparable with the reference results obtained by a dense 2D seismic array of 60 stations (with a spacing of 10-17 km), even though our input data are mainly contributed by only 9 stations (with an averaged spacing of 80 km). Our results have better continuity of 3D topography of subsurface, compared with that obtained by traditional method. Our regularization method could significantly improve the spatial resolution of receiver function imaging both for sparse and dense distribution of seismic stations, especially for imaging relatively complex structures with lateral variations.