Cho Zhng ,Zi’ng Xie, *,Jili Shng ,Jingui Liu ,Tifeng Dong ,Min Tng ,Shoyun Feng,*,Hunjie Ci

a College of Hydraulic Science and Engineering,Yangzhou University,Yangzhou 225009,Jiangsu,China

b Agriculture and Agri-Food Canada,Ottawa Research and Development Centre,960 Carling Avenue,Ottawa,ON K1A 0C6,Canada

c Key Laboratory of Agricultural Soil and Water Engineering in Arid Area of Ministry of Education,Northwest A&F University,Yangling 712100,Shaanxi,China

Keywords:Time-series VI Asymmetric Gaussian function Phenological stage Shape model Remote sensing

ABSTRACT Accurate information about phenological stages is essential for canola field management practices such as irrigation,fertilization,and harvesting.Previous studies in canola phenology monitoring focused mainly on the flowering stage,using its apparent structure features and colors.Additional phenological stages have been largely overlooked.The objective of this study was to improve a shape-model method(SMM) for extracting winter canola phenological stages from time-series top-of-canopy reflectance images collected by an unmanned aerial vehicle (UAV).The transformation equation of the SMM was refined to account for the multi-peak features of the temporal dynamics of three vegetation indices(VIs) (NDVI,EVI,and CIred-edge).An experiment with various seeding scenarios was conducted,including four different seeding dates and three seeding densities.Three mathematical functions: asymmetric Gaussian function (AGF),Fourier function,and double logistic function,were employed to fit timeseries vegetation indices to extract information about phenological stages.The refined SMM effectively estimated the phenological stages of canola,with a minimum root mean square error(RMSE)of 3.7 days for all phenological stages.The AGF function provided the best fitting performance,as it captured multiple peaks in the growth dynamics characteristics for all seeding date scenarios using four scaling parameters.For the three selected VIs,CIred-edge achieved the greatest accuracy in estimating the phenological stage dates.This study demonstrates the high potential of the refined SMM for estimating winter canola phenology.

1.Introduction

Canola (Brassica napus L.),the third largest oilseed producing crop in the world,is a source of edible oil and feed protein and a feedstock for biofuel production[1].China is one of the top canola producers,accounting for 15.6% of total sown area and 17.4% of global production (FAOSTAT: https://www.fao.org/faostat/en).Winter canola accounts for about 90% of canola production in China,and is generally planted in rotation with rice,maize,soybean,and other grain crops for maximum farmland use efficiency.In 2019,the winter canola-sown area and yield accounted for 50.9% and 38.6% of the oil crop of China [2].Timely and accurate information about canola development is desirable for making field management decisions to optimize production.

Phenology information is essential for understanding seasonal crop growth dynamics and performing crop management [3].Insufficient water or inadequate light conditions during the flowering stage of canola can result in pollination failure[4,5],and excessive soil water during the rice tillering period can increase the number of ineffective tillers (inability to grow and ear forming)[6].Accurately identifying key crop phenological stages is essential for management activities such as irrigation scheduling,fertilization,and harvesting [7-10].Phenology dates are also necessary input parameters for crop growth models and remote sensing data assimilation [11-13].

Crop phenological monitoring has relied mainly on manual observation or automated recording,both of which are labor and cost-intensive.With rapid advances in remote-sensing technology,detection of regional crop phenological stages using time-series remote sensing data has been shown to be feasible and effective.The advantage of using remote sensing data is the effective capture of spatial-temporal variability of growth dynamics using sensors with varying revisiting frequencies and spatial resolutions [14].Observations can be categorized into frequency,resolution,and target based on sensor type.Numerous remote sensing techniques,including microwave [15],solar-induced chlorophyll fluorescence(SIF)[16],and synthetic aperture radar(SAR)[17],have been used to identify vegetation phenology.

There is no universally accepted method for determining the phenological stages of all crop types.Vegetation index(VI),an indicator by spectral calculation of two or more bands of light,is commonly used to highlight vegetative properties.The main optical remote-sensing-based phenology detection methods can be classified into three groups:1)threshold methods,which assume that a phenological stage starts when the VI reaches a specific value,including ‘‘fixed” and ‘‘dynamic” thresholds (e.g.,ratio,mean,and median of time-series VI) [18-20];2) change-detection methods,which detect the characteristic points on time-series VI curves,such as maximum,minimum,inflection points,and derivatives to determine phenological stage dates [21-23];and 3) geometrical model fitting,which aims not only to smooth time-series VIs but to align the derived VI curves with ground-observed phenological stages.The typical method used is the shape-model method(SMM),a two-step filtering approach initially proposed by Sakamoto et al.[24].The SMM assumes that the positions of crop phenology stage dates on the time-series VI curve are stable and that the VI curve can be linearly scaled to accommodate other growth dynamics.

