Estimation of Alpine Grassland Forage Nitrogen Coupled with Hyperspectral Characteristics during Di ﬀ erent Growth Periods on the Tibetan Plateau

: The applicability of hyperspectral remote sensing models for forage nitrogen (N) retrieval during di ﬀ erent growth periods is limited. This study aims to develop a multivariate model feasible for estimating the forage N for the growth periods (June to November) in an alpine grassland ecosystem. The random forest (RF) algorithm is employed to determine the optimum combinations of 38 spectral variables capable of capturing dynamic variations in forage N. The results show that (1) throughout the growth period, the red-edge ﬁrst shifts toward longer wavelengths and then shifts toward shorter wavelengths, the amplitude (AMP) and absorption depth (AD) gradually decrease, and the absorption position (AP) changes slightly; (2) the importance of spectral variables for forage N estimation di ﬀ ers during the di ﬀ erent growth periods; (3) the multivariate model achieves better results for the ﬁrst four periods (June to October) than for the last period (when the grass is completely senesced) (V-R 2 : 0.58–0.68 versus 0.23); and (4) for the whole growth period (June to November), the prediction accuracy of the general N estimation model validated by the unknown growth period is lower than that validated by the unknown location (V-R 2 is 0.28 and 0.55 for the validation strategies of Leave-Time-Out and Leave-Location-Out, respectively). This study demonstrates that the changes in the spectral features of the red wavelength (red-edge position, AMP and AD) are well coupled with the forage N content. Moreover, the development of a multivariate RF model for estimating alpine grasslands N content during di ﬀ erent growth periods is promising for the improvement of both the stability and accuracy of the model. areas and ﬁve ﬁeld observations from the end of June to the middle of November (GP170624–GP171115), this study ﬁnds that the changes between the forage N and spectral features (REP, AMP, and AD) are similar throughout the grass growth period. This result indicates that the forage N contents during di ﬀ erent growth periods can be estimated based on the aforementioned relationship. Furthermore, the importance of 38 spectral variables for estimating forage N were measured; results show that the applicability of spectral variables are di ﬀ erent during di ﬀ erent growth periods, indicating that a multivariate forage N estimation model for di ﬀ erent growth periods has the potential to improve the stability and accuracy of the results, especially under the background of the large spatial heterogeneity of alpine grasslands. This study provides insights towards enhancing the universality of forage N models and the sustainable utilization and management of natural alpine grassland resources. However, the result presents an unsatisfactory accuracy (V-R 2 = 0.23) when grass is completely senesced (GP171115), possible related to the low N contents and weak absorption features. This result also provides new insights for further developing a novel model for monitoring the forage N content in dead components of grasslands. The grass withering period lasts for ﬁve months (from November to March) in alpine grasslands, and the nutrition status of grass greatly a ﬀ ects the growth and reproduction of grazing livestock on the Tibetan Plateau during this period. In addition, four target-oriented validation strategies (LOO CV, LTO CV, LLO CV, and LLTO CV) are applied to assess the predictive performance of the general N model, which is established using all data sets. The results indicate that the di ﬀ erence of grassland vegetation at di ﬀ erent growth periods is the primary factor restricting the prediction accuracy of the model. However, it is necessary to verify the predictive ability of spatio-temporal model by using di ﬀ erent validation strategies in further studies.


Introduction
Nitrogen (N) is a key nutrient for vegetation growth and reproduction and plays an important role in the sustainability of natural grassland resources and animal husbandry on the Tibetan Plateau [1,2]. As a crucial component of proteins, nucleic acids, phospholipids, and chlorophyll in plants, N affects plant photosynthesis and water absorption [3][4][5]. Insufficient N supply could result in changes in the external morphology and internal metabolism of plants. According to these changes, corresponding agricultural practices should be implemented to ensure the N status of the plants [6]. The accurate and efficient monitoring of the spatial distribution characteristics and seasonal dynamics of forage N are sensitive to forage N based on RF algorithm; and (4) establish multivariate models to estimate forage N during different growth periods. According to this study, we expect to improve the prediction accuracy of the hyperspectral forage N inversion model and provide technical support for the accurate monitoring of vegetation growth.

Study Area
The study area (33.11 • N-35.57 • N, 100.77 • E-102.98 • E) lies in the northeastern margin of the Tibetan Plateau, and it has an average altitude of over 3000 m. As shown in Figure 1, the area includes three counties, from the northeast to southwest: Xiahe, Luqu, and Maqu. The area is characterized by a continental plateau climate with an average annual temperature of 1.6-13.6 • C, with an average annual precipitation of 400-800 mm, and it exhibits significant spatial differences and temporal variations. The study area is rich in natural alpine grassland resources and 87% of its land is covered with grassland (18,978.3 km 2 ); the major grassland types are alpine meadow and alpine steppe. The dominant species mainly include Stipa aliena, Festuca ovina, Poa pratensis var. pratensis, Kobresia capillifolia, and Potentilla chinensis. Fertilization and mowing are rarely carried out in the study area, and grassland utilization is mainly dominated by four-season rotational grazing and grazing blocked by fencing.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 19 are sensitive to forage N based on RF algorithm; and (4) establish multivariate models to estimate forage N during different growth periods. According to this study, we expect to improve the prediction accuracy of the hyperspectral forage N inversion model and provide technical support for the accurate monitoring of vegetation growth.

