Mapping the Forage Nitrogen-Phosphorus Ratio Based on Sentinel-2 MSI Data and a Random Forest Algorithm in an Alpine Grassland Ecosystem of the Tibetan Plateau

Abstract

Nondestructive and accurate estimating of the forage nitrogen–phosphorus (N:P) ratio is conducive to the real-time diagnosis of nutrient limitation and the formulation of a management scheme during the growth and development of forage. New-generation high-resolution remote sensors equipped with strategic red-edge wavebands offer opportunities and challenges for estimating and mapping forage N:P ratio in support of the sustainable utilization of alpine grassland resources. This study aims to detect the forage N:P ratio as an ecological indicator of grassland nutrient content by employing Sentinel-2 multispectral instrument (MSI) data and a random forest (RF) algorithm. The results showed that the estimation accuracy (R²) of the forage N:P ratio model established by combining the optimized spectral bands and vegetation indices (VIs) is 0.49 and 0.59 in the vigorous growth period (July) and the senescing period (November) of forage, respectively. Moreover, Sentinel-2 MSI B9 and B12 bands contributed greatly to the estimation of the forage N:P ratio, and the VIs (RECI2) constructed by B5 and B8A bands performed well in the estimation of the forage N:P ratio. Overall, it is promising to map the spatial distribution of the forage N:P ratio in alpine grassland using Sentinel-2 MSI data at regional scales. This study will be potentially beneficial in implementing precise positioning of vegetation nutrient deficiency and scientific fertilization management of grassland.

1. Introduction

The biochemical parameters of forage in natural alpine grasslands, such as nitrogen (N), phosphorus (P) and chlorophyll, are important aspects of forage nutrition conditions and feed value and play an important role in the life cycle of plants [1,2]. As an essential component of protein, N promotes the development and formation of stems, leaves and fruits and therefore is a crucial nutrient element affecting grassland biomass [3,4]. P promotes the metabolism of individual plants and is an essential nutrient for grassland plants living in the harsh environment of the Tibetan Plateau [5,6]. The plant nitrogen–phosphorus (N:P) ratio is a pivotal ecological indicator representing the relative absorption of N and P that can reflect the growth status of vegetation [7,8,9].

Despite considerable research [7,10,11], a critical value for the N:P ratio of forage is still not a unified standard in different ecosystems. Overall, the division scheme of threshold developed by Tessier and Raynal can be applied to natural alpine grassland ecosystems [10]. In this scheme, an N:P ratio < 14 is likely plant N limitation, an N:P ratio > 16 tends to indicate P deficiency, and an N:P ratio ranging from 14 to 16 indicates that the growth of vegetation is limited by both N and P [10]. The N:P ratio not only is related to species richness, functional traits and productivity but also indicates potential limitations of N and P [12,13]. Therefore, investigating the spatial distribution characteristics of the N:P ratio at a regional scale can provide useful information for the rational utilization of alpine grassland resources and the diagnosis of forage growth status.

Remote sensing technology has provided a more convenient and efficient approach for the dynamic monitoring of grassland vegetation at the landscape and regional scales [14,15]. The methods of estimating forage biochemical parameters can be mainly divided into two types, viz. the physically based method and empirical method [16]. Radiative transfer models (RTMs) explain the interaction of solar radiation in plant leaves according to the optical laws, thus offer many advantages in the generalizability and robustness. The empirical method offers a preliminary understanding of the relationship between forage biochemical parameters and spectral variables (i.e., spectral features and vegetation indices (VIs)), and is relatively easy to perform [17]. Studies on the estimation of forage N content and spatial inversion based on multisource high-resolution spectral data have become increasingly mature. Forage N have been shown to correlate strongly with chlorophyll [18]. The detection of forage N using remote sensing technology is generally associated with the absorption of chlorophyll, and usually includes the spectral regions related to chlorophyll detection (i.e., the red-edge and near-infrared (NIR) regions). N is involved in the protein synthesis, which facilitate the photosynthetic process [19]. The red-edge position is the point where the vegetation reflectance spectrum has the maximum slope in the 680–750 nm [20]. Studies have indicated that the red-edge is less sensitive to atmospheric effects and soil background. The red-edge position has been extensively applied in estimating leaf area index, N, chlorophyll and biomass [21,22]. However, due to the few identifiable spectral absorption features and low content of P in grassland plants, estimating the P content using high-resolution remote sensing data is still a challenging problem. Since the orders of magnitude of N and P are not consistent, error propagation may occur if the forage N and P contents are estimated separately and then calculating the N:P ratio [23,24]. Therefore, to minimize these errors, it is essential to estimate the forage N:P ratio directly using high-resolution remote sensing data. Ramoelo et al. showed that directly retrieving the forage N:P ratio using hyperspectral data was feasible (R² of 0.69–0.85), and the spectral bands were mainly distributed in the shortwave infrared (SWIR) regions are highly sensitive to the leaf N:P ratio [11]. Gókkaya et al. estimated the N:P ratio of forests using VIs calculated by spaceborne (Hyperion EO-1) and airborne (CASI) imaging spectrometers, and the results showed that some VIs were significantly related to the N:P ratio (R² of 0.34–0.70) [25]. In addition, the blue region is very important for the detection of the canopy N:P ratio [26]. Although field portable spectroradiometers and hyperspectral sensors can obtain rich spectral information, spectroradiometers are not easily accessible to a large number of observations, and hyperspectral images are expensive to acquire and have a small swath width [14,27]. Hence, it is difficult to apply such instruments to monitor grassland biochemical parameters continuously at a regional scale.

