Comparison of Methods for Estimating Fractional Cover of Photosynthetic and Non-Photosynthetic Vegetation in the Otindag Sandy Land Using GF-1 Wide-Field View Data

Photosynthetic vegetation (PV) and non-photosynthetic vegetation (NPV) are important ground cover types for desertification monitoring and land management. Hyperspectral remote sensing has been proven effective for separating NPV from bare soil, but few studies determined fractional cover of PV (f pv) and NPV (f npv) using multispectral information. The purpose of this study is to evaluate several spectral unmixing approaches for retrieval of f pv and f npv in the Otindag Sandy Land using GF-1 wide-field view (WFV) data. To deal with endmember variability, pixel-invariant (Spectral Mixture Analysis, SMA) and pixel-variable (Multi-Endmember Spectral Mixture Analysis, MESMA, and Automated Monte Carlo Unmixing Analysis, AutoMCU) endmember selection approaches were applied. Observed fractional cover data from 104 field sites were used for comparison. For f pv, all methods show statistically significant correlations with observed data, among which AutoMCU had the highest performance (R2 = 0.49, RMSE = 0.17), followed by MESMA (R2 = 0.48, RMSE = 0.21), and SMA (R2 = 0.47, RMSE = 0.27). For f npv, MESMA had the lowest performance (R2 = 0.11, RMSE = 0.24) because of coupling effects of the NPV and bare soil endmembers, SMA overestimates f npv (R2 = 0.41, RMSE = 0.20), but is significantly correlated with observed data, and AutoMCU provides the most accurate predictions of f npv (R2 = 0.49, RMSE = 0.09). Thus, the AutoMCU approach is proven to be more effective than SMA and MESMA, and GF-1 WFV data are capable of distinguishing NPV from bare soil in the Otindag Sandy Land.


Introduction
Accurate and timely information about vegetation cover is fundamental for monitoring desertification [1][2][3] and assessing the impacts of land management strategies [4] in the Otindag Sandy Land, one of the four largest sandy lands in China.This area has experienced severe desertification before 2000 [5,6], and was the focus of massive ecological engineering measures since the launch of the "Beijing and Tianjin Sandstorm Source Controlling Program" in 2001 [7].Remote sensing offers a unique opportunity to estimate vegetation cover at large scales that might otherwise be costly, labor-intensive, and spatially discontinuous [8][9][10].Over the past several decades, common multispectral vegetation indices, which exploit the difference between visible and near-infrared (NIR) reflectance caused by the presence of chlorophyll, such as the normalized difference vegetation index (NDVI) [11] and the enhanced vegetation index (EVI) [12], have been widely utilized for assessing vegetation cover and dynamics.However, these indices are only sensitive to the amount of photosynthetic vegetation (PV), as well as its turgidity and greenness [13].
From a functional perspective, vegetation can be categorized as photosynthetic (green leaves) and non-photosynthetic (wood, senescent material, and litter) [14].Undoubtedly, photosynthetic vegetation is a critical component of vegetation, but it is not the only component.Non-photosynthetic vegetation (NPV) also plays a key role in carbon and nutrient uptake, fire risk and frequency, and wind and water erosion [15].In the Otindag Sandy Land, non-photosynthetic vegetation is common due to the area's arid climate, relatively sparse vegetation coverage, and frequent droughts.Thus, simultaneously acquiring the fractional cover of PV and NPV in the Otindag Sandy Land would provide new insights for desertification monitoring and land management.
Unlike PV, it is difficult to discriminate NPV from soils due to the similarity in reflectance between many soils and NPVs, particularly when multispectral sensors are considered [16].The cellulose absorption index (CAI), based on hyperspectral data, has been proven to be an effective method to resolve NPV cover [17][18][19].An alternative approach is linear unmixing of hyperspectral bands affected by cellulose and lignin to retrieve NPV cover [3] and map forest fuel [20], and to monitor degradation and desertification [21].However, hyperspectral images are difficult to acquire and cannot meet the demand of NPV cover retrieval at large scales.To date, there are few methods for retrieval of NPV dynamics from multispectral imagery.In the past, numerous studies focused on developing a new multispectral spectral index sensitive to NPV such as the normalized difference senescent vegetation index (NDSVI) [22], the soil adjusted total vegetation index (SATVI) [23], a ratio of moderate resolution imaging spectrometer (MODIS) bands 7 and 6 [14], and the dead fuel index (DFI) [24].Although NPV estimation results using these spectral indices in different environments showed good agreement with field data, this approach is area-specific, and not well validated in other environments.
Spectral mixture analysis (SMA) provides another promising method for retrieving NPV cover from multispectral imagery [15,25].A crucial factor affecting the performance of SMA for NPV cover retrieval is the knowledge of the soil background spectrum for each pixel.In order to resolve endmember spectral variability, multiple endmember SMA (MESMA) and relative spectral mixture analysis (RSMA) have been utilized for PV and NPV cover estimation in an agricultural area [26].However, the choice of the appropriate SMA techniques and endmember spectra may still cause uncertainties, even for simple vegetation structure.The selection of appropriate endmembers, which are usually selected either from the image data [21] or from spectral libraries acquired from field investigations [25], is also critical to the accurate estimation of fractional cover.Each approach has distinct merits and drawbacks.Image-based endmembers are ideal for their consistency in data to be analyzed.However, determination of image endmembers often requires the existence of pure pixels, which are rarely available at 10-30 m spatial scale, especially in arid and semi-arid regions.Field-measured endmembers can be readily collected and their quality can be easily controlled.However, there are potential problems related to the generality and scalability of endmembers.In addition, most SMA techniques for PV and NPV cover estimation were applied to multispectral sensors with short-wave infrared (SWIR) bands (e.g., Landsat TM, and MODIS), and there is limited literature reporting PV and NPV cover estimation using visible and NIR bands only.The GF-1 wide field view (WFV) cameras, which only include four bands in the visible and NIR regions, are highly valuable data sources for dynamic monitoring of vegetation cover at large scales, due to their high spatial resolution (16 m), wide coverage (800 km), and high revisit frequency (2-3 days).
Therefore, the objective of this study is to examine, for the first time, the performance of different SMA techniques for simultaneously estimating the fractional cover of PV and NPV with GF-1 WFV data (only visible and NIR bands), and to compare the results to observed fractional cover measurements in the Otindag Sandy Land.