Previous canola phenology studies focused mainly on detecting the flowering stage,owing to the apparent variation of canopy spectral signatures at this stage [25-28].Fang et al.[29] used VIs to estimate the flowering ratio of canola.Wan et al.[30]estimated the number of flowers of oilseed rape by combining VI and image classification.d’Andrimont et al.[4] monitored the peak flowering period of canola based on Sentinel-1/2 satellite remote sensing data with a temporal accuracy of 1-4 days.Han et al.[31] developed a new spectral index(Normalized Rapeseed Flowering Index)using Landsat 8 and Sentinel-1/2 data to identify canola peak flowering dates in Germany.However,accurate detection of canola phenological stages is challenging.In contrast to those of other crops,the temporal profiles of VIs of winter canola show shortterm sharp fluctuations with multiple peaks[32].Delaying seeding for winter crops also shortens the growth period [33].These variations may render threshold and change-detection methods inapplicable.

The SMM has been successfully applied to a wide variety of crop types (including wheat,maize,soybean,and barley) by capturing specific phenological events throughout the growth period [34].Sakamoto et al.[24] used the SMM to estimate the phenological stages of maize and soybean using MODIS (Moderate Resolution Imaging Spectroradiometer) remote sensing data,achieving the minimum root mean square error(RMSE)of 2.9 days for the silking stage.Zeng et al.[35] incorporated a temperature response function into time-series VI to estimate the phenological stages of maize and soybean.Sakamoto [36] improved the ability of the SMM to detect and estimate phenological information for various crops by refining the calibration procedure of SMM.Zhou et al.[37] coupled cumulative temperature with a shape model using a double logistic function to predict winter wheat phenology and obtained the RMSE,ranging from 1.0 day at maturity to 10.3 days at tillering stages.Yang et al.[38]used the SMM for estimating rice phenology,and achieved an RMSE of 4.0 days (68.8% of the mean)and a mean absolute error of 0.32 day.The satisfactory accuracy achieved by these studies suggests that the SMM may be generic for detecting crop phenology with simple seasonal dynamics of VIs(e.g.,one peak).Few studies address extracting phenology from complex curves.

To our knowledge,no studies have addressed the estimation of winter canola phenological stages using the SMM.Most previous phenological monitoring studies were performed in the same region,where the seeding date was consistent locally,resulting in similar crop growth dynamics.The versatility of applying SMM to a canola crop,whose growth is highly variable is unclear.

The main objective of this study was to develop a refined shape model for time-series VI to estimate winter canola phenological stages.A modified SMM was used to account for the features of the temporal dynamics of canola(such as bimodality).Time series of three typical VIs and three filtering functions were employed to fit the growth curves represented by the seasonal VIs.

2.Materials and methods

2.1.Study area and experimental design

The study was performed at an ecological experimental station,located in the Jianghuai Plain of eastern China (32°21′N(xiāo),119°24′E,5 m above sea level).This region is characterized by a subtropical monsoon climate,with a mean annual air temperature of 14.8 °C and precipitation of 1063 mm.The soil texture of the experimental field is loamy with an average field capacity of 0.35 cm3cm-3,permanent wilting point of 0.13 cm3cm-3,and dry bulk density of 1.49 g cm-3.The available carbon,nitrogen,phosphorus,and potassium in the 0-20 cm surface layer were 10.2 g kg-1,0.97 g kg-1,16.3 mg kg-1,and 151.2 mg kg-1,respectively.

