Abstract: To enhance the numerical accuracy and stability of the Lattice Boltzmann (LB) method in solving the incompressible Navier-Stokes (N-S) equations, this paper proposes a class of second-order finite-difference LB methods. The method employs a three-stage second-order scheme for discretizing the temporal derivatives, while the spatial derivatives are discretized using a hybrid scheme combining second-order central differences and second-order upwind schemes. The stability of such numerical methods is verified through von Neumann analysis, and the accuracy of the model is tested via numerical experiments. In the numerical experiments, the impact of the Courant–Friedrichs–Lewy (CFL) number on the computational results is quantitatively analyzed by simulating two-dimensional Taylor vortex flow and two-dimensional lid-driven cavity flow. The results indicate that the second scheme proposed in this paper exhibits the best numerical accuracy and stability. Furthermore, the use of the hybrid difference scheme effectively reduces numerical dissipation and oscillations, significantly improving numerical accuracy. The findings also demonstrate that the model achieves second-order convergence accuracy in both time and space, and outperforms the single-relaxation-time Lattice Boltzmann method (BGK-LBM) and the multiple-relaxation-time Lattice Boltzmann method (MRT-LBM) in terms of error reduction rate. The finite-difference Lattice Boltzmann method proposed in this study provides a new numerical approach for high-precision simulation of complex flow problems and can be further extended to the numerical solution of nonlinear convection-diffusion equations in the future.