Drought Risk Assessment of Winter Wheat at Different Growth Stages in Huang-Huai-Hai Plain Based on Nonstationary Standardized Precipitation Evapotranspiration Index and Crop Coefficient

Abstract

Soil moisture plays a crucial role in determining the yield of winter wheat. The Huang-Huai-Hai (HHH) Plain is the main growing area of winter wheat in China, and frequent occurrence of drought seriously restricts regional agricultural development. Hence, a daily-scale Non-stationary Standardized Precipitation Evapotranspiration Index (NSPEI), based on winter wheat crop coefficient (Kc), was developed in the present study to evaluate the impact of drought characteristics on winter wheat in different growth stages. Results showed that the water demand for winter wheat decreased with the increase in latitude, and the water shortage was affected by effective precipitation, showing a decreasing trend from the middle to both sides in the HHH Plain. Water demand and water shortage showed an increasing trend at the jointing stage and heading stage, while other growth stages showed a decreasing trend. The spatial distributions of drought duration and intensity were consistent, which were higher in the northern region than in the southern region. Moreover, the water shortage and drought intensity at the jointing stage and heading stage showed an increasing trend. The drought had the greatest impact on winter wheat yield at the tillering stage, jointing stage, and heading stage, and the proportions of drought risk vulnerability in these three stages accounted for 0.25, 0.21, and 0.19, respectively. The high-value areas of winter wheat loss due to drought were mainly distributed in the northeast and south-central regions.

1. Introduction

With global warming, the water cycle is accelerating and extreme weather events such as droughts and floods are becoming more frequent and intense [1]. Droughts are triggered by prolonged water scarcity, with severe agricultural and socio-economic impacts [2]. China is a country vulnerable to agricultural drought, with an annual loss of about 30 billion kilograms of all grain due to drought [3]. The Huang-Huai-Hai Plain (HHH) is one of the most important agricultural production areas in China, accounting for 36% and 37% of the country’s arable land and grain output [4]. However, the winter wheat growth period in this region is a water-deficit period, and the water expenditure during the wheat growth period is greater than the water input. Therefore, drought monitoring during the growth period of winter wheat is important for drought prevention and resistance for winter wheat.

The shortage of water resources has become one of the key factors restricting the sustainable development of agriculture. Crop coefficient and crop water demand are important parameters in water-saving irrigation of winter wheat in this region [5]. Vicente proposed to estimate the crop water loss due to evapotranspiration, which is of great significance for the rational allocation of irrigation water and the improvement of water use efficiency [6]. Mirgol used evapotranspiration and crop coefficients to estimate crop water use and analyze the impact of climate change on crop yield [7]. Zhang proposed that effective precipitation and crop water demand reflected the supply and demand of agricultural water under precipitation conditions [8]. The joint risk of effective precipitation and crop water demand is the basis for regional water allocation and irrigation planning. Recep found that drought at different growth stages had different effects on crop phenology and yield, and crop yield loss was not only related to stress intensity, but also to the growth stage of crops [9]. Hence, both crop coefficient and crop water requirement were important parameters for evaluating crop yield. However, agricultural drought is caused by the imbalance between crop water demand and supply. Potential evapotranspiration was used instead of crop water demand in previous drought indices, which limited the applicability of evaluating agricultural drought. Therefore, it is necessary to consider the crop coefficient Kc in the calculation of the drought index, and replace potential evapotranspiration with the crop’s water demand, to improve the monitoring of the drought situation for the crop.

In previous studies, several drought indices were constructed to characterize the duration and intensity of drought. The most widely used meteorological drought indicators are the Palmer Drought Index (PDSI) [10], the Standardized Precipitation Index (SPI) [11], and the Standardized Precipitation Evapotranspiration Index (SPEI). Vicente-Serrano et al. (2010) added potential evapotranspiration to the SPI and constructed SPEI for drought monitoring and evaluation that integrated the joint effects of precipitation and temperature changes [12]. However, previous drought indices were calculated under the assumption of stationarity, which is questioned in changing environments. Under global warming, the assumption of stationarity of factors, such as future precipitation and temperature, is being challenged [13]. Many studies have found that hydrometeorological variables change linearly or nonlinearly with time. Therefore, it is necessary to construct a non-stationary drought index to identify droughts under non-stationary conditions [14,15,16,17,18,19,20].

Crops are more sensitive to drought response during critical growing seasons, and several days of drought can cause them to wilt [21,22]. Previous studies generally underestimated the intensity of short-term continuous strong droughts and overestimated the intensity of long-term weak droughts when monitoring drought intensity [23,24]. Therefore, identifying short-duration drought processes is necessary to monitor crop drought. In addition, the daily-scale SPEI captures the drought process in more detail than the monthly scale SPEI [25,26,27,28,29]. The previous meteorological drought index seldom considered the crop coefficient and effective precipitation which have a greater impact on the actual drought process of crops, which made it difficult for drought indicators to monitor the actual drought process of crops.