A field experiment was designed and conducted during the 2020-2021 growing season,with 12 seeding scenarios.Two factors,seeding date and density,were considered in the experiment.For the seeding date,four time periods: early (ES,September 21),normal (NS1,October 6;and NS2,October 21),and late (LS,November 6) seeding were administered according to the local common seeding calendar.Three seeding densities:low density (LD,12.5 × 104plants ha-1),medium density (MD,25 × 104plants ha-1),and high density (HD,37.5 × 104plants ha-1) were adopted for each of the four seeding dates.Thus,the experiment consisted of 12 treatments with three replicates,totaling 36 plots.

The plot size was 5 m × 5 m,consisting of 13 plant rows and a 2-m buffer between adjacent plots.The experiment followed a split-plot design,with seeding date as the main plot and seeding density as the subplot.Basal fertilization at the rates of 180 kg N ha-1,90 kg P2O5ha-1,120 kg K2O ha-1,and 15 kg B ha-1was applied before seeding.All seedlings were thinned manually by the 5-leaf stage to achieve the designed density.Other field management such as weed,pest,and disease control followed local practices.

2.2.Data collection and preprocessing

2.2.1.Ground-based observations of canola phenology

Days after seeding(DAS)and ground-based phenology observations of winter canola for all plots were simultaneously recorded multiple times throughout the growing season.In general,the growth cycle of canola can be divided into several phenological stages,as defined by the Biologische Bundesanstalt,Bundessortenamt and CHemical Industrie (BBCH) scale [39],including leaf development,stem elongation,inflorescence emergence,flowering,development of fruit,and ripening.Crop development is strongly dependent on the environment,with high variation among geographical regions.In China,winter canola is usually seeded in September and October and undergoes a lengthy (>2 months) leaf and stem development period.The growth can be further divided into six phenological stages: seedling,overwintering(withered-leaf stage),bolting,flowering,podding,and ripening stages[32,40].The DAS and phenological stages of the experiment are summarized in Table 1,and the Red,Green,and Blue (RGB)composites for various phenological stages are shown in Fig.1.

Table 1 Observation of winter canola growth development from the early-seeding scenario.

2.2.2.Spectral image acquisition and preprocessing

Aerial images covering all plots during winter canola growth periods were acquired with a multispectral sensor (Rededge-MX,MicaSense Inc.,Seattle,WA,USA) mounted on a DJI Inspire 2 unmanned aerial vehicle (DJI Inc.,Shenzhen,Guangdong,China).The Rededge-MX camera features five sensors with a 1280 × 960 pixel imaging frame size and a field of view of 47.2°.The five sensors acquire images in five spectral bands (blue,green,red,rededge,and near-infrared),with central wavelengths (bandwidths)of 475 nm (32 nm),560 nm (27 nm),668 nm (14 nm),717 nm(12 nm),and 842 nm(57 nm).An image of a calibrated reflectance panel was captured immediately before and after each flight,to be used to convert raw pixel values to reflectance.The flight path was set at 75% overlap for both along-track and cross-track directions to ensure adequate coverage and post-flight mosaicking.The flight altitude was set to 20 m above the ground,with the camera sensors pointing vertically downwards to achieve a spatial resolution of 1.4 cm.The exposure time was adjusted automatically in response to the irradiance detected by a down-welling light sensor (DLS)mounted on top of the UAV.The flights were carried out at solar noon under cloudless and windless conditions.The UAV flights were repeated every 5 to 10 days throughout the entire growing season.

After each flight,calibrated and georeferenced reflectance spectral images were generated by radiometric correction and geometric correction.An orthomosaic image covering the entire experimental area was generated and each plot’s region of interest(ROI) was manually defined.To avoid interference from edge effects,the edge area of each plot was excluded from subsequent analysis.The plot-based reflectance spectra of the five bands were extracted from the ROI for VI calculation.

2.3.Methodology of the shape-model method

2.3.1.Time-series vegetation indices

