Evaluation of Different Methods on the Estimation of the Daily Crop Coefficient of Winter Wheat

Abstract

Various methods have been developed to estimate daily crop coefficients, but their performance varies. In this paper, a comprehensive evaluation was conducted to estimate the crop coefficient of winter wheat in four growth stages based on the observed data of weighing-type lysimeters and the high-precision automatic weather station in the Wudaogou Hydrological Experimental Station from 2018 to 2019. The three methods include the temperature effect method, the cumulative crop coefficient method, and the radiative soil temperature method. Our results suggest that the performance of these methods was different in each individual growth stage. The temperature effect method was better in the emergence-branching (RMSE = 0.06, r = 0.80) and heading-maturity stages (RMSE = 0.16, r = 0.94) because the temperature is suitable for crop growth during most of these two periods. The cumulative crop coefficient method was better in the greening-jointing (RMSE = 0.16, r = 0.88) and heading-maturity stages (RMSE = 0.20, r = 0.91) because this method is closely related to crop growth, which is vigorous during these two stages. The radiative soil temperature method was better in the emergence-branching (RMSE = 0.20, r = 0.35) and branch-overwintering stages (RMSE = 0.25, r = 0.52) as the energy balance can be ensured by the relatively high level of the effective energy during these periods. By comparing the estimation accuracy indices of the three methods, we found that the temperature effect method performed the best during the emergence-branching stage (RMSE = 0.06, MAE = 0.06, r = 0.80, d_IA = 0.88), branch-overwintering stage (RMSE = 0.13, MAE = 0.11, r = 0.44, d_IA = 0.55), and heading-maturity stage (RMSE = 0.16, MAE = 0.13, r = 0.94, d_IA = 0.97), while the cumulative crop coefficient method performed best during the greening-jointing stage (RMSE = 0.16, MAE = 0.13, r = 0.88, d_IA = 0.89). Based on this result, an integrated modelling procedure was proposed by applying the best method in each growth stage, which provides higher simulation precision than any single method. When the best method was adopted in each growth stage, the estimated accuracy of the whole growth process was RMSE = 0.13, MAE = 0.09, r = 0.98, d_IA = 0.99.

1. Introduction

According to data released by the Food and Agriculture Organization (FAO), the FAO Food Price Index in 2021 reached a 10-year high, and the incidence of moderate or severe food insecurity in the world increased from 22.6% in 2014 to 30.4% in 2020. In 2020, an estimated 720 to 811 million people, or 9.5 to 10.7 percent of the global population, will face food insecurity [1]. Agricultural irrigation accounts for 87% of total water use [2], but water resources are increasingly in short supply under the pressure of climate change and population growth [3]. Accurate estimation of crop water demand is of great significance for rational allocation of water resources.

Wheat is one of the main crops grown and eaten worldwide [4]. From 2011 to 2020, the global average annual planting area reached 219 million hectares, with an annual yield of 733 million tons [5]. Wheat is the second main crop in China [6]; the wheat output of China accounts for approximately 18% of the world’s output (ranking first in the world), and the sown area accounts for approximately 10% of that of the world [7]. Winter wheat is sown from October to November and harvested from May to June of the following year. Planting is performed in dry conditions with little rainfall. Accurate estimation of the daily evapotranspiration and crop coefficient of winter wheat is conducive to the formulation of accurate irrigation schemes to improve the yield of winter wheat and save water resources [8].

Evapotranspiration is an important link in the hydrological cycle and is involved in the surface energy balance and water balance [9]. The crop coefficient is the ratio of actual evapotranspiration and reference evapotranspiration of crops, reflecting the influence of soil, vegetation, and hydrometeorological conditions on evapotranspiration and is often used to calculate crop water requirements [10]. Finding a suitable crop coefficient estimation method is of great significance for further estimating actual evapotranspiration, making irrigation plans, and efficiently utilizing water resources.

Domestic and foreign scholars have performed much research on crop coefficient estimation methods. In the crop coefficient estimation method recommended by FAO, the basic crop coefficient constant is obtained by plotting tabulated values and drawing the crop coefficient curve with a simplified straight-line connection at each growth stage. Finally, the daily crop coefficient value is modified according to wind speed and humidity. This method has been widely adopted and is very convenient for practical application. Er-Raki et al. [11] directly used the basic crop coefficient provided by the FAO Irrigation and Drainage Paper No. 56 (FAO-56) to calculate the soil surface vegetation coverage. Vu et al. [12] compared the crop coefficient recommended by the FAO-56 with the field monitoring value and found that the applicability of the recommended value of FAO-56 was affected by crop variety and growth stage. Ali et al. [13] simulated the Kc curve using four crop coefficient estimation methods, including the FAO-56 recommendation method. Many scholars estimate the crop coefficient according to the crop growth character index and calculate the crop coefficient by measuring the leaf area, plant height, vegetation index, etc. The daily scale is usually obtained by linear interpolation. This method of estimation has strong physical significance and high accuracy. Spiliotopoulos et al. [14] estimated crop coefficients based on the vegetation index using mapping ET of high resolution and internalized calibration models. Zhang et al. [15] established two Ks regression models for crop coefficient inversion and found that the model established by TCARI/RDVI had a better correlation with the crop coefficient. Park et al. [8] estimated crop coefficients for cropland and mixed forest based on the normalized vegetation index, leaf area index, and soil moisture.

