Research on Key Influencing Factors of Ecological Environment Quality in Barcelona Metropolitan Region Based on Remote Sensing

1. Introduction

In the past few decades, environmental assessment, as an important part of environmental management, has gradually attracted people’s attention [1]. People have realized the impact of human activities on environmental change and recognized the key role of environmental assessment in studying the process of environmental change [2]. In today’s world, with the rapid development of urbanization and social and economic construction, the ecological environment has been severely damaged by human activities, and ecological problems such as soil erosion, land salinization and desertification, and a sharp decline in plant and animal diversity have also increased [3,4]. At the same time, with the deterioration of the ecological environment, natural disasters and extreme climates such as floods, droughts, global warming, and high temperature heat waves have occurred frequently [5,6,7]. A good ecological environment is the basis for sustainable social development and human survival. The deterioration of the ecological environment has seriously affected the development of human society and people’s health. We must pay attention to solving ecological and environmental problems while maintaining urban construction and economic development, as well as the continuous progress of science and technology. Only in this way can we achieve a harmonious coexistence between man and nature, reduce the frequent occurrence of natural disasters and extreme weather, and protect people’s health and the sustainable use of natural resources [8].

The scientific evaluation of the ecological environment and the rational optimization and regulation of its pattern are not only hot issues in the field of climate and environment research, but also urgent needs for ecological economic development and ecological civilization construction [9]. At present, satellite remote sensing earth observation systems have been widely used in the field of ecological environment due to their advantages of macro, rapid and real-time. The use of various remote sensing indicators to monitor and evaluate natural forests, grasslands, urban artificial ground, water bodies, and even ecosystems in the entire region has become an important part of the field of ecological environment protection [10]. With the development of “3S” technology (remote sensing, geographic information system and global positioning system), many scholars have begun to use RS and GIS technology to obtain ground information, combined with mathematical models and statistical analysis methods, from single factors to multiple factors, to conduct comprehensive evaluation of the ecological environment in different regions and scales [11,12]. In 2009, Ivits et al. used SPOT image data and the vegetation index NDVI as an ecological indicator to evaluate the ecological environment of forests and farmland bird habitats in China [13]. Marllu J et al. used a variety of ecological factors such as vegetation sensitivity index and ecological isolation index to construct an evaluation system to evaluate the quality of urban ecological environment [13]. Studies have shown that the ecological quality status obtained by this method has reference significance for relevant decisions on urban land use planning. In 2013, Xu Hanqiu proposed a new remote sensing ecological index (RSEI) based on remote sensing information technology. The index couples four indicator factors: vegetation index, humidity component, surface temperature and soil index, and then uses principal component transformation to integrate more complete information in the region [15]. This study used the RSEI index and China’s traditional EI index to analyze the ecological environment in areas with severe soil erosion in Fujian Province, providing prediction and dynamic change analysis, and making up for the shortcomings of the EI index. In 2021, Wang Zi selected four phases of Landsat series remote sensing images as data sources, added research indicators that meet the characteristics of saline-alkali land in the research object, and comprehensively evaluated the ecological environment quality of the Chahan Nur Basin from 1992 to 2020, and analyzed the driving factors of changes in ecological environment quality [8]. In 2022, Hou Yifeng et al. analyzed the evolution of land cover changes in the Tarim River Basin over the past 30 years and explored the contribution of transformations between different land types to the ecological environment [16].

Normalized Difference Vegetation Index (NDVI) can be used to detect plant biomass, vegetation coverage, growth status, etc. It is the best indicator factor reflecting the vulnerability of regional ecological environment and is also the most widely used index in vegetation remote sensing monitoring research [17,18]. Vegetation occupies a major position in the earth’s terrestrial ecosystem and is an important medium for the material and energy cycle of the earth’s ecosystem [19]. NDVI changes play a guiding role in the study of regional ecological environment evolution [20]. High NDVI values usually indicate high vegetation density and healthy growth, which helps to improve the resilience and biodiversity of the ecosystem [21]. In 2004, Wessels et al. studied the land degradation effects caused by human activities in the former homeland area of northern South Africa and found that NDVI values can indirectly reflect the health of the soil and thus judge the quality of the local ecosystem. There is a positive correlation between them [22]. In addition, NDVI can effectively reflect the seasonal changes and long-term trends of vegetation and is effectively applicable to the monitoring and evaluation of different ecosystems [23]. Goward and Prince studied the relationship between NDVI and vegetation net primary productivity (NPP) and found that NDVI can be used as an alternative indicator of NPP to assess the health of the ecological environment [24]. De Jong et al. studied the impact of grazing on vegetation and erosion in Paraguay and found that NDVI monitoring can effectively evaluate the effect of ecosystem restoration [25]. In 2020, Ren Jie systematically studied the driving factors and influence of the evolution of the surface ecological environment in the western region of the Yili Basin based on NDV and combined with the ILUCC theory, and found that the lower the NDVI, the more sensitive the ecological quality [26]. However, at present, scholars have not yet formed a unified standard for evaluating the quality of the ecological environment. The commonly used land weight index is only a numerical value [8,11,12,13,14,15,16,23,24,25,26], which can only generally describe the ecological status of a region in a certain period, and cannot analyze the spatial changes of the ecological environment in different periods. At the same time, different scholars consider this weight index from different perspectives. Therefore, trying to find a more scientific and reasonable indicator system to evaluate the ecological environment quality of the metropolis, which can be applied to changes in land use types in different periods and regions, is one of the important tasks of this study.

