The surrogate-based optimization of hierarchical stiffened composite shells against buckling is a typical multimodal and multivariables optimization problem. To improve the computational efficiency and global optimizing ability of the surrogate-based optimization of hierarchical stiffened composite shells, an enhanced variance reduction method based on Latinized partially stratified sampling and multifidelity analysis methods is proposed in this paper and then integrated into the surrogate-based optimization framework. In the offline step of the optimization framework, candidate pairing strategies of design variables are generated by Latinized partially stratified sampling and compared by performing priori optimizations based on the low-fidelity analysis method, and the optimal pairing strategy is therefore determined. On the basis of the optimal pairing strategy, the surrogate-based optimization is carried out using the highfidelity analysis method in the online step. With less computational cost in the offline step, the proposed enhanced variance reduction method overcomes the limitation of Latinized partially stratified sampling that the optimal pairing strategy is not obvious in complex problems. Then, extensive optimization examples are carried out to verify the efficiency and effectiveness of the proposed optimization framework. Given an approximate computational cost, the optimal buckling result of the proposed framework using enhanced variance reduction method increases by 18.2% than that of the traditional framework based on Latin hypercube sampling. In particular, the advantage of enhanced variance reduction method in the space-filling ability is highlighted in comparison to Latin hypercube sampling. When achieving an approximate global optimal solution, the proposed framework reduces the total computational cost by 76.3% than the traditional framework. Finally, the numerical implementation of asymptotic homogenization method is used herein for the accurate prediction of effective stiffness coefficients of the initial design and optimal results. Through comparison, it is concluded that the high axial stiffness and bending stiffness are the main mechanism for the high loadcarrying capacity of optimal results.