Mapping Plastic-Mulched Farmland with Multi-Temporal Landsat-8 Data

Using plastic mulching for farmland is booming around the world. Despite its benefit of protecting crops from unfavorable conditions and increasing crop yield, the massive use of the plastic-mulching technique causes many environmental problems. Therefore, timely and effective mapping of plastic-mulched farmland (PMF) is of great interest to policy-makers to leverage the trade-off between economic profit and adverse environmental impacts. However, it is still challenging to implement remote-sensing-based PMF mapping due to its changing spectral characteristics with the growing seasons of crops and geographic regions. In this study, we examined the potential of multi-temporal Landsat-8 imagery for mapping PMF. To this end, we gathered the information of spectra, textures, indices, and thermal features into random forest (RF) and support vector machine (SVM) algorithms in order to select the common characteristics for distinguishing PMF from other land cover types. The experiment was conducted in Jizhou, Hebei Province. The results demonstrated that the spectral features and indices features of NDVI (normalized difference vegetation index), GI (greenness index), and textural features of mean are more important than the other features for mapping PMF in Jizhou. With that, the optimal period for mapping PMF is in April, followed by May. A combination of these two times (April and May) is better than later in the season. The highest overall, producer’s, and user’s accuracies achieved were 97.01%, 92.48%, and 96.40% in Jizhou, respectively.


