Fractional Snow-Cover Mapping Based on MODIS and UAV Data over the Tibetan Plateau

Abstract

Moderate-resolution imaging spectroradiometer (MODIS) snow-cover products have relatively low accuracy over the Tibetan Plateau because of its complex terrain and shallow, fragmented snow cover. In this study, fractional snow-cover (FSC) mapping algorithms were developed using a linear regression model (LR), a linear spectral mixture analysis model (LSMA) and a back-propagation artificial neural network model (BP-ANN) based on MODIS data (version 006) and unmanned aerial vehicle (UAV) data. The accuracies of the three models were validated against Landsat 8 Operational Land Imager (OLI) snow-cover maps (Landsat 8 FSC) and compared with the MODIS global FSC product (MOD10A1 FSC, version 005) for the purpose of finding the optimal algorithm for FSC extraction for the Tibetan Plateau. The results showed that (1) the overall retrieval results of the LR and BP-ANN models based on MODIS and UAV data were relatively similar to the OLI snow-cover maps; the accuracy and stability were greatly improved, with even some reduction in errors; compared to the Landsat 8 FSC, the correlation coefficients (r) were 0.8222 and 0.8445 respectively and the root-mean-square errors (RMSEs) were 0.2304 and 0.2201, respectively. (2) The accuracy and stability of the fully constrained LSMA model using the pixel purity index (PPI) endmember extraction method based only on MODIS data suffered the worst performance of the three models; r was only 0.7921 and the RMSE was as large as 0.3485. There were some serious omission phenomena in the study area, specifically for the largest mean absolute error (MAE = 0.2755) and positive mean error (PME = 0.3411). (3) The accuracy of the MOD10A1 FSC product was much lower than that of the LR and BP-ANN models, although its accuracy slightly better that of the LSMA based on comprehensive evaluation of six accuracy indices. (4) The optimal model was the BP-ANN model with combined inputs of surface reflectivity data (R1–R7), elevation (DEM) and temperature (LST), which can easily incorporate auxiliary information (DEM and LST) on the basis of (R1–R7) during the relationship training period and can effectively improve the accuracy of snow area monitoring—it is the ideal algorithm for retrieving FSC for the Tibetan Plateau.

1. Introduction

Snow cover is among the most active natural factors on land surfaces and its characteristics (e.g., distribution, area and depth) are important input parameters for global energy balance, climate, hydrological and ecological models [1]. The Tibetan Plateau is an important seasonal snow-capped area in China and it is of great significance to accurately estimate the snow in the area. Remote sensing technology is the only means for large-scale real-time snow-cover monitoring of the area [2]. The new-generation moderate-resolution imaging spectroradiometer (MODIS) is among the most widely used sources of remote sensing data [3,4,5]. The accuracy of the most commonly used global MODIS product, the binary MODIS snow-cover map [6], is unsatisfactory because of mixed pixel problems over the Tibetan Plateau [7]. To date, fractional snow-cover (FSC) mapping algorithms have been the optimal means for solving the mixed pixel problem.

Currently, there are three main types of MODIS fractional snow-cover mapping algorithms: regression analysis algorithms [8,9,10,11,12,13,14,15,16,17], spectral unmixing algorithms [18,19,20,21,22] and machine learning algorithms [23,24,25,26,27]. Barton et al. (2001) [8], Kaufman et al. (2002) [10] and Salomonson et al. (2004) [13] successively established regression models employing FSCs extracted from Landsat data as true values and normalized difference snow index (NDSI) data derived from MODIS data, which they used to estimate FSC at a subpixel scale. Of those regression models, the univariate regression model proposed by Salomonson et al. (2004) [14] was used by the National Snow and Ice Data Center (NSIDC, US) to produce version 005 of the MODIS global FSC product (denoted as MOD10A1 FSC/MYD10A1 FSC). Employing a validation analysis, Zhang et al. (2014) found that MOD10A1 FSC had a low accuracy in the Tibetan Plateau region (with a correlation coefficient of only 0.7400) [28,29]. Cao et al. (2007) [16] and Zhou et al. (2009) [17] developed a binary regression model (based on FSC from Landsat, NDSI and normalized difference vegetation index (NDVI)) and a segmentation model (based on fully and sparsely snow-covered areas, respectively); however, validation analyses have shown that those models are incapable of providing accurate estimates for extreme (e.g., fully and sparsely snow-covered) and transitional areas (the edge of snow cover extend) over the Tibetan Plateau. Regarding the spectral unmixing method, Painter et al. (2009) [19] proposed a MODIS snow-cover area and grain size retrieval algorithm (MODSCAG) to simultaneously retrieve snow grain size and FSC. The MODSCAG model is the most representative spectral unmixing algorithm based on MODIS data and its accuracy is also quite good [21]. However, the MODSCAG model derives FSC by regularizing shadows, which may generate relatively large numbers of meaningless negative values [22]. In addition, several hyperspectral imaging-based spectral unmixing algorithms have been gradually introduced into MODIS FSC mapping. Zhang et al. (2015) produced MODIS FSC maps using a fully constrained least squares spectral mixture analysis method (FCLS) [30], fully constrained scaled gradient method (FCSGM) [31], a sparse regression method (SPARSE) [32] and a polynomial nonlinear method (POLY) [33,34] for the Tibetan Plateau. Although they found that all of these algorithms can improve the accuracy of snow-cover mapping compared with the MOD10A1 FSC, they also showed that complex terrain can greatly reduce the accuracy of these spectral unmixing algorithms [35]. As a type of multi-factor machine learning method, back-propagation artificial neural network (BP-ANN) models have recently been introduced into MODIS FSC mapping. Dobreva et al. [26] (2009) examined the adaptability of a three-layer back propagation (BP) artificial neural network (ANN) for retrieving FSCs in regions with relatively smooth terrain in the Northern Hemisphere and found that the model exhibited relatively good adaptability. Hou et al. [27] (2014) employed a BP-ANN model to retrieve the FSC for the Heihe River Basin and found that the model exhibited a significantly better accuracy than that of the MOD10A1 FSC.