Study Area
The study area (33.11°N-35.57°N, 100.77°E-102.98°E) lies in the northeastern margin of the Tibetan Plateau, and it has an average altitude of over 3000 m. As shown in Figure 1, the area includes three counties, from the northeast to southwest: Xiahe, Luqu, and Maqu. The area is characterized by a continental plateau climate with an average annual temperature of 1.6-13.6 °C, with an average annual precipitation of 400-800 mm, and it exhibits significant spatial differences and temporal variations. The study area is rich in natural alpine grassland resources and 87% of its land is covered with grassland (18,978.3 km 2 ); the major grassland types are alpine meadow and alpine steppe. The dominant species mainly include Stipa aliena, Festuca ovina, Poa pratensis var. pratensis, Kobresia capillifolia, and Potentilla chinensis. Fertilization and mowing are rarely carried out in the study area, and grassland utilization is mainly dominated by four-season rotational grazing and grazing blocked by fencing.

Grassland Observation Data
Four sampling areas (GJ, YLJ, XC, and AZ) were set in the study area from north to south according to the spatial representation, accessibility and management mode, as shown in Figure 1. Among those sampling areas, AZ is located in the Azi yak propagation bases of Maqu County, XC is located in the Xicang Township of Luqu County, and YLJ and GJ are located in the Yaliji Township and Ganjia Township of Xiahe County, respectively. In comparison with other sampling areas, XC is nearest to YLJ, and the distance between them is approximately 30 km. In each area, five sample plots

Grassland Observation Data
Four sampling areas (GJ, YLJ, XC, and AZ) were set in the study area from north to south according to the spatial representation, accessibility and management mode, as shown in Figure 1. Among those sampling areas, AZ is located in the Azi yak propagation bases of Maqu County, XC is located in the Xicang Township of Luqu County, and YLJ and GJ are located in the Yaliji Township and Ganjia Township of Xiahe County, respectively. In comparison with other sampling areas, XC is nearest to YLJ, and the distance between them is approximately 30 km. In each area, five sample plots were selected as permanent observation sites. The distance between the plots is limited to approximately 3 km in consideration of the homogeneity of the plot, and the dimension of each plot is 100 m × 100 m. Five subplots (0.5 m × 0.5 m) were set up within each plot to obtain the plot variability.
During the grassland growth period (GP, June to November), five field investigations began on 24 June, 27 July, 30 August, 27 September, and 15 November 2017. Each fieldwork campaign lasted approximately four to six days; to present the specific time of each field work more intuitively, the beginning date of every investigation was used as a reference for naming the period, namely, GP170624, GP170727, GP170830, GP170927, and GP171115. GP-ALL is used to represent the whole growth period from June to November. A total of 100 sample were collected. In each subplot, the canopy reflectance spectra of the mixed grassland community were measured 10 times using a portable ASD Field Spec Pro FR2500 spectroradiometer (Analytical Spectral Devices Inc., Boulder, CO, USA) with spectral range of 350-2500 nm and a view angle of 25 • between 10:00 and 14:00 on a sunny day. The height of the sensor was approximately 1 m from the plant canopy. The reflectance spectra were then averaged as the final spectrum of each plot for the data. In addition, we collected the conventional observations for each subplot, including the geographic location, grassland community height and coverage, dominant species, dead component percent of grassland and AGB; grass samples were also collected.
After each investigation, all the grass samples were transferred to the lab for further physical processing (oven-dried at 65 • C for 48 h, smashed and sieved) and chemical analysis. After the samples were boiled and digested with concentrated H 2 SO 4 solution, dehydrated and carbonized, and a series of oxidation reactions were performed. A FIAstar 5000 flow injection analyzer was used to quantify the N content. The method uses K 2 SO 4 /CuSO 4 as a catalyst, and the operation steps are simple and fast and do not interfere with the quantification of N elements. Table 1 summarizes the 38 spectral variables, including 20 VIs, 10 absorption bands, 4 red-edge parameters, and 4 absorption features, which are widely used for estimating forage N. This study aims to assess the performance of all these variables in the estimation of N during different growth periods of grassland.