Study Area
The Otindag Sandy Land is located in central-eastern Inner Mongolia, north of Beijing (Figure 1).The area is about 360 km in length from east to west and 30-100 km in width, with an elevation between 1100 m and 1400 m above sea level, decreasing from southeast to northwest.It is characterized by a temperate continental semi-arid climate with strong wind and less precipitation in winter and spring.Annual average precipitation, mainly occurring in summer and fall with some interannual fluctuation, is 350-400 mm in the southeastern part and below 200 mm in the northwestern part.The major soil type is sandy soil formed by Aeolian processes.The vegetation is characterized by open forest steppe and meadow prairie in the east, semi-arid grassland in the central part, and desert prairie in the west (Figure 2).Because of desertification, most grasslands have experienced different degrees of shrub encroachment.Therefore, the vegetation structure in the Otindag Sandy Land is highly heterogeneous.Due to the scarcity and intra-annual variation of precipitation, NPV accounts for a high percentage of the total vegetation and varies both seasonally and inter-annually.

Study Area
The Otindag Sandy Land is located in central-eastern Inner Mongolia, north of Beijing (Figure 1).The area is about 360 km in length from east to west and 30-100 km in width, with an elevation between 1100 m and 1400 m above sea level, decreasing from southeast to northwest.It is characterized by a temperate continental semi-arid climate with strong wind and less precipitation in winter and spring.Annual average precipitation, mainly occurring in summer and fall with some interannual fluctuation, is 350-400 mm in the southeastern part and below 200 mm in the northwestern part.The major soil type is sandy soil formed by Aeolian processes.The vegetation is characterized by open forest steppe and meadow prairie in the east, semi-arid grassland in the central part, and desert prairie in the west (Figure 2).Because of desertification, most grasslands have experienced different degrees of shrub encroachment.Therefore, the vegetation structure in the Otindag Sandy Land is highly heterogeneous.Due to the scarcity and intra-annual variation of precipitation, NPV accounts for a high percentage of the total vegetation and varies both seasonally and inter-annually.

Remote Sensing Data
The GF-1 satellite operates in a sun-synchronous orbit at an altitude of 645 km and carries four wide-field view (WFV) cameras.In this study, five scenes of GF-1 WFV with four bands (blue, green, red, and NIR), provided by China Center for Resource Satellite Data and Application (CRESDA, http://www.cresda.com/),were utilized to cover the entire Otindag Sandy Land (Table 1).
GF-1 WFV preprocessing is crucial for any quantitative application.First, the digital number (DN) value was converted to radiance using calibration coefficients obtained from the CRESDA.Then, radiance was transformed to surface reflectance using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) algorithm provided with the Environment for Visualizing Images Software (ENVI 5.0).Finally, eight Landsat-8 Operational Land Imager (OLI) datasets provided by the United States Geological Survey (USGS), which have been shown to be geometrically consistent with the field GPS values, were selected as the base map for geometric correction.The geometric correction was conducted in ENVI 5.0, and the resultant geometric coregistration error is less than one pixel of the Landsat-8 OLI panchromatic data (15 m).

Remote Sensing Data
The GF-1 satellite operates in a sun-synchronous orbit at an altitude of 645 km and carries four wide-field view (WFV) cameras.In this study, five scenes of GF-1 WFV with four bands (blue, green, red, and NIR), provided by China Center for Resource Satellite Data and Application (CRESDA, http://www.cresda.com/),were utilized to cover the entire Otindag Sandy Land (Table 1).
GF-1 WFV preprocessing is crucial for any quantitative application.First, the digital number (DN) value was converted to radiance using calibration coefficients obtained from the CRESDA.Then, radiance was transformed to surface reflectance using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) algorithm provided with the Environment for Visualizing Images Software (ENVI 5.0).Finally, eight Landsat-8 Operational Land Imager (OLI) datasets provided by the United States Geological Survey (USGS), which have been shown to be geometrically consistent with the field GPS values, were selected as the base map for geometric correction.The geometric correction was conducted in ENVI 5.0, and the resultant geometric co-registration error is less than one pixel of the Landsat-8 OLI panchromatic data (15 m).