Unmanned aerial vehicles (UAVs) equipped with small digital cameras have several advantages, including the ability to acquire data with high spatial resolutions (up to centimeter scale), small volumes, low weights, low operational costs, high flexibility, the ability to transmit images in real time and the capability of being operated in high-risk areas [36,37,38]. They have been developed in recent years and have played unparalleled roles in eco-environmental monitoring, environmental protection, urban construction and planning and have been used to provide accurate estimates of agricultural production and real-time data acquisition during natural disaster (e.g., snowstorms and landslides) monitoring periods [39,40,41]. At present, UAVs equipped with digital cameras are critical for acquiring high-resolution remote sensing data, which serve as beneficial supplements to ground (camera) and high-altitude (satellite) observational data [42,43].

In this study, the Tibetan Plateau, which is a region with complex terrain and an unpredictable climate, was selected as the study region. MODIS linear regression (LR) and BP-ANN FSC retrieval models were developed based on MODIS and UAV data. In addition, linear spectral mixture analysis (LSMA) FSC maps were generated from the MOD09GA surface reflectance product using a pixel purity index (PPI) endmember extraction algorithm and FCLS method. The accuracies of the three commonly used models for retrieving FSC were evaluated and compared using Landsat 8 OLI data and MOD10A1 FSC with the goal of developing a high-accuracy FSC mapping algorithm suitable for the Tibetan Plateau, which would lay the foundation for accurately evaluating the distribution and variation of snow cover on the Tibetan Plateau and provide more accurate inputs for hydrological and climate models.

2. Study Area

The Tibetan Plateau reaches more than 8800 m in elevation at its highest point, has an average elevation exceeding 4000 m and encompasses a total area of approximately 2.5 million km², making it the highest plateau in the world (Figure 1) [44]. In general, the Tibetan Plateau has relatively low temperatures and is the largest region in China in which glaciers and snow cover are distributed [45]. The formation and evolution of the Tibetan Plateau have had profound impacts on the natural environment and human activity on it and in the nearby regions and they are of great significance in global change research. In recent decades, the Tibetan Plateau has drawn widespread attention and has been studied by the scientific communities in China and elsewhere [46]. As an important surface forcing factor, changes in snow cover on the Tibetan Plateau will inevitably exert significant impacts on regional climate, hydrology and ecology [47]. Therefore, accurately monitoring the snow-cover distribution on the Tibetan Plateau and producing high-accuracy FSC datasets are of great importance to global change research and disaster prevention and reduction [48,49].

3. Data

3.1. Unmanned Aerial Vehicle (UAV)

A field survey was conducted from 17 February to 2 March 2017 and the ground truth FSC data for the Tibetan Plateau were primarily collected using a UAV. The surveyed areas included seven counties and districts in the Tibetan Plateau region: Qilian, Chamdo, Darlag, Sogwo, Nyingchi, Yushu and Zekog. A total of fifteen sample plots with an area of 250,000 m² were surveyed (Figure 1). A DJI Inspire 1 Pro quadcopter UAV equipped with a GoPro Hero 3+ camera was used to vertically photograph the sample plots.

3.2. Landsat 8 Operational Land Imager (OLI)

Landsat 8 OLI data with a spatial resolution of 30 m were obtained from the website of the United States Geological Survey (USGS), from which the FSC data of the study region were extracted. Clear-sky Landsat 8 OLI data (Universal Transverse Mercator projections) from various time periods of nine scenes in the main snow-covered areas on the Tibetan Plateau were collected (

Table 1. Information of the Landsat 8 OLI data in this study.

Data Code	Path	Row	Date	Cloud Fraction
V1	144	36	14 October 2016	0.03
V2	134	40	11 December 2016	0.04
V3	132	41	12 January 2016	0.03
V4	133	38	4 December 2016	0.06
V5	136	35	4 December 2016	0.067
V6	146	36	24 November 2014	0.05
V7	140	36	2 October 2016	0.04
T1	141	34	23 October 2015	0.02
T2	134	33	6 October 2015	0.06

). Of the nine scenes, two (T1 and T2) were used to test the BP-ANN models developed in this study to select an optimal model and seven (V1–V7) were used to validate and compare the three FSC retrieval models and the MOD10A1 FSC (Figure 1).