Spectral Variables
Narrow-band VIs has been widely used to qualitatively and quantitatively evaluate the growth status of grasslands. The specific known protein, chlorophyll, and N absorption bands have been successfully used for the estimation of grass N [14,16,19]. The red-edge parameters, such as REP, amplitude (AMP), Slope725 and Slope_mean, have shown to be significantly correlated with the N of vegetation [31]. The red-edge parameters were calculated based on the first derivative spectrum to effectively reduce the interference of soil and the atmospheric background [32]. Studies have also found that the absorption features of the absorption position (AP), absorption depth (AD), normalized band depth index (NDBI) and band depth ratio (BDR) in the red-edge range (550-750 nm) based on continuum-removed spectra reflectance can effectively extract the biochemical parameters of grassland vegetation, especially N and chlorophyll [8,9,33].

Red-edge parameters
Red-edge position (REP) AMP Wavelength of the red-edge peak (maximum slope position) First derivative value at the red-edge peak (maximum slope) Slope725 First derivative value at 725 nm Slope_mean First derivative value obtained from the corresponding mean red-edge position

RF Algorithm
The nonparametric and multivariate RF algorithm is used to establish a forage N estimation model from different spectral variables during different growth periods. RF can improve the performance of classification and regression trees (CARTs) through the use of a multitude of decision trees. The algorithm was developed by and is commonly used to solve complicated multiple regression problems [27]. The advantages of this algorithm include being less prone to overfitting, handling highly nonlinear data, exhibiting strong anti-noise ability, showing high accuracy, and having relatively simple implementation [54,55]. When using the RF algorithm, it is necessary to optimize three pivotal parameters, mtry (number of predictors tested at each node), ntree (number of regression trees), and nodesize (minimal size of the terminal nodes of the trees). In this study, the RF program is processed using MATLAB 2016a software.

Variable Selection
For different growth periods, the sequential backward search (SBS) method, based on the RF algorithm, is used to determine feature variables. SBS relies on the importance of each spectral variable, which is calculated by measuring the percent increase of the root mean squared error (RMSE) when the out-of-bag (OOB) data of each variable are permuted while all others remain unchanged. We rank the variables (importance not equal to 0) by importance (from large to small) and then evaluate the prediction ability of each RF model, which was established by the first n variables. Finally, we choose the feature set with the least number of variables and the optimal prediction ability as a result of feature selection. For the whole growth period, the forward feature selection (FFS) algorithm that works in conjunction with target-oriented validation is used to determine variables to improve the generalization ability of the model. The algorithm first tunes and trains the RF models using all possible combinations of any two spectral variables. The best initial model in view to target oriented performance is kept. Then, the number of variables is iteratively increased. The improvement of the model is tested for each additional predictor using target-oriented cross-validation (CV). The process stops when none of the remaining variables decreases the error of the current best model. Finally, the predictor variables with best modeling performance is determined as the feature variables. The details of this algorithm can be found in Meyer et al. [56].

Validation Strategies
When establishing the RF model of forage N during different growth periods using the corresponding feature variables, we iterate the two core parameters (ntree and mtry) and choose the default values for the other parameters (e.g., nodesize, minimum impurity split, minimum samples leaves, and minimum samples splits). The ntree values increase from 10 to 1000 in intervals of 10 for each iteration, with a total of 100 iterations, and mtry is increased from 1 to m (m is the number of variables) in intervals of 1 each time for a total of m iterations.
To assess the performance and stability of all models (i.e., models of different growth periods and whole growth period) in forage N estimation, Leave-One-Out (LOO) CV is used in this study because of its merits of unbiased estimation and outlier detection [57]. Moreover, to evaluate the prediction performance of the general model on unknown locations (sample area) and unknown points in time (growth period), three "target-oriented" validation strategies are used is this study: Leave-Time-Out (LTO) CV, Leave-Location-Out (LLO) CV and Leave-Location-and-Time-Out (LLTO) CV [56]. The goodness of fit of the measured and simulated forage N is evaluated by the coefficient of determination (R 2 ), mean absolute error (MAE), RMSE, and coefficient of variation in the RMSE (CVRMSE) (Equation (1)).
where y i and y i are the N values of the test set for the simulated and measured; y is the measured mean value; and n equals the sample size of the test dataset. If the CVRMSE < 10%, the predictive capability is considered excellent; if the CVRMSE ranges from 10-20%, the predictive power is good; if the CVRMSE > 20%, the predictive performance is unsatisfactory.

