In this work, a novel ring polymer representation for a multi-level quantum system is proposed for thermal average calculations. The proposed representation keeps the discreteness of the electronic states: besides position and momentum, each bead in the ring polymer is also characterized by a surface index indicating the electronic energy surface. A path integral molecular dynamics with surface hopping (PIMD-SH) dynamics is also developed to sample the equilibrium distribution of the ring polymer configurational space. The PIMD-SH sampling method is validated theoretically and by numerical examples.