Monitoring Hybrid Rice Phenology at Initial Heading Stage Based on Low-Altitude Remote Sensing Data

Abstract

Accurate monitoring of hybrid rice phenology (RP) is crucial for breeding rice cultivars and controlling fertilizing amount. The aim of this study is to monitor the exact date of hybrid rice initial heading stage (IHS_DAS) based on low-altitude remote sensing data and analyze the influence factors of RP. In this study, six field experiments were conducted in Ezhou city and Lingshui city from 2016 to 2019, which involved different rice cultivars and nitrogen rates. Three low-altitude remote sensing platforms were used to collect rice canopy reflectance. Firstly, we compared the performance of normalized difference vegetation index (NDVI) and red edge chlorophyll index (CIred edge) for monitoring RP. Secondly, double logistic function (DLF), asymmetric gauss function (AGF), and symmetric gauss function (SGF) were used to fit time-series CIred edge for acquiring phenological curves (PC), the feature: maximum curvature (MC) of PC was extracted to monitor IHS_DAS. Finally, we analyzed the influence of rice cultivars, N rates, and air temperature on RP. The results indicated that CIred edge was more appropriate than NDVI for monitoring RP without saturation problem. Compared with DLF and AGF, SGF could fit CIred edge without over fitting problem. MC of SGF_CIred edge from all three platforms showed good performance in monitoring IHS_DAS with good robustness, R² varied between 0.82 and 0.95, RMSE ranged from 2.31 to 3.81. In addition, the results demonstrated that high air temperature might cause a decrease of IHS_DAS, and the growth process of rice was delayed when more nitrogen fertilizer was applied before IHS_DAS. This study illustrated that low-altitude remote sensing technology could be used for monitoring field-scale hybrid rice IHS_DAS accurately.

1. Introduction

With the rapid growth of the world population, food security has become an urgent issue, attracting worldwide attention. Rice is a very important crop that provides staple food for almost 50% of the world’s population [1]. Hybrid vigor (heterosis) is a universal phenomenon in crops, and the yield of hybrid rice is usually higher than that of conventional rice [2,3]. Hybrid rice makes a great contribution to the world’s food security, with more than half the growth of rice yield [4,5,6]. Hybrid rice is widely planted in large areas of Asian region, including China, Japan, and Southeast Asian countries. Breeding hybrid rice cultivars with high yield and good quality has been an important part of these countries’ strategies for developing modernized agriculture [7].

Rice phenology (RP) is critical for farm management and yield evaluation [8]. RP is divided into germination, seeding, tillering, jointing, booting, heading and flowering, grain filling, and milk stages [9]. Breeding hybrid rice is a time-consuming and laborious physical work. It is difficult for breeders to walk in the rice fields, which are fragmented by roads and water networks [10]. The heading and flowering stage is the essential time for breeders to choose superior parents and execute pollination operation. In general, when hybrid rice is at heading and flowering stage, flowers bloom daily between 9:00 am to 2:00 pm lasting 7–12 days [11,12,13], which means that breeders should finish a series of breeding operations accurately and timely, such as selecting appropriate rice spike, pollination, bagging, and labeling. Moreover, the flowering habit of rice is closely related to rice yield [14]. Thus, monitoring field-scale rice initial heading stage (IHS) accurately is not only important for breeding hybrid rice, but also for rice yield estimation [15].

To date, there are three main methods for monitoring vegetation phenology: ground observation [16], vegetation growth models [17] and remote sensing estimation [18]. Traditionally, ground observation is a time-consuming and laborious work, including sampling, visual observation, taking photos and writing reports, which is suitable to acquire accurate in situ phenology information in small areas [19]. Rice growth models (e.g., CERES-Rice, ORYZA2000, RCM, Beta Model, and SIMRIW) can precisely simulate the growth process of rice [20], but these models consist of many parameters (e.g., rainfall, illumination, air temperature, vegetation cultivars, fertilization and irrigation) that are complicated for monitoring region scale vegetation phenology [21,22,23]. Compared with the aforementioned two methods, many studies have proved that remote sensing is a time-saving and cost-effective method for monitoring vegetation phenology at both local and regional scales, avoiding sampling bias and strenuous labor [24,25,26,27,28,29].

At present, RP is mainly monitored by optical satellites and Synthetic Aperture Radar (SAR) satellites on large scale. Optical satellites have shorter revisiting cycle and lower cost [30]. A variety of vegetation indices (VIs) from optical satellites are widely used to monitor vegetation growth and phenology [31,32]. Sakamoto et al. [33] adopted wavelet and Fourier transforms to filter the time-series enhanced vegetation index (EVI) from moderate resolution imaging spectroradiometer (MODIS) images for monitoring rice planting date, heading date and harvesting date. The root mean square errors (RMSE) of the estimated phenological dates against the statistical data were 12.1 days, 9.0 days and 10.6 days, respectively. Onojeghuo et al. [34] fused MODIS data with Landsat data for acquiring the high spatiotemporal data (30 m spatial resolution and 8-day intervals), and proved that normalized difference vegetation index (NDVI) extracted from the fused data could improve the accuracy of RP. Sun et al. [35] proposed the red edge chlorophyll index (CIred edge) from Sentinel-2 images to evaluate cadmium stress in rice at different stages. However, the optical satellites were affected by atmospheric disturbances, solar radiation and cloud cover [36,37]. The estimation error of phenology could be as high as nine days or even more. SAR satellites have the advantage to work at all-weather situation and all-day, and are very sensitive to vegetation morphological structure information [38,39]. Some studies demonstrated that the backscattering information of SAR had good relationship with RP on large scale. Fikriyah et al. [40] discriminated the transplanted rice and direct seeded rice by temporal Sentinel-1 SAR. TerraSAR-X was applied to monitor RP on Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie (BBCH) scale [41,42]. He et al. [43] used the backscattering characteristics of multi-temporal RADARSAT-2 datasets to monitor the four rice phenological phases (transplanting, vegetative, reproductive and maturity), the overall accuracy (OA) is 86.2%. Although these studies indicated that optical satellites and SAR satellites could monitor large-scale RP successfully, it is still difficult to acquire precise field-scale RP due to the lack of high spatiotemporal satellite data.