The above methods use the same method to estimate the crop coefficient at each growth stage, but the growth of crops at different growth stages and the hydrometeorological factors closely related to their growth are different. It is difficult to obtain satisfactory results with a single method, so it is necessary to select the best method to estimate the crop coefficient at different growth stages. Based on the temperature effect, Wang et al. [16] established a daily crop coefficient estimation model of winter wheat and summer maize with a high fitting degree and found that the estimation effect of each growth stage was different. Tang et al. [17] conducted a remote sensing study by UAV and found that there was a good linear correlation between biomass and the cumulative crop coefficient and cumulative transpiration in different treatment areas. Kool et al. [18] evaluated the dual source energy balance (TSEB) model using net radiation, soil heat flux, and surface temperature. Paulino et al. [19] proposed an empirical model based on linear multiple regression to estimate the number of fruits per plant of two sweet oranges and concluded that the estimates of the model varied greatly in the three growth stages. At present, most related studies from both China and abroad focus on crop coefficient estimation methods. Few studies have compared the differences in crop coefficient estimation methods in different growth stages and provided the best estimation methods in each growth stage. In this paper, we attempt to select three methods, the temperature effect method, cumulative crop coefficient method, and radiative soil temperature method, to compare and analyse their estimation results to provide the best estimation method of crop coefficients in different growth stages and provide a basis for the accurate calculation of evapotranspiration in each growth stage.

2. Materials and Methods

2.1. Experimental Area Profile

The measured experimental data of the Wudaogou hydrological experimental station were used. The station is located in Guzhen County (33°09′ N, 117°21′ E), Bengbu City, Anhui Province, in the southern Huaibei Plain, covering an area of 27,000 square meters. It is affected by a subtropical humid monsoon climate and temperate subhumid monsoon climate, with rain and heat at the same time. The rainfall in the area where the experimental station is located varies greatly from year to year and is unevenly distributed within the year. A total of 61.8% of the annual rainfall is concentrated in the flood season (from June to September). The main soil types in the experimental area are sandy ginger black soil and yellow tidal soil, and the main crops are wheat, corn, and soybean. There are 62 sets of nonweighing ground lysimeters and 10 sets of large-scale weighing lysimeters in the station, which can record submersible evaporation, soil moisture, and evapotranspiration. A high-precision weather station is set up about 50 m south of the weighing lysimeters to automatically monitor net radiation, soil heat flux, air temperature, wind speed, and other hydrometeorological elements every 10 min.

According to the measured hydrometeorological data from 1986 to 2021, the average annual temperature is 15.2 °C, the average annual rainfall is 929.5 mm, the average annual flood season rainfall is 582.4 mm, the average annual evaporation is 931.9 mm, the average annual relative humidity is 79.28%, the average annual sunshine duration is 1723.6 h, the average annual wind speed is 1.6 m/s, the annual average surface temperature is 18.4 °C, the annual average soil temperature at 10 cm is 16.6 °C, the maximum temperature is 41.5 °C, and the minimum temperature is −22.7 °C.

2.2. Experimental Facilities and Data Selection

The difference in the three crop coefficient estimation methods at different growth stages and daily scales was studied by using a large-scale lysimeter. The winter wheat relies on rainfall. The lysimeter model is FR101A, the resolution is 0.025 mm, the soil column height is 4.0 m, the diameter area is 2.0 square meters, and the weight data are automatically collected every 10 min. The depth of the groundwater level in the shallow buried area of the Huaibei Plain is 1 to 3 m [20], and the root system of winter wheat is concentrated within 1 m [21]. Therefore, experimental data from a lysimeter with a burial depth of 1.0 m and measured data from high-precision weather stations were selected from 11 November 2018 to 4 June 2019.

2.3. Division of Growth Stages

According to the actual growth status of winter wheat, the whole growth process can be divided into four stages: emergence-branching stage, branch-overwintering stage, greening-jointing stage, and heading-maturity stage. The classification of the stages is based on the characteristics of the crops in each stage. The parameters of each stage are shown in

Table 1. Growth stage division of winter wheat.

