Accurate channel estimation in orthogonal time frequency space (OTFS) systems with massive multiple-input multiple-output (MIMO) configurations is challenging due to high-dimensional sparse representation (SR). Existing methods often face performance degradation and/or high computational complexity. To address these issues and exploit intricate channel sparsity structure, this letter first leverages a novel hybrid burst-sparsity prior to capture the burst/common sparse structure in the angle/delay domain, and then utilizes an independent variational Bayesian inference (VBI) factorization technique to efficiently solve the high-dimensional SR problem. Additionally, an angle/Doppler refinement approach is incorporated into the proposed method to automatically mitigate off-grid mismatches.