A sound understanding of subsurface geological conditions is crucial for the digitalization of underground infrastructure, as well as for smart city applications. The building and updating of underground digital twins heavily rely on sparse geotechnical measurements (e.g., boreholes) retrieved from the ground, and an efficient sampling strategy can facilitate the interpretation of subsurface heterogeneities. Geotechnical sampling design can be viewed as a constrained optimization process that aims to obtain as much geological information as possible from a limited number of sampling locations within a given site boundary. The conventional sampling strategy often specifies equal sampling spacing with regular patterns without considering subsurface stratigraphic variations and irregular site geometries. In this study, a data-driven intelligent sampling strategy is proposed to optimize borehole locations for a multi-stage site investigation of a three-dimensional (3D) geological domain. The initial sampling plan is determined using weighted centroidal Voronoi tessellation, which assigns uniform sampling densities to zones of different importance. Measurements obtained from the initial stage are combined with prior geological knowledge to build underground digital twins using an image-based stochastic modeling method known as IC-XGBoost3D. Multiple geological domains can be developed under the framework of Monte Carlo simulation, and stratigraphic uncertainties associated with multiple random geological domains can be quantified using information entropy. The location with the maximum entropy is adaptively selected as the next optimal sampling location. The proposed method is the first sampling strategy that can explicitly consider 3D subsurface stratigraphic variations. The performance of the proposed multi-stage sampling strategy is demonstrated using a simulation example. The results indicate that the proposed method can efficiently identify the optimal sampling locations while accounting for irregular site geometries and 3D subsurface stratigraphic uncertainties.