VIs are commonly used to describe vegetation growth conditions,and an algebraic combination of two or more spectral bands can enhance signatures of target vegetation properties while suppressing signals from other sources [41].Three representatives of VIs were chosen.The normalized difference vegetation index(NDVI) reveals the greenness of vegetation and is commonly used for retrieving leaf area index (LAI),photosynthetic capacity,chlorophyll,and biomass [42].The enhanced vegetation index(EVI) is similar to the NDVI but has been devised to improve sensitivity at high plant biomass and to suppress canopy background signal and atmospheric influences [43].Red-edge chlorophyll index (CIred-edge),using red-edge reflectance,was developed primarily for estimating chlorophyll content and indicating photosynthesis [44].

2.3.2.Building shape models

Seasonal time-series VI profiles of plants reveal their phenological development.In contrast to those of wheat,maize,and soybean [24,37],the VI curves of winter canola present multiple peaks throughout the growing season (Fig.3).The shape of timeseries VI is directly related to the seeding date.The flowering stage presents a unique feature on the time-series VI curve: the VI continuously decreases until the canola reaches the flowering stage and then drops to a minimum at the blooming stage,behavior that can be easily detected using the threshold method[27].In the present study,the phenological stages from seeding to flowering were the targets.As shown in Fig.2,time-series VI curves before the full flowering period exhibit double peaks that have not been reported in the literature(description and explanation are presented in section 3.1 and 4.1).Based on this signature,the two mathematical models of asymmetric gauss function(AGF)[45]and Fourier series function (FF) [46] were used to match the seasonal VIs dynamics.Because the double logistic function (DLF) has been widely used in fitting single-peak profiles[47],it was also applied.The formulas are as follows:

where DAS is days after seeding.For the AGF model,j and k determine the amplitude of a VI curve,f and h are the central locations of peaks,and g and i are the widths of the peaks.For the FF model,l is a constant associated with the level of base bare soil VI value and p is the fundamental frequency of the signal.For the DLF model,b and d are related to the maximal variation rate of VI for the vegetative stage (increase) and reproductive stage (decrease),and a and c are related to the times of the maximum rates.

To generate a shape model able to represent the majority of the cases of the time-series VI,the data from the early seeding to flowering stage with three density treatments contained the most intensive observations and were thus used as a calibration dataset for building the shape model.Thus,a total of nine plots with observation data collected during 26 field surveys,from seeding to flowering,were used to construct the canola shape model.The remaining three seeding dates of the 27 plots were used for model evaluation.

2.3.3.Transforming the shape model to target VI curves

For the shape model,the shape of the VI curve is assumed to be crop-specific and unaffected by environmental conditions and management practices [38].The shape model is not identical to the fitting function,but rather a linear transformation,which allows the model to be scalable and flexible to fit any time-series VI profile.Two aspects are addressed by the linear transformation.First,as mentioned above,the time-series VI exhibits a double peak pattern;however,the shape of the two peaks is variable(described in Section 3.1),necessitating the modulation of each peak independently.The canola growth cycle is shortened when the seeding date is delayed.Consequently,the x-direction (x: DAS) of the time-series VI needs to be stretched or compressed by the linear transformation.

Fig.1.RGB images showing phenological stages of winter canola in the experiment.

Fig.2.Simplified schematic description of the procedures of the SMM.The horizontal black arrows represent the modulated direction of the parameters xscale and tshift.The vertical red arrows represent the controlled direction of the parameters yscale1 and yscale2.

To transform the shape model into a target VI curve,a refined transformation equation was proposed as shown in Eq.(4).In contrast to the original function with three parameters (xscale,yscale,and tshift)suggested by Sakamoto et al.[24]and Zeng et al.[35],a parameter in the y-direction was added in this study to enable the transformation equation to regulate the shape of each peak.Fig.2 illustrates the schematic procedure of the SMM and the functions of each parameter.The xscale and tshift parameters compress and translate the shape model in the × direction.The yscale1and yscale2parameters adjust the magnitude of the VI.The h (x)denotes the shape model,while the subscript h1,2(x) represents two components of the model.The function g (x) refers to the shape model h (x) geometrically transformed in the × (DAS) and y (VIs) directions using these four scaling parameters.Scaling parameters can be determined by matching the established shape model to the target VI curves(fitted by DLF,AGF,and FF)based on the RMSE (Eq.(5)).The search ranges for each parameter were defined as follows: 0.1

