Induced polarization tomography offers the potential to better characterize the subsurface structures by considering spectral content from the data acquisition over a broad frequency range. Spectral induced polarization tomography is generally defined as a non-linear inverse problem commonly solved through deterministic gradient-based methods. To this end, the spectral parameters, i.e., DC resistivity, chargeability, relaxation time, and frequency exponent, are resolved by individually or simultaneously inverting all frequency data followed by fitting a generalized Cole-Cole model to the inverted complex resistivities. Due to the high correlation between Cole-Cole model parameters and a lack of knowledge about the initial approximation of the spectral parameters, using the classical least-square methods may lead to inaccurate solutions and impede reliable uncertainty analysis. To cope with these limitations, we introduce a new approach based on a hybrid application of a globally convergent homotopic continuation method and Bayesian inference to reconstruct the distribution of the subsurface spectral parameters. The homotopic optimization, owing to its fast and global convergence, is first implemented to invert multi-frequency spectral induced polarization datasets aimed at retrieving the complex-valued resistivity. Then, Bayesian inversion based on a Markov-chain Monte Carlo (McMC) sampling method along with a priori information including lower and upper bounds of the prior distributions is utilized to invert the complex resistivity for Cole-Cole model parameters. By applying the McMC inversion algorithm a full nonlinear uncertainty appraisal can be provided. We numerically evaluate the performance of the proposed method using synthetic and real data examples in the presence of topographical effects. Numerical results prove that the homotopic continuation method outperforms the classic smooth inversion algorithm in the sense of approximation accuracy and computational efficiency. we demonstrate that the proposed hybrid inversion strategy provides reliable representations of the main features and structure of the Earths subsurface in terms of the spectral parameters.