Multi-Temporal Remote Sensing Inversion of Evapotranspiration in the Lower Yangtze River Based on Landsat 8 Remote Sensing Data and Analysis of Driving Factors

: Analysis of the spatial and temporal variation patterns of surface evapotranspiration is important for understanding global climate change, promoting scientiﬁc deployment of regional water resources, and improving crop yield and water productivity. Based on Landsat 8 OIL_TIRS data and remote sensing image data of the lower Yangtze River urban cluster for the same period of 2016–2021, combined with soil and meteorological data of the study area, this paper constructed a multiple linear regression (MLR) model and an extreme learning machine (ELM) inversion model with evapotranspiration as the target and, based on the model inversion, quantitatively and qualitatively analyzed the spatial and temporal variability in surface evapotranspiration in the study area in the past six years. The results show that both models based on feature factors and spectral indices obtained a good inversion accuracy, with the fusion of feature factors effectively improving the inversion ability of the model for ET. The best model for ET in 2016, 2017, and 2021 was MLR, with an R 2 greater than 0.8; the best model for ET in 2018–2019 was ELM, with an R 2 of 0.83 and 0.62, respectively. The inter-annual ET in the study area showed a “double-peak” dynamic variation, with peaks in 2018 and 2020; the intra-annual ET showed a single-peak cycle, with peaks in July–August. Seasonal differences were obvious, and spatially high-ET areas were mainly found in rural areas north of the Yangtze River and central and western China where agricultural land is concentrated. The net solar radiation, soil heat ﬂux, soil temperature and humidity, and fractional vegetation cover all had signiﬁcant positive effects on ET, with correlation coefﬁcients ranging from 0.39 to 0.94. This study can provide methodological and scientiﬁc support for the quantitative and qualitative estimation of regional ET.


Introduction
Evapotranspiration (ET) is the process by which surface water is transferred to the atmosphere, including evaporation of water vapor captured from water bodies, soils, and vegetation surfaces, and transpiration by plants [1]. Evapotranspiration is a major component of the global water cycle and is a key link between terrestrial water, carbon, and surface energy exchange, with strong links to agriculture, hydrology, meteorology, and ecology [2]. Accurate estimates of evapotranspiration are important for understanding global climate change, ecological issues, the water cycle, and hydrological processes, as well as for improving water use efficiency, crop yields, and water productivity [3][4][5].
The conventional estimation method for ET is to analyze meteorological data through the information obtained from meteorological stations and ground flux monitoring stations and 10 m. Reyadh Albarakat [25] and Lida Andalibi et al. [26] found that although the inversion results of the three sensors were close in the same area, the lower spatial resolution and wider sensor bandwidth of AVHRR and MODIS resulted in more significant water vapor absorption and showed noisy behavior in the dataset. For the Sentinel-2 satellite, although it has a spatial resolution of 10 m, its large data volume and limited image storage services make it difficult to perform multi-year time series monitoring. In addition, the 10 m resolution results in a huge amount of computation and requires a high computing power, which is not conducive to application promotion. In contrast, Landsat technology is more mature and easier to apply, and it has considerable typicality and ubiquity, which is suitable for promotion and application.
Vegetation cover is an important factor influencing the amount of surface evapotranspiration. With the extensive use of remotely sensed vegetation indices, it is possible to monitor changes in vegetation cover on a large scale [27]. The normalized difference vegetation index (NDVI), which compares the intensity of reflectance in the visible red and near-infrared bands, is the most commonly used VI to quantify the presence of live green vegetation. By relying on a ratio of band intensities, NDVI removes a large proportion of noise caused by cloud shadows, topographic and solar angle variations, and atmospheric attenuations existing in the visible red and infrared bands, which makes NDVI less susceptible to illumination conditions [27][28][29].
The current inversion of ET using remote sensing mainly focuses on the evaluation of the accuracy of different models or the mechanism, but there are fewer studies on the analysis of ET changes and driving factors in the same area under different years (multitemporal phases). Therefore, this paper took the lower Yangtze River urban agglomeration (Nanjing, Changzhou, and Zhenjiang, hereafter referred to as the three cities) as the study area, and used the Landsat 8 OIL_TIRS series to acquire multispectral remote sensing images of the study area for the same period from 2016 to 2021. In addition, the Food and Agriculture Organization of the United Nations (FAO) Harmonized World Soil Database v 1.2 platform was used to obtain ET, soil temperature and humidity, solar short-wave radiation, and heat flux data of the study area for the same period. On the basis of the reduced dimensional analysis of the spectral data, a multiple linear regression (MLR) model and an extreme learning machine (ELM) inversion model were established, and the spatial distribution of ET and the characteristic factors driving it in different time periods in 2016-2021 and 2017 were analyzed on the basis of the model inversion. This study provides a theoretical basis for modeling the spatial distribution of ET evapotranspiration over different time periods.

