Monitoring Cropland Dynamics of the Yellow River Delta based on Multi-Temporal Landsat Imagery over 1986 to 2015

Abstract

Natural deltas can provide human beings with flat and fertile land to be cultivated. It is important to monitor cropland dynamics to provide policy-relevant information for regional sustainable development. This paper utilized Landsat imagery to study the cropland dynamics of the Yellow River Delta during the last three decades. Multi-temporal Landsat data were used to account for the phenological variations of different plants. Several spectral and textural features were adopted to increase the between-class separability. The robust random forest classifier was used to generate the land cover maps of the Yellow River Delta for 1986, 1995, 2005 and 2015. Experimental results indicated that the proposed methodology showed good performance with an average classification accuracy of 89.44%. The spatial-temporal analysis indicated that the cropland area increased from 467.6 km² in 1986 to 718.5 km² in 2015 with an average growth rate of 8.65 km²/year. The newly created croplands were mainly due to the reclamation of grassland and bare soil while the losses of croplands were due to abandoned cultivation and urban sprawl. The results demonstrate that a sustainable perspective should be adopted by the decision makers in order to simultaneously maintain food security, industrial development and ecosystem safety.

1. Introduction

Natural river deltas are home to many valuable wetland ecosystems around the globe, which provide a variety of functions beneficial to the sustainability of human communities, including flooding control, water quality protection, and biodiversity conservation, to name a few. [1,2,3,4,5]. Meanwhile, the flat topography, sufficient water supply and fertile arable land have made river deltas preferential areas for human settlements and intensified agriculture practices [3,4,5], resulting in plenty of croplands co-existing with coastal wetlands. However, due to intensive anthropogenic activities such as industrialization and urban sprawl, the croplands in natural deltas have experienced great dynamics [6,7,8,9,10], e.g., some croplands were converted to built up areas while some wetlands were reclaimed as the cultivated land in some natural deltas of China. Therefore, it is of great significance to monitor cropland dynamics appropriately and frequently to assist in the sustainable development of natural river deltas.

Considering the dynamic nature of river deltas [3,4,5], it is time-consuming and impractical to rely on traditional field surveys to monitor the land use and land cover change (LUCC). Nevertheless, satellite remote sensing has yielded promising results in mapping and monitoring river deltas due to its synoptic view, multi-temporal coverage, multi-sensor data collection and cost-effectiveness [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]. The Landsat program has been considered as one of the most successful and important missions for Earth Observation in the last forty years [5]. Archived Landsat data are free of charge online [17] and have demonstrated great capabilities in resource investigation, LUCC, disaster monitoring, urban heat island and global change research [11,12,13,14,15,16,18,19,20]. In this study, multi-temporal Landsat imagery was utilized to monitor the cropland dynamics of the Yellow River Delta of China.

Yellow River Delta is China’s largest natural delta and one of the fastest growing deltas in the world [21]. It is located at the northeast of Shandong Province and is home to many valuable coastal wetlands together with China’s second largest oil field, Shengli Oil Field. Due to intensive human activities and global climate change, the Yellow River Delta has been experiencing many challenges, including sea level rise, water discharge decrease, soil salinization and fertile decline, wetland loss and fragmentation, to name a few. [5,21,22,23,24,25,26,27,28], which calls for sustainable actions to achieve the balance between socio-economic development and ecosystem protection. Previous studies [5,21,22,23,24,25,26,27,28] have been conducted on the LUCC of the Yellow River Delta to investigate its spatial-temporal patterns. Yue et al. [22] researched the landscape change detection of the newly created wetlands in the Yellow River Delta using Landsat Thematic Mapper (TM) data of 1984, 1991 and 1996. Significant landscape changes were reported and several measures were suggested to protect the newly created wetlands. Fang et al. [26] studied the land cover and vegetation composition changes of the Yellow River Delta Nature Reserve (YRDNR) based on the Landsat TM data of 1995 and found that the YRDNR has experienced great landscape changes and environmental improvement after 1992 [26]. Liu et al. [27] monitored the land cover changes also in the YRDNR from 1999 to 2008 based on remote sensing, geographical information system (GIS) and Markov chain transition probability matrices. Ottinger et al. [5] utilized Landsat TM data of 1995, 2004 and 2010 to investigate the land cover changes of the entire Yellow River Delta (Dongying City), which showed that land cover changes were mainly caused by intensified farming and urban sprawl. There are several studies whose research area is also the Yellow River Delta; however, their focuses were not on land use and land cover changes. Therefore, they are not discussed in detail in this introduction [21,23,24,25,28].

However, previous studies [5,21,22,23,24,25,26,27,28] mainly discussed the general geomorphic features and changes of the Yellow River Delta, very few studies used Earth Observation techniques to focus on and research into the cropland dynamics and the driving forces behind [29]. The only exception is Zhao’s study [29], which utilized Landsat TM data to detect the cultivated land dynamics of Kenli County. Considering that Zhao only studied a short period from 1987 to 1998 and only part of the whole delta, it is necessary to extend Zhao’s study to a longer period while covering the whole Yellow River Delta at the same time.

Inspired by Zhao et al. [29], we studied the cropland dynamics of the Yellow River Delta from 1986 to 2015 (present) using satellite remote sensing techniques. Meanwhile, the previous studies [5,21,22,23,24,25,26,27,28] only used single-date remote sensing imagery, which could not well account for the seasonal and phenological variations of wetland ecosystems, resulting in the classification errors between cropland and other vegetated land cover types such as grassland and forest. Considering that multi-temporal remote sensing data provide the potential to improve classification accuracy when compared to single-date classification [30,31,32,33], we utilized multi-temporal Landsat imageries to aid in the discrimination between different land cover types. Moreover, given the landscape complexity of the delta, the classifier involved should be robust, efficient and insensitive to outliers and require no assumptions of data distribution [33]. The random forest (RF) classifier [33], which satisfies the above requirements, was selected in this study for the classification of the multi-temporal Landsat data. Additionally, random forest has been rarely documented in the LUCC of the Yellow River Delta and other ecological studies (e.g., forest biomass estimation, chlorophyll-a concentration of the water body), and its performance should be justified.

Overall, the objective of this paper was to monitor the cropland dynamics of the Yellow River Delta using multi-temporal Landsat imagery. Specifically, this paper aimed to (i) study the spatial-temporal dynamics of cropland to provide scientific basis for the sustainable development of the Yellow River Delta; (ii) justify whether the integration of multi-temporal images can increase the classification accuracy; (iii) verify the performance of the random forest classifier in the land cover mapping of the Yellow River Delta.

2. Study Area and Data

The Yellow River Delta is located at the estuary of the Yellow River, which is the second longest river in China with a total length of 5,464 km (Figure 1). In this study, the Yellow River Delta refers to the modern Yellow River Delta which has been formed since 1855, when the river’s lower course translocated the river towards the coast of the Bohai Sea [5]. The modern Yellow River Delta starts from Yuwa village and extends to the mouth of Tiaohe River in the northwest and Song Chunrong Channel in the southeast [34].

