In response to the insufficient accuracy of traditional methods for obtaining viscoelastic parameters of epoxy asphalt mixtures, based on dynamic creep response characteristics, a modified generalized second-order Kelvin model was utilized, and analytical expressions for strain response under various loading modes were constructed. By integrating the creep mechanism and the physical significance of viscoelastic parameters, the fitting method was analyzed from three perspectives: loading mode, selection of fitting period, and choice of fitting parameters. The viscoelastic parameters were further transformed into a Prony series for static penetration finite element numerical simulation. Through deviation rate analysis and residual sum of squares comparison, a fitting method incorporating the Levenberg-Marquardt and global optimization methods was ultimately established. The results indicated that the optimized fitting method significantly enhanced the fitting accuracy during the initial creep stage, with an average error of 4.4% between numerical simulation and experimental results, and a residual sum of squares reduced to 0.003. This work markedly improves the characterization precision of viscoelastic parameters for epoxy asphalt mixtures and offers a novel methodological approach for developing constitutive models of such materials.