Global surrogate modeling aims to build a surrogate model with high accuracy covering the entire design domain. A major challenge to achieve this objective is how to reduce the number of function evaluations of the original computer simulation model. To date, the most widely used approach for global surrogate modeling is the adaptive surrogate modeling method. It starts with an initial surrogate model, which is then refined adaptively using the mean square error or by maximizing the minimum distance criteria. It is observed that current methods may not be able to effectively construct a global surrogate model when the underlying black box function is locally highly nonlinear. A new surrogate modeling method which can allocate more training points in regions with high nonlinearity is needed to overcome this challenge. This paper proposes an efficient global surrogate modeling method based on a multi-layer sampling scheme. An initial surrogate model is constructed first with a small group of training points. After that, samples of input variables are generated in multiple layers over the design domain. The generated samples are space-filling over the design domain and across layers. From the generated samples, new training points and candidate training points are identified in different layers by directly minimizing the prediction bias of the surrogate model. This enables us to improve the accuracy of the surrogate model more effectively than current available methods. Several mathematical examples and a non-linear vibratory system are used to demonstrate the effectiveness of the proposed global surrogate modeling method. The results show that the proposed method performs better than the most widely used variance minimization method.