Introduction
Since the mid-20th century, the practice of plastic-mulched farmland (PMF) has revolutionized agricultural production all over the world [1].With the transparent and energy-saving covering materials, plastic mulch can protect crops from unfavorable growing conditions, such as droughts as well as extreme cold and heat.This manner can largely benefit the crop growth in addition to expanding agricultural production to some unfavorable planting areas.PMF plays an increasing role in modern agriculture.On the other hand, it is also criticized for posing environmental concerns, such as "white pollution" and soil degradation [2,3].Furthermore, plastic mulch also cannot be easily decomposed in the natural environment.This might eventually result in the decrease of crop yields.In addition, the plastic mulch used in PMF holds the special physical characteristics of high reflectance and gas tightness.These features can alter the material and energy exchange between the land surface and the atmosphere in the following three aspects: (1) More sunlight is reflected back to the atmosphere, leading to an increase in the atmosphere temperature.(2) The transparency of plastic mulch allows the transmission of shortwave radiation while blocking the emission of longwave radiation, thus increasing the temperature of the soil.(3) Plastic mulch prevents the water vapors produced from evaporating.As a result, PMF might have an impact on the regional climate [4].The environmental problems caused by PMF have been exacerbated over recent years, creating a pressing demand to monitor and optimize the use of PMF.For this, accurate information of the PMF distribution over a large area is needed for policy makers and researchers in the field of land surface temperature, evapotranspiration, crop phenology, and climate change, with this information being a prerequisite for updating the existing land use database.Field surveys are part of the traditional approach to obtaining this information.Nevertheless, this manner is time-consuming and labor-intensive.A long time is needed to update the land use database.More importantly, the data reported at fixed observation locations lack spatial details.The experience of the investigator also influences the accuracy of the field data.
Remote sensing is considered to be a promising technique for acquiring up-to-date information of land use with multiple temporal and spatial resolutions [5,6].During the past years, land cover mapping with remote-sensing data has drawn increasing attention, especially the one-class land cover type extraction, such as a water body [7,8], impervious surface [4,9], agricultural landscape [10][11][12], snow and ice [13,14], special vegetation, crop type identification [15,16], and so on.
In recent years, increasing attention has been paid to plasticulture.This includes an increased focus on mapping the plastic greenhouses with remote sensing compared to plastic-mulched farmland.For example, Levin et al. explored the feasibility of plasticulture landscape mapping using hyperspectral data [17].Agüera et al. proposed a classification scheme using the best band combination of QuickBird images for plastic-greenhouse mapping [18].After that, Agüera et al. developed a pixel-based method for mapping a plastic greenhouse using texture features from high-resolution images [19].Carvajal et al. mapped the plastic greenhouse using QuickBird and Ikonos images [20].Arcidiacono et al. presented a pixel-based classification of high-resolution satellite images for mapping crop-shelter coverage [21].After that, Arcidiacono et al. proposed a per-pixel classification for improving the crop-shelter coverage mapping by texture analyses of high-resolution satellite images [22].Tarantino et al. mapped the plastic covered vineyards using true color aerial data [23].Koc-San evaluated the performance of different classifiers for differentiating between glass and plastic greenhouses using WorldView-2 images [24].Recently, studies [25][26][27] have developed an object-based approach for mapping plastic greenhouses using high spatial resolution images (GeoEye-1 and WorldView-2).All these proposed methods mostly used high spatial resolution images.Although they efficiently mapped plastic greenhouses in their study area, it will be limited by large spatial extent, large data storage, and costly data procurement, in addition to having a time-consuming process.More recently, Wu et al. [28] and Novelli et al. [29] developed an object-based approach for mapping plastic greenhouses using medium spatial resolution images (Sentinel-2 and Landsat-8).Yang et al. [30] presented a new spectral index for mapping plastic greenhouses with medium spatial resolution satellite data.
However, the remote-sensing characteristics and the spatial distribution pattern of plastic-mulched farmland is different from plastic greenhouses.The distribution area of plastic-mulched farmland is larger than that of plastic greenhouses in China, and the spectral response of the plastic-mulched farmland is changing more quickly than that of plastic greenhouses.Mapping PMF with remote sensing began in the last few years.One of the earliest attempts refers to Wang et al., who probed the feasibility of using multi-angle polarization information for to distinguish plastic-mulched paddy fields from the water background [31].They found that plastic mulch has distinct characteristics of polarization reflection, with these characteristics differing with the type of plastic mulch and wavelength channels.Lu et al. developed a decision-tree classifier for the extraction of plastic-mulched cotton farmland in Xinjiang, China using Landsat-5 TM (Thematic Mapper) images [4].In their work, the effectiveness of the decision-tree classifier for PMF mapping was examined.However, the low pace of revisiting of Landsat sensors limits its use for PMF mapping, as it is difficult to collect cloud-free images in the "planting stage" (the time period that plastic mulch is visible).To overcome this problem, they performed PMF mapping utilizing the time series data of normalized difference vegetation index (NDVI) from the moderate-resolution imaging spectrometer (MODIS) [32].Although they identify PMF in an efficient manner, these two methods are limited by several issues.The plastic mulch used in different regions varies in color, thickness, and transparency, which alters the spectral characteristics of PMF.These differences were rarely considered in the existing threshold methods, limiting their performance when applied to other regions.Most importantly, the spectral information of PMF can be constantly altered by the growth cycle of crops.Using the same standard to identify PMF lacks the flexibility to capture seasonal variations, and may obfuscate the true pixel information.In addition, the low spatial resolution imagery is not suitable for small-patched and fragmented agricultural regions, as some smaller PMF are lost and mixed pixels become more serious.Accordingly, more comprehensive consideration of these issues is needed to improve the rigor of PMF mapping with remote sensing.For this, Hasituya et al. extracted the PMF information by using the spectral and textural features from a single-temporal Landsat-8 OLI (Operational Land Imager) imagery with a support vector machine (SVM) classifier.However, they found that the textural features from a single Landsat-8 OLI imagery provide a limited or negligible improvement in mapping accuracy [33].After that, they used the spectral and textural feature of the high spatial resolution GaoFen-1 satellite image to map plastic-mulched farmland with an SVM classifier, and found that the textural features performed better than spectral features when the spatial resolution of images was high enough [34].However, these two studies both used the single temporal remote-sensing data.
In the field of land cover classification with remote sensing, the use of multi-temporal images and multi-types of features is becoming mainstream.A large number of studies have examined the performance of multi-temporal and multi-types of features-including optical and microwave remote sensing-in mapping and change detection.For example, Zoungrana et al. detected the land use/cover change in the Southwest of Burkina Faso using multi-temporal Landsat images and ancillary data [35].Yamazaki et al. developed a global water body map using multi-temporal Landsat images [36].Gao et al. monitored impervious surface expansion using time series medium-resolution satellite images [37].Fisher et al. mapped the large-area tree cover using multi-temporal SPOT-5 images [38].Karale et al. classified crop types based on multi-temporal satellite remote sensing data [39].Wang et al. applied multi-temporal ENVISAT data for mapping agricultural areas in the Pearl River Delta [40].Larrañaga et al. completed a multi-temporal crop classification framework using quad-polarization Radarsat-2 images [41].The results reported that the multi-temporal images capture the specific temporal characteristics of land cover types or the objects of special interest better than the single temporal images.
The specific techniques used in remote-sensing image classifiers are broadly divided into: (1) supervised and unsupervised classifiers by whether or not prior knowledge is needed; (2) parametric supervised classifier and non-parametric supervised classifier (machine learning classifier) by whether or not parameters are needed; (3) sub-pixel based, per-pixel based, and object-based classifiers according to the basic operating unit; and (4) single classifier and ensemble classifier algorithms based on the number of classifiers.The parametric supervised classifier relies on statistical descriptions of training samples and assumes that the data are normally distributed, which is unsuitable for processing high-dimensional and additional data.Furthermore, when using the non-parametric supervised machine learning classifier, there is no need to make assumptions for data distribution, which can usually result in higher accuracy.The non-parametric classifiers, including the neural network classifier, decision tree classifier, support vector machines, and random forest, have great potential in remote-sensing image classification [8,[42][43][44].The methods for mapping plastic greenhouses include conventional supervised classification, object-based methods, machine learning classifiers (neural network, support vector machine, random forest), index-based threshold methods, and so on.However, the methods for mapping plastic-mulched farmland only include index-based threshold methods, decision tree classifiers, support vector machine classifiers, and so on.
From the review, we found that only five papers discussed mapping PMF with remote sensing.The features include spectral reflectance, spectral indexes, and textural features.Other features (e.g., thermal and temporal features) have not been considered until now.To this end, this paper mainly discusses the effectiveness of multi-temporal and multi-types of features for PMF mapping by machine learning classifiers.The goal relies on the combined use of multi-temporal Landsat-8 OLI and Thermal Infrared Sensor (TIRS) imagery.More specific objectives were (1) to examine the optimal period for PMF mapping; (2) to determine the most effective feature or feature set for PMF mapping; (3) to evaluate whether the temporal combination improves the mapping accuracy of PMF; and (4) to compare the performance of two different machine learning classifiers.

Study Area
China is considered to be one of the largest PMF planting countries in the world.By 2013, PMF had a total coverage of 25 × 10 6 ha in China [45,46] (Figure 1).According to Liu et al., PMF has increased the yield of economic food crops by 20%-30% and the yield of grain crops by 20%-60% [47], contributing significantly to the sustainable development of agriculture in China.
In this study, two PMF areas (Figure 2) with different mulching practices were selected for the development and examination of the new scheme.The first study area is located in Jizhou (latitude and longitude range: 37 • 18 40"N-37 • 44 25 N and 115 • 09 57 E-115 • 41 07 E), Hebei Province, China.Covering an area of 9.22 × 10 4 ha, Jizhou is one of the major agricultural production areas (5.93 × 10 4 ha) located on the North China Plain [48].During the past year, it has experienced a considerable increase in PMF.The area is under a temperate monsoon climate, characterized by a hot and rainy summer that favors agricultural production.The rest of the year is relatively unfavorable for farming.The crops here mainly include cotton, winter wheat, corn, and vegetables.According to our field survey, cotton fields dominate the use of plastic mulch in this area.Other land use types consist of woodland, grassland, water body, and impervious surface (traffic land, residential land, and industrial land).
The second study area is seated in Guyuan, Ningxia Hui Autonomous Region, China.It covers an area of 5.19 × 10 6 ha, 21% of which is dominated by farmland [49,50].As one of the ecological barriers in the western part of China, Ningxia is characterized by a temperate semi-arid climate.The annual average temperature is around 5.4-10.0• C and the annual precipitation is about 169.5-611.8mm.The area is covered by frost for most days in a year, leaving a frost-free time of around 148 days [51].Irrigation is needed to facilitate the cultivation of crops.For the purpose of water-saving in farmland, PMF has been widely used here for planting corn, winter wheat, and vegetables.It was observed that early spring and autumn are the main times for crop mulching.machine learning classifiers.The goal relies on the combined use of multi-temporal Landsat-8 OLI and Thermal Infrared Sensor (TIRS) imagery.More specific objectives were (1) to examine the optimal period for PMF mapping; (2) to determine the most effective feature or feature set for PMF mapping; (3) to evaluate whether the temporal combination improves the mapping accuracy of PMF; and (4) to compare the performance of two different machine learning classifiers.

Study Area
China is considered to be one of the largest PMF planting countries in the world.By 2013, PMF had a total coverage of 25 × 10 6 ha in China [45,46] (Figure 1).According to Liu et al., PMF has increased the yield of economic food crops by 20%-30% and the yield of grain crops by 20%-60% [47], contributing significantly to the sustainable development of agriculture in China.
In this study, two PMF areas (Figure 2) with different mulching practices were selected for the development and examination of the new scheme.The first study area is located in Jizhou (latitude and longitude range: 37°18′40″N-37°44′25′′N and 115°09′57′′E-115°41′07′′E), Hebei Province, China.Covering an area of 9.22 × 10 4 ha, Jizhou is one of the major agricultural production areas (5.93 × 10 4 ha) located on the North China Plain [48].During the past year, it has experienced a considerable increase in PMF.The area is under a temperate monsoon climate, characterized by a hot and rainy summer that favors agricultural production.The rest of the year is relatively unfavorable for farming.The crops here mainly include cotton, winter wheat, corn, and vegetables.According to our field survey, cotton fields dominate the use of plastic mulch in this area.Other land use types consist of woodland, grassland, water body, and impervious surface (traffic land, residential land, and industrial land).
The second study area is seated in Guyuan, Ningxia Hui Autonomous Region, China.It covers an area of 5.19 × 10 6 ha, 21% of which is dominated by farmland [49,50].As one of the ecological barriers in the western part of China, Ningxia is characterized by a temperate semi-arid climate.The annual average temperature is around 5.4-10.0°C and the annual precipitation is about 169.5-611.8mm.The area is covered by frost for most days in a year, leaving a frost-free time of around 148 days [51].Irrigation is needed to facilitate the cultivation of crops.For the purpose of water-saving in farmland, PMF has been widely used here for planting corn, winter wheat, and vegetables.It was observed that early spring and autumn are the main times for crop mulching.

Remote-Sensing Data
Ten Landsat-8 OLI (Operational Land Imager) images (five images for each study area) were collected to develop and validate the new scheme.The five images (path/row: 114/27) for Jizhou were respectively taken on 16 April, 18 May, 3 June, 19 June, and 5 July in 2015.The other five images for the Guyuan area were respectively acquired on 26 April, 12 May, 13 June, 15 July, and 31 July in 2015.Radiometric calibration and atmospheric correction were conducted using the Fast Line of Sight Atmospheric Analysis of Hypercubes (FLAASH) module in the Environment of Visualizing Images (ENVI) software.The brightness temperature was also acquired from the Landsat TIRS (Thermal Infrared Sensor) images by applying the predesigned parameters in the Metadata file.
Eight high spatial resolution GF-1 (GaoFen-1, Chinese satellite) images (2 m for the panchromatic band and 8 m for the multispectral bands with width of 60 km) were also collected as the reference data source for the collection of samples.This dataset contains six images for Jizhou area: four images taken on 5 May 2015 and two images taken on 11 June 2015.Two GF-1 images (acquisition date: 8 April 2015) were collected for the Guyuan area.For the multispectral bands, radiometric calibration and atmospheric correction were carried out using the ENVI software.The multispectral images were then downscaled to a spatial resolution of 2 m using a Gram-Schmidt Pan Sharpening module.These pan-sharpened GF-1 images were then mosaicked to produce the reference images for the two study areas.Geo-referencing between the Landsat-8 OLI and the pansharpened GF-1 images was completed in ENVI with 30 ground control points (GCP).The final registration error was smaller than one GF-1 pixel (i.e., 2 m) evaluated by 20 other independent GCPs.

Remote-Sensing Data
Ten Landsat-8 OLI (Operational Land Imager) images (five images for each study area) were collected to develop and validate the new scheme.The five images (path/row: 114/27) for Jizhou were respectively taken on 16 April, 18 May, 3 June, 19 June, and 5 July in 2015.The other five images for the Guyuan area were respectively acquired on 26 April, 12 May, 13 June, 15 July, and 31 July in 2015.Radiometric calibration and atmospheric correction were conducted using the Fast Line of Sight Atmospheric Analysis of Hypercubes (FLAASH) module in the Environment of Visualizing Images (ENVI) software.The brightness temperature was also acquired from the Landsat TIRS (Thermal Infrared Sensor) images by applying the predesigned parameters in the Metadata file.
Eight high spatial resolution GF-1 (GaoFen-1, Chinese satellite) images (2 m for the panchromatic band and 8 m for the multispectral bands with width of 60 km) were also collected as the reference data source for the collection of samples.This dataset contains six images for Jizhou area: four images taken on 5 May 2015 and two images taken on 11 June 2015.Two GF-1 images (acquisition date: 8 April 2015) were collected for the Guyuan area.For the multispectral bands, radiometric calibration and atmospheric correction were carried out using the ENVI software.The multispectral images were then downscaled to a spatial resolution of 2 m using a Gram-Schmidt Pan Sharpening module.These pan-sharpened GF-1 images were then mosaicked to produce the reference images for the two study areas.Geo-referencing between the Landsat-8 OLI and the pan-sharpened GF-1 images was completed in ENVI with 30 ground control points (GCP).The final registration error was smaller than one GF-1 pixel (i.e., 2 m) evaluated by 20 other independent GCPs.

Field Survey Data
To collect the real-field information of land use, a field survey was carried out during 25-30 April 2015.In the trial, we recorded both the land use types and the coordinates of the sampling points.The land use types mainly consist of PMF, bare soil, impervious surface, vegetation, and water in Jizhou, while the land use types mainly consist of PMF, bare soil, impervious surface, vegetation, water, plastic greenhouse, and mountain areas in Guyuan.After that, we digitized these sampling points into the GF-1 and Landsat-8 OLI images to produce the region samples.Instead of using the real-field sampling points, we enlarged the sampling region to a window of 60 m × 60 m on the GF-1 image (2 × 2 Landsat pixels).It is important to note that we slightly relocated the windows to guarantee the purity of each window.The samples were then modified by comparison with the Landsat images in each season to guarantee the consistency of samples across growing seasons.As shown in Table 1 and Figure 2, a total of 1136 samples were eventually collected for the Jizhou area and another 648 samples were collected for the Guyuan area.The samples of each study area were then equally divided into two groups: training samples and testing samples.

Methodology
The features play a key role in the extraction of PMF.A better combination of features can significantly improve the mapping performance.Therefore, the main idea of this present study is to optimize the selection of features for PMF mapping by integrating the temporal and spatial information.Figure 3 shows the framework of the proposed scheme.It mainly consists of four steps: spectral characteristic analysis, feature extraction, as well as feature selection and mapping.

Field Survey Data
To collect the real-field information of land use, a field survey was carried out during 25-30 April 2015.In the trial, we recorded both the land use types and the coordinates of the sampling points.The land use types mainly consist of PMF, bare soil, impervious surface, vegetation, and water in Jizhou, while the land use types mainly consist of PMF, bare soil, impervious surface, vegetation, water, plastic greenhouse, and mountain areas in Guyuan.After that, we digitized these sampling points into the GF-1 and Landsat-8 OLI images to produce the region samples.Instead of using the real-field sampling points, we enlarged the sampling region to a window of 60 m × 60 m on the GF-1 image (2 × 2 Landsat pixels).It is important to note that we slightly relocated the windows to guarantee the purity of each window.The samples were then modified by comparison with the Landsat images in each season to guarantee the consistency of samples across growing seasons.As shown in Table 1 and Figure 2, a total of 1136 samples were eventually collected for the Jizhou area and another 648 samples were collected for the Guyuan area.The samples of each study area were then equally divided into two groups: training samples and testing samples.

Methodology
The features play a key role in the extraction of PMF.A better combination of features can significantly improve the mapping performance.Therefore, the main idea of this present study is to optimize the selection of features for PMF mapping by integrating the temporal and spatial information.Figure 3 shows the framework of the proposed scheme.It mainly consists of four steps: spectral characteristic analysis, feature extraction, as well as feature selection and mapping.

Spectral Response Analysis
Identifying the distinct spectral characteristics is critical for the development of the PMF mapping method.In this study, pixel-based spectral curves of typical land cover types were manually extracted from the multi-temporal Landsat-8 OLI images.As the sample area is relatively homogeneous, we sampled a total of 200 Landsat pixels from each atmospheric-corrected image.These samples included 40 pixels under five classes: PMF, impervious surface (IS), bare soil (BS), vegetation cover (VC), and water body (WB).The mean spectra of each class were then derived.
Figure 4 shows the mean spectra of typical land cover types in Jizhou.It is evident that the spectral characteristics of PMF change over time.During the planting stage, the spectrum of PMF more resembles the BS and IS.The vegetation signal dominates the spectral response of PMF when the plastic mulch is covered by vegetation.Therefore, it is difficult to separate the PMF from other land cover types based solely on a single image.Accordingly, we conducted PMF mapping by combining the multi-temporal Landsat-8 OLI images.

Spectral Response Analysis
Identifying the distinct spectral characteristics is critical for the development of the PMF mapping method.In this study, pixel-based spectral curves of typical land cover types were manually extracted from the multi-temporal Landsat-8 OLI images.As the sample area is relatively homogeneous, we sampled a total of 200 Landsat pixels from each atmospheric-corrected image.These samples included 40 pixels under five classes: PMF, impervious surface (IS), bare soil (BS), vegetation cover (VC), and water body (WB).The mean spectra of each class were then derived.
Figure 4 shows the mean spectra of typical land cover types in Jizhou.It is evident that the spectral characteristics of PMF change over time.During the planting stage, the spectrum of PMF more resembles the BS and IS.The vegetation signal dominates the spectral response of PMF when the plastic mulch is covered by vegetation.Therefore, it is difficult to separate the PMF from other land cover types based solely on a single image.Accordingly, we conducted PMF mapping by combining the multi-temporal Landsat-8 OLI images.

Feature Extraction
In this study, a total of 84 types of features (Table 2) were derived from the multi-temporal Landsat-8 images.These features can be categorized into four types: spectral features, textural features, indices features, and thermal features.We extracted these 84 features for each Landsat-8 image, eventually producing a total of 420 features for each study area.

Feature Extraction
In this study, a total of 84 types of features (Table 2) were derived from the multi-temporal Landsat-8 images.These features can be categorized into four types: spectral features, textural features, indices features, and thermal features.We extracted these 84 features for each Landsat-8 image, eventually producing a total of 420 features for each study area.

Textural Features
The PMF always exhibits distinct textural characteristics compared with other land cover types.For example, the PMF areas are regularly separated by non-plastic-mulched crops to prevent flooding when there is heavy rain.This manner can largely increase the heterogeneity of the PMF area.In addition, the plastic film used in PMF is quite smooth, which can also increase the surface homogeneity.Therefore, it is worth taking textural information into consideration in PMF mapping.
In this study, eight types of textural features were derived from the multi-temporal Landsat images.A gray level co-occurrence matrix (GLCM) method was applied here to obtain the textural features to ensure good performance in texture analysis [52][53][54][55][56].These features include the mean, variance, homogeneity, contrast, dissimilarity, entropy, angular second moment, and correlation calculated in a window of 3 × 3 pixels.We extracted the eight types of textural features for each band, eventually generating a total of 56 textural features.

Indices Features Tasseled Cap Transformation (TCT)
As plastic mulching can change the soil conditions in terms of temperature and moisture, the indices produced by the tasseled cap transformation (TCT; brightness, greenness, and wetness) were taken as a measure of this effect for PMF mapping.
TCT is a linear transformation of spectral information based on experience.It reduces the dimension of an image into three new physical features (brightness, greenness, and wetness) with minimal inter-correlation.The "brightness" is a weighted sum of all the bands.It represents the variability of the imagery.It is generally related to whether the soil is completely or partially bare, the presence of natural and man-made objects, as well as the variations in topography.The "greenness" is a measure of the conditions of the vegetation.It is determined by the coverage, leaf area, and biomass of vegetation.Soil is characterized by a high brightness value and low greenness value, while the forest shows a low brightness value and a low greenness value.The "wetness" is produced in an orthogonal direction to the first two components.It is a measure of the water content in the soil and vegetation [57].The TCT is conducted in the ENVI software using the predesigned parameters for the Landsat-8 OLI sensor.The transformation coefficients (Table 3) recommended by Muhammad et al. [57] were adopted here for the TCT.Equation (1) shows the linear function applying the coefficients.
where X is a vector of satellite image bands; Y is the vector of brightness, greenness, and wetness; c is the transformation coefficient vector (Table 3); and a is an increment value, which can offset the negative values.

NDVI and NDBI
The normalized difference vegetation index (NDVI) [32] and normalized difference built-up index (NDBI) [58] were taken into account for PMF mapping in this study.NDVI is a pixel-based measure of the vegetation, while NDBI is a measure of built-up area within a Landsat pixel.A higher NDVI means better vegetation conditions, while a higher NDBI means a greater built-up area.Hasituya et al. [33] pointed out that the PMF was frequently confused with impervious surface.Therefore, the incorporation of NDBI is expected to solve this problem.They are defined as follows: where b 4 and b 5 and b 6 are, respectively, the red, NIR, and SWIR-1 bands of the Landsat-8 sensor.

Thermal Features
In this study, we probed the potential of the thermal information for PMF mapping.The brightness temperature was provided by the Landsat-8 TIRS images.It can be computed via the following equation: where T b is at-satellite brightness temperature (K); L λ is the TOA (top of atmosphere) spectral radiance (W/(m 2 × srad × µm)); and K 1 and K 2 are the thermal conversion constants from Landsat metadata (K 1−1 = 774.89,K 1−2 = 480.89,K 2−1 = 1321.08,K 2−2 = 1201.14).In order to keep consistency between different indices, the brightness temperature was rescaled under Equation ( 5) [59], stretching to the range of (0,1).
where T br is the rescaled brightness temperature; T b is the brightness temperature; min(T b ) is the minimum value of the brightness temperature; and max(T b ) is the maximum value of the brightness temperature.

Random Forest for Feature Selection
The efficiency and accuracy of mapping with remote sensing will be improved by selecting the optimal feature subset [60,61].To find effective features for PMF mapping, the feature selection approach-which reduces the correlation and redundancy of features-was applied in our experiment.Random forest (RF), first developed by Breiman in 2001 [62], is an ensemble method for supervised classification based on classification and regression trees (CART).RF is well-known because it is efficient to compute, robust to outliers and noise, in addition to being useful for estimating error, strength, correlation, and variable importance [63].RF is also an effective feature selection algorithm [63,64].In this paper, RF was used to optimize the features for PMF mapping by using the importance of features evaluated by RF [65].A detailed process for measuring the importance of features by RF was presented by Guan et al. [66].
In our experiment, RF was performed to measure the importance of all features in each single temporal feature set and each multi-temporal feature set, respectively.Two parameters were set beforehand: the number of trees and the number of variables.A total of 500 trees were grown each time, and the square root of the number of total input features were used as the number of split variables in this paper.Due to the randomness of RF, calculated feature importance varied slightly across different runs.To avoid these differences, we repeated the feature importance measurements ten times and calculated the average importance value.Following this, the features were sorted in descending order of their average importance, and the cumulative average importances were calculated for each single temporal feature.Next, we developed feature sets according to the cumulative percentage of feature importance for mapping PMF.For these features, the cumulative percentages were 30%, 50%, 80%, 90%, and 100%.After this, we chose the feature set with the highest accuracy to develop the multi-temporal combined features.
After the development of a feature set for mapping PMF in each temporal and each multi-temporal combination, two supervised machine learning classifiers, Random Forest (RF) and Support Vector Machine (SVM), were performed for mapping PMF.Finally, we selected the optimal feature set for mapping PMF according to the classification accuracy.

Random Forest and Support Vector Machine for Classification
To provide the baseline for the efficiency of RF, the support vector machine (SVM)-another machine learning supervised classifier-was used for the same features and the same samples for mapping PMF.SVM is designed to find an optimal hyperplane as a decision function in high-dimensional space based on statistical learning theory.The SVM approach uses the principle of structural risk minimization, not the principle of empirical risk minimization [67].Therefore, the performance of SVM has been examined by many studies [68][69][70][71].
In this paper, the RF-supervised classifier was performed to map PMF in each single temporal and multi-temporal combination based on the selected features by RF and the collected training samples.
In this study, the classification was carried out using EnMAP-Box, an IDL based tool for remotely sensed imagery classification [72].

Accuracy Assessment Accuracy Assessment Using Confusion Matrix
The classification accuracies were assessed quantitatively by using the confusion matrix, which is a common used method in remote sensing.This study used overall accuracy (OA), which is computed by dividing correctly classified pixels by the total number of pixels, and the kappa coefficient (K), which considers the whole confusion matrix instead of using only diagonal elements.In addition, producer's accuracy (PA) and user's accuracy (UA) were used to assess the accuracies of individual classes.

Statistical Significance Testing for Classification Accuracies Using Z-Test
The classification accuracies were further confirmed by using the Z-test [73], which is a method for testing the statistical significance of the K statistic and significance differences of different classification schemes.
With this test, it is possible to statistically compare two analysts, the same analyst over different times, two algorithms, two types of imagery, or even two dates of imagery for examining which produces higher accuracy.To verify the effectiveness of different feature sets and different classifiers, the Z-test was performed on the pairwise error matrix of different analysts.The test statistic for testing if two independent error matrices are significantly different is expressed by the following formula: where K 1 and K 2 denote the estimates of the Kappa statistic for error matrix 1 and 2, respectively; var (K 1 ) and var (K 2 ) are the corresponding estimates of the variance of K as computed from the appropriate equations.At the 99% confidence level, the critical value would be 2.58.Therefore, if the absolute value of the test Z statistic is greater than 2.58, the two analysts are significantly different.

Random Forest Feature Selection for Single Temporal PMF Mapping in Jizhou
Figure 5 shows that there were obvious differences in the importance of different features, which changed with the development of the growing season.In the whole growing season, the spectral features took a dominant position, the indices features and textural features (except for the mean) were inferior to spectral features, while the thermal features were inferior to the others.In 16 April, the spectral and indices features were more important than the other features.In 18 May and 3 June, the spectral and indices features and textural feature of mean were more important than the other features.The importance of some textural features in 18 May were slightly higher than that in 16 April, although this importance was reduced after this time period.In 19 June and 5 July, the spectral and indices features were also the more important feature for PMF mapping.
Table 4 shows the 20 most important features selected by RF in different growing seasons.In the early mulching stage (16 April), the first derivative of the red band, greenness index, NDVI, the red band, and the second derivative of red band were sorted as the five most important features.In the early-mid mulching stage (18 May), the mean of NIR, the first derivative of green band, NDVI, mean of coastal aerosol band, and greenness index were sorted as the five most important features.In the middle mulching stage (3 June), the mean of SWIR2, SWIR2, mean of NIR, the first derivative of SWIR1, and the first derivative of the red band were selected as the five most important features.In the mid-later mulching stage (19 June), the first derivative of SWIR2, the second derivative of SWIR1, the first derivative of red band, the second derivative of SWIR2, and the greenness index were chosen as the five most important features.In the later mulching stage (5 July), the second derivative of NIR, the first derivative of SWIR2, the first derivative of the red band, the second derivative of SWIR1, and the mean of NIR were identified as the five most important features for PMF mapping.
Remote Sens. 2017, 9, 557 12 of 27 SWIR1, and the first derivative of the red band were selected as the five most important features.In the mid-later mulching stage (19 June), the first derivative of SWIR2, the second derivative of SWIR1, the first derivative of red band, the second derivative of SWIR2, and the greenness index were chosen as the five most important features.In the later mulching stage (5 July), the second derivative of NIR, the first derivative of SWIR2, the first derivative of the red band, the second derivative of SWIR1, and the mean of NIR were identified as the five most important features for PMF mapping.2).2).Generally, the spectral features (especially the derivative of spectral reflectance) were the most important features for mapping PMF, while the textural feature (mean) and the indices features (NDVI and GI) were also important.In the early mulching stage, the spectral features of visible-near infrared were more important than the other features.From the early-mid to later mulching stage, the spectral and mean textural features of the near and short-wave infrared became more important for mapping PMF in Jizhou.The thermal features did not appear in the ten most important features in our experiment.