3.3. MODIS Data

The MODIS Terra and Aqua data (including the MOD10A1/MYD10A1 daily snow-cover products, MOD09GA surface reflectivity product, MOD11A1 daily land surface temperature (LST) product and MOD13Q1 vegetation index product) [50,51] from the same areas and dates of the data acquired by the UAV were used to develop the three FSC retrieval models. The MOD10A1/MYD10A1 products (spatial resolution: 500 m, version 006) with sinusoidal projections, including, NDSI, snow albedo (SA), quality assessment (QA), orbit and granule data, were obtained from the NSIDC. The data we used were selected from the values whose quality assessment were 0 (best), 1 (good), 2 (OK). The NDSI data (denoted as MODIS NDSI) were used to establish the LR model in addition to being used as backup input variables for the BP-ANN models. In addition, FSCs retrieved from MODIS daily snow products V005 were only used in this study for comparisons with the other models. The MOD09GA, MOD11A1 and MOD13Q1 products (spatial resolutions: 500 m, 1000 m and 250 m) were obtained from the National Aeronautics and Space Administration website (NASA, US). The MOD09GA product provides estimates of surface spectral reflectivity measured in the absence of atmospheric scattering and absorption, including daily surface reflectivity data from MODIS channels 1–7 (R1–R7). Among them, the data of two bands of 4 and 6 are the main bands for calculating NDSI and snow mapping. They were used to extract the endmember spectra of the study region and decompose the MODIS mixed pixels as well as to provide input variables for the BP-ANN models to produce FSC maps. The MOD11A1 and MOD13Q1 products provide daytime LST data (LST_Day_1 km) and NDVI data, respectively, which were used to develop the BP-ANN models.

3.4. Digital Elevation Model (DEM)

Shuttle Radar Topography Mission (SRTM) DEM data with a spatial resolution of 90 m were obtained from the USGS [52]. After being resampled to a spatial resolution of 500 m using ArcGIS 10.2.2 software, the DEM data were used as input variables for the BP-ANN FSC models.

4. Methodology

4.1. Data Preprocessing

The images acquired by the UAV already with coordinates were then first mosaicked and processed into the same geographic coordinate system as the MODIS data and then orthorectified by Pix4Dmapper software. The spatial resolution of orthorectified image is 5 cm. The FSC data were then calculated using supervised classification via ENVI 5.3 software and were used as the true FSCs for the corresponding MODIS grid. The detailed process was as follows. First, the prepared UAV image was input into ENVI 5.3 software and the “ROI” option was clicked to draw the snow as one category and the other objects as another category. The parameters were then set using the “Supervised” tool in the software’s “Classification” module and the operation was carried out as follows. Finally, the classified binary image (snow/non-snow) was input into ArcGIS 10.2.2 software and the total number of pixels (represented by n_tt) and the number of snow pixels (expressed by n_ss) were calculated. n_ss/n_tt is the UAV FSC value.

The MODIS products (MOD10A1/MYD10A1, MOD09GA, MOD11A1 and MOD13Q1) for each scene were first subjected to format and coordinate conversions using the MODIS Reprojection Tool (MRT). The specific operations included the following procedures. The sinusoidal projections were converted to geographical coordinates using the World Geodetic System 84 ellipsoid, the nearest neighbor algorithm was selected as the resampling algorithm, the image files were converted to the GeoTIFF format and then the output data (NDSI, LST, (R1–R7) and NDVI) were obtained. The LST and NDVI data were resampled to spatial resolutions of 500 m in ArcGIS 10.2.2. Binary snow cover maps were then created from the Landsat 8 OLI data. The detailed steps are as follows. The Landsat 8 OLI images were subjected to radiometric calibration and atmospheric correction. The NDSI (denoted as Landsat NDSI) were calculated using Landsat 8 OLI bands 3 and 6 based on the Landsat 8 OLI band parameters and the SNOMAP algorithm [53], the NDSI threshold value was set to 0.4 and band 5 was set to >0.11 to eliminate interference from water bodies [54,55]. The 30-m-resolution binary snow-cover classification map (1: snow cover; 0: non-snow cover) was then generated [22]. Afterward, the binary snow cover maps generated from the Landsat 8 OLI images and the MODIS NDSI map were subjected to an overlay analysis. n_s represents the number of snow-covered pixels with values of 1 in the Landsat 8 OLI classification map contained in each MODIS grid and $n_t$ represents the number of 30-m pixels contained in each MODIS pixel using the “Fishnet” tool in ArcGIS 10.2.2 [17]. In general, 277 was the approximate number of $n_t$ on average. The 500-m-resolution FSC (denoted as the Landsat 8 FSC) was ultimately calculated:

$$\mathrm{FSC}=n_s/n_t,$$

(1)

where $n_s$ represents the number of snow-covered pixels with values of 1 in the Landsat 8 OLI classification map contained in each MODIS grid and $n_t$ represents the total number of 30-m pixels contained in each MODIS grid.

4.2. LR

