In this paper, we present a Krylov subspace based projection method for reduced-order modeling of large scale bilinear multi-input multi-output (MIMO) systems. The reduced-order bilinear system is constructed in such a way that it can match a desired number of moments of multi-variable transfer functions corresponding to the kernels of Volterra series representation of the original system. Numerical examples report the effectiveness of this method.