Purpose: It has been challenging to recover QSM in the presence of phase errors, which could be caused by the noise or strong local susceptibility shifts in cases of brain hemorrhage and calcification. We propose a Bayesian formulation for QSM where a two-component Gaussian-mixture distribution is used to model the long-tailed noise (error) distribution, and design an approximate message passing (AMP) algorithm with automatic and adaptive parameter estimation. Theory: Wavelet coefficients of the susceptibility map follow the Laplace distribution. The measurement noise follows a two-component Gaussian-mixture distribution where the second Gaussian component models the noise outliers. The distribution parameters are treated as unknown variables and jointly recovered with the susceptibility using AMP. Methods: The proposed AMP with parameter estimation (AMP-PE) is compared with the state-of-the-art nonlinear L1-QSM and MEDI approaches that adopt the L1-norm and L2-norm data-fidelity terms respectively. The three approaches are tested on the Sim2Snr1 data from QSM challenge 2.0, the in vivo data from both healthy and hemorrhage scans. Results: On the simulated Sim2Snr1 dataset, AMP-PE achieved the lowest NRMSE and SSIM, MEDI achieved the lowest HFEN, and each approach also has its own strong suit when it comes to various local evaluation metrics. On the in vivo dataset, AMP-PE is better at preserving structural details and removing streaking artifacts than L1-QSM and MEDI. Conclusion: By leveraging a customized Gaussian-mixture noise prior, AMP-PE achieves better performance on the challenging QSM cases involving hemorrhage and calcification. It is equipped with built-in parameter estimation, which avoids subjective bias from the usual visual fine-tuning step of in vivo reconstruction.