Abstract: The incompressible Navier-Stokes equations play a pivotal role across diverse scientific disciplines including environmental science, biomedical engineering, and fluid mechanics. Given the profound scientific significance and practical implications associated with developing robust and efficient numerical methodologies for solving the incompressible Navier-Stokes equations. To address this need, this study proposes a multiple-distribution-function regularized lattice Boltzmann (MDF-RLB) model for simulating the incompressible Navier-Stokes equations. At the heart of this approach lies a transformation of the incompressible Navier-Stokes equations into a coupled convection-diffusion system and then apply the regularized lattice Boltzmann method to model this system. Specifically, an evolution equation for the distribution function is constructed for each convection-diffusion equation (CDE) within the system. Subsequently, through rigorous Chapman-Enskog analysis, we have conclusively demonstrated that the proposed model can accurately recover the incompressible Navier-Stokes equations based on the convection-diffusion system. Furthermore, this research has yielded important theoretical derivations, including formulations for direct computation of velocity and pressure fields through zeroth and first moments of the distribution functions, as well as local computational schemes based on first moments of non-equilibrium distribution functions for calculating velocity gradients, velocity divergence, strain rate tensor, shear stress, and vorticity. Finally, to comprehensively validate the effectiveness and accuracy of the proposed model and the local computational formulations for these physical quantities, we conducted a series of benchmark numerical simulations, including two-dimensional Poiseuille flow, a simplified two-dimensional four-roll mill problem, and two-dimensional lid-driven cavity flow. The numerical tests reveal that the MDF-RLB model and the local formulations for physical quantities have a second-order convergence rate in space. Furthermore, compared to the multiple-distribution multiple-relaxation-time lattice Boltzmann (MDF-MRTLB) model, the MDF-RLB model demonstrates higher accuracy in certain cases and exhibits superior computational efficiency, achieving a reduction in computational time exceeding 7%.