Due to the deposition of large amount of sediments transported by the Yellow River, the newly created wetland of the Yellow River Delta increases by about 30 km² each year, making it one of the fastest sedimentation areas around the world [21].

Figure 1. Study area. (a) China; (b) True color image of the Yellow River Delta on 5 June 2015.

Figure 1. Study area. (a) China; (b) True color image of the Yellow River Delta on 5 June 2015.

The Yellow River Delta belongs to the temperate continental monsoon climate and has four distinct seasons [21]. It enjoys a hot and humid summer and a cold and dry winter. The annual temperature is about 11.9 °C while the annual average precipitation is about 640 mm [21]. The major crops include winter wheat, rice, corn, cotton and soybeans while the main natural vegetation includes reed, tamarisk, seablite and robinia [5,21].

The Yellow River Delta is also home to Shengli Oil Field which is the basis of Dongying City establishment since 1983 [5]. The intensive oil related industrial activities have caused rapid urban expansion and population growth [5], resulting in the conversion of plenty of wetlands, croplands and bare soil into built up areas including oil fields, chemical factories, harbors and residential areas.

The multi-temporal Landsat-5 TM data of 1986, 1995, 2005 and Landsat-8 OLI (Operational Land Imager) data of 2015 were utilized in this study for land cover classification of the Yellow River Delta (Table 1). In fact, the data used belong to the optical multispectral bands of Landsat 5 TM and Landsat 8 OLI. Six of the seven multispectral bands of OLI are nearly consistent with TM, providing for the compatibility with the historical Landsat data. However, the spectral response and radiometry of TM and OLI are not exactly the same, which may introduce errors and differences when using the data of the two sensors at the same time. All data used were acquired under almost cloud free conditions and were provided by the U.S. Geographical Survey (USGS) [17]. In addition, TM and OLI sensors have a spatial resolution of 30 m while the Landsat 8 platform has a revisit time of 16 days. Meanwhile, the actual revisit of Landsat 5 is not 16 days since the acquisition is not systematic [17]. Given that the incorporation of multi-temporal (or multi-seasonal) imagery can provide additional phenological information for differentiating plant species within a single growing season [30,31,32,33], we also employed multi-temporal Landsat data to increase the classification accuracy.

Table 1. Description of Landsat data used in this study.

Year	Platform	Sensor	Date	Path/row	Resolution	Band number
1986	Landsat-5	TM	1986.04.18	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	1986.06.05	121/34	30 m	1, 2, 3, 4, 5, 7
1995	Landsat-5	TM	1995.03.26	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	1995.09.18	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	1995.10.04	121/34	30 m	1, 2, 3, 4, 5, 7
2005	Landsat-5	TM	2005.04.22	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	2005.05.24	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	2005.10.15	121/34	30 m	1, 2, 3, 4, 5, 7
2015	Landsat-8	OLI	2015.03.01	121/34	30 m	1, 2, 3, 4, 5, 6, 7
	Landsat-8	OLI	2015.05.04	121/34	30 m	1, 2, 3, 4, 5, 6, 7
	Landsat-8	OLI	2015.06.05	121/34	30 m	1, 2, 3, 4, 5, 6, 7

Table 1. Description of Landsat data used in this study.

Year	Platform	Sensor	Date	Path/row	Resolution	Band number
1986	Landsat-5	TM	1986.04.18	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	1986.06.05	121/34	30 m	1, 2, 3, 4, 5, 7
1995	Landsat-5	TM	1995.03.26	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	1995.09.18	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	1995.10.04	121/34	30 m	1, 2, 3, 4, 5, 7
2005	Landsat-5	TM	2005.04.22	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	2005.05.24	121/34	30 m	1, 2, 3, 4, 5, 7
	Landsat-5	TM	2005.10.15	121/34	30 m	1, 2, 3, 4, 5, 7
2015	Landsat-8	OLI	2015.03.01	121/34	30 m	1, 2, 3, 4, 5, 6, 7
	Landsat-8	OLI	2015.05.04	121/34	30 m	1, 2, 3, 4, 5, 6, 7
	Landsat-8	OLI	2015.06.05	121/34	30 m	1, 2, 3, 4, 5, 6, 7

3. Methods

3.1. Workflow

The whole workflow of this study is depicted in Figure 2.

Figure 2. The workflow of this study.

Figure 2. The workflow of this study.

The workflow mainly consists of the following steps: (i) remote sensing data preprocessing (Section 3.2), which includes radiometric calibration, geometric rectification and atmospheric correction of the downloaded Landsat data; (ii) feature calculation (Section 3.3), which includes the Tasseled Cap Transformation (TCT) [35,36,37,38], Normalized Difference Vegetation Index (NDVI) [39], Soil Adjusted Vegetation Index (SAVI) [39], Modified Normalized Difference Water Index (MNDWI) [40] and six texture features derived from the gray level co-occurrence matrix (GLCM) [41,42,43]; (iii) random forest classification (Section 3.5), which utilized the RF to classify the multi-temporal Landsat data into land use and land cover (LULC) maps. Additionally, accuracy assessment (Section 3.6) was carried out to justify the performance of the proposed approach in this study. Finally, spatial-temporal pattern analysis was utilized to study the cropland dynamics of the Yellow River Delta based on the LULC maps.

3.2. Image Preprocessing

Remote sensing image preprocessing was utilized to acquire the ground surface reflectance through several steps, including radiometric calibration, geometric and atmospheric correction. Therein, radiometric calibration was to convert the digital number (DN) values of the original downloaded Landsat data into at-satellite reflectance. The following formula was utilized for the radiometric calibration of Landsat data.

$$Ref=Gain\times DN+offset$$

(1)

where Ref refers to the reflectance, DN refers to the digital number, and Gain and offset are calibration coefficients which can be derived from the header file of the Landsat data. Additionally, images acquired on different dates should be co-registered geographically to maintain the positional accuracy. The nearest neighbor resampling method was adopted in geometric correction to preserve the original spectral characteristics of each pixel. As for atmospheric correction, the by-band 6S model [44] was adopted in this study to remove the atmospheric effects from the remotely sensed data.

3.3. Feature Calculation

The integration of different image features can improve the classification results when compared to the original remote sensing imagery, according to many previous studies [11,12,30,31,32,33]. The image features were calculated for each Landsat image, consisting of spectral features (i.e., TCT components, NDVI, SAVI and MNDWI) and texture features (i.e., the six least correlated features derived from GLCM), which will be documented in detail below.

3.3.1. Tasseled Cap Transformation

Tasseled Cap Transformation was first introduced by Kauth and Thomas (1976) [35] and further developed by Crist and Cicone (1984) [36]. It transformed the original multispectral data orthogonally into a new set of axes associated with physical meanings [37]. Generally, three distinctive axes of brightness, greenness and wetness can be derived through the orthogonal rotation coefficients [37].