Field Spectroscopy
Recognizing the multitude of highly variable PV, NPV, and bare soils in the Otindag Sandy Land, their spectral properties were thoroughly investigated in order to acquire the most generalizable endmember dataset (Figure 1).Spectra were collected with an Analytical Spectral Devices (ASD) full-range (350-2500 nm) Fieldspec ® 4 spectroradiometer with a 25 • sensor foreoptic.All measurements were collected within two hours of local solar noon on clear sky days.The sensor was held 1 m above the top of the PV and NPV canopy or bare soil surface in a vertical downward position.Prior to each measurement, the spectroradiometer was calibrated to a Spectralon ® white reference target.
Field spectra were collected during two periods.The first measurement was conducted on 27 and 31 July 2014, representative of maximum PV existence, quasi-synchronous with the GF-1 data acquisition time.Twenty-nine PV spectra, 14 NPV spectra, and 12 bare soil spectra were recorded across the whole region.In order to acquire additional NPV spectra, a second measurement was conducted on 4 and 5 November 2014, representative of maximum NPV existence.Fourteen NPV spectra and three bare soil spectra were acquired.Finally, an endmember library, consisting of 29 PV, 28 NPV, and 15 bare soil spectra, was established.Based on the spectral response function of the GF-1 WFV sensor, the field spectra were resampled to the GF-1 WFV bands (Equation ( 1)), where R i is the simulated GF-1 WFV reflectance, f i (λ) is the spectral response function at λ wavelength of GF-1 WFV band ith, r (λ) is the field observed reflectance at λ wavelength, and λ min and λ max represents the band width of GF-1 WFV band ith.

.3. Fractional Ground Cover Data
Fractional ground cover data were collected in late July and early August, the season of maximum green vegetation cover in 2014.Based on a stratified random design and accessibility, 104 sites were selected (Figure 1), and one 32 m × 32 m plot was set up in each site.For sampling, two transects coinciding with the diagonal lines across the plot were applied, where vegetation cover in the three categories (non-woody ground cover, woody less than 2 m high, and woody greater than 2 m in height) is recorded.Using a tape laid on the ground, the surveyor recorded the presence of PV, NPV, and bare soil at 1-m intervals using a laser pointer.Non-woody vegetation was observed by looking downward, while woody vegetation was observed by looking upward or downwards where vegetation was less than the observer's height.Total cover fractions per category were calculated by dividing the number of counts of a particular cover type by the total number of counts (90).For each plot, an exposed planar fractional cover was calculated by merging all presented categories, based on the approach proposed by Muir et al. (2011) [27].The coordinates of the cross point of the two transects were recorded by a global position system device (Trimble GeoXT 2800-3000) with a positioning accuracy of approximately ±3 m to match the GF-1 WFV data.

Spectral Mixture Analysis
In SMA, the reflectance of a pixel is assumed to be a linear combination of the reflectance of the spectra of the endmembers, weighted by their fractional cover.Endmembers are fundamental physical components that themselves are not mixtures of other components.SMA assumes that the endmembers are spatially and temporally invariant.In this study, the average spectra of PV, NPV, and bare soil were utilized as the endmember spectra.

Multiple Endmember SMA
Multiple Endmember SMA was developed to handle the endmember variability problem through the utilization of iterative mixture analysis cycles [28].In MESMA, endmembers are allowed to vary on a per-pixel basis.When all endmember combinations are calculated, the best-fit model is determined for each pixel.One criterion most often used in MESMA is the lowest Root Mean Square Error (RMSE) of spectral fit [28].In this study, models were calculated by utilizing all possible combinations among 29 PV spectra, 27 NPV spectra, and 15 bare soil spectra resulting in 11,745 total models.

AutoMCU
Another approach dealing with endmember variability, the Automatically Monte Carlo spectral Unmixing model (AutoMCU), was proposed by Asner and Lobell [29], in which, a large number of endmember combinations for each pixel is calculated by randomly selecting spectral data from a spectral library.Studies [25,29] have shown that the fractional cover of different endmembers is distributed normally when the number of endmember combinations is sufficient.Thus, the average value of the fractional cover for different models was used as the final fractional cover of each pixel.We set up an initial value of 300 iterations, followed by a determination of the appropriate times through the analysis of the average fractional cover dynamic.If the average fractional cover does not vary with the increasing unmixing iterations, this time would be considered as the appropriate time.

Unmixing Strategy
The three approaches are different with regard to endmember selection.The Fully Constrained Least Square (FCLS) algorithm was applied for calculating f k at each time [30].Using FCLS, two important constraints on f k (the fraction sum-to-one constraint (ASC) ∑ n k=1 f k = 1 and the fraction nonnegativity constraint (ANC) f k ≥ 0) could be resolved simultaneously.
According to the linear unmixing model, the spectral signature of a pixel vector r can be represented by a linear regression model of ρ as follows: In order to account for ASC, we included ASC in the signature matrix ρ by introducing a new signature matrix, denoted by M and defined as: with 1 = (1, 1, . . . , 1)T n , and a vector R = δr 1 , where δ is a weight for deciding the degree of ASC; the smaller the δ value, the closer the sum to one.A value of 1 was applied in this study.Combining Equations ( 2) and (3), the following equation is obtained: The nonnegatively constrained least squares (NCLS) imposes the nonnegativity constraint (ANC) on the abundance vector while estimating fractional cover.The estimation of abundances subject to nonnegative constraints is difficult to implement since it results in a set of inequalities and can only be solved by numerical methods.The iteration algorithm proposed by Chang and Heinz [31] was adopted by introducing a Lagrange multiplier vector (λ).The algorithms for FCLS, SMA, MESMS, and AutoMCU were transformed into IDL (Interactive Data Language) routines.

