Mapping Forest Stock Volume Using Phenological Features Derived from Time-Serial Sentinel-2 Imagery in Planted Larch

: As one of the important types of forest resources, mapping forest stock volume (FSV) in larch ( Larix decidua ) forests holds significant importance for forest resource management, carbon cycle research, and climate change monitoring. However, the accuracy of FSV mapping using common spectral and texture features is often limited due to their failure in fully capturing seasonal changes and growth cycle characteristics of vegetation. Phenological features can effectively provide essential information regarding the growth status of forests. In this study, multi-temporal Sentinel-2 satellite imagery were initially acquired in the Wangyedian Forest Farm in Chifeng City, Inner Mongolia. Subsequently, various phenological features were extracted from time series variables constructed by Gaussian Process Regression (GPR) using Savitzky–Golay filters, stepwise differentiation, and Fourier transform techniques. The alternative features were further refined through Pearson’s correlation coefficient analysis and the forward selection algorithm, resulting in six groups of optimal subsets. Finally, four models including the Random Forest (RF), K-Nearest Neighbors (KNN), Support Vector Machine (SVM), and Multiple Linear Regression (MLR) algorithms were developed to estimate FSV. The results demonstrated that incorporating phenological features significantly enhanced model performance, with the SVM model exhibiting the best performance—achieving an R 2 value of 0.77 along with an RMSE value of 46.36 m 3 /hm 2 and rRMSE value of 22.78%. Compared to models without phenological features, inclusion of these features led to a 0.25 increase in R 2 value while reducing RMSE by 10.40 m 3 /hm 2 and rRMSE by 5%. Overall, integration of phenological feature variables not only improves the accuracy of larch forest FSV mapping but also has potential implications for delaying saturation phenomena.


Introduction
Larch (Larix decidua) grows mainly in temperate and sub-boreal regions of the Northern Hemisphere, especially in high-altitude mountainous areas and higher northern latitudes [1].Researchers like Duan Beixing have revealed the significant drought resistance and vitality exhibited by larches.This species can grow in poor soil conditions, playing a crucial role in maintaining local ecological balance [2,3].With the escalating global climate change, larch (Larix decidua), an important and widely distributed tree species, shows significant spatiotemporal variations in its distribution, growth, and biomass accumulation.As a species with excellent timber properties and important ecological functions, monitoring the dynamics of larch is essential for understanding its role in the global carbon cycle, guiding forest management strategies and protecting biodiversity [4].In this context, accurate mapping forest stock volume is critical for forest resource management, ecosystem protection, and research on the global carbon cycle.Therefore, precise monitoring of larch forest stock volume is vital for sustainable forest management, providing a scientific basis for forestry harvesting, renewal, and protection strategies, and playing a key role in assessing and maintaining ecosystem diversity and stability.
Currently, the application of remote sensing inversion technology for forest stock volume has reached a state of maturity [5].Optical remote sensing data, characterized by its multi-temporal, multi-resolution, and multi-spectral advantages as well as ease of acquisition, has been extensively employed in mapping forest stock volume [6][7][8].This approach primarily relies on satellite platforms such as Landsat, MODIS, and Sentinel-2 to provide multispectral and hyperspectral images [9,10].By processing these images, key remote sensing features including spectral variables, vegetation indexes, and texture variables are widely extracted to construct models based on statistical relationships with ground-measured data [11][12][13][14].However, traditional features have certain limitations.For instance, atmospheric conditions, land cover types, and observation angles can affect conventional remote sensing features and subsequently impact estimation accuracy.Moreover, such features often fail to fully capture seasonal changes and growth cycle characteristics of vegetation, restricting their effectiveness in applications involving phenological variations.
Phenological features, such as the onset, peak, and decline phases of vegetation growth, can effectively provide crucial information regarding the growth status of vegetation and ecosystem functions [15][16][17][18].The incorporation of phenological features has introduced a novel perspective in the realm of forest resource monitoring.For instance, Zhang Yali et al. have enhanced the accuracy of estimating forest aboveground biomass by integrating phenological features [19].In comparison to traditional remote sensing features, phenological characteristics more accurately depict seasonal changes and growth cycles of vegetation, thereby augmenting the precision of forest resource monitoring.However, several challenges still exist in applying phenological features for mapping FSV, including accurate extraction techniques for these features, handling noise in large-scale spatial applications, and addressing universality issues across different tree species.
To address these challenges, researchers have developed various methodologies for extracting phenological features, encompassing time series analysis, machine learning algorithms, and model-based approaches [20][21][22].Time series analysis leverages remote sensing data collected at consecutive time points to discern crucial phenological stages (SOS, LOS, EOS) by scrutinizing the trends in vegetation indices over time.Machine learning algorithms enable the processing of extensive datasets and facilitate the extraction of phenological features by discerning patterns within the data.Model-based approaches aim to simulate vegetation growth processes based on ecological principles, thereby enabling the extraction of phenological features.
To improve the accuracy of mapping FSV in larch forests, this study first examines the limitations of existing single methods for extracting phenological features, and then presents a method that integrates time series analysis with ecological modeling to predict vegetation growth.Utilizing the Sentinel-2 satellite remote sensing platform, spectral variables, vegetation indices, and texture features are extracted from time series remote sensing images.Subsequently, the Gaussian Process Regression (GPR) algorithm is employed to smoothen and reconstruct the features extracted from the remote sensing data, aiming to eliminate data noise.To assess the potential of phenological features in mapping the Forest Stock Volume (FSV) in planted larch forests, various types of phenological features are extracted from time series images using three methods, namely Savitzky-Golay filtering, stepwise differentiation, and Fourier transform.Additionally, four machine learning models-Random Forest (RF), K-Nearest Neighbors (KNN), Support Vector Machine (SVM), and Multiple Linear Regression (MLR)-are applied to construct models with the optimal subset derived from the Forward Selection algorithm.The objective is to effectively