Classic multispectral satellite data (i.e., MODIS, Landsat, and SPOT data) are usually restricted to the macro-level monitoring of grassland biomass, coverage, and net primary productivity due to their limited spectral channels and discontinuous broadbands [28,29,30]. Some new-generation satellites, such as RapidEye (launched on 29 August 2008), WorldView-2 (launched on 6 October 2009), Sentinel-2A/B (launched on 23 June 2015 (S2A), and 7 May 2017 (S2B)) and Gaofen-6 (launched on 2 June 2018), equipped with sensors containing specialized red-edge wavebands are more sensitive to biochemical parameters of vegetation. These satellites with a higher spatial resolution, wider spatial coverage, and improved revisit periods have irreplaceable advantages in dynamic monitoring of plant growth status and have been widely applied to the quantitative estimation of grassland vegetation physical and chemical parameters [4,31,32]. In particular, Sentinel-2 multispectral instrument (MSI) data are free and open-source, with a large number of archived data available, and are therefore easier to obtain than other satellite data. From the perspectives of cost and feasibility, these data are more suitable than others for the diagnosis and monitoring of vegetation nutrient limitation at a regional scale. Several studies have used Sentinel-2 data to estimate the chlorophyll, N, and biomass of vegetation in the local area [32,33,34], but these studies mainly focused on the biochemical parameters of crops (i.e., coffee and maize). To the best of our knowledge, Sentinel-2 data have not yet been used to estimate and investigate the forage N:P ratio in the alpine grassland ecosystems.

Machine learning algorithms, especially random forest (RF), have been widely used to estimate the biochemical parameters of vegetation. This is because the algorithms not only have excellent performance to overcome the problems of nonlinearity and multicollinearity, but can also be used to explore the internal relationships of specific vegetation biochemical parameters with multiple satellite wavebands and VIs [4,35].

Here a study applying Sentinel-2 MSI data and an RF algorithm combined with the field measurement data from an alpine grassland during the vigorous growth period and the senescence stage is taken up with the following objectives: (1) to explore the potential of Sentinel-2 spectral bands and VIs to detect the forage N:P ratio, (2) to optimize spectral bands and VIs and establish an estimation model of the forage N:P ratio at different growth stages, and (3) to map the spatial distribution of the forage N:P ratio in the study area and analyze the nutrient limitations of the forage.

2. Material and Methods

2.1. Study Area

Three counties (Maqu, Luqu and Xiahe) and one city (Hezuo) in Gansu Province, China, are designated as the study area (100.77° E–103.11° E, 33.11° N–35.57° N), which is located in the mixed grassland ecoregion on the eastern Tibetan Plateau (Figure 1a). Grasslands are an important dominant resource in the region, occupying approximately 85% of the area. The study area belongs to an alpine grassland region with a typical plateau continental climate, and it has an average altitude of 3000 m. The total precipitation and annual mean temperature of this area are 400–800 mm and 1.6–13.6 °C, respectively, and rain and heat occur over the same period. Affected by monsoons and topography, the spatial and temporal distribution of precipitation in the area are uneven, with the precipitation mainly from July to September and tending to be greater in the south and lower in the north. Alpine steppe and alpine meadow are the two main grassland types in the area, and the dominant species are mainly Poa pratensis var. pratensis, Festuca ovina, Stipa aliena, Potentilla chinensis, and Kobresia capillifolia. Seasonal rotational grazing is the most significant and essential way to utilize natural grassland resources. The grassland in the study area can be divided into cold-season grassland (growth from early November to the end of May) and warm-season grassland (growth from early June to the end of October) according to seasonal characteristics. Under the comprehensive influence of human activities (i.e., overgrazing, land reclamation and grassland ecotourism) and climate change (i.e., snow disasters, hailstone and drought), the health and fragmentation of grassland ecosystems in the area have become increasingly serious.

