The invention discloses a digitalized simulation method of any three-dimensional (3D) rough surface, which improves simulation efficiency and simulation accuracy of statistical parameters. Random initial phase angle sequences are obtained through inverse Fourier transform of white noise sequences, and novel white noise and Fourier transform are obtained by utilizing the initial phase angle sequences. Processing including dispersion and Fourier transform and the like is performed through appointed autocorrelation functions to obtain power spectral density and transfer functions of Gaussian rough surface height sequences. The simulation of the Gaussian rough surface height sequences is finished by utilizing frequency domain dot product and a method for obtaining the inverse Fourier transform. On the basis, non-Gaussian rough surfaces are generated through height distribution statistical characteristic parameters including appointed skewness degree, kurtosis and the like and by utilizing combined Pearson and Johnson non-Gaussian transition systems. If the simulation accuracy of the skewness degree and the kurtosis is below standard, the phase angle sequences and the Fourier transform of the white noise are updated to re-perform Gaussian filtering and non-Gaussian transition till the given accuracy requirement is met.