Many of the above studies have shown that NDVI is an extremely important factor that affects and reflects the ecological environment quality of a region. Therefore, it is particularly important to analyze the factors that can have the most direct impact on NDVI and thus indirectly regulate ecological quality. Many studies have shown that temperature, soil properties, topography and human activities are factors that can directly affect NDVI changes [27,28,29,30,31,32], but there is a more intimate and complex relationship between NDVI and precipitation. Generally speaking, sufficient precipitation can promote plant growth, increase biomass and NDVI values. Conversely, drought and insufficient precipitation will reduce NDVI values [21,33]. However, the relationship between precipitation and NDVI still has obvious differences in different regions and times. Mu Shaojie et al. found [34] that there is a significant time lag effect in the response of NDVI to precipitation in grassland areas. Liu et al. studied the vegetation dynamics of the Amazon rainforest and found that even in areas with abundant precipitation, seasonal precipitation changes can still cause significant fluctuations in NDVI [35]. Xu Yue et al. also proposed that the response rate of NDVI to precipitation changes in different seasons is also different [36]. In previous works, the authors also believed that the effect of precipitation on plant cover is a highly complex factor and it is difficult to determine the true relationship between them, but in fact precipitation does provide a great help to the growth of vegetation [27]. From the current research status, the relationship between precipitation and NDVI is still an unresolved topic.

Based on the above, this article will build a more scientific and reasonable evaluation system based on NDVI, which can be applied to dynamic ecological environment quality in different metropolitan areas and periods, and use this system to evaluate the ecological environment quality of Barcelona Metropolitan Region. Secondly, we will comprehensively analyze the various natural and human factors that can affect NDVI. Most importantly, the relationship between precipitation and NDVI must be deeply analyzed. We try to discover some laws and inspirations for this unresolved topic to a certain extent.

2. Materials and Methods

2.1. The Field of Study

The Barcelona Metropolitan Region (BMR, Figure 1) is located in the northeast of the Iberian Peninsula, in the center of the Mediterranean corridor connecting Spain to the rest of continental Europe, and it has a typical Mediterranean climate. It is the largest metropolitan area in Catalonia and serves as its political, economic and cultural core. It covers an area of approximately 3224.7 km² and contains 164 municipalities. With a population of around 5.4 million inhabitants, the BMR is the most densely populated metropolitan area in the European Union. The Llobregat and Besòs rivers are the two main rivers that flow through the metropolitan region of Barcelona. The coastal and pre-coastal mountain ranges delimit the coastal and pre-coastal depressions where the main population centers are located. It includes the more densely populated Barcelona Metropolitan Area (BMA), about 636km², covering the city of Barcelona and 35 surrounding towns, with a permanent population of about 4.7 million. The urban core area of Barcelona is about 100 km², with a population of more than 1.65 million, and a population density of more than 16,500 inhabitants/km². The Collserola “serralada” (the central part of the coastal mountain range) occupies a central place in the metropolitan area, making it the “green lung” around which the built spaces are located. The Barcelona Metropolitan Territorial Plan (10) pays special attention to the protection of the ecological environment and green spaces of the BMR. The green area of the entire Barcelona Metropolitan Region is about 2514km², which includes urban parks, suburban parks, forests and protected areas, accounting for about 77.4% of the total area. Based on the above situation, we chose BMR as the study area to evaluate the evolution of its ecological environment quality and analyze the impact of various factors, especially precipitation, on NDVI.

2.2. Methodology

