In recent years, multi-stable origami structures have garnered increasing attention for their applications in dynamic scenarios such as robotic arm motions, impact energy absorption, and spectrum gap regulation. Understanding the intricate working mechanisms and exploring the rich dynamics of these structures necessitate the development of dynamic models. However, existing dynamic modeling methods for origami structures are often cumbersome, and the resulting dynamic models often lack interpretability. To overcome these limitations, we propose a novel data-driven dynamic modeling approach based on B-spline Galerkin method and subset selection strategy. This approach directly captures the dynamics of multi-stable origami structures using measured data, eliminating the need for reliance on empirical or prior knowledge. To validate the effectiveness of the proposed approach, we first evaluate it on the Duffing system, which has explicit expressions, successfully reconstructing the governing equation. Subsequently, we apply this method to the dynamic modeling of the origami ball structure with tri-stability and the multi-cell stacked Miura-origami (SMO) structure with high dimensionality through simulation, showing favorable results. Finally, using experimental data collected from a bi-stable SMO structure prototype, we employ the proposed method to obtain a global model that can accurately predict different dynamic behaviors over a broad range of excitation frequencies, including intra-well periodic vibrations, inter-well periodic vibrations, and inter-well chaotic vibrations. Overall, our method showcases outstanding efficacy in formulating interpretable, manageable, and comprehensive dynamic models. It plays a pivotal role in delving into the intricate dynamics of multi-stable origami structures.