Accurate channel estimation is critical for realizing the performance gains of massive multiple-input multiple-output (MIMO) systems. Traditional approaches to channel estimation typically assume ideal receiver hardware and linear signal models. However, practical receivers suffer from impairments such as nonlinearities in the low-noise amplifiers and quantization errors, which invalidate standard model assumptions and degrade the estimation accuracy. In this work, we propose a nonlinear channel estimation framework that models the distortion function arising from hardware impairments using Gaussian process (GP) regression while leveraging the inherent sparsity of massive MIMO channels. First, we form a GP-based surrogate of the distortion function, employing pseudo-inputs to reduce the computational complexity. Then, we integrate the GP-based surrogate of the distortion function into newly developed enhanced sparse Bayesian learning (SBL) methods, enabling distortion-aware sparse channel estimation. Specifically, we propose two nonlinear SBL methods based on distinct optimization objectives, each offering a different trade-off between estimation accuracy and computational complexity. Numerical results demonstrate significant gains over the Bussgang linear minimum mean squared error estimator and linear SBL, particularly under strong distortion and at high signal-to-noise ratio.