This contribution aims to enrich the recently introduced kernel-based regularization method for linear system identification. Instead of a single kernel, we use multiple kernels, which can be instances of any existing kernels for the impulse response estimation of linear systems. We also introduce a new class of kernels constructed based on output error (OE) model estimates. In this way, a more flexible and richer representation of the kernel is obtained. Due to this representation the associated hyper-parameter estimation problem has two good features. First, it is a difference of convex functions programming (DCP) problem. While it is still nonconvex, it can be transformed into a sequence of convex optimization problems with majorization minimization (MM) algorithms and a local minima can thus be found iteratively. Second, it leads to sparse hyper-parameters and thus sparse multiple kernels. This feature shows the kernel-based regularization method with multiple kernels has the potential to tackle various problems of finding sparse solutions in linear system identification.