Based on the corrected sum rules and generalized virial identities, we derive an expression for all modes of excitation spectrum of interacting Bose gases at finite atom numbers in axially anisotropic potentials, in terms of the N-body ground state average. Using the variational Gaussian calculation for the ground-state wave function, its explicit analytic formulas are obtained. These results show clearly the dependence of excitation spectrum on the interaction strength parameter $p = \sqrt{2/\pi}({a_{\rm sc}}/{a_{\rm ho}})(N-1)$ and trap geometry parameter $ \lambda $ for the system with N = 1 through $N\rightarrow \infty $ . For $ \lambda = 0$ and 1 the dependences have simple and intuitive physical interpretations. We compare the low-lying excitation spectra with the existing numerical results and make quantitative predications for future experiments and numerical simulation for higher-lying excitation modes.