Mapping PMF in Jizhou using Single Temporal Features
We classified the cover types in Jizhou by RF and SVM based on the developing features.These cumulative percentages were 30%, 50%, 80%, 90%, and 100%.To compare the effectiveness of selected optimal features, we also classified the land cover types using spectral features alone (seven bands) and the textural features alone (56 textural features).All the OA, Kappa, PA, and UA of RF and SVM are presented in Table 5.Table 5 reveals that the mapping accuracy of PMF varied across different features, different periods, and different classifiers.When using the RF, all the accuracies from the optimal feature set were higher than that from the spectral features alone or textural features alone, especially in the early mulching stage (April and May).It is not so obvious and regular when using SVM.In the later mulching stage (June to July), these characteristics are not so obvious for either RF or SVM.The accuracies of April and May were significantly higher than those of June to July.Therefore, the early mulching stages (mid-April and mid-May) are the optimal periods for mapping PMF in Jizhou.
There were also differences between these two periods.When using the RF classifier, all the features provided the highest accuracy in April and May.The OA, PA, UA, and Kappa were 94.02%, 91.67%, 93.18%, and 0.92, respectively in April, while these values were 94.87%, 86.79%, 91.83%, and 0.93, respectively, in May.The second-highest accuracy provided by 50% of the features occurred in April, the OA, PA, UA, and Kappa were 92.99%, 91.06%, 90.69%, and 0.91 respectively.In May, the second-highest accuracy was provided by 80% of the features, the OA, PA, UA, and Kappa were 94.42%, 86.79%, 90.66%, and 0.93, respectively.The highest PA and UA of the mapping PMF were provided by 50% of features in April by using RF.
When using SVM, 80% of the features provided the highest accuracy in April, the OA, PA, UA, and Kappa were 91.43%, 86.59%, 93.01%, and 0.89, respectively.One hundred percent of the features provided the highest accuracy in May, the OA, PA, UA, and Kappa were 95.63%, 88.82%, 92.39%, and 0.94, respectively.Therefore, the remote-sensing data in April is better for PMF mapping than May.However, after the middle mulching stage (June to July), all the PA and UA values of PMF were lower than 80% using both RF and SVM.Furthermore, later periods resulted in lower accuracies being obtained.Therefore, we can get the preliminary conclusion that the optimal phase of PMF mapping is between April and mid-May.While there were some differences between classifiers, the accuracy generated from RF was higher than that from SVM.
The Z-test was performed to confirm that this classification was meaningful and significantly better than a random classification, and to assess the significance of differences in accuracies of classification obtained from different feature sets and classifiers.The Z-test values between pairs of features and classifiers are given in Table 6.Table 6 shows that the Z-test value was 4.69 (higher than 2.58) when comparing the feature set with the highest accuracy and with the poorest accuracy derived from RF.This means that the performance of these two feature sets was significantly different at the 99% confidence level when using RF.Similarly, when using SVM, the performance of different feature sets differed significantly because the Z-test value was 4.64 (higher than 2.58).For the different classifiers, the Z-test value was 4.36 (higher than 2.58) when comparing the highest accuracy of RF and SVM.The Z-test value was 2.83 (higher than 2.58) when comparing the poorest accuracy of RF and SVM.Therefore, the performance of RF and SVM was significantly different at the 99% confidence level.