Comparison with Observed Data
To compare the performance of different SMA techniques on PV/NPV fractional cover estimation, two metrics were calculated against observed data, the RMSE and coefficient of determination (R 2 ) of linear regression: where n is the number of fields, x i is the estimated fractional cover of field i, y i is the measured fractional cover of field i, x is the average value of the estimated fractional cover, and y is the average value of the measured fractional cover.

Field f pv f npv , and f soil Measurements
Totally, 104 sample plots were investigated, 12 in open woodland, 40 in grassland encroached by shrub and 52 in grassland, respectively.Field observed f pv , f npv , and f soil followed the expected temporal and structural patterns (Table 2).Because the acquiring time corresponds to the maximum vegetation growing season, fields were dominated by green vegetation resulting in high f pv .Even so, NPV takes an important proportion of sample plots, especially in open woodland (average f npv reaches 24%).Compared to grassland, the open woodland was characterized by higher f npv and lower f pv , while no obvious difference was observed for f soil .

Endmember Library and Variability
The endmember library for the Otindag Sandy Land, divided into PV, NPV, and bare soil categories and consistent with GF-1 WFV data, was established (Figure 3).It is obvious that the PV spectra are easily distinguishable from NPV and bare soil spectra.However, the NPV and bare soil spectra are similar, both exhibiting a monotonous increasing tendency, which makes separation more difficult.Even so, there are some differences between NPV and bare soil spectra.First, bare soil spectra show a significantly higher reflectance than the NPV spectra, mainly due to extensive bright sandy substrate.Second, a bow-shaped protuberance, which was first recognized by Gao et al. [32], exists from the blue to the red bands in the bare soil spectra, while it is absent in the NPV spectra.These spectral characteristics of NPV and bare soil provide new opportunities for the discrimination of NPV from bare soil in the Otindag Sandy Land.With regard to intra-variability of the three endmember categories, the PV and NPV spectra are relatively concentrated, while the bare soil spectra vary greatly.Therefore, the bare soil endmember selection will have a greater influence on fractional cover estimation, especially on the fractional cover of NPV.In any case, the choice of optimal endmember combinations is crucial for the spectral unmixing process.While MESMA will try all combinations of PV, NPV, and bare soil endmembers, AutoMCU will utilize as many combinations as necessary through random selection.For SMA, the average spectra of three endmember categories were employed as invariant PV, NPV, and bare soil spectra for spectral unmixing.

SMA
Choosing the average PV, NPV, and bare soil spectra as endmembers, fpv, fnpv, and fsoil in the Otindag Sandy Land were estimated using the FCLS algorithm.The RMSE of SMA for the samples ranges from 1.24% to 19.91%, with an average RMSE value of 3.76%.The spatial distribution of fpv and fnpv is shown in Figure 4a.Overall, the spatial distribution characteristics are consistent with prior field investigations.From east to west, the fractional cover of both PV and NPV shows an obvious decreasing tendency, coincident with decreasing precipitation.The high fpv, mainly located in the eastern part of the Otindag Sandy Land, corresponds to forest, farmland, and wetland.The high fnpv, was mainly found in the grassland along the southern and northern edge of the eastern Otindag Sandy Land, and the high fsoil mainly correspond to desert steppe and active sand dunes, respectively.
Based on the 104 field samples, the validation results of SMA-based fpv and fnpv are shown in Figure 5, and their statistical data are summarized in Table 3.It can be seen that the SMA-based fpv and fnpv are highly correlated with observed data, with R 2 values of 0.47 and 0.41, respectively.However, there is an obvious underestimation of fpv (RMSE = 0.27), and an obvious overestimation of fnpv (RMSE = 0.20).The high correlation demonstrates that PV and NPV can be distinguished by SMA using GF-1 WFV data, but the absolute error suggests that using the average spectra of the PV, NPV, and bare soil endmember library as the invariant endmember is problematic.

SMA
Choosing the average PV, NPV, and bare soil spectra as endmembers, f pv , f npv , and f soil in the Otindag Sandy Land were estimated using the FCLS algorithm.The RMSE of SMA for the samples ranges from 1.24% to 19.91%, with an average RMSE value of 3.76%.The spatial distribution of f pv and f npv is shown in Figure 4a.Overall, the spatial distribution characteristics are consistent with prior field investigations.From east to west, the fractional cover of both PV and NPV shows an obvious decreasing tendency, coincident with decreasing precipitation.The high f pv , mainly located in the eastern part of the Otindag Sandy Land, corresponds to forest, farmland, and wetland.The high f npv , was mainly found in the grassland along the southern and northern edge of the eastern Otindag Sandy Land, and the high f soil mainly correspond to desert steppe and active sand dunes, respectively.
Based on the 104 field samples, the validation results of SMA-based f pv and f npv are shown in Figure 5, and their statistical data are summarized in Table 3.It can be seen that the SMA-based f pv and f npv are highly correlated with observed data, with R 2 values of 0.47 and 0.41, respectively.However, there is an obvious underestimation of f pv (RMSE = 0.27), and an obvious overestimation of f npv (RMSE = 0.20).The high correlation demonstrates that PV and NPV can be distinguished by SMA using GF-1 WFV data, but the absolute error suggests that using the average spectra of the PV, NPV, and bare soil endmember library as the invariant endmember is problematic.

