Evaluation of ecological quality in southeast Chongqing based on modified remote sensing ecological index

Study area

Southeastern Chongqing, China (107° 14′–109° 19′ E, 28° 9′–30° 32′ N), has an area of about 19,800 km² (Fig. 1). The study area has a subtropical monsoon climate. And the area has four distinct seasons, with an annual average temperature of 16.2 °C and abundant rainfall, with an average annual rainfall of 1209 mm. This region is located in the central part of the Wuling mountains, which is characterized by medium and low mountainous landforms, with an average altitude of greater than 1000 m. The water system (the Wujiang River system) in the study area is well developed, with a large drainage area and rich groundwater resources. The soil is dominated by yellow soil and limestone soil, and the sensitivity to soil erosion is high. The district exhibits the typical ecological fragility of karst areas, with barren soil, fragmented surfaces, a single community, and a low ecological carrying capacity. The area includes six counties: Qianjiang district, Shizhu Tujia Autonomous county, Xiushan Tujia and Miao Autonomous county, Youyang Tujia and Miao Autonomous county, Wulong district, and Pengshui Miao and Tujia Autonomous county. The coverage rates of the carbonatite layers in these counties are 42.11, 67.77, 25.70, 34.80, 59.70 and 88.46%, respectively³⁸, and the average coverage of the carbonatite layers is 53.09%, making this a representative area of karst rocky desertification.

Data and image pre-processing

In the study, the remote sensing data were obtained from the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/), including landsat-5 thematic mapper (TM) images acquired in 2001, 2006 and 2011 and Landsat-8 operational land imager (OLI) images obtained in 2016 and 2021 (Table 1). The spatial resolution is 30 m. In order to ensure the comparability of spectral characteristics, the data collection was conducted from May to September when the vegetation grew better. In order to meet the usage requirements, the cloud cover of each image used is below 10%. For the images with poor quality, the adjacent years were selected for replacement. The difference in ecological quality between adjacent years in the same region was not particularly large. In order to represent the actual situation of the ecological environment quality in the target year as much as possible, we tried to minimize the replaced part in each target year. A total of 20 images were collected in this study. The images downloaded were all L1T products, which had undergone systematic radiometric correction and geometric correction, so precise geometric correction was no longer performed. Before the subsequent processing, all 20 images were preprocessed by radiometric calibration, atmospheric correction, image mosaicking and cropping. Then these images were calculated to obtain NDVI, WET, NDBSI, LST and RI. And based on the preprocessed Landsat images, support vector machine classification was performed to obtain the land use (LU) status.