2.2. Grassland Observational Data

Four permanent sampling areas (A1, A2, A3 and A4) are designed in the study area from the northeast to southwest to implement continuous monitoring of grassland vegetation in different seasons depending on the accessibility, grassland utilization pattern and spatial heterogeneity of the selected areas (Figure 1a). These sampling areas are all winter pastures, and there is no grazing during the vegetation growth season. Five fixed sample sites are set up in each sample area. Two field observations were conducted in late July 2017 and mid-November 2017 in the study area, with each observational period lasting approximately 8–10 days. In each field investigation, in addition to observing all the fixed sample sites in the four permanent sampling areas, some random sample sites are also established in other areas based on accessibility and representativeness to obtain more grassland observational data. In total, 66 and 57 sites were observed in July 2017 and November 2017, respectively. Five subplots (0.5 m × 0.5 m) within each site (100 m × 100 m) are established to collect available data (i.e., geographical location, community height, and dominant species). The standing biomass and litter of forage on the soil surface in each subplot are collected using traditional agronomic methods.

2.3. Chemical Analysis

After each round of fieldwork, the grass samples are oven-dried at 65 °C for 48 h to a constant weight, and then ground into a powder for further chemical analysis. Total P content in percentage (g 100 g⁻¹, %) is measured by the phosphomolybdate blue spectrophotometry method. Total N content in percentage (g 100 g⁻¹, %) is assayed by employing an elemental analyzer (Euro EA3000-Single, Euro Vector, Milan, Italy). The forage N:P ratio is then calculated as the ratio between the weight-based contents of N and P.

2.4. Sentinel-2 MSI Data and Processing

