Assessment of the Rice Panicle Initiation by Using NDVI-Based Vegetation Indexes

: The assessment of rice panicle initiation is crucial for the management of nitrogen fertilizer application that affects yield and quality of grain. The occurrence of panicle initiation could be determined via either green ring, internode-elongation, or a 1–2 mm panicle, and was observed through manual dissection. The quadratic polynomial regression model was used to construct the model of the trend of normalized difference vegetation index-based vegetation indexes (NDVI-based VIs) between pre-tillering and panicle differentiation stages. The slope of the quadratic polynomial regression model tended to be alleviated in the period in which the panicle initiation stage should occur. The results indicated that the trend of the NDVI-based VIs was correlated with panicle initiation. NDVI-based VIs could be a useful indicator to remotely assess panicle initiation.


Introduction
Rice (Oryza sativa L.) is one of the most important staple foods for more than half of the world population. Due to the rapid growing of food demand and limited arable land, improving yield potential to boost up future rice production is an urgent need. Rice yield is known to be increased by the nitrogen topdressing at the panicle initiation (PI) stage that is the beginning of the reproductive stage [1][2][3]. According to the recommendation of Agricultural Improvement Committee in Taiwan, when the length of the panicle is found to be 2 mm, the nitrogen fertilizer should be applied within two days. Applying a large amount of nitrogen fertilizer before PI would easily cause excessive stem elongation and thus tend to increase lodging risk. On the contrary, applying nitrogen fertilizer after PI is less effective on rice yield improvement. Therefore, accurate determination of the PI stage is crucial for rice production.
The PI is generally considered as the turning point between the vegetative phase and the reproductive phase. When a rice plant has reached maximum tillering, the internodes of the rice stem start elongating and subsequently panicle initiating. The overlapping period between maximum tillering and PI is also termed as vegetative-lag phase. The differences of the rice appearance between the vegetative-lag and PI are obscure to the naked eye. Therefore, they are difficult to be distinguished directly by human observation.
The general methods to assess PI in the farm are identifying the internode-elongation and green-ring [4,5], length of young panicle (1-2 mm) in the cross section of dissected stem [5], and leaf number index/leaf appearance [6]. When 30% of the main culms have panicles 2 mm or longer, it is considered as the panicle differentiation (PD) stage [3] that is late for nitrogen topdressing. Those methods are inconvenient and inefficient for largescale estimation. Some other convenient, large-scale, and non-destructive approaches that help to monitor the plant growth stages, such as modified-calendar days and heat units, are potential candidates to be used for the precise estimation of PI. The calendar days method is the easiest approach but is not reliable enough as it is largely affected by the weather variability during the cultivation period [7][8][9]. Growing degree day (GDD) is an excellent heat unit that has been widely applied in corn production [10] since it was first proposed to describe the timeline of biological development [7]. However, these methods are rudimentary and are not able to distinguish the variability between fields. Nowadays, quantitative assessment, real-time and site-specific management of precision agriculture are the objectives of community.
Spectral remote sensing is another potential approach for the estimation of various variables that are correlated to plant architecture and physiology, providing high-throughput information non-destructively and rapidly for precision agriculture. The arithmetic combinations of vegetation spectral reflectance, which usually termed as vegetation indexes (VIs), became useful indicators for studying plant health and status. The NDVI (normalized difference vegetation index), which is calculated through a normalization procedure [11], seems to be the most popular and long-established VI. The NDVI is sensitive to responses on the green vegetation [12] as it correlates with the biophysical and physiological changes in plants [13][14][15][16][17][18]. Moreover, the NDVI has been subsequently developed to the other NDVI-based VIs, such as NDRE (normalized differential red-edge), GNDVI (green NDVI), NDSI (normalized difference spectral index), etc. [19][20][21]. Recent studies put efforts on the quantification of nitrogen contents in plants through remote sensing to optimize the fertilization efficiency and increased the yield [22][23][24][25][26]. In the field of rice research, NDVI is one of the most frequently used vegetation indexes. For example, preceding studies revealed that the NDVI is useful for in the research of rice breeding, nitrogen use efficiency monitoring, and rice yield prediction [27][28][29]. Furthermore, NDVI has also been used to monitor the rice growth stages, including the panicle development stage [30]. The preceding studies revealed the trend of NDVI changes during the rice growth cycle; however, the evaluation of PI prediction is currently unavailable.
The main objective of this study, therefore, is to investigate the relationship between NDVI and PI occurrence. Meanwhile, a non-destructive and high temporal resolution approach was also expected to be established for PI assessment in this study.

