Should phenological information be applied to predict agronomic traits across growth stages of winter wheat?

Yu Zho,Yng Meng,Shoyu Hn,Hikun Feng,Guijun Yng,*,Zhenhi Li,c,*

a Key Laboratory of Quantitative Remote Sensing in Agriculture of Ministry of Agriculture and Rural Affairs,Information Technology Research Center,Beijing Academy of Agriculture and Forestry Sciences,Beijing 100097,China

b Collaborative Innovation Center for Modern Crop Production Co-sponsored by Province and Ministry,Nanjing Agricultural University,Nanjing 210095,Jiangsu,China

c College of Geomatics,Shandong University of Science and Technology,Qingdao 266590,Shandong,China

Keywords:Agronomic traits Phenological effect Vegetation index Hierarchical linear model Winter wheat

ABSTRACT Most existing agronomic trait models of winter wheat vary across growing seasons,and it is an open question whether a unified statistical model can be developed to predict agronomic traits using a vegetation index(VI)across multiple growing seasons.In this study,we constructed a hierarchical linear model (HLM) to automatically adapt the relationship between VIs and agronomic traits across growing seasons and tested the model’s performance by sensitivity analysis.Results demonstrated that(1)optical VIs give poor performance in predicting AGB and PNC across all growth stages,whereas VIs perform well for LAI,LGB,LNC,and SPAD.(2) The sensitivity indices of the phenological information in the AGB and PNC prediction models were 0.81-0.86 and 0.66-0.73,whereas LAI,LGB,LNC,and SPAD prediction models produced sensitivity indexes of 0.01-0.02,0.01-0.02,0.01-0.02,and 0.02-0.08,respectively.(3) The AGB and PNC prediction models considering ZS were more accurate than the prediction models based on VI.Whether or not phenological information is used,there was no difference in model accuracy for LGB,LNC,SPAD,and LAI.This study may provide a guideline for deciding whether phenological correction is required for estimation of agronomic traits across multiple growing seasons.

1.Introduction

Crop-growth monitoring provides a basis for real-time knowledge of seedling conditions,farmland management,and early-yield estimation and has become an important part of precision agriculture [1].Advances in remote-sensing technology have led to the development of various statistical and physical methods to monitor agronomic traits[2-4].However,inconsistent variation in agronomic traits and vegetation indices(VIs)over entire growth periods suggests that a VI may not follow a unique mapping to agronomic traits throughout the growing season.

The key to constructing a remote sensing monitoring model is to characterize the time-series variation of agronomic traits and VIs in a growing season.Many agronomic traits are subject to allometric accumulation or reduction across growth stages due to assimilate accumulation and transportation [5],whereas the VI varies with different growth period[6,7].One set of VIs may correspond to multiple agronomic traits,such as aboveground dry biomass (AGB) or plant nitrogen concentration (PNC) [8,9].Most of the literature [10,11] agrees that,although VIs may correlate well with agronomic traits during a given growth stage,they do not transfer to other growth stages.Substantially better performance was achieved by a piecewise agronomic trait-VI relationship method based on phenological information acquired throughout the entire growing season [6,12].

Unfortunately,the discontinuity of the piecewise method can lead to uncertainty in the prediction model for uninvestigated periods[13,14].In dealing with multiple growth stages,it is difficult to estimate agronomic traits from VIs for three reasons: (i) the crop canopy structure and background information in the field observation change constantly,causing the sensitivity of VIs to agronomic traits to depend on the growth stage for which the VI is constructed,(ii) the reflectance of optical bands and VIs saturate at medium-to-high canopy cover,and (iii) optical VIs are not sensitive to obscured information,such as information from stems and lower layers of vegetation close to the ground.

Numerous anti-saturation VIs and soil background-eliminating VIs have been constructed,such as the optimized nonlinear vegetation index[15]and the optimized soil-adjusted vegetation index(OSAVI)[16].Several improvements have been made in estimation of AGB,PNC,and leaf area index(LAI)from image textural characteristics,vertical distribution rule,and plant height,or by combining synthetic aperture radar and spectral information,which corrects the phenological difference in the prediction [17-19].The uncertainties of the radiation transfer model were reduced by varying the model parameters over the growth stages [20].However,the robustness of these new techniques and methods is limited by cost,flexibility,data-processing difficulties,and image quality.