TCT features have long been recognized as an efficient tool in understanding the vegetation growth patterns, estimates of forest structure, providing ancillary information for improving the LULC classification, to name a few. [37,38]. In this study, TCT features were employed to increase the between-class separability, such as the differentiation between cropland and other vegetated land cover types. The TCT coefficients for Landsat-5 TM and Landsat-8 OLI were listed in Table 2 according to previous studies [37,38].

Table 2. TCT coefficients for Landsat-5 TM and Landsat-8 OLI at-satellite reflectance.

Landsat-5 TM	Band 1	Band 2	Band 3	Band 4	Band 5	Band 7
Brightness	0.2909	0.2493	0.4806	0.5568	0.4438	0.1706
Greenness	−0.2728	−0.2174	−0.5508	0.7221	0.0733	−0.1648
Wetness	0.1446	0.1761	0.3322	0.3396	−0.6210	−0.4186
Landsat-8 OLI	Band 2	Band 3	Band 4	Band 5	Band 6	Band 7
Brightness	0.3029	0.2786	0.4733	0.5599	0.5080	0.1872
Greenness	−0.2941	−0.2430	−0.5424	0.7276	0.0713	−0.1608
Wetness	0.1511	0.1973	0.3283	0.3407	−0.7117	−0.4559

Table 2. TCT coefficients for Landsat-5 TM and Landsat-8 OLI at-satellite reflectance.

Landsat-5 TM	Band 1	Band 2	Band 3	Band 4	Band 5	Band 7
Brightness	0.2909	0.2493	0.4806	0.5568	0.4438	0.1706
Greenness	−0.2728	−0.2174	−0.5508	0.7221	0.0733	−0.1648
Wetness	0.1446	0.1761	0.3322	0.3396	−0.6210	−0.4186
Landsat-8 OLI	Band 2	Band 3	Band 4	Band 5	Band 6	Band 7
Brightness	0.3029	0.2786	0.4733	0.5599	0.5080	0.1872
Greenness	−0.2941	−0.2430	−0.5424	0.7276	0.0713	−0.1608
Wetness	0.1511	0.1973	0.3283	0.3407	−0.7117	−0.4559

3.3.2. Spectral Indices

Apart from TCT features, several spectral indices of NDVI, SAVI and MNDWI were also included due to their prevalence in remote sensing fields and the capability in describing the characteristics of the Earth’s surfaces [39,40]. These spectral indices are presented as follows.

$$\text{NDVI}=\frac{\rho (\text{NIR})-\rho (\text{Red})}{\rho (\text{NIR})+\rho (\text{Red})}$$

(2)

$$\text{SAVI}=\frac{\rho (\text{NIR})-\rho (\text{Red})}{\rho (\text{NIR})+\rho (\text{Red})+L}\times (1+L)$$

(3)

$$\text{MNDWI}=\frac{\rho (\text{Green})-\rho (\text{MIR})}{\rho (\text{Green})+\rho (\text{MIR})}$$

(4)

where ρ(Green), ρ(Red), ρ(NIR) and ρ(MIR) refer to the reflectance of the green, red, near-infrared and short wave infrared band of Landsat data and L stands for the correction factor [39]. NDVI was used to separate the vegetated from non-vegetated land cover types while SAVI could improve NDVI’s performance in sparsely vegetated areas of the Yellow River Delta. MNDWI refers to the modified version of NDWI [45], which was proposed by Xu (2006) [40] for Landsat-5 TM data. Although NDWI has been widely used in the remote sensing field, it cannot suppress the signal from building up efficiently and the extracted water features are still mixed with built up land noise [40]. To tackle this issue, MNDWI was proposed by substituting the MIR band for the NIR band of NDWI (Formula (4)). A series of experiments were carried out by Xu [40] for lake water, sea water and river water enhancement and the results indicated that MNDWI can enhance open water features effectively and suppress or even remove the built up, vegetation and soil noise at the same time [40]. MNDWI has also been utilized in other studies, e.g., the water body extraction in the coastal zone east of the Nile Delta of Egypt [46] and the coastline detection of Qinhuangdao of China [47].

3.3.3. Texture Features

According to visual interpretation and field survey, we have found that many land cover types of the Yellow River Delta manifest several distinctive geometric and textural characteristics, e.g., most croplands, reservoirs and cultivated aquatic surfaces have obvious rectangular boundaries, while reeds and seablites show long and narrow distribution patterns along the river and tidal flat, respectively. Thus, the inclusion of texture features can help improve the classification results. Additionally, the role of texture in increasing classification accuracy has been recognized in many studies [41,42,43]. Laliberte et al. [48] analyzed texture measures at multiple scales in object-based analysis to differentiate broad functional groups of vegetation in arid grassland with subdecimeter unmanned aerial vehicle (UAV) imagery. The results indicated that the inclusion of textures can help improve the classification accuracy and overcome the drawbacks of low-cost off-the-shelf digital cameras [48]. Similar findings can also be seen in the studies of Feng et al. [49,50], where textures were added to improve the classification accuracy in urban vegetation mapping and urban flood mapping based on UAV optical imagery. Additionally, Szantoi [42] also adopted texture features to improve the accuracy of analyzing fine-scale wetland composition using areal imagery. The six least correlated texture measures according to Szantoi [42] were used in this study, i.e., mean (MEA), standard deviation (STD), homogeneity (HOM), dissimilarity (DIS), entropy (ENT) and angular second moment (ASM) are as follows:

$${\rm M}{\rm E}{\rm A}i={\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{j=0}^{N-1}i\times P(i,j)}}$$

(5)

$$\mathrm{STD}i=\sqrt{{\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{j=0}^{N-1}P(i,j)\times {(i-\mathrm{MEA}i)}^2}}}$$

(6)

$$\mathrm{HOM}={\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{j=0}^{N-1}\frac{P(i,j)}{1+{(i-j)}^2}}}$$

(7)

$$\mathrm{DIS}={\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{j=0}^{N-1}P(i,j)\times |i-j|}}$$

(8)

$$\mathrm{ENT}={\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{j=0}^{N-1}-P(i,j)\times \ln (P(i,j))}}$$

(9)

$$\mathrm{ASM}={\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{j=0}^{N-1}P{(i,j)}^2}}$$

(10)

where N is the number of grey levels; P is the normalized symmetric GLCM of dimension N × N; P(i, j) is the normalized grey level value in the cell i, j of the co-occurrence matrix, such that the sum of P(i, j) is equal to 1 [42]. In order to determine the optimal size of the moving window within which the texture features were calculated, six different moving windows were chosen for Landsat-8 OLI imagery (2015-06-05) in this study: 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11 and 13 × 13. Texture features derived at different window sizes were added separately as additional ancillary bands to the original multispectral bands for classification. Results indicated that textures calculated at a 5 × 5 moving window yielded the highest accuracy.

