## 1 Introduction
[44] proposed an adaptive regression spline curve model to correct the vegetation error, which proved the feasibility of correcting the vegetation error based on the multipath effect. The different undulation of a terrain was also considered to have adverse impacts on the application of GNSS-IR. [PERSON] et al. [45] proposed a DSNR arc editing method to improve the retrieval accuracy of SM in undulating terrain. Meanwhile, based on the SNR of different carriers, many scholars tried retrieving SM using multisatellite, multi-frequency data. [PERSON] et al. [38] realized SM retrieval through GPS L2 SNR, and the retrieval results were in good agreement with the changing trend of SM. Compared with L2C, the SM retrieved using L5, L1, and L2P showed less difference [46; 47]. The B1 and B2 from the BDS also reflected well the trend of SM [48] ([PERSON] et al., 2019). Compared with a single satellite, SM retrieved by GNSS multisatellite fusion has higher stability and accuracy [49]. In addition, the SM retrieved by multisatellite dual-frequency combined multipath error can effectively improve the time resolution of SM [50]. Since low-quality SNR data usually cause abnormal phases, an SM retrieval method considering the detection and repair of abnormal phases was proposed, and the quality of the phases for each satellite was effectively improved [51]. However, the premise of using these methods to retrieve SM is to accurately separate the reflected signals from each satellite. Presently, the low-order polynomial is normally used to separate direct and reflected signals of satellites, but it is prone to local overfitting and has limited adaptability. Further improving the fitting quality of SNR, [PERSON] et al. [52] proposed a semi-empirical SNR model. Subsequently, the wavelet transform was also well verified to effectively reduce the noise information of SNR and showed a better fitting effect than the low-order polynomial [53]. Empirical Mode Decomposition (EMD) was introducedinto GNSS-IR monitoring to achieve the adaptive decomposition of SNR and improve the model stability [54]. Although wavelet transform and EMD have been used for satellite reflected signal separation, verifying their greater advantages over the traditional low-order polynomial in the fitting, the wavelet transform needs to choose appropriate wavelet bases and decomposition layers. The adaptive decomposition process of SNR data by EMD will generate many intrinsic mode functions (IMF), while there are fewer methods explicitly proposed to determine the satellite reflected signal in SM retrieval. In addition, the effective selection of GNSS satellites is the key to multisatellite combinations. The current satellite selection greatly depends on a priori information or empirical values, which further limits the promotion and application of the GNSS-IR technology. Accordingly, this paper proposes an SM retrieval method based on multisatellite fusion, which combines EMD and cross-correlation satellite selection (CCSS). This study aimed to adaptively decompose the SNR of each satellite through EMD and extract the modulation terms of the SNR using the IMF discriminant method. Meanwhile, these modulation terms were fitted by the nonlinear least squares algorithm (LLS) to gain the satellite phases. Then a cross-correlation satellite selection (CCSS) method was established to select the effective satellites. Subsequently, a multisatellite robust estimation regression (MRER) model was established for SM retrieval using the ICG III weight function. Finally, the model performance was assessed by comparing the retrieval results of different schemes, and the effectiveness of the satellite selection method was also verified.