Abstract: Accurate and rapid magnitude estimation is of paramount importance for earthquake early warning systems （EEWs）. Traditional magnitude estimation methods based on a single characteristic parameter of the initial seismic wave are widely used in EEWs. However, these empirical formulae, established by a single characteristic parameter, fail to fully exploit the information related to magnitude contained in the initial seismic wave, significantly limiting the effectiveness of magnitude estimation. To improve the accuracy of magnitude estimation in EEWs, this paper proposes a Gaussian process regression （GPR） based method that can estimate magnitudes in both scenarios: with and without hypocentral distance. The proposed method, GPR-M, uses multiple characteristic parameters from the time domain, frequency domain, and time-frequency domain as inputs, while GPR-M-R incorporates hypocentral distance. Both methods estimate magnitude by integrating various aspects of information from the initial seismic wave. The study utilized 33698 vertical acceleration records from the Japanese Kiban-Kyoshin Network （KiK-net） for training and testing, and 5353 vertical acceleration records from the Chilean Simulation Based Earthquake Risk and Resilience of Interdependent Systems and Networks （SIBER-RISK） for generalization testing. Additionally, the method’s practical application was validated using three typical earthquake cases in China, with M_S5.4, M_S6.4, and M_S8.0. The performance of the GPR method was compared with the widely adopted \tau ^\max _\mathrmp and P _d methods. The test results from the Japanese records indicate that for initial seismic waves of 3 to 10 s, both GPR-M and GPR-M-R outperform the \tau ^\max _\mathrmp and P _d methods in magnitude estimation. Specifically, the standard deviation of estimation errors for the GPR-M method is reduced by approximately 52.53% to 61.20% compared with the \tau ^\max _\mathrmp method, while the GPR-M-R method reduces the standard deviation of estimation errors by about 37.72% to 41.21% compared with the P _d method. For larger earthquakes （M_W≥6.5）, the magnitude saturation phenomenon is less pronounced in the GPR-M and GPR-M-R methods compared with the \tau ^\max _\mathrmp and P _d methods. The accuracy of magnitude estimation for M_W≥6.5 is improved by 1.4 to 1.5 times with the GPR-M method compared with the \tau ^\max _\mathrmp method, and by 1.2 to 1.45 times with the GPR-M-R method compared with the P _d method. The test results from the Chilean data demonstrate that both the GPR-M and GPR-M-R methods can effectively estimate earthquake magnitudes in Chile. The standard deviation of estimation errors for the GPR-M method is reduced by approximately 53.08% to 55.13% compared with the \tau ^\max _\mathrmp method, and the GPR-M-R method reduces the standard deviation of estimation errors by about 35.88% to 36.59% compared with the P _d method, showing excellent generalization capability. The test results from the three Chinese earthquake cases further confirmed that the GPR methods exhibit better accuracy and reliability compared with the \tau ^\max _\mathrmp and P _d methods. The GPR method can significantly improve the accuracy of magnitude estimation in EEWs and is not affected by regional differences. In conclusion, this study presents a novel GPR-based magnitude estimation method that integrates multiple seismic wave features and optionally incorporates hypocentral distance information. The method demonstrates superior performance in terms of accuracy, reliability, and generalization ability compared with traditional single-parameter approaches. By effectively reducing estimation errors and mitigating magnitude saturation issues, particularly for larger earthquakes, the proposed GPR method offers significant potential for improving the effectiveness of EEWs across diverse geographical regions.