where h (x) is the established shape model and h1,2(x) represents two parts of the model,namely h (x)=h1(x)+h2(x).In the AGF model,

3.Results

3.1.Seasonal variation of VIs under different seeding date scenarios

Great differences in the dynamics of the three VI curves were observed among the four seeding date scenarios (Fig.3).The growth period was reduced from 230 to 185 days with the delay of the seeding date (about 15 days in the seeding interval).The number of days from seeding to flowering stage decreased from 175 to 140 days.Canola plants,in general,have three major components during the growth cycle: leaves,flowers,and pods.The growth dynamics,using NDVI as an example,can be described below.

NDVI rose dramatically following seed germination in the early seeding scenario (ES) and peaked around DAS 30.It began to decline and continued until the end of the overwintering stage(around DAS 120).In the spring,winter canola began to green up again,and the NDVI increased to the second peak,a sign showing the beginning of the flowering period(DAS 160).With the progression of the flowering stage,NDVI dropped to a minimum at full bloom,and then increased again as the flowers faded and pods started to grow.This variation pattern is similar to that reported in Shen et al.[27]and Behrens et al.[48].Although all seeding scenarios in our study showed consistent growth dynamics with three peaks in VI curves,apparent differences in geometrical patterns were observed.With the seeding date delayed,VI decreased significantly after the first peak,to levels even lower than those observed at the later stage.Compared with the NDVI,EVI and CIred-edgecan alleviate the saturation issue at high LAI values.The data of the EVI and CIred-edgeare more scattered,which is advantageous for diagnosing multiple seeding densities.

Fig.3.Dynamics of three VIs during the canola growing season with four seeding dates and three density scenarios.The ES,NS1,NS2,and LS represent early seeding,two normal seeding and late seeding,respectively.The LD,MD,and HD represent low density,medium density,and high density,respectively.(Similarly hereinafter).

3.2.The shape models with the three functions

Three functions,AGF,FF,and DLF,were used to capture the variation of three VIs from the seedling to flowering stages.Results of optimal fitting indicated that all three fitting functions were capable of fitting the dynamics of the time-series VIs(Fig.4).For the NDVI and EVI,although all fitting functions resulted in excellent accuracy with RMSE <0.1,only AGF and FF captured the bimodal feature of the VI curves.In comparison to NDVI and EVI,the bimodal profile was more pronounced on the CIred-edge,greatly reducing the interpretation capability of the DLF.For the CIred-edge,the AGF performed admirably,with an R2of 0.92 and an RMSE of 0.1.The FF was ranked second,and the amplitude of the fit curve was slightly less than that of AGF.The DLF model gave the lowest R2and RMSE values of the three models.The parameters of each function [Eqs.(1),(2),and (3)] were determined by optimal fitting,and the shape model was prepared for transformation by Eq.(4) to describe the changes caused by differences in seeding date.

Fig.4.Establishment of the shape models for three VIs based on three fitting functions.The AGF,FF,and DLF represent asymmetric Gaussian function,Fourier function,and double logistic function,respectively.

3.3.Comparison of the functions for transforming VI curves

According to the performance of three shape models for canola growth matching,only the AGF method preserved the shape signature for all the target curves,despite a few cases of suboptimal accuracy (Fig.5).In most cases,AGF has a greater capability of revealing VI variability associated with phenological changes.The inadequate results of AGF occurred mainly in situations where the amplitude difference between the peaks and troughs was small,such as NDVI (NS1 and NS2) and EVI (NS1 and NS2).Compared with the AGF,the FF showed inconsistent performance.The FF method retained shape features only for the second seeding date scenario(NS1),owing to the similarity of crop growth dynamics of the first seeding date scenario(ES).Still,the FF was unable to capture the dynamics of the time-series VIs,especially the second peak.The DLF method was incapable of capturing the VI dynamics.We note that the fitting accuracy was not exactly equal to the estimation accuracy of the DAS and phenological stage dates,as the estimated value was calculated by Eq.(6) and not retrieved from the fitting function.

Fig.5.Performances of shape models for describing time-series VIs under three fitting functions.The NS1,NS2,and LS represent two normal seeding and late seeding,respectively.The AGF,FF,and DLF represent asymmetric Gaussian function,Fourier function,and double logistic function,respectively.

