It　is　of　importance　to　model　the　input　excitations　considering　the　spatial　variation　and　non-Gaussianity　for　the　reliability　evaluation　of　long　span　structures.　The　existing　multivariate　stationary　non-Gaussian　simulation　methods　are　mainly　based　on　the　target　marginal　distribution　and　cross-spectrum　matrix.　Due　to　the　Cholesky　decomposition　of　cross-spectrum　matrix,　the　computational　cost　of　these　methods　increases　rapidly　with　the　increase　of　spatial　points.　Therefore,　this　paper　develops　a　novel　simulation　method　based　on　the　wavenumber-frequency　spectrum　and　unified　Hermite　polynomial　model　(UHPM).　Firstly,　UHPM　is　extended　to　the　homogenous　and　stationary　non-Gaussian　field.　Then,　a　complete　transformation　model　from　homogenous　and　stationary　non-Gaussian　auto-correlation　function　(ACF)　into　underlying　Gaussian　ACF　with　its　applicable　range　is　proposed.　Furthermore,　two　types　of　incompatibility　between　the　first　four　marginal　moments　and　wavenumber-frequency　spectrum　are　discussed,　and　the　corresponding　remedies　are　provided.　To　facilitate　the　calculation,　the　2D-Fast　Fourier　Transform　(FFT)　technique　is　embedded　in　the　WienerKhintchine　transformation　and　spectral　representation　method　(SRM).　Finally,　a　unified　simulation　framework　for　multivariate　stationary　non-Gaussian　process　based　on　its　relationship　to　the　homogenous　and　stationary　non-Gaussian　field　is　presented.　Two　numerical　examples,　involving　the　simulations　of　non-Gaussian　wind　velocities　and　ground　motion　accelerations,　are　investigated　to　verify　the　capabilities　of　the　proposed　method.