Dynamic Harvest Index Estimation of Winter Wheat Based on UAV Hyperspectral Remote Sensing Considering Crop Aboveground Biomass Change and the Grain Filling Process

: The crop harvest index (HI) is of great signiﬁcance for research on the application of crop variety breeding, crop growth simulation, crop management in precision agriculture and crop yield estimation, among other topics. To obtain spatial information on the crop dynamic HI (D-HI), taking winter wheat as the research object and fully considering the changes in crop biomass and the grain ﬁlling process from the ﬂowering period to the maturity period, the dynamic f G (D- f G ) parameter was estimated as the ratio between the aboveground biomass accumulated in different growth periods, from the ﬂowering stage to the maturity stage, and the aboveground biomass in the corresponding periods. Based on the D- f G parameter estimation using unmanned aerial vehicle (UAV) hyperspectral remote sensing data, a technical method for obtaining spatial information on the winter wheat D-HI was proposed and the accuracy of the proposed method was veriﬁed. A correlation analysis was performed between the normalized difference spectral index (NDSI), which was calculated using pairs of any two bands of the UAV hyperspectral spectrum, and the measured D- f G . Based on this correlation analysis, the center of gravity of the local maximum region of R 2 was used to determine the sensitive band center to accurately estimate D- f G . On this basis, remote sensing estimation of the D- f G was realized by using the NDSI constructed by the sensitive hyperspectral band centers. Finally, based on the D- f G remote sensing parameters and the D-HI estimation model, spatial information on the D-HI of winter wheat was accurately obtained. The results revealed ﬁve pairs of sensitive hyperspectral band centers (i.e., λ (476 nm, 508 nm), λ (444 nm, 644 nm), λ (608 nm, 788 nm), λ (724 nm, 784 nm) and λ (816 nm, 908 nm)) for D- f G estimation, and the results of the D- f G remote sensing estimation showed high precision. The root mean square error (RMSE) was between 0.0436 and 0.0604, the normalized RMSE (NRMSE) was between 10.31% and 14.27% and the mean relative error (MRE) was between 8.28% and 12.55%. In addition, the D- f G parameter estimation, using the NDSI constructed by the above ﬁve sensitive remote sensing band centers, yielded highly accurate spatial D-HI information with an RMSE between 0.0429 and 0.0546, an NRMSE between 9.87% and 12.57% and an MRE between 8.33% and 10.90%. The D-HI estimation results based on the hyperspectral sensitive band centers λ (724 nm, 784 nm) had the highest accuracy, with RMSE, NRMSE and MRE values of 0.0429, 9.87% and 8.33%, respectively. The proposed method of acquiring spatial information on the winter wheat D-HI in this study was shown to be feasible, and it might provide a technical reference toward developing satellite-based indices to monitor large-scale crop HI information. Overview of the methods applied in


Introduction
The crop harvest index (HI) is an important biological parameter to evaluate the level of crop yield and cultivation effectiveness, and it is an important determinant of further improvement of crop yield [1]. For grain crops, the HI refers to the ratio of grain yield to the total aboveground biomass, which essentially reflects the distribution ratio of crop assimilation products in grains and vegetative organs [2,3]. Obtaining crop HI information quickly and accurately is of great scientific significance for agronomists and breeding experts to select new crop varieties, improve varieties and optimize crop cultivation technology. At the same time, it is also of great significance for guiding agricultural management departments in their efforts to assess crop growth and estimate crop yield in a timely manner and effectively carry out agricultural production management [4][5][6].
At present, research on the crop HI is mainly based on field-scale agronomic experiments, mostly focusing on crop variety breeding and dynamic simulation of the crop HI and its influencing factors [7][8][9]. There are relatively few studies on the methods for extracting spatial information on the crop HI at the regional scale. The methods of obtaining the crop HI mainly include the direct method and the indirect method. The direct method involves calculation based on field measurements according to the ratio of grain yield to the total aboveground biomass. Although this method is accurate, it is time-consuming and laborious, it is difficult to carry out over a large range and it is impossible to obtain continuous spatial distribution information for the HI. Indirect methods mainly include two types. The first category regards the crop HI as a function of time and achieves accurate estimation of the crop HI by simulating the functional relationship between the crop dynamic HI and time. This method helps guide crop HI simulation and estimation [10,11]. Another approach is to establish the functional relationships between vegetation information (various remote sensing vegetation indexes, biomass, crop water, etc.) and environmental factors (such as temperature, light, soil water and soil nutrients) in relation to crop growth and the HI to complete the estimation of the HI [12,13].
Remote sensing has gradually become an effective technical means for obtaining largescale crop HIs because of its advantages in speed, accuracy and coverage area. Since the growth of grain crops is mainly divided into the vegetative growth stage (emergence stagebefore flowering) and the reproductive growth stage (flowering stage-maturity stage), these two stages have an important impact on the crop aboveground biomass and grain yield. Most remote sensing methods are based on these two growth stages to obtain crop harvest index information. For example, Moriondo et al. [11] obtained normalized difference vegetation index (NDVI) data in the winter wheat growing season by remote sensing and constructed the parameter 1-NDVI post /NDVI pre to estimate the spatial distribution of the HI according to the average value of the NDVI before and after flowering. This method has important reference significance for using remote sensing information to obtain the HI at the regional scale. This method was further applied by Du et al. [14] in the estimation of the winter wheat HI in Yucheng City, Shandong Province, China. According to the formation mechanism and concept of crop harvest index, Ren et al. [15] used the cumulative NDVI value from the flowering stage to the milking stage (∑NDVI post_anthesis ) to reflect the dry matter accumulation process of winter wheat grain, and they used the cumulative NDVI value from the reviving stage to before the flowering stage (∑NDVI pre_anthesis ) to reflect the dry matter accumulation process of the crop stem and leaf, which affected the later yield formation of winter wheat. Then, the ratio parameter ∑NDVI post_anthesis /∑NDVI pre_anthesis was used to characterize the crop harvest index. By establishing the statistical relationship between the ratio parameter and the measured harvest index of winter wheat, the accurate estimation of the winter wheat HI was realized in Hengshui City, Hebei Province, China. In addition, in the estimation of crop HI values using the indirect method, the ratio f G between the cumulative aboveground biomass during the flowering stage and the maturity stage, and the aboveground biomass during the maturity period, had been widely used and good estimates of the crop HI had been obtained [16,17]. However, this parameter was only applied in the calculation of f G at the maturity stage and the dynamic change in the f G parameter between the flowering stage and the maturity stage has not been considered. It is valuable to investigate the development of a large-scale estimate of the crop HI and f G using remotely sensed data because all f G parameters were calculated using actual measurement data in the aforementioned studies and this method has not been applied at the regional level.
In recent years, unmanned aerial vehicle (UAV) remote sensing technology has also developed rapidly and has become a new means of agricultural remote sensing monitoring and a useful supplement to satellite remote sensing platforms [18][19][20]. Because UAV remote sensing technology has significant advantages, such as a low cost, a small size, simple operation, flexibility, a short operation cycle and little influence from cloud cover, and the load sensor can quickly obtain images with high temporal, spatial and hyperspectral resolutions, UAV remote sensing can meet the needs of small-scale crop remote sensing monitoring and high-frequency acquisition of growth information. Currently, UAV remote sensing sensors mainly include digital cameras, multispectral cameras and hyperspectral cameras. Among them, hyperspectral cameras have more bands and can obtain band information closely related to crop growth conditions, which can provide a richer source of information for dynamic crop growth monitoring and reliably collect information on dynamic crop HI change. Overall, UAV hyperspectral remote sensing has been studied and mainly applied in fine-scale crop type identification [21], physical and chemical parameter estimation (such as the leaf area index (LAI), aboveground biomass and chlorophyll content) [22,23], crop nutrition diagnosis and growth monitoring [24], agricultural disease and pest monitoring [25,26] etc., but there are no relevant reports on its use in the crop HI estimation.
To meet the needs of fine-scale crop management by obtaining key parameters of crop growth in precision agriculture, based on the crop f G parameters proposed by predecessors, winter wheat was taken as the research object and dynamic changes in the aboveground biomass, the grain filling process and the HI from the flowering period to the maturity period were fully considered. Research on the estimation of dynamic HI values for winter wheat was carried out by UAV hyperspectral remote sensing. We hope that this method will serve as a fast and effective method for estimating the dynamic HI of winter wheat during the key growth period, as well as provide a technical reference for crop remote sensing monitoring in precision agriculture.

