Estimating Carbon, Nitrogen, and Phosphorus Contents of West–East Grassland Transect in Inner Mongolia Based on Sentinel-2 and Meteorological Data

Abstract

Estimating the carbon (C), nitrogen (N), and phosphorus (P) contents of a large-span grassland transect is essential for evaluating ecosystem functioning and monitoring biogeochemical cycles. However, the field measurements are scattered, such that they cannot indicate the continuous gradient change in the grassland transect. Although remote sensing methods have been applied for the estimation of nutrient elements at the local scale in recent years, few studies have considered the effective estimation of C, N, and P contents over large-span grassland transects with complex environment including a variety of grassland types (i.e., meadow, typical grassland, and desert grassland). In this paper, an information enhancement algorithm (involving spectral enhancement, regional enhancement, and feature enhancement) is used to extract the weak information related to C, N, and P. First, the spectral simulation algorithm is used to enhance the spectral information of Sentinel-2 imagery. Then, the enhanced spectra and meteorological data are fused to express regional characteristics and the fractional differential (FD) algorithm is used to extract sensitive spectral features related to C, N, and P, in order to construct a partial least-squares regression (PLSR) model. Finally, the C, N, and P contents are estimated over a West–East grassland transect in Inner Mongolia, China. The results demonstrate that: (i) the contents of C, N, and P in large-span transects can be effectively estimated through use of the information enhancement method involving spectral enhancement, regional feature enhancement, and information enhancement, for which the estimation accuracies (R²) were 0.88, 0.78, and 0.85, respectively. Compared with the estimation results of raw Sentinel-2 imagery, the RMSE was reduced by 3.42 g/m², 0.14 g/m², and 13.73 mg/m², respectively; and (ii) the continuous change trend and spatial distribution characteristics of C, N, and P contents in the west–east transect of the Inner Mongolia Plateau were obtained, which showed decreasing trends in C, N, and P contents from east to west and the characteristics of meadow > typical grassland > desert grassland. Thus, the information enhancement algorithm can help to improve estimates of C, N, and P contents when considering large-span grassland transects.

1. Introduction

Carbon (C), nitrogen (N), and phosphorus (P), as the essential three elements for various ecosystems, are closely related to plant health status, ecosystem structure, and function [1,2]. Quantitative analysis of C, N, and P contents plays an important role in exploring the interactions between soil and vegetation, as well as soil organic matter storage [3]. For example, C is the main component of carbohydrates and cellulose, accounting for about 50% of the dry matter content [4], and provides energy for basic ecosystem biogeochemical process [2]. Meanwhile, N and P are the basic nutrients for plant growth, and play a key role in leaf growth, photosynthesis, energy storage, and transmission [5,6]. The C:N ratio can reflect the growth rate and nutrient utilization efficiency of vegetation [7], while the N:P ratio can characterize the uptake of N and P nutrients by plants and reflects the nutrient limitation information of primary productivity [8]. As plant growth and the contents of C, N, and P elements are driven by the climate and soil fertility [9], N and P also have important feedback effects on climate change. Therefore, it is of great significance to obtain the spatial distributions of C, N, and P in grassland quickly, accurately, and objectively, in order to monitor grassland health, estimate grassland primary productivity, analyze grassland–soil interaction models, and monitor the energy cycle.

Transects often extend more than 1000 km along specific environmental gradients (considering factors such as temperature or precipitation) or land-use intensity gradients, providing an important survey method in the field of ecology for studying and revealing the principles and processes of regional or global climate change [10]. However, the traditional methods for obtaining C, N, and P contents through field sampling and manual measurement require a lot of manpower and time, and it is difficult to obtain the gradient changes of nutrients in the whole transect continuously, making large-span grassland transect monitoring unsuitable. Remote sensing techniques provide non-destructive and efficient methods to estimate plant nutrient content [11], and have become a new means for acquiring the physical and chemical parameters of plants. In particular, the development of imaging spectroscopy technology integrates spectroscopy and imagery, which can capture the unique spectral characteristics of different elements and their compounds on the Earth, as well as detecting the specific components of various substances, such as the nutrient contents of N, P, and K in grassland, wheat, and nut crops [12,13,14]. Whether the relevant sensitive features can be extracted is closely related to the accuracy of the model.

Although remote sensing can help to obtain nutrient element data over large-span grassland transects, the weak sensitivity of spectra to C, N, and P is still an important factor restricting the associated estimation accuracy. The leaf spectra contain a lot of information: visible light (400–700 nm) is related to the pigment concentration, the near-infrared region (700–1300 nm; NIR) is related to the leaf structure, and the short-wave infrared region (1300–2500 nm; SWIR) is related to the spectral absorption characteristics of water and proteins [15]. SWIR is also a sensitive region for N and P, and their absorption characteristics are easily covered by atmospheric water vapor or leaf water [16]. The detection of C mainly focuses on the absorption characteristics of lignin, cellulose, and carbohydrates, and there is no clear band for the C element. Thus, due to a lack of observable and clear absorption characteristics, the estimation of C, N, and P contents is still a challenge, in terms of plant traits [17]. Although the construction of vegetation indices (VIs) based on spectral characteristics to estimate chlorophyll, biomass, LWC, and so on is a commonly used method, most VIs only take advantage of a small amount of information and cannot make full use of all spectral information [17]. Furthermore, the VIs that can be directly used to estimate C, N, and P are also rare, such that identifying sensitive wavelengths over the entire spectral range, rather than the limited band, provides a more accurate method for nutrient estimation. The reflectance of canopy spectra is affected by comprehensive factors, such as soil background, illumination geometry, canopy structure, leaf water content, and growth stage [18]. Spectral transformation methods can reduce background noise [19], reduce hyperspectral data redundancy, and enhance weakly sensitive information related to C, N, and P elements, such as fractional differential (FD), continuous wavelet transform (CWT), logarithmic transformation (Log(1/R)), and continuous absorption feature removal (CR) [20,21,22,23]. Compared with Log(1/R), CWT, and CR algorithms, FD has the advantage of not requiring parameters to be set artificially and multi-scale transformation of spectra, which has been widely used in agriculture, ecology, soil, and other fields [24,25]; however, at present, multi-scale spectral transformation methods are most suitable for hyperspectral data with rich spectral information. How to take advantage of multi-spectral satellite data for the estimation of large-span grass transects traits still needs to be explored.

