Accurate mapping of forest carbon (C), nitrogen (N), and phosphorus (P) stocks is essential for advancing our understanding of global biogeochemical cycles. However, compared to carbon, large-scale quantification of forest N and P pools remains limited. We developed a two-step machine learning-based framework for estimating forest aboveground C, N, and P stocks by integrating field measurements, synthetic aperture radar (SAR) and optical remote sensing data, and soil and climatic variables. We tested the method in the permafrost region of the Greater Khingan Mountains, located at the southeastern margin of the boreal forest. We first introduced the fractional vegetation cover to adjust SAR backscatter coefficient, which substantially improved aboveground biomass (AGB) estimation compared with the existing China AGB map (R² value increased from 0.22 to 0.53). We then developed a novel triangular index based on time series of vegetation indices to represent vegetation nutrient uptake and accumulation. This index, together with AGB and auxiliary predictors, was used in a Gaussian Process Regression model to estimate aboveground C, N, and P stocks. The resulting estimates demonstrated higher accuracy than existing datasets, with R² values improving from 0.18, 0.01, and 0.44 to 0.83, 0.76, and 0.77 for C, N, and P stocks, respectively. These improvements were largely attributed to the inclusion of both the triangular index and AGB as key predictors in the model. This study presents an effective approach for large-scale mapping of aboveground C, N, and P stocks in boreal forest ecosystems, offering support for assessments of global carbon and nutrient cycles and for climate change research.