Mapping the PMF in Jizhou Using Multi-Temporal Combined Features
In the multi-temporal experiment, we combined the optimized features that produced the highest accuracy, respectively, in 16 April, 18 May, 3 June, and 19 June.The feature combination includes a two-temporal combination, three-temporal combination, and four-temporal combination (Table 7).RF and SVM classifiers were performed based on the selected multi-temporal features.
Table 7 shows that the highest OA, PA, UA, and Kappa of RF were 97.01%, 93.29%, 96.67%, and 0.96.These were higher than the highest accuracy from single temporal features.The two-temporal combined images of 16 April and 18 May provided the highest accuracy, followed by the two-temporal combined images of 16 April and 3 June, and then the three-temporal combined images of 16 April, 18 May, and 3 June.The highest OA, PA, UA, and Kappa of SVM were 95.80%, 89.63%, 94.23%, and 0.95, produced by the two-temporal combined images of 16 April and 18 May.In the multi-temporal combination experiment, the RF also outperforms the SVM.
In this section, we also performed the Z-test.The Z-test values between pairs of features and classifiers are given in Table 8.Table 8 shows that the Z-test value from the comparison of the highest accuracy from single and multi-temporal features using RF was 6.32, which is higher than 2.58.This Z-test value was 10.82 when comparing the poorest accuracy from single and multi-temporal features using RF.Similarly, when using SVM, the Z-test values were 5.95 and 10.52, respectively, for single and multi-temporal features with highest accuracy and the poorest accuracy.This indicates that the performance of single and multi-temporal features was significantly different at the 99% confidence level when using RF and SVM.For comparing the RF and SVM, the Z-test value was 4.75 using the single temporal features, while the Z-test value was 4.36 using the multi-temporal features.Therefore, the performance of RF is significantly better than SVM using both single and multi-temporal features.We evaluated the importance of multi-temporal features using RF to analyze the contribution of features in a multi-temporal combination, in which DOY (Day of Year) 106 is 16 April, DOY 138 is 18 May, DOY 154 is 3 June, and DOY 170 is 19 June. Figure 6 reveals that the features in April and May were more important than the others in all combinations.Among the two-temporal combinations, the five most important features included the 106-GI in the 106-170 combination, 138-G-1 in the 106-138 combination, 154-SWIR2-M in the 154-170 combination, 106-R-1 in the 106-154 combination, and 138-G-1 in the 138-154 combination.Among the three-temporal combinations, the five most important features included the 138-G-1 in the 106-138-154 combination, 106-GI in the 106-154-170 combination, 138-G-1 in the 138-154-170 combination, 138-G-1 in the 106-138-170 combination, and 138-GI in the 138-154-170 combination.Among the four-temporal combinations, the five most important features included the 138-G-1, 138-SWIR, 138-GI, 138-SWIR1-1, and 138-NIR-M.

