Multi-Temporal Arable Land Monitoring in Arid Region of Northwest China Using a New Extraction Index

Abstract

Development of a high-accuracy method to extract arable land using effective data sources is crucial to detect and monitor arable land dynamics, servicing land protection and sustainable development. In this study, a new arable land extraction index (ALEI) based on spectral analysis was proposed, examined by ground truth data, and then applied to the Hexi Corridor in northwest China. The arable land and its change patterns during 1990–2020 were extracted and identified using 40 Landsat TM/OLI images acquired in 1990, 2000, 2010, and 2020. The results demonstrated that the proposed method can distinguish arable land areas accurately, with the User’s (Producer’s) accuracy and overall accuracy (kappa coefficient) exceeding 0.90 (0.88) and 0.89 (0.87), respectively. The mean relative error calculated using field survey data obtained in 2012 and 2020 was 0.169 and 0.191, respectively, indicating the feasibility of the ALEI method in arable land extracting. The study found that arable land area in the Hexi Corridor was 13217.58 km² in 2020, significantly increased by 25.33% compared to that in 1990. At 10-year intervals, the arable land experienced different change patterns. The study results indicate that ALEI index is a promising tool used to effectively extract arable land in the arid area.

1. Introduction

Human-induced land use/cover change has caused land development (e.g., arable land, urban land), expanding and encroaching into natural areas [1,2]. Arable land in arid regions is mainly distributed in artificial oases. Compared to other land use/cover types in an arid region, e.g., woodland, arable land is more vulnerable to extreme weather and human disturbance [3,4]. Accurately detecting arable land areas and monitoring its dynamic change is indispensible for the development of land management policies to ensure sustainable development.

Early land use/cover change studies have utilized historical records and hand-drawn maps to depict arable land in arid regions [5,6]. Since the 1990s, land use/cover maps from remotely sensed images have been employed as the main data source to monitor the arable land’s status and its response to external disturbance [7,8,9], because remotely sensed images can help retrieve the context, direction, and rate of arable land change in a time- and labor-saving way [10,11]. According to the spatial resolution, remotely sensed images can be categorized to coarse-, middle- and high-resolution classes. High-resolution (<10 m) images are expensive and the imaging scan width is relatively small, and thus are generally not feasible for large-scope land use/cover monitoring (e.g., >10 km²). Coarse-resolution images, such as MODIS (250 m/500 m/1 km) and AVHRR (1 km) series data set, have been widely used to study the land use/cover change at global and regional levels [12,13,14]. However, due to their coarse resolution, pixels are more likely to contain different land use/cover types, which are referred to as mixed pixels. Compared to pure pixel (a pixel that contains only one land use/cover type), a mixed pixel in the coarse-resolution images can lead to low accuracy of arable land detection. Given the high demand of precise land use/cover monitoring, the development of cloud computing platforms such as Google Earth Engine, and efficient methods to process images [15,16], middle-resolution Landsat series images, as first launched in 1972, have gradually became one of the most practical data sources to delineate and monitor the dynamics of land surface change at a regional scale (e.g., [17,18,19]). Compared to the coarse-resolution images and other middle-resolution images, Landsat images can identify land use/cover at low cost and acceptable accuracy [20,21,22] and help monitor the dynamic change of land uses/cover in longer time intervals (1972 to date).

There are volumes and varieties of Landsat data available, however, methods used to delineate arable land at a regional scale suffered development lag. For a long period, a high proportion of information from remotely sensed images was distinguished by visual interpretation, which is labor intensive and time consuming [23]. To improve the efficiency of arable land identification and minimize the artificial errors, vegetation indices, such as Normalized Difference Vegetation Index (NDVI), were employed by studies to differentiate arable land areas [24]. For instance, Xie et al. (2018) [25] used the NDVI threshold method for vegetation extraction and monitored the spatio-temporal changes of oases at a regional scale and found that different NDVI values can help identify the boundary of oases. Venkatappa et al. (2019) [26] developed a robust phenology-based classification method to accurately determine the threshold values of vegetation types.

