It is well known that conventional harmonic lattice dynamics cannot be applied to dynamically unstable crystals at 0 K, such as high temperature body centered cubic (BCC) phase of crystalline Zr. Predicting phonon spectra at finite temperature requires the calculation of force constants to the third, fourth and even higher orders, however, it remains challenging to determine to which order the Taylor expansion of the potential energy surface for different materials should be cut off.Molecular dynamics, on the other hand, intrinsically includes arbitrary orders of phonon anharmonicity, however, its accuracy is severely limited by the empirical potential field used.Recently, machine learning algorithms emerge as a promising tool to build accurate potentials for molecular dynamics simulation. In this work, we approach the problem of predicting phonon dispersion at finite temperature by performing molecular dynamics simulations with machine learning-driven potential fields. We developed Gaussian approximation potential models for both the hexagonal closed-packed (HCP) phase and the body centered cubic (BCC) phase of Zirconium crystals. The developed potential field is first validated with static properties including energyvolume relationship, elastic constants and phonon dispersions at 0 K. Molecular dynamics