Estimating the Characteristic Spatiotemporal Variation in Habitat Quality Using the InVEST Model—A Case Study from Guangdong–Hong Kong–Macao Greater Bay Area

Abstract

The intensity of human activity, habitat loss and habitat degradation have significant impacts on biodiversity. Habitat quality plays an important role in spatial dynamics when evaluating fragmented landscapes and the effectiveness of biodiversity conservation. This study aimed to evaluate the status and characteristic variation in habitat quality to analyze the underlying factors affecting habitat quality in the Guangdong–Hong Kong–Macao Greater Bay Area (GBA). Here, we applied Kendall’s rank correlation method to calculate the sensitivity of habitat types to threat factors for the Integrated Valuation of Ecosystem Services and Tradeoffs habitat quality (InVEST-HQ) model. The spatiotemporal variation in habitat quality of the GBA in the period 1995–2015 was estimated based on the InVEST-HQ model. We analyzed the characteristic habitat quality using different ecosystem classifications and at different elevation gradients. Fractional vegetation cover, the proportion of impervious surface, population distribution and gross domestic product were included as the effect factors for habitat quality. The correlation between the effect factors and habitat quality was analyzed using Pearson’s correlation tests. The results showed that the spatial pattern of habitat quality decreased from fringe areas to central areas in the GBA, that the forest ecosystem had the highest value of habitat quality, and that habitat quality increased with elevation. In the period from 1995 to 2015, habitat quality declined markedly and this could be related to vegetation loss, land use change and intensity of human activity. Built-up land expansion and forest land fragmentation were clear markers of land use change. This study has great significance as an operational approach to mitigating the tradeoff between natural environment conservation and rapid economic development.

1. Introduction

With the increasing intensity of human activity (e.g., urbanization and expansion of agricultural areas), species populations and richness have been changing over time in complex ways [1,2]. Habitat loss and degradation are the largest threats to wildlife and biodiversity [3], particularly the process of urbanization and land use change in developing countries [4,5,6]. Biodiversity supports ecosystem functioning and underpins an ecosystem’s stability in the face of global environment change [7,8,9]. However, the number of species under threat is increasing, and populations of vulnerable taxa are declining, resulting in profound changes to ecosystem function [10,11,12].

Habitat quality refers to the capability of an ecosystem to provide the necessary resources and conditions for all its wildlife or specific populations [13,14]. The study of habitat quality plays an important role in extracting the spatial dynamics of fragmented landscapes and assessing the effectiveness of biodiversity conservation measures [15]. Habitat patches in a landscape will influence wildlife population viability and habitat occupancy in a landscape will ultimately manifest in population dynamics, which highlights the importance of evaluating habitat quality on a large scale [16,17,18,19].

There are two conceptual approaches to assessing habitat quality:

• An indirect approach, which reveals variation in habitat quality by measuring variables for certain species and their populations in different habitats [18]. This approach could be implemented by direct field investigation [9,20] or species distribution modeling. Field investigation could acquire veracious species distribution data and population information, but usually consumes significant human and material resources. Species distribution models, which can be generated using maximum entropy model (MaxEnt) [21] and bioclimatic data [22], use known occupancy data and environmental variables to evaluate the habitat suitability of other nonvisited areas to predict the potential species distribution [23]. This type of model can help describe the relationship between habitat selection and the environment variable, but it requires species occurrence data.
• Another approach for assessing habitat quality is to measure attributes of a habitat directly, such as critical resources and the ecological constraints that could limit the use of resources [18]. Common ways to carry out this approach include expert-based models and ecological process models. An expert-based model can reflect the ecological situation of a study area or the necessary resources for specific populations [24]. However, this approach relies on expert knowledge to select an evaluation index to establish ecological indicators with which to evaluate habitat status. The ecological process model, a simplified version of the ecological process, emphasizes the threat of human activities to habitat quality. Examples include the integrated valuation of environmental services and tradeoffs habitat quality (InVEST-HQ) model [25] and the global biodiversity model (GLOBIO) model [26]. These models have a complete evaluation system to assess habitat quality which can reduce randomness in terms of the selection of evaluation indices. With a consideration of ecological process, these models could provide a more scientific theoretical foundation in the assessment of habitat quality.

In particular, the InVEST-HQ model can estimate habitat quality based on land use/cover data and habitat threat data and provides a concise approach to evaluating the status of habitat quality when there are limited available data and an unsampled area [14,25]. It can be used to map and quantify a numeric value for habitat quality and can be used as a proxy for regional biodiversity to some extent [27,28,29]. A limitation of the InVEST-HQ model is that the parameters rely on empirical values, an issue that requires more attention. The InVEST-HQ model has been applied to assess habitat quality in relation to different anthropogenic threats in many areas. For instance, the model has been used for habitat quality studies that assessed the effects of land use intensity and population density in the Nansihu lake basin [30], urbanization and dam construction in the Llobregat River basin [14], urbanization and landscape patterns in Hangzhou, China [27], built-up areas and water pollution in the Pochote River in León, Nicaragua [31], and urbanization expansion in the Beijing–Tianjing–Hebei region [32].