Hence, in this study, the effective precipitation closer to the crop irrigation water deficit was used to calculate the crop water deficit. Then, based on the water requirement and effective precipitation of winter wheat, a daily NSPEI drought index (NSPEI-Kc) considering the crop coefficient was constructed to realize real-time dynamic monitoring of the regional agricultural drought. Besides, the effect of crop coefficients in different growth stages on the drought index was comprehensively considered, and a machine learning model was used to evaluate the quantitative impact of drought on different growth stages of crop. Finally, the drought risk zoning analysis of winter wheat was carried out through the yield disaster loss risk index. This is conducive to improving the refined monitoring and drought risk assessment of winter wheat drought in the HHH Plain.

2. Data and Methods

2.1. Study Area

Huang-Huai-Hai Plain, the study area (Figure 1), is located in the east-central part of China, between 32°00′~40°24′N, 112°48′~120°45′E. The topography of the region is mainly plain, the land is flat, and the soil is mostly loam and sandy loam. Located south of Yanshan Mountain and north of Huaihe River, it is an alluvial plain of Yellow River, Huaihe River, and Haihe River and part of a hilly mountainous area. This region is a semi-arid and semi-humid region, with annual precipitation of 500–900 mm, high annual variability of precipitation, and uneven seasonal distribution. A total of 60% of the precipitation is concentrated in summer [30], but the growth period of wheat (autumn, winter, and spring) is a dry period of water deficit, which aggravates the drought degree of the winter wheat reproductive period, and the water deficit of winter wheat growth period is 180–350 mm.

2.2. Data

The meteorological data used in this paper were obtained from the National Climate Center of the China Meteorological Administration. In order to ensure the integrity and continuity of the data series, stations with missing data greater than 1% were excluded from analysis. A total of 280 meteorological stations in the Huang-Huai-Hai Plain from 1970 to 2019 were selected for analysis. Meteorological variables included precipitation and temperature, wind speed, air pressure, sunshine duration, relative humidity, and other variables which were needed to calculate potential evapotranspiration. The crop coefficient (Kc) was calculated using the single value average method of FAO segment which is a relatively simple method and does not require a large amount of data and can be widely used to calculate the crop water requirement [31]. However, this method requires a longer time period and does not adequately consider the effect of soil. The conditions in most areas of the Huang-Huai-Hai Plain cannot meet certain standard conditions, so the corrected crop coefficient value of Liu Jia et al. [32] was adopted in this study, as shown in

Table 1. Growth time of winter wheat and Kc value of crop coefficient.

	Sowing	Tillering	Jointing	Heading	Maturing
Start-stop time	9.15–11.30	12.01–3.12	3.13–4.10	4.11–5.3	5.4–5.25
Duration	77	102	29	23	22
Kc value	1.164	1.1555	1.115	1.0145	0.865

. Liu Jia’s revised crop coefficient was determined, based on the research results of the wheat water demand contour map in Anhui Province, the crop coefficient and correction formula of 84 crops recommended by FAO, and the irrigation experimental data of Xinmaqiao Agricultural and Water Comprehensive Experiment Station of Anhui Water Resources Research Institute of recent years. The winter wheat yield data used were unit annual yields (kg/ha) from 2000 to 2019 for each province within the Huang-Huai-Hai Plain and were taken from the provincial statistical yearbooks.

2.3. Methods

The framework shown in Figure 2 describes the overall process of constructing the daily NSPEI-Kc index of winter wheat based on weather station data, crop coefficient, and yield data, and evaluating the drought risk of winter wheat at different growth stages.

2.3.1. Calculation of Required Water Shortage

In the first step, the Penman–Monteith method recommended by 1988FAO was used to calculate daily potential evapotranspiration ET₀, which was calculated as follows [31]:

$$E{T}_0=\frac{0.408\mathsf{\Delta }\left(R_n-G\right)+\gamma \frac{900}{T+273}{u}_2\left(e_s-e_a\right)}{\mathsf{\Delta }+\gamma \left(1+0.34U_2\right)}$$

(1)

where ET₀ is the potential evapotranspiration (mm·d⁻¹); R_n is the net radiation (MJ·m⁻²·d⁻¹); G is the soil heat flux; T is the mean daily air temperature (°C); u₂ is the wind speed at 2 a m height (m·s⁻¹); e_s is the saturated water vapor pressure (kPa); e_a is the actual water vapor pressure (kPa); ∆ is the slope of saturated water vapor pressure–temperature curve (kPa); and γ is dry and wet table constant (kPa/°C).

In the second step, the daily potential evapotranspiration was multiplied by the Kc factor for each reproductive period to obtain the actual water requirement:

$$E{T}_c=E{T}_0*Kc$$

(2)

where ET_c is the actual evapotranspiration of the crop (mm/d), and Kc is the crop coefficient.

The third step was to calculate the effective precipitation which is the amount of rainfall that can be consumed to meet crop evaporation in arid areas, excluding runoff, seepage to areas beyond the root zone, and the deep seepage portion required for leaching salts, so it was calculated using the method of effective precipitation recommended by the USDA Soil Conservation Service, whose expression is as follows [33]:

