Phenology Effects on Physically Based Estimation of Paddy Rice Canopy Traits from UAV Hyperspectral Imagery

: Radiation transform models such as PROSAIL are widely used for crop canopy reﬂectance simulation and biophysical parameter inversion. The PROSAIL model basically assumes that the canopy is turbid homogenous media with a bare soil background. However, the canopy structure changes when crop growth stages develop, which is more or less a departure from this assumption. In addition, a paddy rice ﬁeld is inundated most of the time with ﬂooded soil background. In this study, ﬁeld-scale paddy rice leaf area index (LAI), leaf cholorphyll content (LCC), and canopy chlorophyll content (CCC) were retrieved from unmanned-aerial-vehicle-based hyperspectral images by the PROSAIL radiation transform model using a lookup table (LUT) strategy, with a special focus on the effects of growth-stage development and soil-background signature selection. Results show that involving ﬂooded soil reﬂectance as background reﬂectance for PROSAIL could improve estimation accuracy. When using a LUT with the ﬂooded soil reﬂectance signature ( LUT flooded ) the coefﬁcients of determination ( R 2 ) between observed and estimation variables are 0.70, 0.11, and 0.79 for LAI, LCC, and CCC, respectively, for the entire growing season (from tillering to heading growth stages), and the corresponding mean absolute errors ( MAEs ) are 21.87%, 16.27%, and 12.52%. For LAI and LCC, high model bias mainly occurred in tillering growth stages. There is an obvious overestimation of LAI and underestimation of LCC for in the tillering growth stage. The estimation accuracy of CCC is relatively consistent from tillering to heading growth stages.


Introduction
Crop growth status variation is related to meteorological conditions, field management, soil property, genotype diversity and other factors [1]. Spatial mapping of crop growth-status variation could benefit site-specific field management [2]. Precise field management is the key to improve crop production and light-utilization efficiency [3]. Remote sensing provides an avenue to crop growth status monitoring [4]. Leaf area index (LAI), leaf chlorophyll content (LCC), and canopy chlorophyll content (CCC) are three key parameters that can not only characterize crop growth status, but are also retrievable from remote-sensing technology [5]. LAI, which is defined as half of the all-sided green leaf area per unit ground area, reflects biochemical and physiological processes of crops [6]. LAI mapping is important for a wide range of agricultural studies, such as stress evaluation, growth-status monitoring, and yield estimation [7]. LCC and CCC are chlorophyll (including chlorophyll a and b) content per unit leaf area and ground area, respectively. Chlorophyll content is an important indicator for assessing crop nitrogen content, photosynthetic capacity, and senescent and environmental stress [8].
Numerous approaches have been proposed to retrieve these three key parameters. These approaches can be generally classified into three classes: empirical statistical, physical, and a hybrid of both [9]. The empirical statistical approach builds linear/nonlinear regression models on a carefully designed vegetation index (VI) or selected features from the raw reflectance spectrum, or builds non-parametric regression models on full or subset bands of the reflectance spectra. The widely used non-parametric regression models include partial least-squares regression (PLSR), support vector machine (SVM), random forest (RF), Gaussian process regression (GPR), and artificial neural networks (ANNs) [10]. The physical approach combines radiation transform models (RTMs) with difference inversion strategies [11]. The two most widely used inversion strategies are interactively numerical optimization and the lookup table (LUT) method [12,13]. The general idea of hybrid approach is training statistical models on the RTM generated LUTs [14,15]. The widely used model in hybrid approach is the ANN [16,17]. Recent studies have shown that the SVM [18] and GRP methods also exhibit promising performance [11,19,20]. The numerical approach is straightforward and can usually yield good accuracy. However, it is criticized that is site-, sensor-, phenology-, and crop-specific [21,22]. For the interactively numerical optimization approach, the inversion result may be trapped in local minima, which cannot guarantee a stable and global optimal solution. Moreover, it is computationally intensive. The LUT and hybrid approaches reduce the computational demands, but need more carefully designed simulated databases and proper training [23,24].
When an RTM is involved, different leaf/canopy parameter configurations may yield identical simulated reflectance signatures, which is the "ill-posed" inversion problem [25]. Limited bands of spectral data could amplify this problem. Additionally, the model and measurement bias could raise the inversion inaccuracy. Several regularization schemes have been proposed to mitigate these problems, such as adding additional white Gaussian noise, considering a priori information, and using spectral subsets [26]. For the LUT inversion scheme, a multiple-solutions strategy can make the retrieved result more robust [1,27,28].
For crops, the widely used RTM is the PROSAIL model [24,29], which is a combination of the leaf-level RTM PROPECT [30,31] and canopy bidirectional reflectance model SAIL [32,33]. The PROSAIL model is based on the turbid medium assumption that the leaves are randomly distributed. The PROSAIL model simulates canopy bidirectional reflectance within 400-2500 nm with a step of 1 nm as a function of parameters related to the leaf optical properties, canopy structure, background soil signature, and sun-view geometry. The PROSAIL model has been evaluated on several crops, such as winter wheat, maize, paddy rice, sugar beets, oilseed rape, and potato [23]. In particular, for paddy rice, [28] retrieved paddy rice CCC from multispectral satellite imagery using a PROSAIL-LUT approach. The CCC is retrieved with R 2 = 0.65 and root mean-square error (RMSE = 45 µg cm −2 ). They also evaluated how the LUT size and number of solutions affect the retrieval accuracy. [34] tried to map paddy rice LAI and CCC from unmannedaerial-vehicle-based (UAV-based) multispectral and hyperspectral images by coupling PROSAIL and Bayesian network models. [19,35] highlighted the importance of adequately characterizing the background situation when the PROSAIL model involved paddy rice LAI estimation.
The PROSAIL model generally performs well in most situations for the aforementioned crops, while some studies have pointed out that PROSAIL could not characterize the canopy reflectance of row-planted crops, especially for early growth stages when the canopy is not fully closed [36][37][38] because of foliage clumping and shadowing effects. Paddy rice is normally row-planted and is different from other crops in that a paddy field is flooded most of time during the growth season. The canopy structure of a paddy changes while the growth stage develops, which places the paddy rice canopy in accord with the assumption of the PROSAIL model changing. How the growth-stage development affects the biochemical or biophysical parameters of paddy rice (such as LAI, LCC, and CCC) has seldom been evaluated. In addition, the inundated soil background makes simulation of the paddy rice canopy reflectance signature more complex. How to design a proper soil-background signature for the PROSAIL model has also not been well discussed.
Paddy rice is one of the three major staple food crops in the world [34]. Accurately monitoring its growth status will benefit guiding its management. The main objective of the present study is to retrieve field-scale paddy rice LAI, LCC, and CCC from UAV-based hyperspectral imagery with a PROSAIL RTM using a LUT scheme, with a special focus on the effects of growth-stage development and soil-background signature selection.

