Assessment and simulation of eco-environmental quality changes in rapid rural urbanization: Xiong’an New Area, China

RSEI ecological index

RSEI is a regional EEQ assessment index that is completely based on remote sensing technology. Compared to traditional ecological indices, the most notable improvement in RSEI is its selection of four important indicators, namely greenness, wetness, heat, and dryness, which directly reflect the ecological quality¹⁴. These indicators are used for dimensionality reduction calculation, which not only avoid one-sided evaluation of a single indicator, but also integrate the information from multiple indicators into a concise visual index or result. The functional form of RSEI can be expressed as:

NDVI is capable of reliably assessing vegetation quality and productivity, and it has the advantages of being concise and easy to calculate³⁰. NDVI is currently the most widely used remote sensing vegetation index. In this study, greenness indicator in RSEI is represented by NDVI. NDVI is expressed as:

where NIR and Red are the near-infrared (NIR) band and red band of Landsat-8 reflectance images, respectively.

Tasseled cap transformation (TCT) is a valuable tool for compressing spectral data into a few components that correspond to the physical characteristics of a scene, while minimizing information loss³¹. Among them, the wetness component is associated with soil and plant moistures, water and other moist features³². We used wetness component to characterize the wetness indicator in RSEI. The wetness component of TCT is expressed as³¹:

where Blue, Green, Red, NIR, SWIR1 and SWIR2 corresponding to the blue, green, red, NIR, short-infrared 1 (SWIR1), and short-infrared 2 (SWIR2) bands of Landsat-8 reflectance images, respectively.

Dryness indicator is mainly used to quantitatively depict the “dryness” characteristics on the land surface caused by the replacement of natural cover with artificial surface due to human activities in a given region. In studies focusing on urban built-up area, commonly used indices such as the building/built-up index or the impervious surface index were employed to represent this indicator. However, unlike urban built-up areas, Xiong’an New Area is located in rural areas. Apart from the rural settlements, there are currently large areas of bare soil in this region, which are primarily caused by land leveling activities during the early stages of engineering construction and crop harvesting. In this study, therefore, the normalized difference impervious surface index (NDISI)³³ and the normalized difference soil index (NDSI)³⁴ were combined to construct the normalized difference impervious surface and soil index (NDISSI) to represent the Dryness indicator in RSEI. The NDISSI is expressed as:

where NDISI_nor and NDSI_nor are normalized NDISI and NDSI, respectively. NDISI and NDSI are expressed as:

where TIR is the thermal band; WI represents water index-derived band and can be represented by the modification of normalized difference water index (MNDWI)³⁵.

The normalization algorithm can be expressed as:

where I_nor is the normalized index or indicator; I_i represents the value of the index or indicator in pixel i; I_max and I_min represent the maximum and minimum values of the index or indicator, respectively.

Heat indicator in RSEI is represented by LST. In this study, Landsat-8 TIRS band 10 was utilized. The single channel (SC) algorithm proposed by Jiménez-Muñoz and Sobrino³⁶, Jiménez-Muñoz et al.³⁷, and Cristóbal et al.³⁸ was adopted to retrieve LST. The SC algorithm is expressed as follows:

where L is the at-sensor spectral radiance of TIRS band 10; ε is the land surface emissivity; γ and δ are two parameters dependent on Planck’s function; b_γ = 1324 for TIRS band 10; ψ₁, ψ₂, and ψ₃ are the atmospheric functions calculated as:

where τ is the atmospheric transmissivity, and L^↑ and L^↓ are the upwelling and downwelling atmospheric radiance, respectively. T is the at-sensor brightness temperature calculated as follows:

where K₁ and K₂ are the band-specific thermal conversion constants, K₁ = 774.89 (W/m²/sr/µm) and K₂ = 1321.08 K for TIRS band 10.

The PCA algorithm, a multivariate statistical method, was utilized in the construction of RSEI. Through rotating the coordinate axes of feature spectral space to maximize the removal of correlation among different indicators, the key information of multivariate is concentrated into the first few principal components, such as the first principal component (PC1). The advantage of PCA lies in its utilization of covariance matrix for automatic and objective calculation of the weight distribution for each indicator, and in the objective determination of the contribution of each principal component through the corresponding eigenvalue. Any weighting structure is susceptible to criticism, as assigning weights is a subjective value-dependent process³⁹. While, PCA effectively avoids result biases caused by arbitrary subjective weighting of indicators¹⁴. The RSEI function in Eq. (1) can be further divided into the following two steps:

