Cone penetration test (CPT) is commonly used to determine the stratification of underground soil and the mechanical parameters of soils in stratification. Due to time, resources and/or technical constraints, the number of CPT soundings along with a horizontal direction is generally limited. In such cases, spatial interpolation or stochastic simulation methods is a necessary choice to estimate CPT data at un-sampled locations. This paper proposes an efficient method for simulating CPT data at un-sampled locations directly from a limited number of CPT records. The approach couples the framework of 2D Bayesian compressive sensing with Gibbs sampling, where Kronecker product is introduced for facilitating its simulation efficiency. Both numerical simulations and case histories are used to illustrate the presented method.Results show that the proposed method is reasonable, which can not only reflect the non-stationary characteristics of the data, but also significantly reduce the time cost and have reasonable adaptability after using the sequential updating technique. In addition, the accuracy and reliability of interpolation are negatively and positively proportional to the distance from existing CPT soundings and the number of existing CPT soundings, which demonstrates the data-driven nature of the proposed method.