3.4. Land Cover Classification Schemes

After feature calculation, all of the image features and the original multi-temporal Landsat bands were integrated into a single dataset, from which the training samples were derived to train the random forest classifier. According to field survey and previous studies [5,21,22,23,24,25,26,27,28], the land cover of this study included eight categories, i.e., cropland, forest, grassland, shrubs, water, tidal flat, built up area and bare soil. The landscape descriptions of the classification scheme are shown in Table 3.

Table 3. Classification scheme of the Yellow River Delta.

Land cover class	Description
Cropland	Including irrigated and non-irrigated farmland, e.g. winter wheat, rice, corn, cotton and soybean
Forest	Broadleaved trees, mainly robinia, willow
Grassland	Including densely vegetated areas where reed is dominant and sparsely vegetated areas where seablite is dominant
Shrubs	Sparsely vegetated shrubs, mainly tamarisks
Water	Water bodies including rivers, reservoirs, aquaculture and brine ponds
Tidal flat	Non-vegetated foreshore areas
Built up	Artificial surfaces including residential areas, factories, harbors and oil fields
Bare soil	Non-vegetated bare land, mainly saline and alkaline land

Table 3. Classification scheme of the Yellow River Delta.

Land cover class	Description
Cropland	Including irrigated and non-irrigated farmland, e.g. winter wheat, rice, corn, cotton and soybean
Forest	Broadleaved trees, mainly robinia, willow
Grassland	Including densely vegetated areas where reed is dominant and sparsely vegetated areas where seablite is dominant
Shrubs	Sparsely vegetated shrubs, mainly tamarisks
Water	Water bodies including rivers, reservoirs, aquaculture and brine ponds
Tidal flat	Non-vegetated foreshore areas
Built up	Artificial surfaces including residential areas, factories, harbors and oil fields
Bare soil	Non-vegetated bare land, mainly saline and alkaline land

Training samples of each land cover type were selected randomly using a stratified sampling method in a small polygon block. It was assumed that the pixels within the polygon belong to the same land cover class. The number of training samples for each class was equally set to be 500 to ensure sufficient observations and to avoid any under- or overestimation of the spectral patterns of different land cover types. Specifically, before the stratified sampling process, a number of 600 sampling points were selected for each land cover type. During the sampling process, 500 sampling points out of 600 were randomly selected as the training dataset for each land cover type. The remaining 100 points were selected as the testing dataset for accuracy assessment.

3.5. The Random Forest Classifier

The random forest classifier was utilized to classify the input data into thematic land cover types in this study. As an ensemble learning algorithm, random forest was proposed by Breiman in 2001 [51], which can be regarded as an ensemble of decision trees. The result yielded by random forest is determined by the output of all the decision trees involved [51]. As an ensemble learning method, the random forest classifier has several advantages including that it is easy to parameterize, non-sensitive to over-fitting of the training data and good at dealing with outliers [33,49]. Random forest requires no assumption of training data distribution when compared to statistical methods such as the maximum likelihood classifier (MLC) [49,50]. When compared to decision tree, over-fitting is less of an issue for random forest and there is no need for the cumbersome task of pruning the trees [51]. When compared to support vector machine (SVM), random forest has the advantage of easier parameterization and better generalization capability [49,50]. Due to the above advantages, random forest has been increasingly used for image classification in the remote sensing field [33,49,50,52,53,54]. Schneider [54] utilized the random forest classifier (boosted decision trees) to monitor land cover changes in urban and peri-urban areas based on multi-seasonal Landsat satellite data. Results indicated that both random forest and SVM outperformed MLC but random forest was superior to SVM at handling missing data [54]. Corcoran et al. [33] used the random forest classifier for wetland mapping in Northern Minnesota and analyzed the influence of multi-source and multi-temporal remotely sensed and ancillary data on the classification accuracy. Feng et al. [49,50] also justified the role of random forest in urban vegetation mapping and urban flood mapping based on very high resolution UAV imagery. However, random forest has rarely been applied to the land cover classification of the Yellow River Delta and we are motivated to justify its performance.

Two random selection steps are used to increase the generalization capability of random forest when generating each decision tree [51]. First, only about 2/3 of the training samples are selected with replacement to build every decision tree [51]. The left 1/3 training samples are called out-of-bag (OOB) data [51], which can be utilized for the inner cross validation of random forest. Generally, the plot of OOB error versus the number of trees (ntree) is needed to justify whether there exist sufficient trees in the grown forest [52]. Second, only a subset of the predictable variables is randomly selected to split every decision tree [51]. The number of picked up variables (mty) is set to be 1/3 or the square root of the number of all input variables [51].

The package randomForest in R language [55,56] was adopted for the parameterization of the random forest classifier used in this study. ntree was set to be 200 according to our previous studies [49,50] and the plot of OOB error of the 2015 Landsat dataset versus ntree is depicted in Figure 3 to verify whether 200 decision trees can make the random forest classifier convergent.

Figure 3. OOB error of random forest.

Figure 3. OOB error of random forest.

Figure 3 depicts that the OOB error decreased quickly from 10% to less than 2% when ntree increased from 0 to 50. Afterwards, the OOB error continued to decline slightly and converged to a stable state (1.4%) when ntree reached 200. Therefore, 200 decision trees proved to be feasible for the random forest classifier.

3.6. Accuracy Assessment

In order to assess the performance of the proposed classification procedure, both confusion matrix [49] and visual evaluation were employed. Visual evaluation was used to check the visual effects of the classification results to see whether there existed obvious omission and commission errors. Confusion matrix was derived from reference samples to quantitatively assess the quality of the generated LULC maps through the following indicators: overall accuracy, producer accuracy, user accuracy and Kappa index [49]. Testing points were selected in the stratified sampling process mentioned at the end of Section 3.4. A number of 100 points were selected for each land cover type as the testing dataset for accuracy assessment. The overall accuracy was then weighted by the sampling probability of each class. The reason why proportional random sampling is not used is that the proportion of each land cover type remains unknown before the classification, therefore, we only utilized the stratified sampled testing points for accuracy assessment. Additionally, cross validation was used to determine the meta-parameters of the random forest classifier; however, it was not enough to characterize the results of the “optimal” classification. In remote sensing field, it is mandatory to carry out the accuracy assessment with an independent validation sample dataset.

4. Results and Discussion

4.1. Results of Land Cover Classification

The spatial distributions of cropland of the Yellow River Delta together with other land cover types are illustrated in Figure 4.

Figure 4. LULC maps of the Yellow River Delta in (a) 1986; (b) 1995; (c) 2005; (d) 2015.

Figure 4. LULC maps of the Yellow River Delta in (a) 1986; (b) 1995; (c) 2005; (d) 2015.

