The permeability of a saturated porous medium is an important parameter in the field of water resources and geotechnical engineering. The geometric characteristics of a porous medium are key factors in the prediction of its permeability. In this paper, particles of different shapes are constructed by Cellular Automata (CA) random growth model, and particles with different surface characteristics are constructed by the spherical harmonic function. Then, porous media of different porosities, shapes, surface features, and particle size distributions are generated on the basis of the constructed particles. Three-dimensional Lattice Boltzmann Method is used for the pore-scale simulation of the seepage flow in a porous medium. The numerical results show that the effects of the particle shape and surface characteristics on the permeability are too obvious to be ignored. Using strict univariate analysis, the sensitivity of the various factors to the permeability, ordered from large to small, is as follows: porosity > particle size distribution > particle surface > particle shape. Based on numerical studies, a modified Kozeny–Carman (KC) formula is proposed by considering all the geometrical factors. All the parameters (the Wadell sphericity Sw, Cox roundness Rc, coefficient of non-uniformity Cu, the curve coefficient of curvature Cc, and effective particle size d10) in it are easily obtained in engineering practice and the accuracy of the formula is verified. Although It has been proven that the KC formula is applicable to multi-dispersed spherical particles and non-spherical particles whose surfaces are not very rough, its applicability to rough particles is limited. The modified KC formula does not have this limitation; therefore, it has a wider scope of application than the conventional KC formula.