Permafrost Presence/Absence Mapping of the Qinghai-Tibet Plateau Based on Multi-Source Remote Sensing Data

The Qinghai-Tibet Plateau (QTP) is known as the Third Polar of the earth and the Water Tower of Asia, with more than 70% of the area on the QTP is covered by permafrost possibly. An accurate permafrost distribution map based on valid and available methods is indispensable for the local environment evaluation and engineering constructions planning. Most of the previous permafrost maps have employed traditional mapping method based on field surveys and borehole investigation data. However their accuracy is limited because it is extremely difficulties in obtaining mass data in the high-altitude and cold regions as the QTP; moreover, the mapping method, which would effectively integrate many factors, is still facing great challenges. With the rapid development of remote sensing technology in permafrost mapping, spatial data derived from the satellite sensors can recognize the permafrost environment features and quantitatively estimate permafrost distribution. Until now there is no map indicated permafrost presence/absence on the QTP that has been generated only by remote sensing data as yet. Therefore, this paper used permafrost-influencing factors and examined distribution features of each factor in permafrost regions and seasonally frozen ground regions. Then, using the Decision Tree method with the environmental factors, the 1 km resolution permafrost map over the QTP was obtained. The result shows higher accuracy compared to the previous published map of permafrost on the QTP and the map of the glaciers, frozen ground and deserts in China, which also demonstrates that making comprehensive use of remote sensing technology in permafrost mapping research is fast, macro and feasible. Furthermore, this result provides a simple and valid method for further permafrost research.


Introduction
Permafrost is defined as ground (soil or rock and included ice and organic material) where temperature have remained at or below 0 • C for a period of least two consecutive years [1,2], which is an key component of cryosphere and assessing the impact of the climatic change. It is not only a crucial component in the soil biogeochemistry, local hydrology, ecosystems, heat and energy exchange, but also the sensitive indicator of global warming [3]. The Qinghai-Tibet Plateau (QTP) known as "the the classification thresholds of remotely sensed variables in decision tree analysis are defined by using the areas of agreement from two previous published permafrost maps (1:3,000,000 Map of Permafrost on the QTP Frozen soil and the permafrost map of the QTP derived from 1:4,000,000 Map of the Glaciers, Frozen Ground and Deserts in China), and remote sensing data from 2003 are used in order to reduce time gap with the two previous published permafrost maps. Finally, the result is verified by Western Kunlun, Wenquan and Gaize three areas. The consistency of permafrost distribution is evaluated by overall accuracy and Kappa coefficients among the different published permafrost maps over the QTP.

Study Area
The QTP (73.19 • E~104.47 • E, 26 • N~39.47 • N) area is about 2.62 × 10 6 km 2 , and most of the land surface is estimated to be within permafrost zones. The QTP is considered to be the youngest and most complex plateau in the world. The elevation of the QTP is between 3000-5000 m asl and the average elevation is more than 4000 m asl. It is composed with plenty of mountains and plateaus; the mountains of the QTP are uneven and have great differences, for example, the average elevation of the Himalayas is about 6000 m asl, while the elevation of the Yarlung Zangbo River valley plain is only 3000 m asl. In addition, the mountain series of the QTP can be divided into east-west and south-north directions, the east-west mountains are the most prevalent and occupy the largest area; the south-north mountain are mainly distributed in the southeast of the plateau and near the Hengduan Mountains. These two groups of mountains form the basic geomorphological framework of the QTP [33]. The MAAT in the hinterland of the plateau is below 0 • C, and the monthly average temperature is less than 10 • C. For example, the southern and central QTP is the warm region, where air temperature is lower than 10 • C for nine months in a year; the Qiangtang Plateau and Hoh Xil is the cold region of the QTP, where air temperature is lower than 0 • C for seven months and lower than −10 • C for three or four months in a year [34]. The unique topography of the QTP determines the distribution characteristics of climate and vegetation; from the northwest to the southeast, there are mainly three ecological systems: alpine meadow, alpine steppe, alpine swamp and wetland [35]. The vegetation type of Northwest Plateau is the alpine, grassland and non-vegetation; the center of the QTP is meadow steppe and meadow with a high vegetation cover in this region. Generally, changes in vegetation can cause changes in ground surface conditions, such as evapotranspiration, albedo, and infiltration rate of precipitation, which can indirectly influence the thermal dynamics of permafrost [36], however, NDVI (Normalized Differential Vegetation Index) as the best indicator of vegetation growth state and vegetation coverage is widely used in vegetation remote sensing, it shows NDVI can be used as one of the parameters reflecting the frozen ground. In addition, the basic permafrost types of the QTP are low-latitude and high-altitude plateau permafrost, and the permafrost in the Qilian Mountains, Himalaya Range, Kunlun Mountains, Hengduan Mountains, belongs to a middle-low latitude and high altitude alpine permafrost, which have obvious zonality, especially vertical zonality (e.g., the lower bound of permafrost lowered with the latitude increasing and the thickness of frozen soil increases with the altitude increasing).

