Subsurface structures generally exhibit strong heterogeneities at multiple spatial scales. In this study, two- and three-scale high-order models are developed to effectively predict nonlinear gas transport in heterogeneous porous media with multiscale configurations. The considered heterogeneous porous media are composed of the matrix and inclusions at the mesoscopic and microscopic scales, respectively. First, using the formal two-scale asymptotic analysis, the homogenized solutions, the two-scale first-order and higher-order solutions are derived, with the first-order and second-order cell functions defined at the mesoscopic cell. Second, by further expanding all the mesoscopic cell functions to the microscopic levels, the second-order expansions of the mesoscopic cell functions are established and the upscaled relationships for the permeability tensor from the microscale to the macroscale are developed accordingly. Finally, the three-scale low-order and high-order solutions are constructed by combining the multiscale expansions of the mesoscopic cell functions and the macro–meso two-scale solutions. Several representative cases are simulated to demonstrate the accuracy and reliability of the proposed multiscale solutions. The results show that the high-order solutions can perfectly capture the locally steep pressure fluctuations and non-equilibrium effects caused by the heterogeneities and large permeability contrast in porous media with two- or three-scale configurations. The strategies to obtain the multiscale high-order solutions follow the reverse thought process of the reiteration homogenization method, and can be easily extended to heterogeneous porous media with arbitrary multiple scales.