As shown in Figure 1b, 12 Sentinel-2 Level-1C MSI images are required to completely cover the entire study area. To make the image data better match the periods of grassland observation (late July 2017 and mid-November 2017), 24 images with less cloud cover within 10 days before and after the observation periods are selected as the remote sensing data in this study. MSI data with the World Geodetic System (WGS84) are obtained from the Copernicus Open Access Hub (https://scihub.copernicus.eu/), and the imaging times are 17 July and 14 November 2017. Moreover, MSI has 13 spectral bands from the visible to NIR to SWIR with a high spatial resolution ranging from 10 to 60 m.

Radiometric calibration, atmospheric correction and resampling of the MSI image are performed using the Sentinel Application Platform (SNAP) desktop software (version: 6.0), which is mainly designed and developed to facilitate the further processing of Sentinel data. After resampling, the image mosaic and mask are implemented in ENVI 5.3 software. Finally, Sentinel-2 mosaic images covering the entire study area with a spatial resolution of 20 m are obtained. Spectral reflectance data from Sentinel-2 image pixels corresponding to five subplots at each sample site are extracted, and then the reflectance is averaged to obtain the representative spectral reflectance of each site.

2.5. Spectral Bands and Vegetation Indices

To validate the performance of Sentinel-2 spectral bands and VIs in estimating the forage N:P ratio in different growing seasons of alpine grassland, a total of 29 spectral variables are used, including 11 raw spectral bands and 18 Sentinel-2-derived VIs (

Table 1. Spectral and spatial resolutions of Sentinel-2 MSI data (B1 and B10 are excluded in this study).

Spectral Bands	Band Center (nm)	Bandwidth (nm)	Spatial Resolution (m)	Spectral Region
B2	490	65	10	Blue
B3	560	35	10	Green
B4	665	30	10	Red
B5	705	15	20	Red edge
B6	740	15	20	Red edge
B7	783	20	20	Red edge
B8	842	115	10	NIR
B8A	865	20	20	NIR
B9	945	20	60	NIR
B11	1375	30	20	SWIR
B12	2190	180	20	SWIR

and

Table 2. Summary of the Sentinel-2-derived spectral vegetation indices (VIs) used in this study.

Index	Name	Formulation	Bands	Reference
NDVI	Normalized difference vegetation index	(R842 − R665)/(R842 + R665)	B8, B4	[38]
NDII	Normalized difference infrared index	(R842 − R1610)/(R842 + R1610)	B8, B11	[39]
NDWI	Normalized difference water index	(R865 − R1610)/(R865 + R1610)	B8A, B11	[40]
NDRE1	Normalized difference red-edge 1	(R740 − R705)/(R740 + R705)	B6, B5	[41]
NDRE2	Normalized difference red-edge 2	(R783 − R705)/(R783 + R705)	B7, B5	[42]
RNDVI	Renormalized normalized difference vegetation index	$\left(R842-R665\right)/\sqrt{R842+R665}$	B8, B4	[43]
GNDVI	Green normalized difference vegetation index	(R865 − R560)/(R865 + R560)	B8A, B3	[44]
EVI	Enhanced vegetation index	2.5 × (R865 − R665)/(1 + R865 + 6 × R665 − 7.5 × R490)	B8A, B4, B2	[45]
REP1	Red-edge position 1	705 + 35 × {[(R842 + R665) × 0.5 − R705]/(R740 − R705)}	B8, B4, B5, B6	[46]
REP2	Red-edge position 2	705 + 35 × {[(R783 + R665) × 0.5 − R705]/(R740 − R705)}	B7, B4, B5, B6	[47]
MTCI	MERRIS terrestrial chlorophyll index	(R842 − R705)/(R705 − R665)	B8, B5, B4	[48]
IRECI	Inverted red-edge chlorophyll index	(R842 − R665)/(R740/R705)	B8, B4, B6, B5	[46]
GCI 1	Green chlorophyll index 1	(R842/R560) − 1	B8, B3	[49]
GCI 2	Green chlorophyll index 2	(R783/R560) − 1	B7, B3	[50]
GCI 3	Green chlorophyll index 3	(R865/R560) − 1	B8A, B3	[50]
RECI 1	Red-edge chlorophyll index 1	(R842/R705) − 1	B8, B5	[49]
RECI 2	Red-edge chlorophyll index 2	(R865/R705) − 1	B8A, B5	[50]
WDRVI	Wide dynamic range vegetation index	(0.1 × R865 − R665)/(0.1 × R865 + R665)	B8A, B4	[51]

). The selection of these VIs is dependent on their performance in previous forage N, P, chlorophyll, and biomass estimation studies [34,36,37].

2.6. Random Forest

In this study, the RF method is employed to estimate the forage N:P ratio from unoptimized and optimized combinations of VIs and spectral bands. RF consists of a mass of decision trees and is widely implemented to handle all kinds of classification and regression problems [52]. The method uses feature randomness and bagging when establishing each individual tree. When bootstrapping, only approximately 67% of the data are used. Approximately 33% of the data (out-of-bag data) are not used in the model and can conveniently be used as a validation set [4]. Moreover, the algorithm is very suitable for analyzing high-dimensional data and robust to nonlinear and unbalanced data [53]. For several key parameters in the RF, the number of predictors (mtry) depends on the square root of the total number of predictors used, while the default values of two for the minimal size of the terminal nodes (nodesize) and 500 for the number of regression trees (ntree) are used. Data operations are performed in MATLAB 2016a, including forage N:P ratio modeling and mapping, variable selection and model accuracy assessment.

2.7. Variable Selection and Cross-Validation

A feature selection method (RFFS) based on the increase in node purity (IncNodePurity) of the RF model is applied to optimize spectral bands and VIs. RFFS ranks the input features based on IncNodePurity, and the sequential backward search method is then used to remove the least important feature with the smallest IncNodePurity value from the feature set every time. The root mean squared error (RMSE) of the validation set of each iteration is calculated to evaluate model prediction accuracy. After many iterations, the feature set with the minimum RMSE and the fewest variables is ultimately selected. To guarantee the stability of the experimental results, 10-fold cross-validation is used [54]. In the validation process, the order of importance of variables generated by each iteration is selected as the basis for eliminating features, and the average RMSE is calculated to evaluate the generalization and prediction ability of the model.

2.8. Accuracy Assessment

The mean absolute error (MAE, Equation (1)), R², RMSE, Akaike information criterion (AIC, Equations (2) and (3)), and Bayesian information criterion (BIC, Equations (2) and (4)) are used to comprehensively evaluate the performance of candidate forage N:P ratio models from multiple perspectives. The R² and RMSE are used to evaluate the accuracy and goodness of fit, respectively, of the simulated and measured forage N:P ratios. MAE is used to assess the modeling errors in estimating the forage N:P ratio [34]. In addition, the AIC and BIC are used to balance the complexity and precision of the estimation model [55,56].

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum }\left|{\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}}}_{\mathrm{i}}\right|$$

(1)

$$\mathrm{RSS}={{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\hat{\mathrm{y}}}_i\right)}^2$$

(2)

$$\mathrm{AIC}=2k+n\ln \left(RSS/n\right)$$

(3)

$$\mathrm{BIC}=n\ln \left(RSS/n\right)+k\ln \left(n\right)$$

(4)