The Spatial Distribution of PMF in Jizhou
Figures 7 and 8 show the spatial distribution of PMF generated from different feature subsets and classifiers in different periods.Figures 7 and 8 reveal that there were large differences between different features and between different times.The PMF from April had a smaller distribution than that from the other periods.The results from June to July had serious commission errors.The results from April and May were more acceptable and reasonable in light of field investigation and higher resolution images.The reasons for these results include the PMF being seriously confused with the bare soil in April, while the remote-sensing signature of PMF was influenced by the well-developed plastic-mulched crops in June to July.When using the spectral features alone, there was a greater saltand-pepper effect in the distribution of PMF-especially from May to July.When using the textural features alone, the salt-and-pepper effect was relieved to some extent, but there were also more commission and omission errors.The optimized features provided better spatial distribution for PMF.The temporal combination provided the best result.The distribution maintained good consistency with the knowledge from field investigation.Generally, the PMFs were concentrated in the middle region in addition to being scattered in the southern and northern regions.The RF and SVM followed a similar general trend to the distribution of PMF.

The Spatial Distribution of PMF in Jizhou
Figures 7 and 8 show the spatial distribution of PMF generated from different feature subsets and classifiers in different periods.Figures 7 and 8 reveal that there were large differences between different features and between different times.The PMF from April had a smaller distribution than that from the other periods.The results from June to July had serious commission errors.The results from April and May were more acceptable and reasonable in light of field investigation and higher resolution images.The reasons for these results include the PMF being seriously confused with the bare soil in April, while the remote-sensing signature of PMF was influenced by the well-developed plastic-mulched crops in June to July.When using the spectral features alone, there was a greater salt-and-pepper effect in the distribution of PMF-especially from May to July.When using the textural features alone, the salt-and-pepper effect was relieved to some extent, but there were also more commission and omission errors.The optimized features provided better spatial distribution for PMF.The temporal combination provided the best result.The distribution maintained good consistency with the knowledge from field investigation.Generally, the PMFs were concentrated in the middle region in addition to being scattered in the southern and northern regions.The RF and SVM followed a similar general trend to the distribution of PMF.