Stage of Growth	Emergence-Branching	Branch-Overwintering	Greening-Jointing	Heading-Maturity
Date	2018/11/11–2018/12/1	2018/12/2–2019/2/21	2019/2/22–2019/4/19	2019/4/20–2019/6/4
Number of days	21	82	57	46
ET0 at this stage	21.80 mm	67.80 mm	216.48 mm	220.67 mm
Proportion of total ET0	4.14%	12.87%	41.10%	41.89%
Average daily ET0	1.04 mm	0.83 mm	3.80 mm	4.80 mm

.

2.4. Crop Coefficient, Actual Evapotranspiration, and Reference Evapotranspiration

The crop coefficient is divided into a single crop coefficient and a double crop coefficient. The single crop coefficient involves fewer factors and has higher estimation accuracy, so the single crop coefficient is selected. According to the definition of the crop coefficient, the calculation method is shown in Equation (1):

$${\mathrm{K}}_{\mathrm{c}}=\frac{\mathrm{ET}}{{\mathrm{ET}}_0}$$

(1)

where K_c is the crop coefficient; ET is the actual evapotranspiration (mm); and ET₀ is the reference evapotranspiration (mm).

The actual evapotranspiration (ET) of winter wheat is automatically recorded by a large weighing lysimeter, and the difference between the total weights of two adjacent collections is the actual evapotranspiration of the stage. The basic principle is shown in Equation (2):

$$\mathrm{P}+\mathrm{I}+{\mathrm{E}}_{\mathrm{g}}={\mathrm{P}}_{\mathrm{a}}+\mathrm{ET}+\mathrm{R}+\Delta \mathrm{S}$$

(2)

where P is the rainfall (mm); I is the irrigation water volume (mm); E_g is the diving evaporation (mm); P_a is the deep leakage rate (mm); ET is the actual evapotranspiration (mm); R is the volume of runoff (mm); and ΔS is the soil storage variable (mm).

Xu et al. [22] studied East China and concluded that the FAO-56 Penman–Monteith formula was the best method to calculate the daily reference evapotranspiration. The reference evapotranspiration is calculated by Equation (3):

$${\mathrm{ET}}_0=\frac{0.408\Delta \left({\mathrm{R}}_{\mathrm{n}}-\mathrm{G}\right)+\mathsf{\gamma }\frac{900}{\mathrm{T}+273}{\mathrm{u}}_2\left({\mathrm{e}}_{\mathrm{s}}-{\mathrm{e}}_{\mathrm{a}}\right)}{\Delta +\mathsf{\gamma }\left(1+0.34{\mathrm{u}}_2\right)}$$

(3)

where ET₀ is the reference evapotranspiration (mm × d⁻¹); R_n is the net surface radiation (MJ × m⁻² × d⁻¹); G is the soil heat flux (MJ × m⁻² × d⁻¹); T is the average daily temperature (°C); u₂ is the average wind speed at 2 m above the ground (m × s⁻¹); e_s is the saturated vapour pressure (kPa); e_a is the actual vapour pressure (kPa); Δ is the slope of the saturated vapour pressure and temperature curve (kPa × °C⁻¹); and γ is the dry and wet table constant (kPa × °C⁻¹). Data are from a high precision weather station.

The change curves of actual evapotranspiration and reference evapotranspiration are shown in Figure 1.

2.5. Crop Coefficient Estimation Method and Evaluation Indices

2.5.1. Temperature Effect Method

The temperature effect model proposed by Huang et al. [23] effectively simulated the dynamic process of crop growth and development. Wang et al. [16] used the model structure to construct the calculation formula of the crop coefficient considering the temperature of three basis points, as shown in Equation (4):

$${\mathrm{K}}_{\mathrm{c}}={\mathrm{K}}_0{\mathrm{e}}^{-{\left(\frac{\mathrm{T}-{\mathrm{T}}_0}{\mathsf{\beta }}\right)}^2}$$

(4)

where K_c is the crop coefficient; K₀ is the crop coefficient at the optimum temperature; T is the average temperature (°C); T₀ is the optimum temperature for physiological and ecological processes such as crop growth and photosynthesis (°C); and β is the parameter to be estimated.

Since an invalid temperature is not conducive to crop growth, the average temperature T is calculated using the eliminating invalid temperature method [24,25], where invalid temperature values above the upper limit and below the lower limit are discarded, as shown in Equation (5):

$$\mathrm{T}=\frac{{\mathrm{T}}_{\mathrm{x}}^{\prime }+{\mathrm{T}}_{\mathrm{n}}^{\prime }}{2}$$

(5)

The specific calculation method of each parameter in Equation (5) is