Study Area and Experimental Setup
Experiments were conducted during the 2020 late rice growing season at the rice research field of Guangdong Academy of Agricultural Sciences (23°23 38 N, 113°25 37 E). This region is characterized by the sub-tropical monsoon climate. The average annual precipitation and average temperature are 1700 (mm) and 22.5°C, respectively. The parent material of the soil in the test area is lateritic sandy loam formed by river alluvium. The paddy rice variety Huangruanxiuzhan was used as the test material. The paddy rice was sowed in a nursery bed and transplanted on 9 August 2020.
To acquire necessary reference observations, a paddy rice field with nitrogen treatment was set up and used to conduct field-data collection campaigns (detailed in Figure 1). Nitrogen from 0 to 180 (kg ha −1 ) with 30 (kg ha −1 ) increase were applied to plots J0 to J6 correspondingly 7 days after transplantation (DAT). Nitrogen from 0 to 90 (kg ha −1 ) with 15 (kg ha −1 ) increase were applied to plots F0 to F6 correspondingly at DAT 33. Additionally, 90 (kg ha −1 ) nitrogen were applied to each plots from F0 to F6 at DAT 7. The paddy rice was at tillering and jointing growth stage at DAT 7 and DAT 33, correspondingly. Phosphates and potash fertilizers were applied as basal dressing with P 2 O 5 and K 2 O at 45 and 120 (kg ha −1 ), respectively, for each plot. For both J0 to J6 and F0 to F6, each nitrogen level was comprised of two replicates. Thus, there were 28 plots in total. Each plot was (5 × 4.5 m) in size. The plant density was 400 plants/plot.

