Gaussian Process (GP), as a probabilistic non-linear multi-variable regression model, has been widely used in nonparametric Bayesian framework for the data based modelling of complex processes. The noise dynamics in standard GP regression is assumed to follow a Gaussian distribution. In this setting, the point estimation of the model parameters can be obtained analytically using the maximum likelihood (ML) approach in a straight forward fashion. However, in practical scenarios, processes may have been corrupted by the outliers and other disturbances or have multiple modes of operation, resulting a non-Gaussian data likelihood. In this work, to model such scenarios, we propose to employ a mixture of two Gaussian distributions as the noise model to capture both regular noise and irregular noise, thereby enhancing the robustness of the regression model. Further, we present an Expectation Maximization (EM) algorithm-based approach to obtain the optimal parameters set of the proposed GP regression model. The predictive distribution can then be found according to the estimated hyperparameters from the EM algorithm. The efficacy and practicality of the proposed method are illustrated with two sets of synthetic data, a simulated example, as well as an industrial dataset. (C) 2019 Elsevier Ltd. All rights reserved.