MESMA
The spatial distribution of fpv, fnpv and fsoil is shown in Figure 4b.The overall spatial pattern is consistent with the SMA-based results, but with obvious difference of fpv, fnpv in the east Otindag Sandy Land.The minimal RMSE was selected as the criterion to choose the appropriate endmember from the endmember library.The RMSEs of MESMA for the sample plots range from 0% to 8.46%, with an average RMSE value of 0.53%, significantly lower than that of the SMA (3.76%), indicating that MESMA fits the GF-1 WFV data better than SMA.
In Figure 4, we can see that MESMA-based fpv correlate with observed fpv well (with R 2 of 0.48), but with the obvious underestimation of fpv (RMSE = 0.21).However, the accuracy of MESMA-based fnpv and fsoil are poor, with a weak relationship with observed data (with R 2 of 0.11 and 0.15, respectively).Therefore, MESMA performs better for fpv, because the PV endmember is clearly different from the NPV and bare soil endmembers, but would lead to confusion with regard to fnpv and fsoil because of the similarities of the NPV and bare soil endmember spectra.

MESMA
The spatial distribution of f pv , f npv and f soil is shown in Figure 4b.The overall spatial pattern is consistent with the SMA-based results, but with obvious difference of f pv , f npv in the east Otindag Sandy Land.The minimal RMSE was selected as the criterion to choose the appropriate endmember from the endmember library.The RMSEs of MESMA for the sample plots range from 0% to 8.46%, with an average RMSE value of 0.53%, significantly lower than that of the SMA (3.76%), indicating that MESMA fits the GF-1 WFV data better than SMA.
In Figure 4, we can see that MESMA-based f pv correlate with observed f pv well (with R 2 of 0.48), but with the obvious underestimation of f pv (RMSE = 0.21).However, the accuracy of MESMA-based f npv and f soil are poor, with a weak relationship with observed data (with R 2 of 0.11 and 0.15, respectively).Therefore, MESMA performs better for f pv , because the PV endmember is clearly different from the NPV and bare soil endmembers, but would lead to confusion with regard to f npv and f soil because of the similarities of the NPV and bare soil endmember spectra.

AutoMCU
The AutoMCU was firstly utilized for the GF-1 WFV pixels corresponding to the 104 field samples with 300 iterations.The f pv and f npv distribution and their average values along unmixing times obey the same change rules, one of which is illustrated in Figure 6.The f pv and f npv values are highly variable during the first 50 iterations, indicating that endmember selection evidently influences f pv and f npv estimation results.When the number of iterations exceeds 150, the mean f pv and f npv values are essentially stable.In addition, the f pv and f npv distributions follow Gaussian distribution, indicating that their expected (average) values represent the true values.Therefore, 150 iterations were used for unmixing of the entire GF-1 WFV dataset by AutoMCU.
The validation results of AutoMCU-based f pv , f npv , and f soil show that AutoMCU performs well for f pv and f npv .The R 2 of the regression of f pv and f npv on observed data is 0.49, indicating a significant correlation with observed data.The RMSEs for the f pv and f npv estimation are 0.17 and 0.09, respectively.The AutoMCU was firstly utilized for the GF-1 WFV pixels corresponding to the 104 field samples with 300 iterations.The fpv and fnpv distribution and their average values along unmixing times obey the same change rules, one of which is illustrated in Figure 6.The fpv and fnpv values are highly variable during the first 50 iterations, indicating that endmember selection evidently influences fpv and fnpv estimation results.When the number of iterations exceeds 150, the mean fpv and fnpv values are essentially stable.In addition, the fpv and fnpv distributions follow Gaussian distribution, indicating that their expected (average) values represent the true values.Therefore, 150 iterations were used for unmixing of the entire GF-1 WFV dataset by AutoMCU.
The validation results of AutoMCU-based fpv, fnpv, and fsoil show that AutoMCU performs well for fpv and fnpv.The R 2 of the regression of fpv and fnpv on observed data is 0.49, indicating a significant correlation with observed data.The RMSEs for the fpv and fnpv estimation are 0.17 and 0.09, respectively.

Comparisons of Different Methods
The computation efficiency of the three approaches varies widely, depending on the total unmixing times for each pixel.SMA is the fastest (only one time), followed by AutoMCU (150 times), and MESMA is the slowest since all endmember combinations would be used (11,745 times).
The accuracy of these approaches is also obviously different (Table 3).The accuracy of MESMAbased fpv was improved compared to SMA (increase in R 2 from 0.47 to 0.48, and decrease in RMSE from 0.27 to 0.21), and the underestimation of fpv in SMA was alleviated slightly.However, the

