<p>This study investigates the seismic response of multi-support structures, such as cable-stayed bridges, subjected to spatially and randomly varying ground motions, with a particular focus on uncertainties in soil properties, namely porosity, shear modulus, and layer depth. A refined Improved Descriptive Sampling (IDS) method, validated through 500 realizations, outperforms conventional Monte Carlo and Latin Hypercube techniques in representing soil uncertainty. A refined state space formulation combined with a ground motion simulation approach, benchmarked against Wavelet Packet Transform and Expansion Optimal Linear Estimation-based methods, effectively models non-stationary excitations. The adopted poroviscoelastic soil model incorporates fluid–solid interactions, revealing substantial reductions in the soil amplification function by 52.9%, 56.3%, and 85% for porosity, shear modulus, and layer depth (Coefficient of variation, equal to 30%, 20%, and 50%) along with frequency shifts of the soil layer of 6.7%, 3.5%, and 24.5%, respectively. Ground acceleration time histories at bridge supports indicate peak reductions from 3–3.5 to 2–2.5 m/s². The bridge’s displacement and acceleration responses in the horizontal direction (Degrees of Freedom (DOFs) 1 and 2, corresponding to the tops of the towers) and in the vertical direction (DOF 3, corresponding to the center of the bridge) exhibit reduced amplitudes, broader envelopes, and smoother oscillations, reflecting enhanced damping and period lengthening. These results underscore the importance of probabilistic modeling and soil randomness representation in improving the accuracy of seismic response analyses for multi-support structures founded on porous media.</p>