Variations in Forage N Contents and Reflectance Spectra during Different Growth Periods
The highest N content was observed at site GJ and the lowest at site YLJ during the different growth periods, as shown in Table 2. The overall forage N contents of different sampling areas during different growth periods show a clear decreasing trend with plant growth (P < 0.05). The forage N content in AZ is significantly different in the different growth periods (P < 0.05), indicating that there is significant spatial heterogeneity of grassland vegetation in this sampling area. Moreover, significant differences are not observed in the N contents in YLJ and GJ among GP170727, GP170830, and GP170927. In addition, the mean forage N content of the four sampling areas ranges from 0.89% to 2.16% throughout the growth period and presents a coefficient of variation of 7.71-15.81%. In particular, the coefficient of variation of GP170927 and GP171115 (15.81% and 15.22%) are much greater than those of the other periods, indicating that a large spatial difference occurs in the forage N content in the last two periods. From the typical photographs of samples during different growth periods, some information about the status of the grass can be obtained (i.e., color, morphology, and blossom).
The changes in canopy spectral reflectance at the four sampling areas during different growth periods, as shown in Figure 2. Overall, from the early period of forage growth (GP170624) to the vigorous growth period (GP170830), the plant mainly undergoes vegetative and reproductive growth. The vegetative stage refers to the developmental period comprising leaf growth and development, and the grass grows rapidly and the vegetation coverage increases gradually. The absorption of the visible region by the grass canopy is obviously enhanced, and the reflectance in the visible region decreases gradually. In addition, the absorption in the red valley between 650 and 710 nm becomes more obvious, and the reflectance in the near-infrared bands also increases significantly due to multiple scattering. After the vigorous grass growth period, the vegetation coverage decreases gradually and the forage stops growing and begins to senesce. As the grass gradually senesces, the absorption of the visible region by the grass gradually weakens, the reflectance in the visible region increases while the near-infrared region decreases significantly (GP170830 to GP170927). When the grass is completely senesced in November (GP171115), the spectrum is similar to the spectral characteristics of the soil, showing a slowly increasing trend within the wavelengths from 350-1350 nm without obvious absorption or reflection features. In addition, the higher shortwave near-infrared (1420-1800 nm and 1930-2300 nm) reflectance might be largely due to less water absorption in senesced grass. The absorption features on the spectral curve (i.e., obvious absorption valleys and reflection peaks) are gradually weakened with the advancement of the growth period, which may be related to the biochemical parameters of the forage.      Typical photographs of samples

Spectral Absorption Features and Red-Edge Shift
The changes rule of the spectrum absorption features (AP and AD) and red-edge parameters (REP and AMP) throughout the alpine grassland growth period (June to November) are shown in Figure 3. Overall, the AD of the red region (550-750 nm) in the different sampling areas present a gradually decreasing trend with the advancement of the growth period. The AP (between 676 and 679 nm) of the different sampling areas and growth periods are not significantly different. The AD of the AZ, XC and YLJ change little from GP170624 to GP170830 but significantly change from GP170830 to GP171115. Especially in GP171115 (when the grass is completely senesced), the absorption features become very weak and the AP simultaneously tends to be stable at 677 nm. In addition, the AD in GJ fluctuates slightly (first decreases and then increases) in the first three periods (GP17064-GP170830), which may be related to the large spatial heterogeneity of the sampling area, as shown in Figure 3. According to Figure 3, the REP first shifts toward longer wavelengths and then shifts toward shorter wavelengths with the advancement of the growth period as a whole. Similar to the trend of

Spectral Absorption Features and Red-Edge Shift
The changes rule of the spectrum absorption features (AP and AD) and red-edge parameters (REP and AMP) throughout the alpine grassland growth period (June to November) are shown in Figure 3. Overall, the AD of the red region (550-750 nm) in the different sampling areas present a gradually decreasing trend with the advancement of the growth period. The AP (between 676 and 679 nm) of the different sampling areas and growth periods are not significantly different. The AD of the AZ, XC and YLJ change little from GP170624 to GP170830 but significantly change from GP170830 to GP171115. Especially in GP171115 (when the grass is completely senesced), the absorption features become very weak and the AP simultaneously tends to be stable at 677 nm. In addition, the AD in GJ fluctuates slightly (first decreases and then increases) in the first three periods (GP17064-GP170830), which may be related to the large spatial heterogeneity of the sampling area, as shown in Figure 3.

