To improve the accuracy and robustness of modulation format identification (MFI) in elastic optical network (EON), this paper proposes an MFI method based on Stokes space and a Stacking model. The method extracts one-dimensional probability distribution features of the three axes in Stokes space using kernel density estimation to construct a 240-dimensional feature vector. A genetic algorithm is employed to optimize the combination of base models and meta-models in the Stacking model, while Bayesian optimization is used for hyperparameter tuning, enhancing classification performance under low signal-to-noise ratios. Simulation results show that, within an optical signal-to-noise ratio (OSNR) range of 5~30 dB, the model achieves a macro-average area under the receiver operating characteristic curve (AUC) of 1. The identification accuracy exceeds 98.5% for modulation formats such as polarization-division multiplexing binary phase-shift keying (PDM-BPSK) and polarization-division multiplexing quadrature phase-shift keying (PDM-QPSK), with an average accuracy improvement of 2.05%~5.63% compared to benchmark models like XGBoost and TabNet. Additionally, 100% identification precision is achieved at an OSNR of 18 dB.