It can be observed that the cropland showed distinct spatial patterns. Most cultivated lands were distributed along the corridors of the present and the old courses of the Yellow River. This is mainly due to the high soil moisture content and low soil salinity in those areas [29], which can provide suitable and arable lands for agriculture. Specifically, in 1986, the cropland landscape was rather fragmented and large amount of croplands were surrounded by bare soil and grassland. The year 1995 witnessed a similar cropland distribution pattern to that of 1986. However, the year 2005 manifested a great increase of cropland with numerous large patches. Lots of bare soil and grassland were converted to cropland during the period of 1995 to 2005. Finally, the cropland distribution of 2015 showed a similar spatial pattern to that of 2005.

Besides, shrubs (mainly tamariks) were mainly distributed along the coastal regions with high soil salinity. Grasslands (mainly reed), as an important component of the wetland ecosystem, were distributed between shrubs and cropland. The water bodies were mainly located in the coastal regions where aquaculture and brine ponds were prevalent. Additionally, the water bodies of the inland regions were scattered including reservoirs and small fish ponds.

4.2. Temporal Changes of Cropland and Other Land Cover Types

The temporal changes of cropland in the Yellow River Delta are depicted in Figure 5. During the past three decades, croplands have increased from 467.6 km² in 1986 to 718.5 km² in 2015 with an average growth rate of 8.65 km²/year. During 1986 to 1995, the cropland percentage slightly decreased from 17.4% to 16.6% at a descent rate of 1.27 km²/year. However, the percentage increased dramatically to 24.8% in the year 2005 with a growth rate of 22.38 km²/year. The next decade witnessed a slight increase of the cropland area while the percentage grew to 26% from 2005 to 2015 at a relatively lower rate of 3.85 km²/year.

Figure 5. Temporal changes of cropland in the Yellow River Delta.

Figure 5. Temporal changes of cropland in the Yellow River Delta.

In order to get a better understanding of the LULC dynamics of the whole Yellow River Delta, the overall LULC changes from 1986 to 2005 were summarized in Table 4. It indicates that cropland, grassland, bare soil and tidal flat were the four dominant land cover types, totally occupying 81.4%, 74.8%, 78.1% and 65.6% for the years 1986, 1995, 2005 and 2015, respectively. Therein, the grassland percentage dropped from 26.7% in 1986 to 16.6% in 2015 with a descent of 260.1 km². The bare soil experienced the same trend and decreased from 22.6% to 9.2% during the past three decades. However, the tidal flat did not change much while the built up areas increased from 52.1 km² in 1986 to 221.3 km² in 2015.

Table 4. Statistics of LULC areas of the Yellow River Delta from 1986 to 2015.

Class	1986		1995		2005		2015
	Area (km2)	%	Area (km2)	%	Area (km2)	%	Area (km2)	%
Cropland	467.6	17.4	456.2	16.6	680.0	24.8	718.5	26.0
Forest	38.8	1.4	58.8	2.1	39.9	1.5	52.8	1.9
Grass	717.5	26.7	727.1	26.5	525.9	19.2	457.4	16.6
Shrub	264.5	9.8	333.4	12.1	169.6	6.2	152.7	5.5
Built up	52.1	1.9	92.9	3.4	140.2	5.1	221.3	8.0
Tidal flat	394.7	14.7	471.7	17.2	545.2	19.9	381.9	13.8
Bare soil	608.2	22.6	399.3	14.5	390.2	14.2	254.5	9.2
Water	145.3	5.4	208.9	7.6	249.4	9.1	523.1	18.9
Total	2688.7	100	2748.4	100	2740.3	100	2762.1	100

Table 4. Statistics of LULC areas of the Yellow River Delta from 1986 to 2015.

Class	1986		1995		2005		2015
	Area (km2)	%	Area (km2)	%	Area (km2)	%	Area (km2)	%
Cropland	467.6	17.4	456.2	16.6	680.0	24.8	718.5	26.0
Forest	38.8	1.4	58.8	2.1	39.9	1.5	52.8	1.9
Grass	717.5	26.7	727.1	26.5	525.9	19.2	457.4	16.6
Shrub	264.5	9.8	333.4	12.1	169.6	6.2	152.7	5.5
Built up	52.1	1.9	92.9	3.4	140.2	5.1	221.3	8.0
Tidal flat	394.7	14.7	471.7	17.2	545.2	19.9	381.9	13.8
Bare soil	608.2	22.6	399.3	14.5	390.2	14.2	254.5	9.2
Water	145.3	5.4	208.9	7.6	249.4	9.1	523.1	18.9
Total	2688.7	100	2748.4	100	2740.3	100	2762.1	100

Table 4 also indicates that the water bodies have experienced great changes during the past three decades. To be specific, water bodies have increased from 145.3 km² in 1986 to 523.1 km² in 2015 with a mean increasing rate of 13.02 km²/year. The last decade of 2005 to 2015 has witnessed the fastest growing of water bodies at the speed of 22.87 km²/year. Figure 4 indicates that the increase of water bodies mainly happened in the southeastern, northeastern and northwestern coastal regions. Actually, the above mentioned regions consisted of aquaculture and brine ponds whose establishments were due to the market demand and economic stimulus. Besides, a number of reservoirs have been constructed in the inland regions to meet the increasing demand of fresh water supply due to the population growth, which also contributed to the expansion of water areas.

Meanwhile, the built up areas have increased from 52.1 km² in 1986 to 221.3 km² in 2015 with a percentage increase of 324.7%. This is mainly due to the development of oil related industries resulting in the establishments of a number of oil fields in the east and chemical factories and harbors in the north. The growth of industrial workers also contributed to the expansion of residential areas and the urbanization process.

The past three decades have also witnessed a decline of shrubs from 264.5 km² in 1986 to 152.7 km² in 2015. One reason for this decline is the establishment of the dyke along the coastline. The dyke hindered the exchange of mass and energy between the terrestrial and marine ecosystems, which resulted in the habitat degradation of the shrubs and accounted for the decrease of shrub areas. As well, the area of forest is relatively small and has not changed much during the study period.

In addition, statistical analysis should be carried out to determine whether significant differences exist among the four different dates. A one-way analysis of variance (ANOVA) has been done between different dates using SPSS 20 and the significance (Sig.) is described in Table 5.

Table 5. Significance derived from ANOVA between different dates.

	1986 & 1995	1986 & 2005	1986 & 2015	1995 & 2005	1995 & 2015	2005 & 2015
Sig.	0.951	0.958	0.939	0.993	0.988	0.981

Table 5. Significance derived from ANOVA between different dates.

	1986 & 1995	1986 & 2005	1986 & 2015	1995 & 2005	1995 & 2015	2005 & 2015
Sig.	0.951	0.958	0.939	0.993	0.988	0.981

Table 5 indicates that all Sig. have a value more than 0.05 (Sig. > 0.05), which demonstrates that although the Yellow River Delta experienced great land cover changes during the past three decades, significant differences of land cover types between these dates do not exist.