Field Data Collection
Four field-data collection campaigns were carried out at DAT 18, 32, 47 and 62. The paddy rice was in tillering, jointing, booting, and heading growth stages in the respective field-data collection campaigns. The growth stage was determined according to the rules described by [39].
A destructive sampling routine was used to collect field LAI, LCC, and CCC data. For each field-data collection campaign, three to five rice plants of each plot were random destructively sampled. For each sample, 30 random SPAD-502 (Minolta Camera Co., Osaka, Japan) readings were recorded from paddy rice leaf. Then, all green leaves were scanned by a portable scanner (PERFECTION V39, Seiko Epson Corp., Nagano, Japan). The SPAD readings were converted to LCC ((mg cm −2 )) by relationship (Equation (1)) proposed by [40] and averaged as the LCC of corresponding plot. The leaf pixels in scanned images were extracted by a threshold method. The LAIs of each plot were determined by Equation (2), where ρ green lea f is the green-leaf pixel ratio of each scanned image, area img = 216 × 297 × 10 −4 m is the maximum scan area of the scanner, and density = 400/(5 × 4.5) m −2 is the planting density. CCC was approximated by multiplying LAI and LCC (Equation (3)). During data collection campaigns at DAT 18 and DAT 32, destructive sample only conducted in plots from J0 to J6 (14 samples in total for each data collection campaign). Plots from F0 to F6 were not sampled because they received same fertilization as J3 before DAT 33. During data collection campaigns at DAT 47, destructive sample conducted in all plots (28 samples in total). During data collection campaigns at DAT 62, destructive sample only conducted in plots from F0 to F6 (14 samples in total).
Canopy reflectance data were acquired by a hyperspectral camera Cubert S185 (Cubert GmbH, Ulm, Germany) mounted on a hexacopter (DJI M600Pro, SZ DJI Technology Co., Shenzhen, Guangzhou, China). The hyperspectral camera has a spectra range of 450-950 nm and a spectral resolution of 4 nm. All UAV flight missions were taken at approximately 11:00 am. The flight height was 30 m, which resulted in 16 cm spatial resolution for hyperspectral images. Figure 2 shows the extent and spatial resolution of the S185 imagery. The hyperspectral images were mosaiced and orthographized for further analysis. The bands beyond 850 nm were dropped because of the low spectral quality [41]. Demonstration of S185 image frame and spatial resolution. Right-hand panel is visualization of S185 at original spatial resolution. Left-hand panel is same view with higher spatial resolution. Red rectangle in left-hand panel is a plot boundary. Table Generation Leaf-level RTMs PROPECT-5B [30] were coupled with a canopy bidirectional reflectance model 4SAIL [33] (referred to as PROSAIL hereafter) to simulate paddy rice canopy reflectance. The PROSPECT model simulates leaf directional-hemispherical reflectance and transmittance between 400 and 2500 nm as a function of six parameters: leaf-structure parameter N, leaf chlorophyll content LCC, leaf carotenoid content Car, leaf brown-pigment content Cbrown, leaf equivalent water thickness Cw, and leaf drymatter content Cm. The 4SAIL model simulates canopy bidirectional reflectance by leaf area index LAI, average leaf angle ALA, hotspot parameter hspot, solar zenith angle tts, observer zenith angle tto, relative azimuth angle psi, background soil reflectance rsoil, and leaf bidirectional reflectance and transmittance, which come from PROSPECT-5B. Table 1 shows the parameters used in the PROSAIL model for paddy rice canopy simulation. N, LCC, Car, LAI and ALA were generated by uniform distribution and CCC by multiplying LAI and LCC. The bounds or fixed values of N, LCC, Car, Cw, Cm, LAI, ALA, and hspot were taken from field-collected data and other studies [19,20,28,34,35]. Cw and Cm were set to fix values because water and dry-matter absorption have only marginally effect on the visible-light and near-infrared reflectance. Cbrown was set to zero because no obvious senescent leaves were observed during field-data collection campaigns.

Inversion Method
Three types of soil background signature were evaluated in this study. First, one group of bare soil reflectances was used as rsoil. A simple multiplicative soil brightness factor α soil ranging from 0 to 1 with step 0.1 was applied to a moist and dry soil spectrum to mimic different soil reflectances caused by soil water content and surface roughness (Equation (4), where R moist , R dry , and R Bare are the moist, dry, mimicked spectrum, respectively). Second, one group of flooded soil reflectances was used as rsoil. A mean spectrum of flooded soil was extracted from the UAV images. Then, a factor β ranging from 0.5 to 3.0 with a step of 0.3 was multiplied with the mean spectrum to mimic differences in flooded soil reflectance. Third, a combination of the two above-mentioned ground reflectances were used as rsoil.
Each pixel is assumed to be linearly composed of a fraction vCover of pure vegetation and (1 − vCover) bare/flooded soil signature. The LAI value will be fixed to LAI × vCover. In this study, the vCover is assumed to have a uniform distribution (ranging from 0.6 to 1 with a step of 0.1) and not dependent on the LAI.
A wavelength-dependent random Gaussian noise of 0.4% observed reflectance standard deviation was added to each band of each synthetic spectrum to simulate the instrument and model noise. Equation (5) shows how the Gaussian noise was added, where R * λ and R λ denote noise added and pure simulated spectrum and λ nm, respectively, sd λ denotes the standard deviation of observed reflectance at λ nm.
Finally, the spectra were resampled to the field-collected spectral resolution using a Gaussian response function with a full width at half-maximum at 4 nm. Considering that three types of background signature were evaluated (Bare, Flooded, and Bare+Flooded), the generated LUTs are referred to as LUT Bare , LUT Flooded , and LUT Bare+Flooed , respectively, hereafter.