Vegetation Index-based method is easy to delineate the boundary of arable land and desert, for they show apparently different patterns in the Infrared and Red spectral bands. Recently, the open, shared data sets of Globeland30 that fulfilled the 10-year-interval land use/cover classification were released for free download [27]. However, their classification accuracy of arable land in northwest China is 0.8315 [28], which can be further improved. Besides, in the agro-pastoral zones of arid region, vegetation index values of arable land and other vegetation types (woodland, shaded vegetation, etc.) share range overlap due to their similar spectral change patterns, and thus the vegetation index threshold method may not be sufficient to distinguish arable land in the complicated agro-pastoral transition zones. Currently used classification methods may misclassify the woodland and shaded vegetation to arable land, resulting in more bias and uncertainties of arable land extraction in arid region. Therefore, they need a feasible algorithm to extract arable land to ensure accurate detection of arable land areas in an arid region.

To this end, this study proposed a framework to efficiently distinguish arable land area and its dynamic change patterns. Multi-temporal Landsat images and different sources of auxiliary data were collected and processed to examine the performance of the proposed framework. The method can be used for arable land monitoring in an arid region, and the results can be used as a reference for the agricultural development in the Hexi Corridor region.

2. Materials and Methods

2.1. Study Area

The Hexi Corridor (90°11′–101°29′ E, 37°14′–42°15′ N) is a belt region located in northwest China. The altitude of the study area ranges from 790 to 5780 m and the climate is arid, temperate monsoon, characterized by a hot summer and cold winter. The annual average temperature during the study period ranged from 4 to 10 °C, and the annual average precipitation was from 100 to 400 mm in the region, showing a decreasing trend from southeast to northwest. The region is comprised of three inland rivers, all originated from the Qinghai-Tibet Plateau, including the Shiyanghe river basin, the Heihe river basin, and the Shulehe river basin, distributing from southeast to northwest (Figure 1). Arable land in the study area is mainly distributed in the oases of the Jinchang, Minqin, Sunan, Wuwei, and Yongchang administrative counties in the Shiyanghe river basin; in Gaotai, Jiayuguan, Jinta, Jiuquan, Linze, Minle, Shandan, and Zhangye administrative counties in the Heihe river basin; and in Akesai, Dunhuang, Guazhou, Subei, and Yumen administrative counties in the Shulehe river basin. By thousands of years of cultivation, the arable land has differed clearly from other land use types, so we selected this region to examine the method for arable land extracting.

2.2. General Prodcedure of Data Processing

The arable land was extracted in four successive steps (Figure 2): (1) data selection and processing, (2) pure pixels’ selection, (3) index framework setup, and (4) arable land extracting and accuracy estimation.

2.2.1. Data Sources and Preprocessing