Investigated Regions and Its Permafrost Maps
In order to evaluate the accuracy of this study, three investigated regions ( Figure 1)-Western Kunlun (WKL), Wenquan (WQ), and Gaize (GZ) conducted from 'the comprehensive investigation of permafrost and its environments on the QTP' are selected as verification data. WKL (78.8 • E~79.5 • E, 35.6 • N~36.1 • N) (3508.1 km 2 ) is located in the northwestern margin of the QTP [37], and belongs to the typical mid-latitude alpine permafrost zone. GZ (84 • E~86.03 • E, 32.27 • N~34 • N) is located in the south the QTP and northeast of the Gaize county [38], with altitude is between 4400 m and 6200 m asl and with area is 41,256 km 2 . WQ (99 • 60'E~99 • 42'E, 35 • 60'N~35 • 42'N) is located in the east of the QTP, area is about 2250 km 2 , the permafrost area accounted for the 66.7% of entire region [39]. Three investigated regions are located in the transition zones between permafrost and seasonally frozen ground with different climatic and geographic conditions. Permafrost maps of three investigated regions (IRs) are based on field survey data by employing geophysical exploration (Ground Penetrating Radar, GPR; Time-domain ElectroMagnetic, TEM), mechanical excavation, ground-based observations, and drilling method. Due to the elevation and terrain factors have much greater influence on permafrost presence/absence than that of longitude and latitude in the local region, those IRs permafrost maps were generated using the criteria of lower limits of permafrost in different conditions combined the DEM data and terrain factors (slope and aspect), and a portion of the observed results of geophysical methods and boreholes was reserved to validate the maps [37,38,40]. The validation result of WQ investigated regions permafrost map shows an accuracy of 85.6% [39]. and 6200 m asl and with area is 41,256 km 2 . WQ (99°60'E~99°42'E, 35°60'N~35°42'N) is located in the east of the QTP, area is about 2250 km 2 , the permafrost area accounted for the 66.7% of entire region [39]. Three investigated regions are located in the transition zones between permafrost and seasonally frozen ground with different climatic and geographic conditions. Permafrost maps of three investigated regions (IRs) are based on field survey data by employing geophysical exploration (Ground Penetrating Radar, GPR; Time-domain ElectroMagnetic, TEM), mechanical excavation, ground-based observations, and drilling method. Due to the elevation and terrain factors have much greater influence on permafrost presence/absence than that of longitude and latitude in the local region, those IRs permafrost maps were generated using the criteria of lower limits of permafrost in different conditions combined the DEM data and terrain factors (slope and aspect), and a portion of the observed results of geophysical methods and boreholes was reserved to validate the maps [37,38,40]. The validation result of WQ investigated regions permafrost map shows an accuracy of 85.6% [39].

Elevation
Elevation is the main controlling factor of vertical zoning in the QTP, which cannot only reflect morphological features of terrain, but also extract much surface information, including slope, aspect, and topography. Recently, digital elevation model (DEM) as a conventional terrain tool is usually used in the research of earth's surface topography, geomorphology, surface structure analog, therefore this study uses DEM dataset ( Figure 2a) with a resolution of 1 km × 1 km that is projected to Albers Equal Area Conic in order to ensure the area is not deformed. DEM datasets of the QTP are from "Environmental and Ecological Science Data Center for West China" (http://westdc.westgis.ac.cn).

Elevation
Elevation is the main controlling factor of vertical zoning in the QTP, which cannot only reflect morphological features of terrain, but also extract much surface information, including slope, aspect, and topography. Recently, digital elevation model (DEM) as a conventional terrain tool is usually used in the research of earth's surface topography, geomorphology, surface structure analog, therefore this study uses DEM dataset ( Figure 2a) with a resolution of 1 km × 1 km that is projected to Albers Equal Area Conic in order to ensure the area is not deformed. DEM datasets of the QTP are from "Environmental and Ecological Science Data Center for West China" (http://westdc.westgis.ac.cn).

