This paper develops a regularized lattice Boltzmann method for efficiently simulating the flow of N-phase immiscible incompressible fluids based on the phase field model that satisfies conservation and compatibility. By designing auxiliary moments, this method can accurately recover the second-order Allen-Cahn equation and the modified momentum equation. The correctness and effectiveness of the developed N-phase regularized lattice Boltzmann method are validated through numerical simulations of three-phase liquid lens spreading and Kelvin-Helmholtz instability phenomena. Finally, numerical simulations and analyses of three-phase Rayleigh-Taylor instabilities (RTI) are conducted, focusing on the evolution of the phase interface within the Reynolds number range of 500 \leqslant Re \leqslant 20000 (particularly under high Reynolds number condition of Re = 20000 ). Quantitative analyses are performed on the amplitude variations of bubbles and spikes at the two interfaces, as well as the changes in dimensionless velocity. We find that as the Reynolds number increases, the phase interface curls up at multiple locations due to Kelvin-Helmholtz instability, making the fluid more prone to dispersion and fragmentation. This study also simulates the evolutionary processes of RTI under different interface perturbations. These results demonstrate that RTI first develops at the perturbed interface, with its subsequent evolution inducing instability at a secondary interface.