Limited spectral information is an important restricting factor for the estimation accuracy of C, N, and P using multi-spectral imagery, and many studies have shown that hyperspectral images are generally more advantageous than multi-spectral images for the investigation of plant traits [26], and hyperspectral images with more spectral information have been used to detect changes physical and chemical parameters of the plant [27]. Especially in satellite remote sensing, they can provide large-scale, long-term sequential, and consistent observation data. Compared with multispectral satellites, hyperspectral satellites have relatively few missions (e.g., the EO-1 Hyperion imaging spectrometer, Compact High-Resolution Imaging Spectrometer, CHRIS, on board of PROBA−1, the Chinese HJ-1A HSI sensor, and so on [28,29]), narrow swaths, and high prices, and airborne hyperspectral capture is limited by the flight distance, making it only suitable for small areas. Thus, it is difficult to obtain the hyperspectral data required for the acquisition of transect survey data, and it is imperative to build a high-precision large-scale grassland transect estimation model based on multi-spectral data, considering their economic, spatial, and temporal advantages [30]. The fusion data of hyperspectral and multispectral images can be used to obtain high-spatial and -spectral resolution images, which is the main method used to improve the image information [31]; however, it is generally necessary to obtain the hyperspectral and multispectral images of the same region at the same time. Some scholars have used the physical model PROSAIL to simulate hyperspectral information [32], but the input parameters of the physical model (e.g., chlorophyll a + b content, carotenoids, leaf area index, leaf angle) are difficult to obtain. Pang et al. [30] have proposed a satellite-scale spectral simulation algorithm based on the decomposition of mixed pixels, which can increase the spectral information of multi-spectral satellites, to a certain extent, and obtained high estimation accuracy in re-biomass estimation by combining it with the multi-scale spectral transformation method, which has reference significance for the large-scale estimation of C, N, and P elements.

In previous studies, the estimation of C, N, and P has mostly focused on hyperspectral data measured by spectrometers or UAVs, the range of the study area was relatively small (e.g., sample scale [17] or local scale [1]), and only one grassland type was considered. Thus, the climate and grassland types can be ignored for C, N, and P estimation within a single time scale. However, the transect of the Inner Mongolia Plateau (109° E–125° E, 43.5° N–45.3° N, Figure 1), about 920 km from east to west, including meadow, typical grassland, and desert grassland, is an area of sensitive grassland ecological change in northern China. There are features related to the C, N, and P content variables that are subtle and, in such a complex environment, how to extract effective spectral features based on multispectral images and construct an estimation model suitable for the whole transect still requires further study. To solve this problem, this paper is mainly divided into three parts: (i) The sensitive features related to C, N, and P are extracted by spectral enhancement, regional characteristic enhancement, and information enhancement methods; (ii) an elemental content estimation model for the Inner Mongolia transect is constructed using a partial least-squares regression model; and (iii) we evaluate the estimation accuracy, draw the spatial distribution map, and analyze the changing trend of C, N, and P over the whole transect.

2. Data and Methods

2.1. Study Area

The study area (Figure 1; 109° E–125° E, 43.5° N–45.3° N) is located in the central part of the Inner Mongolia Autonomous Region and western Jilin Province, and crosses Songyuan City of Liaoning Province, Tongliao City of Inner Mongolia Autonomous Region, Xilinguole League, and other cities from east to west. The transects cover a temperate continental climate with an annual average temperature of −1 °C to 10 °C and annual precipitation of 50 mm to 450 mm [33]. Previous studies have shown that water is the main factor restricting the growth of vegetation and vegetation type changes in the Mongolia Plateau [34]. Therefore, a precipitation gradient belt from east to west is set up, covering meadows, typical grasslands, and desert grasslands [34]. The transect is an economical and scientific survey method following the rules of the International Geosphere-Biosphere Program. In recent years, the grassland ecosystem in the Mongolian Plateau has been degraded, due to global climate change and human activities, making it a sensitive area which is vulnerable to ecological damage [35]. Timely monitoring of forage nutrient content in the transect is of great significance for monitoring global climate change and grassland ecosystem protection.

2.2. Data

2.2.1. Field Data

Field sampling was conducted from 3–14 July 2018, at the growth peak of the grasslands, and data of aboveground biomass and geographic information (latitude, longitude, and altitude) were obtained. According to the grassland type, ten sampling areas were arranged at intervals of 90–200 km from east to west, and eight 1 × 1 m quadrats at intervals of 15–55 m were arranged in each sample area, with a total of 80 sampling points (Figure 1). Among them, Z1–Z3 were classified as meadow, Z4–Z7 as grassland, and Z8–Z10 as desert grassland.

Under clear-sky conditions between 10:00 and 14:00, at a height of 1.3 m above the ground, the canopy hyperspectral reflectance (325–1075 nm) of samples was measured using a field spectrometer (ASD Field Spec 2 spectrometer, Analytical Spectral Devices, Boulder, CO, USA). The pure spectra of bare soil and grass which had a coverage of more than 95% were measured, and all measurements were repeated three times at each point, with the average value taken as the final reflectance. The pure end-member spectra were used to construct the enhanced spectra, and the measured canopy spectra of the sample were used to verify the accuracy of the spectral enhancement algorithm and the estimation accuracy for C, N, and P contents.

All of the grasslands in the quadrats were mowed neatly, marked, and refrigerated. After the end of spectral acquisition, all samples were mowed neatly, marked, and refrigerated until further analysis. The plant samples were washed and dried in a 60 °C oven until a constant weight was reached, used for water content analysis, and then powdered using a grinding machine (MM400 ball mill, Retsch, Haan, Germany) for C, N, and P measurements. Parts of samples were acidified with trace metal-grade 68% HNO₃ overnight, then samples were digested using a microwave digestion system (Mars X press Microwave Digestion system, CEM, Matthews, NC, USA) before elemental analysis. Finally, total contents of C and N were measured using an elemental analyzer (VarioMAX CN Elemental Analyzer, Elementar, Hanau, Germany), and total P content was analyzed using an inductively coupled plasma optical emission spectrometer (ICP-OES, Optima 5300 DV, Perkin Elmer, Waltham, MA, USA) [34,36]. The sampling statistics for Z1–Z10 are provided in

Table 1. The statistics of mean, standard deviation (Std), maximum (Max), and minimum (Min) for the considered sample sites.

	Sites	Mean	Std	Min	Max	Sites	Mean	Std	Min	Max
C	Z1	51.67	12.82	40.23	67.60	Z6	31.58	11.07	20.78	55.73
	Z2	41.75	11.89	22.82	61.30	Z7	18.13	4.24	14.76	24.95
	Z3	65.54	15.06	44.42	91.43	Z8	14.52	3.57	10.78	20.34
	Z4	35.88	14.42	20.56	60.01	Z9	5.65	1.59	3.96	8.95
	Z5	14.27	4.84	5.40	19.62	Z10	4.99	1.17	3.32	7.02
N	Z1	2.83	1.01	1.99	5.03	Z6	1.68	0.63	1.14	3.08
	Z2	2.43	0.69	1.31	3.54	Z7	1.18	0.30	0.92	1.69
	Z3	3.54	0.78	2.41	4.90	Z8	1.21	0.24	0.99	1.60
	Z4	1.92	0.74	1.14	3.05	Z9	0.57	0.12	0.42	0.75
	Z5	0.95	0.30	0.34	1.29	Z10	0.54	0.13	0.34	0.75
P	Z1	230.51	59.45	167.21	323.05	Z6	74.10	27.60	46.35	133.47
	Z2	139.30	42.52	77.25	215.57	Z7	64.61	16.04	51.99	95.50
	Z3	203.97	44.84	142.69	279.45	Z8	52.95	13.95	37.83	74.28
	Z4	112.17	47.93	64.22	181.72	Z9	44.99	9.33	33.33	61.16
	Z5	45.10	14.35	16.09	62.28	Z10	33.21	8.65	19.78	47.57