Study Area
The study area is located in Shenzhou County (37.71 •~3 8.16 • N, 115.36 •~1 15.80 • E), Hengshui City, Hebei Province, China, in the main grain-producing region of the North China Plain. The main planted crops in the region are winter wheat and summer corn. Winter wheat is one of the most important summer harvest crops in this area. The region has a warm, temperate, semiarid, monsoon climate with an average annual temperature of approximately 13.4 • C and an average annual precipitation of approximately 480 mm. The planting time of winter wheat in the study area is in early to mid-October, the reviving stage is in early to mid-March of the next year, the jointing stage is from late March to early to mid-April, the heading-flowering period is from late April to early May, the grain filling stage is from mid-May to late May and the maturity stage is in early June. The specific location of the study area is shown in Figure 1.
UAV remote sensing data acquisition was mainly carried out in Yuke Town and Hujiachi Town of Shenzhou County. After the field investigation of winter wheat distribution and crop growth conditions in the study area, six representative plots were selected as the experimental area for UAV remote sensing data acquisition. Each UAV flight plot was an approximate 100 m × 100 m square, with an area of approximately 1 hectare (10,000 m 2 Hujiachi Town of Shenzhou County. After the field investigation of winter wheat distribution and crop growth conditions in the study area, six representative plots were selected as the experimental area for UAV remote sensing data acquisition. Each UAV flight plot was an approximate 100 m × 100 m square, with an area of approximately 1 hectare (10,000 m 2 ). The six plots were sequentially numbered Plot I, Plot II, Plot III, Plot IV, Plot V and Plot VI. In this study, three UAV flight experiments were conducted on 3 May 2021 (flowering period), 25 May 2021 (filling period) and 4 June 2021 (maturity period). At the same time, the aboveground biomass and grain yield data of winter wheat were obtained on the ground. The UAV flight experimental plots and ground sampling points are shown in Figures 1 and 2.   Hujiachi Town of Shenzhou County. After the field investigation of winter wheat distribution and crop growth conditions in the study area, six representative plots were selected as the experimental area for UAV remote sensing data acquisition. Each UAV flight plot was an approximate 100 m × 100 m square, with an area of approximately 1 hectare (10,000 m 2 ). The six plots were sequentially numbered Plot I, Plot II, Plot III, Plot IV, Plot V and Plot VI. In this study, three UAV flight experiments were conducted on 3 May 2021 (flowering period), 25 May 2021 (filling period) and 4 June 2021 (maturity period). At the same time, the aboveground biomass and grain yield data of winter wheat were obtained on the ground. The UAV flight experimental plots and ground sampling points are shown in Figures 1 and 2.

Data Acquisition and Processing
In this study, the data acquired mainly included the crop aboveground biomass, crop grain yield and UAV hyperspectral data. The HI and f G parameters were calculated by the crop aboveground biomass and grain yield data. UAV remote sensing data acquisition and ground data collection were mainly carried out with sunny, windless and cloudless weather. Finally, data acquisition was conducted for three dates: 3 May 2021 (flowering period), 25 May 2021 (filling period), and 4 June 2021 (maturity period). The UAV flew over a total of six plots, each plot included nine evenly distributed ground sampling points, and the position information of each sampling point was recorded by a differential global positioning system (DGPS) for the collection of the aboveground biomass during the flowering-maturity period of winter wheat. The specific point distribution is shown Remote Sens. 2022, 14, 1955 5 of 21 in Figure 2. The collected ground data were mainly divided into two parts. One part was used to construct the crop HI and f G parameter estimation model and the remaining part was used as verification data to evaluate the accuracy of the f G parameter and crop HI estimation model. Finally, a total of 108 ground data sampling points were obtained from six plots, which were divided into a modeling data set (N = 72) and a verification data set (N = 36) at a ratio of 2:1.

Acquisition of Aboveground Biomass Data
During ground data collection in the key growth period of winter wheat, the position information of each sampling point was recorded by the DGPS. The aboveground part of winter wheat with a row length of 20 cm was taken as the sample at each sampling point. Next, after the separation of the stems, leaves and ears, the winter wheat was put into an oven and dried at 105 • C for 30 min. Then, the samples were dried at 85 • C to a constant weight and the dry matter of the stems, leaves and ears of winter wheat at the sampling point was weighed and recorded. Finally, the dry biomass weight of winter wheat at each sampling point was obtained [27][28][29].

Acquisition of the Dynamic Harvest Index (D-HI)
For grain crops (such as wheat and corn), during the period from flowering to maturity, the HI increases with continuous grain filling. As the percentage of crop grain yield to aboveground dry matter increases gradually, the final HI reaches the maximum. Therefore, in this study, the HI obtained at each ground observation time during the grain filling process of winter wheat was called the dynamic HI (D-HI). To obtain the dry weight of the stems, leaves and ears of winter wheat at each sampling point, the wheat ears at each sampling point were threshed and the grain weight at each sampling point was recorded. Finally, the winter wheat D-HI at each ground observation time from grain filling to maturity was calculated. The calculation formula is as follows: where t is the ground sampling time from grain filling to maturity, W G,t is the grain dry weight of winter wheat during the period from filling to maturity and W A,t is the total aboveground dry weight of winter wheat at sampling time t.

UAV Hyperspectral Data Acquisition and Processing
(1) UAV data acquisition and preprocessing In this study, a DJI M600 Pro UAV (SZ DJI Technology Co., Shenzhen, China) remote sensing platform equipped with a Resonon Pika L imaging hyperspectrometer (Resonon Inc., Bozeman, MT, USA) was used for hyperspectral image acquisition. The main technical parameters of the Pika L imaging system are listed in Table 1. For the convenience of data use, the data interval after resampling was 4 nm and 150 spectral bands were obtained. The UAV hyperspectral data acquisition time was local 10:00-14:00, the flight altitude of the UAV was 100 m and the ground resolution of the obtained hyperspectral image was 0.05 m.
The preprocessing of UAV hyperspectral data mainly included image correction, image mosaicking and reflectivity extraction of the sampling points. Image correction included radiometric calibration, atmospheric correction and orthorectification, which were all performed in the UAV image processing software Spectronon Pro (Resonon Inc., Bozeman, MT, USA). Among them, the radiometric calibration mainly depended on the reflectivity curve data of the target, pre-arranged in each UAV flight plot as the correction standard file for radiation correction processing. Then, based on the geographical registration of images, image mosaicking was carried out. Image mosaicking included heading mosaicking and side mosaicking, both of which were completed in ENVI. Finally, from the hyperspectral data, smoothed using the Savitzky-Golay (S-G) filter, the hyperspectral Remote Sens. 2022, 14, 1955 6 of 21 reflectance data corresponding to the sampling points in the UAV hyperspectral image were extracted, according to the position of the ground sampling points. (2) Calculation of the normalized difference spectral index (NDSI) The spectral index can indicate the pigment content, moisture change and nutritional status of green vegetation through the combination of specific bands [23]. In this study, to better use the information in each wavelength of hyperspectral data, pairs of bands of hyperspectral reflectance data from UAVs were combined to construct the NDSI [30,31]. The construction form is as follows: where R λ 1 and R λ 2 are the spectral reflectance values that correspond to the wavelengths of λ 1 and λ 2 , respectively. λ 1 and λ 2 are any two wavelengths in the range of 400~1000 nm. The NDSI range is [−1, 1]. Considering that the crop spectrum at 1350~1415 nm and 1800~1950 nm is greatly affected by the atmosphere and water vapor [32], sensitive band screening for D-f G estimation and remote sensing estimation of the D-HI were carried out in the wavelength range of 400~1000 nm in this study.

The Proposed Dynamic f G (D-f G ) Parameter
From the published research literature [16,17,33,34], the f G parameter (i.e., the ratio between the cumulative aboveground biomass in the flowering-maturity period and the aboveground biomass in the maturity period) was generally calculated in the maturity period, but the dynamic change in f G from the flowering period to the maturity period has not been considered, which may reduce the accuracy of the crop HI estimation using f G . To improve the accuracy of the crop HI estimation model, considering the dynamic change in the f G parameter between the flowering period and the maturity period, the dynamic parameter D-f G (dynamic f G ) was proposed in this study, which is the ratio of the aboveground biomass accumulated from the flowering period to the sampling time to the whole aboveground biomass at the corresponding sampling time. The calculation formula of the D-f G parameter is as follows: where ΣW post is the accumulated aboveground biomass of winter wheat from the flowering period to the sampling time (kg/ha); ΣW whole is the whole aboveground biomass corresponding to the sampling time (kg/ha); t is the sampling time from the flowering period to the maturity period; W A,t is the aboveground dry matter weight of winter wheat at sampling time t (kg/ha); and W F is the aboveground dry matter weight of winter wheat at the flowering stage (kg/ha).

Overall Technical Process
The main process in this study included three parts. The first part was the acquisition and preprocessing of the UAV remote sensing data and ground observation data, as well as the calculation of relevant parameters, such as the NDSI, measured D-f G and measured D-HI, at sampling points. The second part was mainly to determine the sensitive band centers for the D-f G remote sensing estimation. Specifically, the correlation between the NDSI and measured D-f G at sampling points was determined, and a two-dimensional map of fit accuracy (R 2 ) between the NDSI and the D-f G of winter wheat was obtained; then, the sensitive band centers for winter wheat D-f G estimation were identified by determining the local maximum region of R 2 and its center of gravity. The third part was the estimation and verification of the D-HI based on the D-f G parameters obtained from UAV remote sensing data. First, after determining the sensitive band centers for winter wheat D-f G estimation, the corresponding NDSI was constructed by using the reflectivity of the sensitive band of the UAV hyperspectral image, the winter wheat D-f G estimation model was completed and spatial information extraction and accuracy verification of the D-f G parameters were performed. Then, the D-HI estimation model was established based on the measured D-f G and measured D-HI values. Finally, based on the D-f G remote sensing information obtained, spatial information on the D-HI of winter wheat was extracted based on the UAV hyperspectral image and the accuracy was verified. The details of the entire procedure are shown in Figure 3.
where ΣWpost is the accumulated aboveground biomass of winter wheat from the flowering period to the sampling time (kg/ha); ΣWwhole is the whole aboveground biomass corresponding to the sampling time (kg/ha); t is the sampling time from the flowering period to the maturity period; WA,t is the aboveground dry matter weight of winter wheat at sampling time t (kg/ha); and WF is the aboveground dry matter weight of winter wheat at the flowering stage (kg/ha).

Overall Technical Process
The main process in this study included three parts. The first part was the acquisition and preprocessing of the UAV remote sensing data and ground observation data, as well as the calculation of relevant parameters, such as the NDSI, measured D-fG and measured D-HI, at sampling points. The second part was mainly to determine the sensitive band centers for the D-fG remote sensing estimation. Specifically, the correlation between the NDSI and measured D-fG at sampling points was determined, and a two-dimensional map of fit accuracy (R 2 ) between the NDSI and the D-fG of winter wheat was obtained; then, the sensitive band centers for winter wheat D-fG estimation were identified by determining the local maximum region of R 2 and its center of gravity. The third part was the estimation and verification of the D-HI based on the D-fG parameters obtained from UAV remote sensing data. First, after determining the sensitive band centers for winter wheat D-fG estimation, the corresponding NDSI was constructed by using the reflectivity of the sensitive band of the UAV hyperspectral image, the winter wheat D-fG estimation model was completed and spatial information extraction and accuracy verification of the D-fG parameters were performed. Then, the D-HI estimation model was established based on the measured D-fG and measured D-HI values. Finally, based on the D-fG remote sensing information obtained, spatial information on the D-HI of winter wheat was extracted based on the UAV hyperspectral image and the accuracy was verified. The details of the entire procedure are shown in Figure 3.

Determination of the Sensitive Band Centers for D-f G Estimation
Since a correlation between the hyperspectral data bands will increase the redundancy of spectral information, to improve the accuracy of the D-f G remote sensing estimation model, the NDSI was used to screen the sensitive band centers for the D-f G parameter estimation in this study. During the research, based on MATLAB software, a correlation analysis between winter wheat D-f G and NDSI, based on UAV hyperspectral data in the wavelength range of 400~1000 nm, was performed and a two-dimensional diagram of the fit (R 2 ) between the NDSI and winter wheat D-f G data was created. Because the distribution of the local maximum region of R 2 values might be uneven, the local maximum of R 2 and its center of gravity might not be identical; that is, the corresponding band of the local maximum of R 2 might not coincide with the center of the optimal band. Therefore, to Remote Sens. 2022, 14, 1955 8 of 21 ensure the stability of the crop D-f G estimation results, the center of gravity of the local maximum region of R 2 proposed by Liu et al. [35] was used to obtain the sensitive band centers for winter wheat D-f G estimation. The specific process is as follows: First, based on a two-dimensional R 2 map between the NDSI and winter wheat D-f G , the band wavelength range with a high correlation between these two parameters was determined. Second, the local maximum of R 2 was found and all significant points in the eight neighborhoods of the local maximum of R 2 value were traversed. This set of points was marked Ω as the local maximum region of R 2 . Finally, the center of gravity of the local maximum region of R 2 was calculated and used as the sensitive band center of each local maximum region. The specific calculation formula for the center of gravity of the local maximum region of R 2 is as follows: where f (u,v) is the R 2 value of the band coordinate of (u,v); Ω is the local maximum region of R 2 ; and (u, v) is the center coordinate of the sensitive band.

Establishment of the D-f G Estimation Model Based on Hyperspectral NDSI Data
Based on the sensitive band centers for winter wheat D-f G estimation, the NDSI was determined and used for the D-f G remote sensing estimation in this study to provide accurate D-f G remote sensing estimation parameters to acquire the winter wheat D-HI. The linear model between the NDSI and D-f G was established, as shown in Formula (5).
where λ 1 and λ 2 are the UAV hyperspectral bands corresponding to the selected sensitive band centers; a is the coefficient of the first-order term; and b is a constant. D-f G is the ratio between the aboveground biomass accumulated from the flowering period to sampling time and the aboveground biomass accumulated from the growth period to sampling time.

Establishment of the D-HI Estimation Model Based on D-f G Remote Sensing Parameter
Based on the crop HI estimation method including f G in the maturity period, which was proposed by Kemanian et al. [16], a UAV remote sensing-based method for estimating the D-HI with the D-f G information was proposed in this study. The statistical model between the D-f G and the D-HI is shown in Formula (6).
where HI 0 is the intercept, or the value of D-HI when the biomass does not change after the flowering period of crops, that is, when D-f G is 0; and k is the slope of the linear relationship between the D-HI and D-f G .

Accuracy Evaluation of the D-f G Estimation Model and D-HI Estimation Model
In this study, the coefficient of determination (R 2 ), root mean square error (RMSE), normalized RMSE (NRMSE) and mean relative error (MRE) were selected to verify the accuracy of the D-f G and D-HI estimation models. The closer R 2 is to 1, the better the quality of the model. The closer the RMSE is to 0, the closer the simulations are to the observations. When the NRMSE is ≤10% and MRE is ≤10%, the simulation result accuracy is considered excellent; when 10% < NRMSE ≤ 20% and 10% < MRE ≤ 20%, the simulation result accuracy is considered good; when 20% < NRMSE ≤ 30% and 20% < MRE ≤ 30%, the simulation result is considered average; and when NRMSE > 30% and MRE > 30%, the Remote Sens. 2022, 14,1955 9 of 21 accuracy of the simulation result is considered poor. The NRMSE value is given priority in the assessment of model accuracy [36][37][38]. The specific formulas are as follows: where x i is the measured value of D-f G or the D-HI of winter wheat; y i is the estimated value of D-f G or the D-HI of winter wheat; and x is the mean values of x i ; n is the number of samples. In the study, with the UAV hyperspectral reflectance data corresponding to each field sampling point, according to Formula (2), the NDSI of the field sampling points in each plot was calculated and graphed using MATLAB software and the NDSI distribution of winter wheat was obtained. There were 150 × 150 combinations of spectral bands and related NDSI values between any two bands in the 400~1000 nm hyperspectral range. Due to space limitations, only the UAV hyperspectral NDSI calculation results for different growth periods (including the flowering stage, filling stage and maturity stage of winter wheat) at a sampling point (37 • 53 55"N, 115 • 37 43"E) in Plot I are shown, and the specific results are shown in Figure 4. The abscissa (λ 1 ) and ordinate (λ 2 ) are both hyperspectral wavelengths in the range of 400~1000 nm. The points corresponding to the two-dimensional space formed by the horizontal and vertical axes are the NDSI(λ 1 , λ 2 ) values calculated by the corresponding reflectance of any two bands λ 1 and λ 2 .

Determination of Sensitive UAV Hyperspectral Band Centers for D-fG Estimation
MATLAB software was used to analyze the correlation between the D-fG and the NDSI constructed by any two bands of the UAV hyperspectral spectrum. Finally, a two-dimensional diagram of the fit accuracy (R 2 ) of the NDSI and the D-fG was obtained, as shown in Figure 5a. The horizontal and vertical axes constitute a total of 150 × 150 points corresponding to the R 2 two-dimensional space, with each point representing the R 2 value between the corresponding NDSI value and measured D-fG. Figure 5a shows that R 2 is axisymmetric with a connecting line between (400, 400) and (1000, 1000), from which the region and relevant band information with an obvious correlation between the NDSI and the D-fG can be obtained. On this basis, the centers of the sensitive hyperspectral band for the D-fG estimation could be determined. In this paper, only the area on the upper side of the symmetry axis in the R 2 two-dimensional graph was studied.

Determination of Sensitive UAV Hyperspectral Band Centers for D-f G Estimation
MATLAB software was used to analyze the correlation between the D-f G and the NDSI constructed by any two bands of the UAV hyperspectral spectrum. Finally, a twodimensional diagram of the fit accuracy (R 2 ) of the NDSI and the D-f G was obtained, as shown in Figure 5a. The horizontal and vertical axes constitute a total of 150 × 150 points corresponding to the R 2 two-dimensional space, with each point representing the R 2 value between the corresponding NDSI value and measured D-f G . Figure 5a shows that R 2 is axisymmetric with a connecting line between (400, 400) and (1000, 1000), from which the region and relevant band information with an obvious correlation between the NDSI and the D-f G can be obtained. On this basis, the centers of the sensitive hyperspectral band for the D-f G estimation could be determined. In this paper, only the area on the upper side of the symmetry axis in the R 2 two-dimensional graph was studied.
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 (a) Two-dimensional graph of R 2 ( b) Contour map of R 2 Figure 5. Two-dimensional graph and contour map of the R 2 values between the NDSI and D

Remotely Sensed NDSI-Based D-fG Estimation and Its Verification
Based on the 5 pairs of selected sensitive hyperspectral band centers, namely, λ nm, 508 nm), λ(444 nm, 644 nm), λ(608 nm, 788 nm), λ(724 nm, 784 nm) and λ(816 908 nm), the sensitive UAV hyperspectral bands for the D-fG estimation were determ and the NDSI of the six plots in the grain filling period and maturity period were c lated according to Formula (2) with the different sensitive band centers. Then, accor to Formula (5), the linear models of the NDSI and the D-fG based on different sen band centers were constructed, the D-fG spatial estimation results of six plots in the filling stage and maturity stage were obtained based on different sensitive band cen and the verification data set (N = 36) was used for D-fG accuracy verification. Fin based on the centers of the sensitive bands, the statistical relationship between the and the NDSI and the estimation accuracy of the D-fG were obtained. The specific re are shown in Table 2 and Figure 6.
As shown in Table 2, when the NDSI values constructed by the centers of the pairs of selected sensitive hyperspectral bands were fitted with the D-fG, all the rela ships reached a very significant level (p < 0.01), and the determination coefficient R between 0.6047 and 0.6843. The verification results showed that the statistical mod tablished between the NDSI and the D-fG had good estimation effects, with the Dtimation showing a high level of precision. In the estimation accuracy results of the the fit R 2 was between 0.6521 and 0.8595, RMSE was between 0.0436 and 0.0604, NR was between 10.31% and 14.27%, and MRE was between 8.28% and 12.55%. More the NDSI (724 nm, 784 nm) had the highest accuracy in estimating the winter wheat with an RMSE, NRMSE and MRE of 0.0436, 10.31%, and 8.28%, respectively.
Due to space limitations, only the UAV hyperspectral NDSI result and estim D-fG values with the highest estimation accuracy, which correspond to the sensitive centers λ(724 nm, 784 nm), are shown in this paper. The NDSI values calculated for I-VI in the grain filling period and maturity period are shown in Figures 7 and 8, an D-fG remote sensing estimates in the grain filling period and maturity period are sh in Figures 9 and 10, respectively. According to the analysis of six plots, the NDSI du the grain filling period ranged from 0.29 to 0.70, and the average NDSI value in the s area was approximately 0.47. The NDSI at the maturity stage ranged from 0.08 to and the average NDSI value in the study area was approximately 0.27. The estimat the D-fG during the grain filling period were between 0.21 and 0.50, and the average v of the D-fG in the study area was approximately 0.37; the estimates of the D-fG at the turity stage were between 0.32 and 0.65, and the average value of the D-fG in the s area was approximately 0.51. The statistical information of the NDSI and the D-fG in plot is shown in Table 3. According to the significance table for correlation tests, when the number of samples N = 72 and the significance level is 0.01, when R 2 > 0.0910, the NDSI and the D-f G show a very significant correlation. Therefore, the two-dimensional region with R 2 > 0.0910 was selected for research in this paper and the maximum points with R 2 > 0.0910 in Figure 5a were found. Then, all the points with R 2 > 0.0910 in the eight neighborhoods of this point were traversed, and this set of points was marked as the local maximum region of R 2 , which displayed R 2 = 0.1 as the gradient. To display the distribution range of the sensitive band more intuitively, only the results with R 2 ≥ 0.3 are displayed here, as shown in Figure 5b.
To improve the accuracy of the D-f G estimation by the selected sensitive hyperspectral band centers, the two-dimensional region with R 2 > 0.4 was selected to determine the sensitive band centers, as shown in the red and yellow regions in Figure 5b. Finally, Ω A~ΩE were determined as the local maximum regions of R 2 that met the screening conditions. The specific results were as follows: Ω A was 460~500 nm on the horizontal axis and 480~530 nm on the vertical axis; Ω B was 410~480 nm on the horizontal axis and 560~720 nm on the vertical axis; Ω C was 550~660 nm on the horizontal axis and 730~860 nm on the vertical axis; Ω D was 690~760 nm on the horizontal axis and 710~860 nm on the vertical axis; and Ω E was 750~890 nm on the horizontal axis and 820~980 nm on the vertical axis. After determining the sensitive area for the NDSI and the D-f G , the centers λ(λ 1 , λ 2 ) of each local maximum region of R 2 were calculated according to Formula (4), where λ 1 and λ 2 are the abscissa and ordinate of the center of the local maximum region. Finally, the centers of gravity of the local maximum region of Ω A~ΩE were determined as λ(476 nm, 508 nm), λ(444 nm, 644 nm), λ(608 nm, 788 nm), λ(724 nm, 784 nm) and λ(816 nm, 908 nm), respectively.