Contribution of Feature Optimizing
The contribution of feature optimizing was explained by comparing the factors in the confusion matrix of the classification of land cover types with the poorest (textural features alone) and the highest accuracy (two-temporal combined features from April and May) before and after the feature optimizing.Table 9 shows that the factors on the diagonal were smaller and the other factors were relatively greater in the confusion matrix with the poorest accuracy.These trends were found to be opposite in the confusion matrix of the highest accuracy.This means that the feature optimization

Contribution of Feature Optimizing
The contribution of feature optimizing was explained by comparing the factors in the confusion matrix of the classification of land cover types with the poorest (textural features alone) and the highest accuracy (two-temporal combined features from April and May) before and after the feature optimizing.Table 9 shows that the factors on the diagonal were smaller and the other factors were relatively greater in the confusion matrix with the poorest accuracy.These trends were found to be opposite in the confusion matrix of the highest accuracy.This means that the feature optimization can relieve serious confusion between land cover types, in particular the confusion between PMF and BS or IS.From Table 9, we can see that the confusion between PMF and BS was reduced from 18.30% to 7.47% and then to 4.38%, while the confusion between BS and PMF was reduced from 26.02% to 7.11% and then to 6.5%.These three values correspond, respectively, to the progression from the single-temporal with the poorest accuracy to the single-temporal with the highest accuracy and then to the multi-temporal with the highest accuracy by optimizing features.Therefore, the spectral feature of visible-near infrared band (red, green, NIR), indices features (NDVI, GI), and the textural feature (mean) and multi-temporal combination of these features can improve the mapping accuracy significantly.The plastic mulching in Jizhou was completed between April to May.During this period, the remote-sensing characteristics of PMF are influenced by the dust, rain, and by the phenology of the plastic-mulched crop.In addition, the plastic mulch currently used in China is a very thin (0.006-0.008 mm) transparent plastic mulch, and it is placed tightly over the soil surface, so its remote-sensing characteristics are considerably influenced by the characteristics of the covered soil.After the completion of the plastic film mulching, the remote-sensing characteristics of PMF change into the mixed characteristics of plastic mulch and soil if there is wind or rain.Furthermore, with the emergence of the mulched crop, the remote-sensing characteristics of PMF show the characteristics of vegetation, plastic mulch, and soil.In this situation, it makes it difficult to identify PMF.Therefore, the optimal temporal period of mapping PMF is well within the period from completing the mulching to the emergence stage of mulched crop, which occurs before the stage of being covered by well-developed crops.
Due to the specific characteristics of plastic mulch and the mulching mode, the PMF can also be confused with bare land and fallow land in the mulching period (April).Furthermore, in the emergence stage of the mulched crop, PMF contains little green vegetation, while the bare land and fallow land still has no vegetation.Therefore, the combination of the sowing stage (April) and the emergence stage (May) generated the highest accuracy.