The uncertainties of temporal changes of land cover types should also be discussed here. The uncertainties mainly come from human factors and natural factors. As for the Yellow River Delta, the uncertainties with regard to land cover changes mainly stem from shifts in the land use policy, the negative effects of rapid economic development (e.g., the illegal urban encroachment on cropland) and the impact of global climate change (e.g., the decline of tidal flats due to sea level rise). With regard to the above uncertainties, remote sensing can be recognized as a cost-effective tool in the detection and description of temporal and spatial process of land cover changes due to its synoptic view and continuous coverage of the study area. Specifically, through visual interpretation and image classification, remote sensing can provide high quality land cover datasets with good accuracy to detect land cover changes due to the above uncertainties. For instance, the illegal use of land can be monitored and the coastline changes due to climate change can be mapped through remote sensing to aid for the decision making of the regional sustainable development. However, due to the scarcity of land use statistics of the Yellow River Delta, the comparison between the research results and historical statistical data has not been carried out, which is a drawback of the presented study.

To better understand the relationship between cropland dynamics and other LULC changes, the cropland mutual conversion analysis [14] should be implemented and will be documented in detail in the next section.

4.3. Cropland Mutual Conversion Analysis

Undeniably, LULC change is a mutual conversion process [14], which means there exists a conversion between a certain land cover type and the other land cover types within a certain period. In this study, “cropland mutual conversion” mainly refers to two processes, i.e., the conversion of cropland into other land cover types and the conversion of other land cover types into cropland. Through cropland mutual conversion analysis, better understandings of the changing patterns of cropland and the reasons behind these changes emerge. The cropland conversion to other LULC is demonstrated in Figure 6.

Figure 6. (a) Area of cropland converted to other LULC; (b) Percentage of cropland converted to other LULC during 1986 to 2015.

Figure 6. (a) Area of cropland converted to other LULC; (b) Percentage of cropland converted to other LULC during 1986 to 2015.

Figure 6a depicts that the period of 1986 to 1995 witnessed a total loss of 275.3 km² cropland to other LULC types, which ranks the first in the past three decades. The next decade of 1995 to 2005 experienced the least cropland conversion to other LULC types with an area of 130.7 km². The amount of cropland conversion increased again to 243.3 km² during the last decade. Figure 6b illustrates the percentage of other LULC types to which the cropland had been converted. Grassland, bare soil and built up areas were the three major land cover types which were converted from croplands during the past three decades. Specifically, cropland conversion to grassland was 53.4%, 49.3%, 43.1% in the three study periods, respectively. The percentage of cropland converted to bare soil was 20%, 21.2% and 26.2% while the percentage of cropland converted to built up area was 6.3%, 14.7% and 12.3% during the corresponding periods, respectively. The large conversion of cropland to grassland and bare soil indicated the existence of the phenomenon of abandoned cultivation in the Yellow River Delta. Besides, the urbanization and industrialization of the delta had extended the built up areas to rural regions, which in turn resulted in the loss of cropland simultaneously.

The conversion of other LULC types into cropland during different periods is depicted in Figure 7. From 1986 to 1995, about 263.3 km² of cropland were reclaimed from other LULC types. The area of the newly created cropland was 354.5 km² and 281.8 km² for the period of 1995 to 2005 and 2005 to 2015, respectively. Figure 7b demonstrates that grassland and bare soil were the major contributors for the cropland increase. The percentage of grassland converted to cropland was 38.9%, 67.7% and 40.3% for the three periods, respectively, making grassland the dominant contributor. Meanwhile, the percentage of bare soil converted to cropland in the corresponding periods was 47.5%, 14.1% and 37.7%. The above results indicate that the increase of cropland mainly resulted from the reclamation of grassland and bare soil.

Additionally, many land cover studies assume that once land is converted to built up area or urban, it does not get converted back to cropland. However, Figure 7 indicates that there are a number of built up areas that changed back into cropland in the Yellow River Delta during the study periods. Actually, two reasons may account for this phenomenon. First, due to rapid urbanization, several villages were pulled down while high-rise buildings were established in the adjacent areas. After the residents moved to the new buildings, the original villages were transformed into cropland. Second, due to the spectral mixture between built up areas and bare soil, some bare soil areas have been misclassified as built up areas and those bare soil area might be reclaimed into cropland afterwards.

Figure 7. (a) Area of other LULC converted to cropland; (b) Percentage of other LULC converted to cropland during 1986 to 2015.

Figure 7. (a) Area of other LULC converted to cropland; (b) Percentage of other LULC converted to cropland during 1986 to 2015.

In addition, the accuracy of the change maps is dependent on the accuracy of each individual classification and is subject to the error propagation over time [30,57]. According to Congalton and Green [57], for a map with a chosen accuracy of 90% (10% error) and using a 95% confidence level, the minimum number of samples required for the change/no change assessment is 298. Therefore, we randomly selected 298 samples each for change and no change classes for the change/no change validation for the three periods of this study, respectively. The resulted overall accuracy for the three study periods were then listed as follows.

Table 6 indicates that the overall accuracy of change/no change validation is 73.24%, 78.62% and 80.61%, respectively, for the three study periods, which justifies a high change detection accuracy. The period of 2005 to 2015 has the highest change detection accuracy while the period of 1986 to 1995 has the lowest accuracy.

Table 6. Overall accuracy of change/no change validation.

	From 1986 to 1995	From 1995 to 2005	From 2005 to 2015
OA	73.24%	78.62%	80.61%

Table 6. Overall accuracy of change/no change validation.

	From 1986 to 1995	From 1995 to 2005	From 2005 to 2015
OA	73.24%	78.62%	80.61%

4.4. Results of Accuracy Assessment

To quantitatively assess the accuracy of the proposed approach in this study, confusion matrix was utilized to calculate the overall classification accuracy and Kappa index, which are shown in Table 7.

Table 7. Classification accuracy of the Yellow River Delta from 1986 to 2015.

	1986	1995	2005	2015
OA	85.37%	90.80%	90.34%	91.26%
Kappa	0.8328	0.8948	0.8895	0.9000

Table 7. Classification accuracy of the Yellow River Delta from 1986 to 2015.

	1986	1995	2005	2015
OA	85.37%	90.80%	90.34%	91.26%
Kappa	0.8328	0.8948	0.8895	0.9000

It can be observed that the proposed methodology showed good performance with an average classification accuracy of 89.44% and an average Kappa index of 0.8793. This is mainly due to the adoption of multi-temporal dataset and a series of spectral and textural features, through which can increase the separability between different land cover types [30,31,32,33]. Additionally, the robustness and effectiveness of random forest classifier [49,50,52,53] also contributes to the high accuracy here.

