Tracking Reforestation in the Loess Plateau, China after the “Grain for Green” Project through Integrating PALSAR and Landsat Imagery

Abstract

An unprecedented reforestation process happened in the Loess Plateau, China due to the ecological restoration project ‘Grain for Green Project’, which has affected regional carbon and water cycles as well as brought climate feedbacks. Accurately mapping the area and spatial distribution of emerged forests in the Loess Plateau over time is essential for forest management but a very challenging task. Here we investigated the changes of forests in the Loess Plateau after the forest reconstruction project. First, we used a pixel and rule-based algorithm to identify and map the annual forests from 2007 to 2017 in the Loess Plateau by integrating 30 m Landsat data and 25 m resolution PALSAR data in this study. Then, we carried out the accuracy assessment and comparison with several existing forest products. The overall accuracy (OA) and Kappa coefficient of the resultant map, were about 91% and 0.77 in 2010, higher than those of the other forest products (FROM-GLC, GlobeLand30, GLCF-VCF, JAXA, and OU-FDL) with OA ranging from 83.57% to 87.96% and Kappa coefficients from 0.52 to 0.68. Based on the annual forest maps, we found forest area in the Loess Plateau has increased by around 15,000 km² from 2007 to 2017. This study clearly demonstrates the advantages of data fusion between PALSAR and Landsat images for monitoring forest cover dynamics in the Loess Plateau, and the resultant forest maps with lower uncertainty would contribute to the regional forest management.

1. Introduction

The Loess Plateau in northern China has been a fragile zone with severe land degradation issues due to intensive agricultural and grazing activities. In modern times, due to the increasing population, excessive human economic activities, long-term extensive land reclamation have caused large amounts of sediments to be lost every year on the Loess Plateau, and soil erosion was extremely serious. The environment of the Loess Plateau is very fragile with severe soil erosion [1,2] and agricultural production is low and unstable, which to a large degree affects the local economy and human wellbeing [3,4]. An unprecedented ecological restoration project, so-called Grain for Green Project (GGP), was initiated by the Chinese government in 1999. The project returned the marginal croplands to forests, for alleviating soil and water loss, and preventing land degradation through afforestation and reforestation [5]. The project has achieved remarkable effectiveness such as increasing forest cover, enhancing carbon sequestration, controlling soil and water loss, and reducing dust to other regions and countries. However, there have been limited efforts on generating spatially explicit maps of forest changes during the last decade.

Remote sensing has been an essential tool for assessing forest cover dynamics at regional to global scales. Such applications are quite valuable in those regions that are remote with limited human access. Previous studies have generated coarse resolution forest maps primarily from optical sensors like Advanced Very High Resolution Radiometer (AVHRR) [6,7], and Moderate Resolution Imaging Spectroradiometer (MODIS) [8,9,10]. Despite the coarse resolution, those forest maps can adequately describe the general patterns of the spatial distribution of forests and can track the annual dynamics of large intact forests. However, those coarse resolution forest maps have a considerable uncertainty in fragmented landscapes like Loess Plateau [11].

The open access of the archived Landsat imagery enabled the extensive application for forest monitoring [6,12,13,14,15]. These 30-m forest products [6,12] can capture more detailed information in the spatial distribution and dynamic changes of forests than very coarse resolution imagery like AVHRR and MODIS. Hansen et al. [8,9,12] used the Landsat data source for global forest change monitoring. However, there are two main problems in optical remote sensing images. First, the utility of optical remote sensing images is often limited by the quality of observations and restricted by persistent clouds and cloud shadows [16]. Although images from multiple years can increase the quality of observations, it will need high monetary and computational costs. Second, the canopy level characteristics cause the accuracy of the Landsat data to have a considerable commission error because of the similar structural and spectral attributes in shrubs or vegetation with high biomass [17].

Active microwave remote sensing images have the ability for land surfaces mapping and have been widely used to map the forest because of the advantage of all-weather use and have been more applicable especially as wall-to-wall global observations become available. The L-band synthetic aperture radar (SAR) data has a relatively longer band wavelength (23.6 cm), which has greater penetration capacity into forest biomass under the canopy and more sensitive to forest structure interacts with large trunk and branch components [18,19,20,21]. That is particularly suitable for forest monitoring with stronger volume scattering signals than that of X-band sensor and C-band sensor. Therefore, it is preferred to map the forests at regional scale. However, the accuracy of SAR data will be influenced by the underlying soil conditions, moisture content, and high biomass croplands with backscatter coefficients like forests [22,23,24].

To map forests with a high accuracy and low uncertainty is a very challenging task using either optical or SAR remote sensing images. Recently, several studies [2,25,26,27] have focused the integration of both optical and SAR images to map forests. Qin et al. [28] utilized the PALSAR data with 25 m spatial resolution and optical remote sensing data (Landsat) with 30 m spatial resolution to extract the forest information of Sub-Humid and Semi-Arid Regions from Oklahoma, in 2010, which has the highest overall accuracy compared with other forest products, with an overall accuracy of 88.2%. Yang et al. [29] utilized the same method to produce a forest product with 30 m spatial resolution to analyze the agroforests in North China Plain. Considering the above successful cases and complementary advantages between optical and SAR remote sensing data, the integration of SAR and optical remote sensing images may enable improvement of forest mapping [30].