where RSS represents the residual sum of squares between measured and simulated N:P ratios; ${\hat{\mathrm{y}}}_i$ and $y_i$ are the N:P ratios of the validation set for the simulated and measured values, respectively; and n and k are the numbers of observations and variables in the model, respectively.

3. Results

3.1. Variation in the Forage N:P Ratio

The results indicate that the average forage N and P contents in the vigorous growth period (July) are significantly higher than those in the senescing period (November), which indicates that forage N and P gradually decrease in the late growing season, while the average forage N:P ratio shows a pattern contrary to the variation in N and P contents (

Table 3. Descriptive statistics of forage N (%), P (%) and N:P ratio for various periods.

Nutrient	Data Sets	Min	Max	Mean	Median	STDEV	CV (%)	SE	No. of Samples
N	July 2017	1.15	2.78	1.87	1.85	0.30	16	0.04	66
	Nov. 2017	0.32	1.41	0.80	0.80	0.19	23	0.02	57
P	July 2017	0.09	0.29	0.17	0.17	0.04	26	0.01	66
	Nov. 2017	0.03	0.13	0.06	0.06	0.02	32	0.00	57
N:P	July 2017	6.3	21.8	12.0	12.0	3.3	27	0.4	66
	Nov. 2017	5.2	24.8	14.1	13.6	4.4	31	0.6	57

). On the whole, the coefficients of variation in forage N, P and N:P ratio in the senescing period are greater than those in the vigorous growth period, which also means that the grassland vegetation has greater spatial heterogeneity in the senescing period.

The forage N:P ratio is significantly correlated with the N and P contents at different growth stages (

Table 4. Matrix of Pearson correlations between forage N, P and N:P ratio for various periods.

Periods	Nutrient	N	P	N:P
July 2017	N	/	p = 0.1647	p = 0.0005
	P	r = 0.17	/	p = 0.0000
	N:P	r = 0.42	r = −0.74	/
Nov. 2017	N	/	p = 0.0154	p = 0.0000
	P	r = 0.32	/	p = 0.0000
	N:P	r = 0.51	r = −0.59	/

), which further supports using the VIs related to N, P and chlorophyll to estimate the forage N:P ratio. The correlation between the forage N:P ratio and P is strong (r = −0.74 in July 2017 and r = −0.59 in November 2017), indicating that the forage N:P ratio mainly depends upon the forage P content. The correlation between forage N and P is not significant (r = 0.17, p = 0.1647) in the vigorous growth period, but there is a weak correlation between them (r = 0.32) in the senescing period, showing that the correlation between forage N and P changes with the growth stage. Compared with the senescing period, the vigorous growth period appears to be a more suitable period for retrieving the forage N:P ratio.

3.2. Predicting the Forage N:P Ratio with Spectral Bands

The results of the forage N:P ratio estimation based on 11 spectral bands at different growth stages are shown in

Table 5. Performance in estimating the forage N:P ratio with all and optimized Sentinel-2 spectral bands for various periods.

Accuracy Assessment	July 2017		Nov. 2017
	All Bands	Optimized Bands	All Bands	Optimized Bands
No. of variables	11	3	11	5
MAE	2.19	2.1	3.08	2.67
R2	0.41	0.43	0.54	0.55
RMSE	2.5433	2.4783	3.5299	3.1352
AIC	138.36	121.29	151.19	127.09
BIC	161.39	127.58	172.23	136.65

. Using the RFFS algorithm to optimize the spectral bands slightly increases the performance of modeling the forage N:P ratio and reduces the complexity of the model (i.e., the model has a lower AIC and BIC) compared with those observed for the forage N:P ratio model established with entire bands. As shown in Figure 2 and Figure 3, three bands (B8-NIR, B9-NIR and B12-SWIR) and five bands (B3-Green, B4-Red, B5-Red edge, B9-NIR, and B12-SWIR) are selected as important bands for estimating the forage N:P ratio in the vigorous growth period (July) and the senescing period (November), with R² values of 0.43 and 0.55, respectively.

Figure 3a shows the analysis results for the relative importance of different bands in the estimation of the forage N:P ratio. Overall, B9 and B12 play an important part in the estimation of the forage N:P ratio during the vigorous growth period and the senescing period, especially B9, which greatly contributes to forage N:P ratio estimation. B11 and B12 are more important for estimating the forage N:P ratio in the vigorous growth period, but they are less sensitive to the forage N:P ratio detection in the senescing period. In addition, B3 and B5 have a good potential in retrieving the forage N:P ratio in the vigorous growth period, but they have no significant effect on the estimation of the N:P ratio in the senescing period. It can be concluded that the bands sensitive to the estimation of the forage N:P ratio in the vigorous growth period are mainly distributed in the NIR and SWIR regions. Moreover, some bands located in the red (B4) and red-edge (B5) regions also significantly influence the estimation of the forage N:P ratio in the senescing period, in addition to the NIR and SWIR regions.

