The shield method is commonly employed in the construction of urban subways, and the stability of the excavation surface is extremely important for the safety of tunnel construction. When constructing tunnels in submerged or water-rich strata, the seepage force is a key factor affecting the ultimate support force of excavation faces. Uncertainty of natural geotechnical material hydraulic parameters significantly impacts the seepage force near the face. Based on the classic limit analysis and the random field theory, the principle of effective stress in saturated soil is introduced, and the stochastic numerical limit analysis method considering the influence of pore water pressure is established by combining finite element spatial discretization and second-order cone programming. On this basis, the influence of spatial variability of permeability coefficient on the stability of shield tunnel face in water-rich strata is studied. The results show that the heterogeneity of the permeability coefficient leads to an increase in the pore pressure gradient at the face. The increase in the variation of permeability coefficient and the ratio of vertical to horizontal permeability coefficient significantly increases the support force required to maintain the face stability, and the effect of cross correlation coefficient and autocorrelation distance is relatively small. The present research can provide theoretical guidance for evaluating the face stability of tunnels in water-rich strata.