Flash Radiography inspections stand to gain from inversion to infer density distribution of object based on X-ray transmission image. It is indispensable to be able to reliably provide uncertainties associated with the inversions. Although many inversion algorithms have been devised, they often perform poorly due to either their sensitivity to regularization parameter chosen in variational optimization or prohibitive computation and noisy results in stochastic simulation. In this paper, we present a gradual reconstruction algorithm, called TLE-Gibbs (two-level efficient Gibbs sampling), for flash radiography. At its core, TLE-Gibbs is a stochastic approach based on efficient Gibbs sampling and reconstruction refinement. A two-level scheme is proposed that enables high-resolution image to be constrained with uncertainty estimation from high-level reconstruction. Furthermore, a splitting variant that increases flexibility and precision is considered in the two-level scheme. An efficient Markov chain Monte Carlo (MCMC) endowed with first-order truncated conjugate gradient (CG) optimizer is developed to achieve minimal cost per sample and to approximate the posterior distribution. Finally, we adopt an effective refinement method to remove noises remained in the sample meanwhile maintaining sharp edges. For performance evaluation, TLE-Gibbs is applied on both synthetic data in which the influence of system blur is specially investigated and real data, and comparison with state-of-the-art reconstruction methods demonstrates the superiority of the proposed method.