Compared with satellites, low-altitude (below 500 m) remote sensing sensors have obvious advantages: less atmospheric interference, short visiting cycle, high spatial resolution and low cost. Low-altitude sensors, such as unmanned aerial vehicles (UAV) [44], fixed towers [45], and hand-held sensors [46], are gradually applied to monitor field-scale vegetation phenology [47]. Yang et al. [48] adopted a deep learning approach to detect RP based on UAV images, with an accuracy rate of 83.9% and mean absolute error (MAE) of 0.18. Bai et al. [49] utilized the fixed towers mounted with RGB cameras to achieve rice heading stage automatic observation. Lin et al. [50] calculated NDVI from a hand-held radiometer to monitor RP, the results showed that the maximum NDVI had a good correlation ship with the heading and flowering stage. Han et al. [51] proposed a deep convolutional neural network (DCNN) approach to analyze hand-held RGB images for real-time RP detection, with a high overall accuracy of 0.913 and low mean absolute error of 0.090. In addition, Zheng et al. [9] compared the ability of two ground-based spectral instruments for monitoring RP, there were 2.3–9.3 days difference between the observed date and estimated date at different growth stages. In a word, these studies proved that low-altitude remote sensing sensors had good performance in monitoring field-scale RP.

To our knowledge, there are three key steps to monitor vegetation phenology based on remote sensing technology: (1) cleaning and flagging remote sensing data; (2) reconstructing vegetation phenological curves; (3) extracting features for monitoring vegetation phenology [52]. However, previous studies mainly focused on studying the phenological characteristics of rice and exploring appropriate methods for monitoring RP. Few studies compared the advantages and disadvantages of different low-altitude remote sensing sensors for field-scale RP detection, and analyzed the influence factors of RP. Thus, the aim of this study is to (1) analyze the change regularity of typical time-series VIs (NDVI and CIred edge) during the whole growth stage of hybrid rice; (2) evaluate the performance of three functions for reconstructing rice phenological curves based on three low-altitude remote sensing data, and extract features of VIs for monitoring IHS_DAS; (3) analyze the influential factors (rice cultivars, N rates, and air temperature) of RP.

2. Materials and Methods

2.1. Experimental Sites

Six hybrid rice (Oryza sativa L. subsp. indica) experiments were conducted in this study (Figure 1). Experiments 1, 2, 5, and 6 were conducted in Ezhou city (30.3752°N, 114.7450°E), Hubei province, central China from May to October in 2016, 2017, 2018 and 2019. Ezhou city is located in a subtropical monsoon climate. The experiment 3 and 4 were conducted in Lingshui city (18.5302°N, 110.0545°E), Hainan province, southern China from December 2017 to May 2018. Lingshui city is located in a tropical marine monsoon climate.

2.2. Experimental Design

Experiments 1, 2, and 3 were conducted in three growing seasons spanning from 2016 to 2018. Rice seeds were sown in seedling beds. 30 days after sowing (DAS), seedlings were transplanted into rice field with a unified density 225,000 plants/hm². The plot area of each rice cultivar was 40 m² (5.0 m × 8.0 m). For all treatments, N fertilizer (180 kg/hm²) was applied in the form of urea: 50% as basal fertilizer before transplanting, 25% at the tillering stage and 25% at the booting stage. Monocalcium phosphate (P₂O₅ 90 kg/hm²) and potassium chloride (K₂O 180 kg/hm²) were applied before transplanting. Other management decisions (pest control and herbicide application) followed the local standard practices.

Experiment 4 was a randomized complete block design with three replications for each treatment, involving two rice cultivars and four N rates. Rice cultivars were Fengliangyou 4 (FLY4) and Luoyou 9348 (LY9348). Four N rates: 0 (N0), 120 (N1), 180 (N2) and 240 (N3) kg/hm² were applied in the form of urea: 50% as basal fertilizer before transplanting, 50% at the booting stage. The 24 individual plot area was 15 m² (5.0 m × 3.0 m), respectively. In addition, other factors were the same with experiments 1, 2, and 3.

Experiment 5 and 6 were conducted from May to October of 2018 and 2019, involving 1014 and 289 rice cultivars, each rice cultivar was planted in an individual plot, and the area of plot was 1 m² (1.0 m × 1.0 m), respectively. In addition, other factors were the same with experiments 1, 2, and 3.

Six experiments involved different hybrid rice cultivars and N rates, as detailed in Figure 2 and

Table 1. Details of six experiments.

Experiment	Year and Study Site	Number of Plots	Plots Area (m2)	Sowing Time (Month/Day/Year)	Number of Rice Cultivars	N Rates (kg/hm2)
1	2016 Ezhou	16	40	5/7/2016	16	180
2	2017 Ezhou	52	40	5/7/2017	52	180
3	2017–2018 Lingshui	40	40	12/10/2017	40	180
4	2017–2018 Lingshui	24	15	12/10/2017	2	0, 120, 180 and 240
5	2018 Ezhou	1014	1	5/25/2018	1014	180
6	2019 Ezhou	289	1	5/15/2019	289	180

.

2.3. Experimental Data

2.3.1. Meteorological Data