Salomonson et al. (2004) [13] asserted that there a linear relationship between MODIS NDSI and true land surface FSC. By establishing an LR equation between the two factors, FSC can be calculated directly using NDSI. This method is easy to implement and understand. Based on this principle, we used the fifteen FSCs with an area of 250,000 m² acquired by the UAV as the true land surface FSC values (denoted as the UAV FSC) and the corresponding NDSI values were determined using MOD10A1 and MYD10A1 products. Finally, an LR model between UAV FSC and corresponding MODIS NDSI was established using a statistical method that was based on the fifteen samples.

4.3. Fully Constrained LSMA

The FCLS LSMA model is a relatively common SMA algorithm [17]. In that model, the reflectance of a pixel in a certain spectral band is a linear combination of the reflectance and proportionality coefficients of the compositional endmembers [30]. The unmixing process involves two steps. The first step is to extract the spectra of “pure” ground objects (i.e., endmember extraction) and the second step is to use a linear combination of the endmembers to represent the mixed pixels (i.e., SMA). Selecting suitable endmembers is the key to a successful SMA. The existing endmember extraction algorithms are mostly based on hyperspectral images and are mainly classified into convex analysis-based and statistical analysis-based algorithms [7]. Among the main constraint conditions of the FCLS algorithm are that (1) the sum of the components is 1 and that (2) each coefficient is non-negative. This results in the following equations:

$$Y_i={{\displaystyle \sum }}_{j=1}^m{r}_{ij}{f}_j+n_i,$$

(2)

$${{\displaystyle \sum }}_{j=1}^m{r}_{ij}=1,$$

(3)

$$0\le r_{ij}\le 1,$$

(4)

where $m$ represents the number of endmembers in a pixel, $Y_i$ represents the overall reflectance of the mixed pixel in band i, $r_{ij}$ represents the reflectance of endmember j in band i, $f_j$ represents the proportionality coefficient of endmember j in the mixed pixel and $n_i$ represents the spectral error in band i.

Based on the above principle, the fully constrained LSMA was established based on the R1–R7 from MOD09GA. The endmember spectra were first extracted using the PPI and an n-dimensional visualization tool. The FCLS algorithm was then used to decompose the MODIS mixed pixels and obtain the FSC extracted from endmember abundance images in the study area.

4.4. BP-ANN

A BP-ANN uses a gradient descent method to learn; it continuously adjusts its weight and threshold by backward propagation such that the sum of the squared errors reaches a preset value. A BP network adds several layers of neurons (one or more layers) between the input layer and the output layer. These neurons are called hidden units and have no direct connection with the outside world. However, their state changes can affect the relationship and between the inputs and outputs and each layer can have several nodes [56,57].

Due to the complex topography and climate in the study area, the distribution and attribute characteristics of snow cover result from various factors such as (R1–R7), DEM, LST, slope, aspect and NDVI [52]. This research was based on UAV and MODIS data to determine the relationships between UAV FSC and the different factors using a three-layered BP-ANN model. A total of 11 influencing factors, including surface reflectivity data (R1–R7) and four auxiliary sets of data (NDSI, NDVI, LST and DEM) were first selected as the alternative inputs for the BP-ANN models. Moreover, the processed UAV FSC data were used as the outputs for the BP-ANN models. The independence of the input variables was analyzed. The surface reflectivity data (R1–R7) represent the spectral reflectances of different ground objects in different wavebands of the sensor and are important symbols for distinguishing these features [58]. Therefore, to ensure the integrity of the data, (R1–R7) should all be used as the basis data of BP-ANN models for snow identification. Finally, we chose (R1–R7), DEM and LST as the final input variables by analyzing the inputs’ source and significance [59,60]. DEM and LST were used as auxiliary factors to explore the improvement of model performance degree on the basis of (R1–R7).

To obtain an optimal BP-ANN model for the study area, the inputs ((R1–R7), DEM and LST) were first normalized using the deviation standardization method (Equation (5)) for a uniform dimension. We then combined (R1–R7) with the other two auxiliary data sets (DEM and LST) in different ways to serve as inputs. Simulation experiments were performed with four combination schemes of nine input variables that excluded NDSI and NDVI. In Exp. 1, we used only (R1–R7) as inputs. In Exp. 2 and Exp. 3, DEM and LST were used as additional inputs to (R1–R7), respectively. In Exp. 4, DEM and LST were integrated as additional inputs to (R1–R7). The relationship models between the various inputs and true values (i.e., UAV FSCs) were trained using 60% of the samples (nine) and the remaining 40% (six) were used to validate the relationships. All fifteen samples were then used to build the BP-ANN model. The detailed network parameters were set as follows. The learning rate was set to 0.05 (usually between 0.01 and 0.1), the momentum factor was 0.95 (general values of 0.95), the maximum number of failures was 100, the number of iterations was 1000, the minimum value of the error was 0.0001 and the hidden-layer neuron number was 2n + 3 (n is the number of input nodes). The weights between the input and the hidden layers and the hidden and the output layers were all initialized to zero [27].

Table 2. Information on the parameters of the BP-ANN model.