$${\mathrm{T}}_{\mathrm{x}}^{\prime }=\max \left({\mathrm{T}}_{\mathrm{x}},{\mathrm{T}}_{\mathrm{upper}}\right),{\mathrm{T}}_{\mathrm{n}}^{\prime }=\min \left({\mathrm{T}}_{\mathrm{n}},{\mathrm{T}}_{\mathrm{upper}}\right),{\mathrm{T}}_{\mathrm{x}}=\max \left({\mathrm{T}}_{\max },{\mathrm{T}}_{\mathrm{base}}\right),{\mathrm{T}}_{\mathrm{n}}=\max \left({\mathrm{T}}_{\min },{\mathrm{T}}_{\mathrm{base}}\right)$$

(6)

where T_upper is the upper limit temperature (°C), which is 30 °C and T_base indicates the lower limit temperature (°C), which is 3 °C [26,27,28].

The optimal values of unknown parameters K₀, T₀, and β were determined by SPSS software combined with the least squares method and sequential quadratic programming. First, the logarithm of Equation (4) can be obtained:

$$\ln \left({\mathrm{K}}_{\mathrm{c}}\right)=\ln \left({\mathrm{K}}_0\right)-\frac{{\mathrm{T}}_0^2}{{\mathsf{\beta }}^2}+\frac{2{\mathrm{T}}_0}{{\mathsf{\beta }}^2}\times \mathrm{T}-\frac{1}{{\mathsf{\beta }}^2}\times {\mathrm{T}}^2$$

(7)

Setting $\mathrm{y}=\ln \left({\mathrm{K}}_{\mathrm{c}}\right)$, $\mathrm{x}=\mathrm{T}$, ${\mathrm{a}}_1=\ln \left({\mathrm{K}}_0\right)-\frac{{\mathrm{T}}_0^2}{{\mathsf{\beta }}^2}$, ${\mathrm{a}}_2=\frac{2{\mathrm{T}}_0}{{\mathsf{\beta }}^2}$, and ${\mathrm{a}}_3=-\frac{1}{{\mathsf{\beta }}^2}$, Equation (7) is converted to

$$\mathrm{y}={\mathrm{a}}_1+{\mathrm{a}}_2\times \mathrm{x}+{\mathrm{a}}_3\times {\mathrm{x}}^2$$

(8)

Let $\mathrm{Y}=\left[\begin{array}{c}{\mathrm{y}}_1\\ ⋮\\ {\mathrm{y}}_{\mathrm{n}}\end{array}\right]$, $\mathrm{A}=\left[\begin{array}{c}{\mathrm{a}}_1\\ {\mathrm{a}}_2\\ {\mathrm{a}}_3\end{array}\right]$, and $\mathrm{X}=\left[\begin{array}{ccc}1& {\mathrm{x}}_{11}& {\mathrm{x}}_{21}\\ ⋮& ⋮& ⋮\\ 1& {\mathrm{x}}_{1\mathrm{n}}& {\mathrm{x}}_{2\mathrm{n}}\end{array}\right]$; if there are n samples in the whole growth process, then Equation (8) can be expressed as

$$\mathrm{Y}=\mathrm{AX}$$

(9)

The unknown parameter values in the whole growth process were estimated by combining the least square method:

$$\widehat{\mathrm{A}}={\left({\mathrm{X}}^{\mathrm{T}}\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{Y}$$

(10)

The inverse solution is:

$$\widehat{{\mathrm{K}}_0}={\mathrm{e}}^{\widehat{{\mathrm{a}}_0}-\frac{{\widehat{{\mathrm{a}}_1}}^2}{4-\widehat{{\mathrm{a}}_2}}},\widehat{{\mathrm{T}}_0}=-\frac{\widehat{{\mathrm{a}}_1}}{2\times \widehat{{\mathrm{a}}_2}},\widehat{\mathsf{\beta }}=\sqrt{-\frac{1}{\widehat{{\mathrm{a}}_2}}}$$

(11)

The obtained value of Equation (11) is set as the initial value, the sequential quadratic programming method is used to solve the optimal value of each parameter by SPSS software, and the objective function is set as

$$\min {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{\left(\widehat{{\mathrm{K}}_{\mathrm{c}}}-{\mathrm{K}}_{\mathrm{ci}}\right)}^2$$

(12)

The constraint conditions are ${\mathrm{T}}_{\mathrm{base}}<{\mathrm{T}}_0<{\mathrm{T}}_{\mathrm{upper}}$, $0<{\mathrm{K}}_0<3$, and $\mathsf{\beta }>0$.

2.5.2. Cumulative Crop Coefficient Method

Many scholars have used cumulative evapotranspiration [29] and cumulative growth days [30] when estimating crop coefficients and evapotranspiration; in this study, we attempted to use cumulative crop coefficients. The cumulative crop coefficient refers to the cumulative value of the crop coefficient since the day of planting and only considers the days after seeding. This method does not require data on hydrometeorological elements and crop growth traits and is convenient for use in areas where observation conditions and crop growth data are scarce.

The curves of the cumulative crop coefficient and days after seeding of winter wheat are shown in Figure 2, showing an approximately elongated “S” shape. The daily growth value of the cumulative crop coefficient is that of the daily crop coefficient, which increases first and then decreases with time after seeding. The cumulative crop coefficient model was constructed as follows:

$${\mathrm{K}}_{\mathrm{c}\mathrm{cumulative}}=\mathrm{a}+\mathrm{b}\times \cos \left(\mathrm{c}\times \mathrm{D}\right)+\mathrm{d}\times \sin \left(\mathrm{c}\times \mathrm{D}\right)$$

(13)

where K_{c cumulative} is the cumulative crop coefficient value, D is the number of days after seeding, and a, b, c, and d are unknown parameters.

For each growth stage, the least square method was used to estimate the values of unknown parameters a, b, c, and d. According to the definition of the function, its derivative is the daily value function of the crop coefficient, and the value of the crop coefficient on day D can be obtained by substituting D into it. The derivative of Equation (13) with respect to D can be obtained as follows:

$${\mathrm{K}}_{\mathrm{c}}=-\mathrm{b}\times \mathrm{c}\times \sin \left(\mathrm{c}\times \mathrm{D}\right)+\mathrm{d}\times \mathrm{c}\times \cos \left(\mathrm{c}\times \mathrm{D}\right)$$

(14)

2.5.3. Radiative Soil Temperature Method

According to the study of Zhao et al. [31], there is a strong positive correlation between the crop coefficient and evapotranspiration and a significant positive correlation between evapotranspiration and effective energy. Effective energy (R_n-G) refers to the difference between net radiation and soil heat flux, showing a trend of being smaller at the emergence-branching and branch-overwintering stages and gradually increasing at the greening-jointing and heading-maturity stages (Figure 3). During the whole process of winter wheat growth, the correlation coefficient between the crop coefficient and effective energy reached 0.75, showing a strong correlation. The correlation coefficient between the crop coefficient and soil temperature from 0 to 160 cm decreased from 0.66 to 0.70 with increasing depth, as shown in Figure 4.

Because FAO-56 notes that the penetration depth of temperature waves in soil is 0.1 to 0.2 m at a daily scale, soil temperatures at 0 and 10 cm were selected. The following mathematical model was constructed involving radiation and soil temperature:

$${\mathrm{K}}_{\mathrm{c}}={\mathrm{me}}^{\left({\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-\mathrm{n}\left|\frac{{\mathrm{D}}_0}{{\mathrm{D}}_{10}}\right|\right)}$$

(15)

where D₀ is the surface temperature (°C); D₁₀ is the soil temperature of 10 cm (°C); and m and n are unknown parameters. The meanings of the other symbols are the same as described previously.

Parameter calibration was realized by the particle swarm optimization algorithm. After finding the unknown parameter values corresponding to the daily scale data of each group in the whole growth process, the best unknown parameter of each growth stage was obtained by taking the average value. The particle swarm optimization algorithm is a random optimization method based on swarm intelligence proposed by Kennedy and Eberhart [32]. Each particle represents a candidate solution, and the problem is solved through assessment of the simple behaviour of individual particles and information interaction within the group.

When solving the optimization problem, each particle has two state quantities, position and velocity, and the fitness value determined by the objective function. The flight process of the particle is the search process of the individual. In each iteration, individual particles record the best solution found as the current individual extreme value and share it with other particles. The best of all individual extreme values is the extreme value of the current group. All particles adjust their speed and position according to their current individual extreme value and the current group extreme value of the whole particle swarm [33]. Optimization was achieved using MATLAB software programming.

According to the properties of the model, in the two-dimensional solution search space, there is a particle swarm composed of 500 particles. The position of the particle swarm of the k-th iteration in the solution space is assumed to be expressed as Equation (16), and the flight velocity is expressed as Equation (17):

$${\mathrm{X}}_{\mathrm{i}}^{\mathrm{k}}=\left({\mathrm{X}}_{\mathrm{i}1}^{\mathrm{k}},{\mathrm{X}}_{\mathrm{i}2}^{\mathrm{k}}\right),\mathrm{i}=1,2,\dots ,500$$

(16)

$${\mathrm{V}}_{\mathrm{i}}^{\mathrm{k}}=\left({\mathrm{V}}_{\mathrm{i}1}^{\mathrm{k}},{\mathrm{V}}_{\mathrm{i}2}^{\mathrm{k}}\right),\mathrm{i}=1,2,\dots ,500$$

(17)

At this time, the best individual historical position of the jth (j < 500) particle is $\left({\mathrm{X}}_{\mathrm{mj}1}^{\mathrm{k}},{\mathrm{X}}_{\mathrm{mj}2}^{\mathrm{k}}\right)$, and the best historical position of the particle swarm is $\left({\mathrm{Y}}_{\mathrm{m}1}^{\mathrm{k}},{\mathrm{Y}}_{\mathrm{m}2}^{\mathrm{k}}\right)$. Then, at the k + 1 iteration, the velocity and position update formula of the particle is as follows:

$$\left\{\begin{array}{c}{\mathrm{V}}_{\mathrm{j}1}^{\mathrm{k}+1}=\mathrm{w}\times {\mathrm{V}}_{\mathrm{j}1}^{\mathrm{k}}+{\mathrm{c}}_1\times {\mathrm{r}}_1\times \left({\mathrm{X}}_{\mathrm{mj}1}^{\mathrm{k}}-{\mathrm{X}}_{\mathrm{j}1}^{\mathrm{k}}\right)+{\mathrm{c}}_2\times {\mathrm{r}}_2\times \left({\mathrm{Y}}_{\mathrm{m}1}^{\mathrm{k}}-{\mathrm{X}}_{\mathrm{j}1}^{\mathrm{k}}\right)\\ {\mathrm{V}}_{\mathrm{j}2}^{\mathrm{k}+1}=\mathrm{w}\times {\mathrm{V}}_{\mathrm{j}2}^{\mathrm{k}}+{\mathrm{c}}_1\times {\mathrm{r}}_1\times \left({\mathrm{X}}_{\mathrm{mj}2}^{\mathrm{k}}-{\mathrm{X}}_{\mathrm{j}2}^{\mathrm{k}}\right)+{\mathrm{c}}_2\times {\mathrm{r}}_2\times \left({\mathrm{Y}}_{\mathrm{m}2}^{\mathrm{k}}-{\mathrm{X}}_{\mathrm{j}2}^{\mathrm{k}}\right)\end{array}\right.$$

(18)

$$\left\{\begin{array}{c}{\mathrm{X}}_{\mathrm{i}1}^{\mathrm{k}+1}={\mathrm{X}}_{\mathrm{i}1}^{\mathrm{k}}+{\mathrm{V}}_{\mathrm{j}1}^{\mathrm{k}+1}\\ {\mathrm{X}}_{\mathrm{i}2}^{\mathrm{k}+1}={\mathrm{X}}_{\mathrm{i}2}^{\mathrm{k}}+{\mathrm{V}}_{\mathrm{j}2}^{\mathrm{k}+1}\end{array}\right.$$

(19)

$$\mathrm{w}={\mathrm{w}}_{\min }+\left({\mathrm{w}}_{\max }-{\mathrm{w}}_{\min }\right)\times \frac{\mathrm{ger}-\mathrm{times}}{\mathrm{ger}}$$

(20)

where $\mathrm{w}$ is the inertia weight, representing the influence of the current speed on the next movement; ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are the learning factors, which measure the influence of the current individual historical optimal position and group historical optimal position on the next movement; ${\mathrm{r}}_1$ and ${\mathrm{r}}_2$ are random numbers from 0 to 1; ${\mathrm{w}}_{\min }$ is the initial inertia weight with a typical value of 0.4; ${\mathrm{w}}_{\max }$ is the inertia weight under the maximum number of iterations with a typical value of 0.9; $\mathrm{ger}$ is the maximum number of iterations; and $\mathrm{times}$ indicates the current iteration times.

2.5.4. Indices of Evaluation

The evaluation indices include the root mean square error (RMSE), mean absolute error (MAE), correlation coefficient (r), and consistency index (d_IA), which are used to evaluate the error and consistency between the estimated value and the measured value of each estimation method. See Equation (21) to Equation (24) for the calculation formula of each index.

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{x}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{\mathrm{n}}}$$

(21)

$$\mathrm{MAE}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{x}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}\right|}{\mathrm{n}}$$

(22)

$$\mathrm{r}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}{\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)}^2}\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}^2}}$$