Applicability in Another Region
Due to the regional division of natural conditions and the different agricultural production modes, the applicability of methods and techniques in different regions is uncertain.Therefore, we chose Guyuan as an experimental area to test the applicability of the presented technique scheme.The PMFs in Guyuan were mapped following the same scheme as Jizhou.
Table 10 reveals that the PMF mapping accuracy in Guyuan was similar with Jizhou in regularity.The accuracies derived from the optimal features were higher than that from the original features or a single feature.The accuracies of April and May were significantly higher than later in the season.The OA, PA, UA, and Kappa were 85.02%, 74.07%, 81.40%, and 0.79, respectively, in April, while they were 88.13%, 73.60%, 79.53%, and 0.83, respectively, in May.Therefore, the early stages of mulching (April-May) are the optimal periods for mapping PMF in Guyuan.However, the highest accuracies were lower than that in Jizhou.The Z-test value for comparing the differences of the single temporal features with the highest accuracy and the poorest accuracy was computed, with an obtained value of 8.45 (higher than 2.58).Therefore, the feature addition significantly improved the single temporal mapping accuracy.
We used the optimal features selected by RF to develop the multi-temporal combined features for mapping PMF.For each multi-temporal combined feature, we used the RF to classify the land cover types.RF classifier was performed for mapping PMF based on the selected multi-temporal features in Guyuan.Table 11 shows that the highest OA, PA, UA, and Kappa were 90.89%, 80.90%, 87.86%, and 0.87.This was higher than that from single temporal features generally.The three-temporal combination of May, June, and July provided the highest accuracy, followed by the two-temporal combination of May and July.
The Z-test value for comparing the single and multi-temporal features with the highest accuracy and poorest accuracy were calculated, with the multi-temporal mapping with the highest accuracy in Jizhou and Guyuan also being computed.The Z-test results were presented in Table 12.Table 12 shows that the multi-temporal combination significantly improved the mapping accuracy in Guyuan.Furthermore, there were also significant differences in the two different study regions (Z-test value is 18.50).Figure 9 (depicting the spatial distribution of PMF in Guyuan) reveals that large differences exist among different feature subsets and among different seasons.The later mapping period results in more serious commission errors.The PMF distribution generated from the multi-features were better than that from single-features, the distribution from optimized features was better than that from non-optimized features, while the distribution from multi-temporal optimized features was more acceptable than that from single-temporal features.
By comparing the results from the two study areas, we found that the general trends were consistent, although there were also differences, such as the accuracy, temporal combination type, and so on.This inconsistency may be attributed to the land cover types and the agricultural production mode.This should be discussed in further studies.
Generally, the multi-temporal optimized features provided more acceptable results than single-temporal features.Among the single-temporal features, the earlier stages were better than the later stages.The results from April and May were more acceptable and reasonable.This is because the remote-sensing signature of PMF is influenced by the well-developed plastic-mulched crops at later stages.
The difference between these two study areas can be attributed to the land cover types and their distribution pattern.The land cover types in Jizhou are simpler and distributed more uniformly, while the land cover types in Guyuan are more complex and distributed unevenly.In addition, the mulching mode in Guyuan includes mulching in autumn, mulching in early spring, and mulching before sowing.The data used in this paper covered the time period from April to July.Therefore, there may be some discrepancies between the data acquisition time and the mulching time.As a result, the mapping accuracy in Guyuan was lower than that in Jizhou.

Differences Between Classifiers
In this paper, two machine learning classifiers-RF and SVM-were used for mapping PMF in Jizhou.The result indicated that RF provided a higher accuracy for mapping PMF than the SVM.As the SVM classifier is a single machine learning algorithm and the RF is an ensemble machine learning algorithm, RF is better than the SVM in the aspects of noise resistance, robustness, and computational efficiency.Furthermore, RF is more effective in the high-dimensional feature space and can provide higher accuracy.

Comparison with Existing Results
From the review of related studies, we discovered that little research has examined PMF mapping with remote sensing.Two of these studies used the threshold model based on the spectral or the indices features from Landsat TM and MODIS-NDVI.However, due to the changing spectral response of PMF during the crop growth cycle and over regions as well as the very fragmented and small-patched farmland in China, there are seriously mixed pixels in the low spatial resolution imagery.
Table 13 shows that in the case of the same data source and the same study area, the scheme presented in this paper had significantly improved accuracy.The OA improved from 94.14% to 97.01%.Furthermore, in different study areas or from different data sources, the OA of this paper was higher or at a comparable level to the results from other research.In this paper, we considered the multi-temporal and multi-types of features simultaneously to improve the accuracy and efficiency of PMF mapping.Thus, this paper provided more accurate results than that from the existing research.Multi-temporal and multi-types of features can perform better than single temporal features and single-types of features.