Field Experiment Design and Management
Pot experiments were conducted in 2019 at the experimental field (22 • 38 N, 120 • 36 E) of National Pingtung University of Science and Technology, Pingtung, Taiwan, to examine the canopy reflectance behavior as a function of panicle initiation on rice (Oryza sativa L.). A japonica cultivar, Kaohsiung147 (KH147), was planted into 4 groups with 5 replications. Ammonium sulphate was applied at a rate of 150 kg/ha N in each group (20% for basal, 20% for 1st tillering topdressing, 30% for 2nd tillering topdressing, and 30% for panicle initiation topdressing).

Determination of PI through Dissection
The PI stage was determined through manual dissection and features observation. The entrance of the PI stage was verified via either the green ring ( Figure 1A), internodeelongation, or when 1-2 mm panicle was observed ( Figure 1B). The ideal timing for nitrogen topdressing is the interval of Figure 1A,B. When the length of the panicle was over 2 mm, it was considered as the panicle differentiation (PD) stage ( Figure 1C) that was not included in PI determination. The PI occurrence of each group was determined if more than 50% of dissected samples were verified to PI.

Spectral Measurement and NDVI Calculations
The spectral data of the rice canopy was collected by using SpectraPen SP-100 (range 640-1050 nm, 2 nm scan-range, Photon Systems Instruments, Czech), which is a non-imaging, handheld hyperspectral sensor. The data collection period was between 10 DAT (days after transplanting) and 80 DAT. The acquisition time was between 10:00 a.m. and 2:00 p.m., and the sensor was held horizontally at nadir view in the position about 3-5 cm above the highest leaf of the plants. The integration time of spectrum collection was set to auto-sensitivity to minimize the interference of sunlight intensity variability. The collected spectrum data were used to calculate NDVI-based VIs, which can be expressed as follows: where the λa and λb respectively denote the reflectance of near-infrared (NIR) and red wavelengths.

Reference Wavelengths and Recombined NDVI-Based VIs
The key wavelengths were selected by some physiologically related reference wavelengths with the equation of NDVI that have been used in previous studies (Table 1). Those reference wavelengths were recombined and recalculated as the normalized difference procedure with all samples.

Spectral Measurement and NDVI Calculations
The spectral data of the rice canopy was collected by using SpectraPen SP-100 (range 640-1050 nm, 2 nm scan-range, Photon Systems Instruments, Czech), which is a nonimaging, handheld hyperspectral sensor. The data collection period was between 10 DAT (days after transplanting) and 80 DAT. The acquisition time was between 10:00 a.m. and 2:00 p.m., and the sensor was held horizontally at nadir view in the position about 3-5 cm above the highest leaf of the plants. The integration time of spectrum collection was set to auto-sensitivity to minimize the interference of sunlight intensity variability. The collected spectrum data were used to calculate NDVI-based VIs, which can be expressed as follows: where the λ a and λ b respectively denote the reflectance of near-infrared (NIR) and red wavelengths.

Reference Wavelengths and Recombined NDVI-Based VIs
The key wavelengths were selected by some physiologically related reference wavelengths with the equation of NDVI that have been used in previous studies (Table 1). Those reference wavelengths were recombined and recalculated as the normalized difference procedure with all samples.

Estimation of PI Occurrence through First-Order Differentiation
The trends of NDVI-based VIs scatter plots were changed along the rice plant development, and the quadratic slope always gradually decreased after reaching the reproductive phase. We assumed that the NDVI might have reached maxima during PI, and therefore, the slope would be zero ( Figure 2).

Estimation of PI Occurrence through First-Order Differentiation
The trends of NDVI-based VIs scatter plots were changed along the rice plant development, and the quadratic slope always gradually decreased after reaching the reproductive phase. We assumed that the NDVI might have reached maxima during PI, and therefore, the slope would be zero ( Figure 2).

Figure 2.
Estimation of PI occurrence. When the quadratic slope is equal to zero, the actual determination of PI was GDDA = 771.3 through dissection, and the estimation of PI was GDDE = 800.

