Abstract: The slip surface obtained by the limit equilibrium method is not reasonable when the spatial variability of soil parameters is large. The finite element limit analysis method, FELA in simple, based on the generalized variational principle is introduced. By using the mixed finite element technology, the linear stress triangular elements and linear velocity triangular elements are proposed. By use of the strength reduction method and the optimization algorithm, the strict upper bound and lower bound of the safety factor of slopes are calculated. The influences of spatial variability of soil shear strength parameters on slope stability are analyzed, and the results are compared with those of the limit equilibrium methods. The results show that compared to the limit equilibrium method, the proposed FELA can search the critical slip surfaces effectively under spatial variability of soil strengths and provide the strict upper and lower bounds of safety factors of slopes. But the limit equilibrium method will search for unreasonable slip surfaces, and often overestimate the risk of landslides. The spatial variability of soil strengths results in the significant change in the distribution of safety factors. But the single safety factor can not reflect the change. According to the distribution of safety factors, the design value of soil parameters is suggested.