Table 2 summarizes the optimal scaling parameters of the transformation function (Eq.(4)).As mentioned previously,the parameters xscale and tshift primarily modulated the shape model in the time domain,whereas yscale1,2focused on the amplitudes of VI.Taking the result of CIred-edge-AGF as an example,the values of xscale were slightly >1,indicating that the growth cycles of the three later seeding date scenarios were shorter than those of the first seeding date scenario,so that an elongation was required to maintain the shape model’s time span.Positive values for tshift indicated that all shape models must be shifted to the right to approach the target curve.The values of tshift increased gradually for the three later seeding scenarios,denoting that the total growth periods were shortened further with the delay of seeding date.The values of yscale1were all <1 and decreased successively,suggesting that crop growth was influenced by delayed seeding date,so that the canopy could not reach the level of the first seeding date scenario before overwintering.By contrast,the values of yscale2were all >1,indicating that the crop developed more for the later seeding scenarios after the overwintering period.

Table 2 Optimum scaling parameters used in fitting functions.

Table 3 Accuracy assessment of the estimated phenological stage dates against the ground-based observations by the root mean square error (RMSE) for three VIs in three fitting functions.

3.4.Evaluation of phenological stages and DAS estimation

According to the mean accuracy of the estimated phenology,the shape model based on the AGF achieved the highest accuracy,with RMSE ranging from 2.4 to 6.8 days at different growing stages(Table 3).Thus,maintaining the shape features of the growth curve is important for phenology matching and detection even though the parameters of yscale1,2were not involved in the phenological calculation (Eq.(6)).As shown in Fig.6,the AGF model performed well throughout the growth period.The FF model displayed considerable variability in estimating DAS for different seeding scenarios,with RMSE ranging from 3.8 to 142.9 days.Adequate accuracy was obtained only for the NS1 seeding scenario,owing to similar dynamics of the ES scenario.In contrast,poor performances were observed in NS2 and LS scenarios,with large underestimation of the simulated DAS.The DLG model produced unacceptable results for all seeding scenarios,with the maximum RMSE 158.6 days.The results indicated that the AGF model,which is good at revealing crop growth variation,has a greater capability of estimating crop DAS and phenological stage dates than that of the FF and DLF model.

Fig.6.Comparison of the DAP and periods of canola in three VIs and the fitting functions used between ground-based observation and UAV-derived estimation.The AGF,FF,and DLF represent asymmetric Gaussian function,Fourier function,and double logistic function,respectively.

Results of the global estimation accuracy evaluation showed no significant differences were observed among the three VIs when using the AGF-based model(Fig.7).The greatest estimation errors for all three VIs appeared at the bolting stage.Comparable performance was observed in NDVI and EVI,with the respective RMSEs of all stages being 4.6 and 4.4 days.CIred-edgeshowed slightly higher estimation accuracy for the phenology stages (seedling,RMSE 2.8 days;overwintering,RMSE 3.7 days;flowering,RMSE 3.2 days;all stages,RMSE 3.7 days)but not the bolting stage(RMSE 5.4 days).

Fig.7.Comparison of the global estimation accuracy of three VIs and fitting models.SD,seedling stage;OW,overwintering stage;BT,bolting stage;FL,flowering stage;Whole,entire growing period of canola.The AGF,FF,and DLF represent asymmetric Gaussian function,Fourier function,and double logistic function,respectively.

4.Discussion

4.1.Effects of seeding date on crop phenology

Seeding date is the most important factor determining a crop’s growth process if abiotic and biotic stresses are not considered.Our study revealed that seeding date exerts a large effect on crop growth dynamics(Fig.3).For the time-series VIs curves,crops like wheat,maize,soybean,and rice show a typical single peak and plateau [22,38,49].Ma et al.[50] and Zhou et al.[37] suggested that multiple peaks were abnormal fluctuations and noise that should be eliminated or abandoned.However,with the similar growth dynamics proposed by Behrens et al.[48],all the time-series VIs curves for canola in this study displayed typical multi-peak features.The multi-peak situation is attributed to the effect of air temperature on crop growth and spectral response to crop canopy change.

