A method for scattered wave imaging in 3-D with both teleseismic P and S wave receiver function data is introduced. The approach relies on body-wave scattering kernels that are derived from the adjoint data sensitivity kernels which are typically used for full waveform tomography. The forward problem is approximated using ray theory, yielding a computationally efficient imaging algorithm that can resolve dipping and discontinuous interfaces using both P and S wave receiver functions. Traveltime fields for the incident teleseismic arrivals and the receiver point-sources are obtained by solving the Eikonal equation using a Fast Marching code which can handle a 3-D reference velocity model. An energy stable finite-difference method is used to simulate elastic wave propagation in a 2-D hypothetical subduction zone model. The resulting synthetic P and S wave receiver function data sets are used to validate the imaging method. The kernel images are compared with those generated by the Generalized Radon Transform and Common Conversion Point Stacking methods. These results demonstrate the potential of the kernel imaging approach for constraining lithospheric structure in complex geologic environments with sufficiently dense recordings of teleseismic data. Potential imaging targets include short-wavelength compositional variations in the mantle wedge and the slab's lithosphere-asthenosphere boundary.