*Article* **Analysis of the Spatiotemporal Changes in Watershed Landscape Pattern and Its Influencing Factors in Rapidly Urbanizing Areas Using Satellite Data**

**Zhenjie Zhu <sup>1</sup> , Bingjun Liu 2,3,\* , Hailong Wang 2,3 and Maochuan Hu 2,3**


**Abstract:** Analyzing the spatiotemporal characteristics and causes of landscape pattern changes in watersheds around big cities is essential for understanding the ecological consequence of urbanization and provides a basic reference for the watershed management. This study used a land-use transition matrix and landscape indices to explore the spatiotemporal change of land use and landscape pattern over Liuxihe River basin of Guangzhou in the southeast of China from 1980 to 2015 with multitemporal Landsat satellite data in response to the rapid urbanization process. Primary temporal and spatial influencing factors were first quantitatively identified through grey relation analysis (calculating correlation degree between land use changes and influencing factors) and Geodetector (detecting landscape spatial heterogeneity and its driving factors), respectively. Considerable spatial and temporal differences in land use and landscape pattern changes were observed herein, thus determining the influencing factors of these differences in the Liuxihe River basin. These changes were characterized by a large increase in construction land converted from cropland, particularly in the middle and lower reaches of the basin from 2000 to 2010, causing dramatic fragmentation and homogenization of the landscape pattern there. Meanwhile, the landscape pattern gradually transitioned from an agricultural land use dominant landscape to a construction land use dominant landscape in these regions. Furthermore, the rapid growth of a nonagricultural population and the transformation of industry primarily caused the temporal changes of landscape pattern, and the landscape spatial heterogeneity was mainly caused by the interaction of complicated geomorphology and anthropogenic activities in different spatial locations, particularly after 2000. This study not only provides an improved approach to quantifying the main spatiotemporal influencing factors of landscape pattern changes during different time periods, but also offers a reference for decisionmakers to formulate optimal strategies on ecological protection and urban sustainable development of different regions in this study area.

**Keywords:** landscape pattern; spatiotemporal changes; influencing factors; watershed; China SE; satellite data

#### **1. Introduction**

The increasing expansion of big cities has been a common social and economic phenomenon taking place all around the world, especially in the developing countries since the 21st century [1,2]. This process, with no sign of slowing down, may be the most critical anthropogenic force that has brought about dramatic changes in land use and landscape pattern at local, regional, and global scales [3–5]. Numerous studies have found that these immense changes can not only contribute to various environmental issues [6,7], but

**Citation:** Zhu, Z.; Liu, B.; Wang, H.; Hu, M. Analysis of the Spatiotemporal Changes in Watershed Landscape Pattern and Its Influencing Factors in Rapidly Urbanizing Areas Using Satellite Data. *Remote Sens.* **2021**, *13*, 1168. https://doi.org/10.3390/rs13061168

Academic Editor: Francesca Cigna

Received: 14 January 2021 Accepted: 16 March 2021 Published: 18 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

also affect the structure, function, and health of the ecosystem [8,9], and further threaten the sustainable development of big cities [10,11]. Therein, watersheds around rapidly urbanizing areas are more sensitive to these changes due to its richer and fragile natural ecosystems [12]. Moreover, a watershed is a complete natural and unnatural circulation unit, which is more conducive to conduct the ecological protection and restoration. The problem related to landscape pattern changes in these kinds of watersheds has been receiving more and more attention from international scholars in recent decades [13,14]. Therefore, gaining a deep understanding of the processes and causes of landscape pattern changes is crucial for protection, management, and sustainable planning of these areas under rapid urbanization [15,16].

Previous studies have illustrated that the analysis of land use changes is usually regarded as the basis for studying the landscape patterns change [17], because the landscape pattern is usually defined as the spatial arrangement of various landscape patches of different types, sizes, and shapes, which are classified by different land use types [18]. Changes in the landscape pattern were proved to be the results of changes in various land use types [19]. Most scholars choose a land use transition matrix to reflect the mutual transformation characteristics between any two different land use types [20,21], and use landscape metrics to detect the characteristics of spatial-structural composition and configuration in different landscape patches [22,23]. Therein, the former emphasizes changes of land surface properties in different periods [24], and the latter stresses the changes of potential ecological pattern [25]. When studying the changing characteristics of landscape pattern, it is necessary to analyze land use changes first, and emphasize both the temporal and spatial changes of them.