Assessing and analyzing the quality of green ecological environment in metropolitan areas can help protect and restore urban ecosystems, balance species diversity, provide a scientific basis for urban planning, and help formulate reasonable protection and expansion strategies. Studying the changing mechanisms that affect NDVI can enhance the resilience of ecosystems, improve their ability to adapt to extreme climates, and improve the well-being of residents. Corine Land Cover (CLC) provides us with relatively accurate European land cover data with a resolution of 100 meters in 2006, 2012, and 2018 (https://land.copernicus.eu/en/products/corine-land-cover), which records 22 types of green space systems in detail. Based on this dataset, we carried out this study in the following steps.

• First of all, evaluate the Ecological Environment Quality Index (EQI) of BMR. We try to find a more scientific and reasonable indicator system to evaluate the ecological environment quality of metropolises, which can adapt to changes in time and the evaluated regions. To complete this evaluation, we need to achieve the following aspects.

• For the convenience of analysis, we classified the land use coverage of BMR according to the land type classification of CLC.
• Then, we extracted the difference (NDVIaverage-NDBIaverage) of the annual average values of NDVI and NDBI for each type of land in 2006, 2012 and 2018 after reclassification. The MODIS_MOD13Q1 (https://lpdaac.usgs.gov/products/mod13q1v006/) dataset provides relevant remote sensing image data with a resolution of 250 meters.
• Set EQI evaluation parameters. Normalized Difference Vegetation Index (NDVI) is the best indicator of regional ecological environmental vulnerability. In previous studies, the authors believed that the tall tree canopy in summer would make the NDVI estimate higher than the actual value, while the long-wave radiation from the ground in winter was not significantly affected by the tree canopy [37]. Therefore, this study hopes to determine the ecological quality weights of different land uses by using the NDVI in summer and winter.

In addition, we also need to consider the impact of Normalized Difference Built-up Index (NDBI) in ecological benefits. Negative NDBI values represent water bodies, high values represent urban built-up areas, and the NDBI value of vegetation is relatively low. Therefore, we can assume that the ecological quality of a region is inversely proportional to the NDBI index. If we subtract NDBI from NDVI, the result can more truly reflect the ecological effect of the region.

For each year, we homogenize the difference results obtained in the previous step using Formula (1), setting the minimum value of each year to 0 and the maximum value to 1 as the evaluation parameter (E_i) of ecological environmental quality.

$${\mathrm{E}}_{\mathrm{i}}={\displaystyle \frac{{\mathrm{X}}_0-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{\mathrm{X}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}},$$

(1)

where

E_i is the dimensionless ecological environment quality weight index;

X₀ is the NDVIaverage-NDBIaverage of each type of land in that year;

X_max is the maximum value of NDVIaverage-NDBIaverage of various types of land in that year;

X_min is the minimum value of NDVIaverage-NDBIaverage of various types of land in that year.

• Finally, I evaluated the BMR’s annual EQI index based on formula (2).

$${\mathrm{E}\mathrm{Q}\mathrm{I}}_{\mathrm{t}}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{i}}\ast {\displaystyle \frac{{\mathrm{E}}_{\mathrm{i}}}{\mathrm{T}\mathrm{A}}},$$

(2)

where:

EQI_t is the ecological environment quality index in year t;

i is the land type;

A_i is the area of that type of land use;

E_i is the weight of the ecological quality index of this type of land use;

and TA is the total area of the land.

• It is important to pay attention to the land use types that have a greater impact or contribution to the evaluation results when evaluating the ecological environment quality of urban areas. It is the key land use to control the ecological quality of BMR.
• In addition, we also need to draw the distribution map of the annual evaluation results of the ecological environment quality index of BMR and the change map of the results from 2006 to 2018 to analyze the distribution and evolution direction of the ecological quality of BMR from a spatial perspective, and determine the areas where the ecological quality of BMR continues to improve.

2.

Then, we need to study the relationship between the key land use distribution and various possible influencing factors that affect the assessment results. The various types of data listed in the table need to be collected. MODIS can provide land surface temperature (LST). The urban heat island effect (UHI) will be expressed as the urban-rural temperature difference, that is, the difference in the average surface temperature between the urban built-up area and the rural area based on land cover data. Subtracting the average surface temperature of the daytime and nighttime in the rural area each year can get the approximate intensity distribution of the urban heat island effect. The larger the value, the greater the intensity of the urban heat island effect. DEM terrain data comes from SRTM (https://www.earthdata.nasa.gov/sensors/srtm), with a resolution of 30m. Impervious ground data comes from GlobeLand 30 (https://www.webmap.cn/commres.do?method=globeIndex), with a resolution of 30 meters. E-obs can provide annual European precipitation grid data, as well as daily maximum temperature (T_max) and minimum temperature (T_min), but the resolution is 1°, which is very large (https://surfobs.climate.copernicus.eu/dataaccess/access_eobs.php). Therefore, we used the kriging interpolation method based on the E-obs data to reconstruct the relevant climate map of BMR with a resolution of 1km. The night light data is divided into two parts. The data from 1992 to 2013 come from DMSP data with a resolution of 30 arc seconds and a numerical range of 0-63. The night light data after 2013 are provided by VIIRS with a resolution of 15 arc seconds, which is more accurate and very different from DMSP. In the previous work of the authors, the two datasets were calibrated and unified using the stepwise calibration method, which will not be repeated here [38]. After obtaining the above data, we established a 1km grid within the BMR range, extracted the proportion of key land use area in the grid each year as the dependent variable, and established an OLS model to analyze their importance.

3.

Once it is clear that there is a clear interaction between NDVI and the distribution of key land uses that affect EQI assessment results, we must analyze the climate factors that can affect NDVI to discover the indirect impact of climate change on ecological quality. We established a climate database from 2000 to 2023, using the annual average NDVI as the dependent variable and multiple climate data as explanatory variables to analyze the strength of the correlation between them.

4.

In the previous study, the authors found that precipitation is a highly complex control factor [27], and there is an unresolved and strong interaction between it and NDVI, but this complex relationship is difficult to explain through ordinary regression models. Therefore, we try to use the precipitation-NDVI characteristic triangle space to analyze this mechanism.

• Firstly, for each year, a precipitation-NDVI scatter plot was established. In order to facilitate the analysis of the physical meaning of each side of the triangle, we used NDVI as the X-axis and precipitation as the Y-axis. Since the numerical ranges of the two values are quite different, the spatial characteristics of the established scatter plot are not obvious, so we calculated the Ln log function for both NDVI and precipitation.
• Then, the density characteristics of the scatter plot are analyzed using nonparametric kernel density estimation (KDE). Kernel density estimation is a nonparametric method used in probability theory to estimate the probability density kernel function of an unknown random variable, proposed by Rosenblatt [39] and Emanuel [40]. Kernel density estimation can infer the overall data distribution based on a limited sample and analyze the regional location of data aggregation. It does not rely on the overall distribution and its parameters, but instead obtains structural relationships through direct estimation and derivation based on sample data [41]. This study uses Qrigin 2021 for kernel density estimation analysis.
• If it is found that the scattered points are unevenly distributed and there are scattered points with low density distribution around the cluster, we need to perform local outlier factor detection (LOF). This detection method is a density-based outlier detection algorithm, which mainly determines whether it is abnormal by calculating the outlier factor of the sample and comparing whether it is far away from the dense data [42]. Outliers usually appear relatively isolated, so samples with large LOF values have low density and are considered abnormal. In the recognition process, it is considered that when the ratio between the distribution density of the point and the overall average density is much greater than 1, it is more likely to be an outlier [43].
• After removing the abnormal outliers, we can determine the sides of the characteristic triangle. The maximum NDVI value among the retained scattered points will be used as the straight edge. For the hypotenuse, we use 0.01 as a minimum interval and find the maximum and minimum precipitation values in each interval. Then use the least squares method to perform a linear fit with Equation (3) to obtain the hypotenuse of the triangle.

$${\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{p}}_{\mathrm{i}}=\mathrm{a}\ast {\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}}+\mathrm{b},$$

(3)

where:

Precpi is the LN value of precipitation corresponding to the i-th NDVI;

NDVIi is the LN value of NDVI corresponding to the i-th NDVI;

a is the slope;

b is a constant.

• Finally, after completing the most critical triangle fitting, we can analyze the physical meaning of each side of the precipitation-NDVI characteristic triangle and the space, and further analyze the relationship between NDVI and precipitation. This allows us to indirectly understand the impact of precipitation on the green space system.

The research method is summarized in the following idea map.

Figure 2. Research idea map.

Figure 2. Research idea map.

3. Results

3.1. BMR Ecological Environment Quality Assessment

3.1.1. Green Area Is the Key Land Use that Provides Ecological Benefits

In order to establish our envisioned ecological quality assessment system, we first had to reclassify the land use types of BMR according to the CLC’s description of land cover. In addition, in order to more intuitively identify the key land that provides more ecological benefits, we summarized the 11 subdivided land types into 5 main categories (Table 1).

The evaluation parameters of the ecological environment quality of each type of land use in each period can directly reflect the contribution of land to the local ecological environment. Analyzing their changes can understand the evolution of the ecological environment benefits of land use. Figure 3 (a) shows the evaluation parameters of the ecological quality of each type of land use in BMR from 2006 to 2018. The most striking is forest land. In all periods, it provides the highest ecological benefit and has always remained at 1. The second is grassland, whose parameters have remained above 8.2 in three years. The land that provides the least benefit to the ecological environment is not fixed every year. In 2006, the ecological parameter of continuous urban land was 0, but in the second year, it suddenly increased to 0.22, and the wasteland that had performed better before became the worst. In 2018, the ecological benefits of wasteland recovered rapidly, and the parameters of industrial land dropped sharply from 0.3 to 0. The parameters of water bodies were above 0.4 in the first two years, but they were reduced by more than half in the last year. The performance of land use types with lower ecological benefits was more unstable in the three periods, but in general, the ecological parameters were the highest in 2012 and lower in 2018 than in the other two years. Of course, there were a few exceptions: transportation land and water bodies performed worse in 2012 than in other years, and grassland had higher parameters in 2018 than in 2006.

Figure 3 (b) illustrates the changes in the ecological and environmental parameters of the main categories of BMR and the overall process. It is obvious that the green area has the largest ecological benefit and far exceeds other categories. In 2012, the parameter of the green area reached a peak of about 3.33. In 2018, this parameter was the lowest, but it also exceeded 2.9. The second is construction land, whose parameter also reached the highest in 2018, about 0.95. The ecological benefit of urban land is the most unexpected, only 0.1 in the first year, but then soared to 0.43, and then quickly fell to 0.19. Overall, the ecological benefit of land was the best in 2018, reaching 5.15, and the worst in 2018, only 3.79. The evolution of parameters of various categories is also consistent with this development, except for wasteland.

3.1.2. The Green Area Is Shrinking but Becoming More Fragmented

Analyzing the development process of land use distribution that is most beneficial to ecological and environmental protection helps to analyze the attribution of changes in ecological and environmental quality and clarify the direction of environmental protection. Figure 4 seems to tell us a positive message. Between 2006 and 2018, green areas in BMR occupied the vast majority of land, and land with the highest ecological benefits became the main form of land use in the metropolitan region. Large green areas are distributed in the north, east and southwest of BMR, and the central urban area is sparsely distributed.

Figure 5 shows the evolution of the area and number of plaques of the green area of BMR. It is obvious that the area of green area of BMR continues to decrease, but the number of plaques gradually increases. In 2006, the area of green area of BMR reached 2528km², accounting for more than 77.8% of the total area. After 12 years of decline, it has decreased by about 33.3km². Overall, it is satisfactory. In the last year, BMR retained more than 76.8% of the land as green area, even if it is less than before. However, the fragmentation of green areas is getting more and more serious, from 1337 plaques in 2006 to 1377 in the last year. The density of plaques also increased, from 52.9% to 55.2%. The area is reduced, but the fragmentation increases, which means that the complexity of the spatial distribution structure of green areas in BMR has also deepened, which may be unfavorable from the perspective of ecological environmental protection.

3.1.3. The Ecological Quality Results of BMR Are Generally Appreciable, with Green Area Making the Largest Contribution

After the ecological environment quality assessment system of the metropolitan region is established, we can obtain the overall assessment results of each type of land and BMR (Table 3). Overall, the ecological quality results of BMR during the study period are quite impressive, all of which remain above 0.67, especially in 2012, as high as 0.74, which is the only time in these three years that it exceeds 0.7. The assessment results in 2018 are not as good as those in 2006, but only 0.02 less than the former. Among all land categories, as expected, the EQI assessment results of green areas with the highest ecological benefits are also the highest, which almost determines the overall results of BMR ecological quality assessment. In 2006 and 2018, the EQI of green areas was 0.68 and 0.65 respectively, accounting for more than 97% of the total results. In 2012, which was also the year with the highest EQI of green areas, its contribution to the assessment results of BMR decreased to about 92.8%, even though this value is still much higher than other lands. In other years, urban land, which has always accounted for a small proportion, became the second in 2012, accounting for about 4.2% of the total results. The proportion of construction land in the total results is relatively significant compared with the other two types of land. The ecological proportion of water bodies has been decreasing. The EQI of wasteland is the worst among all categories. The specific EQI evaluation results of various land types are detailed in Appendix B.

3.1.4. The EQI of Forests Cannot Be Ignored

Green areas contribute the most to the EQI assessment results and are the decisive factor in the ecological quality of BMR. Therefore, it is extremely necessary to analyze the EQI of various types of land contained in green areas. Figure 6 shows the EQI assessment results of various types of green land in BMR from 2006 to 2018. It is not difficult to find that forests have made outstanding contributions. In the three years, the EQI of forests has always ranked first, but it has shown a downward trend, gradually decreasing from 0.426 to 0.417. The second is farmland, which increased from 0.146 to 0.157 in 2012, reaching a peak, but fell to 0.126 in 2018. The performance of grassland is relatively stable, remaining at around 0.1 during the study period and rising slightly. The worst is recreational land, which is a green space within the city. The EQI was around 0.01 in the three years and fell to a low point in 2018.

3.1.5. BMR Ecological Environment Quality Distribution Needs to Be Improved

The spatial distribution map of EQI assessment results can help us analyze the distribution of environmental quality in the BMR territory and identify areas with good ecological conditions and areas in urgent need of improvement. From Figure 7 (a, b and c), even though BMR has a large area of green areas with high ecological parameters, almost half of its land is in a low EQI situation in three years. The situation in the core Barcelona is the most pessimistic, with only a small area in the north with good ecological quality, and the rest of the EQI is very low. The situation in BMA is slightly better than Barcelona, but it is also dominated by poor ecological quality surfaces. The assessment results of small areas in the west, north and northeast of BMA are good. For BMR, the ecological quality of the coastal area, a large area in the east, and a continuous area in the north is very poor, especially Barcelona and its surrounding areas. Its northern and eastern ends have very high EQIs, and high assessment results appear in the north, central south and northwest, but the distribution is relatively scattered.

In order to more intuitively analyze the evolution of BMR’s EQI between 2006 and 2018, we drew a map of changes in ecological environment quality (d). From the map, we can see that most of BMR’s ecological quality is in a state of unchanged or improved. Barcelona’s situation is still the most unsatisfactory, with most of it deteriorating and unchanged, and only a few areas improving. The areas where the ecological quality of BMA has deteriorated are mostly concentrated in the southwest, and the surface in the north and east has improved significantly. The western end, south-central, north-central and northeastern end of BMR can be regarded as four areas where the deterioration of ecological quality is more concentrated. Its northern performance is outstanding, and the EQI has improved significantly. The rest are mostly in an improved and unchanged state.

3.2. Analysis of Factors Affecting Green Area Distribution

In order to explore the influence of various factors on the distribution of green areas, which plays the most critical role in controlling the EQI results, we established a grid with a side length of 1 km. With the area ratio of green areas in each grid as the dependent variable, longitude, latitude, distance from the sea, aspect, altitude, slope, average NDVI in winter and summer, annual precipitation, daytime LST, nighttime LST, difference between daytime and nighttime LST, maximum temperature (T_max), minimum temperature (T_min), difference between maximum and minimum temperature (T_max-T_min), NDBI, daytime and nighttime urban heat island effect, impervious ground and artificial ground as independent variables, three OLS regression models were established for the data of 2006, 2012 and 2018 (Table 4), with 10 more significant variables, R² reached 0.97. Among the natural factors, the most significant influence on the distribution of green areas is NDVI (+), even though its explanatory power is gradually weakening. Next is annual precipitation (+), the difference between the maximum and minimum temperatures (-) and the maximum temperature (+), especially in the last year, the latter two suddenly increase in influence. Finally, the difference between day and night surface temperature (+). From the perspective of human activities, artificial ground (-) and impervious areas (-) are the main obstacles controlling the distribution of green areas, followed by NDBI (-) and the urban heat island effect occurring at night (-).

In any case, it is undisputed that artificial surfaces and impervious areas are the most fundamental factors controlling the distribution of green areas. The model results show that these two factors alone can explain about 90% of the distribution variation of green areas in BMR. However, compared with human spontaneous activities, it is more meaningful to analyze uncontrollable natural factors in this study. In 2012 and 2016, the average NDVI can explain more than 54% and 67% of the variation of green areas respectively when considering only natural factors, but in 2018, its explanatory power dropped to 35% due to the importance of the difference between the maximum and minimum temperatures and the maximum temperature. If NDVI, artificial surfaces and impervious areas work together, their influence on green areas can reach more than 90%. The reason why diurnal LST, minimum temperature and daytime urban heat island effect are not included in the optimal model is that they have high collinearity with the difference between diurnal LST, maximum temperature and the difference between maximum and minimum temperatures. The relevant Pearson correlation coefficients in the model are shown in Appendix C.

3.3. Analysis of Climate Factors Affecting NDVI Evolution

There is a most significant interaction between the average NDVI and the distribution of green areas, which plays a key role in controlling the ecological environment quality. More importantly, it is necessary to explore the climate factors that are most closely related to the development of NDVI, and then explore the indirect impact of climate change on the ecological environment quality of the metropolitan region. In order to achieve this objective, we used the annual average NDVI from 2000 to 2023 as the dependent variable, and established 5 OLS models with annual precipitation, maximum and minimum temperature, and day and night surface temperature as independent variables to analyze the linear relationship between climate and NDVI change trends. Table 5 presents the results of all model analysis in detail.

The most significant positive correlation with the mean NDVI is precipitation (sig. < 0.001), and the R² of the regression model established between them reaches 0.47, which is the highest fit among the five models. It alone explains about 50% of the annual evolution of NDVI. It is followed by minimum temperature (-) and daytime LST (-), which also have a certain influence on the change of NDVI, but their explanatory power is much worse than that of precipitation. The explanatory variables with almost no effect are nighttime LST (student’s t = 0.62) and maximum temperature (student’s t = 0.25). The Pearson coefficients of the relevant variables in all models are shown in Appendix D.

Figure 8 shows the evolution trend of the annual average NDVI and the statistical averages of the five related climate factors in the metropolitan region from 2000 to 2023. In order to facilitate direct observation and comparison of their change patterns, we normalized the values of all variables so that their annual statistical values are within the range of 0-1. It is obvious that the change pattern of annual precipitation is most similar to that of the average NDVI. They increase and decrease almost in the same period. Even though there were several small surprises in 2001, 2009 and 2016, the overall trend is almost the same. Overall, before 2020, the annual NDVI and precipitation of BMR showed a slow and fluctuating growth trend, after which they fell rapidly and both fell to the lowest value in the past few years in 2023. Although the change patterns of several other temperature-related variables are not as close to those of NDVI as those of precipitation, their growth trends in the last few years have obviously accelerated under the warming climate, reaching a peak in 2022. In the specific research years we selected, namely 2006, 2012 and 2018, the change patterns of NDVI and precipitation were relatively close, and no unexpected situations occurred. Therefore, in the next step of analyzing the specific relationship between the two, we obtained a database with more reference value.

3.4. Spatial Analysis of Precipitation/NDVI Characteristic Triangle

3.4.1. Spatial Definition of Precipitation/NDVI Characteristic Triangle

Using scatter plots to analyze the correlation between variables, identify outliers and trend fitting is a relatively basic reference and research tool [44]. Prior to this, scholars have established scatter plots using NDVI and other various climate factors, confirming that there is a relatively obvious and stable characteristic triangular spatial feature between NDVI and them. In the early days, scholars performed linear fitting on the double edges of the scatter distribution based on the scatter distribution characteristics between LST and NDVI, and obtained the triangular spatial structure of LST/NDVI [45,46]. Later, it was proposed that there is a more stable single-hypotenuse triangle spatial feature between NDVI and night light intensity, which was verified in 5 cities around the world, and the physical characteristics of the three-sided structure of the triangle were used to successfully extract urban areas and classify cities in multiple cities [44]. By analyzing the scatter plot established by using the log function of Ln to process the precipitation and average NDVI in 2006, 2012 and 2018, we can see that there is also a similar double-hypotenuse triangle spatial structure feature between the two (Figure 9).

Inside the triangular space, each scatter point can be represented by precipitation and NDVI, which reflect the degree of humid or dry climate under the given condition of different green vegetation coverage. Among them, the upper hypotenuse represents the positive correlation between precipitation and NDVI, and the bottom hypotenuse reflects the negative relationship between precipitation and NDVI. The vertical edge on the right is the precipitation in the area with high vegetation coverage, and the intersection on the left is the lowest NDVI effective threshold. More meaningful is the vertical line from the vertex of the triangle to the right side, which can help us understand the mechanism of precipitation and NDVI. Above the vertical line is a humid area with abundant precipitation, which promotes the growth of vegetation, but under the relatively dry conditions below the vertical line, the slope of the equation line between vegetation cover and precipitation is negative, that is, the relationship between the two is opposite. Figure 10 more intuitively explains the composition and meaning of each part of the precipitation/NDVI triangle spatial feature.

3.4.2. Construction of the Triangular Spatial Structure of the Precipitation/NDVI Scatter Plot

Firstly, we used kernel density estimation to analyze the scatter point distribution density of three scatter plots of NDVI and precipitation in the metropolitan region in 2006, 2012 and 2018. The results are shown in Figure 11. We can find that the scatter point distribution is very uneven, and the density decreases from the place where the scatter points are highly concentrated to the surrounding area. The red dotted box encircles several relatively obvious low-density distribution scatter points. The presence of a large number of low-density scatter points at the edge of the triangle will definitely cause holes in the triangular spatial structure, which cannot reflect the precipitation and vegetation coverage of the metropolitan region, so we must remove these free scatter points.

Figure 12 shows the results of abnormal outlier identification of scattered point density. The blue bubble represents the LOF detection result of the scattered point. The larger the radius, the greater the possibility of being an outlier, and the red scattered point is the detected outlier that needs to be removed. This result is similar to the prediction in the previous step. A large number of free scattered points with low distribution density are selected, among which the highly concentrated scattered points in the core of the triangular structure are effectively retained. On this basis, we can construct the final precipitation/NDVI triangle structure.

The straight side of each characteristic triangle is obtained by the maximum value of the effectively retained NDVI, and the hypotenuse will be predicted using the segmented extreme value fitting method. We divide Ln (Mean_DNVI) into several small intervals with a width of 0.01. The upper and lower hypotenuses of the characteristic triangle should respectively meet the maximum and minimum precipitation intensity that can be provided under the NDVI conditions in each interval. After the precipitation extreme value and the corresponding NDVI of each interval are found, we use the least squares method for linear fitting. We use the maximum annual precipitation in 2006, the upper hypotenuse of the triangle, as an example to illustrate the hypotenuse fitting method (Figure 13).

The characteristic spatial triangular structure of the scattered point distribution of the three years and the fitting results of each side are finally presented in Figure 14 and Table 6.

3.4.3. Analysis of the Physical Meaning of the Triangular Space of the Precipitation/NDVI Scatter Plot

• Analysis of parameters of spatial structure of precipitation/NDVI characteristic triangle

Table 7 summarizes the parameters of all characteristic triangle structures, intuitively presenting the dynamic changes of the multiple interactions between precipitation and NDVI in these years. Between 2006 and 2018, the maximum precipitation in the metropolitan region first decreased slightly and then increased, with the final 14 mm higher than the initial 12 mm, while the minimum precipitation continued to decrease. The length of the straight edge indicates the threshold of the average annual precipitation. It can be seen that the thresholds in 2006 and 2012 are not much different, only decreasing by about 0.02, and then increasing to 0.48, which is more than twice that of 2012. The precipitation at the intersection represents the boundary of the relationship between precipitation and NDVI. Above this precipitation, there is a positive cooperative relationship between precipitation and NDVI, but if the precipitation is lower than it, there is a mutual inhibition between the two in the opposite direction. Obviously, this critical value continues to increase, from 10.05 in 2006 to 11.51 in 2018. This shows that as time goes by, precipitation must increase in order to significantly and effectively promote the growth of green vegetation on the surface.

The NDVI at the intersection is the minimum effective NDVI in the characteristic triangle space. The values of the three years are not much different, all around 0.14. The straight side of the triangle is the maximum NDVI, which is also the end point of the effective interaction between precipitation and NDVI. The maximum NDVI in 2006 and 2018 are not much different, about 0.83 and 0.85 respectively, and only the minimum in 2012 is around 0.79. The closer to the straight edge, the more scattered points there are, and the more complex the relationship between precipitation and NDVI is.

The strength of the relationship between precipitation and NDVI can be measured by the slope of the hypotenuse. In 2006, the slope of the wet edge (0.1) was much greater than that of the dry edge (0.03), indicating that the main effect of precipitation in that year was to promote the growth of NDVI. In contrast, the dominant relationship mechanism in 2012 was mutual inhibition, where higher precipitation did not effectively promote vegetation growth, and the distribution area of high precipitation on the ground covered by higher vegetation was not as large as the distribution area of low precipitation, and even the average precipitation here was not as good as that on the more exposed surface. The slope of the top hypotenuse in 2018 was 0.11, and the slope of the bottom hypotenuse was slightly larger, with 0.16. However, from the scatter distribution, there was a small blank in the bottom corner of the triangle in 2018, so from the overall consideration, the average precipitation in the area of high NDVI was greater than the model predicted. In the last year, the positive effect between the two may be more significant than the negative effect, or the two relationships are in a dynamic equilibrium mode.

The area of the triangle further confirms the strength and complexity of the two interaction relationships. We divide the precipitation/NDVI characteristic triangle into two triangles, the upper and lower ones, with the vertical line as the dividing line. The top one represents the reaction promotion, and the bottom one represents the reverse effect. As mentioned above, the coexistence pattern of precipitation and NDVI in 2006 was significantly stronger than the mutual interference, while the opposite was true in 2012. The area of the triangle in the dominant state in these two years was much larger than the other one. The areas of the two triangles in 2018 were similar. Considering the blank in the bottom triangle in 2018, the areas of the two triangles would be closer or even have a wider wet hypotenuse. In addition, the larger the area of the triangle, the more scattered points it contains, and the more complex the interaction mechanism between precipitation and NDVI. On the whole, the areas of the characteristic triangles in 2006 and 2012 were similar, with an almost negligible reduction, but in 2018 its area soared to 0.41, more than twice that of previous years. The interaction relationship in 2018 was the most complex.

• Practical analysis of spatial structure of precipitation/NDVI characteristic triangle

In order to find out the actual meaning of each edge of the characteristic triangle structure of precipitation/NDVI in the metropolitan region and the spatial distribution law of the scattered points near them, we selected the three edges of the characteristic triangle and the scattered points near the vertical line from the vertex to the straight edge as samples with a width of 0.01 units, restored the actual pixel positions of these scattered points on the BMR map, and drew Figure 15. Furthermore, we combined the physical meaning of each edge with the EQI results of each land cover type in the metropolitan region evaluated above, providing the possibility of summarizing the relationship between the precipitation/NDVI characteristic triangle and ecological quality.

In 2006, the top edge of the characteristic triangle, where the scatter points with high precipitation are mostly distributed in forests, cultivated land and grasslands in the north of the BMR, where the environmental quality is also very high. In addition, there are a small number of built-up land and urban land with poor ecology. There are also a few such pixels scattered in some places in the north of the central part of the BMR. The NDVI at the location of these pixels increases with the increase of precipitation. Regarding the pixels at the bottom edge, where they appear, the NDVI here is inversely proportional to the precipitation. In this year, they are more concentrated in cultivated land and forests in the south. They can also be seen in the central part of the BMR, especially in the BMA and Barcelona cities with a high degree of human development. Compared with the top edge, almost half of the scatter points at the bottom edge appear in the areas of urban land and built-up land. Where they appear, the NDVI here is inversely proportional to the precipitation. In general, the EQI of the pixels near the bottom edge of the triangle is not as high as that of the top edge.

On the contrary, in 2012, the locations of the pixels near the two hypotenuses of the triangle were reversed. The southern part of the metropolitan region became the main distribution area of the actual location of the wet edge scatter points, and there were also a few such pixels in the northern and central parts of the BMR. Even so, they still tended to appear in green areas with higher EQI. However, there were some exceptions in the BMA, which was covered by urban land and had a lower environmental quality index, but a considerable number of top pixels appeared. There were many dry edge pixels in the forested areas in the north of the metropolitan region, along the coastline, and in the central part and its north and south. These pixels not only preferred to be distributed in green areas with high EQI, but also in urban land and built-up land with poorer EQI, but overall, there was no significant difference in the ecological environment quality of the pixels near the top and bottom edges in this year. During this period, the precipitation in the southern part of the BMR was more abundant, and the vegetation cover was lush due to the wetter vegetation.

The top edge pixels in 2018 appeared in the northern and southern parts of the BMR, and the number of green areas such as forests and cultivated land was obviously greater than that of urban land and construction land. And such scattered points in the east and north of the BMA were also retained to a certain extent. The pixels representing drought appeared in the north-central, central and southern coastal areas of the BMR. Most of them were distributed on artificial surfaces with low ecological index, even though a few were concentrated in the forested areas at the north end.

On the other hand, the scattered points near the straight edge have very high NDVI, and the NDVI and precipitation they reflect are almost unaffected by each other. Their actual geographical location remains at the northern end of the BMR, with the main distribution area being forests, including a small amount of grassland and other green land, and has hardly changed in the three years. From the map, we can roughly understand that these pixels fall on the boundary between undeveloped forests and human activity land at the northern end of the metropolitan region. The pixels of the vertical line of the triangle represent different NDVIs under the same precipitation, which is the critical point of the two functional relationships between precipitation and NDVI. In the map, they are first scattered in the central and southern parts of the BMR, as well as various types of land in the very north. Then in 2012, their distribution area became wider and more scattered, but overall they were mostly in the central part. Finally, they extended to the northern metropolitan region boundary and coastline, forming a pattern of scattered distribution in the central and southern parts, as well as the northern jurisdiction boundary, and concentrated distribution along the northern coastline.

4. Discussion and Conclusions

In this study, we established a more reasonable dynamic ecological environment quality assessment system that is different from previous studies by scholars. Based on remote sensing technology, we used it to evaluate and demonstrate the ecological environment quality and its evolution process in the Barcelona metropolitan region from 2006 to 2018, and analyzed various external factors that may affect the change of land use with the highest ecological benefits. Finally, we studied and summarized the action mode of key climate factors that can indirectly control ecological quality. We found that the scale of green areas with the largest area and the most contribution to the ecological quality index in the metropolitan region is shrinking, but the surface is becoming more and more fragmented, and the spatial distribution structure is more complex. And the overall situation of the ecological environment quality of BMR is difficult to be satisfactory, and there has been no significant improvement in the past few years. Finally, we found that NDVI is the positive natural variable most closely related to the distribution of green space, and precipitation is the climate factor that can most strongly control the ecological environment. There is an extremely complex two-way influence mechanism between them. Generally speaking, the overall EQI of BMR is at a low level, the area of green space that is most beneficial to ecological quality is decreasing, and the degree of fragmentation is deepening. NDVI and precipitation are the two most critical factors that control environmental quality and influence each other in a complex way.

Principally, the overall results of the Ecological Environment Quality Index of BMR are impressive, but the distribution is very uneven, with lower levels covering larger areas than higher levels, especially in the core city of Barcelona, and most of the territory remained unchanged or deteriorated during the period studied, with only a few areas improving. Fortunately, the green areas with the highest ecological benefits are the main land use in BMR, even if they have slightly decreased in area and become more fragmented and complex in structure over the 12 years. Combined with the EQI results of green areas, fragmentation seems to be detrimental to the protection of the ecological environment. In this, the dominant one is forest, which is also the land with the highest EQI of all green areas, three times that of agricultural land or grassland. Urban green areas, perhaps due to their more scattered distribution in the urban layout, have the lowest ecological quality.

We modeled and analyzed 21 man-made and natural factors that may affect the distribution of green space, which plays the most critical role in controlling the EQI results. In addition to artificially built surfaces, which are undoubtedly the most fundamental man-made reasons affecting the layout of green space, NDVI, an uncontrollable natural factor, is the most significant positive explanatory variable associated with the distribution of green space. At the same time, climate change will significantly affect the NDVI value, thereby affecting the ecological quality. The climate factor that most strongly affects NDVI is precipitation. Figure 8 shows that the trend of annual NDVI and precipitation over 23 years is roughly similar, and they increase and decrease almost at the same time, compared with the other four climate variables. In fact, we believe that the relationship between changes in precipitation and NDVI is not as simple as the one-sided positive correlation shown in the model. It is consistent with the research results of most scholars. There is a highly complex and strong relationship between the two, which is difficult to explain by ordinary models [21,27,33,34,35,36].

Based on this, we tried to use the precipitation-NDVI characteristic triangle space to successfully explain this mechanism. The precipitation at the intersection of the two hypotenuses of the characteristic triangle is the critical value, which operates two interactions with NDVI. Above this precipitation, the two are cooperative and mutually promoted, and below this, there is a counter-action relationship between them. And this critical value is increasing every year. The slopes of the two hypotenuses and the size of the area of the two triangles divided by the vertical line both imply the strength of the two relationships. The higher the slope and the larger the area, the greater the intensity of the interaction represented, and the greater the impact of this on the NDVI of that year. As predicted in the article, the dominant position in 2012 is the reverse relationship. NDVI is suppressed by the influence of precipitation, and the related EQI will also decrease. We link it with Table 3 and find that the EQI ratio of green land in this year is the lowest among the three years. In addition, the threshold and total area of the straight side of the characteristic triangle also show that the interaction between BMR precipitation and NDVI is becoming more and more complex. Finally, we analyzed the geographical distribution of pixels corresponding to the scattered points near the three sides and vertical lines of each characteristic triangle and summarized their typical land covers and their EQIs.

According to the results of this study, we can highly confirm that precipitation, which is the dividing value between the two effects of precipitation and NDVI, is increasing, which means that the bottom hypotenuse of the characteristic triangle, that is, the possibility of the reverse effect between the two taking the upper hand, is increasing. In the context of global warming, according to the latest study by Arellano, climate change is leading to a decrease in rainfall, and they firmly believe that the climate in Spain will be significantly drier and warmer in the next 30 years [47]. Based on this, we assume that in the future, the relationship between precipitation and NDVI in BMR may be greatly affected by this, the dividing value between the two effects will continue to increase, the dominance of the reverse relationship will become increasingly dominant, and EQI will also be negatively affected to a certain extent. This will be one of the main goals of our further work.

In this article, a dynamic ecological environment quality assessment system that can be applied to different time and region changes is established through NDVI and NDBI, and the impact of climate change, human activities and other natural factors on the evolution of environmental quality is comprehensively considered, and the most critical direct and indirect influencing variables that control the EQI results are found. Finally, the characteristic triangle spatial model is established to analyze precipitation, an important climate factor that is highly complex and closely related to environmental changes, and to a certain extent reveals the mechanism of action between the two. This study cuts into the perspective of two variables that are crucial to EQI from the perspective of global environmental quality assessment, and deeply analyzes the impact of various factors on environmental quality. However, due to limited access to technology and data and the particularity of precipitation and extreme temperature data, their differences are too large, so we adopted a re-interpolation method, which will affect the accuracy of the experimental results. In future research, we will try to use or establish a more accurate database ourselves.