In this study, we aim to analyze the forest cover changes in the Loess Plateau after GGP based on the integration of PALSAR and Landsat imagery. We tried to answer the following questions: (1) Is the approach integrating Landsat and PALSAR data for mapping forests effective and robust in fragile regions like the Loess Plateau? (2) What is the spatiotemporal pattern of forest cover changes after GGP in the Loess Plateau? To achieve this object, firstly, we produced a new forest map using the 25 m PALSAR and 30 m Landsat data in the Loess Plateau located in the north of central China. Secondly, we compared our forest map with the existing forest maps with 30–50 m resolutions. Thirdly, we analyzed the spatiotemporal pattern of forest dynamics in the Loess Plateau by using the new forest map. This study expects to facilitate an assessment on the effectiveness of GGP on environmental restoration, and to offer a valuable reference for future detection of forest changes.

2. Data and Methods

2.1. Study Area

The Loess Plateau is one of the most corrodible regions in the world which is located at latitude of 34°–40°N and longitude of 103°–114°E and is approximately 1000 km long from west to east, and 700 km wide from south to north. The Loess Plateau stretches across a region north to the Yinshan Mountains in Inner Mongolia, south to the Qinling Mountains, west to the Riyue Mountains, east to the Taihang Mountains which covers a total area of approximately 6.4 × 10⁵ km², 6.5% of the territory in China, with average elevations ranging from 1500 to 2000 m above sea level. This region covers seven provinces, including Gansu, Henan, Ningxia, Inner Mongolia, Qinghai, Shanxi, and Shaanxi (Figure 1). It has complex topography which is full of hills, basins, and alluvial plains with continuous loess of 100–300 m in thickness. It has the continental monsoonal climate which lies in a semi-humid and semi-arid transitional zone with an annual average temperature of about $6$ to $14$°C and annual precipitation of about 200 to 700 mm [31]. Approximately 70% of annual precipitation occurs between June and September during summer. The Loess Plateau is one of the most severely eroded regions globally because of frequent storms, highly erodible soil, low vegetation cover, and steep landscape [32,33]. The major regions of the Loess Plateau are undergoing heavy soil erosion because of human disturbances over a long period of time. Therefore, forests mainly exist in some mountainous areas of the plateau.

Wuqi County was the first place to carry out the GGP, which is located in the northwest of Yan’an County, Shaanxi Province. In 1998, the county took the lead in launching the project. In early 2000, it was listed as a pilot site for forest and grass vegetation restoration on the Loess Plateau of the State Forestry Administration, and implemented conversion from farmland to forests. The county is located in the transitional zone between agriculture and pasture in the southern margin of the Mausu Sandy Land. It covers the middle reaches of the Yellow River, the Beiluo River, and the upper reaches of the Wuding River. It belongs to the Liangzhugou area of the Loess Plateau. The geographical position is 107°38′57′′~108°32′49′′ east longitude, 36°33′33′′~37°24′27′′ north latitude, the total area is 3791 km², and the altitude is 1233~1809 m.

2.2. Data Preprocessing

2.2.1. PALSAR Dataset and Preprocessing

