Image Spectral Resolution Enhancement for Mapping Native Plant Species in a Typical Area of the Three-River Headwaters Region, China

Abstract

Large-scale multispectral remote sensing data are often unavailable for some practical applications. Spectral resolution enhancement for large-scale multispectral remote sensing images by incorporating small-scale hyperspectral remote sensing images is an alternative way to generate remote sensing images with both large spatial range and high spectral resolution. This paper proposes an improved spectral resolution enhancement method (ISREM) using spectral matrix and weighting the spectral angle of the transformation matrix. ISREM is tested in a typical area of the Three-River Headwaters region (TRHR) to produce a synthetic hyperspectral image (HSI). Two existing spectral resolution enhancement methods, the color resolution improvement software package (CRISP) and spectral resolution enhancement method (SREM), are adopted to compare with ISREM. To further test the practicality of the synthetic HSIs generated by the ISREM, CRISP and SREM, they are used to estimate the coverage of native plant species (NPS) using support vector machines (SVM) and random forest (RF) regressions. The experimental results are as follows. (1) For the Pearson correlation coefficient between the synthetic HSI and original image, ISREM yielded the largest value of 0.9582, followed by CRISP and SREM with values of 0.9480 and 0.9514. For spectral similarity, the HSI generated by the ISREM was the closest to the original reference HSI in the spectral curve. It also showed the best cumulative performance with the use of the three quality evaluation indexes. (2) The identification accuracies of native plant species were 93.51%, 90.91%, 89.61% and 89.61% using generated HSIs and original multispectral image (MSI) within a threshold of 20%, respectively. Compared with original MSI, the synthetic HSI showed better ability to identify NPS in the study area, which further illustrated the effectiveness of the ISREM. (3) The ISREM can reduce the strict requirement of pure pixels and maintain the quality of synthetic HSI by spectral angle weighting. Hence, the proposed ISREM outperforms the existing CRISP and SREM methods in image spectral resolution enhancement of multispectral remote sensing images.

1. Introduction

Spectral enhancement for multispectral remote sensing data greatly improves the utilization efficiency of remote sensing images. The methods for enhancing remote sensing data to study grassland degradation in the TRHR mainly include three groups. (1) The first group of enhancement method for identifying grassland degradation is conducting band math of the remote sensing images, i.e., the vegetation index (VI) and biochemical index. The normalized difference vegetation index (NDVI) is one of the most common spectral enhancement methods for vegetation characterization. The VIs and band math expressions method has been widely used to estimate grassland coverage in the TRHR [1,2]. This kind of method can enhance the recognition ability of remote sensing images through the operation of wave bands, especially the ability to distinguish vegetation from the background. However, due to the limited number of bands, this method is weak for identifying vegetation species. (2) The second group of enhancement methods is image sharpening of HSI to obtain a finer spatial resolution for grassland species identification. The remote sensing image sharpening method generally refers to the fusion of panchromatic images and MSI to enhance the spatial resolution of the MSI [3,4,5]. The enhancement of the spatial resolution of remote sensing data has received considerable attention, but it is rarely applied to grassland vegetation species identification. The most commonly used enhancement method is the classic Gram-Schmidt pan-sharpening method [6]. Gram-Schmidt pan-sharpening is usually used to improve the spatial resolution of HSI with a higher spatial resolution MSI in image preprocessing [7]. Although this method enhances the image recognition ability through an improved spatial resolution, it still has the problem of spectral information distortion, which is not conducive to the recognition of grass species. Through field investigation, we found that various types of grassland in the TRHR were mixed-living together, with low grass height and little difference in morphological characteristics. Thus, the improvement of spatial resolution is of little help for grass species identification; it is more important to improve the spectral resolution to identify grass species of the TRHR. (3) The third group of enhancement methods is spectral enhancement for MSI [8]. A model is first built between the MSI and HSI, and then a synthesized HSI is calculated pixel by pixel or using a fixed size filter [9,10]. This kind of method not only maintains the spatial resolution of MSI but also enhances the spectral resolution of MSI.

Spectral enhancement methods are mainly used to improve the spectral resolution of MSI with an HSI and have great potential in some applications requiring high spectral resolution images. The SREM was first proposed in 2014 and applied in mineral alteration information extraction [11]. However, it has not been applied in vegetation identification or grassland degradation as far as we know. Previous studies of grassland degradation have mainly focused on three aspects: (1) total vegetation or biomass based on a long time series to identify grassland degradation; (2) the use of vegetation indexes, such as VIs and Biological indexes (BIs), to assess grassland coverage; and (3) calculation of net primary productivity and local net primary productivity scaling to monitor grassland degradation [12,13]. However, grassland degradation is not only the degradation of biomass or coverage but also the change in grass species, which is manifested as an increase in noxious weeds (NWs) and a decrease in the dominant native grass species. NWs invasion has become serious in recent years and the extremely degraded black soil has attracted the attention of researchers [11]. The coverage of NPS in the TRHR has been proved to be a decreasing trend in recent decades [14]. Traditional methods have difficulty identifying NPS due to the slight spectral differences between NPS and invasive NWs and the low spectral resolution of MSI for large areas. Thus, specific and effective data that can better distinguish between NPS and invasive NWs should be found and utilized. HSI has been widely applied in mineral resource exploration, inland water remote sensing, and vegetation ecology and agriculture, and has been proven to have a good ability to distinguish between similar spectral features [15]. However, to the best of our knowledge, few studies have used hyperspectral data to identify grass species in the alpine grassland of the TRHR because the temporal and spatial consistency of the samples and images is difficult to achieve.

HSI has sensitive spectral discrimination ability and can be applied to the identification of grassland species to realize the study of grassland degradation in the study area. However, satellite HSI is difficult to apply to the identification of grassland components in large regions due to the small width, low spatial resolution and scarce data source. Therefore, it is of great significance to study methods that can maintain the spectral resolution and expand the spatial coverage of HSI. Due to the important symbiosis of grass species in the study area, it was difficult to find evenly distributed grass species over a large area, and pure pixels were difficult to obtain in the image. Therefore, the SREM method could not be used directly for identifying grass species. Thus, a novel spectral enhancement method based on clustering and spectral angular weight was proposed to improve the spectral resolution of the MSI. Compared with SREM, the proposed ISREM makes two main improvements. (1) MSI were pre-clustered with the clustering method, and the result of the clustering was considered a “pure pixel”, which was used in the calculation. (2) The spectral angle was calculated and used as the weight for each cluster result, and then the spectral value of the predicted pixel was calculated. The detailed implementation of ISREM is described in Section 3.1.

