In this contribution, we propose a multi-indicators comprehensive evaluation method of the Global Navigation Satellite System (GNSS) based on entropy weight-grey correlation analysis for reference satellite selection of PPP-AR. The comprehensive evaluation includes the observations index selection, the index normalization, the calculation of index entropy weight and the grey correlation analysis. The position dilution of precision (PDOP), carrier-to-noise density ratio (C/N0) and pseudorange multipath indicators are selected to construct a multi-indicator matrix, normalize the above indicators, and then calculate the entropy weight of each evaluation indicator to determine the indicator weights. Meanwhile, after selecting the optimal indicator set, the matrix is constructed to find the grey correlation coefficient and finally the grey correlation degree. According to the above method, the quality ranking of satellite observations for each epoch can be obtained, and the observation value with the highest ranking is used as the reference satellite during the PPP-AR. One-week observations from 243 MGEX stations are selected to conduct GPS-only, Galileo-only and BDS-3-only kinematic PPP-AR, the reference satellite selection method using the highest-elevation and the proposed method is applied, respectively. The results show that, compared with the traditional method, the positioning accuracies, convergence time and ambiguity fixed rates all can be improved using the proposed method. The positioning accuracies of GPS can be improved by about 10.81%, 5.10%, 2.95% and 5.54% in the east(E), north(N), up(U) and three-dimensional(3D) directions, respectively; 7.94%, 5.11%, 5.04% and 6.02% for Galileo; and 14.94%, 5.80%, 3.98% and 8.81% for BDS-3. The average improvements of convergence time in the east, north and up directions are 4.67%, 2% and 4% for GPS, BDS-3 and Galileo, respectively. The ambiguity fixed rates are improved by 7.31%, 6.89% and 1.45% for GPS, BDS-3 and Galileo between the 80%-100% range, respectively.