Deep-learning algorithms are widely used to predict agronomic traits with excellent results.A well-trained network depends on a wealth of training datasets that are expensive to collect,and their generality for image quality has not been thoroughly tested.To our knowledge,only two studies [13,14] have investigated automatically adapting the relationship between VIs and agronomic traits across growing seasons for AGB and plant area index.No study has yet focused on whether a unified statistical model can be developed to predict agronomic traits using VI,or whether a phenological correction is required for each agronomic trait across multiple growing seasons.

To respond this need,we focused on estimating agronomic traits for winter wheat throughout the entire growing season using six well-known VIs.The research objectives were(1)to investigate the relationship between VI and the agronomic traits during individual growth stages and the entire growth period;(2)to construct a VI-based adaptation approach to predict agronomic traits and assess the accuracy of such predictions;and (3) to evaluate the contributions of phenological information and VIs to AP models.

2.Materials and methods

2.1.Study sites and experimental design

Field experiments were conducted over four consecutive seasons from September 2017 to June 2021 at the Xiaotangshan National Precision Agriculture Research Center in a suburb of Beijing,China (40.17°N,116.43°E,mean altitude: 36 m).

Three experiments were performed in the 2017-2018 (Exp.1),2018-2019(Exp.2),and 2019-2020(Exp.3)winter wheat growing seasons [19] using a completely randomized design.Another experiment conducted in the 2020-2021 growing season (Exp.4)was designed as a three-way factorial arrangement of cultivars(P1: Jinghua 11 and P2: Zhongmai 1062),nitrogen fertilization(N1: 0 kg N ha-1;N2: 90 kg N ha-1;N3: 180 kg N ha-1,and N4:270 kg N ha-1) and irrigation (W0: 25 mm,W1: 148 mm,and W2: 241 mm) treatments in a split-plot design.Detailed designs were used for winter wheat (Table S1).

2.2.Traits measurement

Twenty wheat tillers per experimental plot were randomly sampled at different growth stages to measure agronomic traits(AGB,PNC,LAI,LGB,LNC,and SPAD).The SPAD was the mean value obtained from measuring a leaf of a 20-stem wheat plant using a SPAD-502 instrument.Table S1 presents the number of samples collected in each experiment.The measurement methods and equations for calculating AGB,LGB,PNC,LNC,and LAI followed Zhao et al.[21].

2.3.Field canopy hyperspectral data

Canopy reflectance spectra were collected from 2017 to 2021 at the Xiaotangshan experimental plots with a Field Spec Pro FR2500 spectrometer (Analytical Spectral Devices,Boulder,CO,USA).Measurements were made during cloud-free periods between 10:00 AM and 2:00 PM to minimize variations in illumination conditions.Before measuring the canopy reflectance,a white reflectance panel from Spectrolon (Labsphere,North Sutton,NH,USA)with diffuse reflectance properties was used to calibrate the spectral reflectance.The instrument was held 1.0 m above the canopy,and the view nadir and azimuth angles were both zero.For each plot,hyperspectral data were collected from one representative 1.0 m × 1.0 m sub-plot.To reduce the effects of field conditions,ten scans were acquired and averaged to represent the canopy reflectance of each plot.

2.4.Spectral vegetation indices

Six VIs were selected as good candidates for evaluating agronomic traits (AGB,PNC,LAI,LGB,LNC,and SPAD) (Table 1) and the correlations between agronomic traits and each VI were calculated.These vegetation indices are widely used and can be popularized on the ground,UAV,and satellites.

2.5.Phenological information

The phenological information used in this study was mainly from the digitized growth stage.The Zadoks stages (ZS) [27] days after sowing(DAS),day of the year(DOY),and effective accumulative temperature(EAT)[28]were used to quantify the phenological development process of winter wheat.

2.6.A method based on vegetation indexes to adapt the phenology for predicting agronomic traits

Given that the statistical model parameters for agronomic traits differed among the various phenological periods (Table 3),we bridged the VIs and the phenology and agronomic traits (AGB,PNC,LGB,LNC,LAI,and SPAD) to automatically adjust the VI-toagronomic trait relationships as a function of phenology.A hierarchical linear model (HLM) [29] was introduced to calculate the phenological correction.At level 1,the agronomic traits are linear in phenological stage.The function is characterized by a set of parameters,and the HLM is given by.

Table 1 Summary of spectral vegetation indices used in this study.

where P,β0j,and β1jare the agronomic trait(one of AGB,PNC,LGB,LNC,LAI,or SPAD),the intercept,and the slope of the linear model.

At level 2,these parameters vary across phenological stages,and the parameters of the first layer are automatically corrected by the phenological data of the second layer:

where βj(j=0,1) corresponds to β0and β1from the HLM model,respectively,γm0is the intercept,and γm1is the slope of each PV.PV is the set of phenological variables (ZS,DAS,DOY and EAT).

2.7.Ordinary least squares regression

To assess how phenology affects estimates of the agronomic traits,we constructed and compared conventional ordinary least squares (OLS) regressions using a piecewise model based on multiple growth stages.Regression equations of linear and exponential form were fitted and the best regression model was selected to permit an analysis of the evolving relationships of intercepts and coefficients with increasing growth stages.From these piecewise regression models,we may take a linear regression as an example,such that.

where P and ZS represent agronomic traits and the Zadoks scale[27].The slope coefficients k11,k12,k13,k14,k15,and k16are collectively referred to as k1,and the intercept coefficients k21,k22,k23,k24,k25,and k26are collectively referred to as k2.The exponential regression is k1×eP×k2(LGB and LAI),which is used to replace the linear piecewise function in formula (3).

2.8.Extended Fourier amplitude sensitivity test and sensitivity analysis

To estimate the contribution of VIs and the phenological information of the HLM,we applied an extended Fourier amplitude sensitivity test,which combines a Fourier amplitude sensitivity test with the Sobol algorithm to determine the sensitivity index to quantify how input variables depend on output variables[30].This procedure involved four steps:(1)a set of independent variables is generated by Monte Carlo sampling of distributions,generating a total of 100,000 datasets.(2)The predicted agronomic trait for each sample set of independent variables is calculated by substitution in the HLM model in a MATLAB 2016 (Mathworks,Inc.,Natick,MA,USA)program.(3)The sensitivity index of each independent input variable to the agronomic traits is computed with Simlab (version 2.2.1) software.The total variation is.

where Viis the variance associated with each independent variable of the HLM model,and Vijis the variance associated with the interaction between each pair of independent variables of the HLM model.The sensitivity index is derived from the decomposition of Eq.(4) by dividing the importance measures by var(Y):

where Siis also called the sensitivity index,and Sijis the sensitivity index of the interaction between the VI and PV.The sensitivity index of each independent input variable was computed.

2.9.Model evaluation

Pearson correlations (r) between the VI and agronomic traits was calculated for each growth stage.Differences between regression coefficient were tested by F test.Model performance was evaluated by the adjusted determination coefficient R2and the root mean squared error (RMSE).

3.Results

3.1.Correlation between agronomic traits and vegetation indexes

All VIs correlated weakly with AGB,with OSAVI showing the strongest correlation (r=0.18**),followed by EVI2 (r=0.11**) and PRI (r=-0.10**).The VI with the highest correlation with PNC was GNDVI (r=0.50**),followed by NDVI (r=0.48**),PRI(r=0.48),EVI2(r=0.46**),and OSAVI(r=0.41**)(Table 2).VIs correlated with LGB,LNC,LAI,and SPAD,with the strongest correlations shown by CIrededge(r=0.80**),PRI (r=0.78**),OSAVI(r=0.80**),and OSAVI (r=0.81**).Accordingly,OSAVI,CIrededge,GNDVI,PRI,OSAVI,and OSAVI were used for subsequent analysis.

3.2.Relationship between vegetation indices and agronomic traits in seven phenological stages

The ranges of AGB,LGB,PNC,LNC,LAI,and SPAD from seven growth stages (ZS 33 to ZS 80) in the entire dataset were 0.52-2 2.62 t ha-1,0-3.82 t ha-1,0.86%-3.95%,1.05%-5.07%,0-7.34,and 3.70-56.77,respectively (Fig.1).All vegetation indexes were strongly correlated with the agronomic traits,with correlation coefficients r in each ZS exceeding 0.59 (Table 3).

The prediction-model parameters of the agronomic traits(Table 3) in various growth stages reveal that the AGB prediction parameters,k1and k2,gradually increased with the advance of the phenological stage,ranging from 3.68 to 20.53 and -1.57 to 4.71,respectively.For PNC models in different growth stages,the parameters k1and k2had opposite trends,ranging from 1.01 to 6.61 and -2.51 to 1.42,respectively.Because the parameters of AGB and PNC differed significantly over the growth stages,so do the curves of VI in Fig.1.

Table 2 Correlations between agronomic traits and VIs.

Table 3 Statistical summary of OLS models applied to the various phenological stages.

3.3.Effect of phenological indicators on predicting plant and leaf related traits

