Abstract: The sub-surface shear wave velocity structure of regional sites provides an important basis for the comprehensive consideration of the site seismic effects in the fine simulation of urban ground motion fields and seismic damages of building groups. In this study, a method is proposed to establish the complex staggered sedimentary structure of regional sites. First, the boundary of the staggered soils is confined based on the available information of engineering geology, then the spatial probability distribution of interlaced soils is solved by the sequential indicator simulation method, and then the velocity structure is modeled based on the spatial correlation between shear wave velocity and engineering geological structure quantified by the sequential Gaussian simulation method. Accordingly, the three-dimensional (3D) velocity structure of Harbin urban site is modeled based on the borehole soil information and shear wave velocity data. The feasibility of this method is verified through the qualitative analysis of consistence with engineering geological sections and the quantitative comparison with velocity structures tested in boreholes and estimated by the experimental relation between velocity and depth. The results demonstrate that the 3D model makes a good spatial prediction of engineering geological structure and reflects well on the structure characteristics of regional sites. The shear wave velocity structure of the proposed model shows a good correlation with the engineering geological structure, and even the velocity decreases with the increasing depth. Furthermore, in seismic response analysis, the proposed model is far more conducive to specifying the nonlinear dynamic parameters according to the detailed rock and soil types.