<p>The Bayesian inversion method demonstrates significant potential for solving inverse problems, providing both point estimation and uncertainty quantification (UQ). However, Bayesian maximum a posteriori (MAP) estimation may become unstable when the observed data come from diverse distributions (e.g., solutions of stochastic partial differential equations (SPDEs)). Additionally, Monte Carlo sampling methods are computationally expensive. To address these challenges, we propose a novel two-stage optimization method based on optimal control theory and variational Bayesian methods. This method provides stable solutions to stochastic inverse problems and efficiently quantifies their uncertainty. In the first stage, we introduce a new weighting formulation to ensure the stability of the Bayesian MAP estimation. In the second stage, we derive the necessary condition for efficient uncertainty quantification by combining the weighting formulation with variational inference. Furthermore, we establish an error estimation theorem that relates the optimally estimated solution to the exact solution under different amounts of observed data. Finally, the efficiency of the proposed method is demonstrated through numerical examples.</p>