Assessing spatio-temporal patterns and driving force of ecosystem service value in the main urban area of Guangzhou

This paper uses Landsat TM/OLI images from 1987, 1993, 1999, 2005, 2011 and 2017 by the weight vector AdaBoost (WV AdaBoost) multi-classification algorithm extracting LULC data sets, and the spatiotemporal patterns of LULC over these periods were analysed. The spatiotemporal change patterns and driving force of ESV was estimated. The effect of LULC dynamics on the ESV was evaluated. The flow chart is shown in Fig. 1.

Study area

Guangzhou is located at 112° 57′ ~ 114° 3′ E, 22° 26′ ~ 23° 56′ N in the southeast part of Guangdong Province in the northern margin of the Pearl River Delta. It has a total land area of 7434.40 km². The topography is high in the northeast and low in the southwest, with low mountains and hills in the north and plains in the south and a rich geomorphology. The level of urbanization is high, the land use structure is complex and the land use pattern is changing rapidly. Regional public infrastructure, commercial service land and external transportation land is increasing. Export-oriented industrial agglomeration areas are developing continuously, and the characteristics of export-oriented land use are obvious. Guangzhou has a subtropical oceanic monsoon climate, with an annual average temperature of 20–22 °C and an annual rainfall of about 1720 mm. Guangzhou includes eleven districts (Fig. 2a). This paper focuses on the main urban area of Guangzhou as the research object, a decision based on the city government’s overall urban development strategy for Guangzhou (2010–2020), the main urban area of Guangzhou includes five districts—Liwan, Yuexiu, Haizhu, Baiyun and Tianhe (Fig. 2a). In 2017, the per capita GDP will exceed 136,100 yuan for the first time (Fig. 2b). The population density is large, and the intensity of development is extraordinary. The population has increased from 2.73 million in 1987 to 8.05 million in 2017 (Fig. 2c). The industrial structure has been constantly optimized and improved. The proportion of secondary industries decreased from 31.46% in 1987 to 11.28% in 2017, while the proportion of tertiary industries increased from 61.9% in 1987 to 89.61% in 2017.

Data and data processing

Data

The data used include Landsat series images, an administrative zoning map of Guangzhou, a local historical land use map and social and economic statistics from 1987 to 2017. To explore the socio-economic and natural factors driving the ESV change for different land use types, we examined ten factors—the normalized vegetation index (NDVI), year-end population, elevation, slope, distance from roads, distance from railways, land cover types, GDP, secondary industry GDP and investment in fixed assets.