In addition, only the reference data in 2015 were based on field survey. Reference data of other times were indirectly selected from the LULC maps of the previous studies and Google Earth history image, which may cause some uncertainties in the results of the accuracy assessments. However, the adoption of reference data from previous studies can be acceptable in LUCC study [5,14] since it is sometimes very difficult to get complete history ground truth data. Table 7 also indicates that the land cover classification accuracy for 1986 was notably lower than other periods. Three reasons would account for this issue. First, only two images were employed for the classification of 1986 while three images were employed for the other dates. This led to a decrease of eighteen feature bands and would lead to the lower between-class separability when compared to other dates. Second, the accuracy assessment data for 1986 was sparser since the training-calibration dataset was primarily developed for the year 2015. Third, as can be seen form Figure 4, the land cover configuration for 1986 compromised the inter-class distinction which increased the difficulty of differentiating each land cover category with high accuracy.

The confusion matrix of the year 2015 is presented in Table 8 to better understand the commissions and omissions of different land cover types. It indicates that the classification errors mainly occurred between shrubs and bare soil, and between built up areas and bare soil. For the former case, the shrubs of the Yellow River Delta, mainly tamarisks, were sparsely distributed in the coastal wetlands surrounded by saline land [25], which caused the spectral mixture between shrubs and bare soil. For the latter case, many villages were also sparsely distributed among saline land and showed similar spectral signature to that of bare soil [25], which accounted for the classification errors.

Table 8. Confusion matrix of the classification result of the year 2015.

Class	Validation data
	Crop-land	Forest	Grass	Shrub	Built up	Tidal flat	Bare soil	Water	UA (%)
Cropland	94	1	4	0	1	0	2	0	92.16
Forest	0	95	0	0	0	0	0	0	100.00
Grass	6	4	94	9	2	1	0	6	77.05
Shrub	0	0	2	80	0	0	0	0	97.56
Built up	0	0	0	0	70	0	1	0	98.59
Tidal flat	0	0	0	0	0	99	4	0	96.12
Bare soil	0	0	0	11	27	0	93	0	70.99
Water	0	0	0	0	0	0	0	94	100.00
PA (%)	94.00	95.00	94.00	80.00	70.00	99.00	93.00	94.00
OA (%)	89.88	OA-adj	91.26%		Kappa	0.8843	K-adj	0.9000

Notes: OA-adj, adjusted OA based on the sampling frequency; K-adj, adjusted Kappa index.

Table 8. Confusion matrix of the classification result of the year 2015.

Class	Validation data
	Crop-land	Forest	Grass	Shrub	Built up	Tidal flat	Bare soil	Water	UA (%)
Cropland	94	1	4	0	1	0	2	0	92.16
Forest	0	95	0	0	0	0	0	0	100.00
Grass	6	4	94	9	2	1	0	6	77.05
Shrub	0	0	2	80	0	0	0	0	97.56
Built up	0	0	0	0	70	0	1	0	98.59
Tidal flat	0	0	0	0	0	99	4	0	96.12
Bare soil	0	0	0	11	27	0	93	0	70.99
Water	0	0	0	0	0	0	0	94	100.00
PA (%)	94.00	95.00	94.00	80.00	70.00	99.00	93.00	94.00
OA (%)	89.88	OA-adj	91.26%		Kappa	0.8843	K-adj	0.9000

Notes: OA-adj, adjusted OA based on the sampling frequency; K-adj, adjusted Kappa index.

One of the major objectives of this study was to investigate whether the incorporation of multi-temporal Earth Observation data can increase the LULC classification accuracy of the Yellow River Delta. Here we take the year 2015 as an example. The single-date image classifications were performed in the same manner as the multi-temporal classification, also using the random forest classifier containing 200 decision trees. The accuracy comparison between multi-temporal and single-date classification was illustrated in Table 9.

Table 9. Comparison of the accuracy between multi-temporal and single-date classification.

Accuracy	Muti-temporal		Single-date
	2015	2015-03-01	2015-05-04	2015-06-05
OA	91.26%	78.84%	77.45%	76.84%
Kappa	0.9000	0.7581	0.7423	0.7353
	2005	2005-04-22	2005-05-24	2005-10-15
OA	90.34%	69.02%	71.11%	67.12%
Kappa	0.8895	0.6460	0.6698	0.6242
	1995	1995-03-26	1995-09-18	1995-10-04
OA	90.80%	75.91%	79.88%	77.34%
Kappa	0.8948	0.7247	0.7700	0.7410
	1986	1986-04-18	1986-06-05
OA	85.37%	77.56%	78.93%
Kappa	0.8328	0.7435	0.7592

Table 9. Comparison of the accuracy between multi-temporal and single-date classification.

Accuracy	Muti-temporal		Single-date
	2015	2015-03-01	2015-05-04	2015-06-05
OA	91.26%	78.84%	77.45%	76.84%
Kappa	0.9000	0.7581	0.7423	0.7353
	2005	2005-04-22	2005-05-24	2005-10-15
OA	90.34%	69.02%	71.11%	67.12%
Kappa	0.8895	0.6460	0.6698	0.6242
	1995	1995-03-26	1995-09-18	1995-10-04
OA	90.80%	75.91%	79.88%	77.34%
Kappa	0.8948	0.7247	0.7700	0.7410
	1986	1986-04-18	1986-06-05
OA	85.37%	77.56%	78.93%
Kappa	0.8328	0.7435	0.7592

Table 9 indicates that the inclusion of multi-temporal remote sensing data can increase the classification accuracy significantly for all the four years involved. An average OA increase of 17.44%, 30.77%, 16.84% and 9.11% was observed for the years 2015, 2005, 1995 and 1986, respectively, which verified that there was an accuracy gap between multi-temporal and single-date classification in all the four time steps. This is in accordance with several previous studies [30,31,32,33] and it is mainly due to the fact that the incorporation of multi-temporal Earth Observation data can take into consideration the phenological information [30], which accounts for the increase of between-class separability. For instance, the addition of the spring imagery can help differentiate between the winter wheat and other vegetation, since winter wheat was at its growth peak while other plants just begin to turn green [29]. In addition, the summer and autumn imagery can better separate the bare soil from those vegetated areas.

4.5. Comparison with Other Machine Learning Classifiers

In order to further justify the performance of random forest classifier, it should be compared to other machine learning methods such as support vector machine (SVM) and decision tree (DT) algorithm. Additional comparison experiments were done based on both SVM and DT to classify the multi-temporal remote sensing data of the years 2015, 2005, 1995 and 1986. The same training and validating data were used during all the classifications of RF, SVM and DT to ensure comparability. Specifically, radius basis function (RBF) was utilized as the kernel function of SVM and the DT involved was the classification and regression tree (CART) algorithm. Both the SVM and DT classifications were carried out using the Supervised Classification tool of ENVI 5.1. The accuracy comparison results are listed in Table 10.

Table 10. Accuracy comparison of RF to SVM and DT.

	RF		SVM		DT
	OA	Kappa	OA	Kappa	OA	Kappa
2015	91.26%	0.9000	85.68%	0.8363	82.90%	0.8045
2005	90.34%	0.8895	88.87%	0.8727	85.67%	0.8362
1995	90.80%	0.8948	86.15%	0.8417	84.21%	0.8195
1986	85.37%	0.8328	82.93%	0.8049	82.65%	0.8017