Study Area
The study area is located in the Wangyedian Forest Farm, situated in the southwestern part of Karakin Banner, Chifeng City, Inner Mongolia Autonomous Region (Figure 1) [23].The geographic coordinates of the forest farm range from 118 • 09 ′ to 118 • 30 ′ east longitude and 41 • 21 ′ to 41 • 39 ′ north latitude, respectively.This region predominantly comprises mid-mountain and low-mountain landscapes, with an average altitude ranging from 800 m to 1890 m.The topography exhibits higher elevations towards the southwest, gradually declining towards the northeast.Through a comprehensive investigation and statistical analysis of forest resources within Wangyedian Forest, it was determined that the total forest cover spans approximately 350,000 mu, out of which plantation forests account for around 176,000 mu.Notably, key species such as Larch plantations (Larix principis-rupprechtii), Chinese pine plantations (Pinus tabuliformis), natural secondary Dahurian birch forests (Betula dahurica), and natural secondary silver birch forests (Betula platyphylla) constitute the primary forest resources.(SVM), and Multiple Linear Regression (MLR)-are applied to construct models with the optimal subset derived from the Forward Selection algorithm.The objective is to effectively identify the key phenological stages of vegetation and accurately map the FSV in planted larch forests.

Study Area
The study area is located in the Wangyedian Forest Farm, situated in the southwestern part of Karakin Banner, Chifeng City, Inner Mongolia Autonomous Region (Figure 1) [23].The geographic coordinates of the forest farm range from 118°09′ to 118°30′ east longitude and 41°21′ to 41°39′ north latitude, respectively.This region predominantly comprises mid-mountain and low-mountain landscapes, with an average altitude ranging from 800 m to 1890 m.The topography exhibits higher elevations towards the southwest, gradually declining towards the northeast.Through a comprehensive investigation and statistical analysis of forest resources within Wangyedian Forest, it was determined that the total forest cover spans approximately 350,000 mu, out of which plantation forests account for around 176,000 mu.Notably, key species such as Larch plantations (Larix principis-rupprechtii), Chinese pine plantations (Pinus tabuliformis), natural secondary Dahurian birch forests (Betula dahurica), and natural secondary silver birch forests (Betula platyphylla) constitute the primary forest resources.

Ground Measured Data
In this study, the randomized stratified sampling method [24] was employed to establish ground measured samples based on the distribution of forest volume in planted larch forests.A total of 79 samples were established between October 2017 and September 2019.Each sample has a size of 25 m × 25 m.The position and angles of the center are precisely measured using GNSS, with an error margin of less than 10 cm.Within each sample, various parameters including tree height (H/m), diameter at breast height (D/cm), and crown diameter (C/m) were measured.Additionally, elevation, slope, and aspect were also recorded at the center of each sample plot.Based on the ground measured

Ground Measured Data
In this study, the randomized stratified sampling method [24] was employed to establish ground measured samples based on the distribution of forest volume in planted larch forests.A total of 79 samples were established between October 2017 and September 2019.Each sample has a size of 25 m × 25 m.The position and angles of the center are precisely measured using GNSS, with an error margin of less than 10 cm.Within each sample, various parameters including tree height (H/m), diameter at breast height (D/cm), and crown diameter (C/m) were measured.Additionally, elevation, slope, and aspect were also recorded at the center of each sample plot.Based on the ground measured height and DBH of each tree, the volume of the individual trees was calculated using a binary volume table according to the following formula below, quantifying the relationships among diameter at breast height (DBH), tree height, and volume of larch.The special formula is expressed as follows: where D is the diameter at breast height of a single tree and H is the height of a single tree.
The forest stem volume of the sample was determined by summing the volumes of all trees within the sample.The average diameter at breast height (DBH) of larch ranged from 6.18 to 29.21 cm, while the average height ranged from 5.86 to 23.70 m (Figure 2).
height and DBH of each tree, the volume of the individual trees was calculated using a binary volume table according to the following formula below, quantifying the relationships among diameter at breast height (DBH), tree height, and volume of larch.The special formula is expressed as follows:  = −0.001498+ 0.00007 ×  + 0.000901 ×  + 0.000032 ×  ×  (1) where D is the diameter at breast height of a single tree and H is the height of a single tree.
The forest stem volume of the sample was determined by summing the volumes of all trees within the sample.The average diameter at breast height (DBH) of larch ranged from 6.18 to 29.21 cm, while the average height ranged from 5.86 to 23.70 m (Figure 2).