Currently, previous studies on the change of watershed landscape pattern in rapidly urbanizing areas seldom quantified spatiotemporal processes and causes of landscape pattern changes comprehensively. For example, Su et al. [26] analyzed the land use and landscape pattern in a different period to reflect its spatiotemporal changes characteristics and its causes, but ignored the overall spatial heterogeneity of landscape pattern. Zhang et al. [27] and Shi et al. [28] systematically analyzed the spatiotemporal changes processes of land use and landscape pattern in watersheds. The former study only quantified the temporal influencing factors of land use changes, and the latter study described the temporal and spatial influencing factors respectively, but not quantitatively. In addition, when analyzing the influencing factors of landscape pattern changes quantitatively, many studies failed to solve the problem of insufficient multi-temporal land use data [29–31], nor did they consider the interactive effects of different factors on the spatial landscape heterogeneity [32,33]. Some studies even analyzed the transition driving forces of different land use types in different locations of different period to meet the requirements of large quantities data for commonly used analysis methods [34]. But the causes of the spatial characteristics of landscape pattern were ignored in this case. Considering the fact that the analysis of both spatial and temporal causes were all important for guiding the management and protection of the natural ecosystem in watershed [35]. There is no doubt that researchers should carry out a systematic analysis on the process and cause of watershed landscape pattern spatiotemporal changes quantitatively.

With the continuous improvement of the remote sensing technology, more and more remote sensing data with different sensors, time periods, spectrum and spatial resolution can be acquired [36]. It provides the most stable and accurate multi-temporal data source for land use analysis, thus the land use and landscape pattern changes under the rapid urbanization can be monitored and analyzed spatiotemporally [37]. There existed many studies that using these kinds of data analyzed the land use change of different regions in Guangzhou city [38,39], and compared the relationships between these changes and urban expansion [40,41]. Liu et al. [31] also discussed the land use change and its causes based on the Landsat satellite data with remote sensing technology. These researches all illustrated the feasibility of using land use data provided by remote sensing technology to analyze the changes of landscape patterns in Guangzhou city. Besides, along with the rapidly urban expansion of Guangzhou city, the intensive interaction between natural and human elements in the LXH has brought about the large transition of construction land from lots of forest and cropland [42]. The water quality and natural environment of the LXH were degraded, especially in the downstream [43,44]. The Liuxihe River basin (LXH), as the final ecological barrier on the northwest of Guangzhou city, is an important water resource conservation area. Therefore, compared with areas divided by the administrative units in Guangzhou city, analyzing the landscape pattern changes in the LXH has greater ecological significance and protection value for Guangzhou.

This study mainly analyzes the spatiotemporal differences of process and causes on landscape pattern changes under rapid urbanization of the Guangzhou city. The specific objectives are: (1) To analyze the changes of different land use types in time and space; (2) to characterize the spatial configuration of the landscape pattern in time and space; (3) to establish the relationship between changes in different land use types and different influencing factors temporally; and (4) to determine the impact of different factors on the landscape heterogeneity. Thus, these differences in the spatiotemporal changes and causes of landscape pattern in the LXH under the rapid urbanization since 1980 are revealed. The decision-makers will be more clear about how to formulate an appropriate strategy for planning and management.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Liuxihe River basin is in the rapidly developing and urbanizing city of Guangzhou, southeast of China (Figure 1). The river, about 171 km in length with an area of 2300 km<sup>2</sup> , flows through Guangzhou, and eventually empties into the Beijiang River, a tributary of the Pearl River. The annual precipitation rate over the LXH is 1750 mm and more than 80% of precipitation occurs from April to September. Its daily mean air temperature is about 20 ◦C and annual rate of evaporation is about 1200 mm [45]. The elevation of the LXH falls gradually from northeast to southwest, characterized by mountains in the upstream, hills in the midstream, and plains in the downstream, thus making the upstream more difficult to develop than the middle and lower watershed. At present, the distribution of land use in the LXH is characterized by forests in the upstream, croplands and orchards in the midstream, and construction lands in the downstream. In addition, the speed of the urbanization process in the midstream and downstream of the LXH is faster and stronger than in the upstream [46]. The special location conditions, the unbalancing effect of urbanization, the difference in natural topography, and the land uses/land covers all make the LXH an appropriate case for the research processes and influencing factors for spatiotemporal changes of watershed landscape pattern in rapidly urbanizing areas.

#### *2.2. Data and Data Processing*

Satellite images taken in 1980, 1990, 1995, 2000, 2005, 2008, 2010, and 2015 (Landsat MSS/TM/ETM, Landsat 8) with a spatial resolution of 30 m were provided by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http: //www.resdc.cn, accessed on 30 December 2017). The method of visual interpretation was used to derive thematic land use maps based on the land resource classification system of Chinese academy of sciences [47]. Meanwhile, the accuracy of interpretation was improved through reference data, such as geomorphic maps, vegetation maps, ground truth data at different sample points, and local resident interview data. The calculated results of Kappa coefficient were larger than 0.80, which verified the accuracy and reliability of these land use maps [48]. This paper reclassified these land use maps into nine types including cropland, forest, orchard, grassland, shrub, water, floodplain, construction land, and unused land (Appendix A Table A1).

