In this paper, the regularization approach introduced recently for nonparametric estimation of linear systems is extended to the estimation of nonlinear systems modelled as Volterra series. The kernels of order higher than one, representing higher dimensional impulse responses in the series, are considered to be realizations of multidimensional Gaussian processes. Based on this, prior information about the structure of the Volterra kernel is introduced via an appropriate penalization term in the least squares cost function. It is shown that the proposed method is able to deliver accurate estimates of the Volterra kernels even in the case of a small amount of data points.