Studying the structures, properties and origins of the Earth's internal discontinuities is an important part in the efforts to understand the physical and chemical properties of the layered Earth, as well as to explore the dynamic processes and driving mechanisms of plate tectonics and the whole Earth system. Receiver function imaging is a well-known and widely-adopted seismological method in extracting the structural information of the Earth's internal discontinuities, and has become an indispensable tool to investigate the layering in structure and composition, and the thermal states and deformation behaviors of the crust and upper mantle, lithosphere-asthenosphere, mantle transition zone, and even shallow part of the lower mantle in the deep Earth. Since the receiver function method was proposed about half a century ago, great progress has been made in both methodology and application, targeting to subsurface structures of various spatial scales and from one- to three-dimension. In particular, with more and more seismic arrays being deployed in global and regional scales, and the continuous advancement of computing power and imaging theory during the last two decades, receiver function imaging has only become more powerful to constrain the subsurface structures. In this paper, we first briefly review the development history of the receiver function method. After introducing the basic principles involved, we then outline the major progress made during the last two decades in both methodology and application of this method, including but not limited to receiver function construction and forward modeling, receiver functions analysis for complex media or detailed discontinuity structures (e.g., anisotropy, dipping structures, irregular topography, sharpness of discontinuities), ray and wave-equation based receiver function migration in imaging crustal and upper mantle discontinuities, velocity inversion of receiver functions as well as its combination with other types of data. We focus mainly on the following three aspects: deconvolution techniques to construct receiver functions, imaging of discontinuity structures and inversion of velocity structures using receiver functions, with specific emphasis on the recent advances, challenges, and possible solutions. In the light of the emerging and future trends in seismology, we finally discuss the directions of receiver function studies from the viewpoints of both methodology and application.