The Guangdong–Hong Kong–Macao Greater Bay Area (GBA) is one of the most densely populated and economically active regions in the world [33]. In this region, rapid change in land use/cover and intense human activity has led to a large amount of nature capital consumption and increased pressure on ecosystems [34]. It is vital for the sustainable development of urban agglomerations to maintain the biodiversity and ecosystem services in the GBA during rapid economic and urban development. Thus, evaluating and monitoring the status and dynamic changes in habitat quality under the impact of land use change and human activity has great significance in mitigating the tradeoff between natural environment conservation and rapid economic development. There is still a lack of knowledge regarding the characteristic differentiation of habitat quality, and the relationship between habitat quality and factors affecting it are not clear in the GBA.

In this study, we applied the Kendall’s rank correlation method to calculate the parameters of the InVEST-HQ model, then used the InVEST-HQ model to estimate the spatiotemporal variation in habitat quality over the period 1995–2015. We used different ecosystem classifications and a digital elevation model (DEM) to analyze the characteristics of spatial differentiation of habitat quality in different perspective. Finally, fractional vegetation cover, the proportion of impervious surface, population distribution and gross domestic product were considered as factors affecting habitat quality. The underlying correlation between habitat quality and these factors was analyzed using the Pearson correlation coefficient. The main objectives of this work were to (1) acquire suitable parameters for the sensitivity of land use/cover (LULC) types to the threat factors in the GBA to increase the objectivity of the InVEST-HQ model; (2) assess the spatiotemporal variation in habitat quality in the period 1995–2015; (3) identify the characteristics of spatial differentiation of habitat quality from ecosystem and elevational gradient perspectives; (4) analyze the underlying correlations between the influencing factors and habitat quality.

In this study we assessed habitat quality, which can be used to evaluate the availability of management strategies for biodiversity conservation, and we provide insights for exploring the tradeoffs between natural habitat and economic development in urban agglomerations.

2. Materials and Methods

2.1. Study Area

The Guangdong–Hong Kong–Macao Greater Bay Area (GBA) is located in South China (21–25°N, 111–116°E) with an elevation of between −17 and 1828 m above sea level (Figure 1). The area includes nine cities in Guangdong Province (Guangzhou, Shenzhen, Zhuhai, Jiangmen, Zhaoqing, Foshan, Dongguan, Huizhou and Zhongshan) and two special administrative regions of China (Hong Kong and Macau). The area of the GBA covers approximately 5.6 × 10⁴ km² has several representative geomorphic types including hills, platforms, monadnocks, and plains [35]. The GBA has a subtropical monsoon climate with mild winters and long, hot and humid summers. The average annual precipitation is 1889 mm and the annual average maximum and minimum temperature are 26.4 and 18.9 °C, respectively.

The GBA provides a suitable habitat for a wide diversity of species due to the diversity of its ecosystem types and resources such as abundant wetlands and water [36]. Among the land cover types of the GBA in 2015 were forest land (54% cover), cultivated land (23%), built-up areas (14%), and built-up areas under development from 1995 to 2015 (5%).

2.2. Data Resources and Preparation

For this study we used land use/cover datasets, vector datasets (i.e., road network datasets, residential datasets and industry datasets), socioeconomic datasets (i.e., population distribution (POP) data and gross domestic product (GDP) data), impervious surface data, normalized difference vegetation index (NDVI) data, digital elevation model (DEM) data and terrestrial ecosystem classification data. All datasets in this study are described in

Table 1. The datasets and resources used in this study.

No.	Datasets	Data Description	Data Resources	Data Format
1	Land use/cover datasets	China land use/cover change data	Obtained from Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC), http://www.resdc.cn, DataID = 184 (accessed on 20 January 2021)	TIFF
2	Vector datasets	Road network datasets, residential datasets and industry datasets obtained	Obtained from OpenStreetMap (OSM), https://download.geofabrik.de/ (accessed on 20 January 2021)	Shapefile
3	Socioeconomic datasets	Population distribution (POP) data	Obtained from Oak Ridge National Laborator LandScan global population distribution data, https://landscan.ornl.gov/landscan-datasets (accessed on 20 January 2021)	TIFF
		Gross domestic product (GDP) data	Obtained from RESDC 1 km grid GDP data, http://www.resdc.cn/data.aspx?DATAID=252 (accessed on 20 January 2021)	TIFF
4	Impervious surface data	Global artificial impervious data	Obtained from global artificial impervious area (GAIA) maps, http://data.ess.tsinghua.edu.cn/gaia.html (accessed on 20 January 2021)	TIFF
5	Normalized difference vegetation index data	MOD13Q1 16-day products and RESDC NDVI dataset	Obtained from Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center (LAADS DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 20 January 2021)) and RESDC (http://www.resdc.cn/DOI/doi.aspx?DOIid=49 (accessed on 20 January 2021))	TIFF
6	Supplementary analysis datasets	DEM data	the Shuttle Radar Topography Mission (SRTM), http://srtm.csi.cgiar.org/srtmdata/ (accessed on 20 January 2021)	TIFF
		Terrestrial ecosystem classification data	Obtained from RESDC, http://www.resdc.cn/DOI/doi.aspx?DOIid=103 (accessed on 20 January 2021)	TIFF

.

Land use/cover datasets for the years 1995, 2000, 2005, 2010 and 2015 used China land use/cover change (CNLUCC) data provided by the Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC; http://www.resdc.cn, DataID = 184; accessed on 24 July 2020). CNLUCC data based on Landsat series satellite remote sensing images were processed using the supervised classification method to divide land cover/use categories with comprehensive accuracy exceeding 90% and a fine spatial resolution of 30 m [37]. In this study, CNLUCC data were reclassified into seven categories: cropland, forest, grassland, water bodies, wetland, built-up areas and unused land. The land use/cover datasets are shown in Figure 2.

Vector datasets of the road network, residential areas and industrial areas, were obtained from OpenStreetMap (OSM; https://download.geofabrik.de/, accessed on 12 July 2020). We selected road network datasets with railways, primary roads, trunk roads and secondary roads from the attribution of the road layer (fclass) in OSM data. Industry and residential datasets were extracted from the attribution of the land use layer in OSM data. The vector datasets are shown in Figure 3.

Socioeconomic data included population distribution (POP) data and gross domestic product (GDP) data. Population data for the years 2000–2015 were derived from Oak Ridge National Laboratory (ORNL) LandScan global population distribution data at approximately 1 km (30″ × 30″) spatial resolution (https://landscan.ornl.gov/landscan-datasets, accessed on 20 October 2020). This database uses the best available demographic (census) and geographic data and remote sensing imagery techniques within a multivariate asymmetric model to acquire a reliable population distribution data source; the POP data are shown in Figure A1 in Appendix A. POP data in 1995 were obtained from RESDC (http://www.resdc.cn/data.aspx?DATAID=251, accessed on 21 July 2020) with 1 km spatial resolution. The 1 km grid GDP data we used covered the period from 1995 to 2015 and were provided by RESDC (http://www.resdc.cn/data.aspx?DATAID=252, accessed on 22 July 2020). The GDP data are shown in Figure A2.

Impervious surface data were acquired from global artificial impervious area (GAIA) maps (http://data.ess.tsinghua.edu.cn/gaia.html, accessed on 2 November 2020), which provide highly accurate annual GAIA from 1985 to 2018 with 30 m spatial resolution [38]. We extracted the proportion of impervious surface (PIS) using 300 m grids in the years 1995, 2000, 2005, 2010 and 2015, with values ranging from 0 to 1. The PIS is shown in Figure A3.

NDVI datasets were derived at the spatial resolution of 250 m in MOD13Q1 16-day products from LAADS DAAC (https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 15 October 2020), with data from 2000 onwards. We calculated the mean value of NDVI in each year and then estimated fractional vegetation cover (FVC) based on the pixel dichotomy model. NDVI data for 1995 were obtained from the RESDC website (http://www.resdc.cn/DOI/doi.aspx?DOIid=49, accessed on 25 October 2020), which provides NDVI data from 1998 to 2018. We also estimated FVC and conducted correlation analysis using the 1998 NDVI instead of 1995 NDVI. The fractional vegetation cover (FVC) datasets are shown in Figure A4.

Supplementary analysis datasets in this study included DEM data and terrestrial ecosystem classification data. Terrestrial ecosystem classification data of China were obtained from RESDC (http://www.resdc.cn/data.aspx?DATAID=103, accessed on 28 July 2020) and had seven classifications: farmland, forest, grassland, water and wetland, settlement, desert and other ecosystems. These classification datasets covered the years 1995, 2000, 2005, 2010 and 2015 and were based on 1:100,000 scale land use/cover data. The spatial pattern of the terrestrial ecosystem classification is shown in Figure A5. DEM data at 30 m spatial resolution captured by the Shuttle Radar Topography Mission (SRTM) were downloaded from the Consortium for Spatial Information (CGIAR-CSI; http://srtm.csi.cgiar.org/srtmdata/, accessed on 25 July 2020). The terrestrial ecosystem classification data and DEM data were applied to extract the values of habitat quality, which had been identified using the characteristics of differentiation in habitat quality to extract numerical distribution of habitat quality from different perspectives.

Albers equal-area conic projection was used for all datasets in this study. The central meridian was 105 and the first standard parallel and second standard parallel were 25 and 47, respectively.

2.3. Methods

The analysis workflow of this study consisted of four key steps. First, the original parameter of InVEST-HQ—sensitivity of LULC to a threat factor—was calculated by Kendall’s correlation instead of the empirical value. Second, we estimated the habitat quality of the GBA from 1995 to 2015 based on the optimized InVEST-HQ model. Third, we analyzed spatial differentiation in habitat from different perspectives based on different ecosystem classifications and elevational gradients. Fourth, Pearson’s correlations were used to analyze the relationship between the effect factors and habitat quality. The organizational flow chart of this study is shown in Figure 4.

2.3.1. The InVEST Model

The InVEST model is an open-source software model used to visualize and estimate goods and services from nature on a spatial scale; it was developed by Stanford University, Minnesota University, the World Wide Fund and the Nature Conservancy [25]. Habitat quality is one of the modules of the InVEST model which can rapidly estimate the status of habitat for each LULC scale. The InVEST habitat quality (InVEST-HQ) model has been used to assess historic and dynamic changes in habitat quality in Italy [4], China [29,30] and Nicaragua [31]. We used the InVEST-HQ model (InVEST v3.8.7, Stanford University, Stanford, CA, USA) to estimate the habitat quality of the GBA in 1995, 2000, 2005, 2010 and 2015 (model available at https://naturalcapitalproject.stanford.edu/software/invest, accessed on 25 May 2020). More details can be found in the InVEST user’s guide [25] and relevant literature [14,39], while the parameters of the InVEST-HQ model are described here in Section 3.

2.3.2. Kendall’s Rank Correlation Analysis Methods

Kendall’s rank correlation can be applied to measure the strength of dependence between two classified variables. In the InVEST-HQ model, the original parameter in sensitivity of habitat type to threat factor was determined by empirical values for each LULC and threat factor (both variables were categorical datasets which had been classified as 0 or 1). The absolute value of Kendall’s correlation coefficients served as the sensitivity of habitat types to threat factors in the InVEST-HQ model. All correlation coefficients were calculated using R software in this study.

$$T=\left|\tau \right|=\left|\frac{n_c-n_d}{n\left(n-1\right)/2}\right|$$

(1)

where T is the absolute value of Kendall’s correlation coefficient, $\tau$ is Kendall’s correlation coefficient, n_c is the number of concordant pairs, n_d is the number of discordant pairs, and n is the number of pairs.

2.3.3. Pearson’s Correlation Analysis Methods

Pearson’s correlation analysis was used to investigate the factors affecting habitat quality. Habitat quality is affected by multiple factors. In particular, FVC is an important ecological parameter used for monitoring vegetation cover or forest, POP and GDP are socioeconomic factors that are used to indicate the intensity of human activity, and PIS can indicate the intensity of artificial types of land use. Therefore, POP, PIS, GDP and FVC were selected as the potential effect factors.

The significance of the Pearson’s correlation coefficients were estimated using the t-tests. We calculated the Pearson’s correlation coefficients of all habitat quality datasets and all effect factor datasets from 1995, 2000, 2005, 2010 and 2015.

$$r=\frac{{{\displaystyle \sum }}_{i=0}^{N-1}\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=0}^{N-1}{\left(x_i-\overline{x}\right)}^2{{\displaystyle \sum }}_{i=0}^{N-1}{\left(y_i-\overline{y}\right)}^2}}$$

(2)

where r is the Pearson’s correlation coefficient; x_i is the matrix of habitat quality, y_i is the matrix of each effect factor, N is the number of matrix elements (x, y both have the same matrix element), and $\overline{x}, \overline{y}$ are the mean values of habitat quality and effect factor, respectively.

2.3.4. Variation in Habitat Quality Analysis

A series of habitat quality maps were acquired using the InVEST-HQ model in different years, with the habitat quality scoring between 0 and 1 in each cell. To ascertain the spatiotemporal variation in habitat quality, we referred to the categorized method of Czekajlo et al. [40] on greenness score, which had a similar application environment. This method could express variation at different levels and incorporate the current state with historic trends [40].

In this study, we calculated four comparison groups for a pair-wise comparison sequence between 1995 and 2000, 2000 and 2005, 2005 and 2010, and 2010 and 2015. The difference maps were divided into three classes (decrease (−): difference value less than 0; not significant (0): difference value equal to 0; and increase (+): difference value greater than 0). Then, the habitat quality score (0–1) was categorized into three levels using natural breaks (Jenks) in each year, namely, low (L), moderate (M) and high (H). Last, we acquired 9 variation categories in each difference value group, with each difference value class marked by three habitat quality score levels. The habitat quality variation in this study was categorized by the last year of the comparison (e.g., when comparing a pair-wise sequence between 1995 and 2000, the habitat quality score in 2000 was selected as the categorized level). Natural breaks (Jenks) is an unsupervised classification, which was also successfully applied in another study [40,41]. The categories of variation in habitat quality are shown in Figure 5.

3. Model Parameter Settings

The parameters of the InVEST-HQ model in this study include LULC maps, a raster layer of threat factors, the weight of each threat factor, the maximum distance of each threat factor, the sensitivity score of habitat types to threat factors, and the habitat suitability score. The data resources and descriptions are given in

Table 2. The parameters of the habitat quality model.

Parameter	Parameter Description	Data Resources	Data Format
LULC	Land use/cover map	CNLUCC	TIFF
Threat factor	The layer of threat factor	OSM and CNLUCC	Transferred to TIFF
Threat factor parameters	The weight of threat factor	Based on relevant literature [4,27,30,31] and InVEST user’s guide [25]	CSV
	The maximum distance of each threat factor
	The types of decay over space for the threat
Habitat suitability score	Each LULC type needs to give a habitat suitability score with a value ranging from 0 to 1	Based on relevant literature [4,27] and InVEST user’s guide [25]	CSV
Sensitivity of habitat types to threat factors	The parameter of the sensitivity score for LULC to the threat factor	The absolute value of Kendall’s correlation coefficients	CSV

. The details of the parameter settings will be given in the following sections.

3.1. Threat Factor Parameters

Defining the threat factor parameter for a habitat is a key issue in the InVEST-HQ model [27]. Threat factors include biotic factors (e.g., alien species invasion), abiotic factors (e.g., sunshine duration, temperature, precipitation) and anthropogenic disturbance (e.g., land use change, roads, rural settlements) [30]. Anthropogenic disturbance has been considered in relation to ecology and the environment in previous studies, which have focused on the impact of land use change [30,42,43], and the threat of land use change and roads [31]. In this study, we selected threat factors from anthropogenic disturbance based on the study of Sallustio et al. [4] and the actual situation of the GBA.

The GBA has experienced rapid development and change in land use/cover due to accelerated industrialization and urbanization [44]. Artificial land use/cover types (e.g., cropland, urban areas, plantations) have expanded, undermining ecosystem services along with considerable losses of biodiversity [45]. Artificial land use/cover types should be considered as threat factors, as the most direct manifestation of anthropogenic disturbance to the ecosystem [30]. The transport network has a cumulative and irreversible effect on ecosystem services and biodiversity, and its construction has resulted in habitat fragmentation and has become a severe threat [46,47,48]. We selected cropland, built-up areas and bare land from the LULC map, and road networks (including railways, primary roads, trunk roads and secondary roads), residential and industrial areas from OSM data as they are consistent with these threat factors.

The parameter of threat factors included weight, maximum distance and decay type, which were determined by empirical values or expert knowledge. We determined the weight, maximum distance, and decay type of each threat factor based on relevant studies [4,27,30,31] carried out in similar circumstances, the actual situation of the GBA, and the InVEST model user’s guide [25]. The threat factor parameters of this study are listed in

Table 3. Threat factor parameter settings.

Threat	Weight	Maximum Distance	Decay Type
Bare land	0.2	3	linear
Built-up areas	1	10	exponential
Cropland	0.68	8	linear
Railway	0.9	9	exponential
Trunk road	1	10	exponential
Primary road	1	8	linear
Secondary road	0.75	5	linear
Industry activity	1	12	exponential
Residential	0.5	5	exponential

.

The near edge of the threat factor should be considered which will influence the habitat quality result in the InVEST-HQ model [25]. In this study, we defined the buffer width of road networks, industry and residential activity, then we converted the vector buffer file into a raster surface. The buffer widths were based on Aneseyee et al. [49]: 1.5 km for trunk road, 1.5 km for railway, 1 km for primary road, 0.5 km for secondary road, 0.5 km for residential areas and 1 km for industrial areas. All maps consisted of a value of 1 or 0, where 1 indicated the presence of a threat and 0 indicated the absence of threat.

3.2. Habitat Suitability Score

In the InVEST-HQ model, the habitat suitability score is a numerical parameter based on LULC types regarded as suitable habitat. Each LULC type gives a habitat suitability score with a value ranging from 0 to 1, where a value of 1 indicates the most suitable habitat. We assumed each LULC type was suitable according to the habitat suitability score in the InVEST user’s guide [25] with a consideration for overall biodiversity (γ biodiversity) in the GBA. Built-up areas and unused land that were not suitable habitat were set at a value of 0, while forest and wetland were set the value 1; other LULC values were based on relevant literature [4,27] and the InVEST user’s guide [25]. The habitat suitability scores are shown in

Table 4. Habitat suitability of each land use/land cover (LULC) type.

LULC Types	Cropland	Forest	Grassland	Water Bodies	Wetland	Built-Up Areas	Unused Land
Habitat score	0.35	1	0.4	0.9	1	0	0

.

3.3. Sensitivity of Habitat Types to Threat Factors

The sensitivity of habitat types to threat factors was a numerical parameter given for each LULC type between each threat factor in the InVEST-HQ model, with the value ranging from 0 to 1, with 1 representing high sensitivity to a threat factor. The original parameter of the sensitivity score for LULC to the threat factor within the model were based on the empirical value [50], which generally relies on expert knowledge, meaning that the InVEST-HQ model is implicitly subjective and limited in terms of applicability to special regions.

Therefore, finding an approach which could acquire objective regional representative parameters was deemed necessary. We used Kendall’s correlation analysis to calculate the correlation coefficients between each LULC type and each threat factor in different years, in an attempt to address the limitations of the sensitivity score in the InVEST-HQ model. The method was described in Section 2.3.2 and the results are shown in

Table A1. The sensitivity score of habitat types to threat factors.

Year	LULC Types	Threat Factors
		Cropland	Built-Up Areas	Bare Land	Rail-Way	Trunk Road	Primary Road	Secondary Road	Industry Activity	Residential
1995	Cropland	1	0.17	0.01	0.04	0.08	0.1	0.08	0.04	0.04
	Forest	0.6	0.33	0.01	0.2	0.22	0.24	0.24	0.22	0.19
	Grass land	0.09	0.05	0	0.01	0.02	0.01	0.03	0.01	0.01
	Wetland	0.13	0.07	0	0.05	0.08	0.06	0.06	0.12	0.03
	Water bodies	0.08	0.05	0	0.02	0.04	0.03	0.02	0.07	0.01
	Built-up areas	0.17	1	0	0.25	0.24	0.26	0.29	0.19	0.27
	Unused land	0.01	0.01	0.35	0.01	0	0.01	0	0.02	0
2000	Cropland	1	0.17	0.03	0.05	0.08	0.11	0.09	0.04	0.05
	Forest	0.62	0.31	0.05	0.19	0.22	0.24	0.24	0.22	0.19
	Grass land	0.09	0.04	0.3	0.01	0.01	0.02	0.02	0.02	0.01
	Wetland	0.14	0.07	0.01	0.04	0.07	0.05	0.06	0.11	0.02
	Water bodies	0.09	0.04	0.01	0.02	0.04	0.02	0.01	0.06	0
	Built-up areas	0.17	1	0.01	0.24	0.24	0.25	0.28	0.19	0.27
	Unused land	0.01	0.01	0.21	0.01	0.01	0	0	0.01	0.01
2005	Cropland	1	0.14	0	0	0.05	0.08	0.04	0.02	0.02
	Forest	0.58	0.37	0.01	0.2	0.23	0.25	0.25	0.23	0.19
	Grass land	0.08	0.05	0	0.01	0.01	0.02	0.03	0.02	0.01
	Wetland	0.13	0.08	0	0.03	0.05	0.03	0.04	0.09	0.01
	Water bodies	0.08	0.05	0	0.02	0.04	0.02	0.01	0.06	0
	Built-up areas	0.15	1	0	0.24	0.24	0.25	0.28	0.19	0.24
	Unused land	0.01	0.01	0.77	0.01	0	0	0	0.01	0.01
2010	Cropland	1	0.2	0	0	0.05	0.08	0.04	0.02	0.02
	Forest	0.56	0.4	0.01	0.21	0.23	0.25	0.25	0.23	0.2
	Grass land	0.07	0.05	0	0.01	0.02	0.02	0.03	0.02	0.01
	Wetland	0.12	0.09	0	0.01	0.01	0	0	0.02	0.02
	Water bodies	0.08	0.06	0	0.02	0.04	0.03	0.01	0.06	0.01
	Built-up areas	0.2	1	0	0.3	0.28	0.3	0.35	0.28	0.29
	Unused land	0.01	0.01	0.69	0.01	0	0	0	0.01	0
2015	Cropland	1	0.21	0	0.01	0.05	0.07	0.03	0.07	0.02
	Forest	0.55	0.41	0.01	0.21	0.23	0.25	0.25	0.23	0.2
	Grass land	0.08	0.06	0	0.02	0.02	0.03	0.03	0.02	0.02
	Wetland	0.12	0.09	0	0.01	0.01	0	0	0.02	0.02
	Water bodies	0.08	0.06	0	0.02	0.04	0.03	0.01	0.06	0.01
	Built-up areas	0.21	1	0	0.31	0.29	0.3	0.36	0.29	0.29
	Unused land	0.01	0.01	0.72	0.01	0	0	0	0	0

.

4. Results and Analyses

4.1. Characteristic of Spatiotemporal Pattern in Habitat Quality

4.1.1. Spatial Pattern of Habitat Quality

Figure 6 shows the spatial pattern of habitat quality in the GBA with a range of habitat quality values between 0 and 1 (higher values indicating a more suitable habitat). The spatial pattern displays that the habitat quality index increased from the center to the fringe areas in the period from 1995 to 2015, with some small differences in central areas between years. The overall level of habitat quality was high over the entire study period, with an average score of 0.61, 0.60, 0.59, 0.60 and 0.60 in 1995, 2000, 2005, 2010 and 2015, respectively. The dominant level of high habitat quality (0.8–1) was greater than or equal to 55% and was mainly distributed in the outer ring of the GBA, especially in Zhaoqing, Jiangmen and Huizhou. These areas were dominated by forest in mountain areas with high vegetation fraction cover. The land use intensity was relatively low and the vegetation cover was dense [51], suggesting that the biodiversity value in fringe areas was significantly higher than in central areas. The next most highly distributed levels of habitat quality were the mid-low level (0.2–0.4) and the low level (0–0.2), which accounted for 25% and 10%, respectively, and were concentrated in the central areas of the GBA such as Guangzhou, Shenzhen and Foshan. These areas were dominated by serried industry parks and built-up land with a high intensity of human activities, and thus do not provide much suitable habitat for wildlife. The medium habitat quality level was less than or equal to 5% with a dispersed distribution in Hong Kong, Guangzhou and Dongguan, while the mid-high level had a low distribution of less than or equal to 5% of the total area.

4.1.2. Spatial Pattern of Habitat Quality from Different Perspectives

To evaluate the characteristics of spatial patterns from different ecosystem perspectives, we used the terrestrial ecosystem classification data to extract habitat quality values over the period from 1995 to 2015. As shown in Figure 7, different ecosystem classifications showed a discrepancy in the distribution of habitat quality scores. Farmland and settlements had a consistent distribution with a large value range over the period from 1995 to 2015, which indicates that habitat quality differed within the same ecosystem classification, which could be attributed to differences in climate and elevation. The forest ecosystem had the highest value of habitat quality and was overall close to 1 during the period from 1995 to 2015, which indicates that the forest ecosystem provides a suitable habitat and has high biodiversity value. The grassland ecosystem showed a stable distribution of habitat quality value. Other ecosystems had consistent distributions during the studied period, except for 1995. The desert ecosystem was eliminated from this study because the sample size of this ecosystem was too small to analyze.

To understand the characteristics of spatial patterns from the elevational perspective, a DEM was used to extract habitat quality values at different elevational gradients over the period 1995–2015 (Figure 8). Overall, the habitat quality has an increasing trend with elevation. There was a cluster of habitat quality values in the range of 0.25 to 0.4, which does not change with elevation in the region. It was indicated that this range of habitat quality (0.25 to 0.4) was majorly distributed in plain areas. However, there is a significant trend in the range of 0.5 to 1 that showed that habitat quality increased with altitude, indicating that high habitat quality may be related to elevation. A study in Nanshi Lake basin verified that elevation was significantly correlated with habitat quality [30]. However, this finding differs from a study in the Omo-Gibe Basin, as a different pattern of altitude to habitat quality showed a decreasing trend in the Winike watershed [49].

4.2. Characteristics of Spatiotemporal Variation in Habitat Quality

Spatiotemporal Variation of Habitat Quality

The spatiotemporal variation in habitat quality is shown in Figure 9. In general, the variation in habitat quality in the GBA revealed a conspicuous change in the study years that tended toward habitat degradation overall, except for the period from 1995 to 2000. In the interval value of habitat quality level (low (L), moderate (M) and high (H)) perspective, the high level of habitat quality (e.g., (+H) and (−H)) had a dramatic change with a large distribution in the entire studied period. The moderate level had a small distribution in Guangzhou for the period of 1995 to 2010, though a comparatively larger distribution than that for Foshan, Zhongshan and Dongguan during the period from 2010 to 2015. From the variation perspective, the changing levels of habitat quality showed significant spatial differentiation. The areas with no change (0) were mainly distributed in the cities of central GBA, such as Dongguan, Foshan, Guangzhou and Shenzhen, and had a gradual expansion throughout the course of the studied period. The changing areas in the period from 1995 to 2000 showed that habitat quality improved, with overall positive variation and little habitat degradation in Shenzhen, Foshan and Hong Kong. The variation in the period from 2000 to 2005 showed both positive and negative patterns. However, there was a general decrease both in the period from 2005 to 2010 and from 2010 to 2015.

4.3. Analysis of the Variation of Habitat Quality and Effect Factor

In the period from 1995 to 2000 the habitat quality improved; however, habitat degradation was the general pattern for the GBA during the period from 2000 to 2015. This was associated with multiple factors, such as human activities, land use change and vegetation loss. The effect of land use change was reflected by built-up area expansion and forest land fragmentation. The built-up area increased rapidly from 8.77% to 14.18%, and the area of forest land decreased from 56.42% to 53.94% in the GBA during the period from 1995 to 2015. However, in the period from 1995 to 2000, the percentage of built-up land slightly declined, as shown in Figure 10.

Human activities and vegetation loss were important factors affecting habitat quality in the GBA and, based on the Pearson’s correlation coefficient relating habitat quality to each effect factor in different years, they showed a significant correlation (

Table 5. Pearson’s correlation coefficients of habitat quality in relation to POP, PIS, GDP and FVC.

Year	Effect Factors
	GDP	POP	PIS	FVC
1995	−0.27 **	−0.24 **	−0.24 **	0.47 **
2000	−0.29 **	−0.29 **	−0.23 **	0.51 **
2005	−0.38 **	−0.17 **	−0.28 **	0.48 **
2010	−0.44 **	−0.22 **	−0.25 **	0.5
2015	−0.31 **	−0.35 **	−0.25 **	0.48 **
Mean	−0.34	−0.26	−0.25	0.49

Note: ** Correlation is significant at the 0.01 level (two-tailed); GDP—gross domestic product; POP—population distribution; PIS—proportion of impervious surface; FVC—fractional vegetation cover.

; except for FVC in 2010). FVC, which represented natural vegetation cover, had a strong positive linear relationship with habitat quality. This result indicated that vegetation loss led to habitat degradation. Vegetation loss in urban agglomerations especially decreased the capacity of the urban ecosystem to provide obligatory ecosystem services [49]. POP, PIS and GDP were used to embody human activities, and they had a significant negative linear correlation with habitat quality, which revealed that natural habitat degradation was a result of the growth in population, economic development and the expansion of impervious surface cover. Human activities were verified here as the largest threats to the ecological environment and biodiversity in the GBA. Our findings are similar to those of a previous study of ecological quality in the GBA, which reported that urban expansion and policy variation had a dominant influence on habitat change [52], and another study on the habitat quality of the Pearl River delta found that habitat quality was significantly affected by population and economic density [51].

5. Discussion

5.1. The Role of Sensitivity of Habitat Type to Threat Factors in the InVEST Habitat Quality Model

The sensitivity of habitat types to threat factors is a parameter that uses a value from 0 to 1 to indicate the degree to which the habitat type is affected by each threat factor [25]. A sensitivity value close to 1 signifies that the habitat type is heavily degraded by the threat factor. As mentioned in Section 3.3, the original sensitivity scores were set by an empirical value and relied on an expert knowledge system [25]. However, this parameter is difficult to acquire in some study areas, and gives the HQ model a subjective tendency. Therefore, finding a way to determine sensitivity of habitat types to threat factors using the InVEST-HQ model was considered important in this study. Kendall’s correlation coefficient can describe the strength of association between two categorical variables [53] and has been successfully used in other fields [54,55]. For the InVEST-HQ model we used, the sensitivity of each habitat type to threat factors was estimated by the Kendall’s correlation coefficient (LULC datasets and threat factors were both categorical data), which could calculate an objective sensitivity value specific to the GBA. A limitation of this study is that the results of this strategy to calculate sensitivity were not compared with other strategies, such as the expert-based strategy.

5.2. Exploring the Use of Assessment of Habitat Quality for Biodiversity Conservation

The InVEST-HQ model evaluates habitat quality as an ecological indicator at the landscape scale which can be a proxy for regional biodiversity [27,28]. The model depends on the accessibility and intensity of artificial land use/cover types [56,57]. Recently, the InVEST-HQ model was used to assess the status of biodiversity and habitats worldwide, which could reveal spatial and temporal changes of biodiversity at the largest scale [58].

In this study, the results of the InVEST-HQ model will help to understand the spatial distribution of habitat quality with a consideration for general biodiversity. We obtained a series of habitat quality maps for different years. There was a characteristic spatial pattern which showed increased habitat quality from the center to fringe areas and showed a significant difference among ecosystem classifications and an increase of habitat quality with elevation. This finding can be applied to understand the biodiversity status of the GBA, which could help evaluate the strategies for biodiversity conservation in the region and assist in prioritizing habitat conservation units. A similar mechanism of spatial distribution of habitat quality was observed in other areas, e.g., Nansihu lake basin, China [30], Changli county, China [59], and Italy [4].

The spatiotemporal variation of habitat quality was affected by multiple factors in this study, such as human activities, land use change and vegetation loss. Vegetation loss and intense human activity were the most important factors of habitat degradation, with built-up area expansion and forest land fragmentation as the most clear expressions of land use change. It is necessary in studies of the urban ecosystem to explore the tradeoff areas between nature conservation and economic development, especially in the GBA. Many studies have focused on habitat quality and its influencing factors, such as a study of the habitat quality in Yangtze River Delta Regional, China, which found a significant negative relationship with urbanization intensity [60], and a study of protected areas in Italy found that increased human population density led to a decreasing trend in habitat quality [4].

5.3. Limitations and Future Outlook

The InVEST-HQ model is flexible in terms of parameter settings, which facilitated the integration of a different concept for habitat threat data and its attributes. In particular, the threat factor could be divided into three types: biotic factors, abiotic factors and anthropogenic disturbances [30]. We focused on anthropogenic threats in the InVEST-HQ model in this study. To express the impact of anthropogenic disturbance on habitat quality, we selected cropland, built-up area, bare land, the road network, and residential and industrial areas. However, it is worth noting that we assumed all native species and habitats in the study area were impacted by these threat factors, and we ignored possible alternative consequences. This assumption in relation to threat factors was similar with previous studies, such as a study in Nansihu lake basin, China [30], a study in the state of Minnesota, USA [61], a study in Chaharmahal and Bakhtiari province, Iran [47], and a study in the areas of Hohhot-Baotou-Ordos-Yulin, China [32]. We recognize that this approach could be improved by selecting more threat factors in future studies.

We used Kendall’s correlation analysis to acquire a series of objective sensitivity scores for habitat types to threat factors. Although this study attempted to improve the objectivity by the correlation analysis method, there were some limitations for the other parameters. The weight of the threat factors, the maximum distance of threat factors, and habitat suitability scores were determined by empirical values and expert knowledge, which resulted in some inevitable subjectivity. Previous studies have also optimized some parameters; for example, a study obtained habitat suitability scores by integrating a soil degradation index with the aboveground habitat condition [62], another study obtained habitat suitability scores based on a species distribution model and expert opinion [63], future land use cover types were obtained in a study by cellular automata in different future scenarios, which could predict the future variation in habitat quality [42]. In a future study, we should attempt to find more approaches to add objectivity to the InVEST-HQ model.

Our results of habitat quality in different years were not compared with field site data. The unverified calibration of the field investigation may lead to some uncertainty in the results of the InVEST-HQ model. To overcome this limitation, future studies should obtain a comparison of the performance of our model with other ecological process models, or carry out a field investigation of habitat to calibrate the InVEST-HQ model.

6. Conclusions

In this study we used the Kendall’s correlation coefficient of land use/cover types and threat factors instead of the sensitivity of habitat types to threat factor parameter provided in the InVEST-HQ model. Then, we estimated the spatiotemporal pattern and variation in habitat quality in the GBA during the period from 1995 to 2015 based on the optimized InVEST-HQ model. We analyzed the differentiation of habitat quality across different ecosystem classifications and elevation gradients. Finally, we analyzed the relationship between habitat quality, vegetation cover and anthropogenic threats using Pearson’s correlation coefficients.

Our results indicated that the spatial pattern of habitat quality decreased from fringe areas to central areas, showing significant differences for different ecosystem classifications and elevational gradients. In the period from 1995 to 2015, habitat quality declined as a result of vegetation loss, land use change and intense human activity, and the fraction of vegetation cover was an important effect factor of habitat quality, with a strong positive linear relationship. Human activities were the largest threat to the ecological environment in the GBA, with a significant negative linear relationship. The expansion of built-up areas and forest land fragmentation were clear expressions of land use change. Our method provides an operational approach to evaluating the availability of management strategies for biodiversity conservation and a unique perspective of sustainable development in urban agglomerations that could be used to explore tradeoffs between natural habitats and economic development. The unverified calibration of the field investigation might lead to some uncertainty in the results. Future studies should focus on a comparison of the performance of other ecological process models or design an investigation of habitat to overcome this limitation.

We attempted to develop an optimization method to improve the objectivity of the InVEST-HQ model, and to calculate the sensitivity of habitat types to threat factors by Kendall’s correlation coefficient, which provided an objective and applicable sensitivity value for the GBA. However, a comparison was not made with the original strategy, so future studies should address this limitation.