Temperature and Growing Degree-Days (GDD)
Temperature data were obtained from Chishan meteorological station (22°35′ N, 120°36′ E), which is located approximately 7 km away from experimental site. The GDD of "Method 1" as stated by preceding study was selected as a substitution to represent the cultivation timing, according to which [35]: where the Tmax and Tmin are the daily maximum and minimum air temperature, and the Tbase is the base temperature that was set to 10 °C [8].

Statistical Analysis
The statistical analysis of the data was done by using Microsoft Excel 2013 (Microsoft Corporation, Redmond, WA, USA). The quadratic polynomial regression analysis was performed in SigmaPlot 10.0 (Systat Software Inc., San Jose, CA, USA).
The model was validated by leave-one-out (LOO) cross validation. During the LOO process, each group was progressively and alternately held out for model validation, Figure 2. Estimation of PI occurrence. When the quadratic slope is equal to zero, the actual determination of PI was GDD A = 771.3 through dissection, and the estimation of PI was GDD E = 800.

Temperature and Growing Degree-Days (GDD)
Temperature data were obtained from Chishan meteorological station (22 • 35 N, 120 • 36 E), which is located approximately 7 km away from experimental site. The GDD of "Method 1" as stated by preceding study was selected as a substitution to represent the cultivation timing, according to which [35]: where the T max and T min are the daily maximum and minimum air temperature, and the T base is the base temperature that was set to 10 • C [8].

Statistical Analysis
The statistical analysis of the data was done by using Microsoft Excel 2013 (Microsoft Corporation, Redmond, WA, USA). The quadratic polynomial regression analysis was performed in SigmaPlot 10.0 (Systat Software Inc., San Jose, CA, USA).
The model was validated by leave-one-out (LOO) cross validation. During the LOO process, each group was progressively and alternately held out for model validation, while the remaining groups were used for model construction. The appropriateness of the in which the Ai and Pi are the actual and predicted value of the ith data point and n is the number of data points. Estimations were considered excellent if RE is <10%, good between 10% and 20%, fair between 20% and 30%, and poor if it is >30% [34,36].

PI Determination through Manual Dissection
The PI of each group was observed through manual dissection. The actual determination of PI for group 1 was 60 DAT (GDD = 771.3 • C), group 2 was 56 DAT (GDD = 713.9 • C), group 3 was 57 DAT (GDD = 783.4 • C), and group 4 was 55 DAT (GDD = 791.1 • C) ( Table 2). Group 1 encountered low temperature during the tillering stage. As a result, the required DAT of group 1 for PI was five days longer than group 4, which was thoroughly grown under warm conditions. Besides this, group 2 encountered low temperature before the tillering stage. Consequently, the required GDD of group 2 for PI was approximately 80 • C lower than group 4. These biases indicated that the DAT and GDD are both largely affected by weather variability.

NDVI-Based VIs Selection
All samples were regressed with quadratic polynomial model, as we assumed that the trend of this model (when the slope is equal to zero) might be correlated with the entrance of the PI stage. Considering the ability of explanation, NDVI (700, 720), NDVI (700, 730), NDVI (708, 730), NDVI (660, 760), NDVI (700, 760), NDVI (708, 760), and NDVI (708, 800) were selected by the R 2 value that are above 0.7 of quadratic polynomial regression model (Table 3). On the other hand, all the NDVI-based VIs of groups were stratified, where the distribution of groups 1 and 2 were lower than groups 3 and 4 ( Figure 3A-G). The reason could be that group 1 encountered low temperature during the tillering stage, while group 2 encountered low temperature before the tillering stage. A preceding study indicated that more uniquely expressed proteins were found at 20/12 • C (day/night) and frequently alternating stress/non-stress temperature changes, leading the rice plant to complex stress conditions [37]. This also indicated that cold stress at an early stage would have an irreversible impact to rice plants, since the NDVI-based VIs of groups 1 and 2 remained lower than that of groups 3 and 4.

Estimation of PI Occurrence through First-Order Differentiation
The seven NDVI-based VIs were then tested with first-order differentiation in each group to compare the estimated GDD with the actual GDD of PI occurrence (Table 4). Due to the excellent estimation of GDDs for PI occurrence (RE < 10%), NDVI (700, 720), NDVI (700, 730), NDVI (708, 730), NDVI (700, 760), NDVI (708, 760), and NDVI (708, 800) were considered as the suitable NDVI-based VIs for PI determination. Among these combinations, NDVI (708, 760), which was proposed by Tan et al. (2018), coincidentally had the lowest RE value (RE = 5.63%) [18]. This indicated that the entrance of the reproductive phase might have influenced the photosynthetic capacity of the rice plant, such as a lower rate of crop growth and absorption of photosynthetically active radiation (PAR). The other suitable combinations of NDVI-based VIs have not been researched in previous studies.