The topographical data included the elevation (EV) and slope (SP) data. Among them, the elevation data was provided by the official website of the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/). And the slope data was calculated from the elevation data. The meteorological data, including the monthly average temperature (MT), monthly mean precipitation (PR), monthly even relative humidity (RH), and monthly total sunshine hours (SH) from May to September of the target year, were got from the China Meteorological Data Network (http://data.cma.cn/). In addition, socioeconomic data, including the population density (PD) and gross domestic product (GDP), were obtained from the statistical yearbooks of each district and county in the study area. The nighttime light (NTL) data were obtained from the National Oceanic and Atmospheric Administration (NOAA, https://www.noaa.gov/). The above data and LU were used as the influencing factors of ecological quality to analyze the reasons for the change of local ecological environment quality. The statistical data and monitoring data of each evaluation index used to construct the EI come from the statistical yearbooks, water resources bulletin and soil and water conservation bulletin of each district and county.

Methodology

Study framework

A framework was developed for evaluating the ecological quality in southeastern Chongqing from 2001 to 2021 in the study. And the framework included three parts: data preparation, construction of the MRSEI, and the analysis of the ecological status in the region. Figure 2 presents the detailed information about the framework. The operations of band calculation, normalization and PCA were all carried out using the ENVI 5.3 software (https://www.harrisgeospatial.com).

Indicators used in MRSEI

The greenness, humidity, heat, dryness, and degree of rocky desertification were used to construct the MRSEI. The NDVI³⁹ was chosen to characterize the greenness. The humidity component acquired from the tasseled cap transformation (WET)⁴⁰ was selected to represent the humidity. The LST⁴¹ was used to represent the heat, the normalized difference build-up soil index (NDBSI)⁴² was used to characterize the dryness. The RI was applied to characterize the degree of rocky desertification.

The NDVI is an important indicator for monitoring the physical and chemical properties of vegetation, and it can be employed to calculate the vegetation coverage, leaf area index, and so on¹⁹. In addition, it eliminates some radiation errors and has a stronger response to surface vegetation. It has been widely used in vegetation remote sensing monitoring. The equation for calculating the NDVI is as follows³⁹:

where ({uprho }_{{{text{NIR}}}}) is the reflectance of the near-infrared band and ({uprho }_{{{text{Red}}}}) refers to the reflectance of the red band corresponding to each image.

The WET can effectively reflect the humidity conditions of the surface vegetation, water, and soil, and can reveal the changes in the ecological environment, such as soil degradation. Therefore, it is commonly used in ecological environment monitoring⁴³. The WET can be expressed as^40,43:

where ({uprho }_{{text{i}}} ,) is the reflectance of band i.

The NDBSI is expressed as the average of two indicators, the bare soil index (SI)⁴⁴ and the index-based built-up index (IBI)⁴⁵. It can be applied to characterize the dryness. The calculation formulas are^44,45:

where ({uprho }_{{text{i}}} ,) is the reflectance of band i.

The LST is closely related to natural processes and human phenomena such as crop yield, vegetation growth and distribution, surface water cycle, etc. It can well reflect the state of the surface ecological environment. The atmospheric correction method is used to invert the LST here^46,47, it can be expressed as:

where L is the radiation value in the thermal infrared band, DN is the gray value, gain and bias is the gain value and offset value of the L-band, which was got from the image header file. And T is the temperature value at the sensor; K1 and K2 are calibration parameters respectively (for TM, K1 = 607.76 W/(m² sr μm), K2 = 1260.56 K; for TIRS, K1 = 774.89 W/(m² sr μm), K2 = 1321.08 K); λ is the central wavelength of thermal infrared band; α = 1.438 × 10⁻² m K. ε is the surface emissivity and the value is estimated by the vegetation index mixture model^48,49. It is calculated as follows:

where VFC is the vegetation fractional cover, ({text{NDVI}}_{{{text{Veg}}}}) is the NDVI of the pixel covered by full vegetation and the pixels with NDVI > 0.72 are regarded as pure vegetation pixels; ({text{NDVI}}_{{{text{Soil}}}}) is the NDVI of the bare pixel and the pixels with NDVI < 0.10 are regarded as pure bare soil pixels. Pixels with 0.10 ≤ NDVI ≤ 0.72 are regarded as mixed pixels, Calculated by Eq. (12). ({text{d}}_{{upvarepsilon }}) is the error caused by terrain relief and satellite tilt observations ({upvarepsilon }_{{text{v}}}) is the radiance of pure vegetation pixels, and its value is 0.985, ({upvarepsilon }_{{text{s}}}) is the radiance of pure bare soil pixels, and its value is 0.960, F is the terrain factor, generally 0.55.

At present, there is no definite rocky desertification evaluation system, and remote sensing rocky desertification monitoring needs to consider the characteristics of the remote sensing technology and the application requirements for the rocky desertification evaluation⁵⁰. The bedrock exposure (Fr) and VFC are the most intuitive ground performance features of rocky desertification, and they are the key indicators for rocky desertification evaluation⁵¹. The slope affects the soil erosion, which is one of the reasons for rocky desertification⁵². Therefore, considering the geographical conditions of the study area and the availability of data, combined with previous research results^53,54, VFC, Fr and slope were selected as evaluation indicators. The PCA was performed on these indicators to construct the RI. The slope can be obtained by calculation from the DEM. The vegetation coverage was estimated based on NDVI by the pixel dichotomy method⁵⁵. The normalized difference rock index (NDRI) was used to estimate the exposure of bedrock by pixel dichotomy⁵⁶. The VFC is calculated by Eq. (12), the NDVI with a 5% cumulative frequency was selected as the NDVI for the bare soil pixels, while the NDVI with a 95% cumulative frequency was chosen as the NDVI for the full vegetation coverage pixels here⁵⁵.

The expression for NDRI is⁵⁶:

where ({uprho }_{{{text{SWIR}}}}) and (uprho_{NIR}) represent the reflectivity of the short-wave infrared band and the near-infrared band corresponding to each image, respectively.

And the expression for Fr is⁵⁶:

where ({text{NDRI}}_{{{text{Rock}}}}) is the NDRI of the pixel with full rock coverage, and ({text{NDRI}}_{0}) is the NDVI of the pixel without rock coverage. In the study, the values of ({text{NDRI}}_{{{text{Rock}}}}) and ({text{NDRI}}_{{0}}) are the maximum and minimum values within the confidence interval of 95% confidence in the image, respectively.

Due to the different units of the three indicators used to establish the RI, for the sake of facilitating comparison and calculation and to eliminate the weight imbalance caused by their different units, the three indicators were normalized to a range of 0–1. The expression for the normalized value of the indicators is:

where ({text{I}}) is the value for a pixel, ({text{I}}_{{{text{max}}}}) and ({text{I}}_{{{text{min}}}}) represent the maximum and minimum values of the pixel, respectively.

The band synthesis and PCA of the three indexes after normalization was carried out, and then, the initially RI (RI₀) was obtained as follows:

In order to make the RI correspond to the degree of rocky desertification, for the negative correlations, ({text{RI}}_{0}) was acquired using the following formula:

In addition, RI is obtained by normalizing ({text{RI}}_{{0}}) by formula 15, which is convenient for analysis and comparison between different years.

Calculation of MRSEI

The five component metrics (NDVI, WET, LST, NDBSI, and RI) were normalized to eliminate the dimensional difference, and then band synthesis was carried out. The water system in the study area is well developed. The water area was extracted using the modified normalized difference water index (MNDWI)⁵⁷ to prevent large-scale water from affecting the reflection of the humidity indicators. Then, it was removed by applying the mask. PCA was performed on the synthetic images after removing the water range, and the initial MRSEI (({text{MRSEI}}_{0})) was constructed. And the equation for this index is:

where ({text{PC}}_{{1}}) is the first principal component obtained by PCA. In addition, for the negative correlation results, a positive and negative transposition operation was needed. The initial modified remote sensing ecological index is:

Similarly, it was necessary to normalize the ({text{MRSEI}}_{0}) to obtain the MRSEI for eliminating the influence of dimension. And the quality of the ecological environment was proportional to the size of its value. The larger the MRSEI, the better the ecological quality.

Comprehensiveness verification of MRSEI

The correlation degree was used to determine the comprehensive representative degree of the MRSEI. Its value is between 0 and 1, and the larger the value, the better the comprehensive representation of MRSEI⁵⁸. The formula is:

where ({overline{text{C}}})_p is the average correlation; p, q, r and s are the indicators for the correlation analysis; ({text{C}}_{{text{q}}}), ({text{C}}_{{text{r}}}), and ({text{C}}_{{text{s}}}) are the correlation coefficients between p and q, r and s during the same period, respectively, and n is the number of indicators for the correlation analysis.

Geographical detector method

The geographical detector aims to probing spatial heterogeneity and uncovering its causes. It can not only analyze the relationship between a single independent variable and a dependent variable, but it can also analyze the interactions between two independent variables. In addition, it can accurately determine the strength, direction, and linear or nonlinear relationship between the two variables, and thus, it can reflect the driving force behind the dependent variable to solve practical problems^59,60. The factor detector is mainly applied to detect the spatial differentiation of the dependent variables and the degree of interpretation of the independent variables to the spatial heterogeneity of the dependent variables. The degree of explanation of the independent variable corresponding to the variable is expressed by the q value^61,62,63. In this study, the factor detector was mainly provided to analyze the influences of the various factors on the spatial differentiation of the quality of the ecological environment in southeastern Chongqing. And it can be expressed as:

where q is the influence of the affecting factors on the spatial heterogeneity of the ecological status, ({upsigma }^{{2}}) is the variance of the MRSEI in the entire region, ({upsigma }_{{text{h}}}^{{2}}) is the variance of the MRSEI in sub-level region h, N is the total number of samples in the region, ({text{N}}_{{text{h}}}) is the number of sub-level regional samples, and L is the total number of sub-level regions.

Construction of EI

According to Technical Specifications for Evaluation of Ecological Environment Status (HJ 192-2015)²³, five sub-index and 20 sub-indicators were selected to construct EI. Taking the sub-indicators as the basic data, the sub-index was obtained by calculation, and then the EI is obtained by weighted calculation of the sub-index. The expression for EI is as follows:

where HOI is the Habitat Quality Index, VCI is the Vegetation Cover Index, WND is the Water Network Density Index, LSI is the Land Stress Index, PLI is the Pollution Load Index. The weights of the sub-indicators are determined by AHP, and they are shown in Table 2.

Source: Ecology - nature.com