With delay in seeding date,temperatures before wintering gradually decreased,resulting in a slowing of canola growth rate.Within 0-60 DAS following the first and second seeding date scenarios (ES and NS1),the daily mean temperature was approximately 20 °C,leading to apparent rapid growth of VIs(Fig.3).In contrast,in the third and fourth seeding scenarios,the temperature dropped to about 5 °C,slowing canola growth,and the canopy was unable to reach the same growth level as that of the ES and NS1 scenarios.Similar observations were also reported for other winter canola studies [32,40],noting that winter canola should enter a ‘‘dead leaf” period during the winter.In contrast to other crops (e.g.,wheat,maize,and soybean) that reaching the higher VI values at reproductive growth stages,canola has conspicuous flower petals that dominate the canopy for about 30 days [51].Yellow petals contribute to the reflectance in green and red bands during the flowering stage by absorbing blue light(～450 nm) from carotenoids [26,52,53],resulting in the reduction of VI values.

Comparing seeding scenarios,total growing days were decreased by around 45 days(Fig.3).Studies in grain crops showed that late seeding could shorten the growth period[33,54].This was demonstrated for canola by Lilley et al.[55],as well as the present study.Moreover,the effects of late seeding on crop growth were more pronounced as cold soil further restricted the growth of the small plants,even though comparable values of the three VIs were achieved for the seeding scenarios when entering the reproductive growth stage (flowering stage).

4.2.Performance of the enhanced shape model for fitting canola growth dynamics

Considering the distinctive signatures of temporal VIs,three filtering methods and a refined transformation equation were incorporated to estimate canola phenological stages.The modified SMM effectively established a specific growth curve that can be used to fit the target VI curves with a set of optimal parameters.The average estimation accuracy for all phenological stages based on CIred-edge-AGF was 2.8 to 5.4 days(Fig.7).In comparison,Sakamoto et al.[24]used the SMM to estimate maize and soybean phenological stage,with RMSEs ranging from 2.9 to 7.0 days and 3.2 to 6.9 days,respectively.Zeng et al.[35] modified the SMM estimate maize and soybean phenological stages and obtained a mean RMSE of 1.9 to 4.3 days for maize and 1.9 to 4.9 days for soybean.Zhou et al.[37] coupled cumulative degrees with the shape model of a double logistic function for predicting winter wheat phenology,with minimum RMSE ranging from 1.0 day at maturity to 10.3 days at tillering stage.Compared with the flowering monitoring studies by d’Andrimont et al.[4](RMSE=1-4 days),Han et al.[31](correlation coefficient >0.7),and McNairn et al.[56] (correlation coefficient=0.91-0.96),this study showed comparable accuracy,with RMSE ranging from 3.2 to 4.6 days (AGF-based).The shape model ensured acceptable accuracy not only for the flowering stage but also for other stages.

The accuracy of the SMM in estimating phenological stages was determined by matching the geometrical shape features of the temporal VIs series rather than fitting accuracy.For example,in fitting of the EVI target curve for the NS1 scenario(Fig.5),the R2and RMSE of AGF-based fitting were respectively 0.78 and 0.11,values lower than those of the FF (R2=0.84;RMSE=0.08) and DLF(R2=0.88;RMSE=0.06) functions.Although the AGF-based curve still retained the characteristics of the VI development,the FF and DLF failed.Finally,the global estimation accuracy of the AGF-based function achieved the RMSE of 4.4 days,considerably superior to that of others (RMSE >60 days) (Fig.7).

4.3.Comparison of vegetation indices for estimating phenology

Comparable results were obtained for the three VIs based on the AGF function,corroborating the assertion of Sakamoto et al.[24]that the SMM paid more attention to the geometrical features of crop growth than to the values of the VI curve.Zhou et al.[37]evaluated the capability of a temperature-based shape model in five VIs for identifying wheat phenology and confirmed that the method played a dominant role in prediction accuracy rather than the type of VI used.For the estimation accuracy in this study,each VI gave wins and losses at different phenological stages and under different seeding scenarios(Table 3).The estimation accuracies for all stages using the CIred-edgewere slightly greater than those of NDVI and EVI,as time-series CIred-edgehad more noticeable changes than those of NDVI and EVI,owing to a tendency toward less saturation at high biomass (Fig.3).Gitelson et al.[44],Gong et al.[57],and Ma et al.[50] concluded that the CIred-edgewas appropriate for assessing photosynthesis capability and could be used for monitoring crop growth.Thus,the curve features were highlighted when winter canola growth progressed into the overwintering and flowering period during which the crop photosynthetic rate decreased.