Estimation of PI Occurrence through First-Order Differentiation
The seven NDVI-based VIs were then tested with first-order differentiation in each group to compare the estimated GDD with the actual GDD of PI occurrence (Table 4). Due to the excellent estimation of GDDs for PI occurrence (RE < 10%), NDVI (700, 720), NDVI (700, 730), NDVI (708, 730), NDVI (700, 760), NDVI (708, 760), and NDVI (708, 800) were considered as the suitable NDVI-based VIs for PI determination. Among these combinations, NDVI (708, 760), which was proposed by Tan et al. (2018), coincidentally had the lowest RE value (RE = 5.63%) [18]. This indicated that the entrance of the reproductive phase might have influenced the photosynthetic capacity of the rice plant, such as a lower rate of crop growth and absorption of photosynthetically active radiation (PAR). The other suitable combinations of NDVI-based VIs have not been researched in previous studies. Table 4. Seven selected NDVI-based VIs were tested by first-order differentiation in each group.

Leave-One-Out (LOO) Cross-Validation
The leave-one-out cross-validation procedure was subsequently performed to test the accuracy of the quadratic polynomial regression of NDVI (700, 720), NDVI (700, 730), NDVI (708, 730), NDVI (700, 760), NDVI (708, 760), and NDVI (708, 800). The LOO regression line of groups 1, 3, and 4 showed the highest plateau, while groups 1, 2, and 4 the lowest ( Figure 4A-F). Due to the great variabilities of the plateau distribution, neither setting up an absolute threshold value nor extended the range for the PI occurrence determination seems to be an ideal method. Although the plateau of regression curve between groups was different, the trend of the slopes was relatively similar. This was in agreement with preceding research [38]; however, the plateau of the curve in this study was different, which was at the PI stage rather than the booting-heading stage. The stability between the trend of the slopes and PI occurrence were tested through first-order differentiation in the leave-one-out cross-validation model. The results showed that NDVI (700, 720) was the most reliable NDVI-based VIs for PI determination since it had the lowest RE value (RE = 4.9%). A previous study represented NDVI application on wheat phenology monitoring and obtained satisfying results that achieved the lowest error of 4.61 days during the jointing stage [39]. Our study performed a slightly better approach with the application of NDVI (700, 720) with a quadratic polynomial regression model, which represented a maximum error of 68.9 GDDs (approximately three days) ( Table 5). The others, NDVI (700, 730), NDVI (708, 730), and NDVI (708, 800), are potential candidates that are competent for PI assessment (Table 5). At this point, the relationship between NDVI-based VIs and PI occurrence could be confirmed through the first-order differentiation of quadratic polynomial regression model. 4.9%). A previous study represented NDVI application on wheat phenology monitoring and obtained satisfying results that achieved the lowest error of 4.61 days during the jointing stage [39]. Our study performed a slightly better approach with the application of NDVI (700, 720) with a quadratic polynomial regression model, which represented a maximum error of 68.9 GDDs (approximately three days) ( Table 5). The others, NDVI (700, 730), NDVI (708, 730), and NDVI (708, 800), are potential candidates that are competent for PI assessment (Table 5). At this point, the relationship between NDVI-based VIs and PI occurrence could be confirmed through the first-order differentiation of quadratic polynomial regression model.

Conclusions
This study revealed the relationship between rice panicle initiation and NDVI-based VIs by the first-order differentiation of quadratic polynomial regression. The results showed that NDVI (700, 720), NDVI (700, 730), NDVI (708, 730), NDVI (700, 760), NDVI (708, 760), and NDVI (708, 800) are potential candidates to determine PI stage. Although the observed values had great variabilities between groups, the trend is stable-when the population of the target reached their plateau of NDVI-based VIs distribution, PI occurred. Due to the synchronistic correlation between PI and plateau of the trend, the PI estimation might be achieved through determining the increasing proportional of saturated values of NDVIbased VIs. Cooperating with the UAV-mounted multispectral or hyperspectral camera and pixel-based analysis of the region of interest, NDVI-based VIs could be a useful indicator to assess rice plant PI stage, thus optimizing the management of fertilization practice.