<p>High-precision estimation of regional zenith tropospheric delay (ZTD) is important for global navigation satellite system (GNSS) positioning and monitoring of the regional atmospheric environment. This study introduces a novel regional ZTD modeling method by fusing GNSS with surface meteorological observations, specifically for areas with sparse GNSS stations. First, regional weighted mean temperature and water vapor decrease factor correction models were established to calculate the surface meteorological ZTD. Subsequently, the established regional ZTD vertical correction model and simplified spherical cap harmonic analysis were used to estimate regional ZTD. To validate the accuracy of the proposed fusion method, we used the Yunnan Province in China as a representative area with sparse GNSS stations and conducted experiments using 26 GNSS and 109 surface meteorological stations. Comparative validation with ERA5 surface ZTD data showed that the fusion method improved accuracy by approximately 11.25 and 70.58% compared with using only surface meteorological ZTD or GNSS ZTD, respectively. Validation with continuously operating reference station (CORS) ZTD data further demonstrated that the fusion method enhanced the accuracy by approximately 14.39 and 43.52%, respectively. The proposed fusion method combines the benefits of GNSS and surface meteorological observations to improve the overall accuracy of ZTD estimation in sparse GNSS station areas.</p>