A total of 40 Landsat images (10 images for each year, see Appendix A) were collected for the years 1990, 2000, 2010, and 2020. The Landsat-series images covered seven paths from east to west (131–137) and three rows from north to south (032–034) and contained 30 Landsat-5 TM and 10 Landsat-8 OLI images acquired in the summer season. The Google Earth Engine [21,29] was used to identify the images that were cloud-free and possessed the highest crop coverage rates through the years. Thus, 40 raw format images acquired from 16 June to 6 September during the 30 years were selected and downloaded from the public domain of the Earth Resources Observatory and Science Center (EROS, http://eros.usgs.gov, accessed on 18 October 2020) using a Python bulk download application. The images were then rectified to the Universal Transverse Mercator (UTM) projection using at least 45 well-distributed ground control points (GCPs) and nearest neighborhood resampling method for each image, and the root mean square error (RMSE) was set to be less than 1 pixel. Radiometric calibration, fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) atmospheric correction, and geometry correction were conducted to obtain surface reflectance images using a batch-processing IDL program.

2.2.2. Selection of Pure Pixels

Owing to the high separability characteristics between vegetation (arable land, woodland, etc.) and other land use/cover types (desert, urban and rural areas, glacier, waterbody, etc.), the NDVI threshold method can be used to distinguish and locate vegetation areas in the remotely sensed images. However, due to apparent overlaps in the reflectance and impact of terrain fluctuation, arable land is hard to be accurately separated from woodland and shaded vegetation based on NDVI differences (Figure 3). Here the three land use/cover types (arable land, woodland, and shaded vegetation) were visually interpreted from five randomly selected Landsat images for each year with the aid of Google Earth Pro images, and their NDVI values were calculated using the equation below.

$$NDVI=Float\left(NIR-R\right)/ Float \left(NIR+R\right)$$

(1)

where NIR represents the reflectance value of Near-infrared band, and R is the reflectance value of Red band.

The Minimum Noise Fraction Transform (MNFT) and Pixel Purity Index (PPI) were conducted on the Landsat images to find pure pixels, which were used to detect cues of spectral difference between arable land and other vegetation types [30,31]. At first, the MNFT was used to eliminate the noise in the images by principal component transformations in two steps. The first transformation (based on the estimated noise covariance matrix) was to separate and re-adjust the noise in the data, by minimizing the variance of the transformed noise data without correlation between bands. The second step was to transform the standard principal component of the noise-whitened data [32].

The pixel purity index (PPI) has been widely used in multi-spectral image analysis for pure pixel detection. After the MNFT processing, a set of k unit vectors ${\left\{skewe{r}_j\right\}}_{j=1}^k$ was generated, where k is a sufficiently large, positive integer. For each $skewe{r}_j$, all the data sample vectors were projected onto $skewe{r}_j$ to find sample vectors at its extreme positions to form an extrema set for this particular $skewe{r}_j$, denoted by S_extrema($skewe{r}_j$). Despite the fact that a different $skewe{r}_j$ generates a different extrema set S_extrema($skewe{r}_j$), it is likely that some sample vectors may appear in more than one extrema set. Define an indicator function of a set S, I_S(r) by

$$\begin{gathered} I_S\left(r\right)=\{\begin{array}{c}1, if r\in S\\ 0, if r\notin S\end{array}\mathrm{and} \\ {\mathrm{N}}_{PPI}\left(r\right)={\displaystyle \sum }_j{I}_{Sextrema\left(skewe{r}_j^k\right)}\left(r\right) \end{gathered}$$

(2)

where ${\mathrm{N}}_{PPI}\left(r\right)$ is defined to be the PPI score of sample vector $r$ [33]. A large PPI value (e.g., 1) indicates a high possibility of a pure pixel [34].

The pixels with high PPI values were then exported to establish a new region of interest (ROI) and then linked to the Google Earth Pro images (0.5-m spatial resolution) to further examine the purity of these pixels. Finally, about 400 pure pixel samples were selected in each image for further analysis.

Figure 4 shows the reflectance of the six bands (Blue, Green, Red, Infrared, shortwave 1, and shortwave 2) of the p132r34 image as an example. Comparison analysis found that the spectral change patterns of arable land, woodland, and shaded vegetation were similar to each other. However, the value of the Infrared band was apparently different from other bands and, thus, this band was selected as the basis to design the framework to distinguish arable land areas.

2.2.3. Formulation of Arable Land Extraction Index (ALEI)

Spectral difference is an important foundation to discern arable land. As the wavelength increasing from visible to the shortwave bands, the reflectance of arable land and natural vegetation showed the same change pattern. They increased in the Green band, slightly decreased in the Red band, significantly increased in the Infrared band, and then decreased in the Shortwave bands (Figure 4). The shaded vegetation shared the same change pattern with that of woodland but the reflectance in the Infrared was much less than those of arable land and woodland (Figure 4) and the other images shared the same spectral distribution characteristics. This can help separate shaded vegetation and eliminate its impact on arable land extract. Therefore, the main issue was to separate arable land from woodland in the agro-pastoral zones. After analyzing the spectral change patterns of arable land and woodland, the Infrared and Shortwave 1 were selected to form an arable land extraction index to maximize the contrast between arable land and woodland.

$$ALEI={\rho }_{Infrared}-\omega \ast {\rho }_{Shortwave 1}$$

(3)

where ${\rho }_{Infrared}$ is the reflectance value of the Infrared band and ${\rho }_{Shortwave 1}$ is the reflectance value of the first shortwave band.

For successive identification of arable land by excluding other easily confused land use/cover types, the key is to determine the coefficient ω to maximize the difference between arable land and natural vegetation. So, iteration computations with different ω values were conducted, and the results are presented in Figure 5b. Considering the image acquisition date and land use/cover distribution in the study area, the value cannot be too small or too large, as such may lower the desired difference to separate arable land and woodland. For instance, if ω were 1.4, the index values for the two land use/cover types are positive, and when the ω value increases, e.g., to larger than 2.2, their values become negative. By iterative examinations, the suitable ω value can be set to 1.6 or 2.0 for the formula, but they both fail to fulfill the purpose of maximizing the difference between arable land and woodland. By assigning the coefficient ω to be 1.8, the pixel with positive index value represents arable land and the pixel-processing negative one is woodland. Here, we consented that 0 was a default threshold for separating arable land from woodland, because pixels with positive and negative ALEI values are easier to be classified into different land use/cover categories.

2.2.4. Verification

After the determination of ω parameter, the ALEI was applied to extract the arable land. The arable land extraction framework was built using a two-branch decision tree method. At first, “B1 LT 2500” was used to eliminate the impact of shaded vegetation to arable land extraction, where “B1” is the abbreviation of Band 1, i.e., the Infrared bands of Landsat image, and “LT” represents “Less Than”, i.e., the pixels in the Landsat images with Infrared value of less than 0.25 that were classified to shaded vegetation. Secondly, the remaining pixels were further processed using the “B2 GT 1.8” to distinguish arable land, in which the “B2” (Band 2) is the ALEI-based results calculated using Equation (3) and “GT” is the abbreviation of “Greater Than”. The pixels with ALEI values greater than 1.8 were grouped into the arable land class. In this study, 40 decision trees were established separately and 40 different B1 and B2 sets were used to obtain the computation results.

A 160,000-pixel ground truth data set was randomly selected, and the land use types for each image were visually interpreted based on Google Earth Pro online images (0.5-m spatial resolution). With this ground truth data set, the currency of the ALEI-based results was examined using four coefficients, i.e., user’s accuracy, producer’s accuracy, overall accuracy, and Kappa coefficient [35]. The user’s accuracy is the ratio of the total number of pixels correctly classified into class i (diagonal value in the confusion matrix) and the total number of pixels classified into class i by the classifier (the sum of class i rows in the confusion matrix). The producer’s accuracy is the ratio of the number of pixels (diagonal value in the confusion matrix) that the classifier correctly classifies for the entire image into class i and the true total number of class i (the sum of class i columns in the confusion matrix).

Furthermore, intensive field surveys were conducted for 15 arable land areas in July 2012 and July 2020, for each of which the arable land area boundary was recorded by walking around it using a high-precision tracking GPS (position dilution of precision ≤4, error < 10 m). The recorded boundary for each of the 15 arable land areas was input into ArcGIS 10.6, and then the area was computed. These data were then used for in situ verification of the arable land area obtained from the ALEI method in 2010 and 2020, with the mean relative error (MRE) as below.

$$MRE=\frac{{{\displaystyle \sum }}_{i=1}^n\left(abs(x_i-y_i)/y_i\right)}{n}$$

(4)

where $x_i$ is the ith area obtained from field survey, $y_i$ is the ith area from ALEI method, and $n$ is 15, the total number of the field survey locations.

2.3. Dynamic Change Monitoring of Arable Land

Arable land dynamic change during 1990–2000, 2000–2010, and 2010–2020 in the three river basins was computed by the equations below:

$$\overrightarrow{R}=\left(\sqrt[{\displaystyle t}]{\frac{A_b}{A_a}}-1\right)\ast 100\%$$

(5)

$$T_s=\frac{\mathsf{\Delta }{A}_{gain}-\mathsf{\Delta }{A}_{loss}}{\mathsf{\Delta }{A}_{gain}+\mathsf{\Delta }{A}_{loss}},\mathsf{\Delta }{A}_{gain}+\mathsf{\Delta }{A}_{loss}\ne 0, T_s\in \left[-1,1\right]$$

(6)

where $\overrightarrow{R}$ is the net change rate index and $T_s$ is the status and change trend index. $A_a$ and $A_b$ are the total area of arable land in the initial and end year of the research period, $\mathsf{\Delta }{A}_{gain}$ and $\mathsf{\Delta }{A}_{loss}$ represent the area gained and lost during 1990–2020, respectively, and $t$ is the time duration of the time interval of 10 years. From the Equations (5) and (6), the higher the $\overrightarrow{R}$ is, the more the overall growth is. When $T_s$ is less than 0, the total arable land is reduced; otherwise, it is increased.

3. Results and Discussion

3.1. Verification of Arable Land Extraction

The overall accuracy (kappa coefficient) values of arable land extraction using the ALEI method all exceeded 0.89 (0.87), and both user’s and producer’s accuracies were higher than 0.84, indicating that the classification results were feasible for arable land change analysis (

Table 1. Accuracy assessment of arable land and other vegetation types.

Land Cover	Accuracy Type	Year
		1990	2000	2010	2020
Arable land	User’s Accuracy	0.90–0.92	0.92–0.94	0.92–0.95	0.93–0.94
	Producer’s Accuracy	0.88–0.91	0.91–0.94	0.90–0.92	0.92–0.93
Woodland	User’s Accuracy	0.87–0.90	0.89–0.93	0.89–0.91	0.89–0.93
	Producer’s Accuracy	0.85–0.89	0.86–0.91	0.87–0.90	0.88–0.91
Shaded vegetation	User’s Accuracy	0.90–0.93	0.92–0.95	0.90–0.92	0.92–0.94
	Producer’s Accuracy	0.90–0.93	0.90–0.93	0.89–0.92	0.89–0.91
All	Overall Accuracy	0.89–0.93	0.91–0.94	0.91–0.94	0.90–0.93
All	Kappa Coefficient	0.87–0.90	0.89–0.91	0.90–0.92	0.89–0.91

). Four 160,000-pixel zones were randomly selected in the images of p131r34, p132r34, p133r33, and p135r32 from southeast to northwest to visually show the extraction results. The original Landsat images, the arable land results using the ALEI and the NDVI threshold method, are shown in Figure 6. Both the ALEI and NDVI methods can distinguish vegetation from non-vegetation land use/cover types such as cloud, waterbody, riverbed, urban, and rural areas. However, NDVI-based method confused arable land areas with woodland and shaded vegetation (Figure 6f,l), which resulted in overestimation of arable land areas.

The accuracy of arable land extraction was further verified using areas of arable land obtained in the field survey in the summer season of 2012 and 2020. The correlation between field survey data and ALEI-based data are shown in Figure 7. Both the areas of arable land in 2010 and 2020 showed a high correlation, the coefficients of the years 2010 and 2020 were 0.969 and 0.963, and the MRE values were 0.169 and 0.191, respectively. However, the ALEI results were slightly overestimated compared to the field survey data. This was mainly attributed to the mixing pixels existing in these Landsat images. For instance, a pixel that contained 80% arable land and 20% other land use/cover types was classified as arable land using the ALEI method, resulting in some error. Overall, the classification accuracy using the ALEI method was acceptable and the results can be used to monitor the dynamic change of arable land in the study area.

3.2. Change Pattern of Arable Land in the Study Area

The extracted arable land in the study area during 1990–2020 are presented in Figure 8. In 1990, the total area of arable land was 10,546.33 km², of which 34.45%, 46.42%, and 19.13% was distributed in the Shiyanghe, the Heihe, and the Shulehe river basins, respectively. In 2020, the total arable land area was 13,217.58 km², increased by 25.33% compared to that in 1990, which reached its extreme value. During the three 10-year study intervals, arable land in the Hexi Corridor showed a continually increasing trend with the increase rates of 4.99%, 10.86%, and 7.43%, respectively. The most significant increase occurred during 2000–2010.

In the three inland river basins, arable land showed a significantly expanding trend during the whole study period, while they showed different change patterns in the three 10-year intervals:

1.

In the Shiyanghe river basin, arable land area increased rapidly with an annual rate of 2.36 % from 1990 to 2020, which was the lowest among the three inland river basins (Heihe river basin 2.60%/year, Shulehe river basin 2.68%/year). In this river basin, the most rapid arable land expansion occurred during 1990–2000, with a net change rate $\overrightarrow{R}$ of 0.99%. The trend index $T_s$ of this 10-year interval was 0.51. It had not much difference with the other two 10-year periods (0.53 during 2000–2010 and 0.50 in 2010–2020 interval,

Table 2. Change of arable land in the three 10-year intervals.

River Basin	Net Change Rate ($\overrightarrow{R}$, %)/Status and Trend Index ($T_s$)
	1st Interval	2nd Interval	3rd Interval
Shiyanghe	0.99/0.51	0.51/0.53	0.63/0.50
Heihe	0.17/0.18	1.33/0.82	0.83/0.78
Shulehe	0.33/0.23	1.32/0.69	0.74/0.52

). This indicated that, though arable land in the Shiyanghe basin expanded prominently, the area lost was also significant (Figure 9).

2.

Arable land in the Heihe river basin during 1990–2000 extended slightly in area (0.17 %/year), and the peak cumulative period of arable land occurred in the second 10-year interval (2000–2010, Figure 9). The area of arable land went up to 5681.04 km² with the net change rate $\overrightarrow{R}$ of 1.33 % and trend index $T_s$ of 0.82, both of which were the highest among the three river basins. During 2010–2020, arable land in the Heihe river basin accumulated another 488.15 km², and the net change rate $\overrightarrow{R}$ and trend index $T_s$ showed a decrease trend compared to those in the second 10-year interval (Table 2).

3.

Arable land area of the Shulehe river basin in 1990 (2017.87 km2) only accounted for 55.55% and 41.22 % of that in the Shiyanghe and Heihe basins, respectively. However, arable land in the Shulehe river basin showed a similar change trend with that in the Heihe river basin (Figure 8 and Table 2). During 2000–2010, the area of arable land in the Shulehe basin increased by 292.12 km2 (14.09 %) and during 2010–2020, it further expanded by 7.61 % (Figure 9).

3.3. Discussion

This study established a framework to distinguish arable land and its change patterns in the Hexi Corridor, northwest China. The ALEI method was proposed to distinguish arable land based on 40 Landsat series images. The results demonstrated that the ALEI method is a feasible tool for extracting arable land areas accurately. The arable land in the study area showed an ongoing, expanding trend during 1990–2020, while the three inland river basins experienced different change patterns among the three 10-year intervals.

Compared to the NDVI threshold method, ALEI can help distinguish arable land in the study area more accurately, especially in the agro-pastoral zones. The verification results demonstrated that the producer’s accuracy of arable land extraction using the ALEI method ranged from 0.88 to 0.94, which increased by 6.02% to 13.25% compared to that of the Globeland30 (0.8315) [28]. The study results revealed the status and change patterns of arable land in the Hexi Corridor and can be used as references for decision makers to implement more targeted land protection policies and regulations, to balance the socioeconomic and ecological conflicts.

The ALEI method uses an efficient process framework in which the works can be all or semi-automatically conducted, i.e., ALEI can be computed with a batch-processing module automatically, after the Landsat images are downloaded using a bulk download application. Overall, the proposed methodology is understandable and easy to use. As artificial oases in arid regions share similar vegetation types and Landsat image is one of the most popular remotely sensed data sets covering most of the arid regions worldwide [36], ALEI is a promising tool to accurately extract arable land for land use planning in other arid regions. However, there are some limitations that were noted. Firstly, the field surveys for the verification was limited, covering an area of less than 4 km², because large arable land areas are hard to be measured using the “walking-GPS” way. Recently, an unmanned aerial system (UAS) was applied in the remote sensing studies [37,38] to measure large, arable land areas as the ground truth data. Considering that the field data in 2012 were collected by the “walking-GPS”, so the data in 2020 were also collected by the same way to maintain the consistency and minimize the uncertainties. Secondly, the ALEI method was used and tested in the Hexi Corridor, and more applications in other arid regions could be done to further examine its usability and accuracy.

4. Conclusions

This study developed, examined, and applied a framework and an index of ALEI to identify arable land and its change patterns during past decades in a representative arid region. The results demonstrated that the proposed framework can help identify arable land areas accurately and found that the arable land in the study area extended significantly during 1990–2020. The findings of this study offer essential data to support the formulation and implementation of arable land protection and development policies.