Picard iteration method is commonly used to obtain numerical solution of unsaturated flow in porous media. However, because the system of linear equations derived from Richards equation is seriously ill-conditioned, Picard iteration has slow convergence rate and low computational efficiency, particularly in layered porous media. In this study, control volume method based on non-uniform nodes is used to discrete Richards equation. To improve the convergence rate of Picard iteration, we combine the non-uniform multigrid correction method with the multistep preprocessing technology. Thus, an improved Picard iteration scheme with multistep preconditioner based on non-uniform multigrid correction method (NMG-MPPI(m)) is proposed to model 1D unsaturated flow in layered porous media. Three test cases were used to verify the proposed schemes. The result shows that the condition number of the coefficient matrix has been greatly reduced using the multistep preconditioner. Numerical results indicate that NMG-MPPI(m) can solve Richards equation at a faster convergence rate, with higher calculation accuracy and good robustness. Compared with conventional Picard iteration, NMG-MPPI(m) shows a very high speed-up ratio. As a result, the improved Picard iteration scheme has good application for simulating unsaturated flow in layered porous media.