3.3. Predicting the Forage N:P Ratio with Sentinel-2 Vegetation Indices

The results of the forage N:P ratio estimation based on 18 VIs at different growth stages are shown in

Table 6. Performance in estimating the forage N:P ratio with all and optimized Sentinel-2 vegetation indices (VIs) for various periods.

Accuracy Assessment	July 2017		Nov. 2017
	All VIs	Optimized VIs	All VIs	Optimized VIs
No. of variables	18	4	18	5
MAE	2.06	2.02	3.46	2.94
R2	0.42	0.43	0.41	0.42
RMSE	2.4592	2.4756	4.0742	3.4881
AIC	146.93	118.68	180.48	137.60
BIC	184.63	127.06	214.9	147.16

. Using the RFFS algorithm to optimize the VIs slightly increases the performance of modeling the forage N:P ratio and reduces the complexity of the model (i.e., the model has a lower AIC and BIC) compared with those observed for the forage N:P ratio model established with entire VIs. As shown in Figure 2 and Figure 3, four VIs (enhanced vegetation index (EVI), green chlorophyll index (GCI1), normalized difference infrared index (NDII) and red-edge chlorophyll index (RECI2)) and five VIs (normalized difference red edge 1 (NDRE1), NDRE2, RECI2, normalized difference water index (NDWI) and wide dynamic range vegetation index (WDRVI)) are selected as important VIs for estimating the forage N:P ratio in the vigorous growth period (July) and the senescing period (November), with R² values of 0.43 and 0.42, respectively.

Figure 3b shows the analysis results for the relative importance of different VIs for estimating forage N:P ratio. Overall, RECI2 is important for the forage N:P ratio estimation in both the vigorous growth period and the senescing period. NDII and RECI2 significantly contribute to the forage N:P ratio estimation during the vigorous growth period, and NDWI, RECI2, NDRE1 and NDRE2 are the four most important VIs for the forage N:P ratio estimation in the senescing period. In addition, the seven VIs sensitive to the forage N:P ratio in the vigorous growth period all include the NIR band (B8 or B8A). Among the 5 VIs sensitive to the forage N:P ratio in the senescing period, B5 and B8A in the red-edge and NIR regions are the most frequent bands, respectively. These results indicate that the NIR band has a good potential in estimating the forage N:P ratio as a whole.

3.4. Predicting the Forage N:P Ratio with a Combination of Spectral Bands and Vegetation Indices

The optimized spectral bands and VIs are integrated to establish an RF model for the forage N:P ratio estimation at different growth stages (

Table 7. Performance in estimating the forage N:P ratio with a combination of optimal spectral bands and vegetation indices (VIs).

Accuracy Assessment	Optimized Bands + VIs
	July 2017	Nov. 2017
No. of variables	7	10
MAE	1.93	2.72
R2	0.49	0.59
RMSE	2.2661	3.1095
AIC	117.50	136.58
BIC	132.16	155.70

). Compared with the forage N:P ratio model established by using the optimized spectral bands (n = 3) or VIs (n = 4) alone, the prediction accuracy of the integrated model (R² = 0.49) is improved by 0.06 in the vigorous growth period (July). The accuracy of the integrated N:P ratio model (R² = 0.59) is improved by 0.04 in the senescing period (November) compared with the model that used the optimized spectral bands (n = 5) or VIs (n = 5) alone. This indicates that the integration of optimized spectral bands and VIs increases the estimation accuracy of the forage N:P ratio.

3.5. Mapping of Potential Forage N and P Limitation

According to the integrated model of the forage N:P ratio in the two periods mentioned above, the spatial distribution of the forage N:P ratio is mapped based on Sentinel-2 MSI images. The spatial distribution and nutrient limitation of the forage N:P ratio at different growth stages are shown in Figure 4 (the vigorous growth period) and Figure 5 (the senescing period). The results indicate that the forage N:P ratio ranges from 8.06 to 18.95 in the vigorous growth period. Overall, the growth of forage in most of the study area is limited by N (N:P < 14), the nutrient condition of forage in the local area (i.e., northern Xiahe County, southern Luqu County and northwestern Maqu County) is affected by P deficiency (N:P > 16), and the forage growth in very few areas is restricted by both N and P (14 < N:P < 16) (Figure 4). Moreover, the forage N:P ratio is between 7.75 and 22 in the senescing period (Figure 5). Comparative analysis of the spatial distribution of the forage N:P ratio in the above two periods shows that the forage N:P ratio increases in most areas and decreases slightly in local areas, which suggests that the P content in forage is more deficient than the N content when grassland vegetation stops growing.