Remotely Sensed NDSI-Based D-f G Estimation and Its Verification
Based on the 5 pairs of selected sensitive hyperspectral band centers, namely, λ(476 nm, 508 nm), λ(444 nm, 644 nm), λ(608 nm, 788 nm), λ(724 nm, 784 nm) and λ(816 nm, 908 nm), the sensitive UAV hyperspectral bands for the D-f G estimation were determined, and the NDSI of the six plots in the grain filling period and maturity period were calculated according to Formula (2) with the different sensitive band centers. Then, according to Formula (5), the linear models of the NDSI and the D-f G based on different sensitive band centers were constructed, the D-f G spatial estimation results of six plots in the grain filling stage and maturity stage were obtained based on different sensitive band centers, and the verification data set (N = 36) was used for D-f G accuracy verification. Finally, based on the centers of the sensitive bands, the statistical relationship between the D-f G and the NDSI and the estimation accuracy of the D-f G were obtained. The specific results are shown in Table 2 and Figure 6. Table 2. Relationships between D-f G and the NDSI and their accuracy verification. Centers λ(λ 1 ,  λ 2 ) for D-f G NDSI (λ 1 , λ 2 ) and D-f G (N = 72) Estimation Accuracy of D-f G (N = 36)   λ(λ1, λ2) for D-fG Relationship between the NDSI (λ1, λ2) and D-fG (N = 72 Figure 6. Verification of the D-fG estimates based on the sensitive UAV hyperspectral band centers. Figure 6. Verification of the D-f G estimates based on the sensitive UAV hyperspectral band centers.

Sensitive Band Centers
As shown in Table 2, when the NDSI values constructed by the centers of the five pairs of selected sensitive hyperspectral bands were fitted with the D-f G , all the relationships reached a very significant level (p < 0.01), and the determination coefficient R 2 was between 0.6047 and 0.6843. The verification results showed that the statistical model established between the NDSI and the D-f G had good estimation effects, with the D-f G estimation showing a high level of precision. In the estimation accuracy results of the D-f G , the fit R 2 was between 0.6521 and 0.8595, RMSE was between 0.0436 and 0.0604, NRMSE was between 10.31% and 14.27%, and MRE was between 8.28% and 12.55%. Moreover, the NDSI (724 nm, 784 nm) had the highest accuracy in estimating the winter wheat D-f G , with an RMSE, NRMSE and MRE of 0.0436, 10.31%, and 8.28%, respectively.
Due to space limitations, only the UAV hyperspectral NDSI result and estimated D-f G values with the highest estimation accuracy, which correspond to the sensitive band centers λ(724 nm, 784 nm), are shown in this paper. The NDSI values calculated for Plots I-VI in the grain filling period and maturity period are shown in Figures 7 and 8, and the D-f G remote sensing estimates in the grain filling period and maturity period are shown in Figures 9 and 10, respectively. According to the analysis of six plots, the NDSI during the grain filling period ranged from 0.29 to 0.70, and the average NDSI value in the study area was approximately 0.47. The NDSI at the maturity stage ranged from 0.08 to 0.53, and the average NDSI value in the study area was approximately 0.27. The estimates of the D-f G during the grain filling period were between 0.21 and 0.50, and the average value of the D-f G in the study area was approximately 0.37; the estimates of the D-f G at the maturity stage were between 0.32 and 0.65, and the average value of the D-f G in the study area was approximately 0.51. The statistical information of the NDSI and the D-f G in each plot is shown in Table 3.

Establishment of a D-HI Estimation Model Based on D-fG
In this study, the D-fG and measured D-HI of 72 sampling points were calculated based on the dynamic aboveground biomass and grain yield data during the grain filling period and the maturity period. The estimation model between the D-fG and the D-HI was established according to Formula (6) as follows: The results revealed a significant linear relationship between the measured winter wheat D-fG and the winter wheat D-HI at the level of p < 0.01, and the coefficient of determination (R 2 ) for the linear model was 0.8513 ( Figure 11). Therefore, this model is promising for collecting spatial information on the D-HI based on the D-fG obtained from

Establishment of a D-HI Estimation Model Based on D-f G
In this study, the D-f G and measured D-HI of 72 sampling points were calculated based on the dynamic aboveground biomass and grain yield data during the grain filling period and the maturity period. The estimation model between the D-f G and the D-HI was established according to Formula (6) as follows: The results revealed a significant linear relationship between the measured winter wheat D-f G and the winter wheat D-HI at the level of p < 0.01, and the coefficient of determination (R 2 ) for the linear model was 0.8513 ( Figure 11). Therefore, this model is promising for collecting spatial information on the D-HI based on the D-f G obtained from remote sensing data.

Acquisition and Verification of Spatial D-HI Information via Remote Sensing-Based D-f G Data
In this study, with the D-f G estimates obtained using the NDSI data from five sensitive remote sensing band centers, spatial estimates of the crop D-HI values were calculated in six plots at the filling stage and maturity stage according to Formula (10), and the estimation accuracy of the D-HI remote sensing information was verified by using the reserved verification data set. The overall verification results of the winter wheat D-HI during the filling period and maturity period are shown in Table 4 and Figure 12. . 2022, 14, x FOR PEER REVIEW Figure 11. Establishment of the D-HI estimation model based

Acquisition and Verification of Spatial D-HI Infor D-fG Data
In this study, with the D-fG estimates obtained usi tive remote sensing band centers, spatial estimates of lated in six plots at the filling stage and maturity stage estimation accuracy of the D-HI remote sensing infor reserved verification data set. The overall verification during the filling period and maturity period are show The overall evaluation accuracy results for the Dwhen five pairs of selected sensitive hyperspectral ba D-HI estimates, based on the D-fG parameters obtaine spectral bands, were accurate. The fit accuracy R 2 wa RMSE was between 0.0429 and 0.0546, the NRMSE wa the MRE was between 8.33% and 10.90%. Moreover, band centers λ(724 nm, 784 nm), the accuracy of the sp the highest, with an RMSE, NRMSE and MRE of 0.042 based on hyperspectral sensitive band centers λ(816 spatial D-HI estimation results was relatively low, wit 0.0546, 12.57% and 10.90%, respectively.  Note: N represents the number of samples for D-HI measurement in the validation data set. ** indicates a very significant correlation at the level of p < 0.01. The overall evaluation accuracy results for the D-HI estimation in Table 4 show that when five pairs of selected sensitive hyperspectral band centers were used, the spatial D-HI estimates, based on the D-f G parameters obtained from the sensitive UAV hyperspectral bands, were accurate. The fit accuracy R 2 was between 0.6570 and 0.8573, the RMSE was between 0.0429 and 0.0546, the NRMSE was between 9.87% and 12.57% and the MRE was between 8.33% and 10.90%. Moreover, based on sensitive hyperspectral band centers λ(724 nm, 784 nm), the accuracy of the spatial D-HI estimation results was the highest, with an RMSE, NRMSE and MRE of 0.0429, 9.87% and 8.33%, respectively; based on hyperspectral sensitive band centers λ(816 nm, 908 nm), the accuracy of the spatial D-HI estimation results was relatively low, with an RMSE, NRMSE and MRE of 0.0546, 12.57% and 10.90%, respectively.
Due to space limitations, this paper only shows the D-HI estimation results for Plots I-VI with the highest estimation accuracy, corresponding to the sensitive band centers λ(724 nm, 784 nm) in the grain filling period and maturity period (Figures 13 and 14). The estimates of the D-HI during the grain filling period were between 0.29 and 0.49 among the six plots and the average D-HI in the study area was approximately 0.40. The estimates of the D-HI during the maturity period were between 0.37 and 0.59 and the average value of the D-HI in the study area was 0.49. The statistical results for the D-HI during the grain filling period and maturity period in each plot are shown in Table 5.

The Characteristics and Application Potential of the Proposed Method
(1) Full consideration of dynamic crop growth information Generally, only the crop HI in the maturity period is considered in the common crop HI estimation method and the gradual process of crop HI development is less often considered. Therefore, the dynamic change in the crop HI during grain filling and maturity was considered in this study and the crop HI in the filling period and the maturity period was estimated, with scientific significance and practical value for improving the accuracy of crop growth simulation and effectively improving the level of agricultural production management [5,39,40]. In addition, in most of the existing methods for obtaining the crop HIs based on remote sensing technology, time series of vegetation index changes were

The Characteristics and Application Potential of the Proposed Method
(1) Full consideration of dynamic crop growth information Generally, only the crop HI in the maturity period is considered in the common crop HI estimation method and the gradual process of crop HI development is less often considered. Therefore, the dynamic change in the crop HI during grain filling and maturity was considered in this study and the crop HI in the filling period and the maturity period was estimated, with scientific significance and practical value for improving the accuracy of crop growth simulation and effectively improving the level of agricultural production management [5,39,40]. In addition, in most of the existing methods for obtaining the crop HIs based on remote sensing technology, time series of vegetation index changes were used only to indirectly reflect changes in the crop growth processes, whereas the changes in the crop biomass and grain filling processes were not directly fully considered. In this study, the parameter f G at the maturity stage proposed by Kemanian et al. [16] was developed into a dynamic parameter D-f G , which included crop dynamic growth information and could improve the accuracy and stability of the crop D-HI estimation model.
A comparison with the existing research results shows that R 2 of the statistical model between the HI and the f G at the maturity stage established by Kemanian et al. [16] was 0.38-0.84, while R 2 of the statistical model between the D-HI and the D-f G established in this study, by fully considering the crop dynamic growth information between the flowering stage and the maturity stage was 0.8513. In addition, the RMSE of the crop HI estimation accuracy based on remote sensing technology is generally 0.01-0.06 [11,15,17], while the highest accuracy RMSE of the crop harvest index based on UAV remote sensing information in this study was 0.0429, NRMSE was 9.87% and MRE was 8.33%. These results demonstrate the validity of the proposed crop harvest index estimation method in this study, which fully considered the information of crop dynamic growth.
(2) High-precision acquisition of D-f G based on optimal selection of sensitive bands To transform the general f G parameters into dynamic f G parameters, based on the correlation analysis of the UAV hyperspectral NDSI and the D-f G data, the center of gravity of the local maximum region of R 2 was used in this study to screen sensitive bands for the D-f G remote sensing estimation, and finally, the high-precision D-f G values based on the screening of sensitive hyperspectral bands were acquired. The general field-scale f G parameters need to be obtained by the traditional measurement method. In this research, the accurate f G values were innovatively acquired through UAV remote sensing technology, which not only provided a new idea and technical method for the acquisition of regional f G parameters but also laid a technical foundation for the acquisition of large-scale spatial D-f G information based on remote sensing satellites.
(3) Spatial D-HI information acquisition based on D-f G remote sensing estimation At present, the estimation of crop HI values using the f G parameter at maturity can be carried out only at the field scale by relying on ground-measured data. In view of this limitation, the f G parameter was transformed into the dynamic f G parameter while fully considering the dynamic change in the crop HI. A method for obtaining spatial information on the crop D-HI based on the D-f G remote sensing estimation was proposed in this research, and the proposed method was applied and verified by using UAV hyperspectral information, allowing the acquisition of accurate dynamic spatial information on the crop HI in small areas and in different periods such as the filling period and the maturity period. This proposed method is of great significance for improving the level of small-scale crop management in precision agriculture, laying an important technical foundation for future large-scale crop HI estimation based on satellite remote sensing data.
(4) The potential application of spatial D-HI information acquisition in crops The method of acquiring spatial information on the crop D-HI using UAV hyperspectral data proposed in this research allowed not only the accurate estimation of the HI in the maturity period but also the dynamic estimation of the crop D-HI in different grain filling periods. The above results and technical methods have important application prospects for improving farmland and crop management in precision agriculture. They are also of great significance for accurately simulating and estimating crop yield, estimating the crop biomass, breeding crop varieties, evaluating the crop growth environment, optimizing and implementing evaluations of crop cultivation technology, and researching agricultural responses to climate change. In addition, only the HI of winter wheat was studied in this paper. Analyzed from a theoretical viewpoint, the method proposed in this study also has important reference value for remote sensing-based estimation of the HI in other grain crops (rice, corn, etc.). Although the UAV hyperspectral remote sensing data were used to perform specific research, they can provide guidance and a reference for regional crop HI estimation based on satellite remote sensing data (such as hyperspectral data and multispectral data) in the future and can also provide a theoretical basis for satellite band setting.
Among them, Hyperion hyperspectral satellite can be well matched with the five selected pairs of sensitive band centers in this study, which is very beneficial to the acquisition of the harvest index based on Hyperspectral satellite data. For the RapidEye satellite, the sensitive band centers λ(724 nm, 784 nm) can match the position of the red edge band and near-infrared band of the satellite, and λ(444 nm, 644 nm) can correspond to the blue band and red band of the satellite; For China's GF-6 satellite, the sensitive band λ(724 nm, 784 nm) can match the position of the red edge band 1 and the near-infrared band of the satellite. These in-orbit satellite bands can effectively match the results of the sensitive band centers in this study, which is of great significance to the large-scale application of the proposed harvest index acquisition method in this paper. However, for other broadband multispectral satellites such as MODIS, Sentinel-2, Landsat-8 and Spot-7, the bands cannot be completely individually matched with the sensitive band centers in this study. Thus, it is necessary to further study the impact of the expansion of the sensitive band center widths on the accuracy of the crop harvest index estimation to better select the available satellite remote sensing bands. In addition, the NDSI constructed by the sensitive band centers λ(724 nm, 784 nm) or λ(608 nm, 788 nm) is similar to the general NDVI, which may have better practical operability in estimating the harvest index at the region scale. In this study, the 724 nm red edge band played an important role in estimating the crop harvest index, which will have a certain reference for the setting of the red edge band in satellites in the future [41].

Shortcomings and Improvements of the Proposed Method
In this study, a method for obtaining spatial information on the crop D-HI based on Df G hyperspectral remote sensing estimation was proposed. Although accurate crop dynamic HI estimation results were obtained, this research method still needs to be further improved.
To improve the accuracy and stability of the HI estimation model, multiple plots were selected for the experiment, but the effects of crop water content on the screening of sensitive bands and the D-f G estimation were not considered. In the next step, it is necessary to perform in-depth research on this issue. In the crop HI estimation method applied in this study, only the NDSI was used for remote sensing-based estimation of the D-f G , but for other spectral indices, such as the difference vegetation index (DVI), ratio vegetation index (RVI), enhanced vegetation index (EVI) and green vegetation index (GVI), a detailed comparative study of remote sensing-based D-f G estimation accuracy and HI estimation accuracy needs to be carried out. In the present study, UAV hyperspectral data from only three periods in one year were used to estimate the HI. In the future, it will be necessary to strengthen multiyear UAV hyperspectral data acquisition and corresponding ground observation experiments in different grain filling stages (such as the early grain filling stage, middle grain filling stage and late grain filling stage) to further improve the estimation accuracy of the HI. In this study, the estimation of spatial information for the crop HI was mainly carried out through the ground-measured crop HI and UAV hyperspectral data. In the future, it will be necessary to further carry out crop HI estimation based on the combination of imaging and non-imaging hyperspectral data obtained from UAV imaging hyperspectrometers and ASD spectroradiometers (Analytical Spectral Devices Inc., Boulder, CO, USA) to further optimize the accuracy of the model and improve the reliability of the HI estimates. At the same time, by considering the spatial variation of crop phenology information and crop variety, it will be necessary to further conduct research on large-scale crop HI acquisition with UAV and satellite remote sensing, to apply the crop HI estimation method proposed in this study at the regional scale.

Conclusions
By transforming the static f G parameter into the dynamic D-f G parameter, a method to obtain spatial information for the winter wheat D-HI based on the D-f G remote sensing parameters from UAV hyperspectral data was proposed and verified in this paper. Finally, spatial information for the winter wheat D-HI was accurately estimated. Among them, five pairs of sensitive remote sensing band centers were selected to estimate the D-f G parameter: λ(476 nm, 508 nm), λ(444 nm, 644 nm), λ(608 nm, 788 nm), λ(724 nm, 784 nm) and λ(816 nm, 908 nm). The remote sensing-based D-f G estimates were verified, with an RMSE between 0.0436 and 0.0604, an NRMSE between 10.31% and 14.27% and an MRE between 8.28% and 12.55%. At the same time, for the five pairs of sensitive hyperspectral band centers, the spatial D-HI information estimates were highly accurate, with an RMSE between 0.0429 and 0.0546, an NRMSE between 9.87% and 12.57% and an MRE between 8.33% and 10.90%. The D-HI estimation results based on hyperspectral sensitive band centers λ(724 nm, 784 nm) had the highest accuracy, with RMSE, NRMSE and MRE values of 0.0429, 9.87% and 8.33%, respectively. The estimation results for the D-f G and the D-HI in this study were highly accurate, proving that the proposed method for estimating the spatial information of the winter wheat D-HI based on UAV hyperspectral data was feasible. It has certain reference significance for the estimation of large-scale crop D-HI values based on satellite remote sensing in the future.