Spectral Absorption Features and Red-Edge Shift
The changes rule of the spectrum absorption features (AP and AD) and red-edge parameters (REP and AMP) throughout the alpine grassland growth period (June to November) are shown in Figure 3. Overall, the AD of the red region (550-750 nm) in the different sampling areas present a gradually decreasing trend with the advancement of the growth period. The AP (between 676 and 679 nm) of the different sampling areas and growth periods are not significantly different. The AD of the AZ, XC and YLJ change little from GP170624 to GP170830 but significantly change from GP170830 to GP171115. Especially in GP171115 (when the grass is completely senesced), the absorption features become very weak and the AP simultaneously tends to be stable at 677 nm. In addition, the AD in GJ fluctuates slightly (first decreases and then increases) in the first three periods (GP17064-GP170830), which may be related to the large spatial heterogeneity of the sampling area, as shown in Figure 3. According to Figure 3, the REP first shifts toward longer wavelengths and then shifts toward shorter wavelengths with the advancement of the growth period as a whole. Similar to the trend of According to Figure 3, the REP first shifts toward longer wavelengths and then shifts toward shorter wavelengths with the advancement of the growth period as a whole. Similar to the trend of the variation in AD, the REP of the first three periods (GP170624-GP170830) changes slightly; however, it changes significantly in the last three periods (GP170830-GP171115), which may be closely related to the chlorophyll concentration of the forage. Because the chlorophyll has a close relationship with REP, when the grass gradually senesces in the latter three periods, chlorophyll is gradually decomposed and decreases significantly, resulting in significant REP changes. Moreover, the average REP of the different sampling areas in the GP171115 is 694 nm, which may contribute to identifying the growth period of alpine grassland.
Overall, as the N content decreases, the REP, AMP, and AD present gradually decreasing variation trends throughout the growth period, as shown in Figure 4. Due to the impacts of spatial heterogeneity between different areas, the changes in red region absorption features and red-edge parameters inevitably exhibit a certain degree of fluctuation with changes in the growth period. The above findings indicate that dynamic changes in the spectral characteristics within the red region during different growth periods may be related to the forage N. the variation in AD, the REP of the first three periods (GP170624-GP170830) changes slightly; however, it changes significantly in the last three periods (GP170830-GP171115), which may be closely related to the chlorophyll concentration of the forage. Because the chlorophyll has a close relationship with REP, when the grass gradually senesces in the latter three periods, chlorophyll is gradually decomposed and decreases significantly, resulting in significant REP changes. Moreover, the average REP of the different sampling areas in the GP171115 is 694 nm, which may contribute to identifying the growth period of alpine grassland. Overall, as the N content decreases, the REP, AMP, and AD present gradually decreasing variation trends throughout the growth period, as shown in Figure 4. Due to the impacts of spatial heterogeneity between different areas, the changes in red region absorption features and red-edge parameters inevitably exhibit a certain degree of fluctuation with changes in the growth period. The above findings indicate that dynamic changes in the spectral characteristics within the red region during different growth periods may be related to the forage N.

Relationship between the Forage N Content and Different Variables
Variable importance based on the RF algorithm is applied to rank the importance of the 38 spectral variables. The measured importance of the spectral variables for estimating forage N during different growth periods, as shown in Figure 5. The spectral variables for forage N estimation clearly differs during different growth periods. In GP170727 and GP170830, REP is the most important

Relationship between the Forage N Content and Different Variables
Variable importance based on the RF algorithm is applied to rank the importance of the 38 spectral variables. The measured importance of the spectral variables for estimating forage N during different growth periods, as shown in Figure 5. The spectral variables for forage N estimation clearly differs during different growth periods. In GP170727 and GP170830, REP is the most important variable. In GP170624, NDGI, NDNI, PRI1, and NDCI are the four most important variables for estimating forage N, and the absorption bands, red-edge parameters, and absorption features are not important variables. In GP170727, the Slope_mean and R2060 are the two most important variables for N estimation. In GP171115, the two red region absorption features (NBDI and BDR) are important for N estimation. These findings further indicate that a N model should be developed with a combination of different spectral variables for alpine grassland during different growth periods because of the applicability of different spectral variables during different growth periods.
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 19 variable. In GP170624, NDGI, NDNI, PRI1, and NDCI are the four most important variables for estimating forage N, and the absorption bands, red-edge parameters, and absorption features are not important variables. In GP170727, the Slope_mean and R2060 are the two most important variables for N estimation. In GP171115, the two red region absorption features (NBDI and BDR) are important for N estimation. These findings further indicate that a N model should be developed with a combination of different spectral variables for alpine grassland during different growth periods because of the applicability of different spectral variables during different growth periods.

Estimation Model for Forage N during Different Growth Periods
In this study, the SBS algorithm and V-R 2 are employed to determine the optimal variable combinations that have the least number of variables and the best modeling performance for estimating forage N during different growth periods, as shown in Table 3. From GP170624 to GP171115, the number of selected variables ranges from four to six, thus accounting for 10.5-15.8% of all variables, indicating that the variable selection method based on the RF algorithm yields the desired result (the information contained in a few variables can represent the majority of variables). To further evaluate the predictive ability of the N estimation model at different growth periods, we use LOO CV to assess the accuracy of the model, as shown in Table 3. The results show that the N model for the first four periods (GP170624-GP170927) can explain 58-68% of the variance in the forage N content, except for the last period (GP171115), which is only 23%, as shown in Table 3 and Figure 6. The N model yields higher accuracy in GP170624 (V-R 2 = 0.68, V-RMSE = 0.2046%, 9.46% of the mean) and GP170830 (V-R 2 = 0.67, V-RMSE = 0.2202%, 12.78% of the mean) than in the other periods. The accuracy in GP171115 may be related to the low N content and weak absorption features during the grass senescent period. Overall, the CVRMSE of the model during different growth periods are all less than 20%, indicating that these N models have good predictive capability.