MODIS LST Products
MODerate-resolution Imaging Spectroradiometer (MODIS) sensors equipped with Terra and Aqua satellites provide multi band Land Surface Temperature (LST) data. LST is retrieved by split-window algorithm in infrared band [41,42]. A piecewise function [43] of ground temperature shows that the surface temperature takes on a sinusoidal function in daytime while in nighttime it declines almost linearly. Therefore, daily mean LST can be calculated with four instantaneous observations from the MODIS LST products (MYD11A1 and MOD11A1 with 1 km resolution) every day by integral sum of Sin-Linear piecewise function. Finally, the mean annual surface temperature (MAST) of 2003 ( Figure 2b) is equal to the arithmetic mean value of daily LST in 2003.

NDVI
Vegetation has scattering and reflection effects on solar radiation, which can reduce the amount of surface net radiation and lowered the near surface temperature. Furthermore, it can also intercept snow, preserve heat and increase transpiration, which indirectly affects the spatial distribution of frozen soil [44]. Therefore, in some areas, vegetation was used as an indicator for permafrost classification, for example, Grams [45] deemed that white spruce and birches could indicate permafrost boundary in Alaska valley, Nguyen [26] used vegetation types as an indicator of permafrost or seasonally frozen ground to classify the local permafrost distribution boundary. NDVI is often used to describe the parameters of the quantity and quality of vegetation. Hence, this study select NDVI as one of the parameters related to the permafrost distribution. NDVI data was extracted from the synthesized product data (MOD13A1 and MYD13A1), the image resolution was resampled to 1 km, and a little missing data was filled through the Kriging interpolation method. Finally, the NDVI annual maximum value (Figure 2c) of the QTP was calculated by the band computing tool of ENVI 5.0 (The Environment for Visualizing Images) software.

MODIS LST Products
MODerate-resolution Imaging Spectroradiometer (MODIS) sensors equipped with Terra and Aqua satellites provide multi band Land Surface Temperature (LST) data. LST is retrieved by split-window algorithm in infrared band [41,42]. A piecewise function [43] of ground temperature shows that the surface temperature takes on a sinusoidal function in daytime while in nighttime it declines almost linearly. Therefore, daily mean LST can be calculated with four instantaneous observations from the MODIS LST products (MYD11A1 and MOD11A1 with 1 km resolution) every day by integral sum of Sin-Linear piecewise function. Finally, the mean annual surface temperature (MAST) of 2003 ( Figure 2b) is equal to the arithmetic mean value of daily LST in 2003.

NDVI
Vegetation has scattering and reflection effects on solar radiation, which can reduce the amount of surface net radiation and lowered the near surface temperature. Furthermore, it can also intercept snow, preserve heat and increase transpiration, which indirectly affects the spatial distribution of frozen soil [44]. Therefore, in some areas, vegetation was used as an indicator for permafrost classification, for example, Grams [45] deemed that white spruce and birches could indicate permafrost boundary in Alaska valley, Nguyen [26] used vegetation types as an indicator of permafrost or seasonally frozen ground to classify the local permafrost distribution boundary. NDVI is often used to describe the parameters of the quantity and quality of vegetation. Hence, this study select NDVI as one of the parameters related to the permafrost distribution. NDVI data was extracted from the synthesized product data (MOD13A1 and MYD13A1), the image resolution was resampled to 1 km, and a little missing data was filled through the Kriging interpolation method. Finally, the NDVI annual maximum value (Figure 2c) of the QTP was calculated by the band computing tool of ENVI 5.0 (The Environment for Visualizing Images) software.

Soil Moisture
When the seasons change, the surface wet and dry degree in frozen ground areas have a large difference in time and space. During the freezing and thawing processes of the permafrost active layer near the WuDaoliang of the QTP [46], freezing process usually started in mid-September until completing in late October. In the following year, the freeze-thaw process near the surface occurred in late March, and in early April a stable process of thawing in daytime and freezing in nighttime would form, but the thawing of the active layer of permafrost started about the end of April from the top downwards. In May and June, the surface temperature of the QTP gradually increased, and the active layer began to thaw slowly, resulting in soil volume water increasing. Furthermore, the water content began to significantly change in May, and satellite sensors could easily capture these changes, AMSR-E (Advanced Microwave Scanning Radiometer-EOS) mounted on Aqua satellites provides 1-5 cm soil moisture data (National Snow and Ice Data Center) stored in HDF-EOS format with a resolution of 25 km. hence the average value of AMSR-E soil water content in May 2003 (Figure 2d) was obtained as an auxiliary classification information source due to its resolution was coarse compared to other classification data, such as DEM, LST and NDVI.