Conclusions
In this paper, the potential of multi-temporal Landsat-8 images was evaluated to determine the optimal phase and effective features sets for mapping PMF.The single-temporal multi-features and multi-temporal multi-types of features were developed for RF and SVM classifiers.Our preliminary conclusions were as follows: (1) The Landsat-8 satellite can provide an effective data source for mapping PMF in the northern part of China.The highest accuracy achieved was 97.01% in Jizhou and 90.67% in Guyuan, respectively.(2) In terms of features sets, the spectral features, indices features of NDVI and GI, as well as the textural features of mean are more important than the other features for PMF mapping.
The thermal features have a limited contribution for mapping PMF.The features optimized by RF performed better than the spectral or textural features alone in two different study areas.(3) In terms of the optimal phase for mapping PMF, the features from April provide the highest accuracy, followed by the features from May.The multi-temporal combined features can provide better results than single-temporal features.Among the multi-temporal combinations, the two-temporal combination of April and May provide the best accuracy and a more acceptable spatial pattern of PMF.
(4) In terms of classifiers, the RF performs significantly better than SVM with regards to the mapping accuracy and the efficiency both in sin If their contribution is not too much gle and multi-temporal mapping.
Our study indicated that the presented technique scheme and satellite data can be used to map the spatial distribution of PMF, which will provide basic data for the further study of social and eco-environmental effects of PMF.

Figure 1 .
Figure 1.The statistical results of plastic-mulched farmland (PMF) coverage in China.

Figure 1 .
Figure 1.The statistical results of plastic-mulched farmland (PMF) coverage in China.

Figure 5 .
Figure 5. Feature importance of each period in Jizhou (the abscissa corresponds with the order of the features in Table2).

Figure 5 .
Figure 5. Feature importance of each period in Jizhou (the abscissa corresponds with the order of the features in Table2).

Figure 7 .
Figure 7. Spatial distribution of PMF in Jizhou using RF feature selection and RF classifier.The first line represents the results from spectral reflectance features, the second line represents the result from the textural features, the third line represents the results from the total features, the fourth line represents the results from the optimized features, and the fifth line represents the results from optimal temporal combined features.(a) 2015.04.16;(b) 2015.05.18;(c) 2015.06.03;(d) 2015.06.19;(e) 2015.07.05; (x) two-temporal combination of 2015.04.16 and 2015.05.18; (y) three-temporal combination of 2015.04.16, 2015.05.18, and 2015.06.03; and (z) four-temporal combination of 2015.04.16, 2015.05.18, 2015.06.03, and 2015.06.19.

Figure 7 .
Figure 7. Spatial distribution of PMF in Jizhou using RF feature selection and RF classifier.The first line represents the results from spectral reflectance features, the second line represents the result from the textural features, the third line represents the results from the total features, the fourth line represents the results from the optimized features, and the fifth line represents the results from optimal temporal combined features.(a) 2015.04.16;(b) 2015.05.18;(c) 2015.06.03;(d) 2015.06.19;(e) 2015.07.05; (x) two-temporal combination of 2015.04.16 and 2015.05.18; (y) three-temporal combination of 2015.04.16, 2015.05.18, and 2015.06.03; and (z) four-temporal combination of 2015.04.16, 2015.05.18, 2015.06.03, and 2015.06.19.

Figure 8 .
Figure 8. Spatial distribution of PMF in Jizhou using the SVM feature selection and SVM classifier.The first line represents the results from spectral reflectance features, the second line represents the result from the textural features, the third line represents the results from the total features, the fourth line represents the results from the optimized features, and the fifth line represents the results from optimal temporal combined features.(a) 2015.04.16;(b) 2015.05.18;(c) 2015.06.03;(d) 2015.06.19;(e) 2015.07.05; (x) two-temporal combination of 2015.04.16 and 2015.05.18; (y) three-temporal combination of 2015.04.16, 2015.05.18, and 2015.06.03; and (z) four-temporal combination of 2015.04.16, 2015.05.18, 2015.06.03, and 2015.06.19.

Figure 8 .
Figure 8. Spatial distribution of PMF in Jizhou using the SVM feature selection and SVM classifier.The first line represents the results from spectral reflectance features, the second line represents the result from the textural features, the third line represents the results from the total features, the fourth line represents the results from the optimized features, and the fifth line represents the results from optimal temporal combined features.(a) 2015.04.16;(b) 2015.05.18;(c) 2015.06.03;(d) 2015.06.19;(e) 2015.07.05; (x) two-temporal combination of 2015.04.16 and 2015.05.18; (y) three-temporal combination of 2015.04.16, 2015.05.18, and 2015.06.03; and (z) four-temporal combination of 2015.04.16, 2015.05.18, 2015.06.03, and 2015.06.19.

Figure 9 .
Figure 9. Spatial distribution of PMF in Guyuan using the RF feature selection and RF classifier.The first line represents the results from spectral reflectance features, the second line represents the result from the textural features, the third line represents the results from the total features, the fourth line represents the results from the optimized features, and the fifth line represents the results from optimal temporal combined features.(a) 2015.04.26;(b) 2015.05.12;(c) 2015.06.13;(d) 2015.07.15;(e) 2015.07.31; (x) two-temporal combination of 2015.05.12 and 2015.06.13; (y) three-temporal combination of 2015.05.12, 2015.06.13, and 2015.07.15; and (z) four-temporal combination of 2015.04.26, 2015.05.12, 2015.06.13, and 2015.07.15.

Figure 9 .
Figure 9. Spatial distribution of PMF in Guyuan using the RF feature selection and RF classifier.The first line represents the results from spectral reflectance features, the second line represents the result from the textural features, the third line represents the results from the total features, the fourth line represents the results from the optimized features, and the fifth line represents the results from optimal temporal combined features.(a) 2015.04.26;(b) 2015.05.12;(c) 2015.06.13;(d) 2015.07.15;(e) 2015.07.31; (x) two-temporal combination of 2015.05.12 and 2015.06.13; (y) three-temporal combination of 2015.05.12, 2015.06.13, and 2015.07.15; and (z) four-temporal combination of 2015.04.26, 2015.05.12, 2015.06.13, and 2015.07.15.

Table 1 .
The land cover classification scheme in Jizhou and Guyuan.

Table 1 .
The land cover classification scheme in Jizhou and Guyuan.

Table 2 .
The extracted features from a single temporal Landsat-8 image.NDBI: normalized difference built-up index; NDVI: normalized difference vegetation index.

Table 4 .
Top 20 selected features from RF (random forest) for mapping PMF in different parts of the growing season.

Table 4 .
Top 20 selected features from RF (random forest) for mapping PMF in different parts of the growing season.

Table 5 .
The mapping accuracies generated from RF using different features in Jizhou.OA: overall accuracy; PA: producer's accuracy; UA: user's accuracy.

Table 6 .
Z-test values for the pairwise comparison of the error matrices of single temporal mapping.

Table 7 .
The mapping accuracies of PMF in Jizhou using combined multi-temporal features.

Table 8 .
Z-test values for the pairwise comparison of the error matrices of multi-temporal mapping.

Table 9 .
Confusion matrix of classification with the highest and the poorest accuracies.WB: water body; VC: vegetation cover; PMF: plastic-mulched farmland; BS: bare soil; IS: impervious surface.

Table 10 .
The mapping accuracies generated from RF using different features in Guyuan.

Table 11 .
The mapping accuracy generated from RF using the combined multi-temporal features in Guyuan.

Table 12 .
Z-test values for the pairwise comparison of the error matrices.