For effectively solving the stochastic groundwater flow problems in heterogeneous porous media, a reduced multiscale finite element basis method (RMsBM) was developed by constructing a set of reduced multiscale finite element basis functions independent of the parameter sampling and generating a reduced-order multiscale model. Moreover, in order to improve the online computation of reduced-order model, a matrix discrete empirical interpolation method (MDEIM) was used to affinely decompose the discrete system of the nonaffine parametrized problem. Finally, several numerical experiments were carried out for saturated flow problems in heterogeneous random porous media to evaluate the performance and efficiency of the proposed method. The numerical results showed that the proposed RMsBM can significantly enhance online computation efficiency and maintain comparative accuracy for the flow problems by choosing a suitable coarse mesh and an optimal number of local reduced multiscale basis functions.