Figure S1 shows the correlations between the parameters (k1and k2)of the prediction models for individual growth stages and selected phenological indicators (DAS,DOY,EAT,and ZS).All phenological indicators were strongly correlated with the parameters(k1and k2) of AGB and SPAD prediction models (|r| varied from 0.83 to 0.97).The correlation coefficients |r| of k1for the PNC prediction models with DAS,DOY,GDD,and EAT were 0.61,0.60,0.58,and 0.59,respectively,while k2of the PNC prediction models was weakly correlated with the various phenological factors.For other agronomic traits,such as LGB and LAI,both k1and k2did not show high correlations with the four phenological indicators.

Fig.2 shows the sensitivity of agronomic traits to VI and phenological stage.The sensitivity indices of the phenological information in the AGB and PNC prediction models were 0.81-0.86 and 0.66-0.73,respectively,whereas the LGB,LNC,LAI,and SPAD prediction models produced sensitivity indexes of 0.01-0.02,0.01-0.02,0.01-0.02,and 0.02-0.08,respectively.

3.4.Comparison of agronomic trait models with and without phenological indicators as influencing factors

The predictions of AGB and PNC with the phenological indicators as influencing factors were more accurate than those obtained without the phenological indicators,but no significant difference appeared among the predictions of LGB,LNC,LAI,and SPAD if thephenological indicators were considered(Fig.3).Ignoring the influence of phenological stage,the accuracy of the AGB model based on VI gave R2=0.08 and RMSE=4.70 t ha-1,values ower than those(R2=0.86-0.88 and RMSE=2.92-3.08 t ha-1) obtained with the AGB model combined with phenological information.The model incorporating phenological information gave higher prediction accuracy (R2=0.79-0.87,RMSE=0.42%-0.67%)than the PNC prediction model based on VI(R2=0.45,RMSE=2.16%).The ranges of R2for the prediction model of LGB,LNC,LAI,and SPAD were 0.79-0.83,0.66-0.86,0.75-0.83,and 0.68-0.70,respectively,and the RMSE ranges were 0.27-0.59 t ha-1,0.49%-0.66%,0.60-0.90,and 3.84-4.32,respectively.The HLM model combined with phenological factors showed high accuracy for the remote sensing prediction of agronomic traits.

Fig.1.Regressions of(a)AGB,(b)PNC,(c)LAI,(d)LGB,(e)LNC,and(f)SPAD on VIs at seven phenological stages.VI,vegetation index;AGB,aboveground dry biomass;PNC,plant nitrogen concentration;LAI,leaf area index;LGB,leaf dry biomass;LNC,leaf nitrogen concentration;SPAD,soil and plant analysis development;NDVI,normalized difference vegetation index;OSAVI,optimized soil-adjusted vegetation index;GNDVI,green normalized difference vegetation index;EVI2,enhanced vegetation index 2;PRI,photochemical reflectance index;CIrededge,red edge chlorophyll index;ZS,Zadoks’ stage.

Fig.2.Sensitivity index of the HLM model for four phenological indicators: (a) DAS,(b) DOY,(c) EAT and (d) ZS.AGB,aboveground dry biomass;PNC,plant nitrogen concentration;LAI,leaf area index;LGB,leaf dry biomass;LNC,leaf nitrogen concentration;SPAD,soil and plant analysis development;DAS,days after sowing;DOY,day of the year;EAT,effective accumulative temperature;ZS,Zadoks stage.

In agronomic trait prediction using VIs,HLM automatically adjusted for differences caused by different phenological information.Fig.4 shows the HLM verification results combined with ZS and VI.When 7.5 t ha-1of AGB and 2.5% of PNC were the thresholds for predicting AGB (RMSE=3.41 t ha-1) and PNC(RMSE=0.06%) by the OLS model,the result was a clear overestimate (AGB <7.5 t ha-1or PNC <2.5%) and underestimate(AGB >7.5 t ha-1or PNC >2.5%) as shown in Fig.4a and 4c.The HLM model solved this overestimation and underestimation when predicting AGB (RMSE=0.70 t ha-1)and PNC (RMSE=0.01%).The OLS model and HLM model based on LGB,LNC,LAI,and SPAD did not produce significantly dissimilar results for the verification sets.

Fig.3.Comparing the accuracy of the prediction model for different agronomic traits:(a)R2 and(b)RMSE.The RMSE units of AGB,LGB,LNC and PNC were t ha-1,t ha-1,%,and%,respectively.AGB,aboveground dry biomass;PNC,plant nitrogen concentration;LAI,leaf area index;LGB,leaf dry biomass;LNC,leaf nitrogen concentration;SPAD,soil and plant analysis development.