Time Series Sentinel-2 Data
The remote sensing data in this study were acquired from the Sentinel-2 satellite, launched by the European Space Agency (ESA) in 2015.A total of nine Sentinel-2A L2A images from April to December 2019 was obtained from the Copernicus Open Access Hub (https://scihub.copernicus.eu/accessed on 8 October 2023), excluding those affected by snow cover between January and March.These images underwent specific model processing techniques for acquisition purposes.Subsequently, all band data were resampled at a resolution of 20 m using SNAP 6.0 software for further analysis.Finally, these datasets were imported into Orfeo Toolbox (https://www.orfeo-toolbox.org/,accessed on 10 October 2023) to undergo a series of preprocessing steps including topographic correction and relative radiometric correction [25], ensuring accurate and reliable data for subsequent analysis.Additionally, all these data were divided into four seasons based on the growth characteristics of larch (Figure 3).

Time Series Sentinel-2 Data
The remote sensing data in this study were acquired from the Sentinel-2 satellite, launched by the European Space Agency (ESA) in 2015.A total of nine Sentinel-2A L2A images from April to December 2019 was obtained from the Copernicus Open Access Hub (https://scihub.copernicus.eu/accessed on 8 October 2023), excluding those affected by snow cover between January and March.These images underwent specific model processing techniques for acquisition purposes.Subsequently, all band data were resampled at a resolution of 20 m using SNAP 6.0 software for further analysis.Finally, these datasets were imported into Orfeo Toolbox (https://www.orfeo-toolbox.org/,accessed on 10 October 2023) to undergo a series of preprocessing steps including topographic correction and relative radiometric correction [25], ensuring accurate and reliable data for subsequent analysis.Additionally, all these data were divided into four seasons based on the growth characteristics of larch (Figure 3).
height and DBH of each tree, the volume of the individual trees was calculated using a binary volume table according to the following formula below, quantifying the relationships among diameter at breast height (DBH), tree height, and volume of larch.The special formula is expressed as follows:  = −0.001498+ 0.00007 ×  + 0.000901 ×  + 0.000032 ×  ×  (1) where D is the diameter at breast height of a single tree and H is the height of a single tree.
The forest stem volume of the sample was determined by summing the volumes of all trees within the sample.The average diameter at breast height (DBH) of larch ranged from 6.18 to 29.21 cm, while the average height ranged from 5.86 to 23.70 m (Figure 2).

Time Series Sentinel-2 Data
The remote sensing data in this study were acquired from the Sentinel-2 satellite, launched by the European Space Agency (ESA) in 2015.A total of nine Sentinel-2A L2A images from April to December 2019 was obtained from the Copernicus Open Access Hub (https://scihub.copernicus.eu/accessed on 8 October 2023), excluding those affected by snow cover between January and March.These images underwent specific model processing techniques for acquisition purposes.Subsequently, all band data were resampled at a resolution of 20 m using SNAP 6.0 software for further analysis.Finally, these datasets were imported into Orfeo Toolbox (https://www.orfeo-toolbox.org/,accessed on 10 October 2023) to undergo a series of preprocessing steps including topographic correction and relative radiometric correction [25], ensuring accurate and reliable data for subsequent analysis.Additionally, all these data were divided into four seasons based on the growth characteristics of larch (Figure 3).

Spectral and Texture Feature Extraction
To accurately map FSV using Sentinel-2 images, it is crucial to extract valuable features through various methodologies.Typically, three types of features are considered, spectral features [26], vegetation indices, and texture features.In this study, 12 single-band information (including B1-B12) with a spatial resolution better than 20 m after resampling were selected as independent spectral features.Moreover, the vegetation index (VI) [14][15][16][17] was also employed as an indicator of surface vegetation cover and growth based on optical remote sensing data.After conducting comparative analysis and comprehensive assessment, this study opted for the selection of the top 15 vegetation indices (RVI, IPVI, NDVI, NDVI45, NDVI58a, NDVI56, GNDVI, IRECI, SAVI, ARV I, PSSRa, MTCI, MCARI, REIP and NDII) as key variables for mapping forest stock volume in subsequent modeling and research (Table 1).Additionally, texture elements can provide crucial insights into forest density, tree species distribution, and growth status [14].In this study, the grayscale co-production matrix method was employed to extract eight types of texture information from Sentinel-2A bands (Entropy, Contrast, Second moment, Homogeneity, Dissimilarity, Mean, Variance, and Correlation) [18], with a window size set at 5 × 5 [27][28][29].

Gaussian Process Regression to Reconstruct Time Series
Phenological features are typically derived from time series data in most cases.To mitigate observation errors, it is imperative to initially reconstruct the time series.The Gaussian Process Regression (GPR) algorithm is employed in this study to reconstruct the time series, and subsequently illustrate the extracted candidate information in Figure 4.A Gaussian process can be conceptualized as an ensemble of functions, where any set of function values (corresponding to any set of input values) follows a joint Gaussian distribution.Gaussian processes [30][31][32] are usually defined by the mean function m(x) and the covariance function or kernel function k(x, x ′ ).
where m(x) is the mean value function, which represents the expected function value at a given point x.The kernel function k(x, x ′ ) represents the covariance of the function values between the two points x and x ′ .
mote sensing images from January to March, we employ a methodology that utilizes remote sensing data from October to December of the same year as a surrogate for estimating and supplementing the remote sensing features from March to January of the corresponding year.This approach acknowledges the presence of certain similarities and continuity in remote sensing features between adjacent months within the same geographical region.

Phenological Features (1) Extraction phenology features by S-G filter
Previous studies have demonstrated the significant potential of the Savitzky-Golay (S-G) filtering technique for accurately reconstructing time series data.Consequently, the S-G filter was used to smoothen the time series normalized vegetation index (NDVI) by fitting a low-order polynomial to a subset of selected data points [33][34][35] (Figure 5).Its expression is shown below: where   is the filtered value,    is the original data point,   is the filter coefficient, and  is the polynomial half the size of the window (the window size is 2 + 1).Given a set of observation data points D = {(x1, y1), . . . ,(xn, yn)}, it is possible to predict the output value of the new point x * .In this study, due to the unavailability of remote sensing images from January to March, we employ a methodology that utilizes remote sensing data from October to December of the same year as a surrogate for estimating and supplementing the remote sensing features from March to January of the corresponding year.This approach acknowledges the presence of certain similarities and continuity in remote sensing features between adjacent months within the same geographical region.

Phenological Features (1) Extraction phenology features by S-G filter
Previous studies have demonstrated the significant potential of the Savitzky-Golay (S-G) filtering technique for accurately reconstructing time series data.Consequently, the S-G filter was used to smoothen the time series normalized vegetation index (NDVI) by fitting a low-order polynomial to a subset of selected data points [33][34][35] (Figure 5).Its expression is shown below: where y ′ i is the filtered value, y i+j is the original data point, c j is the filter coefficient, and m is the polynomial half the size of the window (the window size is 2m + 1).In this study, a multi-temporal maximum value synthesis (MVC) method based on the S-G filtering algorithm was used to analyze the NDVI images of the current period.Subsequently, outliers were eliminated through a pixel-wise review of each image, and maximum values were determined by comparing them with the neighboring periods'  In this study, a multi-temporal maximum value synthesis (MVC) method based on the S-G filtering algorithm was used to analyze the NDVI images of the current period.Subsequently, outliers were eliminated through a pixel-wise review of each image, and maximum values were determined by comparing them with the neighboring periods' NDVI images.Interpolation methods were also utilized to recover data from rejected null image elements to enhance analysis accuracy.By analyzing the dynamics of the vegetation growth curves, botanical characteristics such as the start of the season (SOS), the end of the season (EOS), and the length of the season (LOS) [36] were extracted from key time nodes during the growing season (Figure 6).Furthermore, changes in vegetation indices and their textural features during the growing season of joint larch (April to September) could also serve as a phenological variable that reflects the growing season.In this study, a multi-temporal maximum value synthesis (MVC) method based on the S-G filtering algorithm was used to analyze the NDVI images of the current period.Subsequently, outliers were eliminated through a pixel-wise review of each image, and maximum values were determined by comparing them with the neighboring periods' NDVI images.Interpolation methods were also utilized to recover data from rejected null image elements to enhance analysis accuracy.By analyzing the dynamics of the vegetation growth curves, botanical characteristics such as the start of the season (SOS), the end of the season (EOS), and the length of the season (LOS) [36] were extracted from key time nodes during the growing season (Figure 6).Furthermore, changes in vegetation indices and their textural features during the growing season of joint larch (April to September) could also serve as a phenological variable that reflects the growing season.(2) Extraction phenology features by time series seasonal difference Given the deciduous nature of larch, its physiological and ecological features exhibit significant seasonal variations.In this study, the mean values of the featur e variables were initially calculated within four seasons (spring: April and May; s ummer: June, July, and August; autumn: September and October; winter: Novemb er and December) to represent phenological features during spring, summer, autu mn, and winter, respectively.Subsequently, a stepwise differencing method [37,38]  (2) Extraction phenology features by time series seasonal difference Given the deciduous nature of larch, its physiological and ecological features exhibit significant seasonal variations.In this study, the mean values of the feature variables were initially calculated within four seasons (spring: April and May; summer: June, July, and August; autumn: September and October; winter: November and December) to represent phenological features during spring, summer, autumn, and winter, respectively.Subsequently, a stepwise differencing method [37,38] was employed to differentiate the periodic data.The seasonally varying phenological features of larch were extracted by calculating the differences between features extracted from the current season and the previous season.Stepwise differencing is particularly useful for analyzing or predicting time series data or signals as it stabilizes non-stable time series data by eliminating trends and seasonal components.The detailed flow chart is shown in Figure 7.
(3) Phenological features derived by Fourier transform Time series spectral and textural information reconstructed by the Gaussian Process Regression (GPR) was transformed into the frequency domain using the Fourier transform algorithm [39][40][41], as illustrated in Figure 8, enabling the extraction of phenological indicators such as seasonal integrals and growing season amplitudes.The expression for the continuous-time Fourier transform is as follows: where X( f ) is the Fourier transform of the signal x(t), f is the frequency, and e −j2π f t is the complex exponential function used to decompose the signal into sine and cosine waves.Assuming that the frequency range of NDVI index variation is [ f 1 , f 2 ], the frequency components of the growing season by filtering X( f ) can be extracted as follows: where N growth is the number of valid frequencies for k, i.e., f 2 − f 1 + 1.  (3) Phenological features derived by Fourier transform Time series spectral and textural information reconstructed by the Gaussian Process Regression (GPR) was transformed into the frequency domain using the Fourier transform algorithm [39][40][41], as illustrated in Figure 8, enabling the extraction of phenological indicators such as seasonal integrals and growing season amplitudes.The expression for the continuous-time Fourier transform is as follows: where () is the Fourier transform of the signal (),  is the frequency, and  is the complex exponential function used to decompose the signal into sine and cosine waves.Assuming that the frequency range of NDVI index variation is [ ,  ] , the frequency components of the growing season by filtering () can be extracted as follows: where  growth is the number of valid frequencies for k, i.e.,  −  + 1.

Feature Selections
The mapping of forest stock volume is typically achieved by utilizing optimal feature sets extracted from remotely sensed data through feature selection methods [42][43][44].Various feature selection methods have been widely applied to reduce model complexity and

Feature Selections
The mapping of forest stock volume is typically achieved by utilizing optimal feature sets extracted from remotely sensed data through feature selection methods [42][43][44].Various feature selection methods have been widely applied to reduce model complexity and enhance prediction accuracy.In this study, the Forward Selection algorithm (Forward Selection) [45] was employed to select the optimal feature sets in combination with four regression models.To evaluate the capability of phenological features extracted from different methods, all alternative features were divided into six subsets, the subset of spectral and texture features (S + T), the subset of spectral, texture features, and growing season features (S + T + Gs), the subset of spectral, texture, and four-season features (S + T + Fs), the subset of spectral, texture, and seasonal difference features (S + T + Sd), and the subset of spectral, texture, and mean amplitude features fusion (S + T + Fo), and all features (AFC).

Modeling and Accuracy Evaluation
To demonstrate the applicability of various regression models for mapping FSV in Larch forests, the optimal feature sets were employed to construct the models.In this study, there are four models [46][47][48][49], including Multiple Linear Regression (MLR), Random Forest (RF), K-Nearest Neighbors (KNN), and Support Vector Machine (SVM).We utilized wrapper feature selection methods to construct regression models in R 4.2.0 software for mapping the FSV.In this approach, the K-Nearest Neighbors (KNN) model required the tuning of the number of neighbors (K).For the Random Forest (RF) model, the parameters to adjust included the number of trees (ntree) and the number of random variables required to split each node (mtry).The Support Vector Machine (SVM) model necessitated tuning the penalty coefficient (C) and gamma.The Multiple Linear Regression (MLR) model did not require parameter adjustments.Furthermore, the leave-one-out cross-validation (LOOCV) method was also employed to evaluate the accuracy of models.Three accuracy indices, including the coefficient of determination (R 2 ), root mean square error (RMSE), relative root mean square error (rRMSE), were chosen to evaluate the performance of models.The aforementioned methods are summarized.Figure 9 illustrates the technical route for mapping forest stock volume using phenological feature data extraction.
where ŷ denotes the predicted FSV of samples, y i denotes the ground measured FSV of samples, y denotes the mean value of all samples, and n denotes the number of samples.
where  denotes the predicted FSV of samples, y denotes the ground measured FSV of samples,  denotes the mean value of all samples, and n denotes the number of samples.

Sensitivity between Features and FSV
In this study, all alternative features derived from time series Sentinel-2A images, encompassing spectral, texture, and phenological features, were classified into six groups for the purpose of mapping FSV.The sensitivity between these features and FSV was further examined through the absolute values of Pearson's correlation coefficients.Figure 10a-f illustrate the top ten correlation coefficients between each feature group and FSV.Among these features, phenological features extracted from the growing season (Figure 10c) and seasonal differences (Figure 10e) exhibited strong correlations with FSV, as evidenced by their absolute values of Pearson's correlation coefficients, ranging from 0.39 to

Sensitivity between Features and FSV
In this study, all alternative features derived from time series Sentinel-2A images, encompassing spectral, texture, and phenological features, were classified into six groups for the purpose of mapping FSV.The sensitivity between these features and FSV was further examined through the absolute values of Pearson's correlation coefficients.Figure 10a-f illustrate the top ten correlation coefficients between each feature group and FSV.Among these features, phenological features extracted from the growing season (Figure 10c) and seasonal differences (Figure 10e) exhibited strong correlations with FSV, as evidenced by their absolute values of Pearson's correlation coefficients, ranging from 0.39 to 0.50 and 0.33 to 0.54, respectively.Conversely, traditional textural features demonstrated weaker correlations with FSV based on their absolute values of correlation coefficients ranging from 0.23 to 0.39.It is implied that phenological features display higher sensitivity with FSV than spectral and texture features.0.50 and 0.33 to 0.54, respectively.Conversely, traditional textural features demonstrated weaker correlations with FSV based on their absolute values of correlation coefficients ranging from 0.23 to 0.39.It is implied that phenological features display higher sensitivity with FSV than spectral and texture features.

Results of Estimated FSV without Phenological Features
The capability of phenological features was evaluated by initially deriving the estimated FSV using four models.These models utilized a subset of spectral and textural features extracted from Sentinel-2 images acquired between April and September, which corresponds to the larch growing season.The accuracy of the models was assessed using the coefficient of determination (R 2 ), root mean square error (RMSE), and relative root mean square error (rRMSE), as shown in Table 2.It is evident that the accuracy of the mapped FSV varied depending on the choice of model and acquisition date of images.The values for R 2 ranged from 0.14 to 0.52, while rRMSE ranged from 27.89% to 37.47%.Notably, the highest R 2 value and lowest rRMSE were achieved when KNN was employed with the images acquired in July and September, respectively.Furthermore, it was observed that the highest accuracy for mapped FSV was obtained from images acquired in September.

Results of Estimated FSV without Phenological Features
The capability of phenological features was evaluated by initially deriving the estimated FSV using four models.These models utilized a subset of spectral and textural features extracted from Sentinel-2 images acquired between April and September, which corresponds to the larch growing season.The accuracy of the models was assessed using the coefficient of determination (R 2 ), root mean square error (RMSE), and relative root mean square error (rRMSE), as shown in Table 2.It is evident that the accuracy of the mapped FSV varied depending on the choice of model and acquisition date of images.The values for R 2 ranged from 0.14 to 0.52, while rRMSE ranged from 27.89% to 37.47%.Notably, the highest R 2 value and lowest rRMSE were achieved when KNN was employed with the images acquired in July and September, respectively.Furthermore, it was observed that the highest accuracy for mapped FSV was obtained from images acquired in September.To further analyze the results of the estimated FSV, the optimal results (images acquired in September) were chosen to depict the scatter plots between the ground measured and predicted FSV (Figure 11).It is evident that overestimation occurred in young Larch forests, while underestimation was observed in mature ones.Even for samples with FSV ranging from 100 m 3 /hm 2 to 300 m 3 /hm 2 , the accuracy of the results was significantly limited by errors arising from normal spectral and textural features.Therefore, it is imperative to incorporate new features to enhance estimation accuracy.

Results of Estimated FSV after Incorporating Phenological Features
In this study, further explorations were conducted to evaluate the potential of phenological features for mapping forest FSV in Larch forests.Spectral and textural features were extracted from the images acquired in September, while four types of phenological features were derived from time series images using various methods.Subsequently, six

Results of Estimated FSV after Incorporating Phenological Features
In this study, further explorations were conducted to evaluate the potential of phenological features for mapping forest FSV in Larch forests.Spectral and textural features were extracted from the images acquired in September, while four types of phenological features were derived from time series images using various methods.Subsequently, six subsets (S + T, S + T + Gs, S + T + Fs, S + T + Sd, S + T + Fo, AFC) were obtained from an alternative feature set extracted from time series Sentinel-2A images using the Forward Selection algorithm coupled with four corresponding model regressors.A comprehensive overview of the results for each group can be found in Table A1.These selected feature subsets were then utilized to map FSV using four models and the outcomes are presented in Table 3.The results indicate a significant improvement in the accuracy of mapped forest stock volume (FSV) after incorporating various phenological features.Specifically, when compared to subset 1 (S + T), subsets 2 to 5 showed determination coefficients (R 2 ) ranging from 0.38 to 0.77 and relative root mean square error (rRMSE) ranging from 22.78% to 33.03%.This suggests that phenological features have greater potential for capturing changes in forest structure parameters compared to normal spectral and textural features, particularly for Larch forests.Notably, optimal training performance was achieved by including all types of phenological features (AFC) using SVM regressor, resulting in an R 2 value of 0.77 and an rRMSE value of 22.78%.Therefore, incorporating phenological features can improve the accuracy and reliability of mapping FSV, with the degree of improvement closely associated with methods of extracting phenological features.To further explore the response of different feature factors to FSV, scatter plots were generated between predicted and measured FSV for the model with the highest estimation accuracy (Figure 12).Furthermore, it is demonstrated that phenological features hold great potential for improving the accuracy of mapping FSV in Larch forests, resulting in reduced discrepancies between ground measured and predicted FSV values.It is further proven that by incorporating diverse phenological features, the accuracy of FSV estimation can be significantly enhanced, particularly in terms of mitigating underestimation errors.
Based on the accuracy of results (Table 3), the map of FSV in planted Larch forests was generated by Optimal model SVM using spectral and textural features extracted from images acquired in September and four types of phenological features (Figure 13).Within our study area, plantation forests primarily displayed a concentration of FSV ranging from 100 to 300 m 3 /hm 2 , accounting for over 25.28% of the total coniferous plantation area.Specifically, the average volume of coniferous larch within this region was calculated as 208.36 m 3 /hm 2 , with a standard deviation of 33.58 m 3 /hm 2 .
tween predicted and measured FSV for the model with the highest estimation accuracy (Figure 12).Furthermore, it is demonstrated that phenological features hold great potential for improving the accuracy of mapping FSV in Larch forests, resulting in reduced discrepancies between ground measured and predicted FSV values.It is further proven that by incorporating diverse phenological features, the accuracy of FSV estimation can be significantly enhanced, particularly in terms of mitigating underestimation errors.Based on the accuracy of results (Table 3), the map of FSV in planted Larch forests was generated by Optimal model SVM using spectral and textural features extracted from images acquired in September and four types of phenological features (Figure 13).Within our study area, plantation forests primarily displayed a concentration of FSV ranging from 100 to 300 m 3 /hm 2 , accounting for over 25.28% of the total coniferous plantation area.Specifically, the average volume of coniferous larch within this region was calculated as 208.36 m 3 /hm 2 , with a standard deviation of 33.58 m 3 /hm 2 .

Sensitivity Analyses of Phenological Features and FSV
Numerous studies have demonstrated that larch leaves undergo significant phenological changes throughout the annual cycle, reflecting their adaptability to shifts in phenology [2,3].Specifically, from January to March, larches enter a dormant phase characterized by the absence of leaf growth, marking the non-growing season [36].Subsequently, from April to September, larches transition into the leaf growth phase with a rapid development of leaves, showcasing their robust growth vitality.By October, the trees initiate the leaf senescence phase as leaves exhibit signs of aging and yellowing, indicating the onset of autumn.Finally, from November to December, larches re-enter a dormant state.Therefore, larches can be considered as one of the key species for studying

Sensitivity Analyses of Phenological Features and FSV
Numerous studies have demonstrated that larch leaves undergo significant phenological changes throughout the annual cycle, reflecting their adaptability to shifts in phenology [2,3].Specifically, from January to March, larches enter a dormant phase characterized by the absence of leaf growth, marking the non-growing season [36].Subsequently, from April to September, larches transition into the leaf growth phase with a rapid development of leaves, showcasing their robust growth vitality.By October, the trees initiate the leaf senescence phase as leaves exhibit signs of aging and yellowing, indicating the onset of autumn.Finally, from November to December, larches re-enter a dormant state.Therefore, larches can be considered as one of the key species for studying phenological changes.
Previous research suggests that the phenological features reflect the forest's responses to environmental changes and can serve as indicators for assessing the health and productivity of forest [4].Therefore, this study aimed to explore the sensitivity between phenological features and FSV in a planted larch forest.Initially, nine Sentinel-2A remote sensing images were acquired, and textural and spectral features were extracted for each scene.After that, features were subsequently smoothed and reconstructed using the Gaussian Process Regression (GPR) algorithm.Subsequently, the phenological features were extracted from time series images using four methods, such as the Savitzky-Golay filter, stepwise differencing, and Fourier Transform.Additionally, Pearson's correlation analysis was conducted on the spectral variables, vegetation indices, texture features, and phenological features.The results revealed that compared to traditional spectral and textural variables, all phenological features extracted through various methods exhibited significant correlations.Among these features, phenological features extracted from growing season characteristics and seasonal change patterns demonstrated higher sensitivity levels (0.39~0.50).Furthermore, during the process of feature selection using the forward selection method, phenological features also exhibited a substantial contribution within the optimal feature set, thereby further substantiating their significance in estimating forest stock volume.
Existing research has unveiled that the introduction of phenological features brings a new perspective to the domain of forest resource monitoring.Phenological features have demonstrated high sensitivity and practical value in applications such as forest landform identification, monitoring the dynamics of plant communities, assessing the impacts of drought or pestilence, and forecasting the aboveground biomass of forests [50][51][52].Moreover, this study has transcended the approach of extracting phenological features through a single method by employing a variety of techniques, including Savitzky-Golay filters, stepwise differencing, and Fourier transforms.This approach allows for the layered extraction of phenological features from multiple datasets, thereby achieving high precision in the extraction of phenological features.Additionally, the Gaussian Process Regression (GPR) algorithm has effectively addressed the noise issues encountered in the extraction of phenological features in large-scale spatial applications.
In comparison to previous studies, this research provides additional evidence supporting the efficacy of phenological features as functional indicators for assessing forest ecosystem and its associated stock volume.

Contribution of Phenological Features in Mapping FSV
The impact of phenological features on the timing and efficiency of plant growth activities indirectly regulates the growth rate of forest.Previous studies have demonstrated that these features closely associated with forest stock volume can effectively reduce the errors caused by underestimation and delay saturation issues, thereby significantly enhancing the precision of mapping FSV [53,54].In this study, phenological features from Sentinel-2A data were incorporated into the regression model in a planted larch forest, and the results demonstrated that including various phenological features can improve the accuracy of mapping FSV to a certain extent.In particular, integrating all the phenological features yielded the optimal estimation accuracy (R 2 = 0.77, rRMSE = 22.78%).The inclusion of various phenological features resulted in a significant improvement in the determination coefficient (R 2 ), with an increase of 0.25.Additionally, the root mean square error (RMSE) was reduced by 10.40 m 3 /hm 2 and the relative root mean square error (rRMSE) decreased by approximately 5% when compared to traditional spectral and textural variables (Figure 14).logical features yielded the optimal estimation accuracy (R 2 = 0.77, rRMSE = 22.78%).The inclusion of various phenological features resulted in a significant improvement in the determination coefficient (R 2 ), with an increase of 0.25.Additionally, the root mean square error (RMSE) was reduced by 10.40 m 3 /hm 2 and the relative root mean square error (rRMSE) decreased by approximately 5% when compared to traditional spectral and textural variables (Figure 14).Furthermore, previous results have proven that the improvement of accuracy in mapping FSV is attributed to the incorporation of novel features, thereby mitigating the underestimation of predicted FSV [55].Scholars have confirmed that incorporating structural features (such as tree height and stand density) and topographical features (such as slope and aspect), as well as climatic characteristics (including long-term climate data like temperature and precipitation) on the basis of spectral and textural variables, can Furthermore, previous results have proven that the improvement of accuracy in mapping FSV is attributed to the incorporation of novel features, thereby mitigating the underestimation of predicted FSV [55].Scholars have confirmed that incorporating structural features (such as tree height and stand density) and topographical features (such as slope and aspect), as well as climatic characteristics (including long-term climate data like temperature and precipitation) on the basis of spectral and textural variables, can significantly reduce the underestimation of forest stock volume predictions, thereby enhancing the accuracy of the prediction models.This study also demonstrates significant effectiveness in reducing the underestimation of forest stock volume and improving the estimation accuracy by adding phenological features that reflect the seasonal changes in vegetation growth on the basis of spectral and texture variables.Furthermore, studies have shown that phenological features are highly effective in forest land object identification, monitoring plant community growth dynamics, assessing the impacts of drought or pests, and predicting forest aboveground biomass [50][51][52].
This further emphasizes the importance of employing a comprehensive range of phenological features in forest resource estimation while providing novel perspectives and theoretical foundations for future forest resource management and assessment.

Conclusions
To accurately map the forest stock volume of northern larch plantations, this study extracted spectral features, texture features, and phenological features from the Sentinel-2A time series data set.Optimal variables from each variable set were selected in conjunction with Pearson's correlation coefficients and forward feature selection methods; thereafter, four machine learning models were constructed-Multiple Linear Regression (MLR), K-Nearest Neighbors (KNN), Random Forest (RF), and Support Vector Machine (SVM).The results indicate that phenological variables demonstrate a higher sensitivity towards the forest volume of larch compared to spectral and textural variables, with Pearson's correlation coefficient reaching an absolute value of 0.54.Incorporating phenological variables into the forest volume estimation model significantly improved accuracy (R 2 increased by 0.25, RMSE decreased by 10.40 m 3 /hm 2 , and rRMSE decreased by 5%), thereby mitigating the underestimation of samples with larger accumulations.These findings suggest that phenological features extracted from Sentinel-2A time series data hold significant potential for enhancing forest stock estimation.Future studies will further explore the application of phenological features in estimating forest stock for other tree species while optimizing their extraction methods to more accurately reflect the forest stock of each specific tree species.and Z.Y.; resources, H.L. and J.L.; data curation, Q.L.; writing-original draft preparation, Q.L.; writing-review and editing, Q.L. and J.L.; visualization, Q.L.; supervision, H.L. and J.L.; project administration, H.L. and J.L; funding acquisition, H.L. and J.L.All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data supporting the findings of this study, including observed Forest Stock Volume (FSV) data from the sample plots and spatial distribution data of forest resources, are available from the corresponding author upon reasonable request.These data are not publicly accessible due to privacy and confidentiality concerns.The Sentinel-2 images are available from the Copernicus Open Access Hub (https://scihub.copernicus.eu/accessed on 8 October 2023).

Conflicts of Interest:
The authors declare no conflicts of interest.

Forests
phenological stages of vegetation and accurately map the FSV in planted larch forests.

Figure 1 .
Figure 1.The location of the study area and maps of ground samples: (a) is a map of the location of Inner Mongolia in China; (b) is a map of the location of the forest farm in Chifeng City; (c) is a sample point on the ground in the study area; (d,e) is a close-up photograph of a larch plantation forest.

Figure 1 .
Figure 1.The location of the study area and maps of ground samples: (a) is a map of the location of Inner Mongolia in China; (b) is a map of the location of the forest farm in Chifeng City; (c) is a sample point on the ground in the study area; (d,e) is a close-up photograph of a larch plantation forest.

Figure 2 .
Figure 2. The relationships among DBH, Height and FSV: (a) is the scatterplot of FSV and DBH; (b) is the scatterplot of FSV and Height.

Figure 2 .
Figure 2. The relationships among DBH, Height and FSV: (a) is the scatterplot of FSV and DBH; (b) is the scatterplot of FSV and Height.

Figure 2 .
Figure 2. The relationships among DBH, Height and FSV: (a) is the scatterplot of FSV and DBH; (b) is the scatterplot of FSV and Height.

Forests 2024 , 21 Figure 5 .
Figure 5. Schematic diagram of time series data reconstructed by S-G filter.

Figure 5 .
Figure 5. Schematic diagram of time series data reconstructed by S-G filter.

Figure 5 .
Figure 5. Schematic diagram of time series data reconstructed by S-G filter.

Figure 6 .
Figure 6.S-G filtering to reconstruct NDVI and extract phenology features.

Figure 6 .
Figure 6.S-G filtering to reconstruct NDVI and extract phenology features.

Forests 2024 ,
15, x FOR PEER REVIEW 8 of 21 was employed to differentiate the periodic data.The seasonally varying phenologic al features of larch were extracted by calculating the differences between features extracted from the current season and the previous season.Stepwise differencing i s particularly useful for analyzing or predicting time series data or signals as it st abilizes non-stable time series data by eliminating trends and seasonal component s.The detailed flow chart is shown in Figure 7.

Figure 7 .
Figure 7. Extracting phenological features using time series seasonal differences.

Figure 8 .
Figure 8. Fourier extraction of growing season amplitude variables.

Figure 9 .
Figure 9.The technical roadmap for mapping forest stock volume based on phenological feature data extraction.

Figure 9 .
Figure 9.The technical roadmap for mapping forest stock volume based on phenological feature data extraction.

Figure 10 .
Figure 10.Distribution of Pearson's correlation coefficients between different types of features and FSV: (a) Spectral features; (b) Textural features; (c-f) Phenological features by various methods.

Figure 10 .
Figure 10.Distribution of Pearson's correlation coefficients between different types of features and FSV: (a) Spectral features; (b) Textural features; (c-f) Phenological features by various methods.

Forests 2024 , 21 Figure 11 .
Figure 11.The scatterplots between measured and predicted FSV using Sentinel-2A images acquired in September.

Figure 11 .
Figure 11.The scatterplots between measured and predicted FSV using Sentinel-2A images acquired in September.

Figure 12 .
Figure 12.The scatterplots between the measured and predicted FSV using various phenological features.

Figure 12 .
Figure 12.The scatterplots between the measured and predicted FSV using various phenological features.Forests 2024, 15, x FOR PEER REVIEW 15 of 21

Figure 13 .
Figure 13.The map of FSV in the planted larch forest of study area.

Figure 13 .
Figure 13.The map of FSV in the planted larch forest of study area.

Figure 14 .
Figure 14.The histograms of accuracy indices derived from various subsets (a) the histograms of R2 derived from various subsets; (b) the histograms of rRMSE derived from various subsets.

Figure 14 .
Figure 14.The histograms of accuracy indices derived from various subsets (a) the histograms of R2 derived from various subsets; (b) the histograms of rRMSE derived from various subsets.

Funding:
This research was supported by the National Natural Science Foundation of China (Project number: 32171784); the Excellent Youth Project of the Scientific Research Foundation of the Hunan Provincial Department of Education (Project number: 21B0246).

Table 1 .
Sentinel-2A List of Vegetation Index Characteristic Variables.

Table 2 .
The accuracy indices of evaluating estimated FSV using spectral and textural features extracted from Sentinel-2 images acquired between April and September.

Table 3 .
The accuracy indices of evaluating estimated FSV after incorporating phenological features.