Fig.8.Changes in estimated DAS with different observation frequencies.The top and bottom range of the boxes represent percentiles 25 and 75 of samples.The box’s whisker represents the outlier range,and stars are extrema.The line and square within the box represent the median and mean values,respectively.

Other indices have been used in the literature to reduce interference or to highlight specific features.Sakamoto et al.[24],Gitelson [58],and Yang et al.[38] preferred the use of the Wide Dynamic Range Vegetation Index(WDRVI)over the NDVI to detect crop characteristics.Sulik and Long[25]and d’Andrimont et al.[4]suggested that the Normalized Difference Yellow Index (NDYI)could effectively detect canola peak flowering dates owing to its sensitivity to increased yellowness.Han et al.[31] reported that the Normalized Rapeseed Flowering Index (NRFI) outperformed the NDYI in catching peak flowering as soil water contents decreased,and canola during the flowering stage was observed in the short-wave infrared band (2100-2300 nm).

Considering that this study relied on high-temporal-frequency observations,for which data may not be available owing to unfavorable weather conditions and longer revisit cycles of some satellites,e.g.,Sentinel-2 (5 days) and Landsat 8 (16 days),the performance of estimated DAS under different observations was depicted in Fig.8.To capture the uncertainty associated with the observations,200 samples of 15,10,and 5 observations from the original data (ES scenario) were randomly extracted.As observation frequencies decreased,the median of the estimated RMSE for the seeding scenarios increased,and the distributions of the three assumed frequencies expanded gradually.This result was in agreement with Yang et al.[38],who reported that the performance of an SMM degraded when fewer observation data were used.For low frequency (e.g.,5 times),the estimation accuracy of DAS was determined by the observation period.The optimal accuracy of RMSE was as high as 1 day for different seeding scenarios when the shape model could match the observation series.In contrast,when the observation series were concentrated or scattered,with poor matching to the shape characteristics of crop development,errors could exceed 100 days.

5.Conclusions

We introduced a refined shape model method (SMM) for detecting winter canola phenology and tracking days after seeding(DAS) based on time-series UAV remote sensing data.We found that the refined SMM could be used to estimate canola phenology effectively,with the minimum error of RMSE (3.7 days) for all stages.In a comparison of the performances of the asymmetric Gaussian function (AGF),Fourier series function (FF),and double logistic function (DLF) for interpreting time-series VIs of canola,the AGF function provided the highest ability to describe multipeak signatures for all seeding date scenarios with four scaling parameters,whereas the FF and DLF models were not useful.Comparable performance of AGF-based function was observed for NDVI and EVI,with RMSE of 4.5 days.The CIred-edgeshowed a slightly higher estimation accuracy for the phenology stages (seedling,RMSE 2.8 days;overwintering,RMSE 3.7 days;bolting,RMSE 5.4 days;flowering,RMSE 3.2 days;all stages,RMSE 3.7 days).

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.

CRediT authorship contribution statement

Chao Zhang:Investigation,Writing -original draft.Ziang Xie:Conceptualization,Methodology.Jiali shang:Formal analysis,Writing -review &editing.Jiangui Liu:Formal analysis,Writing-review &editing.Taifeng Dong:Methodology,Validation.Min Tang:Writing -review &editing.Shaoyuan Feng:Supervision,Funding acquisition.Huanjie Cai:Conceptualization.

Acknowledgments

This research was supported by the National Natural Science Foundation of China(51909228),the Postdoctoral Science Foundation of China (2020M671623),and the ‘‘Blue Project” of Yangzhou University.The authors would like to acknowledge the Agricultural Water and Hydrological Ecology Experiment Station,Yangzhou University,for providing the experimental facilities.