As the most widely distributed vegetation type in the Qilian Mountains region， subalpine meadows play an important role in maintaining local carbon and water fluxes and responding to climate change. Therefore， accurately detecting their phenological dynamics is crucial for a deeper understanding of mountain ecosystem functioning and its feedback to the climate system. In this study， we conducted a multisource image fusion and land surface phenology extraction experiment in a 15 km × 15 km test area in the northeastern Qilian Mountains， combining ground-level eddy flux data and multisource satellite remote sensing images. We used an Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model （ESTARFM） to fuse multisource images from the ETM+， OLI， and VIIRS sensors， reconstructing a high temporal （minimum 1 day） and high spatial （30 m） resolution time-series image dataset of the 2-band Enhanced Vegetation Index （EVI₂）， Normalized Difference Vegetation Index （NDVI） and Near-Infrared Reflectance of Vegetation （NIR_v） from 2013 to 2020. Based on this， we fitted the growth curves of the GPP flux tower and remote sensing vegetation index images using a Double Hyperbolic Tangent function （DHT） and Global Model Function （GMF）， respectively， and applied a dynamic threshold method to extract the start （SOS）， peak （POS） and end （EOS） of the growing season to evaluate the applicability of different fusion vegetation indices in extracting key phenological parameters of subalpine meadows. The results showed that ESTARFM fused images could accurately reflect the brightness and texture features of the real images， but cloud-contaminated pixels in the input images could also affect the fusion accuracy. At the site scale （without cloud pollution）， NIR_v and EVI₂ exhibited similar fusion accuracy， while at the pixel scale （with cloud pollution， cloud cover < 20%）， the fusion accuracy of NIR_v was significantly higher than that of EVI₂， indicating that NIR_v improved the sensitivity of vegetation partial reflectance in vegetation-bare soil mixed pixels in the algorithm and could maintain high fusion accuracy under cloud pollution conditions. For the growth curve fitting algorithm， DHT + GMF could accurately simulate the seasonal dynamics of the GPP flux tower and remote sensing vegetation indices， with determination coefficients above 0.960 and root mean square errors below 0.062. The comparison of phenology extraction accuracy of the three fused vegetation indices showed that NIR_v had the highest accuracy in extracting SOS and EOS， while NDVI had the highest accuracy in extracting POS， with deviations of 4 d （3 d）， 5 d （5 d）， and 4 d （6 d） at the site （pixel） scales， respectively.