Study Area
The three cities (N 31 • 9 -32 • 36 , E 118 • 21 -120 • 11 ) are located in the southwest of Jiangsu Province, and the monitoring stations are laid out as shown in Figure 1. With a total area of 14,785 km 2 and a population of 16,677,400, the three cities have a gently sloping terrain with an average altitude of 50 m above sea level and are located in a subtropical monsoon climate zone with abundant rainfall of 1200 mm per year and an average annual temperature of 15.4 • C. Due to its location in the lower reaches of the Yangtze River and the Huai River, the total length of the coastline along the river is nearly 360 km, and a large alluvial plain makes up the rich arable land resources, with an arable land area of 406,000 hectares, and an annual grain production of over 35 million tons, making it one of the 13 major grain-producing areas in China.

Satellite Data
The remote sensing data selected for this study were the Landsat 8 OIL_TIRS data (https://www.gscloud.cn/, accessed on 17 February 2023), developed by NASA in collaboration with the United States Geological Survey (USGS). The satellite data are for the years 2016-2021 (the year 2020 is not included in this paper due to the poor results of the 2020 satellite data), and information about the data is shown in Table 1.
The remote sensing data were pre-processed with ENVI5.3 for image fusion (fusion of multispectral bands with panchromatic bands to improve image resolution), radiometric correction, and atmospheric correction. Based on the GPS coordinates of the monitoring points, the average reflectance of the samples in the monitoring area was extracted as the spectral reflectance of the sampling point by plotting the region of interest (ROI), and the spectral reflectance data of different wavelengths were obtained. The remote sensing data selected for this study were the Landsat 8 OIL_TIRS data (https://www.gscloud.cn/, accessed on 17 February 2023), developed by NASA in collaboration with the United States Geological Survey (USGS). The satellite data are for the years 2016-2021 (the year 2020 is not included in this paper due to the poor results of the 2020 satellite data), and information about the data is shown in Table 1.
The remote sensing data were pre-processed with ENVI5.3 for image fusion (fusion of multispectral bands with panchromatic bands to improve image resolution), radiometric correction, and atmospheric correction. Based on the GPS coordinates of the monitoring points, the average reflectance of the samples in the monitoring area was extracted as the spectral reflectance of the sampling point by plotting the region of interest (ROI), and the spectral reflectance data of different wavelengths were obtained. The meteorological information for the study area was obtained from the China Meteorological Data Network (https://data.cma.cn/, accessed on 17 February 2023). The monitoring station layout is consistent with the actual ET measurement points. Zhang et al. [30] pointed out that the soil temperature and surface temperature (LST) are closely related and their temporal variability tends to show similar characteristics. However, this similarity is largely dependent on the influence of meteorological factors and the physical properties of the subsurface (soil heat flux, vegetation cover, snow cover, etc.); therefore, LST data do not truly reflect the true soil temperature, and thus cannot produce more accurate predictions and simulations of soil ET. In summary, the meteorological variables recorded included the air temperature, relative humidity, net solar radiation, and soil temperature and humidity.

Measured ET Data
The actual ET data for this study were obtained from 21 ET monitoring stations in the study area ( Figure 1), with the stations evenly distributed according to the latitude and longitude grid. The data collected in this study included the daily actual ET, monthly actual ET, and annual actual ET at each site in 2016-2021, where the time of collection was aligned with the time of the remote sensing images.

Vegetation Index Construction
Vegetation indices are spectral indices made up of different combinations of wavebands. They are a simple and effective empirical measurement method for different features and their characteristic information, and are widely used in remote sensing monitoring of soil parameters, crop growth, environmental changes, etc. [31][32][33][34]. In this study, 14 vegetation indices were selected as indicators for constructing the ET inversion model of the study area according to the requirements of the evapotranspiration inversion, as shown in Table 2.  [46] Vegetation indices are linear or non-linear combinations of different bands in remote sensing data, which result in a high degree of linear correlation between vegetation counts. In multiple linear regression and multiple regression models, strong correlations between the independent variables can lead to unstable estimates, which can affect the accuracy of the model and its later replication [47]. The correlation between the independent variables can be measured using the variance inflation factor (VIF), which is generally considered to be strong when VIF > 10. The principle of VIF is as follows: where R j represents the negative correlation coefficient of the independent variable on the remaining independent variables for regression analysis, and p is the number of all independent variables.

Cluster Analysis
In the actual modeling process, as most of the vegetation indices reflect the vegetation physiology, soil water content, etc., some of the vegetation indices have the same or a similar meaning for ET interpretation, which may lead to a more redundant model. Therefore, under the condition that the information contained in the vegetation indices is not affected, some of the spectral indices can be clustered and analyzed, and a reasonable number of classifications can be selected for data dimensionality reduction according to the results of the clustering analysis.

Factor Analysis
Factor analysis (FN) is a generalization and development of principal component analysis. FN combines objects (variables or samples) with complex relationships into a small number of factors, retaining the intrinsic links between the factors and the original variables while diluting the linear relationships between the factor variables. The formula for the principle of factor analysis is as follows: where a ij represents the factor loading; A = (a ij ) is the factor loading matrix; ε is the special factor. The factor score function is derived from the factor analysis model as follows: F j = a j1 X 1 + a j2 X 2 + . . . + a j p X p j = 1, 2, . . . , m

Ecosystem Parameters
ET is the total flux of water vapor transported to the atmosphere by vegetation and the ground and consists mainly of vegetation transpiration and soil water evaporation. The intensity of soil water evaporation is mainly influenced by the soil water content, soil temperature, and soil structure. In 1998, the Food and Agriculture Organization of the United Nations (FAO) proposed Penman-Monteith as the only standard for calculating the evaporation of vegetation (ET0). ET0 is calculated as follows: where ET0 is the reference evapotranspiration, mm/day; R n is the net radiation on the crop surface, MJ/(m 2 ·day); G is the soil heat flux, MJ/(m 2 ·day); T is the average daily temperature at an altitude of 2 m, • C; u 2 is the wind speed at a height of 2 m, m/s; e s is the saturated water vapor pressure, kPa; e a is the actual water vapor pressure, kPa; e s − e a is the saturated vapor pressure difference, kPa; ∆ is the slope of the saturated water vapor pressure curve; r is the hygrometer constant, kPa/ • C. In summary, the variables of soil heat flux (Qs), soil temperature (Ts) and moisture (Th), and net solar radiation (D) were introduced as predictive variables to the input model, taking into account the influences of vegetation transpiration and soil evaporation.

Fractional Vegetation Cover
Fractional vegetation cover (FVC) is the percentage of vegetation (including leaves, stems, and branches) per unit area of the vertical projection of vegetation. FVC directly affects soil moisture retention capacity and regional plant water content, and thus surface evapotranspiration. In this study, FVC was calculated based on the two-dimensional pixel model (DPM) with the following equation: where NDVI soil is the NDVI value of bare land, and NDVI veg is the NDVI value of the area covered by vegetation.
2.6. Model and Accuracy Evaluation 2.6.1. ELM ELM is a new learning algorithm for single hidden layer feedforward neural networks (SLFNs) proposed by Huang et al. [48,49]. The essence of ELM is that the learning parameters of the hidden layer are randomly assigned and do not need to be adjusted, while the output weights can be determined via simple generalized inverse matrix operations [50]. Compared to other traditional SLFNs, ELM has a faster learning speed, better generalization performance, and minimal human intervention.
The principle of ELM is shown in Figure 2. The above structure can be expressed by the formula Equation (8) can be rewritten as where T is the full output; L is the number of hidden layers; g(x) are the activation functions; n is the number of input layers; w and β are the connection weights between the layers; H is the implied layer output matrix.

Multiple Linear Regression
Multiple linear regression analysis is a multivariate statistical technique used to examine the relationship between a single dependent variable and a set of independent variables [51]. MLR forecasting models are expressed in the following format: where Y t is the predicted value at time t, X t = (1, x 1t , x 2t , . . . , x kt ) is a vector of k explanatory variables at time t, β = (β0, β1, . . . , βk) is the vector of coefficients, and ε t is a random error term at time t. The errors terms should be independent and have a Gaussian distribution.

Accuracy Evaluation
In order to obtain the optimal model, the coefficient of determination (R 2 ), root mean square error (RMSE), and rate of performance deviation (RPD) were used to evaluate the model and discuss the applicability of each model. R 2 reflects the degree of fit between the predicted and measured values, and the closer its value is to 1, the better the model fits. The RMSE reflects the deviation between the predicted and measured values, and the smaller the value, the better the prediction effect of the model. The principles of each evaluation index are as follows:

Statistical Analysis of Measured ET
The measured ET data at each monitoring site from 2016 to 2021 were analyzed by year for descriptive statistics, and the sample data distribution is shown in Table 3. The data from each year were randomly partitioned using the trian_test_split dataset from the Python Scikit Learn library for the model training set and the test set, i.e., 90% of the samples were used for training and 10% were used for validation.

Variable Correlation Analysis
The correlation analysis of 21 spectral indices (7 band reflectance and 14 vegetation indices) and the measured ET in different years is shown in Table 4.

Data Dimensionality Reduction
The VIF between the spectral indices in this experiment is shown in Figure 3. The result showed that there was a serious problem of multicollinearity between the spectral indices.

Data Dimensionality Reduction
The VIF between the spectral indices in this experiment is shown in Figure 3. The result showed that there was a serious problem of multicollinearity between the spectral indices.
According to Figure 3, it can be seen that there was strong multicollinearity between the spectral indices. In this paper, factor analysis and cluster analysis were used to extract potential collinearity factors to solve the collinearity problem.
Firstly, the spectrum indicators with severe covariance were clustered, and the output of the hierarchical cluster analysis pulse diagram after using the origin is shown in Figure 4.  According to Figure 3, it can be seen that there was strong multicollinearity between the spectral indices. In this paper, factor analysis and cluster analysis were used to extract potential collinearity factors to solve the collinearity problem.
Firstly, the spectrum indicators with severe covariance were clustered, and the output of the hierarchical cluster analysis pulse diagram after using the origin is shown in Figure 4. As can be seen in Figure 4, the spectral indicators were divided into four categories for different years, with band reflectance B1 to B7 grouped together in different years as the different bands are still essentially electromagnetic spectra with similar interpretations As can be seen in Figure 4, the spectral indicators were divided into four categories for different years, with band reflectance B1 to B7 grouped together in different years as the different bands are still essentially electromagnetic spectra with similar interpretations of ground information.
Factor analysis was performed on the classes with a sample size greater than 1 in the clustering results. With a KMO of 0.6 to 0.8 and a p-value of 0.000 *** for the significance of Bartlett's spherical test, the factor analysis was valid. The original spectral indices of the different classes were substituted into the factor analysis score function, and the results obtained were recombined with the spectral indices of the individual classes to form new spectral indices, denoted as F1, F2, F3, and F4. The results of the collinearity test of the new spectral indices are shown in Table 5. By comparing the data and analysis in Figure 4 and Table 4, it can be seen that the VIF of the independent variables was reduced to below 10 after the dimensionality reduction of the cluster factor analysis, and the collinearity between the variables was basically resolved, which improved the stability of the model.

Results of Model Inversions
The four new spectral indices and four characteristic factors affecting ET changes through factor analysis were selected as independent variables, and the measured ET  Table 6. Only spectral parameters were selected to construct the model results, as shown in Table 7. From Table 6, it can be concluded that the top three ranked models in terms of accuracy were 2021-MLR, 2017-ELM, and 2017-MLR, with an R v 2 of 0.87, 0.86, and 0.85, respectively, after combining the three evaluation indicators (R 2 , RMSE, and RPD); the lowest accuracy was found for 2019-MLR, with an R v 2 of 0.59, which showed a more obvious decrease compared with the models of the other years, because the cloud coverage of the Landsat satellite images was larger in 2019 and there were errors in the surface reflectance extracted via remote sensing.
Comparing the results in Table 6 with those in Table 7, it was found that when only spectral parameters were introduced as independent variables, the highest R v 2 was found for 2016-ELM and 2018-MLR, which was only 0.57, much lower than the model accuracy when ET feature parameters were fused under the same model; the reason for this is that, when only spectral parameters were used as independent variables for inversion, the mechanism of interaction between ET and each spectral parameter was unclear, and it was more difficult to perform the inversion through mathematical models. The characteristic parameters, as triggers of ET, are themselves highly correlated with ET, and their introduction as independent variables can effectively explain the correlation between the dependent and independent variables [20].
In order to verify the accuracy of the two models in different years, the predictedmeasured scatter plots of the ELM and MLR validation sets were plotted using the measured ET values of the three cities from 2016 to 2021, as shown in Figure 5. There was a good linear relationship between the estimated and measured values of the two models (R 2 > 0.44, RMSE < 0.08 mm), and the sample points on both sides of the regression line were evenly distributed. The MLR model had the best estimation in 2016, 2017, and 2021, with an R 2 ranging from 0.76 to 0.82, but the estimation in 2018 was relatively poor, with an R 2 of 0.44. Meanwhile, the ELM model had higher accuracy (R 2 > 0.73, RMSE < 0.02 mm) in different years, with a more concentrated sample distribution.

Time-Series Change Characteristics
As shown in Figure 6

Inter-Annual Spatial Distribution Characteristics
Based on the ET remote sensing estimation model constructed in 3.4, the year-byyear data of the study sites were imported into the model. The inversion of the spatial distribution of inter-annual changes from 2016 to 2021 is shown in Figure 8. As can be seen in Figure 8, the areas with higher ET were mainly concentrated in the rural areas north of the Yangtze River and the central and western plains where agricultural land is concentrated, while the urban areas in the study area had lower evapotranspiration due to having less bare soil and cultivated land area and a lower vegetation cover. During the study years, ET ranged from 782 to 979 mm in 2016; in 2018, evapotranspiration increased significantly, with 87% of the study area having inter-annual ET greater than the multiyear average evapotranspiration; in 2019, ET in the study area began to decline slowly, with 51.82% of the total study area having inter-annual ET above the multi-year average The variation in evapotranspiration by month in the study area from 2016 to 2021 is shown in Figure 7. Evapotranspiration rose continuously from January, fluctuated slightly in May, continued to rise to a peak in July-August, and then began to decline, with the minimum value during the year usually occurring in December, with evapotranspiration mainly concentrated in April-August during the year.

Inter-Annual Spatial Distribution Characteristics
Based on the ET remote sensing estimation model constructed in 3.4, the year-byyear data of the study sites were imported into the model. The inversion of the spatial distribution of inter-annual changes from 2016 to 2021 is shown in Figure 8. As can be seen in Figure 8, the areas with higher ET were mainly concentrated in the rural areas north of the Yangtze River and the central and western plains where agricultural land is concentrated, while the urban areas in the study area had lower evapotranspiration due to having less bare soil and cultivated land area and a lower vegetation cover. During the study years, ET ranged from 782 to 979 mm in 2016; in 2018, evapotranspiration increased significantly, with 87% of the study area having inter-annual ET greater than the multiyear average evapotranspiration; in 2019, ET in the study area began to decline slowly,

Inter-Annual Spatial Distribution Characteristics
Based on the ET remote sensing estimation model constructed in Section 3.4, the year-by-year data of the study sites were imported into the model. The inversion of the spatial distribution of inter-annual changes from 2016 to 2021 is shown in Figure 8. As can be seen in Figure 8, the areas with higher ET were mainly concentrated in the rural areas north of the Yangtze River and the central and western plains where agricultural land is concentrated, while the urban areas in the study area had lower evapotranspiration due to having less bare soil and cultivated land area and a lower vegetation cover. During the study years, ET ranged from 782 to 979 mm in 2016; in 2018, evapotranspiration increased significantly, with 87% of the study area having inter-annual ET greater than the multi-year average evapotranspiration; in 2019, ET in the study area began to decline slowly, with 51.82% of the total study area having inter-annual ET above the multi-year average evapotranspiration; in 2021, ET decreased significantly, with an average interannual evapotranspiration of only 845 mm. This is a decrease of 24 mm compared to the multi-year average evapotranspiration.

Seasonal Spatial Distribution Characteristics
Taking 2017 as an example, the intra-annual variation and distribution characteristics of ET based on month-by-month data are shown in Figure 9. The spatial distribution of evapotranspiration in all four seasons was basically consistent with the spatial distribution throughout the year, with high-value areas distributed in areas with high vegetation cover and obvious farmland features, especially in the north and central-west, and lower evapotranspiration along the Yangtze River and near urban clusters. ET changes in the study area all showed obvious seasonal characteristics, with evapotranspiration in February ranging from 14.65 to 46.65 mm, with an average evapotranspiration of 36 mm, and high-evapotranspiration areas mainly concentrated in forests and spring-wheat-growing areas in the central and western parts of the study area. Evapotranspiration in May (summer) was significantly higher, and compared to February, evapotranspiration throughout the study area was increased by more than 10 mm, except for some water bodies that were not included in the study. The evapotranspiration was slightly greater in October (autumn) than in February, with an overall change of 6-18 mm, accounting for 84.83% of the total study area. December (winter) had the least evapotranspiration, with a significant decrease compared to spring, with 90.65% of the study area reduced by 1.6-3.2 mm.

Seasonal Spatial Distribution Characteristics
Taking 2017 as an example, the intra-annual variation and distribution characteristics of ET based on month-by-month data are shown in Figure 9. The spatial distribution of evapotranspiration in all four seasons was basically consistent with the spatial distribution throughout the year, with high-value areas distributed in areas with high vegetation cover and obvious farmland features, especially in the north and central-west, and lower evapotranspiration along the Yangtze River and near urban clusters. ET changes in the study area all showed obvious seasonal characteristics, with evapotranspiration in February ranging from 14.65 to 46.65 mm, with an average evapotranspiration of 36 mm, and highevapotranspiration areas mainly concentrated in forests and spring-wheat-growing areas in the central and western parts of the study area. Evapotranspiration in May (summer) was significantly higher, and compared to February, evapotranspiration throughout the study area was increased by more than 10 mm, except for some water bodies that were not included in the study. The evapotranspiration was slightly greater in October (autumn) than in February, with an overall change of 6-18 mm, accounting for 84.83% of the total study area. December (winter) had the least evapotranspiration, with a significant decrease compared to spring, with 90.65% of the study area reduced by 1.6-3.2 mm.

Environmental Factors
As the most direct receptors of climate change, ecological and environmental factors play an irreplaceable role in all aspects of the ecosystem. To explore the influence of meteorological factors on regional ET, the characteristics of the changes in the mean solar net radiation intensity, soil heat flux, and soil temperature and humidity in the study area from 2016 to 2021 were statistically analyzed, as shown in Figure 10. From Figure 10, it can be seen that, except for the fluctuation in soil moisture, the net solar radiation intensity, soil heat flux, and soil temperature and humidity all showed strong seasonal characteristics, with an increasing trend from January to June and a decreasing trend from August to December, with peaks occurring in July and August, which were basically consistent with the changes in ET during the year. To explore the correlation between the above four environmental factors and the dynamic changes in ET in the study area, the correlation between ET and each of the four environmental factors was calculated for each month from 2016 to 2021, as shown in Figure 11. From Figure 11, it can be seen that the net solar radiation, soil temperature, and soil heat flux had highly significant effects on ET in the study area, with correlation coefficients of 0.94, 0.91, and 0.70 (p < 0.01), respectively. In addition, soil moisture also had a significant effect on ET in the study area, with a correlation coefficient of 0.39 (p < 0.05).

Environmental Factors
As the most direct receptors of climate change, ecological and environmental factors play an irreplaceable role in all aspects of the ecosystem. To explore the influence of meteorological factors on regional ET, the characteristics of the changes in the mean solar net radiation intensity, soil heat flux, and soil temperature and humidity in the study area from 2016 to 2021 were statistically analyzed, as shown in Figure 10. From Figure 10, it can be seen that, except for the fluctuation in soil moisture, the net solar radiation intensity, soil heat flux, and soil temperature and humidity all showed strong seasonal characteristics, with an increasing trend from January to June and a decreasing trend from August to December, with peaks occurring in July and August, which were basically consistent with the changes in ET during the year. To explore the correlation between the above four environmental factors and the dynamic changes in ET in the study area, the correlation between ET and each of the four environmental factors was calculated for each month from 2016 to 2021, as shown in Figure 11. From Figure 11, it can be seen that the net solar radiation, soil temperature, and soil heat flux had highly significant effects on ET in the study area, with correlation coefficients of 0.94, 0.91, and 0.70 (p < 0.01), respectively. In addition, soil moisture also had a significant effect on ET in the study area, with a correlation coefficient of 0.39 (p < 0.05). Remote Sens. 2023, 15, x FOR PEER REVIEW 18 of 25

Factors of FVC
FVC directly affects the soil water holding capacity and the water content within regional plants, thus affecting the amount of surface evapotranspiration [52]. The FVC and ET change curves for the study area from 2016 to 2021 are shown in Figure 12. As can be seen in the figure, FVC had a positive effect on surface evapotranspiration, and both curves are periodic single-peaked curves. From March to August, with the growth of plants, the water in plants increased and the density of FVC increased, providing moisture conditions for surface evapotranspiration. ET and FVC showed a significant upward trend and reached a peak in August-September, and from September to December, due to weather changes, plants started to wilt, and thus FVC and ET declined from September to December, which continued until February of the following year.

Factors of FVC
FVC directly affects the soil water holding capacity and the water content within regional plants, thus affecting the amount of surface evapotranspiration [52]. The FVC and ET change curves for the study area from 2016 to 2021 are shown in Figure 12. As can be seen in the figure, FVC had a positive effect on surface evapotranspiration, and both curves are periodic single-peaked curves. From March to August, with the growth of plants, the water in plants increased and the density of FVC increased, providing moisture conditions for surface evapotranspiration. ET and FVC showed a significant upward trend and reached a peak in August-September, and from September to December, due to weather changes, plants started to wilt, and thus FVC and ET declined from September to December, which continued until February of the following year. The distribution of FVC in the three cities from 2016 to 2021 is shown in Figure 13. Its spatial distribution was basically consistent with the spatial distribution of ET in the same period. The northern, central, and southwestern areas of the study area are mostly The distribution of FVC in the three cities from 2016 to 2021 is shown in Figure 13. Its spatial distribution was basically consistent with the spatial distribution of ET in the same period. The northern, central, and southwestern areas of the study area are mostly cultivated land and natural forest land, where the spatial distribution of vegetation is dense and the soil water content is high, meaning the transpiration of plants and soil evaporation are higher than those in other areas. The spatial distribution of the correlation coefficients Remote Sens. 2023, 15, 2887 19 of 24 between ET and FVC is shown in Figure 14. This is mainly due to the fact that most of the positive-correlation areas are cultivated and forested, and bare soil and vegetation variations have a greater influence on ET.
Remote Sens. 2023, 15, x FOR PEER REVIEW 20 of 25 cultivated land and natural forest land, where the spatial distribution of vegetation is dense and the soil water content is high, meaning the transpiration of plants and soil evaporation are higher than those in other areas. The spatial distribution of the correlation coefficients between ET and FVC is shown in Figure 14. This is mainly due to the fact that most of the positive-correlation areas are cultivated and forested, and bare soil and vegetation variations have a greater influence on ET.  cultivated land and natural forest land, where the spatial distribution of vegetation is dense and the soil water content is high, meaning the transpiration of plants and soil evaporation are higher than those in other areas. The spatial distribution of the correlation coefficients between ET and FVC is shown in Figure 14. This is mainly due to the fact that most of the positive-correlation areas are cultivated and forested, and bare soil and vegetation variations have a greater influence on ET.

Discussion
In this paper, using Landsat 8 remote sensing satellite images of the lower Yangtze River urban agglomeration from 2016 to 2021, we introduced five ecological and environmental factors involved in the Penman-Monteith equation: soil heat flux, soil temperature and humidity, net solar radiation, and vegetation cover, and established a multi-temporal ET prediction model based on ELM and MLR, with refined inversions of both models. The results were good, with an R2 above 0.59. In addition, the spatial and temporal variation characteristics of ET in the study area were analyzed through remote sensing visualization, and the driving forces of soil heat flux and other factors were analyzed and compared. The correlation coefficients showed that the net solar radiation, soil heat flux, soil temperature, and vegetation cover all significantly influenced the ET variation and greatly enhanced the accuracy of the prediction models. This finding indicates that the combination of spectral data and ecological data can be used to carry out the prediction and accuracy optimization of regional remote sensing ET.

Data Optimization Studies
For predictive models, correlations between predictors can have unpredictable and inconsistent effects on parameter estimates and significance, and may lead to biased results [53]. In most remote sensing inversion studies, researchers tend to ignore this issue or simply eliminate variables with strong covariance, which may result in the elimination of variables that contribute more to the model prediction at the expense of the accuracy of the model prediction. In this study, clustering analysis was used to group variables with the same information and perform factor analysis, which can effectively deal with the problem of multiple cointegration between variables and improve the model prediction accuracy and stability. The results showed that the VIF of the variables treated using the cluster factor analysis dropped to below 3, and the models could still maintain a high R 2 .

Model Accuracy Discussion
Many studies have shown that the use of machine learning models can show better predictive power in estimating ET [20]. Compared to the traditional PM equation, the introduction of spectral data, meteorological, and ecological data in this study yielded better predictors (R 2 = 0.59−0.87, RMSE < 0.59) [54]. In addition, when comparing the machine learning model with the linear regression model in this study horizontally, the difference in inversion accuracy between the two was not significant, and even in some year types, the MLR prediction accuracy was higher than that of the ELM model. This is because machine learning techniques rely on a large amount of data to achieve high performance, and when the data samples are the same and have not reached their peak, the machine learning prediction effect does not prevail. In contrast, the traditional linear regression achieved good prediction results based on the theoretical causality of the data [55], while some of the induced factors introduced based on the PM equation, which itself has a good linear relationship with ET, greatly increased the explanatory power of the multiple linear regression for ET.

Selection of Vegetation Indices
NDVI, a widely used vegetation index, plays an important role in monitoring vegetation. However, environmental factors such as atmospheric conditions, topographical effects, and topographic illumination may affect it significantly [56]. Glenn et al. [57] found that the enhanced vegetation index (EVI), compared to NDVI, improved vegetation monitoring through a de-coupling of the canopy background signal and a reduction in atmospheric influences. However, the presence of soil regulators in EVI does not eliminate the influence of topographic effects and has limitations for estimating vegetation cover changes on a large scale. In addition, Moreira et al. [58] found that topographic illumination was negligible in areas where the average watershed slope was below 25 • and the latitude was below 45 • . From a review of the data, the study area (lower Yangtze River plain) has a relatively flat topography, with an average elevation of 20-30 m and an average slope of 2-6 • ; therefore, the choice of NDVI as a parameter for vegetation cover estimation has some applicability.

Spatial-and Temporal-Change-Driven Analysis
Spatially, the distribution of ET in the study area was basically in line with the general trend, mainly concentrated in areas with high vegetation cover such as farmland and forests in the urban periphery and extensive agricultural land. From a temporal perspective, ET showed more obvious seasonal variations within the year, mainly in summer > autumn > spring > winter. According to the analysis of the previous driving factors, ET changes are influenced by the net solar radiation intensity, soil heat flux, and soil temperature and humidity, and these four driving factors themselves are more sensitive to the season, thus leading to seasonal differences in ET. According to the mean ET values measured at each meteorological station, the maximum ET during the year was 912.67 mm in 2018, and the minimum was 845.02 mm in 2021, which is consistent with the results of this paper.

Conclusions
This paper took the lower Yangtze River urban agglomeration as the research object, combined Landsat 8 OIL_TIRS remote sensing data and five types of ecological environment data derived from the Penman-Monteith formula, constructed ET inversion models under different years and seasons based on the ELM and MLR models, and analyzed the distribution characteristics and influence factors of ET. The following conclusions were drawn:

1.
It is feasible to build an ET inversion model using satellite remote sensing spectral data and ecological environment data. Both models could invert the ET in the study area with high accuracy, and the introduction of ecological environment factors significantly improved the model accuracy. Among them, the best model for inversion from 2016 to 2019 was the ELM model, with an R 2 of 0.82, 0.86, 0.83, and 0.62, and an RMSE of 0.31 mm, 0.06 mm, 0.59 mm, and 0.29 mm, respectively; the best model for inversion in 2021 was the MLR model, with an R 2 of 0.87 and an RMSE of 0.29 mm.

2.
The spatial and temporal distributions of ET in the study area showed typical geographical and seasonal characteristics. Spatially, higher-ET areas were concentrated in the north, central, and south-western areas with higher vegetation cover and more extensive agricultural land; ET values were lower in areas close to cities. Temporally, ET peaked in summer, and the values of ET in all seasons, from high to low, followed the order summer > autumn > spring > winter.

3.
Net solar radiation, soil heat flux, soil moisture, and soil temperature significantly affected the ET variation, with correlation coefficients of 0.94 (p < 0.001), 0.70 (p < 0.001), 0.39 (p < 0.05), and 0.91 (p < 0.001), respectively; the average correlation coefficient of vegetation cover reached 0.605, and its spatial and temporal variations were also the main factors in the increase or decrease in ET in the study area.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.