<p>The stochastic variability of nonlinear dynamic properties of soil (such as shear wave velocity <i>V</i>_<i>s</i>, dynamic shear modulus ratio <i>G/G</i>_<i>max</i>, and damping ratio <i>λ</i>) is a critical factor influencing site seismic response characteristics and their associated uncertainties. This study takes a typical alluvial site in Suzhou as an example. Based on limited measured <i>G/G</i>_<i>max</i> and <i>λ</i> data, the Bootstrap resampling technique is employed to generate numerous parameter subsets from the small-sample data. Combined with the inverse Nataf transformation method, this approach constructs stochastic <i>G</i>/<i>G</i>_max and <i>λ</i> curves that vary with depth and exhibit negative correlation characteristics. Simultaneously, a non-stationary random field is used to simulate the spatial variability of the site’s shear wave velocity structure. A refined two-dimensional finite element model is established. The soil dynamic property parameters with representative variability levels are used as input parameters for site seismic response analysis. Based on a limited number of site response calculations, the Bootstrap method is applied to process the corresponding results, thereby revealing the uncertainty characteristics of the seismic response in the study area under the combined influence of bedrock input ground motion characteristics and the variability of soil dynamic parameters.</p>