Comparisons of Different Methods
The computation efficiency of the three approaches varies widely, depending on the total unmixing times for each pixel.SMA is the fastest (only one time), followed by AutoMCU (150 times), and MESMA is the slowest since all endmember combinations would be used (11,745 times).
The accuracy of these approaches is also obviously different (Table 3).The accuracy of MESMA-based f pv was improved compared to SMA (increase in R 2 from 0.47 to 0.48, and decrease in RMSE from 0.27 to 0.21), and the underestimation of f pv in SMA was alleviated slightly.However, the accuracy of MESMA-based f npv deteriorates significantly, with a weak relationship with observed data (R 2 = 0.11) and an increase in the RMSE from 0.22 to 0.24.In addition, the MESMA-based f npv was significantly overestimated, similar to SMA.Therefore, MESMA performs better for estimating f pv , but would lead to serious mistakes for f npv and f soil .Therefore, AutoMCU is an effective approach for f pv and f npv estimation in the Otindag Sandy Land because it performs better than SMA and MESMA.
The accuracies for f pv and f npv estimation are different for different land cover types.For example, in f pv estimation with the AutoMCU approach, grassland has the highest RMSE (0.19), followed by grassland encroached by shrub (RMSE = 0.16) and open woodland (RMSE = 0.14).However, the opposite tendency appears in f npv estimation.Open woodland has the highest RMSE (0.13), followed by grassland encroached by shrub (RMSE = 0.10), and grassland (RMSE = 0.07).Differences in accuracy were consistent with differences in average f pv and f npv values of sample plots in different land cover types (Table 3).

Separability of NPV and Bare Soil in SMA
The separation of PV from others is not difficult, and it is widely used with VIS-NIR bands [11][12][13].Therefore, the method of distinguishing NPV and Bare soil is critical in this study.Unlike other regions, there are two differences between bare soil and NPV spectra in Otindag Sandy Land.First and foremost, the absolute reflectance of bare soil is obviously higher than that of NPV.In addition, a bow-shaped protuberance exists in bare soil when observed from the spectral library, but it is absent in NPV.However, this characteristic is not sufficiently significant because of the coarse spectral resolution.Therefore, the higher reflectance of bare soil relative to NPV due to the extensive sandy substrate is an important prerequisite for separating NPV from bare soil in the Otindag Sandy Land, indicating that this approach cannot be used in regions with a large percentage of dark colored bare soil.
For a given pixel, the key point for successfully estimating f npv and f soil using the GF-1 WFV data is the determination of appropriate endmembers-especially the bare soil endmember, which shows large intra-variability when compared to NPV spectral libraries-and not the absolute difference of reflectance of NPV and bare soil.Thus, identifying the most appropriate NPV and bare soil endmember for unmixing is critical for retrieving ground cover accurately.