ANN Parameters	Description
Input	R1–R7, DEM, LST
Output	FSC
Number of net layers	Three (the input layer, hidden layer and output layer are only one layer)
Number of nodes	Input layer: n; hidden layer: 2n + 3; output layer: 1
Training and learning method	Levenberg-Marquardt, gradient descent method
Performance	RMSE, r
Transfer function	Input layer to hidden layers: tangent hyperbolic
	Hidden layers to output layer: pure linear
Standardization method	Standard deviation method

* Note: n denotes the number of inputs; the number of hidden layers, 2n + 3, is based on empirical values from previous studies [27].

lists the basic parameters used in the models. The representation capabilities of the BP-ANN models were evaluated using RMSE (Equation (6)) and r (Equation (7)). Four of the relatively better models based on the different experiments were preserved and we tested the four BP-ANN models on the test sub-regions (T1 and T2) based on the Landsat 8 FSC of the same date for further screening. Finally, one optimal BP-ANN model suitable for the Tibetan Plateau was eventually selected based on the comprehensive evaluation results with respect to the RMSE and r. To exactly follow how the three types of models (LR, LMSA and BP-ANN models) were developed easily, the development process of three models is shown in

Table 3. The development processes of the LR, LMSA and BP-ANN models for the Tibetan Plateau.

		LR Model	LMSA Model	BP-ANN Model
The procedure of model establishment	Data	MOD10A1/NDSI, MYD10A1/NDSI, UAV/FSC	MOD09GA	R1–R7, DEM, LST, UAV/FSC
	Date	February–March 2017	October 2014–December 2016	February–March 2017
	Tools	Excel 2013	ENVI 5.3	Matlab 2016
	Process	1. Scatter plot (FSC, NDSI) 2. FSC = a × NDSI + b	1. Spectral library building 2. Subpixel unmixing algorithm: FCLS	This ANN was trained with back propagation error to learn the relationship between UAV FSC and R1–R7, DEM and LST
Examples of FSC spatial distribution	Date	14 October 2016
	Region	V1 (randomly selected)
Validation	Data	Binary Landsat 8 OLI snow-cover map
	Date	October 2014–December 2016
	Tools	ArcGIS 10.2.2, Excel 2013
	Region	V1–V7

^* Note: a and b denote undetermined coefficients.

.

$$X=\frac{x-min}{max-min},$$

(5)

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{(y_i-{y\prime }_i)}^2}{n-1}},$$

(6)

$$r=\frac{{{\displaystyle \sum }}_{i=1}^n\left(y_i-\overline{y_i}\right)\left({y\prime }_i-\overline{{y\prime }_i}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{{y\prime }_i}\right)}^2}·{{\displaystyle \sum }}_{i=1}^n{({y\prime }_i-\overline{{y\prime }_i})}^2},$$

(7)

where $max$ represents the maximum value of the sample data, $min$ represents the minimum value of the sample data, $y_i$ represents the true value, ${y\prime }_i$ represents the simulated value, $i$ represents a sample, $n$ represents the number of samples, $\overline{y_i}$ represents the mean hypothetical true value and $\overline{{y\prime }_i}$ represents the mean simulated value.

4.5. Validation

To validate the retrieval accuracies of the three models, a statistical analysis was performed using representative pixels selected from the validation sub-regions (V1–V7) in conjunction with the Landsat 8 FSC and MO10A1 FSC of version 005 of the same date. Six accuracy evaluation indices were used, including the mean absolute error (MAE) (Equation (8)), RMSE, r, mean FSC (MFSC) (Equation (9)), positive mean error (PME) (Equation (10)) and negative mean error (NME) (Equation (11)).

$$\mathrm{MAE}=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left|y_i-{y\prime }_i\right|,$$

(8)

$$\mathrm{MFSC}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{p}_i,$$

(9)

$$\mathrm{PME}=\frac{1}{m}{{\displaystyle \sum }}_{j=1}^m(y_j-{y\prime }_j)(y_j>{y\prime }_j),$$

(10)

$$\mathrm{NME}=\frac{1}{r}{{\displaystyle \sum }}_{k=1}^r(y_k-{y\prime }_k)(y_k<{y\prime }_k),$$

(11)

where $p_i$ represents the FSC in the MODIS pixel, r and m represent the numbers of pixels whose true values are smaller than the simulated values and the number of pixels whose true values are greater than the retrieved values, respectively, $y_j$ and $y_k$ represent the true values (Landsat 8 FSC) of each pixel and ${y\prime }_j$ and ${y\prime }_k$ represent the simulated values.

5. Results

5.1. LR

Figure 2 shows the fifteen snow maps with an area of 250,000 m² obtained from the UAV device and the corresponding FSC and NDSI were also counted and labeled separately. It can be seen that the FSCs varied between 0 and 1 and that the images contained several underlying surface types (including snow cover, vegetation, bare land and rivers). A linear regression analysis between the MODIS NDSI and UAV FSC was performed based on the sampled data and the results are presented as a scatter diagram (Figure 3). There was a significant correlation between NDSI and FSC, the data points were distributed around a straight line, which demonstrates the good goodness-of-fit of the LR model (R² = 0.6864). The LR model (Equation (12)) was then used to retrieve the FSC (denoted as Reg-FSC) for validation sub-region V1 (randomly selected) on 14 October 2016. The results, which were compared with the corresponding retrieval results of the Landsat 8 FSC and MOD10A1 FSC (Figure 4), showed that the LR model results were similar to the Landsat FSC, except for some small errors. Slightly more snow was identified.

