We present a hp-inverse model to estimate a smooth, non-negative source function from a limited number of observations for a two-dimensional linear source inversion problem. A standard least-square inverse model is formulated by using a set of Gaussian radial basis functions (GRBF) on a rectangular mesh system with a uniform grid space. Here, the choice of the mesh system is modeled as a random variable and the generalized polynomial chaos (gPC) expansion is used to represent the random mesh system. It is shown that the convolution of gPC and GRBF provides hierarchical basis functions for the linear source inverse model with the hp-refinement capability. We propose a mixed l1 and l2 regularization to exploit the hierarchical nature of the basis functions to find a sparse solution. The hp-inverse model has an advantage over the standard least-square inverse model when the number of data is limited. It is shown that the hp-inverse model provides a good estimate of the source function even when the number of unknown parameters (m) is much larger the number of data (n), e.g., m∕n>40.