Model Inversion and Evaluation
Root-mean-squared error (RMSE) is used as the cost function with which to measure the agreements between measured and simulated spectra. Equation (6) shows the definition of cost function (L). R Measured λi and R LUT λi are the measured and simulated reflectance, respectively, of i th band. n is the total number of bands and L is the cost value. Multiplesolutions regularization was applied. The median of parameters corresponding to the top 100 spectra that yield lost cost values were considered the final solution.
A pixel-by-pixel inversion was applied to the acquired hyperspectral images. Then, plot-wise mean predicted values were extracted by regions of interest generated by plot boundary (demonstrated in Figure 1) with a 0.5-m inner buffer. The plot-wise mean predicted values were compared with field-collected LAI, LCC, and CCC values. The coefficient of determination (R 2 ), RMSE, and mean relative error (MRE) were calculated to evaluate model performance (Equations (7)-(9), where O i and P i are ith observed and predicted data (LAI, LCC or LAI), respectively, and n is the size of data pairs.).   Figure 3 shows the inversion result from plots of mean spectra. The growth-stagespecific goodness-of-fit is depicted in Table 3. For LAI, the MREs are 27.28%, 21.87%, and 23.26%, respectively, and the R 2 values are 0.53, 0.70, and 0.66 of the entire dataset for LUT bare , LUT Flooded , and LUT bare+ f looded , respectively. Although much variance is explained (R 2 ≥ 0.88), the LAI of tillering growth stage (LAI < 2.0) is always overestimated (MRE ≥ 66.32%) regardless of which rsoil values are used. For the other three growth stages, the growth-stage-specific MRE is generally less than 13.10%, and the growth-stagespecific R 2 is greater than 0.61. For LCC, the MREs are 20.35%, 16.27%, and 17.33%, respectively, and the R 2 values are 0.16, 0.11, and 0.05 of the entire dataset for LUT bare , LUT Flooded , and LUT bare+ f looded , respectively. The LCC of the tillering growth stage (LCC ≥ 40 µg cm −2 ) is generally underestimated (R 2 ≤ 0. 46 Table 3. Growth-stage-specific goodness-of-fit of results of Figure 3. Without considering tillering growth stage, the growth-stage-specific inversion accuracy increases most of the time. At jointing growth stage, using flooded or bare plus flooded soil signatures instead of bare soil signatures in PROSAIL, the MRE decreased from 13.10 % to 12.45% and 12.93%, respectively, for LAI; MRE decreased from 8.78% to 6.01% and 5.97%, respectively, for LCC; and MRE decreased from 15.66% to 11.93% and 12.07%, respectively, for CCC. These results demonstrate that using flooded soil signatures instead of bare soil signatures could increase inversion accuracy at jointing growth stage. For booting and heading growth, the bare-soil-signature-based inversion could obtain the highest R 2 or lowest MRE for LAI, LCC, or CCC. These results demonstrate that there is no obvious advantage of using flooded or bare plus flooded soil signatures at booting and heading growth stages. Figure 4 shows the pixel-wise inversion results with flooded soil signatures for LAI, LCC, and CCC for each tillering, jointing, booting, and heading growth stages. The withinpaddy-rice treatment plot variance is low and the difference between different-paddy-rice treatment plots is obvious. The spatio-temporal patterns of estimated LAI and CCC are consistent with the phenology development of paddy rice.

Discussion
The soil background signature parameter rsoil is an important parameter in the PROSAIL model. Paddy rice is different from other crops in that it is inundated most of the time. This raises the question of whether using flooded soil signatures or a combination of bare and flooded soil signatures instead of bare soil signatures as rsoil could increase model accuracy. Besides, with paddy rice growth stage developing, the canopy structure is changing, whether or not the advantage of either rsoil value has not been widely discussed.
When multiple-growth-stage data were combined, the MRE shows a pattern that LUT f looded < LUT bare+ f looded < LUT bare for paddy rice LAI, LCC, and CCC. This means that using bare soil reflectance signatures is always not an optimal choice. However, examining growth-stage-specific inversion accuracy, the advantage of LUT f looded and LUT bare+ f looded are not obvious at late growth stages (booting and heading growth stages). This is because the canopy is almost closed at this time, and thus the influence of background reflectance signature is weakened.
The SAIL model assumes that the canopy is homogenous, and thus the PROSAIL model is expected to perform better when the crop canopy is closed than when the crop is at early growth stages in which the canopy is not completely closed. For paddy rice, the situation is more complex because at early growth stages paddy rice is inundated most of the time. Several studies have reported the deficiency of the PROSAIL model at early growth stages [34,36,42,43]. It was concluded that the PROSAIL model, which is designed for homogenous canopies, may not be able to characterize the non-homogenous canopy structure and shadow effect caused by row-planted crops, especially at early growth stages. The results of the present study show that, at early growth stage (tillering), although much of the LAI variance (R 2 = 0.92) is explained, an obvious overestimation of LAI (MRE = 66.32% and RMSE = 1.11 with observed mean at 1.50) occurred. Similar results can be found in [44] that the PROSAIL-LUT method overestimated the LAI of grassland when LAI < 2.0. Several studies have shown an underestimation of LAI when LAI > 4.0 or > 4.5 [36,42,44]. However, these phenomena are not obvious in the present results. An obvious underestimation of LCC occurred at tillering growth stage compared to other growth stages. However, the CCC was well retrieved and no obvious over-or underestimation occurred. This supports the results of previous studies that, by converting LAI and LCC into CCC, more regularization of the LUT approach could be obtained, yielding more accurate results.
Previous studies showed the covariance between LAI and ALA makes the retrieval of either difficult [12,27,43]. Both decreasing LAI and ALA may have the same effects on the reflectance spectrum [27]. Furthermore, [45] argued that the ALA cannot be interpreted as physical leaf inclinations, and suggested treating it as a free calibration parameter. The present study is based on a loose assumption regarding the ALA, namely that a wide range of ALA (20-85 • ) was used to generate the LUT. With a multiple-solutions regularization strategy, this setup yielded a good multiple-growth-stage LAI retrieval accuracy and performed well from jointing to heading growth stage.
When multiple-growth-stage data are considered together, the LCC estimation accuracy is not satisfied with R 2 ≤ 0.16 and MRE ≥ 17.33%. However, when considering growth stage specifically, the LCC estimation accuracy is moderate from jointing to heading growth stages, with 0.40 ≤ R 2 ≤ 0.81 and 4.20% ≤ MRE ≤ 19.63% for LUT f looded . LCC estimation accuracy from canopy reflectance has been generally poor in previous studies [27,28,42,44,46,47]. This indicates the poor relationship between canopy reflectance and leaf properties due to poor signal propagation from leaf to canopy scale [28]. Furthermore, the lower R 2 value may be caused by the small LCC variance in the datasets [44,48].

Conclusions
UAV-based hyperspectral images combined with the physical model PROSAIL were used to retrieve paddy leaf area index (LAI), leaf chlorophyll content (LCC), and canopy chlorophyll content (CCC) by a lookup table (LUT) approach with a special focus on the effects on growth-stage development and soil background signature selection. Datacollection campaigns were carried out at each of the following growth stages: tillering, jointing, booting, and heading.
The results suggest that, considering using the flooded soil reflectance signature instead of solely using the bare soil reflectance as the soil background reflectance, could improve retrieval accuracy. When using a LUT with the flooded soil reflectance signature (LUT f looded ), the retrieval accuracies are R 2 = 0.70, 0.11, and 0.79, and MAE = 21.87%, 16.27%, and 12.52% for LAI, LCC, and CCC, respectively, for the entire growing season. When considering growth-stage-specific retrieval accuracy, an obvious overestimation of LAI is apparent, as is an underestimation of LCC for the tillering growth stage. For CCC, the retrieval accuracy lacks the obvious bias of the tillering growth stage with R 2 = 0.90 and MRE = 16.07% (LUT f looded ). For the other growth stages (jointing, booting, and heading), reasonable retrieval accuracies were acquired for all three parameters, i.e., LAI, LCC, and CCC.
These findings could benefit the parameterization of PROSAIL model for paddy rice by designing proper soil background signature, and also give a better understanding of the model performance variance caused by growth stage difference when the LAI, LCC and CCC are inverted by a PROSAIL-LUT approach.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data sharing is not application to this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: