Spatio-Temporal Heterogeneity of Ecological Quality in Hangzhou Greater Bay Area (HGBA) of China and Response to Land Use and Cover Change

Abstract

Human activities have been stressing the ecological environment since we stepped into the Anthropocene Age. It is urgent to formulate a sustainable plan for balancing socioeconomic development and ecological conservation based on a thorough understanding of ecological environment changes. The ecological environment can be evaluated when multiple remote sensing indices are integrated, such as the use of the recently prevalent Remote Sensing-based Ecological Index (RSEI). Currently, most of the RSEI-related studies have focused on the ecological quality evolution in small areas. Less attention was paid to the spatio-temporal heterogeneity of ecological quality in large-scale urban agglomerations and the potential links with Land Use and Cover Change (LUCC). In this study, we monitored the dynamics of the ecological quality in the Hangzhou Greater Bay Area (HGBA) during 1995–2020, using the RSEI as an indicator. During the construction of the RSEI, a percentile de-noising normalization method was proposed to overcome the problem of widespread noises from large-scale regions and make the RSEI-based ecological quality assessment for multiple periods comparable. Combined with the land use data, the quantitative relationship between the ecological quality and the LUCC was revealed. The results demonstrated that: (1) The ecological quality of the HGBA degraded after first improving but was still good (averaged RSEI of 0.638). It was divergent for the prefecture-level cities of the HGBA, presenting degraded, improved, and fluctuant trends for the cities from north to south. (2) For ecological quality, the improved regions have larger area (57.5% vs. 42.5%) but less increment (0.141 vs. −0.195) than the degraded regions. Mountains, downtowns, and coastal wetlands were the hot spots for the improvement and urbanization, and reclamation processes were responsible for the degradation. (3) The ecological quality was improved for forests and urban areas (△RSEI > 0.07) but degraded for farmland (∆RSEI = −0.03). As a result, the ecological cost was reduced among human-dominant environments (e.g., farmland, urban areas) while enlarged for the conversion from nature-(e.g., forests) to human-dominant environments.

1. Introduction

The ecological environment impacts all levels of human life and the sustainability of society [1]. Human activities have been impacting the ecological environment since humankind began populating the planet, and this process is speeding up. Over the past half century, we have confronted numerous ecological environmental problems, such as vegetation diminishment, water pollution, soil fertility decline, biodiversity loss, and urban heat islands, under an unprecedented process of urbanization [2,3,4]. Therefore, a sustainable plan balancing socioeconomic development and ecological conservation based on a thorough understanding of the ecological environment is urgently needed [5,6,7]. For China, the construction of an ecological civilization is a core part of the programs of President Xi, thought of as socialism for a new era. The State Council consistently emphasized the importance of ecological conservation in the 13th and 14th national Five-Year Plans.

Ecological quality evaluation is a fundamental task that must be undertaken before effective measures for ecological conservation can be implemented. Currently, remote sensing technology is widely employed in ecological quality evaluation, as its advantages of synoptic coverage and repeatability enable it to collect the spatial and temporal data on a large scale [8]. Initial studies for ecological quality evaluation mainly used a single remote-sensing-based index to reflect one aspect of the ecological status. For example, Jose and William [9] used a vegetation index to measure vegetation coverage, Marie-Catherine et al. [10] employed a water index to estimate water abundance, and Zeng et al. [11] used surface temperature to evaluate the effect of urban heat islands. In fact, the ecological environment is affected by many factors, and a single remote-sensing-based index can hardly comprehensively represent all aspects of the ecological environment.

Under this condition, some methods integrating multiple remote-sensing-based indices have been proposed. A typical example is the Ecological Index (EI), which integrates five ecological factors of biological abundance, vegetation cover, water network density, land degradation, and environment pollution, as a technical criterion of ecological assessment for China [12]. However, the EI only provides a simple value to measure the ecological condition for a region; thus, it cannot reflect regional heterogeneity of ecological quality. Similar indicators that integrate multiple remote-sensing-based indices include the Forest Disturbance Index (FDI) [13], Ecological Niche Modeling (ENM) [14], and the MODIS Global Disturbance Index (MGDI) [15]. Although these indicators allow spatializing the results, the challenges still lie in weight assignment. The weights are assigned differently when subjective experience is involved, leading to different evaluation results.

Recently, Xu [16] proposed a Remote-Sensing-Based Ecological Index (RSEI), which integrates the indicators of Greenness, Wetness, Dryness, and Heat to most intuitively reflect the ecological conditions according to the Pressure–State–Response (PSR) framework [17]. More importantly, the Principal Component Analysis (PCA) method was used to determine the weights automatically based on the contribution of each involved indicator. This is helpful to avoid the uncertainty from empirical assignments [18,19]. The RSEI, which overcomes the shortcomings of a single index and makes the aggregation of indices more reasonable, has been accepted across a wide range of temperature zones and climatic conditions [20,21,22,23,24].

Until now, the RSEI has been used in many regions worldwide. In China, it was applied in Fuzhou [18], Xiamen [25], Wuhan [26] cities, and Xiong’an new area [27]. In one of the typical studies, Xu et. al. [27] proposed the RSEI and predicted the ecological effects from increasing population, paving the way for forthcoming regional planning of the Xiong’an new area. In the USA, the land surface ecological status from Minneapolis, Dallas, Phoenix, Seattle, and Chicago have been compared using the RSEI [28]. Moreover, urban surface ecological crisis zones can also be recognized by the RSEI, which has been applied in Madrid, Spain; Rome, Italy; Lyon, France; Hamburg, Germany; Ciechanow, Poland; and Budapest, Hungary [29]. Most of the RSEI-related studies focus on the ecological assessment in small areas of cities and districts. However, less attention has been paid to the spatio-temporal heterogeneity of ecological quality for large-scale areas. For small areas, the maximum–minimum normalization was used for normalization before the RSEI was available. However, the performance of the maximum–minimum normalization may not be satisfactory if widespread noises for large-scale regions are not considered. Large-scale regions, such as greater bay areas, have often experienced rapid growth and are looking for ways to implement high-quality development that coordinates economic development with environmental safeguards [30]. In addition, most RSEI-related studies reveal the evolutionary process of regional ecological quality, while the quantitative analysis of potential reasons was less explored [31,32]. In contrast, evaluating land use and cover change (LUCC) is the most intuitive way to measure how humans alter natural environments, quantifying its impacts on regional ecological quality. Such studies benefit future urban development.

This study used Hangzhou Greater Bay Area (HGBA), an important urban agglomeration in the Yangtze River Delta, China, as the study area, (i) to reveal the spatio-temporal heterogeneity during the ecological quality evolution over the past 25 years using an RSEI with a percentile de-noising normalization; and (ii) to explore the quantitative relationship between regional ecological quality and land use types and conversions. The percentile de-noising normalization method we propose can reduce the impact of widespread noise and pave the way for large-scale application of an RSEI-based ecological assessment, and the findings are expected to provide useful guidance for high-quality socioeconomic development for the HGBA, thus being a value to local authorities and policy makers.

2. Materials

2.1. Study Area

Hangzhou Greater Bay Area (HGBA, 118°21′–123°25′E, 28°51′ –31°53′N) is one of the three major bay areas in China. The HGBA is located on the east coast of China and consists of 6 prefecture-level cities—Jiaxing, Hangzhou, Shaoxing, Ningbo, Zhoushan in Zhejiang Province and Shanghai Municipality (Figure 1). The HGBA has a subtropical monsoon climate, with annual precipitation of 1709.0 mm and an average temperature of 17.8 °C [33,34].

Profiting from a series of growth policies (e.g., the Belt and Road Initiative, the Yangtze delta economic plan), socioeconomic development in HGBA has boomed in recent years. During 1990–2019, the GDP for the HGBA increased from CNY 131 thousand billion to CNY 7786 thousand billion, and the proportion of tertiary industry gradually increased and surpassed the secondary industry [35,36]. However, the environmental problems (e.g., deforestation, air and water pollution, biodiversity loss) for the HGBA have become increasingly prominent [37], which, in turn, restricts sustainable socioeconomic development.

2.2. Datasets

In our study, two types of data were used—remote sensing data and auxiliary data. The remote sensing data was Landsat C2L2 products, which were mainly used for generating the RSEI and classifying land use types. The auxiliary data was GlobeLand30 data, which was used for interpreting the class information of training and validation samples to assist with land use classification.

2.2.1. Remote Sensing Data

We used the Landsat Collection 2 Level 2 (C2L2) product, which was downloaded free of charge from the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/ (accessed on 1 August 2022)). This product provides the surface reflectance for each optical band after the pre-procession of geometric correction, radiometric calibration, and atmospheric correction [38]. Unlike the product from previous versions, this product provides a surface temperature band, which is deduced by the combinations of thermal infrared bands, brightness temperature, atmospheric profiles, and emissivity data with the assistance of multi-source remote sensing data [39].

A total of 20 Landsat C2L2 products with low cloud cover percent were collected (

Table 1. Details for Landsat C2L2 products used in our study.

Satellite/ Sensor	Path/ Row	Acquisition Date	Cloud Cover	Satellite/ Sensor	Path/ Row	Acquisition Date	Cloud Cover
Landsat-5 TM	118/38	1995-05-08	0.0%	Landsat-5 TM	118/40	2007-07-28	1.0%
Landsat-5 TM	118/39	1995-08-12	0.0%	Landsat-5 TM	119/39	2008-05-02 *	1.0%
Landsat-5 TM	118/40	1995-05-08	4.0%	Landsat-8 OLI	118/38	2015-08-03	0.0%
Landsat-5 TM	119/39	1994-05-12 *	0.0%	Landsat-8 OLI	118/39	2015-08-03	0.3%
Landsat-7 ETM+	118/38	2000-09-18	0.0%	Landsat-8 OLI	118/40	2015-08-03	0.7%
Landsat-7 ETM+	118/39	2000-09-18	0.0%	Landsat-8 OLI	119/39	2015-05-22	9.3%
Landsat-7 ETM+	118/40	2000-09-18	0.0%	Landsat-8 OLI	118/38	2020-08-16	0.0%
Landsat-5 TM	119/39	2000-09-17	0.0%	Landsat-8 OLI	118/39	2020-08-16	0.3%
Landsat-5 TM	118/38	2007-07-28	0.0%	Landsat-8 OLI	118/40	2020-08-16	12.2% **
Landsat-5 TM	118/39	2007-07-28	0.0%	Landsat-8 OLI	119/39	2020-05-19	4.2%

* Because of the large area in our study, to ensure the image quality and consistency in vegetation phenology, we selected the Landsat products from neighboring years. ** Although the cloud cover seems high, most of the cloud-covered parts are not included in our study.

). These products come from five periods (1995, 2000, 2007, 2015, and 2020) and four swaths (118/038, 118/039, 118/040, and 119/039). To avoid the influence from vegetation phenology, only the products acquired during late spring to early autumn (e.g., May to September) were selected. However, because of frequent cloudy weather, it is not an easy task to collect enough high-quality products for a specific period in the large coastal areas. Therefore, for 1995 and 2007, we selected a product individually from the neighboring year for ensuring the image quality and consistency in vegetation phenology.

2.2.2. Auxiliary Data

The Globeland30 data (http://globeland30.org/ (accessed on 1 August 2022)) provides global land use and cover data with a spatial resolution of 30 m [40]. This data includes 10 land use and cover types, namely, farmland, forest, grassland, bushland, wetland, water, tundra, build-up, bare land, glacier, and permanent snow. The overall accuracy of this data was reported as 85.8% with the Kappa coefficient of 0.82. In our study, the Globeland30 data from 2000, 2010, and 2020 were used for land use and cover classification and accuracy assessment. By overlapping the Globeland30 data and the HGBA area, we found that there was almost no bare land, tundra, glacier, or permanent snow. The grassland and bushland were also rare and commonly appeared surrounding forest areas. Therefore, after integrating grassland and bushland into forest, we determined the final land use types for the HGBA as farmland, forest, wetland, water, and build-up.

The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) was downloaded free of charge from the USGS website (https://earthexplorer.usgs.gov/ (accessed on 1 August 2022)). These data were created from the entire ASTER level 1A archive of scenes acquired between 1 March 2000 and 30 November 2013, with a spatial resolution of 1 arc second. The agreement on both temporal and spatial resolution with the Globeland30 make the ASTER GDEM suitable to delineate the topography of the HGBA after the spatial resampling to 30 m.

3. Methods

3.1. RSEI Construction

The assumption of the RSEI is the PSR framework [41,42], which links pressures on the environment caused by human activities and environmental state (condition) changes. Based on the PSR framework, the indicators used for constructing the RSEI were selected as follows. (1) Pressure is a general awareness that the LUCC has significant effects on the ecological patterns within their boundaries and beyond [43]. The most notable pressure is the LUCC from vacant land to built-up land. Thus, an indicator of Dryness representing the intensity of building and soil bareness needs to be employed. (2) State is the Greenness, which is one of the most successful of many attempts to simply and quickly identify vegetated areas to describe the status quo of the environment [44,45]. (3) Response incorporates Wetness and Heat, usually applied to indicate local climate changes in response to environmental changes and concerns (i.e., humidity and temperature) [46,47].

Therefore, we first constructed remote-sensing-based indices to represent the above four indicators (Greenness, Wetness, Dryness, and Heat). Then, a percentile de-noising and normalization method was used to unify the dimension of each indicator after the removal of outliers. Afterwards, the four indicators were used for the PCA transformation and for obtaining the RSEI. The value of the RSEI ranges from 0 to 1, where 0 represents poor ecological conditions and 1 represents good ecological conditions [16]. For common use, it can be divided into 5 grades linearly with an interval of 0.2, namely Poor (0–0.2), Fair (0.2–0.4), Moderate (0.4–0.6), Good (0.6–0.8), and Excellent (0.8–1.0), to further qualitatively describe the ecological condition [27].

3.1.1. Index Calculation

The formula and description for each index are presented in

Table 2. Formula and description for each remote-sensing index.

Index		Formula	Description
NDVI	$\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}$		${\rho }_{NIR}$ and ${\rho }_{Red}$ respectively denote the surface reflectance from NIR band and red band in the Landsat product [48].
WET	Landsat TM:	$0.0315{\rho }_{Blue}+0.2021{\rho }_{Greem}+0.3102{\rho }_{Red}$ $+0.1594{\rho }_{NIR}-0.6806{\rho }_{SWIR1}-0.6109{\rho }_{SWIR2}$	${\rho }_{Blue}, {\rho }_{Greem}, {\rho }_{Red}, {\rho }_{NIR},{\rho }_{SWIR1},$and ${\rho }_{SWIR2}$respectively denote the surface reflectance from blue band, green band, red band, NIR band, SWIR1 band and SWIR2 band in the Landsat product [49,50,51].
	Landsat ETM+:	$0.2626{\rho }_{Blue}+0.2141{\rho }_{Greem}+0.0926{\rho }_{Red}$ $+0.0656{\rho }_{NIR}-0.7629{\rho }_{SWIR1}-0.5388{\rho }_{SWIR2}$
	Landsat OLI:	$0.1511{\rho }_{Blue}+0.1973{\rho }_{Greem}+0.3283{\rho }_{Red}$ $+0.3407{\rho }_{NIR}-0.7117{\rho }_{SWIR1}-0.4559{\rho }_{SWIR2}$
NDBSI	$NDBSI=\frac{SI+IBI}{2}$ $SI=\frac{\left({\rho }_{SWIR1}+{\rho }_{Red} \right)-\left({\rho }_{NIR}+{\rho }_{Blue}\right)}{\left({\rho }_{SWIR1}+{\rho }_{Red} \right)+\left({\rho }_{NIR}+{\rho }_{Blue}\right)}$ $IBI=\frac{\frac{2{\rho }_{SWIR1}}{{\rho }_{SWIR1}+{\rho }_{NIR} }-\left(\frac{ {\rho }_{NIR}}{ {\rho }_{NIR}+{\rho }_{Red}}+\frac{{\rho }_{Green}}{{\rho }_{Green}+{\rho }_{SWIR1}}\right)}{\frac{2{\rho }_{SWIR1}}{ {\rho }_{SWIR1}+{\rho }_{NIR} }+\left(\frac{ {\rho }_{NIR}}{ {\rho }_{NIR}+{\rho }_{Red}}+\frac{{\rho }_{Green}}{{\rho }_{Green}+{\rho }_{SWIR1}}\right)}$		${\rho }_{Blue}, {\rho }_{Green}, {\rho }_{Red}, {\rho }_{NIR}, {\rho }_{SWIR1}$ respectively denote the surface reflectance from the blue band, green band, red band, NIR band, and SWIR1 band in the Landsat product [52].
LST	$\lambda \rho +\alpha -273.15$		A linear conversion is used to stretch the data to real surface temperature [38]. $\rho$ is the temporal band in the Landsat product; $\lambda$and $a$ are 0.00341802 and 149, respectively.

, specifically:

The Normalized Difference Vegetation Index (NDVI), as one of the common indices to reflect vegetation status, can accurately describe the coverage and growth of the surface vegetation. Therefore, NDVI was used to represent Greenness in our study. Wetness is the ecological indicator that measures the surface water condition, or moisture in soil and vegetation. The WET component from the Tasseled Cap Transform (TCT) based on original Landsat images was used to represent Wetness, since it has been proven to have strong implications for land surface moisture. Artificial impervious surfaces and natural bare soils are the two sources for surface dryness. Therefore, Dryness was measured by the Normalized Difference Build-up and Soil Index (NDBSI), a combination of the Index-Based Built-up Index (IBI) and the Soil Index (SI). Land Surface Temperature (LST) is a useful index for reflecting Heat given it is closely related to the growth and distribution of vegetation and the evaporation of surface water. In our study, the surface temperature band from the Landsat C2L2 product was directly used to represent LST after a simple conversion according to the Landsat user manual.

Given the RSEI is mainly used in terrestrial environments but less suitable for large-scale water areas, the Modification of Normalized Difference Water Index (MNDWI) [53] was employed to mask out water surfaces (e.g., rivers, lakes, seawater). All the indices were calculated using the Band Math tool of the ENVI 5.3 software.

3.1.2. Percentile De-Noising and Normalization

Given the ranges for the original four indices were commonly variant, the weights would be biased if the PCA transformation was directly implemented. Simultaneously, for large areas, such as the HGBA, the noises incurred by poor image quality, tiny clouds/shadows, or omitted water bodies, occurred in the original indices with high possibility. Therefore, before the implementation of the PCA transformation, noise removal and index normalization should be considered.

Commonly, the maximum and minimum values are used for normalization (i.e., maximum–minimum normalization). However, the noises with extremely high or low values limited the validated data in a very narrow range with the most remaining range unused (Figure 2(a1–b2)), which heavily affected the performance of the maximum–minimum normalization. Under this condition, we proposed a percentile de-noising normalization method using the cutoff of α. For a specific index, the values from all the pixels falling in the HGBA were ranked ascendingly, and the values exceeding the range between α% (the lower boundary) and (100-α)% (the upper boundary) percentiles were regarded as noise. The noise below the lower boundary was assigned to the lower boundary value, and the noise over the upper boundary was assigned to the upper boundary. The α was set as 0.5 for almost all periods except 2015 when a certain number of tiny clouds existed. The cooler temperatures in the cloud-covered areas lead to more noise in the LST (Figure 2(c1–d2)). Therefore, we slightly adjusted the α to 1.5 in 2015 for the consistency of the percentile values from adjacent periods. For example, if the α was set 0.5, the lower boundary in 2015 would be 3.9 °C, significantly lower than 24.4 °C in 2007 and 23.0 °C in 2020. When the α was set tp 1.5, the lower boundary in 2015 increased 20.6 °C, which was comparable with that in 2007 and 2020.

After modifying the noises, we used the lower and upper boundary values (also the maximum and minimum values) from all the periods (1995, 2000, 2007, 2015, and 2020) for index normalization (Equation (1)). The data range for each index after the normalization was uniform between [0, 1], paving the way for the PCA transformation. More importantly, the land use with different ecological conditions was clearly discriminated (Figure 2(a3,d3)), in sharp contrast with that from maximum–minimum normalization.

$$NI\left(i\right)=\frac{I\left(i\right)-ILB\left(i\right)}{IUB\left(i\right)-ILB\left(i\right)}$$

(1)

where $NI\left(i\right)$ is the normalized index $i$ ($NDVI, Wet, NDBSI, LST$) after noise removal; $I\left(i\right)$ is the original index$i$; and $I_{LB}\left(i\right)$ and $I_{UB}\left(i\right)$ are the lower and upper boundary for the original index$i$ and $I_{UB}\left(i\right)$ which were determined by all the periods. The whole process of the de-noising and normalization was implemented using the Numpy package in Python 3.6.

3.1.3. PCA Transformation and Analysis

The four indices after the normalization were used for the PCA transformation. The PCA is a multi-dimensional data compression technique that can remove the correlation among the original indices [54]. Moreover, the weights were calculated automatically based on the contributions from each index, avoiding the uncertainty from manual assignment [18]. The first principal component was deemed as the initial RSEI (RSEI₀), which can be calculated using Equation (2).

$${RSEI}_0=PCA1\left(NI\left(NDVI\right),NI\left(WET\right),NI\left(NDBSI\right),NI\left(LST\right)\right)$$

(2)

where PCA1 represents the first principal component from the PCA transformation, and NI(NDVI), NI(WET), NI(NDBSI), and NI(LST) denote the remote sensing indices after the normalization. The percentage de-noising (0.5%) was also implemented on the RSEI₀ to avoid the influence from potential noises. To facilitate the comparison and grading, the RSEI₀ was further normalized to [0, 1] and generated the final one (RSEI) using Equation (3).

$$RSEI=\frac{RSE{I}_0-RSE{I}_{0 min}}{RSE{I}_{0 max}-RSE{I}_{0 min}}$$

(3)

where RSEI_{0 min} and RSEI_{0 max} denote the minimum and maximum value calculated from all the periods from the RSEI₀ after de-noising, for ensuring the temporal comparability. We conducted the PCA transformation in the ENVI 5.3 software using the tool named Forward PCA rotation New Statics and Rotate.

As shown in Figure 3, the loadings were uniformly greater than 0 for the indicators of Greenness and Wetness, indicating that the Greenness and Wetness played a positive role in the RSEI generation, and among the four indicators, the loading for the Greenness always ranked first for each period, highlighting the utmost importance of vegetation in an ecological environment. These results are in high agreement with our common sense. The percentage eigenvalue of the PCA1 to the RSEI exceeded 80% for 2000 and 2020 and was never less than 70% for the rest of periods. That means the PCA1 has integrated the most information from the four indices, and its derived RSEI is likely reliable for comprehensively representing the ecological environment.

3.2. Land Use Classification and Assessment

We used the tool of Create Random Points of the ArcGIS 10.3 software for generating training and validating samples for the land use classification (Figure 4). The number of random points varied from 3000 to 5000 for each scene of the Landsat products based on the coverage in the HGBA and the complexity for land use types. The Globeland30 data was used for determining the land use type. Specifically, when the study period coincided with the acquisition period of the Globeland30 data (e.g., 2000, 2020), the land use type for the random point was directly interpretated from the Globeland30 data; otherwise, the land use type was first interpreted from the Globeland30 data with the closest time to the study period, then modified if the land use change was further detected by interpreting the corresponding Landsat product. In our study, half of the random points (50%) were selected for training the classifier, the performance of which was validated by the rest of the random points (50%).

The Support Vector Machine (SVM) algorithm was selected as the classifier. Based on statistical learning theory, the SVM presents the unique advantages in solving small samples, recognizing non-linear and high-dimensional patterns [55]. This algorithm distinguishes the classes with a decision surface that takes full advantage of the margin between the classes. The surface is termed as optimal hyperplane, and the points of data closest to the hyperplane are named as support vectors. The SVM has been widely applied in classification mapping of land use because of its high accuracy [56,57]. Therefore, the SVM algorithm with the radial basis kernel function was used for land use classification, and this process was completed with the help of ENVI 5.3 software.

The confusion matrix [58] was adopted to evaluate the performance of land use classifications. The overall accuracies of land use classifications in 1995, 2000, 2007, 2015, and 2020 were 82.5%, 83.4%, 85.7%, 83.1%, and 78.2%, with an average of 82.6%, indicating that the land use maps we classified can fairly reflect the real situation of land surface in the HGBA.

Table 3. Accuracy assessments for the land use classifications of the HGBA in 2007 and 2020 (the values in the brackets come from 2007, otherwise the values come from 2020).

	Build-Up	Farmland	Forest	Water	Wetland	User’s Accuracy (%)
Build-up	518 (965)	206 (38)	26 (188)	20 (21)	7 (4)	66.7 (79.4)
Farmland	87 (19)	1133 (1655)	75 (251)	32 (3)	0 (2)	85.4 (85.8)
Forest	13 (249)	78 (210)	1546 (981)	8 (23)	0 (18)	94.0 (66.2)
Water	7 (13)	26 (15)	5 (40)	284 (289)	0 (2)	88.2 (80.5)
Wetland	0 (7)	1 (0)	1 (1)	0 (1)	87 (83)	97.7 (90.2)
Producer’s accuracy (%)	82.9 (77.0)	78.5 (86.3)	93.5 (67.1)	82.6 (85.8)	92.6 (76.1)	66.7 (79.4)
Overall accuracy = 85.7% (78.2%), Kappa coefficient = 0.762 (0.690)

presents the details on accuracy assessment from the best (2007) and worst (2020) classification performance. The classification maps were used for evaluating the variations in ecological conditions among land use and their conversions.

3.3. Exploratory Spatial Data Analysis

Exploratory spatial data analysis is a technique to detect spatial distribution patterns or spatial dependence at different scales [59], mainly including the analyses of global and local spatial autocorrelations.

The global spatial autocorrelation measures the degree and significance of the spatial correlation among variables within the whole area [60], and the Global Moran’s I is a common statistic for this autocorrelation (Equations (4) and (6)). Although facilitating to determine any prominent spatial distribution pattern, the global spatial autocorrelation can hardly tell which objects are similar (or different) to their neighborhood. To overcome this difficulty, the local spatial autocorrelation is proposed to identify specific spatial correlation patterns (High–High, Low–Low, High–Low, and Low–High) as well as the locations [61]. This autocorrelation is usually indicated by Local Moran’s I (Equations (5) and (6)).

$$I=\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j\ne 1}^{n}wij\left(xi-\overline{x}\right)\left(xj-\overline{x}\right)}{S^2\left({{\displaystyle \sum }}_i{{\displaystyle \sum }}_{j}wij\right)}$$

(4)

$$I_i=\frac{\left(xi-\overline{x}\right){{\displaystyle \sum }}_{j=1}^{n}wij\left(xi-\overline{x}\right)}{S^2}$$

(5)

$$S^2=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{(xi-\overline{x})}^2$$

(6)

where $I$ is the Global Moran’s I, $I_i$ is the Local Moran’s I, $n$ is the number of grid units, $x_i\left(x_j\right)$ is the value of the variable at $i\left(j\right)$ grid cells, $w_{ij}$ is the normalized spatial weight matrix, and $S^2$ is the standard variance to measure how much the values at neighborhood locations different from the mean.We selected a K-nearest neighbor model as the spatial weight matrix in the GeoDa software. Both Global Moran’s I and Local Moran’s I were calculated at 53 county-level cities in the HGBA, on the basis of the RSEI change during 1995–2020.

The exploratory spatial data analysis highlighted the hot spots for the changes in ecological quality, where the topography would be further focused. Specifically, the topography was divided into plain and mountainous areas with the GDEM and land use data. The areas with the elevation lower than 200 m and mainly covered by build-up, farmland, and wetland were delimited as the plains, while the areas with the elevation higher than 200 m and predominantly covered by forests were delimited as the mountains. Afterwards, the area of the topography-related hot spots of ecological quality changes were accurately calculated by means of spatial overlay. For example, when referring to the hot spots with ecological quality improvements related in the mountains, we overlapped the mountainous areas with the increased RSEI areas in the High–High clusters.

4. Results

4.1. Spatio-Temporal Heterogeneity in Ecological Quality

The average RSEI in the HGBA during 1995–2020 demonstrated that the overall ecological quality degraded after an initial improvement (Figure 5a). Specifically, during 1995–2000, the ecological quality improved as the average RSEI increased from 0.633 to 0.670. During 2000–2020, the ecological quality degraded, a rapid decline during 2000–2015 in particular, as the average RSEI value decreased from 0.670 to 0.619. During 2015–2020, the ecological quality changed a little with a slight decrease of 0.003 in the average RSEI. In general, the average RSEI in the HGBA was relatively high (0.638) with small fluctuations (<0.1), indicating a good ecological condition. Although the overall ecological quality in the HGBA changed a little numerically, the trend of spatial divergence became increasingly obvious (Figure 6(a1–a5)). This was evidenced by the areas of the poor and excellent grades, which increased 414 km² (30.1%) and 3813 km² (51.7%), respectively (Figure 6b). Moreover, the standard deviation drastically increasing from 0.015 in 1995 to 0.074 in 2020 (Figure 5a) also demonstrated the spatial divergence.

The average RSEI in the prefecture-level city of the HGBA demonstrated that the distinctly different processes for ecological quality changes are mainly related to the locations of the cities (Figure 5b). The trend of ecological quality changes for each prefecture-level city during 1995–2020 can be divided into three categories. (1) Shanghai and Jiaxing are located in the northern part of the HGBA, and their ecological quality degraded significantly. In the past 25 years, the average RSEI decreased 0.087 for Shanghai and 0.140 for Jiaxing (Figure 5b). At the end of 2020, the RSEI for these two cities was lower than 0.6, meaning a moderate ecological condition. Spatially, a gradually radial degradation pattern for the ecological quality from the central urban areas to the neighboring areas can be clearly observed (Figure 6(a1–a5)). (2) Hangzhou is located in the southwest of the HGBA, and the ecological quality of the city improved in spite of some fluctuations. The average RSEI increased from 0.629 to 0.679, demonstrating an increasingly good ecological condition (Figure 5b). The overall rise in the RSEI for the most southwest areas of Hangzhou is the main reason for this improvement (Figure 6(a1–a5)). (3) Shaoxing, Ningbo, and Zhoushan are located in the south of the HGBA. The ecological quality for these cities first improved then degraded, which was similar to that for the HGBA. In the past 25 years, the change of average RSEI was not very large—it increased from 0.623 to 0.670 for Shaoxing, decreased from 0.657 to 0.615 for Ningbo, and slightly decreased from 0.615 to 0.612 for Zhoushan (Figure 5b) and always maintained a good ecological condition. Spatially, there was an opposite trend that the ecological quality radially degraded in the north while uniformly improved in the south for these cites (Figure 6(a1–a5)).

4.2. Hot Spots of Ecological Quality Change

By overlapping the RSEI from 1995 and 2020, we explored the hot spots for the RSEI changes in the HGBA (Figure 7a). In the past 25 years, the area for improved ecological quality reached 20,496.7 km² (57.5%), with an average RSEI increase of 0.141, and the area for degraded ecological quality was 15,164.4 km² (42.5%) with an average RSEI decrease of 0.195. There was an obvious spatial aggregation for the areas with both improved and degraded ecological quality. The areas with degraded ecological quality surrounded the coasts of the HGBA, while the areas with improved ecological quality were mostly distributed in the south or west regions far away from the coasts (Figure 7a).

On the basis of the Global Moran’s I, the scatter points were mainly distributed in the first and third quadrants, indicating positive spatial correlations in the HGBA. Moreover, the Global Moran’s I for each period was greater than 0.654, indicating a high degree of spatial correlation. The Local Moran’s I based on 53 counties demonstrated that most cluster patterns were High–High and Low–Low relations, with few Low–High relations (Figure 7b). Specifically, the High–High relations appeared in parts of Hangzhou (Figure 7b4), central Shanghai (Figure 7b5), and Chongming and Changxing Island (Figure 7b6). Most of the Low–Low relations were widespread in Jiaxing (Figure 7b1,b2), with a few scattered in the north of Shanghai and central Ningbo. For the Low–High relations t were located in the north of Shanghai, intermediated between two High–High relations—central Shanghai and Chongming Island.

The formation of such spatial patterns was likely related to the land use change during 1995–2020. The rapid urbanization process led to the significant degradation in ecological quality for 3397.6 km² of the plain areas surrounding central cities (especially for the suburbs of Jiaxing, Shanghai, and Ningbo), forming the Low–Low clusters. The decreases in the RSEI for the Jiaxing downtown (Figure 7(a1–f1)) and the districts and counties (Tongxiang, Xiuzhou, Jiashan, Jinshan, Pinghu, and Haiyan) (Figure 7(a2–f2)) can reach 0.161, which was evidence of the above process. The natural wetlands disappeared through tidal flat reclamation was also the reason for ecological quality degradation. There were 540 km² natural wetlands replaced by artificial wetlands in the south coasts of Hangzhou Bay, resulting in a decreased RSEI of 0.254 (Figure 7(a3–f3)). In contrast, the most improvements in ecological quality came from mountainous areas of Hangzhou (Yuhang, Gongshu, Xihu, Jianggan Binjiang, and Fuyang), with an area of 2478.9 km². A typical example was the Lin’an District of Hangzhou where the RSEI near-evenly increased 0.119 and formed High–High clusters (Figure 7(a4–f4)). Some central urban areas from prefecture-level cities contributed to the improvement in ecological quality. For Huangpu, Xuhui, and Changning Districts in Shanghai, the gradual configuration of green infrastructures after establishment of built-up areas increased the RSEI of 0.136 (Figure 7(a5–f5)). Moreover, the regeneration of coastal wetlands was another reason for the improvement. Over the last 25 years, the coastal wetland rapidly extended seawards for east coasts of Chongming Island under sufficient sediment from the Yangtze River. As a result, over 290 km² coastal wetlands, which are mainly located on the east coasts of Shanghai and Chongming Island, regrew with an increase of 0.121 in the RSEI (Figure 7(a6–f6)).

4.3. Response of Ecological Quality to Land Use

We quantified the relationship between ecological quality and land use from the perspective of different land use types and their conversions. During the past 25 years, 69.3% of the HGBA has not experienced land use change (Figure 8a). Among them, forest (14,345.2 km²), farmland (9725.7 km²), and built-up (1785.9 km²) were the three predominant land use types, the area of which accounted for more than 90.0% of the total.

Dominated by the natural environment, forests commonly ranked in the Excellent or Good grade of ecological quality (Figure 9a). Except for a small decrease during 2007–2015, the average RSEI for forests increased from 0.72 in 1995 to 0.81 in 2020 with an increment of 14%, indicating a continuous improved ecological quality. Farmland and built-up areas, as two human-dominated ecosystems, presented divergent trends for ecological quality. The ecological quality of farmland mostly included the three grades of Fair, Moderate, and Good (Figure 9b). Despite an initial increase from 0.62 to 0.71 during 1995–2000, the average RSEI substantially decreased to 0.59 in 2020, indicating an overall degradation in ecological quality. On the contrary, the ecological quality for the built-up areas was improved, with the most prominent period during 2000–2007. The average RSEI increased from 0.22 in 1995 to 0.29 in 2020, with a huge increment more than 30% (Figure 9c). Note that the overall ecological quality for the human-dominated ecosystem still degraded due to the overwhelming proportion of farmland. As a result, the gap in ecological quality between nature- and human-dominated environments became enlarged, as reflected by the average RSEI increasing from 0.14 to 0.34.

Simultaneously, the area for land use change in the HGBA in the past 25 years also reached 12,875.9 km², among which, the conversion between farmland and build-up (48.0%) and the conversion between forest and farmland (28.1%) were the two typical modes (Figure 8b). The ecological quality changes for the conversion between cropland and built-up areas fluctuated, especially for the period of 2000–2007 when the range for the RSEI changes reached 0.332 (Figure 9d). In spite of the fluctuations, the ecological quality degradation for the conversion from farmland to built-up areas continued to slow, as the average of the RSEI changes decreased from 0.207 to 0.066. That means the boundary effect between farmland and build-up was weakening, and the ecological cost for the conversion within human-dominant ecosystems became lower, due to the improvement for built-up areas and the degradation for farmland (Figure 9b,c). In contrast, the ecological quality degradation for the conversion from forest to farmland enlarged significantly, with relatively small fluctuations (Figure 9e). The average of the RSEI changes increased from 0.050 to 0.091, with the maximum of 0.131 during the period of 2007–2015. This mainly resulted from the ecological improvement for forest and the degradation for farmland (Figure 9a,b), which pre-warned that the further development of forest resources should be more cautious because the ecological cost for the conversion from forest to a human-dominant environment, farmland in particular, has been becoming higher.

5. Discussion

5.1. Potential Reasons for Ecological Quality Change

As a core part of the Yangtze River Delta, the HGBA experienced rapid urbanization during 1995–2020. Based on our classification maps, the built-up areas increased more than three times from 2555.13 km² to 8346.95 km² in the past 25 years. Under such rapid urbanization, it is interesting to find that the overall ecological quality for the HGBA still maintained a high level (RSEI ≈ 0.62, a good grade), in spite of slight deterioration (Figure 5a).

The ecological quality for forest improved continuously as the average RSEI increased from 0.72 to 0.81 during 25 years, which was an important factor for ecological quality maintenance of the HGBA (Figure 8a). This is likely related to global warming and carbon dioxide (CO₂) concentration. With the development of industrialization and urbanization, regional carbon emission and its derived greenhouse effect gradually increase. Previous studies have reported that the fertilization effect of CO₂ plays an important role in large-scale “greening” of forests in China, especially for southeastern China [62]. This phenomenon emerges owing to the fact that the increase of CO₂ concentration will promote the growth of plant roots, leaves, and stomata and secrete hormones, thus increasing vegetation biomass [63]. Therefore, the ecological quality for the southern and western HGBA where forests were widespread improved globally in the past 25 years (Figure 7(a4–e4)).

Moreover, the ecological quality for the central urban areas of each prefecture-level city also observed improved, Shanghai in particular (the RSEI increased 0.136) (Figure 6(a5–e5) and Figure 8c). With the change of development concepts, urban planning has paid more attention to building a favorable ecological environment for not only new towns but also old towns through renovation [64]. As an important manner, the industrial structure for central urban areas has been optimizing to gradually change the function from traditional industry to modern service [65]. Recently, the ecological conservation and restoration measurement such as urban parks and green spaces was another effective way to improve the ecological quality for central urban areas [66,67,68].

However, the ecological quality for farmland presents an obvious degraded trend, reflected by the RSEI value decreasing from 0.71 in 2000 to 0.59 in 2020 (Figure 7(a2–e2) and Figure 9b), which is related to the non-agriculture and non-food phenomenon during land use change [69]. Our classification maps revealed the sharp shrinkage of farmland during 1995–2020. Urban expansion occupying parts of the high-quality farmland was the main reason for this degradation. Moreover, most of the new reclaimed farmland comes from forests and wetlands, easily leading to problems of land degradation caused by wind erosion and secondary salinization [70]. Moreover, ill usage and poor management for farmland was also the reason for deteriorating ecological quality [71].

Notably, compared to the rapid decrease in the RSEI during 2000–2015, the decrease in the RSEI obviously slowed in the following 5 years with a slight decrease of 0.003 in the average RSEI (Figure 5a). This may be closely related to a series of policies. In 2016, the Chinese government proposed the 13th Five-Year Plan (2016–2020) as a national policy. In this policy, the objectives for ecological conservation and the related countermeasures were established [72]. The local policies were subsequently proposed by Zhejiang Province and Shanghai Municipality such as the plan for institutional reform for promoting ecological progress in 2017, and the delimitation of ecological red line in 2018 [73]. It turns out the trend of ecological degradation can be restrained with the introduction of appropriate regulatory measures.

5.2. Contributions and Limitations of Our Study

The contribution of our study embodies the two aspects. (1) We proposed a percentile de-noising normalization method, which is important for large-scale application of the RSEI-based ecological quality assessment. Previously, the maximum–minimum normalization was used. However, for large-scale regions, widespread noises with extremely high or low values limit the validated data in a very narrow range, heavily affecting the normalization performance. The percentile de-noising normalization method can overcome this problem and make the RSEI-based ecological quality assessment from multiple periods comparable. (2) We quantified the relationship between ecological quality and LUCC, which is likely of value to local authorities and policy makers. Most of the previous studies either revealed spatio-temporal changes in regional ecological quality or predicted the future changes in RSEI value using multiple regression models. Quantitative analysis for the reasons of the changes in ecological quality were rare. The analysis on the LUCC-derived ecological quality change as we have carried out may be the most direct and useful way for serving urban future development in spite of its simplicity.

However, we evaluated the ecological quality for the HGBA with a relatively long interval of ~5-years. The uncertainty from the observations with long time intervals mainly comes from the two aspects. (1) The trend of regional ecological quality may be not very accurate, especially for the key points. For example, the ecological quality of the HGBA was deemed to reach the maximum in 2000 (Figure 5a); however, it is possible to achieve this maximum for another year during 1996–2006 if the shorter time-interval observations are used. (2) The relationship between the ecological quality and land use types was not very favorable because of only five points (periods) (Figure 9a–c). If more periods were involved, the performance for the relationship is expected to improve, and a more detailed relationship can be obtained by sub-dividing it into different periods.

Notably, because of frequent cloudy weather, it is not an easy task to collect cloud-free images for coastal areas. Given the area of the HGBA crosses four Landsat scenes, the collection of the four images with near acquisition time is rather difficult. As a result, the time interval of the adjacent observations and the difference of the acquisition time among images is usually a trade-off. If the time interval become shorter (i.e., the number of time series increases), the difference of the acquisition time may be larger. In our study, although we used the 5-year time interval, the difference of the acquisition time in 1995 still reached nearly 3 months (Table 1), and based on our experiment, some of the indices (e.g., NDVI) used for constructing the RSEI are affected by intra-annual phenology, leading the unignorable difference in the RSEI. For example, the difference between the mean values of the RSEI acquired in May and August is close to 10%, and this difference becomes more prominent for the farmland, which is widespread in Jiaxing city. Therefore, after careful consideration, we selected the time interval of ~5-years with the difference of acquisition time less than 3 months, for a consistent and comparable RSEI from each scene.

Moreover, there may be some limitations for the RSEI. First is the sensitivity to the time span of the image acquisition date [74,75]. Since some involved indicators are affected by intra-annual phenology, the differences for the RSEI can hardly be ignored when the images with the time interval are over a season. Second, is the indicative significance. The RSEI, composed of multiple single indices, can comprehensively measure the general trend of overall ecological quality, while the index itself may be more ambiguous compared to other indicators such as area or temperature. Third is the rationale for the selected indices and weights. This is an inevitable issue for ecological assessment. For example, whether only NDVI can thoroughly reflect the Greenness of regional ecological quality is worthy of discussion. Recently, some other remote-sensing-based indices were also proposed to deal with the different surface conditions. For example, Yang et al. [30] substituted the NDVI for the Vegetation Cover (VC) and Vegetative Health Index (VHI) and proposed a Comprehensive Ecological Evaluation Index (CEEI), which may have the advantage in evaluating the Greenness over the RSEI. Sun et al. [76] proposed a time-series RSEI (ts-RSEI) that can be used for conducting comparable ecological quality assessment at any time interval. Moreover, Musse et al. [77] combined remote sensing and census data and developed an Urban Environmental Quality Index (UEQI), which has been widely used in evaluating the quality of urban ecosystems. In the future, we will further explore the suitability of the above indices for ecological assessment in the HGBA.

6. Conclusions

Large-scale ecological monitoring and assessment to reveal the regional heterogeny is important for implementing effective environmental management. Recently, the RSEI has been proven reliable for evaluating the ecological environment by means of integrating multiple remote sensing indices. In this study, we proposed a percentile de-noising method to overcome the problem of widespread noises from large-scale regions. Simultaneously, instead of normalizing remote sensing indices by scene, we normalized them by all the periods involved to ensure the consistent and comparable RSEI for long-term monitoring. On the basis of this, the dynamics of the ecological quality in the Hangzhou Greater Bay Area (HGBA) during 1995–2020 were tracked using the modified RSEI.

(1) The ecological quality of the HGBA was good but divergent for each prefecture-level city. The ecological quality degraded for the northern cities, maintained for the southern cities, and improved for the southwest city. (2) The areas of mountains, downtowns, and coastal wetlands were the hot spots for the improvement of ecological quality, and urbanization and tidal flat reclamation were the main reasons for the degradation of ecological quality. (3) As the ecological quality improved for forests and urban areas but degraded for farmland, the ecological gap has become lower among human-dominant environments but higher between the conversion from nature- to human-dominant environments. These findings can provide useful guidance for high-quality socioeconomic development for the HGBA.

For the broader large-scale RSEI-based applications, a shorter time-interval is suggested to identify the general trend and turning points of ecological quality more accurately, which is expected to be achieved by incorporating more images from different satellites.