The Villain-Lai-Das Sarma (VLDS) equation has achieved significant attention in surface growth dynamics due to its effective description of Molecular Beam Epitaxy (MBE) growth processes. However, the scaling exponents of the VLDS equation driven by long-range correlated noise remain unclear, as different analytical approximation methods have yielded inconsistent results. The nonlinear term in the VLDS equation presents a challenge for numerical simulation methods, often leading to issues of numerical divergence. Existing numerical approaches primarily employ exponential decay techniques to replace nonlinear terms to alleviate the numerical divergence. However, recent studies have shown that these methods may change the scaling exponents and universality class of the growth system. Therefore, we propose a novel deep neural network-based method to address this issue. First, we construct a fully convolutional neural network to characterize the deterministic terms in the VLDS equation. To train the neural network, we generate training data with traditional finite-difference method before numerical divergence occurs. Then, we train the neural network to represent the deterministic terms, and perform simulations of VLDS driven by long-range temporally and spatially correlated noises based on the neural networks. The simulation results demonstrate that the deep neural networks constructed here possess good numerical stability. It can obtain reliable scaling exponents for the VLDS equation driven by different uncorrelated and correlated noises. Furthermore, this work also discovers that the VLDS system with long-range temporal correlation exhibits mound-shaped morphologies when the temporal correlation exponent is large enough, while the growing surface driven by spatially correlated noise still remains self-affine fractal structure independence on the spatial correlation exponent.