It is important to note that due to the different units of measurement among indicators, it is necessary to normalize each indicator using Eq. (7) before performing PCA. The normalization process ensures that all the indicators are on a unified scale within the range of 0 to 1. Additionally, if there are large water bodies in the study area, they may have an impact on the analysis of indicators loading in PCA¹⁴. In such cases, the MNDWI can be used to mask out the water information. The resulting RSEI is referred to as the normalized RSEI₀ (RSEI_{0 nor}), which has value range from 0 to 1. A higher value of RSEI, that is closer to 1, indicates a better EEQ, while a value of 0 denotes an extremely poor one.

XGBoost algorithm

The XGBoost algorithm, initially proposed by Chen et al.⁴⁰, is an optimized distributed gradient boosting library that has found extensive applications in various fields, including the built environment^41,42, natural geography⁴³, and agricultural and biological⁴⁴. It has demonstrated high accuracy in these domains. As a ensemble decision trees, XGBoost initializes with a base prediction and then iteratively constructs trees to fit the residuals (i.e., the differences between the observed and predicted values) from the previous iterations. This process repeats multiple times, allowing XGBoost to achieve highly accurate predictions⁴⁵. One notable advantage of XGBoost is its ability to refine its predictions through multiple rounds of iterations over the residuals⁴².

XGBoost is a highly scalable machine learning system developed from the concept of Gradient Boost Decision Tree (GBDT)⁴³, with multiple adjusted hyperparameters. Configuring these hyperparameters appropriately is a crucial step in achieving model accuracy. In this study, the seven hyperparameters, namely eta, gamma, max_depth, min_child_weight, subsample, colsample_bytree, and nrounds, were iteratively calculated within a preset hyperparameters tuning space (Table 2). A 5-fold cross validation approach⁴⁶ was employed to determine the optimal hyperparameters combination that minimizes the model’s error. Furthermore, the Shapley Additive exPlanation (SHAP) analysis is based on game theory and relies on SHAP values to determine the importance of individual independent variables⁴⁷. It takes into account their marginal contribution to the XGBoost regression model outcome.

Multi-scenario simulation

This study was based on the 2021 conditions of study area as the basis for multi-scenario simulation. The illustration of multi-scenario simulation workflow is shown in Fig. 2. First, the built-up areas and bare soil areas under development within the study area were extracted using the NDISSI threshold, forming the initial region A1. Then, concentric buffer zones of 200 m, 400 m, and 600 m were created around A1, forming buffer regions A2, A3, A4, respectively. Next, relevant urban-rural construction features of the current study area in 2021 (S0) was obtained, including NDVI_nor, NDISI_nor, percentage of the building rooftop (PB), population density (PopD), and building height (BH). These construction features were then progressively changed by 5% starting from region A1 with a step length of 200 m to simulate the possibility of urban construction in multiple scenarios. Four scenarios were simulated as follows: Scenario One (S1) involved a 5% progressive change in the urban-rural construction features of region A1 based on current situation (S0), specifically a decrease of 5% in NDVI_nor and an increase of 5% in NDISI_nor, PB, PopD, and BH, while other non-A1 regions remain unchanged. The modified urban-rural construction features were used in the XGBoost regression model to simulate the RSEI, resulting in S1. Similarly, Scenario Two (S2) was built upon S1 with a progressive change of 5% in the construction features for regions A1 and A2. Likewise, Scenario Three (S3) was based on S2 and involved a change of 5% in the construction features for regions A1, A2, and A3. Finally, Scenario Four (S4) was built upon S3 with a change of 5% in the construction features for regions A1, A2, A3, and A4.

Illustration of multi-scenario simulation workflow.

It is worth noting that during the simulation, if NDVI_nor falls below 0, it is adjusted to 0. If NDISI_nor and PB exceed values of 1 or 100%, they are constrained to 1 or 100% accordingly. Additionally, in the newly added buffer region, if PopD and BH are both 0, the minimum non-zone values within that region are utilized as substitute.