Effects of Endmember Variability
In this study, three unmixing approaches for determining the appropriate endmember combinations, SMA, MESMA, and AutoMCU, were compared to estimate f pv and f npv in the Otindag Sandy Land with GF-1 WFV data.Table 2 shows that all SMA-based fractions are highly correlated with observed data, but that f pv is clearly underestimated, while f npv is overestimated; only f soil shows no directional bias.Thus, the use of the average spectra of PV and NPV as invariant endmembers is not appropriate, but identifying the most representative PV and NPV endmember spectra is extremely difficult.
In order to resolve this problem, MESMA was developed to explore "optimal" pixel-varied specific endmember combinations by setting up a selection criterion (e.g., the lowest RMSE), which has been proven to be effective for estimating f pv and f npv using Landsat [33], MODIS [34] and, Geoeye [35] data.However, MESMA was shown to be incapable for f npv estimation with GF-1 WFV data in the Otindag Sandy Land, because it would lead to serious confusion between NPV and bare soil, with R 2 values of 0.11 and 0.15, and RMSEs of 0.24 and 0.21, respectively.The endmembers most frequently utilized in MESMA are shown in Figure 7. Three PV (including mean spectra), seven NPV, and four bare soil endmembers were utilized most frequently, accounting for 70%, 68%, and 93%, respectively.
From the spectral curve, we can see that the bare soil endmember spectra differ more widely than those for PV and NPV.Thus, the selection of the bare soil endmember will have a significant effect on the unmixing results.For example, the use of the bare soil endmember with obviously higher reflectance would lead to a significant underestimation of f soil .Therefore, the appropriate bare soil endmember determination would be essential for estimating f npv and f soil accurately.Since the spectral similarity between NPV and bare soil, the actual NPV and bare soil endmember combinations could be replaced by another combinations in MESMA in order to minimize the RMSE of spectral fit.This finding is consistent with Okin et al. [26], who found that MESMA did not improve f soil estimates, because the similarity between soil and NPV spectra can actually lead to errors in MESMA as some bare soil/NPV combinations may be mistaken as combinations of other bare soil/NPV [36].Therefore, it cannot be assumed that MESMA always leads to better f npv and f soil estimates.
Remote Sens. 2016, 8, x FOR PEER REVIEW 13 of 19 endmember determination would be essential for estimating fnpv and fsoil accurately.Since the spectral similarity between NPV and bare soil, the actual NPV and bare soil endmember combinations could be replaced by another combinations in MESMA in order to minimize the RMSE of spectral fit.This finding is consistent with Okin et al. [26], who found that MESMA did not improve fsoil estimates, because the similarity between soil and NPV spectra can actually lead to errors in MESMA as some bare soil/NPV combinations may be mistaken as combinations of other bare soil/NPV [36].Therefore, it cannot be assumed that MESMA always leads to better fnpv and fsoil estimates.The AutoMCU method, which randomly selects endmembers from the spectral library to unmix the image over several iterations, rather than doing an exhaustive search, has successfully been applied to arid and semi-arid rangeland [25,[37][38][39] for fpv and fnpv estimation.Like in previous studies, the AutoMCU is shown to be effective in semi-arid sandy lands.In addition to improving computational efficiency, AutoMCU is able to explicitly quantify the proportion indeterminacy [40].It is worth noting that the AutoMCU was originally designed for hyperspectral sensors [25], and that only few studies applied it to multispectral sensors, mainly Landsat [37,41], but with a lack of quantitative validation of fpv and fnpv.The distinguishing feature of this study is that we were able to show that GF-1 WFV data, lacking the short-wave infrared bands, are capable of distinguishing NPV from bare soil in the Otindag Sandy Land.
Differences in the spatial distribution of fpv and fnpv are shown in Figure 8.Compared to SMA, AutoMCU-based fpv became higher (up to 10%) in east Otindag Sandy Land, except for cropland and wetland with the highest green vegetation cover, and decreased (up to 16%) in the entire west Otindag Sandy Land, corresponding to desert steppe; AutoMCU-based fnpv decreased for most vegetated areas (up to 20%), resolving the overestimation problem of fnpv in the SMA (Figure 5).However, there are extensive regions with increasing fnpv (up to 10%), which were verified to correspond to regions with a bright background of sandy soil.More importantly, these increases do not lead to more serious overestimation.Thus, fnpv estimation accuracy was improved effectively.Compared to the MESMA, AutoMCU-based fpv and fnpv changed dramatically, because MESMA led to a serious confusion between NPV and bare soil.fpv increased for vegetated areas and decreased for sandy lands, which alleviated the underestimation problem of fpv in MESMA.fnpv decreased significantly for salt marshes in west regions, which was misclassified as NPV in MESMA.The AutoMCU method, which randomly selects endmembers from the spectral library to unmix the image over several iterations, rather than doing an exhaustive search, has successfully been applied to arid and semi-arid rangeland [25,[37][38][39] for f pv and f npv estimation.Like in previous studies, the AutoMCU is shown to be effective in semi-arid sandy lands.In addition to improving computational efficiency, AutoMCU is able to explicitly quantify the proportion indeterminacy [40].It is worth noting that the AutoMCU was originally designed for hyperspectral sensors [25], and that only few studies applied it to multispectral sensors, mainly Landsat [37,41], but with a lack of quantitative validation of f pv and f npv .The distinguishing feature of this study is that we were able to show that GF-1 WFV data, lacking the short-wave infrared bands, are capable of distinguishing NPV from bare soil in the Otindag Sandy Land.
Differences in the spatial distribution of f pv and f npv are shown in Figure 8.Compared to SMA, AutoMCU-based f pv became higher (up to 10%) in east Otindag Sandy Land, except for cropland and wetland with the highest green vegetation cover, and decreased (up to 16%) in the entire west Otindag Sandy Land, corresponding to desert steppe; AutoMCU-based f npv decreased for most vegetated areas (up to 20%), resolving the overestimation problem of f npv in the SMA (Figure 5).However, there are extensive regions with increasing f npv (up to 10%), which were verified to correspond to regions with a bright background of sandy soil.More importantly, these increases do not lead to more serious overestimation.Thus, f npv estimation accuracy was improved effectively.Compared to the MESMA, AutoMCU-based f pv and f npv changed dramatically, because MESMA led to a serious confusion between NPV and bare soil.f pv increased for vegetated areas and decreased for sandy lands, which alleviated the underestimation problem of f pv in MESMA.f npv decreased significantly for salt marshes in west regions, which was misclassified as NPV in MESMA.

Cross-Multispectral Sensor Comparison
Based on GF-1 WFV data, the unmixing results indicate that the f pv estimations are relatively stable, regardless of the method used.Estimated f pv is significantly correlated with field observed f pv (p < 0.0001), with R 2 values ranging from 0.47 to 0.49.In addition, the obvious underestimation problem gradually improved from SMA to MESMA and AutoMCU, which shows that the per-pixel variable endmember combinations will increase the accuracy of the f pv estimation.However, the f npv estimations are highly variable: the SMA-and AutoMCU-based f npv estimates are significantly correlated with field observed f npv , with R 2 values of 0.41 and 0.49, respectively, but, the correlation between MESMA-based f npv and in situ f npv is poor (R 2 = 0.11).
In general, the AutoMCU approach results in an optimal estimation of f pv and f npv .Here, f npv is estimated most accurately, with an RMSE of 0.09 (R 2 = 0.49), while f pv is estimated less accurately, with an RMSE of 0.17 (R 2 = 0.49), and f soil shows the worst estimation results, with an RMSE of 0.20 (R 2 = 0.48).An important aspect of this study is the field verification.Few studies compare f pv and f npv from multispectral remote sensing imagery with field measurements.Table 4 [26,42,43] lists recently published results for f pv and f npv estimation based on different multispectral sensors.Since source data, study region, and study period differ, the comparison to our study is limited.For f pv estimation, this study shows lower accuracy, with the lowest RMSE acquired by AutoMCU of 17%, while the RMSE numbers of previous studies range from 7% to 14.7%.However, this study demonstrates improved f npv estimation: the lowest RMSE acquired by AutoMCU is 9%, compared to the previous studies' 12%-20.5%.An important difference is that the accuracy for f npv estimation is higher than that for f pv estimation, which might partly be caused by the lower field observed f npv coverage (≤45%) during the peak growing season in our study area, while field observed f npv values ranged from 0% to 100% when multi-temporal data were considered.Our main goal was to validate the ability of the four-band GF-1 WFV (not including any SWIR bands) to distinguish PV and NPV from bare soil.Based on f pv and f npv estimation accuracy comparisons, we conclude that both f pv and f npv can be determined with modest accuracy using GF-1 WFV data with the AutoMCU approach.

Uncertainties, and Sources of Error
The field patch was usually larger than 2 × 2 GF-1 WFV pixels with a homogeneous vegetation growth status.Therefore, the plot scale (32 m) may represent the scale of a GF-1 WFV pixel (16 m), and the geometric co-registration error from field measurements to GF-1 WFV data was effectively reduced.
Time differences may occur between field observations and the time of GF-1 WFV data acquisition.This effect was minimized because the time difference did not exceed one week.We identified some particular cases of large PV underestimation for grassland sites when the observation was made before grass cutting, but the GF-1 WFV image was acquired after grass cutting.However, grass cutting in early August is relatively rare, and therefore, it has no significant effects on the accuracy assessment.
The f pv and f npv estimation errors could be partially related to errors in the field measurements.Sampling procedure along the diagonal lines has been widely acknowledged to be effective for ground cover measurement of sample plots, wherein vegetation is not distributed in parallel rows [44,45].In particular, field surveyors were well trained and measurements of each transect for the first 20 sample plots were cross-validated to guarantee reliability of the field reference data.However, there were some errors caused by the observer's bias.It is most difficult to acquire consistent field data for NPV and estimates can vary greatly between observers.This can result in the confusion between PV and NPV.

Conclusions
We investigated the use of different spectral unmixing approaches (SMA, MESMA, and AutoMCU) for simultaneously estimating f pv and f npv based on GF-1 WFV data without any SWIR bands.The unmixing results were evaluated against observed data collected concurrently with image acquisition.The main findings are: (1) Despite spectral similarity of NPV and bare soil, there are some differences in the GF-1 WFV bands.First and foremost, the bare soil spectra are significantly higher than the NPV spectra, mainly due to extensive bright sandy substrate.In addition, a bow-shaped protuberance exists from blue to red bands in the bare soil spectra, which is not present in the NPV spectra.(2) Due to the complex and bright soil background of the Otindag Sandy Land, the bare soil endmember libraries show large intra-variability.Therefore, determining the appropriate endmember combinations, especially the bare soil endmember, is a key process for successfully estimating f pv and f npv .(3) Invariant endmember combinations should be used with caution, because they can lead to serious over-or underestimation problems (SMA).The MESMA cannot be assumed to always perform better than SMA, due to the coupling of the NPV and bare soil endmembers.AutoMCU was shown to be effective for dealing with endmember variability while acquiring accurate f pv and f npv estimation.Compared with SMA, both R 2 and RMSE are improved.(4) Compared to other relevant multispectral applications, the GF-1 WFV data were shown to be capable for f pv and f npv estimation in the Otindag Sandy Land, despite a lack of the important SWIR bands, which are considered important for separation of NPV from bare soil.With GF-1 WFV's unique advantage of high spatial resolution (16 m), wide coverage (800 km), and high revisit frequency (2-3 days), there is great potential for future analyses.

Figure 1 .
Figure 1.Location of the study area in North China.Figure 1. Location of the study area in North China.

Figure 1 .Figure 2 .
Figure 1.Location of the study area in North China.Figure 1. Location of the study area in North China.

Figure 2 .
Figure 2. Photos of major vegetation structure and characteristics in the Otindag Sandy Land: (a) open woodland; (b) grassland encroached by shrub; (c) degraded grassland; and (d) high percentage of non-photosynthetic vegetation (NPV).

Figure 3 .
Figure 3. Endmember spectra of (a) Photosynthetic vegetation (PV); (b) Non-photosynthetic vegetation (NPV) and (c) Bare soil used in Spectral Mixture Analysis (SMA), the filled squares are the average values of different types.

Figure 3 .
Figure 3. Endmember spectra of (a) Photosynthetic vegetation (PV); (b) Non-photosynthetic vegetation (NPV) and (c) Bare soil used in Spectral Mixture Analysis (SMA), the filled squares are the average values of different types.

Figure 5 .
Figure 5. f pv , f npv and f soil estimation accuracy validation against field samples.

Figure 6 .
Figure 6.f pv (a); and f npv (b) distribution and their average values relative to unmixing iterations using Automated Monte Carlo Unmixing Analysis (AutoMCU).

Table 1 .
Details of the GF-1 wide-field view (WFV) data used and its properties.

Table 1 .
Details of the GF-1 wide-field view (WFV) data used and its properties.

Table 2 .
Characteristics of sample plots.Values shown for f pv , f npv , and f soil are average ± standard deviations.

Table 3 .
Regression analysis and validation of remote sensing results against observed data.

Table 3 .
Regression analysis and validation of remote sensing results against observed data.

R 2 RMSE T RMSE OW RMSE GS RMSE G R 2 RMSE T RMSE OW RMSE GS RMSE G R 2 RMSE T RMSE OW RMSE GS RMSE G
stands for bare soil.RMSE T stands for total RMSE, RMSE OW stands for the RMSE of open woodland samples, RMSE GS stands for the RMSE of grassland with shrub samples, RMSE G stands for the RMSE of pure grassland samples.* p < 0.0001. BS

Table 4 .
Comparison of multispectral remote sensing-based f pv and f npv estimates.