$$FSC=0.8286\times NDSI+0.3941,$$

(12)

where $FSC$ represents the UAV FSC and $NDSI$ represents the corresponding MODIS NDSI.

5.2. LSMA

The results for validation sub-region V1 were taken as an example to illustrate the endmember spectra that were extracted using the PPI and n-dimensional visualization tool (Figure 5). The most representative six endmembers (snow, cloud, bare land, soil, vegetation and water), compared with theoretical wavelength dependent reflectances of those land cover types, were screened [35]. To verify that the established spectral library can be used to produce a more accurate fractional snow-cover map, the endmember abundance images were obtained using the FCLS algorithm. Figure 6 shows the snow abundance image and RMS error in validation sub-region V1. The Landsat 8 FSC was used as the ground truth for comparison with the FSC generated by the LMSA model (denoted as Mix-FSC) and with the MOD10A1 FSC (Figure 7). It can be seen that the LSMA model produced omission errors in some areas, i.e., the LSMA model identified some snow-covered pixels as snow-free pixels or as pixels with lower FSCs. However, compared with the MOD10A1 FSC, the retrieval results of the LSMA model were closer to the Landsat 8 FSC in validation sub-region V1.

5.3. BP-ANN

Table 4. r and RMSE statistical values between the simulated FSCs of the four BP-ANN experimental schemes and Landsat 8 FSCs in the test sub-regions (T1–T2).

Experiments	Combinations	RMSE	r
Exp. 1	(R1–R7)	0.2606	0.7884
Exp. 2	(R1–R7) + DEM	0.2071	0.8053
Exp. 3	(R1–R7) + LST	0.2267	0.8216
Exp. 4	(R1–R7) + DEM + LST	0.2006	0.8453

shows the r and RMSE values between the simulated FSCs of the four BP-ANN experimental schemes and corresponding Landsat 8 FSCs in the test sub-regions (T1–T2). The r values ranged from 0.7884 to 0.8453 (mean of 0.8152) and the RMSEs ranged from 0.2006 to 0.2606 (mean of 0.2238). It can be seen that the performance of adding a single auxiliary factor was better than the model with only surface reflectivity data. The effects of DEM and LST on snow cover were relatively close based on comprehensive evaluation of r and RMSE. The BP-ANN model trained using Exp. 4 (input variables: (R1–R7) + DEM + LST) had a better r (0.8453) and a smaller RMSE (0.2006) than those of the models trained using the other experimental schemes. Thus, the BP-ANN model using the input variable combination of Exp. 4 was the optimal BP-ANN FSC prediction model for the study area. Figure 8 shows the fitting result between the simulated values of the optimal BP-ANN model and true values (UAV FSCs) based on all fifteen samples. It was found that the true (UAV FSC) and simulated values exhibited a good correlation coefficient (0.9001) and a small RMSE value (0.0707). The simulation results combined with Table 4 showed that the optimal BP-ANN model was a good FSC retrieval method, as expected. Overall, the nine key variables were suitable to be used to build a reasonable and effective BP-ANN model for the study area. According to Figure 9, the optimal BP-ANN model retrieved more snow than the MOD10A1 FSC, although its spatial distribution (denoted as BP-FSC) was closer to that of the Landsat 8 FSC in the snow-covered region.

5.4. Validation

Table 5. MFSC statistics for the three FSC models and the MOD10A1 FSCs in the validation sub-regions (V1–V7).

	MFSC/%	Reg-FSC	Mix-FSC	BP-FSC	MOD10A1-FSC	Landsat 8-FSC
Code
V1		48.08	33.92	51.73	29.67	40.05
V2		43.61	8.21	6.63	16.16	12.94
V3		31.80	9.67	9.14	30.68	26.25
V4		17.48	8.23	26.04	7.58	19.60
V5		36.05	38.08	65.60	39.02	62.29
V6		15.29	2.39	18.71	4.38	7.34
V7		44.38	19.71	26.20	28.11	33.36
Average		33.81	17.17	29.23	22.23	28.83

summarizes the MFSC using the three models, Landsat 8 FSC and MOD10A1 FSC for all the validation sub-regions (V1–V7). The mean FSC values of the three FSC models (Reg-FSC, Mix-FSC and BP-FSC) and MOD10A1 were 33.81%, 17.17%, 29.23% and 22.23%, respectively. The mean FSC extracted from the Landsat 8 OLI data was 28.83%. Thus, it can be concluded that the FSCs retrieved using the LR model and the BP-ANN model were similar to the Landsat 8 FSCs.

Further accuracy analysis results are shown in

Table 6. r, RMSE, MEA, PME and NME statistics for the three FSC models and the MOD10A1 FSCs in the validation sub-regions (V1–V7).

	r	RMSE	MAE	PME	NME
MOD10A1	0.7838	0.2833	0.1877	0.1867	−0.2081
Reg-FSC	0.8222	0.2304	0.1765	0.1431	−0.2570
Mix-FSC	0.7921	0.3485	0.2755	0.3411	−0.0902
BP-FSC	0.8445	0.2201	0.1039	0.1732	−0.1055