(23)

$${\mathrm{d}}_{\mathrm{IA}}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{x}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left(\left|{\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{y}}\right|+\left|{\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right|\right)}^2}$$

(24)

where ${\mathrm{x}}_{\mathrm{i}}$ is the estimated value of ${\mathrm{K}}_{\mathrm{c}}$; ${\mathrm{y}}_{\mathrm{i}}$ is the actual value of ${\mathrm{K}}_{\mathrm{c}}$; $\mathrm{i}$ is the sample ordinal number, $\mathrm{i}=1,2,\dots ,\mathrm{n}$; $\overline{\mathrm{x}}$ is the mean value of ${\mathrm{K}}_{\mathrm{c}}$’s estimation; $\overline{\mathrm{y}}$ is the mean value of ${\mathrm{K}}_{\mathrm{c}}$; and n is the number of samples of the estimated value.

Generally, the closer the root mean square error RMSE and mean absolute error MAE are to 0, the smaller the error and the greater the accuracy. The closer the correlation coefficient r and the consistency index d_IA are to 1, the closer the estimated value of the model is to the actual value, and the stronger its estimation ability is.

3. Results

Based on the observation data of the large weighing lysimeter and the data of the high-precision weather station from 2018 to 2019, the parameters of the three crop coefficient estimation models were calibrated. The unknown parameter values obtained by each method under different growth stages are shown in

