Global Navigation Satellite System interferometric reflectometry (GNSS-IR), as a new remote sensing detection technology, can retrieve surface soil moisture (SM) by separating the modulation terms from the effective signal-to-noise ratio (SNR) data. However, traditional low-order polynomials are prone to over-fitting when separating modulation terms. Moreover, the existing research mainly relies on prior information to select satellites for SM retrieval. Accordingly, this study proposes a method based on empirical modal decomposition (EMD) and cross-correlation satellite selection (CCSS) for SM retrieval. This method intended to adaptively separate the modulation terms of SNR through the combination of EMD and an intrinsic mode functions (IMF) discriminant method, then construct a CCSS method to select available satellites, and finally establish a multisatellite robust estimation regression (MRER) model to retrieve SM. The results indicated that with EMD, the different feature components implied in the SNR data of different satellites could be adaptively decomposed, and the trend and modulation terms of the SNR could more accurately be acquired by the IMF discriminant method. The available satellites could be efficiently selected through CCSS, and the SNR quality of different satellites could also be classified at different accuracy levels. Furthermore, MRER could fuse the multisatellite phases well, which enhanced the accuracy of SM retrieval and further verified the feasibility and effectiveness of combining EMD and CCSS. When $r_m=0.600$ and $r_n=0.700$, the correlation coefficient (r) of the multisatellite combination reached 0.918, an improvement of at least 40% relative to the correlation coefficient of a single satellite. Therefore, this method can improve the adaptive ability of SNR decomposition, and the selection of satellites has high flexibility, which is helpful for the application and popularization of the GNSS-IR technology.