Notes: RF, random forest; SVM, support vector machine; DT, decision tree.

Table 10. Accuracy comparison of RF to SVM and DT.

	RF		SVM		DT
	OA	Kappa	OA	Kappa	OA	Kappa
2015	91.26%	0.9000	85.68%	0.8363	82.90%	0.8045
2005	90.34%	0.8895	88.87%	0.8727	85.67%	0.8362
1995	90.80%	0.8948	86.15%	0.8417	84.21%	0.8195
1986	85.37%	0.8328	82.93%	0.8049	82.65%	0.8017

Notes: RF, random forest; SVM, support vector machine; DT, decision tree.

Table 10 indicates that the random forest classifier outperforms both SVM and DT. Specifically, the OA of RF is slightly higher than that of SVM with an increase of 1.47% to 5.58%. One possible reason is that RF is more robust than SVM in handling high dimensional datasets with some redundancies. The datasets used in this study have a large number of input features due to the layer stacking of multi-temporal images (i.e., 57 input bands for 2015, 54 for 2005 and 1995, and 36 for 1986), and undeniably, some of the input features have some redundancies. However, RF is based on the boosting strategy, which is more robust than SVM in dealing with the collinearity of the training data.

Meanwhile, the OA of RF surpasses that of DT by an increase of 2.72% to 8.36%. This is predictable mainly because RF is the ensemble of CART classifiers and the boosting strategy can improve the accuracy of RF when compared to a single decision tree. In addition, over-fitting of the training data is less an issue for RF than that of DT, which also makes the former outperform the latter.

4.6. Driving Forces and Implications

The research of the driving forces behind the LUCC [21,22] has been recognized as very important for interpreting the reasons and rationales and providing policy-relevant information. In general, the driving forces mainly consist of two interrelated parts, i.e., natural and anthropogenic factors [21,22]. At the regional scale, anthropogenic factors are considered to be more important in characterizing the reasons behind the LULC changes [22]. In this section, we analyzed the driving forces of the cropland changes to provide decision support for the sustainability of the Yellow River Delta.

The development of oil related industries has resulted in the extension of built up areas [5,21,22,34]. It is observed that built up areas increased from 52.1 km² in 1986 to 221.3 km² in 2015 with a growth rate of 5.83 km²/year. Meanwhile, the area of cropland converted into built up areas was 17.4 km², 19.2 km² and 29.9 km² for the three study periods, which verifies that the urban sprawl is one of the driving forces behind the losses of croplands. Moreover, the fact that large croplands being occupied by the oil fields, chemistry factories and residential areas would induce the reclamation of grassland and bare soil to maintain the stable food production and food security. However, the oil related industries have caused several environmental problems, including the oil leakage, the point pollution due to poor or non-treated sewage, to name a few, which could impair the soil fertility of the croplands. In order to protect the Yellow River Delta’s fragile environment from further destruction, strict measures should be taken by the local government. For instance, relative laws and regulations of wetland management should be improved and perfected [22]. Regular monitoring of oil leakage and industrial sewage should be done to protect the aquatic ecosystem. In addition, the irrational economic activities (e.g., establishing oil well in the nature reserves) should be stopped.

The development of both agriculture technology and food marketing in China [34] has spurred the local people and several companies to reclaim large amounts of uncultivated land into farmland to increase their profit. It is observed that grassland and bare soil are the two major contributors to the increase of cropland in the Yellow River Delta during the past three decades. However, the cropland mutual conversion analysis also revealed that there existed the phenomenon of abandoned cultivation due to the decline of soil fertility and soil salinization. The frequent transformation between cropland, grassland and bare soil is the major pattern of LULC changes in the Yellow River Delta. To maintain the sustainable development of croplands, future concerns should be laid on the improvement of cropland’s soil quality instead of reclaiming the uncultivated lands to help protect the fragile local ecosystem.

In addition, the fact that the Yellow River is the main contributor to the fresh water of the delta [5,34] means that the reclamation of cultivated lands will decrease the flux of the Yellow River, which in turn will affect the natural evolution of the entire delta’s ecosystem. The usage of chemical fertilizer and pesticides can also lead to the degradation of soil fertility and non-point pollution to the water resource.

Due to the fact that the value of wetlands has been recognized recently in China, two nature reserves have been established according to the Ramsar Convention on Wetlands [26,27]. Farming is prohibited in those reserves for ecological restoration. Meanwhile, wetland tourism has been more and more popular in recent years, which can compensate the economic loss due to the prohibition of agricultural activities and can improve the sustainability of the entire Yellow River Delta.

Above all, a sustainable perspective should be taken into account seriously for the coordinated development of agriculture, industry and tourism. The improvement of soil quality is important to avoid the phenomenon of abandoned cultivation and the reclamation of grassland and bare soil. The usage of agriculture chemicals should also be controlled to mitigate the degeneration of croplands. In the meantime, pollution from the oil related industry should be eradicated.

5. Conclusions

This paper is aimed at monitoring the cropland dynamics of the Yellow River Delta using multi-temporal Landsat data for the last three decades. Results indicated that the proposed methodology showed good performance in producing the LULC maps of the Yellow River Delta. The inclusion of multi-temporal data has improved the classification accuracy significantly with an overall accuracy increase of 18.54% and a Kappa increase of 0.1572. The adoption of the random forest classifier has yielded good performance with an average accuracy of 89.44% and an average Kappa of 0.8793.

Spatial-temporal analysis demonstrated that the Yellow River Delta experienced severe cropland changes during the past three decades. The cropland area increased from 467.6 km² in 1986 to 718.5 km² in 2015 with an average growth rate of 8.65 km²/year. In addition, the period of 1995 to 2005 witnessed the highest rate of cropland increase (22.38 km²/year) while the period of 1986 to 1995 and 2005 to 2015 experienced rather little change. Due to technical progress and business stimulus, large amount of grassland and bare soil have been reclaimed as cultivated land, resulting in the steady increase of cropland area. Meanwhile, the losses of cropland were mainly from abandoned cultivation due to the decline of soil fertility. As well, the urban sprawl caused by the fast development of oil related industry has also occupied certain amount of cropland. Based on these facts, we suggest that decision makers should adopt a sustainable perspective in the future to keep the balance between food security, industrial development and ecosystem protection. Relative laws and regulations for wetland management should be improved and the monitoring of oil leakage and industrial sewage should be done regularly. All the irrational economic activities should be ceased at the same time.

Above all, this paper discussed the cropland dynamics of the Yellow River Delta based on Earth Observation data and provided some policy-relevant information for decision makers. Future studies should focus on the dynamics of the overall LUCC of the Yellow River Delta and research into the quantitative relationship between the landscape changing patterns and the socio-economic indicators to achieve a further understanding for the sustainability of the Yellow River Delta.