Evapotranspiration of Winter Wheat in the Semi-Arid Southeastern Loess Plateau Based on Multi-Source Satellite Data

Abstract

Continuous monitoring of evapotranspiration (ET) at high spatio-temporal resolutions is vital for managing agricultural water resources in arid and semi-arid regions. This study used the enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM) to calculate the ET of winter wheat between the green-up and milk stages in Linfen Basin, a typical, semi-arid area of the Loess Plateau, at temporal and spatial resolutions of 30 m and 8 d, respectively. We then analyzed the impact of meteorological factors on ET and its variation during the main growth period of winter wheat. The fused ET data displayed the spatial details of the OLI ET data better and could accurately reflect ET variation and local sudden variations during the main growth period of winter wheat. Moreover, winter wheat ET in rain-fed areas is more heavily influenced by meteorological factors, and the effect is more direct. Affected by the synergistic effect of wind velocity, precipitation, and temperature, the ET of winter wheat in rain-fed area was lower in the green-up stage. Then, ET gradually increased, reaching its maximum in the heading–grain filling stage. At the jointing stage, temperature had a significant effect on ET. A combination of precipitation and temperature had the greatest impact on the ET of winter wheat in the heading–filling stage. In the milk stage, meteorological factors had a minor impact on ET. This study serves as a reference for ET in winter wheat in semi-arid areas and its influencing meteorological factors, which can assist in drought mitigation and regional food security strategies.

1. Introduction

ET is a key link in water and energy cycles. It is vital for maintaining the surface water and energy balance, and it plays an important role in regulating regional climates and even the global climate [1,2]. In arid and semi-arid regions, ET has a dominant role in the water budget, with approximately 90% of water outputs occurring in the form of ET [3]. In farmland ecosystems, ET accounts for as much as 99% of water outputs [4]. As a result, accurately estimating farmland ET in arid and semi-arid regions and studying the change laws and dynamic processes of farmland ET, and the influencing factors, can help to improve farmland water-use efficiency and promote sustainable regional agricultural development. Farmland ET includes plant transpiration and soil evaporation [5]. Using traditional methods, such as large lysimeters, eddy correlation, and flux meters, to obtain the single-point farmland ET is costly, time-consuming, and constrained by observation conditions, making significant errors likely [6,7]. Furthermore, single-point ET does not reveal spatial patterns [8]. Remote sensing imaging has the unique advantage of facilitating the monitoring of a large area. It can transmit and record surface conditions as the dynamic process of ET occurs and calculate surface ET over a large area, thereby overcoming the shortcomings of traditional methods [9,10].

Many recent studies successfully estimated regional surface ET based on remote sensing data. A common method is to calculate sensible heat flux using land-surface temperatures from remote sensing and atmospheric temperatures from meteorological observations, and then calculate ET using a surface energy balance method, such as SEBAL [11], the surface energy balance system (SEBS) [12], and the simplified surface energy balance index (S-SEBI) [13]. The calculation of turbulent heat flux in the SEBS model is complex and requires many atmospheric parameters [14]. The S-SEBI model requires temperature extremes at specific locations [15]. The advantages of the SEBAL model are its straightforward data acquisition, clear physical concepts, high degree of precision, and good adaptability, and it can be used to calculate surface ET at large spatial and temporal scales [16]. Additionally, when calculating ET, the SEBAL model considers the spatial heterogeneity of arid and semi-arid regions. The model has been used widely [17,18].

Because the ET of farmland crops changes rapidly during the growing season, multiple time-series data during the critical growth period of crops and remote sensing images with high spatial resolution at specific time periods are required to ensure accuracy [19], but existing remote sensing data lack the necessary temporal or spatial resolution. To overcome this obstacle, scholars have used the spatio-temporal fusion of multi-source remote sensing data, i.e., obtaining both the time phase information of high temporal resolution remote sensing data and the spatial textural information of relatively high-spatial-resolution data to simulate the spatio-temporal dynamics of the ET of farmland crops [20,21]. The spatio-temporal fusion of images, therefore, fuses known “temporally dense” low-spatial-resolution images with “temporally sparse” high-spatial-resolution images corresponding to time points to generate a “temporally dense” sequence of higher spatial-resolution images corresponding to the sequence of low-spatial-resolution images [22].

Current spatio-temporal fusion models can be roughly divided into the following three categories. The first is transformation-based, which process data by transforming them. To obtain a high-resolution image at the moment of prediction, wavelet transformation is used to represent its inverse transformation. The advantage of this method is that detailed spectrum information is preserved, but the disadvantage is the unsolvable problem of mixed pixels and low fusion accuracy [23,24]. The second type of spatio-temporal fusion model is machine learning-based, which uses sparse representation to process images. These models use a two-layer spatio-temporal fusion model, which accounts for differences in data resolution and can capture changing surface reflectance in images. However, they have high data requirements and poor efficiency, so they are not suitable for large-scale studies [25,26]. The third type of model is a reconstruction-based fusion model, which assumes that the reflectance of low-spatial-resolution pixels can be regarded as a linear combination of the reflectance of high-spatial-resolution pixels. Most models use the pixel decomposition method. These models consider both spatial heterogeneity as well as non-linear features of reflectance change. In areas with high heterogeneity, the algorithm has a substantially improved fusion accuracy while retaining more spatial details [27,28]. Typical models include the spatial and temporal adaptive reflectance fusion model (STARFM) [29], the enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM) [30], and the spatial and temporal non-local Filter-based fusion model (STNLFFM) [31]. ESTARFM improves the method for determining the spectral weight coefficient in the STARFM model by fully accounting for the heterogeneity of pixels. It introduces a conversion coefficient for the reflectance change of pure and mixed pixels to improve the accuracy of the model’s reflectance fusion in areas with high heterogeneity [22]. Its fusion results are relatively good, and it is highly effective in remote sensing applications over large areas [32]. It should be noted that existing spatio-temporal data fusion models focus on the fusion of the main inversion parameters [33], such as surface albedo, spectral index, and surface temperature. This does not account for the influence on ET of remote sensing and atmospheric factors, such as soil moisture [34] and vegetation distribution [35], which produces certain limitations.