The Landsat images were taken from the United States Geological Survey website (http://glovis.usgs.gov/), and we downloaded cloud-free Landsat 5 TM and Landsat 8 OLI images (path 122, row 44). Images taken in the dry season (October, November and December) were selected because there was less cloud cover, the change in surface reflectance was much smaller than in other seasons and the image quality was higher. The spatial resolution of imaging is 30 m, and the data of Landsat 5 TM (8 December 1987, 24 December 1993, 25 December 1999, 23 November 2005, 19 November 2011) and Landsat 8 OLI (23 October 2017) were collected.

The elevation data are provided by the Japan Aerospace Exploration Agency (http://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm); the horizontal resolution is 30 m, and the elevation accuracy is 5 m. The road and railway data are provided by Openstreetmap (download.geofabrik.de/), and the distance between roads and railways were calculated. The land use types are provided by Tsinghua University (http://data.ess.tsinghua.edu.cn/) with a resolution of 10 m. GDP, secondary industry GDP and investment in fixed assets were provided by Guangzhou Statistics Bureau²⁸.

Considering the size of the study area and the resolution of the data, we used Fishnet in ArcGIS10.5 to establish a 500 × 500 m grid covering the study area. The values of ten factors were extracted to the grid centre points, and 3915 points were obtained.

Remote sensing extraction of LULC types

To obtain high-precision LULC data, the surface features based on the original remote sensing image were enhanced using the index model (the Index-based Built-up Index (IBI)²⁹, the modified normalized difference water index (MNDWI)³⁰ and the soil-adjusted vegetation index (SAVI)). The humidity, brightness and green chroma indices were transformed using the tasselled-cap method^23. Using the six indices, six images were calculated and superposed into a new multi-band image in the order of IBI, SAVI, MNDWI, brightness, green chroma and humidity (referred to as the six-index image). This was superposed with the original image (band 1, 2, 3, 4, 5, 7) to create an enhanced 12-band image as the data source for LULC classification. The image stretching function was then used to stretch the composite image to better distinguish the characteristics of land use types. After feature enhancement, the LULC types in the main urban area of Guangzhou could be divided into seven types—forest, water body, grassland, cultivated land, high reflectivity building, low reflectivity building and bare land. Of these, high reflectance buildings and low reflectance buildings are sub-categories of the built-up area category. The spectrum of high reflectivity buildings is similar to that of bare land; therefore, they are classified separately to avoid erroneous classification³¹. Based on Google Earth and TM/OLI images, the interpretation marks of LULC types were established (Table 1). To obtain pure and representative samples, the K-means method of unsupervised classification was utilized to cluster the pixels of the ISM image into seven unlabelled classes and to eliminate abnormal clustering. We then selected uniform points and drew polygons on them. Each polygon contains more than 100 pixels, and each class contains between five and seven polygons. All polygons were saved in a shape file as a mask file, which is utilized to extract pure pixels as training samples.

Dou and Chen³¹ suggest using the weight vector AdaBoost (WV AdaBoost) multi-classification algorithm for LULC classification. Compared with AdaBoost, it provides higher classification accuracy and stability. WV AdaBoost includes a C4.5 decision tree, a Naïve Bayes neural network and an artificial neural network. In this study, we use Naïve Bayes-based WV AdaBoost to classify LULC based on Landsat-enhanced images. To obtain a highly accurate LULC classification, it is necessary to post-process the image after classification, overlay the early land use map during the research period, check the incorrectly classified area, and filter and correct some segments. After image classification, the random selection of 200 pixels in each class was checked, and the correctness of the classification and the evaluated accuracy were confirmed.

The average classification accuracy of WV AdaBoost based on Naïve Bayes on the original image is 85.02%, and the kappa coefficient is 0.839. The WV AdaBoost algorithm is based on Naïve Bayes processes for the 12-band images of the new combination (enhanced 12-band image). The average classification accuracy is 88.86%, and the kappa coefficient is 0.871. The classification results still need to be strengthened by post-processing, which achieves good classification accuracy. The final average classification accuracy is 91.97%. The kappa coefficient is 0.907. These classification results agree with the ensuing analysis of LULC.

Methods

Based on Landsat TM/OLI data from 1987 to 2017, a transfer matrix, a land use change index, an ESV evaluation index, a sensitivity model (CS) and a geo-detector (P) were used to analyze the response of ESVs to LULC evolution.

Transfer matrix

The transfer matrix reflects the dynamic process information about mutual transformation LULC types at the beginning (T) and the end (T + 1) of a specified period of time in a certain region (Fig. 3). The general form is:

where S stands for area, and i,j (i,j = 1,2,…, n) represents LULC types before and after the transfer.

The LULC analyzing indices

1. (1)

   Land use dynamics (RS)

   Land use dynamics describe the rate and magnitude of LUCC. The general form is:

   $$RS=frac{{U}_{b}-{U}_{a}}{{U}_{a}}times frac{1}{T}times 100mathrm{%},$$

   (2)

   where ({U}_{a}) is the area of a certain land class at the beginning (km²), ({U}_{b}) is the area of a certain land class at the end (km²), and T is the study period.

2. (2)

   Spatial dynamics of land use (RSS)

   Spatial dynamics of land use describe the degree of spatial change in a certain land use type. The general form is:

   $$RSS=frac{{U}_{in}-{U}_{out}}{{U}_{a}}times frac{1}{T}times 100mathrm{%},$$

   (3)

   where ({U}_{in}) is the sum area of other types converted to this type in study period T, and ({U}_{out}) is the sum area of a certain type converted to other types in study period T. ({U}_{a}) is the area of a certain type at the beginning.
3. (3)

   Land use change state index (PS)

The land use change state index represents the trend and state of LUCC. The general form is:

where ({U}_{in}) is the sum area of other types converted to this type in study period T, and ({U}_{out}) is the sum area of a certain type converted to other types in study period T.

Calculation of ecosystem service value

In this paper, the improved ecosystem services valuation method based on the unit area value equivalent factor proposed by Xie et al.² is employed to evaluate ESV. Therefore, the LULC types of the main urban area of Guangzhou were associated to the corresponding representative biomes (Table 2). The most representative biomes used as a proxy for each LULC type are: (1) cropland for cultivated land, (2) tropical forest for forest, (3) grassland for grassland, and (4) water system and wetland for water body. (5) bare land for bare land.

The LULC types are not exactly identical with the representative biomes. For example, cultivated land in this study accounts areas used for paddy fields and dry land. Therefore, the ESV index of cultivated land is calculated by weighting the area ratio of paddy field and dry land in the statistical yearbook. The average value of the ESV index of the water system and wetland is adopted for water body. The ESV index of build-up area is 0. the ESV index of land use types is shown in Table 2. The value of the universal equivalent factor ESV (D value) of the improved ESV method is 340.65 thousand yuan/km².

In this study, the ESV of each land use unit area is based on the research methods of Costanza³ and Xie et al.². The calculation formula is:

where ESV is the total value (yuan) of ecosystem services, ({A}_{i}) (km²) is the area of class i, and VC_i is the ESV coefficient (yuan/km· year) corresponding to class i. ESV_f is the single ESV, and VC_fi is the value coefficient of the single service function.

Sensitivity analysis model

The sensitivity model is used to calculate the response of ESV to the change of value coefficient (VC)¹⁴ by adjusting the 50% of the ESV coefficient of each land use type up and down and determining the change in ESV over time and the degree of dependence on the value coefficient. The calculation formula is as follows:

where ESV is the estimated ESV, VC is the value coefficient, i and j are the initial value and adjusted value (50% up and down adjustment) and K is the land use type. When CS ≥ 1, ESV is elastic relative to VC; when CS ≤ 1, ESV is inelastic. The larger the CS value is, the more critical the accuracy of the ESV index.

Grid analysis method

Grid analysis method is used to divided the study area into regular grid matrixes with the same size and no overlap, and takes grid as the research object to express and statistical unit in geospatial space⁹. It uses regular square grid as ESV’s spatial statistical unit instead of irregular land-use map spots to ensure the capacity invariance within the unit, which not only highlights the spatial distribution characteristics of ESV, but also facilitates the spatial quantitative statistics of ESV.

The premise of calculating ESV spatial differentiation is to determine the size of grid cell. In this study, based on ArcGIS10.5 software, the area of each land use type in four grid units with side length of 0.5 km, 1 km, 2 km and 3 km was extracted respectively, and then compare the degree of area change, namely the coefficient of variation. The grid of this scale is the optimal grid unit size for the spatial differentiation of ESV in the study area. Finally, the grid analysis method is introduced to construct a (0.5 × 0.5) km square grid as a geospatial statistical unit. By using the equal spacing system sampling method, the study area is divided into 2024 square grids that do not overlap each other (0.5 × 0.5) km, and the grid matrix is composed of these grids. Through the intersection operation of grid matrix and land use data of each research year, the area of various land use types in each grid is counted, and the spatial heterogeneity of land use types and ESV is analyzed.

Geo-detector

Combined with GIS spatial superposition technology and set theory, a statistical method proposed by Wang et al.³² can detect spatial heterogeneity and reveal the driving force to identify the interaction between multiple factors. This model is widely used to analyze the influence mechanism of social economic factors and natural environmental factors. The geo-detector consists of four detectors—a differentiation and factor detector, an interaction detector, a risk area detector and an ecological detector. In this paper, factor detection and interaction detection are utilized to detect and analyze the driving force of ESVs in the main urban area of Guangzhou.

(1) Factor detector: Factor detection can identify the explanatory power of each spatial driving factor in landscape type change, and its model is as follows:

where P is the explanatory power index of ESV influencing factors; n_i is the number of samples in the secondary area; n is the total number of samples; m is the number of samples in the secondary area; ({delta }^{2}) is the variance in land use type change in the whole area; and ({{delta }^{2}}_{i}) is the variance in land use type in the secondary area. Thus, the model is established, assuming ({{delta }^{2}}_{i}) ≠ 0.

The range of values for P is [0,1]. When P = 0, it shows that the spatial distribution of ESV changes is random. The larger the P value is, the greater the impact of longitudinal driving factors on ESV changes.

(2) Interactive detector: Interaction detection can be used to identify the interaction between different spatial drivers. When detecting the interaction of X1 and X2, the explanatory power of the dependent variable Y will increase or decrease. The evaluation method is used to calculate the q value of two factors, X1 and X2, for Y, respectively, q (X1) and q (X2); to calculate the q value of their interaction, q (X1 ∩ X2); and to compare q (X1), q (X2) and q (X1 ∩ X2). The five results of the interactive detector are given in Table 3.

Source: Ecology - nature.com