Meteorological data were collected from the China Meteorological Data Service Center (http://www.weather.com.cn/), which included maximum air temperature (T_max), minimum air temperature (T_min) and weather conditions spanning from 2016 to 2019. In order to show weather conditions clearly, we used different numbers to represent different weather conditions: 1 represented sunny, 2 represented cloudy, 3 represented overcast, 4 represented shower, 5 represented light rain, 6 represented moderate rain, 7 represented heavy rain, and 8 represented rainstorm, respectively.

2.3.2. Field Phenological Data

In our study, the phenological time was recorded as days after sowing (DAS). Ground field observations and records were carried out according to the specifications for agrometeorological observation-Rice (http://www.cma.gov.cn/root7/auto13139/201903/t20190329_519169.html) from the China Meteorological Administration. For example, when 10% of investigated rice grew spikes, the exact date was defined as the time of rice phenology at initial heading stage (IHS_DAS). Then, we investigated 40 hills of rice plants in experiments 1, 2, 3, and 4, and 20 hills of rice plants in experiments 5 and 6 were investigated by random method, respectively.

2.3.3. Daily Reflectance

Three sets of SKYE radiometers (SKR 1860, SKYE Instruments Ltd., Llandrindod Wells, UK) were employed to collect four bands reflectance every 30 min. In order to distinguish SKYE radiometers, we marked them as A_SKYE, B_SKYE and C_SKYE. Each set of SKYE radiometers consisted of two sensors: the upward radiometer sensor (URS) and the downward radiometer sensor (DRS). URS measured downwelling (DW) radiance with the aid of a cosine corrector, which collected radiance from sky in 180° hemisphere above the sensor. DRS measured the 25° field of view upwelling (UW) radiance from rice canopy. During the whole growth stage of rice, all SKYE radiometers were positioned above rice canopy at a height of 2 m with a 25° field-of-view of 0.62 m² area (Figure 3). The specifications of SKYE radiometers are listed in

Table 2. The specifications of SKYE radiometers.

SKYE Radiometers	URS		DRS
	Central Wavelength (nm)	Band Width (nm)	Central Wavelength (nm)	Band Width (nm)
A_SKYE	550.3	39.9	549.7	39.8
	655.1	34.6	655.3	34.7
	717.2	23.7	716.8	23.6
	865.6	28.5	865.0	31.5
B_SKYE	550.0	39.9	549.9	38.7
	655.1	36.7	654.9	36.0
	716.6	23.6	717.7	23.9
	865.3	27.9	865.7	29.7
C_SKYE	551.8	39.6	549.5	39.7
	652.9	35.7	653.4	35.9
	717.0	26.2	715.8	25.2
	857.4	38.0	858.3	37.8

.

SKYE radiometric calibration: the relative reflectance of four bands could be calculated using the following equation.

$$\begin{gathered} R_{\lambda }=\frac{U{W}_{\lambda }}{D{W}_{\lambda }}\times \mathsf{\pi } \\ \left(\lambda =550, 655, 717 \mathrm{and} 865 \mathrm{nm}\right) \end{gathered}$$

(1)

where $U{W}_{\lambda }$ and $D{W}_{\lambda }$ are calculated (mol m⁻² s⁻¹), respectively.

In four growing seasons spanning from 2016 to 2019, SKYE radiometers collected continuously reflectance of three rice cultivars: LY 9348, FLY 4, and Yi Liang You Hua Zhan (YLYHZ), as detailed in

Table 3. Details of SKYE measurements.

Experiment	Year and Study Site	SKYE Radiometers	Rice Cultivars	IHSDAS	N Rates (kg/hm2)
1	2016 Ezhou	C_SKYE	LY 9348	97	180
2	2017 Ezhou	A_SKYE	FLY 4	85	180
		B_SKYE	YLYHZ	77	180
		C_SKYE	LY 9348	87	180
4	2017–2018 Lingshui	A_SKYE	LY 9348	106	120
		B_SKYE	FLY 4	109	120
		C_SKYE	LY 9348	109	240
6	2019 Ezhou	A_SKYE	FLY 4	89	180
		B_SKYE	LY 9348	91	180

. The median reflectance between 10:00 am to 2:00 pm was calculated as daily reflectance to reduce the impact of weather drastic change.

2.3.4. The Hyperspectral Reflectance

The hyper-spectrometer ASD Field Spec 4 (ASD Inc., Boulder, CO, USA) was employed to collect rice canopy reflectance in experiment 1, 2, 3 and 4, with the spectrum ranging from 350–2500 nm. The measurements were conducted on sunny days from 10:00 am to 2:00 pm. A white reference panel (Chemical composition: BaSO₄, Size: 25.4 × 25.4 cm) was used to make relative radiometric correction before each measurement. The optical fiber probe of ASD was positioned above rice canopy at a height of 1 m with a 25° field-of-view of 0.15 m² area. Five fixed points’ spectral measurements were averaged to reduce uncertainty error of spectral measurement in each plot. The hyperspectral reflectance involved different rice cultivars, as detailed in

Table 4. Details of ASD measurements.

Experiment	Year and Study Site	Measurement Dates (Month/Day/Year)	Days after Sowing (DAS)	Number of Rice Cultivars
1	2016 Ezhou	7/28/2016	82	16
		8/11/2016	96	16
		9/2/2016	118	16
		9/21/2016	137	16
2	2017 Ezhou	6/20/2017	44	52
		7/7/2017	61	52
		7/18/2017	72	52
		8/6/2017	91	52
		8/23/2017	108	52
		9/12/2017	128	52
3	2017–2018 Lingshui	2/4/2018	56	40
		2/25/2018	77	40
		3/9/2018	89	40
		3/19/2018	99	40
		3/31/2018	111	40
		4/28/2018	138	40
4	2017–2018 Lingshui	2/4/2018	56	24
		2/22/2018	74	24
		3/4/2018	84	24
		3/14/2018	94	24
		3/24/2018	104	24
		4/4/2018	115	24
		4/29/2018	140	24

.

2.3.5. The UAV Multispectral Reflectance

In this study, a multispectral camera (MCA, Tetracam Inc., Chatsworth, CA, USA) was mounted on an eight-propeller UAV (S1000, DJI Technology Co., Ltd., Shenzhen, China) to acquire the panoramic multispectral images in experiment 5 and 6. The UAV flight was conducted on sunny days between 10:00 am to 2:00 pm. The flight duration varied between 10 and 15 min. The flight height was 120 m with a field of view 83.25 m × 66.50 m, with the ground sampling distance 6.5 cm/pixel.

MCA consisted of twelve individual digital miniature cameras equipped with different band pass filters centered at the wavelength of 490 nm, 520 nm, 550 nm, 570 nm, 670 nm, 680 nm, 700 nm, 720 nm, 800 nm, 850 nm, 900 nm, and 950 nm respectively (

Table 5. Center wavelength and bandwidth of MCA.

Band Number	Center Wavelength (nm)	Bandwidth (nm)	Band Number	Center Wavelength (nm)	Bandwidth (nm)
1	490	10	7	700	10
2	520	10	8	720	10
3	550	10	9	800	10
4	570	10	10	850	10
5	670	10	11	900	20
6	680	10	12	950	40

). The multispectral images are recorded as digital number (DN) values. Prior to the flight, six calibration blankets were placed on the flat road as radiometric calibration targets (Figure 4). One panoramic multispectral image could simultaneously contain the six calibration blankets and the total field plots within the field of view of MCA. The calibration blankets were made with highly durable woven polyester fabric at the size of 1.2 m × 1.2 m, having the relatively constant reflectance of 0.03, 0.12, 0.24, 0.36, 0.56, and 0.80, respectively. As a linear relationship was assumed between surface reflectance and at-sensor radiance, the empirical line method (ELM) was used to convert DN to surface reflectance [53,54]. $Gai{n}_{\mathrm{k}}$ and $offse{t}_{\mathrm{k}}$ were calculated using the least-square method from the relatively constant reflectance (0.03, 0.12, 0.24, 0.36, 0.56, and 0.80) and the corresponding DN values of calibration targets. Then $Gai{n}_{\mathrm{k}}$ and $offse{t}_{\mathrm{k}}$ were applied to calculate the canopy surface of each pixel in multispectral images. The specific equations are as follows.

$$\left(\begin{array}{c}0.03\\ 0.12\\ 0.24\\ 0.36\\ 0.56\\ 0.80\end{array}\right)=\left(\begin{array}{c}D{N}_{\left(1,\mathrm{k}\right)}\\ D{N}_{\left(2,\mathrm{k}\right)}\\ D{N}_{\left(3,\mathrm{k}\right)}\\ D{N}_{\left(4,\mathrm{k}\right)}\\ D{N}_{\left(5,\mathrm{k}\right)}\\ D{N}_{\left(6,\mathrm{k}\right)}\end{array}\right)\times Gai{n}_{\mathrm{k}}+offse{t}_{\mathrm{k}}$$

(2)

$$\begin{gathered} R_{\mathrm{k}}=D{N}_{\mathrm{k}}\times Gai{n}_{\mathrm{k}}+offse{t}_k \\ \left(k=490, 520, 550, 570, 670, 680, 700, 720, 800, 850, 900 \mathrm{and} 950 \mathrm{nm}\right) \end{gathered}$$

(3)

where k is the center wavelengths of 12 bands, DN(i,k) are the digital number values of the calibration target i (i = 1, 2, 3, 4, 5, 6, respectively) in band k. $Gai{n}_{\mathrm{k}}$ and $offse{t}_{\mathrm{k}}$ represent the slope and intercept. $R_{\mathrm{k}}$ and $D{N}_{\mathrm{k}}$ are the canopy surface reflectance and the digital number values of multispectral images.

The UAV system collected multispectral images from May to October of 2018 and 2019, involving 1303 rice cultivars in total (

Table 6. Details of MCA measurements.

Experiment	Year and Study Site	Measurement Dates (Month/Day/Year)	Days after Sowing (DAS)	Number of Rice Cultivars
5	2018 Ezhou	6/27/2018	33	1014
		7/6/2018	42	1014
		7/11/2018	47	1014
		7/16/2018	52	1014
		7/27/2018	63	1014
		8/2/2018	69	1014
		8/9/2018	76	1014
		8/15/2018	82	1014
		8/21/2018	88	1014
		8/27/2018	94	1014
6	2019 Ezhou	2/7/2019	42	289
		6/7/2019	48	289
		14/7/2019	52	289
		22/7/2019	60	289
		27/7/2019	68	289
		1/8/2019	73	289
		6/8/2019	78	289
		11/8/2019	83	289
		16/8/2019	88	289
		22/8/2019	93	289
		29/8/2019	99	289
		3/9/2019	106	289
		9/9/2019	111	289
		17/9/2019	117	289

). PixelWrench2 software (Tetracam Inc., Chatsworth, CA, USA) was used to perform images preprocessing, including noise correction, vignetting correction, lens distortion correction, band-to-band registration and band stacking [55,56]. The radiometric calibration and reflectance extraction processing were performed using the ENVI software (Exelis Visual Information Solutions, Boulder, CO, USA).

After image preprocessing, we selected a region of interest (ROI) in the multispectral image, which included 100 pixels of 0.42 m² area in each plot. The plot-level reflectance was calculated as the average reflectance of all pixels within ROI.

2.4. Research Methods

2.4.1. Vegetation Indices (VIs)

Vegetation indices (VIs) refer to the spectral reflectance transformations of two or more bands, which are related to vegetation photosynthesis and nutrition status [57]. NDVI is the most commonly used VIs for quantitative remote sensing [58]. CIred edge has good performance in the estimation of chlorophyll content, leaf area index and biomass [59]. In this study, NDVI and CIred edge are employed to monitor RP. They are calculated as:

$$NDVI=({\rho }_{NIR}-{\rho }_{Red})/({\rho }_{NIR}+{\rho }_{Red})$$

(4)

$$CIred edge=({\rho }_{NIR}/{\rho }_{Red edge})-1$$

(5)

where ${\rho }_{NIR}$ is the reflectance in the wavelength of 865 nm in SKYE (800 nm in ASD and MCA), ${\rho }_{Red edge}$ is the reflectance in the wavelength of 717 nm in SKYE (720 nm in ASD and MCA), ${\rho }_{Red}$ is the reflectance in the wavelength of 655 nm in SKYE (680 nm in ASD and MCA).

2.4.2. Fitting Functions

There were several methods to smooth and fit time-series VIs curves, such as dynamic threshold approach [60], derivative method [61], savitzky-golay smooth [62], double logistic function (DLF) [25], asymmetric gauss function (AGF) [63] and symmetric gauss function (SGF) [64]. Previous studies reported that DLF and AGF had good performance in fitting time-series VIs for monitoring vegetation phenology [65,66,67,68,69]. In our study, to obtain high quality time-series VIs, three fitting functions (DLF, AGF, and SGF) were applied to fit time-series VIs for reconstructing three phenological curves (PC): DLF_VIs, AGF_VIs, and SGF_VIs. The feature: the maximum curvature (MC) of PC was extracted to monitor IHS_DAS.

The calculation formulas are as follows:

$$\mathrm{DLF}\_\mathrm{VIs}=\frac{c}{\left(1+e^{\left(a+b\ast DAS\right)}\right)}+\frac{g}{\left(1+e^{\left(d+f\ast DAS\right)}\right)}$$

(6)

$$\mathrm{AGF}\_\mathrm{VIs}=h{e}^{-{\left(\frac{DAS-i}{j}\right)}^2}+k{e}^{-{\left(\frac{DAS-l}{m}\right)}^2}$$

(7)

$$\mathrm{SGF}\_\mathrm{VIs}=p{e}^{-{\left(\frac{DAS-q}{r}\right)}^2}$$

(8)

where $DAS$ are the days after rice sowing, corresponding the time of VIs acquisition, $e$ is the natural constant, a, b, c, d, f, g, h, i, j, k, l, m, p, q, and r are the best parameters of three fitting functions (DLF, AGF, and SGF), respectively.

$$\mathrm{MC}=\max (\frac{\left|P{C}^{″}\right|}{{(1+{(P{C}^{\prime })}^2)}^{\frac{3}{2}}})$$

(9)

where PC are the phenological curves: DLF_VIs, AGF_VIs, and SGF_VIs, PC′ is the first derivative of PC, and PC″ is the second derivative of PC. MC is the maximum curvature of PC (Figure 5).

2.4.3. The Process of Modeling

Figure 6 shows the flowchart. Firstly, we calculated the reflectance of different remote sensing platforms. The statistical analysis of meteorological data and field phenological data were carried out. Secondly, daily NDVI and CIred edge were computed using the SKYE reflectance. We analyzed the phenological change of rice during the whole growth stage, and compared the characteristics of NDVI and CIred edge for monitoring RP. Thirdly, the VIs (NDVI or CIred edge) were computed using the reflectance of ASD and MCA. Three fitting functions (DLF, AGF and SGF) were applied to fit the suitable VIs for reconstructing PC: DLF_VIs, AGF_Vis, and SGF_VIs. We evaluated the performance of three fitting functions. Finally, we extracted the MC of PC for monitoring IHS_DAS and evaluated the accuracy.

2.4.4. The Assessment of Models

The coefficient of determination (R²) and root mean square error (RMSE) were used to evaluate the performances of three fitting functions and the accuracy of monitoring IHS_DAS. Mathematically, a higher R² and smaller RMSE represent better model accuracy.

$${\mathrm{R}}^2=1-\frac{{\sum }_{i=1}^n{\left(y_i-x_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(10)

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^n{\left(y_i-x_i\right)}^2}{n}}$$

(11)

where $x_i$ and $y_i$ are the estimated and measured values, $\overline{y}$ is the average of the measured values, and n is the sample number, respectively.

The mathematical models were performed using MATLAB R2018a (MathWorks, Inc., Natick, MA, USA). Graphics were prepared using Origin 9.1 software program (OriginLab Corporation, Northampton, MA, USA).

3. Results

3.1. Statistical Analysis of Meteorological Data

Figure 7 shows that the weather conditions have obvious difference between two climatic zones, the climate of Lingshui city is more stable than that of Ezhou city.

In Ezhou city, the weather conditions in experiment 1 and 2 fluctuate greatly. It often rains during the whole growth period of rice (Figure 7a,c). Sometimes, there is even a regional heavy rain or rainstorm. Further, the continuous high temperature (T_max ≥ 35 °C) appears between DAS 76 to 88 and DAS 64 to 82 (Figure 7b,d). The weather conditions in experiment 5 and 6 are mainly sunny to cloudy with light rain occasionally (Figure 7g,i), however, the continuous high temperature (T_max ≥ 35 °C) appears between DAS 49 to 79 and DAS 70 to 105 (Figure 7h,j).

In Lingshui city, the weather conditions in experiments 3 and 4 are very stable. It is mainly sunny to cloudy with small air temperature change (Figure 7,f), there is not continuous high air temperature appearing during the whole growth period of rice. Good weather conditions in Lingshui city are helpful to obtain accurate remote sensing data.

3.2. Statistical Analysis of IHS_DAS

We carried out statistical analysis of IHS_DAS in six experiments (

Table 7. The statistical characteristics of IHS_DAS.

Experiment	Number of Plots	Minimum	Mean	Maximum	Standard Deviation	Variance	Coefficient of Variation (%)
1	16	92	98.31	104	3.45	11.96	3.51
2	52	76	84.00	94	5.59	31.29	6.66
3	40	96	110.60	119	3.71	13.73	3.35
4	24	103	106.58	110	2.15	4.60	2.01
5	1014	59	76.38	103	7.44	55.40	9.74
6	289	60	87.42	122	9.27	85.84	10.60

). The mean values of IHS_DAS are 98.31, 84.00, 110.60, 76.38, and 87.42 in experiment 1, 2, 3, 5 and 6, while the coefficients of variation are 3.51%, 6.66%, 3.35%, 9.74%, and 10.60%. IHS_DAS ranges from 59 to 122 in six experiments. The statistical indicators of IHS_DAS show the significant differences in IHS_DAS among different rice cultivars. The maximum and minimum value of IHS_DAS in experiment 4 are 103 and 110 involving four N rates, it implies that N rates have obvious effects on IHS_DAS. Therefore, we take into account rice cultivars and nitrogen treatments in this study and focus on building robust and accurate models for monitoring IHS_DAS.

3.3. Comparative Analysis of Daily NDVI and CIred Edge

The SKYE radiometers were used to acquire daily reflectance data of 9 plots in four rice growing seasons. Figure 8 shows that the change regularity of NDVI is different from that of CIred edge. NDVI increases rapidly and achieves the maximum value quickly. Before DAS 70, NDVI increases rapidly; between DAS 80 to 100, NDVI remains stable and is close to 0.9. When rice is at maturity stage (after DAS 100), NDVI declines slowly with the senescence of leaves and stems, and it is close to 0.5 before harvesting. Especially, as Figure 8b–d present, NDVI achieves the maximum value before DAS 55 corresponding to the middle stage of tillering, which indicates that rice grows rapidly after transplanting in field. There is an obvious saturation problem occurring before IHS_DAS (the green lines in Figure 8). Thus, it is difficult to extract features from NDVI for monitoring IHS_DAS.

The curves of CIred edge are different in four rice growing seasons, but the change regularity of CIred edge is similar in the same growing season. CIred edge gradually increases with an exponential growth at tillering and jointing stage. CIred edge achieves the maximum value before IHS_DAS. The maximum value of CIred edge does not exceed 4 without saturation problem. CIred edge remains stable and is close to the maximum value before IHS_DAS. As shown in Figure 8, there are obvious inflection points of CIred edge appearing at IHS_DAS. It suggests that we can take advantage of the feature: the MC of time-series CIred edge for monitoring IHS_DAS. After IHS_DAS, CIred edge decreases linearly with the senescence of leaves and stems.

In addition, there is some noise in time-series CIred edge. For example, Figure 8b–d display abnormal fluctuations: two peaks and one trough appear between DAS 60 to 95. Thus, it is desirable to use some smoothing and fitting methods to eliminate noise of time-series CIred edge.

The comparative analysis indicates that CIred edge presents no saturation problem during the whole growth period of rice. CIred edge is better than NDVI for field-scale RP detection. Meanwhile, CIred edge has the obvious feature (MC) for monitoring IHS_DAS.

3.4. Fitting Several Source CIred Edge

3.4.1. Fitting CIred Edge of SKYE

Three fitting functions (DLF, AGF and SGF) are applied to fit CIred edge of SKYE in four growing seasons. Figure 9 illustrates the relationship between the original CIred edge and the fitted phenological curves. The results show that DLF, AGF and SGF can eliminate noise of the original CIred edge and help to present the change regularity of time-series CIred edge accurately. However, the performance of fitting models varies at different growth stages. DLF has good performance in fitting CIred edge at the early stages (before DAS 80). AGF is suitable to fit CIred edge between DAS 80 to 100, and SGF shows the highest accuracy after rice begins to grow into maturity (after DAS 100).

Table 8. The performance of the models in fitting CIred edge of SKYE.

Models	Fitting Function	R2	RMSE	Fitting Function	R2	RMSE	Fitting Function	R2	RMSE
a	DLF	0.97	0.07	AGF	0.97	0.07	SGF	0.92	0.13
b		0.94	0.19		0.96	0.14		0.88	0.27
c		0.96	0.13		0.96	0.14		0.94	0.16
d		0.94	0.17		0.95	0.16		0.91	0.22
e		0.98	0.13		0.98	0.14		0.95	0.23
f		0.98	0.14		0.98	0.13		0.92	0.36
g		0.98	0.14		0.98	0.15		0.95	0.28
h		0.97	0.14		0.96	0.15		0.89	0.26
i		0.96	0.15		0.97	0.14		0.89	0.25

summarizes the statistics (R² and RMSE) of three fitting models. For DLF models, R² varies between 0.94 and 0.98, RMSE ranges from 0.07 to 0.19. The model b shows the worst performance among DLF models, whereas the model e displays the best performance among DLF models. For AGF models, R² varies between 0.95 and 0.98, and RMSE ranges from 0.07 to 0.16. The model d shows the worst performance among AGF models, whereas the model f displays the best performance among AGF models. For SGF models, R² varies between 0.88 and 0.95, RMSE ranges from 0.13 to 0.36. The model b shows the worst performance among AGF models, whereas the model e displays the best performance among SGF models.

Overall, compared with AGF and SGF, DLF proves to have good performance in fitting CIred edge of SKYE during the whole growth period of rice.

3.4.2. Fitting CIred Edge of ASD

The performances of fitting models in four experiments are shown in

Table 9. The performance of the models in fitting CIred edge of ASD.

Experiment	Number of Plots	Fitting Functions	R2			RMSE
			Minimum	Mean	Maximum	Minimum	Mean	Maximum
1	n = 16	DLF	0.95	0.99	1	0	0.02	0.09
		AGF	1	1	1	0	0	0
		SGF	0.91	0.97	0.99	0.004	0.06	0.11
2	n = 52	DLF	0.84	0.98	1	0	0.05	0.18
		AGF	0.25	0.90	0.99	0.005	0.12	0.44
		SGF	0.62	0.89	0.98	0.08	0.16	0.29
3	n = 40	DLF	0.98	0.99	1	0	0.02	0.07
		AGF	0.92	0.98	0.99	0.01	0.08	0.19
		SGF	0.92	0.97	0.99	0.03	0.09	0.20
4	n = 24	DLF	0.97	0.99	1	0	0.03	0.12
		AGF	0.96	0.98	0.99	0.008	0.10	0.16
		SGF	0.95	0.97	0.99	0.07	0.12	0.17

. In experiment 1, R² and RMSE vary from 0.91–1 and 0–0.11. In experiment 2, R² and RMSE vary from 0.25–1 and 0–0.44. In experiment 3, R² and RMSE vary from 0.92–1 and 0–0.20. In experiment 4, R² and RMSE vary from 0.95–1 and 0–0.17.

We acquire 4, 6, 6, and 7 batches of CIred edge in experiments 1, 2, 3, and 4 (Table 4). There are six, six, and three parameters that need to be optimized in DLF, AGF, and SGF models. When the frequency of times for acquiring CIred edge is low, the over fitting problem of DLF and AGF models appears (R² = 1, RMSE = 0). Compared with DLF and AGF models, the accuracy of SGF models is lower in four experiments. However, there are no obvious over fitting problem appearing in SGF models. The results suggest that the fitting functions with fewer parameters, such as SGF, maybe suitable to avoid over fitting problem. SGF can reconstruct vegetation phenological curves using low-frequency CIred edge.

3.4.3. Fitting CIred Edge of MCA

Table 10. The performance of the models in fitting CIred edge of MCA.

Experiment	Number of Plots	Fitting Functions	R2			RMSE
			Minimum	Mean	Maximum	Minimum	Mean	Maximum
5	n = 1014	DLF	0.18	0.95	0.99	0.015	0.13	0.86
		AGF	0.32	0.92	0.99	0.02	0.19	0.90
		SGF	0.32	0.90	0.99	0.04	0.23	0.96
6	n = 289	DLF	0.69	0.96	0.99	0.05	0.16	0.34
		AGF	0.01	0.94	0.99	0.08	0.20	1.27
		SGF	0.24	0.91	0.98	0.15	0.27	0.47

summarizes the performance of DLF, AGF and SGF. For DLF models, R² varies between 0.18 and 0.99, RMSE ranges from 0.015 to 0.86. For AGF models, R² varies between 0.01 and 0.99, RMSE ranges from 0.02 to 1.27. For SGF models, R² varies between 0.24 and 0.99, RMSE ranges from 0.04 to 0.96. In addition, the performance evaluation implies that the fitting phenological curves (DLF_CIred edge, AGF_CIred edge and SGF_CIred edge) not only have a good correlation ship with the original CIred edge, but also can avoid over fitting problem (R² = 1, RMSE = 0).

In summary, compared with AGF and SGF, DLF is more suitable to fit time-series CIred edge of MCA, with higher R² and lower RMSE. This result is consistent with the result of SKYE data.

3.5. Monitoring IHS_DAS Based on Several Source CIred Edge

3.5.1. Monitoring IHS_DAS Based on CIred Edge of SKYE

We extracted the MC of DLF_CIred edge, AGF_CIred edge and SGF_CIred edge for monitoring IHS_DAS. Figure 10 shows the relationship between MC and IHS_DAS. MC of DLF_CIred edge has relatively lower correlation with IHS_DAS, with the values of R² and RMSE being 0.85 and 4.16. Compared with MC of DLF_CIred edge, MC of AGF_CIred edge and MC of SGF_CIred edge have good correlation with IHS_DAS, with R² values no less than 0.90, with RMSE values no more than 3.00. The results indicate that MC of CIred edge can be used to monitor IHS_DAS effectively.

3.5.2. Monitoring IHS_DAS Based on CIred Edge of ASD

Figure 11 shows the relationship between MC and IHS_DAS in experiment 1, 2, 3 and 4. MC of DLF_CIred edge and MC of AGF_CIred edge have the lowest correlation with IHS_DAS, with R² values of 0.46 and 0.46, RMSE values of 9.18 and 9.23. MC of SGF_CIred edge proves to have good fitting characteristics for monitoring IHS_DAS, with the values of R² and RMSE being 0.95 and 2.86.

3.5.3. Monitoring IHS_DAS Based on CIred Edge of MCA

Figure 12 shows the relationship between MC and IHS_DAS in experiment 5 and 6. MC of DLF_CIred edge and MC of AGF_CIred edge have the lowest correlation with IHS_DAS, with R² values of 0.43 and 0.45, RMSE values of 6.90 and 6.76. MC of SGF_CIred edge has good correlation with IHS_DAS, with the values of R² and RMSE being 0.82 and 3.81. The results are consistent with the results based on CIred edge of ASD.

3.6. Effects of Rice Cultivars and N Rates on IHS_DAS

In order to understand the influence of rice cultivars and N rates on IHS_DAS, experiment 4 was conducted in a randomized complete block designed with three replications for each treatment, involving two rice cultivars (FLY4 and LY9348) and four N rates (N0, N1, N2 and N3). The analysis of variance (ANOVA) indicates that rice cultivars and N rates have obvious effects on IHS_DAS (

Table 11. Analysis of variance (ANOVA) of IHS_DAS.

Experiment	Variable	Factor	Degree of Freedom	Sum of Squares	Mean Square	F Value
4	IHSDAS	Rice Cultivars	1	1.50	1.50	5.14 *
	IHSDAS	N rates	3	95.17	31.72	108.76 **
	IHSDAS	Rice Cultivars × N rates	3	4.50	1.50	5.14 *

Note: * denotes the significant level of α = 0.05, ** denotes the extremely significant level of α = 0.01.

). There is a significant difference between LY9348 and FLY4. The difference of IHS_DAS reaches an extremely significant level under different N rates. The interaction of different rice cultivars and N rates on IHS_DAS reaches a significant level. The ANOVA suggests that we would better take into account rice cultivars and N rates when monitoring IHS_DAS.

Figure 13 shows that IHS_DAS of FLY4 are larger than that of LY9348 under same N rates. Under N1, the average IHS_DAS of FLY4 and LY9348 are 104.66 and 103.66. While under N3, IHS_DAS of FLY 4 and LY 9348 are 110 and 109. With the increase of nitrogen fertilizer application, IHS_DAS of LY9348 and FLY4 increase correspondingly, the results are consistent with previous study [70]. In general, when using more nitrogen fertilizer before heading and flowering stage, the rice growth process might be delayed, leaves would remain green, and spikes appeared later, resulting in a longer life cycle of rice.

4. Discussion

4.1. Daily CIred Edge for Monitoring RP

Previous studies showed that time-series NDVI was suitable for monitoring vegetation phenology [71,72]. However, the major problem of NDVI was that it became saturated for vegetation with moderate to high density [73]. Our study confirmed this fact: daily NDVI of SKYE became significantly saturated before flowers and spikes appeared in rice. Thus, it was difficult to extract features from NDVI for monitoring IHS_DAS. In comparison with NDVI, daily CIred edge of SKYE had no saturation problem during the whole rice growth stage, and was more suitable for monitoring IHS_DAS.

Figure 8 showed that the change regularity of daily CIred edge was closely related to the actual growth process of rice. To our knowledge, CIred edge was originally proposed for high accuracy chlorophyll estimation [74], some studies have demonstrated that the canopy chlorophyll content was closely related to vegetation growth period [75,76]. Therefore CIred edge was better than NDVI for monitoring RP, and the result was consistent with the previous study [9].

In addition, the maximum air temperature (Tmax) has a significant effect on daily CIred edge. Figure 14 displays daily ${\rho }_{NIR}$ and ${\rho }_{Red edge}$ of C_SKYE and Tmax in Ezhou city and Lingshui city. In Lingshui city, Tmax is stable and always below 35 °C. However, the continuous high air temperature (Tmax ≥ 35 °C) often appears with drastic change in Ezhou city. The abnormal change of ${\rho }_{NIR}$ often occurs between DAS 60 to 95 when Tmax has drastic change (Figure 14a). Meanwhile, the abnormal CIred edge also appeared between DAS 60 to 95 (Figure 8d). In general, rice canopy reflectance presents a spectral curve of typical green vegetation before IHS_DAS, ${\rho }_{Red edge}$ is much lower than ${\rho }_{NIR}$ during different growth stage [77,78]. In our study, compared with ${\rho }_{NIR}$, ${\rho }_{Red edge}$ slowly increases with the progress of rice growth and is not sensitive to high air temperature. However, the drastic change of Tmax could cause the sudden change of ${\rho }_{NIR}$, resulting in the abnormal fluctuation of CIred edge appearing.

Rice is a crop that is sensitive to high temperature [79] and water [80]. The continuous high air temperature or additional rain had great influence on the growth process of rice, IHS_DAS of rice may be advanced or postponed [81]. Moreover, strong sunlight can burn rice leaves and change the rice canopy structure. Thus, it is better to take measures (smooth and fitting methods) to eliminate noise before time-series CIred edge was applied to monitor RP. VIs were indeed an effective tool for monitoring vegetation phenology [82,83]. However, if we can make better use of other information, for example, texture, three-dimensional information and coverage, the accuracy of monitoring RP will be improved in the future [84].

4.2. Comparative Analysis of DLF, AGF and SGF for Fitting Several Source CIred Edge

In order to obtain high quality time-series CIred edge, we used DLF, AGF, and SGF to fit CIred edge. The results indicated that DLF and AGF proved to have good performance in fitting CIred edge of SKYE and MCA with higher R² and lower RMSE. However, when the frequency of times for acquiring CIred edge is low, Table 9 shows that the over fitting problem (R² = 1, RMSE = 0) of DLF and AGF models appears. Compared with DLF and AGF models, SGF models have a lower correlation ship with the original CIred edge of ASD, but without obvious over fitting problem.

There are six, six, and three parameters that need to be optimized in DLF, AGF, and SGF models. DLF and AGF are relatively complex, more parameters mean that more CIred edge data should be put into DLF and AGF models. Low-frequency measured data was more likely to cause over fitting problem when more parameters need to be optimized in fitting functions. As we know, the over fitting problem may misinterpret important phenological features and directly affect the accuracy of monitoring vegetation phenology [52]. Thus, DLF and AGF are not suitable to fit low-frequency CIred edge. Acquiring high-frequency CIred edge will help to highlight the phenological features for monitoring RP.

Overall, DLF and AGF could be effectively applied to fit dense CIred edge with higher R² and lower RMSE, and to describe more detailed changes of rice canopy at different growth stages. SGF had good performance in fitting low-frequency CIred edge.

4.3. Monitoring IHS_DAS Based on Several Source CIred Edge

Time-series CIred edge from SKYE, ASD and MCA were reconstructed by DLF, AGF and SGF. We extracted the features: MC of DLF_CIred edge, AGF_CIred edge and SGF_CIred edge for monitoring IHS_DAS.

MC of DLF_CIred edge and AGF_CIred edge from SKYE had good correlation with IHS_DAS. However, MC of DLF_CIred edge and AGF_CIred edge from ASD and MCA did not perform well for monitoring IHS_DAS. MC of SGF_CIred edge from three platforms all have good performance in monitoring IHS_DAS, R² varies between 0.82 and 0.95, RMSE ranges from 2.31 to 3.81, which means that MC of SGF_CIred edge could be applied to monitor IHS_DAS with good robustness. The rice flowers and spikes were usually bright yellow with contrasting colour to the rice leaves [49]. The chlorophyll content of flowers and spikes were also lower than that of green leaves. After flowers and spikes came out, the canopy chlorophyll of rice would decrease dramatically, resulting in the obvious inflection points of CIred edge appearing at IHS_DAS. Thereby, the MC of time-series CIred edge was an essential feature for monitoring field-scale IHS_DAS.

In this study, the accuracy of monitoring IHS_DAS was different using MC extracted from three low-altitude remote sensing data. Three different platforms have the advantages and disadvantages. SKYE radiometers could continuously collect in situ field-scale rice canopy reflectance at very high temporal resolution, which helped analyze the change regularity of phenological curves and extract features for monitoring RP. However, they were mounted in the fixed positions for monitoring several fields, lacking mobility and flexibility [85]. ASD had wider spectrum range and high radiation resolution for extracting accurate VIs, but it was time-consuming to acquire the hyperspectral reflectance of different rice cultivars. The multispectral UAV system had the advantages of mobility, it could collect multispectral images for monitoring real-time large-scale vegetation growth [56,86,87,88]. However, the processing flow of multispectral image was relatively complicated, including radiometric calibration, geometric correction, and image mosaic. It is necessary to acquire the high spatiotemporal and hyperspectral data for improving the accuracy of monitoring field-scale IHS_DAS.

4.4. The Influence Factors on IHS_DAS

There were many factors influencing vegetation phenology, such as illumination, air temperature, vegetation varieties and fertilizer application [89,90]. In our study, we mainly took into account rice cultivars, N rates and air temperature for monitoring IHS_DAS. Table 11 indicates that rice cultivars and N rates have obvious effects on IHS_DAS. The interaction of different rice cultivars and N rates on the IHS_DAS reaches significant level. With the increase of N rates, IHS_DAS of LY 9348 and FLY 4 both increase, which indicates that the growth process of rice is delayed when more nitrogen fertilizer is applied before heading stage, and leaves would remain green, flowers and spikes appear later [70]. In a word, rice cultivars and nitrogen fertilizer application have great influence on RP [91,92].

Air temperature could not only affect the near infrared reflectance, but also disturb RP. Compared with Lingshui city, the weather of Ezhou city fluctuates greatly, it often remains continuous high temperature (T_max ≥ 35 °C). Table 3 shows that IHS_DAS of LY 9348 and FLY 4 planted in Ezhou city are lower than that in Lingshui city. The results confirm that high air temperature might accelerate the whole growth process of rice and shorten growth duration [93,94,95]. The IHS_DAS of hybrid rice might decrease when the same rice cultivars are planted in high temperature area.

5. Conclusions

In this study, we presented the processing steps for monitoring hybrid rice IHS_DAS based on three low-altitude remote sensing data and analyzed the influential factors (rice cultivars, N rates, and air temperature) of RP.

The results indicated that CIred edge could show hybrid rice phenological change accurately without saturation problem, and was more appropriate than NDVI for monitoring IHS_DAS. SGF could eliminate the noise of time series CIred edge without over fitting problem. MC of SGF_CIred edge from three low-altitude remote sensing platforms had good performance in monitoring IHS_DAS with good robustness, as R² varied between 0.82 and 0.95 and RMSE ranged from 2.31 to 3.81.

This study demonstrated that high air temperature might cause a decrease of IHS_DAS when the same rice cultivars were planted in high temperature area. In addition, the interaction of different rice cultivars and N rates on IHS_DAS reached a significant level. The growth process of rice was delayed when more nitrogen fertilizer was applied before IHS_DAS.

Future work should evaluate these results for diverse rice cultivars under different environmental conditions, by combing low-altitude remote sensing data and satellite data. In addition, new methods for reconstructing rice phenological curves will be explored for monitoring RP in the future. We expect that our study will provide useful guidance for breeding rice and fertilizer management with remote sensing technology.