4.Discussion

4.1.How winter wheat phenological stage affects agronomic traits:The VI model

Winter wheat is a typical crop with a long growth cycle,large differences in AGB and PNC during the different growth stages,and nitrogen ‘‘dilution” effects,which in turn produce large differences in the time series of AGB and PNC [21,31].Figure S2 shows that the relationships between LGB and AGB or LNC and PNC were affected for multiple stages.A similar tendency appears in the correlations between VIs and AGB or PNC(Fig.1a,c).During the growth and development of wheat,its external morphology and internal distribution mechanism undergo a series of changes.From the tillering stage to the jointing stage to the flowering stage,the distribution centers of assimilation products are leaves,stems,and spikes [32,33].However,optical information cannot capture the process of nutrient accumulation of winter wheat during its growth period.

4.2.Response of vegetation indices to wheat growth stages and agronomic traits

This study expands efforts to widen the application of VIs for monitoring the agronomic traits of winter wheat with multiple phenological stages.In the VI-based AGB and PNC prediction models,k1varied over a large range (from 3.41 to 21.53 and 1.45 to 6.61,respectively).Yue et al.[34]and Zheng et al.[9]described the phenomenon of AGB and PN at different growth stages.Classical VIs are constructed mainly to reduce the influence of soil background,aerosol,and light,to reduce other noise,and to limit saturation in the satellite image,and less attention is devoted to accurately monitoring phenological effects.Optical data can be used to detect leaves,but their ability to detect plant stems and panicles is limited.As the growth stage advances,the contribution of leaves to canopy information of winter wheat gradually decreases.Therefore,the influence of the phenological stage must be understood to construct a VI that overcomes the phenological effect and thereby accurately estimates AGB and PNC.

4.3.Should we rely on phenological indicators for the agronomic trait model?

A comparison of the sensitivity indices shows that phenological information contributed to both AGB and PNC (Fig.2).In a universal model,VIs has no ability to estimate AGB and PNC over the entire growth period (from tillering to filling) because optical data cannot fully detect the vertical structure of the winter wheat canopy [35].Four more-accessible indicators,DAS,DOY,EAT,and ZS,were used in the present study to estimate the contribution of phenological stage to agronomic trait models.The use of the ZS scale as a continuous variable was highly effective for predicting AGB and PNC,regardless of region.Temperature is considered[19]to drive the evolution of growth stages,but differ in temperature for crop growth and development differ between regions.

DAS,DOY,and EAT are readily accessible[7,36,37]and are suitable for analysis in specific regions,whereas AGB and PNC prediction combined with ZS stage has greater potential for use in larger areas.It is thus feasible to use the phenological phase simulation module of a mature crop growth model to predict the regional ZS and then predict the regional AGB,PNC,and even yield.Information such as plant height and texture features that can represent phenological differences at various growth stages was extracted and good results were obtained for predicting AGB and PNC [8,9].However,on a regional scale,the current resolution of optical remote sensing restricts the extraction of wheat plant height and texture characteristics.

HLM offers potential advantages in solving the problem of nested data[38,39],but these advantages do not extend to the prediction of leaf-scale traits(Fig.4).In the HLM models with LGB,LAI,and SPAD,which were constructed by considering phenology,the phenology information exerted a negative effect (R2decreased and RMSE increased).Because optical remote sensing provides mainly leaf information,VIs could lead to acceptable results when analyzing the leaf aspects in the multi-phenology stage.However,both HLM and OLS overestimated SPAD for small SPAD values,indicating that the phenological correction did not affect SPAD monitoring.As Fig.1 shows,all small SPAD values were obtained from the filling period.With the transfer and transportation of materials,the SPAD at the filling stage is low,but the top spike of the canopy is still green.Optical remote sensing is a kind of passive remote sensing,which cannot completely penetrate the whole canopy of crops,and this characteristic will affect the accuracy of prediction of agronomic traits [35].Each layer of HLM is a linear model,and the nonlinear forms of the model need further investigation(Fig.4).In the future,the structure of the hierarchical linear model needs further changes to predict agronomic traits.

CRediT authorship contribution statement

Yu Zhao:Visualization,Writing -original draft.Zhenhai Li:Writing -review &editing.Guijun Yang:Writing-review &editing.Shaoyu Han:Data curation.Haikuan Feng:Visualization.Yang Meng:Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This study was supported by the National Key Research and Development Program of China (2019YFE0125300),the Shandong Provincial Key R&D Plan(2021LZGC026),and the China Agriculture Research System (CARS-03).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.08.003.