The main contributions of this study are: (1) the proposal of the ISREM for the enhancement of MSI and generation of a synthetic HSI of the study area; (2) features extraction to map the distribution of NPS using the synthetic HSI, and the identification abilities of original MSI and synthetic HSI were compared; and finally, (3) a grassland degradation map was produced based on the distribution map of NPS. As a whole, this study provided a method for identifying the grass types and degradation of plateau grasslands using MSI in the case of HSI absence.

2. Materials and Methods

2.1. Study Area

This study area (Figure 1) is located in the Tongtian River reserve of Yushu Tibetan Autonomous Prefecture, which covers part of southeastern Zhiduo County in Qinghai, China. This area was selected from the hinterland of the TRHR, which is the source region of the Yangtze, Yellow, and Lancang Rivers [11]. The annual average temperature is −0.4°C, with an annual rainfall of ~394 mm/year in the TRHR. The average altitude of the study area is above 4000 m. The climate in this region is dry, cold in winter and cool in summer, with strong solar radiation and thin air. The oxygen content of the air is about 14% in the study area. The cold season lasts nearly 10 months (January–June & September–December), with a large temperature difference between day and night. Approximately 70–90% of the study area is covered by vegetation, and the alpine meadow is the main grassland type in this region [16].

2.2. Grass Species and Samples Collection

In the study area, the grass species mainly include two kinds of NPS (Kobresia and Stipa capillata) and five kinds of NWs (Leontopodium nanum, Potentilla chinensis Ser, Bryophyte, Polygonum sibiricum Laxm, and Lamiophlomis rotata) [17]. NPS were the native plant community, while NWs have become the secondary vegetation due to overgrazing and climate changes in the TRHR. The invasion of NWs has been accelerating due to their adaptive and reproductive capacity for decades. Because of the pungent taste and toxic of NWs, yak and Tibetan sheep prefer NPS to NWs. In the past, for a long time, the coverage of NPS, suitable for grazing, was decreasing, which is an important indicator of grassland degradation [18]. Therefore, mapping the coverage of NPS is critical for investigating the grassland degradation in the TRHR.

The main objective of the sample collection is to collect the spectra of grasses and samples of grassland coverage. Sampling was carried out from 10 to 19 August 2017 for the first time. The second time, a field survey was conducted from 12 to 23 August 2019. Despite a two-year interval of the two surveys, the changes in spectra and coverage of grass species are negligible [19]. The Grass spectral data were collected using the spectrum-measuring instrument HR-512i, a portable object spectrometer produced by the SVC Company, USA. Spectrum acquisition was performed at 11:00–14:00 on sunny days (Figure 2a). For each kind of grass, at least 10 spectrums were collected from 5 different high-density growing areas and combined to 1 average spectrum. The spectrums of Kobresia, Stipa capillata, Leontopodium nanum, Potentilla chinensis Ser, Bryophyte, Polygonum sibiricum Laxm, and Lamiophlomis rotata were finally generated by SVC HR-512i PC software.

For the first filed sampling, 113 30 m × 30 m, the same size of the pixel of HJ-1A MSI, grassland coverage samples were collected. A gentle slope and relatively homogeneous of grass species were necessary to select the sampling area. For each 30 m × 30 m sample, nine quadrats with a size of 1 m × 1 m were established and distributed in an “X” shape, with one quadrat at the center and the other eight quadrats in the “X” directions (Figure 2b). Each quadrat unit was comprised of bare soil, PLS and NWs in the research area. The relationship of soil, PLS and NWs can be described as: $C_{soil}+C_{PLS}+C_{NWs}=1$. $C_{soil}$ is the proportion of soil coverage. $C_{PLS}$ respected the proportion of PLS coverage and $C_{NWs}$ respected the proportion of NWs. For each quadrat unit, coverage information was recorded as Figure 2d. The location information was acquired from Trimble Geo 7X hand-held submeter GPS units with an accuracy of ~0.5 m in the center of each quadrat. Finally, the average coverage $C_{soil}$, $C_{PLS}$ and $C_{NWs}$ of the nine quadrats were calculated, respectively, as the $C_{soil}$, $C_{PLS}$ and $C_{NWs}$ values of the 30 m × 30 m sample. The location information of the 30 m × 30 m sample was the same as center quadrat (NO.1 in Figure 2b).

The second sample collection was implemented using the same equipment and methods. A total of 142 new 30 m × 30 m grassland coverage samples were collected in the second survey. Altogether, 255 30 m × 30 m grassland coverage samples were obtained in the study area and used in this paper.

2.3. Image Data Acquisition and Processing