Table 2. Values of unknown parameters in different growth stages and methods of winter wheat.

Method	Temperature Effect			Cumulative Crop Coefficient				Radiative Soil Temperature
Parameter	${\mathit{K}}_{\mathbf{0}}$	${\mathit{T}}_{\mathbf{0}}$	$\mathit{\beta }$	$\mathit{a}$	$\mathit{b}$	$\mathit{c}$	$\mathit{d}$	$\mathit{m}$	$\mathit{n}$
Emergence-branching stage	1.24	3.00	18.84	158.62	−7.17	0.07	6.54	4.19	5.87
Branch-overwintering stage	1.96	3.00	12.42	216.22	−79.90	0.01	60.00	0.26	3.21
Greening-jointing stage	2.14	20.95	6.61	180.23	−382.33	−0.01	26.25	0.41	5.62
Heading-maturity stage	2.39	20.37	5.16	227.91	0.02	−0.05	43.52	6.16	8.81

.

3.1. The Differences and Causes of Crop Coefficient Estimation by Different Methods

The results of the three methods for estimating winter wheat crop coefficients are shown in Figure 5 and Figure 6. The results of the temperature effect method were better in the emergence-branching and heading-maturity stages, followed by the greening-jointing and branch-overwintering stages. The reason was that there were fewer ineffective temperatures below 3 °C and above 30 °C during the emergence-branching and heading-maturity stages, and crop growth was not inhibited by ineffective temperatures. The cumulative crop coefficient method was better in the greening-jointing and heading-maturity stages but worse in the emergence-branching and branch-overwintering stages. The reason was that the cumulative crop coefficient method is closely related to crop growth [17], the leaf index in the greening-jointing and heading-maturity stages was larger [34,35], and crop growth and development were vigorous. The radiative soil temperature method was better in the emergence-branching and branch-overwintering stages but worse in the greening-jointing and heading-maturity stages. The reason was that the effective energy cannot maintain a high level during the greening-jointing and heading-maturity stages, and the balance between energy and soil temperature was difficult to ensure.

3.2. Determination of The Best Estimation Method for Each Growth Stage

The accuracy of the results of the three methods in each growth stage are shown in

Table 3. The estimated precision index value of the results of three methods in each growth stage.

Stage	Emergence-Branching Stage				Branch-Overwintering Stage				Greening-Jointing Stage				Heading-Maturity Stage
Index	RMSE	MAE	r	dIA	RMSE	MAE	r	dIA	RMSE	MAE	r	dIA	RMSE	MAE	r	dIA
TE	0.06	0.06	0.80	0.88	0.13	0.11	0.44	0.55	0.23	0.18	0.70	0.83	0.16	0.13	0.94	0.97
CCC	0.08	0.07	0.57	0.69	0.13	0.12	0.36	0.51	0.16	0.13	0.88	0.89	0.20	0.16	0.91	0.94
RST	0.20	0.19	0.35	0.51	0.25	0.22	0.52	0.61	0.93	0.79	0.70	0.49	1.10	0.91	0.43	0.49

(in this table, “TE” means the temperature effect method, “CCC” means the cumulative crop coefficient method, and “RST” means the radiative soil temperature method).

The estimation results were compared, and the best estimation method was selected according to the following steps. In the first step, the root mean square error, RMSE, and mean absolute error, MAE, were compared, and the method with the smallest error was selected. In the second step, when the difference in the root mean square error and mean absolute error between multiple methods are no more than 0.03, the correlation coefficient r and consistency index d_IA are compared in turn to select the method with the strongest estimation ability. In the third step, if the difference between the correlation coefficient and consistency index is less than 0.03, then the cumulative crop coefficient method should be directly adopted according to the principle of minimum required observations if it is still in the option list. Otherwise, the most appropriate method should be adopted in combination with the graph.

