Simulation of multivariate non-Gaussian stochastic processes is still a challenging task, where the efficiency and accuracy may not be balanced. In this paper, an efficient methodology is presented to address this challenge. This methodology makes use of a probability weighted moments-based Hermite polynomial model and a wavenumber–frequency spectrum density-based spectral representation approach. First, a probability weighted moments-based Hermite polynomial model is employed to establish a transformation relationship between Gaussian and non-Gaussian samples. Importantly, this model only requires the solution of a solvable system of equations. Second, the target non-Gaussian wavenumber–frequency spectrum density can be promptly transformed into the Gaussian one, where a fast root matching strategy and fast Fourier transform are utilized. Finally, the underlying Gaussian samples are generated by the spectral representation method, and the non-Gaussian samples are obtained via the proposed transformation model. The computational results from numerical examples demonstrate that the proposed methodology is considered to be an accurate, applicable, and efficient approach for simulating multivariate non-Gaussian stochastic processes.