This paper presents a model reduction method for the class of linear quantum stochastic systems often encountered in quantum optics and related fields. The approach is proposed on the basis of an interpolatory projection ensuring that specific input-output responses of the original and the reduced-order systems are matched at multiple selected points (or frequencies). Importantly, the physical realizability property of the original quantum system imposed by the laws of quantum mechanics is preserved under our tangential interpolatory projection. Error bounds are established for the proposed model reduction method and an avenue to select interpolation points is proposed. A passivity preserving model reduction method is also presented. Examples of both active and passive systems are provided to illustrate the merits of our proposed approach.