. The calculated values for MAE, RMSE, r, PME and NME (accuracy evaluation indices) between the three FSC simulated values, MOD10A1 and Landsat 8 FSC of the representative pixels in all validation sub-regions (V1–V7) are listed in Table 6. The results show that the performance of the optimal BP-ANN model was better than that of the other FSC models and MOD10A1 FSC based on comprehensive evaluation of accuracy indices, especially for its best r (0.8445), lowest RMSE (0.2201) and far lower MAE (0.1039), which indicates that the BP-ANN model had the best accuracy and stability. However, these values were slightly lower than the results of Hou’s algorithm used for the Heihe River Basin due to small sample limitations [27]. The LR model, which had a slightly larger NME (−0.2570), had the second-best accuracy (r = 0.8222 and RMSE = 0.2304) among the three models and MODIS FSC. The optimal BP-ANN and LR have nearly identical r and RMSE values, having much better MAE and NME values makes BP-ANN model better than LR model. The retrieval accuracy of the LSMA model was the worst. The r was only 0.7921 and the RMSE was as high as 0.3485; there were some serious omission phenomena in the study area, specifically given its values of MAE (0.2755) and PME (0.3411), which were the largest. Based on the six accuracy evaluation indices, the accuracy of the MOD10A1 FSC product was much lower than those of the LR and BP-ANN models, although its accuracy was slightly better than those of the LSMA and previous similar research results (r = 0.7400) in this area [22]. In short, based on a comprehensive evaluation of the six indices, the optimal model was the BP-ANN model, which can effectively improve the accuracy of snow area monitoring; it is the ideal algorithm for retrieving FSC on the Tibetan Plateau.

6. Discussion

The Tibet Plateau is an important seasonal snow-capped area in China and is of great significance for accurately developing estimates of snow in the area [2,7,9]. In view of the current poor monitoring accuracy of global MODIS snow products for the Tibetan Plateau [35], we developed the LR and BP-ANN models based on MODIS and UAV data, respectively and the LMSA model using PPI, an n-dimensional visualization tool and the FCLS algorithm, based on MODIS data only for this study. The FSC derived from higher spatial resolution Landsat 8 OLI images (Landsat 8 FSC) were then used as ground truth data to validate the accuracy of the three models for the purpose of finding the optimal algorithm for FSC extraction for the Tibetan Plateau. In addition, MOD10A1 FSCs (version 005) were used in this study for comparison with the three models. The correlation analysis of alternative input variables for the BP-ANN model, NDVI and NDSI were calculated separately from bands 1, 2, 4 and 6 of the surface reflectivity data, which may have resulted in some redundancies associated with (R1–R7). As a consequence, we finally chose (R1–R7), DEM and LST as the final input variables for the BP-ANN model [58,59,60]. The results showed that the spectral curve extracted using the PPI and the n-dimensional visualization tool for “bare land” was very different from that of “soil” in validation sub-region V1 (Figure 5) because soil typically contains much more water than bare land [22,35]. A BP-ANN network, which was based on input combination of (R1–R7) + DEM + LST, was trained to establish a fractional snow-cover retrieval model. The network had ideal predictive capability and generalization capacity to meet the requirements of FSC monitoring among the three models and MOD10A1 FSC for the study area. It can easily incorporate auxiliary information (elevation and temperature) on the basis of on the surface reflectivity data during the relationship training period, which can effectively improve the accuracy of snow cover mapping. In short, because the use of UAV data near the ground truth and the ability of the BP-ANN model to integrate various input factors make it good enough. Compared to the similar algorithm of Hou that was employed for the Heihe River Basin [27], the BP-ANN FSC retrieval model based on fifteen samples from the UAV data had a slightly lower accuracy based on a comprehensive evaluation of the accuracy indices, the main reason for which may be that Hou’s algorithm was based on tens of thousands of samples from Landsat data during the training process. Due to the cold, high elevations and fragmented snow distribution on the Tibetan Plateau, it is very difficult to obtain enough ground samples by UAV (a total of fifteen in this study) and to validate the model, we first set the training samples to a lesser extent, herein 60% (corresponding to nine samples) and the remaining 40% were used to verify the relationship samples (the corresponding number of samples was six). All of the samples were then used to build the BP-ANN model in this study. Once BP-ANN is trained, it only needs to reference the saved connection weights, the FSC is retrieved through multiplication and summation. So, it is possible to quickly expand the training of the ANN for better regional applicability by adding training scenes from other snow-covered regions and other periods [27]. Which as a contribution makes the BP-ANN overall better than the LR model for the FSC mapping of the Tibetan Plateau.

The accuracy of the empirical LR model is better than that of the LSMA model based on the spectral analysis in this study, which is an opposite conclusion to those of previous studies [7,19]. The main reason may have been the introduction of UAV aerial data (with spatial resolutions up to one centimeter) during LR model building, unlike the previous studies, which were based on 30-m-resolution Landsat data [8,17]. In addition, other research has shown that grain size and snow pollution are among the important factors that cause albedo changes in snow cover; grain size affects the albedo of snow near the infrared band, snow pollution affects the albedo of snow in visible light [61,62,63] and bidirectional reflectance factors (BRFs) play a key role in snow cover [64]. All of these factors may affect the LMSA model spectral analysis process, which caused the errors of LMSA model to intensify in this study and the accuracy was slightly lower than in previous research results [35].