Estimation Model for Forage N during Different Growth Periods
In this study, the SBS algorithm and V-R 2 are employed to determine the optimal variable combinations that have the least number of variables and the best modeling performance for estimating forage N during different growth periods, as shown in Table 3. From GP170624 to GP171115, the number of selected variables ranges from four to six, thus accounting for 10.5-15.8% of all variables, indicating that the variable selection method based on the RF algorithm yields the desired result (the information contained in a few variables can represent the majority of variables). To further evaluate the predictive ability of the N estimation model at different growth periods, we use LOO CV to assess the accuracy of the model, as shown in Table 3. The results show that the N model for the first four periods (GP170624-GP170927) can explain 58-68% of the variance in the forage N content, except for the last period (GP171115), which is only 23%, as shown in Table 3 and Figure 6. The N model yields higher accuracy in GP170624 (V-R 2 = 0.68, V-RMSE = 0.2046%, 9.46% of the mean) and GP170830 (V-R 2 = 0.67, V-RMSE = 0.2202%, 12.78% of the mean) than in the other periods. The accuracy in GP171115 may be related to the low N content and weak absorption features during the grass senescent period. Overall, the CVRMSE of the model during different growth periods are all less than 20%, indicating that these N models have good predictive capability.

Estimation Model for Forage N throughout the Growth Periods
To further develop a N estimation model applicable for each growth period, the FFS and four CV strategies (LOO, LTO, LLO, and LLTO) are employed to select feature variables and evaluate the predictive ability of the model to data of unknown growth periods and unknown sample areas, as shown in Table 4. For the whole growth period (entire datasets), the NDVI, PRI, REP, Slope_mean, and R640 are determined to detect variations in forage N. Among the four validation strategies, LOO CV shows a good performance between the measured and simulated data (V-R 2 = 0.51, V-RMSE = 0.3741). In addition, the ability of the model to predict the N of the unknown locations (LLO CV) within the four sample areas remains high (V-R 2 = 0.55). However, for an unknown growth period (LTO CV) or faced with unknown location and period simultaneously (LLTO CV), the prediction ability of the model is lower than the other validation strategies (V-R 2 is 0.28 and 0.29 for LTO CV and LLTO CV, respectively). These findings indicate that the model based on the above five variables and the entire datasets can overcome the spatial differences between different sample areas in the study area; however, predicting the temporal changes of forage N in the unknown growth period is difficult. A significant reason may be that the forage N content varies greatly on a temporal scale (different growth period). The spectral variables available in this study cannot easily capture the variations in N at each growth period, which indicates that the sensitivity of spectral variables to N are not consistent in different growth periods.

Estimation Model for Forage N throughout the Growth Periods
To further develop a N estimation model applicable for each growth period, the FFS and four CV strategies (LOO, LTO, LLO, and LLTO) are employed to select feature variables and evaluate the predictive ability of the model to data of unknown growth periods and unknown sample areas, as shown in Table 4. For the whole growth period (entire datasets), the NDVI, PRI, REP, Slope_mean, and R640 are determined to detect variations in forage N. Among the four validation strategies, LOO CV shows a good performance between the measured and simulated data (V-R 2 = 0.51, V-RMSE = 0.3741). In addition, the ability of the model to predict the N of the unknown locations (LLO CV) within the four sample areas remains high (V-R 2 = 0.55). However, for an unknown growth period (LTO CV) or faced with unknown location and period simultaneously (LLTO CV), the prediction ability of the model is lower than the other validation strategies (V-R 2 is 0.28 and 0.29 for LTO CV and LLTO CV, respectively). These findings indicate that the model based on the above five variables and the entire datasets can overcome the spatial differences between different sample areas in the study area; however, predicting the temporal changes of forage N in the unknown growth period is difficult. A significant reason may be that the forage N content varies greatly on a temporal scale (different growth period). The spectral variables available in this study cannot easily capture the variations in N at each growth period, which indicates that the sensitivity of spectral variables to N are not consistent in different growth periods.