According to the above criteria, the temperature effect method during the emergence-branching stage was the best (RMSE = 0.06, MAE = 0.06, r = 0.80, d_IA = 0.88), followed by the cumulative crop coefficient method, and the radiative soil temperature method was inferior. The temperature effect method during the branch-overwintering stage was the best (RMSE = 0.13, MAE = 0.11, r = 0.44, d_IA = 0.55), followed by the cumulative crop coefficient method, and the radiative soil temperature method was inferior. The cumulative crop coefficient method during the greening-jointing stage was the best (RMSE = 0.16, MAE = 0.13, r = 0.88, d_IA = 0.89), followed by the temperature effect method, and the radiative soil temperature method was inferior. The temperature effect method during the heading-maturity stage was the best (RMSE = 0.16, MAE = 0.13, r = 0.94, d_IA = 0.97), followed by the cumulative crop coefficient method, and the radiative soil temperature method was inferior.

The estimated accuracy index of the whole growth process when a single method or the best method is used in each growth stage is shown in

Table 4. The estimation precision index values of the whole growth process.

Method	Root Mean Square Error	Mean Absolute Error	Correlation Coefficient	Consistency Index
TE	0.34	0.25	0.87	0.93
CCC	0.25	0.20	0.93	0.96
RST	0.79	0.58	0.50	0.69
The best	0.13	0.09	0.98	0.99

. When a single method was used in each growth stage, the cumulative crop coefficient method was the best method for estimating the whole growth process (RMSE = 0.25, MAE = 0.20, r = 0.93, d_IA = 0.96), followed by the temperature effect method, and the radiative soil temperature method was inferior. When the cumulative crop coefficient method and temperature effect method were used, the correlation coefficient and consistency index of the whole growth process estimation results were both greater than 0.80, which met the accuracy requirements of estimation and could be used for crop coefficient estimation. When the radiative soil temperature method was used, the correlation coefficient and consistency index of the whole growth process estimation results were only 0.50 and 0.69, respectively, which could not meet the requirement of estimation accuracy.

When the best method was used in each growth stage, the four precision indices were better than when a single method was used; the root mean square error and mean absolute error were 0.13 and 0.09, and the correlation coefficient and consistency index reached 0.98 and 0.99, far greater than 0.80. The best method has higher estimation ability and accuracy than the single method and can be used for crop coefficient estimation.

4. Conclusions and Discussion

1. The results of the temperature effect method were better in the emergence-branching and heading-maturity stages, followed by the greening-jointing and branch-overwintering stages. The reason was that the ineffective temperature was lower during the emergence-branching and heading-maturity stages, and crop growth was not inhibited. The cumulative crop coefficient method was better in the greening-jointing and heading-maturity stages but worse in the emergence-branching and branch-overwintering stages. The reason was that the cumulative crop coefficient method is closely related to crop growth, and the crops grew vigorously during the greening-jointing and heading-maturity stages. The radiative soil temperature method was better in the emergence-branching and branch-overwintering stages but worse in the greening-jointing and heading-maturity stages. The reason was that the effective energy could not maintain a high level during the greening-jointing and heading-maturity stages, and the energy balance was difficult to ensure.

2. The temperature effect method during the emergence-branching stage was the best, followed by the cumulative crop coefficient method, and the radiative soil temperature method was inferior. The temperature effect method during the branch-overwintering stage was the best, followed by the cumulative crop coefficient method, and the radiative soil temperature method was inferior. The cumulative crop coefficient method during the greening-jointing stage was the best, followed by the temperature effect method, and the radiative soil temperature method was inferior. The temperature effect method during the heading-maturity stage was the best, followed by the cumulative crop coefficient method, and the radiative soil temperature method was inferior.

3. When a single method was used in each growth stage, the cumulative crop coefficient method was the best, followed by the temperature effect method, and the radiative soil temperature method was inferior. The cumulative crop coefficient method and temperature effect method meet the accuracy requirements of estimation, but the radiative soil temperature method could not meet the accuracy requirements of estimation. The root mean square error, RMSE = 0.13; mean absolute error, MAE = 0.09; correlation coefficient, r = 0.98; and consistency index, dIA = 0.99 were all better than the single method when the best method was used in each growth stage. The best method had higher estimation ability and accuracy than the single method.

To improve upon the results presented in this paper, the adaptability of different crop coefficient estimation methods in other areas needs to be further explored. In addition, the daily crop coefficient estimation is only discussed under the condition of 1 m burial depth. The difference in the estimation effect of different methods under different burial depths needs to be further studied. In addition to crop growth traits, there is also the standard crop coefficient stage division table recommended by the FAO. Different division methods also have a certain impact on the accuracy of estimation. In the future, more crop coefficient estimation methods can be considered to provide more options for crop coefficient estimation.