The manta ray is a large marine species, which has the ability of gliding efficiently and flapping rapidly. It can autonomously switch between various motion modes, such as gliding, flapping, and group swimming, based on ocean currents and seabed conditions. To address the computational resource and time constraints of traditional numerical simulation methods in modeling the manta ray’s three-dimensional (3D) large-deformation flow field, this study proposes a novel generative artificial intelligence approach based on a denoising probabilistic diffusion model (surf-DDPM). This method predicts the surface flow field of the manta ray by inputting a set of motion parameter variables. Initially, we establish a numerical simulation method for the manta ray’s flapping mode by using the immersed boundary method and the spherical function gas kinetic scheme (IB-SGKS), generating an unsteady flow dataset comprising 180 sets under frequency conditions of 0.3–0.9 Hz and amplitude conditions of 0.1–0.6 body lengths. Data augmentation is then performed. Subsequently, a Markov chain for the noise diffusion process and a neural network model for the denoising generation process are constructed. A pretrained neural network embeds the motion parameters and diffusion time step labels into the flow field data, which are then fed into a U-Net for model training. Notably, a transformer network is incorporated into the U-Net architecture to enable the handling of long-sequence data. Finally, we examine the influence of neural network hyperparameters on model performance and visualize the predicted pressure and velocity fields for multi-flapping postures that were not included in the training set, followed by a quantitative analysis of prediction accuracy, uncertainty, and efficiency. The results demonstrate that the proposed model achieves fast and accurate predictions of the manta ray’s surface flow field, characterized by extensive high-dimensional upsampling. The minimum PSNR value and SSIM value of the predictions are 35.931 dB and 0.9524, respectively, with all data falling within the 95% prediction interval. Compared with CFD simulations, the single-condition simulations by using AI model show that the prediction efficiency is enhanced by 99.97%.