4. Discussion

4.1. Potential of Sentinel-2 Spectral Bands and Vegetation Indices in Estimating the Forage N:P Ratio

This study shows that the Sentinel-2 spectral bands are mainly distributed in the NIR and SWIR regions have an important contribution to the forage N:P estimation of alpine grassland, and some spectral bands located in the red and red-edge regions have a great influence on estimating the forage N:P ratio in the senescing period. In addition, some VIs (such as RECI2) constructed by the Sentinel-2 spectral bands in the red and NIR regions have a good potential in the inversion of the forage N:P ratio. Previous studies have shown that there is a significant correlation between the spectral reflectance of the red-edge and NIR regions and the N content in plant leaves [57,58]. The bands that are sensitive to the P content in vegetation are mainly distributed in the NIR and SWIR regions [59], and VIs based on the red-edge are important variables for predicting the N content in leaves [4]. At different growth stages, the N:P ratio is significantly correlated with the forage N and P contents (Table 4). Therefore, the bands sensitive to N and P may also be applicable to the estimation of N:P ratio in forage. Romoelo et al. [11] estimated the N:P of savanna forages based on hyperspectral data and found that the wavebands were distributed in the SWIR region had a higher sensitivity to N:P in leaves, which is similar to our finding. In different growth stages of grassland, the optimization results of spectral bands based on the RF algorithm show that both B9 and B12 are important bands for estimating the forage N:P ratio. The potential of B12 in estimating the forage N:P ratio can be understood by the absorption feature of N near 2180 nm.

The detection of the forage N using remote sensing technology is generally associated with the absorption of chlorophyll, and usually includes the spectral regions related to chlorophyll detection (i.e., the red-edge and NIR regions) [33]. Previous studies have shown that some known specific absorption bands for chlorophyll, N, and proteins (i.e., 640 nm, 660 nm, 910 nm, and 1510 nm) can be successfully employed to retrieve the forage N content [11,60,61]. The specific absorption features of P are mainly distributed in the NIR and SWIR regions, having broad application prospect in estimating forage P [6,62]. The detection of forage P applying hyperspectral measurements could be related to the stoichiometry (i.e., indirect estimation according to the relationship between other biochemical parameters and P in forage) [63,64]. Preliminary results show the forage N:P ratio is significant correlated with the N and P contents (Table 4), which further supports the feasibility of using the spectral variables related to N, P and chlorophyll to indirectly estimate the forage N:P ratio.

This study confirms that it is feasible to use Sentinel-2 MSI spectral bands and VIs to estimate the forage N:P ratio, and combining optimized spectral bands and VIs further increases the estimation accuracy of the forage N:P ratio. Previous studies have also shown that combining spectral bands and VIs improves the estimation accuracy of vegetation biochemical parameters [4,34], which benefits from enhancing some of the most important variables and integrating effective information to obtain the best results.

4.2. Effects of Different Seasons on Forage N:P Inversion in Natural Alpine Grasslands

The growth and developmental stages of forage directly affect its nutrient content, thereby affecting the feed intake and production performance of animals. This study shows that the average forage N and P contents in the vigorous growth period (July) are significantly higher than those in the senescing period (November) at different growth stages. Many studies have also found that the contents of N, P and protein in forage gradually decrease as plant grows [65]. For example, the P content of eight forages in the Alxa Desert steppe show regular seasonal dynamics; that is, the P content gradually decreases with the advancement of the growth period [66], and the N content in the aboveground part of plants decreases during the plant growth season (May to September) [67]. The forage N content gradually decreases as the forage grows in alpine grassland, especially the forage N that decreases rapidly during the senescing period [68]. In the late growth season, with the expansion of plants, the forage gradually withers, cell senescence occurs, and fiber material increases, resulting in the dilution effect of elements [69]. In addition, plants gradually stop growing after entering the senescing period, and some of the aboveground nutrients are transferred underground for storage to supply the consumption of plants in non-growing period and nutrients needed for plant germination in the coming year [70]. This study reveals that the growth of forage is mainly limited by N during the vigorous growth period in alpine grasslands. With the gradual senescing of forage, the N:P ratio increases significantly, and the vegetation gradually changes from being N limited to P deficient. Especially at the end of the growing season, the limitation effect of P on vegetation growth is greater than that of N, which is potentially related to P being fixed in plant roots during the senescence of forage.