Besides, data from 1980 to 2015 about demographic factors, socioeconomic factors, and urbanized activities of the LXH were collected from the Guangzhou Statistical Yearbook (http://tjj.gz.gov.cn/, accessed on 1 April 2018) and the Outline of Guangzhou Urban

Construction Overall Strategic Concept Plan (http://ghzyj.gz.gov.cn/, accessed on 30 October 2018), which is provided by the Chinese government. The detailed factors are total population (TP), proportion of non-agricultural population (PNAP), gross domestic product (GDP), proportion of primary industry (PPI), proportion of secondary industry (PSI), proportion of tertiary industry (PTI), annual per capital income (APCI), and total investment in real estate development (IRE).

**Figure 1.** The location and digital elevation map (DEM) of the Liuxihe River basin (LXH).

Spatial data such as topographical elements (digital elevation map (DEM) and SLOPE) were obtained from NASA's Earth Observing System Data and Information System (https://search.earthdata.nasa.gov/search, accessed on 30 December 2018). And other spatial data were obtained from Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (Source: China, Satellite images, http://www.resdc.cn, accessed on 30 December 2018), including population density (TP), socioeconomic level (GDP), urbanized activities (NLD), and land use intensity (LUIN). These main parameters of the LXH with a spatial resolution of 1000 m (except for the DEM and SLOPE data with a spatial resolution of 30 m) were all cut out by the ArcGIS software.

#### *2.3. Methods*

2.3.1. Selection of Landscape Metrics and Influencing Factors

• Landscape Metrics

The factor analysis method [49] was taken to screen the notable landscape metrics from 23 frequently used metrics (Appendix B Figure A1). First, if the absolute value of the correlation coefficient between two indices was more than 0.9, only one was used; second, indices representing different aspects of landscape characteristics were selected to reduce the information redundancy among them [50]. Five representative indices were used in the final analysis, including the patch density (PD), aggregation index (AI), largest patch index (LPI), area-weighted mean patch fractal dimension (AWMPFD), and Shannon's diversity index (SHDI). These landscape metrics were calculated for 1980, 1990, 2000, 2010, and 2015 using the public domain software FRAGSTATS 4.2 at both the class-level and landscape-level in 14 different granularities (30 m, 50 m, 100 m, 200 m, 300 m, 400 m, 500 m, 600 m, 700 m, 800 m, 900 m, 1000 m, 1100 m, 1200 m). FRAGSTATS is a computer software program designed to compute a wide variety of landscape metrics for

categorical map patterns at different levels (the detailed introduction of this software can be found at https://www.umass.edu/landeco/research/fragstats/fragstats.html, accessed on 1 February 2018). Thus, the characteristic scale interval which is appropriate for spatial analysis of the Liuxihe River basin was determined (Appendix B Figure A2). The calculation formula and ecological significance of these indices are given in Table 1.



Note: i = 1 . . . m patch types (classes); j = 1 . . . n patches; A = total area of each landscape type (m<sup>2</sup> ); aij = area (m<sup>2</sup> ) of patch ij; Pij = perimeter (m) of patch ij; N<sup>i</sup> = number of patches in the landscape of patch type (class) i; m = number of patch types (classes) present in the landscape, excluding the landscape border if present; P<sup>i</sup> = proportion of the landscape occupied by patch type (class) i.

• Land Use Transition Matrix

The changes in landscapes were detected by calculating each land use type transition matrix of any two adjacent periods from 1980 to 2015. The following equation (Equation (1)) was applied to calculate the matrix:

$$\mathbf{p} = \begin{bmatrix} \mathbf{p}\_{11} & \mathbf{p}\_{12} & \cdots & \mathbf{p}\_{1j} \\ \mathbf{p}\_{21} & \mathbf{p}\_{22} & \cdots & \mathbf{p}\_{2j} \\ \vdots & \vdots & \vdots & \vdots \\ \mathbf{p}\_{i1} & \mathbf{p}\_{i2} & \cdots & \mathbf{p}\_{ij} \end{bmatrix} \tag{1}$$

j=1

where pij indicates the area in transition from landscape i to j. Each element of the transition matrix meets two standards: (1) pij is non-negative, and (2) n ∑ pij = 1.

For a better characterization of landscape changes, the transition matrix between any two adjacent periods was displayed in a two-dimensional table by many researchers [51]. Therein, the diagonal entries of the table reflect the total size of persistent land use types whereas the off-diagonal entries show the transition size of one land use type to another. Besides, the gross gain and gross loss of each land use type were also displayed in the table, which could be used to calculate the total exchange (sum of the gross gain and the gross loss) and net change (the gross gain minus the gross loss).

• Selection of Landscape Change Influencing Factors

The influencing factors on the development of landscape pattern vary at different spatial and temporal scales [52], which were categorized into natural and anthropogenic factors in most of the related researches [53,54]. This study focused on the small-scale watershed around rapidly urbanization city in a short research period. In these areas, anthropogenic factors such as population growth, socioeconomic activities, urbanization activities, and related policies usually play a major role [29]. In terms of natural factors, they are relatively stable and unchanged in short term compared with anthropogenic factors [55], thus their impact on landscape changes can be ignored in this case [56–58]. Finally, only anthropogenic factors were selected to reflect the temporal influencing factors in this study, including the total population (TP), the proportion of the non-agricultural population (PNAP), gross domestic product (GDP), the proportion of the primary industry (PPI), the proportion of the secondary industry (PSI), the proportion of the tertiary industry (PTI), annual per capita income (APCI), and total investment in real estate development (IRE). These factors are about the demographic, socioeconomic, and urbanized activities of a region. Besides these, the policies related to the rapid urbanization were also considered.

In addition, researchers found that climate conditions had no significant influence on the spatial distribution of landscape in a small-scale catchment [33], and topographical and anthropogenic factors were commonly used to interpret the landscape spatial heterogeneity in this case [28,59]. As we all know, the spatial difference of soil and hydrological conditions always depend on topographical factors. Therefore, this study mainly analyzed the spatial difference of topographical and anthropogenic factors on landscape spatial heterogeneity. Herein, the spatially distributed data of GDP and TP were selected to represent the spatial difference between socioeconomic level and population density. The nighttime light (NLD) spatial data were used to represent the urbanization degree. Moreover, the spatial distribution of the LUIN was calculated using the land use comprehensive degree index, representing the urbanization activities related to land use changes. Additional details on the calculation of land use comprehensive degree index can be found in Yu et al. [60]. Besides, the topographical variables (DEM and SLOPE) were also considered.

#### 2.3.2. Quantifying the Influence on the Change of Landscape Pattern

Considering the spatiotemporal characteristics of land use and landscape pattern changes and their interactions, this study mainly analyzed the factors that influence the temporal change of land use types and spatial differences of landscape patterns (i.e., the spatial heterogeneity of landscape). Since the data for each year could not be obtained, it made the commonly used quantitative statistical analysis impractical. This study applied grey correlation analysis to explore the impact of different influencing factors on land use changes in an attempt to effectively solve the problem of insufficient sample data. Meanwhile, spatial correlation analysis using the Geodetector was carried out to assess the influence degree and interaction of each influencing factor toward different spatial characteristics of landscape index.

• Grey Correlation Analysis

Grey correlation analysis is an impacting factor measurement method in the grey system theory proposed by Deng [61] in 1982, which analyzes the uncertain relationship between a main factor and all other influencing factors in a given system. It can complement the defects of statistical analysis methods and can work with small amounts of irregular data; it also negates the inconsistency between quantitative and qualitative results. The analysis method mainly compares the time series of each influencing factor to determine which one is dominant. That is, when the trend of changes between an independent variable and a dependent variable is consistent or the degree of synchronization change is high, a strong correlation results [62]. The relationship is often expressed by grey correlation degree (Equation (2)). The greater the degree of grey correlation is, the more the influence degree of the factor will be and vice versa.

$$\gamma\_{\rm ij} = \frac{1}{n} \sum\_{\mathbf{k}=1}^{n} \frac{\min\_{\mathbf{i}} \min\_{\mathbf{k}} \Delta\_{\mathbf{i}}(\mathbf{k}) + \xi \max\_{\mathbf{i}} \max\_{\mathbf{k}} \Delta\_{\mathbf{i}}(\mathbf{k})}{\Delta\_{\mathbf{i}}(\mathbf{k}) + \xi \max\_{\mathbf{i}} \max\_{\mathbf{k}} \Delta\_{\mathbf{i}}(\mathbf{k})}, \Delta\_{\mathbf{i}}(\mathbf{k}) = \left| \mathbf{x}\_{\mathbf{j}}(\mathbf{k}) - \mathbf{x}\_{\mathbf{i}}(\mathbf{k}) \right| \tag{2}$$

where x<sup>i</sup> and x<sup>j</sup> are the independent variable series and the dependent variable series, respectively; γij is the grey relational degree between independent variables x<sup>i</sup> and dependent variables x<sup>j</sup> ; ξ is the resolution coefficient, ξ ∈ (0, 1), and usually ξ takes a value of 0.5; k = 1, 2, · · · , n is the time series.
