The growth of data recorded by dense seismic arrays has stimulated the development of new array-based receiver function (RF) imaging techniques. This study examines the feasibility and performance of the least-squares migration (LSM) method, a state-of-art technique used in exploration seismology, to lithospheric imaging using teleseismic RFs. Taking advantage of a pair of forward (de-migration) and adjoint (migration) operators, the LSM casts migration as a regularized least-squares optimization problem. We employ the Split-step Fourier method to design the two operators and conduct wavefield propagation in heterogeneous media. Synthetic tests with models containing various Moho geometries demonstrate that LSM enables resolving interfaces at a higher resolution than conventional migration. Then LSM is applied to teleseismic data recorded by the Hi-CLIMB array deployed on the Tibetan Plateau. Considering the irregular and noisy recordings from field acquisition, we adopt signal processing algorithms, including the Radon transform and Singular Spectrum Analysis filter, to regularize the wavefields and precondition the RFs. The proposed workflow produces a significantly improved subsurface image than conventional methods, revealing new observations of 1) two well-defined interfaces at the base of the crust and 2) gently dipping mantle discontinuities extending continuously from the Lhasa Block to the Qiangtang Block. These structures could represent the imbricated Indian and Tibetan crust underlain by the underthrusting Indian lithosphere, implying that the Indian collisional front extends as far north as the Bangong-Nujiang suture. Overall, our study offers a new high-resolution RF imaging tool and inspires the future development of advanced array processing workflows.