This study shows that the integrated forage N:P ratio model achieves ideal estimation accuracy in both the vigorous growth period (R² = 0.49) and the senescing period (R² = 0.59) (Table 7). This finding indicates that some spectral bands and VIs from Sentinel-2 can capture nutrient limitation information for forage growth to a certain extent. Using multispectral data to estimate the forage canopy N:P ratio provides an indirect measurement of the deficiency status of N and P in grassland vegetation, and some spectral bands and VIs that are sensitive to forage N, P and chlorophyll play an important role in N:P ratio estimation [11,26]. The biomass of alpine grasslands reaches its maximum and grazing intensity is relatively high during the vigorous growth period. In terms of practical applications, it is of more practical significance to use remote sensing data to monitor forage growth in this stage than in other stages at a regional scale. In general, the nutrient content of forage decreases and the quality deteriorates during the senescing period, and the canopy spectrum is easily affected by soil and other background factors, which will reduce the sensitivity of the spectral bands and VIs to forage N and P [68,71]. According to our research results, the integrated forage N:P ratio model performs better in the senescing period than in the vigorous growth period, which is potentially attributable to the limited number of sample sites, the optimization of model parameters, and large spatial heterogeneity of grassland. Subsequent research will also explore this problem more deeply based on multisource satellite data and multiple machine learning algorithms.

4.3. Future Perspectives

Although this study maps the spatial distribution of the forage N:P ratio in the vigorous growth period and the senescing period on a regional scale, it is also necessary to understand the potential limitations of this method. The empirical estimation model of vegetation biochemical parameters is susceptible to the spatial heterogeneity of the study area, which reduces the generalization performance of the model in other regions, and therefore still faces many challenges in practical application [19,72]. Due to the cloudy and rainy weather in summer, the grass canopy spectrum and the cloud cover of the optical satellite images are seriously limited by the weather conditions, which confines the remote sensing monitoring of vegetation growth in the study area to a certain extent.

For vegetation applications, hyperspectral sensors (i.e., HYMAP, Hyperion, AISA Eagle and Gaofen-5) can provide abundant spectral information and detect minute spectral features that are masked by the broad bands of satellites [73]. In addition, spectral bands are favorable to the calculation of the red-edge position and narrow band indices, important for vegetation biochemical parameters assessments [74]. Moreover, hyperspectral data have received extensive attention for species discrimination and mapping subtle variations in vegetation biochemical information. However, the high cost and the small swath width involved limits its usage to small spatial extents [14,27]. The recent trend is that hyperspectral data is gradually shifting towards to affordable multispectral sensors with strategic bands (i.e., RapidEye, WorldView-2, Sentinel-2A/B and Gaofen-6), which capable of retrieving plant macronutrients [75,76,77]. Near-surface unmanned aerial vehicle (UAV) remote sensing systems have the advantages of a small size, wide applications and high flexibility, which can overcome the limitations of conventional ground observation-based methods [78,79]. Therefore, exploring integrated monitoring technology for grass biochemical parameters based on near-surface UAV images and multisource, multitemporal and high-resolution satellite data will be the focus of future research.

As RTMs are capable of explaining the interaction and transfer of radiation inside the canopy according to the laws of physics, they offer a specific connection between the canopy reflectance and the vegetation biochemical parameters [80]. A physically based models (i.e., PROSPECT and LIBERTY) have the advantages that a large simulated spectra database can be generated by changing the input parameters [81], which may be applied to other areas. Many studies show that the model inversion has certain potential for estimating vegetation biophysical parameters in heterogeneous grasslands applying hyperspectral data [17,82]. The combination of empirical and physically based models is a reliable approach for estimating forage N and P, although its generalization performance needs to be verified by other data sets.

5. Conclusions

Based on high-resolution multispectral satellite data and an RF machine learning algorithm, this study evaluates the feasibility of using Sentinel-2 MSI spectral bands and VIs to retrieve the forage N:P ratio of alpine grassland on a regional scale. The results indicate that combining optimized spectral bands and VIs increases the accuracy of the forage N:P ratio estimation in comparison with that of models using spectral bands or VIs alone. In addition, the Sentinel-2 spectral bands are mainly distributed in the NIR and SWIR regions have a valuable contribution toward estimating the forage N:P ratio, and some spectral bands located in the red and red-edge regions greatly influence the forage N:P ratio in the senescing period. This study demonstrates the practicality and feasibility of using Sentinel-2 spectral bands and VIs to directly estimate the forage N:P ratio at regional scales. Moreover, airborne and satellite data can be integrated to diagnose the seasonal variations of nutrition limitation at a landscape level in the alpine grassland ecosystem of the Tibetan Plateau and thus represents a promising future development trend.