Linfen Basin, which is in the southeast part of the Loess Plateau in north-central China, is typical of semi-arid regions and a major production area of winter wheat in Shanxi Province [36]. It has been affected by a decline in the winter monsoon, i.e., weak, warm, humid air flows from the Pacific Ocean, which have resulted in severe precipitation and soil water shortages in spring [37]. ET is greater than precipitation, so droughts are a frequent occurrence and have a severe impact on winter wheat yields. It is of great significance to study the influence of climatic conditions on crop yield and ET in this region for ensuring food security and alleviating water shortage. Studies have shown that various meteorological factors, such as temperature, water vapor pressure, sunshine hours, and wind velocity, have an important impact on the growth and ET of winter wheat [38,39]. However, influences such as data accuracy, spatio-temporal scale, underlying surface coverage, and the landscape environment cause substantial differences in winter wheat ET and its main driving factors, with notable regional characteristics [40]. Furthermore, few previous studies have considered the synergistic effects of meteorological factors. These effects have complex responses, and there is no linear relationship with strict statistical standards, making them difficult to convey quantitatively using only the analysis of partial and multiple correlations [41]. Interactions can ensure immunity to multi-collinearity. The synergistic effects of meteorological factors can be used to judge and demonstrate the explanatory power of combined and independent meteorological factors on winter wheat ET, which has extensive prospects [42].

This study looks at the winter wheat-growing area of Linfen Basin and explores the use of the spatio-temporal fusion of multi-source remote sensing data to obtain ET data for winter wheat during the growth stage in the region at a spatial resolution of 30 m and a temporal resolution of 8 d. On this basis, we analyzed the spatio-temporal distribution features of ET during the main growth stage of winter wheat and the influence of meteorological factors to provide a scientific basis for the assessment of changes in ET of farmland in semi-arid regions and the impact of meteorological factors on crop ET, to inform strategies for the sustainable development of agriculture.

2. Materials and Methods

2.1. Study Area

Linfen Basin is located in the southern temperate region in the southeastern part of the Loess Plateau, spanning Shanxi, Shaanxi, and Gansu provinces (35°20′–36°35′N, 110°30′–111°55′E). It borders the Taiyue Mountains to the north and is surrounded by the Taihang Mountains and Lvliang Mountains to the southeast and northwest, respectively. The Fen River watershed runs through the region and meets the Yellow River in the southwest of Linfen Basin. From north to south, it contains the eight county-level administrative areas of Hongtong County, Yaodu District, Xiangfen County, Quwo County, Houma City, Xinjiang County, Jishan County, and Hejin City. This region has a climate with latitudinal zonality and notable vertical variation [37]. It is suitable for the planting and cultivation of crops, such as wheat and corn. The region’s planted area of winter wheat accounts for 58.36% of cultivated land, making it an important production base in north China [36]. The spatial distribution of winter wheat in rain-fed and irrigated areas is shown in Figure 1. The region is a typical, semi-arid region of the Loess Plateau, with frequent droughts in spring and summer, high levels of soil evaporation, and varied topography and water resources, which restrict the development of the region’s agricultural areas.

2.2. Data Sources