Variation in Forage N Contents during Different Growth Periods
This study shows that the forage N content presents a gradual decrease with the advancement of the growth period in the alpine grassland (from the end of June to the middle of November), as shown in Table 2. Especially in the latter three periods (GP170830-GP171115), the forage N decreases rapidly. The N content is the highest in GP170624 (2.16%), when it is approximately 2.5 times higher than that in GP171115. Many previous studies have also found that the crude protein of grass in alpine grasslands on the Tibetan Plateau gradually decreases with plant growth, and the content is the highest in June [58][59][60]. These findings also indirectly confirming our result about the variation trend in forage N.
In the present study, the reason for the changes in forage N can be explained by dividing the whole growth period into two stages. The first stage is from the end of June to the end of August (GP170624-GP170830), when grass mainly undergoes vegetative and reproductive growth. As dry matter accumulates in the plant, the biomass also increases, the mechanical tissues of plants grow, and the proportion of physiologically active, nonmechanical tissues declines gradually [61,62]. Because there are large amounts of crude fiber, such as lignin and cellulose but little N in mechanical tissues, crude protein eventually shows a decreasing trend as the forage quantity increase [63]. The second stage is from the end of August to the middle of November (GP170830-GP171115); when the temperature decreases, the frost-free period becomes shorter, and the forage stops growing and begins to senesce. Meanwhile, the resorption of plants has weakened, resulting in the continuous decrease in forage N content [64,65].

Spectral Characteristics of the Forage Canopy during Different Growth Periods
The spectral reflectance of grassland vegetation is characterized by the canopy morphology and structure, the internal structure of leaves, and the compounds of different components, which are affected by various factors such as water content, chlorophyll content, vegetation coverage, growth period, soil condition, and atmospheric condition [23,24]. In the visible region, the chlorophyll content in leaves is the most important factor that affects the characteristics of leaf reflectance [66][67][68]. Especially in the red-edge region (680-750 nm), the spectral features are closely related to the biochemical parameters, such as N, protein, and chlorophyll [21,31,69,70]. This study has shown that, as the growth period advances, the red-edge first shifts toward longer wavelengths and then toward shorter wavelengths; AMP and AD gradually decrease, and AP changes slightly (most at 678 nm), as shown in Figure 3. In a related study, Zhang et al. also indicated that the REP first shifts to the near-infrared direction then to the green direction from the late stage of returning to green to the vigorous growth period and then to a yellowing and withering period of vegetation [71,72]. The general variation trends of AD and AMP are similar to those observed in this study.
The regular changes in these spectral feature parameters during the growth period can also reflect the growth and development of grassland, which may be related to the morphological changes in vegetation and some physicochemical parameters. For instance, the REP and AMP has a significant relationship with the N content, and the red-edge shift towards longer wavelengths when the N supply increases [31]. For the entire growth period (GP170624-GP171115), the changes of the forage N are similar to that of the spectral features (REP, AMP, and AD) as a whole, as shown in Figure 4; this suggests that the dynamic changes in the spectral feature parameters during different growth periods in the red region may be coupled with the variations in N content. Furthermore, the results of this study also indicated that those features (REP, Slope725, Slope_mean, NBDI, etc.) contribute to the estimation of forage N during different growth periods, as shown in Tables 3 and 4.

Applicability of Spectral Variables for Estimating Forage N during Different Growth Periods
The study shows that the importance of the spectral parameters differs during different growth periods, as shown in Figure 5, and the general model of forage N developed with all available datasets has a lower estimation accuracy (LOO CV, V-R 2 : 0.51 versus V-R 2 : 0.58-0.68) compared with the single growth period (GP170624-GP170927), as shown in Figure 6 and in Tables 3 and 4. The reason for this trend may be that the morphological structure and physicochemical characteristics of the plants change significantly as the growth and development of vegetation; thus, the applicability of different spectral parameters will inevitably differ during different growth periods. According to our results, the N model yields higher accuracy in GP170624 (V-R 2 = 0.68) and GP170830 (V-R 2 = 0.67) than in the other periods. In comparison with the results for estimating the concentration of N by using the reflectance of a tropical grass (Cenchrus ciliaris) (R 2 = 0.73) [31], the estimation results of our model are satisfactory. However, most studies only detect forage N content during a particular period without considering changes in different growth periods. This study attempts to estimate forage N content in alpine grassland during different growth stages for covering the shortcoming due to using a single period.
As shown in Figure 6, the model overestimates N at lower levels and underestimates it at higher ones. This may be attributable to the prediction results of RF model, which are determined by multiple decision trees through voting, and the simulation results tend to be the average of the training data. We analyzed the reasons for model results overestimating and underestimating in this study, and found that (1) the small sample size is one of the fundamental reasons for this phenomenon, and (2) the parameters of the RF model and small correlations between predictors also affect overestimation or underestimation of the model. The overestimation is due to the fact that the OOB observations used to derive predictions from the trees might not be representative. Using stratified subsampling for both tuning parameter selection and error estimation in random forests might therefore be a solution to reduce the bias in the estimated value [73]. In addition, the estimation result of the N model shows a large bias in the last period (GP171115), and the phenomenon of overestimation and underestimation is more obvious. This may be caused by the low N content and vegetation coverage of grassland during this period when most variables are less sensitive to forage N. Moreover, bare soil also has a significant influence on N estimates.
During vegetative growth and reproductive growth, the chlorophyll in plants increases and organic matter accumulates continuously. Some spectral variables (i.e., NDNI, VOG, PRI, SIPI, and REP) that are sensitive to chlorophyll and N play a significant role in the estimation of forage N [31,42,74,75], and the N model in the corresponding period also exhibits acceptable estimation accuracy, as shown in Table 3. In addition, when the leaf area index (LAI) of vegetation is very high (i.e., the vegetation is very dense), the sensitivity of SR and NDVI to N will be decreased [31,76,77]. As the plant undergoes senescence, chlorophyll decomposes gradually and the coverage and canopy density decrease rapidly; accordingly, variables that were previously sensitive to N and chlorophyll previously may become insensitive. However, some spectral variables, such as SAVI, OSAVI, BDR, and NBDI, may contribute significantly to the detection of N during this period because SAVI and OSAVI have the potential to effectively eliminate the soil background effect [47,78]; BDR and NBDI are calculated from the continuum-removed spectrum in which the absorption features of forage nutrients are strengthened [8,26]. These results support the findings of this study that the NBDI, OSAVI, and BDR are three of the first five most important variables for N estimation in GP170927, as shown in Table 3. In particular, when the grassland is completely senesced (GP171115), the majority of the spectral variables are no longer sensitive to the weak absorption features of N, which directly results in the unsatisfactory accuracy of the forage N model (V-R 2 = 0.23).
According to the preliminary results of this study, multivariate models should be developed for the estimation of forage biochemical parameters during different growth periods in the alpine grasslands. However, developing a spectral variable suitable for estimating forage N in alpine grasslands over the whole growth period is challenging. In subsequent studies, it will be meaningful to focus on estimating the forage N during the grass senescing period while further optimizing the multivariate models during different growth periods.