Methods
The regression relationship between the permafrost distribution and influencing factors (elevation, LST, NDVI, SM, albedo et al.) were analyzed by Wang [47], the results showed that the correlation coefficients of four factors (elevation, LST, NDVI and SM) with permafrost are the highest, and elevation is the governing factor. Therefore, elevation, LST, NDVI and SM, are selected as input parameters of decision tree.

Decision Tree Establishment Process
Decision tree algorithm is a typical classification method, it performs a recursive binary partitioning of the feature space [48,49]. Each decision tree has only one root node, but there can be several internal nodes and two or more terminal nodes. Generally, discriminant functions derived from training sample are indispensable because establishment of the next branch or sub-node need to take different discriminant functions as prerequisite until expected decision tree is generated. The following (Sections 3.1.1 and 3.1.2) presents the generation of the discriminant function and the establishment of the decision tree of this study.

Classification Threshold Value
In order to understand distribution features of each influencing factor in permafrost regions and seasonally frozen ground regions, this study selected the two published and widely used permafrost distribution maps [7,10,[50][51][52] as prior information to gain classification threshold of decision tree, which are: (1) The 1:3,000,000 Map of Permafrost on the QTP (Figure 3a), which was compiled by State Key Laboratory of Frozen Soil Engineering, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences [15] (hereafter referred to as "QTP96_map") based on permafrost investigated data, literature, aerial photographs, satellite and others relevant maps; it shows the permafrost area of TP is 1.29 × 10 6 km 2 ; (2) 1:4,000,000 Map of the Glaciers, Frozen Ground and Deserts in China [53], was compiled by following material or means: TOPO30 DEM with 1 km resolution, air temperature and ground temperature, field investigation of frozen soil, borehole data exploration, and interpretation of aerial photographs and satellite images. The permafrost distribution on the QTP used the research result of mean annual ground model, firstly, Nan [54] made regression analysis between mean annual ground temperature and latitude, elevation with the field survey data of 76 boreholes along the Qinghai-Tibet Highway; secondly, the mean average ground temperature over the entire QTP was calculated with elevation and latitude data; lastly, 0.5 • C was taken as the classified value between permafrost and seasonally frozen ground, and finished the permafrost map over the QTP (Figure 3b, hereafter referred to as the "QTP06_map"). Based on the two published permafrost maps (1:3,000,000 Map of Permafrost on the QTP and the permafrost map of the QTP derived from 1:4,000,000 Map of the Glaciers, Frozen Ground and Deserts in China), the permafrost agreeing in both maps is defined "both permafrost", the seasonally frozen ground agreeing in both maps as "both seasonal frozen ground", and the rest as uncertain class (Figure 4). The relationship [55] between lower limit and latitude of permafrost in the northern hemisphere shows that the lower limit of permafrost reaches the maximum value at 25°22′N, and then decreases with the increasing of latitude. In view of latitude affecting the distribution of permafrost, the QTP is equally divided into five zones from north to south according to the latitude. Subsequent classification is completed on each zone.
The elevation of all the pixels on each latitude zone with both permafrost or both seasonally frozen ground is counted, taking the first zone as an example: if 95% elevation value of both permafrost is greater than 3594 m, then the probability that those corresponding pixels will be Based on the two published permafrost maps (1:3,000,000 Map of Permafrost on the QTP and the permafrost map of the QTP derived from 1:4,000,000 Map of the Glaciers, Frozen Ground and Deserts in China), the permafrost agreeing in both maps is defined "both permafrost", the seasonally frozen ground agreeing in both maps as "both seasonal frozen ground", and the rest as uncertain class ( Figure 4). Based on the two published permafrost maps (1:3,000,000 Map of Permafrost on the QTP and the permafrost map of the QTP derived from 1:4,000,000 Map of the Glaciers, Frozen Ground and Deserts in China), the permafrost agreeing in both maps is defined "both permafrost", the seasonally frozen ground agreeing in both maps as "both seasonal frozen ground", and the rest as uncertain class (Figure 4). The relationship [55] between lower limit and latitude of permafrost in the northern hemisphere shows that the lower limit of permafrost reaches the maximum value at 25°22′N, and then decreases with the increasing of latitude. In view of latitude affecting the distribution of permafrost, the QTP is equally divided into five zones from north to south according to the latitude. Subsequent classification is completed on each zone.
The elevation of all the pixels on each latitude zone with both permafrost or both seasonally frozen ground is counted, taking the first zone as an example: if 95% elevation value of both permafrost is greater than 3594 m, then the probability that those corresponding pixels will be with the increasing of latitude. In view of latitude affecting the distribution of permafrost, the QTP is equally divided into five zones from north to south according to the latitude. Subsequent classification is completed on each zone.
The elevation of all the pixels on each latitude zone with both permafrost or both seasonally frozen ground is counted, taking the first zone as an example: if 95% elevation value of both permafrost is greater than 3594 m, then the probability that those corresponding pixels will be seasonally frozen ground is not greater than 5%, so the pixels which elevation is lower than 3594 m are considered as seasonally frozen ground; likewise, if 95% elevation value of both seasonally frozen ground is less than 3938 m, then the probability that those corresponding pixels will be permafrost is not greater than 5%, so the pixels which elevation is greater than 3938 m are considered as permafrost; the pixels which elevation is between 3594 m and 3938 m are still uncertain class. The elevation threshold values (Table 1) of other zones are similar to the first zone. The LST threshold method for each latitude zone is similar to elevation threshold method. For example, in the first zone, if 95% LST value of both permafrost is less than −2.01 • C, then the probability that those corresponding pixels will be seasonally frozen ground is not greater than 5%, so the pixels which LST value is greater than −2.01 • C are considered as seasonally frozen ground; if 95% LST value of both seasonally frozen ground is greater than −4.56 • C, then the probability that those corresponding pixels will be permafrost is not greater than 5%, so the pixels which LST value is less than −4.56 • C are considered as permafrost; the pixels which LST value is between −2.01 • C and −4.56 • C are still uncertain class. The LST threshold values (Table 1) of other zones are similar to the first zone.
Although many literatures pointed out that vegetation can be used as an indicator of frozen soil types [44], the permafrost or seasonally frozen ground surface vegetation types and coverage is not uniform. For example, the vegetation features of permafrost or seasonally frozen ground for each latitude zone are different, and even there are both permafrost and seasonally frozen soil in a same range of NDVI of a latitude zone. Because vegetation as a type indicator of frozen soil is complicated and diverse, NDVI threshold value is different with elevation and LST. NDVI value range of both permafrost and seasonally frozen ground in each latitude zone is counted as Table 2. Obviously, there is a partial intersection in the NDVI range of values, taking the first latitude zone as an example, the NDVI of both permafrost is 0.054-0.133 and the NDVI of both seasonally frozen ground is 0.026-0.219, which means the intersection interval (0.026-0.133) has both permafrost and seasonally frozen ground, so this intersection part is still divided into uncertain class in order to maintain mapping accuracy. After re-division of intersection interval in other latitude zones, NDVI classified threshold is shown in Table 1.
Compared with the elevation, LST and vegetation three indicators, AMSR-E soil moisture with 25 km resolution can only be used as auxiliary classification parameter. The method of obtaining the soil moisture threshold value is the same as elevation threshold method.

Classification Decision Tree
In this paper, a classification decision tree ( Figure 5) was established based on the correlation between elevation, LST, NDVI and soil moisture and permafrost distribution. The numbers (a)-(e) represent five internal branches, which also represent the classification decision tree of the first, second, third, fourth and fifth latitude zone respectively. represent five internal branches, which also represent the classification decision tree of the first, second, third, fourth and fifth latitude zone respectively.

Post Classification Processing
After the gradual classification by four indicators (elevation, LST, NDVI and soil moisture) based on decision tree method, there are still few pixels (uncertain class) (Figure 6d) that cannot be classified and these pixels are in the transition zone between permafrost and seasonally frozen ground.

Post Classification Processing
After the gradual classification by four indicators (elevation, LST, NDVI and soil moisture) based on decision tree method, there are still few pixels (uncertain class) (Figure 6d) that cannot be classified and these pixels are in the transition zone between permafrost and seasonally frozen ground. The First Law of Geography, according to Waldo Tobler [56], reads "everything is related to everything else, but near things are more related than distant things." This theory is also applicable to permafrost, if the number of permafrost pixels (Cp) adjacent to an uncertain pixels is greater than the number of seasonally frozen ground pixels (Cs), then those uncertain pixels are classified permafrost, otherwise-as seasonally frozen ground. If the number of permafrost pixels adjacent to an uncertain pixel is equal to the number of seasonally frozen soil pixels, for example, this case occurs in the 3 × 3 range, then we can extend the neighborhood to the 5 × 5 range and further compare until all the uncertain pixels can be classified (Figure 7). Finally, the permafrost distribution map on the QTP based on multi-source remote sensing data is obtained and shown in Figure 8. The First Law of Geography, according to Waldo Tobler [56], reads "everything is related to everything else, but near things are more related than distant things." This theory is also applicable to permafrost, if the number of permafrost pixels (Cp) adjacent to an uncertain pixels is greater than the number of seasonally frozen ground pixels (Cs), then those uncertain pixels are classified permafrost, otherwise-as seasonally frozen ground. If the number of permafrost pixels adjacent to an uncertain pixel is equal to the number of seasonally frozen soil pixels, for example, this case occurs in the the number of seasonally frozen ground pixels (Cs), then those uncertain pixels are classified permafrost, otherwise-as seasonally frozen ground. If the number of permafrost pixels adjacent to an uncertain pixel is equal to the number of seasonally frozen soil pixels, for example, this case occurs in the 3 × 3 range, then we can extend the neighborhood to the 5 × 5 range and further compare until all the uncertain pixels can be classified (Figure 7). Finally, the permafrost distribution map on the QTP based on multi-source remote sensing data is obtained and shown in Figure 8.

Accuracy Evaluation Method
The overall accuracy and Kappa coefficients are usually used as indicators of qualitative assessment accuracy. The overall accuracy [57] is the proportion of the number of pixels correctly classified to the total number of sampled pixels (1), which shows the consistent degree between the classification results of each random sample and the actual conditions of the corresponding region. Kappa coefficient (2) [58] called KHAT statistics, is used to measure the agreement between two maps. In comparison, the overall accuracy indicator just considers the correctly classified pixels, but kappa coefficient (Ka) takes into account the inconsistent classified pixels.
Overall accuracy = ∑ (1) where n is the number of total pixels, nkk is the correctly classified pixels of the k type frozen soil, nk+ is the total pixels of the k type frozen soil in investigate regions, n+k is the total pixels of the k type frozen soil in the corresponding regions of map (k = 1 represent permafrost, k = 2 represent seasonally frozen soil in the above definition). Cohen [59] has been proposed quality values for K as follow: Ka ≥ 0.8 signifies excellent agreement, 0.6 ≤ Ka < 0.8 represents substantial agreement, 0.4 ≤ Ka < 0.6 represents moderate agreement, 0.2 ≤ Ka < 0.4 represents fair agreement, and a lack of agreement corresponds to Ka < 0.2.

Mapping Result Based on Decision Tree
From the process of classification by decision tree, the area of the certain classes (permafrost

Accuracy Evaluation Method
The overall accuracy and Kappa coefficients are usually used as indicators of qualitative assessment accuracy. The overall accuracy [57] is the proportion of the number of pixels correctly classified to the total number of sampled pixels (1), which shows the consistent degree between the classification results of each random sample and the actual conditions of the corresponding region. Kappa coefficient (2) [58] called KHAT statistics, is used to measure the agreement between two maps. In comparison, the overall accuracy indicator just considers the correctly classified pixels, but kappa coefficient (Ka) takes into account the inconsistent classified pixels.
Overall accuracy = ∑ 2 k=1 n kk n (1) Kappa coefficient = ∑ 2 k=1 n kk − ∑ 2 k=1 n k+ n +k n 2 − ∑ 2 k=1 n k+ n +k where n is the number of total pixels, n kk is the correctly classified pixels of the k type frozen soil, n k+ is the total pixels of the k type frozen soil in investigate regions, n +k is the total pixels of the k type frozen soil in the corresponding regions of map (k = 1 represent permafrost, k = 2 represent seasonally frozen soil in the above definition). Cohen [59] has been proposed quality values for K as follow: Ka ≥ 0.8 signifies excellent agreement, 0.6 ≤ Ka < 0.8 represents substantial agreement, 0.4 ≤ Ka < 0.6 represents moderate agreement, 0.2 ≤ Ka < 0.4 represents fair agreement, and a lack of agreement corresponds to Ka < 0.2.

Mapping Result Based on Decision Tree
From the process of classification by decision tree, the area of the certain classes (permafrost and seasonally frozen ground) is gradually expanded, and the number of unidentified pixels is reduced during the execution of the decision tree. Figure 6 shows the QTP permafrost map changes after the gradual classification of four indicators (elevation, LST, NDVI and soil moisture). Table 3 shows the percentage of frozen ground types to the total area of the corresponding classification map. The results directly reflect elevation and LST are most important factors for the permafrost presence/absence on the QTP, and soil moisture is just an auxiliary index.

Validation at Investigated Regions
In this study, the accuracy of permafrost map over the QTP based on multi-source remote sensing data is evaluated by overall accuracy and Kappa coefficients (Table 4) using the permafrost map of the three investigated regions (WKL, WQ and GZ). From the comparative analysis ( Figure 9) and results of verification, with the improving of the resolution or scale, the total precision of the QTP map has been exactly improved. Especially, the WKL region on QTP06_map was completely divided into permafrost and Kappa coefficient was 0 because the permafrost area of this frozen map was overestimated, which was consistent with the other verification conclusions [60]. However, seasonally frozen soil over the WKL of the QTP permafrost map based on multi-source remote sensing data accounts for 9.5% of the total area, and its spatial distribution of permafrost is in agreement with investigated region. Besides, overall accuracy of WQ and GZ on the QTP permafrost map based on multi-source remote sensing data is more than 85%, and Kappa coefficients of WQ and GZ are 0.71 and 0.69 respectively. However, the Kappa coefficients of WQ and GZ of two previous published maps are both less than 0.5. Compared with the QTP96_map and the QTP06_map, the accuracy of the QTP permafrost map based on multi-source remote sensing data has shown better result. sensing data accounts for 9.5% of the total area, and its spatial distribution of permafrost is in agreement with investigated region. Besides, overall accuracy of WQ and GZ on the QTP permafrost map based on multi-source remote sensing data is more than 85%, and Kappa coefficients of WQ and GZ are 0.71 and 0.69 respectively. However, the Kappa coefficients of WQ and GZ of two previous published maps are both less than 0.5. Compared with the QTP96_map and the QTP06_map, the accuracy of the QTP permafrost map based on multi-source remote sensing data has shown better result.

Discussion
According to the verification result over investigated regions and permafrost area comparison with other permafrost maps of the QTP, permafrost map over the QTP based on Multi-Source remote sensing data is better than two published permafrost maps of QTP. It also demonstrates the macroscopic, dynamic and convenient advantages of remote sensing data in permafrost research.
So far, a number of models (elevation model, Mean Annual Ground Temperature Model (MAGT) [54], Emissivity Model [61], and Annual Average Ground Surface Temperature Model [62]) have been developed and they can basically reflect the permafrost distribution on the QTP. However, they still have several uncertainties, such as: (1) The input parameter of elevation model is not enough and the resolution is coarser, which affect the accuracy of the result; (2) The MAGT is short of consecutive ground temperature observations, so there is much error in Southeast Tibet permafrost regions, North Gangdise and near Himalaya Range; (3) The inversion accuracy of the emissivity ratio in Emissivity Model is limited because of the influence of elevation and MAST Model cannot reflect the distribution of frozen soil in south and southeast Tibet due to scarcity of ground surface temperature observation data.
Besides, although the Frost Number [63] and Temperature at the Top of Permafrost (TTOP) [64] have been transferred to the permafrost distribution on the QTP [65], the simulation results were not satisfactory. On the one hand, high-altitude features of the QTP was not considered because frost number was more suitable for high latitude permafrost distribution, on the other hand, topographic conditions of QTP are so complex that it affected accuracy. However, compared with the previous permafrost distribution models of the QTP, permafrost map based on multi-source remote sensing data does not use single factor like the Elevation Model or the Surface Frost Number Model, but comprehensively considers the governing factors (latitude, elevation, LST, NDVI and soil moisture) related to permafrost distribution. Then, on the local scale, the gradual classification is achieved by the approach of decision tree based on the importance of each factor. Furthermore, the permafrost map of the QTP derived from 1:4,000,000 map of the Glaciers, Frozen Ground and Deserts in China compiled in 2002 and the Map of Permafrost on the Qinghai-Tibetan Plateau compiled in 1996, but Aqua satellite was launched on 4 May 2002. In order to reduce time gap with the prior information (two previous published maps), and keep data integrity and accuracy, this study selects multi-source remote sensing data of 2003 year.
When it comes to a series of published permafrost maps of the QTP simulated using different methods, the results of 1:3,000,000 map of permafrost distribution on the QTP by Li and Cheng (129.8 × 10 4 km 2 of permafrost, 122.4 × 10 4 km 2 of seasonally frozen ground) has been internationally recognized. In fact, the permafrost area of this map and the permafrost map over the QTP derived from the 1:4,000,000 map of snow, ice, and frozen ground in China (hereafter referred to as the "QTP88_map") are overestimated [36,60,66]. The QTP06_map (the frozen soil map of QTP based on MAGT) derived from 1:4,000,000 Map of the Glaciers, Frozen Ground and Deserts in China also was used benchmark (Wang, Rinke, et al., 2016), which permafrost area was 111.8 × 10 4 km 2 . However, the permafrost map over the QTP based on multi-source remote sensing data (hereafter referred to as the "MRSD_QTP map") in this study shows that the permafrost area and seasonally frozen ground area are 111.3 × 10 4 km 2 and 140.9 × 10 4 km 2 respectively, its permafrost area is most similar to the QTP06_map and the new map of the permafrost distribution on the Tibetan Plateau (hereafter referred to as the "QTP16_map") provided by Zou et al. [16]; the overall accuracy are 84.82%, 87.52% respectively, and kappa coefficients are 0.70 and 0.74 (Table 5). In addition, the permafrost area from the permafrost stability type distribution map over the QTP in the 2000s provided by Ran et al. (2017) is approximately 107.9 × 10 4 km 2 when the extremely unstable type of permafrost is neglected because this kind of permafrost mainly belongs to frozen gravel and cave ice that are distributed below the lower limit of permafrost and usually is not added to the total area of permafrost [55], this result is also closer to the MRSD_QTP map. The consistency between the two maps is 87.52%, and the Kappa coefficient is approximately 0.74, which illustrates this study generates a better result. However, this study has rather much difference with other frozen soil maps over the QTP (Table 6), such as Regional statistical survey, HY/CLM4, elevation model and TTOP, which may be caused by permafrost degradation due to climate warming in different times, method variability and historic climate data uncertainties.  Although the permafrost map over the QTP based on multi-source remote sensing data performed better than the two published maps (QTP96_map and QTP06_map) and had substantial agreement with the investigated regions, it still had some uncertainties as follows: (1) The prior information was relatively imperative, although this study selected the two most representative permafrost maps at different stages as the reference basis, there are few uncertainties due to the own errors of those two maps; (2) The final uncertain pixels after the classification with decision tree method were reclassified by a regional growth search algorithm; although this had certain theoretical basis, it was still necessary to find a more scientific and reasonable approach to classify these pixels as permafrost or seasonally frozen ground; (3) AMSR-E soil moisture products used in this study is 25 km in spatial resolution; the resolution had not been substantially improved though it was resampled to 1 km, thereby soil moisture was used as ancillary data. With the development of microwave remote sensing, if the accuracy and resolution of soil moisture are improved, the accuracy of permafrost map will also be improved accordingly; (4) Compared with other prediction models, especially the physical model, the result of this study cannot reflect the temporal dynamic changes of permafrost distribution.
The accuracy verification results of permafrost maps show that the area of seasonally frozen ground of the investigated regions is more than that shown by the other permafrost maps in general and discrepancies exist in transition zone between permafrost and seasonally frozen soil. On the one hand, snow (thickness, duration, time and cover) solar radiation, topographic factor, climate change and other factors, except governing factors, will also affect permafrost spatial distribution. For example, the temperature is extremely low in the eastern, southern and abdomen mountains of the QTP where snow exists all year round, thick and long snow cover on the permafrost. It has the effect of heat preservation because the thermal conductivity of snow is very poor so that solar radiation during the daytime is hard to get into the soil and the heat is not easy to lose in the nighttime. However, in the high plains, valleys and basins of high-altitude areas on the QTP, snow thickness is thin and duration is short, so snow might plays a cooling role [71]; On the other hand, time difference among the different permafrost maps of the QTP in which climate warming caused permafrost degradation, will lead to the area of permafrost reducing and the area of seasonally frozen ground increasing gradually. Therefore, it is possible to generate a more accurate and reasonable permafrost map of the QTP integrating with passive microwave snow products.

Conclusions
This study indirectly obtained the permafrost map on the QTP by exploiting the relationship between environmental factors and frozen ground rather than limited investigated data, it is a pragmatic approach for regional scale permafrost mapping. Compared to traditional mapping method, spatial data derived from the satellite sensors has the macro, dynamic and comprehensive advantages, and remedy the lack of investigated data.
Statistical result of this study shows that permafrost area is 42.5 percent of the QTP area (111.3 × 10 4 km 2 ), and the seasonally frozen ground area is 53.8 percent of the QTP area (140.9 × 10 4 km 2 ), which has the highly consistent with benchmark map (The Map of Permafrost on the Qinghai-Tibetan Plateau). Moreover, the validation result shows that permafrost map of QTP based on multi-source spatial data has high overall accuracy and kappa coefficient with three regions (West Kunlun, Wenquan and Gaize) compared other two existing permafrost maps.
Remote sensing techniques have considerable potential for use as a tool in regional scale permafrost mapping, which might provide a simple and promising approach for further permafrost research. Furthermore, we believe that combined with passive microwave snow products, soil types, energy exchange process model between ground and atmosphere, terrain factors and physical mechanism of permafrost, the mapping model would be improved and might yield still more accurate results.