This study used MOD16A2 ET product data calculated using the surface energy balance method (https://ladsweb.modaps.eosdis.nasa.gov (accessed on 27 January 2022)). Actual ET estimates included three parts: wet-canopy evaporation, canopy transpiration, and soil evaporation. The spatial resolution of the data was 500 m and the temporal resolution was 8 d. Using the MODIS Reprojection Tool (MRT), MOD16A2 data (orbit numbers: H26V05 and H27V05) for Linfen Basin from March to May 2019 were spliced and converted into the GeoTiff format with WGS-1984 and UTM data.

Landsat-8 Operational Land Imager (OLI) remote sensing images for 15 March, 31 March, 16 April, and 2 May 2019 were downloaded from the Geospatial Data Cloud (http://www.gscloud.cn (accessed on 7 October 2021)), and surface reflectance was calculated using radiometric calibration and the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) atmospheric correction tool. The normalized difference vegetation index (NDVI) was calculated using near-infrared and red-light bands. Surface emissivity was derived from the NDVI threshold, and land-surface temperature was calculated using the radiative energy transfer equation [43]. The albedo-to-nadir (AN) coefficient was calculated according to high-quality (QA = 0) MODIS BRDF data in the MCD43A1, and Landsat-8 OLI spectral albedo was determined by multiplying the AN coefficient by the surface albedo. The Landsat-8 OLI surface shortwave broadband albedo was calculated according to the Landsat-8 OLI narrowband-to-broadband conversion coefficient [44].

Observation data of daily atmospheric pressure, sunshine hours, relative humidity, water vapor pressure, temperature, precipitation, wind velocity, and mean maximum temperature were obtained from the China Meteorological Data Service Center (http://data.cma.cn (accessed on 7 October 2021)). Study area topographical data used Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Map (30 m) products, with data obtained from the Geospatial Data Cloud. We selected 23 meteorological stations in and around the study area in 2019 and used the ANUSPLIN tool (with latitude and longitude as independent variables and elevation as a covariate) to conduct spatial interpolation on the meteorological data. Regional-scale meteorological data corresponding to the green-up, jointing, heading–filling, and milk stages of winter wheat, and the interpolated results of wind velocity and temperature were then matched to the dates of Landsat-8 OLI data, with a pixel size of 30 m. Based on the daily observation data from E601B evaporator pans (ET_pan), the measured ET during the main growth period of winter wheat at the dry-farming and irrigated sample points in Linfen Basin in 2019 was calculated using the conversion coefficient of small evaporation pans into the E601B evaporator pan (K) [45], pan coefficient (K_p) [46,47], and crop coefficient (K_c) [48].

2.3. Analysis Methodology

2.3.1. SEBAL Model

The SEBAL model is a single-layer model based on energy balance. It can estimate latent heat flux and other energy balance components with a higher accuracy and better practical value, without needing information on soil, crops, and management practices [18]. The energy balance equation of the SEBAL model is:

$$\lambda ET=R_n-G-H$$

(1)

where $\lambda ET$ is the latent heat flux ($\mathrm{W}/{\mathrm{m}}^2$), and $\lambda$ is the latent heat of vaporization ($\mathrm{J}/\mathrm{kg}$), which can be obtained from the relationship with surface temperature $T_s(K)$:

$$\lambda =[2.501-0.002361\times (T_{\mathrm{s}}-273.15)]\times {10}^6$$

(2)

$R_n$ is the surface net radiation flux ($\mathrm{W}/{\mathrm{m}}^2$), which represents the net balance of surface shortwave and longwave radiation, which are the main sources of surface energy. They are calculated as follows:

$$R_n=\left(1-\alpha \right)\times R_s+\left(\epsilon \times R_{L\downarrow }-R_{L\uparrow }\right)$$

(3)

where $\alpha$ is the surface albedo, $\epsilon$ is the surface emissivity calculated using the NDVI threshold method, $R_s$ is the total solar radiation $\left(\mathrm{W}/{\mathrm{m}}^2\right)$, $R_{L\downarrow }$ is the incoming longwave radiation $\left(\mathrm{W}/{\mathrm{m}}^2\right)$, and $R_{L\uparrow }$ is the reflected longwave radiation $\left(\mathrm{W}/{\mathrm{m}}^2\right)$.

$G$ is the soil heat flux ($\mathrm{W}/{\mathrm{m}}^2$), and $H$ is the sensible heat flux ($\mathrm{W}/{\mathrm{m}}^2$), calculated as follows:

$$G=\left[\frac{T_s}{\alpha }\left(0.0038\alpha +0.0074{\alpha }^2\right)\times \left(1-0.98NDV{I}^4\right)\right]\times R_n$$

(4)

$$H=\frac{{\rho }_{air}\times C_p\times d_T}{r_{ah}}$$

(5)

where $NDVI$ is NDVI, ${\rho }_{air}$ is the air density (1.293 $\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$), $C_p$ is the specific heat capacity of the air at a constant pressure (1004$/\mathrm{k}\mathrm{g}\times K$), $T_s$ is temperature ($K$), and $r_{ah}$ is the aerodynamic impedance $\left(\mathrm{s}/\mathrm{m}\right)$, calculated as follows:

$$r_{ah}=\ln (Z_2/Z_1)/U_{\ast }\times k$$

(6)

where $Z_1$= 0.01 m, $Z_2$= 2 m; $k$ is the Karman constant (0.41); $U_{\ast }$ is the frictional wind velocity under neutral stability conditions, calculated as follows:

$$U_{\ast }=k\times U_r/\ln (Z_r/Z_{\mathrm{om}})$$

(7)

where $U_r$ is the wind velocity at height, $Z_r$. In general, it is assumed that it is stable above 200 m, i.e., Z_r = 200, and it is obtained using the linear regression of the wind velocity from the meteorological data at 2 m; Z_om is the surface roughness, calculated as follows:

$$Z_{\mathrm{om}}=\exp (5.65\times NDVI-6.32)$$

(8)

$d_T$ is the temperature difference between reference height $Z_1$= 0.01 m and reference height $Z_2$= 2 m, calculated as follows:

$$d_T=T_{Z1}-T_{Z2}=a{T}_s+b$$

(9)

where a and b are the regression coefficients. Cold and hot pixel points were selected, correction factors were introduced, and iterative operations were repeated until the value of H was stable, and the values of a and b were determined.

Using the above equations, we calculated R_n, G, and H in the model. Instantaneous ET was calculated according to Equation (1). The model assumes that the evaporation ratio is constant throughout a day, so the instantaneous ET becomes the daily ET. The daily ET of the Landsat-8 OLI data was calculated for 15 March, 31 March, 16 April, and 2 May 2019 according to the above model. We denoted the calculated results of the SEBAL model OLI ET.

2.3.2. ESTARFM Model

The ESTARFM model assumes that the remote sensing data of a certain period from various sources are correlated and that the systematic deviation between them is stable. Based on at least two pairs of high- and low-resolution images before and after the forecast date, and one low-resolution image on the forecast date, the weights and conversion coefficients between corresponding pixels were calculated, thereby simulating high spatio-temporal resolution data for the forecast date [30]. In this study, we used OLI ET (30 m) and MODIS ET (500 m) product data inverted (based on the SEBAL model) to simulate 30 m ET data for the forecast dates, calculated as follows:

$$ET(X_{\omega /2},Y_{\omega /2},t_p)=T_m\times E{T}_m(X_{\omega /2},Y_{\omega /2},t_p)+T_n\times E{T}_n(X_{\omega /2},Y_{\omega /2},t_p)$$

(10)

$$E{T}_{\mathrm{k}}(X_{\omega /2},Y_{\omega /2},t_p)=E{T}_O(X_{\omega /2},Y_{\omega /2},t_k)+{\displaystyle \sum _{\mathrm{i}=1}^N{W}_i{V}_i\left[MET(X_i,Y_i,t_p)-MET(X_i,Y_i,t_k)\right]}$$

(11)

where $ET$ is the ET of the final prediction date; $t_p$ is the prediction date; $E{T}_m$ is the predicted ET at time, $m$; $E{T}_n$ is the predicted ET at time, $n$; $E{T}_k$ is the high-resolution ET based on time prediction, $k$; $E{T}_O$ is the OLI ET inverted according to the SEBAL model at time $k$; MET is the MODIS ET; $\omega$ is the size of the sliding window, with $\left(X_{\omega /2},Y_{\omega /2}\right)$ indicating the position of the central pixel of the sliding window; N is the number of similar pixels in the sliding window; $\left(X_i,Y_i\right)$ is the position of similar pixel i; $V_i$ is the conversion coefficient, i.e., the regression coefficient of two groups of high- and low-resolution ET at the corresponding date; $W_i$ is the weight of the similar pixel i, calculated as follows:

$$W_i=\frac{1/D_i}{{\displaystyle {\sum }_{i=1}^{N}1/D_i}}$$

(12)

$$D_i=\left(1-R_i\right)\times d_i$$

(13)

$$R_i=\frac{E\left[\left(E{T}_i-E\left(E{T}_i\right)\right)\left(ME{T}_i-E\left(ME{T}_i\right)\right)\right]}{\sqrt{D\left(E{T}_i\right)}\times \sqrt{D\left(ME{T}_i\right)}}$$

(14)

$$d_i=1+\frac{\sqrt{\left(X_{\omega /2}-X_i\right)+\left(Y_{\omega /2}-Y_i\right)}}{\omega /2}$$

(15)

where $R_i$ is the spectral weight of the similar pixel, i; X is the distance weight of the similar pixel, i.

T_m and T_n are the weight factors at times m and n, respectively, calculated as follows:

$$T_k=\frac{1/\left|{\displaystyle {\sum }_{j=1}^{\omega }{\displaystyle {\sum }_{i=1}^{\omega }MET\left(X_i,Y_i,t_k\right)}}-{\displaystyle {\sum }_{j=1}^{\omega }{\displaystyle {\sum }_{i=1}^{\omega }MET\left(X_i,Y_i,t_p\right)}}\right|}{\left(1/\left|{\displaystyle {\sum }_{j=1}^{\omega }{\displaystyle {\sum }_{i=1}^{\omega }MET\left(X_i,Y_i,t_k\right)}}-{\displaystyle {\sum }_{j=1}^{\omega }{\displaystyle {\sum }_{i=1}^{\omega }MET\left(X_i,Y_i,t_p\right)}}\right|\right)}$$

(16)

In addition, relative errors (RE) and absolute error (AE) were calculated to test the accuracy of fusion ET. The RE is the ratio of the difference between the measured ET and ESTARFM ET to the measured ET. The AE is the difference between OLI ET and ESTARFM ET.

2.3.3. GeoDetector

GeoDetector is a spatial statistical method for detecting spatial divergence and explaining its driving factors. It is used to explore the spatial heterogeneity between dependent and independent variables and measure the explanatory power (i.e., q value) of independent variables with respect to dependent variables [41]. It is calculated as follows:

$$PD=1-{\displaystyle \sum _{h=1}^L{N}_h{\sigma }_h^2}/N{\sigma }^2$$

(17)

where PD is the explanatory power of winter wheat ET with a value range of 0–1; the larger the value, the stronger the explanatory power of winter wheat ET; h (1, 2, …, L) is the variable classification or partition; N_h and N are the number of layers, h, and regional units, respectively; σ_h² and σ² are the variance of the Y value in the h layer and region, respectively.

Having calculated the PD value of a single meteorological factor, it was necessary to calculate PD (x₁), PD (x₂), and PD (x₁ꓵx₂) to judge and explain the explanatory power of meteorological factors, in combination and independently (increasing or decreasing), on winter wheat ET. Based on the above equation, it was possible to calculate the effects of various meteorological factors (atmospheric pressure, sunshine hours, relative humidity, water vapor pressure, temperature, precipitation, wind velocity, and mean maximum temperature), cultivation methods (rain-fed and irrigated), and growth stages (green-up, jointing, heading–filling, and milk stages) on the ET of winter wheat (Figure 2).

3. Results and Analysis

3.1. Analysis of ESTFARM ET Fusion Results

Based on the ESTARFM algorithm, we fused OLI ET and MODIS ET to obtain ESTARFM ET, with a spatial resolution of 30 m and a temporal resolution of 8 d. Figure 3 shows MODIS ET and OLI ET inverted by the SEBAL model, and ESTARFM ET on 15 March, 31 March, 16 April, and 2 May 2019 in the experimental area, and the scatter plot presents a comparison of OLI ET and ESTARFM ET in Linfen Basin. A comparative analysis of MODIS ET and OLI ET inverted by the SEBAL model, and ESTARFM ET revealed that, compared with the original MODIS ET, ESTARFM ET substantially improved the spatial resolution and provided good spatial detail information. In addition, it was better at showing spatial differences between small surface objects. The correlation coefficient of ESTARFM ET and OLI ET using the SEBAL model was 0.368 on 15 March (which was affected by clouds and rain), and the correlation coefficients on the other dates were all above 0.65, with the scattered points distributed near the 1:1 line, indicating a relatively high linear correlation between ESTARFM ET and OLI ET and highly similar spatial distribution trends. The deviation between ESTARFM ET and OLI ET was 0.003–0.633, and the root mean square error (RMSE) was 0.182–0.754. The values of MODIS ET and OLI ET were quite different on 16 April, resulting in the highest deviation (0.633) between ESTARFM ET and OLI ET, with an RMSE of 0.754 on that date, while the deviation and RMSE on the rest of the dates were lower. The absolute errors between OLI ET and ESTARFM ET maintained a good normal distribution trend on 15 March and 31 March, and the main distribution intervals of the absolute errors were −0.4–0.4, −0.2 and 0.2, respectively. On 16 April, most absolute errors were concentrated in 0.8–0.9, and the frequency was 500 from −0.5 to 0.7. It still showed a certain normal distribution trend. The main distribution range of the absolute error between OLI ET and ESTARFM ET was −0.6 to 0.6 on 2 May. This indicated that the fused ET was consistent with the OLI ET inverted by the SEBAL model and that the data accuracy was good.

We extracted ESTARFM ET and MODIS ET of the sample points in Hongdong County, Xiangfen County, Xinjiang County, Quwo County, Hejin City, and Houma City. Then, we used the Savitzky–Golay (S-G) filter to reconstruct the ESTARFM ET and MODIS ET for the main growth period (March to May) of winter wheat and compared them with the measured ET (Figure 4). It can be seen that MODIS ET was affected by the low spatial resolution, with no marked response to sudden changes in winter wheat ET and a lagging effect. Furthermore, the reconstructed MODIS ET change curve was very low in the simulation, indicating poor suitability. ESTARFM ET reflected the variation in the measured ET better. During the green-up and jointing stages, ET mainly fluctuated upward, with the maximum value occurring at the heading–filling stage (DOY120), with ET falling sharply in the milk stage (DOY140). ESTARFM ET was consistent with measured ET in terms of the variation range and direction of change. Due to the effect of the temporal resolution, ESTARFM ET responded slowly to sudden substantial changes in the measured ET. The RE between the measured ET and ESTARFM ET during the main growth period of winter wheat was 2.11%, whereas the relative error of dry-farmed sampling points was 1.91%, and that of irrigated sampling points was 2.27%. This indicates that ESTARFM ET is better than MODIS ET at representing ET variation during the main growth period of winter wheat in Linfen Basin and that it has a good suitability.

3.2. Temporal and Spatial Characteristics of Winter Wheat ET in Linfen Basin

ESTARFM ET was calculated for the main growth period of winter wheat in Linfen Basin (Figure 5). Spatially, ESTARFM ET was generally higher in central and southern regions, with high values concentrated on both sides of the Fen River watershed that runs through the basin, reaching its maximum value at the boundary between Houma City and Quwo County in the hinterland of the basin. ESTARFM ET gradually decreased from the south-central area to the northeast and southwest. ET at the confluence of the Yellow River and the Fen River in the southwest was substantially higher than that in the Taiyue Mountains in the northeast. In the time series, ESTARFM ET trended upward from the green-up stage, with ET peaking in the jointing stage when ET intensity was greater. The maximum ET value was reached at the end of the jointing stage and the beginning of the heading–filling stage (DOY106). ET intensity decreased notably between the heading–filling stage and the milk stage. In the southwest of Linfen Basin, at the confluence of the Yellow River and Fen River watershed, the variation in ESTARFM ET during the main growth period of winter wheat was substantially higher than that in other areas.

The ET pixels of the winter wheat planting area in Linfen Basin were classified into rain-fed and irrigated areas (Figure 6). Changes in winter wheat ET in the rain-fed and irrigated areas were relatively synchronous, with both increasing from the green-up stage. The ET growth rate was highest from DOY98 to DOY106, and ET reached its maximum value at the end of the jointing stage and the beginning of the heading–filling stage, after which it decreased. The mean ET of winter wheat in the green-up stage, jointing stage, heading–filling stage, and milk stage in irrigated areas was 0.901/mm, 2.217/mm, 3.507/mm, and 2.058/mm, respectively. The total amount of ET during the main growth period of winter wheat was distributed from 200.12 mm to 256.24 mm, which was consistent with the ET inversion results of winter wheat in the North China Plain [49]. The mean ET of winter wheat in the various stages of growth was lower in rain-fed areas (0.718/mm, 1.944/mm, 3.056/mm, and 2.01/mm, respectively), and the increase in ET between the jointing stage and the heading–filling stage was substantially larger in the irrigated-farming samples than in the dry-farming samples. The mean winter wheat ET in irrigated areas was greater than the median value except at the green-up stage, indicating that there was a higher proportion of high ET values in the jointing, heading–filling, and milk stages. In rain-fed areas, the difference in median and mean values of winter wheat ET during different growth stages was not significant.

3.3. Winter Wheat ET Response to Meteorological Factors during the Main Growth Stages

Table 1. PD values of various meteorological factors (X1–X8 represent atmospheric pressure, sunshine, relative humidity, water vapor pressure, temperature, precipitation, wind velocity, and mean maximum temperature, respectively), winter wheat growth stage (X9), and winter wheat cultivation method (X10).

	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10
PD	0.305 **	0.335 **	0.506 **	0.804 **	0.753 **	0.686 **	0.491 **	0.754 **	0.813 **	0.924 **

** p < 0.01.

shows the explanatory power (PD value) of winter wheat ET for various growth stages, meteorological factors, and cultivation methods in Linfen Basin. The order of ET PD values was cultivation method, growth period, water vapor pressure, mean maximum temperature, temperature, precipitation, relative humidity, wind velocity, sunshine hours, and atmospheric pressure, all of which passed the significance test. Of these, the explanatory power values of the cultivation method, growth period, water vapor pressure, mean maximum temperature, temperature, and precipitation were above 0.600, whereas the explanatory power values of the relative humidity, wind velocity, sunshine hours, and atmospheric pressure were lower. The explanatory power of growth period with regard to ET was 0.813, indicating that the water requirements of winter wheat differed substantially at various stages of growth. The explanatory power of the farming method was 0.924, indicating that a change in farming method had a significant impact on ET. Further analysis and discussion of the meteorological factors that affect winter wheat ET is required with a distinction between rain-fed and irrigated areas.

We calculated the explanatory power values of the various growth stages and meteorological factors on winter wheat ET (

Table 2. PD values of meteorological factors (X1–X8 represent atmospheric pressure, sunshine, relative humidity, water vapor pressure, temperature, precipitation, wind velocity, and mean maximum temperature) and various winter wheat growth stages (X9) in rain-fed and irrigated areas.

	X1	X2	X3	X4	X5	X6	X7	X8	X9
rain-fed areas	0.287 **	0.505 **	0.714 **	0.937 **	0.931 **	0.967 **	0.825 **	0.967 **	0.966 **
irrigated areas	0.356 **	0.307 **	0.559 **	0.719 **	0.724 **	0.565 **	0.282 **	0.691 **	0.551 **

** p < 0.01.

), which showed substantial differences in PD values between rain-fed and irrigated areas. The PD value of the explanatory power of the growth stage factor in rain-fed areas reached 0.966, which was 0.415 greater than that in irrigated areas (0.551), indicating that winter wheat ET in rain-fed areas was greatly affected by differences in growth stage. Except for the atmospheric pressure factor, the explanatory power values of meteorological factors for ET in rain-fed areas were higher than those in irrigated areas, suggesting that winter wheat ET in rain-fed areas was affected more by natural meteorological factors. The explanatory power values of precipitation and wind velocity were markedly different. In rain-fed areas, they were 0.967 and 0.825, respectively, and the PD values for irrigated areas were 0.565 and 0.282, respectively. These results demonstrate that, due to the influence of irrigation and other human factors, meteorological factors had less impact on winter wheat ET in irrigated areas than in rain-fed areas. Meteorological factors had a more direct impact on winter wheat ET in rain-fed areas, which is the major reason why the fusion-based ET results were more accurate for rain-fed areas than for irrigated areas.

The PD values of various meteorological factors were also calculated for winter wheat in rain-fed areas in the green-up, jointing, heading–filling, and milk stages (

Table 3. PD values of atmospheric pressure (X1), sunshine (X2), relative humidity (X3), water vapor pressure (X4), temperature (X5), precipitation (X6), wind velocity (X7), and mean maximum temperature (X8) in different growth stages of winter wheat in rain-fed areas.

	X1	X2	X3	X4	X5	X6	X7	X8
green-up stage	0.098	0.348 **	0.546 **	0.524 **	0.564 **	0.475 **	0.579 **	0.609 **
jointing stage	0.333 **	0.311 **	0.648 **	0.888 **	0.949 **	0.845 **	0.701 **	0.932 **
heading-filling stage	0.180 *	0.745 **	0.755 **	0.903 **	0.964 **	0.976 **	0.839 **	0.780 **
milk stage	0.030	0.103 *	0.405 **	0.338 **	0.331 **	0.334 **	0.153 *	0.211 *

* p < 0.05; ** p < 0.01.

). For dry-farmed winter wheat, the PD values of the various growth stages and meteorological factors were markedly different. The explanatory power of the mean maximum temperature in the green-up period was highest (PD = 0.609). The explanatory power values of the atmospheric pressure and sunshine were lower, and the explanatory power of atmospheric pressure did not pass the significance test. In the jointing stage, except for sunshine, the explanatory power values of the other meteorological factors were higher than in the green-up stage. The PD values of the temperature and mean maximum temperature reached 0.949 and 0.932, respectively, and those of water vapor pressure and precipitation both exceeded 0.840. In the heading-filling stage, precipitation replaced temperature as the factor with the highest PD value. Except for atmospheric pressure, the explanatory power values of the other meteorological factors all exceeded 0.700. The explanatory power of sunshine was substantially higher in this stage than in the jointing stage. The milk stage had the lowest PD values for meteorological factors over all stages. The explanatory power of relative humidity was the highest, and that of atmospheric pressure was the lowest, which also failed the significance test.

We calculated interactions between meteorological factors at different growth stages for winter wheat in rain-fed areas (Figure 7). The interactive PD values of meteorological factors were all greater than the PD values of any one factor, indicating that factor interactions produced a mutual and non-linear enhancement, with no independent or attenuating relationship. Atmospheric pressure ꓵ mean maximum temperature in the green-up and heading–filling stages, atmospheric pressure ꓵ sunshine hours, and atmospheric pressure ꓵ relative humidity in the heading–filling stage displayed non-linear enhancement. The PD values of interactions between wind velocity and other meteorological factors increased substantially during the green-up period. Except for the PD value of the interaction with atmospheric pressure (0.589), the PD values of the other interactions were all above 0.770, while the PD value of the wind velocity was just 0.579. This indicates that this factor had a greater synergistic influence, and it had less of a direct influence. Since temperature had a considerable influence on ET in the jointing stage (PD = 0.949), the PD values of its interactions with other meteorological factors were all maximum values, but the increases were not significant, indicating that temperature alone had a marked, direct impact on winter wheat ET in the jointing stage. In the heading–filling stage, the interaction PD value of precipitation and temperature was 0.997, which was the highest interaction PD value across all the growth stages. In the milk stage, the meteorological factor interaction PD values were not substantially higher than those of single factors, indicating that synergistic effects had a negligible impact on ET during this period.

4. Discussion

4.1. Improvements to and Limitations of ESTARFM ET

The spatial resolution of MODIS ET data is low (500 m), making it difficult to clearly make out the details of spatial textures of various surface objects. In the winter wheat planting area of Linfen Basin, there is considerable farmland fragmentation, with many mixed pixels containing both vegetation cover and bare land [36]. As a result, the upper range value of MODIS ET is low during the growth period of winter wheat. This is consistent with the experimental results in the Pampas Region [50]. Furthermore, mixed pixels were the primary reason for ET curve lag and inaccurate simulations [20]. Based on the ESTARFM spatio-temporal fusion model, Landsat-8 OLI (30 m) data could be spatially combined at corresponding times. Most of the spatial details of the latter largely solved the problem of mixed pixels in the MODIS ET data and effectively resolved inaccurate ET variations [30]. It also resolved the scattered time series of Landsat-8 OLI data, facilitating the analysis of the spatial details and textures of winter wheat during the main growth period on a continuous time scale. Additionally, compared with the fusion of the main inversion parameters for ET in previous studies, MODIS ET and OLI ET were both calculated using surface energy balance, which considers the impact on ET of remote sensing and atmospheric features, giving good homogeneity and reducing uncertainty in the data fusion process [22].

In the heading–filling stage of winter wheat, Landsat-8 OLI data were greatly affected by clouds and rain, which caused errors in the determination of similar pixels and the calculation of conversion coefficients and time weights [44]. As a result, the RMSE and deviation of fused ESTARFM ET and OLI ET were too high. Moreover, ET changes rapidly over time; however, the ESTARFM ET change curve for the main growth period of winter wheat in Linfen Basin based on S-G filtering was influenced by the time scale effect, but did not respond sufficiently to sudden changes in ET [51]. In future research, remote sensing images with a higher temporal resolution could be used, such as advanced Himawari imager (AHI) products [52], which have a temporal resolution of 10 min, to improve the speed of response to sudden ET changes. Sentinel-2A and Sentinel-2B data (10 m) can also be used, which have a shorter revisiting cycle and higher spatial resolution, thereby reducing the impact of cloud and rain on image quality and improving ET fusion accuracy. Moreover, deep convolutional networks to image fusion problems reduce the need manually set the fusion rules [53], realize data fusion by learning the implicit rules from a large number of data, make full use of the advantages of automatic feature extraction of convolution operations to extract key information from images, and allow for the fusion of multiple data sources and image reconstruction [54].

4.2. Variation Characteristics of ET and Its Effects on Growth and Development of Winter Wheat

Affected by the difference in the time division of the growth period, regional climatic conditions, and spatial scale, the total amount of ET based on the ESTARFM model (March to May) was different from the total amount of single-point ET [55,56]. It is necessary to calculate ET during the whole growth period of winter wheat, and combine point and surface data based on the assimilation method [57] in future studies. The temporal and spatial variation in ET in the winter wheat planting areas correspond to irrigation conditions and are affected by agricultural management practices [58]. The northeastern part of Linfen Basin has an undulating topography due to the Taiyue Mountains and a highly varied terrain [37]. The winter wheat in this area is mostly rain-fed and is susceptible to spring droughts, with little water for ET; thus, ET is low at each stage of growth. The central and southeastern parts of Linfen Basin are largely composed of plains, with the Fen River running through them, forming an irrigation area with favorable water and heat conditions [59], such that ET is relatively high at each stage of growth.

ET is an important constituent of the water balance of farmland, with a considerable influence on plant growth and yield [3,4]. Before and after the green-up stage, the growth of the leaves, stems, and roots of winter wheat is relatively slow, with less dry matter accumulation, and water demand is lower. ET is mainly soil ET, and temperature is the main factor limiting ET during this period [60]. Due to crop compensatory growth, water stress occurs before and after the green-up stage, with normal irrigation resuming during subsequent growth stages. The water deficit has little effect on plant height, final tillering, leaf area, and dry matter accumulation, so there is little sensitivity to water stress in the green-up stage. From the rapid to middle growth stages, sufficient sunshine ensures that winter wheat leaves, stems, and roots grow simultaneously, and leaf and stem length and volume increase exponentially, while the root system extends rapidly, notably increasing the absorption and utilization of water from the deep soil. Transpiration increases and soil evaporation increases due to the rise in temperature [61]. The daily mean ET of winter wheat is highest in the heading–filling stage. If there is insufficient soil water during this stage of growth, it affects the efficiency and speed of photosynthesis, reduces the synthesis of starch, protein, and organic matter, and greatly reduces grain weight, which will affect final growth and the yield [55,62]. Winter wheat is most sensitive to a water deficit at this stage, so it is necessary to ensure it has sufficient water for its growth requirements. In Linfen Basin, winter wheat is usually irrigated during this period, but winter wheat in rain-fed areas relies on precipitation, resulting in much higher winter wheat ET in irrigated areas than in rain-fed areas. In the milk stage, once nutrients are transported to the grain, leaves turn from green to yellow, photosynthetic capacity decreases, plants gradually wither, roots begin to die, transpiration decreases, and the impact of meteorological factors on ET is attenuated. The variation trend of winter wheat ET in the main growth stages was consistent with the simulation results of winter wheat ET based on the SBES model in Henan Province, with the estimation accuracy in this paper (RE = 2.11%) being higher than the latter (RE = 6.74%) [6]. Irrigation in the late milk stage will cause new tillers to grow from the base of old stems, which consumes nutrients and decreases the thousand-kernel weight [63]. Thus, field irrigation should be reduced during this stage of growth to increase crop yield. We can determine the sensitivity of winter wheat to water demand based on ET at different growth stages, and then improve the water-use efficiency of farmland, i.e., we can ensure adequate irrigation is carried out in the jointing stage and heading–filling stage when ET is high. Appropriate water stress at the early and late growth stages contributes to yield improvement.

4.3. Relationship between Winter Wheat ET and Meteorological Factors in Various Growth Stages

ET plays a key role in the soil–plant–atmosphere continuum and is a key factor affecting water and heat transfer in the land–atmosphere system [2]. Previous studies have shown that the waveform variation in winter wheat ET during the main growth period is greatly affected by meteorological factors [64]. The variation in winter wheat ET and its main driving factors is significantly different in different regions, showing obvious regional characteristics [38,39,65,66,67,68]. Understanding meteorological factor mechanisms can lead to the better utilization of agricultural water resources, thereby increasing yield. ET is a constituent of both the surface energy and water balance [1]. In the ET process, energy mainly originates from solar radiation. Temperature and sunshine hours are closely related to solar radiation and directly determined by net solar radiation. Temperature can increase the stomatal aperture of leaves, thereby increasing their physiological activity [61]. Moreover, crops absorb water from the soil, which requires sufficient root pressure, and temperature plays a decisive role in root pressure. Water vapor transport largely depends on actual water vapor pressure and wind velocity. There is a quadratic polynomial relationship between actual water vapor pressure and precipitation in the Loess Plateau region, reflecting the synergistic effect of temperature and relative humidity [65]. Thus, the meteorological factors that influence ET include sunshine hours, temperature, wind velocity, water vapor pressure, precipitation, and relative humidity. Water vapor pressure is the result of the synergistic effect of precipitation, temperature, and relative humidity. It directly influences the saturation deficit, and it has a dominant role and the greatest explanatory power in relation to ET rate variation [42]. The influence of meteorological factors on ET differs substantially, depending on irrigation methods and the water requirements of winter wheat at the various stages of growth [40]. In irrigated areas, winter wheat has sufficient water, so crops grow better in the jointing and heading–filling stages, the leaf area index increases substantially, and more precipitation is intercepted by the canopy and infiltrates the soil, which is a reason for the higher proportion of ET in this period [66]. Moreover, the increase in temperature accelerates ET, but the slow rise in temperature in the initial period is a major limiting factor, indicating the strong explanatory power of temperature and mean maximum temperature with regard to ET [61].

Because they are unaffected by human factors such as irrigation, meteorological factors have a greater impact on winter wheat ET in rain-fed areas. Cyclonic activity occurs frequently during the green-up stage, and higher wind velocity produces greater turbulence, causing water vapor to spread over a large area, resulting in strong ET and water stress due to a decrease in soil water content [67]. Additionally, due to wind velocity, the temperature rises slowly and precipitation is lower [68]. During this period, there is a marked synergistic effect from wind velocity, precipitation, and temperature. In the jointing stage, winter wheat grows quickly and is very sensitive to temperature changes. Low temperatures, especially in late spring, hinder wheat growth and inhibit transpiration [69]. Low temperature is also the primary factor in curbing soil evaporation. Thus, temperature has a marked, direct influence on ET in the jointing stage, and it has weak synergistic effects [61]. In the heading–filling stage, rain-fed areas of winter wheat are substantially affected by water stress, and precipitation is the primary factor affecting ET. Precipitation in rain-fed areas is mainly from decreasing temperature, which promotes the condensation of water vapor in the air and causes rain, so there is substantial synergy between precipitation and temperature during this stage. Furthermore, sunlight has a notable influence on stomatal apertures at this stage [42]. On sunny days, there is sufficient sunlight and strong solar radiation, leaves are physiologically active, and photosynthetic and transpiration rates are high [70]. In the milk stage, winter wheat growth almost stops, and leaves dry out and turn yellow [63]. The growth status of the crop determines its ET contribution, and the contribution of meteorological factors to ET is relatively low. The above results are helpful to understand the influence mechanism of meteorological factors on ET, and they are of great significance for the rational utilization of farmland water resources, the development of drought prevention measures, and the guarantee of regional food security.

In addition, in future studies, we will use machine learning models such as convolutional neural networks (CvNet), random forests (RF), and deep learning networks to complete the inspection of meteorological data, improve the classification accuracy, and solve the problem of inaccurate discretization in the GeoDetector.

5. Conclusions

The ESTARFM fusion algorithm preserved the temporal variation information of MODIS ET data and conveyed the spatial details of OLI ET data inverted using the SEBAL model with low data uncertainty. It can accurately reflect ET variation and localized sudden changes in the main growth stages of winter wheat, making it better suited to the region’s winter wheat growth characteristics.

ET spatio-temporal characteristics are closely related to the region’s topographical conditions, the growth and development of winter wheat, and farmland management methods. ET in the central and southeastern parts of Linfen Basin was higher than that in the northeast Linfen Basin. The lowest daily mean ET is in the green-up stage, after which ET gradually increases, ET is highest in the heading–filling stage, and ET decreases substantially in the milk stage.

Unaffected by human factors such as irrigation, winter wheat ET in rain-fed areas is more heavily influenced by meteorological factors, and the effect is more direct. The meteorological factors with the greatest explanatory power with regard to ET in the main growth stages of winter wheat in Linfen Basin are (in descending order) water vapor pressure, mean maximum temperature, temperature, precipitation, relative humidity, wind velocity, sunshine hours, and atmospheric pressure. Most have synergistic effects, with no independent or attenuating relationships.

In the green-up stage, winter wheat ET in rain-fed areas was greatly affected by the synergistic effect of wind velocity, precipitation, and temperature; in the jointing stage, the direct effect of air temperature on ET was substantial, and synergistic effects were weaker; in the heading-filling stage, the synergistic effect of precipitation and air temperature had the greatest influence on winter wheat ET; in the milk stage, meteorological factors had a minor effect on ET.