We propose a series-expansion thermal tensor network (SETTN) approach for efficient simulations of quantum lattice models. This continuous-time SETTN method is based on the numerically exact Taylor series expansion of the equilibrium density operator e-$β$H (with H the total Hamiltonian and $β$ the imaginary time), and is thus Trotter-error free. We discover, through simulating XXZ spin chain and square-lattice quantum Ising models, that not only the Hamiltonian H, but also its powers Hn, can be efficiently expressed as matrix product operators, which enables us to calculate with high precision the equilibrium and dynamical properties of quantum lattice models at finite temperatures. Our SETTN method provides an alternative to conventional Trotter-Suzuki renormalization-group (RG) approaches, and achieves a very high standard of thermal RG simulations in terms of accuracy and flexibility.