$${\mathrm{P}}_{\mathrm{e}}=\{\begin{array}{l}\mathrm{P}\left(4.17-0.2\mathrm{P}\right)/4.17\hspace{1em}\left(\mathrm{P}<8.3\mathrm{mm}/\mathrm{d}\right)\\ 4.17+0.1\mathrm{P}\hspace{1em}\left(\mathrm{P}>=8.3\mathrm{mm}/\mathrm{d}\right)\end{array}$$

(3)

where P_e is the effective daily precipitation (mm/d) and P is the total daily precipitation (mm/d).

Step four is to calculate the water deficit, also known as irrigation water demand, which is the difference between water demand and effective precipitation during each reproductive period of winter wheat and was used to indicate the surplus or deficit of water supply. A positive deficit indicated a water deficit, while a negative deficit indicated a water surplus.

$$W=E{T}_c-P_{e }$$

(4)

where W is the crop water deficit during the reproductive period, ET_c is the crop water requirement during the reproductive period, and P_e is the effective precipitation during the reproductive period.

2.3.2. Construction of Non-Stationary Precipitation Evapotranspiration Index (NSPEI-Kc) Based on Kc

In this study, the R language was used to construct the daily NSPEI index. The improved daily NSPEI drought index calculation process was similar to the monthly NSPEI. The specific calculation process is as follows:

The first step was to calculate the difference, D_t, between daily effective precipitation, P_e, and actual daily evapotranspiration (PET_t):

$$D_t=P_e-PE{T}_t$$

(5)

D_t was fitted by Smoothing Splines to determine the linear or nonlinear trend fitting of D.

$${\mathrm{SS}}_{D_t}(h)={\displaystyle \sum }_{\mathrm{t}=1}^n{\left[D_t-f\left(D_t\right)\right]}^2+h{{\displaystyle \int }}_{t_{\min }}^{t_{\max }}{\left[f^n\left(D_t\right)\right]}^2\mathrm{d}t$$

(6)

where f is the linear fitting function of time series D_t; h is a smooth parameter, similar to the amplitude estimated by a local polynomial; t is the time series, t = 1, …, n, t_max is the highest temperature; and t_min indicates the minimum temperature.

The position parameters of log-logistic distribution were replaced by the trend values fitted with smooth splines. The shape and scale parameters of the distribution were unchanged, and only the influence of time changes on time series D was considered. Therefore, information such as extreme values and data distribution of time series D was not lost. The advantage of the smooth spline function was that it combined the characteristics of the data itself, and it did not need to choose a linear or nonlinear model to fit the D time series, nor did it need to judge the stationarity of the D time series.

According to Formulas (5) and (6), the parameters of a time-varying position were obtained as follows:

$${\gamma }_{D_t}=loess\left(D_t\right)$$

(7)

Time-varying D_t-based distribution function:

$$f\left(D_t\mid \alpha ,\beta ,{\gamma }_{D_t}\right)=\frac{a}{\beta }{\left(\frac{D_t-{\gamma }_{D_t}}{\alpha }\right)}^{\beta -1}{\left[1+{\left(\frac{D_t-{\gamma }_{D_t}}{\alpha }\right)}^{\beta }\right]}^{-2}$$

(8)

where α, β, γ are the scale, shape, and position parameters, respectively. In this paper, the probability weighted moments method (PWM_s) with empirical frequencies was used to estimate the parameters α, β, γ. The PWM_s were calculated as follows:

$$w_s=\frac{1}{N}{\displaystyle \sum }_{t=1}^N\frac{{(N-t+0.35)}^s{D}_t}{N}$$

(9)

where w_s is the PWN of order s, where s = 4; N is the length of the data.

$$F(x)={{\displaystyle \int }}_0^{x}f\left(D_t\mid \alpha ,\beta ,{\gamma }_{D_t}\right)dt{\left[1+{\left(\frac{\alpha }{D_t-{\gamma }_{D_t}}\right)}^{\beta }\right]}^{-1}$$

(10)

$$\mathrm{NSPEI}=W-\frac{C_0+C_1+C_2{W}^2}{1+d_{1}W+d_2{W}^2+d_3{W}^3},W=-2\ln (P)$$

(11)

where F is the frequency estimate, when P ≤ 0.5, P is the cumulative probability, P = 1 − F(x); when P ≥ 0.5, then P = 1 − P. The other parameters were as follows: C₀ = 2.515517, C₁ = 0.80285, C₂ = 0.01028, d₁ = 1.432788, d₂ = 0.189269, d₃ = 0.001308.

Table 2. Classification of the NSPEI-Kc index corresponding to drought grade and yield reduction rate.

Drought Grade	NSPEI-Kc Index	Yield Reduction Rate
Mild drought	−0.99 ≤ NSPEI-Kc < 0	10 ≤ yd < 0
Moderate drought	−1.49 ≤ NSPEI-Kc < −1	20 ≤ yd < 10
Severe drought	−1.99 ≤ NSPEI-Kc < −1.5	30 ≤ yd < 20
Extreme drought	NSPEI-Kc ≤ −2	yd ≤ 30

presents the classification table for NSPEI-Kc.

2.3.3. Calculation of Drought Characteristics

The drought process was considered to have occurred when NSPEI-Kc had been in the mild drought category or above for 10 consecutive days. The beginning date of the drought process was the date when the NSPEI-Kc index reached mild drought or above on the first day. During the drought period, the drought was relieved when the NSPEI-Kc index reached no drought grade for 10 consecutive days, while the end date of the drought process was the date when the drought grade reached for the last time.

The run course theory is usually used to identify drought events and calculate drought characteristic variables [34]. The run course theory was used to identify drought duration (D) and drought intensity (I) from the NSPEI-Kc index of the daily scale. As shown in Figure 3, the time from the beginning to the end of the drought process was the duration of the drought (D), and the drought intensity (I) was the sum of drought categories whose NSPEI-Kc index was mild drought or above during all days of the drought process.

NSPEI-Kc at the jointing stage, heading stage, and maturity stage was used at the 30-days (1 month) time scale [35]. Due to the effects of precipitation and evapotranspiration in the first 90 days (about 3 months) during the planting stage in Huang-Huai-Hai Plain [36], the NSPEI-Kc at the 90 days (about 3 months) scale was used during the sowing stage, tillering stage, and whole growth stage.

2.3.4. Mann–Kendall (M-K) Test

The M-K test is a climate diagnostic and prediction technique that is widely used as a statistical test for analyzing trends in climate and hydrologic time series, which avoids local maxima in the data series [37]. The M-K test is also often used to detect trends in precipitation and drought frequency under the influence of climate change. A positive value of the M-K statistic Z A positive value indicates an upward trend and vice versa indicates a downward trend. If Z > 1.96, it indicates a significant trend of change.

2.3.5. Winter Wheat Yield Reduction Rate

There are many factors affecting the formation of the crop yield, and the relationship between the constraints is also very complex [38,39]. Crop yields are usually made up of three components: trend, climate fluctuation, and randomness. The variation trend of output reflects the contribution of factors affecting product development, and the climatic fluctuation component of output is the contribution of climate fluctuation, which is mainly related to meteorological disasters. The key to studying the relationship between drought intensity and winter wheat yield is to isolate this trend from the effects of climate fluctuations. The calculation was as follows:

$${\mathrm{y}=\mathrm{y}}_{\mathrm{t}}{+\mathrm{y}}_{\mathrm{w}}+\Delta \mathrm{y}$$

(12)

where y is the crop yield, y_t is the trend yield, and y_w is the meteorological yield. The “noise” ∆y term accounts for a small percentage and is often ignored and disregarded in practical calculations. This study used HP filtering to eliminate trends in winter wheat yield [40]. HP filtration has been widely used in hydrology, meteorology, and other fields. HP filtering is a complex signal separation method that captures not only the trend of a time series, but also the time variation of the series [41].

The winter wheat yield reduction rate was assigned using year-by-year deviations of actual yield from trend yield relative to meteorological yield, and the final formula was calculated as follows:

$${\mathrm{y}}_{\mathrm{d}}=\frac{{\mathrm{y} - \mathrm{y}}_{\mathrm{t}}}{{\mathrm{y}}_{\mathrm{t}}} \times 100\%$$

(13)

In Equation (13), y_d is the crop yield reduction rate; y is the actual grain yield; y_t is the trend grain yield in kg/hm². The yield reduction rate was classified as mild, moderate, severe, and very severe drought according to the drought level in Table 2.

2.3.6. Drought Risk Index

Due to significant geographical differences, the extent of drought impact on yield varies from region to region, with more significant differences in the vulnerability, resilience, and proportional impact of drought on yield across fertility periods. Disaster risk is the result of the interaction between the intensity and frequency of the causative event and the carrier sensitive to that type of causative factor. In this paper, a risk index reflecting four risk elements was used as a partitioning indicator, which was related to the vulnerability index, disaster loss rate, resilience, and weight of drought impact on yield in different fertility periods. The risk assessment from the perspective of winter wheat yield disaster loss was carried out with the following equation:

$$R=\left(S_P/S_C\right)\times \left(1-X/X_{\mathrm{m}}\right)\times \left({\displaystyle \sum _{i=1}^{\mathrm{m}}{Y}_i/P_i}\right)\times \left(1/\left(K/S_P\right)\right)\times G$$

(14)

In Equation (14), where R is the yield disaster risk index; S_P is the sown area; S_c is the total cultivated area; X is the average yield of wheat; X_m is the maximum yield; Y_i is the multi-year average yield reduction rate; P_i is the frequency of drought; K is the effective irrigated area; and G is the weight of drought affecting yield in different growth stages.

2.3.7. Random Forest (RF)

RF was based on Classification and Regression Trees (CART), which generated several independent trees by two stochastic methods of selecting training samples and selecting variables at each node of the tree to reach the final decision [42]. This randomness alleviated the typical drawbacks of CART, such as overfitting problems and sensitivity to training sample configurations [43].

2.3.8. Support Vector Machine (SVM)

Support vector machine (SVM) is a machine learning method, which is suitable for problems with small sample sizes, nonlinear relationship between variables, and multi-dimensional pattern recognition [44]. Support vector classification and regression relies on the theory of statistical learning and classifies based on the principle of structural risk minimization. Because of its powerful classification and regression ability, SVM has been widely used in image classification, handwriting recognition and remote sensing.

2.3.9. Linear Regression (LR)

Linear regression (LR) can be divided into simple linear regression (single linear regression) or multiple linear regression (multiple linear regression), depending on the number of independent variables. A linear regression model is a common linear regression model used to predict the linear relationship between a continuous target variable and one or more independent variables. The goal of the model is to find the best line so that the error between the predicted result and the actual result is minimal.

3. Result

3.1. Applicability Evaluation of Daily NSPEI-Kc

In this study, the daily NSPEI and daily NSPEI-Kc indices of the Huang-Huai-Hai Plain were input into each machine model and evaluated in combination with the actual yield reduction rate. Drought intensity and trend yield in different growth periods of the Huang-Huai-Hai Plain from 2000 to 2019 were taken as input variables, and the grade of yield reduction rate was taken as output variables. All samples were randomly divided into training sets and test sets, in which the training set accounted for 80% and the test set accounted for 20%. Five machine learning models, including random forest (RF), linear regression (LR), stepwise regression (STEP), Ridge regression (RIDGE), and support vector machine (SVM), were selected. The prediction accuracy of all models was mainly judged by the following indicators: Root mean square error (RMSE), mean square error (MSE), mean absolute percentage error (MAPE), and coefficient of determination (R²). It can be seen from Figure 4 that the daily NSPEI-Kc index was better than the daily NSPEI index, which proved that the daily NSPEI-Kc index was closer to the actual drought process, and the random forest model performed the best among the evaluation indexes. Therefore, the random forest had a better learning ability than other models in this study.

3.2. Temporal and Spatial Characteristics of Water Demand (Deficiency)

From the water demand of the whole growing period (Figure 5), the average annual water demand was 233.6~465.5 mm, which decreased with the increase in latitude. The annual average water demand was 101.8~189.5 mm, decreasing gradually from southeast to northwest. The tillering, jointing, heading, and mature stages showed a decreasing trend from south to north.

As can be seen from Figure 6, the spatial distribution of effective precipitation in each growth period decreased from south to north. The effective precipitation was larger in the central and southern regions, but less in the northern part of the Huang-Huai-Hai Plain. This indicated that the soil of the crop root layer absorbed more water in the central and southern regions, which was one of the reasons for the greater water deficiency in the central region.

It can be seen in Figure 7 that the water deficiency range in the sowing period was 73.9~133.5 mm, and the water deficiency decreased from southeast to northwest. The tillering stage 15.2~34.2 mm, jointing stage 19.4~38.7 mm, heading stage 22.2~39.0 mm, and maturity stage 21.0~37.8 mm were decreasing from middle to south. The whole growth period was 157.6~264.1 mm, decreasing from the central part to the northwest and northeast, showing a relative water deficiency in the central part.

As can be seen from Figure 8, the water deficiency in the sowing period manifested a trend of fluctuation and decline. It can be seen from

Table 3. Average annual water demand and water deficit by growth stage (ETc is the water requirement of winter wheat, and W is the water shortage of winter wheat).

	ETc (mm)				W(mm)
	Daily Avg (mm/d)	Annual Avg	Max	Min	Daily Avg (mm/d)	Annual Avg	Max	Min
Sowing	1.96	146.66	165.00	130.58	1.43	106.88	136.38	58.67
Tillering	0.74	76.67	89.14	57.50	0.51	52.16	79.44	14.82
Jointing	1.50	43.63	56.03	30.91	1.08	31.41	51.62	13.90
Heading	2.06	47.27	54.35	38.45	1.44	33.02	47.51	13.83
Maturing	2.15	47.37	56.12	38.40	1.46	32.10	54.00	8.40
The whole	0.99	361.62	401.59	330.92	0.70	255.58	338.94	187.99

that the average annual water requirement during the sowing period was 146.66 mm, the maximum was 165.00 mm, and the minimum was 130.58 mm. The average water deficiency decreased to 76.67 mm during the tillering stage (Figure 8b), fluctuated sharply from 1998 to 2002, and reached the lowest level in 2001 at 14.82 mm. In the jointing stage, the water demand showed a fluctuating upward trend, and the highest water demand in 2014 was 51.62 mm. The water deficiency in the heading stage fluctuated and increased, and decreased significantly in 1975, 1990, and 2014. The water deficiency at maturity showed a fluctuating trend and decreased significantly in 1997 and 2019. The water demand manifested a decreasing trend in the whole growth period (Figure 8f), and the average change rate of the water demand was −7.9 mm/10a.

The variation trend of water deficiency in domestic demand was roughly similar in each growth period, but the variation range of water deficiency was large. The variation trend of the sowing stage was similar to the tillering stage, except that the jointing stage and heading stage showed an uptrend, and the other growth stages manifested a downtrend. The jointing stage and heading stage are the key growth periods of winter wheat, and irrigation in this period can improve wheat yield and water use efficiency.

3.3. Spatiotemporal Patterns of Drought Characteristic Variables

It can be seen from Figure 9 that the duration of drought during the planting period ranged from 32.2 to 44.5 days. The duration of drought ranged from 45.6 to 61.7 days to tillering stage. The duration of the drought was 8.1–17.5 days in the period of pluck-ripening (the key growth period). The duration of drought during the whole growth period ranged from 171.7 to 211.5 days.

It can be seen from Figure 10 that drought intensity was relatively strong in the sowing period, and from a spatial perspective, the drought intensity in the west was higher than that in the east. In the mature period, the drought intensity in the east and the west was stronger, while the drought intensity in the middle was weaker. The drought intensity in the sowing period ranged from 28.3 to 42.9. The drought intensity ranged from 39.4 to 65.1 during the tillering stage. The drought intensity ranged from 7.6 to 20.7 in the period of plucking to ripening. The drought intensity of the whole growth period ranged from 28.3 to 42.9, decreasing from northwest to southeast.

From the above spatial distribution of drought characteristic variables, the spatial distribution characteristics of drought duration and drought intensity had a certain similar rule, the sowing stage decreased from west to east, and the tillering stage, jointing stage, heading stage and whole growth stage decreased from north to south, and the maturing stage decreased from east to west to middle. Generally speaking, the duration and intensity of drought in the north were higher than those in the south.

Figure 11 shows the temporal trend of drought duration from 1970 to 2019. Drought duration at the sowing stage showed an increasing (decreasing) trend in the south (north) (Figure 11a), while the spatial distribution of drought duration at the heading–tillering–maturity stage was opposite to that at the sowing stage and manifested a decreasing (increasing) trend in the south (north), and the increase and decrease in the number of sites accounted for 29% of all sites, respectively. The drought duration at the tillering stage decreased in the central region and increased in the northern and southern regions. The proportion of decreasing sites increased to 67% throughout the growth stage, and sites with increasing trends were mainly located in the southern part of the Anhui Province and the central part of the Huang-Huai-Hai Plain.

The temporal trends of drought intensity in different fertility periods of winter wheat during 1970–2019 can be seen from Figure 12. The rising sites accounted for 34% of all sites and the declining sites accounted for 30% during the sowing period south of the Qinling and Huaihe Rivers, with the significantly declining sites concentrated in the northern part of the Huang-Huai-Hai Plain. Stations with rising tillering accounted for 65% of the total, with declining stations concentrated in the central region. Rising sites accounted for 58% of the total number of sites in the jointing–heading–maturity period, with decreasing trends concentrated in the southeast and northwest. The sites with a decreasing trend in the whole reproductive period accounted for 77% of all sites, of which the obvious decreasing sites were concentrated in the southwestern part of the Huang-Huai-Hai Plain.

The temporal variation trend of the above drought characteristic variables was similar. The drought characteristic variables in the sowing period manifested a decreasing trend in the north of the Qinling Mountains and Huaihe River, and an increasing trend in the south. Drought characteristics showed an upward trend on the north and south sides, but a downward trend in the middle part. Only a few stations in the southeast manifested a decreasing trend in drought characteristic variables during the maturing period. The duration and intensity of drought during the whole growth period indicated a decreasing trend.

3.4. Drought Risk Assessment of Winter Wheat at Different Growth Stages

3.4.1. Ranking and Weighting of Drought Vulnerability at Different Growth Stages

According to the relative importance score of drought intensity in each growth period in yield prediction in the model, it was standardized to a value of 0–1, and then the change in the contribution of drought intensity in different growth periods to yield was evaluated.

Rain-fed agriculture in the Huang-Huai-Hai Plain is mainly distributed in Anhui, Jiangsu, and southern Henan. As can be seen from Figure 13, the tillering stage had a greater effect on yield, mainly due to the longer duration of the tiller stage. This was followed by the jointing and heading periods. Natural precipitation had a greater impact on the growth of farmland in rainfed agricultural areas compared to non-rainfed agricultural areas. The soil was in a dry state during the period when overwintering and rewetting watering were required at the nodulation and tasseling stages, and no artificial irrigation treatment was performed. Winter wheat was affected by climatic factors, lack of precipitation replenishment, and soil moisture was reduced, but this was the critical period for the formation of spike numbers of winter wheat, so it led to a reduction in the final yield.

Shandong Province and Hebei Province belong to non-rain-fed agricultural areas, and the weight of drought vulnerability in different growth periods was similar. This region was treated with artificial irrigation at the tillering stage, jointing stage, and heading stage, which resulted in less drought vulnerability than in rain-fed farming areas. Combined with the results from different provinces and their mean values, it was found that drought had a great effect on tillering, jointing, and heading.

3.4.2. Spatial and Temporal Distribution of Drought Risk Index

Figure 14 shows that the yield reduction rate of winter wheat in the Huang-Huai-Hai Plain gradually decreased from south to north, with high-value areas appearing in the south of Anhui Province, southeast of Henan Province, west and southwest of Jiangsu Province, and low-value areas appearing in Beijing, Tianjin, northern Hebei Province, Shandong and eastern coastal areas of Jiangsu Province and northern Henan Province, with low yield.

It can be seen from Figure 15 that the spatial distribution of the annual mean value of drought risk growth period in the Huang-Huai-Hai Plain was as a whole high in the south and low in the north. The high-value area included northern Anhui and part of prefecture-level cities in eastern Henan, northwest Shandong, and northeast Hebei. The risk index of the sowing stage and the tillering stage was higher than that of other growth stages.

Since the growth time of each growth period is quite different, the daily mean value of each growth period was obtained to compare the influence of short-duration drought on the drought risk index. It can be seen from Figure 16 that the daily mean value of the drought risk index was significantly different from the total value of the growth stage. The results showed that the average value of the drought risk index at the jointing, heading, and maturing stage was significantly higher than that at the sowing and tillering stage, indicating that the disasters at the jointing, heading, and maturing stage were more likely to increase drought risk and had a greater effect on yield.

The areas with high drought risk were mainly distributed in three regions: (1) the border areas of eastern Henan Province and southwestern Shandong Province; (2) the border areas of western Shandong Province and northwestern Shandong Province and Northeast of Hebei Province. The natural conditions of these two regions were relatively poor, and the natural precipitation was less. (3) Northern Anhui is in the transitional zone of the north–south climate, and the soil is mainly sand and ginger black soil, which has poor cultivability and is prone to waterlogging and drought, which is not conducive to stable and high agricultural yield. The low-value areas are mainly distributed in: (1) The low-value areas along the east coast of Shandong and the east coast of Jiangsu: the average yield reduction rate of these two areas is low, and they are near the sea area with sufficient natural precipitation. (2) Northwestern Hebei–Northern Henan low-value region: this region starts from Handan, Hebei in the south, passes Xingtai and Taihang Mountains, and reaches the southern foot of Yanshan Mountain in the north. This region is located on the windward side of Yanshan Mountain in Taihang Mountains, and has sufficient precipitation in summer, providing sufficient water for winter wheat growth. (3) Low-risk belt in south-central Henan and south-western Henan: the average yield reduction rate in this area is low, and drought damage occasionally occurs in wheat sowing.

In terms of the time variation trend of the drought risk index, the annual mean and the daily mean were consistent. The overall sowing time manifested a downward trend, but it showed an upward trend in the southwest and northeast. There was an obvious increase in the risk index of the southeast region during the tillering stage, which was higher than that in the sowing stage. The risk index of the whole growth stage manifested a downward trend, except that in the southern part of the Huang-Huai-Hai Plain at the maturity stage, and the disaster risk index of most prefecture-level cities manifested no great change. The whole growth period indicated a downward trend, except the central and southern parts of the Huang-Huai-Hai Plain manifested an upward trend.

4. Discussions

The water deficit at different growth stages of crops also had different effects on the components of winter wheat yield, so a comprehensive understanding of crop water demand changes is crucial for effective agricultural water management [45,46]. Some studies have shown that the influence of meteorological elements on irrigation water demand is more complex, and there is a negative correlation with the intensity of the growth period [47]. The soil’s relative water content of about 80% during the sowing period is the best for water absorption, germination, and seedling growth. If the water content is lower than 70%, the emergence of wheat seedlings will be delayed, and there will be a lack of seedlings, which is not conducive to the growth of wheat seedlings. The tillering stage was the longest-growing stage of wheat. It is also a crucial period to determine the number of ears per mu and lay the big ears. During this period, wheat gradually stopped growing and ETc was the lowest (0.74 mm/d). After overwintering, with the increase in wheat evaporation demand and the continuous development of canopy, crop evapotranspiration increased [48], and ETc reached 1.5 mm/d at the jointing stage. From the jointing stage to the maturity stage, wheat grows rapidly, which is the key period to determine wheat grain weight and yield, and also the peak of wheat water demand. ETc gradually rose to 2.15 mm/d. The sowing quantity and sowing time also had a great relationship with wheat yield and wheat quality [49]. The variation trend of ETc in different growth periods was consistent with the water requirement of wheat crops.

The Huang-Huai-Hai Plain was chosen as the study area because it is a paramount base for high-quality wheat production in China. However, due to its wide area and large span between north and south, the planting period of winter wheat in the plain may be quite different. In this paper, a set of crop planting periods and crop coefficient Kc were used to calculate the daily NSPEI-Kc index of winter wheat in different growth periods, which may have an impact on the identification of daily NSPEI-Kc drought process in different regions of the Huang-Huai-Hai Plain. However, one of the focuses of this study was to explore the recognition of drought processes in different growth stages of winter wheat with the addition of daily drought index and crop coefficient. Through comparison of machine learning model results of daily NSPEI-Kc and daily NSPEI drought index combined with yield reduction, the results showed that daily NSPEI-Kc was superior to daily NSPEI in drought recognition. In the current calculation process of the agricultural drought index, the actual evaporation based on crop coefficient was used instead of potential evaporation, which can increase the accuracy of actual crop evaporation.

The drought process was affected by drought duration in each growth stage, and there were longer drought durations and more drought processes in the sowing stage and tillering stage, which led to greater drought intensity. The results indicated that the drought process during jointing–heading–maturity had a great influence on yield, which was similar to the previous research results and the actual situation [50]. The calculation results of the yield disaster risk index also showed that drought risk was more likely to increase in the disaster period of a shorter jointing–heading–maturity period. From the relationship between drought risk and yield reduction, the spatial distribution of drought risk decreased from south to north. Drought risk was affected by drought vulnerability, and drought vulnerability accounted for a large number of growth stages including the tillering stage, jointing stage, and heading stage. From the combination of water demand and drought variables, the water shortage and drought intensity in the jointing stage and heading stage both indicated an upward trend, and the increase in water shortage increased drought intensity. Therefore, it is necessary to strengthen irrigation measures in tillering, jointing, and heading stages to reduce the risk of crop water shortage.

5. Conclusions

In this study, the daily NSPEI-Kc index based on crop coefficient was established, and the daily NSPEI index and daily NSPEI-Kc index combined with yield reduction were compared and analyzed, which verified the superiority of daily NSPEI-Kc in drought monitoring and could better describe the drought process of crops. At the same time, the daily NSPEI-Kc index and water demand were combined to monitor the water deficit and drought process of winter wheat in Huang-Huai-Hai Plain. On this basis, the contribution of drought intensity in different growth periods to winter wheat yield in different provinces was evaluated, the influence weights of drought process in different growth periods on the final yield were analyzed, the drought risk index was constructed, the temporal and spatial distribution of drought risk in different growth periods of winter wheat in the Huang-Huai-Hai Plain was analyzed, and the drought process of winter wheat was refined to each growth period. It is of great significance to ensure the sustainable and stable development of the agricultural economy in Huang-Huai-Hai Plain. From the analysis, the following findings and conclusions were drawn:

(1) In terms of the spatial distribution of water demand, water demand, and effective precipitation generally decreased with the increase in latitude, while water demand decreased from north to south under the influence of effective precipitation. In terms of time variation trend, only the jointing stage and heading stage of the key water demand growth period manifested an upward trend, while the other growth periods indicated a downward trend. The fluctuation trend of domestic water deficiency in each growth period was roughly similar, but the fluctuation of water deficiency was relatively obvious. Due to the influence of effective precipitation, the central and southern parts of the Huang-Huai-Hai Plain were short of water, so artificial agricultural measures should be taken to supplement meteorological conditions, such as temperature and light in different growth periods, to improve the transpiration rate of crops, increase the effective precipitation, and improve the water supply capacity of crops.

(2) In terms of the temporal and spatial characteristics of drought characteristic variables, the spatial distribution characteristics of drought duration and intensity had certain similar rules. The sowing period decreased from west to east, the maturity period decreased from east to west to central, and the other growth periods decreased from north to south. In general, the duration and intensity of winter wheat drought in the north of Huang-Huai-Hai Plain were higher than those in the south. The time variation trend of the drought variable also had a similar law. The sowing time indicated a decreasing trend in the north and an increasing trend in the south of the Qinling-Huaihe River. The tillering stage manifested up on both sides of the north and south, but down in the middle. Only a few sites in the southeast indicated a downward trend during the jointing–heading–maturity period. The duration and intensity of drought in more than 65% of the whole growth period manifested a decreasing trend.

(3) On the Huang-Huai-Hai Plain, the main growth stages that affected drought risk were the tillering stage, jointing stage, and heading stage, the percentage of drought vulnerability was 0.25, 0.21, 0.19, and the risk of drought was greater during these periods. From the perspective of the time variation trend, the whole growth period indicated a downward trend except for the central and southern parts of the Huang-Huai-Hai Plain. From the perspective of spatial distribution, the high drought risk areas were mainly distributed in the northeast and central, and southern regions, mainly due to poor natural conditions, less precipitation, and insufficient water resources utilization capacity. The low-value areas were mainly distributed in the eastern coastal areas, the northern part of Northwest China, and the southern part of Henan, mainly because of low yield reduction, more natural precipitation, and low disaster risk.

The disaster risk of yield in different growth periods can guide the drought resistance in the later period to reduce the loss caused by drought. The combination of the daily NSPEI index and crop coefficient can provide a theoretical basis and technical support for drought monitoring and drought assessment and prediction of winter wheat in the Huang-Huai-Hai Plain.