## 2 Materials and Methods
### Methods  
In the data processing stage of this study, spatial aggregation processing was applied to the 30 m resolution GAIA binary data, converting it into coverage proportion data in 1 km grid units. The specific methodology is as follows (Figure 2): First, in the GAIA binary data, \"1\" represents impervious surface and \"0\" represents non-impervious surface. We used a sliding window with 1000 m sides to aggregate the original raster, counting the number of pixels with a value of 1 within each 1 km grid and dividing by the total number of pixels contained in that grid, thereby obtaining the coverage proportion of that land type within the 1 km grid (value range 0-1). Additionally, to ensure the accuracy of aggregation results, mask processing was applied to invalid values or edge areas in the original data before processing to avoid interference with statistical results. The final output of 1 km resolution raster data can better reflect the spatial distribution density of target land types at the macro scale, providing fundamental data support for regional-scale land use analysis and modeling. In the impact analysis and vulnerability assessment stage, ecological sensitivity factors such as HQI, FVC, LST, population density, and slope are comprehensively considered to construct an ecological environment change weight factor system, and based on weighted calculations, the Ecological Vulnerability Index  
\begin{table}
\begin{tabular}{c c c c} \hline
**Data Name** & **Description** & **Attributes** & **Usage** \\ \hline \multirow{8}{*}{Global Artificial Impervious Area (GAIA)} & \multicolumn{2}{c}{Mapped global urban impervious surface data} & \multirow{8}{*}{30 m} & \multirow{8}{*}{Impervious Surface Coverage (ISC) construction} \\  & (1985–2020) as annual maps. & & \\ \cline{1-1}  & Through spatial interpolation & & \\ \cline{1-1}  & and exclusion–inclusion & & \\ \cline{1-1}  & algorithms to achieve rapid & & \\ \cline{1-1}  & annual mapping of impervious & & \\ \cline{1-1}  & surfaces, which can reasonably & & \\ \cline{1-1}  & reflect the spatiotemporal & & \\ \cline{1-1}  & changes in impervious & & \\ \cline{1-1}  & surface sequences. & & \\ \hline \multirow{8}{*}{Land Use Data} & \multicolumn{2}{c}{CLCD data products from 1985 to 2023 are derived from Landsat-based Chinese annual land use products. The dataset achieved an overall accuracy of 79.31\%, averaged across all years based on 5463 validation samples.} & \multirow{8}{*}{30 m} & \multirow{8}{*}{Calculate HQI Calculate the landscape index} \\  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \hline \multirow{8}{*}{MOD11A1 V6} & \multicolumn{2}{c}{Based on Google Earth Engine, the MODIS-provided Quality Control (QC) bands are applied to filter out low-quality pixels affected by clouds, view angle anomalies, and sensor issues.} & \multirow{8}{*}{1000 m} & \multirow{8}{*}{Calculate FVC and LST} \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \hline \multirow{2}{*}{DEM} & \multicolumn{2}{c}{NASA global 30 m SRTM elevation DEM data} & \multirow{2}{*}{30 m} & \multirow{2}{*}{Calculate the slope} \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \cline{1-1}  & \multicolumn{2}{c}{} & & \\ \end{tabular}
\end{table}
Table 1: The details of the used data.  
(EVI) is computed to assess the impact of urban expansion on ecosystems and regional vulnerability patterns.