Target-Oriented Validation
Our study shows that the predictive performance of the N estimation model based on an unknown growth period is lower than that based on an unknown sample area (V-R 2 = 0.28 of LTO CV vs. V-R 2 = 0.55 of LLO CV), as shown in Table 4, which indicated that the difference in grassland vegetation at different growth periods is the primary factor restricting the prediction accuracy of the model. The reasons for this finding may related to (1) the considerable variations in the vegetation canopy morphology, biochemical substance concentration, coverage, and height, as well as soil background during different growth periods; and (2) limitations of the available sample size in this study. An ideal RF model relies on the analysis of large amounts of data to establish robust models [27]. Most validation strategies for spatial and temporal models in previous studies are based on random k-fold CV or a random validation subset of the entire dataset [79,80]. However, some studies have concluded that the spatio-temporal machine learning models are prone to temporal or spatial over-fitting and need to be evaluated by LLO CV and LTO CV [56,81]. This study attempted to develop a general forage N estimation model suitable for each growth period using a target-oriented FFS algorithm and validation strategies; however, a low predictive ability was observed for the model based on LTO CV. For this, subsequent studies will need to integrate more sample data and a wider variety of spectral variables to further explore this issue to develop a more practical and applicable forage N estimation model across the whole growth period of alpine grassland.

Conclusions
On the basis of the 20 permanent sample plots in four sampling areas and five field observations from the end of June to the middle of November (GP170624-GP171115), this study finds that the changes between the forage N and spectral features (REP, AMP, and AD) are similar throughout the grass growth period. This result indicates that the forage N contents during different growth periods can be estimated based on the aforementioned relationship. Furthermore, the importance of 38 spectral variables for estimating forage N were measured; results show that the applicability of spectral variables are different during different growth periods, indicating that a multivariate forage N estimation model for different growth periods has the potential to improve the stability and accuracy of the results, especially under the background of the large spatial heterogeneity of alpine grasslands.
This study provides insights towards enhancing the universality of forage N models and the sustainable utilization and management of natural alpine grassland resources. However, the result presents an unsatisfactory accuracy (V-R 2 = 0.23) when grass is completely senesced (GP171115), possible related to the low N contents and weak absorption features. This result also provides new insights for further developing a novel model for monitoring the forage N content in dead components of grasslands. The grass withering period lasts for five months (from November to March) in alpine grasslands, and the nutrition status of grass greatly affects the growth and reproduction of grazing livestock on the Tibetan Plateau during this period. In addition, four target-oriented validation strategies (LOO CV, LTO CV, LLO CV, and LLTO CV) are applied to assess the predictive performance of the general N model, which is established using all data sets. The results indicate that the difference of grassland vegetation at different growth periods is the primary factor restricting the prediction accuracy of the model. However, it is necessary to verify the predictive ability of spatio-temporal model by using different validation strategies in further studies.