The 25 m ALOS PALSAR Fine Beam Dual polarization (FBD) products from 2007 to 2017 generated with the selected image acquired between the months of June and October [22] can be downloaded from the Earth Observation Research Center, Japan Aerospace Exploration (JAXA) (http://www.eorc.jaxa.jp/ALOS/en/palsarfnf/data/). The PALSAR (Table 1) 25 m orthorectified mosaic data, provided by JAXA, is created globally using images in the ascending path. The original PALSAR images have an off-nadir angle of 34.3° and cover wet and dry seasons, or winter and summer seasons, and were resampled into the 25 × 25 m mosaic to create one composite per year. Polarized signals can be transmitted and received in the horizontal or vertical dimensions by the L-band PALSAR, allowing for horizontal transmission and reception (HH) and horizontal transmission and vertical reception (HV). While the ALOS has launched, the mosaic of HH and HV polarization data, local incidence, and the mask information (effective area, layover, ocean flag, shadowing, and void area) were included in the datasets. PALSAR HH and HV backscatter data were slope corrected and ortho-rectified with a geometric accuracy of about 12 m to guarantee the geometric accuracy, using the 90 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) [34], before being radiometrically calibrated. The Digital Number (DN) values (amplitude values) were converted into gamma-naught backscattering coefficients in decibels (γ°) using a calibration coefficient from JAXA [35,36]

(1)${\mathsf{\gamma }°=10\times \log }_{10}{<\mathrm{DN}}^2>+\mathrm{CF},$

where $\mathsf{\gamma }°$ is the backscattering coefficient, DN is the digital number value in HH or HV, and CF is the absolute calibration factor of $-83$ [36]. Further, the Ratio and Difference values were calculated using the resultant HH and HV backscattering coefficient in decibels. In which, HH and HV are raster data with 25 m resolution because of their effective classification of the land cover [25,37,38]. The calculation of backscatter coefficients and generation of ratio and difference bands was conducted following the below equations.

(2)$\mathrm{Difference}=\mathrm{HH}-\mathrm{HV},$

(3)$\mathrm{Ratio}=\mathrm{HH}/\mathrm{HV}.$

We also resampled the 25 m PALSAR images (HH, HV, Difference and Ration layers) to 30 m PALSAR images using the neighbor sampling method to match Landsat images at 30 m spatial resolution using the nearest neighbor method. The 30 m PALSAR resultant images were then combined with Landsat images to generate a PALSAR/Landsat forest map.

2.2.2. Landsat Images and Preprocessing

We processed all available surface reflectance products of Landsat (Table 1) accessible from the Google Earth Engine (GEE) which is an open and powerful platform for using Landsat image processing (https://developers.google.com/earth-engine/datasets/catalog) such as atmospheric correction. The international mainstream satellite remote sensing data, including Landsat, Sentine and MODIS, was stored on the Google Earth Engine Cloud Platform. The Google Earth Engine can store, manage, and process massive amounts of remote sensing data, providing a technical means for global or regional land use classification. Some built-up lands, high biomass corn fields, and barren lands which have a complex structure and rough land surface that may bring about high PALSAR backscattering coefficients similar to the forests will produce the commission error. Therefore, we generated an annual maximum NDVI (NDVImax, hereinafter) map to exclude the potential disturbance from buildings which also have backscattering coefficients [26,39]. Poor quality of observations which were caused by the shadows and clouds were removed based on CFmask [40,41] and metadata because the CFmask algorithm can provide cloud, cloud shadow, and snow mask for Landsat TM/ETM+ images. Pixels that are cloudless and not in ETM + scan line zoning were selected as good quality observations.

(4)$\mathrm{NDVI}=(\mathrm{NIR}-\mathrm{RED})/(\mathrm{NIR}+\mathrm{RED}),$

(5)$\mathrm{NDVImax}=\mathrm{Max}(\mathrm{NDVI}1,\mathrm{NDVI}2,\mathrm{NDVI}3,\dots \mathrm{NDVIn}),$

where NIR and RED are the land surface reflectance values of red (RED, 630–690 nm) and near-infrared (NIR, 760–900 nm) spectral bands of Landsat TM/ETM+/OLI image, respectively.

2.3. Mapping Algorithms

The workflow of using PALSAR and Landsat TM/ETM+ images to map the forests is illustrated in Figure 2. We produced the resultant PALSAR/Landsat forest map in two steps. Firstly, we generated the PALSAR-based forest map using the HV, HH, difference, and ratio of the PALSAR data. Secondly, we generated PALSAR/Landsat forest maps with high accuracy and low uncertainties using the integration of PALSAR and Landsat. In summary, forest was first identified using structure-based metrics from PALSAR images, then reduced commission errors using time series of Landsat.

According to the United Nations Food and Agriculture Organization (UN FAO), the forest definition includes both forest structure and canopy information using a unit of land (>0.5 ha) with tree cover > 10% and the minimum height > 5 m. The previous studies found that the threshold-based approach for forest mapping is robust and extendable in different regions, such as the monsoon region in Asia, Oklahoma in USA [27,28], and the North China Plain, in China [29].

A detailed workflow was developed to map the forest and multi-source forest dataset comparison in Figure 2. We used the L-band PALSAR data which can penetrate into forests with substantial volume scattering through the incident energy interaction with trunks and branch components to capture the structure and above ground grass (AGB) of forests. The previous studies show that forest and forest AGB exhibit a certain range of PALSAR backscattering coefficients, respectively [25]. In this study, we tried to use the same decision classification algorithm used in our previous study in the North China Plain, in China [29]. Therefore, we used the HH, HV (Table 1), the difference ($\mathrm{HH}-\mathrm{HV}$), and the ratio (HH/HV) to map the initial PALSAR forest baseline map.

The threshold values were determined based on the workflow addressed in Figure 2 and the final threshold values (Figure 3) were obtained by signature analyses of training samples. In total, we collected 2489 ROIs which include 660 forest ROIs, 410 cropland ROIs, 253 water ROIs, and 1160 built-up land and other ROIs to determine the thresholds in this study. The histogram of the HH image (Figure 3a) shows that water has much lower backscatter than forest and partly overlapped with cropland, as water bodies tend to be relatively calm and smooth, and can reflect most of the backscatter through specular reflection. This means water is separable from forests, built-up lands, and most cropland. Both forests and built-up lands have higher HH backscatter values and some built-up lands such as large building complexes have even higher backscatter values than forests due to the strong reflectance environment caused by their complex structure and corner reflectors. Forests, cropland, and built-up lands also have some overlap of HH backscatter values. The histogram of the HV image (Figure 3b) shows that forests have higher HV values than water and cropland, and have relatively more overlap with built-up lands due to large patches of forests distributed in urban areas. Therefore, HV is an additional and effective indicator to separate forests from water and cropland, but not sufficient for complete separation of forests and built-up lands. The histogram of the difference image (Figure 3c) shows that forests have low Difference values, while water has high Difference values. The histogram of the Ratio image (Figure 3d) shows that forests overlap with the other three land cover types. Based on the 95% confidence intervals of the HH, HV, Ratio, and Difference images, the threshold values for the decision tree algorithm were determined. As some uncertainty may exist in the ROIs of different land cover types, we excluded 2.5% pixels with the lowest and highest values. The threshold values were further rounded to integer numbers for HH and HV images, and to 0.5 for Ratio and Difference images.

The above exploratory data analysis provides useful information for constructing a decision tree algorithm based on the threshold values from the HH, HV, $\mathrm{HH}-\mathrm{HV}$, and HH/HV images. First, water which reflects most of the backscatter through specular reflection can be easily identified because of low values in both HH and HV. Second, forests have high values in both HH and HV, and low values in Difference ($\mathrm{HH}-\mathrm{HV}$) because the L-band of PALSAR has great penetration into forests causing huge volume scattering when its incident energy interacts with large trunks and branch components. Third, the most croplands can be easily identified because the Difference (HH − HV) indicates that the forests have the low backscatter values, while the croplands have relatively high backscatter values. The thresholds of the PALSAR-based forest mapping algorithm were determined as −16 < HV < −8 and 2 < Difference < 8 and 0.3 < Ratio < 0.85.

The PALSAR-based forest maps (Figure 4a) often have lower NDVI causing a commission error due to some natural features like rock lands or manmade structures such as urban, buildings. The NDVI_max values of forests which are more than 0.7 are greater than that of rocky lands and urban lands. The PALSAR data above can use Landsat NDVI_max (>0.7) to reduce commission error of forest maps. Finally, in this study, we combined PALSAR and Landsat NDVI_max (Figure 4b) to generate annual maps of PALSAR/Landsat forest maps using the thresholds: −16 < HV < −8, 2 < Difference < 8, 0.3 < Ratio < 0.85, and NDVI_max > 0.7.

2.4. Validation Samples for Accuracy Assessment of PALSAR/Landsat Forest Maps

The result validation is vital to land cover classification studies. Accurate surface coverage samples are fundamental for ensuring the reliability of the result and effectiveness of classification algorithm. The validation data used in this study was from the high resolution satellite images taken in 2010 from Google Earth which has been widely used in validation of land cover products because of its high spatial accuracy [42,43].

In this study, we used the stratified random sampling method to acquire the validation samples. According to the FAO forest definition above [44], we attained validation samples from the very high resolution (VHR) remote sensing images of Google Earth circa 2010 to guarantee equal probabilities of obtaining the validation samples. To collect samples, we performed the following steps. First, the sampling was stratified by class (forest layer and non-forest layer) based on the previous PALSAR-based [45], and then the random points were generated by using the ArcGIS software (Figure 5a). Second, according to these random points, we established circles with a radius of 25 m as Regions of Interest (ROIs) which can guarantee the homogeneity of landscape and the accuracy of validation result [46,47]. If the size of ROI is too small, the location error inherent in satellite images will affect the verification results and create more uncertainty in landscape heterogeneity. At the same time, we abandoned the ambiguous ROIs to ensure the quality of ROI. Finally, the ROIs were put into the Google Earth and interpreted according to the very high resolution images from around 2010 (Figure 5b). Furthermore, we used the Collect Earth (Figure 5c), which is a tool integrating very high resolution imagery from multiple sources (e.g., Google Earth, Being Maps) and also customized imagery in Google Earth Engine. In total, we acquired 2898 ROIs which include 683 forest ROIs and 2215 non-forest ROIs to validate the forest products in this study.

2.5. Accuracy Assessment and Comparison with Other Forest Maps

To analyze the consistency and uncertainty of the available forest maps at relatively high spatial resolution in the Loess Plateau, we compared five popular existing moderate resolution forest maps in 2010 (Table 2): (1) FROM-GLG Global Land Cover Data, which was generated by Tsinghua University [48]. (2) GlobalLand30 product, which was generated by the National Geomatics Center of China (NGCC) [49]. (3) The University of Maryland’s Global Land Cover Facility Data Center (GLCF) produced a 30 m resolution continuous tree cover product [50] called the Vegetation Continuous Fields (VCF) dataset which was derived from all seven bands of Landsat 5 and 7 images based on rescaling the 250 m MODIS VCF Tree Cover layer(http://glcf.umd.edu/data/landsatTreecover/). (4) The JAXA forest map generated by Earth Observation Research Center, JAXA, using L-band PALSAR data. (5) Through the integration of PALSAR and MODIS/Landsat data, the University of Oklahoma (OU) [26] generated an updated forest data layer (FDL), referred to as OU-FDL in this paper.

The mutual ROIs were used to evaluate the accuracy of the PL-based and the other five forest maps in the Loess Plateau, using the confusion matrix method for the estimates of producer accuracy (PA), user accuracy (UA), and overall accuracy (OA). Then, all the ROIs were divided into two parts to evaluate the accuracy of each forest map in the forest region and non-forest region. Note that we also adjusted the accuracy assessment based on the area of each stratum following previous studies [51].

3. Results

3.1. Accuracy of the PL-Based and the Other Five Forest Maps

The spatial distribution of the PL-based and the other five forest maps are shown in Figure 6. The overall accuracy and the Kappa coefficient of the PL-based forest map are 91.27% and 0.77, respectively (Table 3). The user accuracy and producer accuracy are 91.07% and 76.41%, respectively, in forest area. In the non-forest areas, the user accuracy and producer accuracy are 91.33% and 97.07%, respectively.

In addition, the forest map of FROM-GLC showed the accuracy second to that of the PL forest map in the Loess Plateau, with an overall accuracy of 87.96% and Kappa coefficient of 0.68, and the user accuracy and producer accuracy of forest areas are 81.26% and 71.52%, respectively. In non-forest areas, the user accuracy and producer accuracy are 90.02% and 93.40%, respectively. In non-forest areas, there is a large commission error. The JAXA forest map relies on the PALSAR data. The overall accuracy and Kappa coefficient are 83.57% and 0.5236, respectively. The user accuracy and producer accuracy in forest areas are 59.01% and 67.28%, while those in non-forest areas are 91.15% and 87.82%, respectively. JAXA forest map only approximates the producer accuracy of the PL forest map in the non-forest areas, and other precisions are obviously lower than the PL forest map. The OU-FDL forest map integrated PALSAR and MODIS data sources. The overall accuracy and Kappa coefficient are 84.40% and 0.5622, respectively. The user accuracy and producer accuracy in forest areas are 65.30% and 67.47%, while those in non-forest areas are 90.29% and 89.41%. Compared with JAXA, the accuracy in the forest areas has been improved.

On the pixel scale, these forest maps have a good consistency in the areas with non-forest coverage and greater uncertainties in the areas with forest coverage (Figure 7). Throughout the entire state area, about 1.5% and 80.6% of the pixels were identified by these forest maps as consistent forest and non-forest.

Specifically, the consistent forest between PALSAR/Landsat forest map and FROM-GLC forest map, JAXA forest map and OU-FDL forest map accounted for 10.48%, 7.33%, and 9.22% of the total pixels respectively, as well as 80.87%, 84.90%, and 84.37% of the consistent non-forest pixels (Figure 8). Compared to the PALSAR/Landsat forest map, JAXA forest map clearly underestimated forests in the Loess Plateau and overestimated the non-forest areas while the FROM-GLC forest map underestimated the forest areas in the Loess Plateau.

Compared with the PALSAR/Landsat forest map, the JAXA forest map obviously underestimated the fragile forests in the Loess Plateau, and there were obvious omission errors in the forest areas (Figure 9a,d,g,j). The OU-FDL forest map improved the accuracy of forest area on the basis of JAXA forest map, but the accuracy of non-forest area decreased (Figure 9b,e,h,k). The FROM-GLC forest map underestimated the non-forest areas, resulting in high commission errors (Figure 9c,f,i,l).

By comparing the accuracy of PL forest map with that of the other five forest maps, the PL forest map has a higher accuracy and consistency which lay a good foundation for the analysis of forest dynamic change monitoring in the Loess Plateau.

Figure 10 showcases the forest changes on the Loess Plateau between 2007 and 2017. First, the grey part is the unchanged non-forest areas (Figure 10a), specifically, the non-forest areas in 2007 (Figure 10b) are still in the non-forest areas in 2017 (Figure 10c). Second, the green part means unchanged forest areas (Figure 10d) from 2007 (Figure 10e) to 2017 (Figure 10f). Third, the blue part means the conversion from the non-forest areas to the forest areas where the main change areas are in the black and red circles (Figure 10g), the two major areas have changed from forests in 2007 (Figure 10h) to urban lands in 2017 (Figure 10i). Fourth, the red parts are related to the conversion from the non-forest areas to the forest areas where the main change areas are concentrated in yellow, red, and black circles (Figure 10j), these three major areas have changed from desert areas in 2007 (Figure 10k) to the forests (Figure 10l) in 2017.

3.2. Changes in Area of Forest Cover in LP from 2007 to 2017

According to the PL forest maps, the area of forest covered region increased from 1.05 × 10⁵ km² to 1.20 × 10⁵ km² from 2007 to 2017 (Figure 11), with an increase by 0.15 × 10⁵ km² (~14.3% increase). The annual average forest area was about 1500 km². However, between 2007 and 2010, the area of forests declined from 1.05 × 10⁵ to 0.91 × 10⁵ km², with a decrease of 0.14 × 10⁵ km²; after 2010 the area of forest increased, which may have been caused by different stages of GGP implementation. The first phase of GGP project is from 1999 to 2010, and the second phase started after 2010.

The overall forest area shows a growing trend during the period of the GGP. The trend can be described by the equation Area = 0.06Year + 0.87 (Figure 11). According to this equation, the forest area of the Loess Plateau from 2007 to 2017 grew at an average rate of 600 km² per year.

3.3. Changes in Spatial Pattern of the Forest Cover in LP from 2007 to 2017

The spatial distribution of forest cover in the Loess Plateau from 2007 to 2017 is shown in Figure 12. We can see that the forest and non-forest areas have changed from 2007 to 2017, specifically, the forest area grew while the non-forest area decreased.

The land cover changes can be divided into four types as shown in Figure 13: unchanged non-forest areas in grey, non-forest areas to forest areas in red, forest areas to non-forest areas in blue, and unchanged forest areas in green (Figure 13a). The unchanged non-forest areas mainly concentrated in the northwest of the Loess Plateau, including Inner Mongolia, Gansu, Qinghai, Ningxia, the northern regions of Shaanxi and Shanxi which are adjacent to Inner Mongolia (Figure 13b). These areas with high altitudes and slopes did not change much and remained in non-forest areas.

The unchanged forest areas (Figure 13c) were mainly distributed in the southeastern part of the Loess Plateau, concentrated in the central part of Shanxi Province and Shaanxi Province, and partially distributed in the southern edge of the Loess Plateau. The forests in the central region of Shaanxi and Shanxi grew because of the protection policy while the forests in the southern edge of the Loess Plateau remained stable due to the suitable environment and relatively low human disturbances.

The conversion from forest to non-forest areas (Figure 13d) was mainly concentrated in the southern part of the Shaanxi Province, like Hu County, Xingping County, and Liquan County. The conversion from non-forest to forest areas (Figure 13e) was mainly distributed in the central part of Shaanxi Province, such as Zhidan County, Ansai County, Ganquan County, and the central part of Shanxi Province, like Linxian County, Fangshan County, and Lishi County.

3.4. Difference among Different Counties in the Loess Plateau

Here we analyzed the forest changes in the Wuqi County, which was the first county for the GGP. We found that the forest area of Wuqi County increased from 66.35 km² in 2007 (Figure 14a) into 420.10 km² in 2017 (Figure 14b), with an increase of 353.72 km² or 533%. It can be seen that the conversion from non-forest to forest areas was significantly larger than the conversion from forest to non-forest areas in Wuqi County and the afforestation was mainly distributed in the southeastern part of Wuqi County and parts of the northeast (Figure 14c). It is interesting to notice a unique pattern with remarkably greener landscape in Wuqi County comparing to the neighboring counties (Figure 14d), which also suggests the important effects of the GGP for changing the regional landscape and environment.

4. Discussion

4.1. Advantage of the PALSAR/Landsat Forest Map Compared to Other Data Forest Maps

Managing the forest resources in the Loess Plateau is a challenging task because the area and spatial distribution of forest are being affected by a variety of factors, including complex terrain, changing forest markets, unstable climate, and the economic situation of land owners [1,52]. Therefore, forest management and policy formulation [53] urgently require timely and detailed spatial resolution forest maps. Compared with the other forest maps, the PL-based forest map with the higher accuracy and lower uncertainty in this study could discriminate more forests in the Loess Plateau, especially in the fragmented and fragile regions.

As cloud and cloud shadows limit the availability and practicability of data, Landsat imagery is usually selected from many years to increase the number of high-quality observations and to generate land cover maps, including large-scale forest maps [6,15,48]. Although, these maps cannot detect forest changes and may introduce considerable uncertainty in the estimation of the forest fragmentation, forest biomass and carbon cycle [2,53]. The extensive forest cover loss and benefits were founded in the Loess Plateau according to the resultant PALSAR/Landsat forest maps from 2007 to 2017, especially in the fragile regions.

The integration of PALSAR and Landsat images provided complementary information for other datasets. The PALSAR/Landsat forest maps could provide the forest extension and its temporal and spatial changes for policy makers in this region. Meanwhile, the PL-based forest map integrated the advantage of forest spatial structural information and spectral information which can identify more detailed biomass information and help decision makers to better formulate policies.

The JAXA forest map was purely generated based on the ALOS PALSAR fine dual-polarized data which is vulnerable to the terrain and soil water content. Therefore, the product includes commission errors from the buildings in the urban area with high backscatter values. The PALSAR data was collected from June to October, and the seasonal land change from June to October had some influence on the classification accuracy. According to our previous survey [2,27,28], PALSAR data is one of the most suitable remote sensing images for forest mapping. However, its own advantage did not show in the fragile areas of the Loess Plateau which indicates the PALSAR remote sensing image is suitable for grasping the dense forest and large area of forest coverage, and not suitable for the fragile areas like the Loess Plateau. The OU-FDL forest product improved the accuracy by involving the PALSAR and the MODIS to capture more information of forest and reduced the commission errors compared to the JAXA forest map. However, MODIS data has a coarse resolution and the fusion is not much improved due to the mixed pixel issue. The other forest maps are generated only from Landsat data, but different factors are taken into account in the selection of data sources. Specifically, GlobalLand30 forest product uses a static scenario with the least amount of clouds, supplemented by HJ-1 satellite data, so the impact of clouds and cloud shadows still exists. FROM-GLC and GLCF-VCF forest products use multiple Landsat images to eliminate the effect of cloud and cloud shadows.

The PALSAR/Landsat forest maps were generated by the integration of ALOS PALSAR and Landsat images, which contain more information like structure, biomass, and spectral information of forests to reduce some commission errors in the forest maps compared to other PALSAR-based forest maps and Landsat-based forest maps. This study suggests the PL forest map is suitable for forest mapping in fragmented areas like the Loess Plateau.

4.2. Uncertainty of the Study

The definition of forests plays a vital role in the assessment of the forest areas. Tree density and tree height make up the two main aspects of forest definition. However, different forest products have different requirements for these two standards, which could be the source of uncertainty. The main difference in the spatial distribution of forest maps lies in the areas with low forest coverage, which was considered to be the main reason for the definition of forests [1]. In this study, the PALSAR/Landsat forests and the JAXA forests both use the FAO forest definition: forest canopy coverage > 10% and tree heights > 5 m with a minimum area > 0.5 ha, while the GLCF-VCF forest products and the OU-FDL forest map use the tree density (>10%) and the tree height threshold (5 m). The GlobeLand30 forest maps used the looser forest definition (only tree cover 10%, no criteria for tree height). Forest in the FROM-GLC is defined as having a land cover of > 15% forest canopy coverage and tree heights > 3 m. The uncertainty of the forest maps could be also from the varied algorithms. The PL-based forest map was generated based on a robust algorithm which adopted the improved thresholds for the algorithm and the involvement of higher resolution Landsat images. The reliability of the unified threshold has been proved in our previous study [26]. The JAXA forest map [22] used a decision tree algorithm, whose HV threshold values are slightly different from that of the PL-based forest map in this study. The object-based analysis after segmentation process may lead to the loss of sparse forest patches [22]. As a result, the JAXA forest maps often tend to miss those areas with low forest coverage. The GlobeLand30 adopted a pixel-object-knowledge-based (POK-based) classification approach which also leads to the loss of some potential information in the sparse forests. The FROM-GLC used automatic classification algorithms and the GLCF-VCF used the supervised classification. Both algorithms require the collection of forest and non-forest training samples through on-site survey and/or very high-resolution images for algorithm calibration. Therefore, the algorithm was affected by the richness, spatial distribution, and spatial representation of training samples.

Under the GGP, both restoration and afforestation have contributed to the increase of forest area. However, some studies have shown that due to the increased transpiration, afforestation in unsuitable areas can lead to water shortage, decrease of land surface temperature, and deterioration of soil ecosystems [54,55,56]. In order to more accurately understanding the ecological and environmental effects of GGP, more studies are needed on the biophysical mechanism and ecosystem services due to the ecological restoration project and climate change in the future.

5. Conclusions

There are relatively large uncertainties and low temporal accuracy in fragmented areas when estimating forest distribution based on either synthetic aperture radar (SAR) or optical remote sensing images alone. We proposed a straightforward and robust rule-based algorithm to map the forest based on the Advanced Land Observation Satellite (ALOS) Phased Array type L-band Synthetic Aperture Radar (PALSAR) Fine Beam Dual polarization (FBD) mosaic images and the annual maximum Normalized Difference Vegetation Index of land satellite (NDVI_max). The PALSAR/Landsat forest map in 2010 with the relatively high overall accuracy and lower uncertainty showed improvements against selected forest maps. This indicated that the integration of PALSAR and Landsat data is more advantageous in estimating forest area than using the Landsat or PALSAR alone in traditional fragile regions like the Loess Plateau.

We investigated the temporal and spatial variations of forest cover after GGP in the Loess Plateau from 2007 to 2017, and found the forested areas increased about 0.15 × 10⁵ km² from 2007 to 2017 in the Loess Plateau, suggesting the effectiveness and sustainability of ecological restoration projects like GGP. The resultant PALSAR/Landsat forest maps would provide a timely and accurate data source for forest management and policy formulation under the future GPP in fragile areas with complex climates.

Author Contributions

Conceptualization, F.X. and J.D.; Data curation, F.X. and H.Z.; Formal analysis, Z.Y. and H.Z.; Funding acquisition, H.Z. and G.Z.; Investigation, H.Z. and Z.Y.; Methodology, J.D. and J.Z.; Project administration, G.Z. and Y.Q.; Resources, Z.Y. and H.Z.; Software, Z.Y. and H.Z.; Supervision, F.X. and J.D.; Validation, H.Z.; Visualization, Z.Y. and H.Z.; Writing—original draft, H.Z. and Z.Y.; Writing—review and editing, J.D., Y.Q., and X.X.

Funding

This study is funded by the Strategic Priority Research Program (XDA19040301) and Key Research Program of Frontier Sciences (QYZDB-SSW-DQC005) of Chinese Academy of Sciences (CAS), National Natural Science Foundation of China (61772078), Key R&D Program of Beijing (D171100001817003).

Acknowledgments

We thank Yanyan Pei, Jiadi Yin, Rui Zhao, Arifuzzaman Khondakar, and Mrinal Singha for their valuable comments during the study.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. The study area and location of the Loess Plateau (LP) in China.

Figure 2. The workflow of forest mapping based on the ALOS PALSAR and Landsat TM/ETM+ images.

Figure 3. The backscatter signatures of different land covers in four images. (a) Horizontal transmission and reception (HH) polarization, (b) horizontal transmission and vertical reception (HV) polarization, (c) Difference image (HH − HV), and (d) Ratio image (HH/HV).

Figure 4. The PALSAR-based forest map (a) and the Landsat-based vegetation/non-vegetation map (b) in the Loess Plateau in 2010.

Figure 5. Spatial distribution of the randomly generated validation Regions of Interest (ROIs) in ArcGIS (a). Ground truth samples derived from Google Earth. Two ROIs include the forest region and non-forest region from the very high resolution (VHR) images of Google Earth from around 2010 (b). Using the Collect Earth software to verify the results of ROIs (c).

Figure 6. Spatial distributions of the PL-based and the other five forest maps.

Figure 7. The spatial distribution map of the forest frequency in the Loess Plateau in 2010. The input forest maps were from the PALSAR/Landsat forest map, JAXA forest map, FROM-GLC forest map, and OU-FDL forest map, respectively.

Figure 8. Spatial comparison between the PALSAR/Landsat forest map and (a) FORM-GLC forest map, (b) JAXA forest map, and (c) OU-FDL forest map.

Figure 9. Zoom-in windows of the spatial comparison between the PALSAR/Landsat forest map and FROM-GLC forest map (a,d,g,j); OU-FDL forest map (b,e,h,k); and JAXA forest map (c,f,i,l) in 2010.

Figure 10. Verification of changed areas and unchanged areas using Google Earth on the Loess Plateau between 2007 and 2017. Nonforest areas to Nonforest areas (a–c), Forest areas to Forest areas (d–f), Forest areas to Nonforest areas (g–i), and Nonforest areas to Forest areas (j–l).

Figure 11. Forest variation in the forest cover in the Loess Plateau and the trend of the forest variation cover from 2007 to 2017.

Figure 12. (a–g) Spatial distribution of forest cover in the Loess Plateau from 2007 to 2017 excluding 2011–2014 when the PALSAR data were unavailable.

Figure 13. Spatial pattern of forest changes in the Loess Plateau (a). Unchanged areas of Loess Plateau including unchanged non-forest area (b) and unchanged forest area (c). The forest changes including conversion from forests to non-forest (deforestation) (d) and conversion from non-forest to forests (afforestation) (e).

Figure 14. Reforestation in the pilot county (Wuqi County) of the Grain for Green Project (GGP). The forest map of Wuqi County in 2007 (a) and 2017 (b) the spatial changes of Wuqi County from 2007 to 2017 (c), and the image of Wuqi County and the neighboring counties in 2017 from Google Earth (d).

Table 1

Data specification sheet of the PALSAR and Landsat.

Main Scale of the PALSAR
Mode	High Resolution	Scanning Synthetic Aperture		Polarization
Center frequency	1270 MHz (L-band)
Linear frequency modulation width	28 MHz	14 MHz	14 MHz, 28 MHz	14 MHz
Polarization mode	HH or VV	HH + HV or VV + VH	HH or VV	HH + HV + VH + VV
Incident angle	8°–60°	8°–60°	18°–43°	8°–30°
Spatial resolution	7–44 m	14–88 m	100 m	24–89 m
Width	40–70 km	40–70 km	250–350 km	20–65 km
Satellite parameters of the Landsat
Satellite parameters		TM sensor
Country	United States	Band	Wavelength range (μm)	Resolution (m)
Design life	5	1	0.45~0.53	30
Launch time	1994-04-15	2	0.52~0.60	30
Estimated time to failure		3	0.63~0.69	30
Track type	Near-polar sun synchronization orbit	4	0.76~0.90	30
Track height	705 km	5	1.55~1.75	30
Orbital inclination	98.2°	6	10.40~12.50	60
Operation cycle	98.9 min	7	2.09~2.35	30
Around the earth every day	15	8	0.52~0.90	15

Table 2

Summary of forest products for the comparison in the Loess Plateau in 2010.

Forest Maps	Forest Cover Types	Resolution	Algorithms	Data Sources	Major References
FROM-GLC	Tree cover ≥ 15%, Tree height ≥ 3 m	30 m	Automatic classification algorithms	Landsat	[48]
GlobeLand30	Tree cover ≥ 10%	30 m	POK(pixel-object-konwledge)-based method	Landsat, HJ-1	[49]
GLCF-VCF	Tree cover ≥ 10%, Tree height ≥ 5 m	30 m	Supervised classification	Landsat	[50]
JAXA	Tree cover ≥ 10%, Tree height ≥ 5 m	25 m	Rule-based	PALSAR	[22]
OU-FDL	Tree cover ≥ 10%, Tree height ≥ 5 m	30 m	Decision Tree	PALSAR, MODIS	[26]
PALSAR/Landsat	Tree cover ≥ 10%, Tree height ≥ 5 m	30 m	Rule-based	PALSAR, Landsat	This study

Table 3

Accuracy comparison between PL-based forest map and existing five forest maps with 30–50 m resolutions.

Forest Maps	Type	UA (%)	PA (%)	OA (%)	Kappa
FROM-GLC	Forest	81.26	71.52	87.96	0.68
	Non-Forest	90.02	93.40
GlobeLand30	Forest	78.62	71.22	87.47	0.66
	Non-Forest	90.20	93.19
GLCF-VCF	Forest	72.33	67.30	85.20	0.60
	Non-Forest	89.16	91.27
JAXA	Forest	59.01	67.28	83.57	0.52
	Non-Forest	91.15	87.82
OU-FDL	Forest	65.30	67.47	84.40	0.56
	Non-Forest	90.29	89.41
PALSAR/Landsat	Forest	91.07	76.41	91.27	0.77
	Non-Forest	91.33	97.07