The MSI and HSI used in this paper were provided by the Chinese Environmental Satellite HJ-1A, which is equipped with a CCD camera and an HSI camera. According to the growth period of the vegetation in the study area, one MSI and two HSIs were acquired from the China Centre for Resources Satellite Data and Application (http://www.cresda.com/CN/) on 10 August 2019. The HJ-1A CCD MSI provides 4 visible spectral bands with a spatial resolution of 30 m. The HJ-1A HIS provides 115 bands from 460.04 nm to 951.54 nm with a spatial resolution of 100 m (Table 1).

Atmospheric correction for the HJ-1A MSI/HSI 2A class imagery was carried out with the FLASH module in ENVI 5.3. A vertical stripe removal method [20] based on image rotation was utilized to remove the stripe in band 1 to band 25(410 nm to 490 nm). Terrain correction [21] was then carried out using ASTER GDEM in ENVI 5.3. A G-S fusion was used to improve the spatial resolution of HJ-1A HSIs. Finally, HSI and MSI were corrected to the same spatial reference as the sample data.

2.4. Spectral Resolution Enhancement of HJ-1A MSI Imagery

2.4.1. Basic Model

Spectral enhancement of the MSI using HSI should be based on the following assumptions: (1) Because they use the same physical imaging mechanism, both hyperspectral sensors and multispectral sensors have the same spectral response trend in specific wavelengths for imaging the same ground object. (2) The spectral information of ground objects has spatial scale invariance, that is, the spectral characteristics of the same object in different spatial positions are basically the same. (3) The HSI overlaps to some extent with the MSI to be enhanced, and the features in the overlapped region are consistent with the features in the nonoverlapped region to be enhanced [22,23].

MSI and HSI are extensions of two-dimensional images in spectral latitude and can be regarded as three-dimensional images [24]. To facilitate calculation, we converted the three-dimensional images into two-dimensional matrices, each column of which represented pixels, and the rows represented bands. H was assumed to be an HSI with s samples, c columns and k bands. A k × n matrix was converted for H, where N is calculated by s × c. H can be expressed as:

(1)$H=\left(\begin{array}{cccc}{H}_{\left(1\right)}& H_{\left(2\right)}& \dots & H_{\left(n\right)}\end{array}\right)=\left[\begin{array}{cccc}{h}_{\left(1,1\right)}& h_{\left(1,2\right)}& \dots & h_{\left(1,n\right)}\\ h_{\left(2,1\right)}& h_{\left(2,2\right)}& \dots & h_{\left(2,n\right)}\\ \dots & \dots & \dots & \dots \\ h_{\left(k,1\right)}& h_{\left(k,2\right)}& \dots & h_{\left(k,n\right)}\end{array}\right]$

where k represents the band numbers of the HSI and n represents the pixel number in the HSI. H(i) (i is in the range from 1 to N) represents the column vector of the spectral values, and h(i,j) (j is in the range from 1 to k) represents the value of the i-th pixel in the j band. In the same way, the MSI M can be expressed as

(2)$M=\left(\begin{array}{cccc}{M}_{\left(1\right)}& M_{\left(2\right)}& \dots & M_{\left(n\right)}\end{array}\right)=\left[\begin{array}{cccc}{m}_{\left(1,1\right)}& m_{\left(1,2\right)}& \dots & m_{\left(1,n\right)}\\ m_{\left(2,1\right)}& m_{\left(2,2\right)}& \dots & m_{\left(2,n\right)}\\ \dots & \dots & \dots & \dots \\ m_{\left(l,1\right)}& m_{\left(l,2\right)}& \dots & m_{\left(l,n\right)}\end{array}\right]$

where l represents the band numbers of the MSI and N represents the pixel number in the MSI. M(i) (i is in the range from 1 to n) represents the column vector of the spectral values, and m(i,j) (j is in the range from 1 to k) represents the value of the i-th pixel in the j band. For the HSI and MSI in the same space and under the same sensor, MSI can be regarded as the linear representation of HSI [10]. The HSI H is related to the MSI M through the use of a simple filtering matrix F as

(3)$M=FH+\epsilon$

where F is a l × k matrix, which is the coefficient of the linear combination of hyperspectral bands into MSI. ε is the Gaussian random error. This inverse process can also be written as

(4)$GM=H+\gamma$

where G is a k × l matrix that also acts to approximately reconstruct the HSI. γ is the Gaussian random error. According to the least square principle, Formula (4) can be used to solve the expression of transformation matrix G for constructing HSI as

(5)$G=H{M}^T{\left(M{M}^T\right)}^{-1}$

To calculate G in (5), there must be more pixels than the number of multispectral bands, which is an easy assumption to satisfy. With a transformation coefficient matrix, the artificial hyperspectral value of the ith pixel is then obtained with the following equation:

(6)$H_{\left(i\right)}^{\prime }=G{M}_{\left(i\right)}$

where $H_{\left(i\right)}^{\prime }$ is the ith pixel of the estimated hyperspectral vector and $M_{\left(i\right)}$ is then combined to construct an estimated HSI as

(7)$H^{\prime }=\left[\begin{array}{cccc}{H}_{\left(1\right)}^{\prime }& H_{\left(2\right)}^{\prime }& \begin{array}{ccc}\cdots & H_{\left(i\right)}^{\prime }& \cdots \end{array}& H_{\left(n\right)}^{\prime }\end{array}\right]$

2.4.2. Clustering for Extraction of Spectral Matrix of HSI and MSI

Theoretically, pure pixels are the most ideal for constructing a transformation matrix G of MSI and HSI. However, limited by the low spatial resolution of the images and irregular distribution of ground objects, endmember extraction of pure pixels is a difficult and complicated task, and the results are easily influenced by the image data [25]. Especially in the grassland area of the plateau, because of the low grassland coverage and mixed grassland communities, it is more difficult to extract pure pixels. Thus, clustering, a less-than-ideal alternative method, was then adopted, and k-means clustering was utilized for the pre-classification of the MSI in this study. Supposing that the MSI was clustered into x classes, no less than l (number of multispectral bands) pixels were extracted from both MSI and HSI. For each clustering class, two groups of pixels were used to construct the spectral matrices $Mq$ and $Hq$ (q ranges from 1 to x) in Formulas (8) and (9). Altogether, x pairs of spectral matrices of multispectral and HSI were finally extracted by clustering.

2.4.3. Weighted Spectral Angle for Transformation Matrix G

As the estimated HSI $H^{\prime }$ was created from the MSI using a transformation matrix, the quality of the final product image $H^{\prime }$ depended on transformation matrix G. A global and local model for the generation of G was proposed in previous literature [10]. This model could be useful for rebuilding the spatial resolution of HSI by using Butterworth filters. The image scope of the final product was limited to the overlap area between the HSI and MSI. Beyond the range of the overlap area, this model will no longer be useful.

In fact, G is the mapping from the spectra of MSI pixels to the spectra of HSI pixels. According to the second assumption in Section 3.1.1, the spectral characteristics of the same object in different spatial positions were basically similar. The target region could be divided into x group objects through the visual interpretation of the HSI and MSI with expert experience. The same transformation matrix G can be used jointly by the same ground object. Therefore, x pairs of hyperspectral and multispectral matrices can be constructed according to the x group ground object, and the x group G-matrices can be obtained.

In most HSIs, mixed pixels are ubiquitous, and it is difficult to find pure pixels [26]. Therefore, the direct application of pure pixel objects for constructing HSI and MSI (such as SREM) is limited by the image itself, which makes it difficult to use widely. Therefore, the HSI and MSI in the study area can be first studied by the pre-classification. Then, x pairs of hyperspectral and multispectral matrices can be extracted from the pre-classification result. In this paper, the clustering method (i.e., k-means clustering) was used to obtain ‘pure pixel’ pairs. In the overlapping region of the HSI and MSI, the region is clustered into x classes. For each clustered class q, the H and M matrices are constructed for n groups of objects as shown in Formulas (8) and (9), respectively.

(8)$Hq=\left(\begin{array}{cccc}H{q}_{\left(1\right)}& H{q}_{\left(2\right)}& \dots & H{q}_{\left(n\right)}\end{array}\right)=\left[\begin{array}{cccc}h{q}_{\left(1,1\right)}& h{q}_{\left(1,2\right)}& \dots & h{q}_{\left(1,n\right)}\\ h{q}_{\left(2,1\right)}& h{q}_{\left(2,2\right)}& \dots & h{q}_{\left(2,n\right)}\\ \dots & \dots & \dots & \dots \\ h{q}_{\left(k,1\right)}& h{q}_{\left(k,2\right)}& \dots & h{q}_{\left(k,n\right)}\end{array}\right]$

(9)$Mq=\left(\begin{array}{cccc}M{q}_{\left(1\right)}& M{q}_{\left(2\right)}& \dots & M{q}_{\left(n\right)}\end{array}\right)=\left[\begin{array}{cccc}m{q}_{\left(1,1\right)}& m{q}_{\left(1,2\right)}& \dots & m{q}_{\left(1,n\right)}\\ m{q}_{\left(2,1\right)}& m{q}_{\left(2,2\right)}& \dots & m{q}_{\left(2,n\right)}\\ \dots & \dots & \dots & \dots \\ m{q}_{\left(l,1\right)}& m{q}_{\left(l,2\right)}& \dots & m{q}_{\left(l,n\right)}\end{array}\right]$

where $Hq$ is a k × n matrix that represents the hyperspectral matrix of class q (q represents one class in x, in the range from 1 to x) and $Mq$ is a l × n matrix that represents the multispectral matrix of class q. Note that different classes can have different values of n, but n has to be greater than or equal to l. According to Formula (5), the transformation matrix G of class q can be denoted as

(10)$Gq=HqM{q}^T{\left(HqM{q}^T\right)}^{-1}$

Then, a collection of transformation matrix G can be expressed as $\left\{G1,G2,\dots Gq,\dots ,Gx\right\}$. The G matrix collection could be regarded as all kinds of spectral interpolation from the MSI to the HSI. These G matrices were used to enhance the MSI pixel by pixel, but the best G matrix differs for different pixels due to spectral differences. For each pixel in the MSI, a weighted spectral angle (WSA) for obtaining the optimal transformation matrix G was calculated as

(11)$G{\left(i\right)}^{\prime }={\displaystyle \sum }_{q=1}^{x}WCV\left(q\right)·Gq$

where

(12)$WCV\left(q\right)=\frac{SAM{\left(q\right)}^2}{{{\displaystyle \sum }}_{q=1}^{x}SAM{\left(q\right)}^2}$

(13)$SAM\left(q\right)=arccos\left[\frac{⟨H{q}_{\left(i\right)}{}^{\prime },\overline{Hq}⟩}{{\left|\left|H{q}_{\left(i\right)}{}^{\prime }\right|\right|}_2·{\left|\left|\overline{Hq}\right|\right|}_2}\right]$

where $G{\left(i\right)}^{\prime }$ represents the optimal transformation matrix of each pixel in the MSI to estimate the hyperspectral vector $H_{\left(i\right)}^{\prime }$. WCV(q) is the ratio between the square of $SAM\left(q\right)$ and the summary of the square of all $SAM\left(q\right)$, which represents the weight coefficient value of the q-th G matrix in the G matrix collection. $SAM\left(q\right)$, calculated by Formula (13), is the spectral angle between the average hyperspectral vector and the q-th predicted hyperspectral vector, where $\overline{Hq}$ and $Hq\left(i\right)$ was calculated as

(14)$\overline{Hq}=\raisebox{1ex}{$\left[{{\displaystyle \sum }}_{i=1}^{n}Hq\left(i\right)\right]$}\!\left/ \!\raisebox{-1ex}{$n$}\right.$

(15)$Hq\left(i\right)=Gq{M}_{\left(i\right)}$

where $\overline{Hq}$ is the average spectral value of the q-th clustering class of the HSI metrics, which will not change with i. $Hq\left(i\right)$ is the temporary predicted hyperspectral vector using the matrix $Gq$ and the i-th pixel spectral vector $M_{\left(i\right)}$ of the MSI. For each pixel, the final predicted hyperspectral vector was calculated as

(16)$H_{\left(i\right)}^{\prime }=G{\left(i\right)}^{\prime }{M}_{\left(i\right)}$

Step by step, each optimal transformation matrix $G{\left(i\right)}^{\prime }$ can be calculated using this weighted spectral angle method, and then a hyperspectral vector $H_{\left(i\right)}^{\prime }$ can be calculated for each pixel of the MSI. Finally, a HSI with the same spatial extent as the MSI was combined by {$\begin{array}{cccc}{H}_{\left(1\right)}^{\prime }& H_{\left(2\right)}^{\prime }& \begin{array}{ccc}\cdots & H_{\left(i\right)}^{\prime }& \cdots \end{array}& H_{\left(n\right)}^{\prime }\end{array}$}.

2.4.4. Model Application

To produce an HSI image of the sampling area, five steps were executed as follows:

Step 1:K-means clustering: the shared region of the multispectral and hyperspectral data was extracted by the clipping method, and four classes (grass, water, soil and cloud/snow) of ground were clustered according to field investigation using K-means clustering of the HJ-1A MSI image.

Step 2:Extraction of spectral value matrices from the same pixel of HSI1 and MSI: the spectral value matrices were extracted from the same pixel of HSI 1 and MSI. For each class, sixty pairs of pixels were randomly extracted from the same pixel position on HSI 1 and MSI, and each set of ten groups were merged into a one-pixel pair by means of averaging to ensure that there were at least four (the band number of MSI) pairs of pixel spectral vectors. To reduce the running time of the program and ensure that there was adequate information, we chose 6 averaged pixel spectral pairs for each clustered class after several attempts.

Step 3:Creation of HSI and MSI matrix pairs of each cluster: HSI and MSI matrix pairs were created for the model input. For each clustered class, six pairs of spectra were merged into a matrix. Thus, the spectra of HSI formed a 115 × 6 matrix, and the spectral value of MSI was then merged into a 4 × 6 matrix. At this point, there were 4 groups of matrix pairs: H_(grass), H_(water), H_(soil), H_(cloud/snow) and M_(grass), M_(water), M_(soil) and M_{(cloud/snow).} Thus, the average hyperspectral vector $\overline{Hq}$ was calculated as the average value of H_(grass), H_(water), H_(soil), and H_(cloud/snow).

Step 4:Input Matrixes: Input the matrix pairs to the ISREM of Section 2.4.2, and four weighted G matrixes were generated as G_(grass), G_(water), G_(soil) and G_(cloud/snow) in the process. For each pixel of the MSI image, four temporary hyperspectral vectors were then calculated as H_(grass)(i), H_(water)(i), H_(soil)(i) and H_{(cloud/snow) (i),} and an optimal transformation matrix $G{\left(i\right)}^{\prime }$ was then calculated with Formula (11). A hyperspectral vector was generated using $G{\left(i\right)}^{\prime }$ and the multispectral vector of the i-th pixel of MSI. Finally, a synthetic HSI image with the same extent as the MSI image was predicted pixel by pixel according to Formula (7).

Step 5:Output of the predicted HSI image: to verify the results of the spectral fitting, HSI 2 was used for spectral comparison analysis with the synthetic images. The spectral curve of the same pixel in HSI 2 and MSI was first extracted, and the Pearson correlation coefficient between the overlapping area of HSI 2 and MSI was calculated for each band. Then, three fusion indexes were calculated to evaluate the ISREM results.

Finally, an estimated HSI of the study area, which was used to estimate the coverage of NPS and to investigate grassland degradation, were obtained. The workflow of ISREM in the research is shown in Figure 3.

2.5. Mapping Coverage of NPS

2.5.1. Feature Extraction

When the HJ-1A CCD MSI was enhanced in spectral resolution by the ISREM in this paper, the band numbers were increased from 4 to 115, and the spectral resolution was enhanced to ~5 nm. However, this does not mean that adequate spectral information can achieve better recognition ability. Previous studies have shown that the classification accuracy first increases and then decreases with increasing feature dimensions in the case of limited samples, which is called the Hughes phenomenon [27,28]. To obtain better regression results for NPS, the classification features were extracted from the original wave band, first-order spectral derivative, and VIs.

The spectral difference was essential for the classification of the vegetation [29]. To map the distribution of NPS, it is necessary to find bands that can distinguish NPS from NWs. Thus, spectral differences between NPS and NWs should be determined before classifying them. In our field investigation, an SVC hr-512i ground spectrometer with a spectral range from 350 nm to 1050 nm was used to measure the spectral curves of NPS and NWs. The spectral curves of NPS and NWs were shown in Figure 4a. The spectral reflectance of the NWs was generally higher than that of NPS, which might be due to the withering period of NPS in late August. From the spectral curves, we can see intuitively that NPS and NWs could be distinguished from each other.

In this study, the Mahalanobis distance analysis method [30] was used to select optimal bands for distinguishing NPS and NWs. The spectral reflectances of Kobresia and Stipa capillata were merged for NPS, and those of Leontopodium nanum, Potentilla chinensis Ser, Bryophyte, Polygonum sibiricum Laxm & Lamiophlomis rotate were merged for the NWs. The first-order derivative spectra were then calculated using the spectral math tool in ENVI 5.3. Then, the Mahalanobis distances of the original spectra and first-order derivative spectra between NPS and NWs were calculated.

As shown in Figure 4b,d, the average distance was chosen as the threshold of the selected band. For the original spectra, wavelengths ranging from 340 nm to 523 nm and 645 nm to 693 nm were the optimal band ranges. According to the band characteristics of the HJ-1A HSI, bands 1 to 28 and bands 65 to 75 were selected as classification features. Additionally, for first-order derivative spectra, spectral difference bands ranged from 689 nm to 755 nm. Bands 75 to 87 of the first-order derivative spectra were selected based on the band setting of the HJ-1A HSI. All selected bands of original spectra and first-order derivative spectra are shown in Table 2.

In addition to the selected original spectral and first derivative spectral bands, NDVI is widely used in vegetation growth monitoring, vegetation classification, and coverage modeling, and has a good ability to distinguish vegetation canopies and backgrounds [31,32]. The narrowband green normalized difference vegetation index (nGNDVI) and narrowband normalized difference vegetation index (nNDVI) are very sensitive to the greenness of vegetation [33]. Thus, the NDVI, nGNDVI and nNDVI were selected for the identification.

The biochemical and biophysical parameters within the plant leaves and the canopy, such as chlorophyll a and b, carotenoids and anthocyanin, also affected the spectral reflectance of vegetations. Biochemical indexes, such as narrowband leaf chlorophyll index (nLCI), anthocyanin reflectance index 2 (ARI2), and narrowband photochemical reflectance index (nPRI), has shown obvious differences in different grass species [34,35]. Biochemical indexes, nLCI, ARI2 and nPRI, were also selected. The final selected VIs are shown in Table 3.

2.5.2. SVM and RF Regression

The SVM classifier has been widely used in land use/land cover for classification [41,42,43,44,45]. Compared with other machine learning algorithms, it has the characteristics of high precision and strong generalization ability and is suitable for remote sensing image classification of high-dimensional samples [46]. All features selected in Section 2.5.1 were combined as input images for the SVM regression in EnMAP-Box2.0 [47]. Sample image for NPS was converted from vector to raster format. The sample images must maintain the same extent and spatial resolution as the feature images. The value of the sample images was the coverage ratio of NPS. A random sample tool was used to generate the training samples and validation samples, which were in proportions of 70% and 30%, respectively. The grid search method was utilized to parameterize the SVM regression. Gaussian RBF kernel parameters (g) were initialized with a minimum value of 0.01, a maximum value of 1000, and a multiplier of 10. The regularization parameter (C) was initialized with a minimum value of 0.01, a maximum value of 1000, and a multiplier of 10. Cross-validation was used to determine the optimal parameters. The final determined parameter values were a g value of 0.1 and a C value of 100 with an epsilon (ε) insensitive loss function parameter of 12.33.

Random forest (RF) was employed as a comparison method to estimate the coverage of NPS. Previous studies have shown that RF has been widely used in remote sensing classification [48,49,50]. In recent years, RF has also performed well in grassland classification and identification [12]. For the RF regression, the number of trees was 100, and the number of features was 54. The stop criteria to stop splitting were as follows: (1) the minimum number of samples in a node is 1.0 and (2) the minimum impurity in a node is 0. Then, the selected features and training sample raster files with these parameters were used to parameterize the RF regression model, which was used to retrieve the coverage of NPS.

2.6. Evaluation Indexes of Spectral Enhancement Methods

The Pearson correlation coefficient (Pearson’s r), spectral angle mapper (SAM), erreur relative globale adimensionnelle de synthèse (ERGAS), and universal image quality index (UIQI) were adopted to assess the accuracy of the spectral enhancement methods. Pearson’s r between the same band from original HSI and synthetic HSI reflects the degree of consistency between them. The value range of the Pearson’s r was 0 to 1. The closer the value was to 1, the closer the synthesized image band was to the original band. The SAM is related to spectral distortion [51]. A lower SAM value represents better algorithm performance (the ideal value is 0). ERGAS [52] is mainly related to radiometric distortion. The closer the value is to 0, the better the algorithm performs (the ideal value is 0). UIQI [53] comprehensively considers three factors: related loss, brightness distortion and contrast distortion, and the values range from −1 to 1 (ideally 1).

3. Results

3.1. Predicted Spectra of the ISREM Compared with Those of the CRISP and SREM

3.1.1. Predicted Spectral Curve Comparison of the Three Algorithms

The spectra of the four types of ground surfaces (grassland, bare land, water and snow) were extracted from the synthetic images created by the ISREM compared with the CRISP, SREM and original images. The spectral values of each kind of ground object were randomly selected from the same position of the original image and the synthetic images generated by the CRISP, SREM and ISREM. The final spectra of each kind of ground surface were calculated by the average value and organized in the four different images shown in Figure 5. Generally, the spectral curve vibrated considerably at wavelengths of 460 nm to 520 nm and gradually smoothed out at wavelengths of 520 nm to 950 nm. In Figure 5a,c,d, the spectral curve of the ISREM was closer to the original image than those of the CRISP and SREM, especially at wavelengths of 750 nm to 950 nm. In Figure 5b, the spectral curves of the ISREM and SREM performed better than those of the CRISP.

3.1.2. Statistical Comparison of the Three Algorithms

To further compare the performance of the three algorithms, four central bands of NIR (band 107,862 nm), red (band 61, 647 nm), green (band 39, 585 nm) and blue (band 02,434 nm) bands, were extracted, and the Pearson’s r between the central bands and original bands from the ISREM, CRISP and SREM were calculated. As shown in Table 4, the Pearson’s r of the three algorithms were all above 0.9. The Pearson’s r of the ISREM were significantly higher than those of the other two methods in bands 2 (blue) and 107 (NIR) and slightly higher in bands 39 (green) and 61 (red). The average Pearson’s r indicated that the order of the performance of the three methods was ISREM > SREM > CRISP. In terms of the Pearson’s r, the ISREM performed the best in both the four selected bands and in the average value.

As shown in Table 5, the ISREM performed better than the SREM and CRISP when considering SAM, ERGAS and UIQI. The value of ERGAS that were obtained from the ISREM and SREM reached a satisfactory level. The SAM value obtained from the ISREM also performed very well compared with the other two algorithms. Although the UIQI values of the three algorithms were not as high as those in the existing literature [54], the ISREM still performed the best.

3.2. Coverage of NPS Provided by SVM and RF Regression

The HSI produced by ISREM was used as the original data source for the SVM and RF regressions. First, we combined 39 original bands, 12 first-order bands, and 6 VIs to create an input parameter from the HSI. Then, 255 samples were randomly divided into two categories: 178 training samples and 77 verification samples. To evaluate the spectral enhancement ability of the ISREM, the original multispectral data were also used for regression. Finally, the regression results of NPS were calculated by the SVM and RF, respectively.

3.2.1. Result of SVM and RF Regression

According to the result of species coverage percentage by SVM regression, the coverage percentage of NPS was ranged from 0% to 91% and then divided into 9 grades: 0–10, 10–20, 20–30, 30–40, 40–50, 50–60, 60–70, 70–80, and >80, as shown in Figure 6a. The RF regression of NPS was also carried out. The coverage of NPS ranged from 4% to 87%. A spatial distribution coverage map of NPS with the same coverage grades of RF regression is shown in Figure 6b. Generally, the coverage grade maps of NPS estimated by SVM and RF regressions are consistent. As shown in Figure 6, the low coverage areas of NPS were mainly distributed in the gathering places of the county seat, river valleys, and mountains, while high coverage areas are mainly distributed in the pasture around the town and the sunlit slope of some mountains. This is quite consistent with the results of our field investigation.

3.2.2. Accuracy Evaluation

To quantitatively explore the ability of the ISREM to identify grassland in the study area, the other synthetic HSIs (SREM and CRISP) and the original MSI were used as the comparative data for the regression of the coverage of NPS. For the three generated HSIs, the regression methods in Section 2.5.2 were applied to obtain the abundance information. For the original MSI, the four original bands, VIs and ARI2 in Table 3, were used as input parameters for the SVM and RF regressions.

All 77 verification samples, which were independent of the training samples, were reserved for accuracy verification. The coverage value of the sample pixels was used as a reference to estimate the values. For each validation sample, the difference between the predicted value and actual value was calculated. The quantity and proportion with the difference of less than 10% (0–10%) and less than 20% (0–20%) were used to evaluate the recognition ability of the four images to NPS. The larger number in the threshold (0–10%, 0–20%), the better image performed. As shown in Table 6, ISREM made the largest number of both 0–10% and 0–20%, which means ISREM has shown a recognition ability of NPS.

3.2.3. Degradation Map of Research Area

A degradation map (Shown in Figure 7) was produced according to grassland degradation criterion in the literature [55]. Five grassland degradation levels were described as non-degradation (≥75%), mild degradation (55%~75%), moderate degradation (35%~55%), severe degradation (20~35%), and extreme degradation (≦20%) by the coverage of NPS. As shown in Table 7, moderate and severe degradation areas account for more than 60% of the study area. These areas are mainly distributed in the valley and around the town. Non and mild degradation account for only 19.63% and are mainly distributed in the high mountain slopes, which are difficult to reach. This is consistent with the coverage of NPS.

4. Discussion

4.1. Improvements of ISREM

The SREM was first proposed in 2014 and successfully applied in mineral alteration information extraction. Although the SREM can achieve better spectral and spatial accuracies, it relies too much on pure spectra. These pure spectra were used as spectral predictions for predicting the HSI, which were determined by the minimum spectral angle. However, mixed pixels are common in satellite HSI due to their low spatial resolution, especially in the mixed vegetation areas of the TRHR. Thus, in this study, ISREM made three improvements: (1) Clustering of the overlapping image area was carried out, and the overlapping image was clustered into several types. The average spectra were then extracted for the different types of clustered areas. (2) These average spectra were taken as “pure spectra”, and the weighted spectral angle was used for each pixel of the predicted image. (3) The ISREM can be used in nonoverlapping regions of MSI and HSI; therefore, the coverage of hyperspectral data was expanded and given a higher spatial resolution.

As shown in Section 3.1, the proposed ISREM indeed performed better than both the SREM and CRISP in the spectral and statistical comparison using the HJ-1A MSI and HSI. This might have been caused by two reasons. (1) The difficulty of obtaining pure spectra makes it difficult for the SREM to achieve a better spectral enhancement. (2) To some extent, the weighted spectral angle of “pure spectra” might improve the stability of the ISREM and reduce the effect of mixed spectra on spectral enhancement. Despite the good performance of the ISREM, it focused on mapping grassland degradation in the TRHR, and the data were restricted to HJ-1A MSI and HSI. The ISREM might have the potential for other image data and applications, and this will be explored in the future.

4.2. Potential of Spectral Enhancement in Monitoring Grassland Degradation of the TRHR

Compared to the original MSI in Section 4.2, the predicted HSI had a better performance in mapping NPS in a typical area of the TRHR. HSI can distinguish ground features with similar spectra because of their high spectral resolutions. This also proves that HSI have more potential in mapping grassland degradation in the TRHR [7].

However, limited by the design characteristics of hyperspectral remote sensing sensors, HSI usually has disadvantages, such as low spatial resolution and small spatial coverage. Previous studies have focused on the enhancement of the spectral resolution and have achieved good results [56], while there is very little effort on expanding the spatial coverage of HSI. In this study, the coverage of the HSI was improved by enhancing the spectral resolution of the MSI. Therefore, theoretically, as long as there are small overlaps between the HSI and MSI in the study area, the MSI can be simulated as HSI. This is of great significance for the research areas where remote sensing images, especially HSIs, are scarce.

Although global satellite image data are rapidly growing, relevant image data for the TRHR area, especially effective image data, are still scarce. When collecting remote sensing images to study the degradation of grassland components in a typical area of the TRHR, we found that there were few available HSIs of the study area. Taking HJ-1A imagery data as an example, the only available HSI in the study area was captured in 2013, which was six years before our latest sampling date, making it difficult to guarantee the timeliness of the data. There were two reasons for this: (1) affected by the plateau climate, the grassland growing season in the TRHR is mainly from May to September each year. This short “time window” of less than five months resulted in a significant reduction in the opportunity to obtain useful images; (2) during the growing season, there are many more convective weather and high cloud cover events in the TRHR area, leading most of the images to be unusable. The scarcity of data brings great challenges to grassland degradation research in the TRHR area.

Therefore, the ISREM proposed in this paper can expand the spatial coverage of HSI, enhance the spectral resolution of MSI, and provide a method for obtaining useful data sources. This approach also provided a technical reference for related applications for which it is difficult to directly obtain pure pixels. In addition, by comparing the synthetic HSI with the original MSI, the application potential of HSI in NWs identification was further verified.

5. Conclusions

In this paper, an improved spectral resolution enhancement method for HSI was presented, with the aim of generating a wider HSI located in a typical area of the TRHR. A clustering and weighted spectral angle method was first used as a substitute for pure pixels in generating the transformation matrix. Then, the application of the synthetic HSI to identify NPS in typical areas of the TRHR was carried out using two widely used regression algorithms (RF & SVM). Finally, the coverage maps of NPS were produced. The experiments demonstrated that the clustering and weighted spectral angle provided a better and more stable transformation matrix for predicting the HSI and that the ISREM outperformed two existing spectral resolution enhancement methods in both spectral similarity and correlation with the original reference HSI. Additionally, due to the lack of dependence on pure pixels, the ISREM has more advantages in application scenarios where it is difficult to obtain pure pixels. This alternative clustering and weight spectral angle method may make it universally applicable. The experiments on NPS recognition further proved that HSI generated by ISREM are more advantageous than MSI for identifying NPS. The original purpose of the ISREM was to solve the difficulty in obtaining time- and location-appropriate hyperspectral remote sensing images of the TRHR. Although the method has been shown to work in grassland identification, its application in other scenarios still needs to be further verified. More broadly, the generation of the transformation matrix in the ISREM could be further discussed.

Author Contributions

Conceptualization, B.W. and R.A.; methodology, B.W.; software, B.W. and T.J.; validation, F.X. and F.J.; formal analysis B.W.; investigation, B.W., T.J. and F.X.; resources, F.X.; data curation, T.J.; writing—original draft preparation, B.W.; writing—review and editing, R.A.; visualization, T.J.; supervision, F.J.; project administration, R.A.; funding acquisition, R.A. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by the National Nature Science Foundation of China (No. 41871326; 41271361); Jiangsu Province Key R&D Plan (No. BE2017115).

Acknowledgments

We express our heartfelt gratitude for Zetian Ai, Lijun Huang, Shengyin Zhao, Meiru Zhu, Yina Hu and Xiaoling Zhou for their work of field samples collection. We would like to thank professor Yuehong Chen for his English language modification for this article. We are thankful to anonymous reviewers for their constructive and thoughtful inputs in the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Location of the study area in the TRHR of the Tibetan Plateau. (a) Location of the TRHR of China with the imagery boundary boxes. (b) False-color image clipped by the sample distribution boundary and synthesized by HJ-1A MSI (band 4, band 3 and band 2).

Figure 2. Sample collection work: (a) Spectrum collection of Stipa capillata; (b) a schematic of “X” sampling method; (c) our team recording the quadrat information; (d) main information attribute table structure of the sample collection at this quadrat.

Figure 3. Workflow of ISREM in the study area.

Figure 4. Spectral characteristic curves: (a,b) are the spectral curves of the seven original grass species and the Mahalanobis distance of spectra curves between NPS and NWs. (c,d) are the first-order derivative spectral curves of the seven original grass species and the Mahalanobis distance of first-order derivative spectral curves between NPS and NWs.

Figure 5. Spectral comparison between the original image and synthetic image created by CRISP, SREM, and ISREM, respectively. (a) represents the spectra of grassland, (b)represents the spectra of bare land; (c) represent the spectra of water and, and (d) represents the spectra of snow.

Figure 6. Grassland coverage maps produced by two regressions using ISREM HSI. (a) is NPS abundance level map provide by SVM regression, (a)(i) and (a)(ii) are locally enlarged views of Qumalai and Zhiduo county; (b) is the NPS abundance level map provided by RF regression, (b)(i) and (b)(ii) are the locally enlarged views of Qumalai and Zhiduo county.

Figure 7. Degradation map of a research area based on the coverage map of NPS generated by SVM regression using the ISREM HSI.

Table 1

Comparison between HJ-1A MSI and HSI.

HJ-1A	Bands (nm)		Spatial Resolution (m)
MSI	Band1 430–520	Band3 630–690	30
	Band2 520–600	Band4 760–900
HSI	Band1-Band26	Band61-Band88	100
	460.04–518.81	631.805–759.39
	Band27-Band60	Band89-Band115
	521.475–627.895	765.11–951.54

Table 2

The features selected by Mahalanobis distance analysis from the original spectra and first-order derivative spectra. Feature bands of the fused HJ-1A/HIS were selected based on the corresponding spectral feature scale.

Feature Type	Spectrum (nm)	Bands of HJ-1A/I	Number of Bands
Original spectra	340–523645–693	band 1 to 28band 65 to 75	39
First-order derivative spectra	689–755	band 75 to 87	13

Table 3

VIs selected for regression. In the table, Rx represents the reflectance at wavelength x nm.

VI	Full Name	Formula	Reference
NDVI	Normalized difference vegetation index	(R900 − R680)/((R900 + R680)	[36]
nGNDVI	Narrowband green normalized difference vegetation index	(R780 − R550)/((R780 + R550)	[37]
nNDVI	Narrowband normalized difference vegetation index	(R800 − R670)/((R800 + R670)	[36]
nLCI	Narrowband leaf chlorophyll index	(R850 − R710)/((R850 + R680)	[38]
ARI2	Anthocyanin reflectance index 2	R800 * [(R550)−1 − (R700)−1]	[39]
nPRI	Narrowband photochemical reflectance index	(R550 − R530)/(R550 + R530)	[40]

Table 4

Pearson’s r between selected bands and original bands using ISREM, SREM, and CRISP.

Bands	Pearson’s r of Three Algorithms
	ISREM	SREM	CRISP
B02	0.9300	0.9229	0.9031
B39	0.9539	0.9511	0.9512
B61	0.9625	0.9635	0.9649
B107	0.9716	0.9410	0.9534
Average Value	0.9582	0.9514	0.9480

Table 5

Average cumulative quality indexes of ISREM, CRISP, and SREM.

	ISREM	SREM	CRISP
SAM	8.5591	17.9642	21.8500
ERGAS	0.0481	0.0492	0.0580
UIQI	0.7901	0.6628	0.6771

Table 6

The statistic results of four images using SVM and RF regression (the number of all validation samples is 77. A count of 0–10% and a count of 0–20% represent the number of samples whose difference between the sample value and the estimated value is less than 10% and 20%. Proportion is the ratio of the count value (count of 0–10%, count of 0–20%) to the number of all validation samples.

Type	SVM Regression				RF Regression
	MSI	SREM	CRISP	ISREM	MSI	SREM	CRISP	ISREM
Count of 0–10%	53	52	53	54	49	49	50	52
Proportion (%)	68.83	67.53	68.83	70.13	63.64	63.64	64.94	67.53
Count of 0–20%	69	69	70	72	69	70	70	71
Proportion (%)	89.61%	89.61%	90.91%	93.51%	89.61%	90.91%	90.91%	92.21%

Table 7

Statistics of five grassland degradation levels based on the coverage of NPS.

Degradation Level	Non	Mild	Moderate	Severe	Extreme	Mask
Area (numbers of pixels)	78,392	1,033,373	1,966,572	1,434,292	979,152	171,263
Proportion of grassland (%)	1.43	18.81	35.81	26.12	17.83	/
Proportion of the study area (%)	1.38	18.25	34.73	25.33	17.29	3.02