Note: The units of C, N, P are g/m², g/m², mg/m², respectively.

.

The climatic environment is an important factor affecting grassland types in the Inner Mongolia transect, and the meteorological data needed in this paper were obtained from the China Meteorological Data Center (http://data.cma.cn, accessed 29 August 2020), which mainly include temperature, precipitation, air pressure, wind speed, sunshine duration, and other data from 80 meteorological stations around the transect, from April to July 2018 (e.g., re-greening to cutting). The positions of these meteorological stations are shown in Figure 1, and all data were used to calculate the relative humidity characteristics of the transect, in order to express the differences between regions.

2.2.2. Remote Sensing Data

When carrying out large-span transect analysis, due to a lack of hyperspectral data sources and the large study area, it is often difficult to obtain suitable hyperspectral images. Thus, multispectral images have become the primary choice for large-span transect monitoring. Sentinel-2 imagery has become one of the main data sources in agricultural and ecological remote sensing, as well as other fields, due to its high temporal resolution (5 d) and high spatial resolution (B2, B3, B4, B8 = 10 m; B5, B6, B7, B8a, B11, B12 = 20 m; B1, B9, B10 = 60 m). All Sentinel-2 bands except for B10 were used. To better match with ground sampling data, we selected 13 scenes of Sentinel-2 A\B images with no or few clouds as data sources, between 15 days before and after the sampling time points (1–30 July 2018), and used the SNAP image processing software provided by the European Space Agency ESA (website: https://scihub.copernicus.eu/dhus/#/home, accessed 27 March 2021) to pre-process the data, including radiometric calibration, atmospheric correction, geometric correction, resampling, cloud masking (threshold B2 > 0.2), and so on. The Sentinel-2 image information and corresponding sampling times are given in

Table 2. The information of Sentinel-2 images and the corresponding sampling time in the study area.

Sampling Regions	Sampling Times	Image Scanning Time	Image Names (Sentinel-2)
Z1	10 July 2018	1 July 2018	S2B_MSIL2A_20180701T024549_N0206_R132_T51TWK
Z2	9 July 2018	10 July 2018	S2A_MSIL2A_20180719T025551_N0206_R032_T51TUK
Z3	11 July 2018	19 July 2018	S2A_MSIL2A_20180709T025551_N0206_R032_T50TQR
Z4	8 July 2018	9 July 2018	S2A_MSIL2A_20180709T025551_N0206_R032_T50TPQ
Z5	6 July 2018	12 July 2018	S2A_MSIL2A_20180712T030541_N0206_R075_T50TMQ
Z6	5 July 2018	12 July 2018	S2A_MSIL2A_20180712T030541_N0206_R075_T50TMQ
Z7	7 July 2018	12 July 2018	S2A_MSIL2A_20180712T030541_N0206_R075_T50TNQ
Z8	12 July 2018	10 July 2018	S2B_MSIL2A_20180710T031539_N0206_R118_T50TLP
Z9	14 July 2018	13 July 2018	S2B_MSIL2A_20180713T032539_N0206_R018_T49TGJ
Z10	13 July 2018	13 July 2018	S2B_MSIL2A_20180713T032539_N0206_R018_T49TEJ

.

2.3. Methods

2.3.1. Method Overview

As shown in the flow chart (Figure 2), the work in this paper mainly includes three parts: data, methods, and results and evaluation. The data collected in this paper mainly include the measured data (C, N, and P contents; end-member spectra; measured canopy spectra at grassland sample scale), Sentinel-2 images, meteorological data (temperature, precipitation, wind speed, light, and so on), and the abundance information of grassland and bare soil in the transect, as calculated by Sentinel-2. For the estimation of C, N, and P in large-span grassland transects, regional characteristic differences and weak spectral sensitivity have always been important factors restricting the estimation accuracy. Thus, we used spectral enhancement, regional enhancement, and information enhancement methods to extract the sensitive features to construct the C, N, and P estimation model. First, the abundance information for grassland and bare soil was calculated from Sentinel-2 images, and the fitting spectra were constructed using the abundance information and the measured end-member hyperspectral data of the main ground objects. Then, the simulated spectra algorithm was used to enhance the spectral information of the raw Sentinel-2 imagery, and the accuracies of enhanced spectra were verified using the measured canopy spectra. Second, the relative humidity information, which is closely related to the grassland type change in the Inner Mongolia transect, was calculated from the meteorological data, and the relative humidity and enhanced spectra were used to construct the regional enhanced fused spectrum (FS). Finally, the weak information related to C, N, and P was enhanced using the fractional differential spectral transformation method, and the sequence forward feature selection algorithm was used to select the sensitive feature bands to construct the C, N, and P estimation model. To verify the estimation accuracy of the method proposed in this paper, the estimation results from the raw spectrum (RS) and the FS-FD model were compared with the measured data, the variation trends of element content over the whole transect were analyzed, and the spatial distributions of C, N, and P were mapped.

2.3.2. Sentinel-2 Spectral Information Enhancement

As shown in Figure 3a, the spatial resolution of each pixel in Sentinel-2 is 10 m, which contains the characteristics of vegetation, bare soil, and other ground objects. The spectral reflectance of the pixel is the comprehensive embodiment of the regional characteristics at the unit scale, which is considered to be the real spectra. However, the spectral resolution of the wide-band Sentinel-2 is weakly sensitive to C, N, and P elements. Pang et al. [30] have proposed that, without changing the real spectral reflectance trend, the spectral details of the main ground objects can be assigned to the raw image using the abundance information of ground object types obtained by mixed pixel decomposition as a weight, which can improve the monitoring ability of Sentinel-2 for grassland and achieved significant results in biomass estimation. In this paper, a spectral enhancement algorithm is used to improve the spectral details of raw Sentinel-2 imagery. The calculation process of the spectral enhancement algorithm is shown in Equations (1)–(7) and schematic diagram in Figure 3b.

$$E{S}_{\lambda }=Fi{t}_{\lambda }+\mathrm{VL}-{\mathrm{TL}}_{\lambda },$$

(1)

where $E{S}_{\lambda }$ is the enhanced spectra at wavelength λ, $Fi{t}_{\lambda }$ is the spectra fitted by vegetation coverage and pure end-members, VL is the broken line of Sentinel-2 raw spectra, and TL is the trend line of $Fi{t}_{\lambda }$. The formula for TL is as follows:

$$T{L}_{\lambda }=\left\{\begin{array}{c}\hspace{1em}{a}_1\lambda +b_1, \hspace{1em}\hspace{1em}{W}_1\le \lambda <W_2\hfill \\ \hspace{1em}{a}_2\lambda +b_2,\hspace{1em}\hspace{1em}{W}_2\le \lambda <W_3\hfill \\ \cdots \cdots \\ a_{k-1}\lambda +b_{k-1}, W_{k-1}\le \lambda \le W_k\hfill \end{array}\right.,$$

(2)

where k is the band number of the multispectral satellite, $W_k$ is the central wavelength of the kth band of the multispectral satellite, and $a_k$ and $b_k$ are the slope and intercept of the trend line in [$W_{k-1}, W_k$], respectively. The associated calculation formula is as follows:

$$\left\{\begin{array}{c}{a}_{k-1}=\frac{Fi{t}_{W_k}-Fi{t}_{W_{k-1}}}{W_k-W_{k-1}}\\ b_{k-1}=Fi{t}_{W_{k-1}}-a_{k-1}{W}_{k-1}\end{array}\right.,$$

(3)

where $Fi{t}_{W_k}$ is a spectral curve, which is fitted according to the type of land. Its calculation formula is as follows:

$$Fi{t}_{\lambda }={\displaystyle \sum }_{j=1}^n{f}_j\times p_{\lambda ,j},$$

(4)

where $f_j,\mathrm{and}{p}_{\lambda ,j}$ are the weight and reflectance of the end-member spectra of class j, respectively.

2.3.3. Regional Characteristic Enhancement

There are many grassland types (meadow, typical grassland, and desert grassland) in the Inner Mongolia transect, thus creating a complex ecological environment, and the aboveground biomass (AGB) is significantly different in different grassland types, which is affected by the climate. The contents of C, N, and P, at unit scale, are closely related to the AGB, such that the differences in characteristics by grassland type between regions must be considered in large-span grassland transect monitoring. Although spectral features can reflect vegetation growth information, to a certain extent, there are two drawbacks in estimating C, N, and P contents when transect monitoring only by spectral features: (i) Reflectance is the comprehensive response of multiple factors, such as vegetation type and canopy structure, in the unit pixel, which are prone to the phenomenon of “same subject with heterospectrum” and “foreign body with spectrum”. Mixed pixels are an important factor affecting the accuracy of ground monitoring. Muhuri [37] and Touzi et al. [38] have proposed decomposition algorithms for image mixed information in optical images and Sar images, respectively, and improved the accuracy of the original image to obtain ground object information. Although the enhanced spectra also decomposed mixed pixels, to a certain extent, the introduction of a meteorological weight factor for regional differences can help to further reduce the influence of the heterospectrum phenomenon; and (ii) the differences between grassland types is ignored. Assuming that the elemental content is estimated by partial least-squares regression (Equation (5)) without considering the grassland type (e.g., meadow, desert steppe) under the same standard (reflectance), overestimation or underestimation will occur, thus restricting the estimation accuracy. To balance the overall trend of the transect and to take into account regional differences, we increased the weight information of spectral reflectance features of different grassland types and expanded the regression model of element content to Equation (6).

$$Y=\alpha \times RS+\delta ,$$

(5)

$$Y=\alpha \times \left(W_i\times RS\right)+\delta ,$$

(6)

where Y, RS, and W_i are the total element content per unit area, the reflection features, and the weights of different grassland types, respectively, and $\alpha$ and δ are the regression coefficient and error, respectively.

Taking the grassland type as the weight can alleviate the spectral drawbacks of “same subject with heterospectrum” and “foreign body with spectrum”, to some extent, and better expresses the characteristic differences between regions. However, it is hard to clearly define the range of grassland types in a complex environment, and the grassland types and AGB are largely affected by the climate environment. Related studies have shown that water is a key climatic factor restricting the growth of grassland vegetation in the Inner Mongolia Plateau [34]. Although biomass is an ideal weight factor, there are two problems when using biomass at present: (i) Biomass, itself, is also an unknown quantity, which needs to be estimated using remote sensing data; and (ii) even in the local area, the plant biomass is uneven, such that the instability of biomass as weight can easily lead to different weights for the same grassland type. Therefore, in this paper, we combine the humidity as a weight with the spectral image in order to comprehensively express the regional characteristics from spectral and climatic dimensions (Figure 4), and $\mathrm{RS}\times {\mathrm{W}}_i$ is represented by Equations (7) and (8), as follows:

$$\mathrm{FS}=\mathrm{RS}\times {\mathrm{W}}_i=\mathrm{ES}\times RH{I}_{nor}\times 10,$$

(7)

$$RH{I}_{nor}=\frac{RHI-RH{I}_{min}}{RH{I}_{max}-RH{I}_{min}},$$

(8)

where FS is the fusion spectra of spectral features and meteorological data, RHI_nor is the normalized value of relative humidity indices (RHI_max and RHI_min are the maximum and minimum RHI, respectively) and, to prevent the FS value from being too low and increasing the difference, the final result is multiplied by 10.

The relative humidity index (RHI) is a comprehensive expression of the climate environment, which integrates air temperature, precipitation, air pressure, wind speed, light, and other factors. Compared with a single meteorological factor, it is more suitable for the expression of grassland plant growth environment. Its calculation formula is as follows:

$$\mathrm{RHI}=\frac{\mathrm{P}-\mathrm{PET}}{\mathrm{PET}},$$

(9)

where P is the precipitation (mm) and PET is the potential evapotranspiration (mm). PET can be calculated by the Penman–Monteith method [39]:

$$\mathrm{PET}=\frac{0.408\Delta \left({\mathrm{R}}_{\mathrm{n}}-\mathrm{G}\right)+\mathsf{\gamma }\frac{900}{{\mathrm{T}}_{\mathrm{mean}}+273}{\mathsf{\mu }}_2\left({\mathrm{e}}_{\mathrm{s}}-{\mathrm{e}}_{\mathrm{a}}\right)}{\Delta +\mathsf{\gamma }\left(1+0.34{\mathsf{\mu }}_2\right)},$$

(10)

where R_n is the surface net radiation (MJ·m⁻¹·d⁻¹), G is the soil heat flux (MJ·m⁻²·d⁻¹), T_mean is the daily average temperature (°C), μ₂ is the wind speed at 2 m height (m/s), e_s is the saturated vapor pressure (kPa), e_a is the actual vapor pressure (kPa), Δ is the saturated vapor pressure curve slope (kPa·°C⁻¹), and γ is the dry and wet surface constant (kPa·°C⁻¹).

2.3.4. Spectral Features Enhancement

Fractional differential (FD) is one of the most commonly used spectral transformation methods in recent years, which extends the traditional first- and second-order derivatives to the non-integer-order differential, thus overcoming the data loss problem caused by complementary integer-order differential transformations [25,40]. The FD method has been applied to the fields of soil chemical element content [41], water nitrogen concentration [42], chlorophyll estimation [43], and other fields, in order to improve the estimation accuracy of the target. However, multi-spectral data have less spectral information and, as such, it is difficult to take advantage of the FD algorithm. After data fusion, FS has more spectral characteristics and can be combined with the FD algorithm, which is used to mine the weak sensitive spectral information related to C, N, and P, and to improve the accuracy of elemental estimation. Based on the Grünwald–Letnikov algorithm, the calculation formula can be expressed as [44]:

$$d^{v}f\left(x\right)=\underset{h\to 0}{\lim }\frac{1}{h^v}{\displaystyle \sum }_{m=0}^{\frac{t-v}{h}}{\left(-1\right)}^m\frac{\mathsf{\Gamma }\left(v+1\right)}{m!\mathsf{\Gamma }\left(v-m+1\right)}f\left(x-mh\right),$$

(11)

$$\mathsf{\Gamma }\left(z\right)={{\displaystyle \int }}_0^{\infty }{e}^{-u}{u}^{z-1}du=\left(z-1\right)!,$$

(12)

where v is the order, h is the differential step size, $\mathsf{\Gamma }\left(·\right)$ is the Gamma function, and t and v are the upper and lower limits of the differential, respectively. If the spectral sampling interval is resampled to 1 nm (h = 1), the formula can be further simplified to:

$$\begin{array}{c}\frac{d^{v}f\left(x\right)}{d{x}^v}\approx f\left(x\right)+\left(-v\right)f\left(x-1\right)+\frac{\left(-v\right)\left(-v+1\right)}{2}f\left(x-2\right)+\cdots \\ +\frac{\mathsf{\Gamma }\left(-v+1\right)}{n!\mathsf{\Gamma }\left(-v+n+1\right)}f\left(x-n\right),\end{array}$$

(13)

where n is the difference between the upper and lower limits of the differential and, if the order is 0, it means that the function itself has not been differentiated.

2.3.5. Model Construction and Accuracy Assessment

To evaluate the estimation accuracy of fused spectra for C, N, and P elemental contents in grassland transects over the Inner Mongolia Plateau, the weak sensitive information related to elemental content was extracted by use of the FD transformation method with raw spectra and fused spectra, and the sensitive bands were screened using the sequence forward band selection method. Then, the partial least-squares regression (PLSR) model was used to estimate C, N, and P contents. The coefficient of determination (R²), the root mean square error (RMSE), and the relative percent difference (RPD) between the measured values and the estimated values were used to evaluate the ability of RS and FS to estimate the contents of plant elements in the transect of Inner Mongolia, and the accuracy of the models was tested using a 10-fold cross-subsampling approach.

3. Results

3.1. Correlation Analysis of Climatic Factors and C, N, and P

The environment is one of the important factors affecting plant growth; in particular, the climate factors of ecological monitoring at the large regional scale cannot be ignored. For example, temperature, precipitation, topography, location, and land type are of great value when estimating the C, N, and P contents of leaves [34]. The precipitation, average temperature, and elevation information for the forage grass-growing season (April–July) at the 10 sampling areas, from east to west according to their longitudinal positions, are shown in Figure 5a,b. The region is characterized by the monsoon climate of medium latitudes. The warm and humid southeast wind from the Pacific Ocean in summer makes the windward slopes Z1, Z2, and Z3 rich in precipitation and causes them to have suitable temperature. Going further inland, the effect of the continental climate is significant and the precipitation decreases gradually. Although the average temperature at Z4–Z10 was slightly lower than that of Z1–Z3, being affected by altitude, the difference was small (<4 °C), indicating that the hydrothermal conditions in the eastern region are better than those in the western region.

Temperature and precipitation are the most important climatic factors affecting plant growth, as well as important input parameters in remote sensing estimation models for quantities such as biomass and elemental content; however, these data are often isolated, and the relative humidity index is the ratio of evaporation to precipitation, calculated by complex physical models based on a variety of meteorological factors. Therefore, humidity can express the climatic characteristics of the region better, and its correlations with C, N, and P were higher than those of temperature and precipitation (see Figure 5c). The relative humidity data of the study area were calculated using the Kriging interpolation method with meteorological stations data in the ARCGIS10.1 software. The data of 10 meteorological stations were randomly selected for verification. The results showed that the correlation between measured and interpolated humidity was 0.88, and the RMSE was 0.11 (

Table 3. Comparison of moisture indices between Kriging interpolation results and calculated result by measured data at 10 meteorological stations.

Station Site	Measured Moisture Index	Kriging Interpolation Moisture Index
1	−0.61	−0.53
2	−0.07	−0.03
3	−0.43	−0.38
4	−0.69	−0.64
5	−0.63	−0.65
6	−0.48	−0.47
7	−0.78	−0.53
8	−0.12	−0.31
9	−0.57	−0.53
10	−0.30	−0.39
Evaluating indicator	R = 0.88	RMSE = 0.11

). Therefore, it was considered feasible to obtain the humidity characteristics of the whole study area by Kriging interpolation. The spatial distribution of sampling points, meteorological stations, and humidity in the Inner Mongolia Plateau transect is shown in Figure 5d. The larger the humidity value, the more humid the region. Over the transect, it becomes gradually drier from east to west, affecting the yield of grassland per unit area and nutrient contents. However, it is clear that climate factors have a significant impact on plant growth at a large scale, while climate characteristics at the local scale are similar and insensitive to plant growth. The spectral characteristics can make up for a lack of climate data, facilitating realization of the physiological and biochemical detection of plants in the region. Therefore, the combination of remote sensing data and meteorological data is helpful in improving the monitoring level over large-span grassland transects.

3.2. Comparison of Raw Spectra, Enhanced Spectra, and Fusion Spectra

Fusion spectra (FS) consist of two parts: an enhanced spectrum (ES) and a weight factor, expressed in terms of the relative humidity index (RHI). The main purpose of the ES is to improve the information of the raw spectra (RS), while the RHI mainly provides the characteristic difference information between regions. Thus, the FS can estimate the C, N, and P contents from both the spectral and environmental dimensions.

The differences between RS and ES, RealHRS (measured real hyperspectral by ASD), and ES-Resample (resampled ES by the spectral response function of Sentinel-2) are shown in Figure 6. RealHRS and RS represent the spectral features at sample scale (1 × 1 m) and pixel scale (10 × 10 m), respectively, and have two main differences: (i) RealHRS has more spectral information than RS, which is limited by its spectral resolution; and (ii) the reflectance of RealHRS is higher than RS in the green and infrared ranges. ES was constructed by attaching the spectral detail information of RealHRS to RS, based on the abundance of ground objects without changing the spectral features of the pixel region. Thus, it can be seen that ES has more spectral information than RS, and the trends of ES and ES-Resample were consistent with that of RS. Further, the spectral angle (SA), spectral distance (SD), and correlation (R) between RS and ES-Resample, ES, and RealHRS were calculated (

Table 4. Comparison of spectral angle (SN), spectral distance (SD), and correlation (R) between raw spectra and resampling spectra, measured hyperspectral data, and enhanced spectra.

		RS-ES_Resample			ES-RealHRS
Site	SN	SD	R	SN	SD	R
Z1-1	0.0145	0.0097	0.9994	0.1282	5.1279	0.9823
Z1-2	0.0137	0.0091	0.9995	0.0925	5.1148	0.9881
Z1-3	0.0174	0.0115	0.9992	0.1006	5.0035	0.9894
Z1-4	0.0103	0.0068	0.9997	0.0699	5.0221	0.9964
Z1-5	0.0165	0.0110	0.9993	0.0740	5.0070	0.9944
Z1-6	0.0043	0.0026	0.9999	0.1353	4.6826	0.9853
Z1-7	0.0100	0.0064	0.9997	0.0873	4.7217	0.9934
Z1-8	0.0072	0.0046	0.9999	0.0591	4.7055	0.9937
Z2-1	0.0128	0.0097	0.9997	0.1313	5.7027	0.9815
Z2-2	0.0090	0.0068	0.9998	0.1305	5.7138	0.9836
Z2-3	0.0323	0.0260	0.9982	0.1268	5.7158	0.9781
Z2-4	0.0121	0.0093	0.9997	0.1148	5.8078	0.9844
Z2-5	0.0101	0.0079	0.9998	0.1075	5.4887	0.9891
Z2-6	0.0116	0.0087	0.9998	0.1376	5.6971	0.9772
Z2-7	0.0174	0.0132	0.9995	0.1398	5.7303	0.9727
Z2-8	0.0084	0.0063	0.9999	0.1254	5.6746	0.9807

). The SA and SD between RS and ES-Resample tended to 0, and the correlation tended to 1, indicating that ES generally conforms to the spectral features of the pixel region of RS. The SA between ES and RealHRS was small, and the correlation was high, indicating that ES has a large number of similar spectral details to RealHRS. Thus, ES can make up for the lack of spectral information of Sentinel-2, in order to improve the sensitivity of RS for the estimation of grassland traits. Meanwhile, ES, with rich spectral details, can give full play to the advantages of the multi-scale spectral transformation method for weak information mining.

FS was further supplemented by meteorological information affecting vegetation growth based on ES, and the characteristic differences between the RS and FS for the 10 sampling areas are shown in Figure 7. The wavelength range of the broadband multispectral Sentinel-2 imagery is shown in the gray background and, due to the spectral resolution limitation, RS lacks the necessary spectral information in the imaginary line part, such that the sensitivity of the broadband spectral information is weaker than that of the narrowband when many plant physiological and biochemical parameters are involved. In the actual natural environment, Z1–Z10 comprise a transition between meadow–typical grassland–desert grassland vegetation types, and the corresponding RS also reflects the characteristic changes, from typical vegetation spectra (Z1, Z2) to partial bare soil spectra (Z9, Z10). However, as shown in Figure 7a, the spectral distinction between different sampling areas was not significant in the range of green light reflection peak (560 nm) and near-infrared reflection band (740–850 nm). As shown in Figure 7b, FS had more spectral information than RS—this information was mainly supplemented by ES and tended to be the main ground object type of unit pixel. In addition, Z1–Z10 in FS were more significant than RS, as FS also integrates meteorological data related to plant growth, which enables FS to better express grassland characteristics at the regional scale.

3.3. The Distribution of Weak Spectral Information in Different FD Transform Scales of Fused Spectra

As the spectral fusion algorithm combines meteorological data (e.g., temperature, precipitation, wind speed, atmospheric pressure), surface type, and end-member hyperspectral characteristics, FS can better express regional characteristics. FS and RS were used to analyze the correlations with C, N, and P elements, as shown in Figure 8. Compared with RS, FS has two advantages: (i) FS has more spectral information than RS, and the amount of information depends on the spectral resolution of pure end-members. In this paper, the hyperspectral reflectance of pure end-members was measured by ASD, and there is a lot of narrowband information in the 325–1075 nm range, which means that more information can be used to construct the estimation model; and (ii) considering the single-band correlation, the best correlations of FS-C, FS-N, and FS-P were 0.84, 0.81, and 0.84, respectively, higher than those of RS-C, RS-N, and RS-P (0.78, 0.75, and 0.72, respectively), indicating that the spectral fusion data were better than the raw spectra data, making it more suitable for grassland transect monitoring.

Although the fused spectra had more spectral information, the high redundancy between adjacent bands limits the accuracy of the element estimation model. The fractional differential spectral transform method can reduce spectral redundancy and background noise, mining the weak spectral information sensitive to C, N, and P. We carried out correlation analysis between FS (325–1075 nm) and the elements C, N, and P with 0–4 fractional-order features (Figure 8b–d). It can be seen that the redundancy of adjacent bands was high in the 0–1.5th-order transforms. When the correlation R was higher than 0.78, the C-sensitive bands were mainly distributed in 709–813 nm and 856–1075 nm, the N-sensitive bands were mainly distributed in 725–805 nm and 858–1075 nm, and the P-sensitive bands were mainly distributed in 710–802 nm and 856–1075 nm. The difference between bands increased after the 2–4th-order transformation. The highly sensitive bands were concentrated in a relatively narrow wavelength range, achieving the effect of redundancy removal.

3.4. C, N, and P Content Estimation Model

To verify the effectiveness of FS in estimating C, N, and P contents, the optimal band combination at a single scale was selected using the sequential forward band selection method after the FD spectral transform, with FS and RS for comparison. The PLSR element estimation models were constructed using the ten-fold cross-validation method. The fitting results of the estimated C, N, and P contents and the actual measured values in the verification set are shown in Figure 9. The horizontal coordinate is the measured value, the ordinate is the corresponding estimated value, the red dashed line is a 1:1 straight line, and the blue real line is the fitting trend line between the measured and estimated values. It can be seen that the estimation results of FS were better than those with RS, where the estimation accuracy (R²) for C, N, and P elements was increased by 0.14, 0.14, and 0.18 over RS, respectively. Among them, the optimal estimation accuracy of FS for C was 0.88, the RMSE was 7.57 g/m², and the estimation of N content was the lowest (R² = 0.78). In addition, in the FS estimation model, the fitting line tended to the straight line of 1:1, but it was likely to be underestimated in the high-content region, due to the fewer sampling points. In the RS estimation results, the highest accuracy for the C element was R² = 0.74, RMSE = 10.99, and the R² values for N and P ware both lower than 0.67. The RPDs of the FS model were all greater than 2.0, while the RPDs of RS were between 1.6 and 2.0, indicating that the estimation accuracy of Sentinel-2 imagery can be effectively improved through use of the spectral fusion algorithm and fractional differential transformation.

To verify the effectiveness of fusion data in improving the estimation accuracy of C, N, and P contents, the models were constructed using RS, ES, and measured hyperspectral data (realHRS) without meteorological data and RS-M, FS, and realHRS-M features, with fused meteorological weight data. The estimation results on the validation set are shown in

Table 5. Accuracy evaluation of C, N, and P content estimation with or without meteorological weight in different characteristic models.

Models	C		N		P
	R2	RMSE (g/m2)	R2	RMSE (g/m2)	R2	RMSE (mg/m2)
RS	0.74	10.99	0.64	0.65	0.67	42.65
RS-M	0.73	11.15	0.65	0.64	0.73	38.25
ES	0.85	8.37	0.73	0.57	0.80	33.00
FS	0.88	7.57	0.78	0.51	0.85	28.92
RealHRS	0.93	5.54	0.93	0.30	0.92	20.49
RealHRS-M	0.96	4.52	0.94	0.26	0.95	15.76

Note: RS, ES, RealHRS are raw spectra, enhanced spectra, measured hyperspectral data, respectively; RS-M, FS, and RealHRS-M are the fusion features of meteorological weights with the raw spectra, enhanced spectra, and measured hyperspectral data, respectively.

. The results indicate that, except for the estimation accuracy of C in the RS model being slightly higher than that in RS-M, the other models had better estimation accuracy after integrating meteorological data. In general, the estimation accuracy, from high to low, was realHRS-M > realHRS > FS > ES > RS-M > RS. Thus, meteorological weight data play an important role in improving the estimation accuracy of elemental content in large-span grassland transects. The estimation accuracy of FS was lower than that of realHRS and realHRS-M, which may be due to two reasons: (i) The realHRS spectra are the measured data at the quadrat scale, which is consistent with the element content scale of the sampling point, while FS is at the pixel scale (10 × 10 m) of Sentinel-2; and (ii) although FS increases spectral detail information through the spectral simulation algorithm, its accuracy is still insufficient compared with the measured hyperspectral data. The RS and meteorological data are important factors for constructing the fused spectra, and the random forest method was used to sort the importance of RS and relative humidity. The importance scores are shown in

Table 6. The importance of relative humidity vs. raw spectral features for estimating C, N, and P contents.

Elements	Relative Humidity	Raw Spectra Features
		B1	B2	B3	B4	B5	B6	B7	B8	B8a	B9	B11	B12
C	0.71	0.34	0.32	0.52	0.41	0.38	0.26	0.36	0.31	0.26	0.17	0.52	0.87
N	0.59	0.27	0.33	0.44	0.35	0.37	0.26	0.15	0.29	0.33	0.33	0.46	0.81
P	0.67	0.16	0.30	0.31	0.40	0.37	0.36	0.23	0.30	0.22	0.07	0.27	0.87

, from which it can be seen that the score of relative humidity was second in all the spectral features (after the B12 band). Therefore, FS combined with meteorological data and spectral imagery realized high-precision estimation of trace elements in the large-span transect.

3.5. Spatial Distribution and Trend of C, N, and P

As local observations cannot accurately reveal the pattern characteristics of an ecosystem, continuous and large-scale grassland observations have always posed difficulty in the field of ecology. Remote sensing techniques have become important tools for global environmental monitoring. Figure 9 proves the feasibility of the use of remote sensing data for the estimation of C, N, and P contents. The estimated contents of C, N, and P for the whole transect of Inner Mongolia are shown in Figure 10. Therefore, with the help of remote sensing technology, the grassland trait information of the whole transect can be quickly obtained, thus overcoming the disadvantages of traditional methods that only rely on discrete sampling to obtain C, N, and P content. For example, the spatial distributions of C, N, and P can be obtained in detail from the sub-graph (Z1). Although the spatial distributions of C, N, and P still differed, the biomass per unit area of the same grassland type (meadow) also varied, thus affecting the distributions of C, N, and P. To further explore the spatial distribution law of C, N, and P contents over the whole transect and to verify the estimation accuracy, quantitative analysis of Z1–Z10 was carried out at the point scale (measured and estimated value of sample points) and the local scale (estimated mean value within 1 km buffer of sample point), respectively. As shown in Figure 11, the point scale mainly includes the mean values of measured and predicted samples, and the local scale is the mean content in the buffer of 1 km around the sampling point. It can be seen that: (i) C, N, and P contents showed a consistent decreasing trend from east to west at both scales; (ii) the elemental content was related to the grassland type (meadow > typical grassland > desert grassland), consistent with the natural environmental characteristics that the hydrothermal conditions in the eastern Inner Mongolia Plateau are better than those in the western region; and (iii) there were some differences in elemental content between point scale and local scale. For example, the contents of C and N at the local scale were slightly lower than those at point scale, while the content of P showed the opposite trend. In addition, the predicted values were consistent with the measured values at the point scale, which again proves the feasibility of combining the fusion spectrum with the fractional differential method to estimate C, N, and P contents.

4. Discussion

4.1. Contribution of Information Enhancement Method to C, N, P Estimation

Remote sensing is an effective means to quickly obtain information about C, N, and P contents, along with other important traits, but there are also many new difficulties in the monitoring of large-span grassland transects. For example, the ability to monitor chemical elements is insufficient, due to limitations related to the spectral resolution, different regional characteristics in the grassland transect, and weak sensitivity of spectra to C, N, and P. As shown in Figure 9, the information enhancement algorithm (including spectral enhancement, regional enhancement, and information enhancement) effectively improved the accuracy of the original Sentinel-2 imagery to estimate the C, N, and P contents over a large-span grassland transect.

Spectral enhancement is necessary for estimating C, N, and P contents, as the spectra are weakly sensitive to C, N, and P elements, and the spectral information of wideband multispectral images is insufficient. Multispectral satellites have the advantages of being easy to obtain, low cost, and high spatial resolution, which are of great significance in the monitoring of large-span transects; thus, it is important to develop transect monitoring methods suitable for multi-spectral images. For example, Zhou et al. [45] have attempted to use Sentinel-2, Sentinel-3, Landsat-8, and other data sources to jointly estimate soil C:N and soil organic matter content, but the obtained accuracy was low (R² < 0.5). Lucie Homolova [15] stated that the biggest shortcoming of multispectral satellite images in monitoring vegetation traits is insufficient information, which may have also been the reason for Zhou’s low estimation accuracy of the soil C:N ratio. However, spectral simulation spectra can increase the spectral details of the raw Sentinel-2 imagery (Figure 7), and have been successfully applied to grassland biomass estimation [30]. Furthermore, transects often contain a variety of environments with different regional characteristics, and the contents of C, N, and P elements in plants are affected by climate (temperature and precipitation) and natural environment factors [46]. Under the influence of temperature and precipitation, the whole considered transect presents a transition from meadow to typical grassland to desert grassland from east to west, where the plant element content differs under different grassland types. If only the spectra are used as the estimation standard, it is easy to ignore the difference in grassland types, thus restricting the estimation accuracy (Figure 4). Therefore, the regional characteristics should be reasonably considered in the monitoring of large-span grassland transects.

FS has more information, following spectral enhancement and regional enhancement (Figure 7), and the estimation accuracy of C, N, and P contents by the model constructed using the FS was higher than that with RS (Figure 9). This was mainly because FS is composed of enhanced spectra and meteorological data fusion, which endows FS with two obvious advantages over RS: (i) Spectral fusion spectra have more spectral detail information, integrating multiple information such as surface type, surface feature abundance, and pure end-members; and (ii) regional meteorological data affecting plant growth were fused at the spatial scale, so FS is more suitable for expressing regional differences. By decomposing the mixed pixels of the raw image to obtain the abundance information of ground objects and taking this as the weight of the hyperspectral pure end-member obtained at the ground scale, the reflectance in the unknown wavelength was supplemented without changing the RS reflection characteristics. Therefore, FS has more abundant spectral information for the main types of ground objects, leading to advantages in the construction of vegetation indices and biomass inversion [30]. Meteorological data have always been an important factor affecting the growth of grasslands; in particular, water is the main factor restricting longitudinal change (meadow–typical grassland–desert grassland) of grassland in the Inner Mongolia Plateau [34]. Temperature, precipitation, sunshine, wind, and other factors are closely related to surface water evaporation [47]. Relative humidity can effectively quantify the relationship between surface rainfall and evaporation. Although at a large scale the correlations between relative humidity and C, N, and P contents were better than that of a single meteorological factor and single-band spectral feature (Figure 5 and Figure 8), the climate in adjacent regions was similar, in terms of local regional characteristics. Therefore, the fusion of spectral features and meteorological characteristics can enhance the quality of data and improve the estimation accuracy.

The fractional differential method is an effective technology to analyze spectral reflectance spectra [24], which can strengthen the weak absorption features of vegetation physiology and biochemistry, with respect to the spectra, and eliminate the noise caused by soil background, water absorption, and atmosphere [48]. The spectral transformation features at different scales can be obtained by the FD algorithm, which greatly improves the information content of spectral data and reduces the collinearity between bands (Figure 8c,d). Therefore, it is convenient to extract sensitive bands related to C, N, and P, and select the optimal band combination to estimate the elemental content. At present, the FD algorithm is mainly used to estimate the physiological and biochemical characteristics of N [49] and P [50], leaf water content [51], and soil element content [52] in hyperspectral data. However, as the main data of grassland transect monitoring, multispectral images show difficulty in utilizing the advantages of the FD algorithm, due to spectral resolution constraints. Therefore, combination of the FS and FD algorithms with more spectral information after spectral enhancement and fusion effectively improved the estimation accuracy of grassland traits. In addition, FD can obtain spectral information at multiple scales, but cannot determine the optimal scale related to plant traits in advance, which also increases the difficulty of the application of the fractional differential method, to a certain extent. How to adaptively select the optimal feature combination in multi-scale FD features to construct multi-scale estimation models remains to be explored.

4.2. Spatial Distribution Characteristics of C, N, and P

Point-scale survey data can be extended to the surface-scale by remote sensing, which breaks the limitations of the traditional sampling surveying, and can continuously obtain ecological environmental data over a region [53]. The estimation error of soil organic carbon (SOC) tends to increase as it is extended to the regional or landscape scale, and the coefficient of variation may reach 30% [54]; however, it is necessary for large-scale grassland environmental monitoring. As shown in Figure 10, the C, N, and P contents were estimated by fusion spectrum and fractional differential algorithm. The contents of C, N, and P elements in the transect of the Inner Mongolia Plateau (44° N) decreased from east to west (Figure 11), which was consistent with the climatic characteristics of rich hydrothermal resources in the eastern region (Figure 5a). With the deepening of the inland continental climate, the precipitation decreases gradually, which limits the growth of vegetation. FS is composed of spectral and meteorological data, which is conducive to the description of grassland characteristics from the two aspects of reflection characteristics and climate environment. Therefore, C, N, and P contents had high estimation accuracy, and the R² between measured and estimated values were 0.88, 0.78, and 0.85, respectively (Figure 9).

At present, the estimation of chemical elements is mainly based on hyperspectral data with abundant information. The measured and estimated results in Figure 11 tended to be consistent at the sample scale, indicating that Sentinel-2 imagery can also be used to estimate C, N, and P contents with high accuracy after data fusion, which can help in studying the relationships between the soil and plants. Soil is the main carrier of carbon storage, and remote sensing has become an important means to quickly obtain large-scale soil organic matter [55]. However, most of the soil is covered by plants, and it is difficult to directly estimate soil organic matter. The structural change of vegetation in terrestrial grazing ecosystems has a close coupling effect with carbon and nutrient flows [56] and, if the interaction between vegetation and soil can be quantified, the estimation accuracy of soil organic matter can be further improved. The ratio of carbon to key nutrients (i.e., nitrogen and phosphorus) can regulate and affect the carbon consumption and decomposition processes in ecosystems, which is a key indicator in grassland monitoring [57]. However, most large-scale environmental surveys still rely on manual data collection [58]. There is a certain correlation between the distribution of C, N, and P contents in grassland ecosystems and vegetation types (meadow > typical grassland > desert grassland) at the unit scale, which indicates that the reasonable distribution of sampling points in the estimation model is helpful in improving the estimation accuracy of nutrient elements. Overall, remote sensing is expected to greatly improve the efficiency of access to information in the ecological and environmental fields.

4.3. Limitations and Future Research

Although the information enhancement method could improve the estimation accuracy when using Sentinel-2 imagery (Figure 9), there were still some shortcomings. For example, in contrast to the farmland ecological environment of a single species, it is difficult to clearly distinguish the species in grassland environments with different vegetation types. Thus, in this paper, we only divided the abundance of grassland and bare soil in the transect and took key climatic factors affecting the grassland types in the Inner Mongolia transect as environmental characteristics, and the estimation model for C, N, and P contents was constructed considering spectral and environmental dimensions.

In addition, the quality of the remote sensing imagery itself and spatial scale problems are also important factors restricting the final estimation accuracy. Cloud occlusion is an important factor that restricts image quality, and there were many areas covered by clouds in the estimation of grassland transect, which can cause many ground data to be missing and abnormal results in model estimation. Thus, we utilized mask processing to reduce the influence of clouds. The time error between imaging time and sampling time also restricts the acquisition accuracy of grassland traits. The scale problem—that the ground survey data do not correspond to the spatial resolution of satellite images—is another key factor restricting the monitoring of grassland traits, and poses an urgent problem to be solved in the field of remote sensing [59]. However, there is no effective solution, and this problem can only be alleviated by reasonably controlling the ground sampling points.

A grassland is a complex ecosystem featuring soil–vegetation–climate interactions, in which soil is an important source of elements needed for plant growth [60]. Thus, the introduction of soil is an important direction to further improve the estimation accuracy of C, N, and P, and other chemical element contents in grasslands at sample scale.

5. Conclusions

In this study, multi-spectral satellite (Sentinel-2) imagery were used to estimate the C, N, and P contents over a large-span complex environment (consisting of meadow, typical steppe, and desert steppe) in the Inner Mongolia grassland transect. We evaluated the associated estimation accuracies. The problems of insufficient information of multispectral images, regional characteristic differences, and weak sensitivity of spectra to C, N, and P elements are important factors restricting the accuracy of element content estimation in large-span grassland transects. However, the proposed spectral enhancement algorithm can combine the spatial advantages of satellite imagery with the spectral details of end-member spectra, improve the multispectral spectral information of original satellite imagery, and make up for the lack of sensitivity of multispectral images for monitoring C, N, and P contents, to a certain extent, and the algorithm can be used to construct enhanced spectral images (denoted ES-Sentinel-2). A regional enhancement algorithm can be used to combine meteorological data with spectral data, in order to quantify regional differences in spectral and climatic dimensions, and to reduce the estimation error caused by grassland type differences. The feature enhancement algorithm can expand the original spectral information to multiple dimensions, which is conducive to extracting the spectral features related to C, N, and P.

The spectral features extracted by the information enhancement algorithm (including spectral enhancement, regional feature enhancement, and feature enhancement) were used to construct the partial least-squares regression model, and the estimation accuracies (R²) for C, N, and P were 0.87, 0.78, and 0.85, respectively. Compared with the estimation results of raw Sentinel-2 images, the RMSE was reduced by 3.42 g/m², 0.14 g/m², and 13.73 mg/m², respectively. The C, N, and P contents showed a decreasing trend from east to west over the transect, with the characteristics of meadow > typical grassland > desert grassland.

Therefore, the information enhancement method can improve the estimates of C, N, and P content over large-span grassland transects. At the same time, many areas can be improved; for example, we considered green vegetation in the flourishing grassland period in our model. As such, the estimation accuracy will be restricted by non-green vegetation in other growth stages. In the next stage, our research group will further study the estimation of typical grassland traits in multi-period and long-term sequences.