In general, the MODIS FSC is defined as the proportion of snow per MODIS pixel. It is usually difficult to accurately calculate the true FSC value on surface owing to the influence of various ground obstacles, even the field investigation. In response to this problem, we introduced the UAV equipment in this study and used the results of its aerial photography as the true value of the field FSC. On the one hand, the spatial resolution of UAV can reach centimeter levels, the photos taken are very close to the actual situation of the ground and can be used to maximize the integrity of the preservation information; and in addition, areas where access by humans is difficult and dangerous can be photographed. However, the UAV data have significantly higher resolutions than the spatial resolution of Landsat data, which caused the LR and BP-ANN models to tend to retrieve more snow-cover information in some areas.

The introduction of UAV data has injected new approaches and methods for MODIS FSC and the BP-ANN method, which can approximate any nonlinear function while providing clear physical and conceptual results based on a flexible and changeable topological structure [56], is a widely effective FSC retrieval model for the study area. There are still some problems that require attention. First, the approximation and generalization abilities of the BP-ANN model are closely related to the learning samples and small datasets may be poorly representative and contain conflicting and redundant samples, which may cause higher failure rates during the training period [57]. At present, some researchers have reused the BP-ANN techniques for practical application. For example, Wu et al. (2007) [65] and Hou et al. (2014) [27] repeatedly used their BP-ANN technologies in snow disaster warning and mountain snow area prediction, respectively; Zheng et al. (2015) [66] used BP-ANN model to predict change trend of snow cover in the Manasi River Basin in 2012 and applied the results to runoff forecast in 2017 [67]. The BP-ANN model is an optimal method and can be used to snow mapping and provided basic information for application and decision-making for environmental departments and government agencies, such as snow disaster warning, water resources management, agricultural and animal husbandry arrangement, etc. Furthermore, there are still a large number of mixed pixels in the binary snow map obtained by the 30-m spatial resolution Landsat 8 OLI images of the Tibetan plateau, where the terrain is complex and the snow is thin. The combination of multiple Landsat pixels represents a MODIS pixel and the largest error comes from its spatial influence when it is used as verification data.

7. Conclusions

The Tibetan Plateau was selected as the study area in this work. Based on mixed pixels, three MODIS fractional snow-cover retrieval models (LR, LSMA and BP-ANN models) were constructed based on MODIS and UAV data. The FSCs derived from higher spatial resolution Landsat 8 OLI images (Landsat 8 FSC) were used as ground truth data to validate the accuracies of the three models for the purpose of finding the optimal algorithm for FSC extraction for the Tibetan Plateau. MOD10A1 FSCs (version 005) were also used in this study for comparisons with the other models. The main conclusions are now summarized.

(1)

The optimal input combination of the BP-ANN model obtained from this study was ((R1–R7) + DEM + LST), based on alternative input variables. The BP-ANN model developed using that input combination exhibited better performance than the models trained using the other combinations because the elevations and temperatures, which affected the FSC distribution, were integrated based on an analysis of surface reflectivity.

(2)

The accuracies of the LR and BP-ANN models based on MODIS and UAV data were much better than those of the fully constrained LSMA model based only on MODIS data and the MOD10A1 FSC. The retrieval accuracy of the fully constrained LSMA model using the PPI endmember extraction algorithm was the worst among the three models (r was only 0.7921 and the RMSEs were as high as 0.3485). There were some serious omission phenomena in the study area, specifically because it had the largest MAE (0.2755) and PME (0.3411). However, it was slightly lower than the retrieval accuracies obtained in the MOD10A1 V005 and previous research [35] based on comprehensive evaluation of six accuracy indices. This finding indicates that the introduction of UAV data can improve the accuracy of FSC mapping.

(3)

In brief, the accuracy and stability of the optimal BP-ANN model were the best based on a comprehensive evaluation of the six indices employed in this study, especially for MAE (0.1039) and NME (−0.1055), which were much better than the LR model’s results (MAE = 0.1765, NME = −0.2570). These results indicated it an ideal prediction algorithm for retrieving FSC for the Tibetan Plateau when compared with the LR model, LMSA model and the MODIS global FSC product. The algorithm lays the foundation for accurately evaluating the distribution and variation of snow cover on the Tibetan Plateau and provides more accurate inputs for hydrological and climate models.

This is the first time that the UAV data has been introduced into MODIS FSC mapping and has achieved good results. In the future, we will continue to pay attention to it and study the advantages of UAV data in snow monitoring and find ways to solve the problem of BP-ANN model encountered during the building process, such as factor selection and screening, small sample, parameter setting problems and so on. Especially for spectral mixing, how to introduce the UAV data into spectral mixing method and the theory about how the grain size, BRF and black carbon effect the snow reflectivity and albedo, there is still additional exploration needed.