A three-dimensional rotational failure mechanism of embankment was built based on the upper limit analysis theorem. Drawing on the shear strength theory of Bishop unsaturated soil and considering the spatial distribution of matrix suction inside the unsaturated embankment and its change with the groundwater level, the energy conservation equation encompassing the external force work rate and the internal energy dissipation rate of the embankment was established. An efficient genetic algorithm was developed to search for the least upper-bound value of unsaturated soil embankment. By degrading the unsaturated soil embankment failure mode into a slope failure mode and comparing critical values with existing unsaturated soil slope stability calculations, the correctness of the proposed upper-bound values and the accuracy of the genetic search algorithm were verified. Furthermore, the influence of the pore size distribution, matrix suction, the inverse of the air entry value, embankment inclination angle and effective internal friction angle on the three-dimensional stability of the embankment was systematically studied. The results demonstrated that the stability of unsaturated soil embankment depends not only on the properties of soil but also on the pore size parameters and the inverse of the air entry value. Additionally, the fluctuation of matric suction resulting from changes in groundwater levels significantly impacts the stability of unsaturated soil embankments. These findings offer an important theoretical basis for granular analysis of embankment stability and the design of ultimate fill height, providing valuable insights for improving embankment stability analysis and facilitating the design of embankments.