Bilinear systems are a class of nonlinear dynamical systems that arise in a variety of applications. In order to obtain a sufficiently accurate representation of the underlying physical phenomenon, these models frequently have state-spaces of very large dimension, resulting in the need for model reduction. In this work, we introduce two new methods for the model reduction of bilinear systems in an interpolation framework. Our first approach is to construct reduced models that satisfy multipoint interpolation constraints defined on the Volterra kernels of the full model. We show that this approach can be used to develop an asymptotically optimal solution to the H_2 model reduction problem for bilinear systems. In our second approach, we construct a solution to a bilinear system realization problem posed in terms of constructing a bilinear realization whose kth-order transfer functions satisfy interpolation conditions in k complex variables. The solution to this realization problem can be used to construct a bilinear system realization directly from sampling data on the kth-order transfer functions, without requiring the formation of the realization matrices for the full bilinear system.