Three-dimensional viscoelastic forward modeling is a crucial tool for understanding seismic wave propagation in regions covered by thick loess, red clay and saturated sandy soils. It also serves as a fundamental component in the theoretical study of three-dimensional passive source surface wave detection. The fractional-order wave equation based on the Kjartansson constant-Q model is widely used for two-dimensional viscoelastic forward modeling, reverse time migration and full waveform inversion, owing to its ability to decouple amplitude attenuation from velocity dispersion. However, its application to three-dimensional viscoelastic media faces significant challenges, including high memory requirements, substantial computational costs, and difficulties in accurately implementing the Viscoelastic Stress Image method (VSI) for free boundary conditions. To overcome these limitations, this study proposes a three-dimensional viscoelastic non-fractional wave equation forward modeling approach based on the Generalized Standard Linear Solid (GSLS) model integrated with VSI. The proposed method utilizes a high-order finite-difference scheme on a staggered grid, combined with Multiaxial Convolution Perfectly Matched Layer (MC-PML) absorbing boundary conditions. This approach efficiently yields numerical solutions to wave equations on structured grids, demonstrating strong practicality in exploration seismology. The method enables the successful implementation of VSI, providing high accuracy in simulating surface wave travel times and amplitudes, while maintaining numerical stability in media with high Poisson's ratios. The validity of the proposed method is verified through comparisons between numerical solutions and analytical solutions for both uniform half-space and layered models. Additionally, forward modeling tests on complex three-dimensional viscoelastic models confirm the method's adaptability to intricate geological conditions. The results indicate that the proposed method accurately captures phase distortion, amplitude attenuation, physical dispersion, wave conversion and energy redistribution in viscoelastic media. Furthermore, the approach demonstrates notable advantages in computational efficiency and memory usage, underscoring its potential for application in seismic exploration.